diff --git a/include/Rivet/Tools/Cmp.hh b/include/Rivet/Tools/Cmp.hh --- a/include/Rivet/Tools/Cmp.hh +++ b/include/Rivet/Tools/Cmp.hh @@ -1,303 +1,307 @@ // -*- C++ -*- #ifndef RIVET_Cmp_HH #define RIVET_Cmp_HH #include "Rivet/Config/RivetCommon.hh" #include "Rivet/Projection.hh" #include "Cmp.fhh" #include namespace Rivet { /// Helper class when checking the ordering of two objects. /// /// Cmp is a helper class to be used when checking the ordering of two /// objects. When implicitly converted to an integer the value will be /// negative if the two objects used in the constructor are ordered and /// positive if they are not. Zero will be returned if they are equal. /// /// The main usage of the Cmp class is if several variables should be /// checked for ordering in which case several Cmp objects can be /// combined as follows: cmp(a1, a2) || cmp(b1, b2) || cmp(c1, /// c2) where cmp is a global function for easy creation of Cmp /// objects. template class Cmp final { public: /// @name Standard constructors etc. //@{ /// The default constructor. Cmp(const T& t1, const T& t2) : _value(CmpState::UNDEF), _objects(&t1, &t2) { } /// The copy constructor. template Cmp(const Cmp& x) : _value(x._value), _objects(nullptr, nullptr) { } /// The assignment operator. template const Cmp& operator=(const Cmp& x) { _value = x; return *this; } //@} public: /// Automatically convert to an enum. operator CmpState() const { _compare(); return _value; } /// If this state is equivalent, set this state to the state of \a c. template const Cmp& operator||(const Cmp& c) const { _compare(); if (_value == CmpState::EQ) _value = c; return *this; } private: /// Perform the actual comparison if necessary. void _compare() const { if (_value == CmpState::UNDEF) { std::less l; if ( l(*_objects.first, *_objects.second) ) _value = CmpState::LT; else if ( l(*_objects.second, *_objects.first) ) _value = CmpState::GT; else _value = CmpState::EQ; } /// @todo NEQ? } /// The state of this object. mutable CmpState _value; /// The objects to be compared. const pair _objects; }; /// @brief Specialization of Cmp for checking the ordering of two @a {Projection}s. /// /// Specialization of the Cmp helper class to be used when checking the /// ordering of two Projection objects. When implicitly converted to an /// integer the value will be negative if the two objects used in the /// constructor are ordered and positive if they are not. Zero will be /// returned if they are equal. This specialization uses directly the /// virtual compare() function in the Projection class. /// /// The main usage of the Cmp class is if several variables should be /// checked for ordering in which case several Cmp objects can be /// combined as follows: cmp(a1, a2) || cmp(b1, b2) || cmp(c1, /// c2) where cmp is a global function for easy creation of Cmp /// objects. template <> class Cmp final { public: /// @name Standard constructors and destructors. //@{ /// The default constructor. Cmp(const Projection& p1, const Projection& p2) : _value(CmpState::UNDEF), _objects(&p1, &p2) { } /// The copy constructor. template Cmp(const Cmp& x) : _value(x), _objects(nullptr, nullptr) { } /// The assignment operator. template const Cmp& operator=(const Cmp& x) { _value = x; return *this; } //@} public: /// Automatically convert to an enum. operator CmpState() const { _compare(); return _value; } /// If this state is equivalent, set this state to the state of \a c. template const Cmp& operator||(const Cmp& c) const { _compare(); if (_value == CmpState::EQ) _value = c; return *this; } private: /// Perform the actual comparison if necessary. void _compare() const { if (_value == CmpState::UNDEF) { const std::type_info& id1 = typeid(*_objects.first); const std::type_info& id2 = typeid(*_objects.second); if (id1.before(id2)) _value = CmpState::LT; else if (id2.before(id1)) _value = CmpState::GT; else { - _value = _objects.first->compare(*_objects.second); + //_value = _objects.first->compare(*_objects.second); + CmpState cmps = _objects.first->compare(*_objects.second); + if (cmps == CmpState::UNDEF) _value = CmpState::LT; + else if (cmps == CmpState::NEQ) _value = CmpState::GT; + else _value = cmps; } } } private: /// The state of this object. mutable CmpState _value; /// The objects to be compared. const pair _objects; }; /// @brief Specialization of Cmp for checking the ordering of two floating point numbers. /// /// When implicitly converted to an integer the value will be negative if the /// two objects used in the constructor are ordered and positive if they are /// not. Zero will be returned if they are equal. This specialization uses the /// Rivet fuzzyEquals function to indicate equivalence protected from /// numerical precision effects. /// /// The main usage of the Cmp class is if several variables should be /// checked for ordering in which case several Cmp objects can be /// combined as follows: cmp(a1, a2) || cmp(b1, b2) || cmp(c1, /// c2) where cmp is a global function for easy creation of Cmp /// objects. template <> class Cmp final { public: /// @name Standard constructors and destructors. //@{ /// The default constructor. Cmp(const double p1, const double p2) : _value(CmpState::UNDEF), _numA(p1), _numB(p2) { } /// The copy constructor. template Cmp(const Cmp& x) : _value(x), _numA(0.0), _numB(0.0) { } /// The assignment operator. template const Cmp& operator=(const Cmp& x) { _value = x; return *this; } //@} public: /// Automatically convert to an enum. operator CmpState() const { _compare(); return _value; } /// If this state is equivalent, set this state to the state of \a c. template const Cmp& operator||(const Cmp& c) const { _compare(); if (_value == CmpState::EQ) _value = c; return *this; } private: /// Perform the actual comparison if necessary. void _compare() const { if (_value == CmpState::UNDEF) { if (fuzzyEquals(_numA,_numB)) _value = CmpState::EQ; else if (_numA < _numB) _value = CmpState::LT; else _value = CmpState::GT; } } private: /// The state of this object. mutable CmpState _value; /// The objects to be compared. const double _numA, _numB; }; /////////////////////////////////////////////////////////////////// /// Global helper function for easy creation of Cmp objects. template inline Cmp cmp(const T& t1, const T& t2) { return Cmp(t1, t2); } /// Typedef for Cmp using PCmp = Cmp; /// Global helper function for easy creation of Cmp objects. inline Cmp pcmp(const Projection& p1, const Projection& p2) { return Cmp(p1, p2); } /// Global helper function for easy creation of Cmp objects from /// two parent projections and their common name for the projection to be compared. inline Cmp pcmp(const Projection& parent1, const Projection& parent2, const string& pname) { return Cmp(parent1.getProjection(pname), parent2.getProjection(pname)); } /// Global helper function for easy creation of Cmp objects from /// two parent projections and their common name for the projection to be compared. /// This version takes one parent as a pointer. inline Cmp pcmp(const Projection* parent1, const Projection& parent2, const string& pname) { assert(parent1); return Cmp(parent1->getProjection(pname), parent2.getProjection(pname)); } /// Global helper function for easy creation of Cmp objects from /// two parent projections and their common name for the projection to be compared. /// This version takes one parent as a pointer. inline Cmp pcmp(const Projection& parent1, const Projection* parent2, const string& pname) { assert(parent2); return Cmp(parent1.getProjection(pname), parent2->getProjection(pname)); } /// Global helper function for easy creation of Cmp objects from /// two parent projections and their common name for the projection to be compared. inline Cmp pcmp(const Projection* parent1, const Projection* parent2, const string& pname) { assert(parent1); assert(parent2); return Cmp(parent1->getProjection(pname), parent2->getProjection(pname)); } } #endif diff --git a/src/Core/Projection.cc b/src/Core/Projection.cc --- a/src/Core/Projection.cc +++ b/src/Core/Projection.cc @@ -1,62 +1,63 @@ // -*- 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") { addPdgIdPair(PID::ANY, PID::ANY); } Projection:: ~Projection() = default; 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 CmpState cmpst = compare(p); - // const bool cmp = (cmpst == CmpState::LT || cmpst == CmpState::NEQ || cmpst == CmpState::UNDEF); - const bool cmp = (cmpst != CmpState::GT); + //const bool cmp = (cmpst == CmpState::LT || cmpst == CmpState::NEQ || cmpst == CmpState::UNDEF); + //const bool cmp = (cmpst != CmpState::GT); + const bool cmp = (cmpst == CmpState::LT || cmpst == CmpState::UNDEF); 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 << '\n'; 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/FinalState.cc b/src/Projections/FinalState.cc --- a/src/Projections/FinalState.cc +++ b/src/Projections/FinalState.cc @@ -1,83 +1,84 @@ // -*- C++ -*- #include "Rivet/Projections/FinalState.hh" namespace Rivet { FinalState::FinalState(const Cut& c) : ParticleFinder(c) { setName("FinalState"); const bool isopen = (c == Cuts::open()); MSG_TRACE("Check for open FS conditions: " << std::boolalpha << isopen); if (!isopen) declare(FinalState(), "OpenFS"); } FinalState::FinalState(const FinalState& fsp, const Cut& c) : ParticleFinder(c) { setName("FinalState"); MSG_TRACE("Registering base FSP as 'PrevFS'"); declare(fsp, "PrevFS"); } CmpState FinalState::compare(const Projection& p) const { const FinalState& other = dynamic_cast(p); // First check if there is a PrevFS and it it matches if (hasProjection("PrevFS") != other.hasProjection("PrevFS")) return CmpState::UNDEF; if (hasProjection("PrevFS")) { const PCmp prevcmp = mkPCmp(other, "PrevFS"); if (prevcmp != CmpState::EQ) return prevcmp; } // Then check the extra cuts const bool cutcmp = _cuts == other._cuts; MSG_TRACE(_cuts << " VS " << other._cuts << " -> EQ == " << std::boolalpha << cutcmp); - if (!cutcmp) return CmpState::NEQ; + if (!cutcmp) return CmpState::UNDEF; + //if (!cutcmp) return CmpState::NEQ; // Checks all passed: these FSes are equivalent return CmpState::EQ; } void FinalState::project(const Event& e) { _theParticles.clear(); // Handle "open FS" special case, which should not/cannot recurse 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())) { if (p->status() == 1) { //MSG_TRACE("FS GV = " << p->production_vertex()); _theParticles.push_back(Particle(*p)); } } MSG_TRACE("Number of open-FS selected particles = " << _theParticles.size()); return; } // Base the calculation on PrevFS if available, otherwise OpenFS /// @todo In general, we'd like to calculate a restrictive FS based on the most restricted superset FS. const Particles& allstable = applyProjection(e, (hasProjection("PrevFS") ? "PrevFS" : "OpenFS")).particles(); MSG_TRACE("Beginning Cuts selection"); 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); } }