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