1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// @brief Helicity angle distributions in excited $D$ meson decays
  9  class LHCB_2013_I1243156 : public Analysis {
 10  public:
 11
 12    /// Constructor
 13    RIVET_DEFAULT_ANALYSIS_CTOR(LHCB_2013_I1243156);
 14
 15
 16    /// @name Analysis methods
 17    /// @{
 18
 19    /// Book histograms and initialise projections before the run
 20    void init() {
 21
 22      // Initialise and register projections
 23      declare(UnstableParticles(), "UFS");
 24
 25      // Book histograms
 26      book(_h_D1_ctheta     , 8,1,1);
 27      book(_h_D2_ctheta     , 8,1,2);
 28      book(_h_DStar2S_ctheta, 9,1,1);
 29      book(_h_D3_ctheta     , 9,1,2);
 30      book(_h_D2S_ctheta    ,10,1,1);
 31      book(_h_D21D_ctheta   ,10,1,2);
 32
 33    }
 34
 35
 36    /// Recursively walk the decay tree to find decay products of @a p
 37    void findDecayProducts(Particle mother, Particles & dstar, Particles & d0, Particles & pi,unsigned int & ncount) {
 38      for(const Particle & p: mother.children()) {
 39        if(p.abspid()==413)
 40          dstar.push_back(p);
 41        else if(p.abspid()==421)
 42          d0.push_back(p);
 43        else if(p.abspid()==211)
 44          pi.push_back(p);
 45        ncount +=1;
 46      }
 47    }
 48
 49
 50    /// Perform the per-event analysis
 51    void analyze(const Event& event) {
 52      const Particles usps = apply<UnstableParticles>(event, "UFS")
 53             .particles(Cuts::abspid==425 || Cuts::abspid==10423 || Cuts::abspid==100423 || Cuts::abspid==427||
 54                        Cuts::abspid==100421 || Cuts::abspid==10425 || Cuts::abspid==20425);
 55      for (const Particle& p : usps) {
 56        // decay products
 57        Particles dstar,d0,pi;
 58        unsigned int ncount=0;
 59        findDecayProducts(p,dstar,d0, pi,ncount);
 60        if (ncount!=2 || dstar.size()!=1 || pi.size()!=1 || d0.size()!=0 ) continue;
 61        if (dstar[0].pid()/p.pid()<0) continue;
 62        Particle p2 = dstar[0];
 63        LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(p2.momentum().betaVec());
 64        Vector3 d1 = boost.transform(pi[0].momentum()).p3().unit();
 65        ncount=0;
 66        dstar.clear();
 67        d0.clear();
 68        pi.clear();
 69        findDecayProducts(p2,dstar,d0, pi,ncount);
 70        if (ncount!=2 || dstar.size()!=0 || pi.size()!=1 || d0.size()!=1 ) continue;
 71        if (pi[0].pid()/p2.pid()<0) continue;
 72        Vector3 d2 = boost.transform(pi[0].momentum()).p3().unit();
 73        double cTheta  = d1.dot(d2);
 74        // decay angles
 75        if (p.abspid()==425)
 76          _h_D2_ctheta->fill(cTheta);
 77        else if (p.abspid()==10423)
 78          _h_D1_ctheta->fill(cTheta);
 79        else if (p.abspid()==100423)
 80          _h_DStar2S_ctheta->fill(cTheta);
 81        else if (p.abspid()==427)
 82          _h_D3_ctheta->fill(cTheta);
 83        else if (p.abspid()==100421)
 84          _h_D2S_ctheta->fill(cTheta);
 85        else if (p.abspid()==10425 || p.abspid()==20425)
 86          _h_D21D_ctheta->fill(cTheta);
 87      }
 88    }
 89
 90
 91    /// Normalise histograms etc., after the run
 92    void finalize() {
 93      normalize(_h_D1_ctheta);
 94      normalize(_h_D2_ctheta);
 95      normalize(_h_DStar2S_ctheta);
 96      normalize(_h_D3_ctheta     );
 97      normalize(_h_D2S_ctheta    );
 98      normalize(_h_D21D_ctheta   );
 99    }
100
101    /// @}
102
103
104    /// @name Histograms
105    /// @{
106    Histo1DPtr _h_D1_ctheta,_h_D2_ctheta;
107    Histo1DPtr _h_DStar2S_ctheta,_h_D3_ctheta;
108    Histo1DPtr _h_D2S_ctheta,_h_D21D_ctheta;
109    /// @}
110
111  };
112
113
114
115  RIVET_DECLARE_PLUGIN(LHCB_2013_I1243156);
116
117}