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