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 psi(2S) -> J/psi gamma gamma
11  class BESIII_2012_I1097066 : public Analysis {
12  public:
13
14    /// Constructor
15    RIVET_DEFAULT_ANALYSIS_CTOR(BESIII_2012_I1097066);
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==100443);
25      declare(ufs, "UFS");
26      DecayedParticles PSI(ufs);
27      PSI.addStable(PID::JPSI);
28      declare(PSI, "PSI");
29      declare(Beam(), "Beams");
30      // book histograms
31      for(unsigned int ix=0;ix<2;++ix)
32	book(_h[ix],1,1,1+ix);
33    }
34
35
36    /// Perform the per-event analysis
37    void analyze(const Event& event) {
38      // get the axis, direction of incoming electron
39      const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
40      Vector3 axis;
41      if(beams.first.pid()>0) axis = beams.first .momentum().p3().unit();
42      else                    axis = beams.second.momentum().p3().unit();
43      // find the J/psi decays
44      static const map<PdgId,unsigned int> & mode = { { 443,1},{ 22,2}};
45      DecayedParticles PSI = apply<DecayedParticles>(event, "PSI");
46      if( PSI.decaying().size()!=1) vetoEvent;
47      if(!PSI.modeMatches(0,3,mode)) vetoEvent;
48      if( PSI.decaying()[0].children().size()!=3) vetoEvent;
49      // particles
50      const Particle  & JPsi = PSI.decayProducts()[0].at(443)[0];
51      const Particles & gam  = PSI.decayProducts()[0].at( 22);
52      Vector3 norm = gam[0].p3().cross(gam[1].p3()).unit();
53      _h[0]->fill(abs(axis.dot(norm)));
54      if(JPsi.children().size()!=2) vetoEvent;
55      if(JPsi.children()[0].pid()!=-JPsi.children()[1].pid()) vetoEvent;
56      if(JPsi.children()[0].abspid()!=PID::EMINUS &&
57	 JPsi.children()[0].abspid()!=PID::MUON) vetoEvent;
58      Particle lm = JPsi.children()[0];
59      Particle lp = JPsi.children()[1];
60      if(lm.pid()<0) swap(lm,lp);
61      Vector3 axis2=JPsi.momentum().p3().unit();
62      LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(JPsi.momentum().betaVec());
63      Vector3 axis3 = boost.transform(lm.momentum()).p3().unit();
64      _h[1]->fill(abs(axis2.dot(axis3)));
65    }
66
67
68    /// Normalise histograms etc., after the run
69    void finalize() {
70      for(unsigned int ix=0;ix<2;++ix)
71	normalize(_h[ix],1.,false);
72    }
73
74    /// @}
75
76
77    /// @name Histograms
78    /// @{
79    Histo1DPtr _h[2];
80    /// @}
81
82
83  };
84
85
86  RIVET_DECLARE_PLUGIN(BESIII_2012_I1097066);
87
88}