1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// @brief Add a short analysis description here
  9  class MC_Onium_PiPi_Decay : public Analysis {
 10  public:
 11
 12    /// Constructor
 13    RIVET_DEFAULT_ANALYSIS_CTOR(MC_Onium_PiPi_Decay);
 14
 15
 16    /// @name Analysis methods
 17    //@{
 18
 19    /// Book histograms and initialise projections before the run
 20    void init() {
 21
 22      // Initialise and register projections
 23      declare(UnstableParticles(),"UFS");
 24      // psi 2S
 25      bookHistos(100443,   443,0.6);
 26      // psi(3770)
 27      bookHistos( 30443,   443,0.7);
 28      // Upsilon (4S)
 29      bookHistos(300553,   553,1.2);
 30      bookHistos(300553,100553,0.6);
 31      // Upsilon (3S)
 32      bookHistos(200553,   553,0.9);
 33      bookHistos(200553,100553,0.4);
 34      // Upsilon (2S)
 35      bookHistos(100553,   553,0.6);
 36    }
 37
 38    void bookHistos(int id1, int id2, double deltaM) {
 39      double twompi = 0.378;
 40      _incoming.push_back(id1);
 41      _outgoing.push_back(id2);
 42      std::ostringstream title;
 43      title << "h_" << id1 << "_" << id2 << "_";
 44      _mpipi.push_back(make_pair(Histo1DPtr(), Histo1DPtr()));
 45      book(_mpipi.back().first, title.str()+"mpippim",100,twompi/GeV,deltaM/GeV);
 46      book(_mpipi.back().second, title.str()+"mpi0pi0",100,twompi/GeV,deltaM/GeV);
 47      _hel.push_back(make_pair(Histo1DPtr(), Histo1DPtr()));
 48      book(_hel.back().first, title.str()+"hpippim",100,-1.,1.);
 49      book(_hel.back().second, title.str()+"hpi0pi0",100, 0.,1.);
 50    }
 51    
 52    void findDecayProducts(const Particle & mother,
 53			   unsigned int & nstable,
 54			   Particles& pip, Particles& pim,
 55			   Particles& pi0, Particles & onium) {
 56      for(const Particle & p : mother.children()) {
 57        int id = p.pid();
 58      	if ( id == PID::PIMINUS) {
 59	  pim.push_back(p);
 60	  ++nstable;
 61	}
 62       	else if (id == PID::PIPLUS) {
 63       	  pip.push_back(p);
 64       	  ++nstable;
 65       	}
 66       	else if (id == PID::PI0) {
 67       	  pi0.push_back(p);
 68       	  ++nstable;
 69       	}
 70	else if (abs(id)%1000==443 || abs(id)%1000==553) {
 71	  onium.push_back(p);
 72	  ++nstable;
 73	}
 74	else if ( !p.children().empty() ) {
 75	  findDecayProducts(p,nstable,pip,pim,pi0,onium);
 76	}
 77	else
 78	  ++nstable;
 79      }
 80    }
 81
 82    /// Perform the per-event analysis
 83    void analyze(const Event& event) {
 84      // loop over unstable particles
 85      for(const Particle& vMeson : apply<UnstableParticles>(event, "UFS").particles()) {
 86	int id = vMeson.pid();
 87	if(id%1000!=443 && id%1000!=553) continue;
 88	unsigned int nstable(0);
 89	Particles pip, pim, pi0, onium;
 90	findDecayProducts(vMeson,nstable,pip,pim,pi0,onium);
 91	// check for onium
 92	if(onium.size() !=1 || nstable !=3) continue;
 93	// check for pipi
 94	if( ! ((pip.size()==1 && pim.size() ==1) || pi0.size()==2)) continue;
 95	// check if histos already made
 96	unsigned int iloc=0; bool found(false);
 97	while(!found&&iloc<_incoming.size()) {
 98	  if(_incoming[iloc]==vMeson.pid()&&_outgoing[iloc]==onium[0].pid()) found=true;
 99	  else ++iloc;
100	}
101	// if histos not made, make them
102	if(!found) {
103	  MSG_WARNING("MC_Onium_PiPi_Decay" << vMeson.pid() << " " << onium[0].pid() << " "
104		      << vMeson.mass()-onium[0].mass() << "\n");
105	  continue;
106	}
107	// boost to rest frame of the pion pair
108	FourMomentum q = vMeson.momentum()-onium[0].momentum();
109	LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(q.betaVec());
110	FourMomentum qp = onium[0].momentum();
111	FourMomentum ppi= pip.size()==1 ? pip[0].momentum() : pi0[0].momentum();
112	qp  = boost.transform(qp);
113	ppi = boost.transform(ppi);
114	double cX=-ppi.p3().unit().dot(qp.p3().unit());
115	if(pi0.size()==2) {
116	  _mpipi[iloc].second->fill(q.mass());
117	  _hel  [iloc].second->fill(abs(cX));
118	}
119	else {
120	  _mpipi[iloc].first->fill(q.mass());
121	  _hel  [iloc].first->fill(cX);
122	}
123      }
124    }
125    
126    
127    /// Normalise histograms etc., after the run
128    void finalize() {
129
130      // normalize to unity
131      for(unsigned int ix=0;ix<_mpipi.size();++ix) {
132	normalize(_mpipi[ix].first );
133	normalize(_mpipi[ix].second);
134	normalize(_hel[ix].first );
135	normalize(_hel[ix].second);
136      }
137    }
138
139    //@}
140
141    /**
142     *  Incoming onium states
143     */
144    vector<long> _incoming;
145    
146    /**
147     *  Outgoing onium states
148     */
149    vector<long> _outgoing;
150    
151    /**
152     *  Histograms for the \f$\pi^+\pi^-\f$ masses
153     */
154    vector<pair<Histo1DPtr,Histo1DPtr> > _mpipi;
155    
156    /**
157     *  Histmgrams for the helicity angles
158     */
159    vector<pair<Histo1DPtr,Histo1DPtr> > _hel;
160
161  };
162
163
164  // The hook for the plugin system
165  RIVET_DECLARE_PLUGIN(MC_Onium_PiPi_Decay);
166
167
168}