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 -> etapi0
 10  class CRYSTAL_BALL_1986_I217547 : public Analysis {
 11  public:
 12
 13    /// Constructor
 14    RIVET_DEFAULT_ANALYSIS_CTOR(CRYSTAL_BALL_1986_I217547);
 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,0.683,2.08)) {
 27        book(_nEtaPi, "TMP/nEta", refData(1, 1, 1));
 28        if (inRange(sqrtS()/GeV,0.9,1.1)) {
 29          book(_h_cTheta,2,1,1);
 30        }
 31        else if (inRange(sqrtS()/GeV,1.1,1.480)) {
 32          book(_h_cTheta,2,1,2);
 33        }
 34      }
 35      else {
 36        throw Error("Invalid CMS energy for CRYSTAL_BALL_1986_I217547");
 37      }
 38    }
 39
 40    void findChildren(const Particle & p,map<long,int> & nRes, int &ncount) {
 41      for (const Particle &child : p.children()) {
 42        if (child.children().empty()) {
 43          nRes[child.pid()]-=1;
 44          --ncount;
 45        }
 46        else {
 47          findChildren(child,nRes,ncount);
 48        }
 49      }
 50    }
 51
 52    /// Perform the per-event analysis
 53    void analyze(const Event& event) {
 54      const FinalState& fs = apply<FinalState>(event, "FS");
 55      // find the final-state particles
 56      map<long,int> nCount;
 57      int ntotal(0);
 58      for (const Particle& p : fs.particles()) {
 59        nCount[p.pid()] += 1;
 60        ++ntotal;
 61      }
 62      const FinalState& ufs = apply<FinalState>(event, "UFS");
 63      for (const Particle& p : ufs.particles(Cuts::pid==PID::ETA)) {
 64        if(p.children().empty()) continue;
 65        map<long,int> nRes=nCount;
 66        int ncount = ntotal;
 67        findChildren(p,nRes,ncount);
 68        if(ncount !=1 ) continue;
 69        bool matched = true;
 70        for(auto const & val : nRes) {
 71          if(val.first==PID::PI0) {
 72            if(val.second!=1) {
 73              matched = false;
 74              break;
 75            }
 76          }
 77          else if(val.second!=0) {
 78            matched = false;
 79            break;
 80          }
 81        }
 82        if(matched) {
 83          double cTheta = abs(p.momentum().z()/p.momentum().p3().mod());
 84          if(cTheta<=0.9)    _nEtaPi->fill(sqrtS()/MeV);
 85          if(_h_cTheta ) _h_cTheta ->fill(cTheta);
 86          break;
 87        }
 88      }
 89    }
 90
 91
 92    /// Normalise histograms etc., after the run
 93    void finalize() {
 94      const double fact = crossSection()/nanobarn/sumOfWeights();
 95      if (_h_cTheta ) scale(_h_cTheta ,fact);
 96      scale(_nEtaPi, fact);
 97      Estimate1DPtr tmp;
 98      book(tmp,1,1,1);
 99      barchart(_nEtaPi,tmp);
100    }
101
102    /// @}
103
104
105    /// @name Histograms
106    /// @{
107    Histo1DPtr _nEtaPi, _h_cTheta;
108    /// @}
109
110
111  };
112
113
114  RIVET_DECLARE_PLUGIN(CRYSTAL_BALL_1986_I217547);
115
116}