1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/Beam.hh"
  4#include "Rivet/Projections/UnstableParticles.hh"
  5
  6namespace Rivet {
  7
  8
  9  /// @brief D_s+ in B or continuum
 10  class ARGUS_1992_I319829 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(ARGUS_1992_I319829);
 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      declare(UnstableParticles(), "UFS");
 24      declare(Beam(), "Beams");
 25      // histos
 26      book(_n_Ds, 1, 1, 1);
 27      for (unsigned int ix=0; ix<2; ++ix) {
 28        book(_h_Ds[ix], 2+ix,1,1);
 29      }
 30    }
 31
 32    /// Recursively walk the decay tree to find decay products of @a p
 33    void findDecayProducts(const Particle& mother, Particles& unstable) {
 34      for(const Particle& p: mother.children()) {
 35        const int id = p.abspid();
 36        if (id == 431) {
 37          unstable.push_back(p);
 38        }
 39        else if (!p.children().empty()) {
 40          findDecayProducts(p, unstable);
 41        }
 42      }
 43    }
 44
 45    /// Perform the per-event analysis
 46    void analyze(const Event& event) {
 47      // Get beams and average beam momentum
 48      const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
 49      const double meanBeamMom = 0.5*(beams.first.p3().mod() + beams.second.p3().mod());
 50      MSG_DEBUG("Avg beam momentum = " << meanBeamMom);
 51      // Find the Upsilon(4S) among the unstables
 52      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
 53      Particles upsilons = ufs.particles(Cuts::pid==300553);
 54      // Continuum
 55      if (upsilons.empty()) {
 56        for (const Particle &p : ufs.particles(Cuts::abspid==431)) {
 57          const double xp = p.p3().mod()/sqrt(sqr(meanBeamMom)-sqr(p.mass()));
 58          _h_Ds[0]->fill(xp);
 59          _n_Ds->fill(10);
 60        }
 61      }
 62      else {
 63        for (const Particle& ups : upsilons) {
 64          Particles unstable;
 65          // Find the decay products we want
 66          findDecayProducts(ups, unstable);
 67          // boost to rest frame (if required)
 68          LorentzTransform cms_boost;
 69          if (ups.p3().mod() > 1*MeV) {
 70            cms_boost = LorentzTransform::mkFrameTransformFromBeta(ups.mom().betaVec());
 71          }
 72          for(const Particle& p : unstable) {
 73            const FourMomentum p2 = cms_boost.transform(p.mom());
 74            double xp = p2.p3().mod()/sqrt(0.25*sqr(ups.mass())-sqr(p.mass()));
 75            _h_Ds[1]->fill(xp);
 76          }
 77        }
 78      }
 79    }
 80
 81
 82    /// Normalise histograms etc., after the run
 83    void finalize() {
 84      // bre d+s -> phi pi+
 85      const double br = 0.045;
 86      normalize(_h_Ds, 1.0, false);
 87      scale(_n_Ds, br/sumOfWeights()*crossSection()/picobarn);
 88    }
 89
 90    /// @}
 91
 92
 93    /// @name Histograms
 94    /// @{
 95    BinnedHistoPtr<int> _n_Ds;
 96    Histo1DPtr _h_Ds[2];
 97    /// @}
 98
 99
100  };
101
102
103  RIVET_DECLARE_PLUGIN(ARGUS_1992_I319829);
104
105}