Rivet analyses

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

namespace Rivet {

  /// @brief Measurement of dielectron production in pp at 7 TeV
  class ALICE_2018_I1672792 : public Analysis {
  public:

    /// Constructor
    RIVET_DEFAULT_ANALYSIS_CTOR(ALICE_2018_I1672792);

    /// @name Analysis methods
    //@{

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

      _photon = getOption<int>("PHOTON_Z", 0);
      _bottom = getOption<int>("BOTTOM", 0);
      _charm = getOption<int>("CHARM", 0);
      _open = getOption<int>("OPEN", 0);
      _apid = getOption<int>("APID", 0);

      Cut cuts = (Cuts::pT > 0.2 * GeV) && (Cuts::abseta < 0.8);

      FinalState fs1 = FinalState(Cuts::pid == PID::ELECTRON && cuts);
      declare(fs1, "Elecs");

      FinalState fs2 = FinalState(Cuts::pid == -PID::ELECTRON && cuts);
      declare(fs2, "Posis");


      // Book histograms
      book(_h["m_ee_pp"], 1, 1, 1); // Fig 9
      book(_h["pT_ee_VLMR_pp"], 2, 1, 1); // Fig 10 left
      book(_h["pT_ee_LMR_1_pp"], 4, 1, 1); // Fig 11 left
      book(_h["pT_ee_LMR_2_pp"], 5, 1, 1); // Fig 11 right
      book(_h["pT_ee_IMR_pp"], 7, 1, 1); // FIg 13 left and right (same data)
      book(_h["pT_ee_HMR_pp"], 9, 1, 1); // Fig 16 left
      //Fig 17? why not on hepdata...
    }

    bool valid(const Particle& pa) const {
      // If we are only iterested in open heavy-flavour we reject electrons coming from quarkonia
      if (_open && pa.hasParentWith([&](const Particle& p){
          return PID::isQuarkonium(p.pid());
        })) {
          return false;
      }

      if (_photon && pa.isDirect()) {
        return true;
      }

      if (_bottom && pa.fromBottom()) {
        return true;
      }

      if (_charm && pa.fromCharm()) {
        return true;
      }

      if (!_charm && !_bottom && !_photon) {
        return true;
      }
      return false;
    }

    void bin_em(const Particle& l1, const Particle& l2, double weight= 1.0) {
      // Get the true dilepton pair
      const FourMomentum pll = l1.mom() + l2.mom();

      if (pll.pt() < 8 * GeV) {
        _h["m_ee_pp"]->fill(pll.mass() / GeV,weight);
      }
      if (pll.mass() < 0.14) {
        _h["pT_ee_VLMR_pp"]->fill(pll.pt() / GeV, weight);
      }
      if (0.14 < pll.mass() && pll.mass() < 0.7) {
        _h["pT_ee_LMR_1_pp"]->fill(pll.pt() / GeV, weight);
      }
      if (0.7 < pll.mass() && pll.mass() < 1.1) {
        _h["pT_ee_LMR_2_pp"]->fill(pll.pt() / GeV, weight);
      }
      if (1.1 < pll.mass() && pll.mass() < 2.7) {
        _h["pT_ee_IMR_pp"]->fill(pll.pt() / GeV, weight);
      }
      if (2.7 < pll.mass() && pll.mass() < 3.3) {
        _h["pT_ee_HMR_pp"]->fill(pll.pt() / GeV, weight);
      }
    }
    /**
     * Check if both leptons are valid and not the same
    */
    bool both_valid(const Particle &l1, const Particle &l2) const {
      if (!l1.isSame(l2) && valid(l1) && valid(l2)) {
        // APID filter, can not be moved to valid() since it applies to either one or the other lepton
        if (!_apid || (
            (l1.hasAncestorWith(Cuts::abspid == _apid)) ||
            (l2.hasAncestorWith(Cuts::abspid == _apid)))) {
            return true;
        }
      }
      return false;
    }

    /// Perform the per-event analysis
    void analyze(const Event &event) {
      Particles elecs, posis;
      elecs = apply<ParticleFinder>(event, "Elecs").particlesByPt();
      posis = apply<ParticleFinder>(event, "Posis").particlesByPt();

      Particle l1;
      Particle l2;

      // Unlike sign contributions
      for (const Particle &l1 : elecs) {
        for (const Particle &l2 : posis) {
          if ( both_valid(l1,l2) ) {
            bin_em(l1,l2);
          }
        }
      }
      // Like sign contributions are subtracted through the negative weight
      for (const Particle &l1 : posis) {
        for (const Particle &l2 : posis) {
          if ( both_valid(l1,l2) ) {
            bin_em(l1,l2,-0.5); // only subbtract half since loop double counts
          }
        }
      }
      for (const Particle &l1 : elecs) {
        for (const Particle &l2 : elecs) {
          if ( both_valid(l1,l2) ) {
             bin_em(l1,l2,-0.5); // only subbtract half since loop double counts
          }
        }
      }
    }

    /// Normalise histograms etc., after the run
    void finalize() {
      scale(_h, crossSectionPerEvent() / (millibarn));
    }

    //@}

    /// @name Histograms
    //@{
    map<string, Histo1DPtr> _h;

    int _photon = 0;
    int _bottom = 0;
    int _charm = 0;
    int _open = 0;
    int _apid = 0;

    //@}
  };

  // The hook for the plugin system
  RIVET_DECLARE_PLUGIN(ALICE_2018_I1672792);
}