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

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

namespace Rivet {


  class ATLAS_2010_S8914702 : public Analysis {
  public:

    /// Constructor
    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2010_S8914702);


    /// Book histograms and initialise projections before the run
    void init() {
      FinalState fs;
      declare(fs, "FS");

      FastJets fj(fs, FastJets::KT, 0.5);
      fj.useJetArea(new fastjet::AreaDefinition(fastjet::VoronoiAreaSpec()));
      declare(fj, "KtJetsD05");

      LeadingParticlesFinalState photonfs(FinalState(Cuts::abseta < 1.81 && Cuts::pT > 15*GeV));
      photonfs.addParticleId(PID::PHOTON);
      declare(photonfs, "LeadingPhoton");

      size_t hist_bin = 0;
      for (size_t i = 0; i < _eta_bins.size()-1; ++i) {
        if (fabs(_eta_bins[i] - 1.37) < .0001) continue;
        book(_h_Et_photon[i] ,1, 1, hist_bin+1);
        hist_bin += 1;
      }
    }


    size_t getEtaBin(double eta, bool area_eta) const {
      return (!area_eta) ? binIndex(fabs(eta), _eta_bins) : binIndex(fabs(eta), _eta_bins_areaoffset);
    }


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

      const Particle& leadingPhoton = photons[0];
      if (inRange(leadingPhoton.abseta(), 1.37, 1.52)) vetoEvent;
      const int eta_bin = getEtaBin(leadingPhoton.abseta(), false);

      const Particles& fs = apply<FinalState>(event, "FS").particles();
      FourMomentum mom_in_EtCone;
      for (const Particle& p : fs) {
        // Check if it's in the cone of .4
        if (deltaR(leadingPhoton, p) >= 0.4) continue;
        // Check if it's in the 5x7 central core
        if (fabs(deltaEta(leadingPhoton, p)) < .025*5.0*0.5 &&
            fabs(deltaPhi(leadingPhoton, p)) < (PI/128.)*7.0*0.5) continue;
        // Increment
        mom_in_EtCone += p.momentum();
      }
      MSG_DEBUG("Done with initial Et cone");

      // Get the jet pT densities
      vector< vector<double> > ptDensities(_eta_bins_areaoffset.size()-1);
      FastJets fastjets = apply<FastJets>(event, "KtJetsD05");
      const shared_ptr<fastjet::ClusterSequenceArea> clust_seq_area = fastjets.clusterSeqArea();
      for (const Jet& jet : fastjets.jets()) {
        const double area = clust_seq_area->area(jet); //< Implicit call to pseudojet()
        if (area > 1e-4 && jet.abseta() < _eta_bins_areaoffset.back()) {
          ptDensities.at(getEtaBin(jet.abseta(), true)) += jet.pT()/area;
        }
      }

      // Now compute the median energy densities
      vector<double> ptDensity;
      for (size_t b = 0; b < _eta_bins_areaoffset.size()-1; ++b) {
        ptDensity += ptDensities[b].empty() ? 0 : median(ptDensities[b]);
      }

      // Now figure out the correction
      const double ETCONE_AREA = PI*.4*.4 - (7.0*.025)*(5.0*PI/128.);
      const double correction = ptDensity[getEtaBin(leadingPhoton.abseta(), true)]*ETCONE_AREA;
      MSG_DEBUG("Jet area correction done");

      // Shouldn't need to subtract photon
      // NB. Using expected cut at hadron/particle level, not cut at reco level
      if (mom_in_EtCone.Et() - correction/*-leadingPhoton.Et()*/ > 4.0*GeV) vetoEvent;
      MSG_DEBUG("Passed isolation cut");

      // Fill histogram
      _h_Et_photon[eta_bin]->fill(leadingPhoton.Et()/GeV);
    }


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


  private:

    Histo1DPtr _h_Et_photon[6];

    const vector<double> _eta_bins = {0.00, 0.60, 1.37, 1.52, 1.81};
    const vector<double> _eta_bins_areaoffset = {0.0, 1.5, 3.0};

  };



  RIVET_DECLARE_ALIASED_PLUGIN(ATLAS_2010_S8914702, ATLAS_2010_I882463);

}