1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4#include "Rivet/Tools/HistoGroup.hh"
  5
  6namespace Rivet {
  7
  8
  9  /// @brief ATLAS J/psi psi2s at 7 and 8 TeV
 10  class ATLAS_2016_I1409298 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2016_I1409298);
 15
 16
 17    /// @name Analysis methods
 18    /// @{
 19
 20    /// Book histograms and initialise projections before the run
 21    void init() {
 22
 23      declare(UnstableParticles(Cuts::pid==443 or Cuts::pid==100443), "UFS");
 24
 25      unsigned int iloc=0;
 26      if (isCompatibleWithSqrtS(7000)) {
 27        iloc = 1;
 28      }
 29      else if (isCompatibleWithSqrtS(8000)) {
 30        iloc = 2;
 31      }
 32      else {
 33        throw UserError("Centre-of-mass energy of the given input is neither 7 or 8 TeV.");
 34      }
 35      // binning in y
 36      const vector<double> yedges = {0.,0.25,0.5,0.75,1.0,1.25,1.5,1.75,2.};
 37      // book histos
 38      for (unsigned int ix=0; ix<3; ++ix) {
 39        book(_h_JPsi[ix], yedges);
 40        if(ix<2) book(_h_JPsi[ix+3], yedges);
 41        book(_h_psi2S[ix], yedges);
 42        for (size_t iy=1; iy < yedges.size(); ++iy) {
 43          if (ix == 2) {
 44            // total no for ratios etc
 45            book(_h_JPsi[ix]->bin(iy),  "TMP/JPsi_2_"+toString(iy),  refData(  iloc,1,iy));
 46            book(_h_psi2S[ix]->bin(iy), "TMP/psi2S_2_"+toString(iy), refData(4+iloc,1,iy));
 47            continue;
 48          }
 49          // prompt and non-prompt Jpsi
 50          book(_h_JPsi[ix]->bin(iy), iloc+2*ix,1,iy);
 51          // prompt and non-prompt psi(2S)
 52          book(_h_psi2S[ix]->bin(iy), 4+iloc+2*ix,1,iy);
 53          // extra Jpsi for ratios with psi2s
 54          if(ix<2) book(_h_JPsi[ix+3]->bin(iy), "TMP/JPsi_"+toString(ix+3)+"_"+toString(iy), refData(4+iloc+2*ix,1,iy));
 55        }
 56      }
 57    }
 58
 59
 60    /// Perform the per-event analysis
 61    void analyze(const Event& event) {
 62      
 63      // Final state of unstable particles to get particle spectra
 64      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
 65
 66      for (const Particle& p : ufs.particles()) {
 67        // prompt/non-prompt
 68        bool nonPrompt = p.fromBottom();
 69        const double absrap = p.absrap();
 70        const double xp = p.perp();
 71        if (p.pid()==443) {
 72          _h_JPsi[nonPrompt  ]->fill(absrap,xp);
 73          _h_JPsi[2          ]->fill(absrap,xp);
 74          _h_JPsi[nonPrompt+3]->fill(absrap,xp);
 75        }
 76        else {
 77          _h_psi2S[nonPrompt]->fill(absrap,xp);
 78          _h_psi2S[2        ]->fill(absrap,xp);
 79        }
 80      }
 81    }
 82
 83
 84    /// Normalise histograms etc., after the run
 85    void finalize() {
 86      // 1/2 due rapidity folding +/-
 87      const double factor = 0.5*crossSection()/nanobarn/sumOfWeights();
 88      // br to muons PDG 2021 (psi2s is e+e- due large errors on mu+mu-)
 89      const vector<double> br = {0.05961,0.00793};
 90      // scale histos
 91      for (unsigned int ix=0; ix<5; ++ix) {
 92        scale(_h_JPsi [ix], factor*br[0]);
 93        divByGroupWidth(_h_JPsi[ix]);
 94        if(ix>2) continue;
 95        scale(_h_psi2S[ix], factor*br[1]);
 96        divByGroupWidth(_h_psi2S[ix]);
 97      }
 98      // ratios, first find CMS energy
 99      unsigned int iloc=0;
100      if (isCompatibleWithSqrtS(7000)) {
101        iloc = 1;
102      }
103      else if  (isCompatibleWithSqrtS(8000)) {
104        iloc = 2;
105      }
106      else {
107        throw UserError("Centre-of-mass energy of the given input is neither 7 or 8 TeV.");
108      }
109
110      for (unsigned int iy=1; iy<=_h_JPsi[0]->numBins(); ++iy) {
111        // non-prompt J/psi percentage
112        Estimate1DPtr tmp;
113        book(tmp,8+iloc,1,iy);
114        efficiency(_h_JPsi[1]->bin(iy),_h_JPsi[2]->bin(iy),tmp);
115        tmp->scale(100.);
116        // non-prompt psi2S percentage
117        book(tmp,10+iloc,1,iy);
118        efficiency(_h_psi2S[1]->bin(iy),_h_psi2S[2]->bin(iy),tmp);
119        tmp->scale(100.);
120        // prompt psi(2s)/J/psi percentage
121        book(tmp,12+iloc,1,iy);
122        divide(_h_psi2S[0]->bin(iy),_h_JPsi[3]->bin(iy),tmp);
123        tmp->scale(100.);
124        // non-prompt psi(2s)/J/psi percentage
125        book(tmp,14+iloc,1,iy);
126        divide(_h_psi2S[1]->bin(iy),_h_JPsi[4]->bin(iy),tmp);
127        tmp->scale(100.);
128      }
129    }
130
131    /// @}
132
133
134    /// @name Histograms
135    /// @{
136    Histo1DGroupPtr _h_JPsi[5],_h_psi2S[3];
137    /// @}
138
139
140  };
141
142
143  RIVET_DECLARE_PLUGIN(ATLAS_2016_I1409298);
144
145}