ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/Morgan/src/TotoAnalyzer.cc
Revision: 1.3
Committed: Wed May 28 16:33:59 2008 UTC (16 years, 11 months ago) by mlethuil
Content type: text/plain
Branch: MAIN
Changes since 1.2: +12 -2 lines
Log Message:
Add conversion tracks collection

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