rivet is hosted by Hepforge, IPPP Durham

Rivet analyses reference

CLEOII_2002_I606309

Momentum Spectra for $J/\psi$ and $\psi(2S)$ production in B decays
Experiment: CLEOII (CESR)
Inspire ID: 606309
Status: VALIDATED
Authors:
  • Peter Richardson
References:
  • Phys.Rev.Lett. 89 (2002) 282001
Beams: * *
Beam energies: ANY
Run details:
  • any process with Upsilon(4S) decay, original e+e-

Momentum Spectra for $J/\psi$ and $\psi(2S)$ production in B decays at the $\Upsilon(4S)$

Source code: CLEOII_2002_I606309.cc
  1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// @brief J/Psi and Psi(2S) spectra at the Upsilon(4S)
  9  class CLEOII_2002_I606309 : public Analysis {
 10  public:
 11
 12    /// Constructor
 13    RIVET_DEFAULT_ANALYSIS_CTOR(CLEOII_2002_I606309);
 14
 15
 16    /// @name Analysis methods
 17    ///@{
 18
 19    /// Book histograms and initialise projections before the run
 20    void init() {
 21      // projections
 22      declare(UnstableParticles(), "UFS");
 23      //histograms
 24      book(_weightSum,"/TMP/weightSum");
 25      book(_h_Jpsi            ,3,1,1);
 26      book(_h_Psi_prime       ,3,1,2);
 27      book(_h_cTheta_Jpsi     ,4,1,1);
 28      book(_h_cTheta_Psi_prime,4,1,2);
 29      double bins[4]={0.,.8,1.4,2.};
 30      for(unsigned int ix=0;ix<3;++ix) {
 31	Histo1DPtr temp;
 32	std::ostringstream title;
 33	title << "/TMP/ctheta_" << ix;
 34	book(temp,title.str(),20,-1.,1.);
 35	_h2_cTheta_Jpsi.add(bins[ix],bins[ix+1],temp);
 36      }
 37    }
 38
 39    void findDecayProducts(Particle parent, Particles & charm) {
 40      for (const Particle & p :parent.children()) {
 41        if (p.pid()==443) {
 42	  charm.push_back(p);
 43	  continue;
 44	}
 45	else if (p.pid()==100443) {
 46	  charm.push_back(p);
 47	}
 48	if(!p.children().empty())
 49	  findDecayProducts(p,charm);
 50      }
 51    }
 52
 53    void findLeptons(const Particle & mother,
 54		     unsigned int & nstable,
 55		     Particles& lp, Particles& lm) {
 56      for(const Particle & p : mother.children()) {
 57        int id = p.pid();
 58      	if ( id == 11 || id == 13 ) {
 59	  lm.push_back(p);
 60	  ++nstable;
 61	}
 62       	else if (id == -11 || id==-13) {
 63       	  lp.push_back(p);
 64       	  ++nstable;
 65       	}
 66	else if ( !p.children().empty() ) {
 67	  findLeptons(p,nstable,lp,lm);
 68	}
 69	else
 70	  ++nstable;
 71      }
 72    }
 73
 74    /// Perform the per-event analysis
 75    void analyze(const Event& event) {    
 76      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
 77      for (const Particle& p : ufs.particles(Cuts::pid==300553)) {
 78        _weightSum->fill();
 79        const LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(p.momentum().betaVec());
 80	Particles charm;
 81	findDecayProducts(p,charm);
 82	for(const Particle & pp : charm) {
 83	  FourMomentum pcm = boost.transform(pp.momentum());
 84	  unsigned int nstable = 0;
 85	  Particles lp, lm;
 86	  findLeptons(pp,nstable,lp,lm);
 87	  double cTheta(0.);
 88	  bool foundLeptons(false);
 89	  if(nstable==2&&lp.size()==1&&lm.size()==1) {
 90	    foundLeptons=true;
 91	    FourMomentum pl = boost.transform(lm[0].momentum());
 92	    const LorentzTransform boost2 = LorentzTransform::mkFrameTransformFromBeta(pcm.betaVec());
 93	    pl = boost2.transform(pl);
 94	    cTheta = pl.p3().unit().dot(pcm.p3().unit());
 95	  }
 96	  if(pp.pid()==443) {
 97	    _h_Jpsi->fill(pcm.p3().mod());
 98	    if(foundLeptons) {
 99	      _h_cTheta_Jpsi->fill(cTheta);
100	      _h2_cTheta_Jpsi.fill(pcm.p3().mod(),cTheta);
101	    }
102	  }
103	  else {
104	    _h_Psi_prime->fill(pcm.p3().mod());
105	    if(foundLeptons) {
106	      _h_cTheta_Psi_prime->fill(cTheta);
107	    }
108	  }
109	}
110      }
111    }
112    
113    pair<double,pair<double,double> > calcAlpha(Histo1DPtr hist) {
114      if(hist->numEntries()==0.) return make_pair(0.,make_pair(0.,0.));
115      double d = 3./(pow(hist->xMax(),3)-pow(hist->xMin(),3));
116      double c = 3.*(hist->xMax()-hist->xMin())/(pow(hist->xMax(),3)-pow(hist->xMin(),3));
117      double sum1(0.),sum2(0.),sum3(0.),sum4(0.),sum5(0.);
118      for (auto bin : hist->bins() ) {
119       	double Oi = bin.area();
120	if(Oi==0.) continue;
121	double a =  d*(bin.xMax() - bin.xMin());
122	double b = d/3.*(pow(bin.xMax(),3) - pow(bin.xMin(),3));
123       	double Ei = bin.areaErr();
124	sum1 +=   a*Oi/sqr(Ei);
125	sum2 +=   b*Oi/sqr(Ei);
126	sum3 += sqr(a)/sqr(Ei);
127	sum4 += sqr(b)/sqr(Ei);
128	sum5 +=    a*b/sqr(Ei);
129      }
130      // calculate alpha
131      double alpha = (-c*sum1 + sqr(c)*sum2 + sum3 - c*sum5)/(sum1 - c*sum2 + c*sum4 - sum5);
132      // and error
133      double cc = -pow((sum3 + sqr(c)*sum4 - 2*c*sum5),3);
134      double bb = -2*sqr(sum3 + sqr(c)*sum4 - 2*c*sum5)*(sum1 - c*sum2 + c*sum4 - sum5);
135      double aa =  sqr(sum1 - c*sum2 + c*sum4 - sum5)*(-sum3 - sqr(c)*sum4 + sqr(sum1 - c*sum2 + c*sum4 - sum5) + 2*c*sum5);      
136      double dis = sqr(bb)-4.*aa*cc;
137      if(dis>0.) {
138	dis = sqrt(dis);
139	return make_pair(alpha,make_pair(0.5*(-bb+dis)/aa,-0.5*(-bb-dis)/aa));
140      }
141      else {
142	return make_pair(alpha,make_pair(0.,0.));
143      }
144    }
145
146    /// Normalise histograms etc., after the run
147    void finalize() {
148      if(_weightSum->val()==0.) return;
149      scale(_h_Jpsi     , 50./ *_weightSum);
150      scale(_h_Psi_prime, 50./ *_weightSum);
151      // polarization J/psi
152      normalize(_h_cTheta_Jpsi);
153      pair<double,pair<double,double> > alpha = calcAlpha(_h_cTheta_Jpsi);
154      Scatter2DPtr _h_alpha;
155      book(_h_alpha,1,1,1);
156      _h_alpha->addPoint(1., alpha.first, make_pair(1.,1.),
157			 make_pair(alpha.second.first,alpha.second.second) );
158      // polarization psi(2S)
159      normalize(_h_cTheta_Psi_prime);
160      alpha = calcAlpha(_h_cTheta_Psi_prime);
161      book(_h_alpha,1,1,2);
162      _h_alpha->addPoint(0.8, alpha.first, make_pair(0.8,0.8),
163			 make_pair(alpha.second.first,alpha.second.second) );
164
165      // J/psi as function of momentum
166      double bins[4]={0.,.8,1.4,2.};
167      book(_h_alpha,2,1,1);
168      for(unsigned int ix=0;ix<3;++ix) {
169	normalize(_h2_cTheta_Jpsi.histos()[ix]);
170	double cen = 0.5*(bins[ix+1]+bins[ix]);
171	double wid = 0.5*(bins[ix+1]-bins[ix]);
172	alpha = calcAlpha(_h2_cTheta_Jpsi.histos()[ix]);
173	_h_alpha->addPoint(cen, alpha.first, make_pair(wid,wid),
174			 make_pair(alpha.second.first,alpha.second.second) );
175      }
176    }
177
178    ///@}
179
180
181    /// @name Histograms
182    ///@{
183    // count of weights
184    CounterPtr _weightSum;
185    // histograms
186    Histo1DPtr _h_Jpsi,_h_Psi_prime;
187    Histo1DPtr _h_cTheta_Jpsi,_h_cTheta_Psi_prime;
188    BinnedHistogram _h2_cTheta_Jpsi;
189    ///@}
190
191
192  };
193
194
195  RIVET_DECLARE_PLUGIN(CLEOII_2002_I606309);
196
197}