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 "fastjet/JetDefinition.hh"
  8
  9namespace fastjet {
 10
 11  class 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
 92        book(_h_chMult_qq, { 5.0, 5.5, 6.5, 7.5, 9.5, 13.0, 16.0, 20.0 },
 93             { "d01-x01-y01", "d01-x01-y02", "d01-x01-y03", "d02-x01-y01",
 94               "d02-x01-y02", "d03-x01-y01", "d03-x01-y02" });
 95
 96        book(_h_chFragFunc_qq, { 13., 16., 20. }, { "d05-x01-y01", "d05-x01-y02" });
 97
 98        _sumWEbin.resize(7);
 99        for (size_t i = 0; i < 7; ++i) {
100          book(_sumWEbin[i], "/TMP/sumWEbin" + to_string(i));
101        }
102      }
103    }
104
105
106    /// Perform the per-event analysis
107    void analyze(const Event& event) {
108      // gg mode
109      if (_mode==0) {
110        // find the initial gluons
111        Particles initial;
112        for (ConstGenParticlePtr p : HepMCUtils::particles(event.genEvent())) {
113          ConstGenVertexPtr pv = p->production_vertex();
114          const PdgId pid = p->pdg_id();
115          if(pid!=21) continue;
116          bool passed = false;
117          for (ConstGenParticlePtr pp : HepMCUtils::particles(pv, Relatives::PARENTS)) {
118            const PdgId ppid = abs(pp->pdg_id());
119            passed = (ppid == PID::ELECTRON || ppid == PID::HIGGS || ppid == PID::ZBOSON || ppid == PID::GAMMA);
120            if (passed) break;
121          }
122          if (passed) initial.push_back(Particle(*p));
123        }
124        if (initial.size()!=2) vetoEvent;
125        // use the direction for the event axis
126        Vector3 axis = initial[0].momentum().p3().unit();
127        // fill histograms
128        const Particles& chps = apply<FinalState>(event, "CFS").particles();
129        unsigned int nMult[2] = {0,0};
130        // distribution
131        for (const Particle& p : chps) {
132          const double xE = 2.*p.E()/sqrtS();
133          if (_h_chFragFunc_gg)  _h_chFragFunc_gg->fill(xE);
134          if (p.momentum().p3().dot(axis)>0.) {
135            ++nMult[0];
136          }
137          else {
138            ++nMult[1];
139          }
140        }
141        // multiplicities in jet
142        _h_chMult_gg->fill(nMult[0]);
143        _h_chMult_gg->fill(nMult[1]);
144        _sumW->fill();
145      }
146      // qqbar mode
147      else {
148        // cut on the number of charged particles
149        const Particles& chParticles = apply<FinalState>(event, "CFS").particles();
150        if (chParticles.size() < 5) vetoEvent;
151        // cluster the jets
152        const Particles& particles = apply<FinalState>(event, "FS").particles();
153        fastjet::JetDefinition ee_kt_def(fastjet::ee_kt_algorithm, &p_scheme);
154        PseudoJets pParticles;
155        for (Particle p : particles) {
156          PseudoJet temp = p.pseudojet();
157          if(p.fromBottom()) {
158            temp.set_user_index(5);
159          }
160          pParticles.push_back(temp);
161        }
162        fastjet::ClusterSequence cluster(pParticles, ee_kt_def);
163        // rescale energys to just keep the directions of the jets
164        // and keep track of b tags
165        PseudoJets pJets = sorted_by_E(cluster.exclusive_jets_up_to(3));
166        if (pJets.size() < 3) vetoEvent;
167        array<Vector3, 3> dirs;
168        for (int i=0; i<3; i++) {
169          dirs[i] = Vector3(pJets[i].px(),pJets[i].py(),pJets[i].pz()).unit();
170        }
171        array<bool, 3> bTagged;
172        Jets jets;
173        for (int i=0; i<3; i++) {
174          double Ejet = sqrtS()*sin(angle(dirs[(i+1)%3],dirs[(i+2)%3])) /
175            (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])));
176          jets.push_back(FourMomentum(Ejet,Ejet*dirs[i].x(),Ejet*dirs[i].y(),Ejet*dirs[i].z()));
177          bTagged[i] = false;
178          for (PseudoJet particle : pJets[i].constituents()) {
179            if (particle.user_index() > 1 and !bTagged[i]) {
180              bTagged[i] = true;
181            }
182          }
183        }
184
185        int QUARK1 = 0, QUARK2 = 1, GLUON = 2;
186
187        if (jets[QUARK2].E() > jets[QUARK1].E()) swap(QUARK1, QUARK2);
188        if (jets[GLUON].E() > jets[QUARK1].E())  swap(QUARK1,  GLUON);
189        if (!bTagged[QUARK2]) {
190          if (!bTagged[GLUON]) vetoEvent;
191          else swap(QUARK2, GLUON);
192        }
193        if (bTagged[GLUON]) vetoEvent;
194
195        // exclude collinear or soft jets
196        double k1 = jets[QUARK1].E()*min(angle(jets[QUARK1].momentum(),jets[QUARK2].momentum()),
197                 angle(jets[QUARK1].momentum(),jets[GLUON].momentum()));
198        double k2 = jets[QUARK2].E()*min(angle(jets[QUARK2].momentum(),jets[QUARK1].momentum()),
199                 angle(jets[QUARK2].momentum(),jets[GLUON].momentum()));
200        if (k1<8.0*GeV || k2<8.0*GeV) vetoEvent;
201
202        double sqg = (jets[QUARK1].momentum()+jets[GLUON].momentum()).mass2();
203        double sgq = (jets[QUARK2].momentum()+jets[GLUON].momentum()).mass2();
204        double s = (jets[QUARK1].momentum()+jets[QUARK2].momentum()+jets[GLUON].momentum()).mass2();
205
206        double Eg = 0.5*sqrt(sqg*sgq/s);
207
208        if (Eg < 5.0 || Eg > 46.) { vetoEvent; }
209        else if (Eg > 9.5) {
210          //requirements for experimental reconstructability raise as energy raises
211          if(!bTagged[QUARK1]) {
212            vetoEvent;
213          }
214        }
215
216        // all cuts applied, increment sum of weights
217        _sumWEbin[getEbin(Eg)]->fill();
218
219
220        // transform to frame with event in y-z and glue jet in z direction
221        Matrix3 glueTOz(jets[GLUON].momentum().vector3(), Vector3(0,0,1));
222        Vector3 transQuark = glueTOz*jets[QUARK2].momentum().vector3();
223        Matrix3 quarksTOyz(Vector3(transQuark.x(), transQuark.y(), 0), Vector3(0,1,0));
224
225        // work out transformation to symmetric frame
226        array<double, 3> x_cm;
227        array<double, 3> x_cm_y;
228        array<double, 3> x_cm_z;
229        array<double, 3> x_pr;
230        for (int i=0; i<3; i++) {
231          x_cm[i] = 2*jets[i].E()/sqrt(s);
232          Vector3 p_transf = quarksTOyz*glueTOz*jets[i].p3();
233          x_cm_y[i] = 2*p_transf.y()/sqrt(s);
234          x_cm_z[i] = 2*p_transf.z()/sqrt(s);
235        }
236        x_pr[GLUON] = sqrt(4*(1-x_cm[QUARK1])*(1-x_cm[QUARK2])/(3+x_cm[GLUON]));
237        x_pr[QUARK1] = x_pr[GLUON]/(1-x_cm[QUARK1]);
238        x_pr[QUARK2] = x_pr[GLUON]/(1-x_cm[QUARK2]);
239        double gamma = (x_pr[QUARK1] + x_pr[GLUON] + x_pr[QUARK2])/2;
240        double beta_z = x_pr[GLUON]/(gamma*x_cm[GLUON]) - 1;
241        double beta_y = (x_pr[QUARK2]/gamma - x_cm[QUARK2] - beta_z*x_cm_z[QUARK2])/x_cm_y[QUARK2];
242
243        LorentzTransform toSymmetric = LorentzTransform::mkObjTransformFromBeta(Vector3(0.,beta_y,beta_z)).
244          postMult(quarksTOyz*glueTOz);
245
246        FourMomentum transGlue = toSymmetric.transform(jets[GLUON].momentum());
247        double cutAngle = angle(toSymmetric.transform(jets[QUARK2].momentum()), transGlue)/2;
248
249        int nCh = 0;
250        for (const Particle& chP : chParticles ) {
251          FourMomentum pSymmFrame = toSymmetric.transform(FourMomentum(chP.p3().mod(), chP.px(), chP.py(), chP.pz()));
252          if(angle(pSymmFrame, transGlue) < cutAngle) {
253            _h_chFragFunc_qq->fill(Eg, pSymmFrame.E()*sin(cutAngle)/Eg);
254            nCh++;
255          }
256        }
257        _h_chMult_qq->fill(Eg, nCh);
258      }
259    }
260
261
262    /// Normalise histograms etc., after the run
263    void finalize() {
264      if (_mode==0) {
265        normalize(_h_chMult_gg);
266        if (_h_chFragFunc_gg)  scale(_h_chFragFunc_gg, 0.5/(*_sumW));
267      }
268      else {
269        normalize(_h_chMult_qq);
270        for (auto& b : _h_chFragFunc_qq->bins()) {
271          const double sf = _sumWEbin[b.index()+4]->val();
272          if (sf)  scale(b, 1./sf);
273        }
274      }
275    }
276
277    /// @}
278
279
280  private:
281
282    int getEbin(double E_glue) {
283      int ih = -1;
284      if (inRange(E_glue/GeV, 5.0, 5.5)) {
285        ih = 0;
286      } else if (inRange(E_glue/GeV, 5.5, 6.5)) {
287        ih = 1;
288      } else if (inRange(E_glue/GeV, 6.5, 7.5)) {
289        ih = 2;
290      } else if (inRange(E_glue/GeV, 7.5, 9.5)) {
291        ih = 3;
292      } else if (inRange(E_glue/GeV, 9.5, 13.0)) {
293        ih = 4;
294      } else if (inRange(E_glue/GeV, 13.0, 16.0)) {
295        ih = 5;
296      } else if (inRange(E_glue/GeV, 16.0, 20.0)) {
297        ih = 6;
298      }
299      assert(ih >= 0);
300      return ih;
301    }
302
303
304  private:
305
306    // The mode
307    unsigned int _mode;
308
309    /// needed to normalize as normalised is integral =  <average no particles>
310    vector<CounterPtr> _sumWEbin;
311    CounterPtr _sumW;
312
313    // p scheme jet definition
314    fastjet::P_scheme p_scheme;
315
316
317    /// @name Histograms
318    /// @{
319    BinnedHistoPtr<int> _h_chMult_gg;
320    Histo1DPtr _h_chFragFunc_gg;
321    HistoGroupPtr<double,int> _h_chMult_qq;
322    Histo1DGroupPtr _h_chFragFunc_qq;
323    /// @}
324
325  };
326
327
328  RIVET_DECLARE_PLUGIN(OPAL_2004_I631361);
329
330}