Rivet analyses referenceTASSO_1980_I153656$\pi^\pm$, $K^\pm$ and $p,\bar{p}$ spectra in $e^+e^-$ at 12 and 30 GeVExperiment: TASSO (Petra) Inspire ID: 153656 Status: VALIDATED Authors:
Beam energies: (6.0, 6.0); (15.0, 15.0) GeV Run details:
Measurement of the $\pi^\pm$, $K^\pm$ and $p,\bar{p}$ spectra in $e^+e^-$ collisions for centre-of-mass energies of 12 and 30 GeV by the TASSO experiment at Petra. Beam energy must be specified as analysis option "ENERGY" when rivet-merging samples. Source code: TASSO_1980_I153656.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 pi, K and proton spectra at 12 and 30 GeV
10 class TASSO_1980_I153656 : public Analysis {
11 public:
12
13 /// Constructor
14 RIVET_DEFAULT_ANALYSIS_CTOR(TASSO_1980_I153656);
15
16
17 /// @name Analysis methods
18 /// @{
19
20 /// Book histograms and initialise projections before the run
21 void init() {
22
23 // Initialise and register projections
24 declare(Beam(), "Beams");
25 declare(ChargedFinalState(), "FS");
26
27 // Book histograms
28 _iHist=-1;
29 sqs = 1.;
30 if (isCompatibleWithSqrtS(12*GeV)) {
31 _iHist = 0;
32 sqs = 12.;
33 }
34 else if (isCompatibleWithSqrtS(30*GeV)) {
35 _iHist = 1;
36 sqs = 30.;
37 }
38 else
39 MSG_ERROR("Beam energy " << sqrtS() << " GeV not supported!");
40
41 book(_h["p_pi"],3*_iHist+2,1,1);
42 book(_h["x_pi"],3*_iHist+2,1,2);
43 book(_h["p_K"] ,3*_iHist+3,1,1);
44 book(_h["x_K"] ,3*_iHist+3,1,2);
45 book(_h["p_p"] ,3*_iHist+4,1,1);
46 book(_h["x_p"] ,3*_iHist+4,1,2);
47
48 tribook("pi", 3*_iHist+ 8, 1, 1);
49 tribook("K", 3*_iHist+ 9, 1, 1);
50 tribook("p", 3*_iHist+10, 1, 1);
51
52 if (_iHist) {
53 _axes["pi"] = YODA::Axis<double>({0.325, 0.375, 0.425, 0.5, 0.575, 0.7, 0.9, 1.1, 1.3, 1.5});
54 _axes["K"] = YODA::Axis<double>({0.4, 0.5, 0.575, 0.7, 0.9, 1.1});
55 _axes["p"] = YODA::Axis<double>({0.5, 0.7, 1.2325, 2.0975});
56 _axes["r"] = YODA::Axis<double>({0.4, 0.5, 0.675, 0.925});
57 }
58 else {
59 _axes["pi"] = YODA::Axis<double>({0.3, 0.4, 0.5, 0.675, 1.05, 1.55});
60 _axes["K"] = YODA::Axis<double>({0.4, 0.5, 0.675, 0.925});
61 _axes["p"] = YODA::Axis<double>({0.475, 0.725, 1.2125, 1.9375});
62 _axes["r"] = YODA::Axis<double>({0.4, 0.5, 0.575, 0.7, 0.9, 1.1});
63 }
64 _axes["rp"] = YODA::Axis<double>({0.475, 0.725, 1.0, 2.2});
65 }
66
67 void tribook(const string& label, unsigned int d, unsigned int x, unsigned int y) {
68 book(_h["n_"+label], "TMP/n_"+label, refData<YODA::BinnedEstimate<string>>(d, x, y));
69 book(_h["d_"+label], "TMP/d_"+label, refData<YODA::BinnedEstimate<string>>(d, x, y));
70 book(_r[label], d, x, y);
71 }
72
73
74 /// Perform the per-event analysis
75 void analyze(const Event& event) {
76 if (_edges.empty()) {
77 for (const auto& item : _h) {
78 _edges[item.first] = item.second->xEdges();
79 }
80 for (const auto& item : _r) {
81 _edges["r"+item.first] = item.second->xEdges();
82 }
83 }
84 // First, veto on leptonic events by requiring at least 4 charged FS particles
85 const ChargedFinalState& fs = apply<ChargedFinalState>(event, "FS");
86 const size_t numParticles = fs.particles().size();
87
88 // Even if we only generate hadronic events, we still need a cut on numCharged >= 2.
89 if (numParticles < 2) {
90 MSG_DEBUG("Failed leptonic event cut");
91 vetoEvent;
92 }
93 MSG_DEBUG("Passed leptonic event cut");
94
95 // Get beams and average beam momentum
96 const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
97 const double meanBeamMom = ( beams.first.p3().mod() +
98 beams.second.p3().mod() ) / 2.0;
99 MSG_DEBUG("Avg beam momentum = " << meanBeamMom);
100
101 for (const Particle& p : fs.particles()) {
102 double modp = p.p3().mod();
103 fillND("d_pi", modp);
104 fillND("d_K", modp);
105 fillND("d_p", modp);
106 double beta = modp/p.E();
107 double xE = p.E()/meanBeamMom;
108 if (abs(p.pid())==211) {
109 fillhist("p_pi", modp);
110 fillhist("x_pi", xE, 1./beta);
111 fillND("n_pi", modp);
112 }
113 else if (abs(p.pid())==321) {
114 fillhist("p_K", modp);
115 fillhist("x_K", xE, 1./beta);
116 fillND("n_K", modp);
117 }
118 else if (abs(p.pid())==2212) {
119 fillhist("p_p", modp);
120 fillhist("x_p", xE, 1./beta);
121 fillND("n_p", modp);
122 }
123 }
124 }
125
126 void fillhist(const string& label, const double value, const double weight = 1.0) {
127 string edge = "OTHER";
128 const string tag = label.substr(2);
129 const size_t idx = _axes[tag].index(value);
130 if (idx && idx <= _edges[label].size()) edge = _edges[label][idx];
131 _h[label]->fill(edge, weight);
132 }
133
134 void fillND(const string& label, const double value) {
135 string edge = "OTHER";
136 const string tag = label.substr(2);
137 const size_t idx = _axes[(tag == "p")? "rp" : "r"].index(value);
138 if (idx && idx <= _edges["r"+tag].size()) edge = _edges[label][idx];
139 _h[label]->fill(edge);
140 }
141
142 /// Normalise histograms etc., after the run
143 void finalize() {
144
145 scale(_h["p_pi"], crossSection()/nanobarn/sumOfWeights());
146 scale(_h["x_pi"], sqr(sqs)*crossSection()/microbarn/sumOfWeights());
147 scale(_h["p_K"], crossSection()/nanobarn/sumOfWeights());
148 scale(_h["x_K"], sqr(sqs)*crossSection()/microbarn/sumOfWeights());
149 scale(_h["p_p"], crossSection()/nanobarn/sumOfWeights());
150 scale(_h["x_p"], sqr(sqs)*crossSection()/microbarn/sumOfWeights());
151
152 for (auto& item : _r) {
153 divide(_h["n_"+item.first], _h["d_"+item.first], item.second);
154 }
155
156 }
157
158 /// @}
159
160
161 /// @name Histograms
162 /// @{
163 map<string,BinnedHistoPtr<string>> _h;
164 map<string,BinnedEstimatePtr<string>> _r;
165 map<string, YODA::Axis<double>> _axes;
166 map<string, vector<string>> _edges;
167 int _iHist;
168 double sqs;
169 /// @}
170
171
172 };
173
174
175 RIVET_DECLARE_PLUGIN(TASSO_1980_I153656);
176
177
178}
|