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