ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Utils/src/MVATools.cc
(Generate patch)

Comparing UserCode/MitPhysics/Utils/src/MVATools.cc (file contents):
Revision 1.4 by mingyang, Mon Dec 5 00:47:46 2011 UTC vs.
Revision 1.13 by mingyang, Tue May 29 21:20:12 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 137 | Line 141 | void MVATools::InitializeMVA(int Variabl
141          readers[i]->AddVariable( "RelPreshowerEnergy", &RelPreshowerEnergy );
142        }
143      }
144 +
145 +    if(VariableType==6){
146 +      readers[i]->AddVariable( "HoE", &HoE );
147 +      readers[i]->AddVariable( "covIEtaIEta", &covIEtaIEta );
148 +      readers[i]->AddVariable( "tIso1abs", &tIso1abs );
149 +      readers[i]->AddVariable( "tIso3abs", &tIso3abs );
150 +      readers[i]->AddVariable( "tIso2abs", &tIso2abs );
151 +      readers[i]->AddVariable( "R9", &R9 );
152 +      readers[i]->AddVariable( "absIsoEcal", &absIsoEcal );
153 +      readers[i]->AddVariable( "absIsoHcal", &absIsoHcal );
154 +      readers[i]->AddVariable( "RelE5x5", &RelE5x5 );
155 +      readers[i]->AddVariable( "EtaWidth", &EtaWidth );
156 +      readers[i]->AddVariable( "PhiWidth", &PhiWidth );
157 +      readers[i]->AddVariable( "CoviEtaiPhi", &CoviEtaiPhi );
158 +      readers[i]->AddVariable( "CoviPhiiPhi", &CoviPhiiPhi );
159 +      readers[i]->AddVariable( "NVertexes", &NVertexes );
160 +      if(i==0){
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 +    }
209 +    
210    }
211    
212    fReaderEndcap->BookMVA("BDT method",EndcapWeights);
# Line 147 | Line 217 | void MVATools::InitializeMVA(int Variabl
217    
218   }
219  
220 < Bool_t MVATools::PassMVASelection(const Photon* p,const Vertex* vtx,const TrackCol* trackCol,const VertexCol* vtxCol,Double_t _tRho,const ElectronCol* els,Float_t bdtCutBarrel, Float_t bdtCutEndcap) {
220 > 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) {
221    
222    //initilize the bool value
223    PassMVA=kFALSE;
224    
225 <  Float_t photon_bdt =  MVATools::GetMVAbdtValue(p,vtx,trackCol,vtxCol, _tRho,els);
225 >  Float_t photon_bdt =  MVATools::GetMVAbdtValue_2011(p,vtx,trackCol,vtxCol, _tRho, els, applyElectronVeto);
226    
227    if (isbarrel) {
228      if(bdt>bdtCutBarrel){
# Line 191 | Line 261 | Int_t MVATools::PassElectronVetoInt(cons
261    isbarrel = (fabs(ScEta_MVA)<1.4442);
262  
263    R9 = p->R9();
264 +  //R9 = p->E33()/p->SCluster()->RawEnergy();
265    
266    // check which category it is ...
267    _tCat = 1;
# Line 208 | Line 279 | Int_t MVATools::PassElectronVetoInt(cons
279  
280   //--------------------------------------------------------------------------------------------------
281  
282 < Float_t MVATools::GetMVAbdtValue(const Photon* p,const Vertex* vtx,const TrackCol* trackCol,const VertexCol* vtxCol,Double_t _tRho,const ElectronCol* els) {
282 > Float_t MVATools::GetMVAbdtValue_2012_globe(const Photon* p,const Vertex* vtx,const TrackCol* trackCol,const VertexCol* vtxCol,Double_t _tRho, const PFCandidateCol *fPFCands,Bool_t applyNewScale,const ElectronCol* els,Bool_t applyElectronVeto) {
283    
284    //get the variables used to compute MVA variables
285    ecalIso3 = p->EcalRecHitIsoDr03();
# Line 217 | Line 288 | Float_t MVATools::GetMVAbdtValue(const P
288    
289    wVtxInd = 0;
290    
291 <  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?
291 >  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?
292 >    
293 >  // track iso worst vtx
294 >  trackIso2 = IsolationTools::CiCTrackIsolation(p,vtx, 0.4, 0.02, 0.0, 0.0, 0.1, 1.0, trackCol, &wVtxInd,vtxCol, (!applyElectronVeto ? els : NULL) );
295 >  
296 >  combIso1 = ecalIso3+hcalIso4+trackIso1 - 0.17*_tRho;
297 >  combIso2 = ecalIso4+hcalIso4+trackIso2 - 0.52*_tRho;
298 >  
299 >  RawEnergy = p->SCluster()->RawEnergy();
300 >  
301 >  ScEta = p->SCluster()->Eta();
302    
303 <  // track iso only
304 <  trackIso3 = IsolationTools::CiCTrackIsolation(p,vtx, 0.3, 0.02, 0.0, 0.0, 0.1, 1.0, trackCol);
303 >  //mva varialbes v1 and v2
304 >  tIso1 = (combIso1) *50./p->Et();
305 >  tIso3 = (trackIso1)*50./p->Et();
306 >  tIso2 = (combIso2) *50./(p->MomVtx(vtxCol->At(wVtxInd)->Position()).Pt());
307 >  RelIsoEcal=(ecalIso3-0.17*_tRho)/p->Et();
308 >  RelIsoHcal=(hcalIso4-0.17*_tRho)/p->Et();
309 >
310 >  //compute mva variables for v3
311 >  HoE = p->HadOverEm();
312 >  covIEtaIEta = p->CoviEtaiEta();
313 >  tIso1abs = combIso1;
314 >  tIso3abs = trackIso1;
315 >  tIso2abs = combIso2;
316 >  R9 = p->R9();
317 >
318 >  absIsoEcal=(ecalIso3-0.17*_tRho);
319 >  absIsoHcal=(hcalIso4-0.17*_tRho);
320 >  RelEMax=p->SCluster()->Seed()->EMax()/RawEnergy;
321 >  RelETop=p->SCluster()->Seed()->ETop()/RawEnergy;
322 >  RelEBottom=p->SCluster()->Seed()->EBottom()/RawEnergy;
323 >  RelELeft=p->SCluster()->Seed()->ELeft()/RawEnergy;
324 >  RelERight=p->SCluster()->Seed()->ERight()/RawEnergy;
325 >  RelE2x5Max=p->SCluster()->Seed()->E2x5Max()/RawEnergy;
326 >  RelE2x5Top=p->SCluster()->Seed()->E2x5Top()/RawEnergy;
327 >  RelE2x5Bottom=p->SCluster()->Seed()->E2x5Bottom()/RawEnergy;
328 >  RelE2x5Left=p->SCluster()->Seed()->E2x5Left()/RawEnergy;
329 >  RelE2x5Right=p->SCluster()->Seed()->E2x5Right()/RawEnergy;
330 >  RelE5x5=p->SCluster()->Seed()->E5x5()/RawEnergy;
331    
332 +  EtaWidth=p->EtaWidth();
333 +  PhiWidth=p->PhiWidth();
334 +  CoviEtaiPhi=p->SCluster()->Seed()->CoviEtaiPhi();
335 +  CoviPhiiPhi=p->SCluster()->Seed()->CoviPhiiPhi();
336 +
337 +  RelPreshowerEnergy=p->SCluster()->PreshowerEnergy()/RawEnergy;
338 +  NVertexes=vtxCol->GetEntries();
339 +  
340 +  //spectator variables
341 +  Pt_MVA=p->Pt();
342 +  ScEta_MVA=p->SCluster()->Eta();
343 +
344 +  //
345 +
346 +  isbarrel = (fabs(ScEta_MVA)<1.4442);
347 +  
348 +  //variable 1201
349 +  myphoton_pfchargedisogood03=IsolationTools::PFChargedIsolation(p,vtx,0.3,0,fPFCands);
350 +  myphoton_pfchargedisobad03=IsolationTools::PFChargedIsolation(p,vtx,0.3,0,fPFCands,&wVtxInd,vtxCol);
351 +  myphoton_pfphotoniso03=IsolationTools::PFGammaIsolation(p,0.3,0,fPFCands);
352 +  myphoton_sieie=covIEtaIEta;
353 +  myphoton_sieip=CoviEtaiPhi;
354 +  myphoton_etawidth=EtaWidth;
355 +  myphoton_phiwidth=PhiWidth;
356 +  myphoton_r9=R9;
357 +  myphoton_s4ratio=p->SCluster()->Seed()->E2x2()/p->SCluster()->Seed()->E5x5();
358 +  myphoton_SCeta=ScEta_MVA;
359 +  event_rho= _tRho;
360 +
361 +  myphoton_ESEffSigmaRR=-99;
362 +
363 +  if(!isbarrel){
364 +    myphoton_ESEffSigmaRR=float(sqrt((p->SCluster()->PsEffWidthSigmaXX())*(p->SCluster()->PsEffWidthSigmaXX())+(p->SCluster()->PsEffWidthSigmaYY())*(p->SCluster()->PsEffWidthSigmaYY())));
365 +  }
366 +
367 +  if(applyNewScale){
368 +    if(isbarrel){
369 +      myphoton_s4ratio=1.0055*myphoton_s4ratio;
370 +    }
371 +    if(!isbarrel){
372 +      myphoton_s4ratio=1.0085*myphoton_s4ratio;
373 +      myphoton_ESEffSigmaRR=1.04*myphoton_ESEffSigmaRR;
374 +    }
375 +  }
376 +  
377 +  if (isbarrel) {
378 +    reader = fReaderBarrel;
379 +  }
380 +  else {
381 +    reader = fReaderEndcap;
382 +  }
383 +  
384 +  assert(reader);
385 +
386 +  bdt = reader->EvaluateMVA("BDT method");
387 +
388 +  /* printf("HoE: %f\n",HoE);
389 +  printf("covIEtaIEta: %f\n",covIEtaIEta);
390 +  printf("tIso1abs: %f\n",tIso1abs);
391 +  printf("tIso3abs: %f\n",tIso3abs);
392 +  printf("tIso2abs: %f\n",tIso2abs);
393 +  
394 +  printf("absIsoEcal: %f\n",absIsoEcal);
395 +  printf("absIsoHcal: %f\n",absIsoHcal);
396 +  printf("RelEMax: %f\n",RelEMax);
397 +  printf("RelETop: %f\n",RelETop);
398 +  printf("RelEBottom: %f\n",RelEBottom);
399 +  printf("RelELeft: %f\n",RelELeft);
400 +  printf("RelERight: %f\n",RelERight);
401 +  printf("RelE2x5Max: %f\n",RelE2x5Max);
402 +  printf("RelE2x5Top: %f\n",RelE2x5Top);
403 +  printf("RelE2x5Bottom: %f\n",RelE2x5Bottom);
404 +  printf("RelE2x5Left: %f\n",RelE2x5Left);
405 +  printf("RelE2x5Right;: %f\n",RelE2x5Right);
406 +  printf("RelE5x5: %f\n",RelE5x5);
407 +  
408 +  printf("EtaWidth: %f\n",EtaWidth);
409 +  printf("PhiWidth: %f\n",PhiWidth);
410 +  printf("CoviEtaiPhi: %f\n",CoviEtaiPhi);
411 +  printf("CoviPhiiPhi: %f\n",CoviPhiiPhi);
412 +  
413 +  if (!isbarrel) {
414 +    printf("RelPreshowerEnergy: %f\n",RelPreshowerEnergy);
415 +    }*/
416 +  
417 +  return bdt;
418 + }
419 +
420 + 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) {
421 +  
422 +  //get the variables used to compute MVA variables
423 +  ecalIso3 = p->EcalRecHitIsoDr03();
424 +  ecalIso4 = p->EcalRecHitIsoDr04();
425 +  hcalIso4 = p->HcalTowerSumEtDr04();
426 +  
427 +  wVtxInd = 0;
428 +  
429 +  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?
430 +    
431    // track iso worst vtx
432 <  trackIso2 = IsolationTools::CiCTrackIsolation(p,vtx, 0.4, 0.02, 0.0, 0.0, 0.1, 1.0, trackCol, &wVtxInd,vtxCol);
432 >  trackIso2 = IsolationTools::CiCTrackIsolation(p,vtx, 0.4, 0.02, 0.0, 0.0, 0.1, 1.0, trackCol, &wVtxInd,vtxCol, (!applyElectronVeto ? els : NULL) );
433    
434    combIso1 = ecalIso3+hcalIso4+trackIso1 - 0.17*_tRho;
435    combIso2 = ecalIso4+hcalIso4+trackIso2 - 0.52*_tRho;
436    
437    RawEnergy = p->SCluster()->RawEnergy();
438    
439 <  dRTrack = PhotonTools::ElectronVetoCiC(p, els);
440 <
439 >  ScEta = p->SCluster()->Eta();
440 >  
441    //mva varialbes v1 and v2
442    tIso1 = (combIso1) *50./p->Et();
443 <  tIso3 = (trackIso3)*50./p->Et();
443 >  tIso3 = (trackIso1)*50./p->Et();
444    tIso2 = (combIso2) *50./(p->MomVtx(vtxCol->At(wVtxInd)->Position()).Pt());
445    RelIsoEcal=(ecalIso3-0.17*_tRho)/p->Et();
446    RelIsoHcal=(hcalIso4-0.17*_tRho)/p->Et();
# Line 243 | Line 449 | Float_t MVATools::GetMVAbdtValue(const P
449    HoE = p->HadOverEm();
450    covIEtaIEta = p->CoviEtaiEta();
451    tIso1abs = combIso1;
452 <  tIso3abs = trackIso3;
452 >  tIso3abs = trackIso1;
453    tIso2abs = combIso2;
454    R9 = p->R9();
455  
# Line 261 | Line 467 | Float_t MVATools::GetMVAbdtValue(const P
467    RelE2x5Right=p->SCluster()->Seed()->E2x5Right()/RawEnergy;
468    RelE5x5=p->SCluster()->Seed()->E5x5()/RawEnergy;
469    
470 <  EtaWidth=p->SCluster()->EtaWidth();
471 <  PhiWidth=p->SCluster()->PhiWidth();
470 >  EtaWidth=p->EtaWidth();
471 >  PhiWidth=p->PhiWidth();
472    CoviEtaiPhi=p->SCluster()->Seed()->CoviEtaiPhi();
473    CoviPhiiPhi=p->SCluster()->Seed()->CoviPhiiPhi();
474  
# Line 272 | Line 478 | Float_t MVATools::GetMVAbdtValue(const P
478    //spectator variables
479    Pt_MVA=p->Pt();
480    ScEta_MVA=p->SCluster()->Eta();
481 <  
481 >
482 >  //
483 >
484    isbarrel = (fabs(ScEta_MVA)<1.4442);
485    
486    if (isbarrel) {
# Line 286 | Line 494 | Float_t MVATools::GetMVAbdtValue(const P
494  
495    bdt = reader->EvaluateMVA("BDT method");
496  
497 <  printf("HoE: %f\n",HoE);
497 >  /* printf("HoE: %f\n",HoE);
498    printf("covIEtaIEta: %f\n",covIEtaIEta);
499    printf("tIso1abs: %f\n",tIso1abs);
500    printf("tIso3abs: %f\n",tIso3abs);
# Line 313 | Line 521 | Float_t MVATools::GetMVAbdtValue(const P
521    
522    if (!isbarrel) {
523      printf("RelPreshowerEnergy: %f\n",RelPreshowerEnergy);
524 <  }
524 >    }*/
525    
526    return bdt;
527   }

Diff Legend

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