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 -> rho0 rho0
 10  class ARGUS_1991_I296187 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(ARGUS_1991_I296187);
 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::abspid==113), "UFS");
 25      // book histos
 26      if (inRange(sqrtS()/GeV,1.2,2.2)) {
 27        for (unsigned int ix=0; ix<3; ++ix) {
 28          book(_c[ix],"TMP/nMeson_"+toString(ix+1));
 29        }
 30      }
 31      else {
 32        throw Error("Invalid CMS energy for ARGUS_1991_I296187");
 33      }
 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
 48    /// Perform the per-event analysis
 49    void analyze(const Event& event) {
 50      const FinalState& fs = apply<FinalState>(event, "FS");
 51      // find the final-state particles
 52      map<long,int> nCount;
 53      int ntotal(0);
 54      for (const Particle& p : fs.particles()) {
 55        nCount[p.pid()] += 1;
 56        ++ntotal;
 57      }
 58      bool foundRes=false;
 59      // find any rho mesons
 60      Particles rho=apply<UnstableParticles>(event, "UFS").particles();
 61      for (unsigned int ix=0;ix<rho.size();++ix) {
 62       	if (rho[ix].children().empty()) continue;
 63       	map<long,int> nRes=nCount;
 64       	int ncount = ntotal;
 65       	findChildren(rho[ix],nRes,ncount);
 66        bool matched = false;
 67        for (unsigned int iy=ix+1; iy<rho.size(); ++iy) {
 68          if (rho[iy].children().empty()) continue;
 69          map<long,int> nRes2=nRes;
 70          int ncount2 = ncount;
 71          findChildren(rho[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) {
 81            break;
 82          }
 83        }
 84        if (matched) {
 85          _c[2]->fill();
 86          foundRes=true;
 87          break;
 88        }
 89        else {
 90          bool matched2=true;
 91          for (const auto& val : nRes) {
 92            if (abs(val.first)==211) {
 93              if (val.second!=1) {
 94                matched2 = false;
 95                break;
 96              }
 97            }
 98            else if (val.second!=0) {
 99              matched2 = false;
100              break;
101            }
102          }
103          if (matched2) {
104            _c[0]->fill();
105            foundRes=true;
106            break;
107          }
108        }
109      }
110      // 4 pion final-state
111      if (ntotal==4) {
112        if (nCount[PID::PIPLUS]==2 && nCount[PID::PIMINUS]==2) {
113          if (!foundRes) _c[1]->fill();
114        }
115      }
116    }
117
118
119    /// Normalise histograms etc., after the run
120    void finalize() {
121      scale(_c, crossSection()/nanobarn/sumOfWeights());
122      // loop over tables in paper
123      for (unsigned int ix=0; ix < 3; ++ix) {
124        unsigned int iy=2+ix;
125        if (ix==2) iy=5;
126        Estimate1DPtr mult;
127        book(mult, iy, 1, 1);
128        for (auto& b : mult->bins()) {
129          if (inRange(sqrtS(), b.xMin(), b.xMax())) {
130            b.set(_c[ix]->val(), _c[ix]->err());
131          }
132        }
133      }
134    }
135
136    /// @}
137
138
139    /// @name Histograms
140    /// @{
141    CounterPtr _c[3];
142    /// @}
143
144
145  };
146
147
148  RIVET_DECLARE_PLUGIN(ARGUS_1991_I296187);
149
150}