1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/UnstableParticles.hh"
  5
  6namespace Rivet {
  7
  8
  9  /// @brief charged multiplicity at 4s and nearby continuum
 10  class ARGUS_1992_I319102 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(ARGUS_1992_I319102);
 15
 16
 17    /// @name Analysis methods
 18    /// @{
 19
 20    /// Book histograms and initialise projections before the run
 21    void init() {
 22
 23      // Initialise and register projections
 24      declare(UnstableParticles(), "UFS");
 25      declare(FinalState(), "FS");
 26
 27      // Book histograms
 28      if (isCompatibleWithSqrtS(10.47*GeV)) {
 29        book(_h_N, 2, 1, 1);
 30        book(_h_tot_N,4,1,1);
 31      }
 32      book(_h_N_Upsilon, 3, 1, 1);
 33      book(_h_N_tot_Upsilon,5,1,1);
 34      // counters for R
 35      book(_c_hadrons, "/TMP/sigma_hadrons");
 36      book(_c_muons, "/TMP/sigma_muons");
 37      book(_w_cont,"/TMP/w_cont");
 38      book(_w_ups ,"/TMP/w_ups" );
 39    }
 40
 41    /// Recursively walk the decay tree to find decay products of @a p
 42    void findDecayProducts(Particle mother, unsigned int& nCharged) {
 43      for (const Particle & p: mother.children()) {
 44        if (!p.children().empty()) {
 45          findDecayProducts(p, nCharged);
 46        }
 47        else if (PID::isCharged(p.pid())) {
 48          ++nCharged;
 49        }
 50      }
 51    }
 52
 53
 54    /// Perform the per-event analysis
 55    void analyze(const Event& event) {
 56      const FinalState& fs = apply<FinalState>(event, "FS");
 57      // Find the Upsilons among the unstables
 58      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
 59      Particles upsilons = ufs.particles(Cuts::pid==300553);
 60      // Continuum
 61      if (upsilons.empty()) {
 62        map<long,int> nCount;
 63        int ntotal(0);
 64        unsigned int nCharged(0);
 65        for (const Particle& p : fs.particles()) {
 66          nCount[p.pid()] += 1;
 67          ++ntotal;
 68          if (PID::isCharged(p.pid())) ++nCharged;
 69        }
 70        // mu+mu- + photons
 71        if (nCount[-13]==1 and nCount[13]==1 &&
 72           ntotal==2+nCount[22])
 73          _c_muons->fill();
 74        // everything else
 75        else {
 76          _c_hadrons->fill();
 77          if (_h_N) {
 78            _h_N->fill(nCharged);
 79            _h_tot_N->fill(Ecm1, nCharged);
 80            _w_cont->fill();
 81          }
 82        }
 83      }
 84      // upsilon 4s
 85      else {
 86        for (const Particle& ups : upsilons) {
 87          unsigned int nCharged(0);
 88          findDecayProducts(ups,nCharged);
 89          _h_N_Upsilon->fill(nCharged);
 90          _h_N_tot_Upsilon->fill(Ecm2, nCharged);
 91          _w_ups->fill();
 92        }
 93      }
 94    }
 95
 96
 97    /// Normalise histograms etc., after the run
 98    void finalize() {
 99      Estimate0D R = *_c_hadrons/ *_c_muons;
100      BinnedEstimatePtr<string> mult;
101      book(mult, 1, 1, 1);
102      if (isCompatibleWithSqrtS(9.36*GeV)) {
103        mult->bin(1).set(R.val(), R.errPos());
104      }
105      if (_h_N) {
106        normalize(_h_N, 100.);
107      if (_w_cont->val()!=0)
108        scale(_h_tot_N,1./ *_w_cont);
109      }
110      normalize(_h_N_Upsilon, 100.);
111      if (_w_ups->val()!=0) {
112        scale(_h_N_tot_Upsilon,1./ *_w_ups);
113      }
114    }
115
116    /// @}
117
118
119    /// @name Histograms
120    /// @{
121    BinnedHistoPtr<int> _h_N, _h_N_Upsilon;
122    BinnedHistoPtr<string> _h_tot_N, _h_N_tot_Upsilon;
123    CounterPtr _c_hadrons, _c_muons;
124    CounterPtr _w_cont,_w_ups;
125    const string Ecm1 = "10.47", Ecm2 = "10.575";
126    /// @}
127
128
129  };
130
131
132  RIVET_DECLARE_PLUGIN(ARGUS_1992_I319102);
133
134
135}