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/Makefile.am b/include/Rivet/Makefile.am --- a/include/Rivet/Makefile.am +++ b/include/Rivet/Makefile.am @@ -1,197 +1,199 @@ ## 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/DISKinematics.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/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/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,123 +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; /// 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(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/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/DISLepton.cc b/src/Projections/DISLepton.cc --- a/src/Projections/DISLepton.cc +++ b/src/Projections/DISLepton.cc @@ -1,70 +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 { 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 { 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,48 +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