Rivet analyses referenceLENA_1981_I164397Charged particle multiplicities and thrust in $\Upsilon(1S,2S)$ decays and nearby continuumExperiment: LENA (DORIS) Inspire ID: 164397 Status: VALIDATED Authors:
Beams: * * Beam energies: ANY Run details:
Measurement of the charged particle multiplicities and a thrust-like variable in $\Upsilon(1S,2S)$ decays and nearby continuum. As LENA had no magnetic field the momenta of the charged particles were not measured and therefore a thrust-like variable $T^\prime$ was constructed using only the directions of the charged particles. Beam energy must be specified as analysis option "ENERGY" when rivet-merging samples. Source code: LENA_1981_I164397.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/UnstableParticles.hh"
4#include "Rivet/Projections/ChargedFinalState.hh"
5#include "Rivet/Projections/Thrust.hh"
6
7namespace Rivet {
8
9
10 /// @brief Thrust like variable at Upsilon(1s,2S)
11 class LENA_1981_I164397 : public Analysis {
12 public:
13
14 /// Constructor
15 RIVET_DEFAULT_ANALYSIS_CTOR(LENA_1981_I164397);
16
17
18 /// @name Analysis methods
19 //@{
20
21 /// Book histograms and initialise projections before the run
22 void init() {
23 declare(UnstableParticles(), "UFS");
24 declare(ChargedFinalState(), "FS");
25
26 book(_weightSum_cont, "TMP/weightSum_cont");
27 book(_weightSum_Ups1, "TMP/weightSum_Ups1");
28 book(_weightSum_Ups2, "TMP/weightSum_Ups2");
29 book(_charge_cont, "TMP/charge_cont");
30 book(_charge_Ups1, "TMP/charge_Ups1");
31 book(_charge_Ups2, "TMP/charge_Ups2");
32
33 if(isCompatibleWithSqrtS(9.5149,1e-2)) {
34 book(_hist_T_cont ,4, 1, 1);
35 }
36 else if(isCompatibleWithSqrtS(9.9903,1e-2)) {
37 book(_hist_T_cont ,4, 1, 2);
38 }
39 book(_hist_T_Ups1 ,4, 1, 3);
40 book(_hist_T_Ups2 ,4, 1, 4);
41 }
42
43 /// Recursively walk the decay tree to find the charged decay products of @a p
44 void findDecayProducts(Particle mother, Particles& charged) {
45 for(const Particle & p: mother.children()) {
46 const int id = p.pid();
47 if(!p.children().empty())
48 findDecayProducts(p, charged);
49 else if(PID::isCharged(id))
50 charged.push_back(p);
51 }
52 }
53
54 // defn of thrust in paper used just the direction
55 double thrustPrime(const LorentzTransform & boost, const Particles & particles) {
56 vector<Vector3> vecs;
57 for(const Particle & p : particles) {
58 vecs.push_back(boost.transform(p.momentum()).p3().unit());
59 }
60 Thrust thrust;
61 thrust.calc(vecs);
62 return thrust.thrust();
63 }
64
65 /// Perform the per-event analysis
66 void analyze(const Event& event) {
67 // Find the Upsilons among the unstables
68 const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
69 Particles upsilons = ufs.particles(Cuts::pid==553 or Cuts::pid==100553);
70 if (upsilons.empty()) {
71 MSG_DEBUG("No Upsilons found => continuum event");
72 _weightSum_cont->fill();
73 Particles cfs = apply<ChargedFinalState>(event, "FS").particles();
74 _charge_cont->fill(cfs.size());
75 if(_hist_T_cont) {
76 LorentzTransform boost;
77 _hist_T_cont->fill(thrustPrime(boost,cfs));
78 }
79 }
80 // Upsilon(s) found
81 else {
82 for (const Particle& ups : upsilons) {
83 const int parentId = ups.pid();
84 Particles charged;
85 // boost to rest frame (if required)
86 LorentzTransform boost;
87 if (ups.p3().mod() > 1*MeV)
88 boost = LorentzTransform::mkFrameTransformFromBeta(ups.momentum().betaVec());
89 // Find the decay products we want
90 findDecayProducts(ups, charged);
91 if(parentId==553) {
92 _weightSum_Ups1->fill();
93 _charge_Ups1->fill(charged.size());
94 _hist_T_Ups1->fill(thrustPrime(boost,charged));
95 }
96 else {
97 _weightSum_Ups2->fill();
98 _charge_Ups2->fill(charged.size());
99 _hist_T_Ups2->fill(thrustPrime(boost,charged));
100 }
101 }
102 }
103
104 }
105
106
107 /// Normalise histograms etc., after the run
108 void finalize() {
109 // charged particle multiplicity
110 if(_weightSum_cont->val()>0. ) {
111 scale(_charge_cont,1./ *_weightSum_cont );
112 if(_hist_T_cont) scale(_hist_T_cont,1./ *_weightSum_cont );
113 }
114 if(_weightSum_Ups1->val()>0. ) {
115 scale(_charge_Ups1,1./ *_weightSum_Ups1 );
116 scale(_hist_T_Ups1,1./ *_weightSum_Ups1 );
117 }
118 if(_weightSum_Ups2->val()>0. ) {
119 scale(_charge_Ups2,1./ *_weightSum_Ups2 );
120 scale(_hist_T_Ups2,1./ *_weightSum_Ups2 );
121 }
122 Scatter2D tempScat(refData(3, 1, 1));
123 Scatter2DPtr _mult;
124 book(_mult, 3, 1, 1);
125 for (size_t b = 0; b < tempScat.numPoints(); b++) {
126 const double x = tempScat.point(b).x();
127 pair<double,double> ex = tempScat.point(b).xErrs();
128 pair<double,double> ex2 = ex;
129 if(ex2.first ==0.) ex2. first=0.02;
130 if(ex2.second==0.) ex2.second=0.02;
131 // Upsilon 1S
132 if(b==3) {
133 if (_weightSum_Ups1->val()>0.) {
134 _mult->addPoint(x, _charge_Ups1->val(), ex, make_pair(_charge_Ups1->err(),_charge_Ups1->err()));
135 }
136 else {
137 _mult->addPoint(x, 0., ex, make_pair(0.,.0));
138 }
139 }
140 // Upsilon 2S
141 else if(b==6) {
142 if (_weightSum_Ups2->val()>0.) {
143 _mult->addPoint(x, _charge_Ups2->val(), ex, make_pair(_charge_Ups2->err(),_charge_Ups2->err()));
144 }
145 else {
146 _mult->addPoint(x, 0., ex, make_pair(0.,.0));
147 }
148 }
149 else {
150 if (inRange(sqrtS()/GeV, x-ex2.first, x+ex2.second) && _weightSum_cont->val()>0.) {
151 _mult->addPoint(x, _charge_cont->val(), ex, make_pair(_charge_cont->err(),_charge_cont->err()));
152 }
153 else {
154 _mult->addPoint(x, 0., ex, make_pair(0.,.0));
155 }
156 }
157 }
158 }
159
160 //@}
161
162
163 /// @name Histograms
164 //@{
165 CounterPtr _weightSum_cont, _weightSum_Ups1, _weightSum_Ups2;
166 CounterPtr _charge_cont, _charge_Ups1, _charge_Ups2;
167 Histo1DPtr _hist_T_cont,_hist_T_Ups1,_hist_T_Ups2;
168 //@}
169
170
171 };
172
173
174 RIVET_DECLARE_PLUGIN(LENA_1981_I164397);
175
176}
|