diff --git a/include/Rivet/Projections/CentralityEstimator.hh b/include/Rivet/Projections/CentralityEstimator.hh --- a/include/Rivet/Projections/CentralityEstimator.hh +++ b/include/Rivet/Projections/CentralityEstimator.hh @@ -1,80 +1,83 @@ // -*- C++ -*- #ifndef RIVET_CentralityEstimator_HH #define RIVET_CentralityEstimator_HH #include "Rivet/Projection.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Event.hh" #include "HepMC/GenEvent.h" namespace Rivet { /** @brief Base class for projections giving the value of an observable sensitive to the centrality of a collision. @author Leif Lönnblad The centrality of a collision is not really an observable, but the concept is anyway often used in the heavy ion community as if it were just that. This base class can be used to provide a an estimator for the centrality by projecting down to a single number which then canbe used by a CentralityGroup object to select a histogram to be filled with another observable depending on centrality percentile. The eztimate() should be a non-negative number with large values indicating a higher overlap than small ones. A negative value indicates that the centrality estimate could not be calculated. In the best of all worlds the centrality estimator should be a proper hadron-level observable corrected for detector effects, however, this base class only returns the inverse of the impact_parameter member of the GenHeavyIon object in an GenEvent if present and zero otherwise. */ class CentralityEstimator : public Projection { public: /// Constructor. CentralityEstimator() : _estimate(-1.0) {} /// Clone on the heap. DEFAULT_RIVET_PROJ_CLONE(CentralityEstimator); protected: /// Perform the projection on the Event void project(const Event& e) { _estimate = -1.0; - const HepMC::HeavyIon * hi = e.genEvent()->heavy_ion(); - if ( hi ) _estimate = hi->impact_parameter() > 0.0? - 1.0/hi->impact_parameter(): numeric_limits::max(); + #if HEPMC_VERSION_CODE >= 3000000 + const HepMC::GenHeavyIonPtr hi = e.genEvent()->heavy_ion(); + if (hi) _estimate = hi->impact_parameter > 0.0 ? 1.0/hi->impact_parameter: DBL_MAX; + #else + const HepMC::HeavyIon* hi = e.genEvent()->heavy_ion(); + if (hi) _estimate = hi->impact_parameter() > 0.0 ? 1.0/hi->impact_parameter(): DBL_MAX; + #endif } /// Compare projections int compare(const Projection& p) const { return mkNamedPCmp(p, "CentEst"); } public: /// The value of the centrality estimate. double estimate() const { return _estimate; } protected: /// The value of the centrality estimate. double _estimate; }; } #endif - diff --git a/src/Projections/Beam.cc b/src/Projections/Beam.cc --- a/src/Projections/Beam.cc +++ b/src/Projections/Beam.cc @@ -1,130 +1,140 @@ // -*- C++ -*- #include "Rivet/Config/RivetCommon.hh" #include "Rivet/Tools/Logging.hh" #include "Rivet/Projections/Beam.hh" namespace Rivet { ParticlePair beams(const Event& e) { + #if HEPMC_VERSION_CODE >= 3000000 + assert(e.genEvent()->particles().size() >= 2); + const vector beams = e.genEvent()->beams(); + if (beams.size() == 2 && beams[0] && beams[1]) { + return ParticlePair{beams[0], beams[0]}; + } else { + return ParticlePair{Particle(PID::ANY, FourMomentum()), Particle(PID::ANY, FourMomentum())}; + } + #else assert(e.genEvent()->particles_size() >= 2); if (e.genEvent()->valid_beam_particles()) { pair beams = e.genEvent()->beam_particles(); assert(beams.first && beams.second); return ParticlePair{beams.first, beams.second}; } else if (e.genEvent()->barcode_to_particle(1) && e.genEvent()->barcode_to_particle(2)) { return ParticlePair{e.genEvent()->barcode_to_particle(1), e.genEvent()->barcode_to_particle(2)}; } return ParticlePair{Particle(PID::ANY, FourMomentum()), Particle(PID::ANY, FourMomentum())}; + #endif } double sqrtS(const FourMomentum& pa, const FourMomentum& pb) { const double mom1 = pa.pz(); const double e1 = pa.E(); const double mom2 = pb.pz(); const double e2 = pb.E(); double sqrts = sqrt( sqr(e1+e2) - sqr(mom1+mom2) ); return sqrts; } double asqrtS(const FourMomentum& pa, const FourMomentum& pb) { const static double MNUCLEON = 939*MeV; //< nominal nucleon mass return sqrtS(pa/(pa.mass()/MNUCLEON), pb/(pb.mass()/MNUCLEON)); } double asqrtS(const ParticlePair& beams) { return sqrtS(beams.first.mom()/nuclA(beams.first), beams.second.mom()/nuclA(beams.second)); } FourMomentum acmsBoostVec(const FourMomentum& pa, const FourMomentum& pb) { const static double MNUCLEON = 939*MeV; //< nominal nucleon mass const double na = pa.mass()/MNUCLEON, nb = pb.mass()/MNUCLEON; return cmsBoostVec(pa/na, pb/nb); } FourMomentum acmsBoostVec(const ParticlePair& beams) { return cmsBoostVec(beams.first.mom()/nuclA(beams.first), beams.second.mom()/nuclA(beams.second)); } Vector3 cmsBetaVec(const FourMomentum& pa, const FourMomentum& pb) { // const Vector3 rtn = (pa.p3() + pb.p3()) / (pa.E() + pb.E()); const Vector3 rtn = (pa + pb).betaVec(); return rtn; } Vector3 acmsBetaVec(const FourMomentum& pa, const FourMomentum& pb) { const static double MNUCLEON = 939*MeV; //< nominal nucleon mass const Vector3 rtn = cmsBetaVec(pa/(pa.mass()/MNUCLEON), pb/(pb.mass()/MNUCLEON)); return rtn; } Vector3 acmsBetaVec(const ParticlePair& beams) { const Vector3 rtn = cmsBetaVec(beams.first.mom()/nuclA(beams.first), beams.second.mom()/nuclA(beams.second)); return rtn; } Vector3 cmsGammaVec(const FourMomentum& pa, const FourMomentum& pb) { // const Vector3 rtn = (pa + pb).gammaVec(); const double gamma = (pa.E() + pb.E()) / sqrt( sqr(pa.mass()) + sqr(pb.mass()) + 2*(pa.E()*pb.E() - dot(pa.p3(), pb.p3())) ); const Vector3 rtn = gamma * (pa.p3() + pb.p3()).unit(); return rtn; } Vector3 acmsGammaVec(const FourMomentum& pa, const FourMomentum& pb) { const static double MNUCLEON = 939*MeV; //< nominal nucleon mass Vector3 rtn = cmsGammaVec(pa/(pa.mass()/MNUCLEON), pb/(pb.mass()/MNUCLEON)); return rtn; } Vector3 acmsGammaVec(const ParticlePair& beams) { Vector3 rtn = cmsGammaVec(beams.first.mom()/nuclA(beams.first), beams.second.mom()/nuclA(beams.second)); return rtn; } LorentzTransform cmsTransform(const FourMomentum& pa, const FourMomentum& pb) { /// @todo Automatically choose to construct from beta or gamma according to which is more precise? return LorentzTransform::mkFrameTransformFromGamma(cmsGammaVec(pa, pb)); } LorentzTransform acmsTransform(const FourMomentum& pa, const FourMomentum& pb) { /// @todo Automatically choose to construct from beta or gamma according to which is more precise? return LorentzTransform::mkFrameTransformFromGamma(acmsGammaVec(pa, pb)); } LorentzTransform acmsTransform(const ParticlePair& beams) { return LorentzTransform::mkFrameTransformFromGamma(acmsGammaVec(beams)); } ///////////////////////////////////////////// void Beam::project(const Event& e) { _theBeams = Rivet::beams(e); MSG_DEBUG("Beam particles = " << _theBeams << " => sqrt(s) = " << sqrtS()/GeV << " GeV"); } FourVector Beam::pv() const { HepMC::FourVector v1, v2; const ParticlePair bpair = beams(); if (bpair.first.genParticle() && bpair.first.genParticle()->end_vertex()) v1 = bpair.first.genParticle()->end_vertex()->position(); if (bpair.second.genParticle() && bpair.second.genParticle()->end_vertex()) v2 = bpair.second.genParticle()->end_vertex()->position(); const FourVector rtn = (v1 == v2) ? FourVector(v1.t(), v1.x(), v1.y(), v1.z()) : FourVector(); MSG_DEBUG("Beam PV 4-position = " << rtn); return rtn; } } diff --git a/src/Projections/ChargedFinalState.cc b/src/Projections/ChargedFinalState.cc --- a/src/Projections/ChargedFinalState.cc +++ b/src/Projections/ChargedFinalState.cc @@ -1,48 +1,43 @@ // -*- C++ -*- #include "Rivet/Projections/ChargedFinalState.hh" namespace Rivet { ChargedFinalState::ChargedFinalState(const FinalState& fsp) { setName("ChargedFinalState"); addProjection(fsp, "FS"); } ChargedFinalState::ChargedFinalState(const Cut& c) { setName("ChargedFinalState"); addProjection(FinalState(c), "FS"); } ChargedFinalState::ChargedFinalState(double mineta, double maxeta, double minpt) { setName("ChargedFinalState"); addProjection(FinalState(mineta, maxeta, minpt), "FS"); } int ChargedFinalState::compare(const Projection& p) const { return mkNamedPCmp(p, "FS"); } -} -namespace { - inline bool chargedParticleFilter(const Rivet::Particle& p) { - return Rivet::PID::threeCharge(p.pdgId()) == 0; - } -} -namespace Rivet { void ChargedFinalState::project(const Event& e) { const FinalState& fs = applyProjection(e, "FS"); - _theParticles.clear(); - std::remove_copy_if(fs.particles().begin(), fs.particles().end(), - std::back_inserter(_theParticles), chargedParticleFilter); + // _theParticles.clear(); + // std::remove_copy_if(fs.particles().begin(), fs.particles().end(), + // std::back_inserter(_theParticles), + // [](const Rivet::Particle& p) { return p.charge3() == 0; }); + _theParticles = filter_select(fs.particles(), isCharged); MSG_DEBUG("Number of charged final-state particles = " << _theParticles.size()); if (getLog().isActive(Log::TRACE)) { - for (vector::iterator p = _theParticles.begin(); p != _theParticles.end(); ++p) { - MSG_TRACE("Selected: " << p->pdgId() << ", charge = " << PID::threeCharge(p->pdgId())/3.0); + for (const Particle& p : _theParticles) { + MSG_TRACE("Selected: " << p.pid() << ", charge = " << p.charge()); } } } } diff --git a/src/Projections/DISFinalState.cc b/src/Projections/DISFinalState.cc --- a/src/Projections/DISFinalState.cc +++ b/src/Projections/DISFinalState.cc @@ -1,28 +1,28 @@ // -*- C++ -*- #include "Rivet/Projections/DISFinalState.hh" namespace Rivet { void DISFinalState::project(const Event& e) { const DISKinematics& diskin = applyProjection(e, "Kinematics"); const LorentzTransform hcmboost = (_boosttype == HCM) ? diskin.boostHCM() : diskin.boostBreit(); const DISLepton& dislep = diskin.applyProjection(e, "Lepton"); const FinalState& fs = dislep.applyProjection(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(); - foreach (const Particle& p, fs.particles()) { + const GenParticlePtr dislepGP = dislep.out().genParticle(); + for (const Particle& p : fs.particles()) { if (p.genParticle() != dislepGP) { ///< Ensure that we skip the DIS lepton Particle temp(p); temp.setMomentum(hcmboost.transform(temp.momentum())); _theParticles.push_back(temp); } } } } diff --git a/src/Projections/FinalPartons.cc b/src/Projections/FinalPartons.cc --- a/src/Projections/FinalPartons.cc +++ b/src/Projections/FinalPartons.cc @@ -1,46 +1,45 @@ // -*- C++ -*- #include "Rivet/Projections/FinalPartons.hh" #include "Rivet/Tools/RivetHepMC.hh" namespace Rivet { - bool FinalPartons::accept(const Particle& p) const { + bool FinalPartons::accept(const Particle& p) const { - // reject if *not* a parton - if (!isParton(p)) - return false; + // Reject if *not* a parton + if (!isParton(p)) + return false; - // accept partons if they end on a standard hadronization vertex - if (p.genParticle()->end_vertex() != NULL && - p.genParticle()->end_vertex()->id() == 5) - return true; + // Accept partons if they end on a standard hadronization vertex + if (p.genParticle()->end_vertex() != NULL && + p.genParticle()->end_vertex()->id() == 5) + return true; - // reject if p has a parton child. - foreach (const Particle& c, p.children()) { - if (isParton(c)) - return false; - } - - // reject if from a hadron or tau decay. - if (p.fromDecay()) - return false; - - return _cuts->accept(p); + // Reject if p has a parton child. + for (const Particle& c : p.children()) { + if (isParton(c)) + return false; } + // Reject if from a hadron or tau decay. + if (p.fromDecay()) + return false; - void FinalPartons::project(const Event& e) { - _theParticles.clear(); + return _cuts->accept(p); + } - foreach (const GenParticle* gp, Rivet::particles(e.genEvent())) { - if (!gp) continue; - const Particle p(gp); - if (accept(p)) _theParticles.push_back(p); - } + void FinalPartons::project(const Event& e) { + _theParticles.clear(); + + for (const GenParticlePtr gp : Rivet::particles(e.genEvent())) { + if (!gp) continue; + const Particle p(gp); + if (accept(p)) _theParticles.push_back(p); } + } } diff --git a/src/Projections/FinalState.cc b/src/Projections/FinalState.cc --- a/src/Projections/FinalState.cc +++ b/src/Projections/FinalState.cc @@ -1,82 +1,82 @@ // -*- C++ -*- #include "Rivet/Projections/FinalState.hh" namespace Rivet { FinalState::FinalState(const Cut & c) : ParticleFinder(c) { setName("FinalState"); const bool open = (c == Cuts::open()); MSG_TRACE("Check for open FS conditions: " << std::boolalpha << open); if (!open) addProjection(FinalState(), "OpenFS"); } /// @deprecated, keep for backwards compatibility for now. FinalState::FinalState(double mineta, double maxeta, double minpt) { setName("FinalState"); const bool openpt = isZero(minpt); const bool openeta = (mineta <= -MAXDOUBLE && maxeta >= MAXDOUBLE); MSG_TRACE("Check for open FS conditions:" << std::boolalpha << " eta=" << openeta << ", pt=" << openpt); if (openpt && openeta) { _cuts = Cuts::open(); } else { addProjection(FinalState(), "OpenFS"); if (openeta) _cuts = (Cuts::pT >= minpt); else if ( openpt ) _cuts = Cuts::etaIn(mineta, maxeta); else _cuts = (Cuts::etaIn(mineta, maxeta) && Cuts::pT >= minpt); } } int FinalState::compare(const Projection& p) const { const FinalState& other = dynamic_cast(p); return _cuts == other._cuts ? EQUIVALENT : UNDEFINED; } void FinalState::project(const Event& e) { _theParticles.clear(); // Handle "open FS" special case if (_cuts == Cuts::OPEN) { MSG_TRACE("Open FS processing: should only see this once per event (" << e.genEvent()->event_number() << ")"); - for (const GenParticle* p : Rivet::particles(e.genEvent())) { + for (const GenParticlePtr p : Rivet::particles(e.genEvent())) { if (p->status() == 1) { - MSG_TRACE("FS GV = " << p->production_vertex()); - _theParticles.push_back(Particle(*p)); + MSG_TRACE("FS GV = " << p->production_vertex()->position()); + _theParticles.push_back(Particle(p)); } } return; } // If this is not itself the "open" FS, base the calculations on the open FS' results /// @todo In general, we'd like to calculate a restrictive FS based on the most restricted superset FS. const Particles allstable = applyProjection(e, "OpenFS").particles(); for (const Particle& p : allstable) { const bool passed = accept(p); MSG_TRACE("Choosing: ID = " << p.pid() << ", pT = " << p.pT()/GeV << " GeV" << ", eta = " << p.eta() << ": result = " << std::boolalpha << passed); if (passed) _theParticles.push_back(p); } MSG_TRACE("Number of final-state particles = " << _theParticles.size()); } /// Decide if a particle is to be accepted or not. bool FinalState::accept(const Particle& p) const { // Not having status == 1 should never happen! assert(p.genParticle() == NULL || p.genParticle()->status() == 1); return _cuts->accept(p); } } diff --git a/src/Projections/HeavyHadrons.cc b/src/Projections/HeavyHadrons.cc --- a/src/Projections/HeavyHadrons.cc +++ b/src/Projections/HeavyHadrons.cc @@ -1,64 +1,64 @@ // -*- C++ -*- #include "Rivet/Projections/HeavyHadrons.hh" namespace Rivet { void HeavyHadrons::project(const Event& e) { _theParticles.clear(); _theBs.clear(); _theCs.clear(); /// @todo Allow user to choose whether primary or final HF hadrons are to be returned const Particles& unstables = applyProjection(e, "UFS").particles(); - foreach (const Particle& p, unstables) { + for (const Particle& p : unstables) { // Exclude non-b/c-hadrons if (!isHadron(p)) continue; if (!hasCharm(p) && !hasBottom(p)) continue; MSG_DEBUG("Found a heavy (b or c) unstable hadron: " << p.pid()); // An unbound, or undecayed status 2 hadron: this is weird, but I guess is allowed... if (!p.genParticle() || !p.genParticle()->end_vertex()) { MSG_DEBUG("Heavy hadron " << p.pid() << " with no GenParticle or decay found"); _theParticles.push_back(p); if (hasBottom(p)) _theBs.push_back(p); else _theCs.push_back(p); continue; } // There are descendants -- check them for b or c content /// @todo What about charm hadrons coming from bottom hadron decays? - const vector children = particles_out(p.genParticle(), HepMC::children); + const vector children = Rivet::particles(p.genParticle(), HepMC::children); if (hasBottom(p)) { bool has_b_child = false; - foreach (const GenParticle* p2, children) { + for (const GenParticlePtr p2 : children) { if (PID::hasBottom(p2->pdg_id())) { has_b_child = true; break; } } if (!has_b_child) { _theParticles.push_back(p); _theBs.push_back(p); } } else if (hasCharm(p)) { bool has_c_child = false; - foreach (const GenParticle* p2, children) { + for (const GenParticlePtr p2 : children) { if (PID::hasCharm(p2->pdg_id())) { has_c_child = true; break; } } if (!has_c_child) { _theParticles.push_back(p); _theCs.push_back(p); } } } MSG_DEBUG("Num b hadrons = " << _theBs.size() << ", num c hadrons = " << _theCs.size() << ", total = " << _theParticles.size()); } } diff --git a/src/Tools/Cuts.cc b/src/Tools/Cuts.cc --- a/src/Tools/Cuts.cc +++ b/src/Tools/Cuts.cc @@ -1,449 +1,449 @@ #include "Rivet/Tools/Cuts.hh" #include "Rivet/Particle.hh" #include "Rivet/Jet.hh" #include "Rivet/Math/Vectors.hh" +#include "Rivet/Tools/RivetHepMC.hh" #include "fastjet/PseudoJet.hh" -#include "HepMC/SimpleVector.h" /// @todo Identify what can go into anonymous namespace namespace Rivet { // Base for all wrapper classes that translate ClassToCheck to Cuttable class CuttableBase { public: virtual double getValue(Cuts::Quantity) const = 0; virtual ~CuttableBase() {} }; // Cuttables can be directly passed to @ref _accept template <> bool CutBase::accept(const CuttableBase& t) const { return _accept(t); } // Open Cut singleton class Open_Cut : public CutBase { public: bool operator==(const Cut& c) const { std::shared_ptr cc = dynamic_pointer_cast(c); return bool(cc); } protected: // open cut accepts everything bool _accept(const CuttableBase&) const { return true; } }; const Cut& Cuts::open() { // Only ever need one static open cut object static const Cut open = std::make_shared(); return open; } // Constants for convenient access const Cut& Cuts::OPEN = Cuts::open(); const Cut& Cuts::NOCUT = Cuts::open(); // Cut constructor for == class Cut_Eq : public CutBase { public: Cut_Eq(const Cuts::Quantity qty, double val) : _qty(qty), _val(val) {} bool operator==(const Cut& c) const { std::shared_ptr cc = dynamic_pointer_cast(c); return cc && _qty == cc->_qty && _val == cc->_val; } protected: bool _accept(const CuttableBase& o) const { return o.getValue(_qty) == _val; } private: Cuts::Quantity _qty; int _val; }; // Cut constructor for != class Cut_NEq : public CutBase { public: Cut_NEq(const Cuts::Quantity qty, double val) : _qty(qty), _val(val) {} bool operator==(const Cut& c) const { std::shared_ptr cc = dynamic_pointer_cast(c); return cc && _qty == cc->_qty && _val == cc->_val; } protected: bool _accept(const CuttableBase& o) const { return o.getValue(_qty) != _val; } private: Cuts::Quantity _qty; int _val; }; // Cut constructor for >= class Cut_GtrEq : public CutBase { public: Cut_GtrEq(const Cuts::Quantity qty, double val) : _qty(qty), _val(val) {} bool operator==(const Cut& c) const { std::shared_ptr cc = dynamic_pointer_cast(c); return cc && _qty == cc->_qty && _val == cc->_val; } protected: bool _accept(const CuttableBase& o) const { return o.getValue(_qty) >= _val; } private: Cuts::Quantity _qty; double _val; }; // Cut constructor for < class Cut_Less : public CutBase { public: Cut_Less(const Cuts::Quantity qty, const double val) : _qty(qty), _val(val) {} bool operator==(const Cut& c) const { std::shared_ptr cc = dynamic_pointer_cast(c); return cc && _qty == cc->_qty && _val == cc->_val; } protected: bool _accept(const CuttableBase& o) const { return o.getValue(_qty) < _val; } private: Cuts::Quantity _qty; double _val; }; // Cut constructor for > class Cut_Gtr : public CutBase { public: Cut_Gtr(const Cuts::Quantity qty, const double val) : _qty(qty), _val(val) {} bool operator==(const Cut& c) const { std::shared_ptr cc = dynamic_pointer_cast(c); return cc && _qty == cc->_qty && _val == cc->_val; } protected: bool _accept(const CuttableBase& o) const { return o.getValue(_qty) > _val; } private: Cuts::Quantity _qty; double _val; }; // Cut constructor for <= class Cut_LessEq : public CutBase { public: Cut_LessEq(const Cuts::Quantity qty, const double val) : _qty(qty), _val(val) {} bool operator==(const Cut& c) const { std::shared_ptr cc = dynamic_pointer_cast(c); return cc && _qty == cc->_qty && _val == cc->_val; } protected: bool _accept(const CuttableBase& o) const { return o.getValue(_qty) <= _val; } private: Cuts::Quantity _qty; double _val; }; template inline Cut make_cut(T t) { return std::make_shared(t); } Cut operator == (Cuts::Quantity qty, double n) { return make_cut(Cut_Eq(qty, n)); } Cut operator != (Cuts::Quantity qty, double n) { return make_cut(Cut_NEq(qty, n)); } Cut operator < (Cuts::Quantity qty, double n) { return make_cut(Cut_Less(qty, n)); } Cut operator >= (Cuts::Quantity qty, double n) { return make_cut(Cut_GtrEq(qty, n)); } Cut operator <= (Cuts::Quantity qty, double n) { return make_cut(Cut_LessEq(qty, n)); } Cut operator > (Cuts::Quantity qty, double n) { return make_cut(Cut_Gtr(qty, n)); } Cut Cuts::range(Cuts::Quantity qty, double m, double n) { if (m > n) std::swap(m,n); return (qty >= m) && (qty < n); } ////////////// /// Combiners /// AND, OR, NOT, and XOR objects for combining cuts class CutsOr : public CutBase { public: CutsOr(const Cut& c1, const Cut& c2) : cut1(c1), cut2(c2) {} bool operator==(const Cut& c) const { std::shared_ptr cc = dynamic_pointer_cast(c); return cc && ( ( cut1 == cc->cut1 && cut2 == cc->cut2 ) || ( cut1 == cc->cut2 && cut2 == cc->cut1 )); } protected: bool _accept(const CuttableBase& o) const { return cut1->accept(o) || cut2->accept(o); } private: const Cut cut1; const Cut cut2; }; class CutsAnd : public CutBase { public: CutsAnd(const Cut& c1, const Cut& c2) : cut1(c1), cut2(c2) {} bool operator==(const Cut& c) const { std::shared_ptr cc = dynamic_pointer_cast(c); return cc && ( ( cut1 == cc->cut1 && cut2 == cc->cut2 ) || ( cut1 == cc->cut2 && cut2 == cc->cut1 )); } protected: bool _accept(const CuttableBase& o) const { return cut1->accept(o) && cut2->accept(o); } private: const Cut cut1; const Cut cut2; }; class CutInvert : public CutBase { public: CutInvert(const Cut& c1) : cut(c1) {} bool operator==(const Cut& c) const { std::shared_ptr cc = dynamic_pointer_cast(c); return cc && cut == cc->cut; } protected: bool _accept(const CuttableBase& o) const { return !cut->accept(o); } private: const Cut cut; }; class CutsXor : public CutBase { public: CutsXor(const Cut& c1, const Cut& c2) : cut1(c1), cut2(c2) {} bool operator==(const Cut& c) const { std::shared_ptr cc = dynamic_pointer_cast(c); return cc && ( ( cut1 == cc->cut1 && cut2 == cc->cut2 ) || ( cut1 == cc->cut2 && cut2 == cc->cut1 )); } protected: bool _accept(const CuttableBase& o) const { bool A_and_B = cut1->accept(o) && cut2->accept(o); bool A_or_B = cut1->accept(o) || cut2->accept(o); return A_or_B && (! A_and_B); } private: const Cut cut1; const Cut cut2; }; //////////// ///Operators Cut operator && (const Cut& aptr, const Cut& bptr) { return make_cut(CutsAnd(aptr, bptr)); } Cut operator || (const Cut& aptr, const Cut& bptr) { return make_cut(CutsOr(aptr, bptr)); } Cut operator ! (const Cut& cptr) { return make_cut(CutInvert(cptr)); } Cut operator & (const Cut& aptr, const Cut& bptr) { return make_cut(CutsAnd(aptr, bptr)); } Cut operator | (const Cut& aptr, const Cut& bptr) { return make_cut(CutsOr(aptr, bptr)); } Cut operator ~ (const Cut& cptr) { return make_cut(CutInvert(cptr)); } Cut operator ^ (const Cut& aptr, const Cut& bptr) { return make_cut(CutsXor(aptr, bptr)); } /////////////////////// /// Cuts // Non-functional Cuttable template class. Must be specialized below. template class Cuttable; // Non-cuttables need to be wrapped into a Cuttable first #define SPECIALISE_ACCEPT(TYPENAME) \ template <> \ bool CutBase::accept(const TYPENAME& t) const { \ return _accept(Cuttable(t)); \ } \ inline void qty_not_found() { throw Exception("Missing implementation for a Cuts::Quantity."); } template<> class Cuttable : public CuttableBase { public: Cuttable(const Particle& p) : p_(p) {} double getValue(Cuts::Quantity qty) const { switch ( qty ) { case Cuts::pT: return p_.pT(); case Cuts::Et: return p_.Et(); case Cuts::mass: return p_.mass(); case Cuts::rap: return p_.rap(); case Cuts::absrap: return p_.absrap(); case Cuts::eta: return p_.eta(); case Cuts::abseta: return p_.abseta(); case Cuts::phi: return p_.phi(); case Cuts::pid: return p_.pid(); case Cuts::abspid: return p_.abspid(); case Cuts::charge: return p_.charge(); case Cuts::abscharge: return p_.abscharge(); case Cuts::charge3: return p_.charge3(); case Cuts::abscharge3: return p_.abscharge3(); default: qty_not_found(); } return -999.; } private: const Particle& p_; }; SPECIALISE_ACCEPT(Particle) template<> class Cuttable : public CuttableBase { public: Cuttable(const FourMomentum& fm) : fm_(fm) {} double getValue(Cuts::Quantity qty) const { switch ( qty ) { case Cuts::pT: return fm_.pT(); case Cuts::Et: return fm_.Et(); case Cuts::mass: return fm_.mass(); case Cuts::rap: return fm_.rap(); case Cuts::absrap: return fm_.absrap(); case Cuts::eta: return fm_.eta(); case Cuts::abseta: return fm_.abseta(); case Cuts::phi: return fm_.phi(); default: qty_not_found(); } return -999.; } private: const FourMomentum& fm_; }; SPECIALISE_ACCEPT(FourMomentum) template<> class Cuttable : public CuttableBase { public: Cuttable(const Jet& jet) : jet_(jet) {} double getValue(Cuts::Quantity qty) const { switch ( qty ) { case Cuts::pT: return jet_.momentum().pT(); case Cuts::Et: return jet_.momentum().Et(); case Cuts::mass: return jet_.momentum().mass(); case Cuts::rap: return jet_.momentum().rapidity(); case Cuts::absrap: return std::abs(jet_.momentum().rapidity()); case Cuts::eta: return jet_.momentum().pseudorapidity(); case Cuts::abseta: return std::abs(jet_.momentum().pseudorapidity()); case Cuts::phi: return jet_.momentum().phi(); default: qty_not_found(); } return -999.; } private: const Jet& jet_; }; SPECIALISE_ACCEPT(Jet) template<> class Cuttable : public CuttableBase { public: Cuttable(const fastjet::PseudoJet& pjet) : pjet_(pjet) {} double getValue(Cuts::Quantity qty) const { switch ( qty ) { case Cuts::pT: return pjet_.perp(); case Cuts::Et: return pjet_.Et(); case Cuts::mass: return pjet_.m(); case Cuts::rap: return pjet_.rap(); case Cuts::absrap: return std::abs(pjet_.rap()); case Cuts::eta: return pjet_.eta(); case Cuts::abseta: return std::abs(pjet_.eta()); case Cuts::phi: return pjet_.phi(); default: qty_not_found(); } return -999.; } private: const fastjet::PseudoJet& pjet_; }; SPECIALISE_ACCEPT(fastjet::PseudoJet) template<> class Cuttable : public CuttableBase { public: Cuttable(const HepMC::FourVector& vec) : vec_(vec) {} double getValue(Cuts::Quantity qty) const { switch ( qty ) { case Cuts::pT: return vec_.perp(); /// @todo case Cuts::Et: return vec_.perp(); case Cuts::mass: return vec_.m(); case Cuts::rap: return 0.5*std::log((vec_.t()+vec_.z())/(vec_.t()-vec_.z())); case Cuts::absrap: return std::abs(getValue(Cuts::rap)); case Cuts::eta: return vec_.pseudoRapidity(); case Cuts::abseta: return std::abs(vec_.pseudoRapidity()); case Cuts::phi: return vec_.phi(); default: qty_not_found(); } return -999.; } private: const HepMC::FourVector& vec_; }; SPECIALISE_ACCEPT(HepMC::FourVector) }