HepMCCandAlgos/plugins/ParticleListDrawer.cc

0001 #include <memory>
0002 #include <string>
0003 #include <iostream>
0004 #include <sstream>
0005
0006 #include "FWCore/Framework/interface/Frameworkfwd.h"
0007 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0008 #include "FWCore/Framework/interface/ESHandle.h"
0009 #include "FWCore/Utilities/interface/ESGetToken.h"
0010 #include "FWCore/Framework/interface/Event.h"
0011 #include "FWCore/Framework/interface/EventSetup.h"
0012 #include "FWCore/Framework/interface/MakerMacros.h"
0013 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0014 #include "FWCore/Utilities/interface/InputTag.h"
0015 #include "SimGeneral/HepPDTRecord/interface/ParticleDataTable.h"
0016 #include "DataFormats/Candidate/interface/Candidate.h"
0017 #include "DataFormats/Candidate/interface/CandidateFwd.h"
0018 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0019 #include "DataFormats/Common/interface/Ref.h"
0020
0021 /**
0022    \class   ParticleListDrawer ParticleListDrawer.h "PhysicsTools/HepMCCandAlgos/plugins/ParticleListDrawer.h"
0023
0024    \brief   Module to analyze the particle listing as provided by common event generators
0025
0026    Module to analyze the particle listing as provided by common event generators equivalent
0027    to PYLIST(1) (from pythia). It is expected to run on vectors of  reo::GenParticles. For
0028    an example of use have a look to:
0029
0030    PhysicsTools/HepMCCandAlgos/test/testParticleTreeDrawer.py
0031
0032    Caveats:
0033    Status 3 particles can have daughters both with status 2 and 3. In pythia this is not
0034    the same mother-daughter. The relations are correct but special care has to be taken
0035    when looking at mother-daughter relation which involve status 2 and 3 particles.
0036 */
0037
0038 using namespace std;
0039 using namespace reco;
0040 using namespace edm;
0041
0042 class ParticleListDrawer : public edm::one::EDAnalyzer<> {
0043 public:
0044   explicit ParticleListDrawer(const edm::ParameterSet&);
0045   ~ParticleListDrawer() override{};
0046   void analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) override;
0047
0048 private:
0049   std::string getParticleName(int id) const;
0050
0051   edm::InputTag src_;
0052   edm::EDGetTokenT<reco::CandidateView> srcToken_;
0053   edm::ESHandle<ParticleDataTable> pdt_;
0054   edm::ESGetToken<ParticleDataTable, edm::DefaultRecord> pdtToken_;
0055   int maxEventsToPrint_;  // Must be signed, because -1 is used for no limit
0056   unsigned int nEventAnalyzed_;
0057   bool printOnlyHardInteraction_;
0058   bool printVertex_;
0059   bool printFlags_;
0060   bool useMessageLogger_;
0061 };
0062
0063 ParticleListDrawer::ParticleListDrawer(const edm::ParameterSet& pset)
0064     : src_(pset.getParameter<InputTag>("src")),
0065       srcToken_(consumes<reco::CandidateView>(src_)),
0066       pdtToken_(esConsumes<ParticleDataTable, edm::DefaultRecord>()),
0067       maxEventsToPrint_(pset.getUntrackedParameter<int>("maxEventsToPrint", 1)),
0068       nEventAnalyzed_(0),
0069       printOnlyHardInteraction_(pset.getUntrackedParameter<bool>("printOnlyHardInteraction", false)),
0070       printVertex_(pset.getUntrackedParameter<bool>("printVertex", false)),
0071       printFlags_(pset.getUntrackedParameter<bool>("printFlags", false)),
0072       useMessageLogger_(pset.getUntrackedParameter<bool>("useMessageLogger", false)) {}
0073
0074 std::string ParticleListDrawer::getParticleName(int id) const {
0075   const ParticleData* pd = pdt_->particle(id);
0076   if (!pd) {
0077     std::ostringstream ss;
0078     ss << "P" << id;
0079     return ss.str();
0080   } else
0081     return pd->name();
0082 }
0083
0084 void ParticleListDrawer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0085   Handle<reco::CandidateView> particles;
0086   iEvent.getByToken(srcToken_, particles);
0087   pdt_ = iSetup.getHandle(pdtToken_);
0088
0089   if (maxEventsToPrint_ < 0 || nEventAnalyzed_ < static_cast<unsigned int>(maxEventsToPrint_)) {
0090     ostringstream out;
0091     char buf[256];
0092
0093     out << endl << "[ParticleListDrawer] analysing particle collection " << src_.label() << endl;
0094
0095     snprintf(buf,
0096              sizeof(buf),
0097              " idx  |    ID -       Name |Stat|  Mo1  Mo2  Da1  Da2 |nMo nDa|    pt       eta     phi   |     px       "
0098              "  py         pz        m     |");
0099     out << buf;
0100     if (printVertex_) {
0101       snprintf(buf, sizeof(buf), "        vx       vy        vz     |");
0102       out << buf;
0103     }
0104     out << endl;
0105
0106     int idx = -1;
0107     int iMo1 = -1;
0108     int iMo2 = -1;
0109     int iDa1 = -1;
0110     int iDa2 = -1;
0111     vector<const reco::Candidate*> cands;
0112     vector<const Candidate*>::const_iterator found = cands.begin();
0113     for (CandidateView::const_iterator p = particles->begin(); p != particles->end(); ++p) {
0114       cands.push_back(&*p);
0115     }
0116
0117     for (CandidateView::const_iterator p = particles->begin(); p != particles->end(); p++) {
0118       if (printOnlyHardInteraction_ && p->status() != 3)
0119         continue;
0120
0121       // Particle Name
0122       int id = p->pdgId();
0123       string particleName = getParticleName(id);
0124
0125       // Particle Index
0126       idx = p - particles->begin();
0127
0128       // Particles Mothers and Daighters
0129       iMo1 = -1;
0130       iMo2 = -1;
0131       iDa1 = -1;
0132       iDa2 = -1;
0133       int nMo = p->numberOfMothers();
0134       int nDa = p->numberOfDaughters();
0135
0136       found = find(cands.begin(), cands.end(), p->mother(0));
0137       if (found != cands.end())
0138         iMo1 = found - cands.begin();
0139
0140       found = find(cands.begin(), cands.end(), p->mother(nMo - 1));
0141       if (found != cands.end())
0142         iMo2 = found - cands.begin();
0143
0144       found = find(cands.begin(), cands.end(), p->daughter(0));
0145       if (found != cands.end())
0146         iDa1 = found - cands.begin();
0147
0148       found = find(cands.begin(), cands.end(), p->daughter(nDa - 1));
0149       if (found != cands.end())
0150         iDa2 = found - cands.begin();
0151
0152       char buf[2400];
0153       snprintf(
0154           buf,
0155           sizeof(buf),
0156           " %4d | %5d - %10s | %2d | %4d %4d %4d %4d | %2d %2d | %7.3f %10.3f %6.3f | %10.3f %10.3f %10.3f %8.3f |",
0157           idx,
0158           p->pdgId(),
0159           particleName.c_str(),
0160           p->status(),
0161           iMo1,
0162           iMo2,
0163           iDa1,
0164           iDa2,
0165           nMo,
0166           nDa,
0167           p->pt(),
0168           p->eta(),
0169           p->phi(),
0170           p->px(),
0171           p->py(),
0172           p->pz(),
0173           p->mass());
0174       out << buf;
0175
0176       if (printVertex_) {
0177         snprintf(buf, sizeof(buf), " %10.3f %10.3f %10.3f |", p->vertex().x(), p->vertex().y(), p->vertex().z());
0178         out << buf;
0179       }
0180
0181       if (printFlags_) {
0182         const reco::GenParticle* gp = dynamic_cast<const reco::GenParticle*>(&*p);
0183         if (!gp)
0184           throw cms::Exception("Unsupported", "Status flags can be printed only for reco::GenParticle objects\n");
0185         if (gp->isPromptFinalState())
0186           out << "  PromptFinalState";
0187         if (gp->isDirectPromptTauDecayProductFinalState())
0188           out << "  DirectPromptTauDecayProductFinalState";
0189         if (gp->isHardProcess())
0190           out << "  HardProcess";
0191         if (gp->fromHardProcessFinalState())
0192           out << "  HardProcessFinalState";
0193         if (gp->fromHardProcessBeforeFSR())
0194           out << "  HardProcessBeforeFSR";
0195         if (gp->statusFlags().isFirstCopy())
0196           out << "  FirstCopy";
0197         if (gp->isLastCopy())
0198           out << "  LastCopy";
0199         if (gp->isLastCopyBeforeFSR())
0200           out << "  LastCopyBeforeFSR";
0201       }
0202
0203       out << endl;
0204     }
0205     nEventAnalyzed_++;
0206
0207     if (useMessageLogger_)
0208       LogVerbatim("ParticleListDrawer") << out.str();
0209     else
0210       cout << out.str();
0211   }
0212 }
0213
0214 DEFINE_FWK_MODULE(ParticleListDrawer);