Rivet Analyses Reference

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131

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

namespace Rivet {


  /// Inclusive isolated prompt photon analysis with 2011 LHC data
  class ATLAS_2013_I1263495 : public Analysis {
  public:

    /// Constructor
    ATLAS_2013_I1263495()
      : Analysis("ATLAS_2013_I1263495"),
        _eta_bins{0.00, 1.37, 1.52, 2.37},
        _eta_bins_areaoffset{0.0, 1.5, 3.0}
    {    }


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

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

      // Consider the final state jets for the energy density calculation
      FastJets fj(fs, FastJets::KT, 0.5);
      fj.useJetArea(new fastjet::AreaDefinition(fastjet::VoronoiAreaSpec()));
      declare(fj, "KtJetsD05");

      // Consider the leading pt photon with |eta| < 2.37 and pT > 100 GeV
      LeadingParticlesFinalState photonfs(FinalState(Cuts::abseta < 2.37 && Cuts::pT > 100*GeV));
      photonfs.addParticleId(PID::PHOTON);
      declare(photonfs, "LeadingPhoton");

      // Book the dsigma/dEt (in eta bins) histograms
      book(_h_Et_photon[0], 1, 1, 1);
      book(_h_Et_photon[2], 2, 1, 1);

      // Book the dsigma/d|eta| histogram
      book(_h_eta_photon, 3, 1, 1);
    }


    /// Return eta bin for either dsigma/dET histogram (area_eta=false) or energy density correction (area_eta=true)
    size_t _getEtaBin(double eta_w, bool area_eta) const {
      const double eta = fabs(eta_w);
      if (!area_eta) {
        return binIndex(eta, _eta_bins);
      } else {
        return binIndex(eta, _eta_bins_areaoffset);
      }
    }


    /// Perform the per-event analysis
    void analyze(const Event& event) {
      // Retrieve leading photon
      Particles photons = apply<LeadingParticlesFinalState>(event, "LeadingPhoton").particles();
      if (photons.size() != 1) vetoEvent;
      const Particle& leadingPhoton = photons[0];

      // Veto events with photon in ECAL crack
      if (inRange(leadingPhoton.abseta(), 1.37, 1.52)) vetoEvent;

      // Compute isolation energy in cone of radius .4 around photon (all particles)
      FourMomentum mom_in_EtCone;
      Particles fs = apply<FinalState>(event, "FS").particles();
      for (const Particle& p : fs) {
        // Check if it's outside the cone of 0.4
        if (deltaR(leadingPhoton, p) >= 0.4) continue;
        // Don't count particles in the 5x7 central core
        if (deltaEta(leadingPhoton, p) < .025*5.0*0.5 &&
            deltaPhi(leadingPhoton, p) < (PI/128.)*7.0*0.5) continue;
        // Increment isolation energy
        mom_in_EtCone += p.momentum();
      }

      // Get the area-filtered jet inputs for computing median energy density, etc.
      vector<double> ptDensity;
      vector< vector<double> > ptDensities(_eta_bins_areaoffset.size()-1);
      FastJets fast_jets =apply<FastJets>(event, "KtJetsD05");
      const shared_ptr<fastjet::ClusterSequenceArea> clust_seq_area = fast_jets.clusterSeqArea();
      for (const Jet& jet : fast_jets.jets()) {
        const double area = clust_seq_area->area(jet);
        if (area > 1e-4 && jet.abseta() < _eta_bins_areaoffset.back())
          ptDensities.at( _getEtaBin(jet.abseta(), true) ).push_back(jet.pT()/area);
      }
      // Compute the median energy density, etc.
      for (size_t b = 0; b < _eta_bins_areaoffset.size()-1; b++) {
        ptDensity += ptDensities[b].empty() ? 0 : median(ptDensities[b]);
      }
      // Compute the isolation energy correction (cone area*energy density)
      const double etCone_area = PI*sqr(0.4) - (7.0*.025)*(5.0*PI/128.);
      const double correction = ptDensity[_getEtaBin(leadingPhoton.abseta(), true)]*etCone_area;

      // Apply isolation cut on area-corrected value
      if (mom_in_EtCone.Et() - correction > 7*GeV) vetoEvent;

      // Fill histograms
      const size_t eta_bin = _getEtaBin(leadingPhoton.abseta(), false);
      _h_Et_photon[eta_bin]->fill(leadingPhoton.Et());
      _h_eta_photon->fill(leadingPhoton.abseta());
    }


    /// Normalise histograms etc., after the run
    void finalize() {
      for (size_t i = 0; i < _eta_bins.size()-1; i++) {
        if (fuzzyEquals(_eta_bins[i], 1.37)) continue;
        scale(_h_Et_photon[i], crossSection()/picobarn/sumOfWeights());
      }
      scale(_h_eta_photon, crossSection()/picobarn/sumOfWeights());
    }


  private:

    Histo1DPtr _h_Et_photon[3];
    Histo1DPtr _h_eta_photon;

    vector<double> _eta_bins, _eta_bins_areaoffset;

  };


  RIVET_DECLARE_PLUGIN(ATLAS_2013_I1263495);

}