1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/ChargedFinalState.hh"
  4
  5namespace Rivet {
  6
  7
  8  /// Analysis looking at various distributions of final state particles
  9  class MC_FSPARTICLES : public Analysis {
 10  public:
 11
 12    /// Constructor
 13    RIVET_DEFAULT_ANALYSIS_CTOR(MC_FSPARTICLES);
 14
 15
 16    /// @name Analysis methods
 17    /// @{
 18
 19    /// Book histograms and initialise projections before the run
 20    void init() {
 21
 22      // Projections
 23      FinalState fs(Cuts::abseta < 5 && Cuts::pT > 500*MeV);
 24      declare(fs, "FS");
 25      declare(ChargedFinalState(fs), "CFS");
 26
 27      // Histograms
 28      /// @todo Choose E/pT ranged based on input energies... can't do anything about kin. cuts, though
 29      book(_histMult   , "Mult", 100, -0.5, 199.5);
 30      book(_histMultCh , "MultCh", 100, -0.5, 199.5);
 31
 32      book(_histPt   , "Pt", 300, 0, 30);
 33      book(_histPtCh , "PtCh", 300, 0, 30);
 34
 35      book(_histE   , "E", 100, 0, 200);
 36      book(_histECh , "ECh", 100, 0, 200);
 37
 38      book(_histEtaSumEt , "EtaSumEt", 25, 0, 5);
 39
 40      book(_histEta    , "Eta", 50, -5, 5);
 41      book(_histEtaCh  , "EtaCh", 50, -5, 5);
 42      book(_tmphistEtaPlus, "TMP/EtaPlus", 25, 0, 5);
 43      book(_tmphistEtaMinus, "TMP/EtaMinus", 25, 0, 5);
 44      book(_tmphistEtaChPlus, "TMP/EtaChPlus", 25, 0, 5);
 45      book(_tmphistEtaChMinus, "TMP/EtaChMinus", 25, 0, 5);
 46
 47      book(_histRapidity    , "Rapidity", 50, -5, 5);
 48      book(_histRapidityCh  , "RapidityCh", 50, -5, 5);
 49      book(_tmphistRapPlus, "TMP/RapPlus", 25, 0, 5);
 50      book(_tmphistRapMinus, "TMP/RapMinus", 25, 0, 5);
 51      book(_tmphistRapChPlus, "TMP/RapChPlus", 25, 0, 5);
 52      book(_tmphistRapChMinus, "TMP/RapChMinus", 25, 0, 5);
 53
 54      book(_histPhi    , "Phi", 50, 0, TWOPI);
 55      book(_histPhiCh  , "PhiCh", 50, 0, TWOPI);
 56
 57      book(_histEtaPMRatio , "EtaPMRatio");
 58      book(_histEtaChPMRatio , "EtaChPMRatio");
 59      book(_histRapidityPMRatio , "RapidityPMRatio");
 60      book(_histRapidityChPMRatio , "RapidityChPMRatio");
 61    }
 62
 63
 64    /// Perform the per-event analysis
 65    void analyze(const Event& event) {
 66      // Charged + neutral final state
 67      const FinalState& fs = apply<FinalState>(event, "FS");
 68      MSG_DEBUG("Total multiplicity = " << fs.size());
 69      _histMult->fill(fs.size());
 70      for (const Particle& p : fs.particles()) {
 71        _histEta->fill(p.eta());
 72        _histEtaSumEt->fill(p.abseta(), p.Et());
 73        (p.eta() > 0 ? _tmphistEtaPlus : _tmphistEtaMinus)->fill(p.abseta());
 74        //
 75        _histRapidity->fill(p.rap());
 76        (p.rap() > 0 ? _tmphistRapPlus : _tmphistRapMinus)->fill(p.absrap());
 77        //
 78        _histPt->fill(p.pT()/GeV);
 79        _histE->fill(p.E()/GeV);
 80        _histPhi->fill(p.phi());
 81      }
 82
 83      // Same for the charged FS particles only
 84      const FinalState& cfs = apply<FinalState>(event, "CFS");
 85      MSG_DEBUG("Total charged multiplicity = " << cfs.size());
 86      _histMultCh->fill(cfs.size());
 87      for (const Particle& p : cfs.particles()) {
 88        _histEtaCh->fill(p.eta());
 89        (p.eta() > 0 ? _tmphistEtaChPlus : _tmphistEtaChMinus)->fill(p.abseta());
 90        //
 91        _histRapidityCh->fill(p.rap());
 92        (p.rap() > 0 ? _tmphistRapChPlus : _tmphistRapChMinus)->fill(p.absrap());
 93        //
 94        _histPtCh->fill(p.pT()/GeV);
 95        _histECh->fill(p.E()/GeV);
 96        _histPhiCh->fill(p.phi());
 97      }
 98
 99    }
100
101
102    /// Finalize
103    void finalize() {
104      normalize(_histMult); normalize(_histEta); normalize(_histRapidity);
105      normalize(_histPt); normalize(_histE); normalize(_histPhi);
106      normalize(_histMultCh); normalize(_histEtaCh); normalize(_histRapidityCh);
107      normalize(_histPtCh); normalize(_histECh); normalize(_histPhiCh);
108      divide(_tmphistEtaPlus, _tmphistEtaMinus, _histEtaPMRatio);
109      divide(_tmphistEtaChPlus, _tmphistEtaChMinus, _histEtaChPMRatio);
110      divide(_tmphistRapPlus, _tmphistRapMinus, _histRapidityPMRatio);
111      divide(_tmphistRapChPlus, _tmphistRapChMinus, _histRapidityChPMRatio);
112    }
113
114    /// @}
115
116
117  private:
118
119    /// @name Histograms
120    /// @{
121    Histo1DPtr _histMult, _histEta, _histRapidity, _histPt, _histE, _histPhi;
122    Histo1DPtr _histMultCh,  _histEtaCh, _histRapidityCh, _histPtCh, _histECh, _histPhiCh;
123    Profile1DPtr _histEtaSumEt;
124    Estimate1DPtr _histEtaPMRatio, _histEtaChPMRatio, _histRapidityPMRatio, _histRapidityChPMRatio;
125    /// @}
126
127    /// @name Temporary histos used to calculate +/- rapidity ratio plots
128    /// @{
129    Histo1DPtr _tmphistEtaPlus, _tmphistEtaMinus, _tmphistEtaChPlus, _tmphistEtaChMinus;
130    Histo1DPtr _tmphistRapPlus, _tmphistRapMinus, _tmphistRapChPlus, _tmphistRapChMinus;
131    /// @}
132
133  };
134
135
136  RIVET_DECLARE_PLUGIN(MC_FSPARTICLES);
137
138}