Rivet analyses referenceDELPHI_2001_I526164Charged and Identified particle spectra in $q\bar{q}$ and $WW$ eventsExperiment: DELPHI (LEP) Inspire ID: 526164 Status: VALIDATED Authors:
Beam energies: (66.5, 66.5); (80.5, 80.5); (86.0, 86.0); (91.5, 91.5); (94.5, 94.5); (100.0, 100.0) GeV Run details:
Measurement of charged and identified particle production in both $q\bar{q}$ and $W^+W^-$ events at LEP. Need $q\bar{q}$ at all the energies and $W^+W^-$ at 183 and 189 GeV. Source code: DELPHI_2001_I526164.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/UnstableParticles.hh"
4#include "Rivet/Projections/FinalState.hh"
5#include "Rivet/Projections/Thrust.hh"
6
7namespace Rivet {
8
9
10 /// @brief DELPHI W decay analysis
11 class DELPHI_2001_I526164 : public Analysis {
12 public:
13
14 /// Constructor
15 RIVET_DEFAULT_ANALYSIS_CTOR(DELPHI_2001_I526164);
16
17
18 /// @name Analysis methods
19 ///@{
20
21 /// Book histograms and initialise projections before the run
22 void init() {
23 // projections
24 declare(UnstableParticles(),"UFS");
25 declare(FinalState(),"FS");
26 // Hisotgram booking
27 // qqbar
28 // counters for multiplicities
29 for(unsigned int ix=0;ix<7;++ix) {
30 std::ostringstream title;
31 title << "_n_qq_" << ix;
32 book(_n_qq[ix],title.str());
33 }
34 // spectra
35 _WW=false;
36 if(isCompatibleWithSqrtS(183.)) {
37 book(_h_qq_K0 ,16,1,1);
38 book(_h_qq_Lam,16,1,3);
39 book(_h_p_charged[0],7,1,2);
40 book(_h_p_charged[1],7,1,1);
41 book(_h_p_chargedB[0],"/TMP/h_7_1_2",refData(8,1,1));
42 book(_h_p_chargedB[1],"/TMP/h_7_1_1",refData(8,1,1));
43 book(_h_xi_charged [0],10,1,2);
44 book(_h_xi_charged [1],10,1,1);
45 book(_h_xi_chargedB[0],"/TMP/h_10_1_2",refData(10,1,3));
46 book(_h_xi_chargedB[1],"/TMP/h_10_1_1",refData(10,1,3));
47 book(_h_pT_charged [0],12,1,2);
48 book(_h_pT_charged [1],12,1,1);
49 book(_h_pT_chargedB[0],"/TMP/h_12_1_2",refData(12,1,3));
50 book(_h_pT_chargedB[1],"/TMP/h_12_1_1",refData(12,1,3));
51 _WW=true;
52 }
53 else if (isCompatibleWithSqrtS(189.)) {
54 book(_h_qq_K0 ,16,1,2);
55 book(_h_qq_Lam,16,1,4);
56 book(_h_p_charged [0],5,1,2);
57 book(_h_p_charged [1],5,1,1);
58 book(_h_p_chargedB[0],"/TMP/h_5_1_2",refData(6,1,1));
59 book(_h_p_chargedB[1],"/TMP/h_5_1_1",refData(6,1,1));
60 book(_h_xi_charged [0],9,1,2);
61 book(_h_xi_charged [1],9,1,1);
62 book(_h_xi_chargedB[0],"/TMP/h_9_1_2",refData(9,1,3));
63 book(_h_xi_chargedB[1],"/TMP/h_9_1_1",refData(9,1,3));
64 book(_h_pT_charged [0],11,1,2);
65 book(_h_pT_charged [1],11,1,1);
66 book(_h_pT_chargedB[0],"/TMP/h_11_1_2",refData(11,1,3));
67 book(_h_pT_chargedB[1],"/TMP/h_11_1_1",refData(11,1,3));
68 for(unsigned int ix=0;ix<3;++ix) {
69 for(unsigned int iy=0;iy<4;++iy) {
70 book(_h_xi_ident[ix][iy],13+ix,1,iy+1);
71 }
72 }
73 _WW=true;
74 }
75 if(_WW) {
76 for(unsigned int ix=0;ix<2;++ix) {
77 for(unsigned int iy=0;iy<5;++iy) {
78 std::ostringstream title;
79 title << "_n_WW_" << ix << "_" << iy << "\n";
80 book(_n_WW[ix][iy],title.str());
81 }
82 }
83 }
84 }
85
86 /// Perform the per-event analysis
87 void analyze(const Event& event) {
88 const double mW=80.379*GeV;
89 struct ParticleOrdering {
90 bool operator()(const Particle &p1, const Particle & p2) const {
91 return p1.pid() > p2.pid();
92 }
93 };
94 multimap<Particle,Particles,ParticleOrdering> Wbosons;
95 // loop over final state particles to find W's
96 Particles finalState = apply<FinalState>(event,"FS").particles();
97 // loop over FS particles
98 for(const Particle & p : finalState) {
99 Particle parent=p;
100 while(!parent.parents().empty()) {
101 parent=parent.parents()[0];
102 if(parent.abspid()==24 && parent.mass()>20.) break;
103 }
104 if(parent.abspid()!=24) continue;
105 // find those which came from W
106 bool found=false;
107 for (auto & W : Wbosons) {
108 // W already in list add particle to its decay products
109 if (fuzzyEquals(W.first.momentum(),parent.momentum())) {
110 W.second.push_back(p);
111 found=true;
112 break;
113 }
114 }
115 if(!found) {
116 // check W not child
117 bool Wchild=false;
118 for(const Particle & child : parent.children()) {
119 if(child.abspid()==24) {
120 Wchild=true;
121 break;
122 }
123 }
124 // add to list
125 if(!Wchild) {
126 Particles temp = {p};
127 Wbosons.insert(make_pair(parent,temp));
128 }
129 }
130 }
131 // no W's => q qbar event
132 if (Wbosons.empty()) {
133 _n_qq[0]->fill();
134 UnstableParticles ufs=apply<UnstableParticles>(event,"UFS");
135 for(const Particle & p : ufs.particles()) {
136 if(p.abspid()==PID::PIPLUS) {
137 _n_qq[1]->fill();
138 _n_qq[2]->fill();
139 }
140 else if(p.abspid()==PID::KPLUS) {
141 _n_qq[1]->fill();
142 _n_qq[3]->fill();
143 }
144 else if(p.abspid()==PID::PROTON) {
145 _n_qq[1]->fill();
146 _n_qq[5]->fill();
147 }
148 else if(p.abspid()==130 || p.abspid()==310 ) {
149 _n_qq[4]->fill();
150 if(_h_qq_K0) _h_qq_K0->fill(-log(2.*p.p3().mod()/sqrtS()));
151 }
152 else if(p.abspid()==PID::LAMBDA ) {
153 _n_qq[6]->fill();
154 if(_h_qq_Lam) _h_qq_Lam->fill(-log(2.*p.p3().mod()/sqrtS()));
155 }
156 }
157 }
158 else if (Wbosons.size()==2) {
159 bool leptonic[2] = {false,false};
160 unsigned int iboson=0;
161 for(auto & W : Wbosons) {
162 for(const Particle & child : W.first.children()) {
163 // find lepton bosons
164 if(child.abspid()==11 || child.abspid()==13) {
165 leptonic[iboson]=true;
166 break;
167 }
168 // veto events with W-> tau decays
169 else if(child.abspid()==15 )
170 vetoEvent;
171 }
172 ++iboson;
173 }
174 // only fully hadronic or semi-leptonic events
175 if(leptonic[0] && leptonic[1]) vetoEvent;
176 // 0 SL 1 hadronic
177 unsigned int itype= leptonic[0]||leptonic[1] ? 0 : 1;
178 _n_WW[itype][0]->fill();
179 Particles particles;
180 iboson=0;
181 for(auto & W : Wbosons) {
182 if (!leptonic[iboson])
183 particles.insert(particles.end(),W.second.begin(),W.second.end());
184 ++iboson;
185 }
186 // calculate thrust
187 Thrust thrust;
188 thrust.calc(particles);
189 // type of event
190 // loop over particles and fill histos
191 for(const Particle & p : particles) {
192 if(!PID::isCharged(p.pid())) continue;
193 // Get momentum of each particle.
194 const Vector3 mom3 = p.p3();
195 // Scaled momenta.
196 const double mom = mom3.mod();
197 const double xiW = -log(2.*mom/mW);
198 const double xiL = -log(2.*mom/sqrtS());
199 // Get momenta components w.r.t. thrust
200 const double pTinT = dot(mom3, thrust.thrustMajorAxis());
201 const double pToutT = dot(mom3, thrust.thrustMinorAxis());
202 double pT = sqrt(sqr(pTinT)+sqr(pToutT));
203 // fill charged particle hists
204 _h_p_charged [itype]->fill(mom);
205 _h_xi_charged [itype]->fill(xiL);
206 _h_pT_charged [itype]->fill(pT );
207 _h_p_chargedB [itype]->fill(mom);
208 _h_xi_chargedB[itype]->fill(xiL);
209 _h_pT_chargedB[itype]->fill(pT );
210 // and identified particles
211 if(_h_xi_ident[itype][0]) _h_xi_ident[itype][0]->fill(xiW);
212 _n_WW[itype][1]->fill();
213 if(p.abspid()==PID::PIPLUS) {
214 if(_h_xi_ident[itype][1]) _h_xi_ident[itype][1]->fill(xiW);
215 _n_WW[itype][2]->fill();
216 }
217 else if(p.abspid()==PID::KPLUS) {
218 if(_h_xi_ident[itype][2]) _h_xi_ident[itype][2]->fill(xiW);
219 _n_WW[itype][3]->fill();
220 }
221 else if(p.abspid()==PID::PROTON){
222 if(_h_xi_ident[itype][3]) _h_xi_ident[itype][3]->fill(xiW);
223 _n_WW[itype][4]->fill();
224 }
225 }
226 // boosted specta in W rest frame
227 if(itype==0) {
228 iboson=0;
229 for(auto & W : Wbosons) {
230 ++iboson;
231 if (leptonic[iboson-1]) continue;
232 // boost to rest frame
233 LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(W.first.momentum().betaVec());
234 FourMomentum psum;
235 for(const Particle & p :W.second) psum+=p.momentum();
236 // spectra
237 for(const Particle & p : W.second) {
238 if(!PID::isCharged(p.pid())) continue;
239 // Get momentum of each particle.
240 FourMomentum phad = boost.transform(p.momentum());
241 const Vector3 mom3 = phad.p3();
242 // Scaled momenta.
243 const double mom = mom3.mod();
244 const double scaledMom = 2.*mom/mW;
245 const double xi = -log(scaledMom);
246 // /identified particle spectra
247 if(_h_xi_ident[2][0]) _h_xi_ident[2][0]->fill(xi);
248 if(p.abspid()==PID::PIPLUS) {
249 if(_h_xi_ident[2][1]) _h_xi_ident[2][1]->fill(xi);
250 }
251 else if(p.abspid()==PID::KPLUS) {
252 if(_h_xi_ident[2][2]) _h_xi_ident[2][2]->fill(xi);
253 }
254 else if(p.abspid()==PID::PROTON){
255 if(_h_xi_ident[2][3]) _h_xi_ident[2][3]->fill(xi);
256 }
257 }
258 }
259 }
260 }
261 }
262
263
264 /// Normalise histograms etc., after the run
265 void finalize() {
266 // qq
267 if(_n_qq[0]->effNumEntries()!=0.) {
268 // scale spectra if needed
269 if(_h_qq_K0 ) scale(_h_qq_K0 , 1./ *_n_qq[0]);
270 if(_h_qq_Lam) scale(_h_qq_Lam, 1./ *_n_qq[0]);
271 for(unsigned int ih=1;ih<7;++ih) {
272 if(_n_qq[ih]->effNumEntries()==0.) continue;
273 unsigned int ix= ih==1 ? 1 : 2, iy= ih==1 ? 1 : ih-1;
274 // ratio
275 std::ostringstream title;
276 title << "/TMP/n_" << ix << "_" << iy;
277 Scatter1D sTemp(title.str());
278 Scatter1DPtr s1d = registerAO(sTemp);
279 divide(_n_qq[ih],_n_qq[0],s1d);
280 // hist for axis
281 Scatter2D temphisto(refData(ix, 1, iy));
282 Scatter2DPtr ratio;
283 book(ratio, ix, 1, iy);
284 for (size_t b = 0; b < temphisto.numPoints(); b++) {
285 const double x = temphisto.point(b).x();
286 pair<double,double> ex = temphisto.point(b).xErrs();
287 pair<double,double> ex2 = ex;
288 if(ex2.first ==0.) ex2. first=0.0001;
289 if(ex2.second==0.) ex2.second=0.0001;
290 if (inRange(sqrtS(), x-ex2.first, x+ex2.second)) {
291 ratio->addPoint(x,s1d->points()[0].x(),ex,s1d->points()[0].xErrs());
292 }
293 else {
294 ratio->addPoint(x, 0., ex, make_pair(0.,.0));
295 }
296 }
297 }
298 }
299 // WW
300 if( _WW &&
301 (_n_WW[0][0]->effNumEntries()!=0. ||
302 _n_WW[1][0]->effNumEntries()!=0. )) {
303 Scatter1D ratios[2];
304 unsigned int iloc = isCompatibleWithSqrtS(183.) ? 1 : 0;
305 for(unsigned int iw=0;iw<2;++iw) {
306 if(_n_WW[iw][0]->effNumEntries()==0.) continue;
307 // scale distributions
308 scale(_h_p_charged[iw] , 1./ *_n_WW[iw][0]);
309 scale(_h_xi_charged[iw], 1./ *_n_WW[iw][0]);
310 scale(_h_pT_charged[iw], 1./ *_n_WW[iw][0]);
311 scale(_h_p_chargedB[iw] , 1./ *_n_WW[iw][0]);
312 scale(_h_xi_chargedB[iw], 1./ *_n_WW[iw][0]);
313 scale(_h_pT_chargedB[iw], 1./ *_n_WW[iw][0]);
314 for(unsigned int iy=0;iy<4;++iy) {
315 if(_h_xi_ident[iw][iy]) {
316 scale(_h_xi_ident[iw][iy], 1./ *_n_WW[iw][0]);
317 if(iw==0)
318 scale(_h_xi_ident[2][iy], 1./ *_n_WW[iw][0]);
319 }
320 }
321 // charge mult
322 ratios[iw] = *_n_WW[iw][1]/ *_n_WW[iw][0];
323 // sort out identified particle multiplicities
324 if(iloc==0) {
325 for(unsigned int iy=2;iy<5;++iy) {
326 Scatter1D R = *_n_WW[iw][iy]/ *_n_WW[iw][0];
327 Scatter2D temphisto(refData(4, 1, 2*(iy-1)-iw));
328 Scatter2DPtr mult;
329 book(mult, 4, 1, 2*(iy-1)-iw);
330 for (size_t b = 0; b < temphisto.numPoints(); b++) {
331 const double x = temphisto.point(b).x();
332 pair<double,double> ex = temphisto.point(b).xErrs();
333 pair<double,double> ex2 = ex;
334 if(ex2.first ==0.) ex2. first=0.0001;
335 if(ex2.second==0.) ex2.second=0.0001;
336 if (inRange(sqrtS()/GeV, x-ex2.first, x+ex2.second)) {
337 mult->addPoint(x, R.point(0).x(), ex, R.point(0).xErrs());
338 }
339 else {
340 mult->addPoint(x, 0., ex, make_pair(0.,.0));
341 }
342 }
343 }
344 }
345 }
346 // charged mults
347 Scatter2D temphisto(refData(3, 1, 1));
348 Scatter2DPtr mult;
349 book(mult,3,1,1);
350 for (size_t b = 0; b < temphisto.numPoints(); b++) {
351 const double x = temphisto.point(b).x();
352 pair<double,double> ex = temphisto.point(b).xErrs();
353 if((iloc==1 && b==0) || (iloc==0 && b==1))
354 mult->addPoint(x, ratios[1].point(0).x(), ex, ratios[1].point(0).xErrs());
355 else if((iloc==1 && b==2) || (iloc==0 && b==3))
356 mult->addPoint(x, ratios[0].point(0).x(), ex, ratios[0].point(0).xErrs());
357 else
358 mult->addPoint(x, 0., ex, make_pair(0.,.0));
359 }
360 // difference histos
361 // momentum
362 Scatter2D htemp = mkScatter(*(_h_p_chargedB[0]));
363 htemp.scaleY(2.);
364 Scatter2DPtr hnew;
365 book(hnew,6+2*iloc,1,1);
366 *hnew = YODA::subtract(*(_h_p_chargedB[1]),htemp);
367 hnew->setPath("/"+name()+"/"+mkAxisCode(6+2*iloc,1,1));
368 // xi
369 Scatter2D htemp2 = mkScatter(*(_h_xi_chargedB[0]));
370 htemp2.scaleY(2.);
371 Scatter2DPtr hnew2;
372 book(hnew2,9+iloc,1,3);
373 *hnew2 = YODA::subtract(*(_h_xi_chargedB[1]),htemp2);
374 hnew2->setPath("/"+name()+"/"+mkAxisCode(9+iloc,1,3));
375 // pt
376 Scatter2D hytemp3 = mkScatter(*(_h_pT_chargedB[0]));
377 hytemp3.scaleY(2.);
378 Scatter2DPtr hnew3;
379 book(hnew3,11+iloc,1,3);
380 *hnew3 = YODA::subtract(*(_h_pT_chargedB[1]),hytemp3);
381 hnew3->setPath("/"+name()+"/"+mkAxisCode(11+iloc,1,3));
382 }
383 }
384
385 ///@}
386
387
388 /// @name Histograms
389 ///@{
390 // qqbar
391 Histo1DPtr _h_qq_K0,_h_qq_Lam;
392 CounterPtr _n_qq[7];
393 // WW
394 Histo1DPtr _h_p_charged [2],_h_xi_charged [2],_h_pT_charged [2],_h_xi_ident[3][4];
395 Histo1DPtr _h_p_chargedB[2],_h_xi_chargedB[2],_h_pT_chargedB[2];
396 CounterPtr _n_WW[2][5];
397 bool _WW;
398 ///@}
399
400
401 };
402
403
404 RIVET_DECLARE_PLUGIN(DELPHI_2001_I526164);
405
406}
|