1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// @brief ARGUS vector meson production
  9  ///
 10  /// @author Peter Richardson
 11  class ARGUS_1993_I356616 : public Analysis {
 12  public:
 13
 14    RIVET_DEFAULT_ANALYSIS_CTOR(ARGUS_1993_I356616);
 15
 16    void init() {
 17
 18      declare(UnstableParticles(), "UFS");
 19      for (unsigned int ix=0; ix<3; ++ix) {
 20        for (unsigned int iy=0; iy<5; ++iy) {
 21          if (ix==2 && iy==0) continue;
 22          book(_mult[ix][iy], ix+1, 1, iy+1);
 23        }
 24      }
 25
 26      book(_hist_cont_KStarPlus, 4, 1, 1);
 27      book(_hist_Ups1_KStarPlus, 5, 1, 1);
 28      book(_hist_Ups4_KStarPlus, 6, 1, 1);
 29
 30      book(_hist_cont_KStar0, 7, 1, 1);
 31      book(_hist_Ups1_KStar0, 8, 1, 1);
 32      book(_hist_Ups4_KStar0, 9, 1, 1);
 33
 34      book(_hist_cont_Rho0, 10, 1, 1);
 35      book(_hist_Ups1_Rho0, 11, 1, 1);
 36      book(_hist_Ups4_Rho0, 12, 1, 1);
 37
 38      book(_hist_cont_Omega, 13, 1, 1);
 39      book(_hist_Ups1_Omega, 14, 1, 1);
 40
 41
 42      book(_weightSum_cont,"TMP/weightSumcont");
 43      book(_weightSum_Ups1,"TMP/weightSumUps1");
 44      book(_weightSum_Ups4,"TMP/weightSumUps4");
 45    }
 46
 47
 48    void analyze(const Event& e) {
 49      // Find the upsilons
 50      // First in unstable final state
 51      const UnstableParticles& ufs = apply<UnstableParticles>(e, "UFS");
 52      Particles upsilons = ufs.particles(Cuts::pid==553 || Cuts::pid==300553);
 53      // continuum
 54      if (upsilons.empty()) {
 55        _weightSum_cont->fill();
 56        for (const Particle& p : ufs.particles()) {
 57          int id = p.abspid();
 58          double xp = 2.*p.E()/sqrtS();
 59          double beta = p.p3().mod()/p.E();
 60          if (id == 113) {
 61            _hist_cont_Rho0->fill(xp, 1./beta);
 62            _mult[0][1]->fill(Ecm1);
 63          }
 64          else if (id == 313) {
 65            _hist_cont_KStar0->fill(xp, 1./beta);
 66            _mult[0][2]->fill(Ecm1);
 67          }
 68          else if (id == 223) {
 69            _hist_cont_Omega->fill(xp, 1./beta);
 70            _mult[0][0]->fill(Ecm1);
 71          }
 72          else if (id == 323) {
 73            _hist_cont_KStarPlus->fill(xp,1./beta);
 74            _mult[0][3]->fill(Ecm1);
 75          }
 76          else if (id == 333) {
 77            _mult[0][4]->fill(Ecm1);
 78          }
 79        }
 80      }
 81      // found an upsilon
 82      else {
 83        for (const Particle& ups : upsilons) {
 84          const int parentId = ups.pid();
 85          if(parentId == 553)
 86            _weightSum_Ups1->fill();
 87          else
 88            _weightSum_Ups4->fill();
 89          Particles unstable;
 90          // Find the decay products we want
 91          findDecayProducts(ups,unstable);
 92          LorentzTransform cms_boost;
 93          if (ups.p3().mod() > 0.001)
 94            cms_boost = LorentzTransform::mkFrameTransformFromBeta(ups.momentum().betaVec());
 95          double mass = ups.mass();
 96          for( const Particle & p : unstable) {
 97            int id = p.abspid();
 98            FourMomentum p2 = cms_boost.transform(p.momentum());
 99            double xp = 2.*p2.E()/mass;
100            double beta = p2.p3().mod()/p2.E();
101            if (id == 113) {
102              if (parentId == 553) {
103                _hist_Ups1_Rho0->fill(xp,1./beta);
104                _mult[1][1]->fill(Ecm2);
105              }
106              else {
107                _hist_Ups4_Rho0->fill(xp,1./beta);
108                _mult[2][1]->fill(Ecm3);
109              }
110            }
111            else if (id == 313) {
112              if (parentId == 553) {
113                _hist_Ups1_KStar0->fill(xp,1./beta);
114                _mult[1][2]->fill(Ecm2);
115              }
116              else {
117                _hist_Ups4_KStar0->fill(xp,1./beta);
118                _mult[2][2]->fill(Ecm3);
119              }
120            }
121            else if (id == 223) {
122              if (parentId == 553) {
123                _hist_Ups1_Omega->fill(xp,1./beta);
124                _mult[1][0]->fill(Ecm2);
125              }
126            }
127            else if (id == 323) {
128              if (parentId == 553) {
129                _hist_Ups1_KStarPlus->fill(xp,1./beta);
130                _mult[1][3]->fill(Ecm2);
131              }
132              else {
133                _hist_Ups4_KStarPlus->fill(xp,1./beta);
134                _mult[2][3]->fill(Ecm3);
135              }
136            }
137            else if (id == 333) {
138              if (parentId == 553) {
139                _mult[1][4]->fill(Ecm2);
140              }
141              else {
142                _mult[2][4]->fill(Ecm3);
143              }
144            }
145          }
146        }
147      }
148    }
149
150
151    void finalize() {
152      // multiplicities
153      const vector<CounterPtr> scales = {_weightSum_cont,_weightSum_Ups1,_weightSum_Ups4};
154      for (size_t i=0; i<scales.size(); ++i) {
155        if (!scales[i]->val())  continue;
156        for (size_t j=0; j<5; ++j) {
157          if (i==2 && j==0) continue;
158          _mult[i][j]->scaleW(1.0/ scales[i]->val());
159        }
160      }
161
162      // spectra
163      if (_weightSum_cont->val() > 0.) {
164        scale(_hist_cont_KStarPlus, 1. / *_weightSum_cont);
165        scale(_hist_cont_KStar0   , 1. / *_weightSum_cont);
166        scale(_hist_cont_Rho0     , 1. / *_weightSum_cont);
167        scale(_hist_cont_Omega    , 1. / *_weightSum_cont);
168      }
169      if (_weightSum_Ups1->val() > 0.) {
170        scale(_hist_Ups1_KStarPlus, 1. / *_weightSum_Ups1);
171        scale(_hist_Ups1_KStar0   , 1. / *_weightSum_Ups1);
172        scale(_hist_Ups1_Rho0     , 1. / *_weightSum_Ups1);
173        scale(_hist_Ups1_Omega    , 1. / *_weightSum_Ups1);
174      }
175      if (_weightSum_Ups4->val() > 0.) {
176        scale(_hist_Ups4_KStarPlus, 1. / *_weightSum_Ups4);
177        scale(_hist_Ups4_KStar0   , 1. / *_weightSum_Ups4);
178        scale(_hist_Ups4_Rho0     , 1. / *_weightSum_Ups4);
179      }
180    }
181
182
183  private:
184
185    Histo1DPtr _hist_cont_KStarPlus, _hist_Ups1_KStarPlus, _hist_Ups4_KStarPlus;
186    Histo1DPtr _hist_cont_KStar0, _hist_Ups1_KStar0, _hist_Ups4_KStar0   ;
187    Histo1DPtr _hist_cont_Rho0, _hist_Ups1_Rho0,  _hist_Ups4_Rho0;
188    Histo1DPtr _hist_cont_Omega, _hist_Ups1_Omega;
189
190    BinnedHistoPtr<string> _mult[3][5];
191    CounterPtr _weightSum_cont,_weightSum_Ups1,_weightSum_Ups4;
192    const string Ecm1 = "10.45", Ecm2 = "9.46", Ecm3 = "10.58";
193
194
195    /// Recursively walk the decay tree to find decay products of @a p
196    void findDecayProducts(Particle mother, Particles& unstable) {
197      for(const Particle & p: mother.children()) {
198        const int id = abs(p.pid());
199        if (id == 113 || id == 313 || id == 323 ||
200            id == 333 || id == 223 ) {
201          unstable.push_back(p);
202        }
203        else if(!p.children().empty())
204          findDecayProducts(p, unstable);
205      }
206    }
207
208  };
209
210
211
212  RIVET_DECLARE_ALIASED_PLUGIN(ARGUS_1993_I356616, ARGUS_1993_S2789213);
213
214}