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 SND_2000_I527752 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(SND_2000_I527752);
 15
 16
 17    /// @name Analysis methods
 18    /// @{
 19
 20    /// Book histograms and initialise projections before the run
 21    void init() {
 22
 23      // Initialise and register projections
 24      declare(FinalState(), "FS");
 25      declare(UnstableParticles(), "UFS");
 26      book(_numOmegaPi, 1, 1, 1);
 27      for (const string& en : _numOmegaPi.binning().edges<0>()) {
 28        const size_t idx = en.find("-");
 29        if(idx!=std::string::npos) {
 30          const double emin = std::stod(en.substr(0,idx));
 31          const double emax = std::stod(en.substr(idx+1,string::npos));
 32          if(inRange(sqrtS()/MeV, emin, emax)) {
 33            _ecms = en;
 34            break;
 35          }
 36        }
 37        else {
 38          const double end = std::stod(en)*MeV;
 39          if (isCompatibleWithSqrtS(end)) {
 40            _ecms = en;
 41            break;
 42          }
 43        }
 44      }
 45      if (_ecms.empty())
 46        MSG_ERROR("Beam energy incompatible with analysis.");
 47    }
 48
 49    void findChildren(const Particle & p,map<long,int> & nRes, int &ncount) {
 50      for (const Particle &child : p.children()) {
 51	if(child.children().empty()) {
 52	  nRes[child.pid()]+=1;
 53	  ++ncount;
 54	}
 55	else
 56	  findChildren(child,nRes,ncount);
 57      }
 58    }
 59
 60
 61    /// Perform the per-event analysis
 62    void analyze(const Event& event) {
 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      // three particles (pi0 pi0 gamma)
 72      if(ntotal!=3) vetoEvent;
 73      const FinalState& ufs = apply<FinalState>(event, "UFS");
 74      for (const Particle& p : ufs.particles()) {
 75	if(p.children().empty()) continue;
 76	// find the omega
 77	if(p.pid()==223) {
 78	  map<long,int> nRes;
 79	  int ncount(0);
 80	  findChildren(p,nRes,ncount);
 81	  // only omega to pi0 gamma mode
 82	  if(ncount!=2) continue;
 83	  if(nRes[111]!=1 || nRes[22]!=1) continue;
 84	  // omega pi0
 85	  if(nCount[111]-nRes[111]==1)
 86	    _numOmegaPi->fill(_ecms);
 87	}
 88      }
 89    }
 90
 91
 92    /// Normalise histograms etc., after the run
 93    void finalize() {
 94      scale(_numOmegaPi, crossSection()/ sumOfWeights() /nanobarn);
 95    }
 96
 97    /// @}
 98
 99    /// @name Histograms
100    /// @{
101    BinnedHistoPtr<string> _numOmegaPi;
102    string _ecms;
103    /// @}
104
105
106  };
107
108
109  RIVET_DECLARE_PLUGIN(SND_2000_I527752);
110
111
112}