Rivet analyses referenceLHCB_2013_I1244315J$/\psi$ polarization at 7 TeVExperiment: LHCB (LHC) Inspire ID: 1244315 Status: VALIDATED Authors:
Beam energies: (3500.0, 3500.0) GeV Run details:
Measurement of the polarization of J/$\psi$ at 7 TeV by LHCb Source code: LHCB_2013_I1244315.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 J/psi polarization at 7 TeV
10 class LHCB_2013_I1244315 : public Analysis {
11 public:
12
13 /// Constructor
14 RIVET_DEFAULT_ANALYSIS_CTOR(LHCB_2013_I1244315);
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 // histograms
26 // cross section
27 _ybins={2.0,2.5,3.0,3.5,4.0,4.5};
28 book(_h_Jpsi,_ybins);
29 for(unsigned int iy=0;iy<5;++iy) {
30 book(_h_Jpsi->bin(iy+1),4,1,1+iy);
31 for(unsigned int ix=0;ix<2;++ix) {
32 for(unsigned int iz=0;iz<3;++iz) {
33 string name="TMP/POL_"+toString(ix)+"_"+toString(iy)+"_"+toString(iz);
34 book(_p_Jpsi[ix][iy][iz],name,refData(1+ix,1,1+iy));
35 }
36 }
37 }
38 }
39
40 void findDecayProducts(const Particle & mother, unsigned int & nstable,
41 Particles & mup, Particles & mum) {
42 for(const Particle & p : mother.children()) {
43 int id = p.pid();
44 if (id == PID::MUON ) {
45 ++nstable;
46 mum.push_back(p);
47 }
48 else if (id == PID::ANTIMUON) {
49 ++nstable;
50 mup.push_back(p);
51 }
52 else if (id == PID::PI0 || id == PID::K0S || id == PID::K0L ) {
53 ++nstable;
54 }
55 else if ( !p.children().empty() ) {
56 findDecayProducts(p, nstable, mup, mum);
57 }
58 else
59 ++nstable;
60 }
61 }
62
63 /// Perform the per-event analysis
64 void analyze(const Event& event) {
65 // find the beams
66 const ParticlePair & beams = apply<Beam>(event, "Beams").beams();
67 // Final state of unstable particles to get particle spectra
68 const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
69 for (const Particle& p : ufs.particles(Cuts::pid==443)) {
70 // prompt
71 if(p.fromBottom()) continue;
72 // pT and rapidity
73 double rapidity = p.rapidity();
74 double xp = p.perp();
75 // cross section
76 _h_Jpsi->fill(rapidity,xp);
77 // polarization
78 unsigned int nstable=0;
79 Particles mup,mum;
80 findDecayProducts(p,nstable,mup,mum);
81 if(mup.size()!=1 || mum.size()!=1 || nstable!=2) continue;
82 // find the rapidity interval
83 if(rapidity<=2. || rapidity>=4.5) continue;
84 unsigned int iy=0;
85 for(iy=0;iy<5;++iy) if(rapidity<_ybins[iy+1]) break;
86 // first the CS frame
87 // first boost so upslion momentum =0 in z direction
88 Vector3 beta = p.momentum().betaVec();
89 beta.setX(0.);beta.setY(0.);
90 LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(beta);
91 FourMomentum pp = boost.transform(p.momentum());
92 // and then transverse so pT=0
93 beta = pp.betaVec();
94 LorentzTransform boost2 = LorentzTransform::mkFrameTransformFromBeta(beta);
95 // get all the momenta in this frame
96 Vector3 muDirn = boost2.transform(boost.transform(mup[0].momentum())).p3().unit();
97 FourMomentum p1 = boost2.transform(boost.transform(beams. first.momentum()));
98 FourMomentum p2 = boost2.transform(boost.transform(beams.second.momentum()));
99 if(beams.first.momentum().z()<0.) swap(p1,p2);
100 if(p.rapidity()<0.) swap(p1,p2);
101 Vector3 axisy = (p1.p3().cross(p2.p3())).unit();
102 Vector3 axisz(0.,0.,1.);
103 Vector3 axisx = axisy.cross(axisz);
104 double cTheta = axisz.dot(muDirn);
105 double cPhi = axisx.dot(muDirn);
106 // fill the moments
107 _p_Jpsi[1][iy][0]->fill(xp, 1.25*(3.*sqr(cTheta)-1.));
108 _p_Jpsi[1][iy][1]->fill(xp, 1.25*(1.-sqr(cTheta))*(2.*sqr(cPhi)-1.));
109 _p_Jpsi[1][iy][2]->fill(xp, 2.5 *cTheta*sqrt(1.-sqr(cTheta))*cPhi);
110 // now for the HX frame
111 beta = p.momentum().betaVec();
112 boost = LorentzTransform::mkFrameTransformFromBeta(beta);
113 axisz = pp.p3().unit();
114 axisx = axisy.cross(axisz);
115 cTheta = axisz.dot(muDirn);
116 cPhi = axisx.dot(muDirn);
117 // fill the moments
118 _p_Jpsi[0][iy][0]->fill(xp, 1.25*(3.*sqr(cTheta)-1.));
119 _p_Jpsi[0][iy][1]->fill(xp, 1.25*(1.-sqr(cTheta))*(2.*sqr(cPhi)-1.));
120 _p_Jpsi[0][iy][2]->fill(xp, 2.5 *cTheta*sqrt(1.-sqr(cTheta))*cPhi);
121 }
122 }
123
124
125 /// Normalise histograms etc., after the run
126 void finalize() {
127 double factor = crossSection()/nanobarn/sumOfWeights();
128 scale(_h_Jpsi,factor);
129 divByGroupWidth(_h_Jpsi);
130 for(unsigned int iy=0;iy<5;++iy) {
131 // Loop over frame definition
132 for(unsigned int iframe=0;iframe<2;++iframe) {
133 // book scatters
134 Estimate1DPtr lTheta,lPhi,lThetaPhi;
135 book(lTheta ,iframe+1,1,1+iy);
136 book(lPhi ,iframe+1,3,1+iy);
137 book(lThetaPhi,iframe+1,2,1+iy);
138 // histos for the moments
139 Profile1DPtr moment[3];
140 for(unsigned int ix=0;ix<3;++ix)
141 moment[ix] = _p_Jpsi[iframe][iy][ix];
142 // loop over bins
143 for(unsigned int ibin=1;ibin<=moment[0]->bins().size();++ibin) {
144 // extract moments and errors
145 double val[3],err[3];
146 // m1 = lTheta/(3+lTheta), m2 = lPhi/(3+lTheta), m3 = lThetaPhi/(3+lTheta)
147 for(unsigned int ix=0;ix<3;++ix) {
148 val[ix] = moment[ix]->bins()[ibin].numEntries()>0 && moment[ix]->bins()[ibin].effNumEntries()>0 ? moment[ix]->bins()[ibin].mean(2) : 0.;
149 err[ix] = moment[ix]->bins()[ibin].numEntries()>1 && moment[ix]->bins()[ibin].effNumEntries()>1 ? moment[ix]->bins()[ibin].stdErr(2) : 0.;
150 }
151 // values of the lambdas and their errors
152 double l1 = 3.*val[0]/(1.-val[0]);
153 double e1 = 3./sqr(1.-val[0])*err[0];
154 double l2 = (3.+l1)*val[1];
155 double e2 = 3./sqr(1.-val[0])*sqrt(sqr(err[0]*val[1])+sqr(err[1]*(1.-val[0])));
156 double l3 = (3.+l1)*val[2];
157 double e3 = 3./sqr(1.-val[0])*sqrt(sqr(err[0]*val[1])+sqr(err[1]*(1.-val[0])));
158 // fill the scatters
159 lTheta ->bin(ibin).setVal(l1);
160 lTheta ->bin(ibin).setErr(make_pair(e1,e1));
161 lPhi ->bin(ibin).setVal(l2);
162 lPhi ->bin(ibin).setErr(make_pair(e2,e2));
163 lThetaPhi->bin(ibin).setVal(l3);
164 lThetaPhi->bin(ibin).setErr(make_pair(e3,e3));
165 }
166 }
167 }
168 }
169 /// @}
170
171
172 /// @name Histograms
173 /// @{
174 Profile1DPtr _p_Jpsi[2][5][3];
175 Histo1DGroupPtr _h_Jpsi;
176 vector<double> _ybins;
177 /// @}
178
179
180 };
181
182
183 RIVET_DECLARE_PLUGIN(LHCB_2013_I1244315);
184
185}
|