1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/UnstableParticles.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// @brief rho0, K*=, K*0 spectra at 29 GeV
  9  class HRS_1989_I276948 : public Analysis {
 10  public:
 11
 12    /// Constructor
 13    RIVET_DEFAULT_ANALYSIS_CTOR(HRS_1989_I276948);
 14
 15
 16    /// @name Analysis methods
 17    ///@{
 18
 19    /// Book histograms and initialise projections before the run
 20    void init() {
 21      declare(UnstableParticles(), "UFS");
 22      book(_h["rho"], 1, 1, 1);
 23      book(_h["Kstar0"], 2, 1, 1);
 24      book(_h["Kstarp"], 2, 1, 2);
 25      book(_h_sig_rho    ,3,1,1);
 26      book(_h_sig_Kstar0 ,4,1,1);
 27      book(_h_mult_rho   ,3,1,3);
 28      book(_h_mult_Kstar0,4,1,3);
 29      _axes["rho"] = YODA::Axis<double>({0.0725, 0.0925, 0.114, 0.136, 0.159, 0.183, 0.207, 0.231,
 30                                         0.2555, 0.28, 0.3045, 0.3295, 0.354, 0.3785, 0.4035, 0.4285,
 31                                         0.4595, 0.5155, 0.6025, 0.702, 0.802, 0.902});
 32      _axes["Kstar0"] = YODA::Axis<double>({0.077, 0.119, 0.163, 0.2095, 0.2575, 0.3065,
 33                                            0.3555, 0.405, 0.455, 0.5, 0.6, 0.7, 1.0});
 34      _axes["Kstarp"] = _axes["Kstar0"];
 35    }
 36
 37
 38    /// Perform the per-event analysis
 39    void analyze(const Event& event) {
 40      if (_edges.empty()) {
 41        for (const auto& item : _h) {
 42          _edges[item.first] = item.second->xEdges();
 43        }
 44      }
 45      UnstableParticles ufs = apply<UnstableParticles>(event,"UFS");
 46      for (const Particle& p : ufs.particles(Cuts::abspid==113 || Cuts::abspid==323 || Cuts::abspid==313)) {
 47      	Vector3 mom3 = p.p3();
 48        const double energy = p.E();
 49      	double modp = mom3.mod();
 50      	double beta = modp/energy;
 51      	double xE = 2.*modp/sqrtS();
 52        if (p.pid()==113) {
 53          discfill("rho", xE, 1./beta);
 54          _h_sig_rho->fill(29);
 55          _h_mult_rho->fill(29);
 56        }
 57        else if (p.pid()==313) {
 58          discfill("Kstar0", xE, 1./beta);
 59          _h_sig_Kstar0->fill(29);
 60          _h_mult_Kstar0->fill(29);
 61        }
 62        else {
 63          discfill("Kstarp", xE, 1./beta);
 64        }
 65      }
 66    }
 67
 68    void discfill(const string& name, const double val, const double weight) {
 69      string edge = "OTHER";
 70      const size_t idx = _axes[name].index(val);
 71      if (idx && idx <= _edges[name].size())  edge = _edges[name][idx-1];
 72      _h[name]->fill(edge, weight);
 73    }
 74
 75
 76    /// Normalise histograms etc., after the run
 77    void finalize() {
 78      scale(_h, sqr(sqrtS())*crossSection()/microbarn/sumOfWeights());
 79      scale({_h_sig_rho, _h_sig_Kstar0}, crossSection()/picobarn/sumOfWeights());
 80      scale({_h_mult_rho, _h_mult_Kstar0}, 1./sumOfWeights());
 81      for( auto & hist : _h) {
 82        for(auto & b: hist.second->bins()) {
 83          const size_t idx = b.index();
 84          b.scaleW(1./_axes[hist.first].width(idx));
 85        }
 86      }
 87    }
 88
 89    ///@}
 90
 91
 92    /// @name Histograms
 93    ///@{
 94    BinnedHistoPtr<int> _h_sig_rho,_h_sig_Kstar0;
 95    BinnedHistoPtr<int> _h_mult_rho,_h_mult_Kstar0;
 96    map<string,BinnedHistoPtr<string>> _h;
 97    map<string,YODA::Axis<double>> _axes;
 98    map<string,vector<string>> _edges;
 99    ///@}
100
101
102  };
103
104
105  RIVET_DECLARE_PLUGIN(HRS_1989_I276948);
106
107}