Rivet analyses referenceLHCB_2012_I1119400Measurement of prompt hadron production ratios in $pp$ collisions at $\sqrt{s} = 0.9$ and 7 TeVExperiment: LHCb (LHC) Inspire ID: 1119400 Status: VALIDATED Authors:
Beam energies: (450.0, 450.0); (3500.0, 3500.0) GeV Run details:
WARNING --- particle map is incomplete for determination of promptness, assume unknown particles to have a life time of 0 --- Measurement of the production ratios of prompt charged particles (protons, pions and kaons). Promptness is defined as originating from the primary interaction, either directly, or through the subsequent decay of a resonance. Beam energy must be specified as analysis option "ENERGY" when rivet-merging samples. Source code: LHCB_2012_I1119400.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/ChargedFinalState.hh"
4
5namespace Rivet {
6
7
8 class LHCB_2012_I1119400 : public Analysis {
9 public:
10
11 /// @name Constructors etc.
12 /// @{
13
14 /// Constructor
15 LHCB_2012_I1119400() : Analysis("LHCB_2012_I1119400"),
16 _p_min(5.0),
17 _pt_min(0.0),_pt1_edge(0.8), _pt2_edge(1.2),
18 //_eta_nbins(4),
19 _eta_min(2.5),
20 _eta_max(4.5)
21 { }
22
23 /// @}
24
25
26 public:
27
28 /// @name Analysis methods
29 /// @{
30
31 /// Book histograms and initialise projections before the run
32 void init() {
33 fillMap(_partLftMap);
34
35 int id_shift = 0;
36 if (isCompatibleWithSqrtS(7000*GeV)) id_shift = 1;
37 // define ratios if second pdgid in pair is -1, it means that is a antiparticle/particle ratio
38
39 _ratiotype["pbarp"] = make_pair(2212, -1);
40 _ratiotype["kminuskplus"] = make_pair(321, -1);
41 _ratiotype["piminuspiplus"] = make_pair(211, -1);
42 _ratiotype["ppi"] = make_pair(2212, 211);
43 _ratiotype["kpi"] = make_pair(321, 211);
44 _ratiotype["pk"] = make_pair(2212, 321);
45
46 std::map<string, int > _hepdataid;
47 _hepdataid["pbarp"] = 1 + id_shift;
48 _hepdataid["kminuskplus"] = 3 + id_shift;
49 _hepdataid["piminuspiplus"] = 5 + id_shift;
50 _hepdataid["ppi"] = 7 + id_shift;
51 _hepdataid["kpi"] = 9 + id_shift;
52 _hepdataid["pk"] = 11 + id_shift;
53
54 std::map<std::string, std::pair<int, int> >::iterator it;
55
56 // booking histograms
57 for (it=_ratiotype.begin(); it!=_ratiotype.end(); it++) {
58 book(_h_ratio_lowpt [it->first], _hepdataid[it->first], 1, 1);
59 book(_h_ratio_midpt [it->first], _hepdataid[it->first], 1, 2);
60 book(_h_ratio_highpt[it->first], _hepdataid[it->first], 1, 3);
61 book(_h_num_lowpt [it->first], "TMP/num_l_"+it->first,refData(_hepdataid[it->first], 1, 1));
62 book(_h_num_midpt [it->first], "TMP/num_m_"+it->first,refData(_hepdataid[it->first], 1, 2));
63 book(_h_num_highpt [it->first], "TMP/num_h_"+it->first,refData(_hepdataid[it->first], 1, 3));
64 book(_h_den_lowpt [it->first], "TMP/den_l_"+it->first,refData(_hepdataid[it->first], 1, 1));
65 book(_h_den_midpt [it->first], "TMP/den_m_"+it->first,refData(_hepdataid[it->first], 1, 2));
66 book(_h_den_highpt [it->first], "TMP/den_h_"+it->first,refData(_hepdataid[it->first], 1, 3));
67 }
68
69 declare(ChargedFinalState(Cuts::etaIn(_eta_min, _eta_max) && Cuts::pT >= _pt_min*GeV), "CFS");
70 }
71
72
73 // Perform the per-event analysis
74 void analyze(const Event& event) {
75 const ChargedFinalState& cfs = apply<ChargedFinalState>(event, "CFS");
76
77 for (const Particle& p : cfs.particles()) {
78 int id = p.pid();
79 // continue if particle not a proton, a kaon or a pion
80 if ( !( (abs(id) == 211) || (abs(id) == 321) || (abs(id) == 2212))) {
81 continue;
82 }
83
84 // cut in momentum
85 const FourMomentum& qmom = p.momentum();
86 if (qmom.p3().mod() < _p_min) continue;
87
88 // Lifetime cut: ctau sum of all particle ancestors < 10^-9 m according to the paper (see eq. 5)
89 const double MAX_CTAU = 1.0e-9; // [m]
90 double ancestor_lftsum = getMotherLifeTimeSum(p);
91 if ( (ancestor_lftsum < 0.0) || (ancestor_lftsum > MAX_CTAU) ) continue;
92
93 double eta = qmom.eta();
94 double pT = qmom.pT();
95
96 std::map<std::string, std::pair<int, int> >::iterator it;
97
98 for (it=_ratiotype.begin(); it!=_ratiotype.end(); it++) {
99 // check what type of ratio is
100 if ((it->second.second)==-1) {
101 // check ptbin
102 if (pT < _pt1_edge) {
103 // filling histos for numerator and denominator
104 if (id == -abs(it->second.first)) _h_num_lowpt[it->first]->fill(eta);
105 if (id == abs(it->second.first)) _h_den_lowpt[it->first]->fill(eta);
106 }
107 else if (pT < _pt2_edge) {
108 // filling histos for numerator and denominator
109 if (id == -abs(it->second.first)) _h_num_midpt[it->first]->fill(eta);
110 if (id == abs(it->second.first)) _h_den_midpt[it->first]->fill(eta);
111 }
112 else {
113 // filling histos for numerator and denominator
114 if (id == -abs(it->second.first)) _h_num_highpt[it->first]->fill(eta);
115 if (id == abs(it->second.first)) _h_den_highpt[it->first]->fill(eta);
116 }
117 }
118 else {
119 // check what type of ratio is
120 if (pT < _pt1_edge) {
121 // filling histos for numerator and denominator
122 if (abs(id) == abs(it->second.first)) _h_num_lowpt[it->first]->fill(eta);
123 if (abs(id) == abs(it->second.second)) _h_den_lowpt[it->first]->fill(eta);
124 }
125 else if (pT < _pt2_edge) {
126 // filling histos for numerator and denominator
127 if (abs(id) == abs(it->second.first)) _h_num_midpt[it->first]->fill(eta);
128 if (abs(id) == abs(it->second.second)) _h_den_midpt[it->first]->fill(eta);
129 }
130 else {
131 // filling histos for numerator and denominator
132 if (abs(id) == abs(it->second.first)) _h_num_highpt[it->first]->fill(eta);
133 if (abs(id) == abs(it->second.second)) _h_den_highpt[it->first]->fill(eta);
134 }
135 }
136 }
137 }
138 }
139
140
141 // Generate the ratio histograms
142 void finalize() {
143 std::map<std::string, std::pair<int, int> >::iterator it;
144
145 // booking histograms
146 for (it=_ratiotype.begin(); it!=_ratiotype.end(); it++) {
147 divide(_h_num_lowpt[it->first], _h_den_lowpt[it->first], _h_ratio_lowpt[it->first]);
148 divide(_h_num_midpt[it->first], _h_den_midpt[it->first], _h_ratio_midpt[it->first]);
149 divide(_h_num_highpt[it->first], _h_den_highpt[it->first], _h_ratio_highpt[it->first]);
150 }
151 }
152
153 /// @}
154
155
156 private:
157
158
159 // Get particle lifetime from hardcoded data
160 double getLifeTime(int pid) {
161 pid = abs(pid);
162 double lft = -1.0;
163 map<int, double>::iterator pPartLft = _partLftMap.find(pid);
164 // search stable particle list
165 if (pPartLft == _partLftMap.end()) {
166 if (pid <= 100) return 0.0;
167 for (size_t i=0; i < sizeof(_stablePDGIds)/sizeof(unsigned int); i++) {
168 if (pid == _stablePDGIds[i]) {
169 lft = 0.0;
170 break;
171 }
172 }
173 } else {
174 lft = (*pPartLft).second;
175 }
176 if (lft < 0.0 && PID::isHadron(pid)) {
177 MSG_WARNING("Lifetime map imcomplete --- " << pid
178 << "... assume zero lifetime");
179 lft = 0.0;
180 }
181 return lft;
182 }
183
184 // Data members like post-cuts event weight counters go here
185 double getMotherLifeTimeSum(const Particle& p) {
186 if (p.genParticle() == nullptr) return -1.;
187 double lftSum = 0.;
188 double plft = 0.;
189 ConstGenParticlePtr part = p.genParticle();
190 ConstGenVertexPtr ivtx = part->production_vertex();
191 while(ivtx){
192
193 vector<ConstGenParticlePtr> part_in = HepMCUtils::particles(ivtx, Relatives::PARENTS);
194
195 if (part_in.size() < 1) { lftSum = -1.; break; };
196 part = part_in.at(0);
197 if ( !(part) ) { lftSum = -1.; break; };
198 ivtx = part->production_vertex();
199 if ( (part->pdg_id() == 2212) || !(ivtx) ) break; // reached beam
200 plft = getLifeTime(part->pdg_id());
201 if (plft < 0.) { lftSum = -1.; break; };
202 lftSum += plft;
203 };
204 return (lftSum * c_light);
205 }
206
207 /// @name Private variables
208 // Momentum threshold
209 double _p_min;
210
211 // The edges of the intervals of transversal momentum
212 double _pt_min;
213 double _pt1_edge;
214 double _pt2_edge;
215
216 // The limits of the pseudorapidity window
217 //int _eta_nbins;
218 double _eta_min;
219 double _eta_max;
220
221 // Map between PDG id and particle lifetimes in seconds
222 std::map<int, double> _partLftMap;
223
224 // Set of PDG Ids for stable particles (PDG Id <= 100 are considered stable)
225 static const int _stablePDGIds[205];
226
227 // Define histograms
228 // ratio
229 std::map<std::string, Estimate1DPtr > _h_ratio_lowpt;
230 std::map<std::string, Estimate1DPtr > _h_ratio_midpt;
231 std::map<std::string, Estimate1DPtr > _h_ratio_highpt;
232 // numerator
233 std::map<std::string, Histo1DPtr > _h_num_lowpt;
234 std::map<std::string, Histo1DPtr > _h_num_midpt;
235 std::map<std::string, Histo1DPtr > _h_num_highpt;
236 // denominator
237 std::map<std::string, Histo1DPtr > _h_den_lowpt;
238 std::map<std::string, Histo1DPtr > _h_den_midpt;
239 std::map<std::string, Histo1DPtr > _h_den_highpt;
240
241 // Map of ratios and IDs of numerator and denominator
242 std::map<string, pair<int,int> > _ratiotype;
243
244 // Fill the PDG Id to Lifetime[seconds] map
245 // Data was extracted from LHCb Particle Table through LHCb::ParticlePropertySvc
246 bool fillMap(map<int, double> &m) {
247 m[6] = 4.707703E-25; m[11] = 1.E+16; m[12] = 1.E+16;
248 m[13] = 2.197019E-06; m[14] = 1.E+16; m[15] = 2.906E-13; m[16] = 1.E+16; m[22] = 1.E+16;
249 m[23] = 2.637914E-25; m[24] = 3.075758E-25; m[25] = 9.4E-26; m[35] = 9.4E-26;
250 m[36] = 9.4E-26; m[37] = 9.4E-26; m[84] = 3.335641E-13; m[85] = 1.290893E-12;
251 m[111] = 8.4E-17; m[113] = 4.405704E-24; m[115] = 6.151516E-24; m[117] = 4.088275E-24;
252 m[119] = 2.102914E-24; m[130] = 5.116E-08; m[150] = 1.525E-12; m[211] = 2.6033E-08;
253 m[213] = 4.405704E-24; m[215] = 6.151516E-24; m[217] = 4.088275E-24; m[219] = 2.102914E-24;
254 m[221] = 5.063171E-19; m[223] = 7.752794E-23; m[225] = 3.555982E-24; m[227] = 3.91793E-24;
255 m[229] = 2.777267E-24; m[310] = 8.953E-11; m[313] = 1.308573E-23; m[315] = 6.038644E-24;
256 m[317] = 4.139699E-24; m[319] = 3.324304E-24; m[321] = 1.238E-08; m[323] = 1.295693E-23;
257 m[325] = 6.682357E-24; m[327] = 4.139699E-24; m[329] = 3.324304E-24; m[331] = 3.210791E-21;
258 m[333] = 1.545099E-22; m[335] = 9.016605E-24; m[337] = 7.565657E-24; m[350] = 1.407125E-12;
259 m[411] = 1.04E-12; m[413] = 6.856377E-21; m[415] = 1.778952E-23; m[421] = 4.101E-13;
260 m[423] = 1.000003E-19; m[425] = 1.530726E-23; m[431] = 5.E-13; m[433] = 1.000003E-19;
261 m[435] = 3.291061E-23; m[441] = 2.465214E-23; m[443] = 7.062363E-21; m[445] = 3.242425E-22;
262 m[510] = 1.525E-12; m[511] = 1.525E-12; m[513] = 1.000019E-19; m[515] = 1.31E-23;
263 m[521] = 1.638E-12; m[523] = 1.000019E-19; m[525] = 1.31E-23; m[530] = 1.536875E-12;
264 m[531] = 1.472E-12; m[533] = 1.E-19; m[535] = 1.31E-23; m[541] = 4.5E-13;
265 m[553] = 1.218911E-20; m[1112] = 4.539394E-24; m[1114] = 5.578069E-24; m[1116] = 1.994582E-24;
266 m[1118] = 2.269697E-24; m[1212] = 4.539394E-24; m[1214] = 5.723584E-24; m[1216] = 1.994582E-24;
267 m[1218] = 1.316424E-24; m[2112] = 8.857E+02; m[2114] = 5.578069E-24; m[2116] = 4.388081E-24;
268 m[2118] = 2.269697E-24; m[2122] = 4.539394E-24; m[2124] = 5.723584E-24; m[2126] = 1.994582E-24;
269 m[2128] = 1.316424E-24; m[2212] = 1.E+16; m[2214] = 5.578069E-24; m[2216] = 4.388081E-24;
270 m[2218] = 2.269697E-24; m[2222] = 4.539394E-24; m[2224] = 5.578069E-24; m[2226] = 1.994582E-24;
271 m[2228] = 2.269697E-24; m[3112] = 1.479E-10; m[3114] = 1.670589E-23; m[3116] = 5.485102E-24;
272 m[3118] = 3.656734E-24; m[3122] = 2.631E-10; m[102134] = 4.219309E-23; m[3126] = 8.227653E-24;
273 m[3128] = 3.291061E-24; m[3212] = 7.4E-20; m[3214] = 1.828367E-23; m[3216] = 5.485102E-24;
274 m[3218] = 3.656734E-24; m[3222] = 8.018E-11; m[3224] = 1.838582E-23; m[3226] = 5.485102E-24;
275 m[3228] = 3.656734E-24; m[3312] = 1.639E-10; m[3314] = 6.648608E-23; m[3322] = 2.9E-10;
276 m[3324] = 7.233101E-23; m[3334] = 8.21E-11; m[4112] = 2.991874E-22; m[4114] = 4.088274E-23;
277 m[4122] = 2.E-13; m[4132] = 1.12E-13; m[4212] = 3.999999E-22; m[4214] = 3.291061E-22;
278 m[4222] = 2.951624E-22; m[4224] = 4.417531E-23; m[4232] = 4.42E-13; m[4332] = 6.9E-14;
279 m[4412] = 3.335641E-13; m[4422] = 3.335641E-13; m[4432] = 3.335641E-13; m[5112] = 1.E-19;
280 m[5122] = 1.38E-12; m[5132] = 1.42E-12; m[5142] = 1.290893E-12; m[5212] = 1.E-19;
281 m[5222] = 1.E-19; m[5232] = 1.42E-12; m[5242] = 1.290893E-12; m[5312] = 1.E-19;
282 m[5322] = 1.E-19; m[5332] = 1.55E-12; m[5342] = 1.290893E-12; m[5442] = 1.290893E-12;
283 m[5512] = 1.290893E-12; m[5522] = 1.290893E-12; m[5532] = 1.290893E-12; m[5542] = 1.290893E-12;
284 m[10111] = 2.48382E-24; m[10113] = 4.635297E-24; m[10115] = 2.54136E-24; m[10211] = 2.48382E-24;
285 m[10213] = 4.635297E-24; m[10215] = 2.54136E-24; m[10223] = 1.828367E-24; m[10225] = 3.636531E-24;
286 m[10311] = 2.437823E-24; m[10313] = 7.313469E-24; m[10315] = 3.538775E-24;
287 m[10321] = 2.437823E-24; m[10323] = 7.313469E-24; m[10325] = 3.538775E-24;
288 m[10331] = 4.804469E-24; m[10411] = 4.38E-24; m[10413] = 3.29E-23; m[10421] = 4.38E-24;
289 m[10423] = 3.22653E-23; m[10431] = 6.5821E-22; m[10433] = 6.5821E-22; m[10441] = 6.453061E-23;
290 m[10511] = 4.39E-24; m[10513] = 1.65E-23; m[10521] = 4.39E-24; m[10523] = 1.65E-23;
291 m[10531] = 4.39E-24; m[10533] = 1.65E-23; m[11114] = 2.194041E-24; m[11116] = 1.828367E-24;
292 m[11212] = 1.880606E-24; m[11216] = 1.828367E-24; m[12112] = 2.194041E-24;
293 m[12114] = 2.194041E-24; m[12116] = 5.063171E-24; m[12126] = 1.828367E-24;
294 m[12212] = 2.194041E-24; m[12214] = 2.194041E-24; m[12216] = 5.063171E-24;
295 m[12224] = 2.194041E-24; m[12226] = 1.828367E-24; m[13112] = 6.582122E-24; m[13114] = 1.09702E-23;
296 m[13116] = 5.485102E-24; m[13122] = 1.316424E-23; m[13124] = 1.09702E-23; m[13126] = 6.928549E-24;
297 m[13212] = 6.582122E-24; m[13214] = 1.09702E-23; m[13216] = 5.485102E-24; m[13222] = 6.582122E-24;
298 m[13224] = 1.09702E-23; m[13226] = 5.485102E-24; m[13314] = 2.742551E-23;
299 m[13324] = 2.742551E-23; m[102142] = 1.828367E-22; m[20022] = 1.E+16; m[20113] = 1.567172E-24;
300 m[20213] = 1.567172E-24; m[20223] = 2.708692E-23; m[20313] = 3.782829E-24;
301 m[20315] = 2.384827E-24; m[20323] = 3.782829E-24; m[20325] = 2.384827E-24;
302 m[20333] = 1.198929E-23; m[20413] = 2.63E-24; m[20423] = 2.63E-24; m[20433] = 6.5821E-22;
303 m[20443] = 7.395643E-22; m[20513] = 2.63E-24; m[20523] = 2.63E-24; m[20533] = 2.63E-24;
304 m[21112] = 2.632849E-24; m[21114] = 3.291061E-24; m[21212] = 2.632849E-24;
305 m[21214] = 6.582122E-24; m[22112] = 4.388081E-24; m[22114] = 3.291061E-24;
306 m[22122] = 2.632849E-24; m[22124] = 6.582122E-24; m[22212] = 4.388081E-24;
307 m[22214] = 3.291061E-24; m[22222] = 2.632849E-24; m[22224] = 3.291061E-24;
308 m[23112] = 7.313469E-24; m[23114] = 2.991874E-24; m[23122] = 4.388081E-24;
309 m[23124] = 6.582122E-24; m[23126] = 3.291061E-24; m[23212] = 7.313469E-24;
310 m[23214] = 2.991874E-24; m[23222] = 7.313469E-24; m[23224] = 2.991874E-24;
311 m[30113] = 2.632849E-24; m[30213] = 2.632849E-24; m[30221] = 1.880606E-24;
312 m[30223] = 2.089563E-24; m[30313] = 2.056913E-24; m[30323] = 2.056913E-24;
313 m[30443] = 2.419898E-23; m[31114] = 1.880606E-24; m[31214] = 3.291061E-24;
314 m[32112] = 3.989164E-24; m[32114] = 1.880606E-24; m[32124] = 3.291061E-24;
315 m[32212] = 3.989164E-24; m[32214] = 1.880606E-24; m[32224] = 1.880606E-24;
316 m[33122] = 1.880606E-23; m[42112] = 6.582122E-24; m[42212] = 6.582122E-24;
317 m[43122] = 2.194041E-24; m[53122] = 4.388081E-24; m[100111] = 1.645531E-24;
318 m[100113] = 1.64553E-24; m[100211] = 1.645531E-24; m[100213] = 1.64553E-24;
319 m[100221] = 1.196749E-23; m[100223] = 3.061452E-24; m[100313] = 2.837122E-24;
320 m[100323] = 2.837122E-24; m[100331] = 4.459432E-25; m[100333] = 4.388081E-24;
321 m[100441] = 4.701516E-23; m[100443] = 2.076379E-21; m[100553] = 2.056913E-20;
322 m[200553] = 3.242425E-20; m[300553] = 3.210791E-23; m[9000111] = 8.776163E-24;
323 m[9000211] = 8.776163E-24; m[9000443] = 8.227652E-24; m[9000553] = 5.983747E-24;
324 m[9010111] = 3.164482E-24; m[9010211] = 3.164482E-24; m[9010221] = 9.403031E-24;
325 m[9010443] = 8.438618E-24; m[9010553] = 8.3318E-24; m[9020443] = 1.061633E-23;
326 m[9030221] = 6.038644E-24; m[9042413] = 2.07634E-21; m[9050225] = 1.394517E-24;
327 m[9060225] = 3.291061E-24; m[9080225] = 4.388081E-24; m[9090225] = 2.056913E-24;
328 m[9910445] = 2.07634E-21; m[9920443] = 2.07634E-21;
329 return true;
330 }
331
332 };
333
334
335 const int LHCB_2012_I1119400::_stablePDGIds[205] = {
336 311, 543, 545, 551, 555, 557, 1103, 2101, 2103, 2203, 3101, 3103, 3201, 3203, 3303,
337 4101, 4103, 102144, 4201, 4203, 4301, 4303, 4312, 4314, 4322, 4324, 4334, 4403, 4414,
338 4424, 4434, 4444, 5101, 5103, 5114, 5201, 5203, 5214, 5224, 5301, 5303, 5314, 5324,
339 5334, 5401, 5403, 5412, 5414, 5422, 5424, 5432, 5434, 5444, 5503, 5514, 5524, 5534,
340 5544, 5554, 10022, 10333, 10335, 10443, 10541, 10543, 10551, 10553, 10555, 11112,
341 12118, 12122, 12218, 12222, 13316, 13326, 20543, 20553, 20555, 23314, 23324, 30343,
342 30353, 30363, 30553, 33314, 33324, 41214, 42124, 52114, 52214, 100311, 100315, 100321,
343 100325, 100411, 100413, 100421, 100423, 100551, 100555, 100557, 110551, 110553, 110555,
344 120553, 120555, 130553, 200551, 200555, 210551, 210553, 220553, 1000001, 1000002,
345 1000003, 1000004, 1000005, 1000006, 1000011, 1000012, 1000013, 1000014, 1000015,
346 1000016, 1000021, 1000022, 1000023, 1000024, 1000025, 1000035, 1000037, 1000039,
347 2000001, 2000002, 2000003, 2000004, 2000005, 2000006, 2000011, 2000012, 2000013,
348 2000014, 2000015, 2000016, 3000111, 3000113, 3000211, 3000213, 3000221, 3000223,
349 3000331, 3100021, 3100111, 3100113, 3200111, 3200113, 3300113, 3400113, 4000001,
350 4000002, 4000011, 4000012, 5000039, 9000221, 9900012, 9900014, 9900016, 9900023,
351 9900024, 9900041, 9900042 };
352
353
354 RIVET_DECLARE_PLUGIN(LHCB_2012_I1119400);
355
356}
|