1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// @brief ARGUS pi+/-, K+/- and proton/antiproton spectrum at Upsilon(4S)
  9  ///
 10  /// @author Peter Richardson
 11  class ARGUS_1993_I340894 : public Analysis {
 12  public:
 13
 14    RIVET_DEFAULT_ANALYSIS_CTOR(ARGUS_1993_I340894);
 15
 16
 17    void analyze(const Event& e) {
 18      // Find the upsilons
 19      const UnstableParticles& ufs = apply<UnstableParticles>(e, "UFS");
 20      for (const Particle& p : ufs.particles(Cuts::pid==300553)) {
 21        _weightSum->fill();
 22        Particles pionsA,pionsB,protonsA,protonsB,kaons;
 23        // Find the decay products we want
 24        findDecayProducts(p, pionsA, pionsB, protonsA, protonsB, kaons);
 25        LorentzTransform cms_boost;
 26        if (p.p3().mod() > 1*MeV)
 27          cms_boost = LorentzTransform::mkFrameTransformFromBeta(p.momentum().betaVec());
 28        for (size_t ix = 0; ix < pionsA.size(); ++ix) {
 29          FourMomentum ptemp(pionsA[ix].momentum());
 30          FourMomentum p2 = cms_boost.transform(ptemp);
 31          double pcm = cms_boost.transform(ptemp).vector3().mod();
 32          _histPiA->fill(pcm);
 33        }
 34        _multPiA->fill(Ecm,double(pionsA.size()));
 35        for (size_t ix = 0; ix < pionsB.size(); ++ix) {
 36          double pcm = cms_boost.transform(pionsB[ix].momentum()).vector3().mod();
 37          _histPiB->fill(pcm);
 38        }
 39        _multPiB->fill(Ecm,double(pionsB.size()));
 40        for (size_t ix = 0; ix < protonsA.size(); ++ix) {
 41          double pcm = cms_boost.transform(protonsA[ix].momentum()).vector3().mod();
 42          _histpA->fill(pcm);
 43        }
 44        _multpA->fill(Ecm,double(protonsA.size()));
 45        for (size_t ix = 0; ix < protonsB.size(); ++ix) {
 46          double pcm = cms_boost.transform(protonsB[ix].momentum()).vector3().mod();
 47          _histpB->fill(pcm);
 48        }
 49        _multpB->fill(Ecm,double(protonsB.size()));
 50        for (size_t ix = 0 ;ix < kaons.size(); ++ix) {
 51          double pcm = cms_boost.transform(kaons[ix].momentum()).vector3().mod();
 52          _histKA->fill(pcm);
 53          _histKB->fill(pcm);
 54        }
 55        _multK->fill(Ecm,double(kaons.size()));
 56      }
 57    }
 58
 59
 60    void finalize() {
 61      if (_weightSum->val() > 0.) {
 62        scale(_histPiA, 1. / *_weightSum);
 63        scale(_histPiB, 1. / *_weightSum);
 64        scale(_histKA , 1. / *_weightSum);
 65        scale(_histKB , 1. / *_weightSum);
 66        scale(_histpA , 1. / *_weightSum);
 67        scale(_histpB , 1. / *_weightSum);
 68        scale(_multPiA, 1. / *_weightSum);
 69        scale(_multPiB, 1. / *_weightSum);
 70        scale(_multK  , 1. / *_weightSum);
 71        scale(_multpA , 1. / *_weightSum);
 72        scale(_multpB , 1. / *_weightSum);
 73      }
 74    }
 75
 76
 77    void init() {
 78      declare(UnstableParticles(), "UFS");
 79
 80      // spectra
 81      book(_histPiA ,1, 1, 1);
 82      book(_histPiB ,2, 1, 1);
 83      book(_histKA  ,3, 1, 1);
 84      book(_histKB  ,6, 1, 1);
 85      book(_histpA  ,4, 1, 1);
 86      book(_histpB  ,5, 1, 1);
 87      // multiplicities
 88      book(_multPiA , 7, 1, 1);
 89      book(_multPiB , 8, 1, 1);
 90      book(_multK   , 9, 1, 1);
 91      book(_multpA  ,10, 1, 1);
 92      book(_multpB  ,11, 1, 1);
 93
 94      book(_weightSum, "TMP/weightSum");
 95    }
 96
 97
 98  private:
 99
100    /// @{
101    /// Count of weights
102    CounterPtr _weightSum;
103    /// Spectra
104    Histo1DPtr _histPiA, _histPiB, _histKA, _histKB, _histpA, _histpB;
105    /// Multiplicities
106    BinnedHistoPtr<string> _multPiA, _multPiB, _multK, _multpA, _multpB;
107    const string Ecm = "10.58";
108    /// @}
109
110
111    void findDecayProducts(Particle parent, Particles & pionsA, Particles & pionsB,
112                           Particles & protonsA, Particles & protonsB, Particles & kaons) {
113      int parentId = parent.pid();
114      for (const Particle & p : parent.children()) {
115        int id = abs(p.pid());
116        if (id == PID::PIPLUS) {
117          if (parentId != PID::LAMBDA && parentId != PID::K0S) {
118            pionsA.push_back(p);
119            pionsB.push_back(p);
120          }
121          else
122            pionsB.push_back(p);
123        }
124        else if (id == PID::PROTON) {
125          if (parentId != PID::LAMBDA && parentId != PID::K0S) {
126            protonsA.push_back(p);
127            protonsB.push_back(p);
128          }
129          else
130            protonsB.push_back(p);
131        }
132        else if (id == PID::KPLUS) {
133          kaons.push_back(p);
134        }
135        else if (!p.children().empty())
136          findDecayProducts(p, pionsA, pionsB, protonsA, protonsB, kaons);
137      }
138    }
139
140
141  };
142
143
144
145  RIVET_DECLARE_ALIASED_PLUGIN(ARGUS_1993_I340894, ARGUS_1993_S2653028);
146
147}