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.2 by lethuill, Thu Oct 30 15:58:27 2008 UTC

#	Line 1 | Line 1
1		#include "UserCode/Morgan/interface/PhotonAnalyzer.h"
2	–	#include "DataFormats/EgammaCandidates/interface/PhotonPi0DiscriminatorAssociation.h"
2
3		using namespace std;
4		using namespace reco;
5		using namespace edm;
6
7	<	PhotonAnalyzer::PhotonAnalyzer():verbosity(0)
7	>	PhotonAnalyzer::PhotonAnalyzer(const edm::ParameterSet& producersNames):verbosity_(0), iconvtrack_(0), doPhotonConversion_(true), doVertexCorrection_(true)
8		{
9	+	photonIDProducer_ = producersNames.getParameter<edm::InputTag>("photonIDProducer");
10		}
11	+
12	+	PhotonAnalyzer::PhotonAnalyzer(const edm::ParameterSet& producersNames, int verbosity, bool doPhotonConversion, bool doVertexCorrection):verbosity_(verbosity), iconvtrack_(0), doPhotonConversion_(doPhotonConversion), doVertexCorrection_(doVertexCorrection)
13	+	{
14	+	photonIDProducer_ = producersNames.getParameter<edm::InputTag>("photonIDProducer");
15	+	}
16	+
17		PhotonAnalyzer::~PhotonAnalyzer()
18		{
19		}
20
21	<	void PhotonAnalyzer::Process(const edm::Event& iEvent, TRootEvent* rootEvent, TClonesArray* rootPhotons, string collectionName)
21	>	void PhotonAnalyzer::Process(const edm::Event& iEvent, const edm::EventSetup& iSetup, TRootEvent* rootEvent, TClonesArray* rootPhotons, string collectionName, TClonesArray* conversionTracks, ConversionLikelihoodCalculator* convLikelihoodCalculator, EcalClusterLazyTools& lazyTools)
22		{
23
18	–	if(verbosity>1) cout << "PhotonAnalyzer Collection name: " << collectionName << endl;
24		edm::Handle< reco::PhotonCollection > recoPhotons;
25		iEvent.getByLabel(collectionName, recoPhotons);
26	<	if(verbosity>1) cout << "Number of photons in the event: " << recoPhotons->size() << endl;
27	<
28	<	edm::Handle<reco::PhotonPi0DiscriminatorAssociationMap>  pi0map;
29	<	iEvent.getByLabel("piZeroDiscriminators","PhotonPi0DiscriminatorAssociationMap",  pi0map);
30	<	reco::PhotonPi0DiscriminatorAssociationMap::const_iterator pi0mapIter;
26	>	if(verbosity_>1) std::cout << "   Producer: " << collectionName <<  " PhtonID  Producer: " << photonIDProducer_ << " - Number of photons = " << recoPhotons->size() << std::endl;
27	>
28	>	edm::Handle<reco::PhotonIDAssociationCollection> photonIDMapColl;
29	>	iEvent.getByLabel(photonIDProducer_, photonIDMapColl);
30	>	const reco::PhotonIDAssociationCollection *photonIDMap = photonIDMapColl.product();
31	>
32	>	/*
33	>	// Cluster Shapes with EcalTools
34
35	<	unsigned int iPhot = 0;
36	<	Double_t pi0nn;
35	>	cout << "ecalBarrelRecHitHandle..."<<endl;
36	>	// Get handle to rec hits ecal barrel
37	>	edm::Handle<EcalRecHitCollection>   ecalBarrelRecHitHandle;
38	>	iEvent.getByLabel(reducedBarrelEcalRecHitCollection_, ecalBarrelRecHitHandle);
39	>	if (!ecalBarrelRecHitHandle.isValid())
40	>	{
41	>	edm::LogError("PhotonAnalyzer.cc") << "Error! Can't get the product " << reducedBarrelEcalRecHitCollection_;
42	>	return;
43	>	}
44
45	<	for (unsigned int j=0; j<recoPhotons->size(); j++)
45	>	cout << "ecalEndcapRecHitHandle..."<<endl;
46	>	// Get handle to rec hits ecal endcap
47	>	edm::Handle<EcalRecHitCollection>   ecalEndcapRecHitHandle;
48	>	iEvent.getByLabel(reducedEndcapEcalRecHitCollection_, ecalEndcapRecHitHandle);
49	>	if (!ecalEndcapRecHitHandle.isValid())
50	>	{
51	>	edm::LogError("PhotonAnalyzer.cc") << "Error! Can't get the product " << reducedEndcapEcalRecHitCollection_;
52	>	return;
53	>	}
54	>
55	>	cout << "CaloTopology..."<<endl;
56	>	// Get calorimeter topology
57	>	edm::ESHandle<CaloTopology> caloTopo;
58	>	iSetup.get<CaloTopologyRecord>().get(caloTopo);
59	>	const CaloTopology *topology = caloTopo.product();
60	>	*/
61	>
62	>
63	>	// TODO - add Pi0Disc... not yet implemented in 2.1.X
64	>	//cout << "Pi0Disc..."<<endl;
65	>	//edm::Handle<reco::PhotonPi0DiscriminatorAssociationMap>  pi0map;
66	>	//iEvent.getByLabel("piZeroDiscriminators","PhotonPi0DiscriminatorAssociationMap",  pi0map);
67	>	//reco::PhotonPi0DiscriminatorAssociationMap::const_iterator pi0mapIter;
68	>	//Double_t pi0nn;
69	>
70	>	// Index in TRootPhoton TCloneArray
71	>	unsigned int iPhoton = 0;
72	>
73	>	for (unsigned int iphot=0; iphot<recoPhotons->size(); iphot++)
74		{
75
76	<	pi0mapIter = pi0map->find(edm::Ref<reco::PhotonCollection>(recoPhotons,iPhot));
77	<	pi0nn = ( pi0mapIter == pi0map->end() ? -1 : pi0mapIter->val );
78	<
79	<	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() );
80	<	localPhoton.setHasPixelSeed( (*recoPhotons)[j].hasPixelSeed() );
81	<	localPhoton.setE5x5( (*recoPhotons)[j].e5x5() );
82	<	localPhoton.setR19( (*recoPhotons)[j].r19() );
83	<	localPhoton.setR9( (*recoPhotons)[j].r9() );
84	<	localPhoton.setPi0nn( pi0nn );
85	<
86	<	// FIXME et_had not initialized
87	<	// localPhoton.setEt_had( ( );
88	<
89	<	localPhoton.setSCPosition( (*recoPhotons)[j].superCluster()->position().X(), (*recoPhotons)[j].superCluster()->position().Y(), (*recoPhotons)[j].superCluster()->position().Z() );
90	<	localPhoton.setUnconvPosition( (*recoPhotons)[j].unconvertedPosition().X(), (*recoPhotons)[j].unconvertedPosition().Y(), (*recoPhotons)[j].unconvertedPosition().Z() );
91	<	if(verbosity>2) cout << "  without vtx correction: " << localPhoton << endl;
76	>	// TODO - add Pi0Disc... not yet implemented in 2.1.X
77	>	//pi0mapIter = pi0map->find(edm::Ref<reco::PhotonCollection>(recoPhotons,iPhoton));
78	>	//pi0nn = ( pi0mapIter == pi0map->end() ? -1 : pi0mapIter->val );
79	>
80	>	edm::Ref<reco::PhotonCollection> photonRef(recoPhotons, iphot);
81	>	reco::PhotonIDAssociationCollection::const_iterator photonIter = photonIDMap->find(photonRef);
82	>	const reco::PhotonIDRef &photonID = photonIter->val;
83	>	const reco::PhotonRef &photon = photonIter->key;
84	>
85	>	TRootPhoton localPhoton( photon->px(), photon->py(), photon->pz(), photon->energy(), photon->vx(), photon->vy(), photon->vz(), 22, photon->charge() );
86	>	localPhoton.setHasPixelSeed( photon->hasPixelSeed() );
87	>	localPhoton.setHoE(photon->hadronicOverEm());
88	>	// TODO - add Pi0Disc... not yet implemented in 2.1.X
89	>	//localPhoton.setPi0nn( pi0nn );
90	>
91	>	localPhoton.setSCRawEnergy( photon->superCluster()->rawEnergy() );
92	>	localPhoton.setCaloPosition( photon->caloPosition().X(), photon->caloPosition().Y(), photon->caloPosition().Z() );
93	>
94	>	reco::BasicClusterRef localSeed = photon->superCluster()->seed();
95	>	localPhoton.setEmax( lazyTools.eMax(*localSeed) );
96	>	localPhoton.setE3x3( lazyTools.e3x3(*localSeed) );
97	>	localPhoton.setE5x5( lazyTools.e5x5(*localSeed) );
98	>	localPhoton.setBitsID( photonID->isLooseEM(), photonID->isLoosePhoton(), photonID->isTightPhoton(), photonID->isEBPho(), photonID->isEEPho(), photonID->isEBGap(), photonID->isEEGap(), photonID->isEBEEGap(), photonID->isAlsoElectron() );
99	>
100	>	if(verbosity_>4) std::cout << "seed E5x5=" << lazyTools.e5x5( *localSeed )
101	>	//<< " pi0NN=" << pi0nn
102	>	<< " photonID->isolationEcalRecHit()=" << photonID->isolationEcalRecHit()
103	>	<< " photonID->r9()=" << photonID->r9()
104	>	<< " e3x3 / SC->energy()=" << ( lazyTools.e3x3( *localSeed ) / photon->superCluster()->energy() )
105	>	<< " e3x3 / SC->rawEnergy()=" << ( lazyTools.e3x3( *localSeed ) / photon->superCluster()->rawEnergy() )
106	>	<< " sc->pos=" << photon->superCluster()->position().X() << " , " <<  photon->superCluster()->position().Y() << " , " << photon->superCluster()->position().Z()
107	>	<< " photon->caloPos=" << photon->caloPosition().X() << " , " <<  photon->caloPosition().Y() << " , " << photon->caloPosition().Z()
108	>	<< std::endl;
109
110	<	// Assume photon is coming from primary vertex
51	<	if( rootEvent->nPrimaryVertices()>0 )
110	>	if( doPhotonConversion_ )
111		{
112	<	TVector3 vertex( rootEvent->primaryVertex_x(), rootEvent->primaryVertex_y(), rootEvent->primaryVertex_z() );
113	<	localPhoton.setVertex(vertex);
114	<	if(verbosity>2) cout << "  with vtx correction   : " << localPhoton << endl;
112	>	// If more than one associated convertedPhotons is found, then
113	>	// if nTracks=2, select the convertedPhoton with greatest Likelihood value using Ted's weights file
114	>	// if nTracks=1, select the convertedPhoton with E/p closest to 1
115	>	double eoverp;
116	>	double best_eoverp = 999999.;
117	>	double likely;
118	>	double best_likely = -1.;
119	>	int best_iconv = -1;
120	>	int  best_iconv_likely = -1;
121	>	int  best_iconv_eoverp = -1;
122	>
123	>	std::vector<reco::ConversionRef> conversions = photon->conversions();
124	>	for (unsigned int iconv=0; iconv<conversions.size(); iconv++)
125	>	{
126	>	if(verbosity_>4) cout << "   ["<< setw(3) << iconv << "] Conversion - "
127	>	<< " conv_vertex (x,y,z)=(" << conversions[iconv]->conversionVertex().x()
128	>	<< "," << conversions[iconv]->conversionVertex().y()
129	>	<< "," << conversions[iconv]->conversionVertex().z()
130	>	<< ")  isConverted=" << conversions[iconv]->isConverted()
131	>	<< "  nTracks=" << conversions[iconv]->nTracks()
132	>	<< "  primary_vtx_z=" << conversions[iconv]->zOfPrimaryVertexFromTracks()
133	>	<< "  inv_mass=" << conversions[iconv]->pairInvariantMass()
134	>	<< "  E/p=" << conversions[iconv]->EoverP()
135	>	<< "  cotanTheta=" << conversions[iconv]->pairCotThetaSeparation()
136	>	<< "  likely=" << convLikelihoodCalculator->calculateLikelihood(conversions[iconv])
137	>	<< endl;
138	>
139	>	likely = convLikelihoodCalculator->calculateLikelihood(conversions[iconv]);
140	>	if ( likely>best_likely )
141	>	{
142	>	best_iconv_likely = iconv;
143	>	best_likely = likely;
144	>	}
145	>
146	>	eoverp = conversions[iconv]->EoverP();
147	>	if ( abs(eoverp-1)<abs(best_eoverp-1) )
148	>	{
149	>	best_iconv_eoverp = iconv;
150	>	best_eoverp = eoverp;
151	>	}
152	>	}
153	>
154	>	best_iconv = (best_iconv_likely==-1 ?  best_iconv_eoverp : best_iconv_likely);
155	>
156	>	// Update Photon object with conversion infos
157	>	if ( best_iconv != -1 )
158	>	{
159	>	if(verbosity_>4) cout
160	>	<< "Photon[" << iphot << "]  associated to Conversion[" << best_iconv
161	>	<< "] with E/p="  << conversions[best_iconv]->EoverP()
162	>	<< "  likely=" << convLikelihoodCalculator->calculateLikelihood(conversions[best_iconv])
163	>	<< endl;
164	>
165	>	localPhoton.setConvNTracks( conversions[best_iconv]->nTracks() );
166	>	localPhoton.setConvEoverP( conversions[best_iconv]->EoverP() );
167	>	localPhoton.setConvMass( conversions[best_iconv]->pairInvariantMass() );
168	>	localPhoton.setConvCotanTheta( conversions[best_iconv]->pairCotThetaSeparation() );
169	>	localPhoton.setConvLikely( convLikelihoodCalculator->calculateLikelihood(conversions[best_iconv]) );
170	>	localPhoton.setConvVertex( conversions[best_iconv]->conversionVertex().x(), conversions[best_iconv]->conversionVertex().y(), conversions[best_iconv]->conversionVertex().z() );
171	>	std::vector<math::XYZPoint>  impactVector = conversions[best_iconv]->ecalImpactPosition();
172	>	if ( impactVector.size()>0 ) localPhoton.setConvEcalImpactPosition1( impactVector.at(0).x(), impactVector.at(0).y(), impactVector.at(0).z() );
173	>	if ( impactVector.size()>1 ) localPhoton.setConvEcalImpactPosition2( impactVector.at(1).x(), impactVector.at(1).y(), impactVector.at(1).z() );
174	>
175	>	if ( conversions[best_iconv]->nTracks()>0 )
176	>	{
177	>	std::vector<reco::TrackRef> tracks = conversions[best_iconv]->tracks();
178	>	reco::TrackRef tk1 = tracks.at(0);
179	>	const reco::HitPattern& hit1 = tk1->hitPattern();
180	>	new( (*conversionTracks)[iconvtrack_] ) TRootTrack( tk1->px(), tk1->py(), tk1->pz(), tk1->p(), tk1->vx(), tk1->vy(), tk1->vz(), 0, tk1->charge()
181	>	,hit1.numberOfValidPixelHits(), hit1.numberOfValidTrackerHits(), tk1->chi2(), tk1->d0(), tk1->d0Error(), tk1->dz(), tk1->dzError() );
182	>	// FIXME - Set TRef
183	>	localPhoton.setConvIndexTrack1(iconvtrack_);
184	>	iconvtrack_++;
185	>
186	>	if ( conversions[best_iconv]->nTracks()>1 )
187	>	{
188	>	reco::TrackRef tk2 = tracks.at(1);
189	>	const reco::HitPattern& hit2 = tk2->hitPattern();
190	>	new( (*conversionTracks)[iconvtrack_] ) TRootTrack( tk2->px(), tk2->py(), tk2->pz(), tk2->p(), tk2->vx(), tk2->vy(), tk2->vz(), 0, tk2->charge()
191	>	,hit2.numberOfValidPixelHits(), hit2.numberOfValidTrackerHits(), tk2->chi2(), tk2->d0(), tk2->d0Error(), tk2->dz(), tk2->dzError() );
192	>	// FIXME - Set TRef
193	>	localPhoton.setConvIndexTrack2(iconvtrack_);
194	>	iconvtrack_++;
195	>	}
196	>	}
197	>	}
198		}
199	<	else
199	>
200	>	if(verbosity_>3) cout << "   ["<< setw(3) << iPhoton << "] " << localPhoton << endl;
201	>
202	>	// Vertex correction to photon - Assume photon is coming from primary vertex
203	>	if( doVertexCorrection_)
204		{
205	<	if(verbosity>2) cout << "  with vtx correction: NO PRIMARY VERTEX FOUND !!!" << endl;
205	>	if( rootEvent->nPrimaryVertices()>0 )
206	>	{
207	>	TVector3 vertex( rootEvent->primaryVertex_x(), rootEvent->primaryVertex_y(), rootEvent->primaryVertex_z() );
208	>	localPhoton.setVertex(vertex);
209	>	if(verbosity_>3) cout << "   after vtx correction - ["<< setw(3) << iPhoton << "] " << localPhoton << endl;
210	>	}
211	>	else
212	>	{
213	>	cout << "  PhotonAnalyzer - NO PRIMARY VERTEX FOUND !!!" << endl;
214	>	}
215		}
216
217	<	// Stock new TRootPhoton in TCloneArray
218	<	new( (*rootPhotons)[j] ) TRootPhoton(localPhoton);
219	<	iPhot++;
217	>	// Stock new TRootPhoton in TCloneArray
218	>	new( (*rootPhotons)[iPhoton] ) TRootPhoton(localPhoton);
219	>	iPhoton++;
220		}
66	–
221		}

Diff Legend

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