1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/Beam.hh"
  4#include "Rivet/Projections/UnstableParticles.hh"
  5
  6namespace Rivet {
  7
  8
  9  /// @brief J/psi -> eta(1405)
 10  class BESIII_2012_I1084539 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(BESIII_2012_I1084539);
 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::abspid==443);
 24      declare(ufs, "UFS");
 25      declare(Beam(), "Beams");
 26      // histos
 27      for(unsigned int ix=0;ix<2;++ix)
 28	book(_h[ix],1,1,1+ix);
 29    }
 30
 31
 32    /// Perform the per-event analysis
 33    void analyze(const Event& event) {
 34      // get the axis, direction of incoming electron
 35      const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
 36      Vector3 axis;
 37      if(beams.first.pid()>0)
 38	axis = beams.first .momentum().p3().unit();
 39      else
 40	axis = beams.second.momentum().p3().unit();
 41      Particles psi = apply<UnstableParticles>(event,"UFS").particles();
 42      if(psi.size()!=1 || psi[0].children().size()!=2) vetoEvent;
 43      Particle gamma,eta;
 44      if(psi[0].children()[0].pid()==PID::GAMMA &&
 45	 psi[0].children()[1].pid()==9020221) {
 46	gamma = psi[0].children()[0];
 47	eta   = psi[0].children()[1];
 48      }
 49      else if(psi[0].children()[1].pid()==PID::GAMMA &&
 50	      psi[0].children()[0].pid()==9020221) {
 51	gamma = psi[0].children()[1];
 52	eta   = psi[0].children()[0];
 53      }
 54      else
 55	vetoEvent;
 56      if(eta.children().size()!=2) vetoEvent;
 57      // eta(1405) -> f0(980) pi0
 58      Particle f0,pi0;
 59      if(eta.children()[0].pid()==PID::PI0 &&
 60	 eta.children()[1].pid()==9010221) {
 61	pi0 = eta.children()[0];
 62	f0  = eta.children()[1];
 63      }
 64      else if(eta.children()[1].pid()==PID::PI0 &&
 65	      eta.children()[0].pid()==9010221) {
 66	pi0 = eta.children()[1];
 67	f0  = eta.children()[0];
 68      }
 69      else
 70	vetoEvent;
 71      double cTheta = axis.dot(gamma.p3().unit());
 72      _h[1]->fill(cTheta);
 73      LorentzTransform boost1 = LorentzTransform::mkFrameTransformFromBeta(psi[0].momentum().betaVec());
 74      FourMomentum pGamma = boost1.transform(gamma.momentum());
 75      Vector3 axis1 = pGamma.p3().unit();
 76      FourMomentum pEta = boost1.transform(eta.momentum());
 77      LorentzTransform boost2 = LorentzTransform::mkFrameTransformFromBeta(pEta.betaVec());
 78      FourMomentum pf0 = boost2.transform(boost1.transform(f0.momentum()));
 79      _h[0]->fill(pf0.p3().unit().dot(axis1));
 80    }
 81
 82
 83    /// Normalise histograms etc., after the run
 84    void finalize() {
 85      for(unsigned int ix=0;ix<2;++ix)
 86	normalize(_h[ix],1.,false);
 87    }
 88
 89    /// @}
 90
 91
 92    /// @name Histograms
 93    /// @{
 94    Histo1DPtr _h[2];
 95    /// @}
 96
 97
 98  };
 99
100
101  RIVET_DECLARE_PLUGIN(BESIII_2012_I1084539);
102
103}