1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/ChargedFinalState.hh"
  4#include "Rivet/Projections/Beam.hh"
  5#include "Rivet/Projections/TriggerUA5.hh"
  6
  7namespace Rivet {
  8
  9
 10  /// UA5 \f$ \eta \f$ distributions at 200 and 900 GeV
 11  class UA5_1986_I233599 : public Analysis {
 12  public:
 13
 14    /// Constructor
 15    RIVET_DEFAULT_ANALYSIS_CTOR(UA5_1986_I233599);
 16
 17
 18    /// @name Analysis methods
 19    /// @{
 20
 21    /// Set up projections and histograms
 22    void init() {
 23      declare(TriggerUA5(), "Trigger");
 24      declare(Beam(), "Beams");
 25      declare(ChargedFinalState((Cuts::etaIn(-5.0, 5.0))), "CFS50");
 26
 27      // Histograms
 28      if (isCompatibleWithSqrtS(200.0*GeV)) {
 29        book(_hist_eta_nsd       ,1,1,1);
 30        book(_hist_eta_inelastic ,1,1,2);
 31        _hists_eta_nsd.resize(6);
 32        for (int i = 1; i <= 6; ++i) {
 33          _sumWn.push_back({});
 34          book(_sumWn.back(), "TMP/sumWn"+to_str(i));
 35          book(_hists_eta_nsd[i-1],2,1,i);
 36        }
 37      } else if (isCompatibleWithSqrtS(900.0*GeV)) {
 38        book(_hist_eta_nsd       ,1,1,3);
 39        book(_hist_eta_inelastic ,1,1,4);
 40        _hists_eta_nsd.resize(9);
 41        for (int i = 1; i <= 9; ++i) {
 42          _sumWn.push_back({});
 43          book(_sumWn.back(), "TMP/sumWn"+to_str(i));
 44          book(_hists_eta_nsd[i-1],3,1,i);
 45        }
 46      }
 47      book(_sumWTrig, "sumWtrig");
 48      book(_sumWTrigNSD, "sumWtrigNSD");
 49
 50    }
 51
 52
 53    /// Fill eta histograms (in Nch bins)
 54    void analyze(const Event& event) {
 55      // Trigger
 56      const TriggerUA5& trigger = apply<TriggerUA5>(event, "Trigger");
 57      if (!trigger.sdDecision()) vetoEvent;
 58      const bool isNSD = trigger.nsdDecision();
 59
 60      // Get the index corresponding to the max Nch range histo/sum(w) vector index
 61      const ChargedFinalState& cfs50 = apply<ChargedFinalState>(event, "CFS50");
 62      const int numP = cfs50.size();
 63      const int ni = (int)floor(static_cast<float>(numP-2)/10.0);
 64      const int num_idx = min(ni, (int)_sumWn.size()-1);
 65      MSG_TRACE("Multiplicity index: " << numP << " charged particles -> #" << num_idx);
 66
 67      // Update weights
 68      _sumWTrig->fill();
 69      if (isNSD) {
 70        _sumWTrigNSD->fill();
 71        if (num_idx >= 0) _sumWn[num_idx]->fill();
 72      }
 73
 74      // Fill histos
 75      for (const Particle& p : cfs50.particles()) {
 76        const double eta = p.abseta();
 77        _hist_eta_inelastic->fill(eta);
 78        if (isNSD) {
 79          _hist_eta_nsd->fill(eta);
 80          if (num_idx >= 0) _hists_eta_nsd[num_idx]->fill(eta);
 81        }
 82      }
 83    }
 84
 85
 86    /// Scale histos
 87    void finalize() {
 88      MSG_DEBUG("sumW_NSD,inel = " << _sumWTrigNSD->val() << ", " << _sumWTrig->val());
 89      scale(_hist_eta_nsd, 0.5 / *_sumWTrigNSD);
 90      scale(_hist_eta_inelastic, 0.5 / *_sumWTrig);
 91      //
 92      for (size_t i = 0; i < _hists_eta_nsd.size(); ++i) {
 93        MSG_DEBUG("sumW[n] = " << _sumWn[i]->val());
 94        scale(_hists_eta_nsd[i], 0.5 / *_sumWn[i]);
 95      }
 96    }
 97
 98
 99  private:
100
101    /// @name Weight counters
102    /// @{
103    CounterPtr _sumWTrig;
104    CounterPtr _sumWTrigNSD;
105    vector<CounterPtr> _sumWn;
106    /// @}
107
108    /// @name Histograms
109    /// @{
110    Histo1DPtr _hist_eta_nsd;
111    Histo1DPtr _hist_eta_inelastic;
112    vector<Histo1DPtr> _hists_eta_nsd;
113    /// @}
114
115  };
116
117
118
119  RIVET_DECLARE_ALIASED_PLUGIN(UA5_1986_I233599, UA5_1986_S1583476);
120
121}