1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/Beam.hh"
  4#include "Rivet/Projections/UnstableParticles.hh"
  5#include "Rivet/Projections/DecayedParticles.hh"
  6
  7namespace Rivet {
  8
  9
 10  /// @brief J/psi -> omega K+K-
 11  class BESII_2004_I658084 : public Analysis {
 12  public:
 13
 14    /// Constructor
 15    RIVET_DEFAULT_ANALYSIS_CTOR(BESII_2004_I658084);
 16
 17
 18    /// @name Analysis methods
 19    /// @{
 20
 21    /// Book histograms and initialise projections before the run
 22    void init() {
 23      // Initialise and register projections
 24      UnstableParticles ufs = UnstableParticles(Cuts::abspid==443);
 25      declare(ufs, "UFS");
 26      DecayedParticles PSI(ufs);
 27      PSI.addStable(223);
 28      declare(PSI, "PSI");
 29      declare(Beam(), "Beams");
 30      // histos
 31      for(unsigned int ix=0;ix<2;++ix)
 32	book(_h_mass[ix],1,1,1+ix);
 33      for(unsigned int ix=0;ix<4;++ix)
 34	book(_h_angle[ix],2,1,1+ix);
 35    }
 36
 37    void findChildren(const Particle& p, Particles& pim, Particles& pip,
 38                      Particles & pi0, unsigned int &ncount) {
 39      for (const Particle& child : p.children()) {
 40        if (child.pid()==PID::PIPLUS) {
 41          pip.push_back(child);
 42          ncount+=1;
 43        }
 44        else if (child.pid()==PID::PIMINUS) {
 45          pim.push_back(child);
 46          ncount+=1;
 47        }
 48        else if (child.pid()==PID::PI0) {
 49          pi0.push_back(child);
 50          ncount+=1;
 51        }
 52        else if (child.children().empty()) {
 53          ncount+=1;
 54        }
 55        else {
 56          findChildren(child,pim,pip,pi0,ncount);
 57        }
 58      }
 59    }
 60
 61    /// Perform the per-event analysis
 62    void analyze(const Event& event) {
 63      // get the axis, direction of incoming electron
 64      const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
 65      Vector3 axis;
 66      if (beams.first.pid()>0) {
 67        axis = beams.first .mom().p3().unit();
 68      }
 69      else {
 70        axis = beams.second.mom().p3().unit();
 71      }
 72      // J/psi decay products
 73      DecayedParticles PSI = apply<DecayedParticles>(event, "PSI");
 74      if ( PSI.decaying().size()!=1) vetoEvent;
 75      if (!PSI.modeMatches(0,3,mode)) vetoEvent;
 76      const Particle& Kp    = PSI.decayProducts()[0].at( 321)[0];
 77      const Particle& Km    = PSI.decayProducts()[0].at(-321)[0];
 78      const Particle& omega = PSI.decayProducts()[0].at( 223)[0];
 79      // mass histograms
 80      FourMomentum pKK = Kp.mom()+Km   .mom();
 81      double mKK = pKK.mass();
 82      _h_mass[0]->fill(mKK/GeV);
 83      _h_mass[1]->fill((Kp.mom()+omega.mom()).mass()/GeV);
 84      _h_mass[1]->fill((Km.mom()+omega.mom()).mass()/GeV);
 85      if (mKK<2.) return;
 86      if (abs(Kp.p3().unit().dot(axis))>0.84) return;
 87      if (abs(Km.p3().unit().dot(axis))>0.84) return;
 88      LorentzTransform boost;
 89      if (PSI.decaying()[0].mom().p3().mod() > 1*MeV) {
 90        boost = LorentzTransform::mkFrameTransformFromBeta(PSI.decaying()[0].mom().betaVec());
 91      }
 92      pKK = boost.transform(pKK);
 93      // omega decay
 94      unsigned int ncount=0;
 95      Particles pip,pim,pi0;
 96      findChildren(omega,pim,pip,pi0,ncount);
 97      if ( ncount!=3 || !(pim.size()==1 && pip.size()==1 && pi0.size()==1)) return;
 98      if (abs(pip[0].p3().unit().dot(axis))>0.84) return;
 99      if (abs(pim[0].p3().unit().dot(axis))>0.84) return;
100      FourMomentum pOmega = boost.transform(omega.mom());
101      Vector3 e1Z = pOmega.p3().unit();
102      Vector3 e1Y = e1Z.cross(axis).unit();
103      Vector3 e1X = e1Y.cross(e1Z).unit();
104      const double cOmega = e1Z.dot(axis);
105      _h_angle[1]->fill(cOmega);
106      // kaon angles
107      LorentzTransform boost2 = LorentzTransform::mkFrameTransformFromBeta(pKK.betaVec());
108      FourMomentum pK = boost2.transform(boost.transform(Kp.mom()));
109      const double cK = pK.p3().unit().dot(e1Z);
110      _h_angle[2]->fill(cK);
111      Vector3 transK = pK.p3().unit()-cK*e1Z;
112      // omega angles
113      LorentzTransform boost3 = LorentzTransform::mkFrameTransformFromBeta(pOmega.betaVec());
114      FourMomentum ppip = boost3.transform(boost.transform(pip[0].mom()));
115      FourMomentum ppim = boost3.transform(boost.transform(pim[0].mom()));
116      Vector3 nW = ppip.p3().cross(ppim.p3()).unit();
117      const double bW = nW.dot(e1X);
118      _h_angle[3]->fill(bW);
119      Vector3 transPi = ppip.p3()-ppip.p3().dot(e1Z)*e1Z;
120      const double chi = abs(atan2(transK.cross(transPi).dot(e1Z),transK.dot(transPi)));
121      _h_angle[0]->fill(chi/M_PI*180.);
122    }
123
124
125    /// Normalise histograms etc., after the run
126    void finalize() {
127      normalize(_h_mass, 1.0, false);
128      normalize(_h_angle, 1.0, false);
129    }
130
131    /// @}
132
133
134    /// @name Histograms
135    /// @{
136    Histo1DPtr _h_mass[2],_h_angle[4];
137    const map<PdgId,unsigned int> mode = { { 223,1}, { 321,1},{-321,1}};
138    /// @}
139
140
141  };
142
143
144  RIVET_DECLARE_PLUGIN(BESII_2004_I658084);
145
146}