diff --git a/analyses/pluginHERA/H1_1994_S2919893.cc b/analyses/pluginHERA/H1_1994_S2919893.cc --- a/analyses/pluginHERA/H1_1994_S2919893.cc +++ b/analyses/pluginHERA/H1_1994_S2919893.cc @@ -1,226 +1,228 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Math/Constants.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/DISKinematics.hh" namespace Rivet { /// @brief H1 energy flow and charged particle spectra /// @author Peter Richardson /// Based on the equivalent HZTool analysis class H1_1994_S2919893 : public Analysis { public: /// Constructor H1_1994_S2919893() : Analysis("H1_1994_S2919893") { // Initialise member variables _w77 = make_pair(0.0, 0.0); _w122 = make_pair(0.0, 0.0); _w169 = make_pair(0.0, 0.0); _w117 = make_pair(0.0, 0.0); _wEnergy = make_pair(0.0, 0.0); } /// @name Analysis methods //@{ /// Initialise projections and histograms void init() { // Projections declare(DISLepton(), "Lepton"); declare(DISKinematics(), "Kinematics"); declare(FinalState(), "FS"); // Histos _histEnergyFlowLowX = bookHisto1D(1, 1, 1); _histEnergyFlowHighX = bookHisto1D(1, 1, 2); _histEECLowX = bookHisto1D(2, 1, 1); _histEECHighX = bookHisto1D(2, 1, 2); _histSpectraW77 = bookHisto1D(3, 1, 1); _histSpectraW122 = bookHisto1D(3, 1, 2); _histSpectraW169 = bookHisto1D(3, 1, 3); _histSpectraW117 = bookHisto1D(3, 1, 4); _histPT2 = bookProfile1D(4, 1, 1); } /// Analyse each event void analyze(const Event& event) { // Get the DIS kinematics const DISKinematics& dk = apply(event, "Kinematics"); + if ( dk.failed() ) vetoEvent; const double x = dk.x(); const double w2 = dk.W2(); const double w = sqrt(w2); // Momentum of the scattered lepton const DISLepton& dl = apply(event,"Lepton"); + if ( dl.failed() ) return; const FourMomentum leptonMom = dl.out(); const double ptel = leptonMom.pT(); const double enel = leptonMom.E(); const double thel = leptonMom.angle(dk.beamHadron().mom())/degree; // Extract the particles other than the lepton const FinalState& fs = apply(event, "FS"); Particles particles; particles.reserve(fs.particles().size()); const GenParticle* dislepGP = dl.out().genParticle(); foreach (const Particle& p, fs.particles()) { const GenParticle* loopGP = p.genParticle(); if (loopGP == dislepGP) continue; particles.push_back(p); } // Cut on the forward energy double efwd = 0.0; foreach (const Particle& p, particles) { const double th = p.angle(dk.beamHadron())/degree; if (inRange(th, 4.4, 15)) efwd += p.E(); } // Apply the cuts // Lepton energy and angle, w2 and forward energy MSG_DEBUG("enel/GeV = " << enel/GeV << ", thel = " << thel << ", w2 = " << w2 << ", efwd/GeV = " << efwd/GeV); bool cut = enel/GeV > 14. && thel > 157. && thel < 172.5 && w2 >= 3000. && efwd/GeV > 0.5; if (!cut) vetoEvent; // Weight of the event const double weight = event.weight(); (x < 1e-3 ? _wEnergy.first : _wEnergy.second) += weight; // Boost to hadronic CM const LorentzTransform hcmboost = dk.boostHCM(); // Loop over the particles long ncharged(0); for (size_t ip1 = 0; ip1 < particles.size(); ++ip1) { const Particle& p = particles[ip1]; const double th = p.angle(dk.beamHadron().momentum()) / degree; // Boost momentum to lab const FourMomentum hcmMom = hcmboost.transform(p.momentum()); // Angular cut if (th <= 4.4) continue; // Energy flow histogram const double et = fabs(hcmMom.Et()); const double eta = hcmMom.eta(); (x < 1e-3 ? _histEnergyFlowLowX : _histEnergyFlowHighX)->fill(eta, et*weight); if (PID::threeCharge(p.pid()) != 0) { /// @todo Use units in w comparisons... what are the units? if (w > 50. && w <= 200.) { double xf= 2 * hcmMom.z() / w; double pt2 = hcmMom.pT2(); if (w > 50. && w <= 100.) { _histSpectraW77 ->fill(xf, weight); } else if (w > 100. && w <= 150.) { _histSpectraW122->fill(xf, weight); } else if (w > 150. && w <= 200.) { _histSpectraW169->fill(xf, weight); } _histSpectraW117->fill(xf, weight); /// @todo Is this profile meant to be filled with 2 weight factors? _histPT2->fill(xf, pt2*weight/GeV2, weight); ++ncharged; } } // Energy-energy correlation if (th <= 8.) continue; double phi1 = p.phi(ZERO_2PI); double eta1 = p.eta(); double et1 = fabs(p.momentum().Et()); for (size_t ip2 = ip1+1; ip2 < particles.size(); ++ip2) { const Particle& p2 = particles[ip2]; //double th2 = beamAngle(p2.momentum(), order); double th2 = p2.angle(dk.beamHadron().momentum()) / degree; if (th2 <= 8.) continue; double phi2 = p2.phi(ZERO_2PI); /// @todo Use angle function double deltaphi = phi1 - phi2; if (fabs(deltaphi) > PI) deltaphi = fabs(fabs(deltaphi) - TWOPI); double eta2 = p2.eta(); double omega = sqrt(sqr(eta1-eta2) + sqr(deltaphi)); double et2 = fabs(p2.momentum().Et()); double wt = et1*et2 / sqr(ptel) * weight; (x < 1e-3 ? _histEECLowX : _histEECHighX)->fill(omega, wt); } } // Factors for normalization if (w > 50. && w <= 200.) { if (w <= 100.) { _w77.first += ncharged*weight; _w77.second += weight; } else if (w <= 150.) { _w122.first += ncharged*weight; _w122.second += weight; } else { _w169.first += ncharged*weight; _w169.second += weight; } _w117.first += ncharged*weight; _w117.second += weight; } } // Normalize inclusive single particle distributions to the average number of charged particles per event. void finalize() { normalize(_histSpectraW77, _w77.first/_w77.second); normalize(_histSpectraW122, _w122.first/_w122.second); normalize(_histSpectraW169, _w169.first/_w169.second); normalize(_histSpectraW117, _w117.first/_w117.second); scale(_histEnergyFlowLowX , 1./_wEnergy.first ); scale(_histEnergyFlowHighX, 1./_wEnergy.second); scale(_histEECLowX , 1./_wEnergy.first ); scale(_histEECHighX, 1./_wEnergy.second); } //@} private: /// Polar angle with right direction of the beam inline double beamAngle(const FourVector& v, bool order) { double thel = v.polarAngle()/degree; if (thel < 0) thel += 180.; if (!order) thel = 180 - thel; return thel; } /// @name Histograms //@{ Histo1DPtr _histEnergyFlowLowX, _histEnergyFlowHighX; Histo1DPtr _histEECLowX, _histEECHighX; Histo1DPtr _histSpectraW77, _histSpectraW122, _histSpectraW169, _histSpectraW117; Profile1DPtr _histPT2; //@} /// @name Storage of weights to calculate averages for normalisation //@{ pair _w77, _w122, _w169, _w117, _wEnergy; //@} }; DECLARE_RIVET_PLUGIN(H1_1994_S2919893); } diff --git a/analyses/pluginHERA/H1_1995_S3167097.cc b/analyses/pluginHERA/H1_1995_S3167097.cc --- a/analyses/pluginHERA/H1_1995_S3167097.cc +++ b/analyses/pluginHERA/H1_1995_S3167097.cc @@ -1,120 +1,122 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/DISFinalState.hh" #include "Rivet/Projections/CentralEtHCM.hh" namespace Rivet { /// H1 energy flow in DIS /// /// @todo Make histograms match those in HepData and use autobooking /// /// @author Leif Lonnblad /// @author Andy Buckley class H1_1995_S3167097 : public Analysis { public: /// Constructor H1_1995_S3167097() : Analysis("H1_1995_S3167097") { } /// @name Analysis methods //@{ void init() { // Projections const DISKinematics& diskin = declare(DISKinematics(), "Kinematics"); const DISFinalState& fshcm = declare(DISFinalState(diskin, DISFinalState::HCM), "FS"); declare(CentralEtHCM(fshcm), "Y1HCM"); // Histograms /// @todo Convert to use autobooking and correspond to HepData data tables _sumw.resize(9); _hEtFlow.resize(9); for (size_t i = 0; i < 9; ++i) _hEtFlow[i] = bookHisto1D(to_str(i), 24, -6, 6); _tmphAvEt = Histo1D(9, 1.0, 10.0); _tmphAvX = Histo1D(9, 1.0, 10.0); _tmphAvQ2 = Histo1D(9, 1.0, 10.0); _tmphN = Histo1D(9, 1.0, 10.0); } /// Calculate the bin number from the DISKinematics projection /// @todo Convert to use a HEPUtils Binning1D size_t _getbin(const DISKinematics& dk) { if (inRange(dk.Q2()/GeV2, 5.0, 10.0)) { if (inRange(dk.x(), 1e-4, 2e-4)) return 0; if (inRange(dk.x(), 2e-4, 5e-4) && dk.Q2() > 6.0*GeV2) return 1; } else if (inRange(dk.Q2()/GeV2, 10.0, 20.0)) { if (inRange(dk.x(), 2e-4, 5e-4)) return 2; if (inRange(dk.x(), 5e-4, 8e-4)) return 3; if (inRange(dk.x(), 8e-4, 1.5e-3)) return 4; if (inRange(dk.x(), 1.5e-3, 4e-3)) return 5; } else if (inRange(dk.Q2()/GeV2, 20.0, 50.0)) { if (inRange(dk.x(), 5e-4, 1.4e-3)) return 6; if (inRange(dk.x(), 1.4e-3, 3e-3)) return 7; if (inRange(dk.x(), 3e-3, 1e-2)) return 8; } return -1; } void analyze(const Event& event) { const FinalState& fs = apply(event, "FS"); const DISKinematics& dk = apply(event, "Kinematics"); + if ( dk.failed() ) vetoEvent; const CentralEtHCM& y1 = apply(event, "Y1HCM"); + if ( y1.failed() ) vetoEvent; const int ibin = _getbin(dk); if (ibin < 0) vetoEvent; const double weight = event.weight(); _sumw[ibin] += weight; for (size_t i = 0, N = fs.particles().size(); i < N; ++i) { const double rap = fs.particles()[i].rapidity(); const double et = fs.particles()[i].Et(); _hEtFlow[ibin]->fill(rap, weight * et/GeV); } /// @todo Use fillBin? _tmphAvEt.fill(ibin + 1.5, weight * y1.sumEt()/GeV); _tmphAvX.fill(ibin + 1.5, weight * dk.x()); _tmphAvQ2.fill(ibin + 1.5, weight * dk.Q2()/GeV2); _tmphN.fill(ibin + 1.5, weight); } void finalize() { for (size_t ibin = 0; ibin < 9; ++ibin) scale(_hEtFlow[ibin], 0.5/_sumw[ibin]); /// @todo Improve this! addAnalysisObject(make_shared(_tmphAvEt/_tmphN, histoPath("21")) ); addAnalysisObject(make_shared(_tmphAvX/_tmphN, histoPath("22")) ); addAnalysisObject(make_shared(_tmphAvQ2/_tmphN, histoPath("23")) ); } //@} private: /// Histograms for the \f$ E_T \f$ flow vector _hEtFlow; /// Temporary histograms for averages in different kinematical bins. Histo1D _tmphAvEt, _tmphAvX, _tmphAvQ2, _tmphN; /// Weights counters for each kinematic bin vector _sumw; }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(H1_1995_S3167097); } diff --git a/analyses/pluginHERA/H1_2000_S4129130.cc b/analyses/pluginHERA/H1_2000_S4129130.cc --- a/analyses/pluginHERA/H1_2000_S4129130.cc +++ b/analyses/pluginHERA/H1_2000_S4129130.cc @@ -1,259 +1,261 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Math/Constants.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/DISKinematics.hh" namespace Rivet { /// @brief H1 energy flow and charged particle spectra /// /// @author Peter Richardson /// /// Based on the HZTOOL analysis HZ99091 class H1_2000_S4129130 : public Analysis { public: /// Constructor H1_2000_S4129130() : Analysis("H1_2000_S4129130") { } /// @name Analysis methods //@{ /// Initialise projections and histograms void init() { // Projections declare(DISLepton(), "Lepton"); declare(DISKinematics(), "Kinematics"); declare(FinalState(), "FS"); // Histos Histo1DPtr h; // Histograms and weight vectors for low Q^2 a for (size_t ix = 0; ix < 17; ++ix) { h = bookHisto1D(ix+1, 1, 1); _histETLowQa.push_back(h); _weightETLowQa.push_back(0.); } // Histograms and weight vectors for high Q^2 a for (size_t ix = 0; ix < 7; ++ix) { h = bookHisto1D(ix+18, 1, 1); _histETHighQa.push_back(h); _weightETHighQa.push_back(0.); } // Histograms and weight vectors for low Q^2 b for (size_t ix = 0; ix < 5; ++ix) { h = bookHisto1D(ix+25, 1, 1); _histETLowQb.push_back(h); _weightETLowQb.push_back(0.); } // Histograms and weight vectors for high Q^2 b for (size_t ix = 0; ix < 3; ++ix) { h = bookHisto1D(30+ix, 1, 1); _histETHighQb.push_back(h); _weightETHighQb.push_back(0.0); } // Histograms for the averages _histAverETCentral = bookProfile1D(33, 1, 1); _histAverETFrag = bookProfile1D(34, 1, 1); } /// Analyze each event void analyze(const Event& event) { // DIS kinematics const DISKinematics& dk = apply(event, "Kinematics"); + if ( dk.failed() ) vetoEvent; double q2 = dk.Q2(); double x = dk.x(); double y = dk.y(); double w2 = dk.W2(); // Kinematics of the scattered lepton const DISLepton& dl = apply(event,"Lepton"); + if ( dl.failed() ) return; const FourMomentum leptonMom = dl.out(); const double enel = leptonMom.E(); const double thel = 180 - leptonMom.angle(dl.in().mom())/degree; // Extract the particles other than the lepton const FinalState& fs = apply(event, "FS"); Particles particles; particles.reserve(fs.size()); const GenParticle* dislepGP = dl.out().genParticle(); ///< @todo Is the GenParticle stuff necessary? (Not included in Particle::==?) foreach (const Particle& p, fs.particles()) { const GenParticle* loopGP = p.genParticle(); if (loopGP == dislepGP) continue; particles.push_back(p); } // Cut on the forward energy double efwd = 0.; foreach (const Particle& p, particles) { const double th = 180 - p.angle(dl.in())/degree; if (inRange(th, 4.4, 15.0)) efwd += p.E(); } // There are four possible selections for events bool evcut[4]; // Low Q2 selection a evcut[0] = enel/GeV > 12. && w2 >= 4400.*GeV2 && efwd/GeV > 0.5 && inRange(thel,157.,176.); // Low Q2 selection b evcut[1] = enel/GeV > 12. && inRange(y,0.3,0.5); // High Q2 selection a evcut[2] = inRange(thel,12.,150.) && inRange(y,0.05,0.6) && w2 >= 4400.*GeV2 && efwd > 0.5; // High Q2 selection b evcut[3] = inRange(thel,12.,150.) && inRange(y,0.05,0.6) && inRange(w2,27110.*GeV2,45182.*GeV2); // Veto if fails all cuts /// @todo Can we use all()? if (! (evcut[0] || evcut[1] || evcut[2] || evcut[3]) ) vetoEvent; // Find the bins int bin[4] = {-1,-1,-1,-1}; // For the low Q2 selection a) if (q2 > 2.5*GeV && q2 <= 5.*GeV) { if (x > 0.00005 && x <= 0.0001 ) bin[0] = 0; if (x > 0.0001 && x <= 0.0002 ) bin[0] = 1; if (x > 0.0002 && x <= 0.00035) bin[0] = 2; if (x > 0.00035 && x <= 0.0010 ) bin[0] = 3; } else if (q2 > 5.*GeV && q2 <= 10.*GeV) { if (x > 0.0001 && x <= 0.0002 ) bin[0] = 4; if (x > 0.0002 && x <= 0.00035) bin[0] = 5; if (x > 0.00035 && x <= 0.0007 ) bin[0] = 6; if (x > 0.0007 && x <= 0.0020 ) bin[0] = 7; } else if (q2 > 10.*GeV && q2 <= 20.*GeV) { if (x > 0.0002 && x <= 0.0005) bin[0] = 8; if (x > 0.0005 && x <= 0.0008) bin[0] = 9; if (x > 0.0008 && x <= 0.0015) bin[0] = 10; if (x > 0.0015 && x <= 0.040 ) bin[0] = 11; } else if (q2 > 20.*GeV && q2 <= 50.*GeV) { if (x > 0.0005 && x <= 0.0014) bin[0] = 12; if (x > 0.0014 && x <= 0.0030) bin[0] = 13; if (x > 0.0030 && x <= 0.0100) bin[0] = 14; } else if (q2 > 50.*GeV && q2 <= 100.*GeV) { if (x >0.0008 && x <= 0.0030) bin[0] = 15; if (x >0.0030 && x <= 0.0200) bin[0] = 16; } // check in one of the bins evcut[0] &= bin[0] >= 0; // For the low Q2 selection b) if (q2 > 2.5*GeV && q2 <= 5. *GeV) bin[1] = 0; if (q2 > 5. *GeV && q2 <= 10. *GeV) bin[1] = 1; if (q2 > 10.*GeV && q2 <= 20. *GeV) bin[1] = 2; if (q2 > 20.*GeV && q2 <= 50. *GeV) bin[1] = 3; if (q2 > 50.*GeV && q2 <= 100.*GeV) bin[1] = 4; // check in one of the bins evcut[1] &= bin[1] >= 0; // for the high Q2 selection a) if (q2 > 100.*GeV && q2 <= 400.*GeV) { if (x > 0.00251 && x <= 0.00631) bin[2] = 0; if (x > 0.00631 && x <= 0.0158 ) bin[2] = 1; if (x > 0.0158 && x <= 0.0398 ) bin[2] = 2; } else if (q2 > 400.*GeV && q2 <= 1100.*GeV) { if (x > 0.00631 && x <= 0.0158 ) bin[2] = 3; if (x > 0.0158 && x <= 0.0398 ) bin[2] = 4; if (x > 0.0398 && x <= 1. ) bin[2] = 5; } else if (q2 > 1100.*GeV && q2 <= 100000.*GeV) { if (x > 0. && x <= 1.) bin[2] = 6; } // check in one of the bins evcut[2] &= bin[2] >= 0; // for the high Q2 selection b) if (q2 > 100.*GeV && q2 <= 220.*GeV) bin[3] = 0; else if (q2 > 220.*GeV && q2 <= 400.*GeV) bin[3] = 1; else if (q2 > 400. ) bin[3] = 2; // check in one of*GeV the bins evcut[3] &= bin[3] >= 0; // Veto if fails all cuts after bin selection /// @todo Can we use all()? if (! (evcut[0] || evcut[1] || evcut[2] || evcut[3])) vetoEvent; // Increment the count for normalisation const double weight = event.weight(); if (evcut[0]) _weightETLowQa [bin[0]] += weight; if (evcut[1]) _weightETLowQb [bin[1]] += weight; if (evcut[2]) _weightETHighQa[bin[2]] += weight; if (evcut[3]) _weightETHighQb[bin[3]] += weight; // Boost to hadronic CoM const LorentzTransform hcmboost = dk.boostHCM(); // Loop over the particles double etcent = 0; double etfrag = 0; foreach (const Particle& p, particles) { // Boost momentum to CMS const FourMomentum hcmMom = hcmboost.transform(p.momentum()); double et = fabs(hcmMom.Et()); double eta = hcmMom.eta(); // Averages in central and forward region if (fabs(eta) < .5 ) etcent += et; if (eta > 2 && eta <= 3.) etfrag += et; // Histograms of Et flow if (evcut[0]) _histETLowQa [bin[0]]->fill(eta, et*weight); if (evcut[1]) _histETLowQb [bin[1]]->fill(eta, et*weight); if (evcut[2]) _histETHighQa[bin[2]]->fill(eta, et*weight); if (evcut[3]) _histETHighQb[bin[3]]->fill(eta, et*weight); } // Fill histograms for the average quantities if (evcut[1] || evcut[3]) { _histAverETCentral->fill(q2, etcent, weight); _histAverETFrag ->fill(q2, etfrag, weight); } } // Finalize void finalize() { // Normalization of the Et distributions /// @todo Simplify by using normalize() instead? Are all these being normalized to area=1? for (size_t ix = 0; ix < 17; ++ix) if (_weightETLowQa[ix] != 0) scale(_histETLowQa[ix], 1/_weightETLowQa[ix]); for (size_t ix = 0; ix < 7; ++ix) if (_weightETHighQa[ix] != 0) scale(_histETHighQa[ix], 1/_weightETHighQa[ix]); for (size_t ix = 0; ix < 5; ++ix) if (_weightETLowQb[ix] != 0) scale(_histETLowQb[ix], 1/_weightETLowQb[ix]); for (size_t ix = 0; ix < 3; ++ix) if (_weightETHighQb[ix] != 0) scale(_histETHighQb[ix], 1/_weightETHighQb[ix]); } //@} private: /// @name Histograms //@{ vector _histETLowQa; vector _histETHighQa; vector _histETLowQb; vector _histETHighQb; Profile1DPtr _histAverETCentral; Profile1DPtr _histAverETFrag; //@} /// @name storage of weights for normalisation //@{ vector _weightETLowQa; vector _weightETHighQa; vector _weightETLowQb; vector _weightETHighQb; //@} }; DECLARE_RIVET_PLUGIN(H1_2000_S4129130); } diff --git a/analyses/pluginHERA/H1_2007_I746380.cc b/analyses/pluginHERA/H1_2007_I746380.cc --- a/analyses/pluginHERA/H1_2007_I746380.cc +++ b/analyses/pluginHERA/H1_2007_I746380.cc @@ -1,608 +1,615 @@ // -*- 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(p); return mkNamedPCmp(other, "DISKIN") || mkNamedPCmp(other, "DISFS"); } virtual void project(const Event& e){ const DISKinematics& dk = apply(e, "DISKIN"); + if ( dk.failed() ) { + fail(); + return; + } const Particles& p = apply(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(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(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(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(e, "RAPGAP"); + if ( rg.failed() ) vetoEvent; // 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(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(event, "RapidityGap"); + if ( rg.failed() ) vetoEvent; const DISKinematics& kin = apply(event, "Kinematics"); + if ( kin.failed() ) vetoEvent; const BoostedXSystem& disfsXcm = apply( 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(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/ZEUS_2001_S4815815.cc b/analyses/pluginHERA/ZEUS_2001_S4815815.cc --- a/analyses/pluginHERA/ZEUS_2001_S4815815.cc +++ b/analyses/pluginHERA/ZEUS_2001_S4815815.cc @@ -1,172 +1,173 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/Beam.hh" #include "Rivet/Projections/DISKinematics.hh" #include "Rivet/Projections/FastJets.hh" namespace Rivet { /// @brief ZEUS dijet photoproduction study used in the ZEUS jets PDF fit /// /// This class is a reproduction of the HZTool routine for the ZEUS /// dijet photoproduction paper which was used in the ZEUS jets PDF fit. /// /// @note Cleaning cuts on event pT/sqrt(Et) and y_e are not needed in MC analysis. /// /// @author Andy Buckley /// @author Ilkka Helenius class ZEUS_2001_S4815815 : public Analysis { public: /// Constructor DEFAULT_RIVET_ANALYSIS_CTOR(ZEUS_2001_S4815815); /// @name Analysis methods //@{ // Book projections and histograms void init() { // Projections /// @todo Acceptance FinalState fs; declare(FastJets(fs, FastJets::KT, 1.0), "Jets"); //< R=1 checked with Matt Wing declare(DISKinematics(), "Kinematics"); // Table 1 _h_costh[0] = bookHisto1D(1, 1, 1); _h_costh[1] = bookHisto1D(1, 1, 2); // Table 2 _h_etjet1[1][0] = bookHisto1D(2, 1, 1); _h_etjet1[1][1] = bookHisto1D(3, 1, 1); _h_etjet1[1][2] = bookHisto1D(4, 1, 1); _h_etjet1[1][3] = bookHisto1D(5, 1, 1); _h_etjet1[1][4] = bookHisto1D(6, 1, 1); _h_etjet1[1][5] = bookHisto1D(7, 1, 1); // Table 3 _h_etjet1[0][0] = bookHisto1D(8, 1, 1); _h_etjet1[0][1] = bookHisto1D(9, 1, 1); _h_etjet1[0][2] = bookHisto1D(10, 1, 1); _h_etjet1[0][3] = bookHisto1D(11, 1, 1); _h_etjet1[0][4] = bookHisto1D(12, 1, 1); _h_etjet1[0][5] = bookHisto1D(13, 1, 1); // Table 4 _h_etajet2[1][0] = bookHisto1D(14, 1, 1); _h_etajet2[1][1] = bookHisto1D(15, 1, 1); _h_etajet2[1][2] = bookHisto1D(16, 1, 1); // Table 5 _h_etajet2[0][0] = bookHisto1D(17, 1, 1); _h_etajet2[0][1] = bookHisto1D(18, 1, 1); _h_etajet2[0][2] = bookHisto1D(19, 1, 1); // Table 6 _h_xobsy[0] = bookHisto1D(20, 1, 1); _h_xobsy[1] = bookHisto1D(21, 1, 1); _h_xobsy[2] = bookHisto1D(22, 1, 1); _h_xobsy[3] = bookHisto1D(23, 1, 1); } // Do the analysis void analyze(const Event& event) { // Determine kinematics, including event orientation since ZEUS coord system is for +z = proton direction const DISKinematics& kin = apply(event, "Kinematics"); + if ( kin.failed() ) vetoEvent; const int orientation = kin.orientation(); // Q2 and inelasticity cuts if (kin.Q2() > 1*GeV2) vetoEvent; if (!inRange(kin.y(), 0.2, 0.85)) vetoEvent; // Jet selection const Jets jets = apply(event, "Jets") \ .jets(Cuts::Et > 11*GeV && Cuts::etaIn(-1*orientation, 2.4*orientation), cmpMomByEt); MSG_DEBUG("Jet multiplicity = " << jets.size()); if (jets.size() < 2) vetoEvent; const Jet& j1 = jets[0]; const Jet& j2 = jets[1]; if (j1.Et() < 14*GeV) vetoEvent; // Jet eta and cos(theta*) computation const double eta1 = orientation*j1.eta(), eta2 = orientation*j2.eta(); const double etabar = (eta1 + eta2)/2; const double etadiff = eta1 - eta2; const double costhetastar = tanh(etadiff/2); // Computation of x_y^obs /// @note Assuming Ee is the lab frame positron momentum, not in proton rest frame cf. the ambiguous phrase in the paper const double xyobs = (j1.Et() * exp(-eta1) + j2.Et() * exp(-eta2)) / (2*kin.y()*kin.beamLepton().E()); const size_t i_xyobs = (xyobs < 0.75) ? 0 : 1; // Calculate the invariant mass of the dijet as in the paper const double mjj = sqrt( 2.*j1.Et()*j2.Et()*( cosh(j1.eta() - j2.eta()) - cos(j1.phi() - j2.phi()) ) ); // Fill histograms const double weight = event.weight(); // T1 if (mjj > 42*GeV && inRange(etabar, 0.1, 1.3)) _h_costh[i_xyobs]->fill(abs(costhetastar), weight); // T2, T3: Symmetrize eta selection, each event contribute twice to the cross section for (size_t isel = 0; isel < 2; ++isel) { double etaJet1 = (isel == 0) ? orientation*j1.eta() : orientation*j2.eta(); double etaJet2 = (isel == 0) ? orientation*j2.eta() : orientation*j1.eta(); if (inRange(etaJet1, -1, 0) && inRange(etaJet2, -1, 0)) _h_etjet1[i_xyobs][0]->fill(j1.Et()/GeV, weight); else if (inRange(etaJet1, 0, 1) && inRange(etaJet2, -1, 0)) _h_etjet1[i_xyobs][1]->fill(j1.Et()/GeV, weight); else if (inRange(etaJet1, 0, 1) && inRange(etaJet2, 0, 1)) _h_etjet1[i_xyobs][2]->fill(j1.Et()/GeV, weight); else if (inRange(etaJet1, 1, 2.4) && inRange(etaJet2, -1, 0)) _h_etjet1[i_xyobs][3]->fill(j1.Et()/GeV, weight); else if (inRange(etaJet1, 1, 2.4) && inRange(etaJet2, 0, 1)) _h_etjet1[i_xyobs][4]->fill(j1.Et()/GeV, weight); else if (inRange(etaJet1, 1, 2.4) && inRange(etaJet2, 1, 2.4)) _h_etjet1[i_xyobs][5]->fill(j1.Et()/GeV, weight); // T4, T5 if (inRange(etaJet1, -1, 0)) _h_etajet2[i_xyobs][0]->fill(etaJet2, weight); else if (inRange(etaJet1, 0, 1)) _h_etajet2[i_xyobs][1]->fill(etaJet2, weight); else if (inRange(etaJet1, 1, 2.4)) _h_etajet2[i_xyobs][2]->fill(etaJet2, weight); } // T6 if (inRange(j1.Et()/GeV, 14, 17)) _h_xobsy[0]->fill(xyobs, weight); else if (inRange(j1.Et()/GeV, 17, 25)) _h_xobsy[1]->fill(xyobs, weight); else if (inRange(j1.Et()/GeV, 25, 35)) _h_xobsy[2]->fill(xyobs, weight); else if (inRange(j1.Et()/GeV, 35, 90)) _h_xobsy[3]->fill(xyobs, weight); } // Finalize void finalize() { const double sf = crossSection()/picobarn/sumOfWeights(); for (size_t ix = 0; ix < 2; ++ix) { scale(_h_costh[ix], sf); for (auto& h : _h_etjet1[ix]) scale(h, sf); for (auto& h : _h_etajet2[ix]) scale(h, sf); } for (auto& h : _h_xobsy) scale(h, sf); } //@} private: /// @name Histograms //@{ Histo1DPtr _h_costh[2], _h_etjet1[2][6], _h_etajet2[2][3], _h_xobsy[4]; //@} }; DECLARE_RIVET_PLUGIN(ZEUS_2001_S4815815); } diff --git a/analyses/pluginMC/MC_DIS_Check.cc b/analyses/pluginMC/MC_DIS_Check.cc new file mode 100644 --- /dev/null +++ b/analyses/pluginMC/MC_DIS_Check.cc @@ -0,0 +1,82 @@ +// -*- C++ -*- +#include "Rivet/Analysis.hh" +#include "Rivet/Projections/FinalState.hh" +#include "Rivet/Projections/FastJets.hh" +#include "Rivet/Projections/DISKinematics.hh" + +namespace Rivet { + + + /// @brief A simple analysis to illustrate how to use the + /// DISKinematics projection together with different options. + class MC_DIS_Check : public Analysis { + public: + + /// Constructor + DEFAULT_RIVET_ANALYSIS_CTOR(MC_DIS_Check); + + + /// @name Analysis methods + //@{ + + /// Book histograms and initialise projections before the run + void init() { + + // Initialise and register projections. Note that the definition + // of the scattered lepton can be influenced by sepcifying + // options as declared in the .info file. + DISLepton lepton(options()); + declare(lepton, "Lepton"); + declare(DISKinematics(lepton), "Kinematics"); + + // Book histograms + + _hist_Q2 = bookHisto1D("Q2",logspace(100,0.1, 1000.0)); + _hist_y = bookHisto1D("y",100,0.,1.); + _hist_x = bookHisto1D("xBj",logspace(100,0.00001, 1.0)); + + } + + + /// Perform the per-event analysis + void analyze(const Event& event) { + // Get the DIS kinematics + const DISKinematics& dk = apply(event, "Kinematics"); + if ( dk.failed() ) return; + double x = dk.x(); + double y = dk.y(); + double Q2 = dk.Q2(); + // Weight of the event + const double weight = event.weight(); + _hist_Q2->fill(Q2,weight); + _hist_y->fill(y,weight); + _hist_x->fill(x,weight); + + /// @todo Do the event by event analysis here + + } + + + /// Normalise histograms etc., after the run + void finalize() { + + normalize(_hist_Q2); // normalize to unity + normalize(_hist_y); // normalize to unity + + } + + //@} + + + /// The histograms. + Histo1DPtr _hist_Q2, _hist_y, _hist_x; + + + }; + + + // The hook for the plugin system + DECLARE_RIVET_PLUGIN(MC_DIS_Check); + + +} diff --git a/analyses/pluginMC/MC_DIS_Check.info b/analyses/pluginMC/MC_DIS_Check.info new file mode 100644 --- /dev/null +++ b/analyses/pluginMC/MC_DIS_Check.info @@ -0,0 +1,27 @@ +Name: MC_DIS_Check + +Summary: A simple analysis using the DISKinematics projection. + +Collider: HERA + +Status: UNVALIDATED +Reentrant: True + +Authors: + - Hannes Jung + - Leif Lönnblad + +Options: + - LMode=prompt,any,dressed + - LSort=ENERGY,ETA,ET + - DressDR=# + - IsolDR=# + - Undress=# + +Beams: [[p+, e-], [p+, e+]] + +Description: + A simple analysis to illustrate how to use the DISKinematics + projection together with different options, and to histogram the + obtained x, y and Q2 variables. + diff --git a/analyses/pluginMC/MC_DIS_Check.plot b/analyses/pluginMC/MC_DIS_Check.plot new file mode 100644 --- /dev/null +++ b/analyses/pluginMC/MC_DIS_Check.plot @@ -0,0 +1,11 @@ +BEGIN PLOT /MC_DIS_Check/Q2 +LogX=1 +END PLOT + +BEGIN PLOT /MC_DIS_Check/xBj +LogX=1 +END PLOT + +BEGIN PLOT /MC_DIS_Check/y +END PLOT + diff --git a/include/Rivet/Event.hh b/include/Rivet/Event.hh --- a/include/Rivet/Event.hh +++ b/include/Rivet/Event.hh @@ -1,188 +1,189 @@ // -*- C++ -*- #ifndef RIVET_Event_HH #define RIVET_Event_HH #include "Rivet/Config/RivetCommon.hh" #include "Rivet/Particle.hh" #include "Rivet/Projection.hh" namespace Rivet { /// Rivet wrapper for HepMC event and Projection references. /// /// Event is a concrete class representing an generated event in Rivet. It is /// constructed given a HepMC::GenEvent, a pointer to which is kept by the /// Event object throughout its lifetime. The user must therefore make sure /// that the corresponding HepMC::GenEvent will persist at least as long as /// the Event object. /// /// In addition to the HepMC::GenEvent object the Event also keeps track of /// all Projection objects which have been applied to the Event so far. class Event { public: /// @name Constructors and destructors. //@{ /// Constructor from a HepMC GenEvent pointer Event(const GenEvent* ge) : _genevent_original(ge), _genevent(*ge) { assert(ge); _init(*ge); } /// Constructor from a HepMC GenEvent reference /// @deprecated HepMC uses pointers, so we should talk to HepMC via pointers Event(const GenEvent& ge) : _genevent_original(&ge), _genevent(ge) { _init(ge); } /// Copy constructor Event(const Event& e) : _genevent_original(e._genevent_original), _genevent(e._genevent) { } //@} /// @name Major event properties //@{ /// The generated event obtained from an external event generator const GenEvent* genEvent() const { return &_genevent; } /// @brief The generation weight associated with the event /// /// @todo This needs to be revisited when we finally add the mechanism to /// support NLO counter-events and weight vectors. double weight() const; /// Get the beam particles ParticlePair beams() const; /// Get the beam centre-of-mass energy double sqrtS() const; /// Get the beam centre-of-mass energy per nucleon double asqrtS() const; /// Get the generator centrality (impact-parameter quantile in [0,1]; or -1 if undefined (usual for non-HI generators)) double centrality() const; // /// Get the boost to the beam centre-of-mass // Vector3 beamCMSBoost() const; // /// Get the boost to the beam centre-of-mass // LorentzTransform beamCMSTransform(); //@} /// @name Access to event particles //@{ /// All the raw GenEvent particles, wrapped in Rivet::Particle objects const Particles& allParticles() const; /// @brief All the raw GenEvent particles, wrapped in Rivet::Particle objects, but with a Cut applied /// /// @note Due to the cut, this returns by value, i.e. involves an expensive copy inline Particles allParticles(const Cut& c) const { return filter_select(allParticles(), c); } /// @brief All the raw GenEvent particles, wrapped in Rivet::Particle objects, but with a selection function applied /// /// @note Due to the cut, this returns by value, i.e. involves an expensive copy template inline Particles allParticles(const FN& f) const { return filter_select(allParticles(), f); } //@} /// @name Projection running //@{ /// @brief Add a projection @a p to this Event. /// /// If an equivalent Projection has been applied before, the /// Projection::project(const Event&) of @a p is not called and a reference /// to the previous equivalent projection is returned. If no previous /// Projection was found, the Projection::project(const Event&) of @a p is /// called and a reference to @a p is returned. template const PROJ& applyProjection(PROJ& p) const { Log& log = Log::getLog("Rivet.Event"); log << Log::TRACE << "Applying projection " << &p << " (" << p.name() << ") -> comparing to projections " << _projections << endl; // First search for this projection *or an equivalent* in the already-executed list const Projection* cpp(&p); std::set::const_iterator old = _projections.find(cpp); if (old != _projections.end()) { log << Log::TRACE << "Equivalent projection found -> returning already-run projection " << *old << endl; const Projection& pRef = **old; return pcast(pRef); } // If this one hasn't been run yet on this event, run it and add to the list log << Log::TRACE << "No equivalent projection in the already-run list -> projecting now" << endl; Projection* pp = const_cast(cpp); + pp->_isValid = true; pp->project(*this); _projections.insert(pp); return p; } /// @brief Add a projection @a p to this Event by pointer. template const PROJ& applyProjection(PROJ* pp) const { if (!pp) throw Error("Event::applyProjection(PROJ*): Projection pointer is null."); return applyProjection(*pp); } //@} private: /// @brief Actual (shared) implementation of the constructors from GenEvents void _init(const GenEvent& ge); // /// @brief Convert the GenEvent to use conventional alignment // /// // /// For example, FHerwig only produces DIS events in the unconventional // /// hadron-lepton orientation and has to be corrected for DIS analysis // /// portability. // void _geNormAlignment(); /// @brief The generated event, as obtained from an external generator. /// /// This is the original GenEvent. In practise the version seen by users /// will often/always be a modified one. /// /// @todo Provide access to this via an Event::originalGenEvent() method? If requested... const GenEvent* _genevent_original; /// @brief The GenEvent used by Rivet analysis projections etc. /// /// This version may be rotated to a "normal" alignment, have /// generator-specific particles stripped out, etc. If an analysis is /// affected by these modifications, it is probably an unphysical analysis! /// /// Stored as a non-pointer since it may get overwritten, and memory for /// copying and cleanup is neater this way. /// @todo Change needed for HepMC3? mutable GenEvent _genevent; /// All the GenEvent particles, wrapped as Rivet::Particles /// @note To be populated lazily, hence mutability mutable Particles _particles; /// The set of Projection objects applied so far mutable std::set _projections; }; } #endif diff --git a/include/Rivet/Makefile.am b/include/Rivet/Makefile.am --- a/include/Rivet/Makefile.am +++ b/include/Rivet/Makefile.am @@ -1,197 +1,200 @@ ## Internal headers - not to be installed nobase_dist_noinst_HEADERS = ## Public headers - to be installed nobase_pkginclude_HEADERS = ## Rivet interface nobase_pkginclude_HEADERS += \ Rivet.hh \ Run.hh \ Event.hh \ ParticleBase.hh \ Particle.fhh Particle.hh \ Jet.fhh Jet.hh \ Projection.fhh Projection.hh \ ProjectionApplier.hh \ ProjectionHandler.hh \ Analysis.hh \ AnalysisHandler.hh \ AnalysisInfo.hh \ AnalysisBuilder.hh \ AnalysisLoader.hh ## Build config stuff nobase_pkginclude_HEADERS += \ Config/RivetCommon.hh \ Config/RivetConfig.hh \ Config/BuildOptions.hh ## Projections nobase_pkginclude_HEADERS += \ Projections/AliceCommon.hh \ Projections/AxesDefinition.hh \ Projections/Beam.hh \ Projections/BeamThrust.hh \ Projections/CentralEtHCM.hh \ Projections/CentralityProjection.hh \ Projections/ChargedFinalState.hh \ Projections/ChargedLeptons.hh \ Projections/ConstLossyFinalState.hh \ Projections/DirectFinalState.hh \ + Projections/DISDiffHadron.hh \ Projections/DISFinalState.hh \ Projections/DISKinematics.hh \ Projections/DISLepton.hh \ + Projections/DISRapidityGap.hh \ Projections/DressedLeptons.hh \ Projections/EventMixingFinalState.hh \ Projections/FastJets.hh \ Projections/PxConePlugin.hh \ Projections/FinalPartons.hh \ Projections/FinalState.hh \ Projections/FoxWolframMoments.hh \ Projections/FParameter.hh \ Projections/GeneratedPercentileProjection.hh \ Projections/HadronicFinalState.hh \ Projections/HeavyHadrons.hh \ Projections/Hemispheres.hh \ Projections/IdentifiedFinalState.hh \ Projections/ImpactParameterProjection.hh \ Projections/IndirectFinalState.hh \ Projections/InitialQuarks.hh \ Projections/InvMassFinalState.hh \ Projections/JetAlg.hh \ Projections/JetShape.hh \ Projections/LeadingParticlesFinalState.hh \ Projections/LossyFinalState.hh \ Projections/MergedFinalState.hh \ Projections/MissingMomentum.hh \ Projections/NeutralFinalState.hh \ Projections/NonHadronicFinalState.hh \ Projections/NonPromptFinalState.hh \ Projections/ParisiTensor.hh \ Projections/ParticleFinder.hh \ Projections/PartonicTops.hh \ Projections/PercentileProjection.hh \ Projections/PrimaryHadrons.hh \ Projections/PrimaryParticles.hh \ Projections/PromptFinalState.hh \ Projections/SingleValueProjection.hh \ Projections/SmearedParticles.hh \ Projections/SmearedJets.hh \ Projections/SmearedMET.hh \ Projections/Sphericity.hh \ Projections/Spherocity.hh \ Projections/TauFinder.hh \ Projections/Thrust.hh \ Projections/TriggerCDFRun0Run1.hh \ Projections/TriggerCDFRun2.hh \ Projections/TriggerProjection.hh \ Projections/TriggerUA5.hh \ Projections/UnstableFinalState.hh \ + Projections/UndressBeamLeptons.hh \ Projections/UnstableParticles.hh \ Projections/UserCentEstimate.hh \ Projections/VetoedFinalState.hh \ Projections/VisibleFinalState.hh \ Projections/WFinder.hh \ Projections/ZFinder.hh ## Meta-projection convenience headers nobase_pkginclude_HEADERS += \ Projections/FinalStates.hh \ Projections/Smearing.hh ## Analysis base class headers # TODO: Move to Rivet/AnalysisTools header dir? nobase_pkginclude_HEADERS += \ Analyses/MC_Cent_pPb.hh \ Analyses/MC_ParticleAnalysis.hh \ Analyses/MC_JetAnalysis.hh \ Analyses/MC_JetSplittings.hh ## Tools nobase_pkginclude_HEADERS += \ Tools/AliceCommon.hh \ Tools/AtlasCommon.hh \ Tools/BeamConstraint.hh \ Tools/BinnedHistogram.hh \ Tools/CentralityBinner.hh \ Tools/Cmp.fhh \ Tools/Cmp.hh \ Tools/Correlators.hh \ Tools/Cutflow.hh \ Tools/Cuts.fhh \ Tools/Cuts.hh \ Tools/Exceptions.hh \ Tools/JetUtils.hh \ Tools/Logging.hh \ Tools/Random.hh \ Tools/ParticleBaseUtils.hh \ Tools/ParticleIdUtils.hh \ Tools/ParticleUtils.hh \ Tools/ParticleName.hh \ Tools/Percentile.hh \ Tools/PrettyPrint.hh \ Tools/RivetPaths.hh \ Tools/RivetSTL.hh \ Tools/RivetFastJet.hh \ Tools/RivetHepMC.hh \ Tools/RivetYODA.hh \ Tools/RivetMT2.hh \ Tools/SmearingFunctions.hh \ Tools/MomentumSmearingFunctions.hh \ Tools/ParticleSmearingFunctions.hh \ Tools/JetSmearingFunctions.hh \ Tools/TypeTraits.hh \ Tools/Utils.hh \ Tools/fjcontrib/AxesDefinition.hh \ Tools/fjcontrib/BottomUpSoftDrop.hh \ Tools/fjcontrib/EnergyCorrelator.hh \ Tools/fjcontrib/ExtraRecombiners.hh \ Tools/fjcontrib/IteratedSoftDrop.hh \ Tools/fjcontrib/MeasureDefinition.hh \ Tools/fjcontrib/ModifiedMassDropTagger.hh \ Tools/fjcontrib/Njettiness.hh \ Tools/fjcontrib/NjettinessPlugin.hh \ Tools/fjcontrib/Nsubjettiness.hh \ Tools/fjcontrib/Recluster.hh \ Tools/fjcontrib/RecursiveSoftDrop.hh \ Tools/fjcontrib/RecursiveSymmetryCutBase.hh \ Tools/fjcontrib/SoftDrop.hh \ Tools/fjcontrib/TauComponents.hh \ Tools/fjcontrib/XConePlugin.hh nobase_dist_noinst_HEADERS += \ Tools/osdir.hh ## Maths nobase_pkginclude_HEADERS += \ Math/Matrices.hh \ Math/Vector3.hh \ Math/VectorN.hh \ Math/MatrixN.hh \ Math/MatrixDiag.hh \ Math/MathHeader.hh \ Math/Vectors.hh \ Math/LorentzTrans.hh \ Math/Matrix3.hh \ Math/MathUtils.hh \ Math/Vector4.hh \ Math/Math.hh \ Math/Units.hh \ Math/Constants.hh \ Math/eigen/util.h \ Math/eigen/regressioninternal.h \ Math/eigen/regression.h \ Math/eigen/vector.h \ Math/eigen/ludecompositionbase.h \ Math/eigen/ludecomposition.h \ Math/eigen/linearsolver.h \ Math/eigen/linearsolverbase.h \ Math/eigen/matrix.h \ Math/eigen/vectorbase.h \ Math/eigen/projective.h \ Math/eigen/matrixbase.h diff --git a/include/Rivet/Projection.hh b/include/Rivet/Projection.hh --- a/include/Rivet/Projection.hh +++ b/include/Rivet/Projection.hh @@ -1,185 +1,202 @@ // -*- C++ -*- #ifndef RIVET_Projection_HH #define RIVET_Projection_HH #include "Rivet/Projection.fhh" #include "Rivet/ProjectionApplier.hh" #include "Rivet/ProjectionHandler.hh" #include "Rivet/Config/RivetCommon.hh" #include "Rivet/Tools/Cuts.hh" // NOTE: Cmp.hh, Event.hh and Particle.hh included at the bottom namespace Rivet { // Forward declaration class Event; /// @brief Base class for all Rivet projections. /// /// Projection is the base class of all Projections to be used by /// Rivet. A Projection object can be assigned to an Event object and /// will then define a processed part of the information available in /// the Event, which then can be used by other Projection objects /// and/or Analysis objects. /// /// The main virtual functions to be overridden by concrete sub-classes /// are project(const Event &) and compare(const Projection &). class Projection : public ProjectionApplier { public: /// Event is a friend. friend class Event; /// The Cmp specialization for Projection is a friend. friend class Cmp; /// @name Standard constructors and destructors. //@{ /// The default constructor. Projection(); /// Clone on the heap. virtual unique_ptr clone() const = 0; /// The destructor. virtual ~Projection(); //@} /// Get the name of the projection. virtual std::string name() const { return _name; } + /// Get the state of the projetion. + bool valid() const { + return _isValid; + } + + /// Get the state of the projetion. + bool failed() const { + return !valid(); + } /// @name Projection operation and comparison //@{ /// Take the information available in the Event and make the /// calculations necessary to obtain the projection. Note that this /// function must never be called except inside the /// Event::applyProjection(Projection *) function. virtual void project(const Event& e) = 0; /// This function is used to define a unique ordering between /// different Projection objects of the same class. If this is /// considered to be equivalent to the Projector object, \a p, in the /// argument the function should return 0. If this object should be /// ordered before \a p a negative value should be returned, /// otherwise a positive value should be returned. This function must /// never be called explicitly, but should only be called from the /// operator<(const Projection &). When implementing the function in /// concrete sub-classes, it is then guaranteed that the Projection /// object \a p in the argument is of the same class as the sub-class /// and can be safely dynamically casted to that class. /// /// When implementing this function in a sub-class, the immediate /// base class version of the function should be called first. If the /// base class function returns a non-zero value, that value should /// be returned immediately. Only if zero is returned should this /// function check the member variables of the sub-class to determine /// whether this should be ordered before or after \a p, or if it is /// equivalent with \a p. virtual int compare(const Projection& p) const = 0; /// Determine whether this object should be ordered before the object /// \a p given as argument. If \a p is of a different class than /// this, the before() function of the corresponding type_info /// objects is used. Otherwise, if the objects are of the same class, /// the virtual compare(const Projection &) will be returned. bool before(const Projection& p) const; //@} /// @name Beam configuration /// @todo Does it really make sense to restrict Projections to particular beam configs? Do we use this in practice? //@{ /// Return the allowed beam pairs on which this projection can operate, not /// including recursion. Derived classes should ensure that all contained /// projections are registered in the @a _projections set for the beam /// constraint chaining to work. /// @todo Remove the beam constraints system from projections. virtual const std::set beamPairs() const; /// Add a colliding beam pair. /// @todo This deserves a better name! Projection& addPdgIdPair(PdgId beam1, PdgId beam2) { _beamPairs.insert(PdgIdPair(beam1, beam2)); return *this; } //@} protected: /// Get a Log object based on the getName() property of the calling projection object. Log& getLog() const { string logname = "Rivet.Projection." + name(); return Log::getLog(logname); } /// Used by derived classes to set their name. void setName(const std::string& name) { _name = name; } + /// Set the projection in an unvalid state. + void fail() { + _isValid = false; + } + /// Shortcut to make a named Cmp comparison with the @c *this /// object automatically passed as one of the parent projections. Cmp mkNamedPCmp(const Projection& otherparent, const std::string& pname) const; /// Shortcut to make a named Cmp comparison with the @c *this /// object automatically passed as one of the parent projections. /// /// @note Alias for mkNamedPCmp Cmp mkPCmp(const Projection& otherparent, const std::string& pname) const; /// Block Projection copying virtual Projection& operator = (const Projection&); private: /// Name variable is used by the base class messages to identify /// which derived class is being handled. string _name; /// Beam-type constraint. /// @todo Remove? set _beamPairs; + /// Flag to tell if this projection is in a valid state. + bool _isValid; + }; } /// Define "less" operator for Projection* containers in terms of the Projection::before virtual method. inline bool std::less::operator()(const Rivet::Projection* x, const Rivet::Projection* y) const { return x->before(*y); } #endif #include "Rivet/Event.hh" #include "Rivet/Particle.hh" #include "Rivet/Tools/Cmp.hh" /// @def DEFAULT_RIVET_PROJ_CLONE /// Preprocessor define to prettify the manky constructor with name string argument #define DEFAULT_RIVET_PROJ_CLONE(clsname) \ virtual unique_ptr clone() const { return unique_ptr(new clsname(*this)); } diff --git a/include/Rivet/Projections/Beam.hh b/include/Rivet/Projections/Beam.hh --- a/include/Rivet/Projections/Beam.hh +++ b/include/Rivet/Projections/Beam.hh @@ -1,206 +1,208 @@ // -*- C++ -*- #ifndef RIVET_Beam_HH #define RIVET_Beam_HH #include "Rivet/Projection.hh" #include "Rivet/Event.hh" #include "Rivet/Particle.hh" #include "Rivet/Math/LorentzTrans.hh" namespace Rivet { /// @name Standalone beam kinematics functions //@{ /// Get beam particles from an event ParticlePair beams(const Event& e); /// Get beam particle IDs from a pair of Particles /// @deprecated Use pids(beams) inline PdgIdPair beamIds(const ParticlePair& beams) { return pids(beams); } /// Get beam particle IDs from an event /// @deprecated Use pids(e.beams()) inline PdgIdPair beamIds(const Event& e) { return pids(beams(e)); } /// Get beam centre-of-mass energy from a pair of beam momenta double sqrtS(const FourMomentum& pa, const FourMomentum& pb); /// Get beam centre-of-mass energy from a pair of Particles inline double sqrtS(const ParticlePair& beams) { return sqrtS(beams.first.momentum(), beams.second.momentum()); } /// Get beam centre-of-mass energy from an Event inline double sqrtS(const Event& e) { return sqrtS(beams(e)); } /// Get per-nucleon beam centre-of-mass energy from a pair of beam momenta /// @note Uses a nominal nucleon mass of 0.939 GeV to convert masses to A double asqrtS(const FourMomentum& pa, const FourMomentum& pb); /// Get per-nucleon beam centre-of-mass energy from a pair of Particles /// @note Uses the sum of nuclear mass numbers A for each beam double asqrtS(const ParticlePair& beams); /// Get per-nucleon beam centre-of-mass energy from an Event /// @note Uses the sum of nuclear mass numbers A for each beam inline double asqrtS(const Event& e) { return asqrtS(beams(e)); } /// Get the Lorentz boost to the beam centre-of-mass system (CMS) from a pair of beam momenta inline FourMomentum cmsBoostVec(const FourMomentum& pa, const FourMomentum& pb) { return pa + pb; } /// Get the Lorentz boost to the beam centre-of-mass system (CMS) from a pair of Particles inline FourMomentum cmsBoostVec(const ParticlePair& beams) { return cmsBoostVec(beams.first, beams.second); } /// Get the Lorentz boost to the beam centre-of-mass system (CMS) from a pair of beam momenta FourMomentum acmsBoostVec(const FourMomentum& pa, const FourMomentum& pb); /// Get the Lorentz boost to the beam centre-of-mass system (CMS) from a pair of Particles FourMomentum acmsBoostVec(const ParticlePair& beams); /// Get the Lorentz boost to the beam centre-of-mass system (CMS) from a pair of beam momenta Vector3 cmsBetaVec(const FourMomentum& pa, const FourMomentum& pb); /// Get the Lorentz boost to the beam centre-of-mass system (CMS) from a pair of Particles inline Vector3 cmsBetaVec(const ParticlePair& beams) { return cmsBetaVec(beams.first, beams.second); } /// Get the Lorentz boost to the per-nucleon beam centre-of-mass system (ACMS) from a pair of beam momenta /// @note Uses a nominal nucleon mass of 0.939 GeV to convert masses to A Vector3 acmsBetaVec(const FourMomentum& pa, const FourMomentum& pb); /// Get the Lorentz boost to the per-nucleon beam centre-of-mass system (ACMS) from a pair of Particles /// @note Uses the sum of nuclear mass numbers A for each beam Vector3 acmsBetaVec(const ParticlePair& beams); /// Get the Lorentz boost to the beam centre-of-mass system (CMS) from a pair of beam momenta Vector3 cmsGammaVec(const FourMomentum& pa, const FourMomentum& pb); /// Get the Lorentz boost to the beam centre-of-mass system (CMS) from a pair of Particles inline Vector3 cmsGammaVec(const ParticlePair& beams) { return cmsGammaVec(beams.first, beams.second); } /// Get the Lorentz boost to the per-nucleon beam centre-of-mass system (ACMS) from a pair of beam momenta /// @note Uses a nominal nucleon mass of 0.939 GeV to convert masses to A Vector3 acmsGammaVec(const FourMomentum& pa, const FourMomentum& pb); /// Get the Lorentz boost to the per-nucleon beam centre-of-mass system (ACMS) from a pair of Particles /// @note Uses the sum of nuclear mass numbers A for each beam Vector3 acmsGammaVec(const ParticlePair& beams); /// Get the Lorentz transformation to the beam centre-of-mass system (CMS) from a pair of beam momenta LorentzTransform cmsTransform(const FourMomentum& pa, const FourMomentum& pb); /// Get the Lorentz transformation to the beam centre-of-mass system (CMS) from a pair of Particles inline LorentzTransform cmsTransform(const ParticlePair& beams) { return cmsTransform(beams.first, beams.second); } /// Get the Lorentz transformation to the per-nucleon beam centre-of-mass system (CMS) from a pair of beam momenta /// @note Uses a nominal nucleon mass of 0.939 GeV to convert masses to A LorentzTransform acmsTransform(const FourMomentum& pa, const FourMomentum& pb); /// Get the Lorentz transformation to the per-nucleon beam centre-of-mass system (CMS) from a pair of Particles /// @note Uses the sum of nuclear mass numbers A for each beam LorentzTransform acmsTransform(const ParticlePair& beams); //@} /// @brief Project out the incoming beams class Beam : public Projection { public: /// Default (and only) constructor Beam() { setName("Beam"); } /// Clone on the heap DEFAULT_RIVET_PROJ_CLONE(Beam); /// @name Beam particles and kinematics //@{ /// The pair of beam particles in the current collision const ParticlePair& beams() const { return _theBeams; } /// The pair of beam particle PDG codes in the current collision /// @deprecated Use pids(beams()) PdgIdPair beamIds() const { return pids(beams()); } /// Get centre of mass energy, \f$ \sqrt{s} \f$ double sqrtS() const { return Rivet::sqrtS(beams()); } /// Get the Lorentz boost to the beam centre-of-mass FourMomentum cmsBoostVec() const { return Rivet::cmsBoostVec(beams()); } /// Get the Lorentz transform to the beam centre-of-mass LorentzTransform cmsTransform() const { return Rivet::cmsTransform(beams()); } /// Get the beta factor vector for the Lorentz boost to the beam centre-of-mass Vector3 cmsBetaVec() const { return Rivet::cmsBetaVec(beams()); } /// Get the gamma factor vector for the Lorentz boost to the beam centre-of-mass Vector3 cmsGammaVec() const { return Rivet::cmsGammaVec(beams()); } //@} /// @name Per-nucleon beam kinematics //@{ /// Get per-nucleon centre of mass energy, \f$ \sqrt{s}/(A_1 + A_2) \f$ double asqrtS() const { return Rivet::asqrtS(beams()); } /// Get the Lorentz boost to the per-nucleon beam centre-of-mass Vector3 acmsBetaVec() const { return Rivet::acmsBetaVec(beams()); } /// Get the Lorentz boost to the per-nucleon beam centre-of-mass Vector3 acmsGammaVec() const { return Rivet::acmsGammaVec(beams()); } /// Get the Lorentz transform to the per-nucleon beam centre-of-mass LorentzTransform acmsTransform() const { return Rivet::acmsTransform(beams()); } //@} /// Get the beam interaction primary vertex (PV) position FourVector pv() const; /// Project on to the Event virtual void project(const Event& e); + protected: + + /// The beam particles in the current collision + ParticlePair _theBeams; + private: /// Compare with other projections -- it's always the same, since there are no params virtual int compare(const Projection&) const { return EQUIVALENT; } - /// The beam particles in the current collision - ParticlePair _theBeams; - }; } #endif diff --git a/include/Rivet/Projections/DISDiffHadron.hh b/include/Rivet/Projections/DISDiffHadron.hh new file mode 100644 --- /dev/null +++ b/include/Rivet/Projections/DISDiffHadron.hh @@ -0,0 +1,64 @@ +// -*- C++ -*- +#ifndef RIVET_DISDiffHadron_HH +#define RIVET_DISDiffHadron_HH + +#include "Rivet/Projections/Beam.hh" +#include "Rivet/Projections/HadronicFinalState.hh" +#include "Rivet/Particle.hh" +#include "Rivet/Event.hh" + +namespace Rivet { + + + /// @brief Get the incoming and outgoing hadron in a diffractive ep + /// event. + class DISDiffHadron : public Projection { + public: + + /// @name Constructors. + //@{ + + /// Default constructor. + DISDiffHadron() { + setName("DISDiffHadron"); + addProjection(Beam(), "Beam"); + addProjection(FinalState(), "FS"); + } + + /// Clone on the heap. + DEFAULT_RIVET_PROJ_CLONE(DISDiffHadron); + + //@} + + + protected: + + /// Perform the projection operation on the supplied event. + virtual void project(const Event& e); + + /// Compare with other projections. + virtual int compare(const Projection& p) const; + + + public: + + /// The incoming lepton + const Particle& in() const { return _incoming; } + + /// The outgoing lepton + const Particle& out() const { return _outgoing; } + + private: + + /// The incoming lepton + Particle _incoming; + + /// The outgoing lepton + Particle _outgoing; + + }; + +} + + +#endif diff --git a/include/Rivet/Projections/DISLepton.hh b/include/Rivet/Projections/DISLepton.hh --- a/include/Rivet/Projections/DISLepton.hh +++ b/include/Rivet/Projections/DISLepton.hh @@ -1,111 +1,138 @@ // -*- C++ -*- #ifndef RIVET_DISLepton_HH #define RIVET_DISLepton_HH #include "Rivet/Projections/Beam.hh" #include "Rivet/Projections/PromptFinalState.hh" #include "Rivet/Projections/HadronicFinalState.hh" #include "Rivet/Projections/DressedLeptons.hh" +#include "Rivet/Projections/UndressBeamLeptons.hh" #include "Rivet/Particle.hh" #include "Rivet/Event.hh" namespace Rivet { /// @brief Get the incoming and outgoing leptons in a DIS event. class DISLepton : public Projection { public: + /// Enum to enable different orderings for selecting scattered + /// leptons in case several were found. + enum SortOrder { ENERGY, ETA, ET }; + /// @name Constructors. //@{ /// Default constructor taking general options. The recognised /// options are: LMODE, taking the options "prompt", "any" and /// "dressed"; DressedDR giving a delta-R cone radius where photon /// momenta are added to the lepton candidates for LMODE=dresses; - /// and IsolDR giving a cone in delta-R where no hadrons are - /// allowed around a lepton candidate. + /// IsolDR giving a cone in delta-R where no hadrons are allowed + /// around a lepton candidate; and Undress giving a cone around + /// the incoming incoming beam in which photons are considered + /// initial state rafiation for which the momentum is subtracted + /// from the beam momentum. DISLepton(const std::map & opts = - std::map()): _isolDR(0.0) { + std::map()) + : _isolDR(0.0), _sort(ENERGY) { setName("DISLepton"); - addProjection(Beam(), "Beam"); addProjection(HadronicFinalState(), "IFS"); + auto sorting = opts.find("LSort"); + if ( sorting != opts.end() && sorting->second == "ETA" ) + _sort = ETA; + else if ( sorting != opts.end() && sorting->second == "ET" ) + _sort = ET; + + double undresstheta = 0.0; + auto undress = opts.find("Undress"); + if ( undress != opts.end() ) + undresstheta = std::stod(undress->second); + if ( undresstheta > 0.0 ) + addProjection(UndressBeamLeptons(undresstheta), "Beam"); + else + addProjection(Beam(), "Beam"); + auto isol = opts.find("IsolDR"); if ( isol != opts.end() ) _isolDR = std::stod(isol->second); double dressdr = 0.0; auto dress = opts.find("DressDR"); if ( dress != opts.end() ) dressdr = std::stod(dress->second); - auto lmode = opts.find("LMODE"); + auto lmode = opts.find("LMode"); if ( lmode != opts.end() && lmode->second == "any" ) addProjection(FinalState(), "LFS"); else if ( lmode != opts.end() && lmode->second == "dressed" ) addProjection(DressedLeptons(dressdr), "LFS"); else addProjection(PromptFinalState(), "LFS"); } /// Constructor taking the following arguments: a final state /// projection defining which lepton candidates to consider; a /// beam projection detining the momenta of the incoming lepton /// beam, and a final state projection defining the particles not /// allowed witin a delta-R of @a isolationcut of a lepton /// candidate. DISLepton(const FinalState & leptoncandidates, const Beam & beamproj = Beam(), const FinalState & isolationfs = FinalState(), - double isolationcut = 0.0): _isolDR(isolationcut) { + double isolationcut = 0.0, SortOrder sorting = ENERGY) + : _isolDR(isolationcut), _sort(sorting) { addProjection(leptoncandidates, "LFS"); addProjection(isolationfs, "IFS"); addProjection(beamproj, "Beam"); } /// Clone on the heap. DEFAULT_RIVET_PROJ_CLONE(DISLepton); //@} protected: /// Perform the projection operation on the supplied event. virtual void project(const Event& e); /// Compare with other projections. virtual int compare(const Projection& p) const; public: /// The incoming lepton const Particle& in() const { return _incoming; } /// The outgoing lepton const Particle& out() const { return _outgoing; } /// Sign of the incoming lepton pz component int pzSign() const { return sign(_incoming.pz()); } private: /// The incoming lepton Particle _incoming; /// The outgoing lepton Particle _outgoing; /// If larger than zerp an isolation cut around the lepton is required. double _isolDR; + /// How to sort leptons + SortOrder _sort; + }; } #endif diff --git a/include/Rivet/Projections/DISRapidityGap.hh b/include/Rivet/Projections/DISRapidityGap.hh new file mode 100644 --- /dev/null +++ b/include/Rivet/Projections/DISRapidityGap.hh @@ -0,0 +1,94 @@ +// -*- C++ -*- +#ifndef RIVET_DISRapidityGap_HH +#define RIVET_DISRapidityGap_HH + +#include "Rivet/Projections/DISKinematics.hh" +#include "Rivet/Projections/DISFinalState.hh" +#include "Rivet/Particle.hh" +#include "Rivet/Event.hh" + +namespace Rivet { + + + /// @brief Get the incoming and outgoing hadron in a diffractive ep + /// event. + class DISRapidityGap : public Projection { + + public: + + /// Type of DIS boost to apply + enum Frame { HCM, LAB, XCM }; + + DISRapidityGap() { + setName("DISRapidityGap"); + addProjection(DISKinematics(), "DISKIN"); + addProjection(DISFinalState(DISFinalState::HCM), "DISFS"); + } + + DEFAULT_RIVET_PROJ_CLONE(DISRapidityGap); + + 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; + else return _ePpzX_HCM; + } + const double EmPzX(Frame f) const { + if (f == LAB) return _eMpzX_LAB; + else if (f == XCM) return _eMpzX_XCM; + else return _eMpzX_HCM; + } + const FourMomentum pX(Frame f) const { + if (f == LAB) return _momX_LAB; + else if (f == XCM) return _momX_XCM; + else return _momX_HCM; + } + const FourMomentum pY(Frame f) const { + if (f == LAB) return _momY_LAB; + else if (f == XCM) return _momY_XCM; + else return _momY_HCM; + } + const Particles& systemX(Frame f) const { + if (f == LAB) return _pX_LAB; + else if (f == XCM) return _pX_XCM; + else return _pX_HCM; + } + const Particles& systemY(Frame f) const { + if (f == LAB) return _pY_LAB; + else if (f == XCM) return _pY_XCM; + else return _pY_HCM; + } + + protected: + + virtual int compare(const Projection& p) const; + + virtual void project(const Event& e); + + void clearAll(); + + void findgap(const Particles& particles, const DISKinematics& diskin); + + private: + + double _M2X, _M2Y, _t; + double _gap, _gapUpp, _gapLow; + double _ePpzX_LAB, _eMpzX_LAB; + double _ePpzX_HCM, _eMpzX_HCM; + double _ePpzX_XCM, _eMpzX_XCM; + FourMomentum _momX_HCM, _momY_HCM; + FourMomentum _momX_LAB, _momY_LAB; + FourMomentum _momX_XCM, _momY_XCM; + Particles _pX_HCM, _pY_HCM, _pX_LAB, _pY_LAB, _pX_XCM, _pY_XCM; + + }; + +} + + +#endif diff --git a/include/Rivet/Projections/UndressBeamLeptons.hh b/include/Rivet/Projections/UndressBeamLeptons.hh new file mode 100644 --- /dev/null +++ b/include/Rivet/Projections/UndressBeamLeptons.hh @@ -0,0 +1,46 @@ +// -*- C++ -*- +#ifndef RIVET_UndressBeamLeptons_HH +#define RIVET_UndressBeamLeptons_HH + +#include "Rivet/Projections/Beam.hh" +#include "Rivet/Projections/FinalState.hh" + +namespace Rivet { + + + /// @brief Project out the incoming beams, but subtract any colinear + /// photons from lepton beams within a given cone. + class UndressBeamLeptons : public Beam { + public: + + /// Default (and only) constructor. Takes an angle as + /// argument. The momentum of any photon within This angle wrt. a + /// charged lepton beam will be subtracted from the beam lepton + /// momentum. + UndressBeamLeptons(double theta = 0.0): _thetamax(theta) { + setName("UndressBeamLeptons"); + addProjection(FinalState(), "FS"); + } + + /// Clone on the heap + DEFAULT_RIVET_PROJ_CLONE(UndressBeamLeptons); + + + /// Project on to the Event + virtual void project(const Event& e); + + + private: + + /// Compare with other projections. + virtual int compare(const Projection & p) const; + + /// The beam particles in the current collision + double _thetamax; + + }; + + +} + +#endif diff --git a/src/Core/Projection.cc b/src/Core/Projection.cc --- a/src/Core/Projection.cc +++ b/src/Core/Projection.cc @@ -1,60 +1,60 @@ // -*- C++ -*- #include "Rivet/Event.hh" #include "Rivet/Projection.hh" #include "Rivet/Tools/Logging.hh" #include "Rivet/Tools/BeamConstraint.hh" #include "Rivet/Tools/Cmp.hh" namespace Rivet { Projection::Projection() - : _name("BaseProjection") + : _name("BaseProjection"), _isValid(true) { addPdgIdPair(PID::ANY, PID::ANY); } Projection:: ~Projection() { } Projection& Projection::operator = (const Projection&) { return *this; } bool Projection::before(const Projection& p) const { const std::type_info& thisid = typeid(*this); const std::type_info& otherid = typeid(p); if (thisid == otherid) { const bool cmp = compare(p) < 0; MSG_TRACE("Comparing projections of same RTTI type: " << this << " < " << &p << " = " << cmp); return cmp; } else { const bool cmp = thisid.before(otherid); MSG_TRACE("Ordering projections of different RTTI type: " << this << " < " << &p << " = " << cmp); return cmp; } } const set Projection::beamPairs() const { set ret = _beamPairs; set projs = getProjections(); for (set::const_iterator ip = projs.begin(); ip != projs.end(); ++ip) { ConstProjectionPtr p = *ip; getLog() << Log::TRACE << "Proj addr = " << p << endl; if (p) ret = intersection(ret, p->beamPairs()); } return ret; } Cmp Projection::mkNamedPCmp(const Projection& otherparent, const string& pname) const { return pcmp(*this, otherparent, pname); } Cmp Projection::mkPCmp(const Projection& otherparent, const string& pname) const { return pcmp(*this, otherparent, pname); } } diff --git a/src/Projections/CentralEtHCM.cc b/src/Projections/CentralEtHCM.cc --- a/src/Projections/CentralEtHCM.cc +++ b/src/Projections/CentralEtHCM.cc @@ -1,17 +1,21 @@ // -*- C++ -*- #include "Rivet/Projections/CentralEtHCM.hh" namespace Rivet { void CentralEtHCM::project(const Event& e) { const DISFinalState& fs = applyProjection(e, "FS"); + if ( fs.failed() ) { + fail(); + return; + } _sumet = 0.0; for (const Particle& p : fs.particles()) { /// @todo Generalise rapidity cut value if (fabs(p.rapidity()) < 0.5) _sumet += p.Et(); } } } diff --git a/src/Projections/DISDiffHadron.cc b/src/Projections/DISDiffHadron.cc new file mode 100644 --- /dev/null +++ b/src/Projections/DISDiffHadron.cc @@ -0,0 +1,49 @@ +// -*- C++ -*- +#include "Rivet/Projections/DISDiffHadron.hh" + +namespace Rivet { + + + int DISDiffHadron::compare(const Projection& p) const { + return mkNamedPCmp(p, "Beam") || mkNamedPCmp(p, "FS"); + } + + + void DISDiffHadron::project(const Event& e) { + + // Find incoming lepton beam + const ParticlePair& inc = applyProjection(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 { + fail(); + return; + } + + const FinalState & fs = applyProjection(e, "FS"); + Particles fshadrons; + if ( _incoming.momentum().pz() >= 0.0 ) + fshadrons = fs.particles(isHadron, cmpMomByDescEta); + else + fshadrons = fs.particles(isHadron, cmpMomByEta); + + Particles sfhadrons = filter_select(fshadrons, + Cuts::pid == _incoming.pid()); + MSG_DEBUG("SF hadrons = " << sfhadrons.size() << + ", all hadrons = " << fshadrons.size()); + if (!sfhadrons.empty()) { + _outgoing = sfhadrons.front(); + } else if (!fshadrons.empty()) { + _outgoing = fshadrons.front(); + } else { + fail(); + } + + } + + +} diff --git a/src/Projections/DISFinalState.cc b/src/Projections/DISFinalState.cc --- a/src/Projections/DISFinalState.cc +++ b/src/Projections/DISFinalState.cc @@ -1,33 +1,41 @@ // -*- C++ -*- #include "Rivet/Projections/DISFinalState.hh" namespace Rivet { void DISFinalState::project(const Event& e) { const DISKinematics& diskin = apply(e, "Kinematics"); + if ( diskin.failed() ) { + fail(); + return; + } LorentzTransform hcmboost; //< Null boost = LAB frame by default if (_boosttype == HCM) hcmboost = diskin.boostHCM(); else if (_boosttype == BREIT) hcmboost = diskin.boostBreit(); const DISLepton& dislep = diskin.apply(e, "Lepton"); + if ( diskin.failed() ) { + fail(); + return; + } const FinalState& fs = apply(e, "FS"); // Fill the particle list with all particles _other_ than the DIS scattered // lepton, with momenta boosted into the appropriate frame. _theParticles.clear(); _theParticles.reserve(fs.particles().size()-1); const GenParticle* dislepGP = dislep.out().genParticle(); // const GenParticle* dislepIN = dislep.in().genParticle(); for (const Particle& p : fs.particles()) { ///< Ensure that we skip the DIS lepton Particle temp = p; if (_boosttype != LAB) temp.setMomentum(hcmboost.transform(temp.momentum())); if (p.genParticle() != dislepGP) _theParticles.push_back(temp); } } } diff --git a/src/Projections/DISKinematics.cc b/src/Projections/DISKinematics.cc --- a/src/Projections/DISKinematics.cc +++ b/src/Projections/DISKinematics.cc @@ -1,80 +1,85 @@ // -*- C++ -*- #include "Rivet/Projections/DISKinematics.hh" #include "Rivet/Math/Constants.hh" namespace Rivet { void DISKinematics::project(const Event& e) { // Find appropriate DIS leptons const DISLepton& dislep = applyProjection(e, "Lepton"); + if ( dislep.failed() ) { + fail(); + return; + } _outLepton = dislep.out(); // Identify beam hadron const ParticlePair& inc = applyProjection(e, "Beam").beams(); const bool firstIsHadron = PID::isHadron(inc.first.pid()); const bool secondIsHadron = PID::isHadron(inc.second.pid()); if (firstIsHadron && !secondIsHadron) { _inHadron = inc.first; _inLepton = dislep.in(); // inc.second; } else if (!firstIsHadron && secondIsHadron) { _inHadron = inc.second; _inLepton = dislep.in(); // inc.first; } else { - throw Error("DISKinematics could not find the correct beam hadron"); + fail(); + return; } // Get the DIS lepton and store some of its properties const FourMomentum pHad = _inHadron.momentum(); const FourMomentum pLepIn = _inLepton.momentum(); const FourMomentum pLepOut = _outLepton.momentum(); const FourMomentum pGamma = pLepIn - pLepOut; const FourMomentum tothad = pGamma + pHad; _theQ2 = -pGamma.mass2(); _theW2 = tothad.mass2(); _theX = Q2()/(2.0 * pGamma * pHad); _theY = (pGamma * pHad) / (pLepIn * pHad); _theS = invariant(pLepIn + pHad); // Calculate boost vector for boost into HCM-system LorentzTransform tmp; tmp.setBetaVec(-tothad.boostVector()); // Rotate so the photon is in x-z plane in HCM rest frame FourMomentum pGammaHCM = tmp.transform(pGamma); tmp.preMult(Matrix3(Vector3::mkZ(), -pGammaHCM.azimuthalAngle())); pGammaHCM = tmp.transform(pGamma); assert(isZero(dot(pGammaHCM.vector3(), Vector3::mkY()))); // Rotate so the photon is along the positive z-axis const double rot_angle = pGammaHCM.polarAngle() * (pGammaHCM.px() >= 0 ? -1 : 1); tmp.preMult(Matrix3(Vector3::mkY(), rot_angle)); // Check that final HCM photon lies along +ve z as expected pGammaHCM = tmp.transform(pGamma); assert(isZero(dot(pGammaHCM.vector3(), Vector3::mkX()), 1e-3)); assert(isZero(dot(pGammaHCM.vector3(), Vector3::mkY()), 1e-3)); assert(isZero(angle(pGammaHCM.vector3(), Vector3::mkZ()), 1e-3)); // Finally rotate so outgoing lepton at phi = 0 FourMomentum pLepOutHCM = tmp.transform(pLepOut); tmp.preMult(Matrix3(Vector3::mkZ(), -pLepOutHCM.azimuthalAngle())); assert(isZero(tmp.transform(pLepOut).azimuthalAngle())); _hcm = tmp; // Boost to Breit frame (use opposite convention for photon --- along *minus* z) tmp.preMult(Matrix3(Vector3::mkX(), PI)); const double bz = 1 - 2*x(); _breit = LorentzTransform::mkObjTransformFromBeta(Vector3::mkZ() * bz).combine(tmp); assert(isZero(angle(_breit.transform(pGamma).vector3(), -Vector3::mkZ()), 1e-3)); assert(isZero(_breit.transform(pLepOut).azimuthalAngle(), 1e-3)); } int DISKinematics::compare(const Projection & p) const { const DISKinematics& other = pcast(p); return mkNamedPCmp(other, "Lepton"); } } diff --git a/src/Projections/DISLepton.cc b/src/Projections/DISLepton.cc --- a/src/Projections/DISLepton.cc +++ b/src/Projections/DISLepton.cc @@ -1,69 +1,78 @@ // -*- C++ -*- #include "Rivet/Projections/DISLepton.hh" namespace Rivet { int DISLepton::compare(const Projection& p) const { const DISLepton& other = pcast(p); return mkNamedPCmp(other, "Beam") || mkNamedPCmp(other, "LFS") || - mkNamedPCmp(other, "IFS"); + mkNamedPCmp(other, "IFS") || cmp(_sort, other._sort); } void DISLepton::project(const Event& e) { // Find incoming lepton beam const ParticlePair& inc = applyProjection(e, "Beam").beams(); bool firstIsLepton = PID::isLepton(inc.first.pid()); bool secondIsLepton = PID::isLepton(inc.second.pid()); if (firstIsLepton && !secondIsLepton) { _incoming = inc.first; } else if (!firstIsLepton && secondIsLepton) { _incoming = inc.second; } else { - throw Error("DISLepton could not find the correct beam"); + fail(); + return; } // If no graph-connected scattered lepton, use the hardest // (preferably same-flavour) prompt FS lepton in the event. - - const Particles fsleptons = - applyProjection(e, "LFS").particles(isLepton, cmpMomByE); + const FinalState & fs = applyProjection(e, "LFS"); + Particles fsleptons; + if ( _sort == ET ) + fsleptons = fs.particles(isLepton, cmpMomByEt); + else if ( _sort == ETA && _incoming.momentum().pz() >= 0.0 ) + fsleptons = fs.particles(isLepton, cmpMomByDescEta); + else if ( _sort == ETA && _incoming.momentum().pz() < 0.0 ) + fsleptons = fs.particles(isLepton, cmpMomByEta); + else + fsleptons = fs.particles(isLepton, cmpMomByE); + Particles sfleptons = filter_select(fsleptons, Cuts::pid == _incoming.pid()); MSG_DEBUG("SF leptons = " << sfleptons.size() << ", all leptons = " << fsleptons.size()); if ( sfleptons.empty() ) sfleptons = fsleptons; if ( _isolDR > 0.0 ) { const Particles & other = applyProjection(e, "IFS").particles(); while (!sfleptons.empty()) { bool skip = false; Particle testlepton = sfleptons.front(); for ( auto p: other ) { if ( skip ) break; if ( deltaR(p, testlepton) < _isolDR ) skip = true; for ( auto c : testlepton.constituents() ) { if ( c.genParticle() == p.genParticle() ) { skip = false; break; } } } if ( !skip ) break; sfleptons.erase(sfleptons.begin()); } } if ( !sfleptons.empty() ) { _outgoing = sfleptons.front(); } else { - throw Error("Could not find the scattered lepton"); + fail(); } } } diff --git a/src/Projections/DISRapidityGap.cc b/src/Projections/DISRapidityGap.cc new file mode 100644 --- /dev/null +++ b/src/Projections/DISRapidityGap.cc @@ -0,0 +1,154 @@ +// -*- C++ -*- +#include "Rivet/Projections/DISRapidityGap.hh" + +namespace Rivet { + + + int DISRapidityGap::compare(const Projection& p) const { + return mkNamedPCmp(p, "DISKIN") || mkNamedPCmp(p, "DISFS"); + } + + + void DISRapidityGap::project(const Event& e) { + const DISKinematics& dk = + apply(e, "DISKIN"); + const Particles& p = + apply(e, "DISFS").particles(cmpMomByEta); + findgap(p, dk); + } + + void DISRapidityGap::clearAll() { + _M2X = _M2Y = _t = _gap = 0.; + _gapUpp = _gapLow = -8.; + _ePpzX_HCM = _eMpzX_HCM =_ePpzX_LAB = + _eMpzX_LAB = _ePpzX_XCM = _eMpzX_XCM = 0.; + _momX_HCM.setPE(0., 0., 0., 0.); + _momY_HCM.setPE(0., 0., 0., 0.); + _momX_XCM.setPE(0., 0., 0., 0.); + _momY_XCM.setPE(0., 0., 0., 0.); + _momX_LAB.setPE(0., 0., 0., 0.); + _momY_LAB.setPE(0., 0., 0., 0.); + _pX_HCM.clear(); + _pY_HCM.clear(); + _pX_XCM.clear(); + _pY_XCM.clear(); + _pX_LAB.clear(); + _pY_LAB.clear(); + } + + void DISRapidityGap::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 = abs(particles[i+1].eta() - particles[i].eta()); + if (tmpGap > _gap) { + _gap = tmpGap; + _gapLow = (dir > 0) ? particles[i].eta() : dir * particles[i+1].eta(); + _gapUpp = (dir > 0) ? particles[i+1].eta() : dir * particles[i].eta(); + } + } + + // Define the two systems X and Y. + Particles tmp_pX, tmp_pY; + foreach (const Particle& ip, particles) { + if (dir * ip.eta() > _gapLow) tmp_pX.push_back(ip); + else tmp_pY.push_back(ip); + } + + Particles pX, pY; + pX = (dir < 0) ? tmp_pY : tmp_pX; + pY = (dir < 0) ? tmp_pX : tmp_pY; + + // 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; + for (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; + for (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); + } + + for (const Particle& jp : pX) { + // Boost from HCM to LAB. + FourMomentum lab = hcminverse.transform(jp.momentum()); + _ePpzX_LAB += lab.E() + dir * lab.pz(); + _eMpzX_LAB += lab.E() - dir * 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); + } + } + + for ( 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); + } + } + + // 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. + FourMomentum proton = hcmboost.transform(diskin.beamHadron().momentum()); + FourMomentum pPom = proton - _momY_HCM; + _t = pPom * pPom; + + } +} diff --git a/src/Projections/Makefile.am b/src/Projections/Makefile.am --- a/src/Projections/Makefile.am +++ b/src/Projections/Makefile.am @@ -1,47 +1,50 @@ noinst_LTLIBRARIES = libRivetProjections.la libRivetProjections_la_SOURCES = \ Beam.cc \ BeamThrust.cc \ ChargedFinalState.cc \ ChargedLeptons.cc \ CentralEtHCM.cc \ + DISDiffHadron.cc \ DISFinalState.cc \ DISKinematics.cc \ DISLepton.cc \ + DISRapidityGap.cc \ DressedLeptons.cc \ FastJets.cc \ PxConePlugin.cc \ FinalPartons.cc \ FinalState.cc \ FParameter.cc \ HadronicFinalState.cc \ HeavyHadrons.cc \ Hemispheres.cc \ IdentifiedFinalState.cc \ InitialQuarks.cc \ InvMassFinalState.cc \ JetAlg.cc \ JetShape.cc \ LeadingParticlesFinalState.cc \ MergedFinalState.cc \ MissingMomentum.cc \ NeutralFinalState.cc \ NonHadronicFinalState.cc \ NonPromptFinalState.cc \ ParisiTensor.cc \ ParticleFinder.cc \ PrimaryHadrons.cc \ PrimaryParticles.cc \ PromptFinalState.cc \ Sphericity.cc \ Spherocity.cc \ TauFinder.cc \ Thrust.cc \ TriggerCDFRun0Run1.cc \ TriggerCDFRun2.cc \ TriggerUA5.cc \ UnstableFinalState.cc \ + UndressBeamLeptons.cc \ VetoedFinalState.cc \ VisibleFinalState.cc \ WFinder.cc \ ZFinder.cc diff --git a/src/Projections/UndressBeamLeptons.cc b/src/Projections/UndressBeamLeptons.cc new file mode 100644 --- /dev/null +++ b/src/Projections/UndressBeamLeptons.cc @@ -0,0 +1,36 @@ +// -*- C++ -*- +#include "Rivet/Projections/UndressBeamLeptons.hh" + +namespace Rivet { + + + int UndressBeamLeptons::compare(const Projection & p) const { + const UndressBeamLeptons & o = + dynamic_cast(p); + return cmp(_thetamax, o._thetamax) || mkNamedPCmp(o, "FS"); + } + + + void UndressBeamLeptons::project(const Event& e) { + Beam::project(e); + if ( _thetamax <= 0.0 ) return; + bool l1 = _theBeams.first.isChargedLepton(); + bool l2 = _theBeams.second.isChargedLepton(); + if ( !l1 && !l2 ) return; + FourMomentum b1 = _theBeams.first.momentum(); + FourMomentum b2 = _theBeams.second.momentum(); + Vector3 b10 = b1.vector3(); + Vector3 b20 = b2.vector3(); + for ( auto p : apply(e, "FS").particles() ) { + if ( p.pid() != PID::PHOTON ) continue; + if ( p.momentum().angle(b10) < _thetamax ) + b1 -= p.momentum(); + if ( p.momentum().angle(b20) < _thetamax ) + b2 -= p.momentum(); + } + _theBeams.first.setMomentum(b1); + _theBeams.second.setMomentum(b2); + } + + +}