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 vector<double> bins = {2.,2.6,3.4,4.9};
37 for(unsigned int ix=0;ix<3;++ix) {
38 Histo1DPtr temp;
39 // book cos theta*
40 if(ix<=1) {
41 std::ostringstream title;
42 title << "/TMP/cThetaStar_" << ix+1;
43 book(temp,title.str(),5,-1.,1.);
44 }
45 else {
46 book(temp,4,1,2);
47 }
48 _h_cThetaStar.add(bins[ix],bins[ix]+1,temp);
49 // book cos theta_H
50 if(ix<=1) {
51 std::ostringstream title;
52 title << "/TMP/cThetaH_" << ix+1;
53 book(temp,title.str(),5,-1.,1.);
54 }
55 else {
56 book(temp,4,2,2);
57 }
58 _h_cThetaH.add(bins[ix],bins[ix]+1,temp);
59 }
60 book(_h_cS_low,4,1,1);
61 book(_h_cS_high,"/TMP/cS_high",5,-1.,1.);
62 book(_h_cH_low,4,2,1);
63 book(_h_cH_high,"/TMP/cH_high",5,-1.,1.);
64 }
65
66 void findLeptons(const Particle & mother,
67 unsigned int & nstable,
68 Particles& lp, Particles& lm) {
69 for(const Particle & p : mother.children()) {
70 int id = p.pid();
71 if ( id == 11 || id == 13 ) {
72 lm.push_back(p);
73 ++nstable;
74 }
75 else if (id == -11 || id==-13) {
76 lp.push_back(p);
77 ++nstable;
78 }
79 else if ( !p.children().empty() ) {
80 findLeptons(p,nstable,lp,lm);
81 }
82 else
83 ++nstable;
84 }
85 }
86
87 /// Perform the per-event analysis
88 void analyze(const Event& event) {
89 for(const Particle & p : apply<UnstableParticles>("UFS",event).particles(Cuts::pid==443 or Cuts::pid==100443 )) {
90 LorentzTransform boost = cmsTransform( beams() );
91 // check if prompt (i.e. not from B decay)
92 if(p.fromBottom()) continue;
93 bool feedDown = false;
94 FourMomentum mom = boost.transform(p.momentum());
95 double pStar = mom.p3().mod();
96 if(p.pid()==443) {
97 Particle parent=p;
98 while(!parent.parents().empty()) {
99 parent=parent.parents()[0];
100 if(p.pid()==parent.pid()) continue;
101 if((parent.abspid()%1000)/10==44) {
102 feedDown=true;
103 break;
104 }
105 }
106 _h_Jpsi->fill(pStar);
107 _h_sig_JPsi_all->fill(10.6);
108 if(pStar>2.) {
109 _h_sig_Jpsi_high->fill(10.6);
110 if(feedDown) {
111 _h_feed->fill(pStar);
112 _h_sig_Jpsi_feed->fill(10.6);
113 }
114 double cThetaS = cos(mom.p3().polarAngle());
115 _h_cThetaStar.fill(pStar,cThetaS);
116 if(pStar<3.4) _h_cS_low ->fill(cThetaS);
117 else _h_cS_high->fill(cThetaS);
118 // leptons from J/psi decay
119 unsigned int nstable = 0;
120 Particles lp, lm;
121 findLeptons(p,nstable,lp,lm);
122 if(nstable==2&&lp.size()==1&&lm.size()==1) {
123 FourMomentum pl = boost.transform(lp[0].momentum());
124 LorentzTransform b2 = LorentzTransform::mkFrameTransformFromBeta(p.momentum().betaVec());
125 pl = b2.transform(pl);
126 double cThetaH = pl.p3().unit().dot(p.p3().unit());
127 _h_cThetaH.fill(pStar,cThetaH);
128 if(pStar<3.4) _h_cH_low ->fill(cThetaH);
129 else _h_cH_high->fill(cThetaH);
130 }
131 }
132 }
133 else {
134 _h_Psi2->fill(pStar);
135 if(pStar>2.) _h_sig_Psi2_high->fill(10.6);
136 }
137 }
138 }
139
140 pair<double,pair<double,double> > calcAlpha(Histo1DPtr hist) {
141 if(hist->numEntries()==0.) return make_pair(0.,make_pair(0.,0.));
142 double sum1(0.),sum2(0.),sum3(0.),sum4(0.),sum5(0.);
143 for (auto bin : hist->bins() ) {
144 double Oi = bin.area();
145 if(Oi==0.) continue;
146 double a = 1.5*(bin.xMax() - bin.xMin());
147 double b = 0.5*(pow(bin.xMax(),3) - pow(bin.xMin(),3));
148 double Ei = bin.areaErr();
149 sum1 += a*Oi/sqr(Ei);
150 sum2 += b*Oi/sqr(Ei);
151 sum3 += sqr(a)/sqr(Ei);
152 sum4 += sqr(b)/sqr(Ei);
153 sum5 += a*b/sqr(Ei);
154 }
155 // calculate alpha
156 double alpha = (-3*sum1 + 9*sum2 + sum3 - 3*sum5)/(sum1 - 3*sum2 + 3*sum4 - sum5);
157 // and error
158 double cc = -pow((sum3 + 9*sum4 - 6*sum5),3);
159 double bb = -2*sqr(sum3 + 9*sum4 - 6*sum5)*(sum1 - 3*sum2 + 3*sum4 - sum5);
160 double aa = sqr(sum1 - 3*sum2 + 3*sum4 - sum5)*(-sum3 - 9*sum4 + sqr(sum1 - 3*sum2 + 3*sum4 - sum5) + 6*sum5);
161 double dis = sqr(bb)-4.*aa*cc;
162 if(dis>0.) {
163 dis = sqrt(dis);
164 return make_pair(alpha,make_pair(0.5*(-bb+dis)/aa,-0.5*(-bb-dis)/aa));
165 }
166 else {
167 return make_pair(alpha,make_pair(0.,0.));
168 }
169 }
170
171 /// Normalise histograms etc., after the run
172 void finalize() {
173 // spectra
174 normalize(_h_Jpsi,1.,false);
175 normalize(_h_feed,1.,false);
176 normalize(_h_Psi2,1.,false);
177 // cross sections
178 double fact = 1./sumOfWeights()*crossSection()/picobarn;
179 scale(_h_sig_JPsi_all ,fact);
180 scale(_h_sig_Jpsi_high,fact);
181 scale(_h_sig_Jpsi_feed,fact);
182 scale(_h_sig_Psi2_high,fact);
183 // angular dists and parameters from them
184 vector<double> bins = {2.,2.6,3.4,4.9};
185 Scatter2DPtr _h_A;
186 book(_h_A ,2,1,1);
187 Scatter2DPtr _h_alpha;
188 book(_h_alpha,2,1,2);
189 for(unsigned int ix=0;ix<3;++ix) {
190 normalize(_h_cThetaStar.histos()[ix]);
191 pair<double,pair<double,double> > alpha = calcAlpha(_h_cThetaStar.histos()[ix]);
192 double centre = 0.5*(bins[ix+1]+bins[ix]);
193 double width = 0.5*(bins[ix+1]-bins[ix]);
194 _h_A->addPoint(centre, alpha.first, make_pair(width,width),
195 make_pair(alpha.second.first,alpha.second.second) );
196 normalize(_h_cThetaH.histos()[ix]);
197 alpha = calcAlpha(_h_cThetaH.histos()[ix]);
198 _h_alpha->addPoint(centre, alpha.first, make_pair(width,width),
199 make_pair(alpha.second.first,alpha.second.second) );
200 }
201 double centre1 = 0.5*(bins[2]+bins[0]);
202 double width1 = 0.5*(bins[2]-bins[0]);
203 double centre2 = 0.5*(bins[3]+bins[0]);
204 double width2 = 0.5*(bins[3]-bins[0]);
205 normalize(_h_cS_low);
206 pair<double,pair<double,double> > alpha = calcAlpha(_h_cS_low);
207 book(_h_A ,2,2,1);
208 _h_A->addPoint(centre1, alpha.first, make_pair(width1,width1),
209 make_pair(alpha.second.first,alpha.second.second) );
210 normalize(_h_cS_high);
211 alpha = calcAlpha(_h_cS_high);
212 book(_h_A ,2,3,1);
213 _h_A->addPoint(centre2, alpha.first, make_pair(width2,width2),
214 make_pair(alpha.second.first,alpha.second.second) );
215 normalize(_h_cH_low);
216 alpha = calcAlpha(_h_cH_low);
217 book(_h_alpha,2,2,2);
218 _h_alpha->addPoint(centre1, alpha.first, make_pair(width1,width1),
219 make_pair(alpha.second.first,alpha.second.second) );
220 normalize(_h_cH_high);
221 alpha = calcAlpha(_h_cH_high);
222 book(_h_alpha,2,3,2);
223 _h_alpha->addPoint(centre2, alpha.first, make_pair(width2,width2),
224 make_pair(alpha.second.first,alpha.second.second) );
225 }
226
227 ///@}
228
229
230 /// @name Histograms
231 ///@{
232 Histo1DPtr _h_Jpsi,_h_Psi2,_h_feed;
233 Histo1DPtr _h_sig_JPsi_all,_h_sig_Jpsi_high,_h_sig_Jpsi_feed,_h_sig_Psi2_high;
234 Histo1DPtr _h_cS_low,_h_cS_high,_h_cH_low,_h_cH_high;
235 BinnedHistogram _h_cThetaStar,_h_cThetaH;
236
237
238 ///@}
239
240
241 };
242
243
244 RIVET_DECLARE_PLUGIN(BELLE_2002_I563840);
245
246}
|