HepMCCandAlgos/plugins/ParticleDecayDrawer.cc

0001 /* class ParticleDecayDrawer
0002  *
0003  * \author Luca Lista, INFN
0004  */
0005 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0006 #include "FWCore/Utilities/interface/InputTag.h"
0007 #include "FWCore/Framework/interface/ESHandle.h"
0008 #include "FWCore/Utilities/interface/ESGetToken.h"
0009 #include "SimGeneral/HepPDTRecord/interface/ParticleDataTable.h"
0010 #include "DataFormats/Common/interface/View.h"
0011 #include "DataFormats/Candidate/interface/Candidate.h"
0012
0013 class ParticleDecayDrawer : public edm::one::EDAnalyzer<> {
0014 public:
0015   ParticleDecayDrawer(const edm::ParameterSet &);
0016
0017 private:
0018   void analyze(const edm::Event &, const edm::EventSetup &) override;
0019   edm::EDGetTokenT<edm::View<reco::Candidate> > srcToken_;
0020   std::string decay(const reco::Candidate &, std::list<const reco::Candidate *> &) const;
0021   edm::ESHandle<ParticleDataTable> pdt_;
0022   edm::ESGetToken<ParticleDataTable, edm::DefaultRecord> pdtToken_;
0023   /// print parameters
0024   bool printP4_, printPtEtaPhi_, printVertex_;
0025   /// print 4 momenta
0026   std::string printP4(const reco::Candidate &) const;
0027   /// accept candidate
0028   bool accept(const reco::Candidate &, const std::list<const reco::Candidate *> &) const;
0029   /// select candidate
0030   bool select(const reco::Candidate &) const;
0031   /// has valid daughters in the chain
0032   bool hasValidDaughters(const reco::Candidate &) const;
0033 };
0034
0035 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0036 #include "FWCore/Framework/interface/Event.h"
0037 #include "DataFormats/Common/interface/Handle.h"
0038 #include "FWCore/Framework/interface/EventSetup.h"
0039 #include "FWCore/Utilities/interface/EDMException.h"
0040 #include <iostream>
0041 #include <sstream>
0042 #include <algorithm>
0043 using namespace std;
0044 using namespace edm;
0045 using namespace reco;
0046
0047 ParticleDecayDrawer::ParticleDecayDrawer(const ParameterSet &cfg)
0048     : srcToken_(consumes<edm::View<reco::Candidate> >(cfg.getParameter<InputTag>("src"))),
0049       pdtToken_(esConsumes<ParticleDataTable, edm::DefaultRecord>()),
0050       printP4_(cfg.getUntrackedParameter<bool>("printP4", false)),
0051       printPtEtaPhi_(cfg.getUntrackedParameter<bool>("printPtEtaPhi", false)),
0052       printVertex_(cfg.getUntrackedParameter<bool>("printVertex", false)) {}
0053
0054 bool ParticleDecayDrawer::accept(const reco::Candidate &c, const list<const Candidate *> &skip) const {
0055   if (find(skip.begin(), skip.end(), &c) != skip.end())
0056     return false;
0057   return select(c);
0058 }
0059
0060 bool ParticleDecayDrawer::select(const reco::Candidate &c) const { return c.status() == 3; }
0061
0062 bool ParticleDecayDrawer::hasValidDaughters(const reco::Candidate &c) const {
0063   size_t ndau = c.numberOfDaughters();
0064   for (size_t i = 0; i < ndau; ++i)
0065     if (select(*c.daughter(i)))
0066       return true;
0067   return false;
0068 }
0069
0070 void ParticleDecayDrawer::analyze(const Event &event, const EventSetup &es) {
0071   pdt_ = es.getHandle(pdtToken_);
0072   Handle<View<Candidate> > particles;
0073   event.getByToken(srcToken_, particles);
0074   list<const Candidate *> skip;
0075   vector<const Candidate *> nodes, moms;
0076   for (View<Candidate>::const_iterator p = particles->begin(); p != particles->end(); ++p) {
0077     if (p->numberOfMothers() > 1) {
0078       if (select(*p)) {
0079         skip.push_back(&*p);
0080         nodes.push_back(&*p);
0081         for (size_t j = 0; j < p->numberOfMothers(); ++j) {
0082           const Candidate *mom = p->mother(j);
0083           const Candidate *grandMom;
0084           while ((grandMom = mom->mother()) != nullptr)
0085             mom = grandMom;
0086           if (select(*mom)) {
0087             moms.push_back(mom);
0088           }
0089         }
0090       }
0091     }
0092   }
0093   cout << "-- decay: --" << endl;
0094   if (!moms.empty()) {
0095     if (moms.size() > 1)
0096       for (size_t m = 0; m < moms.size(); ++m) {
0097         string dec = decay(*moms[m], skip);
0098         if (!dec.empty())
0099           cout << "{ " << dec << " } ";
0100       }
0101     else
0102       cout << decay(*moms[0], skip);
0103   }
0104   if (!nodes.empty()) {
0105     cout << "-> ";
0106     if (nodes.size() > 1) {
0107       for (size_t n = 0; n < nodes.size(); ++n) {
0108         skip.remove(nodes[n]);
0109         string dec = decay(*nodes[n], skip);
0110         if (!dec.empty()) {
0111           if (dec.find("->", 0) != string::npos)
0112             cout << " ( " << dec << " )";
0113           else
0114             cout << " " << dec;
0115         }
0116       }
0117     } else {
0118       skip.remove(nodes[0]);
0119       cout << decay(*nodes[0], skip);
0120     }
0121   }
0122   cout << endl;
0123 }
0124
0125 string ParticleDecayDrawer::printP4(const Candidate &c) const {
0126   ostringstream cout;
0127   if (printP4_)
0128     cout << " (" << c.px() << ", " << c.py() << ", " << c.pz() << "; " << c.energy() << ")";
0129   if (printPtEtaPhi_)
0130     cout << " [" << c.pt() << ": " << c.eta() << ", " << c.phi() << "]";
0131   if (printVertex_)
0132     cout << " {" << c.vx() << ", " << c.vy() << ", " << c.vz() << "}";
0133   return cout.str();
0134 }
0135
0136 string ParticleDecayDrawer::decay(const Candidate &c, list<const Candidate *> &skip) const {
0137   string out;
0138   if (find(skip.begin(), skip.end(), &c) != skip.end())
0139     return out;
0140   skip.push_back(&c);
0141
0142   int id = c.pdgId();
0143   const ParticleData *pd = pdt_->particle(id);
0144   assert(pd != nullptr);
0145   out += (pd->name() + printP4(c));
0146
0147   size_t validDau = 0, ndau = c.numberOfDaughters();
0148   for (size_t i = 0; i < ndau; ++i)
0149     if (accept(*c.daughter(i), skip))
0150       ++validDau;
0151   if (validDau == 0)
0152     return out;
0153
0154   out += " ->";
0155
0156   for (size_t i = 0; i < ndau; ++i) {
0157     const Candidate *d = c.daughter(i);
0158     if (accept(*d, skip)) {
0159       string dec = decay(*d, skip);
0160       if (dec.find("->", 0) != string::npos)
0161         out += (" ( " + dec + " )");
0162       else
0163         out += (" " + dec);
0164     }
0165   }
0166   return out;
0167 }
0168
0169 #include "FWCore/Framework/interface/MakerMacros.h"
0170
0171 DEFINE_FWK_MODULE(ParticleDecayDrawer);