ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/Morgan/src/TotoAnalyzer.cc
Revision: 1.1
Committed: Mon May 19 16:12:29 2008 UTC (16 years, 11 months ago) by mlethuil
Content type: text/plain
Branch: MAIN
Log Message:
*** empty log message ***

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