rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

OPAL_2004_I631361

Gluon jet charged multiplicities and fragmentation functions
Experiment: OPAL (LEP)
Inspire ID: 631361
Status: VALIDATED
Authors:
  • Daniel Reichelt
References:
  • Phys. Rev. D69, 032002,2004
  • hep-ex/0310048
Beams: e+ e-
Beam energies: (5.2, 5.2); (6.0, 6.0); (7.0, 7.0); (8.4, 8.4); (10.9, 10.9); (14.2, 14.2); (17.7, 17.7); (45.6, 45.6) GeV
Run details:
  • The fictional $e^+e^-\to gg$ process

Measurement of gluon jet properties using the jet boost algorithm, a technique to select unbiased samples of gluon jets in $e^+e^-$ annihilation, i.e. gluon jets free of biases introduced by event selection or jet finding criteria. Two modes are provided, the prefer option is to produce the fictional $e^+e^-\to g g $ process to be used due to the corrections applied to the data, PROCESS=GG. The original analysis technique to extract gluon jets from hadronic $e^+e^-$ events using $e^+e^-\to q\bar{q}$ events, PROCESS=QQ, is also provided but cannot be used for tuning as the data has been corrected for impurities, however it is still useful qualitatively in order to check the properties of gluon jets in the original way in which there were measured rather than using a fictitious process.

Source code: OPAL_2004_I631361.cc
  1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/ChargedFinalState.hh"
  5#include "Rivet/Projections/FastJets.hh"
  6#include "Rivet/Projections/HadronicFinalState.hh"
  7#include "Rivet/Tools/BinnedHistogram.hh"
  8#include "fastjet/JetDefinition.hh"
  9namespace fastjet {
 10
 11class P_scheme : public JetDefinition::Recombiner {
 12 public:
 13  std::string description() const {return "";}
 14  void recombine(const PseudoJet & pa, const PseudoJet & pb,
 15                         PseudoJet & pab) const {
 16    PseudoJet tmp = pa + pb;
 17    double E = sqrt(tmp.px()*tmp.px() + tmp.py()*tmp.py() + tmp.pz()*tmp.pz());
 18    pab.reset_momentum(tmp.px(), tmp.py(), tmp.pz(), E);
 19  }
 20  void preprocess(PseudoJet & p) const {
 21    double E = sqrt(p.px()*p.px() + p.py()*p.py() + p.pz()*p.pz());
 22    p.reset_momentum(p.px(), p.py(), p.pz(), E);
 23  }
 24  ~P_scheme() { }
 25};
 26
 27}
 28
 29namespace Rivet {
 30
 31
 32  class OPAL_2004_I631361 : public Analysis {
 33  public:
 34
 35    /// Constructor
 36    RIVET_DEFAULT_ANALYSIS_CTOR(OPAL_2004_I631361);
 37
 38    /// @name Analysis methods
 39    //@{
 40
 41    void init() {
 42
 43      // Get options from the new option system
 44      _mode = 0;
 45      if ( getOption("PROCESS") == "GG" ) _mode = 0;
 46      if ( getOption("PROCESS") == "QQ" ) _mode = 1;
 47      // projections we need for both cases      
 48      const FinalState fs;
 49      declare(fs, "FS");
 50      const ChargedFinalState cfs;
 51      declare(cfs, "CFS");
 52      // additional projections for q qbar
 53      if(_mode==1) {
 54        declare(HadronicFinalState(fs), "HFS");
 55        declare(HadronicFinalState(cfs), "HCFS");
 56      }
 57
 58      // book the histograms
 59      if(_mode==0) {
 60        int ih(0), iy(0);
 61        if (inRange(0.5*sqrtS()/GeV, 5.0, 5.5)) {
 62          ih = 1;
 63          iy = 1;
 64        } else if (inRange(0.5*sqrtS()/GeV, 5.5, 6.5)) {
 65          ih = 1;
 66          iy = 2;
 67        } else if (inRange(0.5*sqrtS()/GeV, 6.5, 7.5)) {
 68          ih = 1;
 69          iy = 3;
 70        } else if (inRange(0.5*sqrtS()/GeV, 7.5, 9.5)) {
 71          ih = 2;
 72          iy = 1;
 73        } else if (inRange(0.5*sqrtS()/GeV, 9.5, 13.0)) {
 74          ih = 2;
 75          iy = 2;
 76        } else if (inRange(0.5*sqrtS()/GeV, 13.0, 16.0)) {
 77          ih = 3;
 78          iy = 1;
 79        } else if (inRange(0.5*sqrtS()/GeV, 16.0, 20.0)) {
 80          ih = 3;
 81          iy = 2;
 82        }
 83        if (!ih)  MSG_WARNING("Option \"PROCESS=GG\" not compatible with this beam energy!");
 84        assert(ih>0);
 85        book(_h_chMult_gg, ih,1,iy);
 86        if(ih==3)  book(_h_chFragFunc_gg, 5,1,iy);
 87        else       _h_chFragFunc_gg = nullptr;
 88	book(_sumW,"/TMP/sumW");
 89      }
 90      else {
 91        Histo1DPtr dummy;
 92        _h_chMult_qq.add( 5.0,  5.5, book(dummy, 1,1,1));
 93        _h_chMult_qq.add( 5.5,  6.5, book(dummy, 1,1,2));
 94        _h_chMult_qq.add( 6.5,  7.5, book(dummy, 1,1,3));
 95        _h_chMult_qq.add( 7.5,  9.5, book(dummy, 2,1,1));
 96        _h_chMult_qq.add( 9.5, 13.0, book(dummy, 2,1,2));
 97        _h_chMult_qq.add(13.0, 16.0, book(dummy, 3,1,1));
 98        _h_chMult_qq.add(16.0, 20.0, book(dummy, 3,1,2));
 99	
100        _h_chFragFunc_qq.add(13.0, 16.0, book(dummy, 5,1,1));
101        _h_chFragFunc_qq.add(16.0, 20.0, book(dummy, 5,1,2));
102
103        _sumWEbin.resize(7);
104        for (size_t i = 0; i < 7; ++i) {
105          book(_sumWEbin[i], "/TMP/sumWEbin" + to_string(i));
106        }
107      }
108    }
109
110
111    /// Perform the per-event analysis
112    void analyze(const Event& event) {
113      // gg mode
114      if(_mode==0) {
115        // find the initial gluons
116        Particles initial;
117        for (ConstGenParticlePtr p : HepMCUtils::particles(event.genEvent())) {
118          ConstGenVertexPtr pv = p->production_vertex();
119          const PdgId pid = p->pdg_id();
120          if(pid!=21) continue;
121          bool passed = false;
122          for (ConstGenParticlePtr pp : HepMCUtils::particles(pv, Relatives::PARENTS)) {
123            const PdgId ppid = abs(pp->pdg_id());
124            passed = (ppid == PID::ELECTRON || ppid == PID::HIGGS ||
125                ppid == PID::ZBOSON   || ppid == PID::GAMMA);
126            if(passed) break;
127          }
128          if(passed) initial.push_back(Particle(*p));
129        }
130        if(initial.size()!=2) vetoEvent;
131        // use the direction for the event axis
132        Vector3 axis = initial[0].momentum().p3().unit();
133        // fill histograms
134        const Particles& chps = applyProjection<FinalState>(event, "CFS").particles();
135        unsigned int nMult[2] = {0,0};
136        // distribution
137        for (const Particle& p : chps) {
138          double xE = 2.*p.E()/sqrtS();
139          if(_h_chFragFunc_gg) _h_chFragFunc_gg->fill(xE);
140          if(p.momentum().p3().dot(axis)>0.)
141            ++nMult[0];
142          else
143            ++nMult[1];
144        }
145        // multiplicities in jet
146        _h_chMult_gg->fill(nMult[0]);
147        _h_chMult_gg->fill(nMult[1]);
148	_sumW->fill();
149      }
150      // qqbar mode
151      else {
152        // cut on the number of charged particles
153        const Particles& chParticles = applyProjection<FinalState>(event, "CFS").particles();
154        if(chParticles.size() < 5) vetoEvent;  
155        // cluster the jets
156        const Particles& particles = applyProjection<FinalState>(event, "FS").particles();
157        fastjet::JetDefinition ee_kt_def(fastjet::ee_kt_algorithm, &p_scheme);
158        PseudoJets pParticles;
159        for (Particle p : particles) {
160          PseudoJet temp = p.pseudojet();
161          if(p.fromBottom()) {
162            temp.set_user_index(5);
163          }
164          pParticles.push_back(temp);
165        }
166        fastjet::ClusterSequence cluster(pParticles, ee_kt_def);
167        // rescale energys to just keep the directions of the jets
168        // and keep track of b tags
169        PseudoJets pJets = sorted_by_E(cluster.exclusive_jets_up_to(3));
170        if(pJets.size() < 3) vetoEvent;
171        array<Vector3, 3> dirs;
172        for(int i=0; i<3; i++) {
173          dirs[i] = Vector3(pJets[i].px(),pJets[i].py(),pJets[i].pz()).unit();
174        }
175        array<bool, 3> bTagged;
176        Jets jets;
177        for(int i=0; i<3; i++) {
178          double Ejet = sqrtS()*sin(angle(dirs[(i+1)%3],dirs[(i+2)%3])) /
179            (sin(angle(dirs[i],dirs[(i+1)%3])) + sin(angle(dirs[i],dirs[(i+2)%3])) + sin(angle(dirs[(i+1)%3],dirs[(i+2)%3])));
180          jets.push_back(FourMomentum(Ejet,Ejet*dirs[i].x(),Ejet*dirs[i].y(),Ejet*dirs[i].z()));
181          bTagged[i] = false;
182          for (PseudoJet particle : pJets[i].constituents()) {
183            if(particle.user_index() > 1 and !bTagged[i]) {
184              bTagged[i] = true;
185            }
186          }
187        }
188        
189        int QUARK1 = 0, QUARK2 = 1, GLUON = 2; 
190        
191        if(jets[QUARK2].E() > jets[QUARK1].E()) swap(QUARK1, QUARK2);
192        if(jets[GLUON].E() > jets[QUARK1].E())  swap(QUARK1,  GLUON);
193        if(!bTagged[QUARK2]) {
194          if(!bTagged[GLUON]) vetoEvent;
195          else swap(QUARK2, GLUON);
196        }
197        if(bTagged[GLUON]) vetoEvent;
198        
199        // exclude collinear or soft jets
200        double k1 = jets[QUARK1].E()*min(angle(jets[QUARK1].momentum(),jets[QUARK2].momentum()),
201                 angle(jets[QUARK1].momentum(),jets[GLUON].momentum())); 
202        double k2 = jets[QUARK2].E()*min(angle(jets[QUARK2].momentum(),jets[QUARK1].momentum()),
203                 angle(jets[QUARK2].momentum(),jets[GLUON].momentum()));
204        if(k1<8.0*GeV || k2<8.0*GeV) vetoEvent;
205        
206        double sqg = (jets[QUARK1].momentum()+jets[GLUON].momentum()).mass2();
207        double sgq = (jets[QUARK2].momentum()+jets[GLUON].momentum()).mass2();
208        double s = (jets[QUARK1].momentum()+jets[QUARK2].momentum()+jets[GLUON].momentum()).mass2();
209        
210        double Eg = 0.5*sqrt(sqg*sgq/s);
211        
212        if(Eg < 5.0 || Eg > 46.) { vetoEvent; }
213        else if(Eg > 9.5) {
214          //requirements for experimental reconstructability raise as energy raises
215          if(!bTagged[QUARK1]) {
216            vetoEvent;
217          }
218        }
219        
220        // all cuts applied, increment sum of weights
221        _sumWEbin[getEbin(Eg)]->fill();
222        
223        
224        // transform to frame with event in y-z and glue jet in z direction
225        Matrix3 glueTOz(jets[GLUON].momentum().vector3(), Vector3(0,0,1));
226        Vector3 transQuark = glueTOz*jets[QUARK2].momentum().vector3();
227        Matrix3 quarksTOyz(Vector3(transQuark.x(), transQuark.y(), 0), Vector3(0,1,0));
228        
229        // work out transformation to symmetric frame
230        array<double, 3> x_cm;
231        array<double, 3> x_cm_y;
232        array<double, 3> x_cm_z;
233        array<double, 3> x_pr;
234        for(int i=0; i<3; i++) {
235          x_cm[i] = 2*jets[i].E()/sqrt(s);
236          Vector3 p_transf = quarksTOyz*glueTOz*jets[i].p3();
237          x_cm_y[i] = 2*p_transf.y()/sqrt(s);
238          x_cm_z[i] = 2*p_transf.z()/sqrt(s);
239        }
240        x_pr[GLUON] = sqrt(4*(1-x_cm[QUARK1])*(1-x_cm[QUARK2])/(3+x_cm[GLUON]));
241        x_pr[QUARK1] = x_pr[GLUON]/(1-x_cm[QUARK1]);
242        x_pr[QUARK2] = x_pr[GLUON]/(1-x_cm[QUARK2]);
243        double gamma = (x_pr[QUARK1] + x_pr[GLUON] + x_pr[QUARK2])/2;
244        double beta_z = x_pr[GLUON]/(gamma*x_cm[GLUON]) - 1;
245        double beta_y = (x_pr[QUARK2]/gamma - x_cm[QUARK2] - beta_z*x_cm_z[QUARK2])/x_cm_y[QUARK2];
246        
247        LorentzTransform toSymmetric = LorentzTransform::mkObjTransformFromBeta(Vector3(0.,beta_y,beta_z)).
248          postMult(quarksTOyz*glueTOz);
249        
250        FourMomentum transGlue = toSymmetric.transform(jets[GLUON].momentum());
251        double cutAngle = angle(toSymmetric.transform(jets[QUARK2].momentum()), transGlue)/2;
252        
253        int nCh = 0;
254        for (const Particle& chP : chParticles ) {
255          FourMomentum pSymmFrame = toSymmetric.transform(FourMomentum(chP.p3().mod(), chP.px(), chP.py(), chP.pz()));
256          if(angle(pSymmFrame, transGlue) < cutAngle) {
257            _h_chFragFunc_qq.fill(Eg, pSymmFrame.E()*sin(cutAngle)/Eg);
258            nCh++;
259          }
260        }
261        _h_chMult_qq.fill(Eg, nCh);
262      }
263    }
264
265
266    /// Normalise histograms etc., after the run
267    void finalize() {
268      if(_mode==0) {
269        normalize(_h_chMult_gg);
270        if(_h_chFragFunc_gg) scale(_h_chFragFunc_gg, 0.5/(*_sumW));
271      }
272      else {
273        for (Histo1DPtr hist : _h_chMult_qq.histos()) {
274          normalize(hist);
275        }
276        for(int i=0; i<2; i++) {
277          if(!isZero(_sumWEbin[i+5]->val())) {
278            scale(_h_chFragFunc_qq.histos()[i], 1./(*_sumWEbin[i+5]));
279          }
280        }
281      }
282    }
283
284    //@}
285
286
287  private:
288
289    int getEbin(double E_glue) {
290      int ih = -1;
291      if (inRange(E_glue/GeV, 5.0, 5.5)) {
292        ih = 0;
293      } else if (inRange(E_glue/GeV, 5.5, 6.5)) {
294        ih = 1;
295      } else if (inRange(E_glue/GeV, 6.5, 7.5)) {
296        ih = 2;
297      } else if (inRange(E_glue/GeV, 7.5, 9.5)) {
298        ih = 3;
299      } else if (inRange(E_glue/GeV, 9.5, 13.0)) {
300        ih = 4;
301      } else if (inRange(E_glue/GeV, 13.0, 16.0)) {
302        ih = 5;
303      } else if (inRange(E_glue/GeV, 16.0, 20.0)) {
304        ih = 6;
305      }
306      assert(ih >= 0);
307      return ih;
308    }
309
310
311    class PScheme : public JetDefinition::Recombiner {
312     public:
313      std::string description() const {return "";}
314      void recombine(const PseudoJet & pa, const PseudoJet & pb, PseudoJet & pab) const {
315        PseudoJet tmp = pa + pb;
316        double E = sqrt(tmp.px()*tmp.px() + tmp.py()*tmp.py() + tmp.pz()*tmp.pz());
317        pab.reset_momentum(tmp.px(), tmp.py(), tmp.pz(), E);
318      }
319      void preprocess(PseudoJet & p) const {
320        double E = sqrt(p.px()*p.px() + p.py()*p.py() + p.pz()*p.pz());
321        p.reset_momentum(p.px(), p.py(), p.pz(), E);
322      }
323      ~PScheme() { }
324    };
325
326
327  private:
328
329    // The mode
330    unsigned int _mode;
331
332    /// needed to normalize as normalised is integral =  <average no particles>
333    vector<CounterPtr> _sumWEbin;
334    CounterPtr _sumW;
335    
336    // p scheme jet definition
337    fastjet::P_scheme p_scheme;
338
339
340    /// @name Histograms
341    //@{
342    Histo1DPtr _h_chMult_gg;
343    Histo1DPtr _h_chFragFunc_gg;
344    BinnedHistogram _h_chMult_qq;
345    BinnedHistogram _h_chFragFunc_qq;
346    //@}
347
348  };
349
350
351  // The hook for the plugin system
352  RIVET_DECLARE_PLUGIN(OPAL_2004_I631361);
353
354}