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  gamma, pi0 and eta spectra at 14, 22.5 and 34.4 GeV
 11  class JADE_1985_I213948 : public Analysis {
 12  public:
 13
 14    /// Constructor
 15    RIVET_DEFAULT_ANALYSIS_CTOR(JADE_1985_I213948);
 16
 17
 18    /// @name Analysis methods
 19    /// @{
 20
 21    /// Book histograms and initialise projections before the run
 22    void init() {
 23      declare(Beam(), "Beams");
 24      declare(FinalState(), "FS");
 25      declare(UnstableParticles(), "UFS");
 26      // find the beam energy
 27      int ioff=-1;
 28      if (isCompatibleWithSqrtS(34.5*GeV))  ioff=0;
 29      else if(isCompatibleWithSqrtS(22.5*GeV)) 	ioff=1;
 30      else if(isCompatibleWithSqrtS(14.0*GeV)) 	ioff=2;
 31      else {
 32      	MSG_ERROR("Beam energy " << sqrtS() << " not supported!");
 33      }
 34      // book histos
 35      book(_h["gamma"], ioff+1, 1, 1);
 36      book(_h["pi0"],   ioff+4, 1, 1);
 37      if (ioff==0) {
 38        book(_h["eta"], 7, 1, 1);
 39        _axes.emplace_back<vector<double>>({0.00, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.10,
 40                                            0.11, 0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 0.19, 0.20,
 41                                            0.22, 0.24, 0.26, 0.28, 0.30, 0.32, 0.34, 0.36, 0.38, 0.40,
 42                                            0.45, 0.50, 0.55, 0.60, 0.70, 0.80, 0.90}); // d01
 43        _axes.emplace_back<vector<double>>({0.011, 0.023, 0.035, 0.047, 0.057, 0.069, 0.079, 0.083,
 44                                            0.093, 0.107, 0.113, 0.119, 0.143, 0.167, 0.191, 0.239}); // d04
 45        _axes.emplace_back<vector<double>>({0.040, 0.056, 0.088, 0.274}); // d07
 46      }
 47      else if (ioff==2) {
 48        _axes.emplace_back<vector<double>>({0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.10,
 49                                            0.12, 0.14, 0.16, 0.18, 0.20, 0.25, 0.40, 0.60, 1.00}); // d02
 50        _axes.emplace_back<vector<double>>({0.018, 0.036, 0.054, 0.072, 0.108, 0.144, 0.216, 0.360}); // d05
 51      }
 52      else {
 53        _axes.emplace_back<vector<double>>({0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.10,
 54                                            0.12, 0.14, 0.16, 0.18, 0.20, 0.25, 0.40, 0.60, 1.00}); //d03
 55        _axes.emplace_back<vector<double>>({0.030, 0.058, 0.086, 0.114, 0.170, 0.226, 0.338, 0.562}); // d06
 56      }
 57    }
 58
 59
 60    /// Perform the per-event analysis
 61    void analyze(const Event& event) {
 62      if (_edges.empty()) {
 63        _edges.resize(_axes.size());
 64        _edges[0] = _h["gamma"]->xEdges();
 65        _edges[1] = _h["pi0"]->xEdges();
 66        if (_h.count("eta"))  _edges[3] = _h["eta"]->xEdges();
 67      }
 68      // Get beams and average beam momentum
 69      const ParticlePair& beams = apply<Beam>(event, "Beams").beams();
 70      const double meanBeamMom = 0.5*(beams.first.p3().mod()+beams.second.p3().mod());
 71      // gamma
 72      const FinalState& fs = apply<FinalState>(event, "FS");
 73      for (const Particle& p : fs.particles(Cuts::pid==22)) {
 74      	const double xE = p.E()/meanBeamMom;
 75        _h["gamma"]->fill(map2string(xE, 0));
 76      }
 77      // pi0, eta
 78      const UnstableParticles& ufs = apply<UnstableParticles>(event, "UFS");
 79      for (const Particle& p : ufs.particles(Cuts::pid==111 or Cuts::pid==221)) {
 80      	double xE = p.E()/meanBeamMom;
 81      	if (p.pid()==111) {
 82          _h["pi0"]->fill(map2string(xE, 1));
 83        }
 84      	else if(_h.count("eta")) {
 85          _h["eta"]->fill(map2string(xE, 2));
 86        }
 87      }
 88    }
 89
 90    string map2string(const double val, const size_t k) const {
 91      const size_t idx = _axes[k].index(val);
 92      if (idx && idx <= _axes[k].numBins())  return _edges[k][idx-1];
 93      return "OTHER";
 94    }
 95
 96
 97    /// Normalise histograms etc., after the run
 98    void finalize() {
 99      const double sf = crossSection()*sqr(sqrtS())/microbarn/sumOfWeights();
100      scale(_h, sf);
101    }
102
103    /// @}
104
105
106    /// @name Histograms
107    /// @{
108    map<string, BinnedHistoPtr<string>> _h;
109    vector<YODA::Axis<double>> _axes;
110    vector<vector<string>> _edges;
111    /// @}
112
113
114  };
115
116
117  RIVET_DECLARE_PLUGIN(JADE_1985_I213948);
118
119}