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