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