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