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