diff --git a/include/Rivet/Event.hh b/include/Rivet/Event.hh --- a/include/Rivet/Event.hh +++ b/include/Rivet/Event.hh @@ -1,174 +1,174 @@ // -*- 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); } + : _genevent_original(ge) + { assert(ge); _genevent = *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 //@{ /// 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 boost to the beam centre-of-mass // Vector3 beamCMSBoost() const; // /// Get the boost to the beam centre-of-mass // LorentzTransform beamCMSTransform(); /// The generated event obtained from an external event generator const GenEvent* genEvent() const { return &_genevent; } /// 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); } /// @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; //@} /// @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 { const Projection* cpp(&p); std::set::const_iterator old = _projections.find(cpp); if (old != _projections.end()) { const Projection& pRef = **old; return pcast(pRef); } // Add the projection via the Projection base class (only // possible because Event is a friend of Projection) Projection* pp = const_cast(cpp); 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/src/Core/Event.cc b/src/Core/Event.cc --- a/src/Core/Event.cc +++ b/src/Core/Event.cc @@ -1,114 +1,114 @@ // -*- C++ -*- #include "Rivet/Event.hh" #include "Rivet/Tools/BeamConstraint.hh" #include "Rivet/Projections/Beam.hh" #include "HepMC/GenEvent.h" namespace Rivet { ParticlePair Event::beams() const { return Rivet::beams(*this); } // PdgIdPair Event::beamIds() const { return pids(beams()); } double Event::sqrtS() const { return Rivet::sqrtS(beams()); } double Event::asqrtS() const { return Rivet::asqrtS(beams()); } // Vector3 Event::beamCMSBoost() const { return Rivet::beamCMSBoost(*this); } // LorentzTransform Event::beamCMSTransform() const { return Rivet::beamCMSTransform(*this); } void Event::_init(const GenEvent& ge) { // Use Rivet's preferred units if possible #ifdef HEPMC_HAS_UNITS _genevent.use_units(HepMC::Units::GEV, HepMC::Units::MM); #endif // Use the conventional alignment // _geNormAlignment(); /// @todo Filter the event to remove generator-specific particles: optional /// behaviour? Maybe disableable in an inconvenient way, e.g. with an env /// var, to communicate the appropriate distaste for this sort of truth /// analysis ;-) // Debug printout to check that copying/mangling has worked /// @todo Enable this when HepMC has been fixed to allow printing to a stream like the Rivet logger. //_genevent.print(); } // namespace { // unnamed namespace for hiding // // void _geRot180x(GenEvent& ge) { // /// @todo Use nicer iterators over HepMC particles // for (HepMC::GenEvent::particle_iterator ip = ge.particles_begin(); ip != ge.particles_end(); ++ip) { // const HepMC::FourVector& mom = (*ip)->momentum(); // (*ip)->set_momentum(HepMC::FourVector(mom.px(), -mom.py(), -mom.pz(), mom.e())); // } // /// @todo Use nicer iterators over HepMC vertices // for (HepMC::GenEvent::vertex_iterator iv = ge.vertices_begin(); iv != ge.vertices_end(); ++iv) { // const HepMC::FourVector& pos = (*iv)->position(); // (*iv)->set_position(HepMC::FourVector(pos.x(), -pos.y(), -pos.z(), pos.t())); // } // } // // } // void Event::_geNormAlignment() { // if (!_genevent.valid_beam_particles()) return; // typedef pair GPPair; // GPPair bps = _genevent.beam_particles(); // // // Rotate e+- p and ppbar to put p along +z // /// @todo Is there an e+ e- convention for longitudinal boosting, e.g. at B-factories? Different from LEP? // // if (compatible(beamids, make_pdgid_pair(ELECTRON, PROTON)) || // // compatible(beamids, make_pdgid_pair(POSITRON, PROTON)) || // // compatible(beamids, make_pdgid_pair(ANTIPROTON, PROTON)) ) { // // Log::getLog("Rivet.Event") << Log::TRACE << "May need to rotate event..." << endl; // bool rot = false; // const HepMC::GenParticle* plusgp = 0; // if (bps.first->pdg_id() != PID::PROTON || bps.second->pdg_id() != PID::PROTON) { // if (bps.first->pdg_id() == PID::PROTON) { // plusgp = bps.first; // } else if (bps.second->pdg_id() == PID::PROTON) { // plusgp = bps.second; // } // if (plusgp && plusgp->momentum().pz() < 0) { // rot = true; // } // } // // // Do the rotation // if (rot) { // if (Log::getLog("Rivet.Event").isActive(Log::TRACE)) { // Log::getLog("Rivet.Event") << Log::TRACE << "Rotating event\n" // << "Before rotation: " // << bps.first->pdg_id() << "@pz=" << bps.first->momentum().pz()/GeV << ", " // << bps.second->pdg_id() << "@pz=" << bps.second->momentum().pz()/GeV << endl; // } // _geRot180x(_genevent); // } // } const Particles& Event::allParticles() const { if (_particles.empty()) { //< assume that empty means no attempt yet made - for (const GenParticle* gp : particles(genEvent())) { + for (const GenParticlePtr gp : particles(genEvent())) { _particles += Particle(gp); } } return _particles; } double Event::weight() const { return (!_genevent.weights().empty()) ? _genevent.weights()[0] : 1.0; } } diff --git a/src/Core/Run.cc b/src/Core/Run.cc --- a/src/Core/Run.cc +++ b/src/Core/Run.cc @@ -1,155 +1,163 @@ // -*- C++ -*- #include "Rivet/Run.hh" #include "Rivet/AnalysisHandler.hh" #include "HepMC/IO_GenEvent.h" #include "Rivet/Math/MathUtils.hh" #include namespace Rivet { Run::Run(AnalysisHandler& ah) : _ah(ah), _fileweight(1.0), _xs(NAN) { } Run::~Run() { } Run& Run::setCrossSection(const double xs) { _xs = xs; return *this; } double Run::crossSection() const { return _ah.crossSection(); } Run& Run::setListAnalyses(const bool dolist) { _listAnalyses = dolist; return *this; } // Fill event and check for a bad read state bool Run::readEvent() { /// @todo Clear rather than new the GenEvent object per-event? _evt.reset(new GenEvent()); if (_io->rdstate() != 0 || !_io->fill_next_event(_evt.get()) ) { Log::getLog("Rivet.Run") << Log::DEBUG << "Read failed. End of file?" << endl; return false; } // Rescale event weights by file-level weight, if scaling is non-trivial if (!fuzzyEquals(_fileweight, 1.0)) { for (size_t i = 0; i < (size_t) _evt->weights().size(); ++i) { _evt->weights()[i] *= _fileweight; } } return true; } - + // Fill event and check for a bad read state --- to skip, maybe HEPMC3 will have a better way bool Run::skipEvent() { if (_io->rdstate() != 0 || !_io->fill_next_event(_evt.get()) ) { Log::getLog("Rivet.Run") << Log::DEBUG << "Read failed. End of file?" << endl; return false; } return true; } bool Run::openFile(const std::string& evtfile, double weight) { // Set current weight-scaling member _fileweight = weight; // Set up HepMC input reader objects if (evtfile == "-") { _io.reset(new HepMC::IO_GenEvent(std::cin)); } else { // Ignore the HepMC::IO_GenEvent(filename, ios) constructor, since it's only available from HepMC 2.4 _istr.reset(new std::fstream(evtfile.c_str(), std::ios::in)); _io.reset(new HepMC::IO_GenEvent(*_istr)); } if (_io->rdstate() != 0) { Log::getLog("Rivet.Run") << Log::ERROR << "Read error on file " << evtfile << endl; return false; } return true; } bool Run::init(const std::string& evtfile, double weight) { if (!openFile(evtfile, weight)) return false; // Read first event to define run conditions bool ok = readEvent(); if (!ok) return false; + #if HEPMC_VERSION_CODE >= 300000 + if (_evt->particles().empty()) { + #else if (_evt->particles_size() == 0) { + #endif Log::getLog("Rivet.Run") << Log::ERROR << "Empty first event." << endl; return false; } // Initialise AnalysisHandler with beam information from first event _ah.init(*_evt); // Set cross-section from command line if (!std::isnan(_xs)) { Log::getLog("Rivet.Run") << Log::DEBUG << "Setting user cross-section = " << _xs << " pb" << endl; _ah.setCrossSection(_xs); } // List the chosen & compatible analyses if requested if (_listAnalyses) { foreach (const std::string& ana, _ah.analysisNames()) { cout << ana << endl; } } return true; } bool Run::processEvent() { // Set cross-section if found in event and not from command line #ifdef HEPMC_HAS_CROSS_SECTION if (std::isnan(_xs) && _evt->cross_section()) { + #if HEPMC_VERSION_CODE >= 300000 + const double xs = _evt->cross_section()->cross_section; ///< in pb + #else const double xs = _evt->cross_section()->cross_section(); ///< in pb + #endif Log::getLog("Rivet.Run") << Log::DEBUG << "Setting cross-section = " << xs << " pb" << endl; _ah.setCrossSection(xs); } #endif // Complain about absence of cross-section if required! if (_ah.needCrossSection() && !_ah.hasCrossSection()) { Log::getLog("Rivet.Run") << Log::ERROR << "Total cross-section needed for at least one of the analyses. " << "Please set it (on the command line with '-x' if using the 'rivet' program)" << endl; return false; } // Analyze event _ah.analyze(*_evt); return true; } bool Run::finalize() { _evt.reset(); _istr.reset(); _io.reset(); if (!std::isnan(_xs)) _ah.setCrossSection(_xs); _ah.finalize(); return true; } }