Rivet analyses referenceBELLE_2020_I1777678Single and dihadron scaled momenta spectra at 10.58 GeVExperiment: BELLE (KEKB) Inspire ID: 1777678 Status: VALIDATED Authors:
Beam energies: (5.3, 5.3) GeV Run details:
Measurement of single and di-hadron spectra by the BELLE collaboration at 10.58 GeV Source code: BELLE_2020_I1777678.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/ChargedFinalState.hh"
4#include "Rivet/Projections/FinalState.hh"
5#include "Rivet/Projections/Thrust.hh"
6
7namespace Rivet {
8
9
10 /// @brief Single and di-hadron spectra
11 class BELLE_2020_I1777678 : public Analysis {
12 public:
13
14 /// Constructor
15 RIVET_DEFAULT_ANALYSIS_CTOR(BELLE_2020_I1777678);
16
17
18 /// @name Analysis methods
19 /// @{
20
21 /// Book histograms and initialise projections before the run
22 void init() {
23
24 // Initialise and register projections
25 declare(ChargedFinalState(Cuts::abspid==211 or Cuts::abspid==321 or Cuts::abspid==2212),"CFS");
26 // projections
27 FinalState fs;
28 declare(fs,"FS");
29 declare(Thrust(fs),"Thrust");
30 // single particle hists
31 vector<int> pdg={211,321,2212};
32 for (size_t ix=0; ix<3; ++ix) {
33 book(_s_all [pdg[ix]], 1, ix+1, 1);
34 book(_s_strong[pdg[ix]], 1, ix+1, 2);
35 }
36 // dihadron histograms
37 const vector<double> edges{0.20,0.25,0.30,0.35,0.40,0.45,0.50,0.55,0.60,
38 0.65,0.70,0.75,0.80,0.85,0.90,0.95,1.00};
39 size_t i0=1;
40 for (size_t defn=0; defn<3; ++defn) {
41 for (size_t hemi=0; hemi<3; ++hemi) {
42 for (size_t ip=0; ip<6; ++ip) {
43 ++i0;
44 size_t ymax=16;
45 if(i0==7 || i0==19) ymax=15;
46 else if(i0>=8 && i0<=12) ymax=14;
47 else if(i0==13) ymax=13;
48 else if(i0==26) ymax=10;
49 else if(i0==27||i0==30) ymax= 9;
50 else if(i0==28) ymax= 7;
51 else if(i0==29) ymax= 8;
52 else if(i0==31||i0==44) ymax= 6;
53 else if(i0==45||i0==48) ymax= 5;
54 else if(i0==46||i0==47) ymax= 4;
55 else if(i0==49 ) ymax= 3;
56 book(_d_all[ip][defn][hemi], edges);
57 book(_d_strong[ip][defn][hemi], edges);
58 for (size_t iy=1; iy < _d_all[ip][defn][hemi]->numBins()+1; ++iy) {
59 if (iy <= ymax) {
60 book(_d_all[ip][defn][hemi]->bin(iy), i0, 1, iy);
61 book(_d_strong[ip][defn][hemi]->bin(iy), i0, 2, iy);
62 }
63 else {
64 _d_all[ip][defn][hemi]->maskBin(iy);
65 _d_strong[ip][defn][hemi]->maskBin(iy);
66 }
67
68 }
69 }
70 }
71 }
72 }
73
74 bool isWeak(const Particle & p) {
75 bool weak = false;
76 if(p.parents().empty()) return weak;
77 Particle parent = p.parents()[0];
78 while (!parent.parents().empty()) {
79 if(parent.abspid()==411 || parent.abspid()==421 || parent.abspid()==431 ||
80 parent.abspid()==4122 || parent.abspid()==4232 || parent.abspid()==4132 ||
81 parent.abspid()==4332) {
82 weak=true;
83 break;
84 }
85 parent = parent.parents()[0];
86 }
87 return weak;
88 }
89
90 void fillHistos(int ip,bool strong,bool same,bool opp,
91 const Particle & p1, const Particle & p2) {
92 for (size_t def=0; def<3; ++def) {
93 double z1 = 0., z2 = 0.;
94 if(def==0) {
95 z1 = 2.*p1.momentum().t()/sqrtS();
96 z2 = 2.*p2.momentum().t()/sqrtS();
97 }
98 else if(def==1) {
99 z1 = 2.*p1.momentum().t()/sqrtS();
100 z2 = (p1.momentum()*p2.momentum())/p1.momentum().t()/sqrtS();
101 }
102 else if(def==2) {
103 double p1p2 = p1.momentum()*p2.momentum();
104 double p1q = p1.momentum().t()*sqrtS();
105 double p2q = p2.momentum().t()*sqrtS();
106 z1 = (p1p2-p1.mass2()*p2.mass2()/p1p2)/(p2q-p2.mass2()*p1q/p1p2);
107 z2 = (p1p2-p1.mass2()*p2.mass2()/p1p2)/(p1q-p1.mass2()*p2q/p1p2);
108 }
109 _d_all[ip][def][0]->fill(z1,z2,0.5);
110 if (strong) _d_strong[ip][def][0]->fill(z1,z2,0.5);
111 if (same) {
112 _d_all[ip][def][1]->fill(z1,z2,0.5);
113 if (strong) _d_strong[ip][def][1]->fill(z1,z2,0.5);
114 }
115 if (opp) {
116 _d_all[ip][def][2]->fill(z1,z2,0.5);
117 if (strong) _d_strong[ip][def][2]->fill(z1,z2,0.5);
118 }
119 }
120 }
121
122 /// Perform the per-event analysis
123 void analyze(const Event& event) {
124 // apply projection
125 const ChargedFinalState& cfs = apply<ChargedFinalState>(event, "CFS");
126 // fill single particle histos
127 for (const Particle& p : cfs.particles()) {
128 const double z = 2.*p.momentum().t()/sqrtS();
129 _s_all[p.abspid()]->fill(z);
130 if (!isWeak(p)) _s_strong[p.abspid()]->fill(z);
131 }
132 // get thrust
133 const Thrust thrust = apply<Thrust>(event,"Thrust");
134 ThreeVector axis = thrust.thrustAxis();
135 Particles piK = cfs.particles(Cuts::abspid==PID::KPLUS or Cuts::abspid==PID::PIPLUS);
136 for (size_t ix=0; ix<piK.size(); ++ix) {
137 double dot1 = axis.dot(piK[ix].momentum().p3());
138 bool weak1 = isWeak(piK[ix]);
139 for (size_t iy=0; iy<piK.size(); ++iy) {
140 if (ix==iy) continue;
141 double dot2 = axis.dot(piK[iy].momentum().p3());
142 bool weak2 = isWeak(piK[iy]);
143 bool strong = !weak1 && !weak2;
144 bool same = thrust.thrust()>0.8 && dot1*dot2>0.;
145 bool opp = thrust.thrust()>0.8 && dot1*dot2<0.;
146 unsigned int ip=0;
147 if (piK[ix].pid()==PID::PIPLUS) {
148 if (piK[iy].pid()==PID::PIPLUS) ip=1;
149 else if (piK[iy].pid()==PID::PIMINUS) ip=0;
150 else if (piK[iy].pid()==PID::KPLUS ) ip=3;
151 else if (piK[iy].pid()==PID::KMINUS) ip=2;
152 }
153 else if (piK[ix].pid()==PID::PIMINUS) {
154 if (piK[iy].pid()==PID::PIPLUS) ip=0;
155 else if (piK[iy].pid()==PID::PIMINUS) ip=1;
156 else if (piK[iy].pid()==PID::KPLUS ) ip=2;
157 else if (piK[iy].pid()==PID::KMINUS ) ip=3;
158 }
159 else if(piK[ix].pid()==PID::KPLUS) {
160 if(piK[iy].pid()==PID::PIPLUS) ip=3;
161 else if(piK[iy].pid()==PID::PIMINUS) ip=2;
162 else if(piK[iy].pid()==PID::KPLUS) ip=5;
163 else if(piK[iy].pid()==PID::KMINUS) ip=4;
164 }
165 else if(piK[ix].pid()==PID::KMINUS) {
166 if(piK[iy].pid()==PID::PIPLUS) ip=2;
167 else if(piK[iy].pid()==PID::PIMINUS) ip=3;
168 else if(piK[iy].pid()==PID::KPLUS) ip=4;
169 else if(piK[iy].pid()==PID::KMINUS) ip=5;
170 }
171 fillHistos(ip,strong,same,opp,piK[ix],piK[iy]);
172 }
173 }
174 }
175
176 /// Normalise histograms etc., after the run
177 void finalize() {
178
179 const double sf = crossSection()/femtobarn/sumOfWeights();
180 scale(_s_all, sf);
181 scale(_s_strong, sf);
182 for (size_t ix=0; ix<6; ++ix) {
183 for (size_t iy=0; iy<3; ++iy) {
184 scale(_d_all[ix][iy], sf);
185 divByGroupWidth(_d_all[ix][iy]);
186 scale(_d_strong[ix][iy], sf);
187 divByGroupWidth(_d_strong[ix][iy]);
188 }
189 }
190 }
191
192 /// @}
193
194
195 /// @name Histograms
196 /// @{
197 map<int,Histo1DPtr> _s_all,_s_strong;
198 Histo1DGroupPtr _d_all[6][3][3], _d_strong[6][3][3];
199 /// @}
200
201
202 };
203
204
205 RIVET_DECLARE_PLUGIN(BELLE_2020_I1777678);
206
207}
|