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  tau polarization at 57 GeV
 10  class AMY_1990_I298238 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(AMY_1990_I298238);
 15
 16
 17    /// @name Analysis methods
 18    /// @{
 19
 20    /// Book histograms and initialise projections before the run
 21    void init() {
 22      // Initialise and register projections
 23      declare(ChargedFinalState(), "FS");
 24      declare(UnstableParticles(), "UFS");
 25      // book histos
 26      book(_h_e  ,"_t_e " , 20,-1,1);
 27      book(_h_mu ,"_t_mu" , 20,-1,1);
 28      book(_h_pi ,"_t_pi" , 20,-1,1);
 29      book(_h_rho,"_t_rho", 20,-1,1);
 30    }
 31
 32    void findTau(const Particle & p, unsigned int & nprod,
 33     		 Particles & piP,Particles & pi0, Particles & ell, Particles & nu_ell,
 34		 Particles & nu_tau) {
 35      for(const Particle & child : p.children()) {
 36	if(child.pid()==PID::ELECTRON || child.pid()==PID::MUON) {
 37	  ++nprod;
 38	  ell.push_back(child);
 39	}
 40	else if(child.pid()==PID::NU_EBAR || child.pid()==PID::NU_MUBAR) {
 41	  ++nprod;
 42	  nu_ell.push_back(child);
 43	}
 44	else if(child.pid()==PID::PIMINUS) {
 45	  ++nprod;
 46	  piP.push_back(child);
 47	}
 48	else if(child.pid()==PID::PI0) {
 49	  ++nprod;
 50	  pi0.push_back(child);
 51	}
 52	else if(child.pid()==PID::NU_TAU) {
 53	  ++nprod;
 54	  nu_tau.push_back(child);
 55	}
 56	else if(child.pid()==PID::GAMMA)
 57	  continue;
 58	else if(child.children().empty() || child.pid()==221 || child.pid()==331) {
 59	  ++nprod;
 60	}
 61	else {
 62	  findTau(child,nprod,piP,pi0,ell,nu_ell,nu_tau);
 63	}
 64      }
 65    }
 66
 67    /// Perform the per-event analysis
 68    void analyze(const Event& event) {
 69      // require 2 chanrged particles to veto hadronic events
 70      if(apply<ChargedFinalState>(event, "FS").particles().size()!=2) vetoEvent;
 71      // loop over tau leptons
 72      for(const Particle& p : apply<UnstableParticles>(event, "UFS").particles(Cuts::pid==15)) {
 73	unsigned int nprod(0);
 74	Particles piP, pi0, ell, nu_ell, nu_tau;
 75	findTau(p,nprod,piP, pi0, ell, nu_ell, nu_tau);
 76	LorentzTransform boost1 = LorentzTransform::mkFrameTransformFromBeta(p.momentum().betaVec());
 77	if(nprod==2 && nu_tau.size()==1 && piP.size()==1) {
 78	  FourMomentum pPi = boost1.transform(piP[0].momentum());
 79	  double cTheta = pPi.p3().unit().dot(p.momentum().p3().unit());
 80	  _h_pi->fill(cTheta);
 81	}
 82	else if(nprod==3 && nu_tau.size()==1 && ell.size()==1 && nu_ell.size()==1) {
 83	  if(ell[0].pid()==PID::ELECTRON) {
 84	    _h_e ->fill(2.*ell[0].momentum().t()/sqrtS());
 85	  }
 86	  else {
 87	    _h_mu->fill(2.*ell[0].momentum().t()/sqrtS());
 88	  }
 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      BinnedEstimatePtr<string> h_P;
125      book(h_P,1,1,1);
126      normalize(_h_e);
127      pair<double,double> P_e  = calcP(_h_e,1);
128      double s1 = P_e.first/sqr(P_e.second);
129      double s2 = 1./sqr(P_e.second);
130      normalize(_h_mu);
131      pair<double,double> P_mu = calcP(_h_mu,1);
132      s1 += P_mu.first/sqr(P_mu.second);
133      s2 += 1./sqr(P_mu.second);
134      normalize(_h_pi);
135      pair<double,double> P_pi = calcP(_h_pi,0);
136      s1 += P_pi.first/sqr(P_pi.second);
137      s2 += 1./sqr(P_pi.second);
138      normalize(_h_rho);
139      pair<double,double> P_rho = calcP(_h_rho,0);
140      P_rho.first  /=0.46;
141      P_rho.second /=0.46;
142      s1 += P_rho.first/sqr(P_rho.second);
143      s2 += 1./sqr(P_rho.second);
144      // average
145      h_P->bin(1).set(s1/s2, sqrt(1./s2));
146    }
147
148    /// @}
149
150
151    /// @name Histograms
152    /// @{
153    Histo1DPtr _h_e,_h_mu,_h_pi,_h_rho;
154    /// @}
155
156
157  };
158
159
160  RIVET_DECLARE_PLUGIN(AMY_1990_I298238);
161
162}