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 gamma gamma -> omega omega
 10  class ARGUS_1996_I403304 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(ARGUS_1996_I403304);
 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 histos
 26      if (inRange(sqrtS()/GeV, 1.6, 3.0)) {
 27        for (unsigned int ix=0;ix<3;++ix)
 28          book(_nMeson[ix],"TMP/nMeson_"+toString(ix+1));
 29      }
 30      else {
 31	     throw Error("Invalid CMS energy for ARGUS_1996_I403304");
 32      }
 33    }
 34
 35    void findChildren(const Particle & p,map<long,int> & nRes, int &ncount) {
 36      for (const Particle &child : p.children()) {
 37        if (child.children().empty()) {
 38          nRes[child.pid()]-=1;
 39          --ncount;
 40        }
 41        else {
 42	       findChildren(child,nRes,ncount);
 43        }
 44      }
 45    }
 46
 47    /// Perform the per-event analysis
 48    void analyze(const Event& event) {
 49      const FinalState& fs = apply<FinalState>(event, "FS");
 50      // find the final-state particles
 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 omega mesons
 58      Particles omega=apply<UnstableParticles>(event, "UFS").particles(Cuts::abspid==223);
 59      bool nonRes=true;
 60      for (unsigned int ix=0; ix<omega.size(); ++ix) {
 61       	if(omega[ix].children().empty()) continue;
 62       	map<long,int> nRes=nCount;
 63       	int ncount = ntotal;
 64       	findChildren(omega[ix],nRes,ncount);
 65        bool matched=false;
 66        // omega omega
 67        for (unsigned int iy=ix+1; iy<omega.size(); ++iy) {
 68          if (omega[iy].children().empty()) continue;
 69          map<long,int> nRes2=nRes;
 70          int ncount2 = ncount;
 71          findChildren(omega[iy],nRes2,ncount2);
 72          if (ncount2 !=0 ) continue;
 73          matched = true;
 74          for (const auto& val : nRes2) {
 75            if (val.second!=0) {
 76              matched = false;
 77              break;
 78            }
 79          }
 80          if (matched)  break;
 81        }
 82        if (matched) {
 83          _nMeson[0]->fill();
 84          nonRes=false;
 85          break;
 86        }
 87        matched=true;
 88        for (const auto& val : nRes) {
 89          if (abs(val.first)==211 || val.first==111) {
 90            if(val.second!=1) {
 91              matched = false;
 92              break;
 93            }
 94          }
 95          else {
 96            if (val.second!=0) {
 97              matched = false;
 98              break;
 99            }
100          }
101        }
102        if (matched) {
103          nonRes=false;
104          _nMeson[1]->fill();
105          break;
106        }
107      }
108      if (nonRes && ntotal==6 && nCount[PID::PI0]==2 &&
109         nCount[PID::PIPLUS]==2 && nCount[PID::PIMINUS]==2 ) {
110        _nMeson[2]->fill();
111      }
112    }
113
114
115    /// Normalise histograms etc., after the run
116    void finalize() {
117      scale(_nMeson, crossSection()/nanobarn/sumOfWeights());
118      // loop over tables in paper
119      for (unsigned int ix=0; ix<3; ++ix) {
120        Estimate1DPtr mult;
121        book(mult, 1, 1, 1+ix);
122        for (auto& b : mult->bins()) {
123          if (inRange(sqrtS(), b.xMin(), b.xMax())) {
124            b.setVal(_nMeson[ix]->val());
125            b.setErr(_nMeson[ix]->err());
126          }
127        }
128      }
129    }
130
131    /// @}
132
133
134    /// @name Histograms
135    /// @{
136    CounterPtr _nMeson[3];
137    /// @}
138
139
140  };
141
142
143  RIVET_DECLARE_PLUGIN(ARGUS_1996_I403304);
144
145}