diff --git a/src/Core/Jet.cc b/src/Core/Jet.cc
--- a/src/Core/Jet.cc
+++ b/src/Core/Jet.cc
@@ -1,191 +1,191 @@
 #include "Rivet/Jet.hh"
 #include "Rivet/Tools/Cuts.hh"
 #include "Rivet/Tools/ParticleName.hh"
 #include "Rivet/Tools/Logging.hh"
 #include "Rivet/Tools/ParticleIdUtils.hh"
 
 namespace Rivet {
 
 
   Jet& Jet::clear() {
     _momentum = FourMomentum();
     _pseudojet.reset(0,0,0,0);
     _particles.clear();
     return *this;
   }
 
 
   Jet& Jet::setState(const FourMomentum& mom, const Particles& particles, const Particles& tags) {
     clear();
     _momentum = mom;
     _pseudojet = fastjet::PseudoJet(mom.px(), mom.py(), mom.pz(), mom.E());
     _particles = particles;
     _tags = tags;
     return *this;
   }
 
 
   Jet& Jet::setState(const fastjet::PseudoJet& pj, const Particles& particles, const Particles& tags) {
     clear();
     _pseudojet = pj;
     _momentum = FourMomentum(pj.e(), pj.px(), pj.py(), pj.pz());
     _particles = particles;
     _tags = tags;
     // if (_particles.empty()) {
     //   for (const fastjet::PseudoJet pjc : _pseudojet.constituents()) {
     //     // If there is no attached user info, we can't create a meaningful particle, so skip
     //     if (!pjc.has_user_info<RivetFJInfo>()) continue;
     //     const RivetFJInfo& fjinfo = pjc.user_info<RivetFJInfo>();
     //     // Don't add ghosts to the particles list
     //     if (fjinfo.isGhost) continue;
     //     // Otherwise construct a Particle from the PseudoJet, preferably from an associated GenParticle
     //     ?if (fjinfo.genParticle != NULL) {
     //       _particles.push_back(Particle(fjinfo.genParticle));
     //     } else {
     //       if (fjinfo.pid == 0) continue; // skip if there is a null PID entry in the FJ info
     //       const FourMomentum pjcmom(pjc.e(), pjc.px(), pjc.py(), pjc.pz());
     //       _particles.push_back(Particle(fjinfo.pid, pjcmom));
     //     }
     //   }
     // }
     return *this;
   }
 
 
   Jet& Jet::setParticles(const Particles& particles) {
     _particles = particles;
     return *this;
   }
 
 
   bool Jet::containsParticle(const Particle& particle) const {
-    #if HEPMC_VERSION_CODE >= 300000
+    #if HEPMC_VERSION_CODE >= 3000000
     const int idcode = particle.genParticle()->id();
     for (const Particle& p : particles()) {
       if (p.genParticle()->id() == idcode) return true;
     }
     #else
     const int barcode = particle.genParticle()->barcode();
     for (const Particle& p : particles()) {
       if (p.genParticle()->barcode() == barcode) return true;
     }
     #endif
     return false;
   }
 
 
   bool Jet::containsParticleId(PdgId pid) const {
     for (const Particle& p : particles()) {
       if (p.pid() == pid) return true;
     }
     return false;
   }
 
 
   bool Jet::containsParticleId(const vector<PdgId>& pids) const {
     for (const Particle& p : particles()) {
       for (PdgId pid : pids) {
         if (p.pid() == pid) return true;
       }
     }
     return false;
   }
 
 
   /// @todo Jet::containsMatch(Matcher m) { ... if m(pid) return true; ... }
 
 
   Jet& Jet::transformBy(const LorentzTransform& lt) {
     _momentum = lt.transform(_momentum);
     for (Particle& p : _particles) p.transformBy(lt);
     for (Particle& t : _tags) t.transformBy(lt);
     _pseudojet.reset(_momentum.px(), _momentum.py(), _momentum.pz(), _momentum.E()); //< lose ClusterSeq etc.
     return *this;
   }
 
 
   double Jet::neutralEnergy() const {
     double e_neutral = 0.0;
     for (const Particle& p : particles()) {
       if (p.charge3() == 0) e_neutral += p.E();
     }
     return e_neutral;
   }
 
 
   double Jet::hadronicEnergy() const {
     double e_hadr = 0.0;
     for (const Particle& p : particles()) {
       if (isHadron(p)) e_hadr += p.E();
     }
     return e_hadr;
   }
 
 
   bool Jet::containsCharm(bool include_decay_products) const {
     for (const Particle& p : particles()) {
       if (p.abspid() == PID::CQUARK) return true;
       if (isHadron(p) && hasCharm(p)) return true;
       if (include_decay_products) {
         const HepMC::GenVertexPtr gv = p.genParticle()->production_vertex();
         if (gv) {
           for (const GenParticlePtr pi : Rivet::particles(gv, HepMC::ancestors)) {
             const PdgId pid2 = pi->pdg_id();
             if (PID::isHadron(pid2) && PID::hasCharm(pid2)) return true;
           }
         }
       }
     }
     return false;
   }
 
 
   bool Jet::containsBottom(bool include_decay_products) const {
     for (const Particle& p : particles()) {
       if (p.abspid() == PID::BQUARK) return true;
       if (isHadron(p) && hasBottom(p)) return true;
       if (include_decay_products) {
         const HepMC::GenVertexPtr gv = p.genParticle()->production_vertex();
         if (gv) {
           for (const GenParticlePtr pi : Rivet::particles(gv, HepMC::ancestors)) {
             const PdgId pid2 = pi->pdg_id();
             if (PID::isHadron(pid2) && PID::hasBottom(pid2)) return true;
           }
         }
       }
     }
     return false;
   }
 
 
   Particles Jet::tags(const Cut& c) const {
     return filter_select(tags(), c);
   }
 
   Particles Jet::bTags(const Cut& c) const {
     Particles rtn;
     for (const Particle& tp : tags()) {
       if (hasBottom(tp) && c->accept(tp)) rtn.push_back(tp);
     }
     return rtn;
   }
 
   Particles Jet::cTags(const Cut& c) const {
     Particles rtn;
     for (const Particle& tp : tags()) {
       /// @todo Is making b and c tags exclusive the right thing to do?
       if (hasCharm(tp) && !hasBottom(tp) && c->accept(tp)) rtn.push_back(tp);
     }
     return rtn;
   }
 
   Particles Jet::tauTags(const Cut& c) const {
     Particles rtn;
     for (const Particle& tp : tags()) {
       if (isTau(tp) && c->accept(tp)) rtn.push_back(tp);
     }
     return rtn;
   }
 
 
 }
diff --git a/src/Core/Particle.cc b/src/Core/Particle.cc
--- a/src/Core/Particle.cc
+++ b/src/Core/Particle.cc
@@ -1,266 +1,258 @@
 #include "Rivet/Particle.hh"
 #include "Rivet/Tools/Cuts.hh"
 #include "Rivet/Tools/ParticleIdUtils.hh"
 
 namespace Rivet {
 
 
   Particle& Particle::transformBy(const LorentzTransform& lt) {
     _momentum = lt.transform(_momentum);
     return *this;
   }
 
 
   bool Particle::isVisible() const {
     // Charged particles are visible
     if ( PID::threeCharge(pid()) != 0 ) return true;
     // Neutral hadrons are visible
     if ( PID::isHadron(pid()) ) return true;
     // Photons are visible
     if ( pid() == PID::PHOTON ) return true;
     // Gluons are visible (for parton level analyses)
     if ( pid() == PID::GLUON ) return true;
     // Everything else is invisible
     return false;
   }
 
 
   bool Particle::isStable() const {
     return genParticle() != NULL &&
       genParticle()->status() == 1 &&
       genParticle()->end_vertex() == NULL;
   }
 
 
   vector<Particle> Particle::parents(const Cut& c) const {
     vector<Particle> rtn;
-    /// @todo Remove this const mess crap when HepMC doesn't suck
+    #if HEPMC_VERSION_CODE >= 3000000
+    for (const GenParticlePtr gp : genParticle()->parents()) {
+    const Particle p(gp);
+    if (c != Cuts::OPEN && !c->accept(p)) continue;
+      rtn += p;
+    }
+    #else
     GenVertexPtr gv = const_cast<GenVertexPtr>( genParticle()->production_vertex() );
     if (gv == NULL) return rtn;
-    /// @todo Would like to do this, but the range objects are broken
-    // foreach (const GenParticlePtr gp, gv->particles(HepMC::children))
-    //   rtn += Particle(gp);
     for (GenVertex::particle_iterator it = gv->particles_begin(HepMC::parents); it != gv->particles_end(HepMC::parents); ++it) {
-      const Particle p(*it);
+    const Particle p(*it);
       if (c != Cuts::OPEN && !c->accept(p)) continue;
       rtn += p;
     }
+    #endif
     return rtn;
   }
 
 
   vector<Particle> Particle::children(const Cut& c) const {
     vector<Particle> rtn;
     if (isStable()) return rtn;
-    /// @todo Remove this const mess crap when HepMC doesn't suck
+    #if HEPMC_VERSION_CODE >= 3000000
+    for (const GenParticlePtr gp : genParticle()->children()) {
+      const Particle p(gp);
+      if (c != Cuts::OPEN && !c->accept(p)) continue;
+      rtn += p;
+    }
+    #else
     GenVertexPtr gv = const_cast<GenVertexPtr>( genParticle()->end_vertex() );
     if (gv == NULL) return rtn;
-    /// @todo Would like to do this, but the range objects are broken
-    // foreach (const GenParticlePtr gp, gv->particles(HepMC::children))
-    //   rtn += Particle(gp);
     for (GenVertex::particle_iterator it = gv->particles_begin(HepMC::children); it != gv->particles_end(HepMC::children); ++it) {
       const Particle p(*it);
       if (c != Cuts::OPEN && !c->accept(p)) continue;
       rtn += p;
     }
+    #endif
     return rtn;
   }
 
 
   /// @todo Insist that the current particle is post-hadronization, otherwise throw an exception?
   /// @todo Use recursion through replica-avoiding functions to avoid bookkeeping duplicates
   vector<Particle> Particle::allDescendants(const Cut& c, bool remove_duplicates) const {
     vector<Particle> rtn;
     if (isStable()) return rtn;
-    /// @todo Remove this const mess crap when HepMC doesn't suck
+    #if HEPMC_VERSION_CODE >= 3000000
+    for (const GenParticlePtr gp : genParticle()->descendants()) {
+      const Particle p(gp);
+      if (c != Cuts::OPEN && !c->accept(p)) continue;
+      if (remove_duplicates) {
+        bool dup = false;
+        for (const GenParticlePtr gp2 : gp->children()) {
+          if (gp->pid() == gp2->pid()) { dup = true; break; }
+        }
+        if (dup) continue;
+      }
+      rtn += p;
+    }
+    #else
     GenVertexPtr gv = const_cast<GenVertexPtr>( genParticle()->end_vertex() );
     if (gv == NULL) return rtn;
     /// @todo Would like to do this, but the range objects are broken
     // foreach (const GenParticlePtr gp, gv->particles(HepMC::descendants))
     for (GenVertex::particle_iterator it = gv->particles_begin(HepMC::descendants); it != gv->particles_end(HepMC::descendants); ++it) {
       const Particle p(*it);
       if (c != Cuts::OPEN && !c->accept(p)) continue;
       if (remove_duplicates && (*it)->end_vertex() != NULL) {
         // size_t n = 0; ///< @todo Only remove 1-to-1 duplicates?
         bool dup = false;
         /// @todo Yuck, HepMC
         for (GenVertex::particle_iterator it2 = (*it)->end_vertex()->particles_begin(HepMC::children); it2 != (*it)->end_vertex()->particles_end(HepMC::children); ++it2) {
           // n += 1; if (n > 1) break;
           if ((*it)->pdg_id() == (*it2)->pdg_id()) { dup = true; break; }
         }
         if (dup) continue;
       }
       rtn += p;
     }
+    #endif
     return rtn;
   }
 
 
   /// @todo Insist that the current particle is post-hadronization, otherwise throw an exception?
   vector<Particle> Particle::stableDescendants(const Cut& c) const {
     vector<Particle> rtn;
     if (isStable()) return rtn;
-    /// @todo Remove this const mess crap when HepMC doesn't suck
+    #if HEPMC_VERSION_CODE >= 3000000
+    for (const GenParticlePtr gp : genParticle()->descendants()) {
+      const Particle p(gp);
+      if (!p.isStable()) continue;
+      if (c != Cuts::OPEN && !c->accept(p)) continue;
+      rtn += p;
+    }
+    #else
     GenVertexPtr gv = const_cast<GenVertexPtr>( genParticle()->end_vertex() );
     if (gv == NULL) return rtn;
-    /// @todo Would like to do this, but the range objects are broken
-    // foreach (const GenParticlePtr gp, gv->particles(HepMC::descendants))
     for (GenVertex::particle_iterator it = gv->particles_begin(HepMC::descendants); it != gv->particles_end(HepMC::descendants); ++it) {
       // if ((*it)->status() != 1 || (*it)->end_vertex() != NULL) continue;
       const Particle p(*it);
       if (!p.isStable()) continue;
       if (c != Cuts::OPEN && !c->accept(p)) continue;
       rtn += p;
     }
+    #endif
     return rtn;
   }
 
 
   double Particle::flightLength() const {
     if (isStable()) return -1;
-    if (genParticle() == NULL) return 0;
-    if (genParticle()->production_vertex() == NULL) return 0;
+    if (!genParticle()) return 0;
+    if (!genParticle()->production_vertex()) return 0;
     const HepMC::FourVector v1 = genParticle()->production_vertex()->position();
     const HepMC::FourVector v2 = genParticle()->end_vertex()->position();
     return sqrt(sqr(v2.x()-v1.x()) + sqr(v2.y()-v1.y()) + sqr(v2.z()-v1.z()));
   }
 
 
   bool Particle::hasParent(PdgId pid) const {
     return _hasRelativeWith(HepMC::parents, hasPID(pid));
   }
 
   bool Particle::hasParentWith(const Cut& c) const {
     return hasParentWith([&](const Particle& p){return c->accept(p);});
   }
 
 
   bool Particle::hasAncestor(PdgId pid) const {
     return _hasRelativeWith(HepMC::ancestors, hasPID(pid));
   }
 
   bool Particle::hasAncestorWith(const Cut& c) const {
     return hasAncestorWith([&](const Particle& p){return c->accept(p);});
   }
 
 
   bool Particle::fromBottom() const {
     return _hasRelativeWith(HepMC::ancestors, [](const Particle& p){
         return p.genParticle()->status() == 2 && p.isHadron() && p.hasBottom();
       });
-    // const GenVertexPtr prodVtx = genParticle()->production_vertex();
-    // if (prodVtx == NULL) return false;
-    // foreach (const GenParticlePtr ancestor, particles(prodVtx, HepMC::ancestors)) {
-    //   const PdgId pid = ancestor->pdg_id();
-    //   if (ancestor->status() == 2 && (PID::isHadron(pid) && PID::hasBottom(pid))) return true;
-    // }
-    // return false;
   }
 
 
   bool Particle::fromCharm() const {
     return _hasRelativeWith(HepMC::ancestors, [](const Particle& p){
         return p.genParticle()->status() == 2 && p.isHadron() && p.hasCharm();
       });
-    // const GenVertexPtr prodVtx = genParticle()->production_vertex();
-    // if (prodVtx == NULL) return false;
-    // foreach (const GenParticlePtr ancestor, particles(prodVtx, HepMC::ancestors)) {
-    //   const PdgId pid = ancestor->pdg_id();
-    //   if (ancestor->status() == 2 && (PID::isHadron(pid) && PID::hasCharm(pid) && !PID::hasBottom(pid))) return true;
-    // }
-    // return false;
   }
 
 
   bool Particle::fromHadron() const {
     return _hasRelativeWith(HepMC::ancestors, [](const Particle& p){
         return p.genParticle()->status() == 2 && p.isHadron();
       });
-    // const GenVertexPtr prodVtx = genParticle()->production_vertex();
-    // if (prodVtx == NULL) return false;
-    // foreach (const GenParticlePtr ancestor, particles(prodVtx, HepMC::ancestors)) {
-    //   const PdgId pid = ancestor->pdg_id();
-    //   if (ancestor->status() == 2 && PID::isHadron(pid)) return true;
-    // }
-    // return false;
   }
 
 
   bool Particle::fromTau(bool prompt_taus_only) const {
     if (prompt_taus_only && fromHadron()) return false;
     return _hasRelativeWith(HepMC::ancestors, [](const Particle& p){
         return p.genParticle()->status() == 2 && isTau(p);
       });
-    // const GenVertexPtr prodVtx = genParticle()->production_vertex();
-    // if (prodVtx == NULL) return false;
-    // foreach (const GenParticlePtr ancestor, particles(prodVtx, HepMC::ancestors)) {
-    //   const PdgId pid = ancestor->pdg_id();
-    //   if (ancestor->status() == 2 && abs(pid) == PID::TAU) return true;
-    // }
-    // return false;
   }
 
 
-  // bool Particle::fromDecay() const {
-  //   const GenVertexPtr prodVtx = genParticle()->production_vertex();
-  //   if (prodVtx == NULL) return false;
-  //   foreach (const GenParticlePtr ancestor, particles(prodVtx, HepMC::ancestors)) {
-  //     const PdgId pid = ancestor->pdg_id();
-  //     if (ancestor->status() == 2 && (PID::isHadron(pid) || abs(pid) == PID::TAU)) return true;
-  //   }
-  //   return false;
-  // }
-
-
   bool Particle::isPrompt(bool allow_from_prompt_tau, bool allow_from_prompt_mu) const {
     if (genParticle() == NULL) return false; // no HepMC connection, give up! Throw UserError exception?
     const GenVertexPtr prodVtx = genParticle()->production_vertex();
     if (prodVtx == NULL) return false; // orphaned particle, has to be assume false
+    #if HEPMC_VERSION_CODE >= 3000000
+    const pair<GenParticlePtr, GenParticlePtr> beams = std::make_pair(prodVtx->parent_event()->beams()[0], prodVtx->parent_event()->beams()[1]);
+    #else
     const pair<GenParticlePtr, GenParticlePtr> beams = prodVtx->parent_event()->beam_particles();
-
+    #endif
     /// @todo Would be nicer to be able to write this recursively up the chain, exiting as soon as a parton or string/cluster is seen
     for (const GenParticlePtr ancestor : Rivet::particles(prodVtx, HepMC::ancestors)) {
       const PdgId pid = ancestor->pdg_id();
       if (ancestor->status() != 2) continue; // no non-standard statuses or beams to be used in decision making
       if (ancestor == beams.first || ancestor == beams.second) continue; // PYTHIA6 uses status 2 for beams, I think... (sigh)
       if (PID::isParton(pid)) continue; // PYTHIA6 also uses status 2 for some partons, I think... (sigh)
       if (PID::isHadron(pid)) return false; // prompt particles can't be from hadron decays
       if (abs(pid) == PID::TAU && abspid() != PID::TAU && !allow_from_prompt_tau) return false; // allow or ban particles from tau decays (permitting tau copies)
       if (abs(pid) == PID::MUON && abspid() != PID::MUON && !allow_from_prompt_mu) return false; // allow or ban particles from muon decays (permitting muon copies)
     }
     return true;
   }
 
 
 
   ///////////////////////
 
 
 
   string to_str(const Particle& p) {
     string pname;
     try {
       pname = PID::toParticleName(p.pid());
     } catch (...) {
       pname = "PID=" + to_str(p.pid());
     }
     stringstream out;
     out << pname << " @ " << p.momentum() << " GeV";
     return out.str();
   }
 
 
   string to_str(const ParticlePair& pair) {
     stringstream out;
     out << "[" << pair.first << ", " << pair.second << "]";
     // out << "["
     //     << PID::toParticleName(pair.first.pid()) << " @ "
     //     << pair.first.momentum().E()/GeV << " GeV, "
     //     << PID::toParticleName(pair.second.pid()) << " @ "
     //     << pair.second.momentum().E()/GeV << " GeV]";
     return out.str();
   }
 
 
 
 }
diff --git a/src/Core/Run.cc b/src/Core/Run.cc
--- a/src/Core/Run.cc
+++ b/src/Core/Run.cc
@@ -1,163 +1,163 @@
 // -*- C++ -*-
 #include "Rivet/Run.hh"
 #include "Rivet/AnalysisHandler.hh"
 #include "HepMC/IO_GenEvent.h"
 #include "Rivet/Math/MathUtils.hh"
 #include <limits>
 
 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 HEPMC_VERSION_CODE >= 3000000
     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
+    #if HEPMC_VERSION_CODE >= 3000000
       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;
   }
 
 
 
 
 }