Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-09-18 09:31:23

0001 // -*- C++ -*-
0002 #ifndef RIVET_Event_HH
0003 #define RIVET_Event_HH
0004 
0005 #include "Rivet/Config/RivetCommon.hh"
0006 #include "Rivet/Particle.hh"
0007 #include "Rivet/Projection.hh"
0008 
0009 namespace Rivet {
0010 
0011 
0012   /// @brief Representation of a HepMC event, and enabler of Projection caching
0013   ///
0014   /// Event is a concrete class representing an generated event in Rivet. It is
0015   /// constructed given a HepMC::GenEvent, a pointer to which is kept by the
0016   /// Event object throughout its lifetime. The user must therefore make sure
0017   /// that the corresponding HepMC::GenEvent will persist at least as long as
0018   /// the Event object.
0019   ///
0020   /// In addition to the HepMC::GenEvent object the Event also keeps track of
0021   /// all Projection objects which have been applied to the Event so far.
0022   class Event {
0023   public:
0024 
0025     /// @name Constructors and destructors.
0026     /// @{
0027 
0028     /// Constructor from a HepMC GenEvent pointer
0029     ///
0030     /// @note Although Rivet is an analysis system and should not modify the
0031     /// *physics* of provided events, it needs to ensure that the events are
0032     /// "sane" for analysis purposes, including fixing units and potentially
0033     /// orientations and unphysical graph structures. The passed event pointer
0034     /// is hence not const.
0035     Event(GenEvent* ge, const vector<size_t>& weightindices={}) //, bool strip=false)
0036       : _weightIndices(weightindices), _genevent(*ge)
0037     {
0038       if (!ge) throw UserError("User provided a null GenEvent pointer for analysis");
0039       // if (strip) _strip(_genevent);
0040       _fixGenEvent();
0041     }
0042 
0043     /// Constructor from a HepMC GenEvent reference
0044     ///
0045     /// @note See the pointer-based constructor for discussion of the
0046     /// non-constness of the passed GenEvent.
0047     ///
0048     /// @deprecated HepMC uses pointers, so we should talk to HepMC via pointers: no need to duplicate.
0049     Event(GenEvent& ge, const vector<size_t>& weightindices={}) //, bool strip=false)
0050       : _weightIndices(weightindices), _genevent(ge)
0051     {
0052       // if (strip) _strip(_genevent);
0053       _fixGenEvent();
0054     }
0055 
0056     /// Copy constructor
0057     Event(const Event& e)
0058       : _weightIndices(e._weightIndices),
0059         _genevent(e._genevent) {  }
0060 
0061     /// @}
0062 
0063 
0064     /// @name Major event properties
0065     /// @{
0066 
0067     /// The generated HepMC event obtained from an external event generator
0068     const GenEvent& hepmcEvent() const { return _genevent; }
0069     // /// @brief The generated HepMC event obtained from an external event generator
0070     // ///
0071     // /// Convenience alias for hepmcEventPtr()
0072     // const GenEvent& hepmc() const { return _genevent; }
0073 
0074     /// The generated HepMC event pointer obtained from an external event generator
0075     const GenEvent* hepmcEventPtr() const { return &_genevent; }
0076     /// @brief The generated HepMC event pointer obtained from an external event generator
0077     ///
0078     /// Backward-compatibility alias for hepmcEventPtr()
0079     const GenEvent* genEvent() const { return &_genevent; }
0080 
0081     /// Get the beam particles
0082     ParticlePair beams() const;
0083 
0084     /// Get the beam centre-of-mass energy
0085     double sqrtS() const;
0086 
0087     /// Get the beam centre-of-mass energy per nucleon
0088     double asqrtS() const;
0089 
0090     /// @}
0091 
0092 
0093     /// @name Access to event particles
0094     /// @{
0095 
0096     /// All the raw GenEvent particles, wrapped in Rivet::Particle objects
0097     const Particles& allParticles() const;
0098 
0099     /// @brief All the raw GenEvent particles, wrapped in Rivet::Particle objects, but with a Cut applied
0100     ///
0101     /// @note Due to the cut, this returns by value, i.e. involves an expensive copy
0102     inline Particles allParticles(const Cut& c) const {
0103       return select(allParticles(), c);
0104     }
0105 
0106     /// @brief All the raw GenEvent particles, wrapped in Rivet::Particle objects, but with a selection function applied
0107     ///
0108     /// @note Due to the cut, this returns by value, i.e. involves an expensive copy
0109     template <typename FN>
0110     inline Particles allParticles(const FN& f) const {
0111       return select(allParticles(), f);
0112     }
0113 
0114     /// @brief The generation weights associated with the event
0115     ///
0116     /// This vector contains only the weights selected by the AnalysisHandler,
0117     /// so is of constant length for all events.
0118     ///
0119     /// @note Excluded generation weights, and their associated weight names
0120     /// if set, can be accessed via the GenEvent pointer hepmcEventPtr().
0121     std::valarray<double> weights() const;
0122 
0123     /// @brief The generation cross-sections associated with the event
0124     ///
0125     /// This vector contains only the cross-sections selected by the AnalysisHandler,
0126     /// so is of constant length for all events.
0127     ///
0128     /// @note Excluded cross-sections, and their associated names if set,
0129     /// can be accessed via the GenEvent pointer hepmcEventPtr().
0130     std::vector<std::pair<double, double>> crossSections() const;
0131 
0132     /// @}
0133 
0134     /// @name Projection running
0135     /// @{
0136 
0137     /// @brief Add a projection @a p to this Event.
0138     ///
0139     /// If an equivalent Projection has been applied before, the
0140     /// Projection::project(const Event&) of @a p is not called and a reference
0141     /// to the previous equivalent projection is returned. If no previous
0142     /// Projection was found, the Projection::project(const Event&) of @a p is
0143     /// called and a reference to @a p is returned.
0144     ///
0145     /// @todo Can make this non-templated, since only cares about ptr to Projection base class
0146     ///
0147     /// @note Comparisons here are by direct pointer comparison, because
0148     /// equivalence is guaranteed if pointers are equal, and inequivalence
0149     /// guaranteed if they aren't, thanks to the ProjectionHandler registry
0150     template <typename PROJ>
0151     const PROJ& applyProjection(PROJ& p) const {
0152       static bool docaching = getEnvParam("RIVET_CACHE_PROJECTIONS", true);
0153       if (docaching) {
0154         MSG_TRACE("Applying projection " << &p << " (" << p.name() << ") -> comparing to projections " << _projections);
0155         // First search for this projection *or an equivalent* in the already-executed list
0156         const Projection* cpp(&p);
0157         /// @note Currently using reint cast to integer type to bypass operator==(Proj*, Proj*)
0158         // std::set<const Projection*>::const_iterator old = _projections.find(cpp);
0159         std::set<const Projection*>::const_iterator old = std::begin(_projections);
0160         std::uintptr_t recpp = reinterpret_cast<std::uintptr_t>(cpp);
0161         for (; old != _projections.end(); ++old)
0162           if (reinterpret_cast<std::uintptr_t>(*old) == recpp) break;
0163         if (old != _projections.end()) {
0164           MSG_TRACE("Equivalent projection found -> returning already-run projection " << *old);
0165           const Projection& pRef = **old;
0166           return pcast<PROJ>(pRef);
0167         }
0168         MSG_TRACE("No equivalent projection in the already-run list -> projecting now");
0169       } else {
0170         MSG_TRACE("Applying projection " << &p << " (" << p.name() << ") WITHOUT projection caching & comparison");
0171       }
0172       // If this one hasn't been run yet on this event, run it and add to the list
0173       Projection* pp = const_cast<Projection*>(&p);
0174       pp->_isValid = true;
0175       pp->project(*this);
0176       if (docaching) _projections.insert(pp);
0177       return p;
0178     }
0179 
0180 
0181     /// @brief Add a projection @a p to this Event by pointer.
0182     template <typename PROJ>
0183     const PROJ& applyProjection(PROJ* pp) const {
0184       if (!pp) throw Error("Event::applyProjection(PROJ*): Projection pointer is null.");
0185       return applyProjection(*pp);
0186     }
0187 
0188     /// @}
0189 
0190 
0191   private:
0192 
0193     /// Tweak the GenEvent to Rivet's expected standards if necessary
0194     void _fixGenEvent();
0195 
0196     /// Get a Log object for Event
0197     Log& getLog() const;
0198 
0199 
0200     // /// @brief Remove uninteresting or unphysical particles in the
0201     // /// GenEvent to speed up searches
0202     // ///
0203     // /// @todo Remove!
0204     // void _strip(GenEvent& ge);
0205 
0206     // /// @brief Convert the GenEvent to use conventional alignment
0207     // ///
0208     // /// For example, FHerwig only produces DIS events in the unconventional
0209     // /// hadron-lepton orientation and has to be corrected for DIS analysis
0210     // /// portability.
0211     // void _geNormAlignment();
0212 
0213     /// @brief The indices of the selected weights, as instructed to the AnalysisHandler.
0214     ///
0215     /// The user can (de-)select weights and the AnalysisHandler knows about the subset
0216     /// that match the specifications.
0217     const std::vector<size_t> _weightIndices;
0218 
0219     /// The GenEvent used by Rivet analysis projections etc.
0220     const GenEvent& _genevent;
0221 
0222     /// All the GenEvent particles, wrapped as Rivet::Particles
0223     ///
0224     /// @note To be populated lazily, hence mutability
0225     mutable Particles _particles;
0226 
0227     /// The set of Projection objects applied so far
0228     mutable std::set<ConstProjectionPtr> _projections;
0229 
0230     /// A cached set of event weights (a single unit weight if the original weight vector was empty)
0231     mutable std::valarray<double> _weights;
0232 
0233     /// A cached set of event cross-section & errors (a pair of zeros if the original weight vector was empty)
0234     mutable std::vector<std::pair<double,double>> _xsecs;
0235   };
0236 
0237 
0238 }
0239 
0240 #endif