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
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189

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

namespace Rivet {


  /// @brief Measurement of isolated gamma + jet + X differential cross-sections
  ///
  /// Inclusive isolated gamma + jet cross-sections, differential in pT(gamma), for
  /// various photon and jet rapidity bins.
  ///
  /// @author Andy Buckley
  /// @author Gavin Hesketh
  class D0_2008_S7719523 : public Analysis {
  public:

    RIVET_DEFAULT_ANALYSIS_CTOR(D0_2008_S7719523);


    /// @name Analysis methods
    /// @{

    /// Set up projections and book histograms
    void init() {
      // General FS
      FinalState fs;
      declare(fs, "FS");

      // Get leading photon
      LeadingParticlesFinalState photonfs(FinalState((Cuts::etaIn(-1.0, 1.0) && Cuts::pT >=  30.0*GeV)));
      photonfs.addParticleId(PID::PHOTON);
      declare(photonfs, "LeadingPhoton");

      // FS excluding the leading photon
      VetoedFinalState vfs(fs);
      vfs.addVetoOnThisFinalState(photonfs);
      declare(vfs, "JetFS");

      // Jets
      FastJets jetpro(vfs, FastJets::D0ILCONE, 0.7);
      declare(jetpro, "Jets");

      // Histograms
      book(_h_central_same_cross_section ,1, 1, 1);
      book(_h_central_opp_cross_section  ,2, 1, 1);
      book(_h_forward_same_cross_section ,3, 1, 1);
      book(_h_forward_opp_cross_section  ,4, 1, 1);

      // Ratio histos to be filled by divide()
      book(_h_cen_opp_same, 5, 1, 1);
      book(_h_fwd_opp_same, 8, 1, 1);
      // Ratio histos to be filled manually, since the num/denom inputs don't match
      book(_h_cen_same_fwd_same, 6, 1, 1, true);
      book(_h_cen_opp_fwd_same, 7, 1, 1, true);
      book(_h_cen_same_fwd_opp, 9, 1, 1, true);
      book(_h_cen_opp_fwd_opp, 10, 1, 1, true);
    }



    /// Do the analysis
    void analyze(const Event& event) {
      // Get the photon
      const FinalState& photonfs = apply<FinalState>(event, "LeadingPhoton");
      if (photonfs.particles().size() != 1) {
        vetoEvent;
      }
      const FourMomentum photon = photonfs.particles().front().momentum();

      // Isolate photon by ensuring that a 0.4 cone around it contains less than 7% of the photon's energy
      double egamma = photon.E();
      double eta_P = photon.eta();
      double phi_P = photon.phi();
      double econe = 0.0;
      for (const Particle& p : apply<FinalState>(event, "JetFS").particles()) {
        if (deltaR(eta_P, phi_P, p.eta(), p.phi()) < 0.4) {
          econe += p.E();
          // Veto as soon as E_cone gets larger
          if (econe/egamma > 0.07) {
            MSG_DEBUG("Vetoing event because photon is insufficiently isolated");
            vetoEvent;
          }
        }
      }

      Jets jets = apply<FastJets>(event, "Jets").jetsByPt(15.0*GeV);
      if (jets.empty()) vetoEvent;

      FourMomentum leadingJet = jets[0].momentum();
      if (deltaR(eta_P, phi_P, leadingJet.eta(), leadingJet.phi()) < 0.7) {
        vetoEvent;
      }

      int photon_jet_sign = sign( leadingJet.rapidity() * photon.rapidity() );

      // Veto if leading jet is outside plotted rapidity regions
      const double abs_y1 = fabs(leadingJet.rapidity());
      if (inRange(abs_y1, 0.8, 1.5) || abs_y1 > 2.5) {
        MSG_DEBUG("Leading jet falls outside acceptance range; |y1| = " << abs_y1);
        vetoEvent;
      }

      // Fill histos
      if (fabs(leadingJet.rapidity()) < 0.8) {
        Histo1DPtr h = (photon_jet_sign >= 1) ? _h_central_same_cross_section : _h_central_opp_cross_section;
        h->fill(photon.pT());
      } else if (inRange( fabs(leadingJet.rapidity()), 1.5, 2.5)) {
        Histo1DPtr h = (photon_jet_sign >= 1) ? _h_forward_same_cross_section : _h_forward_opp_cross_section;
        h->fill(photon.pT());
      }

    }


    /// Finalize
    void finalize() {
      const double lumi_gen = sumOfWeights()/crossSection();
      const double dy_photon = 2.0;
      const double dy_jet_central = 1.6;
      const double dy_jet_forward = 2.0;

      // Cross-section ratios (6 plots)
      // Central/central and forward/forward ratios
      divide(_h_central_opp_cross_section, _h_central_same_cross_section, _h_cen_opp_same);
      divide(_h_forward_opp_cross_section, _h_forward_same_cross_section, _h_fwd_opp_same);
      // Central/forward ratio combinations
      /// @note The central/forward histo binnings are not the same! Hence the need to do these by hand :-(
      for (size_t i = 0; i < _h_cen_same_fwd_same->numPoints(); ++i) {
        const YODA::HistoBin1D& cen_same_bini = _h_central_same_cross_section->bin(i);
        const YODA::HistoBin1D& cen_opp_bini = _h_central_opp_cross_section->bin(i);
        const YODA::HistoBin1D& fwd_same_bini = _h_central_same_cross_section->bin(i);
        const YODA::HistoBin1D& fwd_opp_bini = _h_central_opp_cross_section->bin(i);
        _h_cen_same_fwd_same->point(i).setY(_safediv(cen_same_bini.sumW(), fwd_same_bini.sumW(), 0),
                                            add_quad(cen_same_bini.relErr(), fwd_same_bini.relErr()));
        _h_cen_opp_fwd_same->point(i).setY(_safediv(cen_opp_bini.sumW(), fwd_same_bini.sumW(), 0),
                                           add_quad(cen_opp_bini.relErr(), fwd_same_bini.relErr()));
        _h_cen_same_fwd_opp->point(i).setY(_safediv(cen_same_bini.sumW(), fwd_opp_bini.sumW(), 0),
                                           add_quad(cen_same_bini.relErr(), fwd_opp_bini.relErr()));
        _h_cen_opp_fwd_opp->point(i).setY(_safediv(cen_opp_bini.sumW(), fwd_opp_bini.sumW(), 0),
                                          add_quad(cen_opp_bini.relErr(), fwd_opp_bini.relErr()));
      }

      // Use generator cross section for remaining histograms
      // Each of these needs the additional factor 2 because the
      // y_photon * y_jet requirement reduces the corresponding 2D "bin width"
      // by a factor 1/2.
      scale(_h_central_same_cross_section, 2.0/lumi_gen * 1.0/dy_photon * 1.0/dy_jet_central);
      scale(_h_central_opp_cross_section, 2.0/lumi_gen * 1.0/dy_photon * 1.0/dy_jet_central);
      scale(_h_forward_same_cross_section, 2.0/lumi_gen * 1.0/dy_photon * 1.0/dy_jet_forward);
      scale(_h_forward_opp_cross_section, 2.0/lumi_gen * 1.0/dy_photon * 1.0/dy_jet_forward);
    }

    /// @}


  private:

    // A local scope function for division, handling the div-by-zero case
    /// @todo Why isn't the math divide() function being found?
    double _safediv(double a, double b, double result_if_err) {
      return (b != 0) ? a/b : result_if_err;
    }

    /// @name Histograms
    /// @{
    Histo1DPtr _h_central_same_cross_section;
    Histo1DPtr _h_central_opp_cross_section;
    Histo1DPtr _h_forward_same_cross_section;
    Histo1DPtr _h_forward_opp_cross_section;

    Scatter2DPtr _h_cen_opp_same;
    Scatter2DPtr _h_fwd_opp_same;
    Scatter2DPtr _h_cen_same_fwd_same;
    Scatter2DPtr _h_cen_opp_fwd_same;
    Scatter2DPtr _h_cen_same_fwd_opp;
    Scatter2DPtr _h_cen_opp_fwd_opp;
    /// @}

  };



  RIVET_DECLARE_ALIASED_PLUGIN(D0_2008_S7719523, D0_2008_I782968);

}