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 Upsilon(1,2,3S) production at 7TeV
 10  class ATLAS_2013_I1204994 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2013_I1204994);
 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==553 || Cuts::pid==100553 || Cuts::pid==200553), "UFS");
 24
 25      // branching ratios
 26      _br = {0.0248,0.0193,0.0218};
 27
 28      // total cross section
 29      book(_h_total,1,1,1);
 30
 31      // double differential
 32      const vector<double> ybins = {0.,1.2,2.25};
 33      for(unsigned int ix=0;ix<3;++ix) {
 34      	book(_h_Upsilon_fid_y[ix],5,1,ix+1);
 35      	book(_h_Upsilon_y    [ix],9,1,ix+1);
 36      	book(_h_Upsilon_y_r  [ix],"TMP/Ups"+toString(ix),refData(12,1,1));
 37        book(_h_Upsilon_fid[ix], ybins);
 38        book(_h_Upsilon[ix],     ybins);
 39        book(_h_Upsilon_r[ix],   ybins);
 40        for (size_t iy=1; iy < ybins.size(); ++iy) {
 41      	  book(_h_Upsilon_fid[ix]->bin(iy), ix+2,1,iy);
 42      	  book(_h_Upsilon[ix]->bin(iy), ix+6,1,iy);
 43      	  book(_h_Upsilon_r[ix]->bin(iy), "TMP/Ups"+toString(ix)+"_"+toString(iy), refData(11,1,iy));
 44      	}
 45      }
 46    }
 47
 48    void findDecayProducts(const Particle & mother, unsigned int & nstable,
 49                           Particles & mup, Particles & mum, unsigned int & ngamma) {
 50      for (const Particle & p : mother.children()) {
 51        int id = p.pid();
 52        if (id == PID::MUON ) {
 53          ++nstable;
 54          mum.push_back(p);
 55        }
 56        else if (id == PID::ANTIMUON) {
 57          ++nstable;
 58          mup.push_back(p);
 59        }
 60        else if (id == PID::PI0 || id == PID::K0S || id == PID::K0L ) {
 61          ++nstable;
 62        }
 63        else if (id == PID::GAMMA && p.children().empty() ) {
 64          ++ngamma;
 65          ++nstable;
 66        }
 67        else if ( !p.children().empty() ) {
 68          findDecayProducts(p, nstable, mup, mum,ngamma);
 69        }
 70        else
 71          ++nstable;
 72      }
 73    }
 74
 75    /// Perform the per-event analysis
 76    void analyze(const Event& event) {
 77
 78      // Final state of unstable particles to get particle spectra
 79      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
 80
 81      for (const Particle& p : ufs.particles()) {
 82      	// check if in fiducal region
 83      	unsigned int nstable=0, ngamma=0;
 84      	Particles mup,mum;
 85      	findDecayProducts(p,nstable,mup,mum,ngamma);
 86      	bool fiducial = false;
 87      	if (mup.size()==1 && mum.size()==1 && nstable==ngamma+2) {
 88      	  fiducial = true;
 89      	  if(mup[0].perp()<4. || mup[0].abseta()>2.4) fiducial = false;
 90      	  if(mum[0].perp()<4. || mum[0].abseta()>2.4) fiducial = false;
 91      	}
 92      	// pT and rapidity
 93      	double absrap = p.absrap();
 94      	double xp = p.perp();
 95      	// type of upsilon
 96      	unsigned int iups=p.pid()/100000;
 97      	if (fiducial) {
 98      	  if (xp<70.) _h_Upsilon_fid_y[iups]->fill(absrap);
 99      	  _h_Upsilon_fid[iups]->fill(absrap,xp);
100      	}
101      	if (xp<70.) {
102          _h_Upsilon_y[iups]->fill(absrap);
103          _h_Upsilon_y_r[iups]->fill(absrap);
104        }
105      	_h_Upsilon[iups]->fill(absrap,xp);
106      	_h_Upsilon_r[iups]->fill(absrap,xp);
107      	if (absrap<2.25 && xp<70) _h_total->fill(_edges[iups],_br[iups]);
108      }
109    }
110
111
112    /// Normalise histograms etc., after the run
113    void finalize() {
114
115      scale(_h_total,crossSection()/nanobarn/sumOfWeights());
116
117      for (unsigned int ix=0; ix<3; ++ix) {
118      	scale(_h_Upsilon_fid_y[ix], 0.5*crossSection()/picobarn/sumOfWeights());
119      	scale(_h_Upsilon_y[ix],     0.5*_br[ix]*crossSection()/picobarn/sumOfWeights());
120      	scale(_h_Upsilon_y_r[ix],   0.5*_br[ix]*crossSection()/picobarn/sumOfWeights());
121
122      	scale(_h_Upsilon_fid[ix], 0.5*crossSection()/femtobarn/sumOfWeights());
123      	scale(_h_Upsilon[ix],     0.5*_br[ix]*crossSection()/femtobarn/sumOfWeights());
124      	scale(_h_Upsilon_r[ix],   0.5*_br[ix]*crossSection()/femtobarn/sumOfWeights());
125        divByGroupWidth(_h_Upsilon_fid[ix]);
126        divByGroupWidth(_h_Upsilon[ix]);
127        divByGroupWidth(_h_Upsilon_r[ix]);
128      }
129
130      // ratios
131      for (size_t iy=1; iy<_h_Upsilon[0]->numBins()+1; ++iy) {
132      	Estimate1DPtr tmp;
133      	book(tmp,10,1,iy);
134      	divide(_h_Upsilon_r[1]->bin(iy), _h_Upsilon_r[0]->bin(iy), tmp);
135      	book(tmp,11,1,iy);
136      	divide(_h_Upsilon_r[2]->bin(iy), _h_Upsilon_r[0]->bin(iy), tmp);
137      }
138      for (unsigned int iy=0; iy<2; ++iy) {
139      	Estimate1DPtr tmp;
140      	book(tmp,12,1,iy+1);
141      	divide(_h_Upsilon_y_r[iy+1], _h_Upsilon_y_r[0], tmp);
142      }
143    }
144    /// @}
145
146
147    /// @name Histograms
148    /// @{
149    BinnedHistoPtr<string> _h_total;
150    Histo1DPtr _h_Upsilon_fid_y[3],_h_Upsilon_y[3],_h_Upsilon_y_r[3];
151    Histo1DGroupPtr _h_Upsilon[3],_h_Upsilon_fid[3],_h_Upsilon_r[3];
152    vector<double> _br;
153    vector<string> _edges = { "P P --> UPSILON(1S) < MU+ MU- > X",
154                              "P P --> UPSILON(2S) < MU+ MU- > X",
155                              "P P --> UPSILON(3S) < MU+ MU- > X" };
156    /// @}
157
158
159  };
160
161
162  RIVET_DECLARE_PLUGIN(ATLAS_2013_I1204994);
163
164}