Reference Manual

Go to the documentation of this file.
00001 #include <algorithm>
00002 #include <vector>
00003 
00004 #include <Math/GenVector/PxPyPzE4D.h>
00005 
00006 #include <HepMC/GenEvent.h>
00007 #include <HepMC/GenParticle.h>
00008 #include <HepMC/GenVertex.h>
00009 
00010 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00011 #include "FWCore/Utilities/interface/Exception.h"
00012 
00013 #include "GeneratorInterface/LHEInterface/interface/JetInput.h"
00014 
00015 namespace lhef {
00016 
00017 JetInput::JetInput() :
00018         partonicFinalState(false),
00019         onlyHardProcess(false),
00020         excludeResonances(false),
00021         tausAsJets(false),
00022         ptMin(0.0)
00023 {
00024 }
00025 
00026 JetInput::JetInput(const edm::ParameterSet &params) :
00027         partonicFinalState(params.getParameter<bool>("partonicFinalState")),
00028         onlyHardProcess(params.getParameter<bool>("onlyHardProcess")),
00029         excludeResonances(false),
00030         tausAsJets(params.getParameter<bool>("tausAsJets")),
00031         ptMin(0.0)
00032 {
00033         if (params.exists("ignoreParticleIDs"))
00034                 setIgnoredParticles(
00035                         params.getParameter<std::vector<unsigned int> >(
00036                                                         "ignoreParticleIDs"));
00037         if (params.exists("excludedResonances"))
00038                 setExcludedResonances(
00039                         params.getParameter<std::vector<unsigned int> >(
00040                                                         "excludedResonances"));
00041         if (params.exists("excludedFromResonances"))
00042                 setExcludedFromResonances(
00043                         params.getParameter<std::vector<unsigned int> >(
00044                                                 "excludedFromResonances"));
00045 }
00046 
00047 JetInput::~JetInput()
00048 {
00049 }
00050 
00051 void JetInput::setIgnoredParticles(
00052                         const std::vector<unsigned int> &particleIDs)
00053 {
00054         ignoreParticleIDs = particleIDs;
00055         std::sort(ignoreParticleIDs.begin(), ignoreParticleIDs.end());
00056 }
00057 
00058 void JetInput::setExcludedResonances(
00059                         const std::vector<unsigned int> &particleIDs)
00060 {
00061         setExcludeResonances(true);
00062         excludedResonances = particleIDs;
00063         std::sort(excludedResonances.begin(), excludedResonances.end());
00064 }
00065 
00066 void JetInput::setExcludedFromResonances(
00067                         const std::vector<unsigned int> &particleIDs)
00068 {
00069         excludedFromResonances = particleIDs;
00070         std::sort(excludedFromResonances.begin(), excludedFromResonances.end());
00071 }
00072 
00073 bool JetInput::isParton(int pdgId) const
00074 {
00075         pdgId = (pdgId > 0 ? pdgId : -pdgId) % 10000;
00076         return (pdgId > 0 && pdgId < 6) || pdgId == 7 ||
00077                pdgId == 9 || (tausAsJets && pdgId == 15) || pdgId == 21;
00078         // tops are not considered "regular" partons
00079         // but taus eventually are (since they may hadronize later)
00080 }
00081 
00082 bool JetInput::isHadron(int pdgId)
00083 {
00084         pdgId = (pdgId > 0 ? pdgId : -pdgId) % 10000;
00085         return (pdgId > 100 && pdgId < 900) ||
00086                (pdgId > 1000 && pdgId < 9000);
00087 }
00088 
00089 static inline bool isContained(const std::vector<unsigned int> &list, int id)
00090 {
00091         unsigned int absId = (unsigned int)(id > 0 ? id : -id);
00092         std::vector<unsigned int>::const_iterator pos =
00093                         std::lower_bound(list.begin(), list.end(), absId);
00094         return pos != list.end() && *pos == absId;
00095 }
00096 
00097 bool JetInput::isResonance(int pdgId) const
00098 {
00099         if (excludedResonances.empty()) {
00100                 // gauge bosons and tops
00101                 pdgId = (pdgId > 0 ? pdgId : -pdgId) % 10000;
00102                 return (pdgId > 21 && pdgId <= 39) || pdgId == 6 || pdgId == 8;
00103         }
00104 
00105         return isContained(excludedResonances, pdgId);
00106 }
00107 
00108 bool JetInput::isIgnored(int pdgId) const
00109 {
00110         return isContained(ignoreParticleIDs, pdgId);
00111 }
00112 
00113 bool JetInput::isExcludedFromResonances(int pdgId) const
00114 {
00115         if (excludedFromResonances.empty())
00116                 return true;
00117 
00118         return isContained(excludedFromResonances, pdgId);
00119 }
00120 
00121 bool JetInput::isHardProcess(const HepMC::GenVertex *vertex,
00122                              const VertexVector &hardProcess) const
00123 {
00124         std::vector<const HepMC::GenVertex*>::const_iterator pos =
00125                         std::lower_bound(hardProcess.begin(),
00126                                          hardProcess.end(), vertex);
00127         return pos != hardProcess.end() && *pos == vertex;
00128 }
00129 
00130 static unsigned int partIdx(const JetInput::ParticleVector &p,
00131                             const HepMC::GenParticle *particle)
00132 {
00133         JetInput::ParticleVector::const_iterator pos =
00134                         std::lower_bound(p.begin(), p.end(), particle);
00135         if (pos == p.end() || *pos != particle)
00136                 throw cms::Exception("CorruptedData")
00137                         << "HepMC::GenEvent corrupted: Unlisted particles"
00138                            " in decay tree." << std::endl;
00139 
00140         return pos - p.begin();
00141 }
00142 
00143 static void invalidateTree(JetInput::ParticleBitmap &invalid,
00144                            const JetInput::ParticleVector &p,
00145                            const HepMC::GenVertex *v)
00146 {
00147         if (!v)
00148                 return;
00149 
00150         for(HepMC::GenVertex::particles_out_const_iterator iter =
00151                                         v->particles_out_const_begin();
00152             iter != v->particles_out_const_end(); ++iter) {
00153                 unsigned int idx = partIdx(p, *iter);
00154 
00155                 if (invalid[idx])
00156                         continue;
00157 
00158                 invalid[idx] = true;
00159 
00160                 const HepMC::GenVertex *v = (*iter)->end_vertex();
00161                 invalidateTree(invalid, p, v);
00162         }
00163 }
00164 
00165 int JetInput::testPartonChildren(JetInput::ParticleBitmap &invalid,
00166                                  const JetInput::ParticleVector &p,
00167                                  const HepMC::GenVertex *v) const
00168 {
00169         if (!v)
00170                 return 0;
00171 
00172         for(HepMC::GenVertex::particles_out_const_iterator iter =
00173                                         v->particles_out_const_begin();
00174             iter != v->particles_out_const_end(); ++iter) {
00175                 unsigned int idx = partIdx(p, *iter);
00176 
00177                 if (invalid[idx])
00178                         continue;
00179 
00180                 if (isParton((*iter)->pdg_id()))
00181                         return 1;
00182                 if (isHadron((*iter)->pdg_id()))
00183                         return -1;
00184 
00185                 const HepMC::GenVertex *v = (*iter)->end_vertex();
00186                 int result = testPartonChildren(invalid, p, v);
00187                 if (result)
00188                         return result;
00189         }
00190 
00191         return 0;
00192 }
00193 
00194 bool JetInput::hasPartonChildren(ParticleBitmap &invalid,
00195                                  const ParticleVector &p,
00196                                  const HepMC::GenParticle *particle) const
00197 {
00198         return testPartonChildren(invalid, p, particle->end_vertex()) > 0;
00199 }
00200 
00201 bool JetInput::fromHardProcess(ParticleBitmap &invalid,
00202                                const ParticleVector &p,
00203                                const HepMC::GenParticle *particle,
00204                                const VertexVector &hardProcess) const
00205 {
00206         unsigned int idx = partIdx(p, particle);
00207 
00208         if (invalid[idx])
00209                 return false;
00210 
00211         const HepMC::GenVertex *v = particle->production_vertex();
00212         if (!v)
00213                 return false;
00214         else if (isHardProcess(v, hardProcess))
00215                 return true;
00216 
00217         for(HepMC::GenVertex::particles_in_const_iterator iter =
00218                                         v->particles_in_const_begin();
00219             iter != v->particles_in_const_end(); ++iter)
00220                 if (fromHardProcess(invalid, p, *iter, hardProcess))
00221                         return true;
00222 
00223         return false;
00224 }
00225 
00226 bool JetInput::fromSignalVertex(ParticleBitmap &invalid,
00227                                 const ParticleVector &p,
00228                                 const HepMC::GenParticle *particle,
00229                                 const HepMC::GenVertex *signalVertex) const
00230 {
00231         unsigned int idx = partIdx(p, particle);
00232 
00233         if (invalid[idx])
00234                 return false;
00235 
00236         const HepMC::GenVertex *v = particle->production_vertex();
00237         if (!v)
00238                 return false;
00239         else if (v == signalVertex)
00240                 return true;
00241 
00242         for(HepMC::GenVertex::particles_in_const_iterator iter =
00243                                         v->particles_in_const_begin();
00244             iter != v->particles_in_const_end(); ++iter)
00245                 if (fromSignalVertex(invalid, p, *iter, signalVertex))
00246                         return true;
00247 
00248         return false;
00249 }
00250 
00251 JetInput::ResonanceState
00252 JetInput::fromResonance(ParticleBitmap &invalid,
00253                         const ParticleVector &p,
00254                         const HepMC::GenParticle *particle,
00255                         const HepMC::GenVertex *signalVertex,
00256                         const VertexVector &hardProcess) const
00257 {
00258         unsigned int idx = partIdx(p, particle);
00259         int id = particle->pdg_id();
00260 
00261         if (invalid[idx])
00262                 return kIndirect;
00263 
00264         if (particle->status() == 3 && isResonance(id))
00265                 return kDirect;
00266 
00267         if (!isIgnored(id) && excludedFromResonances.empty() && isParton(id))
00268                 return kNo;
00269 
00270         const HepMC::GenVertex *v = particle->production_vertex();
00271         if (!v)
00272                 return kNo;
00273         else if (v == signalVertex && excludedResonances.empty())
00274                 return kIndirect;
00275         else if (v == signalVertex || isHardProcess(v, hardProcess))
00276                 return kNo;
00277 
00278         for(HepMC::GenVertex::particles_in_const_iterator iter =
00279                                         v->particles_in_const_begin();
00280             iter != v->particles_in_const_end(); ++iter) {
00281                 ResonanceState result =
00282                         fromResonance(invalid, p, *iter,
00283                                       signalVertex, hardProcess);
00284                 switch(result) {
00285                     case kNo:
00286                         break;
00287                     case kDirect:
00288                         if ((*iter)->pdg_id() == id)
00289                                 return kDirect;
00290                         if (!isExcludedFromResonances(id))
00291                                 break;
00292                     case kIndirect:
00293                         return kIndirect;
00294                 }
00295         }
00296 
00297         return kNo;
00298 }
00299 
00300 JetInput::ParticleVector JetInput::operator () (
00301                                 const HepMC::GenEvent *event) const
00302 {
00303         if (!event->signal_process_vertex())
00304                 throw cms::Exception("InvalidHepMCEvent")
00305                         << "HepMC event is lacking signal vertex."
00306                         << std::endl;
00307 
00308         VertexVector hardProcess, toLookAt;
00309         std::pair<HepMC::GenParticle*,HepMC::GenParticle*> beamParticles =
00310                                                 event->beam_particles();
00311         toLookAt.push_back(event->signal_process_vertex());
00312         while(!toLookAt.empty()) {
00313                 std::vector<const HepMC::GenVertex*> next;
00314                 for(std::vector<const HepMC::GenVertex*>::const_iterator v =
00315                         toLookAt.begin(); v != toLookAt.end(); ++v) {
00316                         if (!*v || isHardProcess(*v, hardProcess))
00317                                 continue;
00318 
00319                         bool veto = false;
00320                         for(HepMC::GenVertex::particles_in_const_iterator iter =
00321                                         (*v)->particles_in_const_begin();
00322                             iter != (*v)->particles_in_const_end(); ++iter) {
00323                                 if (*iter == beamParticles.first ||
00324                                     *iter == beamParticles.second) {
00325                                         veto = true;
00326                                         break;
00327                                 }
00328                         }
00329                         if (veto)
00330                                 continue;
00331 
00332                         hardProcess.push_back(*v);
00333                         std::sort(hardProcess.begin(), hardProcess.end());
00334 
00335                         for(HepMC::GenVertex::particles_in_const_iterator iter =
00336                                         (*v)->particles_in_const_begin();
00337                             iter != (*v)->particles_in_const_end(); ++iter) {
00338                                 const HepMC::GenVertex *pv =
00339                                                 (*iter)->production_vertex();
00340                                 if (pv)
00341                                         next.push_back(pv);
00342                         }
00343                 }
00344 
00345                 toLookAt = next;
00346                 std::sort(toLookAt.begin(), toLookAt.end());
00347         }
00348 
00349         ParticleVector particles;
00350         for(HepMC::GenEvent::particle_const_iterator iter = event->particles_begin();
00351             iter != event->particles_end(); ++iter)
00352                 particles.push_back(*iter);
00353 
00354         std::sort(particles.begin(), particles.end());
00355         unsigned int size = particles.size();
00356 
00357         ParticleBitmap selected(size, false);
00358         ParticleBitmap invalid(size, false);
00359 
00360         for(unsigned int i = 0; i < size; i++) {
00361                 const HepMC::GenParticle *particle = particles[i];
00362                 if (invalid[i])
00363                         continue;
00364                 if (particle->status() == 1)
00365                         selected[i] = true;
00366 
00367                 if (partonicFinalState && isParton(particle->pdg_id())) {
00368                         if (!particle->end_vertex() &&
00369                             particle->status() != 1)
00370                                 // some brokenness in event...
00371                                 invalid[i] = true;
00372                         else if (!hasPartonChildren(invalid, particles,
00373                                                     particle)) {
00374                                 selected[i] = true;
00375                                 invalidateTree(invalid, particles,
00376                                                particle->end_vertex());
00377                         }
00378                 }
00379 
00380                 if (onlyHardProcess &&
00381                     !fromHardProcess(invalid, particles,
00382                                      particle, hardProcess) &&
00383                     !isHardProcess(particle->end_vertex(), hardProcess))
00384                         invalid[i] = true;
00385         }
00386 
00387         ParticleVector result;
00388         for(unsigned int i = 0; i < size; i++) {
00389                 const HepMC::GenParticle *particle = particles[i];
00390                 if (!selected[i] || invalid[i])
00391                         continue;
00392 
00393                 if (excludeResonances &&
00394                     fromResonance(invalid, particles, particle,
00395                                   event->signal_process_vertex(),
00396                                   hardProcess)) {
00397                         invalid[i] = true;
00398                         continue;
00399                 }
00400 
00401                 if (isIgnored(particle->pdg_id()))
00402                         continue;
00403 
00404                 if (particle->momentum().perp() >= ptMin)
00405                         result.push_back(particle);
00406         }
00407 
00408         return result;
00409 }
00410 
00411 } // namespace lhef
/data/refman/pasoursint/CMSSW_5_3_0/src/GeneratorInterface/LHEInterface/src/JetInput.cc
