Rivet analyses referenceHRS_1985_I201482Event shapes in $e^+e^-$ collisions at 29 GeVExperiment: HRS (PEP) Inspire ID: 201482 Status: VALIDATED Authors:
Beam energies: (14.5, 14.5) GeV Run details:
Measurement of a range of event shapes at 29 GeV by the HRS experiment. The event are seperate into two ($S\leq0.25$, $A\leq0.1$) and jet three($S>0.25$, $A\leq0.1$) jet regions. The mean values of event shapes are not implemented. Source code: HRS_1985_I201482.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/Beam.hh"
4#include "Rivet/Projections/ChargedFinalState.hh"
5#include "Rivet/Projections/Sphericity.hh"
6#include "Rivet/Projections/Thrust.hh"
7
8namespace Rivet {
9
10
11 /// @brief event shapes at 29 GeV
12 class HRS_1985_I201482 : public Analysis {
13 public:
14
15 /// Constructor
16 RIVET_DEFAULT_ANALYSIS_CTOR(HRS_1985_I201482);
17
18
19 /// @name Analysis methods
20 /// @{
21
22 /// Book histograms and initialise projections before the run
23 void init() {
24
25 // Initialise and register projections
26 declare(Beam(), "Beams");
27 const ChargedFinalState cfs;
28 declare(cfs, "FS");
29 declare(Sphericity(cfs), "Sphericity");
30 const Thrust thrust(cfs);
31 declare(thrust, "Thrust");
32
33 // Book histograms
34 book(_histSphericity, 1, 1, 1);
35 book(_histThrust , 3, 1, 1);
36 book(_histThrust2Jet, 4, 1, 1);
37 book(_histAplanarity, 6, 1, 1);
38 book(_histZ , 10, 1, 1);
39 book(_histZ2Jet , 11, 1, 1);
40 book(_histZScale , 12, 1, 1);
41 book(_histZJet[0] , 13, 1, 1);
42 book(_histZJet[1] , 14, 1, 1);
43 book(_histZJet[2] , 15, 1, 1);
44 book(_histXFeyn , 16, 1, 1);
45 book(_histXFeyn2Jet , 17, 1, 1);
46 book(_histRap[0] , 19, 1, 1);
47 book(_histRap[1] , 20, 1, 1);
48 book(_histPtT , 22, 1, 1);
49 book(_histPtT2Jet , 23, 1, 1);
50 book(_histPtTIn , 24, 1, 1);
51 book(_histPtTOut , 25, 1, 1);
52 book(_wSum ,"TMP/wSum");
53 book(_wSum2,"TMP/wSum2");
54
55 _axes[0] = YODA::Axis<double>({-5.0, -4.0, -3.5, -3.0, -2.75, -2.5, -2.3, -2.0, -1.75, -1.5,
56 -1.25, -1.0, -0.75, -0.5, -0.25, 0.0, 0.25, 0.5, 0.75, 1.0,
57 1.25, 1.5, 1.75, 2.0, 2.3, 2.5, 2.75, 3.0, 3.5, 4.0, 5.0});
58 _axes[1] = YODA::Axis<double>({-4.0, -3.5, -3.0, -2.75, -2.5, -2.25, -2.0, -1.75, -1.5,
59 -1.25, -1.0, -0.75, -0.5, -0.25, 0.0, 0.25, 0.5, 0.75,
60 1.0, 1.25, 1.5, 1.75, 2.0, 2.25, 2.5, 2.75, 3.0, 3.5, 4.0});
61 }
62
63
64 /// Perform the per-event analysis
65 void analyze(const Event& event) {
66 if (_edges[0].empty()) {
67 _edges[0] = _histRap[0]->xEdges();
68 _edges[1] = _histRap[1]->xEdges();
69 }
70 // require 5 charged particles
71 const FinalState& fs = apply<FinalState>(event, "FS");
72 const size_t numParticles = fs.particles().size();
73 if(numParticles<5) vetoEvent;
74 // Get beams and average beam momentum
75 const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
76 const double meanBeamMom = ( beams.first.p3().mod() +
77 beams.second.p3().mod() ) / 2.0;
78 MSG_DEBUG("Avg beam momentum = " << meanBeamMom);
79 // calc thrust and sphericity
80 const Thrust& thrust = apply<Thrust>(event, "Thrust");
81 Vector3 axis = thrust.thrustAxis();
82 const Sphericity& sphericity = apply<Sphericity>(event, "Sphericity");
83 // identify two and three jet regions
84 bool twoJet = sphericity.sphericity()<=0.25 && sphericity.aplanarity()<=0.1;
85 //bool threeJet = sphericity.sphericity() >0.25 && sphericity.aplanarity()<=0.1;
86 _wSum->fill();
87 if (twoJet) _wSum2->fill();
88 // basic event shapes
89 _histSphericity->fill(sphericity.sphericity());
90 _histThrust ->fill(thrust.thrust());
91 _histAplanarity->fill(sphericity.aplanarity());
92 if(twoJet) _histThrust2Jet->fill(thrust.thrust());
93 double pTSqIn = 0.;
94 double pTSqOut = 0.;
95 unsigned int iPlus(0),iMinus(0);
96 // single particle dists
97 for(const Particle & p : sortBy(fs.particles(),cmpMomByP)) {
98 const double z = p.p3().mod()/meanBeamMom;
99 const double momT = axis.dot(p.p3());
100 const double xF = fabs(momT)/meanBeamMom;
101 const double energy = p.E();
102 const double rap = 0.5 * std::log((energy + momT) / (energy - momT));
103 const double pTin = dot(p.p3(), thrust.thrustMajorAxis());
104 const double pTout = dot(p.p3(), thrust.thrustMinorAxis());
105 const double pT2 = sqr(pTin)+sqr(pTout);
106 pTSqIn += sqr(dot(p.p3(), sphericity.sphericityMajorAxis()));
107 pTSqOut += sqr(dot(p.p3(), sphericity.sphericityMinorAxis()));
108 _histZ->fill(z);
109 _histZScale->fill(z);
110 _histXFeyn ->fill(xF, z);
111 _histRap[0]->fill(map2string(rap, 0));
112 _histPtT->fill(pT2);
113 if(twoJet) {
114 _histZ2Jet->fill(z);
115 _histXFeyn2Jet->fill(xF, z);
116 _histRap[1]->fill(map2string(rap, 1));
117 _histPtT2Jet->fill(pT2);
118 if(momT>0.&&iPlus<3) {
119 _histZJet[iPlus]->fill(z);
120 iPlus+=1;
121 }
122 else if(momT<0.&&iMinus<3) {
123 _histZJet[iMinus]->fill(z);
124 iMinus+=1;
125 }
126 }
127 }
128 _histPtTIn ->fill(pTSqIn /numParticles);
129 _histPtTOut->fill(pTSqOut/numParticles);
130 }
131
132
133 /// Normalise histograms etc., after the run
134 void finalize() {
135
136 normalize(_histSphericity);
137 normalize(_histThrust);
138 normalize(_histThrust2Jet);
139 normalize(_histAplanarity);
140 scale(_histZ, 1./ *_wSum);
141 scale(_histZScale, sqr(sqrtS())*crossSection()/microbarn/sumOfWeights());
142 scale(_histXFeyn, 1./M_PI/ *_wSum);
143 scale(_histRap[0], 1./ *_wSum);
144 scale(_histZ2Jet, 1./ *_wSum2);
145 scale(_histXFeyn2Jet, 1./M_PI/ *_wSum2);
146 scale(_histRap[1], 1./ *_wSum2);
147 scale(_histPtT, 1./ *_wSum);
148 scale(_histPtT2Jet, 1./ *_wSum2);
149 scale(_histPtTIn, 1./ *_wSum);
150 scale(_histPtTOut, 1./ *_wSum);
151 for (size_t i=0; i<3; ++i) {
152 scale(_histZJet[i], 0.5/ *_wSum2);
153 }
154 for(unsigned int ix=0;ix<2;++ix) {
155 for(auto & b: _histRap[ix]->bins()) {
156 const size_t idx = b.index();
157 b.scaleW(1./_axes[ix].width(idx));
158 }
159 }
160 }
161
162 /// @}
163
164 string map2string(const double value, const size_t k) const {
165 const size_t idx = _axes[k].index(value);
166 if (idx && idx <= _edges[k].size()) return _edges[k][idx-1];
167 return "OTHER";
168 }
169
170
171 /// @name Histograms
172 /// @{
173 Histo1DPtr _histSphericity, _histThrust, _histThrust2Jet, _histAplanarity,
174 _histZ, _histZ2Jet, _histZScale, _histXFeyn, _histXFeyn2Jet,
175 _histPtT, _histPtT2Jet, _histPtTIn, _histPtTOut, _histZJet[3];
176 BinnedHistoPtr<string> _histRap[2];
177 CounterPtr _wSum,_wSum2;
178
179 YODA::Axis<double> _axes[2];
180 vector<string> _edges[2];
181
182 /// @}
183
184
185 };
186
187
188 RIVET_DECLARE_PLUGIN(HRS_1985_I201482);
189
190
191}
|
