Rivet analyses referenceBELLE_2002_I563840Prompt Charmonium production at 10.6 GeVExperiment: BELLE (KEKB) Inspire ID: 563840 Status: VALIDATED Authors:
Beam energies: (5.3, 5.3) GeV Run details:
Measurement of the spectra for prompt $J/\psi$ and $\psi(2S)$ production at 10.6 GeV by BELLE. Source code: BELLE_2002_I563840.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/UnstableParticles.hh"
4#include "Rivet/Projections/Beam.hh"
5
6
7namespace Rivet {
8
9
10 /// @brief charmonium production at 10.6 GeV
11 class BELLE_2002_I563840 : public Analysis {
12 public:
13
14 /// Constructor
15 RIVET_DEFAULT_ANALYSIS_CTOR(BELLE_2002_I563840);
16
17
18 /// @name Analysis methods
19 /// @{
20
21 /// Book histograms and initialise projections before the run
22 void init() {
23 // projections
24 declare(UnstableParticles(),"UFS");
25 // book the histograms
26 // spectra
27 book(_h_Jpsi,3,1,1);
28 book(_h_feed,3,1,2);
29 book(_h_Psi2,3,2,1);
30 // cross sections
31 book(_h_sig_JPsi_all ,1,1,1);
32 book(_h_sig_Jpsi_high,1,2,1);
33 book(_h_sig_Jpsi_feed,1,2,2);
34 book(_h_sig_Psi2_high,1,2,3);
35 // angular distributions
36 const vector<double> bins = {2.,2.6,3.4,4.9};
37 book(_h_cThetaStar, bins);
38 book(_h_cThetaH, bins);
39 for (size_t ix=1; _h_cThetaH->numBins()+1; ++ix) {
40 if (ix <= 2) {
41 const string suff = to_string(ix+1);
42 book(_h_cThetaStar->bin(ix), "/TMP/cThetaStar_"+suff, 5, -1.0, 1.0);
43 book(_h_cThetaH->bin(ix), "/TMP/cThetaH_"+suff, 5, -1.0, 1.0);
44 }
45 else {
46 book(_h_cThetaStar->bin(ix), 4, 1, 2);
47 book(_h_cThetaH->bin(ix), 4, 2, 2);
48 }
49 }
50 book(_h_cS_low, 4, 1, 1);
51 book(_h_cS_high, "/TMP/cS_high", 5, -1.0, 1.0);
52 book(_h_cH_low, 4, 2, 1);
53 book(_h_cH_high, "/TMP/cH_high", 5, -1.0, 1.0);
54 }
55
56 void findLeptons(const Particle & mother, unsigned int & nstable, Particles& lp, Particles& lm) {
57 for (const Particle &p : mother.children()) {
58 int id = p.pid();
59 if (id == 11 || id == 13) {
60 lm.push_back(p);
61 ++nstable;
62 }
63 else if (id == -11 || id==-13) {
64 lp.push_back(p);
65 ++nstable;
66 }
67 else if (!p.children().empty()) {
68 findLeptons(p,nstable,lp,lm);
69 }
70 else {
71 ++nstable;
72 }
73 }
74 }
75
76 /// Perform the per-event analysis
77 void analyze(const Event& event) {
78 for (const Particle &p : apply<UnstableParticles>("UFS",event).particles(Cuts::pid==443 or Cuts::pid==100443 )) {
79 LorentzTransform boost = cmsTransform( beams() );
80 // check if prompt (i.e. not from B decay)
81 if (p.fromBottom()) continue;
82 bool feedDown = false;
83 FourMomentum mom = boost.transform(p.momentum());
84 double pStar = mom.p3().mod();
85 if (p.pid()==443) {
86 Particle parent=p;
87 while (!parent.parents().empty()) {
88 parent=parent.parents()[0];
89 if(p.pid()==parent.pid()) continue;
90 if((parent.abspid()%1000)/10==44) {
91 feedDown=true;
92 break;
93 }
94 }
95 _h_Jpsi->fill(pStar);
96 _h_sig_JPsi_all->fill(10.6);
97 if (pStar>2.) {
98 _h_sig_Jpsi_high->fill(10.6);
99 if(feedDown) {
100 _h_feed->fill(pStar);
101 _h_sig_Jpsi_feed->fill(10.6);
102 }
103 double cThetaS = cos(mom.p3().polarAngle());
104 _h_cThetaStar->fill(pStar, cThetaS);
105 if(pStar<3.4) _h_cS_low ->fill(cThetaS);
106 else _h_cS_high->fill(cThetaS);
107 // leptons from J/psi decay
108 unsigned int nstable = 0;
109 Particles lp, lm;
110 findLeptons(p,nstable,lp,lm);
111 if (nstable==2&&lp.size()==1&&lm.size()==1) {
112 FourMomentum pl = boost.transform(lp[0].momentum());
113 LorentzTransform b2 = LorentzTransform::mkFrameTransformFromBeta(p.momentum().betaVec());
114 pl = b2.transform(pl);
115 double cThetaH = pl.p3().unit().dot(p.p3().unit());
116 _h_cThetaH->fill(pStar, cThetaH);
117 if (pStar<3.4) _h_cH_low ->fill(cThetaH);
118 else _h_cH_high->fill(cThetaH);
119 }
120 }
121 }
122 else {
123 _h_Psi2->fill(pStar);
124 if (pStar>2.) _h_sig_Psi2_high->fill(10.6);
125 }
126 }
127 }
128
129 pair<double,pair<double,double> > calcAlpha(Histo1DPtr hist) {
130 if (hist->numEntries()==0.) return make_pair(0.,make_pair(0.,0.));
131 double sum1(0.),sum2(0.),sum3(0.),sum4(0.),sum5(0.);
132 for (const auto& bin : hist->bins() ) {
133 double Oi = bin.sumW();
134 if(Oi==0.) continue;
135 double a = 1.5*(bin.xMax() - bin.xMin());
136 double b = 0.5*(pow(bin.xMax(),3) - pow(bin.xMin(),3));
137 double Ei = bin.errW();
138 sum1 += a*Oi/sqr(Ei);
139 sum2 += b*Oi/sqr(Ei);
140 sum3 += sqr(a)/sqr(Ei);
141 sum4 += sqr(b)/sqr(Ei);
142 sum5 += a*b/sqr(Ei);
143 }
144 // calculate alpha
145 double alpha = (-3*sum1 + 9*sum2 + sum3 - 3*sum5)/(sum1 - 3*sum2 + 3*sum4 - sum5);
146 // and error
147 double cc = -pow((sum3 + 9*sum4 - 6*sum5),3);
148 double bb = -2*sqr(sum3 + 9*sum4 - 6*sum5)*(sum1 - 3*sum2 + 3*sum4 - sum5);
149 double aa = sqr(sum1 - 3*sum2 + 3*sum4 - sum5)*(-sum3 - 9*sum4 + sqr(sum1 - 3*sum2 + 3*sum4 - sum5) + 6*sum5);
150 double dis = sqr(bb)-4.*aa*cc;
151 if (dis>0.) {
152 dis = sqrt(dis);
153 return make_pair(alpha,make_pair(0.5*(-bb+dis)/aa,-0.5*(-bb-dis)/aa));
154 }
155 else {
156 return make_pair(alpha,make_pair(0.,0.));
157 }
158 }
159
160 /// Normalise histograms etc., after the run
161 void finalize() {
162 // spectra
163 normalize(_h_Jpsi,1.,false);
164 normalize(_h_feed,1.,false);
165 normalize(_h_Psi2,1.,false);
166 // cross sections
167 double fact = 1./sumOfWeights()*crossSection()/picobarn;
168 scale(_h_sig_JPsi_all ,fact);
169 scale(_h_sig_Jpsi_high,fact);
170 scale(_h_sig_Jpsi_feed,fact);
171 scale(_h_sig_Psi2_high,fact);
172 // angular dists and parameters from them
173 vector<double> bins = {2.,2.6,3.4,4.9};
174 Estimate1DPtr _h_A;
175 book(_h_A ,2,1,1);
176 Estimate1DPtr _h_alpha;
177 book(_h_alpha,2,1,2);
178 for (size_t ix=1; _h_cThetaH->numBins()+1; ++ix) {
179 normalize(_h_cThetaStar->bin(ix));
180 pair<double,pair<double,double> > alpha = calcAlpha(_h_cThetaStar->bin(ix));
181 _h_A->bin(1).set(alpha.first, alpha.second);
182 normalize(_h_cThetaH->bin(ix));
183 alpha = calcAlpha(_h_cThetaH->bin(ix));
184 _h_alpha->bin(1).set(alpha.first, alpha.second);
185 }
186 normalize(_h_cS_low);
187 pair<double,pair<double,double> > alpha = calcAlpha(_h_cS_low);
188 book(_h_A,2,2,1);
189 _h_A->bin(1).set(alpha.first, alpha.second);
190 normalize(_h_cS_high);
191 alpha = calcAlpha(_h_cS_high);
192 book(_h_A,2,3,1);
193 _h_A->bin(1).set(alpha.first, alpha.second);
194 normalize(_h_cH_low);
195 alpha = calcAlpha(_h_cH_low);
196 book(_h_alpha,2,2,2);
197 _h_alpha->bin(1).set(alpha.first, alpha.second);
198 normalize(_h_cH_high);
199 alpha = calcAlpha(_h_cH_high);
200 book(_h_alpha,2,3,2);
201 _h_alpha->bin(1).set(alpha.first, alpha.second);
202 }
203
204 /// @}
205
206
207 /// @name Histograms
208 /// @{
209 Histo1DPtr _h_Jpsi,_h_Psi2,_h_feed;
210 Histo1DPtr _h_sig_JPsi_all,_h_sig_Jpsi_high,_h_sig_Jpsi_feed,_h_sig_Psi2_high;
211 Histo1DPtr _h_cS_low,_h_cS_high,_h_cH_low,_h_cH_high;
212 Histo1DGroupPtr _h_cThetaStar, _h_cThetaH;
213
214
215 /// @}
216
217
218 };
219
220
221 RIVET_DECLARE_PLUGIN(BELLE_2002_I563840);
222
223}
|