Rivet analyses

// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/ChargedFinalState.hh"
#include "Rivet/Projections/FastJets.hh"
#include "Rivet/Projections/FinalState.hh"
#include "Rivet/Projections/IdentifiedFinalState.hh"
#include "Rivet/Projections/VetoedFinalState.hh"
#include "Rivet/Projections/VisibleFinalState.hh"

namespace Rivet {


  /// @author Peter Richardson
  class ATLAS_2012_I946427 : public Analysis {
  public:

    /// @name Constructors etc.
    /// @{

    /// Constructor
    ATLAS_2012_I946427()
        : Analysis("ATLAS_2012_I946427") { }

    /// @}


  public:

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

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

      // photons
      IdentifiedFinalState photonfs(Cuts::abseta < 1.81 && Cuts::pT > 25 * GeV);
      photonfs.acceptId(PID::PHOTON);
      declare(photonfs, "Photon");

      //
      FinalState fs;
      declare(fs, "FS");

      // Used for pTmiss
      declare(VisibleFinalState(Cuts::abseta < 4.9), "vfs");

      // Book histograms
      book(_count_SR, "count_SR", 1, 0., 1.);

      book(_hist_ET_photon, "hist_ET_photon", 48, 20., 500.);
      book(_hist_met, "hist_met", 100, 0., 500.);
    }


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


      // require at least 2 photons in final state
      Particles photons = apply<IdentifiedFinalState>(event, "Photon").particlesByPt();
      if (photons.size() < 2) {
        vetoEvent;
      }

      // Loop over photons and fill vector of isolated ones
      Particles fs = apply<FinalState>(event, "FS").particles();
      Particles isolated_photons;
      for (const Particle& photon : photons) {
        // remove photons in crack
        double eta_P = photon.eta();
        if (fabs(eta_P) >= 1.37 && fabs(eta_P) < 1.52) continue;

        double phi_P = photon.phi();

        FourMomentum mom_in_EtCone = -photon.momentum();
        for (const Particle& p : fs) {
          // check if it's in the cone of .2
          if (deltaR(eta_P, phi_P, p.eta(), p.phi()) >= 0.2) continue;
          mom_in_EtCone += p.momentum();
        }
        // apply isolation
        if (mom_in_EtCone.Et() > 5.) continue;

        // add photon to list of isolated ones
        isolated_photons.push_back(photon);
      }

      // need two isolated photons
      if (isolated_photons.size() < 2) {
        vetoEvent;
      }

      // pTmiss
      Particles vfs_particles = apply<VisibleFinalState>(event, "vfs").particles();
      FourMomentum pTmiss;
      for (const Particle& p : vfs_particles) {
        pTmiss -= p.momentum();
      }
      double eTmiss = pTmiss.pT();

      _hist_ET_photon->fill(isolated_photons[0].Et());
      _hist_met->fill(eTmiss);

      if (eTmiss > 125.) _count_SR->fill(0.5);
    }


    void finalize() {

      double norm = crossSection() / femtobarn * 1.07 / sumOfWeights();
      // these are number of events at 1.07fb^-1 per 10 GeV
      scale(_hist_ET_photon, 10. * norm);
      // these are number of events at 1.07fb^-1 per  5 GeV
      scale(_hist_met, 5. * norm);
      // these are number of events at 1.07fb^-1
      scale(_count_SR, norm);
    }

    /// @}


  private:

    Histo1DPtr _count_SR;
    Histo1DPtr _hist_ET_photon;
    Histo1DPtr _hist_met;
  };


  RIVET_DECLARE_PLUGIN(ATLAS_2012_I946427);

}