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* poliarzation at 29 GeV
 10  class TPC_1991_I316132 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(TPC_1991_I316132);
 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.3, 0.4, 0.5, 0.6, 0.7, 0.8, 1.0});
 29      for (auto& b : _h_ctheta->bins()) {
 30        const string name = "ctheta_" + std::to_string(b.index()-1);
 31        book(b, name, 20, -1.0, 1.0);
 32      }
 33      book(_h_ctheta_all, "ctheta_all",20,-1,1);
 34
 35      book(_h_phi, {0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 1.0});
 36      for (auto& b : _h_phi->bins()) {
 37        const string name = "phi_" + std::to_string(b.index()-1);
 38        book(b, name, 20, -M_PI, M_PI);
 39      }
 40      book(_h_phi_all, "phi_all",20,-M_PI,M_PI);
 41
 42      book(_h_01, {0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 1.0});
 43      for (auto& b : _h_01->bins()) {
 44        const string name = "h_01_" + std::to_string(b.index()-1);
 45        book(b, name, 20, -1.0, 1.0);
 46      }
 47      book(_h_01_all, "h_01_all",20,-1,1);
 48    }
 49
 50
 51    /// Perform the per-event analysis
 52    void analyze(const Event& event) {
 53      // Get beams and average beam momentum
 54      const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
 55      const double Emax = ( beams.first.p3().mod() + beams.second.p3().mod() ) / 2.0;
 56      Vector3 axis;
 57      if (beams.first.pid()>0) {
 58        axis = beams.first .momentum().p3().unit();
 59      }
 60      else {
 61        axis = beams.second.momentum().p3().unit();
 62      }
 63
 64      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
 65      for  (const Particle& p : ufs.particles(Cuts::abspid==413)) {
 66        if (p.children().size()!=2) continue;
 67        int sign = p.pid()/413;
 68        Particle D0;
 69        if (p.children()[0].pid()==sign*421 && p.children()[1].pid()==sign*211) {
 70          D0 = p.children()[0];
 71        }
 72        else if (p.children()[1].pid()==sign*421 && p.children()[0].pid()==sign*211) {
 73          D0 = p.children()[1];
 74        }
 75        else {
 76          continue;
 77        }
 78        LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(p.momentum().betaVec());
 79        const double xE = p.momentum().t()/Emax;
 80        Vector3 e1z = p.momentum().p3().unit();
 81        Vector3 e1y = e1z.cross(axis).unit();
 82        Vector3 e1x = e1y.cross(e1z).unit();
 83        Vector3 axis1 = boost.transform(D0.momentum()).p3().unit();
 84        const double ctheta = e1z.dot(axis1);
 85        _h_ctheta->fill(xE,ctheta);
 86        const double phi = atan2(e1y.dot(axis1),e1x.dot(axis1));
 87        _h_phi->fill(xE, phi);
 88        _h_01->fill(xE, ctheta, cos(phi));
 89
 90        if(xE>0.3) {
 91          _h_ctheta_all->fill(ctheta);
 92          _h_phi_all->fill(phi);
 93          _h_01_all->fill(ctheta, cos(phi));
 94        }
 95      }
 96    }
 97
 98    pair<double,double> calcRho(Histo1DPtr hist,unsigned int mode) {
 99      if(hist->numEntries()==0.) return make_pair(0.,0.);
100      double sum1(0.),sum2(0.);
101      for (const auto& bin : hist->bins() ) {
102        double Oi = bin.sumW();
103        if (Oi==0.) continue;
104        double ai,bi;
105        if (mode==0) {
106          ai = 0.25*(bin.xMax()*(3.-sqr(bin.xMax())) - bin.xMin()*(3.-sqr(bin.xMin())));
107          bi = 0.75*(bin.xMin()*(1.-sqr(bin.xMin())) - bin.xMax()*(1.-sqr(bin.xMax())));
108        }
109        else if (mode==1) {
110          ai = 0.5/M_PI*(bin.xMax()-bin.xMin());
111          bi = 0.5/M_PI*(sin(2.*bin.xMin())-sin(2.*bin.xMax()));
112        }
113        else  {
114          ai=0.;
115          bi= sqrt(0.5)*((1.-sqr(bin.xMax()))*sqrt(1.-sqr(bin.xMax()))-
116             (1.-sqr(bin.xMin()))*sqrt(1.-sqr(bin.xMin())));
117        }
118        double Ei = bin.errW();
119        sum1 += sqr(bi/Ei);
120        sum2 += bi/sqr(Ei)*(Oi-ai);
121      }
122      return make_pair(sum2/sum1,sqrt(1./sum1));
123    }
124
125    pair<double,pair<double,double> > calcAlpha(Histo1DPtr hist) {
126      if(hist->numEntries()==0.) return make_pair(0.,make_pair(0.,0.));
127      double d = 3./(pow(hist->xMax(),3)-pow(hist->xMin(),3));
128      double c = 3.*(hist->xMax()-hist->xMin())/(pow(hist->xMax(),3)-pow(hist->xMin(),3));
129      double sum1(0.),sum2(0.),sum3(0.),sum4(0.),sum5(0.);
130      for (const auto& bin : hist->bins() ) {
131       	double Oi = bin.sumW();
132        if (Oi==0.) continue;
133        double a =  d*(bin.xMax() - bin.xMin());
134        double b = d/3.*(pow(bin.xMax(),3) - pow(bin.xMin(),3));
135       	double Ei = bin.errW();
136        sum1 +=   a*Oi/sqr(Ei);
137        sum2 +=   b*Oi/sqr(Ei);
138        sum3 += sqr(a)/sqr(Ei);
139        sum4 += sqr(b)/sqr(Ei);
140        sum5 +=    a*b/sqr(Ei);
141      }
142      // calculate alpha
143      double alpha = (-c*sum1 + sqr(c)*sum2 + sum3 - c*sum5)/(sum1 - c*sum2 + c*sum4 - sum5);
144      // and error
145      double cc = -pow((sum3 + sqr(c)*sum4 - 2*c*sum5),3);
146      double bb = -2*sqr(sum3 + sqr(c)*sum4 - 2*c*sum5)*(sum1 - c*sum2 + c*sum4 - sum5);
147      double aa =  sqr(sum1 - c*sum2 + c*sum4 - sum5)*(-sum3 - sqr(c)*sum4 + sqr(sum1 - c*sum2 + c*sum4 - sum5) + 2*c*sum5);
148      double dis = sqr(bb)-4.*aa*cc;
149      if (dis>0.) {
150        dis = sqrt(dis);
151        return make_pair(alpha,make_pair(0.5*(-bb+dis)/aa,-0.5*(-bb-dis)/aa));
152      }
153      else {
154        return make_pair(alpha,make_pair(0.,0.));
155      }
156    }
157
158    /// Normalise histograms etc., after the run
159    void finalize() {
160      Estimate1DPtr h_alpha, h_rho00, h_rhooff, h_01;
161      book(h_alpha , 1,1,1);
162      book(h_rho00 , 2,1,1);
163      book(h_rhooff, 2,1,2);
164      book(h_01    , 2,1,3);
165      for (size_t ix = 1; ix < _h_ctheta->numBins()+1; ++ix) {
166        const double integral = _h_ctheta->bin(ix)->integral();
167        scale(_h_ctheta->bin(ix), 1./integral);
168        scale(_h_01->bin(ix), 1./integral);
169        normalize(_h_phi->bin(ix));
170
171        pair<double,double> rho00 = calcRho(_h_ctheta->bin(ix), 0);
172        h_rho00->bin(ix).set(rho00.first, rho00.second);
173
174        pair<double,pair<double,double> > alpha = calcAlpha(_h_ctheta->bin(ix));
175        h_alpha->bin(ix).set(alpha.first, alpha.second);
176
177        pair<double,double> rhooff = calcRho(_h_phi->bin(ix), 1);
178        h_rhooff->bin(ix).set(rhooff.first, rhooff.second);
179
180        pair<double,double> rho01 = calcRho(_h_01->bin(ix), 2);
181        h_01->bin(ix).set(rho01.first, rho01.second);
182      }
183      // integral over z
184      double integral = _h_ctheta_all->integral();
185      scale(_h_ctheta_all, 1./integral);
186      scale(_h_01_all, 1./integral);
187      normalize(_h_phi_all);
188
189      pair<double,double> rho00 = calcRho(_h_ctheta_all,0);
190      h_rho00->bin(7).set(rho00.first, rho00.second);
191
192      pair<double,pair<double,double> > alpha = calcAlpha(_h_ctheta_all);
193      h_alpha->bin(7).set(alpha.first, alpha.second);
194
195      pair<double,double> rhooff = calcRho(_h_phi_all, 1);
196      h_rhooff->bin(7).set(rhooff.first, rhooff.second);
197
198      pair<double,double> rho01 = calcRho(_h_01_all, 2);
199      h_01->bin(7).set(rho01.first, rho01.second);
200    }
201
202    /// @}
203
204
205    /// @name Histograms
206    /// @{
207    Histo1DGroupPtr _h_ctheta,_h_phi, _h_01;
208    Histo1DPtr _h_ctheta_all, _h_phi_all, _h_01_all;
209    /// @}
210
211
212  };
213
214
215  RIVET_DECLARE_PLUGIN(TPC_1991_I316132);
216
217
218}