GenTreeViewer/src/GenTreeViewer.cc

// -*- C++ -*-
//
// Package:    GenTreeViewer
// Class:      GenTreeViewer
// 
/**\class GenTreeViewer GenTreeViewer.cc DataInspection/GenTreeViewer/src/GenTreeViewer.cc

 Description: [one line class summary]

 Implementation:
     [Notes on implementation]
*/
//
// Original Author:  Yutaro Iiyama,512 1-005,+41227670489,
//         Created:  Sun Jul  8 21:52:04 CEST 2012
// $Id$
//
//


// system include files
#include <memory>
#include <algorithm>
#include <iomanip>
#include <sstream>

#include "TString.h"

// user include files
#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/EDAnalyzer.h"

#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/MakerMacros.h"

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

#include "FWCore/Utilities/interface/InputTag.h"

#include "DataFormats/Common/interface/Handle.h"

#include "DataFormats/HepMCCandidate/interface/GenParticle.h"
#include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h"
//
// class declaration
//

struct PNode {
  PNode() : mother(0), daughters(0), pdgId(0), status(0), mass(0.) {}
  PNode& operator=(PNode const& _rhs)
  {
    mother = _rhs.mother;
    daughters = _rhs.daughters;
    pdgId = _rhs.pdgId;
    status = _rhs.status;
    mass = _rhs.mass;
    pt = _rhs.pt;
    return *this;
  }
  PNode* mother;
  std::vector<PNode*> daughters;
  int pdgId;
  int status;
  double mass;
  double pt;
  std::string print(std::vector<bool>&);
  bool hasMother() { return mother != 0; }
};

std::string
PNode::print(std::vector<bool>& last)
{
  using namespace std;

  stringstream ss;
  ss << "+" << setw(5) << pdgId << " (" << setw(5) << fixed << setprecision(1) << pt << ") " << status << " ";
  for(unsigned iD(0); iD < daughters.size(); iD++){
    if(iD > 0){
      ss << endl;
      for(unsigned i(0); i < last.size(); i++)
        ss << (last[i] ? " " : "|") << "                ";
    }
    last.push_back(iD == daughters.size() - 1);
    ss << daughters[iD]->print(last);
    last.pop_back();
  }

  return ss.str();
}

PNode*
setDaughters(reco::GenParticle const* gen, std::map<reco::GenParticle const*, PNode*>& nodeMap, float _minPt)
{
  if(gen->status() == 1 && gen->pt() < _minPt) return 0;
  unsigned nD(gen->numberOfDaughters());
  if(nD == 0 && gen->status() != 1) return 0;

  PNode* node(new PNode);
  node->pdgId = gen->pdgId();
  node->status = gen->status();
  node->mass = gen->mass();
  node->pt = gen->pt();
  nodeMap[gen] = node;


  for(unsigned iD(0); iD < nD; iD++){
    reco::GenParticle const* daughter(static_cast<reco::GenParticle const*>(gen->daughter(iD)));
    if(nodeMap.find(daughter) != nodeMap.end()){
      PNode* dnode(nodeMap[daughter]);
      PNode* mother(dnode->mother);
      if(mother){
        if(mother == node) continue;

        int mpdg(std::abs(mother->pdgId));
        bool mhad((mpdg / 100) % 10 != 0 || mpdg == 21 || (mpdg > 80 && mpdg < 101));
        int pdg(std::abs(node->pdgId));
        bool had((pdg / 100) % 10 != 0 || pdg == 21 || (pdg > 80 && pdg < 101));
        bool takeAway(false);
        if(had && mhad)
          takeAway = node->pt > mother->pt;
        else if(!had && mhad)
          takeAway = true;

        if(takeAway){
          dnode->mother = node;
          node->daughters.push_back(dnode);
          std::vector<PNode*>::iterator dItr(std::find(mother->daughters.begin(), mother->daughters.end(), dnode));
          mother->daughters.erase(dItr);
        }
      }
      else{
        node->daughters.push_back(dnode);
        dnode->mother = node;
      }
    }
    else{
      PNode* dnode(setDaughters(daughter, nodeMap, _minPt));
      if(dnode){
        dnode->mother = node;
        node->daughters.push_back(dnode);
      }
    }
  }

  return node;
}

void
cleanDaughters(PNode* node)
{
  std::vector<PNode*> daughters(node->daughters);
  int motherPdg(std::abs(node->pdgId));
  for(unsigned iD(0); iD < daughters.size(); iD++){
    PNode* dnode(daughters[iD]);
    cleanDaughters(dnode);

    unsigned nGD(dnode->daughters.size());
    bool intermediateTerminal(nGD == 0 && dnode->status != 1);
    bool noDecay(nGD == 1 && dnode->pdgId == dnode->daughters.front()->pdgId);
    int pdg(std::abs(dnode->pdgId));
    bool hadronicIntermediate(dnode->status != 1 &&
                              ((pdg / 100) % 10 != 0 || pdg == 21 || (pdg > 80 && pdg < 101)));
    bool firstHeavyHadron((motherPdg / 1000) % 10 < 4 && (motherPdg / 100) % 10 < 4 &&
                          ((pdg / 1000) % 10 >= 4 || (pdg / 100) % 10 >= 4));
    bool lightDecayingToLight(pdg < 4);
    for(unsigned iGD(0); iGD < nGD; iGD++){
      if(std::abs(dnode->daughters[iGD]->pdgId) > 3){
        lightDecayingToLight = false;
        break;
      }
    }
    if(intermediateTerminal || noDecay || (hadronicIntermediate && !firstHeavyHadron) || lightDecayingToLight){
      node->daughters[iD] = 0;
      for(unsigned iGD(0); iGD < nGD; iGD++){
        node->daughters.push_back(dnode->daughters[iGD]);
        dnode->daughters[iGD]->mother = node;
        daughters.push_back(dnode->daughters[iGD]);
      }
    }
  }

  for(unsigned iD(0); iD < node->daughters.size(); iD++){
    if(!node->daughters[iD]){
      node->daughters.erase(node->daughters.begin() + iD);
      iD--;
    }
  }
}


class GenTreeViewer : public edm::EDAnalyzer {
   public:
      explicit GenTreeViewer(const edm::ParameterSet&);
      ~GenTreeViewer();

      static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);


   private:
      virtual void beginJob() ;
      virtual void analyze(const edm::Event&, const edm::EventSetup&);
      virtual void endJob() ;

      virtual void beginRun(edm::Run const&, edm::EventSetup const&);
      virtual void endRun(edm::Run const&, edm::EventSetup const&);
      virtual void beginLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&);
      virtual void endLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&);

      // ----------member data ---------------------------

  edm::InputTag genParticlesTag_;
  int cleaningMode_;
  float minPt_;
};

//
// constants, enums and typedefs
//

//
// static data member definitions
//

//
// constructors and destructor
//
GenTreeViewer::GenTreeViewer(const edm::ParameterSet& iConfig) :
  genParticlesTag_(iConfig.getParameter<edm::InputTag>("genParticlesTag")),
  cleaningMode_(iConfig.getParameter<int>("cleaningMode")),
  minPt_(iConfig.getParameter<double>("minPt"))
{
   //now do what ever initialization is needed

}


GenTreeViewer::~GenTreeViewer()
{
 
   // do anything here that needs to be done at desctruction time
   // (e.g. close files, deallocate resources etc.)

}


//
// member functions
//

// ------------ method called for each event  ------------
void
GenTreeViewer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
{
   using namespace edm;

   Handle<reco::GenParticleCollection> gpHndl;
   if(!iEvent.getByLabel(genParticlesTag_, gpHndl)) return;

   std::map<reco::GenParticle const*, PNode*> nodeMap;
   std::vector<PNode*> rootNodes;

   for(reco::GenParticleCollection::const_iterator genItr(gpHndl->begin()); genItr != gpHndl->end(); ++genItr){
     reco::GenParticle const& gen(*genItr);

     if(gen.numberOfMothers() == 0){
       PNode* node(setDaughters(&gen, nodeMap, minPt_));
       if(node) rootNodes.push_back(node);
     }
   }

   if(cleaningMode_ == 0 || cleaningMode_ == 2){
     std::cout << "=== FULL DECAY TREE ===" << std::endl << std::endl;
     for(unsigned iN(0); iN < rootNodes.size(); iN++){
       std::vector<bool> last(1, true);
       std::cout << rootNodes[iN]->print(last);
       std::cout << std::endl;
     }
   }

   if(cleaningMode_ == 1 || cleaningMode_ == 2){
     std::cout << "--- CLEANED DECAY TREE ---" << std::endl << std::endl;
     for(unsigned iN(0); iN < rootNodes.size(); iN++){
       cleanDaughters(rootNodes[iN]);
       std::vector<bool> last(1, true);
       std::cout << rootNodes[iN]->print(last);
       std::cout << std::endl;
     }
   }

   for(std::map<reco::GenParticle const*, PNode*>::iterator nItr(nodeMap.begin()); nItr != nodeMap.end(); ++nItr)
     delete nItr->second;

   std::cout << std::endl;
}


// ------------ method called once each job just before starting event loop  ------------
void 
GenTreeViewer::beginJob()
{
}

// ------------ method called once each job just after ending the event loop  ------------
void 
GenTreeViewer::endJob() 
{
}

// ------------ method called when starting to processes a run  ------------
void 
GenTreeViewer::beginRun(edm::Run const&, edm::EventSetup const&)
{
}

// ------------ method called when ending the processing of a run  ------------
void 
GenTreeViewer::endRun(edm::Run const&, edm::EventSetup const&)
{
}

// ------------ method called when starting to processes a luminosity block  ------------
void 
GenTreeViewer::beginLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&)
{
}

// ------------ method called when ending the processing of a luminosity block  ------------
void 
GenTreeViewer::endLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&)
{
}

// ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
void
GenTreeViewer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
  //The following says we do not know what parameters are allowed so do no validation
  // Please change this to state exactly what you do use, even if it is no parameters
  edm::ParameterSetDescription desc;
  desc.setUnknown();
  descriptions.addDefault(desc);
}

//define this as a plug-in
DEFINE_FWK_MODULE(GenTreeViewer);
Revision:	1.1
Committed:	Fri Jul 13 17:16:04 2012 UTC (12 years, 9 months ago) by yiiyama
Content type:	text/plain
Branch:	MAIN
Log Message:	first commit
#	User	Rev	Content
1	yiiyama	1.1	// -- C++ --
2			//
3			// Package: GenTreeViewer
4			// Class: GenTreeViewer
5			//
6			/**\class GenTreeViewer GenTreeViewer.cc DataInspection/GenTreeViewer/src/GenTreeViewer.cc
7
8			Description: [one line class summary]
9
10			Implementation:
11			[Notes on implementation]
12			*/
13			//
14			// Original Author: Yutaro Iiyama,512 1-005,+41227670489,
15			// Created: Sun Jul 8 21:52:04 CEST 2012
16			// $Id$
17			//
18			//
19
20
21			// system include files
22			#include <memory>
23			#include <algorithm>
24			#include <iomanip>
25			#include <sstream>
26
27			#include "TString.h"
28
29			// user include files
30			#include "FWCore/Framework/interface/Frameworkfwd.h"
31			#include "FWCore/Framework/interface/EDAnalyzer.h"
32
33			#include "FWCore/Framework/interface/Event.h"
34			#include "FWCore/Framework/interface/MakerMacros.h"
35
36			#include "FWCore/ParameterSet/interface/ParameterSet.h"
37
38			#include "FWCore/Utilities/interface/InputTag.h"
39
40			#include "DataFormats/Common/interface/Handle.h"
41
42			#include "DataFormats/HepMCCandidate/interface/GenParticle.h"
43			#include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h"
44			//
45			// class declaration
46			//
47
48			struct PNode {
49			PNode() : mother(0), daughters(0), pdgId(0), status(0), mass(0.) {}
50			PNode& operator=(PNode const& _rhs)
51			{
52			mother = _rhs.mother;
53			daughters = _rhs.daughters;
54			pdgId = _rhs.pdgId;
55			status = _rhs.status;
56			mass = _rhs.mass;
57			pt = _rhs.pt;
58			return *this;
59			}
60			PNode* mother;
61			std::vector<PNode*> daughters;
62			int pdgId;
63			int status;
64			double mass;
65			double pt;
66			std::string print(std::vector<bool>&);
67			bool hasMother() { return mother != 0; }
68			};
69
70			std::string
71			PNode::print(std::vector<bool>& last)
72			{
73			using namespace std;
74
75			stringstream ss;
76			ss << "+" << setw(5) << pdgId << " (" << setw(5) << fixed << setprecision(1) << pt << ") " << status << " ";
77			for(unsigned iD(0); iD < daughters.size(); iD++){
78			if(iD > 0){
79			ss << endl;
80			for(unsigned i(0); i < last.size(); i++)
81			ss << (last[i] ? " " : "\|") << " ";
82			}
83			last.push_back(iD == daughters.size() - 1);
84			ss << daughters[iD]->print(last);
85			last.pop_back();
86			}
87
88			return ss.str();
89			}
90
91			PNode*
92			setDaughters(reco::GenParticle const* gen, std::map<reco::GenParticle const, PNode>& nodeMap, float _minPt)
93			{
94			if(gen->status() == 1 && gen->pt() < _minPt) return 0;
95			unsigned nD(gen->numberOfDaughters());
96			if(nD == 0 && gen->status() != 1) return 0;
97
98			PNode* node(new PNode);
99			node->pdgId = gen->pdgId();
100			node->status = gen->status();
101			node->mass = gen->mass();
102			node->pt = gen->pt();
103			nodeMap[gen] = node;
104
105
106			for(unsigned iD(0); iD < nD; iD++){
107			reco::GenParticle const* daughter(static_cast<reco::GenParticle const*>(gen->daughter(iD)));
108			if(nodeMap.find(daughter) != nodeMap.end()){
109			PNode* dnode(nodeMap[daughter]);
110			PNode* mother(dnode->mother);
111			if(mother){
112			if(mother == node) continue;
113
114			int mpdg(std::abs(mother->pdgId));
115			bool mhad((mpdg / 100) % 10 != 0 \|\| mpdg == 21 \|\| (mpdg > 80 && mpdg < 101));
116			int pdg(std::abs(node->pdgId));
117			bool had((pdg / 100) % 10 != 0 \|\| pdg == 21 \|\| (pdg > 80 && pdg < 101));
118			bool takeAway(false);
119			if(had && mhad)
120			takeAway = node->pt > mother->pt;
121			else if(!had && mhad)
122			takeAway = true;
123
124			if(takeAway){
125			dnode->mother = node;
126			node->daughters.push_back(dnode);
127			std::vector<PNode*>::iterator dItr(std::find(mother->daughters.begin(), mother->daughters.end(), dnode));
128			mother->daughters.erase(dItr);
129			}
130			}
131			else{
132			node->daughters.push_back(dnode);
133			dnode->mother = node;
134			}
135			}
136			else{
137			PNode* dnode(setDaughters(daughter, nodeMap, _minPt));
138			if(dnode){
139			dnode->mother = node;
140			node->daughters.push_back(dnode);
141			}
142			}
143			}
144
145			return node;
146			}
147
148			void
149			cleanDaughters(PNode* node)
150			{
151			std::vector<PNode*> daughters(node->daughters);
152			int motherPdg(std::abs(node->pdgId));
153			for(unsigned iD(0); iD < daughters.size(); iD++){
154			PNode* dnode(daughters[iD]);
155			cleanDaughters(dnode);
156
157			unsigned nGD(dnode->daughters.size());
158			bool intermediateTerminal(nGD == 0 && dnode->status != 1);
159			bool noDecay(nGD == 1 && dnode->pdgId == dnode->daughters.front()->pdgId);
160			int pdg(std::abs(dnode->pdgId));
161			bool hadronicIntermediate(dnode->status != 1 &&
162			((pdg / 100) % 10 != 0 \|\| pdg == 21 \|\| (pdg > 80 && pdg < 101)));
163			bool firstHeavyHadron((motherPdg / 1000) % 10 < 4 && (motherPdg / 100) % 10 < 4 &&
164			((pdg / 1000) % 10 >= 4 \|\| (pdg / 100) % 10 >= 4));
165			bool lightDecayingToLight(pdg < 4);
166			for(unsigned iGD(0); iGD < nGD; iGD++){
167			if(std::abs(dnode->daughters[iGD]->pdgId) > 3){
168			lightDecayingToLight = false;
169			break;
170			}
171			}
172			if(intermediateTerminal \|\| noDecay \|\| (hadronicIntermediate && !firstHeavyHadron) \|\| lightDecayingToLight){
173			node->daughters[iD] = 0;
174			for(unsigned iGD(0); iGD < nGD; iGD++){
175			node->daughters.push_back(dnode->daughters[iGD]);
176			dnode->daughters[iGD]->mother = node;
177			daughters.push_back(dnode->daughters[iGD]);
178			}
179			}
180			}
181
182			for(unsigned iD(0); iD < node->daughters.size(); iD++){
183			if(!node->daughters[iD]){
184			node->daughters.erase(node->daughters.begin() + iD);
185			iD--;
186			}
187			}
188			}
189
190
191			class GenTreeViewer : public edm::EDAnalyzer {
192			public:
193			explicit GenTreeViewer(const edm::ParameterSet&);
194			~GenTreeViewer();
195
196			static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
197
198
199			private:
200			virtual void beginJob() ;
201			virtual void analyze(const edm::Event&, const edm::EventSetup&);
202			virtual void endJob() ;
203
204			virtual void beginRun(edm::Run const&, edm::EventSetup const&);
205			virtual void endRun(edm::Run const&, edm::EventSetup const&);
206			virtual void beginLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&);
207			virtual void endLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&);
208
209			// ----------member data ---------------------------
210
211			edm::InputTag genParticlesTag_;
212			int cleaningMode_;
213			float minPt_;
214			};
215
216			//
217			// constants, enums and typedefs
218			//
219
220			//
221			// static data member definitions
222			//
223
224			//
225			// constructors and destructor
226			//
227			GenTreeViewer::GenTreeViewer(const edm::ParameterSet& iConfig) :
228			genParticlesTag_(iConfig.getParameter<edm::InputTag>("genParticlesTag")),
229			cleaningMode_(iConfig.getParameter<int>("cleaningMode")),
230			minPt_(iConfig.getParameter<double>("minPt"))
231			{
232			//now do what ever initialization is needed
233
234			}
235
236
237			GenTreeViewer::~GenTreeViewer()
238			{
239
240			// do anything here that needs to be done at desctruction time
241			// (e.g. close files, deallocate resources etc.)
242
243			}
244
245
246			//
247			// member functions
248			//
249
250			// ------------ method called for each event ------------
251			void
252			GenTreeViewer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
253			{
254			using namespace edm;
255
256			Handle<reco::GenParticleCollection> gpHndl;
257			if(!iEvent.getByLabel(genParticlesTag_, gpHndl)) return;
258
259			std::map<reco::GenParticle const, PNode> nodeMap;
260			std::vector<PNode*> rootNodes;
261
262			for(reco::GenParticleCollection::const_iterator genItr(gpHndl->begin()); genItr != gpHndl->end(); ++genItr){
263			reco::GenParticle const& gen(*genItr);
264
265			if(gen.numberOfMothers() == 0){
266			PNode* node(setDaughters(&gen, nodeMap, minPt_));
267			if(node) rootNodes.push_back(node);
268			}
269			}
270
271			if(cleaningMode_ == 0 \|\| cleaningMode_ == 2){
272			std::cout << "=== FULL DECAY TREE ===" << std::endl << std::endl;
273			for(unsigned iN(0); iN < rootNodes.size(); iN++){
274			std::vector<bool> last(1, true);
275			std::cout << rootNodes[iN]->print(last);
276			std::cout << std::endl;
277			}
278			}
279
280			if(cleaningMode_ == 1 \|\| cleaningMode_ == 2){
281			std::cout << "--- CLEANED DECAY TREE ---" << std::endl << std::endl;
282			for(unsigned iN(0); iN < rootNodes.size(); iN++){
283			cleanDaughters(rootNodes[iN]);
284			std::vector<bool> last(1, true);
285			std::cout << rootNodes[iN]->print(last);
286			std::cout << std::endl;
287			}
288			}
289
290			for(std::map<reco::GenParticle const, PNode>::iterator nItr(nodeMap.begin()); nItr != nodeMap.end(); ++nItr)
291			delete nItr->second;
292
293			std::cout << std::endl;
294			}
295
296
297			// ------------ method called once each job just before starting event loop ------------
298			void
299			GenTreeViewer::beginJob()
300			{
301			}
302
303			// ------------ method called once each job just after ending the event loop ------------
304			void
305			GenTreeViewer::endJob()
306			{
307			}
308
309			// ------------ method called when starting to processes a run ------------
310			void
311			GenTreeViewer::beginRun(edm::Run const&, edm::EventSetup const&)
312			{
313			}
314
315			// ------------ method called when ending the processing of a run ------------
316			void
317			GenTreeViewer::endRun(edm::Run const&, edm::EventSetup const&)
318			{
319			}
320
321			// ------------ method called when starting to processes a luminosity block ------------
322			void
323			GenTreeViewer::beginLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&)
324			{
325			}
326
327			// ------------ method called when ending the processing of a luminosity block ------------
328			void
329			GenTreeViewer::endLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&)
330			{
331			}
332
333			// ------------ method fills 'descriptions' with the allowed parameters for the module ------------
334			void
335			GenTreeViewer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
336			//The following says we do not know what parameters are allowed so do no validation
337			// Please change this to state exactly what you do use, even if it is no parameters
338			edm::ParameterSetDescription desc;
339			desc.setUnknown();
340			descriptions.addDefault(desc);
341			}
342
343			//define this as a plug-in
344			DEFINE_FWK_MODULE(GenTreeViewer);