1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// @brief J/psi and psi(2S) at 7 TeV
  9  class CMS_2012_I944755 : public Analysis {
 10  public:
 11
 12    /// Constructor
 13    RIVET_DEFAULT_ANALYSIS_CTOR(CMS_2012_I944755);
 14
 15
 16    /// @name Analysis methods
 17    /// @{
 18
 19    /// Book histograms and initialise projections before the run
 20    void init() {
 21      declare(UnstableParticles(), "UFS");
 22      // J/psi/psi(2s) histos
 23      for (unsigned int ipsi=0; ipsi<2; ++ipsi) {
 24        for (unsigned int iprompt=0; iprompt<2; ++iprompt) {
 25          book(_h_pT[ipsi][iprompt],
 26               "TMP/pT_"+toString(ipsi)+"_"+toString(iprompt),
 27               refData(20, 1, 1+iprompt));
 28        }
 29        const vector<double> raps = ipsi==0 ? vector<double>({0.,0.9,1.2,1.6,2.1,2.4}) : vector<double>({0.,1.2,1.6,2.4});
 30        for (unsigned int iprompt=0; iprompt<3; ++iprompt) {
 31          book(_h_psi[ipsi][iprompt], raps);
 32          if (ipsi==1 && iprompt<2) {
 33            book(_h_ratio[0][iprompt], raps);
 34            book(_h_ratio[1][iprompt], raps);
 35          }
 36          for (unsigned int iy=1; iy<raps.size(); ++iy) {
 37            if (iprompt<2) {
 38              book(_h_psi[ipsi][iprompt]->bin(iy), iy+5*ipsi, 1, 1+iprompt);
 39            }
 40            else {
 41              book(_h_psi[ipsi][iprompt]->bin(iy),
 42                   "TMP/psi_"+toString(ipsi)+"_"+toString(iy),
 43                   refData(iy+5*ipsi,1,1));
 44            }
 45            if (iprompt<2 && ipsi==1) {
 46              book(_h_ratio[0][iprompt]->bin(iy),
 47                   "TMP/ratio_0_"+toString(iprompt)+"_"+toString(iy),
 48                   refData(16+iy,1,1+iprompt));
 49              book(_h_ratio[1][iprompt]->bin(iy),
 50                   "TMP/ratio_1_"+toString(iprompt)+"_"+toString(iy),
 51                   refData(16+iy,1,1+iprompt));
 52            }
 53          }
 54        }
 55      }
 56    }
 57
 58
 59    /// Perform the per-event analysis
 60    void analyze(const Event& event) {
 61      // Final state of unstable particles to get particle spectra
 62      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
 63      // loop over onium states
 64      for (const Particle& p : ufs.particles(Cuts::pid==443 || Cuts::pid==100443)) {
 65        // cuts on rapidity
 66        const double y = p.absrap();
 67        if (y>2.4) continue;
 68        const double pT = p.perp();
 69        // prompt
 70        unsigned int iprompt = p.fromBottom();
 71        unsigned int ipsi=p.pid()/100000;
 72        _h_psi  [ipsi][iprompt]->fill(y, pT);
 73        _h_psi  [ipsi][   2   ]->fill(y, pT);
 74        _h_ratio[ipsi][iprompt]->fill(y, pT);
 75        _h_pT[ipsi][iprompt]->fill(pT);
 76      }
 77    }
 78
 79
 80    /// Normalise histograms etc., after the run
 81    void finalize() {
 82      // 0.5 due folded rapidity
 83      const double factor = 0.5*crossSection() / nanobarn/ sumOfWeights();
 84      // branching ratios to muons
 85      const vector<double> brs={0.05961,0.00793};
 86      for (unsigned int ipsi=0; ipsi<2; ++ipsi) {
 87        for(unsigned int iprompt=0; iprompt<2; ++iprompt) {
 88          scale(_h_psi[ipsi][iprompt], factor*brs[ipsi]);
 89          divByGroupWidth(_h_psi[ipsi][iprompt]);
 90          scale(_h_ratio[ipsi][iprompt],factor * brs[ipsi]);
 91          divByGroupWidth(_h_ratio[ipsi][iprompt]);
 92          scale(_h_pT[ipsi][iprompt], factor*brs[ipsi]);
 93        }
 94        scale(_h_psi[ipsi][2], factor*brs[ipsi]);
 95        divByGroupWidth(_h_psi[ipsi][2]);
 96      }
 97      // ratios psi(2S) over J/psi
 98      for (unsigned int iprompt=0; iprompt<2; ++iprompt) {
 99        for (unsigned int iy=0; iy<_h_ratio[0][iprompt]->numBins(); ++iy) {
100          Estimate1DPtr tmp;
101          book(tmp, 17+iy, 1, 1+iprompt);
102          divide(_h_ratio[1][iprompt]->bin(1+iy), _h_ratio[0][iprompt]->bin(1+iy), tmp);
103        }
104        Estimate1DPtr tmp;
105        book(tmp, 20, 1, 1+iprompt);
106        divide(_h_pT[1][iprompt], _h_pT[0][iprompt],tmp);
107      }
108      // non-prompt fraction
109      for (unsigned int ipsi=0; ipsi<2; ++ipsi) {
110        for (unsigned int iy=0; iy<_h_psi[ipsi][1]->numBins(); ++iy) {
111          Estimate1DPtr tmp;
112          book(tmp, 21+iy+5*ipsi, 1, 1);
113          efficiency(_h_psi[ipsi][1]->bin(1+iy), _h_psi[ipsi][2]->bin(1+iy), tmp);
114        }
115      }
116    }
117
118    /// @}
119
120
121    /// @name Histograms
122    /// @{
123    Histo1DGroupPtr _h_psi[2][3];
124    Histo1DGroupPtr _h_ratio[2][2];
125    Histo1DPtr _h_pT[2][2];
126    /// @}
127
128
129  };
130
131
132  RIVET_DECLARE_PLUGIN(CMS_2012_I944755);
133
134}