1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/PromptFinalState.hh"
  4#include "Rivet/Projections/VisibleFinalState.hh"
  5#include "Rivet/Projections/VetoedFinalState.hh"
  6#include "Rivet/Math/MathUtils.hh"
  7#include "Rivet/Projections/FastJets.hh"
  8
  9namespace Rivet {
 10
 11
 12  /// @brief Isolated diphoton + X differential cross-sections with full run-2
 13  class ATLAS_2021_I1887997 : public Analysis {
 14  public:
 15
 16    // Constructor
 17    RIVET_DEFAULT_ANALYSIS_CTOR(ATLAS_2021_I1887997);
 18
 19
 20    // Book histograms and initialise projections before the run
 21    void init() {
 22
 23      // Calorimeter particles for photon isolation
 24      VisibleFinalState visFS;
 25      VetoedFinalState calo_fs(visFS);
 26      calo_fs.addVetoPairId(PID::MUON);
 27      declare(calo_fs, "calo_fs");
 28
 29      // Photons
 30      declare(PromptFinalState(Cuts::abspid == PID::PHOTON), "Photons");
 31
 32      // Jets for UE subtraction with jet-area method
 33      FastJets fj(FinalState(), JetAlg::KT, 0.5, JetMuons::NONE, JetInvisibles::NONE);
 34      fj.useJetArea(new fastjet::AreaDefinition(fastjet::VoronoiAreaSpec(0.9)));
 35      declare(fj, "KtJetsD05");
 36
 37      // Histograms
 38      book(_xs, "yy_xs");
 39      _observables = {"ph1_pt","ph2_pt","yy_cosTS","yy_m",
 40                      "yy_phiStar","yy_piMDphi","yy_pT","yy_pTt"};
 41      for (auto name : _observables) {
 42        book(_h[name], name);
 43      }
 44    }
 45
 46    // Perform the per-event analysis
 47    void analyze(const Event& event) {
 48
 49      // Require at least 2 prompt photons in final state
 50      Particles photons = apply<PromptFinalState>(event, "Photons").particlesByPt();
 51      if (photons.size() < 2) vetoEvent;
 52      photons.resize(2);
 53      const FourMomentum ph1 = photons[0];
 54      const FourMomentum ph2 = photons[1];
 55
 56      // Leading photon should have pT > 40 GeV, subleading > 30 GeV
 57      const double ph1_pt = ph1.pT();
 58      const double ph2_pt = ph2.pT();
 59      if (ph1_pt < 40.*GeV || ph2_pt < 30.*GeV) vetoEvent;
 60
 61      // Apply photon eta cuts
 62      iselect(photons, (Cuts::abseta < 2.37) && ( (Cuts::abseta <= 1.37) || (Cuts::abseta >= 1.52) ));
 63      if (photons.size() < 2) vetoEvent;
 64
 65      // Require the two photons to be separated in dR
 66      if (deltaR(ph1,ph2) < 0.4) vetoEvent;
 67
 68      // Get UE pt densities rho for subtraction later
 69      const vector<double> eta_bins = {0.0, 1.5, 3.0};
 70      vector<double> rho(eta_bins.size()-1, 0.0);
 71      FastJets ktjets = apply<FastJets>(event, "KtJetsD05");
 72      for (size_t ieta = 0; ieta < eta_bins.size()-1; ++ieta) {
 73        fastjet::Selector fjselector(fastjet::SelectorAbsRapRange(eta_bins[ieta], eta_bins[ieta+1]));
 74        double sigma, area;
 75        ktjets.clusterSeqArea()->get_median_rho_and_sigma(fjselector, true,
 76                                                          rho[ieta], sigma, area);
 77      }
 78
 79      // Loop over photons and require isolation
 80      const double isoRCone=0.2;
 81      for (const Particle& photon : photons) {
 82        // Compute calo isolation via particles within a cone around the photon
 83        const Particles fs = apply<VetoedFinalState>(event, "calo_fs").particles();
 84        FourMomentum mom_in_EtCone;
 85        for (const Particle& p : fs) {
 86          // Reject if not in cone
 87          if (deltaR(photon.momentum(), p.momentum()) > isoRCone)  continue;
 88          // Sum momentum
 89          mom_in_EtCone += p.momentum();
 90        }
 91        // subtract core photon
 92        mom_in_EtCone -= photon.momentum();
 93        // UE subtraction energy
 94        double UEpT = M_PI*sqr(isoRCone) * rho[binIndex(fabs(photon.eta()), eta_bins)];
 95        // Use photon if energy in isolation cone is low enough
 96        if (mom_in_EtCone.Et() - UEpT > 0.09*photon.momentum().pT()) vetoEvent;
 97      }
 98
 99      map<string, double> obs;
100      obs["ph1_pt"] = ph1_pt;
101      obs["ph2_pt"] = ph2_pt;
102      const FourMomentum yy = ph1 + ph2;
103      obs["yy_m"] = yy.mass();
104      obs["yy_pT"] = yy.pT();
105      obs["yy_piMDphi"] = PI-mapAngle0ToPi(ph1.phi() - ph2.phi());
106
107      obs["yy_cosTS"] = fabs(sinh(( ph1.eta() - ph2.eta() ))*2.0*ph1_pt*ph2_pt/sqrt(sqr(obs["yy_m"])+sqr(obs["yy_pT"]))/obs["yy_m"]); // Collins Soper frame
108      const double yy_cosTSLab = fabs(tanh(( ph1.eta() - ph2.eta() ) / 2.)); // Lab frame
109      const double sinthetastar_ = sqrt(1. - pow(yy_cosTSLab, 2));
110      obs["yy_phiStar"] = tan(0.5 * obs["yy_piMDphi"]) * sinthetastar_;
111
112      // a_t
113      const Vector3 t_hat(ph1.x()-ph2.x(), ph1.y()-ph2.y(), 0.);
114      const double factor = t_hat.mod();
115      const Vector3 t_hatx(t_hat.x()/factor, t_hat.y()/factor, t_hat.z()/factor);
116      const Vector3 At(ph1.x()+ph2.x(), ph1.y()+ph2.y(), 0.);
117      // Compute a_t transverse component with respect to t_hat
118      obs["yy_pTt"] = At.cross(t_hatx).mod();
119
120      // Fill fiducial cross section
121      _xs->fill();
122
123      // Fill histograms
124      for (auto name : _observables) {
125        _h[name]->fill(obs[name]);
126      }
127    }
128
129
130    // Normalise histograms etc., after the run
131    void finalize() {
132      const double sf = crossSection() / (picobarn * sumOfWeights());
133      scale(_xs, sf);
134      scale(_h, sf);
135    }
136
137
138  private:
139
140    CounterPtr _xs;
141    map<string, Histo1DPtr> _h;
142    vector<string> _observables;
143
144  };
145
146
147  RIVET_DECLARE_PLUGIN(ATLAS_2021_I1887997);
148
149}