ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/Morgan/src/PhotonAnalyzer.cc
Revision: 1.14
Committed: Wed Jun 10 11:17:06 2009 UTC (15 years, 10 months ago) by lethuill
Content type: text/plain
Branch: MAIN
CVS Tags: all_2_2_9_03, all_2_2_9_02, all_2_2_9_01
Branch point for: CMSSW_2_2_X_br
Changes since 1.13: +369 -315 lines
Log Message:
Better protection against missing collection / Cleaning data format selection / Last iteration for migration to PAT of Photons

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