1// -*- C++ -*-
  2#include "Rivet/Analysis.hh"
  3#include "Rivet/Projections/FinalState.hh"
  4#include "Rivet/Projections/IdentifiedFinalState.hh"
  5#include "Rivet/Tools/BinnedHistogram.hh"
  6
  7namespace Rivet {
  8
  9
 10  /// D0 direct photon pair production
 11  class D0_2010_S8570965 : public Analysis {
 12  public:
 13
 14    RIVET_DEFAULT_ANALYSIS_CTOR(D0_2010_S8570965);
 15
 16
 17    void init() {
 18      FinalState fs;
 19      declare(fs, "FS");
 20
 21      IdentifiedFinalState ifs(Cuts::abseta < 0.9 && Cuts::pT > 20*GeV);
 22      ifs.acceptId(PID::PHOTON);
 23      declare(ifs, "IFS");
 24
 25      book(_h_M ,1, 1, 1);
 26      book(_h_pT ,2, 1, 1);
 27      book(_h_dPhi ,3, 1, 1);
 28      book(_h_costheta ,4, 1, 1);
 29
 30      std::pair<double, double> M_ranges[] = { std::make_pair(30.0, 50.0),
 31                                               std::make_pair(50.0, 80.0),
 32                                               std::make_pair(80.0, 350.0) };
 33
 34      for (size_t i = 0; i < 3; ++i) {
 35        Histo1DPtr a,b,c;
 36        _h_pT_M.add(M_ranges[i].first, M_ranges[i].second, book(a, 5+3*i, 1, 1));
 37        _h_dPhi_M.add(M_ranges[i].first, M_ranges[i].second, book(b, 6+3*i, 1, 1));
 38        _h_costheta_M.add(M_ranges[i].first, M_ranges[i].second, book(c, 7+3*i, 1, 1));
 39      }
 40    }
 41
 42
 43    /// Perform the per-event analysis
 44    void analyze(const Event& event) {
 45      const double weight = 1.0;
 46
 47      Particles photons = apply<IdentifiedFinalState>(event, "IFS").particlesByPt();
 48      if (photons.size() < 2 ||
 49          (photons[0].pT() < 21.0*GeV)) {
 50        vetoEvent;
 51      }
 52
 53      // Isolate photons with ET_sum in cone
 54      Particles isolated_photons;
 55      Particles fs = apply<FinalState>(event, "FS").particles();
 56      for (const Particle& photon : photons) {
 57        double eta_P = photon.eta();
 58        double phi_P = photon.phi();
 59        double Etsum=0.0;
 60        for (const Particle& p : fs) {
 61          if (HepMCUtils::uniqueId(p.genParticle()) != HepMCUtils::uniqueId(photon.genParticle()) &&
 62              deltaR(eta_P, phi_P, p.eta(), p.phi()) < 0.4) {
 63            Etsum += p.Et();
 64          }
 65        }
 66        if (Etsum < 2.5*GeV) {
 67          isolated_photons.push_back(photon);
 68        }
 69      }
 70
 71      if (isolated_photons.size() != 2) {
 72        vetoEvent;
 73      }
 74      std::sort(isolated_photons.begin(), isolated_photons.end(), cmpMomByPt);
 75
 76      FourMomentum y1=isolated_photons[0].momentum();
 77      FourMomentum y2=isolated_photons[1].momentum();
 78      if (deltaR(y1, y2)<0.4) {
 79        vetoEvent;
 80      }
 81
 82      FourMomentum yy=y1+y2;
 83      double Myy = yy.mass()/GeV;
 84      if (Myy<30.0 || Myy>350.0) {
 85        vetoEvent;
 86      }
 87
 88      double pTyy = yy.pT()/GeV;
 89      if (Myy<pTyy) {
 90        vetoEvent;
 91      }
 92
 93      double dPhiyy = mapAngle0ToPi(y1.phi()-y2.phi());
 94      if (dPhiyy<0.5*M_PI) {
 95        vetoEvent;
 96      }
 97
 98      double costhetayy = fabs(tanh((y1.eta()-y2.eta())/2.0));
 99
100      _h_M->fill(Myy, weight);
101      _h_pT->fill(pTyy, weight);
102      _h_dPhi->fill(dPhiyy, weight);
103      _h_costheta->fill(costhetayy, weight);
104
105      _h_pT_M.fill(Myy, pTyy, weight);
106      _h_dPhi_M.fill(Myy, dPhiyy, weight);
107      _h_costheta_M.fill(Myy, costhetayy, weight);
108    }
109
110
111    void finalize() {
112
113      scale(_h_M, crossSection()/sumOfWeights());
114      scale(_h_pT, crossSection()/sumOfWeights());
115      scale(_h_dPhi, crossSection()/sumOfWeights());
116      scale(_h_costheta, crossSection()/sumOfWeights());
117
118      _h_pT_M.scale(crossSection()/sumOfWeights(), this);
119      _h_dPhi_M.scale(crossSection()/sumOfWeights(), this);
120      _h_costheta_M.scale(crossSection()/sumOfWeights(), this);
121
122    }
123
124
125  private:
126
127    Histo1DPtr _h_M;
128    Histo1DPtr _h_pT;
129    Histo1DPtr _h_dPhi;
130    Histo1DPtr _h_costheta;
131    BinnedHistogram _h_pT_M;
132    BinnedHistogram _h_dPhi_M;
133    BinnedHistogram _h_costheta_M;
134
135  };
136
137
138
139  RIVET_DECLARE_ALIASED_PLUGIN(D0_2010_S8570965, D0_2010_I846997);
140
141}