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

Comparing UserCode/MitPhysics/Utils/src/JetIDMVA.cc (file contents):
Revision 1.9 by ceballos, Sat May 12 06:17:43 2012 UTC vs.
Revision 1.14 by bendavid, Tue May 29 17:09:02 2012 UTC

# Line 35 | Line 35 | JetIDMVA::JetIDMVA() :
35    fFrac02   (0),
36    fFrac03   (0),
37    fFrac04   (0),
38 <  fFrac05   (0)
38 >  fFrac05   (0),
39 >  fDR2Mean  (0)
40   {    
41    fReader = 0;
42   }
43   //--------------------------------------------------------------------------------------------------
44   JetIDMVA::~JetIDMVA() {
45  
46 <  fReader = 0;
46 >  delete fReader;
47 >  delete fLowPtReader;
48   }
49  
50   //--------------------------------------------------------------------------------------------------
# Line 56 | Line 58 | void JetIDMVA::Initialize( JetIDMVA::Cut
58    fIsInitialized = kTRUE;
59    fType          = iType;
60    fCutType       = iCutType;
61 +  
62 +  std::string lCutId = "JetIdParams";
63 +  if(fType == k42)  lCutId = "PuJetIdOptMVA_wp";
64 +  if(fType == k52)  lCutId = "full_5x_wp";
65 +  if(fType == kCut) lCutId = "PuJetIdCutBased_wp";
66 +  //Load Cut Matrix
67 +  edm::ParameterSet lDefConfig = edm::readPSetsFrom(iCutFileName.Data())->getParameter<edm::ParameterSet>("JetIdParams");
68 +  edm::ParameterSet lConfig    = edm::readPSetsFrom(iCutFileName.Data())->getParameter<edm::ParameterSet>(lCutId);
69 +  std::string lCutType = "Tight";
70 +  if(fCutType == kMedium) lCutType = "Medium";
71 +  if(fCutType == kLoose ) lCutType = "Loose";
72 +  if(fCutType == kMET   ) lCutType = "MET";
73 +  if(fType != kCut) {
74 +    std::string lLowPtCut = "MET";
75 +    std::vector<double> lPt010  = lDefConfig.getParameter<std::vector<double> >(("Pt010_" +lLowPtCut).c_str());
76 +    std::vector<double> lPt1020 = lConfig.getParameter<std::vector<double> >(("Pt1020_"+lCutType).c_str());
77 +    std::vector<double> lPt2030 = lConfig.getParameter<std::vector<double> >(("Pt2030_"+lCutType).c_str());
78 +    std::vector<double> lPt3050 = lConfig.getParameter<std::vector<double> >(("Pt3050_"+lCutType).c_str());
79 +    for(int i0 = 0; i0 < 4; i0++) fMVACut[0][i0] = lPt010 [i0];
80 +    for(int i0 = 0; i0 < 4; i0++) fMVACut[1][i0] = lPt1020[i0];
81 +    for(int i0 = 0; i0 < 4; i0++) fMVACut[2][i0] = lPt2030[i0];
82 +    for(int i0 = 0; i0 < 4; i0++) fMVACut[3][i0] = lPt3050[i0];
83 +  }
84 +  if(fType == kCut) {
85 +    for(int i0 = 0; i0 < 2; i0++) {
86 +      std::string lFullCutType = lCutType;
87 +      if(i0 == 0) lFullCutType = "BetaStar"+ lCutType;
88 +      if(i0 == 1) lFullCutType = "RMS"     + lCutType;
89 +      std::vector<double> pt010  = lConfig.getParameter<std::vector<double> >(("Pt010_" +lFullCutType).c_str());
90 +      std::vector<double> pt1020 = lConfig.getParameter<std::vector<double> >(("Pt1020_"+lFullCutType).c_str());
91 +      std::vector<double> pt2030 = lConfig.getParameter<std::vector<double> >(("Pt2030_"+lFullCutType).c_str());
92 +      std::vector<double> pt3050 = lConfig.getParameter<std::vector<double> >(("Pt3050_"+lFullCutType).c_str());
93 +      if(i0 == 0) {
94 +        for(int i2 = 0; i2 < 4; i2++) fBetaStarCut[0][i2] = pt010 [i2];
95 +        for(int i2 = 0; i2 < 4; i2++) fBetaStarCut[1][i2] = pt1020[i2];
96 +        for(int i2 = 0; i2 < 4; i2++) fBetaStarCut[2][i2] = pt2030[i2];
97 +        for(int i2 = 0; i2 < 4; i2++) fBetaStarCut[3][i2] = pt3050[i2];
98 +      }
99 +      if(i0 == 1) {
100 +        for(int i2 = 0; i2 < 4; i2++) fRMSCut[0][i2] = pt010 [i2];
101 +        for(int i2 = 0; i2 < 4; i2++) fRMSCut[1][i2] = pt1020[i2];
102 +        for(int i2 = 0; i2 < 4; i2++) fRMSCut[2][i2] = pt2030[i2];
103 +        for(int i2 = 0; i2 < 4; i2++) fRMSCut[3][i2] = pt3050[i2];
104 +      }
105 +    }
106 +    return;
107 +  }
108 +  
109 +  
110 +  fLowPtReader   = 0;
111 +  fLowPtReader   = new TMVA::Reader( "!Color:!Silent:Error" );  
112 +  fLowPtReader->AddVariable( "nvtx"     , &fNVtx      );
113 +  fLowPtReader->AddVariable( "jetPt"    , &fJPt1      );  
114 +  fLowPtReader->AddVariable( "jetEta"   , &fJEta1     );
115 +  fLowPtReader->AddVariable( "jetPhi"   , &fJPhi1     );            
116 +  fLowPtReader->AddVariable( "dZ"       , &fJDZ1      );
117 +  fLowPtReader->AddVariable( "d0"       , &fJD01      );
118 +  fLowPtReader->AddVariable( "beta"     , &fBeta      );
119 +  fLowPtReader->AddVariable( "betaStar" , &fBetaStar  );
120 +  fLowPtReader->AddVariable( "nCharged" , &fNCharged  );
121 +  fLowPtReader->AddVariable( "nNeutrals", &fNNeutrals );
122 +  fLowPtReader->AddVariable( "dRMean"   , &fDRMean    );
123 +  fLowPtReader->AddVariable( "frac01"   , &fFrac01    );
124 +  fLowPtReader->AddVariable( "frac02"   , &fFrac02    );
125 +  fLowPtReader->AddVariable( "frac03"   , &fFrac03    );
126 +  fLowPtReader->AddVariable( "frac04"   , &fFrac04    );
127 +  fLowPtReader->AddVariable( "frac05"   , &fFrac05    );
128 +  
129    fReader        = 0;
130    fReader        = new TMVA::Reader( "!Color:!Silent:Error" );  
131    if (fType == kBaseline) {
# Line 76 | Line 146 | void JetIDMVA::Initialize( JetIDMVA::Cut
146      fReader->AddVariable( "frac04"   , &fFrac04    );
147      fReader->AddVariable( "frac05"   , &fFrac05    );
148    }
149 <  fReader->BookMVA(fLowPtMethodName  , iLowPtWeights );
149 >  if (fType == k42) {
150 >    fReader->AddVariable( "frac01"   , &fFrac01    );
151 >    fReader->AddVariable( "frac02"   , &fFrac02    );
152 >    fReader->AddVariable( "frac03"   , &fFrac03    );
153 >    fReader->AddVariable( "frac04"   , &fFrac04    );
154 >    fReader->AddVariable( "frac05"   , &fFrac05    );
155 >    fReader->AddVariable( "nvtx"     , &fNVtx      );
156 >    fReader->AddVariable( "nNeutrals", &fNNeutrals );
157 >    fReader->AddVariable( "beta"     , &fBeta      );
158 >    fReader->AddVariable( "betaStar" , &fBetaStar  );
159 >    fReader->AddVariable( "dZ"       , &fJDZ1      );
160 >    fReader->AddVariable( "nCharged" , &fNCharged  );
161 >    fReader->AddSpectator( "jetPt"    , &fJPt1      );  
162 >    fReader->AddSpectator( "jetEta"   , &fJEta1     );
163 >  }
164 >  if (fType == k52) {
165 >    fReader->AddVariable( "frac01"   , &fFrac01    );
166 >    fReader->AddVariable( "frac02"   , &fFrac02    );
167 >    fReader->AddVariable( "frac03"   , &fFrac03    );
168 >    fReader->AddVariable( "frac04"   , &fFrac04    );
169 >    fReader->AddVariable( "frac05"   , &fFrac05    );
170 >    fReader->AddVariable( "dR2Mean"  , &fDR2Mean   );
171 >    fReader->AddVariable( "nvtx"     , &fNVtx      );
172 >    fReader->AddVariable( "nNeutrals", &fNNeutrals );
173 >    fReader->AddVariable( "beta"     , &fBeta      );
174 >    fReader->AddVariable( "betaStar" , &fBetaStar  );
175 >    fReader->AddVariable( "dZ"       , &fJDZ1      );
176 >    fReader->AddVariable( "nCharged" , &fNCharged  );
177 >    fReader->AddSpectator( "jetPt"    , &fJPt1      );  
178 >    fReader->AddSpectator( "jetEta"   , &fJEta1     );
179 >  }
180 >
181 >  fLowPtReader->BookMVA(fLowPtMethodName  , iLowPtWeights );
182    fReader->BookMVA(fHighPtMethodName , iHighPtWeights );
183    std::cout << "Jet ID MVA Initialization\n";
184    std::cout << "MethodName : " << fLowPtMethodName << " , type == " << fType << std::endl;
185  
84 –  //Load Cut Matrix
85 –  edm::ParameterSet lConfig = edm::readPSetsFrom(iCutFileName.Data())->getParameter<edm::ParameterSet>("JetIdParams");
86 –  std::string lCutType = "Tight";
87 –  if(fCutType == kMedium) lCutType = "Medium";
88 –  if(fCutType == kLoose ) lCutType = "Loose";
89 –  if(fCutType == kMET   ) lCutType = "MET";
90 –  std::vector<double> lPt010  = lConfig.getParameter<std::vector<double> >(("Pt010_" +lCutType).c_str());
91 –  std::vector<double> lPt1020 = lConfig.getParameter<std::vector<double> >(("Pt1020_"+lCutType).c_str());
92 –  std::vector<double> lPt2030 = lConfig.getParameter<std::vector<double> >(("Pt2030_"+lCutType).c_str());
93 –  std::vector<double> lPt3050 = lConfig.getParameter<std::vector<double> >(("Pt3050_"+lCutType).c_str());
94 –  for(int i0 = 0; i0 < 4; i0++) fMVACut[0][i0] = lPt010 [i0];
95 –  for(int i0 = 0; i0 < 4; i0++) fMVACut[1][i0] = lPt1020[i0];
96 –  for(int i0 = 0; i0 < 4; i0++) fMVACut[2][i0] = lPt2030[i0];
97 –  for(int i0 = 0; i0 < 4; i0++) fMVACut[3][i0] = lPt3050[i0];
98 –  //std::cout << " Working Points : << " << std::endl;
99 –  //for(int i0 = 0; i0 < 4; i0++) for(int i1 = 0; i1 < 4; i1++)
100 –  //  std::cout << " ==> " << i0 << " -- " << i1 << " -- " << fMVACut[i0][i1] << std::endl;
186   }
187  
188   //--------------------------------------------------------------------------------------------------
# Line 117 | Line 202 | Double_t JetIDMVA::MVAValue(
202                              Float_t iFrac02  ,
203                              Float_t iFrac03  ,
204                              Float_t iFrac04  ,
205 <                            Float_t iFrac05  
205 >                            Float_t iFrac05  ,
206 >                            Float_t iDR2Mean  
207                              ){
208    
209    if(!fIsInitialized) {
# Line 141 | Line 227 | Double_t JetIDMVA::MVAValue(
227    fFrac03    = iFrac03;
228    fFrac04    = iFrac04;
229    fFrac05    = iFrac05;
230 +  fDR2Mean   = iDR2Mean;
231  
232    Double_t lMVA = -9999;  
233 <  if(iJPt1 < 10) lMVA = fReader->EvaluateMVA( fLowPtMethodName  );
233 >  if(iJPt1 < 10) lMVA = fLowPtReader->EvaluateMVA( fLowPtMethodName  );
234    if(iJPt1 > 10) lMVA = fReader->EvaluateMVA( fHighPtMethodName );
235  
236    return lMVA;
# Line 151 | Line 238 | Double_t JetIDMVA::MVAValue(
238   //--------------------------------------------------------------------------------------------------
239   Bool_t JetIDMVA::pass(const PFJet *iJet,const Vertex *iVertex,const VertexCol *iVertices,
240                        FactorizedJetCorrector *iJetCorrector,
241 <                      const PileupEnergyDensityCol *iPileupEnergyDensity) {
241 >                      const PileupEnergyDensityCol *iPileupEnergyDensity,
242 >                      RhoUtilities::RhoType type) {
243    
244    if(!JetTools::passPFLooseId(iJet))                 return false;
245    if(fabs(iJet->Eta()) > 4.99)                       return false;
246    
247    double lMVA = MVAValue   (iJet,iVertex,iVertices,iJetCorrector,iPileupEnergyDensity);
248 <  double lPt  = correctedPt(iJet,                  iJetCorrector,iPileupEnergyDensity);
248 >  double lPt  = correctedPt(iJet,                  iJetCorrector,iPileupEnergyDensity,type);
249    if(lPt < fJetPtMin)                         return false;
250    
251    int lPtId = 0;
# Line 173 | Line 261 | Bool_t JetIDMVA::pass(const PFJet *iJet,
261    double lMVACut = fMVACut[lPtId][lEtaId];
262    if(lMVA < lMVACut) return false;
263    return true;
264 <   //if( fabs(JetTools::impactParameter(iJet,iVertex,true)) < 0.2) return true;
265 <  //if(correctedPt(iJet,iJetCorrector,iPileupEnergyDensity) < fJetPtMin && iJet->TrackCountingHighEffBJetTagsDisc() == -100) return false;
266 <  //This line is a bug in the Met training
267 <  //if(lMVA < -0.8)                            return false;
268 <  //if(lMVA < -0.5 && fabs(iJet->Eta()) > 3.0) return false;
264 > }
265 > //--------------------------------------------------------------------------------------------------
266 > Bool_t JetIDMVA::passCut(const PFJet *iJet,const Vertex *iVertex,const VertexCol *iVertices) {
267 >  if(!JetTools::passPFLooseId(iJet))                 return false;
268 >  if(iJet->Pt()        < fJetPtMin) return false;
269 >  if(fabs(iJet->Eta()) > 4.99)      return false;
270 >  //if(fType == kCut) passCut(iJet,iVertex,iVertices);
271 >
272 >  double lPt = iJet->Pt();
273 >  int lPtId = 0;
274 >  if(lPt > 10 && lPt < 20) lPtId = 1;
275 >  if(lPt > 20 && lPt < 30) lPtId = 2;
276 >  if(lPt > 30                   ) lPtId = 3;
277 >  
278 >  int lEtaId = 0;
279 >  if(fabs(iJet->Eta()) > 2.5  && fabs(iJet->Eta()) < 2.75) lEtaId = 1;
280 >  if(fabs(iJet->Eta()) > 2.75 && fabs(iJet->Eta()) < 3.0 ) lEtaId = 2;
281 >  if(fabs(iJet->Eta()) > 3.0  && fabs(iJet->Eta()) < 5.0 ) lEtaId = 3;
282 >  float betaStarModified = JetTools::betaStarClassic(iJet,iVertex,iVertices)/log(iVertices ->GetEntries()-0.64);
283 >  float dR2Mean          = JetTools::dR2Mean(iJet,-1);
284 >  
285 >  if(betaStarModified < fBetaStarCut[lPtId][lEtaId] &&
286 >     dR2Mean          < fRMSCut     [lPtId][lEtaId]
287 >     ) return true;
288 >  
289 >  return false;
290   }
291   //--------------------------------------------------------------------------------------------------
292   Bool_t JetIDMVA::pass(const PFJet *iJet,const Vertex *iVertex,const VertexCol *iVertices) {
293    if(!JetTools::passPFLooseId(iJet))                 return false;
294    if(iJet->Pt()        < fJetPtMin) return false;
295    if(fabs(iJet->Eta()) > 4.99)      return false;
296 +  if(fType == kCut) return passCut(iJet,iVertex,iVertices);
297    double lMVA = MVAValue(iJet,iVertex,iVertices);
298    
299    int lPtId = 0;
# Line 229 | Line 339 | Double_t JetIDMVA::MVAValue(const PFJet
339    fNNeutrals = iJet->NeutralMultiplicity();
340  
341    fDRMean    = JetTools::dRMean(iJet,-1);
342 +  fDR2Mean   = JetTools::dR2Mean(iJet,-1);
343    fFrac01    = JetTools::frac  (iJet,0.1,0. ,-1);
344    fFrac02    = JetTools::frac  (iJet,0.2,0.1,-1);
345    fFrac03    = JetTools::frac  (iJet,0.3,0.2,-1);
# Line 236 | Line 347 | Double_t JetIDMVA::MVAValue(const PFJet
347    fFrac05    = JetTools::frac  (iJet,0.5,0.4,-1);
348  
349    double lMVA = 0;
350 <  if(fJPt1 < 10) lMVA = fReader->EvaluateMVA( fLowPtMethodName  );
350 >  if(fJPt1 < 10) lMVA = fLowPtReader->EvaluateMVA( fLowPtMethodName  );
351    if(fJPt1 > 10) lMVA = fReader->EvaluateMVA( fHighPtMethodName );
352    if (printDebug == kTRUE) {
353      std::cout << "Debug Jet MVA: "
# Line 255 | Line 366 | Double_t JetIDMVA::MVAValue(const PFJet
366                << fFrac02    << " "
367                << fFrac03    << " "
368                << fFrac04    << " "
369 <              << fFrac05    
369 >              << fFrac05    << " "
370 >              << fDRMean    
371                << " === : === "
372                << lMVA << " "    
373                << std::endl;
# Line 284 | Line 396 | Double_t JetIDMVA::MVAValue(const PFJet
396    fNCharged  = iJet->ChargedMultiplicity();
397    fNNeutrals = iJet->NeutralMultiplicity();
398  
399 <  fDRMean    = JetTools::dRMean(iJet,-1);
400 <  fFrac01    = JetTools::frac  (iJet,0.1,0. ,-1);
401 <  fFrac02    = JetTools::frac  (iJet,0.2,0.1,-1);
402 <  fFrac03    = JetTools::frac  (iJet,0.3,0.2,-1);
403 <  fFrac04    = JetTools::frac  (iJet,0.4,0.3,-1);
404 <  fFrac05    = JetTools::frac  (iJet,0.5,0.4,-1);
399 >  fDRMean    = JetTools::dRMean (iJet,-1);
400 >  fDR2Mean   = JetTools::dR2Mean(iJet,-1);
401 >  fFrac01    = JetTools::frac   (iJet,0.1,0. ,-1);
402 >  fFrac02    = JetTools::frac   (iJet,0.2,0.1,-1);
403 >  fFrac03    = JetTools::frac   (iJet,0.3,0.2,-1);
404 >  fFrac04    = JetTools::frac   (iJet,0.4,0.3,-1);
405 >  fFrac05    = JetTools::frac   (iJet,0.5,0.4,-1);
406  
407    double lMVA = 0;
408 <  if(fJPt1 < 10) lMVA = fReader->EvaluateMVA( fLowPtMethodName  );
409 <  if(fJPt1 > 10) lMVA = fReader->EvaluateMVA( fHighPtMethodName );
408 >  if(fJPt1 < 10) lMVA = fLowPtReader->EvaluateMVA( fLowPtMethodName  );
409 >  if(fJPt1 > 10) lMVA = fReader     ->EvaluateMVA( fHighPtMethodName );
410    
411    if (printDebug == kTRUE) {
412      std::cout << "Debug Jet MVA: "
# Line 312 | Line 425 | Double_t JetIDMVA::MVAValue(const PFJet
425                << fFrac02    << " "
426                << fFrac03    << " "
427                << fFrac04    << " "
428 <              << fFrac05    
428 >              << fFrac05    << " "
429 >              << fDRMean    
430                << " === : === "
431                << lMVA << " "    
432                << std::endl;
# Line 320 | Line 434 | Double_t JetIDMVA::MVAValue(const PFJet
434  
435    return lMVA;
436   }
437 < Double_t JetIDMVA::correctedPt(const PFJet *iJet, FactorizedJetCorrector *iJetCorrector,const PileupEnergyDensityCol *iPUEnergyDensity) {
437 > Double_t JetIDMVA::correctedPt(const PFJet *iJet, FactorizedJetCorrector *iJetCorrector,
438 >                               const PileupEnergyDensityCol *iPUEnergyDensity,
439 >                               RhoUtilities::RhoType type) {
440 >  Double_t Rho = 0.0;
441 >  switch(type) {
442 >  case RhoUtilities::MIT_RHO_VORONOI_HIGH_ETA:
443 >    Rho = iPUEnergyDensity->At(0)->Rho();
444 >    break;
445 >  case RhoUtilities::MIT_RHO_VORONOI_LOW_ETA:
446 >    Rho = iPUEnergyDensity->At(0)->RhoLowEta();
447 >    break;
448 >  case RhoUtilities::MIT_RHO_RANDOM_HIGH_ETA:
449 >    Rho = iPUEnergyDensity->At(0)->RhoRandom();
450 >    break;
451 >  case RhoUtilities::MIT_RHO_RANDOM_LOW_ETA:
452 >    Rho = iPUEnergyDensity->At(0)->RhoRandomLowEta();
453 >    break;
454 >  case RhoUtilities::CMS_RHO_RHOKT6PFJETS:
455 >    Rho = iPUEnergyDensity->At(0)->RhoKt6PFJets();
456 >    break;
457 >  default:
458 >    // use the old default
459 >    Rho = iPUEnergyDensity->At(0)->Rho();
460 >    break;
461 >  }
462 >    
463    const FourVectorM rawMom = iJet->RawMom();
464    iJetCorrector->setJetEta(rawMom.Eta());
465    iJetCorrector->setJetPt (rawMom.Pt());
466    iJetCorrector->setJetPhi(rawMom.Phi());
467    iJetCorrector->setJetE  (rawMom.E());
468 <  iJetCorrector->setRho   (iPUEnergyDensity->At(0)->RhoHighEta());
468 >  iJetCorrector->setRho   (Rho);
469    iJetCorrector->setJetA  (iJet->JetArea());
470    iJetCorrector->setJetEMF(-99.0);    
471    Double_t correction = iJetCorrector->getCorrection();

Diff Legend

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