ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/Morgan/src/PhotonAnalyzer.cc
Revision: 1.15
Committed: Fri Sep 18 14:14:21 2009 UTC (15 years, 7 months ago) by lethuill
Content type: text/plain
Branch: MAIN
CVS Tags: all_3_3_2_01, all_3_2_5_02, all_3_2_5_01, HEAD
Changes since 1.14: +77 -75 lines
Log Message:
Update for 3.2.X

File Contents

# User Rev Content
1 lethuill 1.5 #include "../interface/PhotonAnalyzer.h"
2 mlethuil 1.1
3     using namespace std;
4     using namespace reco;
5     using namespace edm;
6    
7 lethuill 1.2
8 lethuill 1.10 PhotonAnalyzer::PhotonAnalyzer(const edm::ParameterSet& producersNames, const edm::ParameterSet& myConfig, int verbosity):verbosity_(verbosity), iconvtrack_(0)
9 lethuill 1.2 {
10 lethuill 1.14 dataType_ = producersNames.getUntrackedParameter<string>("dataType","unknown");
11     photonProducer_ = producersNames.getParameter<edm::InputTag>("photonProducer");
12     doPhotonConversion_ = myConfig.getUntrackedParameter<bool>("doPhotonConversion", false);
13     doVertexCorrection_ = myConfig.getUntrackedParameter<bool>("doPhotonVertexCorrection", false);
14     useMC_ = myConfig.getUntrackedParameter<bool>("doMuonMC", false);
15     allowMissingCollection_ = producersNames.getUntrackedParameter<bool>("allowMissingCollection", false);
16 lethuill 1.2 }
17    
18 lethuill 1.14
19 mlethuil 1.1 PhotonAnalyzer::~PhotonAnalyzer()
20     {
21     }
22    
23 lethuill 1.14
24     bool PhotonAnalyzer::process(const edm::Event& iEvent, const edm::EventSetup& iSetup, TRootEvent* rootEvent, TClonesArray* rootPhotons, TClonesArray* conversionTracks, ConversionLikelihoodCalculator* convLikelihoodCalculator, EcalClusterLazyTools* lazyTools)
25 mlethuil 1.1 {
26 lethuill 1.14
27     unsigned int nPhotons=0;
28 mlethuil 1.1
29 lethuill 1.14 bool doPhotonID = true;
30     edm::Handle< reco::PhotonCollection > recoPhotons;
31 lethuill 1.15 const edm::ValueMap<Bool_t> *loosePhotonIdMap = 0;
32     const edm::ValueMap<Bool_t> *tightPhotonIdMap = 0;
33    
34     if( dataType_=="RECO" )
35 lethuill 1.14 {
36     try
37     {
38     iEvent.getByLabel(photonProducer_, recoPhotons);
39     nPhotons = recoPhotons->size();
40     }
41     catch (cms::Exception& exception)
42     {
43     if ( !allowMissingCollection_ )
44     {
45     cout << " ##### ERROR IN PhotonAnalyzer::process => reco::Photon collection is missing #####"<<endl;
46     throw exception;
47     }
48     if(verbosity_>1) cout << " ===> No reco::Photon collection, skip photons info" << endl;
49     return false;
50     }
51 lethuill 1.15
52 lethuill 1.14 try
53     {
54 lethuill 1.15 edm::Handle< edm::ValueMap<Bool_t> > loosePhotonId;
55     iEvent.getByLabel("PhotonIDProd", "PhotonCutBasedIDLoose", loosePhotonId);
56     loosePhotonIdMap = loosePhotonId.product();
57 lethuill 1.14 }
58     catch (cms::Exception& exception)
59     {
60     if ( !allowMissingCollection_ )
61     {
62 lethuill 1.15 cout << " ##### ERROR IN PhotonAnalyzer::process => PhotonCutBasedIDLoose is missing #####"<<endl;
63 lethuill 1.14 throw exception;
64     }
65 lethuill 1.15 if(verbosity_>1) cout << " ===> No PhotonCutBasedIDLoose collection, skip photonID info" << endl;
66     doPhotonID = false;
67     }
68    
69     try
70     {
71     edm::Handle< edm::ValueMap<Bool_t> > tightPhotonId;
72     iEvent.getByLabel("PhotonIDProd", "PhotonCutBasedIDTight", tightPhotonId);
73     tightPhotonIdMap = tightPhotonId.product();
74 lethuill 1.14 }
75 lethuill 1.15 catch (cms::Exception& exception)
76     {
77     if ( !allowMissingCollection_ )
78     {
79     cout << " ##### ERROR IN PhotonAnalyzer::process => PhotonCutBasedIDTight is missing #####"<<endl;
80     throw exception;
81     }
82     if(verbosity_>1) cout << " ===> No PhotonCutBasedIDTight collection, skip photonID info" << endl;
83     doPhotonID = false;
84     }
85     }
86 lethuill 1.14
87     edm::Handle < std::vector <pat::Photon> > patPhotons;
88     if( dataType_=="PAT" )
89     {
90     try
91     {
92     iEvent.getByLabel(photonProducer_, patPhotons);
93     nPhotons = patPhotons->size();
94     }
95     catch (cms::Exception& exception)
96     {
97     if ( !allowMissingCollection_ )
98     {
99     cout << " ##### ERROR IN PhotonAnalyzer::process => pat::Photon collection is missing #####"<<endl;
100     throw exception;
101     }
102     if(verbosity_>1) cout << " ===> No pat::Photon collection, skip photons info" << endl;
103     return false;
104     }
105     }
106    
107     if(verbosity_>1) std::cout << " Number of photons = " << nPhotons << " Label: " << photonProducer_.label() << " Instance: " << photonProducer_.instance() << std::endl;
108    
109    
110     // TODO - add Pi0Disc... not implemented in 2.X.X
111     //edm::Handle<reco::PhotonPi0DiscriminatorAssociationMap> pi0map;
112     //iEvent.getByLabel("piZeroDiscriminators","PhotonPi0DiscriminatorAssociationMap", pi0map);
113     //reco::PhotonPi0DiscriminatorAssociationMap::const_iterator pi0mapIter;
114     //Double_t pi0nn;
115    
116    
117     for (unsigned int j=0; j<nPhotons; j++)
118     {
119     const reco::Photon* photon = 0;
120     if( dataType_=="RECO" ) photon = &((*recoPhotons)[j]);
121     if( dataType_=="PAT" ) photon = (const reco::Photon*) ( & ((*patPhotons)[j]) );
122    
123     TRootPhoton localPhoton(
124     photon->px()
125     ,photon->py()
126     ,photon->pz()
127     ,photon->energy()
128     ,photon->vx()
129     ,photon->vy()
130     ,photon->vz()
131     ,photon->pdgId()
132     ,photon->charge()
133     );
134    
135    
136     // Variables from reco::Photon
137     localPhoton.setCaloPosition( photon->caloPosition().X(), photon->caloPosition().Y(), photon->caloPosition().Z() );
138     localPhoton.setHoE(photon->hadronicOverEm());
139     localPhoton.setHasPixelSeed( photon->hasPixelSeed() );
140    
141    
142     // Variables from reco::SuperCluster
143     reco::SuperClusterRef superCluster = photon->superCluster();
144     if ( superCluster.isNonnull() )
145     {
146     localPhoton.setNbClusters(superCluster->clustersSize());
147     localPhoton.setSuperClusterRawEnergy( superCluster->rawEnergy() );
148     localPhoton.setPreshowerEnergy(superCluster->preshowerEnergy());
149     }
150    
151    
152     // Cluster Shape variables
153 lethuill 1.15 // need reco::SuperCluster and reco::CaloCluster
154 lethuill 1.14 if ( superCluster.isNonnull() )
155     {
156 lethuill 1.15 reco::CaloClusterPtr seedCaloCluster = superCluster->seed();
157     if ( seedCaloCluster.isNonnull() ) localPhoton.setClusterAlgo(seedCaloCluster->algo());
158 lethuill 1.14
159     // dR of the cone centered on the reco::Photon and containing all its basic clusters constituents
160     bool atLeastOneBC = false;
161     Float_t caloConeSize = 0.;
162 lethuill 1.15 for (reco::CaloCluster_iterator caloCluster = (*superCluster).clustersBegin(); caloCluster != (*superCluster).clustersEnd(); ++caloCluster )
163 lethuill 1.14 {
164     atLeastOneBC = true;
165     // Take SC position which is not vertex corrected
166     TVector3 scPosition( superCluster->position().x(), superCluster->position().y(), superCluster->position().z() );
167 lethuill 1.15 Float_t dR = scPosition.DeltaR( TVector3( (*caloCluster)->position().x(), (*caloCluster)->position().y(), (*caloCluster)->position().z() ) );
168 lethuill 1.14 if (dR > caloConeSize) caloConeSize = dR;
169     }
170     if (! atLeastOneBC) caloConeSize = -999.;
171     localPhoton.setCaloConeSize(caloConeSize);
172    
173     // need reduced Ecal RecHits Collections for EcalClusterLazyTools
174 lethuill 1.15 if ( seedCaloCluster.isNonnull() && lazyTools != 0 )
175 lethuill 1.14 {
176 lethuill 1.15 localPhoton.setE2x2(lazyTools->e2x2(*seedCaloCluster));
177     localPhoton.setE3x3( lazyTools->e3x3(*seedCaloCluster) );
178     localPhoton.setE5x5( lazyTools->e5x5(*seedCaloCluster) );
179     localPhoton.setEMax( lazyTools->eMax(*seedCaloCluster) );
180 lethuill 1.14 }
181     }
182    
183    
184     // Photon conversions
185     // Variables from reco::Conversion
186     if( doPhotonConversion_ )
187     {
188     // If more than one associated convertedPhotons is found, then
189     // if nTracks=2, select the convertedPhoton with greatest Likelihood value using Ted's weights file
190     // if nTracks=1, select the convertedPhoton with E/p closest to 1
191     double eoverp;
192     double best_eoverp = 999999.;
193     double likely;
194     double best_likely = -1.;
195     int best_iconv = -1;
196     int best_iconv_likely = -1;
197     int best_iconv_eoverp = -1;
198    
199 lethuill 1.15 reco::ConversionRefVector conversions = photon->conversions();
200 lethuill 1.14 for (unsigned int iconv=0; iconv<conversions.size(); iconv++)
201     {
202     if(verbosity_>4) cout << " ["<< setw(3) << iconv << "] Conversion - "
203     << " conv_vertex (x,y,z)=(" << conversions[iconv]->conversionVertex().x()
204     << "," << conversions[iconv]->conversionVertex().y()
205     << "," << conversions[iconv]->conversionVertex().z()
206     << ") isConverted=" << conversions[iconv]->isConverted()
207     << " nTracks=" << conversions[iconv]->nTracks()
208     << " primary_vtx_z=" << conversions[iconv]->zOfPrimaryVertexFromTracks()
209     << " inv_mass=" << conversions[iconv]->pairInvariantMass()
210     << " E/p=" << conversions[iconv]->EoverP()
211     << " cotanTheta=" << conversions[iconv]->pairCotThetaSeparation()
212     << " likely=" << convLikelihoodCalculator->calculateLikelihood(conversions[iconv])
213     << endl;
214    
215     likely = convLikelihoodCalculator->calculateLikelihood(conversions[iconv]);
216     if ( likely>best_likely )
217     {
218     best_iconv_likely = iconv;
219     best_likely = likely;
220     }
221    
222     eoverp = conversions[iconv]->EoverP();
223     if ( abs(eoverp-1)<abs(best_eoverp-1) )
224     {
225     best_iconv_eoverp = iconv;
226     best_eoverp = eoverp;
227     }
228     }
229    
230     best_iconv = (best_iconv_likely==-1 ? best_iconv_eoverp : best_iconv_likely);
231    
232     // Update Photon object with conversion infos
233     if ( best_iconv != -1 )
234     {
235     if(verbosity_>4) cout
236     << " ===> TRootPhoton[" << j << "] associated to Conversion[" << best_iconv
237     << "] with E/p=" << conversions[best_iconv]->EoverP()
238     << " likely=" << convLikelihoodCalculator->calculateLikelihood(conversions[best_iconv])
239     << endl;
240    
241     localPhoton.setConvNTracks( conversions[best_iconv]->nTracks() );
242     localPhoton.setConvEoverP( conversions[best_iconv]->EoverP() );
243     localPhoton.setConvMass( conversions[best_iconv]->pairInvariantMass() );
244     localPhoton.setConvCotanTheta( conversions[best_iconv]->pairCotThetaSeparation() );
245     localPhoton.setConvLikely( convLikelihoodCalculator->calculateLikelihood(conversions[best_iconv]) );
246     localPhoton.setConvVertex( conversions[best_iconv]->conversionVertex().x(), conversions[best_iconv]->conversionVertex().y(), conversions[best_iconv]->conversionVertex().z() );
247     std::vector<math::XYZPoint> impactVector = conversions[best_iconv]->ecalImpactPosition();
248     if ( impactVector.size()>0 ) localPhoton.setConvEcalImpactPosition1( impactVector.at(0).x(), impactVector.at(0).y(), impactVector.at(0).z() );
249     if ( impactVector.size()>1 ) localPhoton.setConvEcalImpactPosition2( impactVector.at(1).x(), impactVector.at(1).y(), impactVector.at(1).z() );
250    
251     if ( conversionTracks!=0 ) // Check branch is activated
252     {
253     if ( conversions[best_iconv]->nTracks()>0 )
254     {
255     std::vector<reco::TrackRef> tracks = conversions[best_iconv]->tracks();
256     reco::TrackRef tk1 = tracks.at(0);
257     const reco::HitPattern& hit1 = tk1->hitPattern();
258     new( (*conversionTracks)[iconvtrack_] ) TRootTrack( tk1->px(), tk1->py(), tk1->pz(), tk1->p(), tk1->vx(), tk1->vy(), tk1->vz(), 0, tk1->charge()
259     ,hit1.pixelLayersWithMeasurement(), hit1.stripLayersWithMeasurement(), tk1->chi2(), tk1->d0(), tk1->d0Error(), tk1->dz(), tk1->dzError() );
260     localPhoton.setConvIndexTrack1(iconvtrack_);
261     localPhoton.setConvTrack1((*conversionTracks)[iconvtrack_]);
262     iconvtrack_++;
263    
264     if ( conversions[best_iconv]->nTracks()>1 )
265     {
266     reco::TrackRef tk2 = tracks.at(1);
267     const reco::HitPattern& hit2 = tk2->hitPattern();
268     new( (*conversionTracks)[iconvtrack_] ) TRootTrack( tk2->px(), tk2->py(), tk2->pz(), tk2->p(), tk2->vx(), tk2->vy(), tk2->vz(), 0, tk2->charge()
269     ,hit2.pixelLayersWithMeasurement(), hit2.stripLayersWithMeasurement(), tk2->chi2(), tk2->d0(), tk2->d0Error(), tk2->dz(), tk2->dzError() );
270     localPhoton.setConvIndexTrack2(iconvtrack_);
271     localPhoton.setConvTrack2((*conversionTracks)[iconvtrack_]);
272     iconvtrack_++;
273     }
274     }
275     }
276     }
277     }
278 lethuill 1.15
279    
280     // Photon isolation calculated by PhotonID
281     // now embeded in reco::Photon
282     localPhoton.setIsolation(
283     photon->ecalRecHitSumEtConeDR03()
284     ,photon->hcalTowerSumEtConeDR03()
285     ,photon->trkSumPtSolidConeDR03()
286     ,photon->trkSumPtHollowConeDR03()
287     ,photon->nTrkSolidConeDR03()
288     ,photon->nTrkHollowConeDR03()
289     );
290 lethuill 1.14
291    
292     if( dataType_=="RECO" )
293     {
294     // Some specific methods requiring RECO / AOD format
295     // Do association to genParticle ?
296    
297     if (doPhotonID)
298     {
299     // Photon ID
300 lethuill 1.15 // need corresponding PhotonID ValueMap
301     edm::Ref<reco::PhotonCollection> photonRef(recoPhotons,j);
302     Bool_t looseID = (*loosePhotonIdMap)[photonRef];
303     Bool_t tightID = (*tightPhotonIdMap)[photonRef];
304     localPhoton.setBitsID(
305     looseID
306     ,tightID
307     ,photon->isEB()
308     ,photon->isEE()
309     ,photon->isEBGap()
310     ,photon->isEEGap()
311     ,photon->isEBEEGap()
312 lethuill 1.14 );
313     }
314    
315     /*
316 lethuill 1.15 if(verbosity_>4) std::cout << "seed E5x5=" << lazyTools->e5x5( *seedCaloCluster )
317 lethuill 1.14 //<< " pi0NN=" << pi0nn
318     << " photonID->isolationEcalRecHit()=" << photonID->isolationEcalRecHit()
319     << " photonID->r9()=" << photonID->r9()
320 lethuill 1.15 << " e3x3 / SC->energy()=" << ( lazyTools->e3x3( *seedCaloCluster ) / photon->superCluster()->energy() )
321     << " e3x3 / SC->rawEnergy()=" << ( lazyTools->e3x3( *seedCaloCluster ) / photon->superCluster()->rawEnergy() )
322 lethuill 1.14 << " sc->pos=" << photon->superCluster()->position().X() << " , " << photon->superCluster()->position().Y() << " , " << photon->superCluster()->position().Z()
323     << " photon->caloPos=" << photon->caloPosition().X() << " , " << photon->caloPosition().Y() << " , " << photon->caloPosition().Z() << std::endl;
324     */
325     }
326    
327    
328     if( dataType_=="PAT" )
329     {
330     // Some specific methods to pat::Photon
331     const pat::Photon *patPhoton = dynamic_cast<const pat::Photon*>(&*photon);
332    
333     // Photon ID
334 lethuill 1.15 // Use the PhotonID pairs embeded in pat::Photon
335     Bool_t looseID = false;
336     Bool_t tightID = false;
337     if (patPhoton->isPhotonIDAvailable("PhotonCutBasedIDLoose")) looseID = patPhoton->photonID("PhotonCutBasedIDLoose");
338     if (patPhoton->isPhotonIDAvailable("PhotonCutBasedIDTight")) tightID = patPhoton->photonID("PhotonCutBasedIDTight");
339     localPhoton.setBitsID(
340     looseID
341     ,tightID
342     ,patPhoton->isEB()
343     ,patPhoton->isEE()
344 lethuill 1.14 ,patPhoton->isEBGap()
345     ,patPhoton->isEEGap()
346     ,patPhoton->isEBEEGap()
347     );
348    
349     if(useMC_)
350     {
351     // MC truth associator index
352     if ((patPhoton->genParticleRef()).isNonnull()) {
353     localPhoton.setGenParticleIndex((patPhoton->genParticleRef()).index());
354     } else {
355     localPhoton.setGenParticleIndex(-1);
356     }
357     }
358     }
359    
360    
361     // FIXME - Add pi0nn
362     //pi0mapIter = pi0map->find(edm::Ref<reco::PhotonCollection>(recoPhotons,iPhoton));
363     //pi0nn = ( pi0mapIter == pi0map->end() ? -1 : pi0mapIter->val );
364     //localPhoton.setPi0nn( pi0nn );
365    
366    
367     // Vertex correction to photon - Assume photon is coming from primary vertex
368     if( doVertexCorrection_)
369     {
370     if(verbosity_>3) cout << " Before vertex correction ["<< setw(3) << j << "] " << localPhoton << endl;
371     if( rootEvent->primaryVertex()!=0 )
372     {
373     TVector3 vertex( rootEvent->primaryVertex()->x(), rootEvent->primaryVertex()->y(), rootEvent->primaryVertex()->z() );
374     localPhoton.setVertex(vertex);
375     }
376     else
377     {
378     cout << " PhotonAnalyzer - NO SELECTED PRIMARY VERTEX FOUND !!!" << endl;
379     }
380     }
381    
382    
383     // Stock new TRootPhoton in TCloneArray
384     new( (*rootPhotons)[j] ) TRootPhoton(localPhoton);
385     if(verbosity_>2) cout << " ["<< setw(3) << j << "] " << localPhoton << endl;
386     }
387    
388     return true;
389 lethuill 1.10 }