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 }
185 ++iboson;
186 }
187 // calculate thrust
188 Thrust thrust;
189 thrust.calc(particles);
190 // type of event
191 // loop over particles and fill histos
192 for (const Particle & p : particles) {
193 if (!PID::isCharged(p.pid())) continue;
194 // Get momentum of each particle.
195 const Vector3 mom3 = p.p3();
196 // Scaled momenta.
197 const double mom = mom3.mod();
198 const double xiW = -log(2.*mom/mW);
199 const double xiL = -log(2.*mom/sqrtS());
200 // Get momenta components w.r.t. thrust
201 const double pTinT = dot(mom3, thrust.thrustMajorAxis());
202 const double pToutT = dot(mom3, thrust.thrustMinorAxis());
203 double pT = sqrt(sqr(pTinT)+sqr(pToutT));
204 // fill charged particle hists
205 _h_p_charged [itype]->fill(mom);
206 _h_xi_charged [itype]->fill(xiL);
207 _h_pT_charged [itype]->fill(pT );
208 _h_p_chargedB [itype]->fill(mom);
209 _h_xi_chargedB[itype]->fill(xiL);
210 _h_pT_chargedB[itype]->fill(pT );
211 // and identified particles
212 if(_h_xi_ident[itype][0]) _h_xi_ident[itype][0]->fill(xiW);
213 _n_WW[itype][1]->fill();
214 if(p.abspid()==PID::PIPLUS) {
215 if(_h_xi_ident[itype][1]) _h_xi_ident[itype][1]->fill(xiW);
216 _n_WW[itype][2]->fill();
217 }
218 else if(p.abspid()==PID::KPLUS) {
219 if(_h_xi_ident[itype][2]) _h_xi_ident[itype][2]->fill(xiW);
220 _n_WW[itype][3]->fill();
221 }
222 else if (p.abspid()==PID::PROTON){
223 if (_h_xi_ident[itype][3]) _h_xi_ident[itype][3]->fill(xiW);
224 _n_WW[itype][4]->fill();
225 }
226 }
227 // boosted specta in W rest frame
228 if (itype==0) {
229 iboson=0;
230 for (auto& W : Wbosons) {
231 ++iboson;
232 if (leptonic[iboson-1]) continue;
233 // boost to rest frame
234 LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(W.first.momentum().betaVec());
235 FourMomentum psum;
236 for (const Particle& p :W.second) psum+=p.momentum();
237 // spectra
238 for (const Particle& p : W.second) {
239 if (!PID::isCharged(p.pid())) continue;
240 // Get momentum of each particle.
241 FourMomentum phad = boost.transform(p.momentum());
242 const Vector3 mom3 = phad.p3();
243 // Scaled momenta.
244 const double mom = mom3.mod();
245 const double scaledMom = 2.*mom/mW;
246 const double xi = -log(scaledMom);
247 // /identified particle spectra
248 if (_h_xi_ident[2][0]) _h_xi_ident[2][0]->fill(xi);
249 if (p.abspid()==PID::PIPLUS) {
250 if(_h_xi_ident[2][1]) _h_xi_ident[2][1]->fill(xi);
251 }
252 else if (p.abspid()==PID::KPLUS) {
253 if (_h_xi_ident[2][2]) _h_xi_ident[2][2]->fill(xi);
254 }
255 else if (p.abspid()==PID::PROTON){
256 if (_h_xi_ident[2][3]) _h_xi_ident[2][3]->fill(xi);
257 }
258 }
259 }
260 }
261 }
262 }
263
264
265 /// Normalise histograms etc., after the run
266 void finalize() {
267 // qq
268 if (_n_qq[0]->effNumEntries()!=0.) {
269 // scale spectra if needed
270 if (_h_qq_K0 ) scale(_h_qq_K0 , 1./ *_n_qq[0]);
271 if (_h_qq_Lam) scale(_h_qq_Lam, 1./ *_n_qq[0]);
272 for (unsigned int ih=1; ih<7; ++ih) {
273 if (_n_qq[ih]->effNumEntries()==0.) continue;
274 unsigned int ix= ih==1 ? 1 : 2, iy= ih==1 ? 1 : ih-1;
275 // hist for axis
276 BinnedEstimatePtr<int> ratio;
277 book(ratio, ix, 1, iy);
278 for (auto& b : ratio->bins()) {
279 if ( isCompatibleWithSqrtS(b.xEdge()) ) {
280 b = *_n_qq[ih] / *_n_qq[0];
281 }
282 }
283 }
284 }
285 // WW
286 if( _WW && (_n_WW[0][0]->effNumEntries()!=0. || _n_WW[1][0]->effNumEntries()!=0. )) {
287 Estimate0D ratios[2];
288 unsigned int iloc = isCompatibleWithSqrtS(183.) ? 1 : 0;
289 for(unsigned int iw=0;iw<2;++iw) {
290 if(_n_WW[iw][0]->effNumEntries()==0.) continue;
291 // scale distributions
292 scale(_h_p_charged[iw] , 1./ *_n_WW[iw][0]);
293 scale(_h_xi_charged[iw], 1./ *_n_WW[iw][0]);
294 scale(_h_pT_charged[iw], 1./ *_n_WW[iw][0]);
295 scale(_h_p_chargedB[iw] , 1./ *_n_WW[iw][0]);
296 scale(_h_xi_chargedB[iw], 1./ *_n_WW[iw][0]);
297 scale(_h_pT_chargedB[iw], 1./ *_n_WW[iw][0]);
298 for(unsigned int iy=0;iy<4;++iy) {
299 if(_h_xi_ident[iw][iy]) {
300 scale(_h_xi_ident[iw][iy], 1./ *_n_WW[iw][0]);
301 if(iw==0)
302 scale(_h_xi_ident[2][iy], 1./ *_n_WW[iw][0]);
303 }
304 }
305 // charge mult
306 ratios[iw] = *_n_WW[iw][1]/ *_n_WW[iw][0];
307 // sort out identified particle multiplicities
308 if(iloc==0) {
309 for(unsigned int iy=2;iy<5;++iy) {
310 Estimate0D R = *_n_WW[iw][iy]/ *_n_WW[iw][0];
311 Estimate1DPtr mult;
312 book(mult, 4, 1, 2*(iy-1)-iw);
313 for (auto& b : mult->bins()) {
314 if ( isCompatibleWithSqrtS(b.xMid()) ) {
315 b.set(R.val(), R.errPos());
316 }
317 }
318 }
319 }
320 }
321 // charged mults
322 Estimate1DPtr mult;
323 book(mult,3,1,1);
324 for (auto& b : mult->bins()) {
325 if ((iloc==1 && b.index()==1) || (iloc==0 && b.index()==2))
326 b.set(ratios[1].val(), ratios[1].errPos());
327 else if ((iloc==1 && b.index()==3) || (iloc==0 && b.index()==4))
328 b.set(ratios[0].val(), ratios[0].errPos());
329 }
330 // difference histos
331 // momentum
332 YODA::Estimate1D htemp = _h_p_chargedB[0]->mkEstimate();
333 htemp.scale(2.);
334 Estimate1DPtr hnew;
335 book(hnew,6+2*iloc,1,1);
336 *hnew = *_h_p_chargedB[1] - htemp;
337 hnew->setPath("/"+name()+"/"+mkAxisCode(6+2*iloc,1,1));
338 // xi
339 YODA::Estimate1D htemp2 = _h_xi_chargedB[0]->mkEstimate();
340 htemp2.scale(2.);
341 Estimate1DPtr hnew2;
342 book(hnew2,9+iloc,1,3);
343 *hnew2 = *_h_xi_chargedB[1] - htemp2;
344 hnew2->setPath("/"+name()+"/"+mkAxisCode(9+iloc,1,3));
345 // pt
346 YODA::Estimate1D hytemp3 = _h_pT_chargedB[0]->mkEstimate();
347 hytemp3.scale(2.);
348 Estimate1DPtr hnew3;
349 book(hnew3,11+iloc,1,3);
350 *hnew3 = *_h_pT_chargedB[1] - hytemp3;
351 hnew3->setPath("/"+name()+"/"+mkAxisCode(11+iloc,1,3));
352 }
353 }
354
355 ///@}
356
357
358 /// @name Histograms
359 ///@{
360 // qqbar
361 Histo1DPtr _h_qq_K0,_h_qq_Lam;
362 CounterPtr _n_qq[7];
363 // WW
364 Histo1DPtr _h_p_charged [2],_h_xi_charged [2],_h_pT_charged [2],_h_xi_ident[3][4];
365 Histo1DPtr _h_p_chargedB[2],_h_xi_chargedB[2],_h_pT_chargedB[2];
366 CounterPtr _n_WW[2][5];
367 bool _WW;
368 ///@}
369
370
371 };
372
373
374 RIVET_DECLARE_PLUGIN(DELPHI_2001_I526164);
375
376}
|