ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/Morgan/src/TotoAnalyzer.cc
Revision: 1.4
Committed: Mon Jun 2 16:31:24 2008 UTC (16 years, 11 months ago) by mlethuil
Content type: text/plain
Branch: MAIN
Changes since 1.3: +47 -17 lines
Log Message:
Add TRootTrack collection with PixelHits info
Tracker isolation for TRootPhoton

File Contents

# User Rev Content
1 mlethuil 1.1 #include "UserCode/Morgan/interface/TotoAnalyzer.h"
2    
3     using namespace std;
4     using namespace reco;
5     using namespace edm;
6    
7    
8     TotoAnalyzer::TotoAnalyzer(const edm::ParameterSet& iConfig)
9     {
10     myConfig = iConfig.getParameter<ParameterSet>("myConfig");
11     }
12    
13    
14     TotoAnalyzer::~TotoAnalyzer()
15     {}
16    
17    
18    
19     // ------------ method called once each job just before starting event loop ------------
20     void TotoAnalyzer::beginJob(const edm::EventSetup&)
21     {
22    
23     // Load Config parameters
24     verbosity = myConfig.getUntrackedParameter<int>("verbosity", 0);
25     rootFileName_ = myConfig.getUntrackedParameter<string>("RootFileName","noname.root");
26     doMC = myConfig.getUntrackedParameter<bool>("doMC",false);
27     doSignalMC = myConfig.getUntrackedParameter<bool>("doSignalMC",false);
28     doTrack = myConfig.getUntrackedParameter<bool>("doTrack",false);
29     doJet = myConfig.getUntrackedParameter<bool>("doJet",false);
30     doMuon = myConfig.getUntrackedParameter<bool>("doMuon",false);
31     doElectron = myConfig.getUntrackedParameter<bool>("doElectron",false);
32     doPhoton = myConfig.getUntrackedParameter<bool>("doPhoton",false);
33     doCluster = myConfig.getUntrackedParameter<bool>("doCluster",false);
34 mlethuil 1.4 doPhotonIsolation = myConfig.getUntrackedParameter<bool>("doPhotonIsolation",false);
35 mlethuil 1.3 doPhotonConversion = myConfig.getUntrackedParameter<bool>("doPhotonConversion",false);
36 mlethuil 1.1
37     nTotEvt_ = 0;
38    
39     // initialize root output file
40     if(verbosity>0) cout << "New RootFile " << rootFileName_.c_str() << " is created" << endl;
41     rootFile_ = new TFile(rootFileName_.c_str(), "recreate");
42     rootFile_->cd();
43     rootTree_ = new TTree("tree", "RootTree a Toto");
44    
45     rootEvent = 0;
46     rootTree_->Branch ("Event", "TRootEvent", &rootEvent);
47    
48     if(doSignalMC)
49     {
50     if(verbosity>0) cout << "Signal generator info will be added to rootuple" << endl;
51     rootMCSignalEvent = 0;
52     cout << "Create MuMuGamma Branch" << endl;
53     rootTree_->Branch ("MuMuGamma", "TRootSignalEvent", &rootMCSignalEvent);
54     }
55    
56     if(doMC)
57     {
58     if(verbosity>0) cout << "MCParticles info will be added to rootuple" << endl;
59     mcParticles = new TClonesArray("TRootParticle", 1000);
60     rootTree_->Branch ("MCParticles", "TClonesArray", &mcParticles);
61     }
62    
63     if(doTrack)
64     {
65     if(verbosity>0) cout << "Tracks info will be added to rootuple" << endl;
66 mlethuil 1.4 tracks = new TClonesArray("TRootTrack", 1000);
67 mlethuil 1.1 rootTree_->Branch ("Tracks", "TClonesArray", &tracks);
68     }
69    
70     if(doJet)
71     {
72     if(verbosity>0) cout << "Jets info will be added to rootuple" << endl;
73     jets = new TClonesArray("TRootJet", 1000);
74     rootTree_->Branch ("Jets", "TClonesArray", &jets);
75     }
76    
77     if(doMuon)
78     {
79     if(verbosity>0) cout << "Muons info will be added to rootuple" << endl;
80     muons = new TClonesArray("TRootMuon", 1000);
81     rootTree_->Branch ("Muons", "TClonesArray", &muons);
82     }
83    
84     if(doElectron)
85     {
86     if(verbosity>0) cout << "Electrons info will be added to rootuple" << endl;
87     electrons = new TClonesArray("TRootElectron", 1000);
88     rootTree_->Branch ("Electron", "TClonesArray", &electrons);
89     }
90    
91     if(doPhoton)
92     {
93     if(verbosity>0) cout << "Photons info will be added to rootuple" << endl;
94     uncorrectedPhotons = new TClonesArray("TRootPhoton", 1000);
95     rootTree_->Branch ("uncorrectedPhotons", "TClonesArray", &uncorrectedPhotons);
96     correctedPhotons = new TClonesArray("TRootPhoton", 1000);
97     rootTree_->Branch ("correctedPhotons", "TClonesArray", &correctedPhotons);
98     }
99    
100     if(doCluster)
101     {
102     if(verbosity>0) cout << "ECAL Clusters info will be added to rootuple" << endl;
103     clusters = new TClonesArray("TRootCluster", 1000);
104     rootTree_->Branch ("BasicClusters", "TClonesArray", &clusters);
105     superClusters = new TClonesArray("TRootSuperCluster", 1000);
106     rootTree_->Branch ("SuperClusters", "TClonesArray", &superClusters);
107     }
108 mlethuil 1.3
109     if(doPhotonConversion)
110     {
111     if(verbosity>0) cout << "Conversion Tracks info will be added to rootuple" << endl;
112     conversionTracks = new TClonesArray("TRootParticle", 1000);
113     rootTree_->Branch ("ConversionTracks", "TClonesArray", &conversionTracks);
114     }
115 mlethuil 1.1 }
116    
117    
118     // ------------ method called once each job just after ending the event loop ------------
119     void TotoAnalyzer::endJob()
120     {
121     std::cout << "Total number of events: " << nTotEvt_ << std::endl;
122     std::cout << "Closing rootfile " << rootFile_->GetName() << std::endl;
123     rootFile_->Write();
124     rootFile_->Close();
125     delete mcParticles;
126     delete tracks;
127     delete jets;
128     delete muons;
129     delete electrons;
130     delete correctedPhotons;
131     delete uncorrectedPhotons;
132     delete clusters;
133     delete superClusters;
134 mlethuil 1.3 delete conversionTracks;
135 mlethuil 1.1 }
136    
137    
138    
139     // ------------ method called to for each event ------------
140     void TotoAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
141     {
142    
143     rootFile_->cd();
144     nTotEvt_++;
145     if( (verbosity>1) || (verbosity>0 && nTotEvt_%10==0 && nTotEvt_<=100) || (verbosity>0 && nTotEvt_%100==0 && nTotEvt_>100) )
146     cout << endl << endl << "####### TotoAnalyzer - Event " << nTotEvt_ << " #######" << endl;
147    
148    
149     // Global Event Infos
150     rootEvent = new TRootEvent();
151     rootEvent->setNb(nTotEvt_);
152    
153    
154     // MC Info
155     if(doMC)
156     {
157     rootMCSignalEvent = new TRootSignalEvent();
158     if(verbosity>1) cout << "Analysing MC info..." << endl;
159     myMCAnalyzer = new MCAnalyzer(&myConfig);
160     //myMCAnalyzer->SetVerbosity(1);
161     myMCAnalyzer->Process(iEvent, iSetup, rootEvent, rootMCSignalEvent, mcParticles);
162     delete myMCAnalyzer;
163     }
164    
165    
166     // Get Primary Vertices
167     if(verbosity>1) cout << "Analysing primary vertices collection..." << endl;
168     edm::Handle< reco::VertexCollection > recoVertex;
169     iEvent.getByLabel("offlinePrimaryVerticesFromCTFTracks", recoVertex);
170     if(verbosity>1) std::cout << "Number of primary vertices in the event: " << recoVertex->size() << std::endl;
171     for (unsigned int j=0; j<recoVertex->size(); j++)
172     {
173     if(verbosity>1) cout << "offlinePrimaryVerticesFromCTFTracks vx=" << (*recoVertex)[j].x() << " vy=" << (*recoVertex)[j].y() << " vz=" << (*recoVertex)[j].z() << endl;
174     rootEvent->addPrimaryVertex( TVector3( (*recoVertex)[j].x(), (*recoVertex)[j].y(), (*recoVertex)[j].z() ) );
175     }
176    
177    
178     // Tracks
179     if(doTrack)
180     {
181     if(verbosity>1) cout << "Analysing tracks collection..." << endl;
182     edm::Handle<reco::TrackCollection> recoTracks;
183     iEvent.getByLabel("ctfWithMaterialTracks", recoTracks);
184     if(verbosity>1) std::cout << "Number of tracks in the event: " << recoTracks->size() << std::endl;
185 mlethuil 1.4
186 mlethuil 1.1 for (unsigned int j=0; j<recoTracks->size(); j++)
187     {
188 mlethuil 1.4 const reco::HitPattern& hit = (*recoTracks)[j].hitPattern();
189    
190     //TRootTrack localTrack( (*recoTracks)[j].px(), (*recoTracks)[j].py(), (*recoTracks)[j].pz(), (*recoTracks)[j].p(), (*recoTracks)[j].vx(), (*recoTracks)[j].vy(), (*recoTracks)[j].vz(), 0, (*recoTracks)[j].charge() );
191    
192     TRootTrack localTrack(
193     (*recoTracks)[j].px()
194     ,(*recoTracks)[j].py()
195     ,(*recoTracks)[j].pz()
196     ,(*recoTracks)[j].p()
197     ,(*recoTracks)[j].vx()
198     ,(*recoTracks)[j].vy()
199     ,(*recoTracks)[j].vz()
200     ,0
201     ,(*recoTracks)[j].charge()
202     ,hit.numberOfValidPixelHits()
203     ,hit.numberOfValidTrackerHits()
204     ,(*recoTracks)[j].chi2()
205     ,(*recoTracks)[j].d0()
206     ,(*recoTracks)[j].d0Error()
207     ,(*recoTracks)[j].dz()
208     ,(*recoTracks)[j].dzError()
209     );
210    
211     new( (*tracks)[j] ) TRootTrack(localTrack);
212 mlethuil 1.1 if(verbosity>2) cout << "Track - "<< localTrack << endl;
213 mlethuil 1.4
214 mlethuil 1.1 }
215     }
216    
217    
218     // Calo Jet
219     if(doJet)
220     {
221     if(verbosity>1) cout << "Analysing jets collection..." << endl;
222     edm::Handle< reco::CaloJetCollection > recoJets;
223     iEvent.getByLabel("iterativeCone5CaloJets", recoJets);
224     if(verbosity>1) std::cout << "Number of jets in the event: " << recoJets->size() << std::endl;
225     for (unsigned int j=0; j<recoJets->size(); j++)
226     {
227     if(verbosity>2) cout << "jet pt="<< (*recoJets)[j].pt() << endl;
228     new( (*jets)[j] ) TRootJet( (*recoJets)[j].px(), (*recoJets)[j].py(), (*recoJets)[j].pz(), (*recoJets)[j].energy() );
229     }
230     }
231    
232    
233     // Muons
234     if(doMuon)
235     {
236     Float_t sintheta = 0.;
237     if(verbosity>1) cout << "Analysing muons collection..." << endl;
238     edm::Handle < reco::MuonCollection > recoMuons;
239     iEvent.getByLabel ("muons", recoMuons);
240     if(verbosity>1) std::cout << "Number of muons in the event: " << recoMuons->size() << std::endl;
241     for (unsigned int j=0; j<recoMuons->size(); j++)
242     {
243     TRootMuon localMuon( (*recoMuons)[j].px(), (*recoMuons)[j].py(), (*recoMuons)[j].pz(), (*recoMuons)[j].energy(), (*recoMuons)[j].vx(), (*recoMuons)[j].vy(), (*recoMuons)[j].vz(), 13, (*recoMuons)[j].charge() );
244     sintheta = sin( localMuon.Theta() );
245     localMuon.setCaloEnergy( (*recoMuons)[j].getCalEnergy().em * sintheta, (*recoMuons)[j].getCalEnergy().emS9 * sintheta,
246     (*recoMuons)[j].getCalEnergy().had * sintheta, (*recoMuons)[j].getCalEnergy().hadS9 * sintheta,
247     (*recoMuons)[j].getCalEnergy().ho * sintheta, (*recoMuons)[j].getCalEnergy().hoS9 * sintheta, (*recoMuons)[j].getCaloCompatibility() );
248     localMuon.setIsoR03( (*recoMuons)[j].getIsolationR03().emEt, (*recoMuons)[j].getIsolationR03().hadEt, (*recoMuons)[j].getIsolationR03().hoEt,
249     (*recoMuons)[j].getIsolationR03().sumPt, (*recoMuons)[j].getIsolationR03().nTracks, (*recoMuons)[j].getIsolationR03().nJets);
250     localMuon.setIsoR05( (*recoMuons)[j].getIsolationR05().emEt, (*recoMuons)[j].getIsolationR05().hadEt, (*recoMuons)[j].getIsolationR05().hoEt,
251     (*recoMuons)[j].getIsolationR05().sumPt, (*recoMuons)[j].getIsolationR05().nTracks, (*recoMuons)[j].getIsolationR05().nJets);
252     new( (*muons)[j] ) TRootMuon(localMuon);
253     if(verbosity>2) cout << localMuon << endl;
254     }
255     }
256    
257    
258     // Electrons
259     if(doElectron)
260     {
261     if(verbosity>1) cout << "Analysing electrons collection..." << endl;
262     edm::Handle< reco::PixelMatchGsfElectronCollection > recoElectrons;
263     iEvent.getByLabel("pixelMatchGsfElectrons", recoElectrons);
264     if(verbosity>1) std::cout << "Number of electrons in the event: " << recoElectrons->size() << std::endl;
265     for (unsigned int j=0; j<recoElectrons->size(); j++)
266     {
267     float elec_e = ((*recoElectrons)[j].superCluster())->energy();
268     float elec_px = (*recoElectrons)[j].px() * elec_e / (*recoElectrons)[j].p();
269     float elec_py = (*recoElectrons)[j].py() * elec_e / (*recoElectrons)[j].p();
270     float elec_pz = (*recoElectrons)[j].pz() * elec_e / (*recoElectrons)[j].p();
271     TRootElectron localElectron( elec_px, elec_py, elec_pz, elec_e, (*recoElectrons)[j].vx(), (*recoElectrons)[j].vy(), (*recoElectrons)[j].vz(), 11, (*recoElectrons)[j].charge() );
272     new( (*electrons)[j] ) TRootElectron(localElectron);
273     if(verbosity>2) cout << localElectron << endl;
274     }
275     }
276    
277    
278     // Photons
279     if(doPhoton)
280     {
281     if(verbosity>1) cout << "Analysing photons collection..." << endl;
282     myPhotonAnalyzer = new PhotonAnalyzer();
283     myPhotonAnalyzer->SetVerbosity(verbosity);
284     myPhotonAnalyzer->Process(iEvent, rootEvent, uncorrectedPhotons, "photons");
285     myPhotonAnalyzer->Process(iEvent, rootEvent, correctedPhotons, "correctedPhotons");
286     delete myPhotonAnalyzer;
287     }
288    
289    
290     // Ecal Clusters
291     if(doCluster)
292     {
293     if(verbosity>1) cout << "Analysing BasicClusters collection..." << endl;
294     myClusterAnalyzer = new ClusterAnalyzer();
295     myClusterAnalyzer->SetVerbosity(verbosity);
296     myClusterAnalyzer->Process(iEvent, rootEvent, clusters, "islandBasicClusters","islandBarrelBasicClusters", "islandBarrelShape", 110);
297     myClusterAnalyzer->Process(iEvent, rootEvent, clusters, "islandBasicClusters", "islandEndcapBasicClusters", "islandEndcapShape", 120);
298     myClusterAnalyzer->Process(iEvent, rootEvent, clusters, "hybridSuperClusters", "", "", 210);
299     delete myClusterAnalyzer;
300    
301     if(verbosity>1) cout << "Analysing SuperClusters collection..." << endl;
302     mySClusterAnalyzer = new SuperClusterAnalyzer();
303     mySClusterAnalyzer->SetVerbosity(verbosity);
304     mySClusterAnalyzer->Process(iEvent, rootEvent, superClusters, "islandSuperClusters","islandBarrelSuperClusters",110);
305     mySClusterAnalyzer->Process(iEvent, rootEvent, superClusters, "islandSuperClusters","islandEndcapSuperClusters",120);
306     mySClusterAnalyzer->Process(iEvent, rootEvent, superClusters, "correctedIslandEndcapSuperClusters","",121);
307     mySClusterAnalyzer->Process(iEvent, rootEvent, superClusters, "correctedEndcapSuperClustersWithPreshower","",122);
308     mySClusterAnalyzer->Process(iEvent, rootEvent, superClusters, "hybridSuperClusters","",210);
309     mySClusterAnalyzer->Process(iEvent, rootEvent, superClusters, "correctedHybridSuperClusters","",211);
310     delete mySClusterAnalyzer;
311     }
312    
313     /*
314     cout << "Nb of BC110=" << rootEvent->nBasicClusters(110) << endl;
315     cout << "Nb of BC120"= << rootEvent->nBasicClusters(120) << endl;
316     cout << "Nb of BC210=" << rootEvent->nBasicClusters(210) << endl; );
317     if(verbosity>2) cout << "electron et="<< (*recoElectrons)[j].pt() <
318     cout << "Nb of SC110=" << rootEvent->nSuperClusters(110) << endl;
319     cout << "Nb of SC120=" << rootEvent->nSuperClusters(120) << endl;
320     cout << "Nb of SC121=" << rootEvent->nSuperClusters(121) << endl;Misc
321     cout << "Nb of SC122=" << rootEvent->nSuperClusters(122) << endl;
322     cout << "Nb of SC210=" << rootEvent->nSuperClusters(210) << endl;
323     cout << "Nb of SC211=" << rootEvent->nSuperClusters(211) << endl;
324     */
325    
326    
327    
328     if(doCluster)
329     {
330     ClusterAssociator* myClusterAssociator = new ClusterAssociator();
331     myClusterAssociator->SetVerbosity(verbosity);
332     myClusterAssociator->Process(superClusters, clusters);
333     //myClusterAssociator->PrintSuperClusters(superClusters, clusters,0); // 0 to print all types SC
334     delete myClusterAssociator;
335     }
336    
337    
338     if(doPhoton && doCluster)
339     {
340     PhotonAssociator* myPhotonAssociator = new PhotonAssociator();
341     myPhotonAssociator->SetVerbosity(verbosity);
342 mlethuil 1.2 myPhotonAssociator->AssociateSuperCluster(correctedPhotons, superClusters);
343     //myPhotonAssociator->AssociateSuperCluster(uncorrectedPhotons, superClusters);
344 mlethuil 1.3 if(doPhotonConversion) myPhotonAssociator->AssociateConversion(iEvent, correctedPhotons, superClusters, conversionTracks);
345     //if(doPhotonConversion) myPhotonAssociator->AssociateConversion(iEvent, uncorrectedPhotons, superClusters, conversionTracks);
346 mlethuil 1.1 //myPhotonAssociator->PrintPhotons(correctedPhotons,superClusters,0); // 0 to print all types associated SC
347     myPhotonAssociator->FullPrintPhotons(correctedPhotons,superClusters,clusters,0); // 0 to print all types associated SC
348     delete myPhotonAssociator;
349     }
350    
351    
352 mlethuil 1.4 if(doPhoton && doPhotonIsolation) //FIXME - Test HCAL rechits availability
353 mlethuil 1.1 {
354     if(verbosity>1) cout << "PhotonIsolator......" << endl;
355    
356     // Get HB/HE rechits
357     edm::Handle<HBHERecHitCollection> hbheHandle;
358     iEvent.getByLabel("hbhereco", hbheHandle);
359     const HBHERecHitCollection* hbheHits = hbheHandle.product();
360    
361     // Get HO rechits
362     edm::Handle<HORecHitCollection> hoHandle;
363     iEvent.getByLabel("horeco", hoHandle);
364     const HORecHitCollection* hoHits = hoHandle.product();
365    
366     // Get HF rechits
367     edm::Handle<HFRecHitCollection> hfHandle;
368     iEvent.getByLabel("hfreco", hfHandle);
369     const HFRecHitCollection* hfHits = hfHandle.product();
370    
371     // Get Calo Geometry
372     edm::ESHandle<CaloGeometry> geomHandle;
373     (iSetup.get<IdealGeometryRecord>()).get(geomHandle);
374     const CaloGeometry* caloGeom = geomHandle.product();
375    
376     Float_t ecalIslandIsolation;
377     Float_t ecalDoubleConeIsolation;
378     Float_t hcalRecHitIsolation;
379 mlethuil 1.4 Float_t isoTracks;
380     Int_t isoNTracks;
381     TRootPhoton* localPhoton;
382    
383 mlethuil 1.1 PhotonIsolator* myPhotonIsolator = new PhotonIsolator(&myConfig);
384     for (int iphoton=0; iphoton<correctedPhotons->GetEntriesFast(); iphoton++)
385     {
386 mlethuil 1.4 localPhoton = (TRootPhoton*)correctedPhotons->At(iphoton);
387    
388 mlethuil 1.1 ecalIslandIsolation=-1.;
389 mlethuil 1.4 if (doCluster) ecalIslandIsolation = myPhotonIsolator->EcalIslandIsolation(localPhoton, superClusters, clusters);
390    
391 mlethuil 1.1 ecalDoubleConeIsolation=-1.;
392 mlethuil 1.4 if (doCluster) ecalDoubleConeIsolation = myPhotonIsolator->EcalDoubleConeIsolation(localPhoton, superClusters, clusters);
393    
394     hcalRecHitIsolation=-1.;
395 mlethuil 1.1 hcalRecHitIsolation = myPhotonIsolator->HcalRecHitIsolation(localPhoton, hbheHits, hoHits, hfHits, caloGeom);
396 mlethuil 1.4
397     isoTracks=-1.;
398     isoNTracks=-1;
399     if(doTrack) isoTracks = myPhotonIsolator->TrackerIsolation(localPhoton, tracks, isoNTracks);
400    
401     localPhoton->setIsolation(ecalIslandIsolation, ecalDoubleConeIsolation, hcalRecHitIsolation, isoTracks, isoNTracks);
402     if(verbosity>2) cout << "PhotonIsolator for [" << iphoton << "] " << *localPhoton << endl;
403 mlethuil 1.1 }
404     delete myPhotonIsolator;
405     }
406    
407    
408    
409     if(verbosity>2) cout << "Filling rootuple..." << endl;
410     rootTree_->Fill();
411     if(verbosity>2) cout << "Deleting objects..." << endl;
412     delete rootEvent;
413     if(doSignalMC) delete rootMCSignalEvent;
414     if(doMC) (*mcParticles).Delete();
415     if(doTrack) (*tracks).Delete();
416     if(doJet) (*jets).Delete();
417     if(doMuon) (*muons).Delete();
418     if(doElectron) (*electrons).Delete();
419     if(doPhoton) (*uncorrectedPhotons).Delete();
420     if(doPhoton) (*correctedPhotons).Delete();
421     if(doCluster) (*clusters).Delete();
422     if(doCluster) (*superClusters).Delete();
423 mlethuil 1.3 if(doPhotonConversion) (*conversionTracks).Delete();
424 mlethuil 1.1
425     }
426    
427    
428    
429     //define this as a plug-in
430     //DEFINE_ANOTHER_FWK_MODULE(TotoAnalyzer);