1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/Beam.hh"
  4#include "Rivet/Projections/UnstableParticles.hh"
  5
  6namespace Rivet {
  7
  8
  9  /// @brief Upsilon polarization at 7 TeV
 10  class CMS_2013_I1185414 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(CMS_2013_I1185414);
 15
 16
 17    /// @name Analysis methods
 18    /// @{
 19
 20    /// Book histograms and initialise projections before the run
 21    void init() {
 22      // projections
 23      declare(Beam(), "Beams");
 24      declare(UnstableParticles(), "UFS");
 25      // book profile hists for the moments
 26      // loop upsilon states
 27      for (unsigned int iups=0; iups<3; ++iups) {
 28        // rapidity intervals
 29        for (unsigned int iy=0; iy<2; ++iy) {
 30          // frame defns
 31          for (unsigned int iframe=0; iframe<3; ++iframe) {
 32            // 3 different moments
 33            for (unsigned int iobs=0; iobs<3; ++iobs) {
 34              const string name = "TMP/POL_"+toString(iups)+"_"+toString(iy)
 35                                 +"_"+toString(iframe)+"_"+toString(iobs);
 36              book(_p_Ups[iups][iy][iframe][iobs], name, refData(1,1,1));
 37            }
 38          }
 39        }
 40      }
 41    }
 42
 43    void findDecayProducts(const Particle& mother, unsigned int& nstable,
 44                           Particles & mup, Particles & mum) {
 45      for (const Particle& p : mother.children()) {
 46        int id = p.pid();
 47        if (id == PID::MUON ) {
 48          ++nstable;
 49          mum.push_back(p);
 50        }
 51        else if (id == PID::ANTIMUON) {
 52          ++nstable;
 53          mup.push_back(p);
 54        }
 55        else if (id == PID::PI0 || id == PID::K0S || id == PID::K0L ) {
 56          ++nstable;
 57        }
 58        else if ( !p.children().empty() ) {
 59          findDecayProducts(p, nstable, mup, mum);
 60        }
 61        else {
 62          ++nstable;
 63        }
 64      }
 65    }
 66
 67
 68    /// Perform the per-event analysis
 69    void analyze(const Event& event) {
 70      // find the beams
 71      const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
 72      // Final state of unstable particles to get particle spectra
 73      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
 74      for (const Particle& p : ufs.particles(Cuts::pid==553 or Cuts::pid==100553 or Cuts::pid==200553)) {
 75        // check mu+mu- decay and find muons
 76      	unsigned int nstable=0;
 77      	Particles mup,mum;
 78      	findDecayProducts(p,nstable,mup,mum);
 79      	if (mup.size()!=1 || mum.size()!=1 || nstable!=2) continue;
 80      	// pT and rapidity
 81      	const double absrap = p.absrap();
 82      	const double xp = p.perp();
 83        if (absrap>1.2 || xp<10. || xp>50.) continue;
 84      	// type of upsilon
 85      	unsigned int iups=p.pid()/100000;
 86        // rapidity interval
 87        unsigned int iy = absrap>0.6;
 88        // first the CS frame
 89        // first boost so upslion momentum =0 in z direction
 90        Vector3 beta = p.mom().betaVec();
 91        beta.setX(0.); beta.setY(0.);
 92        LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(beta);
 93        FourMomentum pp = boost.transform(p.mom());
 94        // and then transverse so pT=0
 95        beta = pp.betaVec();
 96        LorentzTransform boost2 = LorentzTransform::mkFrameTransformFromBeta(beta);
 97        // get all the momenta in this frame
 98        Vector3 muDirn = boost2.transform(boost.transform(mup[0].mom())).p3().unit();
 99        FourMomentum p1 = boost2.transform(boost.transform(beams. first.mom()));
100        FourMomentum p2 = boost2.transform(boost.transform(beams.second.mom()));
101        if (beams.first.mom().z()<0.) swap(p1,p2);
102        if (p.rapidity()<0.) swap(p1,p2);
103        Vector3 axisy = (p1.p3().cross(p2.p3())).unit();
104        Vector3 axisz(0.,0.,1.);
105        Vector3 axisx = axisy.cross(axisz);
106        double cTheta = axisz.dot(muDirn);
107        double cPhi   = axisx.dot(muDirn);
108        // fill the moments
109        _p_Ups[iups][iy][0][0]->fill(xp, 1.25*(3.*sqr(cTheta)-1.));
110        _p_Ups[iups][iy][0][1]->fill(xp, 1.25*(1.-sqr(cTheta))*(2.*sqr(cPhi)-1.));
111        _p_Ups[iups][iy][0][2]->fill(xp, 2.5 *cTheta*sqrt(1.-sqr(cTheta))*cPhi);
112        // now for the HX frame
113        beta = p.mom().betaVec();
114        boost = LorentzTransform::mkFrameTransformFromBeta(beta);
115        axisz = pp.p3().unit();
116        axisx = axisy.cross(axisz);
117        cTheta = axisz.dot(muDirn);
118        cPhi   = axisx.dot(muDirn);
119        // fill the moments
120        _p_Ups[iups][iy][1][0]->fill(xp, 1.25*(3.*sqr(cTheta)-1.));
121        _p_Ups[iups][iy][1][1]->fill(xp, 1.25*(1.-sqr(cTheta))*(2.*sqr(cPhi)-1.));
122        _p_Ups[iups][iy][1][2]->fill(xp, 2.5 *cTheta*sqrt(1.-sqr(cTheta))*cPhi);
123        // then PX
124        axisz = (p1.p3().unit()+p2.p3().unit()).unit();
125        axisx = axisy.cross(axisz);
126        cTheta = axisz.dot(muDirn);
127        cPhi   = axisx.dot(muDirn);
128        // fill the moments
129        _p_Ups[iups][iy][2][0]->fill(xp, 1.25*(3.*sqr(cTheta)-1.));
130        _p_Ups[iups][iy][2][1]->fill(xp, 1.25*(1.-sqr(cTheta))*(2.*sqr(cPhi)-1.));
131        _p_Ups[iups][iy][2][2]->fill(xp, 2.5 *cTheta*sqrt(1.-sqr(cTheta))*cPhi);
132      }
133    }
134
135
136    /// Normalise histograms etc., after the run
137    void finalize() {
138      // Loop over states
139      for (unsigned int iups=0; iups<3; ++iups) {
140        // Loop over rapidity ranges
141        for (unsigned int iy=0; iy<2; ++iy) {
142          // Loop over frame definition
143          for (unsigned int iframe=0; iframe<3; ++iframe) {
144            // base no for the ihistos in rivet
145            unsigned int ibase = 24*iups+iy+8*iframe;
146            // book scatters
147            Estimate1DPtr lTheta, lPhi, lThetaPhi, lTilde;
148            book(lTheta, ibase+1, 1, 1);
149            book(lPhi, ibase+3, 1, 1);
150            book(lThetaPhi, ibase+5, 1, 1);
151            book(lTilde   , ibase+7, 1, 1);
152            // histos for the moments
153            Profile1DPtr moment[3];
154            for (unsigned int ix=0;ix<3;++ix) {
155              moment[ix] = _p_Ups[iups][iy][iframe][ix];
156            }
157            // loop over bins
158            for (unsigned int ibin=1; ibin<=moment[0]->bins().size(); ++ibin) {
159              // extract moments and errors
160              double val[3], err[3];
161              // m1 = lTheta/(3+lTheta), m2 = lPhi/(3+lTheta), m3 = lThetaPhi/(3+lTheta)
162              for (unsigned int ix=0; ix<3; ++ix) {
163                bool has0 = moment[ix]->bin(ibin).numEntries()>0 && moment[ix]->bin(ibin).effNumEntries()>0;
164                bool has1 = moment[ix]->bin(ibin).numEntries()>1 && moment[ix]->bin(ibin).effNumEntries()>1;
165                val[ix] = has0? moment[ix]->bin(ibin).mean(2)   : 0.;
166                err[ix] = has1? moment[ix]->bin(ibin).stdErr(2) : 0.;
167              }
168              // values of the lambdas and their errors
169              double l1 = 3.*val[0]/(1.-val[0]);
170              double l2 = (3.+l1)*val[1];
171              // fill the scatters
172              lTheta->bin(ibin).setVal(l1);
173              lTheta->bin(ibin).setErr(3./sqr(1.-val[0])*err[0]);
174              lPhi  ->bin(ibin).setVal(l2);
175              lPhi  ->bin(ibin).setErr(3./sqr(1.-val[0])*sqrt(sqr(err[0]*val[1])+sqr(err[1]*(1.-val[0]))));
176              lThetaPhi->bin(ibin).setVal((3.+l1)*val[2]);
177              lThetaPhi->bin(ibin).setErr(3./sqr(1.-val[0])*sqrt(sqr(err[0]*val[1])+sqr(err[1]*(1.-val[0]))));
178              lTilde   ->bin(ibin).setVal((l1+3.*l2)/(1.-l2));
179              lTilde   ->bin(ibin).setErr( 3./sqr(1.-val[0]-3*val[1])*sqrt(sqr(err[0])+9.*sqr(err[1])));
180            }
181          }
182        }
183      }
184    }
185
186    /// @}
187
188
189    /// @name Histograms
190    /// @{
191    Profile1DPtr _p_Ups[3][2][3][3];
192    /// @}
193
194
195  };
196
197
198  RIVET_DECLARE_PLUGIN(CMS_2013_I1185414);
199
200}