ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/Morgan/src/TotoAnalyzer.cc
Revision: 1.9
Committed: Wed Sep 3 15:02:32 2008 UTC (16 years, 7 months ago) by mlethuil
Content type: text/plain
Branch: MAIN
Changes since 1.8: +16 -4 lines
Log Message:
Bug correction to use Ted Likelihood for converted photon

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