1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// @brief Upsilon(1,2,3S) at 7 TeV
  9  class CMS_2013_I1225274 : public Analysis {
 10  public:
 11
 12    /// Constructor
 13    RIVET_DEFAULT_ANALYSIS_CTOR(CMS_2013_I1225274);
 14
 15    /// @name Analysis methods
 16    ///@{
 17
 18    /// Book histograms and initialise projections before the run
 19    void init() {
 20      // projection
 21      declare(UnstableParticles(), "UFS");
 22      // histograms
 23      for (unsigned int ix=0; ix<3; ++ix) {
 24        if (ix<2)  book(_h_total[ix],1,1,1+ix);
 25        for (unsigned int iy=0; iy<2; ++iy) {
 26          book(_h_pT_y[ix][iy], {0.,0.4,0.8,1.2,1.6,2.0,2.4});
 27          for(unsigned int iz=0; iz<6; ++iz) {
 28            book(_h_pT_y[ix][iy]->bin(iz+1), 5+3*iz+ix, 1, iy+1);
 29          }
 30          book(_h_pT[ix][iy], 2 +ix, 1, iy+1);
 31          book(_h_y [ix][iy], 23+ix, 1, iy+1);
 32          book(_h_r [ix][iy], "TMP/h_r_"+toString(ix)+"_"+toString(iy), refData(26,1,iy+1));
 33        }
 34        book(_h_pT_acc[ix], 29+ix, 1, 1);
 35      }
 36    }
 37
 38    void findChildren(const Particle& p, Particles& mum, Particles& mup, unsigned int& nstable) {
 39      for (const Particle & child : p.children()) {
 40        if (child.pid()==PID::MUON) {
 41          mum.push_back(child);
 42          ++nstable;
 43        }
 44        else if (child.pid()==PID::ANTIMUON) {
 45          mup.push_back(child);
 46          ++nstable;
 47        }
 48        else if (child.pid()==PID::PHOTON) {
 49          continue;
 50        }
 51        else if (child.children().empty()) {
 52          ++nstable;
 53        }
 54        else {
 55          findChildren(child,mum,mup,nstable);
 56        }
 57      }
 58    }
 59
 60    // from eqn 1 of paper
 61    bool acceptMuon(const Particle & p) const {
 62      const double abseta = p.abseta();
 63      const double xp = p.perp();
 64      if (abseta<0.8)      return xp>3.75;
 65      else if (abseta<1.6) return xp>3.5;
 66      else if (abseta<2.4) return xp>3.0;
 67      else return false;
 68    }
 69
 70    /// Perform the per-event analysis
 71    void analyze(const Event& event) {
 72      if (_edges.empty()) _edges = _h_total[0]->xEdges();
 73      // Final state of unstable particles to get particle spectra
 74      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
 75
 76      for (const Particle& p : ufs.particles(Cuts::pid==553 || Cuts::pid==100553 || Cuts::pid==200553)) {
 77        const double absrap = p.absrap();
 78        // rapidity cut
 79        if (absrap>2.4) continue;
 80        double xp = p.perp();
 81        unsigned int iloc=0;
 82        if      (p.pid()==   553) iloc=0;
 83        else if (p.pid()==100553) iloc=1;
 84        else if (p.pid()==200553) iloc=2;
 85        // acceptance corrected only hist
 86        if (absrap<1.2) {
 87          _h_pT_acc[iloc]->fill(xp);
 88        }
 89        // check if children muons and within acceptance
 90        unsigned int imin=1;
 91        // find the children
 92        Particles mum,mup;
 93        unsigned int nstable(0);
 94        findChildren(p,mum,mup,nstable);
 95        if (mup.size()==1 && mup.size()==1 && nstable==2) {
 96          if (acceptMuon(mup[0]) && acceptMuon(mum[0])) imin=0;
 97        }
 98        // fill the histos
 99        for (unsigned int ix=imin;ix<2;++ix) {
100          _h_pT_y[iloc][ix] ->fill(absrap,xp);
101          _h_pT  [iloc][ix]->fill(xp);
102          _h_r   [iloc][ix]->fill(xp);
103          if (xp<50.) _h_y[iloc][ix]->fill(absrap);
104          if (ix==0) {
105            _h_total[ix]->fill(_edges[iloc]);
106          }
107          else {
108            _h_total[ix]->fill(_edges[iloc], _br[iloc]);
109          }
110        }
111      }
112    }
113
114
115    /// Normalise histograms etc., after the run
116    void finalize() {
117      double factor = crossSection() / nanobarn/ sumOfWeights();
118      for (unsigned int ix=0; ix<2; ++ix) {
119        // total cross sections, just the factor
120        scale(_h_total[ix], factor);
121        for (unsigned int iy=0; iy<3; ++iy) {
122          double factor2 = factor;
123          if (ix==1) factor2*=_br[iy];
124          // pT integrated over y, just the factor
125          scale(_h_pT[iy][ix], factor2);
126          scale(_h_r[iy][ix], factor2);
127          if (ix==1) scale(_h_pT_acc[iy],factor2);
128          // not integrated over y, alsso undo y +/- folding
129          scale(_h_y[iy][ix], 0.5*factor2);
130          scale(_h_pT_y[iy][ix], 0.5*factor2);
131          divByGroupWidth(_h_pT_y[iy][ix]);
132        }
133        // ratios
134        Estimate1DPtr tmp;
135        // ups 3/ ups1
136        book(tmp, 26, 1, ix+1);
137        divide(_h_r[2][ix], _h_r[0][ix], tmp);
138        // ups 2/ ups1
139        book(tmp, 27, 1, ix+1);
140        divide(_h_r[1][ix], _h_r[0][ix], tmp);
141        // ups 3/ ups2
142        book(tmp, 28, 1, ix+1);
143        divide(_h_r[2][ix], _h_r[1][ix], tmp);
144      }
145    }
146
147    ///@}
148
149
150    /// @name Histograms
151    ///@{
152    Histo1DGroupPtr _h_pT_y[3][2];
153    BinnedHistoPtr<string> _h_total[2];
154    Histo1DPtr _h_pT[3][2],_h_pT_acc[3],_h_y[3][2],_h_r[3][2];
155    vector<double> _br{0.0248,0.0193,0.0218};
156    vector<string> _edges;
157    ///@}
158
159
160  };
161
162
163  RIVET_DECLARE_PLUGIN(CMS_2013_I1225274);
164
165}