/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_7/src/GeneratorInterface/LHEInterface/src/JetInput.cc

Go to the documentation of this file.

#include <algorithm>
#include <vector>

#include <Math/GenVector/PxPyPzE4D.h>

#include <HepMC/GenEvent.h>
#include <HepMC/GenParticle.h>
#include <HepMC/GenVertex.h>

#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/Utilities/interface/Exception.h"

#include "GeneratorInterface/LHEInterface/interface/JetInput.h"

namespace lhef {

JetInput::JetInput() :
        partonicFinalState(false),
        onlyHardProcess(false),
        excludeResonances(false),
        tausAsJets(false),
        ptMin(0.0)
{
}

JetInput::JetInput(const edm::ParameterSet &params) :
        partonicFinalState(params.getParameter<bool>("partonicFinalState")),
        onlyHardProcess(params.getParameter<bool>("onlyHardProcess")),
        excludeResonances(false),
        tausAsJets(params.getParameter<bool>("tausAsJets")),
        ptMin(0.0)
{
        if (params.exists("ignoreParticleIDs"))
                setIgnoredParticles(
                        params.getParameter<std::vector<unsigned int> >(
                                                        "ignoreParticleIDs"));
        if (params.exists("excludedResonances"))
                setExcludedResonances(
                        params.getParameter<std::vector<unsigned int> >(
                                                        "excludedResonances"));
        if (params.exists("excludedFromResonances"))
                setExcludedFromResonances(
                        params.getParameter<std::vector<unsigned int> >(
                                                "excludedFromResonances"));
}

JetInput::~JetInput()
{
}

void JetInput::setIgnoredParticles(
                        const std::vector<unsigned int> &particleIDs)
{
        ignoreParticleIDs = particleIDs;
        std::sort(ignoreParticleIDs.begin(), ignoreParticleIDs.end());
}

void JetInput::setExcludedResonances(
                        const std::vector<unsigned int> &particleIDs)
{
        setExcludeResonances(true);
        excludedResonances = particleIDs;
        std::sort(excludedResonances.begin(), excludedResonances.end());
}

void JetInput::setExcludedFromResonances(
                        const std::vector<unsigned int> &particleIDs)
{
        excludedFromResonances = particleIDs;
        std::sort(excludedFromResonances.begin(), excludedFromResonances.end());
}

bool JetInput::isParton(int pdgId) const
{
        pdgId = (pdgId > 0 ? pdgId : -pdgId) % 10000;
        return (pdgId > 0 && pdgId < 6) || pdgId == 7 ||
               pdgId == 9 || (tausAsJets && pdgId == 15) || pdgId == 21;
        // tops are not considered "regular" partons
        // but taus eventually are (since they may hadronize later)
}

bool JetInput::isHadron(int pdgId)
{
        pdgId = (pdgId > 0 ? pdgId : -pdgId) % 10000;
        return (pdgId > 100 && pdgId < 900) ||
               (pdgId > 1000 && pdgId < 9000);
}

static inline bool isContained(const std::vector<unsigned int> &list, int id)
{
        unsigned int absId = (unsigned int)(id > 0 ? id : -id);
        std::vector<unsigned int>::const_iterator pos =
                        std::lower_bound(list.begin(), list.end(), absId);
        return pos != list.end() && *pos == absId;
}

bool JetInput::isResonance(int pdgId) const
{
        if (excludedResonances.empty()) {
                // gauge bosons and tops
                pdgId = (pdgId > 0 ? pdgId : -pdgId) % 10000;
                return (pdgId > 21 && pdgId <= 39) || pdgId == 6 || pdgId == 8;
        }

        return isContained(excludedResonances, pdgId);
}

bool JetInput::isIgnored(int pdgId) const
{
        return isContained(ignoreParticleIDs, pdgId);
}

bool JetInput::isExcludedFromResonances(int pdgId) const
{
        if (excludedFromResonances.empty())
                return true;

        return isContained(excludedFromResonances, pdgId);
}

bool JetInput::isHardProcess(const HepMC::GenVertex *vertex,
                             const VertexVector &hardProcess) const
{
        std::vector<const HepMC::GenVertex*>::const_iterator pos =
                        std::lower_bound(hardProcess.begin(),
                                         hardProcess.end(), vertex);
        return pos != hardProcess.end() && *pos == vertex;
}

static unsigned int partIdx(const JetInput::ParticleVector &p,
                            const HepMC::GenParticle *particle)
{
        JetInput::ParticleVector::const_iterator pos =
                        std::lower_bound(p.begin(), p.end(), particle);
        if (pos == p.end() || *pos != particle)
                throw cms::Exception("CorruptedData")
                        << "HepMC::GenEvent corrupted: Unlisted particles"
                           " in decay tree." << std::endl;

        return pos - p.begin();
}

static void invalidateTree(JetInput::ParticleBitmap &invalid,
                           const JetInput::ParticleVector &p,
                           const HepMC::GenVertex *v)
{
        if (!v)
                return;

        for(HepMC::GenVertex::particles_out_const_iterator iter =
                                        v->particles_out_const_begin();
            iter != v->particles_out_const_end(); ++iter) {
                unsigned int idx = partIdx(p, *iter);

                if (invalid[idx])
                        continue;

                invalid[idx] = true;

                const HepMC::GenVertex *v = (*iter)->end_vertex();
                invalidateTree(invalid, p, v);
        }
}

int JetInput::testPartonChildren(JetInput::ParticleBitmap &invalid,
                                 const JetInput::ParticleVector &p,
                                 const HepMC::GenVertex *v) const
{
        if (!v)
                return 0;

        for(HepMC::GenVertex::particles_out_const_iterator iter =
                                        v->particles_out_const_begin();
            iter != v->particles_out_const_end(); ++iter) {
                unsigned int idx = partIdx(p, *iter);

                if (invalid[idx])
                        continue;

                if (isParton((*iter)->pdg_id()))
                        return 1;
                if (isHadron((*iter)->pdg_id()))
                        return -1;

                const HepMC::GenVertex *v = (*iter)->end_vertex();
                int result = testPartonChildren(invalid, p, v);
                if (result)
                        return result;
        }

        return 0;
}

bool JetInput::hasPartonChildren(ParticleBitmap &invalid,
                                 const ParticleVector &p,
                                 const HepMC::GenParticle *particle) const
{
        return testPartonChildren(invalid, p, particle->end_vertex()) > 0;
}

bool JetInput::fromHardProcess(ParticleBitmap &invalid,
                               const ParticleVector &p,
                               const HepMC::GenParticle *particle,
                               const VertexVector &hardProcess) const
{
        unsigned int idx = partIdx(p, particle);

        if (invalid[idx])
                return false;

        const HepMC::GenVertex *v = particle->production_vertex();
        if (!v)
                return false;
        else if (isHardProcess(v, hardProcess))
                return true;

        for(HepMC::GenVertex::particles_in_const_iterator iter =
                                        v->particles_in_const_begin();
            iter != v->particles_in_const_end(); ++iter)
                if (fromHardProcess(invalid, p, *iter, hardProcess))
                        return true;

        return false;
}

bool JetInput::fromSignalVertex(ParticleBitmap &invalid,
                                const ParticleVector &p,
                                const HepMC::GenParticle *particle,
                                const HepMC::GenVertex *signalVertex) const
{
        unsigned int idx = partIdx(p, particle);

        if (invalid[idx])
                return false;

        const HepMC::GenVertex *v = particle->production_vertex();
        if (!v)
                return false;
        else if (v == signalVertex)
                return true;

        for(HepMC::GenVertex::particles_in_const_iterator iter =
                                        v->particles_in_const_begin();
            iter != v->particles_in_const_end(); ++iter)
                if (fromSignalVertex(invalid, p, *iter, signalVertex))
                        return true;

        return false;
}

JetInput::ResonanceState
JetInput::fromResonance(ParticleBitmap &invalid,
                        const ParticleVector &p,
                        const HepMC::GenParticle *particle,
                        const HepMC::GenVertex *signalVertex,
                        const VertexVector &hardProcess) const
{
        unsigned int idx = partIdx(p, particle);
        int id = particle->pdg_id();

        if (invalid[idx])
                return kIndirect;

        if (particle->status() == 3 && isResonance(id))
                return kDirect;

        if (!isIgnored(id) && excludedFromResonances.empty() && isParton(id))
                return kNo;

        const HepMC::GenVertex *v = particle->production_vertex();
        if (!v)
                return kNo;
        else if (v == signalVertex && excludedResonances.empty())
                return kIndirect;
        else if (v == signalVertex || isHardProcess(v, hardProcess))
                return kNo;

        for(HepMC::GenVertex::particles_in_const_iterator iter =
                                        v->particles_in_const_begin();
            iter != v->particles_in_const_end(); ++iter) {
                ResonanceState result =
                        fromResonance(invalid, p, *iter,
                                      signalVertex, hardProcess);
                switch(result) {
                    case kNo:
                        break;
                    case kDirect:
                        if ((*iter)->pdg_id() == id)
                                return kDirect;
                        if (!isExcludedFromResonances(id))
                                break;
                    case kIndirect:
                        return kIndirect;
                }
        }

        return kNo;
}

JetInput::ParticleVector JetInput::operator () (
                                const HepMC::GenEvent *event) const
{
        if (!event->signal_process_vertex())
                throw cms::Exception("InvalidHepMCEvent")
                        << "HepMC event is lacking signal vertex."
                        << std::endl;

        VertexVector hardProcess, toLookAt;
        std::pair<HepMC::GenParticle*,HepMC::GenParticle*> beamParticles =
                                                event->beam_particles();
        toLookAt.push_back(event->signal_process_vertex());
        while(!toLookAt.empty()) {
                std::vector<const HepMC::GenVertex*> next;
                for(std::vector<const HepMC::GenVertex*>::const_iterator v =
                        toLookAt.begin(); v != toLookAt.end(); ++v) {
                        if (!*v || isHardProcess(*v, hardProcess))
                                continue;

                        bool veto = false;
                        for(HepMC::GenVertex::particles_in_const_iterator iter =
                                        (*v)->particles_in_const_begin();
                            iter != (*v)->particles_in_const_end(); ++iter) {
                                if (*iter == beamParticles.first ||
                                    *iter == beamParticles.second) {
                                        veto = true;
                                        break;
                                }
                        }
                        if (veto)
                                continue;

                        hardProcess.push_back(*v);
                        std::sort(hardProcess.begin(), hardProcess.end());

                        for(HepMC::GenVertex::particles_in_const_iterator iter =
                                        (*v)->particles_in_const_begin();
                            iter != (*v)->particles_in_const_end(); ++iter) {
                                const HepMC::GenVertex *pv =
                                                (*iter)->production_vertex();
                                if (pv)
                                        next.push_back(pv);
                        }
                }

                toLookAt = next;
                std::sort(toLookAt.begin(), toLookAt.end());
        }

        ParticleVector particles;
        for(HepMC::GenEvent::particle_const_iterator iter = event->particles_begin();
            iter != event->particles_end(); ++iter)
                particles.push_back(*iter);

        std::sort(particles.begin(), particles.end());
        unsigned int size = particles.size();

        ParticleBitmap selected(size, false);
        ParticleBitmap invalid(size, false);

        for(unsigned int i = 0; i < size; i++) {
                const HepMC::GenParticle *particle = particles[i];
                if (invalid[i])
                        continue;
                if (particle->status() == 1)
                        selected[i] = true;

                if (partonicFinalState && isParton(particle->pdg_id())) {
                        if (!particle->end_vertex() &&
                            particle->status() != 1)
                                // some brokenness in event...
                                invalid[i] = true;
                        else if (!hasPartonChildren(invalid, particles,
                                                    particle)) {
                                selected[i] = true;
                                invalidateTree(invalid, particles,
                                               particle->end_vertex());
                        }
                }

                if (onlyHardProcess &&
                    !fromHardProcess(invalid, particles,
                                     particle, hardProcess) &&
                    !isHardProcess(particle->end_vertex(), hardProcess))
                        invalid[i] = true;
        }

        ParticleVector result;
        for(unsigned int i = 0; i < size; i++) {
                const HepMC::GenParticle *particle = particles[i];
                if (!selected[i] || invalid[i])
                        continue;

                if (excludeResonances &&
                    fromResonance(invalid, particles, particle,
                                  event->signal_process_vertex(),
                                  hardProcess)) {
                        invalid[i] = true;
                        continue;
                }

                if (isIgnored(particle->pdg_id()))
                        continue;

                if (particle->momentum().perp() >= ptMin)
                        result.push_back(particle);
        }

        return result;
}

} // namespace lhef