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.0465, 0.058, 0.0695, 0.081,
 44                                            0.093, 0.1045, 0.119, 0.1395, 0.1625, 0.1915, 0.2265}); // d04
 45        _axes.emplace_back<vector<double>>({0.040, 0.056, 0.088, 0.274}); // d07
 46      }
 47      else if (ioff==1) {
 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[2] = _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      vector<string> st = {"gamma","pi0","eta"};
102      for(unsigned int ix=0;ix<3;++ix) {
103        if(_h.count(st[ix])==0) continue;
104        for(auto & b: _h[st[ix]]->bins()) {
105          const size_t idx = b.index();
106          b.scaleW(1./_axes[ix].width(idx));
107        }
108      } 
109    }
110
111    /// @}
112
113
114    /// @name Histograms
115    /// @{
116    map<string, BinnedHistoPtr<string>> _h;
117    vector<YODA::Axis<double>> _axes;
118    vector<vector<string>> _edges;
119    /// @}
120
121
122  };
123
124
125  RIVET_DECLARE_PLUGIN(JADE_1985_I213948);
126
127}