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 e+ e- > omega pi0 -> 4pi
 10  class SND_2023_I2693057 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(SND_2023_I2693057);
 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(Cuts::pid==223), "UFS");
 25      book(_nOmegaPi, 1, 1, 1);
 26      for (const string& en : _nOmegaPi.binning().edges<0>()) {
 27        double end = std::stod(en)*MeV;
 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()];
 40          --ncount;
 41        }
 42        else {
 43          findChildren(child,nRes,ncount);
 44        }
 45      }
 46    }
 47
 48    /// Perform the per-event analysis
 49    void analyze(const Event& event) {
 50      // final state particles
 51      const FinalState& fs = apply<FinalState>(event, "FS");
 52
 53      map<long,int> nCount;
 54      int ntotal(0);
 55      for (const Particle& p : fs.particles()) {
 56      	nCount[p.pid()] += 1;
 57      	++ntotal;
 58      }
 59      // require 4 pi final state
 60      if (ntotal!=4) vetoEvent;
 61      if (nCount[211]!=1 || nCount[-211]!=1 || nCount[111]!=2) vetoEvent;
 62      // omega
 63      for (const Particle& p : apply<FinalState>(event, "UFS").particles()) {
 64      	if (p.children().empty()) continue;
 65       	// find the omega
 66      	map<long,int> nRes = nCount;
 67      	int ncount = ntotal;
 68       	findChildren(p,nRes,ncount);
 69      	// omega pi0
 70       	if (ncount!=1) continue;
 71       	bool matched = true;
 72      	for (const auto& val : nRes) {
 73      	  if (val.first==111) {
 74      	    if (val.second !=1) {
 75      	      matched = false;
 76      	      break;
 77      	    }
 78      	  }
 79      	  else if (val.second!=0) {
 80      	    matched = false;
 81      	    break;
 82      	  }
 83      	}
 84      	if (matched) {
 85      	  _nOmegaPi->fill(_ecms);
 86          break;
 87        }
 88      }
 89    }
 90
 91
 92    /// Normalise histograms etc., after the run
 93    void finalize() {
 94      scale(_nOmegaPi, crossSection()/ sumOfWeights() /nanobarn);
 95    }
 96
 97    /// @}
 98
 99
100    /// @name Histograms
101    /// @{
102    BinnedHistoPtr<string> _nOmegaPi;
103    string _ecms;
104    /// @}
105
106
107  };
108
109
110  RIVET_DECLARE_PLUGIN(SND_2023_I2693057);
111
112}