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 Add a short analysis description here
 10  class DM2_1990_I297706 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(DM2_1990_I297706);
 15
 16
 17    /// @name Analysis methods
 18    /// @{
 19
 20    /// Book histograms and initialise projections before the run
 21    void init() {
 22      // Initialise and register projections
 23      declare(FinalState(), "FS");
 24      declare(UnstableParticles(), "UFS");
 25      book(_nProton, 1, 1, 1);
 26      book(_nLambda, 2, 1, 1);
 27    }
 28
 29    void findChildren(const Particle & p,map<long,int> & nRes, int &ncount) {
 30      for (const Particle &child : p.children()) {
 31        if(child.children().empty()) {
 32          nRes[child.pid()]-=1;
 33          --ncount;
 34        }
 35        else {
 36          findChildren(child,nRes,ncount);
 37        }
 38      }
 39    }
 40
 41    /// Perform the per-event analysis
 42    void analyze(const Event& event) {
 43      const FinalState& fs = apply<FinalState>(event, "FS");
 44      // total hadronic and muonic cross sections
 45      map<long,int> nCount;
 46      int ntotal(0);
 47      for (const Particle& p : fs.particles()) {
 48        nCount[p.pid()] += 1;
 49        ++ntotal;
 50      }
 51      if (ntotal==2 && nCount[2212]==1 && nCount[-2212]==1) {
 52        _nProton->fill(Ecm);
 53      }
 54
 55      // find the Lambdas
 56      const FinalState& ufs = apply<UnstableParticles>(event, "UFS");
 57      for (unsigned int ix=0; ix<ufs.particles().size(); ++ix) {
 58        const Particle& p1 = ufs.particles()[ix];
 59        if(abs(p1.pid())!=3122) continue;
 60        bool matched = false;
 61        // check fs
 62        bool fs = true;
 63        for (const Particle& child : p1.children()) {
 64          if(child.pid()==p1.pid()) {
 65            fs = false;
 66            break;
 67          }
 68        }
 69        if(!fs) continue;
 70        // find the children
 71        map<long,int> nRes = nCount;
 72        int ncount = ntotal;
 73        findChildren(p1,nRes,ncount);
 74        for (unsigned int iy=ix+1;iy<ufs.particles().size();++iy) {
 75          const Particle& p2 = ufs.particles()[iy];
 76          if (abs(p2.pid())!=3122) continue;
 77          // check fs
 78          bool fs = true;
 79          for (const Particle& child : p2.children()) {
 80            if(child.pid()==p2.pid()) {
 81              fs = false;
 82              break;
 83            }
 84          }
 85          if(!fs) continue;
 86          map<long,int> nRes2 = nRes;
 87          int ncount2 = ncount;
 88          findChildren(p2,nRes2,ncount2);
 89          if (ncount2!=0) continue;
 90          matched=true;
 91          for (const auto& val : nRes2) {
 92            if(val.second!=0) {
 93              matched = false;
 94              break;
 95            }
 96          }
 97          if (matched) {
 98            _nLambda->fill(Ecm);
 99            break;
100          }
101        }
102        if (matched)  break;
103      }
104    }
105
106
107    /// Normalise histograms etc., after the run
108    void finalize() {
109      const double fact = crossSection()/ sumOfWeights() /picobarn;
110      scale({_nProton, _nLambda}, fact);
111    }
112    /// @}
113
114    /// @name Histograms
115    /// @{
116    BinnedHistoPtr<string> _nProton,_nLambda;
117    const string Ecm = "2.386";
118    /// @}
119
120  };
121
122
123  RIVET_DECLARE_PLUGIN(DM2_1990_I297706);
124
125
126}