1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/Beam.hh"
  4#include "Rivet/Projections/FinalState.hh"
  5#include "Rivet/Projections/UnstableParticles.hh"
  6
  7namespace Rivet {
  8
  9
 10  /// @brief kaon production at low energies
 11  class PLUTO_1981_I165122 : public Analysis {
 12  public:
 13
 14    /// Constructor
 15    RIVET_DEFAULT_ANALYSIS_CTOR(PLUTO_1981_I165122);
 16
 17
 18    /// @name Analysis methods
 19    /// @{
 20
 21    /// Book histograms and initialise projections before the run
 22    void init() {
 23
 24      // Initialise and register projections
 25      declare(Beam(), "Beams");
 26      declare(FinalState(), "FS");
 27      declare(UnstableParticles(), "UFS");
 28      // // Book histograms
 29      book(_c_hadrons , "/TMP/sigma_hadrons", refData<YODA::BinnedEstimate<string>>(1, 1, 1));
 30      for(unsigned int ix=0;ix<2;++ix) {
 31        book(_c_muons[ix], "/TMP/sigma_muons_"+toString(ix), refData<YODA::BinnedEstimate<string>>(1+2*ix, 1, 1));
 32        book(_c_kaons[ix], "/TMP/sigma_kaons_"+toString(ix), refData<YODA::BinnedEstimate<string>>(1+2*ix, 1, 1));
 33        for (const string& en : _c_muons[ix].binning().edges<0>()) {
 34          // these points are really Upsilon(1S) decay
 35          if(en == "9.45 - 9.466"s || en=="9.45 - 9.456"s) continue;
 36          const size_t idx = en.find("-");
 37          if(idx!=std::string::npos) {
 38            const double emin = std::stod(en.substr(0,idx));
 39            const double emax = std::stod(en.substr(idx+1,string::npos));
 40            if(inRange(sqrtS()/GeV, emin, emax)) {
 41              _ecms[ix] = en;
 42              break;
 43            }
 44          }
 45          else {
 46            const double end = std::stod(en)*GeV;
 47            if (isCompatibleWithSqrtS(end)) {
 48              _ecms[ix] = en;
 49              break;
 50            }
 51          }
 52        }
 53      }
 54      if (isCompatibleWithSqrtS(9.4*GeV)) {
 55        book(_h_spectrum1, 5, 1, 1);
 56      }
 57      else if (isCompatibleWithSqrtS(30.0*GeV, 1E-2)) {
 58        book(_h_spectrum1, 4, 1, 1);
 59      }
 60      book(_h_spectrum2, 6, 1, 1);
 61      book(_c_hadronsY,"TMP/nUps");
 62    }
 63
 64    /// Recursively walk the decay tree to find decay products of @a p
 65    void findDecayProducts(Particle mother, Particles &kaons, Particles& stable) {
 66      for(const Particle & p: mother.children()) {
 67        const int id = p.pid();
 68        if(id==130 || id ==310) {
 69          kaons.push_back(p);
 70        }
 71        if (id==111 or p.children().empty())
 72          stable.push_back(p);
 73        else
 74          findDecayProducts(p, kaons, stable);
 75      }
 76    }
 77
 78
 79    /// Perform the per-event analysis
 80    void analyze(const Event& event) {
 81      // Get beams and average beam momentum
 82      const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
 83      const double meanBeamMom = ( beams.first.p3().mod() +
 84                                   beams.second.p3().mod() ) / 2.0;
 85      MSG_DEBUG("Avg beam momentum = " << meanBeamMom);
 86      // Find the Upsilons among the unstables
 87      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
 88      Particles upsilons = ufs.particles(Cuts::pid==553);
 89      // Continuum
 90      if (upsilons.empty()) {
 91        MSG_DEBUG("No Upsilons found => continuum event");
 92        // final state particles
 93        const FinalState& fs = apply<FinalState>(event, "FS");
 94        map<long,int> nCount;
 95        int ntotal(0);
 96        for (const Particle& p : fs.particles()) {
 97          nCount[p.pid()] += 1;
 98          ++ntotal;
 99        }
100        // mu+mu- + photons
101        if(nCount[-13]==1 and nCount[13]==1 &&
102           ntotal==2+nCount[22]) {
103          for(unsigned int ix=0;ix<2;++ix) _c_muons[ix]->fill(_ecms[ix]);
104        }
105        // everything else
106        else
107          _c_hadrons->fill(_ecms[0]);
108        // unstable particles
109        for (const Particle& p : ufs.particles(Cuts::pid==130 or Cuts::pid==310)) {
110          if(_h_spectrum1) {
111            double xp = p.p3().mod()/meanBeamMom;
112            _h_spectrum1->fill(xp);
113          }
114          for(unsigned int ix=0;ix<2;++ix) _c_kaons[ix]->fill(_ecms[ix]);
115        }
116      }
117      else {
118        MSG_DEBUG("Upsilons found => resonance event");
119        for (const Particle& ups : upsilons) {
120          Particles kaons,stable;
121          // Find the decay products we want
122          findDecayProducts(ups, kaons, stable);
123          // boost to rest frame (if required)
124          LorentzTransform cms_boost;
125          if (ups.p3().mod() > 1*MeV)
126            cms_boost = LorentzTransform::mkFrameTransformFromBeta(ups.momentum().betaVec());
127          const double mass = ups.mass();
128
129          map<long,int> nCount;
130          int ntotal(0);
131          for (const Particle& p : stable) {
132            nCount[p.pid()] += 1;
133            ++ntotal;
134          }
135          for( const Particle & kaon : kaons) {
136            const FourMomentum p2 = cms_boost.transform(kaon.momentum());
137            const double xp = 2.*p2.p3().mod()/mass;
138            _h_spectrum2->fill(xp);
139            _c_kaons[0]->fill("9.45 - 9.466"s);
140            _c_kaons[1]->fill("9.45 - 9.456"s);
141          }
142          // mu+mu- + photons
143          if(nCount[-13]==1 and nCount[13]==1 &&
144             ntotal==2+nCount[22]) {
145            _c_muons[0]->fill("9.45 - 9.466"s);
146            _c_muons[1]->fill("9.45 - 9.456"s);
147          }
148          // everything else
149          else {
150            _c_hadrons->fill("9.45 - 9.466"s);
151            _c_hadronsY->fill();
152          }
153        }
154      }
155    }
156
157
158    /// Normalise histograms etc., after the run
159    void finalize() {
160      BinnedEstimatePtr<string> ratio;
161      book(ratio,1,1,1);
162      divide(_c_kaons[0],_c_muons[0],ratio);
163      book(ratio,2,1,1);
164      divide(_c_kaons[0],_c_hadrons,ratio);
165      book(ratio,3,1,1);
166      divide(_c_kaons[1],_c_muons[1],ratio);
167      // normalize the spectra if required
168      if(_h_spectrum1) {
169        scale(_h_spectrum1, sqr(sqrtS())*crossSection()/microbarn/sumOfWeights());
170      }
171      if(_h_spectrum2) {
172        scale(_h_spectrum2, 1./ *_c_hadronsY);
173      }
174    }
175
176    /// @}
177
178
179    /// @name Histograms
180    /// @{
181    BinnedHistoPtr<string> _c_hadrons,_c_muons[2],_c_kaons[2];
182    string _ecms[2];
183    Histo1DPtr _h_spectrum1, _h_spectrum2;
184    CounterPtr _c_hadronsY;
185    /// @}
186
187
188  };
189
190
191  RIVET_DECLARE_PLUGIN(PLUTO_1981_I165122);
192
193
194}