Rivet analyses referenceOPAL_2004_I631361Gluon jet charged multiplicities and fragmentation functionsExperiment: OPAL (LEP) Inspire ID: 631361 Status: VALIDATED Authors:
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:
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}
|