rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

SLD_2004_S5693039

Production of $\pi^+$, $\pi^-$, $K^+$, $K^-$, $p$ and $\bar p$ in Light ($uds$), $c$ and $b$ Jets from Z Decays
Experiment: SLD (SLC)
Inspire ID: 630327
Status: VALIDATED
Authors:
  • Peter Richardson
References: Beams: e+ e-
Beam energies: (45.6, 45.6) GeV
Run details:
  • Hadronic Z decay events generated on the Z pole ($\sqrt{s} = 91.2$ GeV)

Measurements of the differential production rates of stable charged particles in hadronic $Z^0$ decays, and of charged pions, kaons and protons identified over a wide momentum range. In addition to flavour-inclusive $Z^0$ decays, measurements are made for $Z^0$ decays into light ($u$, $d$, $s$), $c$ and $b$ primary flavors.

Source code: SLD_2004_S5693039.cc
  1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/Beam.hh"
  4#include "Rivet/Projections/FinalState.hh"
  5#include "Rivet/Projections/ChargedFinalState.hh"
  6#include "Rivet/Projections/Thrust.hh"
  7
  8#define I_KNOW_THE_INITIAL_QUARKS_PROJECTION_IS_DODGY_BUT_NEED_TO_USE_IT
  9#include "Rivet/Projections/InitialQuarks.hh"
 10
 11namespace Rivet {
 12
 13
 14  /// @brief SLD flavour-dependent fragmentation paper
 15  ///
 16  /// @author Peter Richardson
 17  class SLD_2004_S5693039 : public Analysis {
 18  public:
 19
 20    /// Constructor
 21    RIVET_DEFAULT_ANALYSIS_CTOR(SLD_2004_S5693039);
 22
 23
 24    /// @name Analysis methods
 25    //@{
 26
 27    void analyze(const Event& e) {
 28      // First, veto on leptonic events by requiring at least 2 charged FS particles
 29      const FinalState& fs = apply<FinalState>(e, "FS");
 30      const size_t numParticles = fs.particles().size();
 31
 32      // Even if we only generate hadronic events, we still need a cut on numCharged >= 2.
 33      if (numParticles < 2) {
 34        MSG_DEBUG("Failed ncharged cut");
 35        vetoEvent;
 36      }
 37      MSG_DEBUG("Passed ncharged cut");
 38
 39      // Get beams and average beam momentum
 40      const ParticlePair& beams = apply<Beam>(e, "Beams").beams();
 41      const double meanBeamMom = ( beams.first.p3().mod() +
 42                                   beams.second.p3().mod() ) / 2.0;
 43      MSG_DEBUG("Avg beam momentum = " << meanBeamMom);
 44      int flavour = 0;
 45      const InitialQuarks& iqf = apply<InitialQuarks>(e, "IQF");
 46
 47      // If we only have two quarks (qqbar), just take the flavour.
 48      // If we have more than two quarks, look for the highest energetic q-qbar pair.
 49      Particles quarks;
 50      if (iqf.particles().size() == 2) {
 51        flavour = iqf.particles().front().abspid();
 52        quarks = iqf.particles();
 53      }
 54      else {
 55        map<int, Particle > quarkmap;
 56        for (const Particle& p : iqf.particles()) {
 57          if (quarkmap.find(p.pid())==quarkmap.end())
 58            quarkmap[p.pid()] = p;
 59          else if (quarkmap[p.pid()].E() < p.E())
 60            quarkmap[p.pid()] = p;
 61        }
 62        double maxenergy = 0.;
 63        for (int i = 1; i <= 5; ++i) {
 64          double energy(0.);
 65          if(quarkmap.find( i)!=quarkmap.end())
 66            energy += quarkmap[ i].E();
 67          if(quarkmap.find(-i)!=quarkmap.end())
 68            energy += quarkmap[-i].E();
 69          if (energy > maxenergy) flavour = i;
 70        }
 71        if(quarkmap.find( flavour)!=quarkmap.end())
 72          quarks.push_back(quarkmap[ flavour]);
 73        if(quarkmap.find(-flavour)!=quarkmap.end())
 74          quarks.push_back(quarkmap[-flavour]);
 75      }
 76
 77      // total multiplicities
 78      switch (flavour) {
 79      case PID::DQUARK:
 80      case PID::UQUARK:
 81      case PID::SQUARK:
 82        _weightLight  ->fill();
 83        _weightedTotalChargedPartNumLight  ->fill(numParticles);
 84        break;
 85      case PID::CQUARK:
 86        _weightCharm  ->fill();
 87        _weightedTotalChargedPartNumCharm  ->fill(numParticles);
 88        break;
 89      case PID::BQUARK:
 90        _weightBottom ->fill();
 91        _weightedTotalChargedPartNumBottom ->fill(numParticles);
 92        break;
 93      }
 94      // thrust axis for projections
 95      Vector3 axis = apply<Thrust>(e, "Thrust").thrustAxis();
 96      double dot(0.);
 97      if(!quarks.empty()) {
 98        dot = quarks[0].p3().dot(axis);
 99        if(quarks[0].pid()<0) dot *= -1.;
100      }
101      // spectra and individual multiplicities
102      for (const Particle& p : fs.particles()) {
103        double pcm = p.p3().mod();
104        const double xp = pcm/meanBeamMom;
105
106        // if in quark or antiquark hemisphere
107        bool quark = p.p3().dot(axis)*dot>0.;
108
109        _h_PCharged ->fill(pcm     );
110        // all charged
111        switch (flavour) {
112        case PID::DQUARK:
113        case PID::UQUARK:
114        case PID::SQUARK:
115          _h_XpChargedL->fill(xp);
116          break;
117        case PID::CQUARK:
118          _h_XpChargedC->fill(xp);
119          break;
120        case PID::BQUARK:
121          _h_XpChargedB->fill(xp);
122          break;
123        }
124
125        int id = p.abspid();
126        // charged pions
127        if (id == PID::PIPLUS) {
128          _h_XpPiPlus->fill(xp);
129          _h_XpPiPlusTotal->fill(xp);
130          switch (flavour) {
131          case PID::DQUARK:
132          case PID::UQUARK:
133          case PID::SQUARK:
134            _h_XpPiPlusL->fill(xp);
135            _h_NPiPlusL->fill(sqrtS());
136            if( ( quark && p.pid()>0 ) || ( !quark && p.pid()<0 ))
137              _h_RPiPlus->fill(xp);
138            else
139              _h_RPiMinus->fill(xp);
140            break;
141          case PID::CQUARK:
142            _h_XpPiPlusC->fill(xp);
143            _h_NPiPlusC->fill(sqrtS());
144            break;
145          case PID::BQUARK:
146            _h_XpPiPlusB->fill(xp);
147            _h_NPiPlusB->fill(sqrtS());
148            break;
149          }
150        }
151        else if (id == PID::KPLUS) {
152          _h_XpKPlus->fill(xp);
153          _h_XpKPlusTotal->fill(xp);
154          switch (flavour) {
155          case PID::DQUARK:
156          case PID::UQUARK:
157          case PID::SQUARK:
158            _h_XpKPlusL->fill(xp);
159            _h_NKPlusL->fill(sqrtS());
160            if( ( quark && p.pid()>0 ) || ( !quark && p.pid()<0 ))
161              _h_RKPlus->fill(xp);
162            else
163              _h_RKMinus->fill(xp);
164            break;
165          case PID::CQUARK:
166            _h_XpKPlusC->fill(xp);
167            _h_NKPlusC->fill(sqrtS());
168            break;
169          case PID::BQUARK:
170            _h_XpKPlusB->fill(xp);
171            _h_NKPlusB->fill(sqrtS());
172            break;
173          }
174        }
175        else if (id == PID::PROTON) {
176          _h_XpProton->fill(xp);
177          _h_XpProtonTotal->fill(xp);
178          switch (flavour) {
179          case PID::DQUARK:
180          case PID::UQUARK:
181          case PID::SQUARK:
182            _h_XpProtonL->fill(xp);
183            _h_NProtonL->fill(sqrtS());
184            if( ( quark && p.pid()>0 ) || ( !quark && p.pid()<0 ))
185              _h_RProton->fill(xp);
186            else
187              _h_RPBar  ->fill(xp);
188            break;
189          case PID::CQUARK:
190            _h_XpProtonC->fill(xp);
191            _h_NProtonC->fill(sqrtS());
192            break;
193          case PID::BQUARK:
194            _h_XpProtonB->fill(xp);
195            _h_NProtonB->fill(sqrtS());
196            break;
197          }
198        }
199      }
200    }
201
202
203    void init() {
204      // Projections
205      declare(Beam(), "Beams");
206      declare(ChargedFinalState(), "FS");
207      declare(InitialQuarks(), "IQF");
208      declare(Thrust(FinalState()), "Thrust");
209
210      // Book histograms
211      book(_h_PCharged   , 1, 1, 1);
212      book(_h_XpPiPlus   , 2, 1, 2);
213      book(_h_XpKPlus    , 3, 1, 2);
214      book(_h_XpProton   , 4, 1, 2);
215      book(_h_XpPiPlusTotal , 2, 2, 2);
216      book(_h_XpKPlusTotal  , 3, 2, 2);
217      book(_h_XpProtonTotal , 4, 2, 2);
218      book(_h_XpPiPlusL  , 5, 1, 1);
219      book(_h_XpPiPlusC  , 5, 1, 2);
220      book(_h_XpPiPlusB  , 5, 1, 3);
221      book(_h_XpKPlusL   , 6, 1, 1);
222      book(_h_XpKPlusC   , 6, 1, 2);
223      book(_h_XpKPlusB   , 6, 1, 3);
224      book(_h_XpProtonL  , 7, 1, 1);
225      book(_h_XpProtonC  , 7, 1, 2);
226      book(_h_XpProtonB  , 7, 1, 3);
227      book(_h_XpChargedL , 8, 1, 1);
228      book(_h_XpChargedC , 8, 1, 2);
229      book(_h_XpChargedB , 8, 1, 3);
230
231      book(_h_NPiPlusL  , 5, 2, 1);
232      book(_h_NPiPlusC  , 5, 2, 2);
233      book(_h_NPiPlusB  , 5, 2, 3);
234      book(_h_NKPlusL   , 6, 2, 1);
235      book(_h_NKPlusC   , 6, 2, 2);
236      book(_h_NKPlusB   , 6, 2, 3);
237      book(_h_NProtonL  , 7, 2, 1);
238      book(_h_NProtonC  , 7, 2, 2);
239      book(_h_NProtonB  , 7, 2, 3);
240
241      book(_h_RPiPlus  , 9, 1, 1);
242      book(_h_RPiMinus , 9, 1, 2);
243      book(_h_RKPlus   ,10, 1, 1);
244      book(_h_RKMinus  ,10, 1, 2);
245      book(_h_RProton  ,11, 1, 1);
246      book(_h_RPBar    ,11, 1, 2);
247
248      // Ratios: used as target of divide() later
249      book(_s_PiM_PiP,  9, 1, 3);
250      book(_s_KM_KP  , 10, 1, 3);
251      book(_s_Pr_PBar, 11, 1, 3);
252
253      book(_weightedTotalChargedPartNumLight, "_weightedTotalChargedPartNumLight");
254      book(_weightedTotalChargedPartNumCharm, "_weightedTotalChargedPartNumCharm");
255      book(_weightedTotalChargedPartNumBottom, "_weightedTotalChargedPartNumBottom");
256      book(_weightLight, "_weightLight");
257      book(_weightCharm, "_weightCharm");
258      book(_weightBottom, "_weightBottom");
259
260      book(tmp1, 8, 2, 1, true);
261      book(tmp2, 8, 2, 2, true);
262      book(tmp3, 8, 2, 3, true);
263      book(tmp4, 8, 3, 2, true);
264      book(tmp5, 8, 3, 3, true);
265
266
267    }
268
269
270    /// Finalize
271    void finalize() {
272
273      // Multiplicities
274      /// @todo Include errors
275      const double avgNumPartsLight = _weightedTotalChargedPartNumLight->val() / _weightLight->val();
276      const double avgNumPartsCharm = _weightedTotalChargedPartNumCharm->val() / _weightCharm->val();
277      const double avgNumPartsBottom = _weightedTotalChargedPartNumBottom->val() / _weightBottom->val();
278      tmp1->point(0).setY(avgNumPartsLight);
279      tmp2->point(0).setY(avgNumPartsCharm);
280      tmp3->point(0).setY(avgNumPartsBottom);
281      tmp4->point(0).setY(avgNumPartsCharm - avgNumPartsLight);
282      tmp5->point(0).setY(avgNumPartsBottom - avgNumPartsLight);
283
284      // Do divisions
285      divide(*_h_RPiMinus - *_h_RPiPlus, *_h_RPiMinus + *_h_RPiPlus, _s_PiM_PiP);
286      divide(*_h_RKMinus - *_h_RKPlus, *_h_RKMinus + *_h_RKPlus, _s_KM_KP);
287      divide(*_h_RProton - *_h_RPBar, *_h_RProton + *_h_RPBar, _s_Pr_PBar);
288
289      // Scale histograms
290      scale(_h_PCharged,      1./sumOfWeights());
291      scale(_h_XpPiPlus,      1./sumOfWeights());
292      scale(_h_XpKPlus,       1./sumOfWeights());
293      scale(_h_XpProton,      1./sumOfWeights());
294      scale(_h_XpPiPlusTotal, 1./sumOfWeights());
295      scale(_h_XpKPlusTotal,  1./sumOfWeights());
296      scale(_h_XpProtonTotal, 1./sumOfWeights());
297      scale(_h_XpPiPlusL,     1. / *_weightLight);
298      scale(_h_XpPiPlusC,     1. / *_weightCharm);
299      scale(_h_XpPiPlusB,     1. / *_weightBottom);
300      scale(_h_XpKPlusL,      1. / *_weightLight);
301      scale(_h_XpKPlusC,      1. / *_weightCharm);
302      scale(_h_XpKPlusB,      1. / *_weightBottom);
303      scale(_h_XpProtonL,     1. / *_weightLight);
304      scale(_h_XpProtonC,     1. / *_weightCharm);
305      scale(_h_XpProtonB,     1. / *_weightBottom);
306
307      scale(_h_XpChargedL, 1. / *_weightLight);
308      scale(_h_XpChargedC, 1. / *_weightCharm);
309      scale(_h_XpChargedB, 1. / *_weightBottom);
310
311      scale(_h_NPiPlusL, 1. / *_weightLight);
312      scale(_h_NPiPlusC, 1. / *_weightCharm);
313      scale(_h_NPiPlusB, 1. / *_weightBottom);
314      scale(_h_NKPlusL,  1. / *_weightLight);
315      scale(_h_NKPlusC,  1. / *_weightCharm);
316      scale(_h_NKPlusB,  1. / *_weightBottom);
317      scale(_h_NProtonL, 1. / *_weightLight);
318      scale(_h_NProtonC, 1. / *_weightCharm);
319      scale(_h_NProtonB, 1. / *_weightBottom);
320
321      // Paper suggests this should be 0.5/weight but it has to be 1.0 to get normalisations right...
322      scale(_h_RPiPlus,  1. / *_weightLight);
323      scale(_h_RPiMinus, 1. / *_weightLight);
324      scale(_h_RKPlus,   1. / *_weightLight);
325      scale(_h_RKMinus,  1. / *_weightLight);
326      scale(_h_RProton,  1. / *_weightLight);
327      scale(_h_RPBar,    1. / *_weightLight);
328
329      // convert ratio to %
330      _s_PiM_PiP->scaleY(100.);
331      _s_KM_KP  ->scaleY(100.);
332      _s_Pr_PBar->scaleY(100.);
333    }
334
335    //@}
336
337
338  private:
339
340    Scatter2DPtr tmp1, tmp2, tmp3, tmp4, tmp5;
341
342    /// Multiplicities
343    CounterPtr _weightedTotalChargedPartNumLight, _weightedTotalChargedPartNumCharm, _weightedTotalChargedPartNumBottom;
344
345    /// Weights
346    CounterPtr _weightLight, _weightCharm, _weightBottom;
347
348    // Histograms
349    /// @{
350    Histo1DPtr _h_PCharged;
351    Histo1DPtr _h_XpPiPlus, _h_XpKPlus, _h_XpProton;
352    Histo1DPtr _h_XpPiPlusTotal, _h_XpKPlusTotal, _h_XpProtonTotal;
353    Histo1DPtr _h_XpPiPlusL, _h_XpPiPlusC, _h_XpPiPlusB;
354    Histo1DPtr _h_XpKPlusL, _h_XpKPlusC, _h_XpKPlusB;
355    Histo1DPtr _h_XpProtonL, _h_XpProtonC, _h_XpProtonB;
356    Histo1DPtr _h_XpChargedL, _h_XpChargedC, _h_XpChargedB;
357    Histo1DPtr _h_NPiPlusL, _h_NPiPlusC, _h_NPiPlusB;
358    Histo1DPtr _h_NKPlusL, _h_NKPlusC, _h_NKPlusB;
359    Histo1DPtr _h_NProtonL, _h_NProtonC, _h_NProtonB;
360    Histo1DPtr _h_RPiPlus, _h_RPiMinus, _h_RKPlus;
361    Histo1DPtr _h_RKMinus, _h_RProton, _h_RPBar;
362    Scatter2DPtr _s_PiM_PiP, _s_KM_KP, _s_Pr_PBar;
363    /// @}
364
365  };
366
367
368
369  RIVET_DECLARE_ALIASED_PLUGIN(SLD_2004_S5693039, SLD_2004_I630327);
370
371}