rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

LHCB_2017_I1621596

Measurement of $\Upsilon(1,2,3S)$ polarization at 7 and 8 TeV
Experiment: LHCB (LHC)
Inspire ID: 1621596
Status: VALIDATED
Authors:
  • Peter Richardson
References:
  • JHEP 12 (2017) 110
Beams: p+ p+
Beam energies: (3500.0, 3500.0); (4000.0, 4000.0) GeV
Run details:
  • Upsilon production

Measurement of the polarization of $\Upsilon(1,2,3)$ at 7 and 8 TeV by LHCb. Beam energy must be specified as analysis option "ENERGY" when rivet-merging samples.

Source code: LHCB_2017_I1621596.cc
  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 and 8 TeV
 10  class LHCB_2017_I1621596 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(LHCB_2017_I1621596);
 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      int iloc=-1;
 26      if (isCompatibleWithSqrtS(7000)) {
 27	iloc = 0;
 28      }
 29      else if  (isCompatibleWithSqrtS(8000)) {
 30	iloc = 1;
 31      }
 32      else
 33	throw UserError("Centre-of-mass energy of the given input is neither 7 or 8 TeV.");
 34      // histograms
 35      _ybins={2.2,3.0,3.5,4.5};
 36      for(unsigned int iups=0;iups<3;++iups) {
 37	for(unsigned int iframe=0;iframe<3;++iframe) {
 38	  for(unsigned int imom=0;imom<3;++imom) {
 39	    for(unsigned int iy=0;iy<3;++iy) {
 40	      book(_p_Upsilon[iups][iframe][iy][imom],
 41		   "TMP/UPS_"+toString(iups)+"_"+toString(iframe)+"_"+toString(iy)+"_"+toString(imom),
 42		   refData(32*iups+4*iloc+8*iframe+1,1,iy+1));
 43	    }
 44	    book(_p_Upsilon[iups][iframe][3][imom],
 45		 "TMP/UPS_"+toString(iups)+"_"+toString(iframe)+"_3_"+toString(imom),
 46		 refData(32*iups+4*iloc+25,1,iframe+1));
 47	  }
 48	}
 49      }
 50    }
 51
 52    void findDecayProducts(const Particle & mother, unsigned int & nstable,  
 53                           Particles & mup, Particles & mum) {
 54      for(const Particle & p : mother.children()) {
 55        int id = p.pid();
 56        if (id == PID::MUON ) {
 57          ++nstable;
 58    	  mum.push_back(p);
 59    	}
 60        else if (id == PID::ANTIMUON) {
 61          ++nstable;
 62    	  mup.push_back(p);
 63        }
 64        else if (id == PID::PI0 || id == PID::K0S || id == PID::K0L ) {
 65          ++nstable;
 66        }
 67        else if ( !p.children().empty() ) {
 68          findDecayProducts(p, nstable, mup, mum);
 69        }
 70        else
 71          ++nstable;
 72      }
 73    }
 74
 75    /// Perform the per-event analysis
 76    void analyze(const Event& event) {
 77      // find the beams
 78      const ParticlePair & beams = apply<Beam>(event, "Beams").beams();
 79      // Final state of unstable particles to get particle spectra
 80      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
 81      for (const Particle& p : ufs.particles(Cuts::pid==553 || Cuts::pid==100553 || Cuts::pid==200553)) {
 82      	// pT and rapidity
 83      	double rapidity = p.rapidity();
 84      	double xp = p.perp();
 85	if(rapidity<2.2 || rapidity >4.5) continue;
 86	// which upsilon
 87	unsigned int iups=p.pid()/100000;
 88	// polarization
 89      	unsigned int nstable=0;
 90      	Particles mup,mum;
 91      	findDecayProducts(p,nstable,mup,mum);
 92      	if(mup.size()!=1 || mum.size()!=1 || nstable!=2) continue;
 93	unsigned int iy=0;
 94	for(iy=0;iy<3;++iy) if(rapidity<_ybins[iy+1]) break;
 95	// first the CS frame
 96	// first boost so upslion momentum =0 in z direction
 97	Vector3 beta = p.momentum().betaVec();
 98	beta.setX(0.);beta.setY(0.);
 99	LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(beta);
100	FourMomentum pp = boost.transform(p.momentum());
101	// and then transverse so pT=0
102        beta = pp.betaVec();
103	LorentzTransform boost2 = LorentzTransform::mkFrameTransformFromBeta(beta);
104	// get all the momenta in this frame
105	Vector3 muDirn = boost2.transform(boost.transform(mup[0].momentum())).p3().unit();
106	FourMomentum p1 = boost2.transform(boost.transform(beams. first.momentum()));
107	FourMomentum p2 = boost2.transform(boost.transform(beams.second.momentum()));
108	if(beams.first.momentum().z()<0.) swap(p1,p2);
109	if(p.rapidity()<0.) swap(p1,p2);
110	Vector3 axisy = (p1.p3().cross(p2.p3())).unit();
111	Vector3 axisz(0.,0.,1.);
112	Vector3 axisx = axisy.cross(axisz);
113	double cTheta = axisz.dot(muDirn);
114	double cPhi   = axisx.dot(muDirn);
115	// fill the moments
116	_p_Upsilon[iups][1][iy][0]->fill(xp, 1.25*(3.*sqr(cTheta)-1.));
117	_p_Upsilon[iups][1][iy][1]->fill(xp, 1.25*(1.-sqr(cTheta))*(2.*sqr(cPhi)-1.));
118	_p_Upsilon[iups][1][iy][2]->fill(xp, 2.5 *cTheta*sqrt(1.-sqr(cTheta))*cPhi);
119	_p_Upsilon[iups][1][3 ][0]->fill(xp, 1.25*(3.*sqr(cTheta)-1.));
120	_p_Upsilon[iups][1][3 ][1]->fill(xp, 1.25*(1.-sqr(cTheta))*(2.*sqr(cPhi)-1.));
121	_p_Upsilon[iups][1][3 ][2]->fill(xp, 2.5 *cTheta*sqrt(1.-sqr(cTheta))*cPhi);
122	// Gottfried-Jackson frame
123	axisz = p1.p3().unit();
124	axisx = axisy.cross(axisz);
125	cTheta = axisz.dot(muDirn);
126	cPhi   = axisx.dot(muDirn);
127	// fill the moments
128	_p_Upsilon[iups][2][iy][0]->fill(xp, 1.25*(3.*sqr(cTheta)-1.));
129	_p_Upsilon[iups][2][iy][1]->fill(xp, 1.25*(1.-sqr(cTheta))*(2.*sqr(cPhi)-1.));
130	_p_Upsilon[iups][2][iy][2]->fill(xp, 2.5 *cTheta*sqrt(1.-sqr(cTheta))*cPhi);
131	_p_Upsilon[iups][2][3 ][0]->fill(xp, 1.25*(3.*sqr(cTheta)-1.));
132	_p_Upsilon[iups][2][3 ][1]->fill(xp, 1.25*(1.-sqr(cTheta))*(2.*sqr(cPhi)-1.));
133	_p_Upsilon[iups][2][3 ][2]->fill(xp, 2.5 *cTheta*sqrt(1.-sqr(cTheta))*cPhi);
134	// now for the HX frame
135	beta = p.momentum().betaVec();
136	boost = LorentzTransform::mkFrameTransformFromBeta(beta);
137	axisz = pp.p3().unit();
138	axisx = axisy.cross(axisz);
139	cTheta = axisz.dot(muDirn);
140	cPhi   = axisx.dot(muDirn);
141	// fill the moments
142	_p_Upsilon[iups][0][iy][0]->fill(xp, 1.25*(3.*sqr(cTheta)-1.));
143	_p_Upsilon[iups][0][iy][1]->fill(xp, 1.25*(1.-sqr(cTheta))*(2.*sqr(cPhi)-1.));
144	_p_Upsilon[iups][0][iy][2]->fill(xp, 2.5 *cTheta*sqrt(1.-sqr(cTheta))*cPhi);
145	_p_Upsilon[iups][0][3 ][0]->fill(xp, 1.25*(3.*sqr(cTheta)-1.));
146	_p_Upsilon[iups][0][3 ][1]->fill(xp, 1.25*(1.-sqr(cTheta))*(2.*sqr(cPhi)-1.));
147	_p_Upsilon[iups][0][3 ][2]->fill(xp, 2.5 *cTheta*sqrt(1.-sqr(cTheta))*cPhi);
148      }
149    }
150
151    /// Normalise histograms etc., after the run
152    void finalize() {
153      int iloc=-1;
154      if (isCompatibleWithSqrtS(7000)) {
155	iloc = 0;
156      }
157      else if  (isCompatibleWithSqrtS(8000)) {
158	iloc = 1;
159      }
160      // loop over upslion
161      for(unsigned int iups=0;iups<3;++iups) {
162	// loop over iframe
163	for(unsigned int iframe=0;iframe<3;++iframe) {
164	  unsigned int ibase = 32*iups+4*iloc+8*iframe;
165	  unsigned int ibase2 = 32*iups+4*iloc+24;
166	  // rapidity range
167	  for(unsigned int iy=0;iy<4;++iy) {
168	    // book scatters
169	    Estimate1DPtr lTheta,lPhi,lThetaPhi,lTilde;
170	    if(iy<3) {
171	      book(lTheta   ,ibase+1,1,1+iy);
172	      book(lPhi     ,ibase+3,1,1+iy);
173	      book(lThetaPhi,ibase+2,1,1+iy);
174	      book(lTilde   ,ibase+4,1,1+iy);
175	    }
176	    else {
177	      book(lTheta   ,ibase2+1,1,1+iframe);
178	      book(lPhi     ,ibase2+3,1,1+iframe);
179	      book(lThetaPhi,ibase2+2,1,1+iframe);
180	      book(lTilde   ,ibase2+4,1,1+iframe);
181	    }
182	    // histos for the moments
183	    Profile1DPtr moment[3];
184	    for(unsigned int ix=0;ix<3;++ix)
185	      moment[ix] = _p_Upsilon[iups][iframe][iy][ix];
186	    // loop over bins
187	    for(unsigned int ibin=1;ibin<=moment[0]->numBins();++ibin) {
188	      // extract moments and errors
189	      double val[3],err[3];
190	      // m1 = lTheta/(3+lTheta), m2 = lPhi/(3+lTheta), m3 = lThetaPhi/(3+lTheta)
191	      for(unsigned int ix=0;ix<3;++ix) {
192		val[ix] = moment[ix]->bins()[ibin].numEntries()>0 &&
193                  moment[ix]->bins()[ibin].effNumEntries()>0 ? moment[ix]->bins()[ibin].mean(2)   : 0.;
194		err[ix] = moment[ix]->bins()[ibin].numEntries()>1 && moment[ix]->bins()[ibin].effNumEntries()>1 ? moment[ix]->bins()[ibin].stdErr(2) : 0.;
195	      }
196	      // values of the lambdas and their errors
197              double l1 = 3.*val[0]/(1.-val[0]);
198              double l2 = (3.+l1)*val[1];
199	      lTheta   ->bin(ibin).setVal(l1);
200              lTheta   ->bin(ibin).setErr(3./sqr(1.-val[0])*err[0]);
201	      lPhi     ->bin(ibin).setVal(l2);
202	      lPhi     ->bin(ibin).setErr(3./sqr(1.-val[0])*sqrt(sqr(err[0]*val[1])+sqr(err[1]*(1.-val[0]))));
203	      lThetaPhi->bin(ibin).setVal((3.+l1)*val[2]);
204	      lThetaPhi->bin(ibin).setErr(3./sqr(1.-val[0])*sqrt(sqr(err[0]*val[1])+sqr(err[1]*(1.-val[0]))));
205	      lTilde   ->bin(ibin).setVal((l1+3.*l2)/(1.-l2));
206	      lTilde   ->bin(ibin).setErr(3./sqr(1.-val[0]-3*val[1])*sqrt(sqr(err[0])+9.*sqr(err[1])));
207	    }
208	  }
209	}
210      }
211    }
212
213    /// @}
214
215
216    /// @name Histograms
217    /// @{
218    Profile1DPtr _p_Upsilon[3][3][4][3];
219    vector<double>  _ybins;
220    /// @}
221
222
223  };
224
225
226  RIVET_DECLARE_PLUGIN(LHCB_2017_I1621596);
227
228}