ViewVC Help

View File | Revision Log | Show Annotations | Root Listing

root/cvsroot/UserCode/MitPhysics/Utils/src/MVATools.cc

Comparing UserCode/MitPhysics/Utils/src/MVATools.cc (file contents):
Revision 1.6 by mingyang, Wed Dec 7 00:45:17 2011 UTC vs.
Revision 1.15 by mingyang, Mon Jun 11 09:04:33 2012 UTC

#	Line 4 | Line 4
4		#include "MitPhysics/Utils/interface/MVATools.h"
5		#include "MitPhysics/Utils/interface/ElectronTools.h"
6		#include "MitPhysics/Utils/interface/IsolationTools.h"
7	+	#include "MitAna/DataTree/interface/PFCandidateCol.h"
8		#include "MitAna/DataTree/interface/StableData.h"
9	+	#include <TMath.h>
10		#include <TFile.h>
11		#include <TRandom3.h>
12		#include "TMVA/Tools.h"//MVA
13		#include "TMVA/Reader.h"//MVA
14
15	+
16		ClassImp(mithep::MVATools)
17
18		using namespace mithep;
#	Line 55 | Line 58 | MVATools::MVATools():
58		RelIsoHcal(0),
59		tIso1(0),
60		tIso3(0),
61	<	tIso2(0)
61	>	tIso2(0),
62	>	ScEta(0.)
63		{
64		// Constructor.
65		}
#	Line 157 | Line 161 | void MVATools::InitializeMVA(int Variabl
161		readers[i]->AddVariable( "RelPreshowerEnergy", &RelPreshowerEnergy );
162		}
163		}
164	+
165	+	if(VariableType==7){
166	+	readers[i]->AddVariable( "HoE", &HoE );
167	+	readers[i]->AddVariable( "covIEtaIEta", &covIEtaIEta );
168	+	readers[i]->AddVariable( "tIso1abs", &tIso1abs );
169	+	readers[i]->AddVariable( "tIso3abs", &tIso3abs );
170	+	readers[i]->AddVariable( "tIso2abs", &tIso2abs );
171	+	readers[i]->AddVariable( "R9", &R9 );
172	+	readers[i]->AddVariable( "absIsoEcal", &absIsoEcal );
173	+	readers[i]->AddVariable( "absIsoHcal", &absIsoHcal );
174	+	readers[i]->AddVariable( "NVertexes", &NVertexes );
175	+	readers[i]->AddVariable( "ScEta", &ScEta );
176	+	}
177	+
178	+	if(VariableType==10){
179	+	readers[i]->AddVariable( "HoE", &HoE );
180	+	readers[i]->AddVariable( "covIEtaIEta", &covIEtaIEta );
181	+	readers[i]->AddVariable( "tIso1abs", &tIso1abs );
182	+	readers[i]->AddVariable( "tIso3abs", &tIso3abs );
183	+	readers[i]->AddVariable( "tIso2abs", &tIso2abs );
184	+	readers[i]->AddVariable( "R9", &R9 );
185	+	readers[i]->AddVariable( "absIsoEcal", &absIsoEcal );
186	+	readers[i]->AddVariable( "absIsoHcal", &absIsoHcal );
187	+	readers[i]->AddVariable( "NVertexes", &NVertexes );
188	+	readers[i]->AddVariable( "ScEta", &ScEta );
189	+	readers[i]->AddVariable( "EtaWidth", &EtaWidth );
190	+	readers[i]->AddVariable( "PhiWidth", &PhiWidth );
191	+	}
192
193	+	if(VariableType==1201){
194	+	/*readers[i]->AddVariable( "myphoton_pfchargedisogood03", &myphoton_pfchargedisogood03);
195	+	readers[i]->AddVariable( "myphoton_pfchargedisobad03", &myphoton_pfchargedisobad03);
196	+	readers[i]->AddVariable( "myphoton_pfphotoniso03", &myphoton_pfphotoniso03 );
197	+	readers[i]->AddVariable( "myphoton_sieie", &myphoton_sieie );
198	+	readers[i]->AddVariable( "myphoton_sieip", &myphoton_sieip );
199	+	readers[i]->AddVariable( "myphoton_etawidth", &myphoton_etawidth );
200	+	readers[i]->AddVariable( "myphoton_phiwidth", &myphoton_phiwidth );
201	+	readers[i]->AddVariable( "myphoton_r9", &myphoton_r9 );
202	+	readers[i]->AddVariable( "myphoton_s4ratio", &myphoton_s4ratio );
203	+	readers[i]->AddVariable( "myphoton_SCeta", &myphoton_SCeta );
204	+	readers[i]->AddVariable( "event_rho", &event_rho );
205	+	if(i==0){
206	+	readers[i]->AddVariable( "myphoton_ESEffSigmaRR", &myphoton_ESEffSigmaRR);
207	+	}*/
208	+	readers[i]->AddVariable( "ph.r9", &myphoton_r9 );
209	+	readers[i]->AddVariable( "ph.sigietaieta", &myphoton_sieie );
210	+	readers[i]->AddVariable( "ph.scetawidth", &myphoton_etawidth );
211	+	readers[i]->AddVariable( "ph.scphiwidth", &myphoton_phiwidth );
212	+	readers[i]->AddVariable( "ph.idmva_CoviEtaiPhi", &myphoton_sieip );
213	+	readers[i]->AddVariable( "ph.idmva_s4ratio", &myphoton_s4ratio );
214	+	readers[i]->AddVariable( "ph.idmva_GammaIso", &myphoton_pfphotoniso03 );
215	+	readers[i]->AddVariable( "ph.idmva_ChargedIso_selvtx", &myphoton_pfchargedisogood03);
216	+	readers[i]->AddVariable( "ph.idmva_ChargedIso_worstvtx", &myphoton_pfchargedisobad03);
217	+	readers[i]->AddVariable( "ph.sceta", &myphoton_SCeta );
218	+	readers[i]->AddVariable( "rho", &event_rho );
219	+	if(i==0){
220	+	readers[i]->AddVariable( "1.00023*ph.idmva_PsEffWidthSigmaRR + 0.0913", &myphoton_ESEffSigmaRR);
221	+	}
222	+	}
223	+
224		}
225
226		fReaderEndcap->BookMVA("BDT method",EndcapWeights);
#	Line 168 | Line 231 | void MVATools::InitializeMVA(int Variabl
231
232		}
233
234	<	Bool_t MVATools::PassMVASelection(const Photon* p,const Vertex* vtx,const TrackCol* trackCol,const VertexCol* vtxCol,Double_t _tRho,Float_t bdtCutBarrel, Float_t bdtCutEndcap) {
234	>	Bool_t MVATools::PassMVASelection(const Photon* p,const Vertex* vtx,const TrackCol* trackCol,const VertexCol* vtxCol,Double_t _tRho,Float_t bdtCutBarrel, Float_t bdtCutEndcap, const ElectronCol* els, Bool_t applyElectronVeto) {
235
236		//initilize the bool value
237		PassMVA=kFALSE;
238
239	<	Float_t photon_bdt =  MVATools::GetMVAbdtValue(p,vtx,trackCol,vtxCol, _tRho);
239	>	Float_t photon_bdt =  MVATools::GetMVAbdtValue_2011(p,vtx,trackCol,vtxCol, _tRho, els, applyElectronVeto);
240
241		if (isbarrel) {
242		if(bdt>bdtCutBarrel){
#	Line 212 | Line 275 | Int_t MVATools::PassElectronVetoInt(cons
275		isbarrel = (fabs(ScEta_MVA)<1.4442);
276
277		R9 = p->R9();
278	+	//R9 = p->E33()/p->SCluster()->RawEnergy();
279
280		// check which category it is ...
281		_tCat = 1;
#	Line 229 | Line 293 | Int_t MVATools::PassElectronVetoInt(cons
293
294		//--------------------------------------------------------------------------------------------------
295
296	<	Float_t MVATools::GetMVAbdtValue(const Photon* p,const Vertex* vtx,const TrackCol* trackCol,const VertexCol* vtxCol,Double_t _tRho) {
296	>	Float_t MVATools::GetMVAbdtValue_2012_globe(const Photon* p,const Vertex* vtx,const TrackCol* trackCol,const VertexCol* vtxCol,Double_t _tRho, const PFCandidateCol *fPFCands,const ElectronCol* els,Bool_t applyElectronVeto) {
297
298		//get the variables used to compute MVA variables
299		ecalIso3 = p->EcalRecHitIsoDr03();
#	Line 238 | Line 302 | Float_t MVATools::GetMVAbdtValue(const P
302
303		wVtxInd = 0;
304
305	<	trackIso1 = IsolationTools::CiCTrackIsolation(p,vtx, 0.3, 0.02, 0.0, 0.0, 0.1, 1.0, trackCol);//Question Ming:whyfPV->At(0) instead of selected vertex using ranking method?
305	>	trackIso1 = IsolationTools::CiCTrackIsolation(p,vtx, 0.3, 0.02, 0.0, 0.0, 0.1, 1.0, trackCol, NULL, NULL, (!applyElectronVeto ? els : NULL) );//Question Ming:whyfPV->At(0) instead of selected vertex using ranking method?
306	>
307	>	// track iso worst vtx
308	>	trackIso2 = IsolationTools::CiCTrackIsolation(p,vtx, 0.4, 0.02, 0.0, 0.0, 0.1, 1.0, trackCol, &wVtxInd,vtxCol, (!applyElectronVeto ? els : NULL) );
309
310	<	// track iso only
311	<	trackIso3 = IsolationTools::CiCTrackIsolation(p,vtx, 0.3, 0.02, 0.0, 0.0, 0.1, 1.0, trackCol);
310	>	combIso1 = ecalIso3+hcalIso4+trackIso1 - 0.17*_tRho;
311	>	combIso2 = ecalIso4+hcalIso4+trackIso2 - 0.52*_tRho;
312	>
313	>	RawEnergy = p->SCluster()->RawEnergy();
314
315	+	ScEta = p->SCluster()->Eta();
316	+
317	+	//mva varialbes v1 and v2
318	+	tIso1 = (combIso1) *50./p->Et();
319	+	tIso3 = (trackIso1)*50./p->Et();
320	+	tIso2 = (combIso2) *50./(p->MomVtx(vtxCol->At(wVtxInd)->Position()).Pt());
321	+	RelIsoEcal=(ecalIso3-0.17*_tRho)/p->Et();
322	+	RelIsoHcal=(hcalIso4-0.17*_tRho)/p->Et();
323	+
324	+	//compute mva variables for v3
325	+	HoE = p->HadOverEm();
326	+	covIEtaIEta = p->CoviEtaiEta();
327	+	tIso1abs = combIso1;
328	+	tIso3abs = trackIso1;
329	+	tIso2abs = combIso2;
330	+	R9 = p->R9();
331	+
332	+	absIsoEcal=(ecalIso3-0.17*_tRho);
333	+	absIsoHcal=(hcalIso4-0.17*_tRho);
334	+	RelEMax=p->SCluster()->Seed()->EMax()/RawEnergy;
335	+	RelETop=p->SCluster()->Seed()->ETop()/RawEnergy;
336	+	RelEBottom=p->SCluster()->Seed()->EBottom()/RawEnergy;
337	+	RelELeft=p->SCluster()->Seed()->ELeft()/RawEnergy;
338	+	RelERight=p->SCluster()->Seed()->ERight()/RawEnergy;
339	+	RelE2x5Max=p->SCluster()->Seed()->E2x5Max()/RawEnergy;
340	+	RelE2x5Top=p->SCluster()->Seed()->E2x5Top()/RawEnergy;
341	+	RelE2x5Bottom=p->SCluster()->Seed()->E2x5Bottom()/RawEnergy;
342	+	RelE2x5Left=p->SCluster()->Seed()->E2x5Left()/RawEnergy;
343	+	RelE2x5Right=p->SCluster()->Seed()->E2x5Right()/RawEnergy;
344	+	RelE5x5=p->SCluster()->Seed()->E5x5()/RawEnergy;
345	+
346	+	EtaWidth=p->EtaWidth();
347	+	PhiWidth=p->PhiWidth();
348	+	CoviEtaiPhi=p->SCluster()->Seed()->CoviEtaiPhi();
349	+	CoviPhiiPhi=p->SCluster()->Seed()->CoviPhiiPhi();
350	+
351	+	RelPreshowerEnergy=p->SCluster()->PreshowerEnergy()/RawEnergy;
352	+	NVertexes=vtxCol->GetEntries();
353	+
354	+	//spectator variables
355	+	Pt_MVA=p->Pt();
356	+	ScEta_MVA=p->SCluster()->Eta();
357	+
358	+	//
359	+
360	+	isbarrel = (fabs(ScEta_MVA)<1.4442);
361	+
362	+	//variable 1201
363	+	myphoton_pfchargedisogood03=IsolationTools::PFChargedIsolation(p,vtx,0.3,0,fPFCands);
364	+	myphoton_pfchargedisobad03=IsolationTools::PFChargedIsolation(p,vtx,0.3,0,fPFCands,&wVtxInd,vtxCol);
365	+	myphoton_pfphotoniso03=IsolationTools::PFGammaIsolation(p,0.3,0,fPFCands);
366	+	myphoton_sieie=covIEtaIEta;
367	+	myphoton_sieip=CoviEtaiPhi;
368	+	myphoton_etawidth=EtaWidth;
369	+	myphoton_phiwidth=PhiWidth;
370	+	myphoton_r9=R9;
371	+	myphoton_s4ratio=p->S4Ratio();
372	+	myphoton_SCeta=ScEta_MVA;
373	+	event_rho= _tRho;
374	+
375	+	myphoton_ESEffSigmaRR=-99;
376	+
377	+	if(!isbarrel){
378	+	myphoton_ESEffSigmaRR=p->EffSigmaRR();
379	+	}
380	+
381	+	if (isbarrel) {
382	+	reader = fReaderBarrel;
383	+	}
384	+	else {
385	+	reader = fReaderEndcap;
386	+	}
387	+
388	+	assert(reader);
389	+
390	+	bdt = reader->EvaluateMVA("BDT method");
391	+
392	+	/* printf("HoE: %f\n",HoE);
393	+	printf("covIEtaIEta: %f\n",covIEtaIEta);
394	+	printf("tIso1abs: %f\n",tIso1abs);
395	+	printf("tIso3abs: %f\n",tIso3abs);
396	+	printf("tIso2abs: %f\n",tIso2abs);
397	+
398	+	printf("absIsoEcal: %f\n",absIsoEcal);
399	+	printf("absIsoHcal: %f\n",absIsoHcal);
400	+	printf("RelEMax: %f\n",RelEMax);
401	+	printf("RelETop: %f\n",RelETop);
402	+	printf("RelEBottom: %f\n",RelEBottom);
403	+	printf("RelELeft: %f\n",RelELeft);
404	+	printf("RelERight: %f\n",RelERight);
405	+	printf("RelE2x5Max: %f\n",RelE2x5Max);
406	+	printf("RelE2x5Top: %f\n",RelE2x5Top);
407	+	printf("RelE2x5Bottom: %f\n",RelE2x5Bottom);
408	+	printf("RelE2x5Left: %f\n",RelE2x5Left);
409	+	printf("RelE2x5Right;: %f\n",RelE2x5Right);
410	+	printf("RelE5x5: %f\n",RelE5x5);
411	+
412	+	printf("EtaWidth: %f\n",EtaWidth);
413	+	printf("PhiWidth: %f\n",PhiWidth);
414	+	printf("CoviEtaiPhi: %f\n",CoviEtaiPhi);
415	+	printf("CoviPhiiPhi: %f\n",CoviPhiiPhi);
416	+
417	+	if (!isbarrel) {
418	+	printf("RelPreshowerEnergy: %f\n",RelPreshowerEnergy);
419	+	}*/
420	+
421	+	return bdt;
422	+	}
423	+
424	+	Float_t MVATools::GetMVAbdtValue_2011(const Photon* p,const Vertex* vtx,const TrackCol* trackCol,const VertexCol* vtxCol,Double_t _tRho,const ElectronCol* els,Bool_t applyElectronVeto) {
425	+
426	+	//get the variables used to compute MVA variables
427	+	ecalIso3 = p->EcalRecHitIsoDr03();
428	+	ecalIso4 = p->EcalRecHitIsoDr04();
429	+	hcalIso4 = p->HcalTowerSumEtDr04();
430	+
431	+	wVtxInd = 0;
432	+
433	+	trackIso1 = IsolationTools::CiCTrackIsolation(p,vtx, 0.3, 0.02, 0.0, 0.0, 0.1, 1.0, trackCol, NULL, NULL, (!applyElectronVeto ? els : NULL) );//Question Ming:whyfPV->At(0) instead of selected vertex using ranking method?
434	+
435		// track iso worst vtx
436	<	trackIso2 = IsolationTools::CiCTrackIsolation(p,vtx, 0.4, 0.02, 0.0, 0.0, 0.1, 1.0, trackCol, &wVtxInd,vtxCol);
436	>	trackIso2 = IsolationTools::CiCTrackIsolation(p,vtx, 0.4, 0.02, 0.0, 0.0, 0.1, 1.0, trackCol, &wVtxInd,vtxCol, (!applyElectronVeto ? els : NULL) );
437
438		combIso1 = ecalIso3+hcalIso4+trackIso1 - 0.17*_tRho;
439		combIso2 = ecalIso4+hcalIso4+trackIso2 - 0.52*_tRho;
440
441		RawEnergy = p->SCluster()->RawEnergy();
442
443	+	ScEta = p->SCluster()->Eta();
444	+
445		//mva varialbes v1 and v2
446		tIso1 = (combIso1) *50./p->Et();
447	<	tIso3 = (trackIso3)*50./p->Et();
447	>	tIso3 = (trackIso1)*50./p->Et();
448		tIso2 = (combIso2) *50./(p->MomVtx(vtxCol->At(wVtxInd)->Position()).Pt());
449		RelIsoEcal=(ecalIso3-0.17*_tRho)/p->Et();
450		RelIsoHcal=(hcalIso4-0.17*_tRho)/p->Et();
#	Line 262 | Line 453 | Float_t MVATools::GetMVAbdtValue(const P
453		HoE = p->HadOverEm();
454		covIEtaIEta = p->CoviEtaiEta();
455		tIso1abs = combIso1;
456	<	tIso3abs = trackIso3;
456	>	tIso3abs = trackIso1;
457		tIso2abs = combIso2;
458		R9 = p->R9();
459
#	Line 280 | Line 471 | Float_t MVATools::GetMVAbdtValue(const P
471		RelE2x5Right=p->SCluster()->Seed()->E2x5Right()/RawEnergy;
472		RelE5x5=p->SCluster()->Seed()->E5x5()/RawEnergy;
473
474	<	EtaWidth=p->SCluster()->EtaWidth();
475	<	PhiWidth=p->SCluster()->PhiWidth();
474	>	EtaWidth=p->EtaWidth();
475	>	PhiWidth=p->PhiWidth();
476		CoviEtaiPhi=p->SCluster()->Seed()->CoviEtaiPhi();
477		CoviPhiiPhi=p->SCluster()->Seed()->CoviPhiiPhi();
478
#	Line 291 | Line 482 | Float_t MVATools::GetMVAbdtValue(const P
482		//spectator variables
483		Pt_MVA=p->Pt();
484		ScEta_MVA=p->SCluster()->Eta();
485	<
485	>
486	>	//
487	>
488		isbarrel = (fabs(ScEta_MVA)<1.4442);
489
490		if (isbarrel) {

Diff Legend

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