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}