Rivet analyses referenceTPC_1988_I262143Charged particle spectra at 29 GeVExperiment: TPC (PEP) Inspire ID: 262143 Status: VALIDATED Authors:
Beam energies: (14.5, 14.5) GeV Run details:
Momentum spectra for all charged particles, $\pi^\pm$, $K^\pm$ $p,\bar{p}$ in $e^+e^-$ collisions at 29 GeV measured by the TPC experiment. Source code: TPC_1988_I262143.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/ChargedFinalState.hh"
4#include "Rivet/Projections/Beam.hh"
5
6namespace Rivet {
7
8
9 /// @brief Add a short analysis description here
10 class TPC_1988_I262143 : public Analysis {
11 public:
12
13 /// Constructor
14 RIVET_DEFAULT_ANALYSIS_CTOR(TPC_1988_I262143);
15
16
17 /// @name Analysis methods
18 /// @{
19
20 /// Book histograms and initialise projections before the run
21 void init() {
22 declare(Beam(), "Beams");
23 declare(ChargedFinalState(), "FS");
24
25 // Book histograms
26 book(_h["z_pi"], 1, 1, 1);
27 book(_h["z_K"], 1, 1, 2);
28 book(_h["z_p"], 1, 1, 3);
29 book(_h["z_all"], 1, 1, 4);
30
31 book(_h["z2_pi"], 5, 1, 1);
32 book(_h["z2_K"], 5, 1, 2);
33 book(_h["z2_p"], 5, 1, 3);
34
35 tribook("_pi", 6, 1, 1);
36 tribook("_K", 6, 1, 2);
37 tribook("_p", 6, 1, 3);
38 tribook("2_K", 7, 1, 1);
39 tribook("2_p", 7, 1, 2);
40 tribook("3_p", 7, 1, 3);
41 }
42
43 void tribook(const string& label, unsigned int d, unsigned int x, unsigned int y) {
44 book(_nd["n_"+label], "TMP/n" + label, refData<YODA::BinnedEstimate<string>>(d, x, y));
45 book(_nd["d_"+label], "TMP/d" + label, refData<YODA::BinnedEstimate<string>>(d, x, y));
46 book(_ratio[label], d, x, y);
47 }
48
49 /// Perform the per-event analysis
50 void analyze(const Event& event) {
51
52 if (_edges.empty()) {
53 for (const auto& item : _ratio) {
54 _edges[item.first] = _nd["n_" + item.first]->xEdges();
55 }
56 }
57
58 // First, veto on leptonic events by requiring at least 4 charged FS particles
59 const FinalState& fs = apply<FinalState>(event, "FS");
60 const size_t numParticles = fs.particles().size();
61
62 // Even if we only generate hadronic events, we still need a cut on numCharged >= 2.
63 if (numParticles < 2) {
64 MSG_DEBUG("Failed leptonic event cut");
65 vetoEvent;
66 }
67 MSG_DEBUG("Passed leptonic event cut");
68
69 // Get beams and average beam momentum
70 const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
71 const double meanBeamMom = (beams.first.p3().mod() + beams.second.p3().mod()) / 2.0;
72 MSG_DEBUG("Avg beam momentum = " << meanBeamMom);
73 for (const Particle& p : fs.particles()) {
74 double xP = p.p3().mod()/meanBeamMom;
75 _h["z_all"]->fill(xP);
76 histfill("d_pi", xP);
77 histfill("d_K", xP);
78 histfill("d_p", xP);
79 int id = abs(p.pid());
80 if (id==211) {
81 _h["z_pi"]->fill(xP);
82 _h["z2_pi"]->fill(xP, xP);
83 histfill("n_pi", xP, 100.);
84 histfill("d2_K", xP);
85 histfill("d2_p", xP);
86 histfill("d3_p", xP);
87 }
88 else if (id==321) {
89 _h["z_K"]->fill(xP);
90 _h["z2_K"]->fill(xP, xP);
91 histfill("n_K", xP, 100.);
92 histfill("n2_K", xP);
93 histfill("d3_p", xP);
94 }
95 else if (id==2212) {
96 _h["z_p"]->fill(xP);
97 _h["z2_p"]->fill(xP, xP);
98 histfill("n_p", xP, 100.);
99 histfill("n2_p", xP);
100 histfill("n3_p", xP);
101 }
102 }
103 }
104
105 /// Normalise histograms etc., after the run
106 void finalize() {
107 scale(_h, 1./sumOfWeights());
108 for (auto& item : _ratio) {
109 divide(_nd["n_"+item.first], _nd["d_"+item.first], item.second);
110 }
111 }
112
113 /// @}
114
115 void histfill(const string& label, const double value, const double weight = 1.0) {
116 string edge = "OTHER";
117 const string tag = label.substr(1);
118 const size_t idx = _axis.index(value);
119 if (tag == "_pi") {
120 if (0.035 <= value && value <= 0.07) edge = _edges[label][idx-3];
121 if (0.09 <= value && value <= 0.11) edge = _edges[label][idx-7];
122 if (0.16 <= value && value <= 0.7 ) edge = _edges[label][idx-11];
123 }
124 if (tag == "_K") {
125 if (0.035 <= value && value <= 0.07) edge = _edges[label][idx-3];
126 if (0.09 <= value && value <= 0.14) edge = _edges[label][idx-7];
127 if (0.25 <= value && value <= 0.7 ) edge = _edges[label][idx-12];
128 }
129 if (tag == "_p" || tag == "2_p") {
130 if (0.035 <= value && value <= 0.11) edge = _edges[label][idx-3];
131 if (0.25 <= value && value <= 0.7 ) edge = _edges[label][idx-11];
132 }
133 if (tag == "2_K") {
134 if (0.025 <= value && value <= 0.11) edge = _edges[label][idx];
135 if (0.25 <= value && value <= 0.7 ) edge = _edges[label][idx-8];
136 }
137 if (tag == "3_p") {
138 if (0.035 <= value && value <= 0.07) edge = _edges[label][idx-3];
139 if (0.09 <= value && value <= 0.11) edge = _edges[label][idx-7];
140 if (0.25 <= value && value <= 0.7 ) edge = _edges[label][idx-15];
141 }
142 _nd[label]->fill(edge, weight);
143 }
144
145
146 /// @name Histograms
147 /// @{
148 map<string, Histo1DPtr> _h;
149 map<string, BinnedHistoPtr<string>> _nd;
150 map<string, BinnedEstimatePtr<string>> _ratio;
151 map<string, vector<string>> _edges;
152 YODA::Axis<double> _axis{0.025, 0.03, 0.035, 0.04, 0.045, 0.05, 0.055, 0.06, 0.06, 0.065,
153 0.07, 0.075, 0.08, 0.085, 0.09, 0.1, 0.11, 0.12, 0.13, 0.14,
154 0.16, 0.18, 0.2, 0.22, 0.25, 0.3, 0.35, 0.4, 0.5, 0.6, 0.7};
155
156
157 /// }
158
159
160 };
161
162
163 RIVET_DECLARE_PLUGIN(TPC_1988_I262143);
164
165
166}
|