rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

BELLE_2002_I563840

Prompt Charmonium production at 10.6 GeV
Experiment: BELLE (KEKB)
Inspire ID: 563840
Status: VALIDATED
Authors:
  • Peter Richardson
References:
  • Phys.Rev.Lett. 88 (2002) 052001
Beams: e+ e-
Beam energies: (5.3, 5.3) GeV
Run details:
  • e+e- to hadrons, only the continuuum contributes not Upslion(4S) decays

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}