1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/UnstableParticles.hh"
  5#include "Rivet/Projections/Beam.hh"
  6
  7namespace Rivet {
  8
  9
 10  /// @brief pi0 and gamma spectra at 14, 22 and 34
 11  class CELLO_1983_I191415 : public Analysis {
 12  public:
 13
 14    /// Constructor
 15    RIVET_DEFAULT_ANALYSIS_CTOR(CELLO_1983_I191415);
 16
 17
 18    /// @name Analysis methods
 19    /// @{
 20
 21    /// Book histograms and initialise projections before the run
 22    void init() {
 23
 24      // Initialise and register projections
 25      declare(Beam(), "Beams");
 26      declare(FinalState(), "FS");
 27      declare(UnstableParticles(), "UFS");
 28
 29      unsigned int iloc(0);
 30      if (isCompatibleWithSqrtS(14*GeV)) {
 31        iloc=1;
 32        _axes["gamma"] = YODA::Axis<double>({0.0428, 0.057, 0.07125, 0.0855, 0.09975, 0.114, 0.12825,
 33                                             0.14255, 0.17105, 0.19955, 0.22805, 0.25655, 0.2855,
 34                                             0.35625, 0.42755, 0.4988, 0.57, 0.8552});
 35        _axes["pi0"] = YODA::Axis<double>({0.1325, 0.1635, 0.2505, 0.3375, 0.434, 0.54, 0.791, 1.04});
 36      }
 37      else if (isCompatibleWithSqrtS(22*GeV)) {
 38        iloc=2;
 39        _axes["gamma"] = YODA::Axis<double>({0.02725, 0.03635, 0.04545, 0.05455, 0.06365, 0.07275,
 40                                             0.08185, 0.09095, 0.1091, 0.1273, 0.1455, 0.1637, 0.1819,
 41                                             0.22735, 0.27285, 0.31835, 0.36385, 0.54575, 0.72765});
 42        _axes["pi0"] = YODA::Axis<double>({0.0825, 0.1055, 0.1605, 0.2155, 0.2705, 0.3555, 0.5145, 0.733, 0.961});
 43      }
 44      else if (isCompatibleWithSqrtS(34*GeV)) {
 45        iloc=3;
 46        _axes["gamma"] = YODA::Axis<double>({0.01765, 0.02355, 0.02945, 0.0353, 0.04115, 0.04705, 0.05295,
 47                                             0.05885, 0.0706, 0.08235, 0.0941, 0.106, 0.1177, 0.1471, 0.1765,
 48                                             0.2059, 0.2353, 0.35295, 0.47065});
 49        _axes["pi0"] = YODA::Axis<double>({0.054, 0.068, 0.104, 0.14, 0.176, 0.23, 0.338, 0.45, 0.71, 0.97});
 50      }
 51      else {
 52        MSG_ERROR("Beam energy not supported!");
 53      }
 54
 55
 56      // Book histograms
 57      book(_h["gamma"], iloc,   1, 1);
 58      book(_h["pi0"],   iloc+3, 1, 1);
 59    }
 60
 61
 62    /// Perform the per-event analysis
 63    void analyze(const Event& event) {
 64
 65      if (_edges.empty()) {
 66        for (const auto& item : _h) {
 67          _edges[item.first] = item.second->xEdges();
 68        }
 69      }
 70
 71      // at least 5 charged FS particles
 72      const FinalState& fs = apply<FinalState>(event, "FS");
 73      const size_t numParticles = fs.particles().size();
 74
 75      if (numParticles < 5) {
 76        MSG_DEBUG("Failed leptonic event cut");
 77        vetoEvent;
 78      }
 79      MSG_DEBUG("Passed leptonic event cut");
 80
 81      // Get beams and average beam momentum
 82      const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
 83      const double meanBeamMom = 0.5*(beams.first.p3().mod() + beams.second.p3().mod());
 84      MSG_DEBUG("Avg beam momentum = " << meanBeamMom);
 85
 86      // Final state to get particle spectra
 87      for (const Particle& p : apply<UnstableParticles>(event, "UFS").particles(Cuts::pid==111)) {
 88        double xE = p.E()/meanBeamMom;
 89        discfill("pi0", xE);
 90      }
 91      for (const Particle& p : apply<FinalState>(event, "FS").particles(Cuts::pid==22)) {
 92        double xE = p.E()/meanBeamMom;
 93        discfill("gamma", xE);
 94      }
 95    }
 96
 97    void discfill(const string& name, const double value) {
 98      string edge = "OTHER";
 99      const size_t idx = _axes[name].index(value);
100      if (idx && idx <= _edges[name].size())  edge = _edges[name][idx-1];
101      _h[name]->fill(edge);
102    }
103
104
105    /// Normalise histograms etc., after the run
106    void finalize() {
107
108      const double fact = sqr(sqrtS())/GeV2*crossSection()/microbarn/sumOfWeights();
109      scale(_h, fact);
110      for( auto & hist : _h) {
111        for(auto & b: hist.second->bins()) {
112          const size_t idx = b.index();
113          b.scaleW(1./_axes[hist.first].width(idx));
114        }
115      }
116
117    }
118
119    /// @}
120
121
122    /// @name Histograms
123    /// @{
124    map<string, BinnedHistoPtr<string>> _h;
125    map<string, YODA::Axis<double>> _axes;
126    map<string, vector<string>> _edges;
127    /// @}
128
129
130  };
131
132
133  RIVET_DECLARE_PLUGIN(CELLO_1983_I191415);
134
135
136}