1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4
  5
  6namespace Rivet {
  7
  8
  9  /// @brief D meson production at 5.02 and 7 TeV
 10  class ALICE_2019_I1716440 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(ALICE_2019_I1716440);
 15
 16
 17    /// @name Analysis methods
 18    /// @{
 19
 20    /// Book histograms and initialise projections before the run
 21    void init() {
 22      // projections
 23      declare(UnstableParticles(Cuts::abspid==411 || Cuts::abspid==421 || Cuts::abspid==431 || Cuts::abspid==413), "UFS");
 24      // histograms
 25      if (isCompatibleWithSqrtS(5020.0)) {
 26        for (unsigned int ix=0; ix<4; ++ix) {
 27          book(_h_prompt[ix],1+ix,1,1);
 28          book(_h_ratio_num[ix],"TMP/h_ratio_num_"+toString(ix+1), refData(6+ix,1,1));
 29          book(_h_ratio_den[ix],"TMP/h_ratio_den_"+toString(ix+1), refData(6+ix,1,1));
 30        }
 31        book(_h_incl,   5, 1, 1);
 32        book(_h_total, 18, 1, 1);
 33      }
 34      else if (isCompatibleWithSqrtS(7000.0)){
 35        for (unsigned int ix=0;ix<4;++ix) {
 36          book(_h_ratio_num[ix],"TMP/h_ratio_num_"+toString(ix+1), refData(10+ix,1,1));
 37          book(_h_ratio_den[ix],"TMP/h_ratio_den_"+toString(ix+1), refData(10+ix,1,1));
 38        }
 39      }
 40      else {
 41        throw UserError("Centre-of-mass energy of the given input is neither 5020 nor 7000 GeV.");
 42      }
 43      for (unsigned int ix=0; ix<4; ++ix) {
 44        book(_h_energy[ix],"h_energy_"+toString(ix+1), refData(14+ix,1,1));
 45      }
 46    }
 47
 48
 49    /// Perform the per-event analysis
 50    void analyze(const Event& event) {
 51      // Final state of unstable particles to get particle spectra
 52      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
 53      for (const Particle& p : ufs.particles()) {
 54        // no mixing and |y|<0.5
 55        if(p.children().size()==1 || p.absrap()>0.5) continue;
 56        unsigned int imeson=0;
 57        if      (p.abspid()==411) {
 58          imeson=1;
 59        }
 60        else if (p.abspid()==413) {
 61          imeson=2;
 62        }
 63        else if (p.abspid()==431) {
 64          imeson=3;
 65        }
 66        const double pT=p.perp();
 67        if (p.fromBottom()) {
 68          if (imeson==0 && _h_incl) _h_incl->fill(pT);
 69          continue;
 70        }
 71        // prompt at 5.02 TeV
 72        if (_h_prompt[imeson]) {
 73          _h_prompt[imeson]->fill(pT);
 74          _h_total->fill(_edges[imeson]);
 75        }
 76        if (imeson==0) {
 77          _h_ratio_den[0]->fill(pT);
 78          _h_ratio_den[1]->fill(pT);
 79          _h_ratio_den[2]->fill(pT);
 80          if (_h_incl) _h_incl->fill(pT);
 81        }
 82        else if (imeson==1) {
 83          _h_ratio_num[0]->fill(pT);
 84          _h_ratio_den[3]->fill(pT);
 85        }
 86        else if (imeson==2) {
 87          _h_ratio_num[1]->fill(pT);
 88        }
 89        else if (imeson==3) {
 90          _h_ratio_num[2]->fill(pT);
 91          _h_ratio_num[3]->fill(pT);
 92        }
 93        _h_energy[imeson]->fill(pT);
 94      }
 95    }
 96
 97
 98    /// Normalise histograms etc., after the run
 99    void finalize() {
100      const double factor = crossSection()/microbarn/sumOfWeights();
101      if (_h_prompt[0]) {
102        for(unsigned int ix=0;ix<4;++ix) {
103          scale(_h_prompt[ix],factor);
104        }
105        scale(_h_total,factor);
106        scale(_h_incl,factor);
107      }
108      int ioff = 0;
109      if (isCompatibleWithSqrtS(7000.0)) {
110        ioff = 1;
111      }
112      for (unsigned int ix=0;ix<4;++ix) {
113        Estimate1DPtr tmp;
114        scale(_h_energy   [ix],factor);
115        scale(_h_ratio_num[ix],factor);
116        scale(_h_ratio_den[ix],factor);
117        book(tmp,6+4*ioff+ix,1,1);
118        divide(_h_ratio_num[ix],_h_ratio_den[ix],tmp);
119      }
120    }
121
122    /// @}
123
124
125    /// @name Histograms
126    /// @{
127    Histo1DPtr _h_prompt[4],_h_incl;
128    BinnedHistoPtr<string> _h_total;
129    Histo1DPtr _h_ratio_num[4],_h_ratio_den[4];
130    Histo1DPtr _h_energy[4];
131    vector<string> _edges = { "P P --> D0 (Q=PROMPT) X",
132                              "P P --> D+ (Q=PROMPT) X",
133                              "P P --> D*+ (Q=PROMPT) X",
134                              "P P --> D/s+ (Q=PROMPT) X" };
135    /// @}
136
137
138  };
139
140
141  RIVET_DECLARE_PLUGIN(ALICE_2019_I1716440);
142
143}