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["p_K"] ,3*_iHist+3,1,1);
43 book(_h["p_p"] ,3*_iHist+4,1,1);
44
45 tribook("pi", 3*_iHist+ 8, 1, 1);
46 tribook("K", 3*_iHist+ 9, 1, 1);
47 tribook("p", 3*_iHist+10, 1, 1);
48
49 if (_iHist) {
50 _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});
51 _axes["K"] = YODA::Axis<double>({0.4, 0.5, 0.575, 0.7, 0.9, 1.1});
52 _axes["p"] = YODA::Axis<double>({0.5, 0.7, 1.2325, 2.0975});
53 _axes["r"] = YODA::Axis<double>({0.4, 0.5, 0.575, 0.7, 0.9, 1.1});
54 }
55 else {
56 _axes["pi"] = YODA::Axis<double>({0.3, 0.4, 0.5, 0.6, 1.0, 1.6});
57 _axes["K"] = YODA::Axis<double>({0.4, 0.5, 0.6, 1.0});
58 _axes["p"] = YODA::Axis<double>({0.5, 0.7, 0.9, 2.43});
59 _axes["r"] = YODA::Axis<double>({0.4, 0.5, 0.6, 1.0});
60 }
61 _axes["rp"] = YODA::Axis<double>({0.475, 0.725, 1.0, 2.2});
62 }
63
64 void tribook(const string& label, unsigned int d, unsigned int x, unsigned int y) {
65 book(_h["n_"+label], "TMP/n_"+label, refData<YODA::BinnedEstimate<string>>(d, x, y));
66 book(_h["d_"+label], "TMP/d_"+label, refData<YODA::BinnedEstimate<string>>(d, x, y));
67 book(_r[label], d, x, y);
68 }
69
70
71 /// Perform the per-event analysis
72 void analyze(const Event& event) {
73 if (_edges.empty()) {
74 for (const auto& item : _h) {
75 _edges[item.first] = item.second->xEdges();
76 }
77 for (const auto& item : _r) {
78 _edges["r"+item.first] = _h["n_"+item.first]->xEdges();
79 }
80 }
81 // First, veto on leptonic events by requiring at least 4 charged FS particles
82 const ChargedFinalState& fs = apply<ChargedFinalState>(event, "FS");
83 const size_t numParticles = fs.particles().size();
84
85 // Even if we only generate hadronic events, we still need a cut on numCharged >= 2.
86 if (numParticles < 2) {
87 MSG_DEBUG("Failed leptonic event cut");
88 vetoEvent;
89 }
90 MSG_DEBUG("Passed leptonic event cut");
91
92 // Get beams and average beam momentum
93 const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
94 const double meanBeamMom = ( beams.first.p3().mod() +
95 beams.second.p3().mod() ) / 2.0;
96 MSG_DEBUG("Avg beam momentum = " << meanBeamMom);
97
98 for (const Particle& p : fs.particles()) {
99 double modp = p.p3().mod();
100 fillND("d_pi", modp);
101 fillND("d_K", modp);
102 fillND("d_p", modp);
103 if (abs(p.pid())==211) {
104 fillhist("p_pi", modp);
105 fillND("n_pi", modp);
106 }
107 else if (abs(p.pid())==321) {
108 fillhist("p_K", modp);
109 fillND("n_K", modp);
110 }
111 else if (abs(p.pid())==2212) {
112 fillhist("p_p", modp);
113 fillND("n_p", modp);
114 }
115 }
116 }
117
118 void fillhist(const string& label, const double value, const double weight = 1.0) {
119 string edge = "OTHER";
120 const string tag = label.substr(2);
121 const size_t idx = _axes[tag].index(value);
122 if (idx && idx <= _edges[label].size()) edge = _edges[label][idx-1];
123 _h[label]->fill(edge, weight);
124 }
125
126 void fillND(const string& label, const double value) {
127 string edge = "OTHER";
128 const string tag = label.substr(2);
129 const size_t idx = _axes[(tag == "p")? "rp" : "r"].index(value);
130 if (idx && idx <= _edges["r"+tag].size()) edge = _edges["r"+tag][idx-1];
131 _h[label]->fill(edge);
132 }
133
134 /// Normalise histograms etc., after the run
135 void finalize() {
136
137 scale(_h["p_pi"], crossSection()/nanobarn/sumOfWeights());
138 scale(_h["p_K"], crossSection()/nanobarn/sumOfWeights());
139 scale(_h["p_p"], crossSection()/nanobarn/sumOfWeights());
140 vector<string> s1={"pi","K","p"};
141 for ( auto & key : s1 ) {
142 for(auto & b: _h["p_"+key]->bins()) {
143 const size_t idx = b.index();
144 b.scaleW(1./_axes[key].width(idx));
145 }
146 }
147 for (auto& item : _r) {
148 divide(_h["n_"+item.first], _h["d_"+item.first], item.second);
149 }
150
151 }
152
153 /// @}
154
155
156 /// @name Histograms
157 /// @{
158 map<string,BinnedHistoPtr<string>> _h;
159 map<string,BinnedEstimatePtr<string>> _r;
160 map<string, YODA::Axis<double>> _axes;
161 map<string, vector<string>> _edges;
162 int _iHist;
163 double sqs;
164 /// @}
165
166
167 };
168
169
170 RIVET_DECLARE_PLUGIN(TASSO_1980_I153656);
171
172
173}
|