1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4#include "Rivet/Projections/DecayedParticles.hh"
  5
  6namespace Rivet {
  7
  8
  9  /// @brief Upsilon_2 -> pi+ pi- Upsilon
 10  class BABAR_2010_I850492 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(BABAR_2010_I850492);
 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      UnstableParticles ufs = UnstableParticles(Cuts::pid==20555);
 24      declare(ufs, "UFS");
 25      DecayedParticles Upsilon2(ufs);
 26      Upsilon2.addStable(PID::PI0);
 27      Upsilon2.addStable(553);
 28      declare(Upsilon2, "Upsilon2");
 29      for(unsigned int ix=0;ix<3;++ix)
 30	book(_h[ix],1,1,1+ix);
 31    }
 32
 33
 34    /// Perform the per-event analysis
 35    void analyze(const Event& event) {
 36      static const map<PdgId,unsigned int> & mode = { { 211,1}, {-211,1}, {553,1} };
 37      DecayedParticles Upsilon2 = apply<DecayedParticles>(event, "Upsilon2");
 38      // loop over particles
 39      for(unsigned int ix=0;ix<Upsilon2.decaying().size();++ix) {
 40	if ( !Upsilon2.modeMatches(ix,3,mode) ) continue;
 41       	const Particle & pip= Upsilon2.decayProducts()[ix].at( 211)[0];
 42       	const Particle & pim= Upsilon2.decayProducts()[ix].at(-211)[0];
 43       	const Particle & ups= Upsilon2.decayProducts()[ix].at( 553)[0];
 44	FourMomentum ptot = pip.momentum()+pim.momentum();
 45	_h[0]->fill(ptot.mass());
 46	// boost to Upsilon_2 rest frame
 47	LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(Upsilon2.decaying()[ix].momentum().betaVec());
 48	FourMomentum pDir = boost.transform(ptot);
 49	Matrix3 ptoz(-pDir.p3().unit(), Vector3(0,0,1));
 50	boost.preMult(ptoz);
 51	FourMomentum p2 = boost.transform(ups.momentum());
 52	FourMomentum ppip = boost.transform(pip.momentum());
 53	FourMomentum ppim = boost.transform(pim.momentum());
 54	ptot = ppip+ppim;
 55	// pion angle
 56	LorentzTransform boostPi = LorentzTransform::mkFrameTransformFromBeta(ptot.betaVec());
 57      	Vector3 axisPi = boostPi.transform(ppip).p3().unit();
 58	double cosPi = axisPi.dot(ptot.p3().unit());
 59      	_h[2]->fill(abs(cosPi));
 60	if(ups.children().size()!=2) continue;
 61	Particle ep,em;
 62	if ( ups.children()[0].pid()==-ups.children()[1].pid() &&
 63	     (ups.children()[0].abspid()==11 || ups.children()[0].abspid()==13)) {
 64	  ep = ups.children()[0];
 65	  em = ups.children()[1];
 66	}
 67	else
 68	  continue;
 69	if(em.pid()<0) swap(ep,em);
 70      	LorentzTransform boostUps = LorentzTransform::mkFrameTransformFromBeta(p2.betaVec());
 71	FourMomentum pe = boost.transform(ep .momentum());
 72      	Vector3 axisE = boostUps.transform(pe).p3().unit();
 73      	axisPi.setZ(0.);
 74      	axisE.setZ(0.);
 75      	double chi = abs(atan2(axisE.cross(axisPi).dot(p2.p3().unit()), axisE.dot(axisPi)));
 76	if(chi>M_PI) chi=2.*M_PI-chi;
 77      	_h[1]->fill(chi);
 78      }
 79    }
 80
 81
 82    /// Normalise histograms etc., after the run
 83    void finalize() {
 84      for(unsigned int ix=0;ix<3;++ix)
 85	normalize(_h[ix],1.,false);
 86    }
 87
 88    /// @}
 89
 90
 91    /// @name Histograms
 92    /// @{
 93    Histo1DPtr _h[3];
 94    /// @}
 95
 96
 97  };
 98
 99
100  RIVET_DECLARE_PLUGIN(BABAR_2010_I850492);
101
102}