Rivet analyses

// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/ChargedFinalState.hh"

namespace Rivet {


  /// Analysis looking at various distributions of final state particles
  class MC_FSPARTICLES : public Analysis {
  public:

    /// Constructor
    RIVET_DEFAULT_ANALYSIS_CTOR(MC_FSPARTICLES);


    /// @name Analysis methods
    /// @{

    /// Book histograms and initialise projections before the run
    void init() {

      // Get kinematic cut variables from the analysis options
      const double PTMIN = getOption<double>("PTMIN", 0.0) * GeV;
      const double PTMAX = getOption<double>("PTMAX", max(30.0 * GeV, 2 * PTMIN)) * GeV;

      // Projections
      FinalState fs(Cuts::abseta < 5 && Cuts::pT > PTMIN);
      declare(fs, "FS");
      declare(ChargedFinalState(fs), "CFS");

      // Histograms
      book(_histMult, "Mult", 100, -0.5, 199.5);
      book(_histMultCh, "MultCh", 100, -0.5, 199.5);

      book(_histPt, "Pt", 300, PTMIN / GeV, PTMAX / GeV);
      book(_histPtCh, "PtCh", 300, PTMIN / GeV, PTMAX / GeV);

      book(_histE, "E", logspace(100, max(PTMIN / GeV, 0.1), 10 * PTMAX));
      book(_histECh, "ECh", logspace(100, max(PTMIN / GeV, 0.1), 10 * PTMAX));

      book(_histEtaSumEt, "EtaSumEt", 25, 0, 5);

      book(_histEta, "Eta", 50, -5, 5);
      book(_histEtaCh, "EtaCh", 50, -5, 5);
      book(_tmphistEtaPlus, "TMP/EtaPlus", 25, 0, 5);
      book(_tmphistEtaMinus, "TMP/EtaMinus", 25, 0, 5);
      book(_tmphistEtaChPlus, "TMP/EtaChPlus", 25, 0, 5);
      book(_tmphistEtaChMinus, "TMP/EtaChMinus", 25, 0, 5);

      book(_histRapidity, "Rapidity", 50, -5, 5);
      book(_histRapidityCh, "RapidityCh", 50, -5, 5);
      book(_tmphistRapPlus, "TMP/RapPlus", 25, 0, 5);
      book(_tmphistRapMinus, "TMP/RapMinus", 25, 0, 5);
      book(_tmphistRapChPlus, "TMP/RapChPlus", 25, 0, 5);
      book(_tmphistRapChMinus, "TMP/RapChMinus", 25, 0, 5);

      book(_histPhi, "Phi", 50, 0, TWOPI);
      book(_histPhiCh, "PhiCh", 50, 0, TWOPI);

      book(_histEtaPMRatio, "EtaPMRatio");
      book(_histEtaChPMRatio, "EtaChPMRatio");
      book(_histRapidityPMRatio, "RapidityPMRatio");
      book(_histRapidityChPMRatio, "RapidityChPMRatio");
    }


    /// Perform the per-event analysis
    void analyze(const Event& event) {
      // Charged + neutral final state
      const FinalState& fs = apply<FinalState>(event, "FS");
      MSG_DEBUG("Total multiplicity = " << fs.size());
      _histMult->fill(fs.size());
      for (const Particle& p : fs.particles()) {
        _histEta->fill(p.eta());
        _histEtaSumEt->fill(p.abseta(), p.Et());
        (p.eta() > 0 ? _tmphistEtaPlus : _tmphistEtaMinus)->fill(p.abseta());
        //
        _histRapidity->fill(p.rap());
        (p.rap() > 0 ? _tmphistRapPlus : _tmphistRapMinus)->fill(p.absrap());
        //
        _histPt->fill(p.pT() / GeV);
        _histE->fill(p.E() / GeV);
        _histPhi->fill(p.phi());
      }

      // Same for the charged FS particles only
      const FinalState& cfs = apply<FinalState>(event, "CFS");
      MSG_DEBUG("Total charged multiplicity = " << cfs.size());
      _histMultCh->fill(cfs.size());
      for (const Particle& p : cfs.particles()) {
        _histEtaCh->fill(p.eta());
        (p.eta() > 0 ? _tmphistEtaChPlus : _tmphistEtaChMinus)->fill(p.abseta());
        //
        _histRapidityCh->fill(p.rap());
        (p.rap() > 0 ? _tmphistRapChPlus : _tmphistRapChMinus)->fill(p.absrap());
        //
        _histPtCh->fill(p.pT() / GeV);
        _histECh->fill(p.E() / GeV);
        _histPhiCh->fill(p.phi());
      }
    }


    /// Finalize
    void finalize() {
      normalize(_histMult);
      normalize(_histEta);
      normalize(_histRapidity);
      normalize(_histPt);
      normalize(_histE);
      normalize(_histPhi);
      normalize(_histMultCh);
      normalize(_histEtaCh);
      normalize(_histRapidityCh);
      normalize(_histPtCh);
      normalize(_histECh);
      normalize(_histPhiCh);
      divide(_tmphistEtaPlus, _tmphistEtaMinus, _histEtaPMRatio);
      divide(_tmphistEtaChPlus, _tmphistEtaChMinus, _histEtaChPMRatio);
      divide(_tmphistRapPlus, _tmphistRapMinus, _histRapidityPMRatio);
      divide(_tmphistRapChPlus, _tmphistRapChMinus, _histRapidityChPMRatio);
    }

    /// @}


  private:

    /// @name Histograms
    /// @{
    Histo1DPtr _histMult, _histEta, _histRapidity, _histPt, _histE, _histPhi;
    Histo1DPtr _histMultCh, _histEtaCh, _histRapidityCh, _histPtCh, _histECh, _histPhiCh;
    Profile1DPtr _histEtaSumEt;
    Estimate1DPtr _histEtaPMRatio, _histEtaChPMRatio, _histRapidityPMRatio, _histRapidityChPMRatio;
    /// @}

    /// @name Temporary histos used to calculate +/- rapidity ratio plots
    /// @{
    Histo1DPtr _tmphistEtaPlus, _tmphistEtaMinus, _tmphistEtaChPlus, _tmphistEtaChMinus;
    Histo1DPtr _tmphistRapPlus, _tmphistRapMinus, _tmphistRapChPlus, _tmphistRapChMinus;
    /// @}
  };


  RIVET_DECLARE_PLUGIN(MC_FSPARTICLES);

}