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 Add a short analysis description here
 11  class CELLO_1989_I276764 : public Analysis {
 12  public:
 13
 14    /// Constructor
 15    RIVET_DEFAULT_ANALYSIS_CTOR(CELLO_1989_I276764);
 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(Cuts::pid==111 || Cuts::pid==221), "UFS");
 28      // Book histograms
 29      for(unsigned int ix=0;ix<4;++ix)
 30        book(_h[ix], 2+ix, 1, 1);
 31    }
 32
 33
 34    /// Perform the per-event analysis
 35    void analyze(const Event& event) {
 36      if(_edges[0].empty()) {
 37        for(unsigned int ix=0;ix<4;++ix)
 38          _edges[ix] = _h[ix]->xEdges();
 39      }
 40      // at least 5 charged FS particles
 41      const FinalState& fs = apply<FinalState>(event, "FS");
 42      const size_t numParticles = fs.particles().size();
 43
 44      if (numParticles < 5) {
 45        MSG_DEBUG("Failed leptonic event cut");
 46        vetoEvent;
 47      }
 48      MSG_DEBUG("Passed leptonic event cut");
 49
 50      // Get beams and average beam momentum
 51      const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
 52      const double meanBeamMom = ( beams.first.p3().mod() +
 53                                   beams.second.p3().mod() ) / 2.0;
 54      MSG_DEBUG("Avg beam momentum = " << meanBeamMom);
 55
 56      // Final state to get particle spectra
 57      for (const Particle& p : apply<UnstableParticles>(event, "UFS").particles()) {
 58        double xE = p.E()/meanBeamMom;
 59        if(p.pid()==111) {
 60          _h[1]->fill( map2string(xE,1));
 61          _h[2]->fill( map2string(xE,2));
 62        }
 63        else
 64          _h[3]->fill( map2string(xE,3));
 65      }
 66      for (const Particle& p : apply<FinalState>(event, "FS").particles(Cuts::pid==22)) {
 67        _h[0]->fill(map2string(p.E()/meanBeamMom,0));
 68      }
 69    }
 70
 71    string map2string(const double val, const size_t axis) const {
 72      const size_t idx = _axes[axis].index(val);
 73      if (idx && idx <= _edges[axis].size())  return _edges[axis][idx-1];
 74      return "OTHER";
 75    }
 76
 77    /// Normalise histograms etc., after the run
 78    void finalize() {
 79      double fact = sqr(sqrtS())/GeV2*crossSection()/microbarn/sumOfWeights();
 80      for(unsigned int ix=0;ix<4;++ix) {
 81        scale(_h[ix], fact);
 82        for(auto & b: _h[ix]->bins()) {
 83          const size_t idx = b.index();
 84          b.scaleW(1./_axes[ix].width(idx));
 85        }
 86      }
 87    }
 88    
 89    /// @}
 90
 91
 92    /// @name Histograms
 93    /// @{
 94    BinnedHistoPtr<string> _h[4];
 95    vector<string> _edges[4];
 96    YODA::Axis<double> _axes[4] = { YODA::Axis<double>({0.01,0.02,0.03,0.04,0.05,0.06,0.07,0.08,0.09,0.10,0.11,0.12,0.13,0.14,0.15,0.16,
 97                                                            0.175,0.195,0.215,0.235,0.255,0.275,0.295,0.315,0.3535,0.394,0.434,0.4735,
 98                                                            0.514,0.555,0.5965,0.641,0.691,0.74,0.785,0.839,0.901}),
 99                                    YODA::Axis<double>({0.039,0.059,0.079,0.109,0.139,0.193,0.277,0.439}),
100                                    YODA::Axis<double>({0.2705,0.3415,0.412,0.482,0.552,0.6225,0.6935,0.7805,0.876,0.962}),
101                                    YODA::Axis<double>({0.0225,0.0775,0.1325,0.2875,0.4625,0.9175}) };
102    /// @}
103
104  };
105
106
107  RIVET_DECLARE_PLUGIN(CELLO_1989_I276764);
108
109
110}