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 rho+/- and omega polarization
 11  class OPAL_2000_I502750 : public Analysis {
 12  public:
 13
 14    /// Constructor
 15    RIVET_DEFAULT_ANALYSIS_CTOR(OPAL_2000_I502750);
 16
 17
 18    /// @name Analysis methods
 19    /// @{
 20
 21    /// Book histograms and initialise projections before the run
 22    void init() {
 23
 24      // Initialise and register projections
 25      declare(Beam(), "Beams");
 26      declare(ChargedFinalState(), "FS");
 27      declare(UnstableParticles(), "UFS");
 28
 29      // Book histograms
 30      book(_h_ctheta_rho, {0.025, 0.05, 0.1, 0.15, 0.3, 0.6});
 31      book(_h_ctheta_omega, {0.025, 0.05, 0.1, 0.15, 0.3, 0.6});
 32      for (size_t i = 0; i<_h_ctheta_rho->numBins(); ++i) {
 33        book(_h_ctheta_rho->bin(i+1), "ctheta_rho_"+to_string(i), 20, -1.0, 1.0);
 34        book(_h_ctheta_omega->bin(i+1), "ctheta_omega_"+to_string(i), 20, -1.0, 1.0);
 35      }
 36      book(_h_ctheta_omega_all, "ctheta_omega_all",20,-1.,1.);
 37    }
 38
 39    pair<double,double> calcRho(Histo1DPtr hist) {
 40      if(hist->numEntries()==0.) return make_pair(0.,0.);
 41      double sum1(0.),sum2(0.);
 42      for (const auto& bin : hist->bins() ) {
 43	double Oi = bin.sumW();
 44	if(Oi==0.) continue;
 45	double ai = 0.25*(bin.xMax()*(3.-sqr(bin.xMax())) - bin.xMin()*(3.-sqr(bin.xMin())));
 46	double bi = 0.75*(bin.xMin()*(1.-sqr(bin.xMin())) - bin.xMax()*(1.-sqr(bin.xMax())));
 47	double Ei = bin.errW();
 48	sum1 += sqr(bi/Ei);
 49	sum2 += bi/sqr(Ei)*(Oi-ai);
 50      }
 51      return make_pair(sum2/sum1,sqrt(1./sum1));
 52    }
 53
 54    bool findOmegaDecay(Particle omega,Particles & pi0, Particles & pip, Particles & pim) {
 55      for(const Particle & child : omega.children()) {
 56	if(child.pid()==211)
 57	  pip.push_back(child);
 58	else if(child.pid()==-211)
 59	  pim.push_back(child);
 60	else if(child.pid()==111)
 61	  pi0.push_back(child);
 62	else if(!child.children().empty()) {
 63	  if(!findOmegaDecay(child,pi0,pip,pim)) return false;
 64	}
 65	else
 66	  return false;
 67      }
 68      return true;
 69    }
 70
 71    /// Perform the per-event analysis
 72    void analyze(const Event& event) {
 73      // First, veto on leptonic events by requiring at least 4 charged FS particles
 74      const FinalState& fs = apply<FinalState>(event, "FS");
 75      const size_t numParticles = fs.particles().size();
 76
 77      // Even if we only generate hadronic events, we still need a cut on numCharged >= 2.
 78      if (numParticles < 2) {
 79        MSG_DEBUG("Failed leptonic event cut");
 80        vetoEvent;
 81      }
 82      MSG_DEBUG("Passed leptonic event cut");
 83      // Get beams and average beam momentum
 84      const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
 85      const double meanBeamMom = ( beams.first.p3().mod() +
 86                                   beams.second.p3().mod() ) / 2.0;
 87      MSG_DEBUG("Avg beam momentum = " << meanBeamMom);
 88      // loop over rho and omega mesons
 89      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
 90      for (const Particle& p : ufs.particles(Cuts::abspid==213 || Cuts::abspid==223)) {
 91	double xE = p.momentum().t()/meanBeamMom;
 92	Vector3 e1z = p.momentum().p3().unit();
 93	LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(p.momentum().betaVec());
 94	if(p.abspid()==213) {
 95	  if(p.children().size()!=2) continue;
 96	  int sign = p.pid()/213;
 97	  Particle pion;
 98	  if(p.children()[0].pid()==sign*211 && p.children()[1].pid()==111) {
 99	    pion = p.children()[0];
100	  }
101	  else if(p.children()[1].pid()==sign*211 && p.children()[0].pid()==111) {
102	    pion = p.children()[1];
103	  }
104	  else
105	    continue;
106	  Vector3 axis1 = boost.transform(pion.momentum()).p3().unit();
107	  double ctheta = e1z.dot(axis1);
108	  _h_ctheta_rho->fill(xE,ctheta);
109	}
110	else {
111	  Particles pi0,pip,pim;
112	  bool three_pi = findOmegaDecay(p,pi0,pip,pim);
113	  if(!three_pi || pi0.size()!=1 || pip.size()!=1 || pim.size()!=1)
114	    continue;
115	  Vector3 v1 = boost.transform(pi0[0].momentum()).p3().unit();
116	  Vector3 v2 = boost.transform(pip[0].momentum()).p3().unit();
117	  Vector3 norm = v1.cross(v2).unit();
118	  double ctheta = e1z.dot(norm);
119	  _h_ctheta_omega->fill(xE,ctheta);
120	  if(xE>0.025) _h_ctheta_omega_all->fill(ctheta);
121	}
122      }
123    }
124
125
126    /// Normalise histograms etc., after the run
127    void finalize() {
128      vector<double> x = {0.025,0.05,0.1,0.15,0.3,0.6};
129      Estimate1DPtr h_rho  ;
130      book(h_rho, 1,1,1);
131      Estimate1DPtr h_omega;
132      book(h_omega, 2,1,1);
133      for (size_t ix=1; ix<_h_ctheta_rho->numBins()+1; ++ix) {
134        // rho
135        normalize(_h_ctheta_rho->bin(ix));
136        pair<double,double> rho00 = calcRho(_h_ctheta_rho->bin(ix));
137        h_rho->bin(ix).set(rho00.first, rho00.second);
138        // omega
139        normalize(_h_ctheta_omega->bin(ix));
140        rho00 = calcRho(_h_ctheta_omega->bin(ix));
141        h_omega->bin(ix).set(rho00.first, rho00.second);
142      }
143      // omega over whole range
144      Estimate1DPtr h_omega_all;
145      book(h_omega_all,2,2,1);
146      normalize(_h_ctheta_omega_all);
147      pair<double,double> rho00 = calcRho(_h_ctheta_omega_all);
148      h_omega_all->bin(1).set(rho00.first, rho00.second);
149    }
150
151    /// @}
152
153
154    /// @name Histograms
155    /// @{
156    Histo1DGroupPtr _h_ctheta_rho, _h_ctheta_omega;
157    Histo1DPtr _h_ctheta_omega_all;
158    /// @}
159
160
161  };
162
163
164  RIVET_DECLARE_PLUGIN(OPAL_2000_I502750);
165
166
167}