Rivet analyses

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

namespace Rivet {

    /// @brief Dielectron production in pp and p-Pb collisions at 5.02 TeV
    class ALICE_2020_I1797621 : public Analysis {
  public:

    /// Constructor
    RIVET_DEFAULT_ANALYSIS_CTOR(ALICE_2020_I1797621);

    /// @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);

      const ParticlePair& pp = beams();
      _hasLeadBeam = pp.first.pid() == PID::LEAD || pp.second.pid() == PID::LEAD;

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

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

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

      FinalState afs1 = FinalState(Cuts::pid == PID::ELECTRON && nocuts);
      declare(afs1, "AllElecs");

      FinalState afs2 = FinalState(Cuts::pid == -PID::ELECTRON && nocuts);
      declare(afs2, "AllPosis");

      // Book histograms
      book(_h["m_ee_pp"], 1,1,1);
      book(_h["m_ee_pPb"], 2, 1, 1);
      book(_h["pT_ee_low_pp"], 3, 1, 1);
      book(_h["pT_ee_low_pPb"], 4, 1, 1);
      book(_h["pT_ee_high_pp"], 5, 1, 1);
      book(_h["pT_ee_high_pPb"], 6, 1, 1);
    }

    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;
      }

      // Comes from Drell-Yan or virtual photon
      if (_photon && pa.isDirect()) {
        return true;
      }

      // Comes from bottom hadron decay
      if (_bottom && pa.fromBottom()) {
        return true;
      }

      // Comes from charm hadron decay
      if (_charm && pa.fromCharm()) {
        return true;
      }

      // By default there are no restrictions -> accept all
      if (!_charm && !_bottom && !_photon) {
        return true;
      }

      return false;
    }

    // With analysis cuts
    void bin_em(const Particle &l1, const Particle &l2,double weight, const Event& event) {
      // Combine the dilepton pair
      const FourMomentum pll = l1.mom() + l2.mom();

      if (pll.pt() < 8 * GeV) {
        if (!_hasLeadBeam) _h["m_ee_pp"]->fill(pll.mass() / GeV,weight);
        else               _h["m_ee_pPb"]->fill(pll.mass() / GeV,weight);
      }
      if (0.5 < pll.mass() && pll.mass() < 1.10) {
        if (!_hasLeadBeam)  _h["pT_ee_low_pp"]->fill(pll.pt() / GeV,weight);
        else                _h["pT_ee_low_pPb"]->fill(pll.pt() / GeV,weight);
      }
      if (1.1 < pll.mass() && pll.mass() < 2.70) {
        if (!_hasLeadBeam)  _h["pT_ee_high_pp"]->fill(pll.pt() / GeV,weight);
        else                _h["pT_ee_high_pPb"]->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;

      // Now perform same-sign approximation. 2401.12875 Eq. (3.1)

      // Unlike sign contributions
      for (const Particle &l1 : elecs) {
        for (const Particle &l2 : posis) {
          if ( both_valid(l1,l2) ) {
            bin_em(l1,l2,1.0,event);
          }
        }
      }

      // 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,event); // only subtract 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,event); // only subtract 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;
    bool _hasLeadBeam = false;

    /// @}

    };

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