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
13 namespace H1_2015_I1343110_PROJECTIONS {
14
15 /// Projection to find the largest gaps and the masses of the two
16 /// systems separated by the gap. Based on the HZTools gap-finding
17 /// method (hzhadgap.F). Note that gaps are found in the HCM frame.
18 ///
19 /// @author Christine O. Rasmussen.
20 class RapidityGap : public Projection {
21 public:
22
23 /// Type of DIS boost to apply
24 enum Frame { HCM, LAB, XCM };
25
26 RapidityGap() {
27 setName("RapidityGap");
28 declare(DISKinematics(), "DISKIN");
29 declare(DISFinalState(DISFrame::HCM), "DISFS");
30 }
31
32 RIVET_DEFAULT_PROJ_CLONE(RapidityGap);
33
34 /// Import to avoid warnings about overload-hiding
35 using Projection::operator =;
36
37 double M2X() const { return _M2X; }
38 double M2Y() const { return _M2Y; }
39 double t() const { return _t; }
40 double gap() const { return _gap; }
41 double gapUpp() const { return _gapUpp; }
42 double gapLow() const { return _gapLow; }
43 double EpPzX(Frame f) const {
44 if (f == LAB) return _ePpzX_LAB;
45 else if (f == XCM) return _ePpzX_XCM;
46 else return _ePpzX_HCM;
47 }
48 double EmPzX(Frame f) const {
49 if (f == LAB) return _eMpzX_LAB;
50 else if (f == XCM) return _eMpzX_XCM;
51 else return _eMpzX_HCM;
52 }
53 FourMomentum pX(Frame f) const {
54 if (f == LAB) return _momX_LAB;
55 else if (f == XCM) return _momX_XCM;
56 else return _momX_HCM;
57 }
58 FourMomentum pY(Frame f) const {
59 if (f == LAB) return _momY_LAB;
60 else if (f == XCM) return _momY_XCM;
61 else return _momY_HCM;
62 }
63 const Particles& systemX(Frame f) const {
64 if (f == LAB) return _pX_LAB;
65 else if (f == XCM) return _pX_XCM;
66 else return _pX_HCM;
67 }
68 const Particles& systemY(Frame f) const {
69 if (f == LAB) return _pY_LAB;
70 else if (f == XCM) return _pY_XCM;
71 else return _pY_HCM;
72 }
73
74 protected:
75
76 virtual CmpState compare(const Projection& p) const {
77 const RapidityGap& other = pcast<RapidityGap>(p);
78 return mkNamedPCmp(other, "DISKIN") || mkNamedPCmp(other, "DISFS");
79 }
80
81 virtual void project(const Event& e){
82 const DISKinematics& dk = apply<DISKinematics>(e, "DISKIN");
83 const Particles& p = apply<DISFinalState>(e, "DISFS").particles(cmpMomByEta);
84 findgap(p, dk);
85 }
86
87 void clearAll(){
88 _M2X = _M2Y = _t = _gap = 0.;
89 _gapUpp = _gapLow = -8.;
90 _ePpzX_HCM = _eMpzX_HCM =_ePpzX_LAB = _eMpzX_LAB = _ePpzX_XCM = _eMpzX_XCM = 0.;
91 _momX_HCM.setPE(0., 0., 0., 0.);
92 _momY_HCM.setPE(0., 0., 0., 0.);
93 _momX_XCM.setPE(0., 0., 0., 0.);
94 _momY_XCM.setPE(0., 0., 0., 0.);
95 _momX_LAB.setPE(0., 0., 0., 0.);
96 _momY_LAB.setPE(0., 0., 0., 0.);
97 _pX_HCM.clear();
98 _pY_HCM.clear();
99 _pX_XCM.clear();
100 _pY_XCM.clear();
101 _pX_LAB.clear();
102 _pY_LAB.clear();
103 }
104
105 void findgap(const Particles& particles, const DISKinematics& diskin){
106
107 clearAll();
108
109 // Begin by finding largest gap and gapedges between all final
110 // state particles in HCM frame.
111 int nP = particles.size();
112 int dir = diskin.orientation();
113 for (int i = 0; i < nP-1; ++i){
114 double tmpGap = abs(particles[i+1].eta() - particles[i].eta());
115 if (tmpGap > _gap) {
116 _gap = tmpGap;
117 _gapLow = (dir > 0) ? particles[i].eta() : dir * particles[i+1].eta();
118 _gapUpp = (dir > 0) ? particles[i+1].eta() : dir * particles[i].eta();
119 }
120 }
121
122 // Define the two systems X and Y.
123 Particles tmp_pX, tmp_pY;
124 for (const Particle& ip : particles) {
125 if (dir * ip.eta() > _gapLow) tmp_pX.push_back(ip);
126 else tmp_pY.push_back(ip);
127 }
128
129 Particles pX, pY;
130 pX = (dir < 0) ? tmp_pY : tmp_pX;
131 pY = (dir < 0) ? tmp_pX : tmp_pY;
132
133 // Find variables related to HCM frame.
134 // Note that HCM has photon along +z, as opposed to
135 // H1 where proton is along +z. This results in a sign change
136 // as compared to H1 papers!
137
138 // X - side
139 FourMomentum momX;
140 for (const Particle& jp : pX) {
141 momX += jp.momentum();
142 _ePpzX_HCM += jp.E() - jp.pz(); // Sign + => -
143 _eMpzX_HCM += jp.E() + jp.pz(); // Sign - => +
144 }
145 _momX_HCM = momX;
146 _pX_HCM = pX;
147 _M2X = _momX_HCM.mass2();
148
149 // Y - side
150 FourMomentum momY;
151 for (const Particle& kp : pY) momY += kp.momentum();
152 _momY_HCM = momY;
153 _pY_HCM = pY;
154 _M2Y = _momY_HCM.mass2();
155
156 // Find variables related to LAB frame
157 const LorentzTransform hcmboost = diskin.boostHCM();
158 const LorentzTransform hcminverse = hcmboost.inverse();
159 _momX_LAB = hcminverse.transform(_momX_HCM);
160 _momY_LAB = hcminverse.transform(_momY_HCM);
161
162 // Find momenta in XCM frame. Note that it is HCM frame that is
163 // boosted, resulting in a sign change later!
164 const bool doXCM = (momX.betaVec().mod2() < 1.);
165 if (doXCM) {
166 const LorentzTransform xcmboost =
167 LorentzTransform::mkFrameTransformFromBeta(momX.betaVec());
168 _momX_XCM = xcmboost.transform(momX);
169 _momY_XCM = xcmboost.transform(momY);
170 }
171
172 for (const Particle& jp : pX) {
173 // Boost from HCM to LAB.
174 FourMomentum lab = hcminverse.transform(jp.momentum());
175 _ePpzX_LAB += lab.E() + dir * lab.pz();
176 _eMpzX_LAB += lab.E() - dir * lab.pz();
177 Particle plab = jp;
178 plab.setMomentum(lab);
179 _pX_LAB.push_back(plab);
180 // Set XCM. Note that since HCM frame is boosted to XCM frame,
181 // we have a sign change
182 if (doXCM) {
183 const LorentzTransform xcmboost =
184 LorentzTransform::mkFrameTransformFromBeta(_momX_HCM.betaVec());
185 FourMomentum xcm = xcmboost.transform(jp.momentum());
186 _ePpzX_XCM += xcm.E() - xcm.pz(); // Sign + => -
187 _eMpzX_XCM += xcm.E() + xcm.pz(); // Sign - => +
188 Particle pxcm = jp;
189 pxcm.setMomentum(xcm);
190 _pX_XCM.push_back(pxcm);
191 }
192 }
193
194 for (const Particle& jp : pY) {
195 // Boost from HCM to LAB
196 FourMomentum lab = hcminverse.transform(jp.momentum());
197 Particle plab = jp;
198 plab.setMomentum(lab);
199 _pY_LAB.push_back(plab);
200 // Boost from HCM to XCM
201 if (doXCM) {
202 const LorentzTransform xcmboost =
203 LorentzTransform::mkFrameTransformFromBeta(_momX_HCM.betaVec());
204 FourMomentum xcm = xcmboost.transform(jp.momentum());
205 Particle pxcm = jp;
206 pxcm.setMomentum(xcm);
207 _pY_XCM.push_back(pxcm);
208 }
209 }
210
211 // Boost p from LAB to HCM frame to find t.
212 ///
213 /// @note Currently can only handle gap on proton side.
214 ///
215 /// @todo Expand to also handle gap on photon side
216 const FourMomentum proton = hcmboost.transform(diskin.beamHadron().momentum());
217 FourMomentum pPom = proton - _momY_HCM;
218 _t = pPom * pPom;
219
220 }
221
222 private:
223
224 double _M2X, _M2Y, _t;
225 double _gap, _gapUpp, _gapLow;
226 double _ePpzX_LAB, _eMpzX_LAB, _ePpzX_HCM, _eMpzX_HCM, _ePpzX_XCM, _eMpzX_XCM;
227 FourMomentum _momX_HCM, _momY_HCM,_momX_LAB, _momY_LAB, _momX_XCM, _momY_XCM;
228 Particles _pX_HCM, _pY_HCM, _pX_LAB, _pY_LAB, _pX_XCM, _pY_XCM;
229
230 };
231
232 /// Projection to boost system X (photon+Pomeron) particles into its rest frame.
233 ///
234 /// @author Ilkka Helenius
235 class BoostedXSystem : public FinalState {
236 public:
237
238 BoostedXSystem(const FinalState& fs) {
239 setName("BoostedXSystem");
240 declare(fs,"FS");
241 declare(RapidityGap(), "RAPGAP");
242 }
243
244 // Return the boost to XCM frame.
245 const LorentzTransform& boost() const { return _boost; }
246
247 RIVET_DEFAULT_PROJ_CLONE(BoostedXSystem);
248
249 /// Import to avoid warnings about overload-hiding
250 using Projection::operator =;
251
252 protected:
253
254 // Apply the projection on the supplied event.
255 void project(const Event& e){
256
257 const RapidityGap& rg = apply<RapidityGap>(e, "RAPGAP");
258
259 // Total momentum of the system X.
260 const FourMomentum pX = rg.pX(RapidityGap::HCM);
261
262 // Reset the boost. Is there a separate method for this?
263 _boost = combine(_boost, _boost.inverse());
264
265 // Define boost only when numerically safe, otherwise negligible.
266 if (pX.betaVec().mod2() < 1.)
267 _boost = LorentzTransform::mkFrameTransformFromBeta(pX.betaVec());
268
269 // Boost the particles from system X.
270 _theParticles.clear();
271 _theParticles.reserve(rg.systemX(RapidityGap::HCM).size());
272 for (const Particle& p : rg.systemX(RapidityGap::HCM)) {
273 Particle temp = p;
274 temp.setMomentum(_boost.transform(temp.momentum()));
275 _theParticles.push_back(temp);
276 }
277
278 }
279
280 // Compare projections.
281 CmpState compare(const Projection& p) const {
282 const BoostedXSystem& other = pcast<BoostedXSystem>(p);
283 return mkNamedPCmp(other, "RAPGAP") || mkNamedPCmp(other, "FS");
284 }
285
286 private:
287
288 LorentzTransform _boost;
289
290 };
291
292 }
293
294
295 /// @brief H1 diffractive dijets
296 ///
297 /// Diffractive dijets H1 with 920 GeV p and 27.5 GeV e
298 /// Tagged protons & jets found in gamma*p rest frame.
299 ///
300 /// @author Christine O. Rasmussen
301 class H1_2015_I1343110 : public Analysis {
302 public:
303
304 /// Constructor
305 RIVET_DEFAULT_ANALYSIS_CTOR(H1_2015_I1343110);
306
307 typedef H1_2015_I1343110_PROJECTIONS::RapidityGap RapidityGap;
308 typedef H1_2015_I1343110_PROJECTIONS::BoostedXSystem BoostedXSystem;
309
310 /// @name Analysis methods
311 /// @{
312
313 // Book projections and histograms
314 void init() {
315
316 declare(DISKinematics(), "Kinematics");
317 const DISFinalState& disfs = declare(DISFinalState(DISFrame::HCM), "DISFS");
318 const BoostedXSystem& disfsXcm = declare( BoostedXSystem(disfs), "BoostedXFS");
319 declare(FastJets(disfsXcm, fastjet::JetAlgorithm::kt_algorithm, fastjet::RecombinationScheme::pt_scheme, 1.0,
320 JetMuons::ALL, JetInvisibles::NONE, nullptr), "DISFSJets");
321 declare(DISDiffHadron(), "Hadron");
322 declare(RapidityGap(), "RapidityGap");
323
324 // Book histograms from REF data
325 book(_h_PHO_sig_sqrts, 1, 1, 1);
326 book(_h_DIS_sig_sqrts, 2, 1, 1);
327 book(_h_PHODIS_sqrts, 3, 1, 1);
328
329 book(_h_DIS_dsigdz, 4, 1, 1);
330 book(_h_DIS_dsigdxPom, 5, 1, 1);
331 book(_h_DIS_dsigdy, 6, 1, 1);
332 book(_h_DIS_dsigdQ2, 7, 1, 1);
333 book(_h_DIS_dsigdEtj1, 8, 1, 1);
334 book(_h_DIS_dsigdMX, 9, 1, 1);
335 book(_h_DIS_dsigdDeltaEta, 10, 1, 1);
336 book(_h_DIS_dsigdAvgEta, 11, 1, 1);
337
338 book(_h_PHO_dsigdz, 12, 1, 1);
339 book(_h_PHO_dsigdxPom, 13, 1, 1);
340 book(_h_PHO_dsigdy, 14, 1, 1);
341 book(_h_PHO_dsigdxGam, 15, 1, 1);
342 book(_h_PHO_dsigdEtj1, 16, 1, 1);
343 book(_h_PHO_dsigdMX, 17, 1, 1);
344 book(_h_PHO_dsigdDeltaEta, 18, 1, 1);
345 book(_h_PHO_dsigdAvgEta, 19, 1, 1);
346
347 book(_h_PHODIS_deltaEta, 20, 1, 1);
348 book(_h_PHODIS_y, 21, 1, 1);
349 book(_h_PHODIS_z, 22, 1, 1);
350 book(_h_PHODIS_Etj1, 23, 1, 1);
351
352 isPHO = false;
353 nVeto1 = 0;
354 nVeto2 = 0;
355 nVeto3 = 0;
356 nVeto4 = 0;
357 nVeto5 = 0;
358 nVeto6 = 0;
359 nPHO = 0;
360 nDIS = 0;
361 }
362
363 // Do the analysis
364 void analyze(const Event& event) {
365
366 // Event weight
367 isPHO = false;
368
369 // Projections - special handling of events where no proton found:
370 const RapidityGap& rg = apply<RapidityGap>(event, "RapidityGap");
371 const DISKinematics& kin = apply<DISKinematics>(event, "Kinematics");
372 const BoostedXSystem& disfsXcm = apply<BoostedXSystem>( event, "BoostedXFS");
373 Particle hadronOut;
374 Particle hadronIn;
375 try {
376 const DISDiffHadron& diffhadr = apply<DISDiffHadron>(event, "Hadron");
377 hadronOut = diffhadr.out();
378 hadronIn = diffhadr.in();
379 } catch (const Error& e){
380 vetoEvent;
381 }
382
383 // Determine kinematics: H1 has +z = proton direction
384 int dir = kin.orientation();
385 double y = kin.y();
386 double Q2 = kin.Q2();
387
388 // Separate into DIS and PHO regimes else veto
389 if (Q2 < 2.*GeV2 && inRange(y, 0.2, 0.70)) {
390 isPHO = true;
391 ++nPHO;
392 } else if (inRange(Q2, 4.0*GeV2, 80.*GeV2) && inRange(y, 0.2, 0.7)) {
393 isPHO = false;
394 ++nDIS;
395 } else vetoEvent;
396 ++nVeto1;
397
398 // Find diffractive variables as defined in paper.
399 // Note tagged protons in VFPS => smaller allowed xPom range
400 // xPom = 1 - E'/E, M2X from hadrons, t = (P-P')^2
401 const double M2X = rg.M2X();
402 const double abst = abs(rg.t());
403 const double xPom = 1. - hadronOut.energy() / hadronIn.energy();
404
405 //cout << "\nhadout=" << hadronOut.energy() << ", hadin=" << hadronIn.energy() << endl;
406 //cout << "xPomH1=" << (Q2+M2X) / (y * sqr(sqrtS())) << endl;
407 //cout << "|t|=" << abst << ", xPom=" << xPom << endl;
408 // Veto if outside allowed region
409 if (abst > 0.6 * GeV2) vetoEvent;
410 ++nVeto2;
411 if (!inRange(xPom, 0.010, 0.024)) vetoEvent;
412 ++nVeto3;
413
414 // Jet selection. Note jets are found in XCM frame, but
415 // eta cut is applied in lab frame!
416 Cut jetcuts = Cuts::Et > 4.* GeV;
417 Jets jets = apply<FastJets>(event, "DISFSJets").jets(jetcuts, cmpMomByEt);
418 // Veto if not dijets and if Et_j1 < 5.5
419 if (jets.size() < 2) vetoEvent;
420 if (jets[0].Et() < 5.5 * GeV) vetoEvent;
421 ++nVeto4;
422 // Find Et_jet1 in XCM frame
423 double EtJet1 = jets[0].Et() * GeV;
424
425 //cout << "gamma*p frame:" << endl;
426 //cout << "Et1=" << jets[0].Et() << ", E1=" << jets[0].E() << ", pz1=" << jets[0].pz() << ", m1=" << jets[0].mass() << endl;
427 //cout << "Et2=" << jets[1].Et() << ", E2=" << jets[1].E() << ", pz2=" << jets[1].pz() << ", m2=" << jets[1].mass() << endl;
428
429 // Transform from XCM to HCM
430 const LorentzTransform xcmboost = disfsXcm.boost();
431 for (int i = 0; i < 2; ++i) jets[i].transformBy(xcmboost.inverse());
432
433 // Find mass of jets and EpPz, EmPz of jets in HCM frame.
434 FourMomentum momJets = jets[0].momentum() + jets[1].momentum();
435 double M2jets = momJets.mass2();
436 double EpPzJets = 0.;
437 double EmPzJets = 0.;
438 // Note sign change wrt. H1 because photon is in +z direction
439 for (int i = 0; i < 2; ++i){
440 EpPzJets += jets[i].E() - jets[i].pz(); // Sign: + => -
441 EmPzJets += jets[i].E() + jets[i].pz(); // Sign: - => +
442 }
443
444 // Transform the jets from HCM to LAB frame where eta cut is
445 // applied for photoproduction.
446 const LorentzTransform hcmboost = kin.boostHCM();
447 for (int i = 0; i < 2; ++i) jets[i].transformBy(hcmboost.inverse());
448 double etaLabJet1 = dir * jets[0].eta();
449 double etaLabJet2 = dir * jets[1].eta();
450 if (!inRange(etaLabJet1, -1., 2.5)) vetoEvent;
451 if (!inRange(etaLabJet2, -1., 2.5)) vetoEvent;
452 ++nVeto5;
453
454 // Pseudorapidity distributions are examined in lab frame:
455 double deltaEtaJets = abs(dir * jets[0].eta() - dir * jets[1].eta());
456 double avgEtaJets = 0.5 * (dir * jets[0].eta() + dir * jets[1].eta());
457
458 // Evaluate observables
459 double zPomJets, xGamJets = 0.;
460 if (isPHO){
461 zPomJets = EpPzJets / rg.EpPzX(RapidityGap::HCM);
462 xGamJets = EmPzJets / rg.EmPzX(RapidityGap::HCM);
463 //cout << "xGamJets=" << xGamJets << endl;
464 } else {
465 zPomJets = (Q2 + M2jets) / (Q2 + M2X);
466 }
467
468 //cout << "lab frame:" << endl;
469 //cout << "Et1=" << jets[0].Et() << ", E1=" << jets[0].E() << ", pz1=" << jets[0].pz() << ", m1=" << jets[0].mass() << endl;
470 //cout << "Et2=" << jets[1].Et() << ", E2=" << jets[1].E() << ", pz2=" << jets[1].pz() << ", m2=" << jets[1].mass() << endl;
471 //cout << "EpPzJets=" << EpPzJets << ", EmPzJets=" << EmPzJets << endl;
472 //cout << "Et*exp(eta)=" << jets[0].Et()*exp(etaLabJet1) + jets[1].Et()*exp(etaLabJet2) << endl;
473 //cout << "Et*exp(-eta)=" << jets[0].Et()*exp(-etaLabJet1) + jets[1].Et()*exp(-etaLabJet2) << endl;
474 //cout << "EpPz=" << rg.EpPzX(RapidityGap::HCM) << ", EmPz=" << rg.EmPzX(RapidityGap::HCM) << endl;
475 //cout << "2 xPom Ep=" << 2. * xPom * kin.beamHadron().E() << ", 2 y Ee=" << 2. * y * kin.beamLepton().E() << endl;
476 //cout << "xGam=" << xGamJets << ", zPom=" << zPomJets << endl;
477 //cout << "M12=" << M2jets << ", deltaEta=" << deltaEtaJets << ", avgEta=" << avgEtaJets << endl;
478
479 // Veto events with zPom > 0.8
480 if (zPomJets > 0.8) vetoEvent;
481 ++nVeto6;
482
483 // Now fill histograms
484 if (isPHO){
485 _h_PHO_sig_sqrts ->fill(319);
486 _h_PHO_dsigdz ->fill(zPomJets);
487 _h_PHO_dsigdxPom ->fill(xPom);
488 _h_PHO_dsigdy ->fill(y);
489 _h_PHO_dsigdxGam ->fill(xGamJets);
490 _h_PHO_dsigdEtj1 ->fill(EtJet1);
491 _h_PHO_dsigdMX ->fill(sqrt(M2X)/GeV);
492 _h_PHO_dsigdDeltaEta ->fill(deltaEtaJets);
493 _h_PHO_dsigdAvgEta ->fill(avgEtaJets);
494 } else {
495 _h_DIS_sig_sqrts ->fill(319);
496 _h_DIS_dsigdz ->fill(zPomJets);
497 _h_DIS_dsigdxPom ->fill(xPom);
498 _h_DIS_dsigdy ->fill(y);
499 _h_DIS_dsigdQ2 ->fill(Q2);
500 _h_DIS_dsigdEtj1 ->fill(EtJet1);
501 _h_DIS_dsigdMX ->fill(sqrt(M2X)/GeV);
502 _h_DIS_dsigdDeltaEta ->fill(deltaEtaJets);
503 _h_DIS_dsigdAvgEta ->fill(avgEtaJets);
504 }
505
506 }
507
508 // Finalize
509 void finalize() {
510 // Normalise to cross section
511 // Remember to manually scale the cross section afterwards with
512 // the number of rejected events.
513 const double norm = crossSection()/picobarn/sumOfWeights();
514
515 scale(_h_PHO_sig_sqrts, norm);
516 scale(_h_PHO_dsigdz, norm);
517 scale(_h_PHO_dsigdxPom, norm);
518 scale(_h_PHO_dsigdy, norm);
519 scale(_h_PHO_dsigdxGam, norm);
520 scale(_h_PHO_dsigdEtj1, norm);
521 scale(_h_PHO_dsigdMX, norm);
522 scale(_h_PHO_dsigdDeltaEta, norm);
523 scale(_h_PHO_dsigdAvgEta, norm);
524
525 scale(_h_DIS_sig_sqrts, norm);
526 scale(_h_DIS_dsigdz, norm);
527 scale(_h_DIS_dsigdxPom, norm);
528 scale(_h_DIS_dsigdy, norm);
529 scale(_h_DIS_dsigdQ2, norm);
530 scale(_h_DIS_dsigdEtj1, norm);
531 scale(_h_DIS_dsigdMX, norm);
532 scale(_h_DIS_dsigdDeltaEta, norm);
533 scale(_h_DIS_dsigdAvgEta, norm);
534
535 if (_h_DIS_sig_sqrts->numEntries() != 0)
536 divide(_h_PHO_sig_sqrts, _h_DIS_sig_sqrts, _h_PHODIS_sqrts);
537 if (_h_DIS_dsigdDeltaEta->numEntries() != 0)
538 divide(_h_PHO_dsigdDeltaEta, _h_DIS_dsigdDeltaEta, _h_PHODIS_deltaEta);
539 if (_h_DIS_dsigdy->numEntries() != 0)
540 divide(_h_PHO_dsigdy, _h_DIS_dsigdy, _h_PHODIS_y);
541 if (_h_DIS_dsigdz->numEntries() != 0)
542 divide(_h_PHO_dsigdz, _h_DIS_dsigdz, _h_PHODIS_z);
543 if (_h_DIS_dsigdEtj1->numEntries() != 0)
544 divide(_h_PHO_dsigdEtj1, _h_DIS_dsigdEtj1, _h_PHODIS_Etj1);
545
546 const double dPHO = nPHO;
547 MSG_INFO("H1_2015_I1343110");
548 MSG_INFO("Cross section = " << crossSection()/picobarn << " pb");
549 MSG_INFO("Number of events = " << numEvents() << ", sumW = " << sumOfWeights());
550 MSG_INFO("Number of PHO = " << nPHO << ", number of DIS = " << nDIS);
551 MSG_INFO("Events passing electron veto = " << nVeto1 << " (" << nVeto1/dPHO * 100. << "%)" );
552 MSG_INFO("Events passing t veto = " << nVeto2 << " (" << nVeto2/dPHO * 100. << "%)" );
553 MSG_INFO("Events passing xPom = " << nVeto3 << " (" << nVeto3/dPHO * 100. << "%)" );
554 MSG_INFO("Events passing jet Et veto = " << nVeto4 << " (" << nVeto4/dPHO * 100. << "%)" );
555 MSG_INFO("Events passing jet eta veto = " << nVeto5 << " (" << nVeto5/dPHO * 100. << "%)" );
556 MSG_INFO("Events passing zPom veto = " << nVeto6 << " (" << nVeto6/dPHO * 100. << "%)" );
557
558 }
559
560 /// @}
561
562
563 private:
564
565 /// @name Histograms
566 /// @{
567 // Book histograms from REF data
568 BinnedHistoPtr<int> _h_PHO_sig_sqrts;
569 BinnedHistoPtr<int> _h_DIS_sig_sqrts;
570 BinnedEstimatePtr<int> _h_PHODIS_sqrts;
571
572 Histo1DPtr _h_DIS_dsigdz;
573 Histo1DPtr _h_DIS_dsigdxPom;
574 Histo1DPtr _h_DIS_dsigdy;
575 Histo1DPtr _h_DIS_dsigdQ2;
576 Histo1DPtr _h_DIS_dsigdEtj1;
577 Histo1DPtr _h_DIS_dsigdMX;
578 Histo1DPtr _h_DIS_dsigdDeltaEta;
579 Histo1DPtr _h_DIS_dsigdAvgEta;
580
581 Histo1DPtr _h_PHO_dsigdz;
582 Histo1DPtr _h_PHO_dsigdxPom;
583 Histo1DPtr _h_PHO_dsigdy;
584 Histo1DPtr _h_PHO_dsigdxGam;
585 Histo1DPtr _h_PHO_dsigdEtj1;
586 Histo1DPtr _h_PHO_dsigdMX;
587 Histo1DPtr _h_PHO_dsigdDeltaEta;
588 Histo1DPtr _h_PHO_dsigdAvgEta;
589
590 Estimate1DPtr _h_PHODIS_deltaEta;
591 Estimate1DPtr _h_PHODIS_y;
592 Estimate1DPtr _h_PHODIS_z;
593 Estimate1DPtr _h_PHODIS_Etj1;
594 /// @}
595
596 bool isPHO;
597 int nVeto1, nVeto2, nVeto3, nVeto4, nVeto5, nVeto6;
598 int nPHO, nDIS;
599 };
600
601 RIVET_DECLARE_PLUGIN(H1_2015_I1343110);
602
603}
|