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(unsigned int 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 double bins[17]={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 unsigned int i0=1;
40 for(unsigned int defn=0;defn<3;++defn) {
41 for(unsigned int hemi=0;hemi<3;++hemi) {
42 for(unsigned int ip=0;ip<6;++ip) {
43 i0+=1;
44 unsigned int 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 for(unsigned int iy=0;iy<ymax;++iy) {
57 Histo1DPtr temp;
58 book(temp,i0,1,iy+1);
59 _d_all [ip][defn][hemi].add(bins[iy],bins[iy+1],temp);
60 book(temp,i0,2,iy+1);
61 _d_strong[ip][defn][hemi].add(bins[iy],bins[iy+1],temp);
62 }
63 }
64 }
65 }
66 }
67
68 bool isWeak(const Particle & p) {
69 bool weak = false;
70 if(p.parents().empty()) return weak;
71 Particle parent = p.parents()[0];
72 while (!parent.parents().empty()) {
73 if(parent.abspid()==411 || parent.abspid()==421 || parent.abspid()==431 ||
74 parent.abspid()==4122 || parent.abspid()==4232 || parent.abspid()==4132 ||
75 parent.abspid()==4332) {
76 weak=true;
77 break;
78 }
79 parent = parent.parents()[0];
80 }
81 return weak;
82 }
83
84 void fillHistos(int ip,bool strong,bool same,bool opp,
85 const Particle & p1, const Particle & p2) {
86 for(unsigned int def=0;def<3;++def) {
87 double z1 = 0., z2 = 0.;
88 if(def==0) {
89 z1 = 2.*p1.momentum().t()/sqrtS();
90 z2 = 2.*p2.momentum().t()/sqrtS();
91 }
92 else if(def==1) {
93 z1 = 2.*p1.momentum().t()/sqrtS();
94 z2 = (p1.momentum()*p2.momentum())/p1.momentum().t()/sqrtS();
95 }
96 else if(def==2) {
97 double p1p2 = p1.momentum()*p2.momentum();
98 double p1q = p1.momentum().t()*sqrtS();
99 double p2q = p2.momentum().t()*sqrtS();
100 z1 = (p1p2-p1.mass2()*p2.mass2()/p1p2)/(p2q-p2.mass2()*p1q/p1p2);
101 z2 = (p1p2-p1.mass2()*p2.mass2()/p1p2)/(p1q-p1.mass2()*p2q/p1p2);
102 }
103 _d_all[ip][def][0].fill(z1,z2,0.5);
104 if(strong) _d_strong[ip][def][0].fill(z1,z2,0.5);
105 if(same) {
106 _d_all[ip][def][1].fill(z1,z2,0.5);
107 if(strong) _d_strong[ip][def][1].fill(z1,z2,0.5);
108 }
109 if(opp) {
110 _d_all[ip][def][2].fill(z1,z2,0.5);
111 if(strong) _d_strong[ip][def][2].fill(z1,z2,0.5);
112 }
113 }
114 }
115
116 /// Perform the per-event analysis
117 void analyze(const Event& event) {
118 // apply projection
119 const ChargedFinalState& cfs = apply<ChargedFinalState>(event, "CFS");
120 // fill single particle histos
121 for (const Particle& p : cfs.particles()) {
122 const double z = 2.*p.momentum().t()/sqrtS();
123 _s_all [p.abspid()]->fill(z);
124 if(!isWeak(p)) _s_strong[p.abspid()]->fill(z);
125 }
126 // get thrust
127 const Thrust thrust = apply<Thrust>(event,"Thrust");
128 ThreeVector axis = thrust.thrustAxis();
129 Particles piK = cfs.particles(Cuts::abspid==PID::KPLUS or
130 Cuts::abspid==PID::PIPLUS);
131 for(unsigned int ix=0;ix<piK.size();++ix) {
132 double dot1 = axis.dot(piK[ix].momentum().p3());
133 bool weak1 = isWeak(piK[ix]);
134 for(unsigned int iy=0;iy<piK.size();++iy) {
135 if(ix==iy) continue;
136 double dot2 = axis.dot(piK[iy].momentum().p3());
137 bool weak2 = isWeak(piK[iy]);
138 bool strong = !weak1 && !weak2;
139 bool same = thrust.thrust()>0.8 && dot1*dot2>0.;
140 bool opp = thrust.thrust()>0.8 && dot1*dot2<0.;
141 unsigned int ip=0;
142 if(piK[ix].pid()==PID::PIPLUS) {
143 if(piK[iy].pid()==PID::PIPLUS) ip=1;
144 else if(piK[iy].pid()==PID::PIMINUS) ip=0;
145 else if(piK[iy].pid()==PID::KPLUS ) ip=3;
146 else if(piK[iy].pid()==PID::KMINUS) ip=2;
147 }
148 else if(piK[ix].pid()==PID::PIMINUS) {
149 if(piK[iy].pid()==PID::PIPLUS) ip=0;
150 else if(piK[iy].pid()==PID::PIMINUS) ip=1;
151 else if(piK[iy].pid()==PID::KPLUS ) ip=2;
152 else if(piK[iy].pid()==PID::KMINUS ) ip=3;
153 }
154 else if(piK[ix].pid()==PID::KPLUS) {
155 if(piK[iy].pid()==PID::PIPLUS) ip=3;
156 else if(piK[iy].pid()==PID::PIMINUS) ip=2;
157 else if(piK[iy].pid()==PID::KPLUS) ip=5;
158 else if(piK[iy].pid()==PID::KMINUS) ip=4;
159 }
160 else if(piK[ix].pid()==PID::KMINUS) {
161 if(piK[iy].pid()==PID::PIPLUS) ip=2;
162 else if(piK[iy].pid()==PID::PIMINUS) ip=3;
163 else if(piK[iy].pid()==PID::KPLUS) ip=4;
164 else if(piK[iy].pid()==PID::KMINUS) ip=5;
165 }
166 fillHistos(ip,strong,same,opp,piK[ix],piK[iy]);
167 }
168 }
169 }
170
171 /// Normalise histograms etc., after the run
172 void finalize() {
173 for (const auto& kv : _s_all)
174 scale(kv.second,crossSection()/femtobarn/sumOfWeights());
175 for (const auto& kv : _s_strong)
176 scale(kv.second,crossSection()/femtobarn/sumOfWeights());
177 for(unsigned int ix=0;ix<6;++ix) {
178 for(unsigned int iy=0;iy<3;++iy) {
179 for(unsigned int iz=0;iz<3;++iz) {
180 _d_all [ix][iy][iz].scale(crossSection()/femtobarn/sumOfWeights(),this);
181 _d_strong[ix][iy][iz].scale(crossSection()/femtobarn/sumOfWeights(),this);
182 }
183 }
184 }
185 }
186
187 ///@}
188
189
190 /// @name Histograms
191 ///@{
192 map<int,Histo1DPtr> _s_all,_s_strong;
193 BinnedHistogram _d_all[6][3][3],_d_strong[6][3][3];
194 ///@}
195
196
197 };
198
199
200 RIVET_DECLARE_PLUGIN(BELLE_2020_I1777678);
201
202}
|