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.9 by lethuill, Mon Mar 9 13:44:38 2009 UTC vs.
Revision 1.13 by lethuill, Thu Apr 16 14:50:22 2009 UTC

#	Line 6 | Line 6 | using namespace edm;
6
7		PhotonAnalyzer::PhotonAnalyzer(const edm::ParameterSet& producersNames):verbosity_(0), iconvtrack_(0), doPhotonConversion_(true), doVertexCorrection_(true), useMC_(false)
8		{
9	+	dataType_ = producersNames.getUntrackedParameter<string>("dataType","unknown");
10	+	photonProducer_ = producersNames.getParameter<edm::InputTag>("photonProducer");
11		photonIDProducer_ = producersNames.getParameter<edm::InputTag>("photonIDProducer");
12		}
13
14	<	PhotonAnalyzer::PhotonAnalyzer(const edm::ParameterSet& producersNames, const edm::ParameterSet& myConfig, int verbosity):verbosity_(verbosity)
14	>	PhotonAnalyzer::PhotonAnalyzer(const edm::ParameterSet& producersNames, const edm::ParameterSet& myConfig, int verbosity):verbosity_(verbosity), iconvtrack_(0)
15		{
16		dataType_ = producersNames.getUntrackedParameter<string>("dataType","unknown");
17	+	photonProducer_ = producersNames.getParameter<edm::InputTag>("photonProducer");
18		photonIDProducer_ = producersNames.getParameter<edm::InputTag>("photonIDProducer");
19		doPhotonConversion_ = myConfig.getUntrackedParameter<bool>("doPhotonConversion");
20		doVertexCorrection_ = myConfig.getUntrackedParameter<bool>("doPhotonVertexCorrection");
#	Line 22 | Line 25 | PhotonAnalyzer::~PhotonAnalyzer()
25		{
26		}
27
28	<	void PhotonAnalyzer::Process(const edm::Event& iEvent, const edm::EventSetup& iSetup, TRootEvent* rootEvent, TClonesArray* rootPhotons, string collectionName, TClonesArray* conversionTracks, ConversionLikelihoodCalculator* convLikelihoodCalculator, EcalClusterLazyTools& lazyTools)
28	>	void PhotonAnalyzer::Process(const edm::Event& iEvent, const edm::EventSetup& iSetup, TRootEvent* rootEvent, TClonesArray* rootPhotons, TClonesArray* conversionTracks, ConversionLikelihoodCalculator* convLikelihoodCalculator, EcalClusterLazyTools& lazyTools)
29		{
30
31	+	unsigned int nPhotons=0;
32	+
33		edm::Handle< reco::PhotonCollection > recoPhotons;
34	<	iEvent.getByLabel(collectionName, recoPhotons);
35	<	if(verbosity_>1) std::cout << "   Producer: " << collectionName <<  " PhotonID  Producer: " << photonIDProducer_.label() << " - Number of photons = " << recoPhotons->size() << std::endl;
34	>	const reco::PhotonIDAssociationCollection *photonIDMap = 0;
35	>	if( dataType_=="RECO" || dataType_=="AOD" )
36	>	{
37	>	iEvent.getByLabel(photonProducer_, recoPhotons);
38	>	nPhotons = recoPhotons->size();
39	>	// Photon identification
40	>	edm::Handle<reco::PhotonIDAssociationCollection> photonIDMapColl;
41	>	iEvent.getByLabel(photonIDProducer_, photonIDMapColl);
42	>	photonIDMap = photonIDMapColl.product();
43	>	}
44
45	<	edm::Handle<reco::PhotonIDAssociationCollection> photonIDMapColl;
46	<	iEvent.getByLabel(photonIDProducer_, photonIDMapColl);
47	<	const reco::PhotonIDAssociationCollection *photonIDMap = photonIDMapColl.product();
45	>	edm::Handle < std::vector <pat::Photon> > patPhotons;
46	>	if( dataType_=="PAT" || dataType_=="PATAOD" )
47	>	{
48	>	iEvent.getByLabel(photonProducer_, patPhotons);
49	>	nPhotons = patPhotons->size();
50	>	}
51
52	+	if(verbosity_>1) std::cout << "   Number of photons = " << nPhotons << "   Label: " << photonProducer_.label() << "   Instance: " << photonProducer_.instance() << std::endl;
53
54	<	// TODO - add Pi0Disc... not yet implemented in 2.1.X
54	>
55	>	// TODO - add Pi0Disc... not yet implemented in 2.X.X
56		//cout << "Pi0Disc..."<<endl;
57		//edm::Handle<reco::PhotonPi0DiscriminatorAssociationMap>  pi0map;
58		//iEvent.getByLabel("piZeroDiscriminators","PhotonPi0DiscriminatorAssociationMap",  pi0map);
59		//reco::PhotonPi0DiscriminatorAssociationMap::const_iterator pi0mapIter;
60		//Double_t pi0nn;
61
44	–	// Index in TRootPhoton TCloneArray
45	–	unsigned int iPhoton = 0;
62
63	<	for (unsigned int iphot=0; iphot<recoPhotons->size(); iphot++)
63	>	for (unsigned int j=0; j<nPhotons; j++)
64		{
65	+	const reco::Photon* photon = 0;
66	+	if( dataType_=="RECO" || dataType_=="AOD" ) photon =  &((*recoPhotons)[j]);
67	+	if( dataType_=="PAT" || dataType_=="PATAOD" ) photon = (const reco::Photon*) ( & ((*patPhotons)[j]) );
68	+
69	+	TRootPhoton localPhoton(
70	+	photon->px()
71	+	,photon->py()
72	+	,photon->pz()
73	+	,photon->energy()
74	+	,photon->vx()
75	+	,photon->vy()
76	+	,photon->vz()
77	+	,photon->pdgId()
78	+	,photon->charge()
79	+	);
80
81	<	// TODO - add Pi0Disc... not yet implemented in 2.1.X
82	<	//pi0mapIter = pi0map->find(edm::Ref<reco::PhotonCollection>(recoPhotons,iPhoton));
83	<	//pi0nn = ( pi0mapIter == pi0map->end() ? -1 : pi0mapIter->val );
81	>	// Variables from reco::Photon
82	>	localPhoton.setCaloPosition( photon->caloPosition().X(), photon->caloPosition().Y(), photon->caloPosition().Z() );
83	>	localPhoton.setHoE(photon->hadronicOverEm());
84	>	localPhoton.setHasPixelSeed( photon->hasPixelSeed() );
85
86	<	edm::Ref<reco::PhotonCollection> photonRef(recoPhotons, iphot);
87	<	reco::PhotonIDAssociationCollection::const_iterator photonIter = photonIDMap->find(photonRef);
88	<	const reco::PhotonIDRef &photonID = photonIter->val;
89	<	const reco::PhotonRef &photon = photonIter->key;
86	>	// Variables from reco::SuperCluster
87	>	reco::SuperClusterRef superCluster = photon->superCluster();
88	>	//const reco::SuperCluster* sc = superCluster.get();
89	>	if ( superCluster.isNonnull() )
90	>	{
91	>	localPhoton.setSCRawEnergy( superCluster->rawEnergy() );
92	>	//localPhoton.setNbClusters(electron->numberOfClusters());
93	>	//localPhoton.setPreshowerEnergy(superCluster->preshowerEnergy());
94	>	}
95
96	<	TRootPhoton localPhoton( photon->px(), photon->py(), photon->pz(), photon->energy(), photon->vx(), photon->vy(), photon->vz(), photon->pdgId(), photon->charge() );
97	<	localPhoton.setHasPixelSeed( photon->hasPixelSeed() );
98	<	localPhoton.setHoE(photon->hadronicOverEm());
99	<	// TODO - add Pi0Disc... not yet implemented in 2.1.X
100	<	//localPhoton.setPi0nn( pi0nn );
96	>	// Cluster Shape variables
97	>	// need reco::SuperCluster and reco::BasicCluster
98	>	if ( superCluster.isNonnull() )
99	>	{
100	>	reco::BasicClusterRef seedBasicCluster = superCluster->seed();
101	>	//if ( seedBasicCluster.isNonnull() ) localPhoton.setClusterAlgo(seedBasicCluster->algo());
102
103	<	localPhoton.setSCRawEnergy( photon->superCluster()->rawEnergy() );
104	<	localPhoton.setCaloPosition( photon->caloPosition().X(), photon->caloPosition().Y(), photon->caloPosition().Z() );
103	>	/*
104	>	// dR of the cone centered on the reco::GsfElectron and containing all its basic clusters constituents
105	>	Float_t caloConeSize = 0;
106	>	for (reco::basicCluster_iterator basicCluster = (*superCluster).clustersBegin(); basicCluster != (*superCluster).clustersEnd(); ++basicCluster )
107	>	{
108	>	Float_t dR = localElectron.DeltaR(TLorentzVector( (*basicCluster)->position().x(), (*basicCluster)->position().y(), (*basicCluster)->position().z(), 0. ) );
109	>	if (dR > caloConeSize) caloConeSize = dR;
110	>	}
111	>	// FIXME - if no BasicCluster collection, init to -999.
112	>	localElectron.setCaloConeSize(caloConeSize);
113	>	*/
114	>
115	>	// need reduced Ecal RecHits Collections for EcalClusterLazyTools
116	>	// FIXME - check lazyTools is initialized ==> lazyTools has to be a pointer
117	>	if ( seedBasicCluster.isNonnull() )
118	>	{
119	>	//localPhoton.setE2x2(lazyTools.e2x2(*seedBasicCluster));
120	>	localPhoton.setE3x3( lazyTools.e3x3(*seedBasicCluster) );
121	>	localPhoton.setE5x5( lazyTools.e5x5(*seedBasicCluster) );
122	>	localPhoton.setEmax( lazyTools.eMax(*seedBasicCluster) );
123	>	}
124	>	}
125
68	–	reco::BasicClusterRef localSeed = photon->superCluster()->seed();
69	–	localPhoton.setEmax( lazyTools.eMax(*localSeed) );
70	–	localPhoton.setE3x3( lazyTools.e3x3(*localSeed) );
71	–	localPhoton.setE5x5( lazyTools.e5x5(*localSeed) );
72	–	localPhoton.setBitsID( photonID->isLooseEM(), photonID->isLoosePhoton(), photonID->isTightPhoton(), photonID->isEBPho(), photonID->isEEPho(), photonID->isEBGap(), photonID->isEEGap(), photonID->isEBEEGap(), photonID->isAlsoElectron() );
126
74	–	if(verbosity_>4) std::cout << "seed E5x5=" << lazyTools.e5x5( *localSeed )
75	–	//<< " pi0NN=" << pi0nn
76	–	<< " photonID->isolationEcalRecHit()=" << photonID->isolationEcalRecHit()
77	–	<< " photonID->r9()=" << photonID->r9()
78	–	<< " e3x3 / SC->energy()=" << ( lazyTools.e3x3( *localSeed ) / photon->superCluster()->energy() )
79	–	<< " e3x3 / SC->rawEnergy()=" << ( lazyTools.e3x3( *localSeed ) / photon->superCluster()->rawEnergy() )
80	–	<< " sc->pos=" << photon->superCluster()->position().X() << " , " <<  photon->superCluster()->position().Y() << " , " << photon->superCluster()->position().Z()
81	–	<< " photon->caloPos=" << photon->caloPosition().X() << " , " <<  photon->caloPosition().Y() << " , " << photon->caloPosition().Z()
82	–	<< std::endl;
127
128	+	// FIXME - Add pi0nn
129	+	//pi0mapIter = pi0map->find(edm::Ref<reco::PhotonCollection>(recoPhotons,iPhoton));
130	+	//pi0nn = ( pi0mapIter == pi0map->end() ? -1 : pi0mapIter->val );
131	+	//localPhoton.setPi0nn( pi0nn );
132	+
133	+
134	+
135	+	// Photon conversions
136	+	// Variables from reco::Conversion
137		if( doPhotonConversion_ )
138		{
139		// If more than one associated convertedPhotons is found, then
#	Line 90 | Line 143 | void PhotonAnalyzer::Process(const edm::
143		double best_eoverp = 999999.;
144		double likely;
145		double best_likely = -1.;
146	<	int best_iconv = -1;
147	<	int  best_iconv_likely = -1;
148	<	int  best_iconv_eoverp = -1;
146	>	int best_iconv = -1;
147	>	int  best_iconv_likely = -1;
148	>	int  best_iconv_eoverp = -1;
149
150		std::vector<reco::ConversionRef> conversions = photon->conversions();
151		for (unsigned int iconv=0; iconv<conversions.size(); iconv++)
#	Line 121 | Line 174 | void PhotonAnalyzer::Process(const edm::
174		if ( abs(eoverp-1)<abs(best_eoverp-1) )
175		{
176		best_iconv_eoverp = iconv;
177	<	best_eoverp = eoverp;
177	>	best_eoverp = eoverp;
178		}
179		}
180
#	Line 131 | Line 184 | void PhotonAnalyzer::Process(const edm::
184		if ( best_iconv != -1 )
185		{
186		if(verbosity_>4) cout
187	<	<< "Photon[" << iphot << "]  associated to Conversion[" << best_iconv
187	>	<< "   ===> Photon[" << j << "] associated to Conversion[" << best_iconv
188		<< "] with E/p="  << conversions[best_iconv]->EoverP()
189		<< "  likely=" << convLikelihoodCalculator->calculateLikelihood(conversions[best_iconv])
190		<< endl;
191	<
191	>
192		localPhoton.setConvNTracks( conversions[best_iconv]->nTracks() );
193		localPhoton.setConvEoverP( conversions[best_iconv]->EoverP() );
194		localPhoton.setConvMass( conversions[best_iconv]->pairInvariantMass() );
#	Line 146 | Line 199 | void PhotonAnalyzer::Process(const edm::
199		if ( impactVector.size()>0 ) localPhoton.setConvEcalImpactPosition1( impactVector.at(0).x(), impactVector.at(0).y(), impactVector.at(0).z() );
200		if ( impactVector.size()>1 ) localPhoton.setConvEcalImpactPosition2( impactVector.at(1).x(), impactVector.at(1).y(), impactVector.at(1).z() );
201
202	<	if ( conversions[best_iconv]->nTracks()>0 )
203	<	{
204	<	std::vector<reco::TrackRef> tracks = conversions[best_iconv]->tracks();
205	<	reco::TrackRef tk1 = tracks.at(0);
206	<	const reco::HitPattern& hit1 = tk1->hitPattern();
207	<	new( (*conversionTracks)[iconvtrack_] ) TRootTrack( tk1->px(), tk1->py(), tk1->pz(), tk1->p(), tk1->vx(), tk1->vy(), tk1->vz(), 0, tk1->charge()
208	<	,hit1.numberOfValidPixelHits(), hit1.numberOfValidTrackerHits(), tk1->chi2(), tk1->d0(), tk1->d0Error(), tk1->dz(), tk1->dzError() );
209	<	localPhoton.setConvIndexTrack1(iconvtrack_);
210	<	localPhoton.setConvTrack1((*conversionTracks)[iconvtrack_]);
211	<	iconvtrack_++;
212	<
160	<	if ( conversions[best_iconv]->nTracks()>1 )
161	<	{
162	<	reco::TrackRef tk2 = tracks.at(1);
163	<	const reco::HitPattern& hit2 = tk2->hitPattern();
164	<	new( (*conversionTracks)[iconvtrack_] ) TRootTrack( tk2->px(), tk2->py(), tk2->pz(), tk2->p(), tk2->vx(), tk2->vy(), tk2->vz(), 0, tk2->charge()
165	<	,hit2.numberOfValidPixelHits(), hit2.numberOfValidTrackerHits(), tk2->chi2(), tk2->d0(), tk2->d0Error(), tk2->dz(), tk2->dzError() );
166	<	localPhoton.setConvIndexTrack2(iconvtrack_);
167	<	localPhoton.setConvTrack2((*conversionTracks)[iconvtrack_]);
202	>	if ( conversionTracks!=0 )  // Check branch is activated
203	>	{
204	>	if ( conversions[best_iconv]->nTracks()>0 )
205	>	{
206	>	std::vector<reco::TrackRef> tracks = conversions[best_iconv]->tracks();
207	>	reco::TrackRef tk1 = tracks.at(0);
208	>	const reco::HitPattern& hit1 = tk1->hitPattern();
209	>	new( (*conversionTracks)[iconvtrack_] ) TRootTrack( tk1->px(), tk1->py(), tk1->pz(), tk1->p(), tk1->vx(), tk1->vy(), tk1->vz(), 0, tk1->charge()
210	>	,hit1.pixelLayersWithMeasurement(), hit1.stripLayersWithMeasurement(), tk1->chi2(), tk1->d0(), tk1->d0Error(), tk1->dz(), tk1->dzError() );
211	>	localPhoton.setConvIndexTrack1(iconvtrack_);
212	>	localPhoton.setConvTrack1((*conversionTracks)[iconvtrack_]);
213		iconvtrack_++;
214	+
215	+	if ( conversions[best_iconv]->nTracks()>1 )
216	+	{
217	+	reco::TrackRef tk2 = tracks.at(1);
218	+	const reco::HitPattern& hit2 = tk2->hitPattern();
219	+	new( (*conversionTracks)[iconvtrack_] ) TRootTrack( tk2->px(), tk2->py(), tk2->pz(), tk2->p(), tk2->vx(), tk2->vy(), tk2->vz(), 0, tk2->charge()
220	+	,hit2.pixelLayersWithMeasurement(), hit2.stripLayersWithMeasurement(), tk2->chi2(), tk2->d0(), tk2->d0Error(), tk2->dz(), tk2->dzError() );
221	+	localPhoton.setConvIndexTrack2(iconvtrack_);
222	+	localPhoton.setConvTrack2((*conversionTracks)[iconvtrack_]);
223	+	iconvtrack_++;
224	+	}
225		}
226		}
227		}
228		}
229
230	<	if(verbosity_>3) cout << "   ["<< setw(3) << iPhoton << "] " << localPhoton << endl;
230	>
231	>	if( dataType_=="RECO" || dataType_=="AOD" )
232	>	{
233	>	// Some specific methods requiring  RECO / AOD format
234	>	// Do association to genParticle ?
235	>
236	>	// Photon ID
237	>	edm::Ref<reco::PhotonCollection> photonRef(recoPhotons, j);
238	>	reco::PhotonIDAssociationCollection::const_iterator photonIter = photonIDMap->find(photonRef);
239	>	const reco::PhotonIDRef &photonID = photonIter->val;
240	>	//const reco::PhotonRef &photon = photonIter->key;
241	>	localPhoton.setBitsID(
242	>	photonID->isLooseEM()
243	>	,photonID->isLoosePhoton()
244	>	,photonID->isTightPhoton()
245	>	,photonID->isEBPho()
246	>	,photonID->isEEPho()
247	>	,photonID->isEBGap()
248	>	,photonID->isEEGap()
249	>	,photonID->isEBEEGap()
250	>	,photonID->isAlsoElectron()
251	>	);
252	>
253	>	/*
254	>	if(verbosity_>4) std::cout << "seed E5x5=" << lazyTools.e5x5( *seedBasicCluster )
255	>	//<< " pi0NN=" << pi0nn
256	>	<< " photonID->isolationEcalRecHit()=" << photonID->isolationEcalRecHit()
257	>	<< " photonID->r9()=" << photonID->r9()
258	>	<< " e3x3 / SC->energy()=" << ( lazyTools.e3x3( *seedBasicCluster ) / photon->superCluster()->energy() )
259	>	<< " e3x3 / SC->rawEnergy()=" << ( lazyTools.e3x3( *seedBasicCluster ) / photon->superCluster()->rawEnergy() )
260	>	<< " sc->pos=" << photon->superCluster()->position().X() << " , " <<  photon->superCluster()->position().Y() << " , " << photon->superCluster()->position().Z()
261	>	<< " photon->caloPos=" << photon->caloPosition().X() << " , " <<  photon->caloPosition().Y() << " , " << photon->caloPosition().Z()                                     << std::endl;
262	>	*/
263	>	}
264	>
265	>
266	>	if( dataType_=="PATAOD" || dataType_=="PAT" )
267	>	{
268	>	// Some specific methods to pat::Photon
269	>	const pat::Photon *patPhoton = dynamic_cast<const pat::Photon*>(&*photon);
270	>
271	>	// Photon ID
272	>	localPhoton.setBitsID(
273	>	patPhoton->isLooseEM()
274	>	,patPhoton->isLoosePhoton()
275	>	,patPhoton->isTightPhoton()
276	>	,patPhoton->isEBPho()
277	>	,patPhoton->isEEPho()
278	>	,patPhoton->isEBGap()
279	>	,patPhoton->isEEGap()
280	>	,patPhoton->isEBEEGap()
281	>	,patPhoton->isAlsoElectron()
282	>	);
283	>
284	>	// Isolation ?
285	>	/*
286	>	pair < Float_t, Int_t > trackerIso;
287	>	trackerIso = patElectron->trackerIsoDeposit()->depositAndCountWithin(0.1);
288	>	localElectron.setIsoR01_sumPt(trackerIso.first);
289	>	localElectron.setIsoR01_nTracks(trackerIso.second);
290	>	*/
291	>
292	>
293	>	if(useMC_)
294	>	{
295	>	// MC truth associator index
296	>	if ((patPhoton->genParticleRef()).isNonnull()) {
297	>	localPhoton.setGenParticleIndex((patPhoton->genParticleRef()).index());
298	>	} else {
299	>	localPhoton.setGenParticleIndex(-1);
300	>	}
301	>	}
302	>	}
303	>
304	>
305	>
306	>
307	>	// FIXME - associator supercluster <-> photon
308	>	//localElectron.setSCRef(superCluster->toto());
309	>
310	>
311
312		// Vertex correction to photon - Assume photon is coming from primary vertex
313		if( doVertexCorrection_)
314		{
315	<	if( rootEvent->nPrimaryVertices()>0 )
315	>	if(verbosity_>3) cout << "   Before vertex correction  ["<< setw(3) << j << "] " << localPhoton << endl;
316	>	if( rootEvent->primaryVertex()!=0 )
317		{
318	<	TVector3 vertex( rootEvent->primaryVertex_x(), rootEvent->primaryVertex_y(), rootEvent->primaryVertex_z() );
318	>	TVector3 vertex( rootEvent->primaryVertex()->x(), rootEvent->primaryVertex()->y(), rootEvent->primaryVertex()->z() );
319		localPhoton.setVertex(vertex);
183	–	if(verbosity_>3) cout << "   after vtx correction - ["<< setw(3) << iPhoton << "] " << localPhoton << endl;
320		}
321		else
322		{
323	<	cout << "  PhotonAnalyzer - NO PRIMARY VERTEX FOUND !!!" << endl;
323	>	cout << "  PhotonAnalyzer - NO SELECTED PRIMARY VERTEX FOUND !!!" << endl;
324		}
325		}
326
191	–	/*if(useMC_)
192	–	{
193	–	if( dataType_=="PAT" || dataType_=="PATAOD" )
194	–	// MC truth associator index
195	–	if ((photon.genParticleRef()).isNonnull()) {
196	–	localPhoton.setGenParticleIndex((photon.genParticleRef()).index());
197	–	} else {
198	–	localPhoton.setGenParticleIndex(-1);
199	–	}
200	–	}
201	–	}*/
327
328		// Stock new TRootPhoton in TCloneArray
329	<	new( (*rootPhotons)[iPhoton] ) TRootPhoton(localPhoton);
330	<	iPhoton++;
329	>	new( (*rootPhotons)[j] ) TRootPhoton(localPhoton);
330	>	//if(verbosity_>2) cout << "   ["<< setw(3) << j << "] " << localPhoton << endl;
331		}
332		}
333	+

Diff Legend

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