1// -*- C++ -*-
 2#include "Rivet/Analysis.hh"
 3#include "Rivet/Projections/UnstableParticles.hh"
 4
 5namespace Rivet {
 6
 7
 8  /// @brief D*2 spectrum and decay angle
 9  class ARGUS_1989_I268577 : public Analysis {
10  public:
11
12    /// Constructor
13    RIVET_DEFAULT_ANALYSIS_CTOR(ARGUS_1989_I268577);
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_rate     , 1, 1, 1);
27      book(_h_D2_x     , 2, 1, 1);
28      book(_h_D2_ctheta, 3, 1, 1);
29    }
30
31
32   /// Recursively walk the decay tree to find decay products of @a p
33   void findDecayProducts(Particle mother, Particles & d, Particles & pi, unsigned int & ncount) {
34     for (const Particle & p: mother.children()) {
35       if (p.abspid()==411)
36         d.push_back(p);
37       else if (p.abspid()==211)
38         pi.push_back(p);
39       ncount +=1;
40     }
41   }
42
43    /// Perform the per-event analysis
44    void analyze(const Event& event) {
45
46      if (sedges.empty()) {
47        sedges = _h_D2_x->xEdges();
48      }
49
50      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
51      for (const Particle& p : ufs.particles(Cuts::abspid==425)) {
52        const double xp = 2.*p.p3().mod()/sqrtS();
53        const size_t idx = axisMap.index(xp)-1;
54        _h_D2_x->fill(sedges[idx]);
55        // decay products
56        Particles d,pi;
57        unsigned int ncount=0;
58        findDecayProducts(p,d,pi,ncount);
59        if(ncount!=2 || d.size()!=1 || pi.size()!=1) continue;
60        if(d[0].pid()/p.pid()<0) continue;
61        LorentzTransform boost = LorentzTransform::mkFrameTransformFromBeta(p.momentum().betaVec());
62        Vector3 axis = boost.transform(pi[0].momentum()).p3().unit();
63        double cosL  = axis.dot(p.momentum().p3().unit());
64        // decay angles
65        _h_D2_ctheta->fill(cosL);
66        _h_rate->fill(10);
67      }
68    }
69
70
71    /// Normalise histograms etc., after the run
72    void finalize() {
73      normalize(_h_D2_x);
74      normalize(_h_D2_ctheta);
75      // br of D mode used from PDG2018
76      const double br = 0.0898;
77      scale(_h_rate,br/sumOfWeights()*crossSection()/picobarn);
78    }
79
80    /// @}
81
82
83    /// @name Histograms
84    /// @{
85    BinnedHistoPtr<int> _h_rate;
86    BinnedHistoPtr<string> _h_D2_x;
87    Histo1DPtr _h_D2_ctheta;
88    vector<string> sedges;
89    YODA::Axis<double> axisMap{0.0, 0.4, 0.6, 0.7125, 0.85, 1.0};
90    /// @}
91
92
93  };
94
95
96  RIVET_DECLARE_PLUGIN(ARGUS_1989_I268577);
97
98
99}