Rivet analyses referenceD0_2008_I782968Isolated γ + jet cross-sections, differential in pT(γ) for various y binsExperiment: D0 (Tevatron Run 2) Inspire ID: 782968 Status: VALIDATED Authors:
Beam energies: (980.0, 980.0) GeV Run details:
The process pˉp -> photon + jet + X as studied by the D0 detector at the Fermilab Tevatron collider at center-of-mass energy √s=1.96TeV. Photons are reconstructed in the central rapidity region |yγ|<1.0 with transverse momenta in the range 30--400 GeV, while jets are reconstructed in either the central |yjet|<0.8 or forward 1.5<|yjet|<2.5 rapidity intervals with pjet⊥>15GeV. The differential cross section d3σ/dpγ⊥dyγdyjet is measured as a function of pγ⊥ in four regions, differing by the relative orientations of the photon and the jet. MC predictions have trouble with simultaneously describing the measured normalization and pγ⊥ dependence of the cross section in any of the four measured regions. Source code: D0_2008_I782968.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/FinalState.hh"
4#include "Rivet/Projections/LeadingParticlesFinalState.hh"
5#include "Rivet/Projections/VetoedFinalState.hh"
6#include "Rivet/Projections/FastJets.hh"
7
8namespace Rivet {
9
10
11 /// @brief Measurement of isolated gamma + jet + X differential cross-sections
12 ///
13 /// Inclusive isolated gamma + jet cross-sections, differential in pT(gamma), for
14 /// various photon and jet rapidity bins.
15 ///
16 /// @author Andy Buckley
17 /// @author Gavin Hesketh
18 class D0_2008_I782968 : public Analysis {
19 public:
20
21 RIVET_DEFAULT_ANALYSIS_CTOR(D0_2008_I782968);
22
23
24 /// @name Analysis methods
25 /// @{
26
27 /// Set up projections and book histograms
28 void init() {
29 // General FS
30 FinalState fs;
31 declare(fs, "FS");
32
33 // Get leading photon
34 LeadingParticlesFinalState photonfs(FinalState((Cuts::etaIn(-1.0, 1.0) && Cuts::pT >= 30.0*GeV)));
35 photonfs.addParticleId(PID::PHOTON);
36 declare(photonfs, "LeadingPhoton");
37
38 // FS excluding the leading photon
39 VetoedFinalState vfs(fs);
40 vfs.addVetoOnThisFinalState(photonfs);
41 declare(vfs, "JetFS");
42
43 // Jets
44 FastJets jetpro(vfs, JetAlg::D0ILCONE, 0.7);
45 declare(jetpro, "Jets");
46
47 // Histograms
48 book(_h_central_same_cross_section ,1, 1, 1);
49 book(_h_central_opp_cross_section ,2, 1, 1);
50 book(_h_forward_same_cross_section ,3, 1, 1);
51 book(_h_forward_opp_cross_section ,4, 1, 1);
52
53 // Ratio histos to be filled by divide()
54 book(_h_cen_opp_same, 5, 1, 1);
55 book(_h_fwd_opp_same, 8, 1, 1);
56 // Ratio histos to be filled manually, since the num/denom inputs don't match
57 book(_h_cen_same_fwd_same, 6, 1, 1);
58 book(_h_cen_opp_fwd_same, 7, 1, 1);
59 book(_h_cen_same_fwd_opp, 9, 1, 1);
60 book(_h_cen_opp_fwd_opp, 10, 1, 1);
61 }
62
63
64
65 /// Do the analysis
66 void analyze(const Event& event) {
67 // Get the photon
68 const FinalState& photonfs = apply<FinalState>(event, "LeadingPhoton");
69 if (photonfs.particles().size() != 1) {
70 vetoEvent;
71 }
72 const FourMomentum photon = photonfs.particles().front().momentum();
73
74 // Isolate photon by ensuring that a 0.4 cone around it contains less than 7% of the photon's energy
75 double egamma = photon.E();
76 double eta_P = photon.eta();
77 double phi_P = photon.phi();
78 double econe = 0.0;
79 for (const Particle& p : apply<FinalState>(event, "JetFS").particles()) {
80 if (deltaR(eta_P, phi_P, p.eta(), p.phi()) < 0.4) {
81 econe += p.E();
82 // Veto as soon as E_cone gets larger
83 if (econe/egamma > 0.07) {
84 MSG_DEBUG("Vetoing event because photon is insufficiently isolated");
85 vetoEvent;
86 }
87 }
88 }
89
90 Jets jets = apply<FastJets>(event, "Jets").jetsByPt(Cuts::pT > 15.0*GeV);
91 if (jets.empty()) vetoEvent;
92
93 FourMomentum leadingJet = jets[0].momentum();
94 if (deltaR(eta_P, phi_P, leadingJet.eta(), leadingJet.phi()) < 0.7) {
95 vetoEvent;
96 }
97
98 int photon_jet_sign = sign( leadingJet.rapidity() * photon.rapidity() );
99
100 // Veto if leading jet is outside plotted rapidity regions
101 const double abs_y1 = fabs(leadingJet.rapidity());
102 if (inRange(abs_y1, 0.8, 1.5) || abs_y1 > 2.5) {
103 MSG_DEBUG("Leading jet falls outside acceptance range; |y1| = " << abs_y1);
104 vetoEvent;
105 }
106
107 // Fill histos
108 if (fabs(leadingJet.rapidity()) < 0.8) {
109 Histo1DPtr h = (photon_jet_sign >= 1) ? _h_central_same_cross_section : _h_central_opp_cross_section;
110 h->fill(photon.pT());
111 } else if (inRange( fabs(leadingJet.rapidity()), 1.5, 2.5)) {
112 Histo1DPtr h = (photon_jet_sign >= 1) ? _h_forward_same_cross_section : _h_forward_opp_cross_section;
113 h->fill(photon.pT());
114 }
115
116 }
117
118
119 /// Finalize
120 void finalize() {
121 const double lumi_gen = sumOfWeights()/crossSection()/picobarn;
122 const double dy_photon = 2.0;
123 const double dy_jet_central = 1.6;
124 const double dy_jet_forward = 2.0;
125
126 // Cross-section ratios (6 plots)
127 // Central/central and forward/forward ratios
128 divide(_h_central_opp_cross_section, _h_central_same_cross_section, _h_cen_opp_same);
129 divide(_h_forward_opp_cross_section, _h_forward_same_cross_section, _h_fwd_opp_same);
130 // Central/forward ratio combinations
131 /// @note The central/forward histo binnings are not the same! Hence the need to do these by hand :-(
132 for (size_t i = 1; i < _h_cen_same_fwd_same->numBins()+1; ++i) {
133 const auto& cen_same_bini = _h_central_same_cross_section->bin(i);
134 const auto& cen_opp_bini = _h_central_opp_cross_section->bin(i);
135 const auto& fwd_same_bini = _h_central_same_cross_section->bin(i);
136 const auto& fwd_opp_bini = _h_central_opp_cross_section->bin(i);
137 _h_cen_same_fwd_same->bin(i).set(_safediv(cen_same_bini.sumW(), fwd_same_bini.sumW(), 0),
138 add_quad(cen_same_bini.relErrW(), fwd_same_bini.relErrW()));
139 _h_cen_opp_fwd_same->bin(i).set(_safediv(cen_opp_bini.sumW(), fwd_same_bini.sumW(), 0),
140 add_quad(cen_opp_bini.relErrW(), fwd_same_bini.relErrW()));
141 _h_cen_same_fwd_opp->bin(i).set(_safediv(cen_same_bini.sumW(), fwd_opp_bini.sumW(), 0),
142 add_quad(cen_same_bini.relErrW(), fwd_opp_bini.relErrW()));
143 _h_cen_opp_fwd_opp->bin(i).set(_safediv(cen_opp_bini.sumW(), fwd_opp_bini.sumW(), 0),
144 add_quad(cen_opp_bini.relErrW(), fwd_opp_bini.relErrW()));
145 }
146
147 // Use generator cross section for remaining histograms
148 // Each of these needs the additional factor 2 because the
149 // y_photon * y_jet requirement reduces the corresponding 2D "bin width"
150 // by a factor 1/2.
151 scale(_h_central_same_cross_section, 2.0/lumi_gen * 1.0/dy_photon * 1.0/dy_jet_central);
152 scale(_h_central_opp_cross_section, 2.0/lumi_gen * 1.0/dy_photon * 1.0/dy_jet_central);
153 scale(_h_forward_same_cross_section, 2.0/lumi_gen * 1.0/dy_photon * 1.0/dy_jet_forward);
154 scale(_h_forward_opp_cross_section, 2.0/lumi_gen * 1.0/dy_photon * 1.0/dy_jet_forward);
155 }
156
157 /// @}
158
159
160 private:
161
162 // A local scope function for division, handling the div-by-zero case
163 /// @todo Why isn't the math divide() function being found?
164 double _safediv(double a, double b, double result_if_err) {
165 return (b != 0) ? a/b : result_if_err;
166 }
167
168 /// @name Histograms
169 /// @{
170 Histo1DPtr _h_central_same_cross_section;
171 Histo1DPtr _h_central_opp_cross_section;
172 Histo1DPtr _h_forward_same_cross_section;
173 Histo1DPtr _h_forward_opp_cross_section;
174
175 Estimate1DPtr _h_cen_opp_same;
176 Estimate1DPtr _h_fwd_opp_same;
177 Estimate1DPtr _h_cen_same_fwd_same;
178 Estimate1DPtr _h_cen_opp_fwd_same;
179 Estimate1DPtr _h_cen_same_fwd_opp;
180 Estimate1DPtr _h_cen_opp_fwd_opp;
181 /// @}
182
183 };
184
185
186
187 RIVET_DECLARE_ALIASED_PLUGIN(D0_2008_I782968, D0_2008_S7719523);
188
189}
|
