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