rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

DELPHI_2001_I526164

Charged and Identified particle spectra in $q\bar{q}$ and $WW$ events
Experiment: DELPHI (LEP)
Inspire ID: 526164
Status: VALIDATED
Authors:
  • Peter Richardson
References:
  • Eur.Phys.J.C 18 (2000) 203-228, 2000
Beams: e+ e-
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:
  • e+ e- -> q qbar amd W+W- events

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}