1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/Beam.hh"
  4#include "Rivet/Projections/UnstableParticles.hh"
  5
  6namespace Rivet {
  7
  8
  9  /// @brief D* helicity in e+e- at 10.5 GeV
 10  class CLEO_1998_I467595 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(CLEO_1998_I467595);
 15
 16
 17    /// @name Analysis methods
 18    /// @{
 19
 20    /// Book histograms and initialise projections before the run
 21    void init() {
 22
 23      // Initialise and register projections
 24      declare(Beam(), "Beams");
 25      declare(UnstableParticles(), "UFS");
 26
 27      // Book histograms
 28      book(_h_ctheta, {0.25, 0.45, 0.55, 0.65, 0.75, 0.85, 1.});
 29      for (auto& b : _h_ctheta->bins()) {
 30        book(b, 5, 1, b.index());
 31      }
 32
 33    }
 34
 35
 36    /// Perform the per-event analysis
 37    void analyze(const Event& event) {
 38      // Get beams and average beam momentum
 39      const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
 40      const double Emax = ( beams.first.p3().mod() + beams.second.p3().mod() ) / 2.0;
 41      const double Pmax = sqrt(sqr(Emax)-sqr(2.01026));
 42
 43      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
 44      for (const Particle& p : ufs.particles(Cuts::abspid==413)) {
 45        if (p.children().size()!=2) continue;
 46        int sign = p.pid()/413;
 47        Particle D0;
 48        if(p.children()[0].pid()==sign*421 && p.children()[1].pid()==sign*211) {
 49          D0 = p.children()[0];
 50        }
 51        else if(p.children()[1].pid()==sign*421 && p.children()[0].pid()==sign*211) {
 52          D0 = p.children()[1];
 53        }
 54        else {
 55          continue;
 56        }
 57        LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(p.momentum().betaVec());
 58        double xp = (p.momentum().p3().mod()+p.momentum().t())/(Pmax+Emax);
 59        Vector3 e1z = p.momentum().p3().unit();
 60        Vector3 axis1 = boost.transform(D0.momentum()).p3().unit();
 61        double ctheta = e1z.dot(axis1);
 62        _h_ctheta->fill(xp,ctheta);
 63      }
 64    }
 65
 66    pair<double,double> calcRho(const Histo1DPtr& hist) {
 67      if(hist->numEntries()==0.) return make_pair(0.,0.);
 68      double sum1(0.),sum2(0.);
 69      for (const auto& bin : hist->bins() ) {
 70        double Oi = bin.sumW();
 71        if(Oi==0.) continue;
 72        double ai = 0.25*(bin.xMax()*(3.-sqr(bin.xMax())) - bin.xMin()*(3.-sqr(bin.xMin())));
 73        double bi = 0.75*(bin.xMin()*(1.-sqr(bin.xMin())) - bin.xMax()*(1.-sqr(bin.xMax())));
 74        double Ei = bin.errW();
 75        sum1 += sqr(bi/Ei);
 76        sum2 += bi/sqr(Ei)*(Oi-ai);
 77      }
 78      return make_pair(sum2/sum1,sqrt(1./sum1));
 79    }
 80
 81    pair<double,pair<double,double> > calcAlpha(const Histo1DPtr& hist) {
 82      if(hist->numEntries()==0.) return make_pair(0.,make_pair(0.,0.));
 83      double d = 3./(pow(hist->xMax(),3)-pow(hist->xMin(),3));
 84      double c = 3.*(hist->xMax()-hist->xMin())/(pow(hist->xMax(),3)-pow(hist->xMin(),3));
 85      double sum1(0.),sum2(0.),sum3(0.),sum4(0.),sum5(0.);
 86      for (const auto& bin : hist->bins() ) {
 87       	double Oi = bin.sumW();
 88        if(Oi==0.) continue;
 89        double a =  d*(bin.xMax() - bin.xMin());
 90        double b = d/3.*(pow(bin.xMax(),3) - pow(bin.xMin(),3));
 91       	double Ei = bin.errW();
 92        sum1 +=   a*Oi/sqr(Ei);
 93        sum2 +=   b*Oi/sqr(Ei);
 94        sum3 += sqr(a)/sqr(Ei);
 95        sum4 += sqr(b)/sqr(Ei);
 96        sum5 +=    a*b/sqr(Ei);
 97      }
 98      // calculate alpha
 99      double alpha = (-c*sum1 + sqr(c)*sum2 + sum3 - c*sum5)/(sum1 - c*sum2 + c*sum4 - sum5);
100      // and error
101      double cc = -pow((sum3 + sqr(c)*sum4 - 2*c*sum5),3);
102      double bb = -2*sqr(sum3 + sqr(c)*sum4 - 2*c*sum5)*(sum1 - c*sum2 + c*sum4 - sum5);
103      double aa =  sqr(sum1 - c*sum2 + c*sum4 - sum5)*(-sum3 - sqr(c)*sum4 + sqr(sum1 - c*sum2 + c*sum4 - sum5) + 2*c*sum5);
104      double dis = sqr(bb)-4.*aa*cc;
105      if(dis>0.) {
106        dis = sqrt(dis);
107        return make_pair(alpha,make_pair(0.5*(-bb+dis)/aa,-0.5*(-bb-dis)/aa));
108      }
109      else {
110        return make_pair(alpha,make_pair(0.,0.));
111      }
112    }
113
114    /// Normalise histograms etc., after the run
115    void finalize() {
116      vector<double> x = {0.25,0.45,0.55,0.65,0.75,0.85,1.};
117      Estimate1DPtr h_alpha;
118      book(h_alpha, 3,1,1);
119      Estimate1DPtr h_rho;
120      book(h_rho, 4,1,1);
121      for (auto& b : _h_ctheta->bins()) {
122        normalize(b);
123        pair<double,double> rho00 = calcRho(b);
124        h_rho->bin(b.index()).set(rho00.first, rho00.second);
125        pair<double,pair<double,double> > alpha = calcAlpha(b);
126        h_alpha->bin(b.index()).set(alpha.first, alpha.second);
127      }
128    }
129
130    /// @}
131
132
133    /// @name Histograms
134    /// @{
135    Histo1DGroupPtr _h_ctheta;
136    /// @}
137
138
139  };
140
141
142  RIVET_DECLARE_PLUGIN(CLEO_1998_I467595);
143
144
145}