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