1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// @brief Production of the $\eta'(958)$ and $f_0(980)$ in $e^+e^-$ annihilation in the Upsilon region
  9  ///
 10  /// @author Peter Richardson
 11  class ARGUS_1993_I342061 : public Analysis {
 12  public:
 13
 14    RIVET_DEFAULT_ANALYSIS_CTOR(ARGUS_1993_I342061);
 15
 16
 17    void init() {
 18      declare(UnstableParticles(), "UFS");
 19
 20      book(_weightSum_cont, "TMP/weightSum_cont");
 21      book(_weightSum_Ups1, "TMP/weightSum_Ups1");
 22      book(_weightSum_Ups2, "TMP/weightSum_Ups2");
 23
 24      book(_count_etaPrime_highZ, 1,1,1);
 25      book(_count_etaPrime_allZ , 1,1,2);
 26      book(_count_f0            , 5,1,1);
 27
 28      book(_hist_cont_f0 ,2, 1, 1);
 29      book(_hist_Ups1_f0 ,3, 1, 1);
 30      book(_hist_Ups2_f0 ,4, 1, 1);
 31    }
 32
 33
 34    void analyze(const Event& e) {
 35      // Find the Upsilons among the unstables
 36      const UnstableParticles& ufs = apply<UnstableParticles>(e, "UFS");
 37      Particles upsilons = ufs.particles(Cuts::pid==553 or Cuts::pid==100553);
 38      // Continuum
 39      if (upsilons.empty()) {
 40        MSG_DEBUG("No Upsilons found => continuum event");
 41        _weightSum_cont->fill();
 42        for (const Particle& p : ufs.particles()) {
 43          const int id = p.pid();
 44          const double xp = 2.*p.E()/sqrtS();
 45          const double beta = p.p3().mod() / p.E();
 46          if (id == 9010221) {
 47            _hist_cont_f0->fill(xp, 1./beta);
 48            _count_f0->fill(string("10.45"));
 49          } else if (id == 331) {
 50            if (xp > 0.35) _count_etaPrime_highZ->fill(string("9.36 - 10.45"));
 51            _count_etaPrime_allZ->fill(string("9.36 - 10.45"));
 52          }
 53        }
 54      }
 55      // Upsilon(s) found
 56      else {
 57        for (const Particle& ups : upsilons) {
 58          const int parentId = ups.pid();
 59          if(parentId==553) {
 60            _weightSum_Ups1->fill();
 61          }
 62          else {
 63            _weightSum_Ups2->fill();
 64          }
 65          Particles unstable;
 66          // Find the decay products we want
 67          findDecayProducts(ups, unstable);
 68          // boost to rest frame (if required)
 69          LorentzTransform cms_boost;
 70          if (ups.p3().mod() > 1*MeV)
 71            cms_boost = LorentzTransform::mkFrameTransformFromBeta(ups.momentum().betaVec());
 72          const double mass = ups.mass();
 73          // loop over decay products
 74          for(const Particle& p : unstable) {
 75            const int id = p.pid();
 76            const FourMomentum p2 = cms_boost.transform(p.momentum());
 77            const double xp = 2.*p2.E()/mass;
 78            const double beta = p2.p3().mod()/p2.E();
 79            if (id == 9010221) {
 80              if(parentId == 553 ) {
 81        	_hist_Ups1_f0->fill(xp, 1./beta);
 82        	_count_f0->fill(string("9.46"));
 83              }
 84              else {
 85        	_hist_Ups2_f0->fill(xp, 1./beta);
 86        	_count_f0->fill(string("10.02"));
 87              }
 88            }
 89            else if ( id == 331 ) {
 90              if (parentId == 553) {
 91        	if (xp > 0.35) _count_etaPrime_highZ->fill(string("9.46"));
 92        	_count_etaPrime_allZ->fill(string("9.46"));
 93              }
 94              else {
 95        	_count_etaPrime_allZ->fill(string("10.02"));
 96              }
 97            }
 98          }
 99        }
100      }
101    }
102
103    void finalize() {
104      // Point at 9.905
105      if (_weightSum_cont->val() > 0) { 
106        _count_etaPrime_highZ->binAt(string("9.36 - 10.45")).scaleW(1./_weightSum_cont->val());
107        _count_etaPrime_allZ ->binAt(string("9.36 - 10.45")).scaleW(1./_weightSum_cont->val());
108        _count_f0            ->binAt(       string("10.45")).scaleW(1./_weightSum_cont->val());
109      }
110      // Point at 9.460
111      if (_weightSum_Ups1->val() > 0) {  
112        _count_etaPrime_highZ->binAt(string("9.46")).scaleW(1./_weightSum_Ups1->val());
113        _count_etaPrime_allZ ->binAt(string("9.46")).scaleW(1./_weightSum_Ups1->val());
114        _count_f0            ->binAt(string("9.46")).scaleW(1./_weightSum_Ups1->val());
115      }
116      // Point at 10.02
117      if (_weightSum_Ups2->val() > 0) {
118        _count_etaPrime_allZ->binAt(string("10.02")).scaleW(1./_weightSum_Ups2->val());
119        _count_f0           ->binAt(string("10.02")).scaleW(1./_weightSum_Ups2->val());
120      }
121      // Scale histos
122      if (_weightSum_cont->val() > 0.) scale(_hist_cont_f0, 1./ *_weightSum_cont);
123      if (_weightSum_Ups1->val() > 0.) scale(_hist_Ups1_f0, 1./ *_weightSum_Ups1);
124      if (_weightSum_Ups2->val() > 0.) scale(_hist_Ups2_f0, 1./ *_weightSum_Ups2);
125    }
126
127
128  private:
129
130    /// @name Counters
131    /// @{
132    BinnedHistoPtr<string> _count_etaPrime_highZ, _count_etaPrime_allZ, _count_f0;
133    CounterPtr _weightSum_cont,_weightSum_Ups1,_weightSum_Ups2;
134    /// @}
135
136
137    /// Histos
138    Histo1DPtr _hist_cont_f0, _hist_Ups1_f0, _hist_Ups2_f0;
139
140
141    /// Recursively walk the decay tree to find decay products of @a p
142    void findDecayProducts(Particle mother, Particles& unstable) {
143      for (const Particle & p: mother.children()) {
144        const int id = p.pid();
145        if (id == 331 || id == 9010221) {
146          unstable.push_back(p);
147        }
148        else if(!p.children().empty()) {
149          findDecayProducts(p, unstable);
150        }
151      }
152    }
153
154  };
155
156
157
158  RIVET_DECLARE_ALIASED_PLUGIN(ARGUS_1993_I342061, ARGUS_1993_S2669951);
159
160}