ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/Morgan/src/TotoAnalyzer.cc
Revision: 1.11
Committed: Thu Oct 30 16:25:34 2008 UTC (16 years, 6 months ago) by lethuill
Content type: text/plain
Branch: MAIN
Changes since 1.10: +130 -150 lines
Log Message:
Updated for 2.1.X

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 lethuill 1.11 myConfig_ = iConfig.getParameter<ParameterSet>("myConfig");
11     producersNames_ = iConfig.getParameter<ParameterSet>("producersNames");
12 mlethuil 1.1 }
13    
14    
15     TotoAnalyzer::~TotoAnalyzer()
16 lethuill 1.11 {
17     }
18 mlethuil 1.1
19    
20    
21     // ------------ method called once each job just before starting event loop ------------
22     void TotoAnalyzer::beginJob(const edm::EventSetup&)
23     {
24    
25     // Load Config parameters
26 lethuill 1.11 verbosity = myConfig_.getUntrackedParameter<int>("verbosity", 0);
27     rootFileName_ = myConfig_.getUntrackedParameter<string>("RootFileName","noname.root");
28     isCSA07Soup = myConfig_.getUntrackedParameter<bool>("isCSA07Soup",false);
29     doHLT = myConfig_.getUntrackedParameter<bool>("doHLT",false);
30     doMC = myConfig_.getUntrackedParameter<bool>("doMC",false);
31     doSignalMC = myConfig_.getUntrackedParameter<bool>("doSignalMC",false);
32     doTrack = myConfig_.getUntrackedParameter<bool>("doTrack",false);
33     doJet = myConfig_.getUntrackedParameter<bool>("doJet",false);
34     doMuon = myConfig_.getUntrackedParameter<bool>("doMuon",false);
35     doElectron = myConfig_.getUntrackedParameter<bool>("doElectron",false);
36     doPhoton = myConfig_.getUntrackedParameter<bool>("doPhoton",false);
37     doCluster = myConfig_.getUntrackedParameter<bool>("doCluster",false);
38     doPhotonIsolation = myConfig_.getUntrackedParameter<bool>("doPhotonIsolation",false);
39     doPhotonConversion = myConfig_.getUntrackedParameter<bool>("doPhotonConversion",false);
40     doPhotonConversionMC = myConfig_.getUntrackedParameter<bool>("doPhotonConversionMC",false);
41 mlethuil 1.7
42 mlethuil 1.1 nTotEvt_ = 0;
43    
44     // initialize root output file
45     if(verbosity>0) cout << "New RootFile " << rootFileName_.c_str() << " is created" << endl;
46     rootFile_ = new TFile(rootFileName_.c_str(), "recreate");
47     rootFile_->cd();
48 mlethuil 1.5
49     runInfos_ = new TRootRun();
50     runTree_ = new TTree("runTree", "Global Run Infos");
51     runTree_->Branch ("runInfos", "TRootRun", &runInfos_);
52    
53 mlethuil 1.1 rootEvent = 0;
54 mlethuil 1.5 eventTree_ = new TTree("eventTree", "Event Infos");
55     eventTree_->Branch ("Event", "TRootEvent", &rootEvent);
56    
57     if(doHLT)
58     {
59 lethuill 1.11 hltAnalyzer_ = new HLTAnalyzer(&myConfig_);
60 mlethuil 1.5 hltAnalyzer_->setVerbosity(verbosity);
61     }
62    
63 mlethuil 1.1 if(doSignalMC)
64     {
65     if(verbosity>0) cout << "Signal generator info will be added to rootuple" << endl;
66     rootMCSignalEvent = 0;
67     cout << "Create MuMuGamma Branch" << endl;
68 mlethuil 1.5 eventTree_->Branch ("MuMuGamma", "TRootSignalEvent", &rootMCSignalEvent);
69 mlethuil 1.1 }
70    
71     if(doMC)
72     {
73     if(verbosity>0) cout << "MCParticles info will be added to rootuple" << endl;
74     mcParticles = new TClonesArray("TRootParticle", 1000);
75 mlethuil 1.5 eventTree_->Branch ("MCParticles", "TClonesArray", &mcParticles);
76 mlethuil 1.1 }
77    
78 mlethuil 1.7 if(doPhotonConversionMC)
79     {
80     if(verbosity>0) cout << "Converted MC Photons info will be added to rootuple" << endl;
81     mcPhotons = new TClonesArray("TRootMCPhoton", 1000);
82     eventTree_->Branch ("MCPhotons", "TClonesArray", &mcPhotons);
83     }
84    
85 mlethuil 1.1 if(doTrack)
86     {
87     if(verbosity>0) cout << "Tracks info will be added to rootuple" << endl;
88 mlethuil 1.4 tracks = new TClonesArray("TRootTrack", 1000);
89 mlethuil 1.5 eventTree_->Branch ("Tracks", "TClonesArray", &tracks);
90 mlethuil 1.1 }
91    
92     if(doJet)
93     {
94     if(verbosity>0) cout << "Jets info will be added to rootuple" << endl;
95     jets = new TClonesArray("TRootJet", 1000);
96 mlethuil 1.5 eventTree_->Branch ("Jets", "TClonesArray", &jets);
97 mlethuil 1.1 }
98    
99     if(doMuon)
100     {
101     if(verbosity>0) cout << "Muons info will be added to rootuple" << endl;
102     muons = new TClonesArray("TRootMuon", 1000);
103 mlethuil 1.5 eventTree_->Branch ("Muons", "TClonesArray", &muons);
104 mlethuil 1.1 }
105    
106     if(doElectron)
107     {
108     if(verbosity>0) cout << "Electrons info will be added to rootuple" << endl;
109     electrons = new TClonesArray("TRootElectron", 1000);
110 mlethuil 1.5 eventTree_->Branch ("Electrons", "TClonesArray", &electrons);
111 mlethuil 1.1 }
112    
113     if(doPhoton)
114     {
115     if(verbosity>0) cout << "Photons info will be added to rootuple" << endl;
116 lethuill 1.11 photons = new TClonesArray("TRootPhoton", 1000);
117     eventTree_->Branch ("Photons", "TClonesArray", &photons);
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 lethuill 1.11 std::string weightsString = myConfig_.getUntrackedParameter<string>("conversionLikelihoodWeightsFile","RecoEgamma/EgammaTools/data/TMVAnalysis_Likelihood.weights.txt");
137 mlethuil 1.9 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.9 if(doPhotonConversion) delete conversionLikelihoodCalculator_;
164 mlethuil 1.1 }
165    
166    
167    
168     // ------------ method called to for each event ------------
169     void TotoAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
170     {
171    
172     rootFile_->cd();
173     nTotEvt_++;
174     if( (verbosity>1) || (verbosity>0 && nTotEvt_%10==0 && nTotEvt_<=100) || (verbosity>0 && nTotEvt_%100==0 && nTotEvt_>100) )
175 lethuill 1.10 cout << endl << endl
176     << "####### TotoAnalyzer - Cumulated Events " << nTotEvt_
177     << " - Run " << iEvent.id().run()
178     << " - Event " << iEvent.id().event()
179     << " #######" << endl;
180    
181 mlethuil 1.1 // Global Event Infos
182     rootEvent = new TRootEvent();
183     rootEvent->setNb(nTotEvt_);
184    
185 lethuill 1.11 /*
186 mlethuil 1.6 // CSA07 Process Id and Event Weight
187     if(isCSA07Soup)
188     {
189     edm::Handle< double> weightHandle;
190     iEvent.getByLabel ("csaweightproducer","weight", weightHandle);
191     double weight = * weightHandle;
192     int processId = csa07::csa07ProcessId(iEvent, 1000., "csaweightproducer");
193     char * pname= csa07::csa07ProcessName(processId);
194     if(verbosity>1) cout << "CSA07 Soup - id=" << processId << " - "<< pname << " - weight=" << weight << endl;
195     rootEvent->setCsa07id(processId);
196     rootEvent->setCsa07weight(weight);
197     rootEvent->setCsa07process(pname);
198     }
199 mlethuil 1.5
200 mlethuil 1.6
201 mlethuil 1.5 // Trigger
202 lethuill 1.10 //rootEvent->setGlobalHLT(true);
203 mlethuil 1.5 if(doHLT)
204     {
205     cout << "Get TriggerResults..." << endl;
206     if (nTotEvt_==1) hltAnalyzer_->init(iEvent, rootEvent);
207     hltAnalyzer_->process(iEvent, rootEvent);
208     }
209 lethuill 1.10 //if ( ! rootEvent->passGlobalHLT() ) return;
210 mlethuil 1.5
211    
212 mlethuil 1.1 // MC Info
213     if(doMC)
214     {
215     rootMCSignalEvent = new TRootSignalEvent();
216     if(verbosity>1) cout << "Analysing MC info..." << endl;
217 lethuill 1.11 myMCAnalyzer = new MCAnalyzer(&myConfig_);
218 mlethuil 1.1 //myMCAnalyzer->SetVerbosity(1);
219 mlethuil 1.7 myMCAnalyzer->Process(iEvent, iSetup, rootEvent, rootMCSignalEvent, mcParticles, mcPhotons);
220 mlethuil 1.1 delete myMCAnalyzer;
221     }
222 lethuill 1.11 */
223 mlethuil 1.1
224     // Get Primary Vertices
225 lethuill 1.11 if(verbosity>1) cout << endl << "Analysing primary vertices collection..." << endl;
226     std::string vertexProducer = producersNames_.getParameter<string>("primaryVertexProducer");
227 mlethuil 1.1 edm::Handle< reco::VertexCollection > recoVertex;
228 lethuill 1.11 iEvent.getByLabel(vertexProducer, recoVertex);
229     if(verbosity>1) std::cout << " Producer: " << vertexProducer << " - Number of primary vertices = " << recoVertex->size() << std::endl;
230 mlethuil 1.1 for (unsigned int j=0; j<recoVertex->size(); j++)
231     {
232 lethuill 1.11 if(verbosity>2) cout << " ["<< setw(3) << j << "] Vertex - (vx,vy,vz)=(" << setw(12) << (*recoVertex)[j].x() << ", " << setw(12) << (*recoVertex)[j].y() << ", " << setw(12) << (*recoVertex)[j].z() << ")" << endl;
233 mlethuil 1.1 rootEvent->addPrimaryVertex( TVector3( (*recoVertex)[j].x(), (*recoVertex)[j].y(), (*recoVertex)[j].z() ) );
234     }
235    
236 lethuill 1.11
237 mlethuil 1.1 // Tracks
238     if(doTrack)
239     {
240 lethuill 1.11 if(verbosity>1) cout << endl << "Analysing tracks collection..." << endl;
241     std::string trackProducer = producersNames_.getParameter<string>("trackProducer");
242 mlethuil 1.1 edm::Handle<reco::TrackCollection> recoTracks;
243 lethuill 1.11 iEvent.getByLabel(trackProducer, recoTracks);
244     if(verbosity>1) std::cout << " Producer: " << trackProducer << " - Number of tracks = " << recoTracks->size() << std::endl;
245 mlethuil 1.4
246 mlethuil 1.1 for (unsigned int j=0; j<recoTracks->size(); j++)
247     {
248 mlethuil 1.4 const reco::HitPattern& hit = (*recoTracks)[j].hitPattern();
249    
250     TRootTrack localTrack(
251     (*recoTracks)[j].px()
252     ,(*recoTracks)[j].py()
253     ,(*recoTracks)[j].pz()
254     ,(*recoTracks)[j].p()
255     ,(*recoTracks)[j].vx()
256     ,(*recoTracks)[j].vy()
257     ,(*recoTracks)[j].vz()
258     ,0
259     ,(*recoTracks)[j].charge()
260     ,hit.numberOfValidPixelHits()
261     ,hit.numberOfValidTrackerHits()
262     ,(*recoTracks)[j].chi2()
263     ,(*recoTracks)[j].d0()
264     ,(*recoTracks)[j].d0Error()
265     ,(*recoTracks)[j].dz()
266     ,(*recoTracks)[j].dzError()
267     );
268    
269     new( (*tracks)[j] ) TRootTrack(localTrack);
270 lethuill 1.11 if(verbosity>2) cout << " ["<< setw(3) << j << "] " << localTrack << endl;
271 mlethuil 1.4
272 mlethuil 1.1 }
273     }
274    
275 lethuill 1.11
276 mlethuil 1.1 // Calo Jet
277     if(doJet)
278     {
279 lethuill 1.11 if(verbosity>1) cout << endl << "Analysing jets collection..." << endl;
280     std::string jetProducer = producersNames_.getParameter<string>("jetProducer");
281 mlethuil 1.1 edm::Handle< reco::CaloJetCollection > recoJets;
282 lethuill 1.11 iEvent.getByLabel(jetProducer, recoJets);
283     if(verbosity>1) std::cout << " Producer: " << jetProducer << " - Number of jets = " << recoJets->size() << std::endl;
284 mlethuil 1.1 for (unsigned int j=0; j<recoJets->size(); j++)
285     {
286 lethuill 1.11 TRootJet localJet(
287     (*recoJets)[j].px()
288     ,(*recoJets)[j].py()
289     ,(*recoJets)[j].pz()
290     ,(*recoJets)[j].energy()
291     );
292     new( (*jets)[j] ) TRootJet(localJet);
293     if(verbosity>2) cout << " ["<< setw(3) << j << "] " << localJet << endl;
294 mlethuil 1.1 }
295     }
296    
297    
298     // Muons
299     if(doMuon)
300     {
301     Float_t sintheta = 0.;
302 lethuill 1.11 if(verbosity>1) cout << endl << "Analysing muons collection..." << endl;
303     std::string muonProducer = producersNames_.getParameter<string>("muonProducer");
304 mlethuil 1.1 edm::Handle < reco::MuonCollection > recoMuons;
305 lethuill 1.11 iEvent.getByLabel(muonProducer, recoMuons);
306     if(verbosity>1) std::cout << " Producer: " << muonProducer << " - Number of muons = " << recoMuons->size() << std::endl;
307 mlethuil 1.1 for (unsigned int j=0; j<recoMuons->size(); j++)
308     {
309     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() );
310     sintheta = sin( localMuon.Theta() );
311 lethuill 1.11 localMuon.setCaloEnergy( (*recoMuons)[j].calEnergy().em * sintheta, (*recoMuons)[j].calEnergy().emS9 * sintheta,
312     (*recoMuons)[j].calEnergy().had * sintheta, (*recoMuons)[j].calEnergy().hadS9 * sintheta,
313     (*recoMuons)[j].calEnergy().ho * sintheta, (*recoMuons)[j].calEnergy().hoS9 * sintheta, (*recoMuons)[j].caloCompatibility() );
314     localMuon.setIsoR03( (*recoMuons)[j].isolationR03().emEt, (*recoMuons)[j].isolationR03().hadEt, (*recoMuons)[j].isolationR03().hoEt,
315     (*recoMuons)[j].isolationR03().sumPt, (*recoMuons)[j].isolationR03().nTracks, (*recoMuons)[j].isolationR03().nJets);
316     localMuon.setIsoR05( (*recoMuons)[j].isolationR05().emEt, (*recoMuons)[j].isolationR05().hadEt, (*recoMuons)[j].isolationR05().hoEt,
317     (*recoMuons)[j].isolationR05().sumPt, (*recoMuons)[j].isolationR05().nTracks, (*recoMuons)[j].isolationR05().nJets);
318 mlethuil 1.1 new( (*muons)[j] ) TRootMuon(localMuon);
319 lethuill 1.11 if(verbosity>2) cout << " ["<< setw(3) << j << "] " << localMuon << endl;
320 mlethuil 1.1 }
321     }
322    
323    
324     // Electrons
325     if(doElectron)
326     {
327 lethuill 1.11 if(verbosity>1) cout << endl << "Analysing electrons collection..." << endl;
328     std::string electronProducer = producersNames_.getParameter<string>("electronProducer");
329     edm::Handle< reco::GsfElectronCollection > recoElectrons;
330     iEvent.getByLabel(electronProducer, recoElectrons);
331     if(verbosity>1) std::cout << " Producer: " << electronProducer << " - Number of electrons = " << recoElectrons->size() << std::endl;
332 mlethuil 1.1 for (unsigned int j=0; j<recoElectrons->size(); j++)
333     {
334     float elec_e = ((*recoElectrons)[j].superCluster())->energy();
335     float elec_px = (*recoElectrons)[j].px() * elec_e / (*recoElectrons)[j].p();
336     float elec_py = (*recoElectrons)[j].py() * elec_e / (*recoElectrons)[j].p();
337     float elec_pz = (*recoElectrons)[j].pz() * elec_e / (*recoElectrons)[j].p();
338     TRootElectron localElectron( elec_px, elec_py, elec_pz, elec_e, (*recoElectrons)[j].vx(), (*recoElectrons)[j].vy(), (*recoElectrons)[j].vz(), 11, (*recoElectrons)[j].charge() );
339     new( (*electrons)[j] ) TRootElectron(localElectron);
340 lethuill 1.11 if(verbosity>2) cout << " ["<< setw(3) << j << "] " << localElectron << endl;
341 mlethuil 1.1 }
342     }
343    
344 lethuill 1.11
345     // Lazy Tools to calculate Cluster shape variables
346     if(verbosity>1) cout << endl << "Loading egamma LazyTools..." << endl;
347     edm::InputTag reducedBarrelEcalRecHitCollection_ = producersNames_.getParameter<edm::InputTag>("reducedBarrelEcalRecHitCollection");
348     edm::InputTag reducedEndcapEcalRecHitCollection_ = producersNames_.getParameter<edm::InputTag>("reducedEndcapEcalRecHitCollection");
349     EcalClusterLazyTools lazyTools( iEvent, iSetup, reducedBarrelEcalRecHitCollection_, reducedEndcapEcalRecHitCollection_ );
350    
351    
352 mlethuil 1.1 // Photons
353     if(doPhoton)
354     {
355 lethuill 1.11 if(verbosity>1) cout << endl << "Analysing photons collection..." << endl;
356     bool doPhotonVertexCorrection = false;
357     myPhotonAnalyzer = new PhotonAnalyzer(producersNames_, verbosity, doPhotonConversion, doPhotonVertexCorrection);
358     std::string photonProducer = producersNames_.getParameter<string>("photonProducer");
359     myPhotonAnalyzer->Process(iEvent, iSetup, rootEvent, photons, photonProducer, conversionTracks, conversionLikelihoodCalculator_, lazyTools);
360 mlethuil 1.1 delete myPhotonAnalyzer;
361     }
362 lethuill 1.11
363    
364 mlethuil 1.1 // Ecal Clusters
365     if(doCluster)
366     {
367 lethuill 1.11 if(verbosity>1) cout << endl << "Analysing BasicClusters collection..." << endl;
368 mlethuil 1.1 myClusterAnalyzer = new ClusterAnalyzer();
369     myClusterAnalyzer->SetVerbosity(verbosity);
370 lethuill 1.11 myClusterAnalyzer->Process(iEvent, rootEvent, lazyTools, clusters, "hybridSuperClusters","hybridBarrelBasicClusters", 110);
371     myClusterAnalyzer->Process(iEvent, rootEvent, lazyTools, clusters, "multi5x5BasicClusters", "multi5x5EndcapBasicClusters", 320);
372 mlethuil 1.1 delete myClusterAnalyzer;
373 lethuill 1.11
374     if(verbosity>1) cout << endl << "Analysing SuperClusters collection..." << endl;
375 mlethuil 1.1 mySClusterAnalyzer = new SuperClusterAnalyzer();
376     mySClusterAnalyzer->SetVerbosity(verbosity);
377     mySClusterAnalyzer->Process(iEvent, rootEvent, superClusters, "hybridSuperClusters","",210);
378     mySClusterAnalyzer->Process(iEvent, rootEvent, superClusters, "correctedHybridSuperClusters","",211);
379 lethuill 1.11 mySClusterAnalyzer->Process(iEvent, rootEvent, superClusters, "multi5x5SuperClusters","multi5x5EndcapSuperClusters",320);
380     mySClusterAnalyzer->Process(iEvent, rootEvent, superClusters, "multi5x5SuperClustersWithPreshower","",323);
381     mySClusterAnalyzer->Process(iEvent, rootEvent, superClusters, "correctedMulti5x5SuperClustersWithPreshower","",322);
382 mlethuil 1.1 delete mySClusterAnalyzer;
383     }
384    
385    
386 lethuill 1.11 //cout << "Nb of BC110=" << rootEvent->nBasicClusters(110) << endl;
387     //cout << "Nb of BC120"= << rootEvent->nBasicClusters(120) << endl;
388     //cout << "Nb of BC210=" << rootEvent->nBasicClusters(210) << endl; );
389     //if(verbosity>2) cout << "electron et="<< (*recoElectrons)[j].pt() <
390     //cout << "Nb of SC110=" << rootEvent->nSuperClusters(110) << endl;
391     //cout << "Nb of SC120=" << rootEvent->nSuperClusters(120) << endl;
392     //cout << "Nb of SC121=" << rootEvent->nSuperClusters(121) << endl;Misc
393     //cout << "Nb of SC122=" << rootEvent->nSuperClusters(122) << endl;
394     //cout << "Nb of SC210=" << rootEvent->nSuperClusters(210) << endl;
395     //cout << "Nb of SC211=" << rootEvent->nSuperClusters(211) << endl;
396 mlethuil 1.1
397    
398     if(doCluster)
399     {
400     ClusterAssociator* myClusterAssociator = new ClusterAssociator();
401     myClusterAssociator->SetVerbosity(verbosity);
402     myClusterAssociator->Process(superClusters, clusters);
403     //myClusterAssociator->PrintSuperClusters(superClusters, clusters,0); // 0 to print all types SC
404     delete myClusterAssociator;
405     }
406 lethuill 1.11
407    
408 mlethuil 1.1 if(doPhoton && doCluster)
409     {
410     PhotonAssociator* myPhotonAssociator = new PhotonAssociator();
411     myPhotonAssociator->SetVerbosity(verbosity);
412 lethuill 1.11 myPhotonAssociator->AssociateSuperCluster(photons, superClusters);
413     // Obsolete... converted photon now contained in reco::Photon object. Done in PhotonAnalyzer
414     //if(doPhotonConversion) myPhotonAssociator->AssociateConversion(iEvent, photons, superClusters, conversionTracks, conversionLikelihoodCalculator_);
415 mlethuil 1.1 delete myPhotonAssociator;
416     }
417    
418    
419 lethuill 1.11 if(doPhoton && doPhotonIsolation)
420 mlethuil 1.1 {
421 lethuill 1.11 if(verbosity>1) cout << endl << "Calculating Photon isolation..." << endl;
422 mlethuil 1.1
423     Float_t ecalIslandIsolation;
424     Float_t ecalDoubleConeIsolation;
425     Float_t hcalRecHitIsolation;
426 mlethuil 1.4 Float_t isoTracks;
427     Int_t isoNTracks;
428     TRootPhoton* localPhoton;
429    
430 lethuill 1.11 PhotonIsolator* myPhotonIsolator = new PhotonIsolator(&myConfig_, &producersNames_);
431     for (int iphoton=0; iphoton<photons->GetEntriesFast(); iphoton++)
432 mlethuil 1.1 {
433 lethuill 1.11 localPhoton = (TRootPhoton*)photons->At(iphoton);
434 mlethuil 1.4
435 mlethuil 1.1 ecalIslandIsolation=-1.;
436 mlethuil 1.4 if (doCluster) ecalIslandIsolation = myPhotonIsolator->EcalIslandIsolation(localPhoton, superClusters, clusters);
437    
438 mlethuil 1.1 ecalDoubleConeIsolation=-1.;
439 mlethuil 1.4 if (doCluster) ecalDoubleConeIsolation = myPhotonIsolator->EcalDoubleConeIsolation(localPhoton, superClusters, clusters);
440    
441     hcalRecHitIsolation=-1.;
442 lethuill 1.11 //FIXME - Test HCAL rechits availability - LoadHcalRecHits() returning bool
443     myPhotonIsolator->LoadHcalRecHits(iEvent, iSetup);
444     hcalRecHitIsolation = myPhotonIsolator->HcalRecHitIsolation(localPhoton);
445 mlethuil 1.4
446     isoTracks=-1.;
447     isoNTracks=-1;
448     if(doTrack) isoTracks = myPhotonIsolator->TrackerIsolation(localPhoton, tracks, isoNTracks);
449    
450     localPhoton->setIsolation(ecalIslandIsolation, ecalDoubleConeIsolation, hcalRecHitIsolation, isoTracks, isoNTracks);
451 lethuill 1.11 if(verbosity>4) cout << " After isolation - ["<< setw(3) << iphoton << "] " << *localPhoton << endl;
452 mlethuil 1.1 }
453     delete myPhotonIsolator;
454     }
455    
456 lethuill 1.11 // Print Photons / SuperClusters / BasicClusters arborescence
457     if(doPhoton)
458     {
459     PhotonAssociator* myPhotonAssociator = new PhotonAssociator();
460     if(verbosity>2) myPhotonAssociator->FullPrintPhotons(photons,superClusters,clusters,0);
461     delete myPhotonAssociator;
462     }
463 mlethuil 1.1
464    
465 lethuill 1.11 if(verbosity>1) cout << endl << "Filling rootuple..." << endl;
466 mlethuil 1.5 eventTree_->Fill();
467 lethuill 1.11 if(verbosity>1) cout << endl << "Deleting objects..." << endl;
468 mlethuil 1.1 delete rootEvent;
469 lethuill 1.11 //if(doSignalMC) delete rootMCSignalEvent;
470     //if(doMC) (*mcParticles).Delete();
471     //if(doPhotonConversionMC) (*mcPhotons).Delete();
472 mlethuil 1.1 if(doTrack) (*tracks).Delete();
473     if(doJet) (*jets).Delete();
474     if(doMuon) (*muons).Delete();
475     if(doElectron) (*electrons).Delete();
476 lethuill 1.11 if(doPhoton) (*photons).Delete();
477 mlethuil 1.1 if(doCluster) (*clusters).Delete();
478     if(doCluster) (*superClusters).Delete();
479 lethuill 1.11 //if(doPhotonConversion) (*conversionTracks).Delete();
480     if(verbosity>0) cout << endl;
481 mlethuil 1.1 }
482    
483    
484    
485     //define this as a plug-in
486     //DEFINE_ANOTHER_FWK_MODULE(TotoAnalyzer);