ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitProd/TreeFiller/src/FillGenParts.cc
Revision: 1.3
Committed: Tue Jun 3 07:21:45 2008 UTC (16 years, 11 months ago) by paus
Content type: text/plain
Branch: MAIN
Changes since 1.2: +5 -3 lines
Log Message:
Minor cleanup and adding examples (see MitSoft TWiki).

File Contents

# User Rev Content
1 paus 1.3 // $Id: FillGenParts.cc,v 1.2 2008/06/02 04:52:56 loizides Exp $
2 loizides 1.1
3     #include "MitProd/TreeFiller/interface/FillGenParts.h"
4     #include "FWCore/MessageLogger/interface/MessageLogger.h"
5     #include "FWCore/Framework/interface/ESHandle.h"
6     #include "DataFormats/Common/interface/Handle.h"
7     #include "FWCore/ServiceRegistry/interface/Service.h"
8    
9     #include "SimDataFormats/HepMCProduct/interface/HepMCProduct.h"
10     #include "MitAna/DataTree/interface/Particle.h"
11    
12     #include "TLorentzVector.h"
13    
14     using namespace std;
15     using namespace edm;
16     using namespace mithep;
17    
18     //-------------------------------------------------------------------------------------------------
19     FillGenParts::FillGenParts(const edm::ParameterSet &iConfig)
20     : anaMCSource_(iConfig.getUntrackedParameter<string>("anaMCSource" , "source"))
21     {
22     genParticles_ = new mithep::Vector<mithep::Particle>();
23     }
24    
25     //-------------------------------------------------------------------------------------------------
26     FillGenParts::~FillGenParts()
27     {
28     cout << " Fillgenpars done " <<endl;
29     }
30    
31     //-------------------------------------------------------------------------------------------------
32     void FillGenParts::analyze(const edm::Event &theEvent,
33     const edm::EventSetup &iSetup)
34     {
35     genParticles_->Reset();
36    
37     Handle<HepMCProduct> theMCProduct;
38     try {
39     theEvent.getByLabel(anaMCSource_, theMCProduct);
40     } catch (cms::Exception& ex) {
41     printf("Error! can't get collection with label %s\n",anaMCSource_.c_str());
42     }
43    
44     const HepMC::GenEvent GenEvent = theMCProduct->getHepMCData();
45    
46     int nGen = 0;
47     for (HepMC::GenEvent::particle_const_iterator pgen =
48     GenEvent.particles_begin();
49     pgen != GenEvent.particles_end(); ++pgen) {
50    
51     HepMC::GenParticle* mcPart = (*pgen);
52     if(!mcPart) continue;
53    
54 loizides 1.2 mithep::Particle genParticle(mcPart->momentum().x(),mcPart->momentum().y(),
55     mcPart->momentum().z(),mcPart->momentum().e());
56 loizides 1.1
57     //printf("ngen %d\n",nGen);
58 paus 1.3
59 loizides 1.1 #if 0
60     genParticle.setPdgID(MCPart->pdg_id());
61     genParticle.setBarCode(MCPart->barcode()-1);
62     genParticle.setStatus(MCPart->status());
63    
64     HepMC::GenVertex * momVert = MCPart->production_vertex();
65     //genParticle.setMother(-1);
66 paus 1.3 if (momVert) {
67     if (momVert->particles_in_size() == 1) {
68 loizides 1.1 HepMC::GenVertex::particles_in_const_iterator mom = momVert->particles_in_const_begin();
69     genParticle.setMother((*mom)->barcode() - 1);
70     } else {genParticle.setMother(-1);}
71     } else {genParticle.setMother(-1);}
72     //if(MCPart->mother()) genParticle.setMother(MCPart->mother()->barcode()-1);
73     //else genParticle.setMother(-1);
74     genParticle.setPtr(nGen);
75     #endif
76 paus 1.3
77 loizides 1.1 genParticles_->Add(genParticle);
78     nGen++;
79     }
80    
81     //genParticles_->resize(nGen);
82     }
83    
84     //-------------------------------------------------------------------------------------------------
85     void FillGenParts::beginJob(edm::EventSetup const &iEvent)
86     {
87     Service<TreeService> ts;
88     TreeWriter *tws = ts->get();
89     if(!tws) {
90     //todo error
91     return;
92     }
93    
94     tws->AddBranch("Particles","mithep::Vector<mithep::Particle>",&genParticles_,32000,99);
95     }
96    
97     //-------------------------------------------------------------------------------------------------
98     void FillGenParts::endJob()
99     {
100     edm::LogInfo("FillGenParts::endJob") << "Ending Job" << endl;
101     }