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 Transversality angles for $\eta_c\to f_2f^{(\prime)}_2$
10  class BELLE_2008_I754089 : public Analysis {
11  public:
12
13    /// Constructor
14    RIVET_DEFAULT_ANALYSIS_CTOR(BELLE_2008_I754089);
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==441);
24      declare(ufs, "UFS");
25      DecayedParticles ETA(ufs);
26      ETA.addStable(225);
27      ETA.addStable(335);
28      declare(ETA, "ETA");
29      // histos
30      for (unsigned int ix=0; ix<2; ++ix) {
31        book(_h[ix],1,1,1+ix);
32      }
33    }
34
35
36    /// Perform the per-event analysis
37    void analyze(const Event& event) {
38      DecayedParticles ETA = apply<DecayedParticles>(event, "ETA");
39      // loop over particles
40      for (unsigned int ix=0; ix<ETA.decaying().size(); ++ix) {
41        int imode=-1;
42      	if      (ETA.modeMatches(ix,2,mode1)) imode=0;
43      	else if (ETA.modeMatches(ix,2,mode2)) imode=1;
44      	else {
45          continue;
46        }
47        Vector3 trans[2];
48        bool found=true;
49        LorentzTransform boost1 = LorentzTransform::mkFrameTransformFromBeta(ETA.decaying()[ix].mom().betaVec());
50        Vector3 axis = boost1.transform(ETA.decaying()[ix].children()[0].mom()).p3().unit();
51        for (unsigned int iy=0;iy<ETA.decaying()[ix].children().size();++iy) {
52          int iMeson = ETA.decaying()[ix].children()[iy].pid()==225 ? 211 : 321;
53          if (ETA.decaying()[ix].children()[iy].children().size()!=2) found=false;
54          if (ETA.decaying()[ix].children()[iy].children()[0].pid()!=
55             -ETA.decaying()[ix].children()[iy].children()[1].pid()) found=false;
56          if (ETA.decaying()[ix].children()[iy].children()[0].abspid()!=iMeson) found=false;
57          if (!found) break;
58          FourMomentum pF2 = boost1.transform(ETA.decaying()[ix].children()[iy].mom());
59          LorentzTransform boost2 = LorentzTransform::mkFrameTransformFromBeta(pF2.betaVec());
60          FourMomentum pMeson = boost2.transform(boost1.transform(ETA.decaying()[ix].children()[iy].children()[0].mom()));
61          trans[iy] = pMeson.p3().unit()-pMeson.p3().unit().dot(axis)*axis;
62        }
63        if (!found) continue;
64        const double chi = atan2(trans[0].cross(trans[1]).dot(axis),trans[0].dot(trans[1]));
65        _h[imode]->fill(chi);
66      }
67    }
68
69
70    /// Normalise histograms etc., after the run
71    void finalize() {
72      normalize(_h, 1.0, false);
73    }
74
75    /// @}
76
77
78    /// @name Histograms
79    /// @{
80    Histo1DPtr _h[2];
81    const map<PdgId,unsigned int> mode1 = { { 225,2}};
82    const map<PdgId,unsigned int> mode2 = { { 225,1},{ 335,1}};
83    /// @}
84
85
86  };
87
88
89  RIVET_DECLARE_PLUGIN(BELLE_2008_I754089);
90
91}