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 J/psi -> K*0 Kbar*0
 10  class BES_1999_I505287 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(BES_1999_I505287);
 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      DecayedParticles PSI(ufs);
 26      PSI.addStable(PID::K0S);
 27      PSI.addStable(PID::PI0);
 28      PSI.addStable( 313);
 29      PSI.addStable(-313);
 30      declare(PSI, "PSI");
 31      // histograms
 32      book(_h_mass,1,1,1);
 33      for (unsigned int ix=0; ix<2; ++ix) {
 34        book(_h_angle[ix],2,1,1+ix);
 35      }
 36    }
 37
 38
 39    /// Perform the per-event analysis
 40    void analyze(const Event& event) {
 41      // find the J/psi decays
 42      DecayedParticles PSI = apply<DecayedParticles>(event, "PSI");
 43      // loop over particles
 44      for (unsigned int ix=0;ix<PSI.decaying().size();++ix) {
 45        if (!PSI.modeMatches(ix,3,mode)) continue;
 46        const Particle  & Kstar0  = PSI.decayProducts()[0].at( 313)[0];
 47        const Particle  & KstarB0 = PSI.decayProducts()[0].at(-313)[0];
 48        FourMomentum phad = Kstar0.mom()+KstarB0.mom();
 49        _h_mass->fill(phad.mass());
 50        Particle Kp;
 51        if (Kstar0.children()[0].pid()== 321 &&
 52            Kstar0.children()[1].pid()==-211) {
 53          Kp = Kstar0.children()[0];
 54        }
 55        else if (Kstar0.children()[1].pid()== 321 &&
 56                 Kstar0.children()[0].pid()==-211) {
 57          Kp = Kstar0.children()[1];
 58        }
 59        else {
 60          continue;
 61        }
 62        Particle Km;
 63        if (KstarB0.children()[0].pid()==-321 &&
 64            KstarB0.children()[1].pid()== 211) {
 65          Km = KstarB0.children()[0];
 66        }
 67        else if (KstarB0.children()[1].pid()==-321 &&
 68                 KstarB0.children()[0].pid()== 211) {
 69          Km = KstarB0.children()[1];
 70        }
 71        else {
 72          continue;
 73        }
 74        const LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(phad.betaVec());
 75        // K*0 decay
 76        FourMomentum pKstar0 = boost.transform(Kstar0 .mom());
 77        Vector3 axis1 = pKstar0.p3().unit();
 78        const LorentzTransform boost1 = LorentzTransform::mkFrameTransformFromBeta(pKstar0 .betaVec());
 79        FourMomentum pKp = boost1.transform(boost.transform(Kp.mom()));
 80        double cTheta1 = pKp.p3().unit().dot(axis1);
 81        _h_angle[0]->fill(cTheta1);
 82        Vector3 trans1 = pKp.p3() - cTheta1*pKp.p3().mod()*axis1;
 83        // Kbar*0 decay
 84        FourMomentum pKstarB0 = boost.transform(KstarB0.mom());
 85        Vector3 axis2 = pKstar0.p3().unit();
 86        const LorentzTransform boost2 = LorentzTransform::mkFrameTransformFromBeta(pKstarB0.betaVec());
 87        FourMomentum pKm = boost2.transform(boost.transform(Km.mom()));
 88        double cTheta2 = pKm.p3().unit().dot(axis2);
 89        _h_angle[0]->fill(cTheta2);
 90        Vector3 trans2 = pKm.p3() - cTheta1*pKm.p3().mod()*axis2;
 91        double chi = abs(atan2(trans1.cross(trans2).dot(axis1),trans1.dot(trans2)))/M_PI*180.;
 92        _h_angle[1]->fill(chi);
 93      }
 94    }
 95
 96
 97    /// Normalise histograms etc., after the run
 98    void finalize() {
 99      normalize(_h_mass, 1.0, false);
100      normalize(_h_angle, 1.0, false);
101    }
102
103    /// @}
104
105
106    /// @name Histograms
107    /// @{
108    Histo1DPtr _h_mass, _h_angle[2];
109    const map<PdgId,unsigned int> mode = { { 313,1}, {-313,1}, { 22,1}};
110    /// @}
111
112
113  };
114
115
116  RIVET_DECLARE_PLUGIN(BES_1999_I505287);
117
118}