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)
{    
  fReader = 0;
}
//--------------------------------------------------------------------------------------------------
JetIDMVA::~JetIDMVA() {

  fReader = 0;
}

//--------------------------------------------------------------------------------------------------
void JetIDMVA::Initialize( JetIDMVA::CutType iCutType,
                           TString iLowPtWeights,
                           TString iHighPtWeights, 
                           JetIDMVA::MVAType iType,
                           TString iCutFileName) 
{ 
  
  fIsInitialized = kTRUE;
  fType          = iType;
  fCutType       = iCutType;
  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    );
  }
  fReader->BookMVA(fLowPtMethodName  , iLowPtWeights );
  fReader->BookMVA(fHighPtMethodName , iHighPtWeights );
  std::cout << "Jet ID MVA Initialization\n";
  std::cout << "MethodName : " << fLowPtMethodName << " , type == " << fType << std::endl;

  //Load Cut Matrix
  edm::ParameterSet lConfig = edm::readPSetsFrom(iCutFileName.Data())->getParameter<edm::ParameterSet>("JetIdParams");
  std::string lCutType = "Tight";
  if(fCutType == kMedium) lCutType = "Medium";
  if(fCutType == kLoose ) lCutType = "Loose";
  if(fCutType == kMET   ) lCutType = "MET";
  std::vector<double> lPt010  = lConfig.getParameter<std::vector<double> >(("Pt010_" +lCutType).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  
                            ){
  
  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;

  Double_t lMVA = -9999;  
  if(iJPt1 < 10) lMVA = fReader->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) { 
  
  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);
  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;
   //if( fabs(JetTools::impactParameter(iJet,iVertex,true)) < 0.2) return true;
  //if(correctedPt(iJet,iJetCorrector,iPileupEnergyDensity) < fJetPtMin && iJet->TrackCountingHighEffBJetTagsDisc() == -100) return false; 
  //This line is a bug in the Met training
  //if(lMVA < -0.8)                            return false;
  //if(lMVA < -0.5 && fabs(iJet->Eta()) > 3.0) return false;
}
//--------------------------------------------------------------------------------------------------
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);
  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 = fReader->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    
              << " === : === "
              << 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);
  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 = fReader->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    
              << " === : === "
              << lMVA << " "    
              << std::endl;
  }

  return lMVA;
}
Double_t JetIDMVA::correctedPt(const PFJet *iJet, FactorizedJetCorrector *iJetCorrector,const PileupEnergyDensityCol *iPUEnergyDensity) { 
  const FourVectorM rawMom = iJet->RawMom();
  iJetCorrector->setJetEta(rawMom.Eta());
  iJetCorrector->setJetPt (rawMom.Pt());
  iJetCorrector->setJetPhi(rawMom.Phi());
  iJetCorrector->setJetE  (rawMom.E());
  iJetCorrector->setRho   (iPUEnergyDensity->At(0)->RhoHighEta());
  iJetCorrector->setJetA  (iJet->JetArea());
  iJetCorrector->setJetEMF(-99.0);     
  Double_t correction = iJetCorrector->getCorrection();
  return rawMom.Pt()*correction;
}
Revision:	1.9
Committed:	Sat May 12 06:17:43 2012 UTC (12 years, 11 months ago) by ceballos
Content type:	text/plain
Branch:	MAIN
Changes since 1.8:	+0 -2 lines
Log Message:	remove jpt>50 requiment
#	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			fFrac05 (0)
39	pharris	1.1	{
40			fReader = 0;
41			}
42			//--------------------------------------------------------------------------------------------------
43	pharris	1.4	JetIDMVA::~JetIDMVA() {
44
45	pharris	1.1	fReader = 0;
46			}
47
48			//--------------------------------------------------------------------------------------------------
49	pharris	1.4	void JetIDMVA::Initialize( JetIDMVA::CutType iCutType,
50			TString iLowPtWeights,
51			TString iHighPtWeights,
52			JetIDMVA::MVAType iType,
53			TString iCutFileName)
54	pharris	1.1	{
55
56			fIsInitialized = kTRUE;
57			fType = iType;
58	pharris	1.4	fCutType = iCutType;
59	pharris	1.1	fReader = 0;
60			fReader = new TMVA::Reader( "!Color:!Silent:Error" );
61			if (fType == kBaseline) {
62	pharris	1.4	fReader->AddVariable( "nvtx" , &fNVtx );
63			fReader->AddVariable( "jetPt" , &fJPt1 );
64			fReader->AddVariable( "jetEta" , &fJEta1 );
65			fReader->AddVariable( "jetPhi" , &fJPhi1 );
66			fReader->AddVariable( "dZ" , &fJDZ1 );
67			fReader->AddVariable( "d0" , &fJD01 );
68			fReader->AddVariable( "beta" , &fBeta );
69			fReader->AddVariable( "betaStar" , &fBetaStar );
70			fReader->AddVariable( "nCharged" , &fNCharged );
71			fReader->AddVariable( "nNeutrals", &fNNeutrals );
72			fReader->AddVariable( "dRMean" , &fDRMean );
73			fReader->AddVariable( "frac01" , &fFrac01 );
74			fReader->AddVariable( "frac02" , &fFrac02 );
75			fReader->AddVariable( "frac03" , &fFrac03 );
76			fReader->AddVariable( "frac04" , &fFrac04 );
77			fReader->AddVariable( "frac05" , &fFrac05 );
78	pharris	1.1	}
79	pharris	1.4	fReader->BookMVA(fLowPtMethodName , iLowPtWeights );
80			fReader->BookMVA(fHighPtMethodName , iHighPtWeights );
81	pharris	1.1	std::cout << "Jet ID MVA Initialization\n";
82	pharris	1.4	std::cout << "MethodName : " << fLowPtMethodName << " , type == " << fType << std::endl;
83
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	pharris	1.6	if(fCutType == kMET ) lCutType = "MET";
90	pharris	1.4	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	pharris	1.5	//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;
101	pharris	1.1	}
102
103			//--------------------------------------------------------------------------------------------------
104			Double_t JetIDMVA::MVAValue(
105			Float_t iNPV ,
106			Float_t iJPt1 ,
107			Float_t iJEta1 ,
108			Float_t iJPhi1 ,
109			Float_t iJD01 ,
110			Float_t iJDZ1 ,
111	pharris	1.4	Float_t iBeta ,
112			Float_t iBetaStar,
113			Float_t iNCharged,
114			Float_t iNNeutrals,
115			Float_t iDRMean ,
116			Float_t iFrac01 ,
117			Float_t iFrac02 ,
118			Float_t iFrac03 ,
119			Float_t iFrac04 ,
120			Float_t iFrac05
121	pharris	1.1	){
122
123	pharris	1.4	if(!fIsInitialized) {
124	pharris	1.1	std::cout << "Error: JetIDMVA not properly initialized.\n";
125			return -9999;
126			}
127
128	pharris	1.4	fNVtx = iNPV;
129			fJPt1 = iJPt1;
130			fJEta1 = iJEta1;
131			fJPhi1 = fJPhi1;
132			fJD01 = iJD01;
133			fJDZ1 = iJDZ1;
134			fBeta = iBeta;
135			fBetaStar = iBetaStar;
136			fNCharged = iNCharged;
137			fNNeutrals = iNNeutrals;
138			fDRMean = iDRMean;
139			fFrac01 = iFrac01;
140			fFrac02 = iFrac02;
141			fFrac03 = iFrac03;
142			fFrac04 = iFrac04;
143			fFrac05 = iFrac05;
144	pharris	1.1
145			Double_t lMVA = -9999;
146	pharris	1.4	if(iJPt1 < 10) lMVA = fReader->EvaluateMVA( fLowPtMethodName );
147			if(iJPt1 > 10) lMVA = fReader->EvaluateMVA( fHighPtMethodName );
148	pharris	1.1
149			return lMVA;
150			}
151			//--------------------------------------------------------------------------------------------------
152	pharris	1.4	Bool_t JetIDMVA::pass(const PFJet iJet,const Vertex iVertex,const VertexCol *iVertices,
153	pharris	1.1	FactorizedJetCorrector *iJetCorrector,
154			const PileupEnergyDensityCol *iPileupEnergyDensity) {
155	pharris	1.4
156	pharris	1.1	if(!JetTools::passPFLooseId(iJet)) return false;
157	pharris	1.7	if(fabs(iJet->Eta()) > 4.99) return false;
158	pharris	1.5
159	pharris	1.4	double lMVA = MVAValue (iJet,iVertex,iVertices,iJetCorrector,iPileupEnergyDensity);
160			double lPt = correctedPt(iJet, iJetCorrector,iPileupEnergyDensity);
161	pharris	1.8	if(lPt < fJetPtMin) return false;
162	pharris	1.7
163	pharris	1.4	int lPtId = 0;
164	pharris	1.7	if(lPt > 10 && lPt < 20) lPtId = 1;
165			if(lPt > 20 && lPt < 30) lPtId = 2;
166	pharris	1.4	if(lPt > 30 ) lPtId = 3;
167
168			int lEtaId = 0;
169			if(fabs(iJet->Eta()) > 2.5 && fabs(iJet->Eta()) < 2.75) lEtaId = 1;
170			if(fabs(iJet->Eta()) > 2.75 && fabs(iJet->Eta()) < 3.0 ) lEtaId = 2;
171			if(fabs(iJet->Eta()) > 3.0 && fabs(iJet->Eta()) < 5.0 ) lEtaId = 3;
172
173			double lMVACut = fMVACut[lPtId][lEtaId];
174			if(lMVA < lMVACut) return false;
175	pharris	1.1	return true;
176	pharris	1.4	//if( fabs(JetTools::impactParameter(iJet,iVertex,true)) < 0.2) return true;
177			//if(correctedPt(iJet,iJetCorrector,iPileupEnergyDensity) < fJetPtMin && iJet->TrackCountingHighEffBJetTagsDisc() == -100) return false;
178			//This line is a bug in the Met training
179			//if(lMVA < -0.8) return false;
180			//if(lMVA < -0.5 && fabs(iJet->Eta()) > 3.0) return false;
181	pharris	1.1	}
182			//--------------------------------------------------------------------------------------------------
183	pharris	1.4	Bool_t JetIDMVA::pass(const PFJet iJet,const Vertex iVertex,const VertexCol *iVertices) {
184	pharris	1.2	if(!JetTools::passPFLooseId(iJet)) return false;
185	pharris	1.5	if(iJet->Pt() < fJetPtMin) return false;
186	pharris	1.7	if(fabs(iJet->Eta()) > 4.99) return false;
187	pharris	1.4	double lMVA = MVAValue(iJet,iVertex,iVertices);
188
189			int lPtId = 0;
190			if(iJet->Pt() > 10 && iJet->Pt() < 20) lPtId = 1;
191			if(iJet->Pt() > 20 && iJet->Pt() < 30) lPtId = 2;
192			if(iJet->Pt() > 30 ) lPtId = 3;
193
194			int lEtaId = 0;
195			if(fabs(iJet->Eta()) > 2.5 && fabs(iJet->Eta()) < 2.75) lEtaId = 1;
196			if(fabs(iJet->Eta()) > 2.75 && fabs(iJet->Eta()) < 3.0 ) lEtaId = 2;
197			if(fabs(iJet->Eta()) > 3.0 && fabs(iJet->Eta()) < 5.0 ) lEtaId = 3;
198
199			double lMVACut = fMVACut[lPtId][lEtaId];
200			if(lMVA < lMVACut) return false;
201	pharris	1.2	return true;
202	pharris	1.4	//if(lMVA < -0.8) return false;
203			//if(lMVA < -0.5 && fabs(iJet->Eta()) > 3.0) return false;
204			//if(iJet->Pt() < fJetPtMin && iJet->TrackCountingHighEffBJetTagsDisc() == -100) return false; //This line is a bug in the Met training
205			//if( fabs(JetTools::impactParameter(iJet,iVertex,true)) < 0.2) return true;
206	pharris	1.2	}
207			//--------------------------------------------------------------------------------------------------
208	pharris	1.4	Double_t JetIDMVA::MVAValue(const PFJet iJet,const Vertex iVertex,const VertexCol *iVertices, //Vertex here is the PV
209	pharris	1.1	FactorizedJetCorrector *iJetCorrector,
210			const PileupEnergyDensityCol *iPileupEnergyDensity,
211			Bool_t printDebug) {
212
213			if (!fIsInitialized) {
214			std::cout << "Error: JetIDMVA not properly initialized.\n";
215			return -9999;
216			}
217			if(!JetTools::passPFLooseId(iJet)) return -2.;
218
219			//set all input variables
220	pharris	1.4	fNVtx = iVertices->GetEntries();
221	pharris	1.1	fJPt1 = correctedPt(iJet,iJetCorrector,iPileupEnergyDensity);
222			fJEta1 = iJet->RawMom().Eta();
223			fJPhi1 = iJet->RawMom().Phi();
224	pharris	1.4	fJD01 = JetTools::impactParameter(iJet,iVertex);
225			fJDZ1 = JetTools::impactParameter(iJet,iVertex,true);
226			fBeta = JetTools::Beta(iJet,iVertex,fDZCut);
227			fBetaStar = JetTools::betaStar(iJet,iVertex,iVertices,fDZCut);
228			fNCharged = iJet->ChargedMultiplicity();
229			fNNeutrals = iJet->NeutralMultiplicity();
230
231			fDRMean = JetTools::dRMean(iJet,-1);
232			fFrac01 = JetTools::frac (iJet,0.1,0. ,-1);
233			fFrac02 = JetTools::frac (iJet,0.2,0.1,-1);
234			fFrac03 = JetTools::frac (iJet,0.3,0.2,-1);
235			fFrac04 = JetTools::frac (iJet,0.4,0.3,-1);
236			fFrac05 = JetTools::frac (iJet,0.5,0.4,-1);
237
238			double lMVA = 0;
239			if(fJPt1 < 10) lMVA = fReader->EvaluateMVA( fLowPtMethodName );
240			if(fJPt1 > 10) lMVA = fReader->EvaluateMVA( fHighPtMethodName );
241	pharris	1.1	if (printDebug == kTRUE) {
242			std::cout << "Debug Jet MVA: "
243	pharris	1.4	<< fNVtx << " "
244			<< fJPt1 << " "
245			<< fJEta1 << " "
246			<< fJPhi1 << " "
247			<< fJD01 << " "
248			<< fJDZ1 << " "
249			<< fBeta << " "
250			<< fBetaStar << " "
251			<< fNCharged << " "
252			<< fNNeutrals << " "
253			<< fDRMean << " "
254			<< fFrac01 << " "
255			<< fFrac02 << " "
256			<< fFrac03 << " "
257			<< fFrac04 << " "
258			<< fFrac05
259	pharris	1.1	<< " === : === "
260			<< lMVA << " "
261			<< std::endl;
262			}
263
264			return lMVA;
265			}
266	pharris	1.4	Double_t JetIDMVA::MVAValue(const PFJet iJet,const Vertex iVertex, const VertexCol *iVertices,//Vertex here is the PV
267	pharris	1.2	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.2	fJPt1 = iJet->Pt();
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			fFrac01 = JetTools::frac (iJet,0.1,0. ,-1);
289			fFrac02 = JetTools::frac (iJet,0.2,0.1,-1);
290			fFrac03 = JetTools::frac (iJet,0.3,0.2,-1);
291			fFrac04 = JetTools::frac (iJet,0.4,0.3,-1);
292			fFrac05 = JetTools::frac (iJet,0.5,0.4,-1);
293
294			double lMVA = 0;
295			if(fJPt1 < 10) lMVA = fReader->EvaluateMVA( fLowPtMethodName );
296			if(fJPt1 > 10) lMVA = fReader->EvaluateMVA( fHighPtMethodName );
297	pharris	1.5
298	pharris	1.2	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			<< fFrac05
316	pharris	1.2	<< " === : === "
317			<< lMVA << " "
318			<< std::endl;
319			}
320
321			return lMVA;
322			}
323	pharris	1.1	Double_t JetIDMVA::correctedPt(const PFJet iJet, FactorizedJetCorrector iJetCorrector,const PileupEnergyDensityCol *iPUEnergyDensity) {
324			const FourVectorM rawMom = iJet->RawMom();
325			iJetCorrector->setJetEta(rawMom.Eta());
326			iJetCorrector->setJetPt (rawMom.Pt());
327			iJetCorrector->setJetPhi(rawMom.Phi());
328			iJetCorrector->setJetE (rawMom.E());
329			iJetCorrector->setRho (iPUEnergyDensity->At(0)->RhoHighEta());
330			iJetCorrector->setJetA (iJet->JetArea());
331			iJetCorrector->setJetEMF(-99.0);
332			Double_t correction = iJetCorrector->getCorrection();
333			return rawMom.Pt()*correction;
334			}