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: GenTreeViewer.cc,v 1.1 2012/07/13 17:16:04 yiiyama Exp $
//
//


// 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(7) << 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.2
Committed:	Mon Oct 8 12:31:54 2012 UTC (12 years, 6 months ago) by yiiyama
Content type:	text/plain
Branch:	MAIN
CVS Tags:	HEAD
Changes since 1.1:	+3 -3 lines
Log Message:	Correct offset Added RA3 version
#	Content
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: GenTreeViewer.cc,v 1.1 2012/07/13 17:16:04 yiiyama Exp $
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(7) << 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);