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