diff --git a/src/Projections/MergedFinalState.cc b/src/Projections/MergedFinalState.cc --- a/src/Projections/MergedFinalState.cc +++ b/src/Projections/MergedFinalState.cc @@ -1,41 +1,40 @@ // -*- C++ -*- #include "Rivet/Projections/MergedFinalState.hh" namespace Rivet { int MergedFinalState::compare(const Projection& p) const { /// @todo: Currently A+B is not recognised to be the same as B+A. return mkNamedPCmp(p, "FSA") || mkNamedPCmp(p, "FSB"); } void MergedFinalState::project(const Event& e) { const FinalState& fsa = applyProjection(e, "FSA"); const FinalState& fsb = applyProjection(e, "FSB"); _theParticles.clear(); - foreach (const Particle& pa, fsa.particles()){ + for (const Particle& pa : fsa.particles()){ _theParticles.push_back(pa); } - foreach (const Particle& pb, fsb.particles()){ - const GenParticle* originalb = pb.genParticle(); + for (const Particle& pb : fsb.particles()){ + const GenParticlePtr originalb = pb.genParticle(); bool notfound = true; - foreach (const Particle& pa, fsa.particles()){ - const GenParticle* originala = pa.genParticle(); + for (const Particle& pa : fsa.particles()){ + const GenParticlePtr originala = pa.genParticle(); if (originala == originalb) { notfound = false; break; } } - if (notfound) { - _theParticles.push_back(pb); - } + if (notfound) _theParticles.push_back(pb); + } MSG_DEBUG("Number of particles in the two final states to be merged: = \n" << " 1st final state = " << fsa.particles().size() << "\n" << " 2nd final state = " << fsb.particles().size()); MSG_DEBUG("Number of merged final-state particles = " << _theParticles.size()); } } diff --git a/src/Projections/PrimaryHadrons.cc b/src/Projections/PrimaryHadrons.cc --- a/src/Projections/PrimaryHadrons.cc +++ b/src/Projections/PrimaryHadrons.cc @@ -1,39 +1,39 @@ // -*- C++ -*- #include "Rivet/Projections/PrimaryHadrons.hh" namespace Rivet { void PrimaryHadrons::project(const Event& e) { _theParticles.clear(); const Particles& unstables = applyProjection(e, "UFS").particles(); - foreach (const Particle& p, unstables) { + for (const Particle& p : unstables) { // Exclude taus etc. if (!isHadron(p)) continue; // A spontaneously appearing hadron: this is weird, but I guess is allowed... and is primary if (!p.genParticle() || !p.genParticle()->production_vertex()) { MSG_DEBUG("Hadron " << p.pid() << " with no GenParticle or parent found: treating as primary"); _theParticles.push_back(p); continue; } // There are ancestors -- check them for status=2 hadronic content - const vector ancestors = particles_in(p.genParticle(), HepMC::ancestors); + const vector ancestors = Rivet::particles(p.genParticle(), HepMC::ancestors); bool has_hadron_parent = false; - foreach (const GenParticle* pa, ancestors) { + for (const GenParticlePtr pa : ancestors) { if (pa->status() != 2) continue; /// @todo Are hadrons from tau decays "primary hadrons"? I guess not if (PID::isHadron(pa->pdg_id()) || abs(pa->pdg_id()) == PID::TAU) { has_hadron_parent = true; break; } } // If the particle seems to be a primary hadron, add it to the list if (!has_hadron_parent) _theParticles.push_back(p); } MSG_DEBUG("Number of primary hadrons = " << _theParticles.size()); } } diff --git a/src/Projections/UnstableFinalState.cc b/src/Projections/UnstableFinalState.cc --- a/src/Projections/UnstableFinalState.cc +++ b/src/Projections/UnstableFinalState.cc @@ -1,77 +1,63 @@ // -*- C++ -*- #include "Rivet/Projections/UnstableFinalState.hh" -/// @todo Replace with PID::isParton() -#define IS_PARTON_PDGID(id) ( abs(id) <= 100 && abs(id) != 22 && (abs(id) < 11 || abs(id) > 18) ) - namespace Rivet { /// @todo Add a FIRST/LAST/ANY enum to specify the mode for uniquifying replica chains (default = LAST) void UnstableFinalState::project(const Event& e) { _theParticles.clear(); /// @todo Replace PID veto list with PID:: functions? vector vetoIds; vetoIds += 22; // status 2 photons don't count! vetoIds += 110; vetoIds += 990; vetoIds += 9990; // Reggeons //vetoIds += 9902210; // something weird from PYTHIA6 - foreach (const GenParticle* p, Rivet::particles(e.genEvent())) { + for (const GenParticlePtr p : Rivet::particles(e.genEvent())) { const int st = p->status(); bool passed = (st == 1 || (st == 2 && !contains(vetoIds, abs(p->pdg_id())))) && - !IS_PARTON_PDGID(p->pdg_id()) && ///< Always veto partons - !p->is_beam() && // Filter beam particles - _cuts->accept(p->momentum()); + !PID::isParton(p->pdg_id()) && ///< Always veto partons + p != e.beams().first.genParticle() && p != e.beams().second.genParticle() && // Filter beam particles + _cuts->accept(p); // Avoid double counting by re-marking as unpassed if ID == (any) parent ID - const GenVertex* pv = p->production_vertex(); - // if (passed && pv) { - // for (GenVertex::particles_in_const_iterator pp = pv->particles_in_const_begin(); pp != pv->particles_in_const_end(); ++pp) { - // if (p->pdg_id() == (*pp)->pdg_id() && (*pp)->status() == 2) { - // passed = false; - // break; - // } - // } - // } - // + const GenVertexPtr pv = p->production_vertex(); // Avoid double counting by re-marking as unpassed if ID == any child ID - const GenVertex* dv = p->end_vertex(); + const GenVertexPtr dv = p->end_vertex(); if (passed && dv) { - foreach (GenParticle* pp, particles_out(dv)) { + for (GenParticlePtr pp : particles_out(dv)) { if (p->pdg_id() == pp->pdg_id() && pp->status() == 2) { passed = false; break; } } } // Add to output particles collection - if (passed) _theParticles.push_back(Particle(*p)); + if (passed) _theParticles.push_back(Particle(p)); // Log parents and children if (getLog().isActive(Log::TRACE)) { MSG_TRACE("ID = " << p->pdg_id() << ", status = " << st << ", pT = " << p->momentum().perp() << ", eta = " << p->momentum().eta() << ": result = " << std::boolalpha << passed); if (pv) { - for (GenVertex::particles_in_const_iterator pp = pv->particles_in_const_begin(); pp != pv->particles_in_const_end(); ++pp) { - MSG_TRACE(" parent ID = " << (*pp)->pdg_id()); - } + for (const GenParticlePtr pp : Rivet::particles(pv, HepMC::parents)) + MSG_TRACE(" parent ID = " << pp->pdg_id()); } if (dv) { - for (GenVertex::particles_out_const_iterator pp = dv->particles_out_const_begin(); pp != dv->particles_out_const_end(); ++pp) { - MSG_TRACE(" child ID = " << (*pp)->pdg_id()); - } + for (const GenParticlePtr pp : Rivet::particles(dv, HepMC::children)) + MSG_TRACE(" child ID = " << pp->pdg_id()); } } } MSG_DEBUG("Number of unstable final-state particles = " << _theParticles.size()); } } diff --git a/src/Projections/VetoedFinalState.cc b/src/Projections/VetoedFinalState.cc --- a/src/Projections/VetoedFinalState.cc +++ b/src/Projections/VetoedFinalState.cc @@ -1,147 +1,128 @@ // -*- C++ -*- #include "Rivet/Projections/VetoedFinalState.hh" namespace Rivet { int VetoedFinalState::compare(const Projection& p) const { const PCmp fscmp = mkNamedPCmp(p, "FS"); if (fscmp != EQUIVALENT) return fscmp; /// @todo We can do better than this... if (_vetofsnames.size() != 0) return UNDEFINED; const VetoedFinalState& other = dynamic_cast(p); return \ cmp(_vetoCodes, other._vetoCodes) || cmp(_compositeVetoes, other._compositeVetoes) || cmp(_parentVetoes, other._parentVetoes); } void VetoedFinalState::project(const Event& e) { const FinalState& fs = applyProjection(e, "FS"); _theParticles.clear(); _theParticles.reserve(fs.particles().size()); - foreach (const Particle& p, fs.particles()) { + + for (const Particle& p : fs.particles()) { if (getLog().isActive(Log::TRACE)) { vector codes; for (VetoDetails::const_iterator code = _vetoCodes.begin(); code != _vetoCodes.end(); ++code) { codes.push_back(code->first); } const string codestr = "{ " + join(codes) + " }"; MSG_TRACE(p.pid() << " vs. veto codes = " << codestr << " (" << codes.size() << ")"); } VetoDetails::iterator iter = _vetoCodes.find(p.pid()); if (iter == _vetoCodes.end()) { MSG_TRACE("Storing with PDG code = " << p.pid() << ", pT = " << p.pT()); _theParticles.push_back(p); } else { // This particle code is listed as a possible veto... check pT. // Make sure that the pT range is sensible: BinaryCut ptrange = iter->second; assert(ptrange.first <= ptrange.second); stringstream rangess; if (ptrange.first < numeric_limits::max()) rangess << ptrange.second; rangess << " - "; if (ptrange.second < numeric_limits::max()) rangess << ptrange.second; MSG_TRACE("ID = " << p.pid() << ", pT range = " << rangess.str()); stringstream debugline; debugline << "with PDG code = " << p.pid() << " pT = " << p.pT(); if (p.pT() < ptrange.first || p.pT() > ptrange.second) { MSG_TRACE("Storing " << debugline.str()); _theParticles.push_back(p); } else { MSG_TRACE("Vetoing " << debugline.str()); } } } set toErase; for (set::iterator nIt = _nCompositeDecays.begin(); nIt != _nCompositeDecays.end() && !_theParticles.empty(); ++nIt) { map, FourMomentum> oldMasses; map, FourMomentum> newMasses; set start; start.insert(_theParticles.begin()); oldMasses.insert(pair, FourMomentum> (start, _theParticles.begin()->momentum())); for (int nParts = 1; nParts != *nIt; ++nParts) { for (map, FourMomentum>::iterator mIt = oldMasses.begin(); mIt != oldMasses.end(); ++mIt) { Particles::iterator pStart = *(mIt->first.rbegin()); for (Particles::iterator pIt = pStart + 1; pIt != _theParticles.end(); ++pIt) { FourMomentum cMom = mIt->second + pIt->momentum(); set pList(mIt->first); pList.insert(pIt); newMasses[pList] = cMom; } } oldMasses = newMasses; newMasses.clear(); } for (map, FourMomentum>::iterator mIt = oldMasses.begin(); mIt != oldMasses.end(); ++mIt) { double mass2 = mIt->second.mass2(); if (mass2 >= 0.0) { double mass = sqrt(mass2); for (CompositeVeto::iterator cIt = _compositeVetoes.lower_bound(*nIt); cIt != _compositeVetoes.upper_bound(*nIt); ++cIt) { BinaryCut massRange = cIt->second; if (mass < massRange.second && mass > massRange.first) { for (set::iterator lIt = mIt->first.begin(); lIt != mIt->first.end(); ++lIt) { toErase.insert(*lIt); } } } } } } for (set::reverse_iterator p = toErase.rbegin(); p != toErase.rend(); ++p) { _theParticles.erase(*p); } // Remove particles whose parents match entries in the parent veto PDG ID codes list - /// @todo There must be a nice way to do this -- an STL algorithm (or we provide a nicer wrapper) - foreach (PdgId vetoid, _parentVetoes) { - for (Particles::iterator ip = _theParticles.begin(); ip != _theParticles.end(); ++ip) { - const GenVertex* startVtx = ip->genParticle()->production_vertex(); - if (startVtx == NULL) continue; - // Loop over parents and test their IDs - /// @todo Could use any() here? - foreach (const GenParticle* parent, Rivet::particles(startVtx, HepMC::ancestors)) { - if (vetoid == parent->pdg_id()) { - ip = _theParticles.erase(ip); --ip; //< Erase this _theParticles entry - break; + for (PdgId vetoid : _parentVetoes) + ifilter_discard(_theParticles, [&](const Particle& p) { + // Loop over parents and test their IDs + for (const GenParticlePtr parent : Rivet::particles(p.genParticle(), HepMC::ancestors)) + if (parent->pdg_id() == vetoid) return true; + return false; }); + + // Now veto on FS particles from other projections + for (const string& ifs : _vetofsnames) { + const Particles& vfsp = applyProjection(e, ifs).particles(); + ifilter_discard(_theParticles, [&](const Particle& p) { + for (const Particle& pcheck : vfsp) { + MSG_TRACE("Comparing with veto particle: " << p << " vs. " << pcheck); + if (p.genParticle() == pcheck.genParticle()) return true; } - } - } + return false; }); } - // Now veto on the FS - foreach (const string& ifs, _vetofsnames) { - const FinalState& vfs = applyProjection(e, ifs); - const Particles& vfsp = vfs.particles(); - for (Particles::iterator icheck = _theParticles.begin(); icheck != _theParticles.end(); ++icheck) { - if (icheck->genParticle() == NULL) continue; - bool found = false; - for (Particles::const_iterator ipart = vfsp.begin(); ipart != vfsp.end(); ++ipart){ - if (ipart->genParticle() == NULL) continue; - MSG_TRACE("Comparing barcode " << icheck->genParticle()->barcode() - << " with veto particle " << ipart->genParticle()->barcode()); - if (ipart->genParticle()->barcode() == icheck->genParticle()->barcode()){ - found = true; - break; - } - } - if (found) { - _theParticles.erase(icheck); - --icheck; - } - } - } } }