Rivet analyses

// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/DISKinematics.hh"
#include "Rivet/Projections/FastJets.hh"
#include "Rivet/Projections/FinalState.hh"

namespace Rivet {


  /// A configurable reconstruction of a DIS final state and calculation of standard variables
  class MC_DIS : public Analysis {
  public:

    /// Constructor
    RIVET_DEFAULT_ANALYSIS_CTOR(MC_DIS);


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

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

      // Convert analysis options to params
      ObjOrdering lsort = ObjOrdering::ENERGY;
      const string lsortopt = toUpper(retrieve(options(), string("LSort"), string("ENERGY")));
      if (toUpper(lsortopt) == "ETA") lsort = ObjOrdering::ETA;
      if (toUpper(lsortopt) == "ET") lsort = ObjOrdering::ET;

      LeptonReco lreco = LeptonReco::PROMPT_DRESSED;
      const string lrecoopt = toUpper(retrieve(options(), string("LReco"), string("PROMPT_DRESSED")));
      if (lrecoopt == "ALL") {
        lreco = LeptonReco::ALL;
      }
      else if (lrecoopt == "ALL_DRESSED") {
        lreco = LeptonReco::ALL_DRESSED;
      }
      else if (lrecoopt == "BARE" || lrecoopt == "PROMPT_BARE") {
        lreco = LeptonReco::PROMPT_BARE;
      }

      double beamundresstheta = 0.0;
      const string lundressopt = retrieve(options(), string("UndressTheta"), string("0.0"));
      beamundresstheta = stod(lundressopt);

      double isoldr = 0.0;
      const string isoldropt = retrieve(options(), "IsolDR", "0.0");
      isoldr = stod(isoldropt);

      double dressdr = 0.0;
      const string dressdropt = retrieve(options(), "DressDR", "0.0");
      dressdr = stod(dressdropt);


      // Initialise and register projections. The definition of the
      // scattered lepton can be influenced by the analysis options.
      declare(FinalState(), "FS");
      DISLepton lepton(lreco, lsort, beamundresstheta, isoldr, dressdr);
      declare(lepton, "Lepton");
      declare(DISKinematics(lepton), "Kinematics");


      // Histograms
      book(_h_charge_electron, "chargeelectron", 2, -1.0, 1.0);
      book(_h_x, "x", logspace(100, 1e-6, 1.0)); // powspace(100, 1e-6, 1.0, 0.01);
      //
      book(_h_eminuspz, "eminuspz", 240, 0.0, 60.0);
      book(_h_etot_remnant, "etotremnant", 100, 0.0, 1000.0);
      book(_h_pt_remnant, "ptremnant", 50, 0.0, 5.0);
      //
      book(_h_pttot, "pttot", 200, 0.0, 200.0);
      book(_h_pttot_leptons, "pttotleptons", 200, 0.0, 200.0);
      book(_h_pttot_hadrons, "pttothadrons", 200, 0.0, 200.0);
      book(_h_pttot_gamma, "pttotgamma", 200, 0.0, 200.0);
      //
      book(_h_e_electron, "eelectron", 240, 0.0, 60.0);
      book(_h_pt_electron, "ptelectron", 240, 0.0, 60.0);
      book(_h_y, "y", 100, 0.0, 1.0);
      book(_h_W2, "W2", 100, 0.0, 100000);
      book(_h_Q2, "Q2", logspace(100, 0.1, 1e5)); // powspace(100, 0.1, 1e5, 0.01);
      book(_h_gammahad, "gammahad", 180, 0.0, 180);
      book(_h_theta_electron, "thetaelectron", 180, 0.0, 180);
    }


    /// Perform the per-event analysis
    void analyze(const Event& event) {

      /// We analyze event an extract DIS kinematics
      const DISKinematics& dk = apply<DISKinematics>(event, "Kinematics");
      const DISLepton& dl = apply<DISLepton>(event, "Lepton");

      const double q2 = dk.Q2();
      const double x = dk.x();
      const double y = dk.y();
      const double W2 = dk.W2();
      const double gammahad = dk.gammahad() / degree;

      _h_x->fill(x);
      _h_y->fill(y);
      _h_W2->fill(W2);
      _h_Q2->fill(q2);
      _h_gammahad->fill(gammahad);
      _h_theta_electron->fill(dl.out().angle(dk.beamHadron().mom()) / degree);
      _h_e_electron->fill(dl.out().E());
      _h_pt_electron->fill(dl.out().pT());
      _h_charge_electron->fill(0.5 * (dl.in().charge() > 0 ? 1. : -1));

      double eminuspz = 0;
      double etot_remnant = 0;
      double pttot = 0;         /// transverse momentum of all particles but the scattered lepton
      double pttot_leptons = 0; /// transverse momentum of all leptons but the scattered one
      double pttot_hadrons = 0; /// transverse momentum of all hadrons
      double pttot_gamma = 0;   /// transverse momentum of all gammas
      const FinalState& fs = apply<FinalState>(event, "FS");
      for (const Particle& p : fs.particles()) {
        eminuspz += (p.E() + p.pz() * dl.pzSign());
        if (p.genParticle() == dl.out().genParticle()) continue;
        pttot += p.pT();
        if (p.isLepton()) pttot_leptons += p.pT();
        if (p.abspid() == PID::PHOTON) pttot_gamma += p.pT();
        if (p.isVisible() && !p.isLepton() && !(p.abspid() == PID::PHOTON)) pttot_hadrons += p.pT();

        if (p.abseta() < 6) continue;
        etot_remnant += p.E();
        _h_pt_remnant->fill(p.pT());
      }
      _h_eminuspz->fill(eminuspz);
      _h_etot_remnant->fill(etot_remnant);
      _h_pttot->fill(pttot);
      _h_pttot_leptons->fill(pttot_leptons);
      _h_pttot_hadrons->fill(pttot_hadrons);
      _h_pttot_gamma->fill(pttot_gamma);
    }


    /// Normalise histograms etc., after the run
    void finalize() {
      scale(_h_charge_electron, crossSection() / picobarn / sumOfWeights());
      normalize({_h_y, _h_W2, _h_x, _h_Q2, _h_gammahad, _h_eminuspz, _h_pt_remnant, _h_etot_remnant, _h_pttot,
                 _h_pttot_leptons, _h_pttot_hadrons, _h_pttot_gamma, _h_e_electron, _h_pt_electron,
                 _h_theta_electron});
    }

    /// @}


  private:

    /// @name Histograms
    /// @{
    Histo1DPtr _h_charge_electron;
    Histo1DPtr _h_y, _h_W2, _h_x, _h_Q2, _h_gammahad;
    Histo1DPtr _h_eminuspz;
    Histo1DPtr _h_pt_remnant;
    Histo1DPtr _h_etot_remnant;
    Histo1DPtr _h_pttot, _h_pttot_leptons, _h_pttot_hadrons, _h_pttot_gamma;
    Histo1DPtr _h_e_electron, _h_pt_electron, _h_theta_electron;
    /// @}
  };


  RIVET_DECLARE_PLUGIN(MC_DIS);

}