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
132
133
134
135
136
137
138
139
140
141
142

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

namespace Rivet {


  /// @brief Measurement of isolated diphoton + X differential cross-sections
  ///
  /// Inclusive isolated gamma gamma cross-sections, differential in M(gg), pT(gg), dphi(gg)
  ///
  /// @author Giovanni Marchiori
  class ATLAS_2011_S9120807 : public Analysis {
  public:

    /// Constructor
    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2011_S9120807);


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

      IdentifiedFinalState photonfs(Cuts::abseta < 2.37 && Cuts::pT > 16*GeV);
      photonfs.acceptId(PID::PHOTON);
      declare(photonfs, "Photon");

      book(_h_M    ,1, 1, 1);
      book(_h_pT   ,2, 1, 1);
      book(_h_dPhi ,3, 1, 1);
    }


    size_t getEtaBin(double eta) const {
      const double aeta = fabs(eta);
      return binIndex(aeta, _eta_bins_areaoffset);
    }


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

      // Compute jet pT densities
      vector<double> ptDensity;
      vector< vector<double> > ptDensities(_eta_bins_areaoffset.size()-1);
      const shared_ptr<fastjet::ClusterSequenceArea> clust_seq_area = apply<FastJets>(event, "KtJetsD05").clusterSeqArea();
      for (const Jet& jet : apply<FastJets>(event, "KtJetsD05").jets()) {
        const double area = clust_seq_area->area(jet); // .pseudojet() called implicitly
        /// @todo Should be 1e-4?
        if (area > 10e-4 && jet.abseta() < _eta_bins_areaoffset.back()) {
          ptDensities.at(getEtaBin(jet.abseta())).push_back(jet.pT()/area);
        }
      }

      // Compute the median energy density
      for (size_t b = 0; b < _eta_bins_areaoffset.size()-1; ++b) {
        ptDensity += ptDensities[b].empty() ? 0 : median(ptDensities[b]);
      }

      // Loop over photons and fill vector of isolated ones
      Particles isolated_photons;
      for (const Particle& photon : photons) {

        // Remove photons in crack
        if (inRange(photon.abseta(), 1.37, 1.52)) continue;

        // Standard ET cone isolation
        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(photon, p) >= 0.4) continue;
          // Veto if it's in the 5x7 central core
          if (fabs(deltaEta(photon, p)) < 0.025*5.0*0.5 &&
              fabs(deltaPhi(photon, p)) < (M_PI/128.)*7.0*0.5) continue;
          // Increment isolation cone ET sum
          mom_in_EtCone += p.momentum();
        }

        // Now figure out the correction (area*density)
        const double ETCONE_AREA = M_PI*.4*.4 - (7.0*.025)*(5.0*M_PI/128.);
        const double correction = ptDensity[getEtaBin(photon.abseta())] * ETCONE_AREA;

        // 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 > 4.0*GeV) continue;

        // Add to list of isolated photons
        isolated_photons.push_back(photon);
      }

      // Require at least two isolated photons
      if (isolated_photons.size() < 2) vetoEvent;

      // Select leading pT pair
      std::sort(isolated_photons.begin(), isolated_photons.end(), cmpMomByPt);
      const FourMomentum y1 = isolated_photons[0].momentum();
      const FourMomentum y2 = isolated_photons[1].momentum();

      // Require the two photons to be separated (dR>0.4)
      if (deltaR(y1, y2) < 0.4) vetoEvent;

      const FourMomentum yy = y1 + y2;
      const double Myy = yy.mass()/GeV;
      const double pTyy = yy.pT()/GeV;
      const double dPhiyy = deltaPhi(y1.phi(), y2.phi());

      _h_M->fill(Myy);
      _h_pT->fill(pTyy);
      _h_dPhi->fill(dPhiyy);
    }


    /// Normalise histograms etc., after the run
    void finalize() {
      scale(_h_M, crossSection()/sumOfWeights());
      scale(_h_pT, crossSection()/sumOfWeights());
      scale(_h_dPhi, crossSection()/sumOfWeights());
    }


  private:

    Histo1DPtr _h_M, _h_pT, _h_dPhi;
    const vector<double> _eta_bins_areaoffset = {0.0, 1.5, 3.0};

  };


  RIVET_DECLARE_ALIASED_PLUGIN(ATLAS_2011_S9120807, ATLAS_2011_I916832);

}