Utils/src/JetIDMVA.cc

#include "FWCore/PythonParameterSet/interface/MakeParameterSets.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "MitPhysics/Utils/interface/JetIDMVA.h"
#include "MitPhysics/Utils/interface/JetTools.h"
#include "MitAna/DataTree/interface/StableData.h"
#include <TFile.h>
#include <TRandom3.h>
#include "TMVA/Tools.h"
#include "TMVA/Reader.h"


ClassImp(mithep::JetIDMVA)

using namespace mithep;

//--------------------------------------------------------------------------------------------------
JetIDMVA::JetIDMVA() :
  fJetPtMin(0.)  , //We need to lower this
  fDZCut   (0.2),
  fLowPtMethodName ("JetIDMVALowPt" ),
  fHighPtMethodName("JetIDMVAHighPt"),
  fIsInitialized(kFALSE),
  fNVtx     (0),
  fJPt1     (0),
  fJEta1    (0),
  fJPhi1    (0),
  fJD01     (0),
  fJDZ1     (0),
  fBeta     (0),
  fBetaStar (0),
  fNCharged (0),
  fNNeutrals(0),
  fDRMean   (0),
  fFrac01   (0),
  fFrac02   (0),
  fFrac03   (0),
  fFrac04   (0),
  fFrac05   (0),
  fDR2Mean  (0)
{    
  fReader = 0;
}
//--------------------------------------------------------------------------------------------------
JetIDMVA::~JetIDMVA() {

  delete fReader;
  delete fLowPtReader;
}

//--------------------------------------------------------------------------------------------------
void JetIDMVA::Initialize( JetIDMVA::CutType iCutType,
                           TString iLowPtWeights,
                           TString iHighPtWeights, 
                           JetIDMVA::MVAType iType,
                           TString iCutFileName) 
{ 
  
  fIsInitialized = kTRUE;
  fType          = iType;
  fCutType       = iCutType;
  fLowPtReader   = 0;
  fLowPtReader   = new TMVA::Reader( "!Color:!Silent:Error" );  
  fLowPtReader->AddVariable( "nvtx"     , &fNVtx      ); 
  fLowPtReader->AddVariable( "jetPt"    , &fJPt1      );  
  fLowPtReader->AddVariable( "jetEta"   , &fJEta1     );
  fLowPtReader->AddVariable( "jetPhi"   , &fJPhi1     );             
  fLowPtReader->AddVariable( "dZ"       , &fJDZ1      );
  fLowPtReader->AddVariable( "d0"       , &fJD01      );
  fLowPtReader->AddVariable( "beta"     , &fBeta      );
  fLowPtReader->AddVariable( "betaStar" , &fBetaStar  );
  fLowPtReader->AddVariable( "nCharged" , &fNCharged  );
  fLowPtReader->AddVariable( "nNeutrals", &fNNeutrals );
  fLowPtReader->AddVariable( "dRMean"   , &fDRMean    );
  fLowPtReader->AddVariable( "frac01"   , &fFrac01    );
  fLowPtReader->AddVariable( "frac02"   , &fFrac02    );
  fLowPtReader->AddVariable( "frac03"   , &fFrac03    );
  fLowPtReader->AddVariable( "frac04"   , &fFrac04    );
  fLowPtReader->AddVariable( "frac05"   , &fFrac05    );
  
  fReader        = 0;
  fReader        = new TMVA::Reader( "!Color:!Silent:Error" );  
  if (fType == kBaseline) {
    fReader->AddVariable( "nvtx"     , &fNVtx      ); 
    fReader->AddVariable( "jetPt"    , &fJPt1      );  
    fReader->AddVariable( "jetEta"   , &fJEta1     );
    fReader->AddVariable( "jetPhi"   , &fJPhi1     );             
    fReader->AddVariable( "dZ"       , &fJDZ1      );
    fReader->AddVariable( "d0"       , &fJD01      );
    fReader->AddVariable( "beta"     , &fBeta      );
    fReader->AddVariable( "betaStar" , &fBetaStar  );
    fReader->AddVariable( "nCharged" , &fNCharged  );
    fReader->AddVariable( "nNeutrals", &fNNeutrals );
    fReader->AddVariable( "dRMean"   , &fDRMean    );
    fReader->AddVariable( "frac01"   , &fFrac01    );
    fReader->AddVariable( "frac02"   , &fFrac02    );
    fReader->AddVariable( "frac03"   , &fFrac03    );
    fReader->AddVariable( "frac04"   , &fFrac04    );
    fReader->AddVariable( "frac05"   , &fFrac05    );
  }
  if (fType == k42) {
    fReader->AddVariable( "frac01"   , &fFrac01    );
    fReader->AddVariable( "frac02"   , &fFrac02    );
    fReader->AddVariable( "frac03"   , &fFrac03    );
    fReader->AddVariable( "frac04"   , &fFrac04    );
    fReader->AddVariable( "frac05"   , &fFrac05    );
    fReader->AddVariable( "nvtx"     , &fNVtx      ); 
    fReader->AddVariable( "nNeutrals", &fNNeutrals );
    fReader->AddVariable( "beta"     , &fBeta      );
    fReader->AddVariable( "betaStar" , &fBetaStar  );
    fReader->AddVariable( "dZ"       , &fJDZ1      );
    fReader->AddVariable( "nCharged" , &fNCharged  );
    fReader->AddSpectator( "jetPt"    , &fJPt1      );  
    fReader->AddSpectator( "jetEta"   , &fJEta1     );
  }
  if (fType == k52) {
    fReader->AddVariable( "frac01"   , &fFrac01    );
    fReader->AddVariable( "frac02"   , &fFrac02    );
    fReader->AddVariable( "frac03"   , &fFrac03    );
    fReader->AddVariable( "frac04"   , &fFrac04    );
    fReader->AddVariable( "frac05"   , &fFrac05    );
    fReader->AddVariable( "dR2Mean"  , &fDR2Mean   );
    fReader->AddVariable( "nvtx"     , &fNVtx      ); 
    fReader->AddVariable( "nNeutrals", &fNNeutrals );
    fReader->AddVariable( "beta"     , &fBeta      );
    fReader->AddVariable( "betaStar" , &fBetaStar  );
    fReader->AddVariable( "dZ"       , &fJDZ1      );
    fReader->AddVariable( "nCharged" , &fNCharged  );
    fReader->AddSpectator( "jetPt"    , &fJPt1      );  
    fReader->AddSpectator( "jetEta"   , &fJEta1     );
  }

  fLowPtReader->BookMVA(fLowPtMethodName  , iLowPtWeights );
  fReader->BookMVA(fHighPtMethodName , iHighPtWeights );
  std::cout << "Jet ID MVA Initialization\n";
  std::cout << "MethodName : " << fLowPtMethodName << " , type == " << fType << std::endl;

  std::string lCutId = "JetIdParams";
  if(fType == k42) lCutId = "PuJetIdOptMVA_wp";
  if(fType == k52) lCutId = "full_5x_wp";
  //Load Cut Matrix
  edm::ParameterSet lDefConfig = edm::readPSetsFrom(iCutFileName.Data())->getParameter<edm::ParameterSet>("JetIdParams");
  edm::ParameterSet lConfig    = edm::readPSetsFrom(iCutFileName.Data())->getParameter<edm::ParameterSet>(lCutId);
  std::string lCutType = "Tight";
  if(fCutType == kMedium) lCutType = "Medium";
  if(fCutType == kLoose ) lCutType = "Loose";
  if(fCutType == kMET   ) lCutType = "MET";
  std::string lLowPtCut = "MET";
  std::vector<double> lPt010  = lDefConfig.getParameter<std::vector<double> >(("Pt010_" +lLowPtCut).c_str());
  std::vector<double> lPt1020 = lConfig.getParameter<std::vector<double> >(("Pt1020_"+lCutType).c_str());
  std::vector<double> lPt2030 = lConfig.getParameter<std::vector<double> >(("Pt2030_"+lCutType).c_str());
  std::vector<double> lPt3050 = lConfig.getParameter<std::vector<double> >(("Pt3050_"+lCutType).c_str());
  for(int i0 = 0; i0 < 4; i0++) fMVACut[0][i0] = lPt010 [i0];
  for(int i0 = 0; i0 < 4; i0++) fMVACut[1][i0] = lPt1020[i0];
  for(int i0 = 0; i0 < 4; i0++) fMVACut[2][i0] = lPt2030[i0];
  for(int i0 = 0; i0 < 4; i0++) fMVACut[3][i0] = lPt3050[i0];
  //std::cout << " Working Points : << " << std::endl;
  //for(int i0 = 0; i0 < 4; i0++) for(int i1 = 0; i1 < 4; i1++) 
  //  std::cout << " ==> " << i0 << " -- " << i1 << " -- " << fMVACut[i0][i1] << std::endl;
}

//--------------------------------------------------------------------------------------------------
Double_t JetIDMVA::MVAValue(    
                            Float_t iNPV    ,
                            Float_t iJPt1   ,
                            Float_t iJEta1  ,
                            Float_t iJPhi1  ,
                            Float_t iJD01   ,
                            Float_t iJDZ1   ,
                            Float_t iBeta   ,
                            Float_t iBetaStar,
                            Float_t iNCharged,
                            Float_t iNNeutrals,
                            Float_t iDRMean  ,
                            Float_t iFrac01  ,
                            Float_t iFrac02  ,
                            Float_t iFrac03  ,
                            Float_t iFrac04  ,
                            Float_t iFrac05  ,
                            Float_t iDR2Mean  
                            ){
  
  if(!fIsInitialized) { 
    std::cout << "Error: JetIDMVA not properly initialized.\n"; 
    return -9999;
  }
  
  fNVtx      = iNPV;
  fJPt1      = iJPt1;
  fJEta1     = iJEta1;
  fJPhi1     = fJPhi1;
  fJD01      = iJD01;
  fJDZ1      = iJDZ1;
  fBeta      = iBeta;
  fBetaStar  = iBetaStar;
  fNCharged  = iNCharged;
  fNNeutrals = iNNeutrals;
  fDRMean    = iDRMean;
  fFrac01    = iFrac01;
  fFrac02    = iFrac02;
  fFrac03    = iFrac03;
  fFrac04    = iFrac04;
  fFrac05    = iFrac05;
  fDR2Mean   = iDR2Mean;

  Double_t lMVA = -9999;  
  if(iJPt1 < 10) lMVA = fLowPtReader->EvaluateMVA( fLowPtMethodName  );
  if(iJPt1 > 10) lMVA = fReader->EvaluateMVA( fHighPtMethodName );

  return lMVA;
}
//--------------------------------------------------------------------------------------------------
Bool_t JetIDMVA::pass(const PFJet *iJet,const Vertex *iVertex,const VertexCol *iVertices,
                      FactorizedJetCorrector *iJetCorrector,
                      const PileupEnergyDensityCol *iPileupEnergyDensity,
                      RhoUtilities::RhoType type) { 
  
  if(!JetTools::passPFLooseId(iJet))                 return false;
  if(fabs(iJet->Eta()) > 4.99)                       return false;
  
  double lMVA = MVAValue   (iJet,iVertex,iVertices,iJetCorrector,iPileupEnergyDensity);
  double lPt  = correctedPt(iJet,                  iJetCorrector,iPileupEnergyDensity,type);
  if(lPt < fJetPtMin)                         return false; 
  
  int lPtId = 0; 
  if(lPt > 10 && lPt < 20) lPtId = 1;
  if(lPt > 20 && lPt < 30) lPtId = 2;
  if(lPt > 30                   ) lPtId = 3;
  
  int lEtaId = 0;
  if(fabs(iJet->Eta()) > 2.5  && fabs(iJet->Eta()) < 2.75) lEtaId = 1; 
  if(fabs(iJet->Eta()) > 2.75 && fabs(iJet->Eta()) < 3.0 ) lEtaId = 2; 
  if(fabs(iJet->Eta()) > 3.0  && fabs(iJet->Eta()) < 5.0 ) lEtaId = 3; 
  
  double lMVACut = fMVACut[lPtId][lEtaId];
  if(lMVA < lMVACut) return false;
  return true;
}
//--------------------------------------------------------------------------------------------------
Bool_t JetIDMVA::pass(const PFJet *iJet,const Vertex *iVertex,const VertexCol *iVertices) { 
  if(!JetTools::passPFLooseId(iJet))                 return false;
  if(iJet->Pt()        < fJetPtMin) return false; 
  if(fabs(iJet->Eta()) > 4.99)      return false;
  double lMVA = MVAValue(iJet,iVertex,iVertices);
  
  int lPtId = 0; 
  if(iJet->Pt() > 10 && iJet->Pt() < 20) lPtId = 1;
  if(iJet->Pt() > 20 && iJet->Pt() < 30) lPtId = 2;
  if(iJet->Pt() > 30                   ) lPtId = 3;
  
  int lEtaId = 0;
  if(fabs(iJet->Eta()) > 2.5  && fabs(iJet->Eta()) < 2.75) lEtaId = 1; 
  if(fabs(iJet->Eta()) > 2.75 && fabs(iJet->Eta()) < 3.0 ) lEtaId = 2; 
  if(fabs(iJet->Eta()) > 3.0  && fabs(iJet->Eta()) < 5.0 ) lEtaId = 3; 
  
  double lMVACut = fMVACut[lPtId][lEtaId];
  if(lMVA < lMVACut) return false;
  return true;
  //if(lMVA < -0.8)                            return false;
  //if(lMVA < -0.5 && fabs(iJet->Eta()) > 3.0) return false;
  //if(iJet->Pt() < fJetPtMin && iJet->TrackCountingHighEffBJetTagsDisc() == -100) return false; //This line is a bug in the Met training
  //if( fabs(JetTools::impactParameter(iJet,iVertex,true)) < 0.2) return true;
}
//--------------------------------------------------------------------------------------------------
Double_t JetIDMVA::MVAValue(const PFJet *iJet,const Vertex *iVertex,const VertexCol *iVertices, //Vertex here is the PV
                            FactorizedJetCorrector *iJetCorrector,
                            const PileupEnergyDensityCol *iPileupEnergyDensity,
                            Bool_t printDebug) {
  
  if (!fIsInitialized) { 
    std::cout << "Error: JetIDMVA not properly initialized.\n"; 
    return -9999;
  }
  if(!JetTools::passPFLooseId(iJet)) return -2.;

  //set all input variables
  fNVtx      = iVertices->GetEntries();
  fJPt1      = correctedPt(iJet,iJetCorrector,iPileupEnergyDensity);
  fJEta1     = iJet->RawMom().Eta();
  fJPhi1     = iJet->RawMom().Phi();
  fJD01      = JetTools::impactParameter(iJet,iVertex);  
  fJDZ1      = JetTools::impactParameter(iJet,iVertex,true);
  fBeta      = JetTools::Beta(iJet,iVertex,fDZCut);
  fBetaStar  = JetTools::betaStar(iJet,iVertex,iVertices,fDZCut);
  fNCharged  = iJet->ChargedMultiplicity();
  fNNeutrals = iJet->NeutralMultiplicity();

  fDRMean    = JetTools::dRMean(iJet,-1);
  fDR2Mean   = JetTools::dR2Mean(iJet,-1);
  fFrac01    = JetTools::frac  (iJet,0.1,0. ,-1);
  fFrac02    = JetTools::frac  (iJet,0.2,0.1,-1);
  fFrac03    = JetTools::frac  (iJet,0.3,0.2,-1);
  fFrac04    = JetTools::frac  (iJet,0.4,0.3,-1);
  fFrac05    = JetTools::frac  (iJet,0.5,0.4,-1);

  double lMVA = 0;
  if(fJPt1 < 10) lMVA = fLowPtReader->EvaluateMVA( fLowPtMethodName  );
  if(fJPt1 > 10) lMVA = fReader->EvaluateMVA( fHighPtMethodName );
  if (printDebug == kTRUE) {
    std::cout << "Debug Jet MVA: "
              << fNVtx      << " "
              << fJPt1      << " "
              << fJEta1     << " "
              << fJPhi1     << " "
              << fJD01      << " "
              << fJDZ1      << " "
              << fBeta      << " "
              << fBetaStar  << " "
              << fNCharged  << " "
              << fNNeutrals << " "
              << fDRMean    << " "
              << fFrac01    << " "
              << fFrac02    << " "
              << fFrac03    << " "
              << fFrac04    << " "
              << fFrac05    << " "
              << fDRMean    
              << " === : === "
              << lMVA << " "    
              << std::endl;
  }

  return lMVA;
}
Double_t JetIDMVA::MVAValue(const PFJet *iJet,const Vertex *iVertex, const VertexCol *iVertices,//Vertex here is the PV
                            Bool_t printDebug) {
  
  if (!fIsInitialized) { 
    std::cout << "Error: JetIDMVA not properly initialized.\n"; 
    return -9999;
  }
  if(!JetTools::passPFLooseId(iJet)) return -2.;

  //set all input variables
  fNVtx      = iVertices->GetEntries();
  fJPt1      = iJet->Pt();
  fJEta1     = iJet->RawMom().Eta();
  fJPhi1     = iJet->RawMom().Phi();
  fJD01      = JetTools::impactParameter(iJet,iVertex);  
  fJDZ1      = JetTools::impactParameter(iJet,iVertex,true);
  fBeta      = JetTools::Beta(iJet,iVertex,fDZCut);
  fBetaStar  = JetTools::betaStar(iJet,iVertex,iVertices,fDZCut);
  fNCharged  = iJet->ChargedMultiplicity();
  fNNeutrals = iJet->NeutralMultiplicity();

  fDRMean    = JetTools::dRMean (iJet,-1);
  fDR2Mean   = JetTools::dR2Mean(iJet,-1);
  fFrac01    = JetTools::frac   (iJet,0.1,0. ,-1);
  fFrac02    = JetTools::frac   (iJet,0.2,0.1,-1);
  fFrac03    = JetTools::frac   (iJet,0.3,0.2,-1);
  fFrac04    = JetTools::frac   (iJet,0.4,0.3,-1);
  fFrac05    = JetTools::frac   (iJet,0.5,0.4,-1);

  double lMVA = 0;
  if(fJPt1 < 10) lMVA = fLowPtReader->EvaluateMVA( fLowPtMethodName  );
  if(fJPt1 > 10) lMVA = fReader     ->EvaluateMVA( fHighPtMethodName );
  
  if (printDebug == kTRUE) {
    std::cout << "Debug Jet MVA: "
              << fNVtx      << " "
              << fJPt1      << " "
              << fJEta1     << " "
              << fJPhi1     << " "
              << fJD01      << " "
              << fJDZ1      << " "
              << fBeta      << " "
              << fBetaStar  << " "
              << fNCharged  << " "
              << fNNeutrals << " "
              << fDRMean    << " "
              << fFrac01    << " "
              << fFrac02    << " "
              << fFrac03    << " "
              << fFrac04    << " "
              << fFrac05    << " "
              << fDRMean    
              << " === : === "
              << lMVA << " "    
              << std::endl;
  }

  return lMVA;
}
Double_t JetIDMVA::correctedPt(const PFJet *iJet, FactorizedJetCorrector *iJetCorrector,
                               const PileupEnergyDensityCol *iPUEnergyDensity,
                               RhoUtilities::RhoType type) { 
  Double_t Rho = 0.0;
  switch(type) {
  case RhoUtilities::MIT_RHO_VORONOI_HIGH_ETA:
    Rho = iPUEnergyDensity->At(0)->Rho();
    break;
  case RhoUtilities::MIT_RHO_VORONOI_LOW_ETA:
    Rho = iPUEnergyDensity->At(0)->RhoLowEta();
    break;
  case RhoUtilities::MIT_RHO_RANDOM_HIGH_ETA:
    Rho = iPUEnergyDensity->At(0)->RhoRandom();
    break;
  case RhoUtilities::MIT_RHO_RANDOM_LOW_ETA:
    Rho = iPUEnergyDensity->At(0)->RhoRandomLowEta();
    break;
  case RhoUtilities::CMS_RHO_RHOKT6PFJETS:
    Rho = iPUEnergyDensity->At(0)->RhoKt6PFJets();
    break;
  default:
    // use the old default
    Rho = iPUEnergyDensity->At(0)->Rho();
    break;
  }
    
  const FourVectorM rawMom = iJet->RawMom();
  iJetCorrector->setJetEta(rawMom.Eta());
  iJetCorrector->setJetPt (rawMom.Pt());
  iJetCorrector->setJetPhi(rawMom.Phi());
  iJetCorrector->setJetE  (rawMom.E());
  iJetCorrector->setRho   (Rho);
  iJetCorrector->setJetA  (iJet->JetArea());
  iJetCorrector->setJetEMF(-99.0);     
  Double_t correction = iJetCorrector->getCorrection();
  return rawMom.Pt()*correction;
}
Revision:	1.11
Committed:	Mon May 14 09:42:27 2012 UTC (12 years, 11 months ago) by pharris
Content type:	text/plain
Branch:	MAIN
Changes since 1.10:	+82 -18 lines
Log Message:	Added final MVA
#	User	Rev	Content
1	pharris	1.4	#include "FWCore/PythonParameterSet/interface/MakeParameterSets.h"
2			#include "FWCore/ParameterSet/interface/ParameterSet.h"
3	pharris	1.1	#include "MitPhysics/Utils/interface/JetIDMVA.h"
4			#include "MitPhysics/Utils/interface/JetTools.h"
5			#include "MitAna/DataTree/interface/StableData.h"
6			#include <TFile.h>
7			#include <TRandom3.h>
8			#include "TMVA/Tools.h"
9			#include "TMVA/Reader.h"
10
11
12			ClassImp(mithep::JetIDMVA)
13
14			using namespace mithep;
15
16			//--------------------------------------------------------------------------------------------------
17			JetIDMVA::JetIDMVA() :
18	pharris	1.8	fJetPtMin(0.) , //We need to lower this
19	pharris	1.4	fDZCut (0.2),
20			fLowPtMethodName ("JetIDMVALowPt" ),
21			fHighPtMethodName("JetIDMVAHighPt"),
22	pharris	1.1	fIsInitialized(kFALSE),
23	pharris	1.4	fNVtx (0),
24			fJPt1 (0),
25			fJEta1 (0),
26			fJPhi1 (0),
27			fJD01 (0),
28			fJDZ1 (0),
29			fBeta (0),
30			fBetaStar (0),
31			fNCharged (0),
32			fNNeutrals(0),
33			fDRMean (0),
34			fFrac01 (0),
35			fFrac02 (0),
36			fFrac03 (0),
37			fFrac04 (0),
38	pharris	1.11	fFrac05 (0),
39			fDR2Mean (0)
40	pharris	1.1	{
41			fReader = 0;
42			}
43			//--------------------------------------------------------------------------------------------------
44	pharris	1.4	JetIDMVA::~JetIDMVA() {
45
46	pharris	1.11	delete fReader;
47			delete fLowPtReader;
48	pharris	1.1	}
49
50			//--------------------------------------------------------------------------------------------------
51	pharris	1.4	void JetIDMVA::Initialize( JetIDMVA::CutType iCutType,
52			TString iLowPtWeights,
53			TString iHighPtWeights,
54			JetIDMVA::MVAType iType,
55			TString iCutFileName)
56	pharris	1.1	{
57
58			fIsInitialized = kTRUE;
59			fType = iType;
60	pharris	1.4	fCutType = iCutType;
61	pharris	1.11	fLowPtReader = 0;
62			fLowPtReader = new TMVA::Reader( "!Color:!Silent:Error" );
63			fLowPtReader->AddVariable( "nvtx" , &fNVtx );
64			fLowPtReader->AddVariable( "jetPt" , &fJPt1 );
65			fLowPtReader->AddVariable( "jetEta" , &fJEta1 );
66			fLowPtReader->AddVariable( "jetPhi" , &fJPhi1 );
67			fLowPtReader->AddVariable( "dZ" , &fJDZ1 );
68			fLowPtReader->AddVariable( "d0" , &fJD01 );
69			fLowPtReader->AddVariable( "beta" , &fBeta );
70			fLowPtReader->AddVariable( "betaStar" , &fBetaStar );
71			fLowPtReader->AddVariable( "nCharged" , &fNCharged );
72			fLowPtReader->AddVariable( "nNeutrals", &fNNeutrals );
73			fLowPtReader->AddVariable( "dRMean" , &fDRMean );
74			fLowPtReader->AddVariable( "frac01" , &fFrac01 );
75			fLowPtReader->AddVariable( "frac02" , &fFrac02 );
76			fLowPtReader->AddVariable( "frac03" , &fFrac03 );
77			fLowPtReader->AddVariable( "frac04" , &fFrac04 );
78			fLowPtReader->AddVariable( "frac05" , &fFrac05 );
79
80	pharris	1.1	fReader = 0;
81			fReader = new TMVA::Reader( "!Color:!Silent:Error" );
82			if (fType == kBaseline) {
83	pharris	1.4	fReader->AddVariable( "nvtx" , &fNVtx );
84			fReader->AddVariable( "jetPt" , &fJPt1 );
85			fReader->AddVariable( "jetEta" , &fJEta1 );
86			fReader->AddVariable( "jetPhi" , &fJPhi1 );
87			fReader->AddVariable( "dZ" , &fJDZ1 );
88			fReader->AddVariable( "d0" , &fJD01 );
89			fReader->AddVariable( "beta" , &fBeta );
90			fReader->AddVariable( "betaStar" , &fBetaStar );
91			fReader->AddVariable( "nCharged" , &fNCharged );
92			fReader->AddVariable( "nNeutrals", &fNNeutrals );
93			fReader->AddVariable( "dRMean" , &fDRMean );
94			fReader->AddVariable( "frac01" , &fFrac01 );
95			fReader->AddVariable( "frac02" , &fFrac02 );
96			fReader->AddVariable( "frac03" , &fFrac03 );
97			fReader->AddVariable( "frac04" , &fFrac04 );
98			fReader->AddVariable( "frac05" , &fFrac05 );
99	pharris	1.1	}
100	pharris	1.11	if (fType == k42) {
101			fReader->AddVariable( "frac01" , &fFrac01 );
102			fReader->AddVariable( "frac02" , &fFrac02 );
103			fReader->AddVariable( "frac03" , &fFrac03 );
104			fReader->AddVariable( "frac04" , &fFrac04 );
105			fReader->AddVariable( "frac05" , &fFrac05 );
106			fReader->AddVariable( "nvtx" , &fNVtx );
107			fReader->AddVariable( "nNeutrals", &fNNeutrals );
108			fReader->AddVariable( "beta" , &fBeta );
109			fReader->AddVariable( "betaStar" , &fBetaStar );
110			fReader->AddVariable( "dZ" , &fJDZ1 );
111			fReader->AddVariable( "nCharged" , &fNCharged );
112			fReader->AddSpectator( "jetPt" , &fJPt1 );
113			fReader->AddSpectator( "jetEta" , &fJEta1 );
114			}
115			if (fType == k52) {
116			fReader->AddVariable( "frac01" , &fFrac01 );
117			fReader->AddVariable( "frac02" , &fFrac02 );
118			fReader->AddVariable( "frac03" , &fFrac03 );
119			fReader->AddVariable( "frac04" , &fFrac04 );
120			fReader->AddVariable( "frac05" , &fFrac05 );
121			fReader->AddVariable( "dR2Mean" , &fDR2Mean );
122			fReader->AddVariable( "nvtx" , &fNVtx );
123			fReader->AddVariable( "nNeutrals", &fNNeutrals );
124			fReader->AddVariable( "beta" , &fBeta );
125			fReader->AddVariable( "betaStar" , &fBetaStar );
126			fReader->AddVariable( "dZ" , &fJDZ1 );
127			fReader->AddVariable( "nCharged" , &fNCharged );
128			fReader->AddSpectator( "jetPt" , &fJPt1 );
129			fReader->AddSpectator( "jetEta" , &fJEta1 );
130			}
131
132			fLowPtReader->BookMVA(fLowPtMethodName , iLowPtWeights );
133	pharris	1.4	fReader->BookMVA(fHighPtMethodName , iHighPtWeights );
134	pharris	1.1	std::cout << "Jet ID MVA Initialization\n";
135	pharris	1.4	std::cout << "MethodName : " << fLowPtMethodName << " , type == " << fType << std::endl;
136
137	pharris	1.11	std::string lCutId = "JetIdParams";
138			if(fType == k42) lCutId = "PuJetIdOptMVA_wp";
139			if(fType == k52) lCutId = "full_5x_wp";
140	pharris	1.4	//Load Cut Matrix
141	pharris	1.11	edm::ParameterSet lDefConfig = edm::readPSetsFrom(iCutFileName.Data())->getParameter<edm::ParameterSet>("JetIdParams");
142			edm::ParameterSet lConfig = edm::readPSetsFrom(iCutFileName.Data())->getParameter<edm::ParameterSet>(lCutId);
143	pharris	1.4	std::string lCutType = "Tight";
144			if(fCutType == kMedium) lCutType = "Medium";
145			if(fCutType == kLoose ) lCutType = "Loose";
146	pharris	1.6	if(fCutType == kMET ) lCutType = "MET";
147	pharris	1.11	std::string lLowPtCut = "MET";
148			std::vector<double> lPt010 = lDefConfig.getParameter<std::vector<double> >(("Pt010_" +lLowPtCut).c_str());
149	pharris	1.4	std::vector<double> lPt1020 = lConfig.getParameter<std::vector<double> >(("Pt1020_"+lCutType).c_str());
150			std::vector<double> lPt2030 = lConfig.getParameter<std::vector<double> >(("Pt2030_"+lCutType).c_str());
151			std::vector<double> lPt3050 = lConfig.getParameter<std::vector<double> >(("Pt3050_"+lCutType).c_str());
152			for(int i0 = 0; i0 < 4; i0++) fMVACut[0][i0] = lPt010 [i0];
153			for(int i0 = 0; i0 < 4; i0++) fMVACut[1][i0] = lPt1020[i0];
154			for(int i0 = 0; i0 < 4; i0++) fMVACut[2][i0] = lPt2030[i0];
155			for(int i0 = 0; i0 < 4; i0++) fMVACut[3][i0] = lPt3050[i0];
156	pharris	1.5	//std::cout << " Working Points : << " << std::endl;
157			//for(int i0 = 0; i0 < 4; i0++) for(int i1 = 0; i1 < 4; i1++)
158			// std::cout << " ==> " << i0 << " -- " << i1 << " -- " << fMVACut[i0][i1] << std::endl;
159	pharris	1.1	}
160
161			//--------------------------------------------------------------------------------------------------
162			Double_t JetIDMVA::MVAValue(
163			Float_t iNPV ,
164			Float_t iJPt1 ,
165			Float_t iJEta1 ,
166			Float_t iJPhi1 ,
167			Float_t iJD01 ,
168			Float_t iJDZ1 ,
169	pharris	1.4	Float_t iBeta ,
170			Float_t iBetaStar,
171			Float_t iNCharged,
172			Float_t iNNeutrals,
173			Float_t iDRMean ,
174			Float_t iFrac01 ,
175			Float_t iFrac02 ,
176			Float_t iFrac03 ,
177			Float_t iFrac04 ,
178	pharris	1.11	Float_t iFrac05 ,
179			Float_t iDR2Mean
180	pharris	1.1	){
181
182	pharris	1.4	if(!fIsInitialized) {
183	pharris	1.1	std::cout << "Error: JetIDMVA not properly initialized.\n";
184			return -9999;
185			}
186
187	pharris	1.4	fNVtx = iNPV;
188			fJPt1 = iJPt1;
189			fJEta1 = iJEta1;
190			fJPhi1 = fJPhi1;
191			fJD01 = iJD01;
192			fJDZ1 = iJDZ1;
193			fBeta = iBeta;
194			fBetaStar = iBetaStar;
195			fNCharged = iNCharged;
196			fNNeutrals = iNNeutrals;
197			fDRMean = iDRMean;
198			fFrac01 = iFrac01;
199			fFrac02 = iFrac02;
200			fFrac03 = iFrac03;
201			fFrac04 = iFrac04;
202			fFrac05 = iFrac05;
203	pharris	1.11	fDR2Mean = iDR2Mean;
204	pharris	1.1
205			Double_t lMVA = -9999;
206	pharris	1.11	if(iJPt1 < 10) lMVA = fLowPtReader->EvaluateMVA( fLowPtMethodName );
207	pharris	1.4	if(iJPt1 > 10) lMVA = fReader->EvaluateMVA( fHighPtMethodName );
208	pharris	1.1
209			return lMVA;
210			}
211			//--------------------------------------------------------------------------------------------------
212	pharris	1.4	Bool_t JetIDMVA::pass(const PFJet iJet,const Vertex iVertex,const VertexCol *iVertices,
213	pharris	1.1	FactorizedJetCorrector *iJetCorrector,
214	ceballos	1.10	const PileupEnergyDensityCol *iPileupEnergyDensity,
215			RhoUtilities::RhoType type) {
216	pharris	1.4
217	pharris	1.1	if(!JetTools::passPFLooseId(iJet)) return false;
218	pharris	1.7	if(fabs(iJet->Eta()) > 4.99) return false;
219	pharris	1.5
220	pharris	1.4	double lMVA = MVAValue (iJet,iVertex,iVertices,iJetCorrector,iPileupEnergyDensity);
221	ceballos	1.10	double lPt = correctedPt(iJet, iJetCorrector,iPileupEnergyDensity,type);
222	pharris	1.8	if(lPt < fJetPtMin) return false;
223	pharris	1.7
224	pharris	1.4	int lPtId = 0;
225	pharris	1.7	if(lPt > 10 && lPt < 20) lPtId = 1;
226			if(lPt > 20 && lPt < 30) lPtId = 2;
227	pharris	1.4	if(lPt > 30 ) lPtId = 3;
228
229			int lEtaId = 0;
230			if(fabs(iJet->Eta()) > 2.5 && fabs(iJet->Eta()) < 2.75) lEtaId = 1;
231			if(fabs(iJet->Eta()) > 2.75 && fabs(iJet->Eta()) < 3.0 ) lEtaId = 2;
232			if(fabs(iJet->Eta()) > 3.0 && fabs(iJet->Eta()) < 5.0 ) lEtaId = 3;
233
234			double lMVACut = fMVACut[lPtId][lEtaId];
235			if(lMVA < lMVACut) return false;
236	pharris	1.1	return true;
237			}
238			//--------------------------------------------------------------------------------------------------
239	pharris	1.4	Bool_t JetIDMVA::pass(const PFJet iJet,const Vertex iVertex,const VertexCol *iVertices) {
240	pharris	1.2	if(!JetTools::passPFLooseId(iJet)) return false;
241	pharris	1.5	if(iJet->Pt() < fJetPtMin) return false;
242	pharris	1.7	if(fabs(iJet->Eta()) > 4.99) return false;
243	pharris	1.4	double lMVA = MVAValue(iJet,iVertex,iVertices);
244
245			int lPtId = 0;
246			if(iJet->Pt() > 10 && iJet->Pt() < 20) lPtId = 1;
247			if(iJet->Pt() > 20 && iJet->Pt() < 30) lPtId = 2;
248			if(iJet->Pt() > 30 ) lPtId = 3;
249
250			int lEtaId = 0;
251			if(fabs(iJet->Eta()) > 2.5 && fabs(iJet->Eta()) < 2.75) lEtaId = 1;
252			if(fabs(iJet->Eta()) > 2.75 && fabs(iJet->Eta()) < 3.0 ) lEtaId = 2;
253			if(fabs(iJet->Eta()) > 3.0 && fabs(iJet->Eta()) < 5.0 ) lEtaId = 3;
254
255			double lMVACut = fMVACut[lPtId][lEtaId];
256			if(lMVA < lMVACut) return false;
257	pharris	1.2	return true;
258	pharris	1.4	//if(lMVA < -0.8) return false;
259			//if(lMVA < -0.5 && fabs(iJet->Eta()) > 3.0) return false;
260			//if(iJet->Pt() < fJetPtMin && iJet->TrackCountingHighEffBJetTagsDisc() == -100) return false; //This line is a bug in the Met training
261			//if( fabs(JetTools::impactParameter(iJet,iVertex,true)) < 0.2) return true;
262	pharris	1.2	}
263			//--------------------------------------------------------------------------------------------------
264	pharris	1.4	Double_t JetIDMVA::MVAValue(const PFJet iJet,const Vertex iVertex,const VertexCol *iVertices, //Vertex here is the PV
265	pharris	1.1	FactorizedJetCorrector *iJetCorrector,
266			const PileupEnergyDensityCol *iPileupEnergyDensity,
267			Bool_t printDebug) {
268
269			if (!fIsInitialized) {
270			std::cout << "Error: JetIDMVA not properly initialized.\n";
271			return -9999;
272			}
273			if(!JetTools::passPFLooseId(iJet)) return -2.;
274
275			//set all input variables
276	pharris	1.4	fNVtx = iVertices->GetEntries();
277	pharris	1.1	fJPt1 = correctedPt(iJet,iJetCorrector,iPileupEnergyDensity);
278			fJEta1 = iJet->RawMom().Eta();
279			fJPhi1 = iJet->RawMom().Phi();
280	pharris	1.4	fJD01 = JetTools::impactParameter(iJet,iVertex);
281			fJDZ1 = JetTools::impactParameter(iJet,iVertex,true);
282			fBeta = JetTools::Beta(iJet,iVertex,fDZCut);
283			fBetaStar = JetTools::betaStar(iJet,iVertex,iVertices,fDZCut);
284			fNCharged = iJet->ChargedMultiplicity();
285			fNNeutrals = iJet->NeutralMultiplicity();
286
287			fDRMean = JetTools::dRMean(iJet,-1);
288	pharris	1.11	fDR2Mean = JetTools::dR2Mean(iJet,-1);
289	pharris	1.4	fFrac01 = JetTools::frac (iJet,0.1,0. ,-1);
290			fFrac02 = JetTools::frac (iJet,0.2,0.1,-1);
291			fFrac03 = JetTools::frac (iJet,0.3,0.2,-1);
292			fFrac04 = JetTools::frac (iJet,0.4,0.3,-1);
293			fFrac05 = JetTools::frac (iJet,0.5,0.4,-1);
294
295			double lMVA = 0;
296	pharris	1.11	if(fJPt1 < 10) lMVA = fLowPtReader->EvaluateMVA( fLowPtMethodName );
297	pharris	1.4	if(fJPt1 > 10) lMVA = fReader->EvaluateMVA( fHighPtMethodName );
298	pharris	1.1	if (printDebug == kTRUE) {
299			std::cout << "Debug Jet MVA: "
300	pharris	1.4	<< fNVtx << " "
301			<< fJPt1 << " "
302			<< fJEta1 << " "
303			<< fJPhi1 << " "
304			<< fJD01 << " "
305			<< fJDZ1 << " "
306			<< fBeta << " "
307			<< fBetaStar << " "
308			<< fNCharged << " "
309			<< fNNeutrals << " "
310			<< fDRMean << " "
311			<< fFrac01 << " "
312			<< fFrac02 << " "
313			<< fFrac03 << " "
314			<< fFrac04 << " "
315	pharris	1.11	<< fFrac05 << " "
316			<< fDRMean
317	pharris	1.1	<< " === : === "
318			<< lMVA << " "
319			<< std::endl;
320			}
321
322			return lMVA;
323			}
324	pharris	1.4	Double_t JetIDMVA::MVAValue(const PFJet iJet,const Vertex iVertex, const VertexCol *iVertices,//Vertex here is the PV
325	pharris	1.2	Bool_t printDebug) {
326
327			if (!fIsInitialized) {
328			std::cout << "Error: JetIDMVA not properly initialized.\n";
329			return -9999;
330			}
331			if(!JetTools::passPFLooseId(iJet)) return -2.;
332
333			//set all input variables
334	pharris	1.4	fNVtx = iVertices->GetEntries();
335	pharris	1.2	fJPt1 = iJet->Pt();
336			fJEta1 = iJet->RawMom().Eta();
337			fJPhi1 = iJet->RawMom().Phi();
338	pharris	1.4	fJD01 = JetTools::impactParameter(iJet,iVertex);
339			fJDZ1 = JetTools::impactParameter(iJet,iVertex,true);
340			fBeta = JetTools::Beta(iJet,iVertex,fDZCut);
341			fBetaStar = JetTools::betaStar(iJet,iVertex,iVertices,fDZCut);
342			fNCharged = iJet->ChargedMultiplicity();
343			fNNeutrals = iJet->NeutralMultiplicity();
344
345	pharris	1.11	fDRMean = JetTools::dRMean (iJet,-1);
346			fDR2Mean = JetTools::dR2Mean(iJet,-1);
347			fFrac01 = JetTools::frac (iJet,0.1,0. ,-1);
348			fFrac02 = JetTools::frac (iJet,0.2,0.1,-1);
349			fFrac03 = JetTools::frac (iJet,0.3,0.2,-1);
350			fFrac04 = JetTools::frac (iJet,0.4,0.3,-1);
351			fFrac05 = JetTools::frac (iJet,0.5,0.4,-1);
352	pharris	1.4
353			double lMVA = 0;
354	pharris	1.11	if(fJPt1 < 10) lMVA = fLowPtReader->EvaluateMVA( fLowPtMethodName );
355			if(fJPt1 > 10) lMVA = fReader ->EvaluateMVA( fHighPtMethodName );
356	pharris	1.5
357	pharris	1.2	if (printDebug == kTRUE) {
358			std::cout << "Debug Jet MVA: "
359	pharris	1.4	<< fNVtx << " "
360			<< fJPt1 << " "
361			<< fJEta1 << " "
362			<< fJPhi1 << " "
363			<< fJD01 << " "
364			<< fJDZ1 << " "
365			<< fBeta << " "
366			<< fBetaStar << " "
367			<< fNCharged << " "
368			<< fNNeutrals << " "
369			<< fDRMean << " "
370			<< fFrac01 << " "
371			<< fFrac02 << " "
372			<< fFrac03 << " "
373			<< fFrac04 << " "
374	pharris	1.11	<< fFrac05 << " "
375			<< fDRMean
376	pharris	1.2	<< " === : === "
377			<< lMVA << " "
378			<< std::endl;
379			}
380
381			return lMVA;
382			}
383	ceballos	1.10	Double_t JetIDMVA::correctedPt(const PFJet iJet, FactorizedJetCorrector iJetCorrector,
384			const PileupEnergyDensityCol *iPUEnergyDensity,
385			RhoUtilities::RhoType type) {
386			Double_t Rho = 0.0;
387			switch(type) {
388			case RhoUtilities::MIT_RHO_VORONOI_HIGH_ETA:
389			Rho = iPUEnergyDensity->At(0)->Rho();
390			break;
391			case RhoUtilities::MIT_RHO_VORONOI_LOW_ETA:
392			Rho = iPUEnergyDensity->At(0)->RhoLowEta();
393			break;
394			case RhoUtilities::MIT_RHO_RANDOM_HIGH_ETA:
395			Rho = iPUEnergyDensity->At(0)->RhoRandom();
396			break;
397			case RhoUtilities::MIT_RHO_RANDOM_LOW_ETA:
398			Rho = iPUEnergyDensity->At(0)->RhoRandomLowEta();
399			break;
400			case RhoUtilities::CMS_RHO_RHOKT6PFJETS:
401			Rho = iPUEnergyDensity->At(0)->RhoKt6PFJets();
402			break;
403			default:
404			// use the old default
405			Rho = iPUEnergyDensity->At(0)->Rho();
406			break;
407			}
408
409	pharris	1.1	const FourVectorM rawMom = iJet->RawMom();
410			iJetCorrector->setJetEta(rawMom.Eta());
411			iJetCorrector->setJetPt (rawMom.Pt());
412			iJetCorrector->setJetPhi(rawMom.Phi());
413			iJetCorrector->setJetE (rawMom.E());
414	ceballos	1.10	iJetCorrector->setRho (Rho);
415	pharris	1.1	iJetCorrector->setJetA (iJet->JetArea());
416			iJetCorrector->setJetEMF(-99.0);
417			Double_t correction = iJetCorrector->getCorrection();
418			return rawMom.Pt()*correction;
419			}