diff --git a/.hgignore b/.hgignore --- a/.hgignore +++ b/.hgignore @@ -1,72 +1,80 @@ ~$ ^(anas|test|run).*$ \.orig$ \.(o|Po|lo|Plo|la|a|so|dylib|pyc|tar\.bz2|tar\.gz|fifo|hepmc)$ ^config\.guess$ ^config\.status$ ^config\.sub$ ^config\.log$ ^configure$ ^depcomp$ ^compile$ ^ar-lib$ ^install-sh$ ^INSTALL$ ^libtool$ ^test-driver$ ^ltmain\.sh$ ^m4/ltoptions\.m4$ ^m4/ltsugar\.m4$ ^m4/ltversion\.m4$ ^m4/lt~obsolete\.m4$ ^missing$ ^autom4te\.cache$ ^include/Rivet/Config/stamp-h.$ ^tmp$ ^rivet.pc$ Makefile$ Makefile\.in$ \.(pdf|toc|bbl|blg|sty|bib|html|tex)$ Doxyfile$ \.libs$ \.deps$ #Rivet.yoda #.*/Rivet\.yoda ^bin/rivet-buildplugin$ ^bin/rivet-config$ ^bin/rivet-nopy$ ^aclocal\.m4$ ^pyext/build$ ^pyext/rivet/core.cpp$ ^pyext/rivet/rivetwrap\.py$ ^pyext/rivet/rivetwrap_wrap\.cc$ ^pyext/rivet/fix-out-of-source$ ^pyext/setup\.py$ ^rivetenv\.(sh|csh)$ ^test/test(Api|Cmp|NaN)$ ^include/Rivet/Config/DummyConfig\.hh\.in$ ^include/Rivet/Config/BuildOptions\.hh$ ^include/Rivet/Config/DummyConfig\.hh$ ^include/Rivet/Config/RivetConfig\.hh$ ^doc/.*\.(log|aux|out)$ ^doc/ana ^test/log$ ^test/.*\.log$ ^test/out\.yoda$ ^test/NaN.aida$ ^test/Rivet.yoda$ ^test/testBoost$ ^test/.*\.trs$ ^test/testMatVec$ ^test/testMath$ ^test/testBeams$ ^dev$ ^devval$ ^newanalyses$ ^.*plots$ ^a\.out$ ^a\.out\.dSYM$ ^Rivet-.\..\..$ ^local/.*$ ^(for|analyses)\d\d\d$ ^analyses/data$ + +# added by aHg on Tue Mar 5 14:30:12 2019 +syntax: glob +doc/rivet-manual.dvi + +# added by aHg on Tue Mar 5 14:30:21 2019 +syntax: glob +doc/rivet-manual.ps diff --git a/analyses/pluginCMS/CMS_2013_I1224539_ZJET.cc b/analyses/pluginCMS/CMS_2013_I1224539_ZJET.cc --- a/analyses/pluginCMS/CMS_2013_I1224539_ZJET.cc +++ b/analyses/pluginCMS/CMS_2013_I1224539_ZJET.cc @@ -1,195 +1,201 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/FastJets.hh" #include "Rivet/Projections/WFinder.hh" #include "Rivet/Projections/ZFinder.hh" #include "fastjet/tools/Filter.hh" #include "fastjet/tools/Pruner.hh" namespace Rivet { class CMS_2013_I1224539_ZJET : public Analysis { public: /// @name Constructors etc. //@{ /// Constructor CMS_2013_I1224539_ZJET() : Analysis("CMS_2013_I1224539_ZJET"), _filter(fastjet::Filter(fastjet::JetDefinition(fastjet::cambridge_algorithm, 0.3), fastjet::SelectorNHardest(3))), _trimmer(fastjet::Filter(fastjet::JetDefinition(fastjet::kt_algorithm, 0.2), fastjet::SelectorPtFractionMin(0.03))), _pruner(fastjet::Pruner(fastjet::cambridge_algorithm, 0.1, 0.5)) { } //@} public: /// @name Analysis methods //@{ /// Book histograms and initialise projections before the run void init() { FinalState fs(Cuts::abseta < 2.4); declare(fs, "FS"); // Find Zs with pT > 120 GeV ZFinder zfinder(fs, Cuts::abseta < 2.4 && Cuts::pT > 30*GeV, PID::ELECTRON, 80*GeV, 100*GeV, 0.2, ZFinder::CLUSTERNODECAY, ZFinder::TRACK); declare(zfinder, "ZFinder"); // Z+jet jet collections declare(FastJets(zfinder.remainingFinalState(), FastJets::ANTIKT, 0.7), "JetsAK7_zj"); declare(FastJets(zfinder.remainingFinalState(), FastJets::CAM, 0.8), "JetsCA8_zj"); declare(FastJets(zfinder.remainingFinalState(), FastJets::CAM, 1.2), "JetsCA12_zj"); // Histograms /// @note These are 2D histos rendered into slices const int zjetsOffset = 28; for (size_t i = 0; i < N_PT_BINS_vj; ++i ) { _h_ungroomedJetMass_AK7_zj[i] = bookHisto1D(zjetsOffset+i+1+0*N_PT_BINS_vj, 1, 1); _h_filteredJetMass_AK7_zj[i] = bookHisto1D(zjetsOffset+i+1+1*N_PT_BINS_vj,1,1); _h_trimmedJetMass_AK7_zj[i] = bookHisto1D(zjetsOffset+i+1+2*N_PT_BINS_vj,1,1); _h_prunedJetMass_AK7_zj[i] = bookHisto1D(zjetsOffset+i+1+3*N_PT_BINS_vj,1,1); _h_prunedJetMass_CA8_zj[i] = bookHisto1D(zjetsOffset+i+1+4*N_PT_BINS_vj,1,1); if (i > 0) _h_filteredJetMass_CA12_zj[i] = bookHisto1D(zjetsOffset+i+5*N_PT_BINS_vj,1,1); } } bool isBackToBack_zj(const ZFinder& zf, const fastjet::PseudoJet& psjet) { const FourMomentum& z = zf.bosons()[0].momentum(); const FourMomentum& l1 = zf.constituents()[0].momentum(); const FourMomentum& l2 = zf.constituents()[1].momentum(); /// @todo We should make FourMomentum know how to construct itself from a PseudoJet const FourMomentum jmom(psjet.e(), psjet.px(), psjet.py(), psjet.pz()); return (deltaPhi(z, jmom) > 2.0 && deltaR(l1, jmom) > 1.0 && deltaR(l2, jmom) > 1.0); } // Find the pT histogram bin index for value pt (in GeV), to hack a 2D histogram equivalent /// @todo Use a YODA axis/finder alg when available size_t findPtBin(double ptJ) { const double ptBins_vj[N_PT_BINS_vj+1] = { 125.0, 150.0, 220.0, 300.0, 450.0 }; for (size_t ibin = 0; ibin < N_PT_BINS_vj; ++ibin) { if (inRange(ptJ, ptBins_vj[ibin], ptBins_vj[ibin+1])) return ibin; } return N_PT_BINS_vj; } /// Perform the per-event analysis void analyze(const Event& event) { const double weight = event.weight(); // Get the Z const ZFinder& zfinder = apply(event, "ZFinder"); if (zfinder.bosons().size() != 1) vetoEvent; const Particle& z = zfinder.bosons()[0]; - const Particle& l1 = zfinder.constituents()[0]; - const Particle& l2 = zfinder.constituents()[1]; + if ( zfinder.constituents().size() < 2 ) { + cerr << "ZFinder problems: #constinuents: " + << zfinder.constituents().size() << endl + << " #Z-constituents: " + << z.constituents().size() << endl; + } + const Particle l1 = zfinder.constituents()[0]; + const Particle l2 = zfinder.constituents()[1]; // Require a high-pT Z (and constituents) if (l1.pT() < 30*GeV || l2.pT() < 30*GeV || z.pT() < 120*GeV) vetoEvent; // AK7 jets const PseudoJets& psjetsAK7_zj = apply(event, "JetsAK7_zj").pseudoJetsByPt(50.0*GeV); if (!psjetsAK7_zj.empty()) { // Get the leading jet and make sure it's back-to-back with the Z const fastjet::PseudoJet& j0 = psjetsAK7_zj[0]; if (isBackToBack_zj(zfinder, j0)) { const size_t njetBin = findPtBin(j0.pt()/GeV); if (njetBin < N_PT_BINS_vj) { fastjet::PseudoJet filtered0 = _filter(j0); fastjet::PseudoJet trimmed0 = _trimmer(j0); fastjet::PseudoJet pruned0 = _pruner(j0); _h_ungroomedJetMass_AK7_zj[njetBin]->fill(j0.m()/GeV, weight); _h_filteredJetMass_AK7_zj[njetBin]->fill(filtered0.m()/GeV, weight); _h_trimmedJetMass_AK7_zj[njetBin]->fill(trimmed0.m()/GeV, weight); _h_prunedJetMass_AK7_zj[njetBin]->fill(pruned0.m()/GeV, weight); } } } // CA8 jets const PseudoJets& psjetsCA8_zj = apply(event, "JetsCA8_zj").pseudoJetsByPt(50.0*GeV); if (!psjetsCA8_zj.empty()) { // Get the leading jet and make sure it's back-to-back with the Z const fastjet::PseudoJet& j0 = psjetsCA8_zj[0]; if (isBackToBack_zj(zfinder, j0)) { const size_t njetBin = findPtBin(j0.pt()/GeV); if (njetBin < N_PT_BINS_vj) { fastjet::PseudoJet pruned0 = _pruner(j0); _h_prunedJetMass_CA8_zj[njetBin]->fill(pruned0.m()/GeV, weight); } } } // CA12 jets const PseudoJets& psjetsCA12_zj = apply(event, "JetsCA12_zj").pseudoJetsByPt(50.0*GeV); if (!psjetsCA12_zj.empty()) { // Get the leading jet and make sure it's back-to-back with the Z const fastjet::PseudoJet& j0 = psjetsCA12_zj[0]; if (isBackToBack_zj(zfinder, j0)) { const size_t njetBin = findPtBin(j0.pt()/GeV); if (njetBin>0 && njetBin < N_PT_BINS_vj) { fastjet::PseudoJet filtered0 = _filter(j0); _h_filteredJetMass_CA12_zj[njetBin]->fill( filtered0.m() / GeV, weight); } } } } /// Normalise histograms etc., after the run void finalize() { const double normalizationVal = 1000; for (size_t i = 0; i < N_PT_BINS_vj; ++i ) { normalize( _h_ungroomedJetMass_AK7_zj[i], normalizationVal); normalize( _h_filteredJetMass_AK7_zj[i], normalizationVal); normalize( _h_trimmedJetMass_AK7_zj[i], normalizationVal); normalize( _h_prunedJetMass_AK7_zj[i], normalizationVal); normalize( _h_prunedJetMass_CA8_zj[i], normalizationVal); if (i > 0) normalize( _h_filteredJetMass_CA12_zj[i], normalizationVal); } } //@} private: /// @name FastJet grooming tools (configured in constructor init list) //@{ const fastjet::Filter _filter; const fastjet::Filter _trimmer; const fastjet::Pruner _pruner; //@} /// @name Histograms //@{ enum BINS_vj { PT_125_150_vj=0, PT_150_220_vj, PT_220_300_vj, PT_300_450_vj, N_PT_BINS_vj }; Histo1DPtr _h_ungroomedJetMass_AK7_zj[N_PT_BINS_vj]; Histo1DPtr _h_filteredJetMass_AK7_zj[N_PT_BINS_vj]; Histo1DPtr _h_trimmedJetMass_AK7_zj[N_PT_BINS_vj]; Histo1DPtr _h_prunedJetMass_AK7_zj[N_PT_BINS_vj]; Histo1DPtr _h_prunedJetMass_CA8_zj[N_PT_BINS_vj]; Histo1DPtr _h_filteredJetMass_CA12_zj[N_PT_BINS_vj]; //@} }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(CMS_2013_I1224539_ZJET); } diff --git a/include/Rivet/Event.hh b/include/Rivet/Event.hh --- a/include/Rivet/Event.hh +++ b/include/Rivet/Event.hh @@ -1,183 +1,186 @@ // -*- 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) { 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 //@{ /// The generated event obtained from an external event generator const GenEvent* genEvent() const { return &_genevent; } + /// The generated event obtained from an external event generator + const GenEvent* originalGenEvent() const { return _genevent_original; } + /// @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; /// 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 generator centrality (impact-parameter quantile in [0,1]; or -1 if undefined (usual for non-HI generators)) //double centrality() const; // /// Get the boost to the beam centre-of-mass // Vector3 beamCMSBoost() const; // /// Get the boost to the beam centre-of-mass // LorentzTransform beamCMSTransform(); //@} /// @name Access to event particles //@{ /// 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); } //@} /// @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/Projections/Beam.cc b/src/Projections/Beam.cc --- a/src/Projections/Beam.cc +++ b/src/Projections/Beam.cc @@ -1,147 +1,147 @@ // -*- C++ -*- #include "Rivet/Config/RivetCommon.hh" #include "Rivet/Tools/Logging.hh" #include "Rivet/Projections/Beam.hh" namespace Rivet { ParticlePair beams(const Event& e) { // First try the official way: ask the GenEvent for the beam pointers assert(HepMCUtils::particles_size(e.genEvent()) >= 2); - std::pair thebeams = HepMCUtils::beams(e.genEvent()); + std::pair thebeams = HepMCUtils::beams(e.originalGenEvent()); if ( thebeams.first && thebeams.second ) { return ParticlePair{thebeams.first, thebeams.second}; } // Ok, that failed: let's find the status = 4 particles by hand const vector pstat4s = e.allParticles([](const Particle& p){ return p.genParticle()->status() == 4; }); if (pstat4s.size() >= 2) { return ParticlePair{pstat4s[0], pstat4s[1]}; } /// There are no barcodes in HepMC3 /// @todo implement some other fallback rubric? #ifndef ENABLE_HEPMC_3 // Hmm, this sucks. Last guess is that barcodes 1 and 2 are the beams 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)}; } #endif // Give up: return null beams return ParticlePair{Particle(), Particle()}; } 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 { RivetHepMC::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/Tools/RivetHepMC_3.cc b/src/Tools/RivetHepMC_3.cc --- a/src/Tools/RivetHepMC_3.cc +++ b/src/Tools/RivetHepMC_3.cc @@ -1,108 +1,111 @@ // -*- C++ -*- #include "Rivet/Tools/RivetHepMC.hh" #include "HepMC3/ReaderAscii.h" #include "HepMC3/ReaderAsciiHepMC2.h" namespace Rivet{ namespace HepMCUtils{ std::vector particles(ConstGenEventPtr ge){ return ge->particles(); } std::vector particles(const GenEvent *ge){ assert(ge); return ge->particles(); } std::vector vertices(ConstGenEventPtr ge){ return ge->vertices(); } std::vector vertices(const GenEvent *ge){ assert(ge); return ge->vertices(); } std::vector particles(ConstGenVertexPtr gv, const Relatives &relo){ return relo(gv); } std::vector particles(ConstGenParticlePtr gp, const Relatives &relo){ return relo(gp); } int particles_size(ConstGenEventPtr ge){ return particles(ge).size(); } int particles_size(const GenEvent *ge){ return particles(ge).size(); } int uniqueId(ConstGenParticlePtr gp){ return gp->id(); } std::pair beams(const GenEvent *ge) { std::vector beamlist = ge->beams(); - if ( beamlist.size() < 2 ) return std::pair(); + if ( beamlist.size() < 2 ) { + std::cerr << "CANNOT FIND ANY BEAMS!" << std::endl; + return std::pair(); + } return std::make_pair(beamlist[0], beamlist[1]); } bool readEvent(std::shared_ptr io, std::shared_ptr evt){ io->read_event(*evt); return !io->failed(); } shared_ptr makeReader(std::istream & istr, std::string * errm) { shared_ptr ret; // First scan forward and check if there is some hint as to what // kind of file we are looking att. int nchar = 256; std::string header; while ( nchar-- && !istr.eof() ) header += char(istr.get()); // If this stream was too short to contain any reasonable number // of events, just give up. if ( !istr ) { if ( errm ) *errm = "Could not find HepMC header information"; return shared_ptr(); } // Now reset the stream to its original state ... for ( int i = header.length() - 1; i >= 0; --i ) istr.putback(header[i]); // ... and check which kind of format it was and create the // corresponding reader. First try the HepM3 ascii format. if ( header.find("HepMC::Asciiv3-START_EVENT_LISTING") != std::string::npos ) ret = make_shared(istr); // Check if the file is written by WriterRoot or WriterRootTree. else if ( header.substr(0, 4) == "root" ) { if ( errm ) *errm = "Rivet cancurrently not read HepMC root files."; return ret; } // The default is to assume it is a good old HepMC2 ascii file. else { ret = make_shared(istr); } // Check that everything was ok. if ( ret->failed() ) { if ( errm ) *errm = "Problems reading from HepMC file."; ret = shared_ptr(); } return ret; } } }