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      for(auto & b : _h_D2_x->bins()) {
 75        b.scaleW(1./axisMap.width(b.index()));
 76      }
 77      normalize(_h_D2_ctheta);
 78      // br of D mode used from PDG2018
 79      const double br = 0.0898;
 80      scale(_h_rate,br/sumOfWeights()*crossSection()/picobarn);
 81    }
 82
 83    /// @}
 84
 85
 86    /// @name Histograms
 87    /// @{
 88    BinnedHistoPtr<int> _h_rate;
 89    BinnedHistoPtr<string> _h_D2_x;
 90    Histo1DPtr _h_D2_ctheta;
 91    vector<string> sedges;
 92    YODA::Axis<double> axisMap{0.0, 0.5, 0.6, 0.7, 0.85, 1.0};
 93    /// @}
 94
 95
 96  };
 97
 98
 99  RIVET_DECLARE_PLUGIN(ARGUS_1989_I268577);
100
101
102}