ViewVC Help

View File | Revision Log | Show Annotations | Root Listing

root/cvsroot/UserCode/MitPhysics/Mods/src/PhotonCiCMod.cc

Comparing UserCode/MitPhysics/Mods/src/PhotonCiCMod.cc (file contents):
Revision 1.1 by fabstoec, Wed Jun 1 18:11:52 2011 UTC vs.
Revision 1.3 by fabstoec, Wed Jun 15 10:46:51 2011 UTC

#	Line 18 | Line 18 | PhotonCiCMod::PhotonCiCMod(const char *n
18		fTrackBranchName   (Names::gkTrackBrn),
19		fPileUpDenName     (Names::gkPileupEnergyDensityBrn),
20		fElectronName      ("Electrons"),
21	<	fPhotonPtMin(15.0),
21	>	fPhotonPtMin(20.0),
22		fApplySpikeRemoval(kFALSE),
23		fAbsEtaMax(999.99),
24		fPhotons(0),
#	Line 40 | Line 40 | void PhotonCiCMod::Process()
40		// Process entries of the tree.
41
42		LoadEventObject(fPhotonBranchName,   fPhotons);
43	+
44	+	PhotonOArr *GoodPhotons = new PhotonOArr;
45	+	GoodPhotons->SetName(fGoodPhotonsName);
46
47		Double_t _tRho = -1.;
48		if (fPhotons->GetEntries()>0) {
#	Line 51 | Line 54 | void PhotonCiCMod::Process()
54		if(fPileUpDen->GetEntries() > 0)
55		_tRho = (Double_t) fPileUpDen->At(0)->Rho();
56
57	<	} else return;
57	>	} else {
58	>	AddObjThisEvt(GoodPhotons);
59	>	return;
60	>	}
61	>
62	>	int numVertices = fPV->GetEntries();
63	>	if(numVertices == 0) {
64	>	AddObjThisEvt(GoodPhotons);
65	>	return;
66	>	}
67
68	<	PhotonOArr *GoodPhotons = new PhotonOArr;
69	<	GoodPhotons->SetName(fGoodPhotonsName);
68	>	PhotonOArr* preselPh  = new PhotonOArr;
69	>	PhotonOArr* failingPh = new PhotonOArr;
70	>
71	>	const Vertex** vtxRanked = new const Vertex*[numVertices];
72	>
73	>	// 1. we do the pre-selection; but keep the non-passing photons in a secont container...
74	>	for (UInt_t i=0; i<fPhotons->GetEntries(); ++i) {
75	>	const Photon *ph = fPhotons->At(i);
76	>	if(ph->AbsEta() > 2.5) {
77	>	failingPh->Add(ph);
78	>	continue;
79	>	}
80	>	if(ph->Et() < 30.) {
81	>	failingPh->Add(ph);
82	>	continue;
83	>	}
84	>	if(ph->HadOverEm() > 0.15) {
85	>	failingPh->Add(ph);
86	>	continue;
87	>	}
88	>	if(ph->AbsEta() < 1.5) {
89	>	if(ph->CoviEtaiEta() > 0.013) {
90	>	failingPh->Add(ph);
91	>	continue;
92	>	}
93	>	if(ph->EcalRecHitIsoDr03() > 10.) {
94	>	failingPh->Add(ph);
95	>	continue;
96	>	}
97	>	} else {
98	>	if(ph->CoviEtaiEta() > 0.03) {
99	>	failingPh->Add(ph);
100	>	continue;
101	>	}
102	>	if(ph->EcalRecHitIsoDr03() > 10.) {
103	>	failingPh->Add(ph);
104	>	continue;
105	>	}
106	>	}
107	>	preselPh->Add(ph);
108	>	}
109
110	<	std::vector<const Photon*> phCat1;
111	<	std::vector<const Photon*> phCat2;
112	<	std::vector<const Photon*> phCat3;
113	<	std::vector<const Photon*> phCat4;
110	>	// sort both by pt
111	>	preselPh->Sort();
112	>	failingPh->Sort();
113	>
114	>	// if we don't have two photons in total, nothing to do...
115	>	if( (preselPh->GetEntries() + failingPh->GetEntries() ) < 2) {
116	>	AddObjThisEvt(GoodPhotons);
117	>	delete preselPh;
118	>	delete failingPh;
119	>	return;
120	>	}
121
122	+	// now we assign the two hardest photons to be the ones to select the Vertex...
123	+	double ptggV[2];
124	+	ptggV[0] = 0.;
125	+	ptggV[1] = 0.;
126	+	for(unsigned int iPh =0; iPh < preselPh->GetEntries() && iPh < 2 ; ++iPh) {
127	+	ptggV[0] += preselPh->At(iPh)->Px();
128	+	ptggV[1] += preselPh->At(iPh)->Px();
129	+	}
130	+	for(unsigned int iPh =preselPh->GetEntries(); iPh < (preselPh->GetEntries()+failingPh->GetEntries()) && iPh < 2 ; ++iPh) {
131	+	ptggV[0] += failingPh->At(iPh)->Px();
132	+	ptggV[1] += failingPh->At(iPh)->Px();
133	+	}
134	+	double ptgg = TMath::Sqrt(ptggV[0]*ptggV[0]+ptggV[1]+ptggV[1]);
135	+
136	+	// loop over all vertices and assigne the ranks
137	+	int* ptbal_rank  = new int[fPV->GetEntries()];
138	+	int* ptasym_rank = new int[fPV->GetEntries()];
139	+	int* total_rank  = new int[fPV->GetEntries()];
140	+	double* ptbal = new double[fPV->GetEntries()];
141	+	double* ptasym = new double[fPV->GetEntries()];
142	+
143	+	for(int iVtx = 0; iVtx < numVertices; ++iVtx) {
144	+	const Vertex* _tVtx = fPV->At(iVtx);
145	+	ptbal [iVtx] = 0.0;
146	+	ptasym[iVtx] = 0.0;
147	+	ptbal_rank [iVtx] = 1;
148	+	ptasym_rank[iVtx] = 1;
149	+	double ptvtx = 0.;
150	+	for(unsigned int iTrk = 0; iTrk < _tVtx->NTracks(); ++iTrk) {
151	+	const Track* _tTrk = _tVtx->Trk(iTrk);
152	+	ptvtx += _tTrk->Pt();
153	+	ptbal[iVtx]+= (_tTrk->Px()*ptggV[0])+(_tTrk->Py()*ptggV[1]);
154	+	}
155	+	ptbal [iVtx] = -ptbal[iVtx]/ptgg;
156	+	ptasym[iVtx] = (ptvtx - ptgg)*(ptvtx+ptgg);
157	+	for(int cVtx =0; cVtx < iVtx; ++cVtx) {
158	+	if(ptbal [iVtx] > ptbal [cVtx])
159	+	ptbal_rank [cVtx]++;
160	+	else
161	+	ptbal_rank [iVtx]++;
162	+	if(ptasym [iVtx] > ptasym [cVtx])
163	+	ptasym_rank [cVtx]++;
164	+	else
165	+	ptasym_rank [iVtx]++;
166	+	}
167	+	}
168
169	<	const BaseVertex* theVtx = NULL;
170	<	if(fPV->GetEntries() > 0)
171	<	theVtx = fPV->At(0);
172	<	else return;
169	>	// compute the total rank
170	>	for(int iVtx = 0; iVtx < numVertices; ++iVtx) {
171	>	ptbal_rank [iVtx] = ptbal_rank [iVtx]*ptasym_rank [iVtx]*(iVtx+1);
172	>	total_rank [iVtx] = 0;
173	>	for(int cVtx =0; cVtx < iVtx; ++cVtx) {
174	>	if(ptbal_rank [iVtx] > ptbal_rank [cVtx])
175	>	total_rank[iVtx]++;
176	>	else
177	>	total_rank[cVtx]++;
178	>	}
179	>	}
180
181	<	float cic4_allcuts_temp_sublead[] = {
182	<	3.77459,     2.18305,     1.76811,     1.30029,
183	<	11.5519,     3.48306,     3.87394,     1.89038,
184	<	3.47103,      2.1822,     2.26931,     1.43769,
185	<	0.0105631,  0.00974116,   0.0282572,   0.0265096,
186	<	0.0834648,   0.0821447,   0.0648449,   0.0476437,
187	<	0.94,    0.355935,    0.94,    0.316358,
188	<	98.9979,     1.94497,     96.2292,     96.2294,
189	<	1.5,         1.5,         1.5,         1.5 };   // THIS IS ????
181	>	// sort the Vtx by total rank
182	>	for(int iVtx = 0; iVtx < numVertices; ++iVtx)
183	>	vtxRanked[total_rank[iVtx]] = fPV->At(iVtx);
184	>
185	>	// pick the higest ranked vertex
186	>	const BaseVertex* theVtx = vtxRanked[0];
187	>
188	>	// delete all utility arrays
189	>	delete ptbal_rank  ;
190	>	delete ptasym_rank ;
191	>	delete total_rank  ;
192	>	delete ptbal       ;
193	>	delete ptasym      ;
194
195	+	// fix all the photon kinematics using the new Vtx
196	+	PhotonOArr *fixedPhotons = new PhotonOArr;
197		for (UInt_t i=0; i<fPhotons->GetEntries(); ++i) {
198	<	const Photon *ph = fPhotons->At(i);
198	>	const Photon *ph = fPhotons->At(i);
199	>	// first copy this guy
200	>	Photon* fPh = new Photon(*ph);
201	>	const SuperCluster* sc = ph->SCluster();
202	>	// compute the new Photon momentum (assuming zero mass)
203	>	ThreeVectorC tPh_3mom  = sc->Point() - theVtx->Position();
204	>	tPh_3mom = tPh_3mom/tPh_3mom.R();
205	>	fPh->SetMom(sc->Energy()*tPh_3mom.X(), sc->Energy()*tPh_3mom.Y(), sc->Energy()*tPh_3mom.Z(), sc->Energy());
206	>	fixedPhotons->Add(fPh);
207	>	}
208	>
209	>	// these values are taken from the H2GGlobe code... (actually from Marco/s mail)
210	>	float cic4_allcuts_temp_sublead[] = {
211	>	3.8,         2.2,         1.77,        1.29,
212	>	11.7,        3.4,         3.9,         1.84,
213	>	3.5,         2.2,         2.3,         1.45,
214	>	0.0106,      0.0097,      0.028,       0.027,
215	>	0.082,       0.062,       0.065,       0.048,
216	>	0.94,        0.36,        0.94,        0.32,
217	>	1.,          0.062,       0.97,        0.97,
218	>	1.5,         1.5,         1.5,         1.5 };  // the last line is POixelmatchVeto and un-used
219	>
220	>	for (UInt_t i=0; i<fixedPhotons->GetEntries(); ++i) {
221	>	const Photon *ph = fixedPhotons->At(i);
222
223	<	if (ph->Pt() <= fPhotonPtMin) continue;
224	<	if (ph->AbsEta() >= fAbsEtaMax) continue;
223	>	// cut on Et instead of Pt???
224	>	if ( ph->Pt() <= fPhotonPtMin   ) continue;
225	>	if ( ph->AbsEta() >= fAbsEtaMax ) continue;
226	>
227	>	if (  ph->SCluster()->AbsEta()>=2.5 || (ph->SCluster()->AbsEta()>=1.4442 && ph->SCluster()->AbsEta()<=1.566)  )
228	>	continue;
229
230		Bool_t isbarrel = ph->SCluster()->AbsEta()<1.5;
231
#	Line 101 | Line 245 | void PhotonCiCMod::Process()
245		double trackIso1 = IsolationTools::CiCTrackIsolation(ph, theVtx, 0.3, 0.02, 0.0, 0.0, 0.1, 1.0, fTracks, NULL);
246		double trackIso2 = IsolationTools::CiCTrackIsolation(ph, theVtx, 0.3, 0.02, 0.0, 0.0, 0.1, 1.0, fTracks, fPV);
247		double trackIso3 = IsolationTools::CiCTrackIsolation(ph, theVtx, 0.4, 0.02, 0.0, 0.0, 0.1, 1.0, fTracks, NULL);
248	<
248	>
249		double combIso1 = ecalIso+hcalIso+trackIso1 - 0.17*_tRho;
250		double combIso2 = ecalIso+hcalIso+trackIso2 - 0.52*_tRho;
251
252	<	double tIso1 = (combIso1)*50./ph->Et();
253	<	double tIso2 = (combIso2)*50./ph->Et();
252	>	double tIso1 = (combIso1) *50./ph->Et();
253	>	double tIso2 = (combIso2) *50./ph->Et();
254		double tIso3 = (trackIso3)*50./ph->Et();
255
256		double covIEtaIEta  =ph->CoviEtaiEta();
#	Line 116 | Line 260 | void PhotonCiCMod::Process()
260
261		double dRTrack = PhotonTools::ElectronVetoCiC(ph, fElectrons);
262
263	<	// stupid hard coded, will update later...
264	<	if( tIso1                       < cic4_allcuts_temp_sublead[0+0*4]  &&
265	<	tIso2                       < cic4_allcuts_temp_sublead[0+1*4]  &&
266	<	tIso3                       < cic4_allcuts_temp_sublead[0+2*4]  &&
267	<	covIEtaIEta                 < cic4_allcuts_temp_sublead[0+3*4]  &&
268	<	HoE                         < cic4_allcuts_temp_sublead[0+4*4]  &&
269	<	R9                          > cic4_allcuts_temp_sublead[0+5*4]  &&
270	<	dRTrack*100.                > cic4_allcuts_temp_sublead[0+6*4]  &&
271	<	isbarrel ) phCat1.push_back(ph);
272	<
273	<	if( tIso1                       < cic4_allcuts_temp_sublead[1+0*4]  &&
274	<	tIso2                       < cic4_allcuts_temp_sublead[1+1*4]  &&
275	<	tIso3                       < cic4_allcuts_temp_sublead[1+2*4]  &&
132	<	covIEtaIEta                 < cic4_allcuts_temp_sublead[1+3*4]  &&
133	<	HoE                         < cic4_allcuts_temp_sublead[1+4*4]  &&
134	<	R9                          > cic4_allcuts_temp_sublead[1+5*4]  &&
135	<	dRTrack*100.                > cic4_allcuts_temp_sublead[1+6*4]  &&
136	<	isbarrel ) phCat2.push_back(ph);
137	<
138	<	if( tIso1                       < cic4_allcuts_temp_sublead[2+0*4]  &&
139	<	tIso2                       < cic4_allcuts_temp_sublead[2+1*4]  &&
140	<	tIso3                       < cic4_allcuts_temp_sublead[2+2*4]  &&
141	<	covIEtaIEta                 < cic4_allcuts_temp_sublead[2+3*4]  &&
142	<	HoE                         < cic4_allcuts_temp_sublead[2+4*4]  &&
143	<	R9                          > cic4_allcuts_temp_sublead[2+5*4]  &&
144	<	dRTrack*100.                > cic4_allcuts_temp_sublead[2+6*4]     )  phCat3.push_back(ph);
145	<
146	<	if( tIso1                       < cic4_allcuts_temp_sublead[3+0*4]  &&
147	<	tIso2                       < cic4_allcuts_temp_sublead[3+1*4]  &&
148	<	tIso3                       < cic4_allcuts_temp_sublead[3+2*4]  &&
149	<	covIEtaIEta                 < cic4_allcuts_temp_sublead[3+3*4]  &&
150	<	HoE                         < cic4_allcuts_temp_sublead[3+4*4]  &&
151	<	R9                          > cic4_allcuts_temp_sublead[3+5*4]  &&
152	<	dRTrack*100.                > cic4_allcuts_temp_sublead[3+6*4]     )  phCat4.push_back(ph);
153	<
154	<	}
155	<
156	<	// check which category has two vgalid passes
157	<	bool foundCat=false;
158	<	if(phCat1.size() > 1) {
159	<	// pick the first two guys, assuming they are sorted...
160	<	GoodPhotons->Add(phCat1[0]);
161	<	GoodPhotons->Add(phCat1[1]);
162	<	} else if(phCat2.size() > 1) {
163	<	GoodPhotons->Add(phCat2[0]);
164	<	GoodPhotons->Add(phCat2[1]);
165	<	} else {
166	<	if(phCat3.size() > 1) {
167	<	// one must be a Endcap photons... (stupid, but that's how it is in CiC)
168	<	if(phCat3[0]->SCluster()->AbsEta() > 1.5 || phCat3[1]->SCluster()->AbsEta()) {
169	<	GoodPhotons->Add(phCat3[0]);
170	<	GoodPhotons->Add(phCat3[1]);
171	<	foundCat=true;
172	<	}
173	<	if (!foundCat) {
174	<	if(phCat4.size() > 1) {
175	<	if(phCat4[0]->SCluster()->AbsEta() > 1.5 || phCat4[1]->SCluster()->AbsEta()) {
176	<	GoodPhotons->Add(phCat4[0]);
177	<	GoodPhotons->Add(phCat4[1]);
178	<	}
179	<	}
180	<	}
181	<	}
263	>	// check which category it is ...
264	>	int _tCat = 1;
265	>	if ( !isbarrel ) _tCat = 3;
266	>	if ( R9 < 0.94 ) _tCat++;
267	>
268	>	if( tIso1                       < cic4_allcuts_temp_sublead[_tCat-1+0*4]  &&
269	>	tIso2                       < cic4_allcuts_temp_sublead[_tCat-1+1*4]  &&
270	>	tIso3                       < cic4_allcuts_temp_sublead[_tCat-1+2*4]  &&
271	>	covIEtaIEta                 < cic4_allcuts_temp_sublead[_tCat-1+3*4]  &&
272	>	HoE                         < cic4_allcuts_temp_sublead[_tCat-1+4*4]  &&
273	>	R9                          > cic4_allcuts_temp_sublead[_tCat-1+5*4]  &&
274	>	dRTrack*100.                > cic4_allcuts_temp_sublead[_tCat-1+6*4]  &&
275	>	isbarrel ) GoodPhotons->Add(ph);
276		}
277
278		// sort according to pt
279		GoodPhotons->Sort();
280	<
280	>
281		// add to event for other modules to use
282		AddObjThisEvt(GoodPhotons);
283		}
#	Line 201 | Line 295 | void PhotonCiCMod::SlaveBegin()
295		ReqEventObject(fPVName,             fPV,        fPVFromBranch);
296
297		}
298	+

Diff Legend

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