1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/Beam.hh"
  4#include "Rivet/Projections/ChargedFinalState.hh"
  5#include "Rivet/Projections/UnstableParticles.hh"
  6
  7namespace Rivet {
  8
  9
 10  /// @brief  e+e- > tau+ tau-
 11  class ALEPH_2001_I555653 : public Analysis {
 12  public:
 13
 14    /// Constructor
 15    RIVET_DEFAULT_ANALYSIS_CTOR(ALEPH_2001_I555653);
 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(ChargedFinalState(), "FS");
 26      declare(UnstableParticles(), "UFS");
 27
 28      // book hists
 29      book(_h_e,   {-0.9, -0.7, -0.5, -0.3, -0.1, 0.1, 0.3, 0.5, 0.7, 0.9 });
 30      book(_h_mu,  {-0.9, -0.7, -0.5, -0.3, -0.1, 0.1, 0.3, 0.5, 0.7, 0.9 });
 31      book(_h_pi,  {-0.9, -0.7, -0.5, -0.3, -0.1, 0.1, 0.3, 0.5, 0.7, 0.9 });
 32      book(_h_rho, {-0.9, -0.7, -0.5, -0.3, -0.1, 0.1, 0.3, 0.5, 0.7, 0.9 });
 33      for (size_t ix=0; ix < _h_e->numBins(); ++ix) {
 34        const string suff = to_string(ix);
 35        book(_h_e->bin(ix+1),   "_h_e_"+suff,   20, -1.0, 1.0);
 36        book(_h_mu->bin(ix+1),  "_h_mu_"+suff,  20, -1.0, 1.0);
 37        book(_h_pi->bin(ix+1),  "_h_pi_"+suff,  20, -1.0, 1.0);
 38        book(_h_rho->bin(ix+1), "_h_rho_"+suff, 20, -1.0, 1.0);
 39      }
 40    }
 41
 42    void findTau(const Particle & p, unsigned int & nprod, Particles& piP,
 43     		         Particles& pi0, Particles& ell, Particles& nu_ell, Particles& nu_tau) {
 44      for (const Particle& child : p.children()) {
 45        if (child.pid()==PID::ELECTRON || child.pid()==PID::MUON) {
 46          ++nprod;
 47          ell.push_back(child);
 48        }
 49        else if (child.pid()==PID::NU_EBAR || child.pid()==PID::NU_MUBAR) {
 50          ++nprod;
 51          nu_ell.push_back(child);
 52        }
 53        else if (child.pid()==PID::PIMINUS) {
 54          ++nprod;
 55          piP.push_back(child);
 56        }
 57        else if (child.pid()==PID::PI0) {
 58          ++nprod;
 59          pi0.push_back(child);
 60        }
 61        else if (child.pid()==PID::NU_TAU) {
 62          ++nprod;
 63          nu_tau.push_back(child);
 64        }
 65        else if (child.pid()==PID::GAMMA) {
 66          continue;
 67        }
 68        else if (child.children().empty() || child.pid()==221 || child.pid()==331) {
 69          ++nprod;
 70        }
 71        else {
 72          findTau(child,nprod,piP,pi0,ell,nu_ell,nu_tau);
 73        }
 74      }
 75    }
 76
 77    /// Perform the per-event analysis
 78    void analyze(const Event& event) {
 79      // require 2 chanrged particles to veto hadronic events
 80      if (apply<ChargedFinalState>(event, "FS").particles().size()!=2) vetoEvent;
 81      // Get beams and average beam momentum
 82      const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
 83      Vector3 axis;
 84      if (beams.first.pid()>0) {
 85        axis = beams.first .momentum().p3().unit();
 86      }
 87      else {
 88        axis = beams.second.momentum().p3().unit();
 89      }
 90      // loop over tau leptons
 91      for (const Particle& p : apply<UnstableParticles>(event, "UFS").particles(Cuts::pid==15)) {
 92        unsigned int nprod(0);
 93        Particles piP, pi0, ell, nu_ell, nu_tau;
 94        findTau(p,nprod,piP, pi0, ell, nu_ell, nu_tau);
 95        LorentzTransform boost1 = LorentzTransform::mkFrameTransformFromBeta(p.momentum().betaVec());
 96        double cBeam = axis.dot(p.momentum().p3().unit());
 97        if (nprod==2 && nu_tau.size()==1 && piP.size()==1) {
 98          FourMomentum pPi = boost1.transform(piP[0].momentum());
 99          double cTheta = pPi.p3().unit().dot(p.momentum().p3().unit());
100          _h_pi->fill(cBeam,cTheta);
101        }
102        else if (nprod==3 && nu_tau.size()==1 && ell.size()==1 && nu_ell.size()==1) {
103          if(ell[0].pid()==PID::ELECTRON) {
104            _h_e->fill(cBeam,2.*ell[0].momentum().t()/sqrtS());
105          }
106          else {
107            _h_mu->fill(cBeam,2.*ell[0].momentum().t()/sqrtS());
108          }
109        }
110        else if (nprod==3 && nu_tau.size()==1 && piP.size()==1&& pi0.size()==1) {
111          FourMomentum pRho = boost1.transform(piP[0].momentum()+pi0[0].momentum());
112          double cTheta = pRho.p3().unit().dot(p.momentum().p3().unit());
113          _h_rho->fill(cBeam, cTheta);
114        }
115      }
116    }
117
118    pair<double,double> calcP(Histo1DPtr hist,unsigned int imode) {
119      if (hist->numEntries()==0.) return make_pair(0.,0.);
120      double sum1(0.),sum2(0.);
121      for (const auto& bin : hist->bins()) {
122        double Oi = bin.sumW();
123        if (Oi==0.) continue;
124        double ai(0.),bi(0.);
125        // tau -> pi/rho nu
126        if (imode==0) {
127          ai = 0.5*(bin.xMax()-bin.xMin());
128          bi = 0.5*ai*(bin.xMax()+bin.xMin());
129        }
130        // lepton mode
131        else {
132          ai = (-5*bin.xMin() + 3*pow(bin.xMin(),3) -   pow(bin.xMin(),4) + 5*bin.xMax() - 3*pow(bin.xMax(),3) +   pow(bin.xMax(),4))/3.;
133          bi = (  -bin.xMin() + 3*pow(bin.xMin(),3) - 2*pow(bin.xMin(),4) +   bin.xMax() - 3*pow(bin.xMax(),3) + 2*pow(bin.xMax(),4))/3.;
134        }
135        double Ei = bin.errW();
136        sum1 += sqr(bi/Ei);
137        sum2 += bi/sqr(Ei)*(Oi-ai);
138      }
139      return make_pair(sum2/sum1,sqrt(1./sum1));
140    }
141
142    /// Normalise histograms etc., after the run
143    void finalize() {
144      Estimate1DPtr _h_P;
145      book(_h_P,2,1,1);
146      for (size_t ix=1; ix < _h_e->numBins()+1; ++ix) {
147        normalize(_h_e->bin(ix));
148        pair<double,double> P_e  = calcP(_h_e->bin(ix), 1);
149        double s1 = P_e.first/sqr(P_e.second);
150        double s2 = 1./sqr(P_e.second);
151
152        normalize(_h_mu->bin(ix));
153        pair<double,double> P_mu = calcP(_h_mu->bin(ix), 1);
154        s1 += P_mu.first/sqr(P_mu.second);
155        s2 += 1./sqr(P_mu.second);
156
157        normalize(_h_pi->bin(ix));
158        pair<double,double> P_pi = calcP(_h_pi->bin(ix), 0);
159        s1 += P_pi.first/sqr(P_pi.second);
160        s2 += 1./sqr(P_pi.second);
161
162        normalize(_h_rho->bin(ix));
163        pair<double,double> P_rho = calcP(_h_rho->bin(ix), 0);
164        P_rho.first  /=0.46;
165        P_rho.second /=0.46;
166        s1 += P_rho.first/sqr(P_rho.second);
167        s2 += 1./sqr(P_rho.second);
168        // average
169        _h_P->bin(ix).set(s1/s2, sqrt(1./s2));
170      }
171    }
172
173    /// @}
174
175
176    /// @name Histograms
177    /// @{
178    Histo1DGroupPtr _h_e, _h_mu, _h_pi, _h_rho;
179    /// @}
180
181
182  };
183
184
185  RIVET_DECLARE_PLUGIN(ALEPH_2001_I555653);
186
187}