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      book(_h2_cTheta_Jpsi, {0.,0.8,1.4,2.});
 30      for (auto& b : _h2_cTheta_Jpsi->bins()) {
 31        book(b, "/TMP/ctheta_"+to_string(b.index()-1), 20, -1.0, 1.0);
 32      }
 33    }
 34
 35    void findDecayProducts(Particle parent, Particles& charm) {
 36      for (const Particle & p :parent.children()) {
 37        if (p.pid()==443) {
 38          charm.push_back(p);
 39          continue;
 40        }
 41        else if (p.pid()==100443) {
 42          charm.push_back(p);
 43        }
 44        if(!p.children().empty()) {
 45          findDecayProducts(p,charm);
 46        }
 47      }
 48    }
 49
 50    void findLeptons(const Particle & mother, unsigned int& nstable, Particles& lp, Particles& lm) {
 51      for (const Particle& p : mother.children()) {
 52        int id = p.pid();
 53      	if ( id == 11 || id == 13 ) {
 54          lm.push_back(p);
 55          ++nstable;
 56        }
 57       	else if (id == -11 || id==-13) {
 58       	  lp.push_back(p);
 59       	  ++nstable;
 60       	}
 61        else if ( !p.children().empty() ) {
 62          findLeptons(p,nstable,lp,lm);
 63        }
 64        else {
 65          ++nstable;
 66        }
 67      }
 68    }
 69
 70    /// Perform the per-event analysis
 71    void analyze(const Event& event) {
 72      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
 73      for (const Particle& p : ufs.particles(Cuts::pid==300553)) {
 74        _weightSum->fill();
 75        const LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(p.momentum().betaVec());
 76        Particles charm;
 77        findDecayProducts(p,charm);
 78        for (const Particle& pp : charm) {
 79          FourMomentum pcm = boost.transform(pp.momentum());
 80          unsigned int nstable = 0;
 81          Particles lp, lm;
 82          findLeptons(pp,nstable,lp,lm);
 83          double cTheta(0.);
 84          bool foundLeptons(false);
 85          if (nstable==2&&lp.size()==1&&lm.size()==1) {
 86            foundLeptons=true;
 87            FourMomentum pl = boost.transform(lm[0].momentum());
 88            const LorentzTransform boost2 = LorentzTransform::mkFrameTransformFromBeta(pcm.betaVec());
 89            pl = boost2.transform(pl);
 90            cTheta = pl.p3().unit().dot(pcm.p3().unit());
 91          }
 92          if(pp.pid()==443) {
 93            _h_Jpsi->fill(pcm.p3().mod());
 94            if(foundLeptons) {
 95              _h_cTheta_Jpsi->fill(cTheta);
 96              _h2_cTheta_Jpsi->fill(pcm.p3().mod(), cTheta);
 97            }
 98          }
 99          else {
100            _h_Psi_prime->fill(pcm.p3().mod());
101            if(foundLeptons) {
102              _h_cTheta_Psi_prime->fill(cTheta);
103            }
104          }
105        }
106      }
107    }
108
109    pair<double,pair<double,double> > calcAlpha(const Histo1DPtr& hist) {
110      if(hist->numEntries()==0.) return make_pair(0.,make_pair(0.,0.));
111      double d = 3./(pow(hist->xMax(),3)-pow(hist->xMin(),3));
112      double c = 3.*(hist->xMax()-hist->xMin())/(pow(hist->xMax(),3)-pow(hist->xMin(),3));
113      double sum1(0.),sum2(0.),sum3(0.),sum4(0.),sum5(0.);
114      for (const auto& bin : hist->bins() ) {
115       	double Oi = bin.sumW();
116        if(Oi==0.) continue;
117        double a =  d*(bin.xMax() - bin.xMin());
118        double b = d/3.*(pow(bin.xMax(),3) - pow(bin.xMin(),3));
119       	double Ei = bin.errW();
120        sum1 +=   a*Oi/sqr(Ei);
121        sum2 +=   b*Oi/sqr(Ei);
122        sum3 += sqr(a)/sqr(Ei);
123        sum4 += sqr(b)/sqr(Ei);
124        sum5 +=    a*b/sqr(Ei);
125      }
126      // calculate alpha
127      double alpha = (-c*sum1 + sqr(c)*sum2 + sum3 - c*sum5)/(sum1 - c*sum2 + c*sum4 - sum5);
128      // and error
129      double cc = -pow((sum3 + sqr(c)*sum4 - 2*c*sum5),3);
130      double bb = -2*sqr(sum3 + sqr(c)*sum4 - 2*c*sum5)*(sum1 - c*sum2 + c*sum4 - sum5);
131      double aa =  sqr(sum1 - c*sum2 + c*sum4 - sum5)*(-sum3 - sqr(c)*sum4 + sqr(sum1 - c*sum2 + c*sum4 - sum5) + 2*c*sum5);
132      double dis = sqr(bb)-4.*aa*cc;
133      if (dis>0.) {
134        dis = sqrt(dis);
135        return make_pair(alpha,make_pair(0.5*(-bb+dis)/aa,-0.5*(-bb-dis)/aa));
136      }
137      else {
138        return make_pair(alpha,make_pair(0.,0.));
139      }
140    }
141
142    /// Normalise histograms etc., after the run
143    void finalize() {
144      if(_weightSum->val()==0.) return;
145      scale(_h_Jpsi     , 50./ *_weightSum);
146      scale(_h_Psi_prime, 50./ *_weightSum);
147      // polarization J/psi
148      normalize(_h_cTheta_Jpsi);
149      pair<double,pair<double,double> > alpha = calcAlpha(_h_cTheta_Jpsi);
150      Estimate1DPtr _h_alpha;
151      book(_h_alpha,1,1,1);
152      _h_alpha->bin(1).set(alpha.first, alpha.second);
153      // polarization psi(2S)
154      normalize(_h_cTheta_Psi_prime);
155      alpha = calcAlpha(_h_cTheta_Psi_prime);
156      book(_h_alpha,1,1,2);
157      _h_alpha->bin(1).set(alpha.first, alpha.second);
158
159      // J/psi as function of momentum
160      book(_h_alpha,2,1,1);
161      for (auto& b : _h2_cTheta_Jpsi->bins()) {
162        normalize(b);
163        alpha = calcAlpha(b);
164        _h_alpha->bin(b.index()).set(alpha.first, alpha.second);
165      }
166    }
167
168    /// @}
169
170
171    /// @name Histograms
172    /// @{
173    // count of weights
174    CounterPtr _weightSum;
175    // histograms
176    Histo1DPtr _h_Jpsi,_h_Psi_prime;
177    Histo1DPtr _h_cTheta_Jpsi,_h_cTheta_Psi_prime;
178    Histo1DGroupPtr _h2_cTheta_Jpsi;
179    /// @}
180
181
182  };
183
184
185  RIVET_DECLARE_PLUGIN(CLEOII_2002_I606309);
186
187}