ViewVC Help

View File | Revision Log | Show Annotations | Root Listing

root/cvsroot/UserCode/Morgan/src/PhotonAnalyzer.cc

Comparing UserCode/Morgan/src/PhotonAnalyzer.cc (file contents):
Revision 1.1 by mlethuil, Mon May 19 16:12:28 2008 UTC vs.
Revision 1.14 by lethuill, Wed Jun 10 11:17:06 2009 UTC

#	Line 1 | Line 1
1	<	#include "UserCode/Morgan/interface/PhotonAnalyzer.h"
2	<	#include "DataFormats/EgammaCandidates/interface/PhotonPi0DiscriminatorAssociation.h"
1	>	#include "../interface/PhotonAnalyzer.h"
2
3		using namespace std;
4		using namespace reco;
5		using namespace edm;
6
7	<	PhotonAnalyzer::PhotonAnalyzer():verbosity(0)
7	>
8	>	PhotonAnalyzer::PhotonAnalyzer(const edm::ParameterSet& producersNames, const edm::ParameterSet& myConfig, int verbosity):verbosity_(verbosity), iconvtrack_(0)
9		{
10	+	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		}
18	+
19	+
20		PhotonAnalyzer::~PhotonAnalyzer()
21		{
22		}
23
15	–	void PhotonAnalyzer::Process(const edm::Event& iEvent, TRootEvent* rootEvent, TClonesArray* rootPhotons, string collectionName)
16	–	{
24
25	<	if(verbosity>1) cout << "PhotonAnalyzer Collection name: " << collectionName << endl;
26	<	edm::Handle< reco::PhotonCollection > recoPhotons;
27	<	iEvent.getByLabel(collectionName, recoPhotons);
28	<	if(verbosity>1) cout << "Number of photons in the event: " << recoPhotons->size() << endl;
22	<
23	<	edm::Handle<reco::PhotonPi0DiscriminatorAssociationMap>  pi0map;
24	<	iEvent.getByLabel("piZeroDiscriminators","PhotonPi0DiscriminatorAssociationMap",  pi0map);
25	<	reco::PhotonPi0DiscriminatorAssociationMap::const_iterator pi0mapIter;
26	<
27	<	unsigned int iPhot = 0;
28	<	Double_t pi0nn;
29	<
30	<	for (unsigned int j=0; j<recoPhotons->size(); j++)
31	<	{
32	<
33	<	pi0mapIter = pi0map->find(edm::Ref<reco::PhotonCollection>(recoPhotons,iPhot));
34	<	pi0nn = ( pi0mapIter == pi0map->end() ? -1 : pi0mapIter->val );
35	<
36	<	TRootPhoton localPhoton( (*recoPhotons)[j].px(), (*recoPhotons)[j].py(), (*recoPhotons)[j].pz(), (*recoPhotons)[j].energy(), (*recoPhotons)[j].vx(), (*recoPhotons)[j].vy(), (*recoPhotons)[j].vz(), 22, (*recoPhotons)[j].charge() );
37	<	localPhoton.setHasPixelSeed( (*recoPhotons)[j].hasPixelSeed() );
38	<	localPhoton.setE5x5( (*recoPhotons)[j].e5x5() );
39	<	localPhoton.setR19( (*recoPhotons)[j].r19() );
40	<	localPhoton.setR9( (*recoPhotons)[j].r9() );
41	<	localPhoton.setPi0nn( pi0nn );
42	<
43	<	// FIXME et_had not initialized
44	<	// localPhoton.setEt_had( ( );
45	<
46	<	localPhoton.setSCPosition( (*recoPhotons)[j].superCluster()->position().X(), (*recoPhotons)[j].superCluster()->position().Y(), (*recoPhotons)[j].superCluster()->position().Z() );
47	<	localPhoton.setUnconvPosition( (*recoPhotons)[j].unconvertedPosition().X(), (*recoPhotons)[j].unconvertedPosition().Y(), (*recoPhotons)[j].unconvertedPosition().Z() );
48	<	if(verbosity>2) cout << "  without vtx correction: " << localPhoton << endl;
49	<
50	<	// Assume photon is coming from primary vertex
51	<	if( rootEvent->nPrimaryVertices()>0 )
52	<	{
53	<	TVector3 vertex( rootEvent->primaryVertex_x(), rootEvent->primaryVertex_y(), rootEvent->primaryVertex_z() );
54	<	localPhoton.setVertex(vertex);
55	<	if(verbosity>2) cout << "  with vtx correction   : " << localPhoton << endl;
56	<	}
57	<	else
58	<	{
59	<	if(verbosity>2) cout << "  with vtx correction: NO PRIMARY VERTEX FOUND !!!" << endl;
60	<	}
61	<
62	<	// Stock new TRootPhoton in TCloneArray
63	<	new( (*rootPhotons)[j] ) TRootPhoton(localPhoton);
64	<	iPhot++;
65	<	}
25	>	bool PhotonAnalyzer::process(const edm::Event& iEvent, const edm::EventSetup& iSetup, TRootEvent* rootEvent, TClonesArray* rootPhotons, TClonesArray* conversionTracks, ConversionLikelihoodCalculator* convLikelihoodCalculator, EcalClusterLazyTools* lazyTools)
26	>	{
27	>
28	>	unsigned int nPhotons=0;
29
30	+	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		}

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines
>	Changed lines