1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// @brief psi(2S) -> J/Psi pi+pi-
  9  class BES_1999_I507637 : public Analysis {
 10  public:
 11
 12    /// Constructor
 13    RIVET_DEFAULT_ANALYSIS_CTOR(BES_1999_I507637);
 14
 15
 16    /// @name Analysis methods
 17    /// @{
 18
 19    /// Book histograms and initialise projections before the run
 20    void init() {
 21      // Initialise and register projections
 22      declare(UnstableParticles(),"UFS");
 23      book(_h_mpipi,1,1,1);
 24      book(_h_cosl ,1,1,2);
 25      book(_h_cosX ,1,1,3);
 26      book(_h_cospi,1,1,4);
 27    }
 28
 29    void findDecayProducts(const Particle & mother,
 30                           unsigned int & nstable,
 31                           Particles& pip, Particles& pim,
 32                           Particles & onium) {
 33      for (const Particle & p : mother.children()) {
 34        int id = p.pid();
 35        if ( id == PID::PIMINUS) {
 36          pim.push_back(p);
 37          ++nstable;
 38        }
 39        else if (id == PID::PIPLUS) {
 40          pip.push_back(p);
 41          ++nstable;
 42        }
 43        else if (id==443) {
 44          onium.push_back(p);
 45          ++nstable;
 46        }
 47        else if ( !p.children().empty() ) {
 48          findDecayProducts(p,nstable,pip,pim,onium);
 49        }
 50        else
 51          ++nstable;
 52      }
 53    }
 54
 55    void findLeptons(const Particle & mother,
 56                     unsigned int & nstable,
 57                     Particles& lp, Particles& lm) {
 58      for (const Particle & p : mother.children()) {
 59        int id = p.pid();
 60        if ( id == 11 || id == 13 ) {
 61          lm.push_back(p);
 62          ++nstable;
 63        }
 64        else if (id == -11 || id==-13) {
 65          lp.push_back(p);
 66          ++nstable;
 67        }
 68        else if ( !p.children().empty() ) {
 69          findLeptons(p,nstable,lp,lm);
 70        }
 71        else
 72          ++nstable;
 73      }
 74    }
 75
 76    /// Perform the per-event analysis
 77    void analyze(const Event& event) {
 78      // loop over unstable particles
 79      for (const Particle& psi : apply<UnstableParticles>(event, "UFS").particles(Cuts::pid==100443)) {
 80        unsigned int nstable(0);
 81        Particles pip, pim, onium;
 82        findDecayProducts(psi,nstable,pip,pim,onium);
 83        // check for onium
 84        if (onium.size() !=1 || nstable !=3 ||
 85           pip.size()!=1 || pim.size() !=1 ) continue;
 86        FourMomentum q = pip[0].momentum()+pim[0].momentum();
 87        _h_mpipi->fill(q.mass());
 88        _h_cosX->fill(cos(q.polarAngle()));
 89        // leptons from J/psi decay
 90        nstable = 0;
 91        Particles lp, lm;
 92        findLeptons(onium[0],nstable,lp,lm);
 93        if (nstable==2&&lp.size()==1&&lm.size()==1) {
 94          LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(onium[0].momentum().betaVec());
 95          FourMomentum pl = boost.transform(lp[0].momentum());
 96          _h_cosl->fill(cos(pl.polarAngle()));
 97        }
 98        // pions in rest frame
 99        LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(q.betaVec());
100        FourMomentum ppi = boost.transform(pip[0].momentum());
101        Vector3 axis1 = q.p3().unit();
102        _h_cospi->fill(axis1.dot(ppi.p3().unit()));
103      }
104
105    }
106
107
108    /// Normalise histograms etc., after the run
109    void finalize() {
110      normalize(_h_mpipi);
111      normalize(_h_cosl ,1.,false);
112      normalize(_h_cosX );
113      normalize(_h_cospi);
114    }
115
116    /// @}
117
118
119    /// @name Histograms
120    /// @{
121    Histo1DPtr _h_mpipi,_h_cosl,_h_cosX,_h_cospi;
122    /// @}
123
124
125  };
126
127
128  RIVET_DECLARE_PLUGIN(BES_1999_I507637);
129
130}