Rivet analyses referenceH1_2015_I1343110Diffractive dijets in DIS and photoproductionExperiment: H1 (HERA) Inspire ID: 1343110 Status: UNVALIDATED Authors:
Beam energies: (920.0, 27.5) GeV Run details:
H1 diffractive jets from proton--positron collisions at beam energies of 920 GeV and 27.5 GeV. The cross section of the diffractive process e+p -> e+X+p is measured at a centre-of-mass energy of 318 GeV, where the system X contains at least two jets and the leading final state proton p is detected in the H1 Very Forward Proton Spectrometer. The measurement is performed in photoproduction with photon virtualities Q^2 < 2 GeV^2 and in deep-inelastic scattering with 4 GeV^2 < Q2 < 80 GeV^2. Source code: H1_2015_I1343110.cc 1// -*- C++ -*-
2#include "Rivet/Analysis.hh"
3#include "Rivet/Projections/Beam.hh"
4#include "Rivet/Projections/FinalState.hh"
5#include "Rivet/Projections/DISKinematics.hh"
6#include "Rivet/Projections/DISFinalState.hh"
7#include "Rivet/Projections/DISDiffHadron.hh"
8#include "Rivet/Projections/FastJets.hh"
9
10namespace Rivet {
11
12 /// @brief H1 diffractive dijets
13 ///
14 /// Diffractive dijets H1 with 920 GeV p and 27.5 GeV e
15 /// Tagged protons & jets found in gamma*p rest frame.
16 ///
17 /// @author Christine O. Rasmussen
18 class H1_2015_I1343110 : public Analysis {
19 public:
20
21 /// Constructor
22 RIVET_DEFAULT_ANALYSIS_CTOR(H1_2015_I1343110);
23
24 /// @name Analysis methods
25 /// @{
26
27 // Book projections and histograms
28 void init() {
29
30 declare(DISKinematics(), "Kinematics");
31 declare(DISFinalState(DISFrame::LAB), "FSLab");
32 declare(FastJets(DISFinalState(DISFrame::HCM),
33 fastjet::JetAlgorithm::kt_algorithm, fastjet::RecombinationScheme::pt_scheme, 1.0,
34 JetMuons::ALL, JetInvisibles::NONE, nullptr), "DISFSJets");
35 declare(DISDiffHadron(), "Hadron");
36
37 // Book histograms from REF data
38 book(_h_PHO_sig_sqrts, 1, 1, 1);
39 book(_h_DIS_sig_sqrts, 2, 1, 1);
40 book(_h_PHODIS_sqrts, 3, 1, 1);
41
42 book(_h_DIS_dsigdz, 4, 1, 1);
43 book(_h_DIS_dsigdxPom, 5, 1, 1);
44 book(_h_DIS_dsigdy, 6, 1, 1);
45 book(_h_DIS_dsigdQ2, 7, 1, 1);
46 book(_h_DIS_dsigdEtj1, 8, 1, 1);
47 book(_h_DIS_dsigdMX, 9, 1, 1);
48 book(_h_DIS_dsigdDeltaEta, 10, 1, 1);
49 book(_h_DIS_dsigdAvgEta, 11, 1, 1);
50
51 book(_h_PHO_dsigdz, 12, 1, 1);
52 book(_h_PHO_dsigdxPom, 13, 1, 1);
53 book(_h_PHO_dsigdy, 14, 1, 1);
54 book(_h_PHO_dsigdxGam, 15, 1, 1);
55 book(_h_PHO_dsigdEtj1, 16, 1, 1);
56 book(_h_PHO_dsigdMX, 17, 1, 1);
57 book(_h_PHO_dsigdDeltaEta, 18, 1, 1);
58 book(_h_PHO_dsigdAvgEta, 19, 1, 1);
59
60 book(_h_PHODIS_deltaEta, 20, 1, 1);
61 book(_h_PHODIS_y, 21, 1, 1);
62 book(_h_PHODIS_z, 22, 1, 1);
63 book(_h_PHODIS_Etj1, 23, 1, 1);
64
65 isPHO = false;
66 nVeto0 = 0;
67 nVeto1 = 0;
68 nVeto2 = 0;
69 nVeto3 = 0;
70 nVeto4 = 0;
71 nVeto5 = 0;
72 nVeto6 = 0;
73 nPHO = 0;
74 nDIS = 0;
75 }
76
77 // Do the analysis
78 void analyze(const Event& event) {
79
80 // Event weight
81 isPHO = false;
82
83 // Projections - special handling of events where no proton found:
84 const DISKinematics& kin = apply<DISKinematics>(event, "Kinematics");
85 if (kin.failed()) vetoEvent;
86 Particle hadronOut;
87 Particle hadronIn;
88 try {
89 const DISDiffHadron& diffhadr = apply<DISDiffHadron>(event, "Hadron");
90 hadronOut = diffhadr.out();
91 hadronIn = diffhadr.in();
92 } catch (const Error& e){
93 vetoEvent;
94 }
95 ++nVeto0;
96
97 // Determine kinematics: H1 has +z = proton direction
98 int dir = kin.orientation();
99 double y = kin.y();
100 double Q2 = kin.Q2();
101
102 // Separate into DIS and PHO regimes else veto
103 if (Q2 < 2.*GeV2 && inRange(y, 0.2, 0.70)) {
104 isPHO = true;
105 ++nPHO;
106 } else if (inRange(Q2, 4.0*GeV2, 80.*GeV2) && inRange(y, 0.2, 0.7)) {
107 isPHO = false;
108 ++nDIS;
109 } else vetoEvent;
110 ++nVeto1;
111
112 // Find diffractive variables as defined in paper.
113 // Note tagged protons in VFPS => smaller allowed xPom range
114 // xPom = 1 - E'/E, t = (P-P')^2
115 const double abst = abs((hadronIn.mom() - hadronOut.mom()).mass2());
116 const double xPom = 1. - hadronOut.energy() / hadronIn.energy();
117 // Veto if outside allowed region
118 if (abst > 0.6 * GeV2) vetoEvent;
119 ++nVeto2;
120 if (!inRange(xPom, 0.010, 0.024)) vetoEvent;
121 ++nVeto3;
122
123 // Jet selection. Note jets are found in HCM frame, but
124 // eta cut is applied in lab frame!
125 Jets jets = apply<FastJets>(event, "DISFSJets").jets(Cuts::Et > 4.* GeV, cmpMomByEt);
126 // Veto if not dijets and if Et_j1 < 5.5
127 if (jets.size() < 2) vetoEvent;
128 if (jets[0].Et() < 5.5 * GeV) vetoEvent;
129 ++nVeto4;
130
131 const double EtJet1 = jets[0].Et();
132
133 // Transform from HCM to LAB and apply eta cut
134 const LorentzTransform boost = kin.boostHCM();
135 for (int i = 0; i < 2; ++i) {
136 jets[i].transformBy(boost.inverse());
137 if (!inRange(dir * jets[i].eta(),-1,2.5)) vetoEvent;
138 }
139 ++nVeto5;
140
141 double sumE(0.), sumPz(0.);
142 FourMomentum tot_mom;
143 for(const Particle & part : apply<FinalState>(event,"FSLab").particles()) {
144 if (part.genParticle() != hadronOut.genParticle()) {
145 sumE += part.energy();
146 sumPz += part.pz();
147 tot_mom += part.mom();
148 }
149 }
150 double EplusPz = sumE + dir * sumPz;
151 double EminusPz = sumE - dir * sumPz;
152 double M2X = tot_mom.mass2();
153
154 // Pseudorapidity distributions are examined in lab frame:
155 double deltaEtaJets = abs(dir * jets[0].eta() - dir * jets[1].eta());
156 double avgEtaJets = 0.5 * (dir * jets[0].eta() + dir * jets[1].eta());
157
158 // Evaluate observables
159 double zPomJets(0.), xGamJets(0.);
160 FourMomentum pj = jets[0].momentum()+jets[1].momentum();
161 if (isPHO){
162 zPomJets = (pj.E() + dir * pj.pz()) / EplusPz;
163 xGamJets = (pj.E() - dir * pj.pz()) / EminusPz;
164 } else {
165 zPomJets = (Q2 + pj.mass2()) / (Q2 + M2X);
166 }
167 // Veto events with zPom > 0.8
168 if (zPomJets > 0.8) vetoEvent;
169 ++nVeto6;
170
171 // Now fill histograms
172 if (isPHO){
173 _h_PHO_sig_sqrts ->fill(sqrtS()/GeV);
174 _h_PHO_dsigdz ->fill(zPomJets);
175 _h_PHO_dsigdxPom ->fill(xPom);
176 _h_PHO_dsigdy ->fill(y);
177 _h_PHO_dsigdxGam ->fill(xGamJets);
178 _h_PHO_dsigdEtj1 ->fill(EtJet1);
179 _h_PHO_dsigdMX ->fill(sqrt(M2X)/GeV);
180 _h_PHO_dsigdDeltaEta ->fill(deltaEtaJets);
181 _h_PHO_dsigdAvgEta ->fill(avgEtaJets);
182 } else {
183 _h_DIS_sig_sqrts ->fill(sqrtS()/GeV);
184 _h_DIS_dsigdz ->fill(zPomJets);
185 _h_DIS_dsigdxPom ->fill(xPom);
186 _h_DIS_dsigdy ->fill(y);
187 _h_DIS_dsigdQ2 ->fill(Q2);
188 _h_DIS_dsigdEtj1 ->fill(EtJet1);
189 _h_DIS_dsigdMX ->fill(sqrt(M2X)/GeV);
190 _h_DIS_dsigdDeltaEta ->fill(deltaEtaJets);
191 _h_DIS_dsigdAvgEta ->fill(avgEtaJets);
192 }
193
194 }
195
196 // Finalize
197 void finalize() {
198 // Normalise to cross section
199 // Remember to manually scale the cross section afterwards with
200 // the number of rejected events.
201 const double norm = crossSection()/picobarn/sumOfWeights();
202
203 scale(_h_PHO_sig_sqrts, norm);
204 scale(_h_PHO_dsigdz, norm);
205 scale(_h_PHO_dsigdxPom, norm);
206 scale(_h_PHO_dsigdy, norm);
207 scale(_h_PHO_dsigdxGam, norm);
208 scale(_h_PHO_dsigdEtj1, norm);
209 scale(_h_PHO_dsigdMX, norm);
210 scale(_h_PHO_dsigdDeltaEta, norm);
211 scale(_h_PHO_dsigdAvgEta, norm);
212
213 scale(_h_DIS_sig_sqrts, norm);
214 scale(_h_DIS_dsigdz, norm);
215 scale(_h_DIS_dsigdxPom, norm);
216 scale(_h_DIS_dsigdy, norm);
217 scale(_h_DIS_dsigdQ2, norm);
218 scale(_h_DIS_dsigdEtj1, norm);
219 scale(_h_DIS_dsigdMX, norm);
220 scale(_h_DIS_dsigdDeltaEta, norm);
221 scale(_h_DIS_dsigdAvgEta, norm);
222
223 if (_h_DIS_sig_sqrts->numEntries() != 0)
224 divide(_h_PHO_sig_sqrts, _h_DIS_sig_sqrts, _h_PHODIS_sqrts);
225 if (_h_DIS_dsigdDeltaEta->numEntries() != 0)
226 divide(_h_PHO_dsigdDeltaEta, _h_DIS_dsigdDeltaEta, _h_PHODIS_deltaEta);
227 if (_h_DIS_dsigdy->numEntries() != 0)
228 divide(_h_PHO_dsigdy, _h_DIS_dsigdy, _h_PHODIS_y);
229 if (_h_DIS_dsigdz->numEntries() != 0)
230 divide(_h_PHO_dsigdz, _h_DIS_dsigdz, _h_PHODIS_z);
231 if (_h_DIS_dsigdEtj1->numEntries() != 0)
232 divide(_h_PHO_dsigdEtj1, _h_DIS_dsigdEtj1, _h_PHODIS_Etj1);
233
234 const double dPHO = numEvents();
235 MSG_INFO("H1_2015_I1343110");
236 MSG_INFO("Cross section = " << crossSection()/picobarn << " pb");
237 MSG_INFO("Number of events = " << numEvents() << ", sumW = " << sumOfWeights());
238 MSG_INFO("Number of PHO = " << nPHO << ", number of DIS = " << nDIS);
239 MSG_INFO("Events passing event setup = " << nVeto0 << " (" << nVeto0/dPHO * 100. << "%)" );
240 MSG_INFO("Events passing electron veto = " << nVeto1 << " (" << nVeto1/dPHO * 100. << "%)" );
241 MSG_INFO("Events passing t veto = " << nVeto2 << " (" << nVeto2/dPHO * 100. << "%)" );
242 MSG_INFO("Events passing xPom = " << nVeto3 << " (" << nVeto3/dPHO * 100. << "%)" );
243 MSG_INFO("Events passing jet Et veto = " << nVeto4 << " (" << nVeto4/dPHO * 100. << "%)" );
244 MSG_INFO("Events passing jet eta veto = " << nVeto5 << " (" << nVeto5/dPHO * 100. << "%)" );
245 MSG_INFO("Events passing zPom veto = " << nVeto6 << " (" << nVeto6/dPHO * 100. << "%)" );
246
247 }
248
249 /// @}
250
251
252 private:
253
254 /// @name Histograms
255 /// @{
256 // Book histograms from REF data
257 BinnedHistoPtr<int> _h_PHO_sig_sqrts;
258 BinnedHistoPtr<int> _h_DIS_sig_sqrts;
259 BinnedEstimatePtr<int> _h_PHODIS_sqrts;
260
261 Histo1DPtr _h_DIS_dsigdz;
262 Histo1DPtr _h_DIS_dsigdxPom;
263 Histo1DPtr _h_DIS_dsigdy;
264 Histo1DPtr _h_DIS_dsigdQ2;
265 Histo1DPtr _h_DIS_dsigdEtj1;
266 Histo1DPtr _h_DIS_dsigdMX;
267 Histo1DPtr _h_DIS_dsigdDeltaEta;
268 Histo1DPtr _h_DIS_dsigdAvgEta;
269
270 Histo1DPtr _h_PHO_dsigdz;
271 Histo1DPtr _h_PHO_dsigdxPom;
272 Histo1DPtr _h_PHO_dsigdy;
273 Histo1DPtr _h_PHO_dsigdxGam;
274 Histo1DPtr _h_PHO_dsigdEtj1;
275 Histo1DPtr _h_PHO_dsigdMX;
276 Histo1DPtr _h_PHO_dsigdDeltaEta;
277 Histo1DPtr _h_PHO_dsigdAvgEta;
278
279 Estimate1DPtr _h_PHODIS_deltaEta;
280 Estimate1DPtr _h_PHODIS_y;
281 Estimate1DPtr _h_PHODIS_z;
282 Estimate1DPtr _h_PHODIS_Etj1;
283 /// @}
284
285 bool isPHO;
286 int nVeto0, nVeto1, nVeto2, nVeto3, nVeto4, nVeto5, nVeto6;
287 int nPHO, nDIS;
288 };
289
290 RIVET_DECLARE_PLUGIN(H1_2015_I1343110);
291
292}
|
