Rivet analyses

title: ATLAS_2011_I916832

Inclusive isolated diphoton analysis

Experiment: ATLAS (LHC 7TeV)

Inspire ID: 916832

Status: VALIDATED

Authors: - Giovanni Marchiori

References: - Expt page: ATLAS-STDM-2011-05 - arXiv: 1107.0581

Beams: p+ p+

Beam energies: (3500.0, 3500.0)GeV

Run details: - Inclusive diphoton + $X$ events at $\sqrt{s} = 7$~TeV.

A measurement of the cross section for inclusive isolated photon production at $\sqrt{s} = 7$ TeV. The measurement is done in bins of $M_{\gamma\gamma}$, $p_{T\gamma\gamma}$, and $\Delta\phi_{\gamma\gamma}$, for isolated photons with $|\eta|<2.37$ and $E_T^\gamma>16$~GeV. The measurement uses 37 pb$^{-1}$ of integrated luminosity collected with the ATLAS detector.

Source code:ATLAS_2011_I916832.cc

```c++ // -- C++ --

include "Rivet/Analysis.hh"

include "Rivet/Projections/FastJets.hh"

include "Rivet/Projections/FinalState.hh"

include "Rivet/Projections/IdentifiedFinalState.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_I916832 : public Analysis { public:

/// Constructor
RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2011_I916832);


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

  FastJets fj(fs, JetAlg::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() / picobarn / sumOfWeights());
  scale(_h_pT, crossSection() / picobarn / sumOfWeights());
  scale(_h_dPhi, crossSection() / picobarn / 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_I916832, ATLAS_2011_S9120807);

} ```

Aliases: - ATLAS_2011_S9120807