Rivet analyses referenceLHCB_2017_I1621596Measurement of $\Upsilon(1,2,3S)$ polarization at 7 and 8 TeVExperiment: LHCB (LHC) Inspire ID: 1621596 Status: VALIDATED Authors:
Beam energies: (3500.0, 3500.0); (4000.0, 4000.0) GeV Run details:
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}
|