diff --git a/analyses/pluginHERA/H1_2007_I746380.cc b/analyses/pluginHERA/H1_2007_I746380.cc
new file mode 100644
--- /dev/null
+++ b/analyses/pluginHERA/H1_2007_I746380.cc
@@ -0,0 +1,608 @@
+// -*- C++ -*-
+#include "Rivet/Analysis.hh"
+#include "Rivet/Projections/Beam.hh"
+#include "Rivet/Projections/FinalState.hh"
+#include "Rivet/Projections/DISKinematics.hh"
+#include "Rivet/Projections/DISFinalState.hh"
+#include "Rivet/Projections/FastJets.hh"
+
+namespace Rivet {
+
+ // Projection to find the largest gaps and the masses of the two
+ // systems separated by the gap. Based on the HZTools gap-finding
+ // method (hzhadgap.F). Note that gaps are found in the HCM frame.
+ // Author Christine O. Rasmussen.
+ class RapidityGap : public Projection {
+
+ public:
+
+ /// Type of DIS boost to apply
+ enum Frame { HCM, LAB, XCM };
+
+ RapidityGap() {
+ setName("RapidityGap");
+ addProjection(DISKinematics(), "DISKIN");
+ addProjection(DISFinalState(DISFinalState::HCM), "DISFS");
+ }
+
+ DEFAULT_RIVET_PROJ_CLONE(RapidityGap);
+
+ const double M2X() const {return _M2X;}
+ const double M2Y() const {return _M2Y;}
+ const double t() const {return _t;}
+ const double gap() const {return _gap;}
+ const double gapUpp() const {return _gapUpp;}
+ const double gapLow() const {return _gapLow;}
+ const double EpPzX(Frame f) const {
+ if (f == LAB) return _ePpzX_LAB;
+ else if (f == XCM) return _ePpzX_XCM;
+ return _ePpzX_HCM;
+ }
+ const double EmPzX(Frame f) const {
+ if (f == LAB) return _eMpzX_LAB;
+ else if (f == XCM) return _eMpzX_XCM;
+ return _eMpzX_HCM;
+ }
+ const FourMomentum pX(Frame f) const {
+ if (f == LAB) return _momX_LAB;
+ else if (f == XCM) return _momX_XCM;
+ return _momX_HCM;
+ }
+ const FourMomentum pY(Frame f) const {
+ if (f == LAB) return _momY_LAB;
+ else if (f == XCM) return _momY_XCM;
+ return _momY_HCM;
+ }
+ const Particles& systemX(Frame f) const {
+ if (f == LAB) return _pX_LAB;
+ else if (f == XCM) return _pX_XCM;
+ return _pX_HCM;
+ }
+ const Particles& systemY(Frame f) const {
+ if (f == LAB) return _pY_LAB;
+ else if (f == XCM) return _pY_XCM;
+ return _pY_HCM;
+ }
+
+ protected:
+
+ virtual int compare(const Projection& p) const {
+ const RapidityGap& other = pcast<RapidityGap>(p);
+ return mkNamedPCmp(other, "DISKIN") || mkNamedPCmp(other, "DISFS");
+ }
+
+ virtual void project(const Event& e){
+ const DISKinematics& dk = apply<DISKinematics>(e, "DISKIN");
+ const Particles& p = apply<DISFinalState>(e, "DISFS").particles(cmpMomByEta);
+ findgap(p, dk);
+ }
+
+ void clearAll(){
+ _M2X = _M2Y = _t = 0.;
+ _gap = 0.;
+ _gapUpp = _gapLow = -8.;
+ _ePpzX_HCM = _eMpzX_HCM =_ePpzX_LAB = _eMpzX_LAB = _ePpzX_XCM = _eMpzX_XCM = 0.;
+ }
+
+ void findgap(const Particles& particles, const DISKinematics& diskin){
+
+ clearAll();
+
+ // Begin by finding largest gap and gapedges between all final
+ // state particles in HCM frame.
+ int nP = particles.size();
+ int dir = diskin.orientation();
+ for (int i = 0; i < nP-1; ++i){
+ double tmpGap = dir * particles[i+1].eta()
+ - dir * particles[i].eta();
+ if (tmpGap > _gap) {
+ _gap = tmpGap;
+ _gapLow = dir * particles[i].eta();
+ _gapUpp = dir * particles[i+1].eta();
+ }
+ }
+
+ // Define the two systems X and Y.
+ Particles pX, pY;
+ foreach (const Particle& ip, particles) {
+ if (dir * ip.eta() > _gapLow) pX.push_back(ip);
+ else pY.push_back(ip);
+ }
+
+ // Find variables related to HCM frame.
+ // Note that HCM has photon along +z, as opposed to
+ // H1 where proton is along +z. This results in a sign change
+ // as compared to H1 papers!
+ // X - side
+ FourMomentum momX;
+ foreach (const Particle& jp, pX) {
+ momX += jp.momentum();
+ _ePpzX_HCM += jp.E() - jp.pz(); // Sign + => -
+ _eMpzX_HCM += jp.E() + jp.pz(); // Sign - => +
+ }
+ _momX_HCM = momX;
+ _pX_HCM = pX;
+ _M2X = _momX_HCM.mass2();
+ // Y - side
+ FourMomentum momY;
+ foreach (const Particle& kp, pY) momY += kp.momentum();
+ _momY_HCM = momY;
+ _pY_HCM = pY;
+ _M2Y = _momY_HCM.mass2();
+
+ // Find variables related to LAB frame
+ const LorentzTransform hcmboost = diskin.boostHCM();
+ const LorentzTransform hcminverse = hcmboost.inverse();
+ _momX_LAB = hcminverse.transform(_momX_HCM);
+ _momY_LAB = hcminverse.transform(_momY_HCM);
+
+ // Find momenta in XCM frame. Note that it is HCM frame that is
+ // boosted, resulting in a sign change later!
+ const bool doXCM = (momX.betaVec().mod2() < 1.);
+ if (doXCM) {
+ const LorentzTransform xcmboost =
+ LorentzTransform::mkFrameTransformFromBeta(momX.betaVec());
+ _momX_XCM = xcmboost.transform(momX);
+ _momY_XCM = xcmboost.transform(momY);
+ } else {
+ _momX_XCM = momX;
+ _momY_XCM = momY;
+ }
+
+ double pInPlusZ = (diskin.beamHadron().pz() > 0.) ? 1. : -1.;
+ foreach (const Particle& jp, pX) {
+ // Boost from HCM to LAB. Here check if proton is along +z
+ // and adjust sign wrt. this.
+ FourMomentum lab = hcminverse.transform(jp.momentum());
+ _ePpzX_LAB += lab.E() + pInPlusZ * lab.pz();
+ _eMpzX_LAB += lab.E() - pInPlusZ * lab.pz();
+ Particle plab = jp;
+ plab.setMomentum(lab);
+ _pX_LAB.push_back(plab);
+ // Set XCM. Note that since HCM frame is boosted to XCM frame,
+ // we have a sign change
+ if (doXCM) {
+ const LorentzTransform xcmboost =
+ LorentzTransform::mkFrameTransformFromBeta(_momX_HCM.betaVec());
+ FourMomentum xcm = xcmboost.transform(jp.momentum());
+ _ePpzX_XCM += xcm.E() - xcm.pz(); // Sign + => -
+ _eMpzX_XCM += xcm.E() + xcm.pz(); // Sign - => +
+ Particle pxcm = jp;
+ pxcm.setMomentum(xcm);
+ _pX_XCM.push_back(pxcm);
+ }
+ }
+
+ foreach (const Particle& jp, pY) {
+ // Boost from HCM to LAB
+ FourMomentum lab = hcminverse.transform(jp.momentum());
+ Particle plab = jp;
+ plab.setMomentum(lab);
+ _pY_LAB.push_back(plab);
+ // Boost from HCM to XCM
+ if (doXCM) {
+ const LorentzTransform xcmboost =
+ LorentzTransform::mkFrameTransformFromBeta(_momX_HCM.betaVec());
+ FourMomentum xcm = xcmboost.transform(jp.momentum());
+ Particle pxcm = jp;
+ pxcm.setMomentum(xcm);
+ _pY_XCM.push_back(pxcm);
+ }
+ }
+ if (!doXCM) {
+ _ePpzX_XCM = _ePpzX_HCM;
+ _eMpzX_XCM = _eMpzX_HCM;
+ _pX_XCM = _pX_HCM;
+ _pY_XCM = _pY_HCM;
+ }
+
+ // Find t: Currently can only handle gap on proton side.
+ // @TODO: Expand to also handle gap on photon side
+ // Boost p from LAB to HCM frame to find t.
+ const FourMomentum proton = hcmboost.transform(diskin.beamHadron().momentum());
+ FourMomentum pPom = proton - _momY_HCM;
+ _t = pPom * pPom;
+
+ }
+
+ private:
+
+ double _M2X, _M2Y, _t;
+ double _gap, _gapUpp, _gapLow;
+ double _ePpzX_LAB, _eMpzX_LAB, _ePpzX_HCM, _eMpzX_HCM, _ePpzX_XCM, _eMpzX_XCM;
+ FourMomentum _momX_HCM, _momY_HCM,_momX_LAB, _momY_LAB, _momX_XCM, _momY_XCM;
+ Particles _pX_HCM, _pY_HCM, _pX_LAB, _pY_LAB, _pX_XCM, _pY_XCM;
+
+ };
+
+ // Projection to pick out the scattered beam proton with diffraction.
+ // Author Ilkka Helenius
+ class DiffHadron : public Projection {
+
+ public:
+
+ /// @name Constructors.
+ //@{
+
+ DiffHadron(){
+ setName("DiffHadron");
+ addProjection(Beam(), "Beam");
+ addProjection(PromptFinalState(), "PromptFS");
+ }
+
+ /// Clone on the heap.
+ DEFAULT_RIVET_PROJ_CLONE(DiffHadron);
+
+ //@}
+
+ protected:
+
+ void project(const Event& e) {
+
+ // Find incoming hadron beam
+ const ParticlePair& inc = applyProjection<Beam>(e, "Beam").beams();
+ bool firstIsHadron = PID::isHadron(inc.first.pid());
+ bool secondIsHadron = PID::isHadron(inc.second.pid());
+ if (firstIsHadron && !secondIsHadron) {
+ _incoming = inc.first;
+ } else if (!firstIsHadron && secondIsHadron) {
+ _incoming = inc.second;
+ } else {
+ throw Error("DiffHadron could not find the correct beam");
+ }
+
+ const Particles fshadrons = applyProjection<FinalState>(e, "PromptFS").particles(isHadron, cmpMomByE);
+ const Particles sfhadrons = filter_select(fshadrons, Cuts::pid == _incoming.pid());
+ if (!sfhadrons.empty()) {
+ _outgoing = sfhadrons.front();
+ } else if (!fshadrons.empty()) {
+ _outgoing = fshadrons.front();
+ } else {
+ throw Error("Could not find the scattered hadron");
+ }
+
+ }
+
+ int compare(const Projection& p) const {
+ const DiffHadron& other = pcast<DiffHadron>(p);
+ return mkNamedPCmp(other, "Beam") || mkNamedPCmp(other, "PromptFS");
+ }
+
+
+ public:
+
+ /// The incoming hadron
+ const Particle& in() const { return _incoming; }
+
+ /// The outgoing hadron
+ const Particle& out() const { return _outgoing; }
+
+ /// Sign of the incoming lepton pz component
+ int pzSign() const { return sign(_incoming.pz()); }
+
+
+ private:
+
+ /// The incoming hadron
+ Particle _incoming;
+
+ /// The outgoing hadron
+ Particle _outgoing;
+
+ // /// The charge sign of the DIS current
+ // double _charge;
+
+ };
+
+ // Projection to boost system X (photon+Pomeron) particles into its rest frame.
+ // Author Ilkka Helenius
+ class BoostedXSystem : public FinalState {
+ public:
+
+ BoostedXSystem(const FinalState& fs) {
+ setName("BoostedXSystem");
+ declare(fs,"FS");
+ addProjection(RapidityGap(), "RAPGAP");
+ }
+
+ // Return the boost to XCM frame.
+ const LorentzTransform& boost() const { return _boost; }
+
+ DEFAULT_RIVET_PROJ_CLONE(BoostedXSystem);
+
+ protected:
+
+ // Apply the projection on the supplied event.
+ void project(const Event& e){
+
+ const RapidityGap& rg = apply<RapidityGap>(e, "RAPGAP");
+
+ // Total momentum of the system X.
+ const FourMomentum pX = rg.pX(RapidityGap::HCM);
+
+ // Reset the boost. Is there a separate method for this?
+ _boost = combine(_boost, _boost.inverse());
+
+ // Define boost only when numerically safe, otherwise negligible.
+ if (pX.betaVec().mod2() < 1.)
+ _boost = LorentzTransform::mkFrameTransformFromBeta(pX.betaVec());
+
+ // Boost the particles from system X.
+ _theParticles.clear();
+ _theParticles.reserve(rg.systemX(RapidityGap::HCM).size());
+ for (const Particle& p : rg.systemX(RapidityGap::HCM)) {
+ Particle temp = p;
+ temp.setMomentum(_boost.transform(temp.momentum()));
+ _theParticles.push_back(temp);
+ }
+
+ }
+
+ // Compare projections.
+ int compare(const Projection& p) const {
+ const BoostedXSystem& other = pcast<BoostedXSystem>(p);
+ return mkNamedPCmp(other, "RAPGAP") || mkNamedPCmp(other, "FS");
+ }
+
+ private:
+
+ LorentzTransform _boost;
+
+ };
+
+ /// @brief H1 diffractive dijets
+ ///
+ /// Diffractive dijets H1 with 920 GeV p and 27.5 GeV e
+ /// Note tagged protons!
+ ///
+ /// @author Christine O. Rasmussen
+ class H1_2007_I746380 : public Analysis {
+ public:
+
+ /// Constructor
+ DEFAULT_RIVET_ANALYSIS_CTOR(H1_2007_I746380);
+
+ /// @name Analysis methods
+ //@{
+
+ // Book projections and histograms
+ void init() {
+
+ declare(DISKinematics(), "Kinematics");
+ const DISFinalState& disfs = declare(DISFinalState(DISFinalState::HCM), "DISFS");
+ const BoostedXSystem& disfsXcm = declare( BoostedXSystem(disfs), "BoostedXFS");
+ declare(FastJets(disfsXcm, fastjet::JetAlgorithm::kt_algorithm, fastjet::RecombinationScheme::pt_scheme, 1.0,
+ JetAlg::ALL_MUONS, JetAlg::NO_INVISIBLES, nullptr), "DISFSJets");
+ declare(DiffHadron(), "Hadron");
+ declare(RapidityGap(), "RapidityGap");
+
+ // Book histograms from REF data
+ _h_DIS_dsigdzPom = bookHisto1D(1, 1, 1);
+ _h_DIS_dsigdlogXpom = bookHisto1D(2, 1, 1);
+ _h_DIS_dsigdW = bookHisto1D(3, 1, 1);
+ _h_DIS_dsigdQ2 = bookHisto1D(4, 1, 1);
+ _h_DIS_dsigdEtJet1 = bookHisto1D(5, 1, 1);
+ _h_DIS_dsigdAvgEta = bookHisto1D(6, 1, 1);
+ _h_DIS_dsigdDeltaEta = bookHisto1D(7, 1, 1);
+
+ _h_PHO_dsigdzPom = bookHisto1D(8, 1, 1);
+ _h_PHO_dsigdxGam = bookHisto1D(9, 1, 1);
+ _h_PHO_dsigdlogXpom = bookHisto1D(10, 1, 1);
+ _h_PHO_dsigdW = bookHisto1D(11, 1, 1);
+ _h_PHO_dsigdEtJet1 = bookHisto1D(12, 1, 1);
+ _h_PHO_dsigdAvgEta = bookHisto1D(13, 1, 1);
+ _h_PHO_dsigdDeltaEta = bookHisto1D(14, 1, 1);
+ _h_PHO_dsigdMjets = bookHisto1D(15, 1, 1);
+
+ isDIS = false;
+ nVeto0 = 0;
+ nVeto1 = 0;
+ nVeto2 = 0;
+ nVeto3 = 0;
+ nVeto4 = 0;
+ nVeto5 = 0;
+ nPHO = 0;
+ nDIS = 0;
+ }
+
+ // Do the analysis
+ void analyze(const Event& event) {
+
+ // Event weight
+ const double weight = event.weight();
+ isDIS = false;
+
+ // Projections - special handling of events where no proton found:
+ const RapidityGap& rg = apply<RapidityGap>(event, "RapidityGap");
+ const DISKinematics& kin = apply<DISKinematics>(event, "Kinematics");
+ const BoostedXSystem& disfsXcm = apply<BoostedXSystem>( event, "BoostedXFS");
+
+ // Determine kinematics: H1 has +z = proton direction
+ int dir = kin.orientation();
+ double W2 = kin.W2();
+ double W = sqrt(W2);
+ double y = kin.y();
+ double Q2 = kin.Q2();
+
+ // Separate into DIS and PHO regimes else veto
+ if (!inRange(W, 165.*GeV, 242.*GeV)) vetoEvent;
+ if (Q2 < 0.01*GeV2) {
+ isDIS = false;
+ ++nPHO;
+ } else if (inRange(Q2, 4.0*GeV2, 80.*GeV2)) {
+ isDIS = true;
+ ++nDIS;
+ } else {
+ vetoEvent;
+ }
+ ++nVeto0;
+
+ // Find diffractive variables as defined in paper.
+ const double M2Y = rg.M2Y();
+ const double M2X = rg.M2X();
+ const double abst = abs(rg.t());
+ const double xPom = (isDIS) ? (Q2 + M2X) / (Q2 + W2) :
+ rg.EpPzX(RapidityGap::LAB) / (2. * kin.beamHadron().E());
+
+ // Veto if outside allowed region
+ if (sqrt(M2Y) > 1.6*GeV) vetoEvent;
+ ++nVeto1;
+ if (abst > 1.0*GeV2) vetoEvent;
+ ++nVeto2;
+ if (xPom > 0.03) vetoEvent;
+ ++nVeto3;
+
+ // Jet selection. Note jets are found in photon-proton (XCM)
+ // frame, but eta cut is applied in lab frame!
+ Cut jetcuts = Cuts::Et > 4.* GeV;
+ Jets jets = apply<FastJets>(event, "DISFSJets").jets(jetcuts, cmpMomByEt);
+ // Veto if not dijets and if Et_j1 < 5.0
+ if (jets.size() < 2) vetoEvent;
+ if (jets[0].Et() < 5.*GeV) vetoEvent;
+ ++nVeto4;
+ // Find Et_jet1 and deltaEta* in XCM frame
+ double EtJet1 = jets[0].Et() * GeV;
+ double etaXCMJet1 = dir * jets[0].eta();
+ double etaXCMJet2 = dir * jets[1].eta();
+ double deltaEtaJets = abs(etaXCMJet1 - etaXCMJet2);
+
+ // Transform from XCM to HCM
+ const LorentzTransform xcmboost = disfsXcm.boost();
+ for (int i = 0; i < 2; ++i) jets[i].transformBy(xcmboost.inverse());
+ // Find mass of jets and EpPz, EmPz of jets
+ FourMomentum momJets = jets[0].momentum() + jets[1].momentum();
+ double M2jets = momJets.mass2();
+ double EpPzJets = 0.;
+ double EmPzJets = 0.;
+ // DIS variables are found in XCM frame, so boost back again
+ if (isDIS){
+ for (int i = 0; i < 2; ++i) jets[i].transformBy(xcmboost);
+ }
+ // Note sign change wrt. H1 because photon is in +z direction
+ for (int i = 0; i < 2; ++i){
+ EpPzJets += jets[i].E() - jets[i].pz(); // Sign: + => -
+ EmPzJets += jets[i].E() + jets[i].pz(); // Sign: - => +
+ }
+
+ // Transform the jets from HCM to LAB frame where eta cut is
+ // applied for photoproduction.
+ const LorentzTransform hcmboost = kin.boostHCM();
+ for (int i = 0; i < 2; ++i) jets[i].transformBy(hcmboost.inverse());
+ double etaLabJet1 = dir * jets[0].eta();
+ double etaLabJet2 = dir * jets[1].eta();
+ double etaMin = (isDIS) ? -3. : -1.;
+ double etaMax = (isDIS) ? 0. : 2.;
+ double eta1 = (isDIS) ? etaXCMJet1 : etaLabJet1;
+ double eta2 = (isDIS) ? etaXCMJet2 : etaLabJet2;
+ if (!inRange(eta1, etaMin, etaMax)) vetoEvent;
+ if (!inRange(eta2, etaMin, etaMax)) vetoEvent;
+ ++nVeto5;
+
+ // Pseudorapidity distributions are examined in lab frame:
+ double avgEtaJets = 0.5 * (etaLabJet1 + etaLabJet2);
+
+ double zPomJets, xGamJets;
+ if (isDIS) {
+ zPomJets = (Q2 + M2jets) / (Q2 + M2X);
+ xGamJets = EmPzJets / rg.EmPzX(RapidityGap::XCM);
+ } else {
+ // Boost E_p, E_e to HCM frame
+ const LorentzTransform hcminverse = hcmboost.inverse();
+ FourMomentum lep = hcminverse.transform(kin.beamLepton().momentum());
+ FourMomentum had = hcminverse.transform(kin.beamHadron().momentum());
+ zPomJets = EpPzJets / (2. * xPom * had.E());
+ xGamJets = EmPzJets / (2. * y * lep.E());
+ }
+
+ // Now fill histograms
+ if (isDIS){
+ _h_DIS_dsigdzPom ->fill(zPomJets, weight);
+ _h_DIS_dsigdlogXpom ->fill(log10(xPom), weight);
+ _h_DIS_dsigdW ->fill(W, weight);
+ _h_DIS_dsigdQ2 ->fill(Q2, weight);
+ _h_DIS_dsigdEtJet1 ->fill(EtJet1, weight);
+ _h_DIS_dsigdAvgEta ->fill(avgEtaJets, weight);
+ _h_DIS_dsigdDeltaEta ->fill(deltaEtaJets, weight);
+ } else {
+ _h_PHO_dsigdzPom ->fill(zPomJets, weight);
+ _h_PHO_dsigdxGam ->fill(xGamJets, weight);
+ _h_PHO_dsigdlogXpom ->fill(log10(xPom), weight);
+ _h_PHO_dsigdW ->fill(W, weight);
+ _h_PHO_dsigdEtJet1 ->fill(EtJet1, weight);
+ _h_PHO_dsigdAvgEta ->fill(avgEtaJets, weight);
+ _h_PHO_dsigdDeltaEta ->fill(deltaEtaJets, weight);
+ _h_PHO_dsigdMjets ->fill(sqrt(M2jets), weight);
+ }
+
+ }
+
+ // Finalize
+ void finalize() {
+ // Normalise to cross section
+ const double norm = crossSection()/picobarn/sumOfWeights();
+
+ scale( _h_DIS_dsigdzPom , norm);
+ scale( _h_DIS_dsigdlogXpom , norm);
+ scale( _h_DIS_dsigdW , norm);
+ scale( _h_DIS_dsigdQ2 , norm);
+ scale( _h_DIS_dsigdEtJet1 , norm);
+ scale( _h_DIS_dsigdAvgEta , norm);
+ scale( _h_DIS_dsigdDeltaEta, norm);
+
+ scale( _h_PHO_dsigdzPom , norm);
+ scale( _h_PHO_dsigdxGam , norm);
+ scale( _h_PHO_dsigdlogXpom , norm);
+ scale( _h_PHO_dsigdW , norm);
+ scale( _h_PHO_dsigdEtJet1 , norm);
+ scale( _h_PHO_dsigdAvgEta , norm);
+ scale( _h_PHO_dsigdDeltaEta, norm);
+ scale( _h_PHO_dsigdMjets , norm);
+
+ const double dPHO = nPHO;
+ MSG_INFO("H1_2007_I746380");
+ MSG_INFO("Cross section = " << crossSection()/picobarn << " pb");
+ MSG_INFO("Number of events = " << numEvents() << ", sumW = " << sumOfWeights());
+ MSG_INFO("Number of PHO = " << nPHO << ", number of DIS = " << nDIS);
+ MSG_INFO("Events passing electron veto = " << nVeto0 << " (" << nVeto0/dPHO * 100. << "%)" );
+ MSG_INFO("Events passing MY = " << nVeto1 << " (" << nVeto1/dPHO * 100. << "%)" );
+ MSG_INFO("Events passing t veto = " << nVeto2 << " (" << nVeto2/dPHO * 100. << "%)" );
+ MSG_INFO("Events passing xPom = " << nVeto3 << " (" << nVeto3/dPHO * 100. << "%)" );
+ MSG_INFO("Events passing jet Et veto = " << nVeto4 << " (" << nVeto4/dPHO * 100. << "%)" );
+ MSG_INFO("Events passing jet eta veto = " << nVeto5 << " (" << nVeto5/dPHO * 100. << "%)" );
+
+ }
+
+ //@}
+
+
+ private:
+
+ /// @name Histograms
+ //@{
+ // Book histograms from REF data
+ Histo1DPtr _h_DIS_dsigdzPom ;
+ Histo1DPtr _h_DIS_dsigdlogXpom ;
+ Histo1DPtr _h_DIS_dsigdW ;
+ Histo1DPtr _h_DIS_dsigdQ2 ;
+ Histo1DPtr _h_DIS_dsigdEtJet1 ;
+ Histo1DPtr _h_DIS_dsigdAvgEta ;
+ Histo1DPtr _h_DIS_dsigdDeltaEta;
+
+ Histo1DPtr _h_PHO_dsigdzPom ;
+ Histo1DPtr _h_PHO_dsigdxGam ;
+ Histo1DPtr _h_PHO_dsigdlogXpom ;
+ Histo1DPtr _h_PHO_dsigdW ;
+ Histo1DPtr _h_PHO_dsigdEtJet1 ;
+ Histo1DPtr _h_PHO_dsigdAvgEta ;
+ Histo1DPtr _h_PHO_dsigdDeltaEta;
+ Histo1DPtr _h_PHO_dsigdMjets ;
+ //@}
+
+ bool isDIS;
+ int nVeto0, nVeto1, nVeto2, nVeto3, nVeto4, nVeto5;
+ int nPHO, nDIS;
+ };
+
+ DECLARE_RIVET_PLUGIN(H1_2007_I746380);
+
+}
diff --git a/analyses/pluginHERA/H1_2007_I746380.info b/analyses/pluginHERA/H1_2007_I746380.info
new file mode 100644
--- /dev/null
+++ b/analyses/pluginHERA/H1_2007_I746380.info
@@ -0,0 +1,43 @@
+Name: H1_2007_I746380
+Year: 2007
+Summary: Tests of QCD Factorisation in the diffractive production of dijets in deep-inelastic scattering and photoproduction at HERA
+Experiment: H1
+Collider: HERA
+InspireID: 746380
+Status: UNVALIDATED
+Authors:
+ - Christine O. Rasmussen <christine.rasmussen@thep.lu.se>
+ - Ilkka Helenius <ilkka.helenius@uni-tuebingen.de>
+References:
+- Eur. Phys. J. C51 (2007) 549
+- arXiv:hep-ex/0703022
+RunInfo:
+ 820 GeV protons colliding with 27.5 GeV positrons;
+ Diffractive photoproduction of dijets;
+ Jet $E_T > 4$ GeV;
+NumEvents: 1000000
+Beams: [p+, e+]
+Energies: [[820, 27.5]]
+PtCuts: [0]
+Description:
+ H1 diffractive jets from proton--positron
+ collisions at beam energies of 820~GeV on 27.5~GeV.
+BibKey: Andreev:2015cwa
+BibTeX: '@article{Aktas:2007hn,
+ author = "Aktas, A. and others",
+ title = "{Tests of QCD factorisation in the diffractive production
+ of dijets in deep-inelastic scattering and photoproduction
+ at HERA}",
+ collaboration = "H1",
+ journal = "Eur. Phys. J.",
+ volume = "C51",
+ year = "2007",
+ pages = "549-568",
+ doi = "10.1140/epjc/s10052-007-0325-4",
+ eprint = "hep-ex/0703022",
+ archivePrefix = "arXiv",
+ primaryClass = "hep-ex",
+ reportNumber = "DESY-07-018",
+ SLACcitation = "%%CITATION = HEP-EX/0703022;%%"
+}'
+
diff --git a/analyses/pluginHERA/H1_2007_I746380.plot b/analyses/pluginHERA/H1_2007_I746380.plot
new file mode 100644
--- /dev/null
+++ b/analyses/pluginHERA/H1_2007_I746380.plot
@@ -0,0 +1,174 @@
+# BEGIN PLOT /H1_2007_I746380/d01-x01-y01
+Title= DIS
+YLabel=$\mathrm{d}\sigma / \mathrm{d} z_{\mathbb{P}}$
+XLabel=$z_{\mathbb{P}}$
+LogY=0
+LegendXPos=0.1
+LegendYPos=0.35
+RatioPlotYMin=0.1
+RatioPlotYMax=2.2
+# END PLOT
+
+# BEGIN PLOT /H1_2007_I746380/d02-x01-y01
+Title= DIS
+YLabel=$\mathrm{d}\sigma / \mathrm{d} \mathrm{log}_{10}x_{\mathbb{P}}$
+XLabel=$\mathrm{log}_{10}x_{\mathbb{P}}$
+LogY=0
+LegendXPos=0.1
+LegendYPos=0.95
+RatioPlotYMin=0.1
+RatioPlotYMax=2.2
+# END PLOT
+
+# BEGIN PLOT /H1_2007_I746380/d03-x01-y01
+Title= DIS
+YLabel=$\mathrm{d}\sigma / \mathrm{d} W$
+XLabel=$W$
+LogY=0
+LegendXPos=0.1
+LegendYPos=0.35
+RatioPlotYMin=0.1
+RatioPlotYMax=2.2
+# END PLOT
+
+# BEGIN PLOT /H1_2007_I746380/d04-x01-y01
+Title= DIS
+YLabel=$\mathrm{d}\sigma / \mathrm{d} Q^2$
+XLabel=$Q^2$
+LogY=0
+LegendXPos=0.6
+LegendYPos=0.95
+RatioPlotYMin=0.1
+RatioPlotYMax=2.2
+# END PLOT
+
+# BEGIN PLOT /H1_2007_I746380/d05-x01-y01
+Title= DIS
+YLabel=$\mathrm{d}\sigma / \mathrm{d} E_T^{*\,\mathrm{jet\,1}}$
+XLabel=$E_T^{*\,\mathrm{jet\,1}}$
+LogY=0
+LegendXPos=0.6
+LegendYPos=0.95
+RatioPlotYMin=0.1
+RatioPlotYMax=2.2
+# END PLOT
+
+# BEGIN PLOT /H1_2007_I746380/d06-x01-y01
+Title=DIS
+YLabel=$\mathrm{d}\sigma / \mathrm{d}<\eta_{\mathrm{jet}}^{\mathrm{lab}}>$
+XLabel=$<\eta_{\mathrm{jet}}^{\mathrm{lab}}>$
+LogY=0
+LegendXPos=0.1
+LegendYPos=0.95
+RatioPlotYMin=0.1
+RatioPlotYMax=2.2
+# END PLOT
+
+# BEGIN PLOT /H1_2007_I746380/d07-x01-y01
+Title=DIS
+YLabel=$\mathrm{d}\sigma / \mathrm{d}|\Delta\eta_{\mathrm{jet}}^{*}|$
+XLabel=$|\Delta\eta_{\mathrm{jet}}^{*}|$
+LogY=0
+LegendXPos=0.6
+LegendYPos=0.95
+RatioPlotYMin=0.1
+RatioPlotYMax=2.2
+# END PLOT
+
+# BEGIN PLOT /H1_2007_I746380/d08-x01-y01
+Title=Photoproduction
+YLabel=$\mathrm{d}\sigma / \mathrm{d} z_{\mathbb{P}}$
+XLabel=$z_{\mathbb{P}}$
+LogY=0
+LegendXPos=0.1
+LegendYPos=0.35
+RatioPlotYMin=0.1
+RatioPlotYMax=2.2
+# END PLOT
+
+# BEGIN PLOT /H1_2007_I746380/d09-x01-y01
+Title=Photoproduction
+YLabel=$\mathrm{d}\sigma / \mathrm{d} x_{\gamma}$
+XLabel=$x_{\gamma}$
+LogY=0
+LegendXPos=0.1
+LegendYPos=0.95
+RatioPlotYMin=0.1
+RatioPlotYMax=2.2
+# END PLOT
+
+# BEGIN PLOT /H1_2007_I746380/d10-x01-y01
+Title=Photoproduction
+YLabel=$\mathrm{d}\sigma / \mathrm{d} \mathrm{log}_{10}x_{\mathbb{P}}$
+XLabel=$\mathrm{log}_{10}x_{\mathbb{P}}$
+LogY=0
+LegendXPos=0.1
+LegendYPos=0.95
+RatioPlotYMin=0.1
+RatioPlotYMax=2.2
+# END PLOT
+
+# BEGIN PLOT /H1_2007_I746380/d11-x01-y01
+Title= Photoproduction
+YLabel=$\mathrm{d}\sigma / \mathrm{d} W$
+XLabel=$W$
+LogY=0
+LegendXPos=0.1
+LegendYPos=0.35
+RatioPlotYMin=0.1
+RatioPlotYMax=2.2
+# END PLOT
+
+# BEGIN PLOT /H1_2007_I746380/d12-x01-y01
+Title= Photoproduction
+YLabel=$\mathrm{d}\sigma / \mathrm{d} E_T^{*\,\mathrm{jet\,1}}$
+XLabel=$E_T^{*\,\mathrm{jet\,1}}$
+LogY=1
+LegendXPos=0.1
+LegendYPos=0.35
+RatioPlotYMin=0.1
+RatioPlotYMax=2.2
+# END PLOT
+
+# BEGIN PLOT /H1_2007_I746380/d13-x01-y01
+Title= Photoproduction
+YLabel=$\mathrm{d}\sigma / \mathrm{d}<\eta_{\mathrm{jet}}^{\mathrm{lab}}>$
+XLabel=$<\eta_{\mathrm{jet}}^{\mathrm{lab}}>$
+LogY=0
+LegendXPos=0.6
+LegendYPos=0.95
+RatioPlotYMin=0.1
+RatioPlotYMax=2.2
+# END PLOT
+
+# BEGIN PLOT /H1_2007_I746380/d14-x01-y01
+Title=Photoproduction
+YLabel=$\mathrm{d}\sigma / \mathrm{d}|\Delta\eta_{\mathrm{jet}}^{*}|$
+XLabel=$|\Delta\eta_{\mathrm{jet}}^{*}|$
+LogY=0
+LegendXPos=0.6
+LegendYPos=0.95
+RatioPlotYMin=0.1
+RatioPlotYMax=2.2
+# END PLOT
+
+# BEGIN PLOT /H1_2007_I746380/d15-x01-y01
+Title=Photoproduction
+YLabel=$\mathrm{d}\sigma / \mathrm{d}M_{12}$
+XLabel=$M_{12}$
+LogY=0
+LegendXPos=0.6
+LegendYPos=0.95
+RatioPlotYMin=0.1
+RatioPlotYMax=2.2
+# END PLOT
+
+# ... add more histograms as you need them ...
+# BEGIN PLOT /H1_2007_I746380/d01-x01-y01
+#Title=[Uncomment and insert title for histogram d01-x01-y01 here]
+#XLabel=[Uncomment and insert x-axis label for histogram d01-x01-y01 here]
+#YLabel=[Uncomment and insert y-axis label for histogram d01-x01-y01 here]
+# + any additional plot settings you might like, see make-plots documentation
+# END PLOT
+
+# ... add more histograms as you need them ...
diff --git a/analyses/pluginHERA/H1_2007_I746380.yoda b/analyses/pluginHERA/H1_2007_I746380.yoda
new file mode 100644
--- /dev/null
+++ b/analyses/pluginHERA/H1_2007_I746380.yoda
@@ -0,0 +1,176 @@
+BEGIN YODA_SCATTER2D_V2 /REF/H1_2007_I746380/d01-x01-y01
+IsRef: 1
+Path: /REF/H1_2007_I746380/d01-x01-y01
+Title: doi:10.17182/hepdata.45555.v1/t1
+Type: Scatter2D
+---
+# xval xerr- xerr+ yval yerr- yerr+
+2.500000e-01 1.500000e-01 1.500000e-01 5.900000e+01 1.118034e+01 1.118034e+01
+5.000000e-01 1.000000e-01 1.000000e-01 3.400000e+01 8.944272e+00 8.944272e+00
+7.000000e-01 1.000000e-01 1.000000e-01 1.600000e+01 5.000000e+00 5.000000e+00
+9.000000e-01 1.000000e-01 1.000000e-01 5.200000e+00 2.332381e+00 2.332381e+00
+END YODA_SCATTER2D_V2
+BEGIN YODA_SCATTER2D_V2 /REF/H1_2007_I746380/d02-x01-y01
+IsRef: 1
+Path: /REF/H1_2007_I746380/d02-x01-y01
+Title: doi:10.17182/hepdata.45555.v1/t2
+Type: Scatter2D
+---
+# xval xerr- xerr+ yval yerr- yerr+
+-2.100000e+00 2.000000e-01 2.000000e-01 2.030000e+01 5.024938e+00 5.024938e+00
+-1.800000e+00 1.000000e-01 1.000000e-01 4.300000e+01 8.602325e+00 8.602325e+00
+-1.600000e+00 1.000000e-01 1.000000e-01 6.200000e+01 1.476482e+01 1.476482e+01
+END YODA_SCATTER2D_V2
+BEGIN YODA_SCATTER2D_V2 /REF/H1_2007_I746380/d03-x01-y01
+IsRef: 1
+Path: /REF/H1_2007_I746380/d03-x01-y01
+Title: doi:10.17182/hepdata.45555.v1/t3
+Type: Scatter2D
+---
+# xval xerr- xerr+ yval yerr- yerr+
+1.750000e+02 1.000000e+01 1.000000e+01 4.500000e-01 1.166190e-01 1.166190e-01
+1.950000e+02 1.000000e+01 1.000000e+01 4.000000e-01 9.433981e-02 9.433981e-02
+2.150000e+02 1.000000e+01 1.000000e+01 3.700000e-01 8.602325e-02 8.602325e-02
+2.335000e+02 8.500000e+00 8.500000e+00 3.300000e-01 8.602325e-02 8.602325e-02
+END YODA_SCATTER2D_V2
+BEGIN YODA_SCATTER2D_V2 /REF/H1_2007_I746380/d04-x01-y01
+IsRef: 1
+Path: /REF/H1_2007_I746380/d04-x01-y01
+Title: doi:10.17182/hepdata.45555.v1/t4
+Type: Scatter2D
+---
+# xval xerr- xerr+ yval yerr- yerr+
+5.000000e+00 1.000000e+00 1.000000e+00 3.400000e+00 8.602325e-01 8.602325e-01
+7.000000e+00 1.000000e+00 1.000000e+00 2.400000e+00 6.403124e-01 6.403124e-01
+1.000000e+01 2.000000e+00 2.000000e+00 8.300000e-01 2.701851e-01 2.701851e-01
+1.600000e+01 4.000000e+00 4.000000e+00 6.300000e-01 1.720465e-01 1.720465e-01
+2.500000e+01 5.000000e+00 5.000000e+00 4.700000e-01 1.442221e-01 1.442221e-01
+3.500000e+01 5.000000e+00 5.000000e+00 1.800000e-01 7.810250e-02 7.810250e-02
+6.000000e+01 2.000000e+01 2.000000e+01 8.100000e-02 3.342155e-02 3.342155e-02
+END YODA_SCATTER2D_V2
+BEGIN YODA_SCATTER2D_V2 /REF/H1_2007_I746380/d05-x01-y01
+IsRef: 1
+Path: /REF/H1_2007_I746380/d05-x01-y01
+Title: doi:10.17182/hepdata.45555.v1/t5
+Type: Scatter2D
+---
+# xval xerr- xerr+ yval yerr- yerr+
+6.000000e+00 1.000000e+00 1.000000e+00 9.300000e+00 1.700000e+00 1.700000e+00
+8.000000e+00 1.000000e+00 1.000000e+00 4.100000e+00 1.029563e+00 1.029563e+00
+1.000000e+01 1.000000e+00 1.000000e+00 1.000000e+00 4.472136e-01 4.472136e-01
+END YODA_SCATTER2D_V2
+BEGIN YODA_SCATTER2D_V2 /REF/H1_2007_I746380/d06-x01-y01
+IsRef: 1
+Path: /REF/H1_2007_I746380/d06-x01-y01
+Title: doi:10.17182/hepdata.45555.v1/t6
+Type: Scatter2D
+---
+# xval xerr- xerr+ yval yerr- yerr+
+-5.500000e-01 1.500000e-01 1.500000e-01 2.100000e+01 5.830952e+00 5.830952e+00
+-2.500000e-01 1.500000e-01 1.500000e-01 2.900000e+01 7.211103e+00 7.211103e+00
+5.000000e-02 1.500000e-01 1.500000e-01 2.100000e+01 5.830952e+00 5.830952e+00
+END YODA_SCATTER2D_V2
+BEGIN YODA_SCATTER2D_V2 /REF/H1_2007_I746380/d07-x01-y01
+IsRef: 1
+Path: /REF/H1_2007_I746380/d07-x01-y01
+Title: doi:10.17182/hepdata.45555.v1/t7
+Type: Scatter2D
+---
+# xval xerr- xerr+ yval yerr- yerr+
+2.500000e-01 2.500000e-01 2.500000e-01 3.000000e+01 5.830952e+00 5.830952e+00
+7.500000e-01 2.500000e-01 2.500000e-01 1.500000e+01 3.605551e+00 3.605551e+00
+1.500000e+00 5.000000e-01 5.000000e-01 5.900000e+00 1.526434e+00 1.526434e+00
+END YODA_SCATTER2D_V2
+BEGIN YODA_SCATTER2D_V2 /REF/H1_2007_I746380/d08-x01-y01
+IsRef: 1
+Path: /REF/H1_2007_I746380/d08-x01-y01
+Title: doi:10.17182/hepdata.45555.v1/t8
+Type: Scatter2D
+---
+# xval xerr- xerr+ yval yerr- yerr+
+2.500000e-01 1.500000e-01 1.500000e-01 2.900000e+02 5.385165e+01 5.385165e+01
+5.000000e-01 1.000000e-01 1.000000e-01 3.400000e+02 7.280110e+01 7.280110e+01
+7.000000e-01 1.000000e-01 1.000000e-01 3.100000e+02 5.385165e+01 5.385165e+01
+9.000000e-01 1.000000e-01 1.000000e-01 1.500000e+02 3.162278e+01 3.162278e+01
+END YODA_SCATTER2D_V2
+BEGIN YODA_SCATTER2D_V2 /REF/H1_2007_I746380/d09-x01-y01
+IsRef: 1
+Path: /REF/H1_2007_I746380/d09-x01-y01
+Title: doi:10.17182/hepdata.45555.v1/t9
+Type: Scatter2D
+---
+# xval xerr- xerr+ yval yerr- yerr+
+2.500000e-01 1.500000e-01 1.500000e-01 1.500000e+02 3.162278e+01 3.162278e+01
+5.000000e-01 1.000000e-01 1.000000e-01 2.500000e+02 5.385165e+01 5.385165e+01
+7.000000e-01 1.000000e-01 1.000000e-01 3.700000e+02 4.472136e+01 4.472136e+01
+9.000000e-01 1.000000e-01 1.000000e-01 3.500000e+02 7.280110e+01 7.280110e+01
+END YODA_SCATTER2D_V2
+BEGIN YODA_SCATTER2D_V2 /REF/H1_2007_I746380/d10-x01-y01
+IsRef: 1
+Path: /REF/H1_2007_I746380/d10-x01-y01
+Title: doi:10.17182/hepdata.45555.v1/t10
+Type: Scatter2D
+---
+# xval xerr- xerr+ yval yerr- yerr+
+-2.200000e+00 1.000000e-01 1.000000e-01 1.100000e+02 3.162278e+01 3.162278e+01
+-2.000000e+00 1.000000e-01 1.000000e-01 2.000000e+02 3.162278e+01 3.162278e+01
+-1.800000e+00 1.000000e-01 1.000000e-01 3.500000e+02 5.385165e+01 5.385165e+01
+-1.600000e+00 1.000000e-01 1.000000e-01 5.500000e+02 1.044031e+02 1.044031e+02
+END YODA_SCATTER2D_V2
+BEGIN YODA_SCATTER2D_V2 /REF/H1_2007_I746380/d11-x01-y01
+IsRef: 1
+Path: /REF/H1_2007_I746380/d11-x01-y01
+Title: doi:10.17182/hepdata.45555.v1/t11
+Type: Scatter2D
+---
+# xval xerr- xerr+ yval yerr- yerr+
+1.750000e+02 1.000000e+01 1.000000e+01 3.800000e+00 6.324555e-01 6.324555e-01
+1.950000e+02 1.000000e+01 1.000000e+01 3.500000e+00 5.385165e-01 5.385165e-01
+2.150000e+02 1.000000e+01 1.000000e+01 3.000000e+00 4.472136e-01 4.472136e-01
+2.335000e+02 8.500000e+00 8.500000e+00 2.400000e+00 4.472136e-01 4.472136e-01
+END YODA_SCATTER2D_V2
+BEGIN YODA_SCATTER2D_V2 /REF/H1_2007_I746380/d12-x01-y01
+IsRef: 1
+Path: /REF/H1_2007_I746380/d12-x01-y01
+Title: doi:10.17182/hepdata.45555.v1/t12
+Type: Scatter2D
+---
+# xval xerr- xerr+ yval yerr- yerr+
+6.000000e+00 1.000000e+00 1.000000e+00 8.500000e+01 1.236932e+01 1.236932e+01
+8.000000e+00 1.000000e+00 1.000000e+00 2.800000e+01 4.472136e+00 4.472136e+00
+1.000000e+01 1.000000e+00 1.000000e+00 7.300000e+00 2.061553e+00 2.061553e+00
+END YODA_SCATTER2D_V2
+BEGIN YODA_SCATTER2D_V2 /REF/H1_2007_I746380/d13-x01-y01
+IsRef: 1
+Path: /REF/H1_2007_I746380/d13-x01-y01
+Title: doi:10.17182/hepdata.45555.v1/t13
+Type: Scatter2D
+---
+# xval xerr- xerr+ yval yerr- yerr+
+-5.500000e-01 1.500000e-01 1.500000e-01 1.400000e+02 3.162278e+01 3.162278e+01
+-2.500000e-01 1.500000e-01 1.500000e-01 2.300000e+02 3.162278e+01 3.162278e+01
+5.000000e-02 1.500000e-01 1.500000e-01 1.900000e+02 3.162278e+01 3.162278e+01
+5.000000e-01 3.000000e-01 3.000000e-01 8.900000e+01 1.431782e+01 1.431782e+01
+END YODA_SCATTER2D_V2
+BEGIN YODA_SCATTER2D_V2 /REF/H1_2007_I746380/d14-x01-y01
+IsRef: 1
+Path: /REF/H1_2007_I746380/d14-x01-y01
+Title: doi:10.17182/hepdata.45555.v1/t14
+Type: Scatter2D
+---
+# xval xerr- xerr+ yval yerr- yerr+
+2.500000e-01 2.500000e-01 2.500000e-01 1.790000e+02 3.036445e+01 3.036445e+01
+7.500000e-01 2.500000e-01 2.500000e-01 1.570000e+02 2.284732e+01 2.284732e+01
+1.250000e+00 2.500000e-01 2.500000e-01 8.800000e+01 1.431782e+01 1.431782e+01
+1.750000e+00 2.500000e-01 2.500000e-01 5.500000e+01 1.029563e+01 1.029563e+01
+END YODA_SCATTER2D_V2
+BEGIN YODA_SCATTER2D_V2 /REF/H1_2007_I746380/d15-x01-y01
+IsRef: 1
+Path: /REF/H1_2007_I746380/d15-x01-y01
+Title: doi:10.17182/hepdata.45555.v1/t15
+Type: Scatter2D
+---
+# xval xerr- xerr+ yval yerr- yerr+
+1.300000e+01 4.000000e+00 4.000000e+00 2.560000e+01 4.002499e+00 4.002499e+00
+2.200000e+01 5.000000e+00 5.000000e+00 3.800000e+00 7.615773e-01 7.615773e-01
+END YODA_SCATTER2D_V2