Rivet analyses referencePLUTO_1983_I191161Measurement of average transverse momentum w.r.t thurst and Sphericity Axes for $E_{\text{CMS}}=7.7\to31.6$ GeVExperiment: PLUTO (DORIS/Petra) Inspire ID: 191161 Status: VALIDATED Authors:
Beams: e- e+ Beam energies: (3.9, 3.9); (4.7, 4.7); (6.0, 6.0); (6.5, 6.5); (8.5, 8.5); (11.0, 11.0); (13.8, 13.8); (15.4, 15.4) GeV Run details:
Measurements of the average transverse momentum w.r.t thurst and Sphericity Axes for $E_{\text{CMS}}=7.7\to31.6$ GeV. Beam energy must be specified as analysis option "ENERGY" when rivet-merging samples. Source code: PLUTO_1983_I191161.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/Beam.hh"
4#include "Rivet/Projections/Thrust.hh"
5#include "Rivet/Projections/Sphericity.hh"
6#include "Rivet/Projections/FinalState.hh"
7#include "Rivet/Projections/ChargedFinalState.hh"
8
9namespace Rivet {
10
11
12 /// @brief average pT w.r.t. thrust and sphericity axes
13 class PLUTO_1983_I191161 : public Analysis {
14 public:
15
16 /// Constructor
17 RIVET_DEFAULT_ANALYSIS_CTOR(PLUTO_1983_I191161);
18
19
20 /// @name Analysis methods
21 /// @{
22
23 /// Book histograms and initialise projections before the run
24 void init() {
25
26 // Initialise and register projections
27 const FinalState fs;
28 declare(fs, "FS");
29 declare(ChargedFinalState(), "CFS");
30 // Thrust
31 const Thrust thrust(fs);
32 declare(thrust, "Thrust");
33 const Sphericity sphericity(fs);
34 declare(sphericity, "Sphericity");
35
36 sqs = 1.0;
37 if (isCompatibleWithSqrtS(7.7)) sqs = 7.7;
38 else if (isCompatibleWithSqrtS(9.4)) sqs = 9.4;
39 else if (isCompatibleWithSqrtS(12.)) sqs = 12.;
40 else if (isCompatibleWithSqrtS(13.)) sqs = 13.;
41 else if (isCompatibleWithSqrtS(17.)) sqs = 17.;
42 else if (isCompatibleWithSqrtS(22.)) sqs = 22.;
43 else if (isCompatibleWithSqrtS(27.6)) sqs = 27.6;
44 else if (isCompatibleWithSqrtS(30.8)) sqs = 30.8;
45 else MSG_ERROR("Beam energy " << sqrtS() << " not supported!");
46 }
47
48
49 /// Perform the per-event analysis
50 void analyze(const Event& event) {
51 // at least 4 charged particles
52 if (apply<ChargedFinalState>(event, "CFS").particles().size()<4) vetoEvent;
53 // Sphericities
54 MSG_DEBUG("Calculating sphericity");
55 const Sphericity& sphericity = apply<Sphericity>(event, "Sphericity");
56 MSG_DEBUG("Calculating thrust");
57 const Thrust& thrust = apply<Thrust>(event, "Thrust");
58 const FinalState & fs = apply<FinalState>(event, "FS");
59 // remove pi0->gammagamma decay products and replace with the pi0s
60 // needed to get the average pTs right
61 Particles fsParticles;
62 set<ConstGenParticlePtr> pi0;
63 for(const Particle & p : fs.particles()) {
64 if ((p.pid()!=PID::PHOTON && p.abspid()!=PID::ELECTRON)|| p.parents().empty() || p.parents()[0].pid()!=PID::PI0)
65 fsParticles.push_back(p);
66 else {
67 if (pi0.find(p.parents()[0].genParticle())==pi0.end())
68 fsParticles.push_back(p.parents()[0]);
69 }
70 }
71 double pT_T_sum(0.),pT2_T_sum(0.);
72 double pT_S_sum(0.),pT2_S_sum(0.);
73 for(const Particle & p : fsParticles) {
74 const Vector3 mom3 = p.p3();
75 const double pTinT = dot(mom3, thrust.thrustMajorAxis());
76 const double pToutT = dot(mom3, thrust.thrustMinorAxis());
77 const double pTinS = dot(mom3, sphericity.sphericityMajorAxis());
78 const double pToutS = dot(mom3, sphericity.sphericityMinorAxis());
79 const double pT2_T = sqr(pTinT) + sqr(pToutT);
80 const double pT2_S = sqr(pTinS) + sqr(pToutS);
81 const double pT_T = sqrt(pT2_T);
82 const double pT_S = sqrt(pT2_S);
83 pT_T_sum += sqrt(pT2_T);
84 pT2_T_sum += pT2_T ;
85 pT_S_sum += sqrt(pT2_S);
86 pT2_S_sum += pT2_S ;
87 _p_thrust_pt .fill(pT_T /MeV );
88 _p_thrust_pt2 .fill(pT2_T/1e3/sqr(MeV));
89 _p_sphere_pt .fill(pT_S /MeV );
90 _p_sphere_pt2 .fill(pT2_S/1e3/sqr(MeV));
91 }
92 _p_thrust_sum_pt .fill(pT_T_sum /GeV );
93 _p_thrust_sum_pt2 .fill(pT2_T_sum/GeV );
94 _p_sphere_sum_pt .fill(pT_S_sum /GeV );
95 _p_sphere_sum_pt2 .fill(pT2_S_sum/GeV );
96 }
97
98
99 /// Normalise histograms etc., after the run
100 void finalize() {
101 for(unsigned int ix=1;ix<3;++ix) {
102 for(unsigned int iy=1;iy<5;++iy) {
103 double value = 0.0, error = 0.0;
104 if (ix==1) {
105 if (iy==1) {
106 value = _p_thrust_pt.xMean();
107 error = _p_thrust_pt.xStdErr();
108 }
109 else if (iy==2) {
110 value = _p_thrust_pt2.xMean();
111 error = _p_thrust_pt2.xStdErr();
112 }
113 else if (iy==3) {
114 value = _p_thrust_sum_pt.xMean();
115 error = _p_thrust_sum_pt.xStdErr();
116 }
117 else if (iy==4) {
118 value = _p_thrust_sum_pt2.xMean();
119 error = _p_thrust_sum_pt2.xStdErr();
120 }
121 }
122 else {
123 if (iy==1) {
124 value = _p_sphere_pt.xMean();
125 error = _p_sphere_pt.xStdErr();
126 }
127 else if (iy==2) {
128 value = _p_sphere_pt2.xMean();
129 error = _p_sphere_pt2.xStdErr();
130 }
131 else if (iy==3) {
132 value = _p_sphere_sum_pt.xMean();
133 error = _p_sphere_sum_pt.xStdErr();
134 }
135 else if (iy==4) {
136 value = _p_sphere_sum_pt2.xMean();
137 error = _p_sphere_sum_pt2.xStdErr();
138 }
139 }
140 Scatter2D temphisto(refData(ix, 1, iy));
141 Scatter2DPtr mult;
142 book(mult, ix, 1, iy);
143 for (size_t b = 0; b < temphisto.numPoints(); b++) {
144 const double x = temphisto.point(b).x();
145 pair<double,double> ex = temphisto.point(b).xErrs();
146 pair<double,double> ex2 = ex;
147 if (ex2.first ==0.) ex2. first=0.0001;
148 if (ex2.second==0.) ex2.second=0.0001;
149 if (inRange(sqs, x-ex2.first, x+ex2.second)) {
150 mult->addPoint(x, value, ex, make_pair(error,error));
151 }
152 else {
153 mult->addPoint(x, 0., ex, make_pair(0.,.0));
154 }
155 }
156 }
157 }
158 }
159
160 /// @}
161
162
163 /// @name Histograms
164 //@{
165 YODA::Dbn1D _p_thrust_pt, _p_thrust_pt2, _p_thrust_sum_pt, _p_thrust_sum_pt2;
166 YODA::Dbn1D _p_sphere_pt, _p_sphere_pt2, _p_sphere_sum_pt, _p_sphere_sum_pt2;
167 double sqs;
168 //@}
169
170
171 };
172
173
174 // The hook for the plugin system
175 RIVET_DECLARE_PLUGIN(PLUTO_1983_I191161);
176
177
178}
|