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*+ polarization
 10  class CLEO_1991_I314060 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(CLEO_1991_I314060);
 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(_h_ctheta, {0.25, 0.45, 0.55, 0.65, 0.75, 0.85, 1.00});
 28      for (auto& b : _h_ctheta->bins()) {
 29        book(b, 2, 1, b.index());
 30      }
 31
 32    }
 33
 34    /// Recursively walk the decay tree to find decay products of @a p
 35    void findDecayProducts(Particle mother, Particles & d0, Particles & pi,unsigned int & ncount) {
 36      for (const Particle& p: mother.children()) {
 37        if(p.abspid()==421) {
 38          d0.push_back(p);
 39        }
 40        else if(p.abspid()==211) {
 41          pi.push_back(p);
 42        }
 43        ncount +=1;
 44      }
 45    }
 46
 47    /// Perform the per-event analysis
 48    void analyze(const Event& event) {
 49      // Get beams and average beam momentum
 50      const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
 51      const double meanBeamMom = ( beams.first.p3().mod() +
 52                                   beams.second.p3().mod() ) / 2.0;
 53      MSG_DEBUG("Avg beam momentum = " << meanBeamMom);
 54      // loop over D*+ mesons
 55      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
 56      for (const Particle& p: ufs.particles(Cuts::abspid==413)) {
 57        // calc x+
 58        double x = (p.momentum().E()+p.momentum().z())/(meanBeamMom + sqrt(sqr(meanBeamMom)+p.mass2()));
 59        // checck decay products
 60        Particles d0,pi;
 61        unsigned int ncount=0;
 62        findDecayProducts(p,d0, pi,ncount);
 63        if (ncount!=2 || pi.size()!=1 || d0.size()!=1 ) continue;
 64        if (pi[0].pid()/p.pid()<0) continue;
 65        LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(p.momentum().betaVec());
 66        Vector3 d1 = boost.transform(pi[0].momentum()).p3().unit();
 67        double ctheta  = d1.dot(p.momentum().p3().unit());
 68        _h_ctheta->fill(x, ctheta);
 69      }
 70    }
 71
 72    pair<double,pair<double,double> > calcAlpha(const Histo1DPtr& hist) {
 73      if (hist->numEntries()==0.) return make_pair(0.,make_pair(0.,0.));
 74      double sum1(0.),sum2(0.),sum3(0.),sum4(0.),sum5(0.);
 75      for (const auto& bin : hist->bins()) {
 76       	double Oi = bin.sumW();
 77        if(Oi==0.) continue;
 78        double a =  1.5*(bin.xMax() - bin.xMin());
 79        double b = 0.5*(pow(bin.xMax(),3) - pow(bin.xMin(),3));
 80        double Ei = bin.errW();
 81        sum1 +=   a*Oi/sqr(Ei);
 82        sum2 +=   b*Oi/sqr(Ei);
 83        sum3 += sqr(a)/sqr(Ei);
 84        sum4 += sqr(b)/sqr(Ei);
 85        sum5 +=    a*b/sqr(Ei);
 86      }
 87      // calculate alpha
 88      double alpha = (-3*sum1 + 9*sum2 + sum3 - 3*sum5)/(sum1 - 3*sum2 + 3*sum4 - sum5);
 89      // and error
 90      double cc = -pow((sum3 + 9*sum4 - 6*sum5),3);
 91      double bb = -2*sqr(sum3 + 9*sum4 - 6*sum5)*(sum1 - 3*sum2 + 3*sum4 - sum5);
 92      double aa =  sqr(sum1 - 3*sum2 + 3*sum4 - sum5)*(-sum3 - 9*sum4 + sqr(sum1 - 3*sum2 + 3*sum4 - sum5) + 6*sum5);
 93      double dis = sqr(bb)-4.*aa*cc;
 94      if (dis>0.) {
 95        dis = sqrt(dis);
 96        return make_pair(alpha,make_pair(0.5*(-bb+dis)/aa,-0.5*(-bb-dis)/aa));
 97      }
 98      else {
 99        return make_pair(alpha,make_pair(0.,0.));
100      }
101    }
102
103    /// Normalise histograms etc., after the run
104    void finalize() {
105      //vector<double> x   = {0.35,0.5 ,0.6 ,0.7 ,0.8 ,0.925};
106      //vector<double> wid = {0.10,0.05,0.05,0.05,0.05,0.075};
107      Estimate1DPtr h_alpha;
108      book(h_alpha, 1,1,1);
109      Estimate1DPtr h_rho;
110      book(h_rho,   1,1,2);
111      Estimate1DPtr h_eta;
112      book(h_eta,   1,1,3);
113      unsigned int ix=0;
114      for (auto& b : _h_ctheta->bins()) {
115        // normalize
116        normalize(b);
117        // alpha
118        pair<double,pair<double,double> > alpha = calcAlpha(b);
119        h_alpha->bin(1+ix).set(alpha.first, alpha.second);
120        // rho
121        const double rho = (1.+alpha.first)/(3.+alpha.first);
122        pair<double,double> rho_error;
123        rho_error.first  = 2.*alpha.second.first /sqr(3.+alpha.first);
124        rho_error.second = 2.*alpha.second.second/sqr(3.+alpha.first);
125        h_rho->bin(1+ix).set(rho, rho_error);
126        // eta
127        const double eta = alpha.first/(3.+alpha.first);
128        pair<double,double> eta_error;
129        eta_error.first  = 3.*alpha.second.first /sqr(3.+alpha.first);
130        eta_error.second = 3.*alpha.second.second/sqr(3.+alpha.first);
131        h_eta->bin(1+ix).set(eta, eta_error);
132        ++ix;
133      }
134    }
135
136    /// @}
137
138
139    /// @name Histograms
140    /// @{
141    Histo1DGroupPtr _h_ctheta;
142    /// @}
143
144
145  };
146
147
148  RIVET_DECLARE_PLUGIN(CLEO_1991_I314060);
149
150
151}