Rivet analyses referenceCELLO_1982_I12010Energy-Energy correlation at 22 and 34 GeVExperiment: CELLO (Petra) Inspire ID: 12010 Status: VALIDATED Authors:
Beam energies: (11.0, 11.0); (17.0, 17.0) GeV Run details:
Measurement of the energy-energy correlation, and its assymetry, at centre-of-mass energies 22 and 34 GeV Beam energy must be specified as analysis option "ENERGY" when rivet-merging samples. Source code: CELLO_1982_I12010.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/Beam.hh"
4#include "Rivet/Projections/FinalState.hh"
5
6namespace Rivet {
7
8
9 /// @brief Add a short analysis description here
10 class CELLO_1982_I12010 : public Analysis {
11 public:
12
13 /// Constructor
14 RIVET_DEFAULT_ANALYSIS_CTOR(CELLO_1982_I12010);
15
16
17 /// @name Analysis methods
18 /// @{
19
20 /// Book histograms and initialise projections before the run
21 void init() {
22 // Initialise and register projections
23 declare(FinalState(), "FS");
24
25 // Book histograms
26 unsigned int iloc(0);
27 if (isCompatibleWithSqrtS(22*GeV)) {
28 iloc=1;
29 }
30 else if (isCompatibleWithSqrtS(34*GeV)) {
31 iloc=2;
32 }
33 else {
34 MSG_ERROR("Beam energy not supported!");
35 }
36 book(_histEEC, 1, 1, iloc);
37 book(_histAEEC, 3, 1, iloc);
38 book(_weightSum, "TMP/weightSum");
39 }
40
41
42 /// Perform the per-event analysis
43 void analyze(const Event& event) {
44
45 if (edges[0].empty()) edges[0] = _histEEC->xEdges();
46 if (edges[1].empty()) edges[1] = _histAEEC->xEdges();
47
48 // First, veto on leptonic events by requiring at least 4 charged FS particles
49 const FinalState& fs = apply<FinalState>(event, "FS");
50 // Even if we only generate hadronic events, we still need a cut on numCharged >= 2.
51 if ( fs.particles().size() < 2) {
52 MSG_DEBUG("Failed leptonic event cut");
53 vetoEvent;
54 }
55 MSG_DEBUG("Passed leptonic event cut");
56 _weightSum->fill();
57
58 double Evis = 0.0;
59 for (const Particle& p : fs.particles()) {
60 Evis += p.E();
61 }
62 double Evis2 = sqr(Evis);
63 // (A)EEC
64 // Need iterators since second loop starts at current outer loop iterator, i.e. no "foreach" here!
65 for (Particles::const_iterator p_i = fs.particles().begin(); p_i != fs.particles().end(); ++p_i) {
66 for (Particles::const_iterator p_j = p_i; p_j != fs.particles().end(); ++p_j) {
67 const Vector3 mom3_i = p_i->momentum().p3();
68 const Vector3 mom3_j = p_j->momentum().p3();
69 const double energy_i = p_i->momentum().E();
70 const double energy_j = p_j->momentum().E();
71 const double thetaij = mom3_i.unit().angle(mom3_j.unit());
72 double eec = (energy_i*energy_j) / Evis2;
73 if (p_i != p_j) eec *= 2.;
74 _histEEC->fill(map2string(thetaij, 0), eec);
75 if (thetaij < 0.5*M_PI) {
76 _histAEEC->fill(map2string(thetaij, 1), -eec);
77 }
78 else {
79 _histAEEC->fill(map2string(M_PI-thetaij, 1), eec);
80 }
81 }
82 }
83 }
84
85 string map2string(const double val, const size_t axis) const {
86 const size_t idx = axes[axis].index(val);
87 if (idx && idx <= edges[axis].size()) return edges[axis][idx-1];
88 return "OTHER";
89 }
90
91 /// Normalise histograms etc., after the run
92 void finalize() {
93 scale(_histEEC, 1.0/ *_weightSum);
94 for(auto & b: _histEEC->bins()) {
95 const size_t idx = b.index();
96 b.scaleW(1./axes[0].width(idx));
97 }
98 scale(_histAEEC, 1.0/ *_weightSum);
99 for(auto & b: _histAEEC->bins()) {
100 const size_t idx = b.index();
101 b.scaleW(1./axes[1].width(idx));
102 }
103 }
104
105 /// @}
106
107
108 /// @name Histograms
109 /// @{
110 BinnedHistoPtr<string> _histEEC, _histAEEC;
111 vector<string> edges[2];
112 CounterPtr _weightSum;
113 YODA::Axis<double> axes[2] = { YODA::Axis<double>(50, 0.0, M_PI),
114 YODA::Axis<double>(24, 0.0628, 0.5*M_PI) };
115 /// @}
116
117
118 };
119
120 RIVET_DECLARE_PLUGIN(CELLO_1982_I12010);
121}
|