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(iJet->Pt() < fJetPtMin)                         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 > 50)                                return true; //==> we can raise this 
  
  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(iJet->Pt()        > 50)        return true; //==> we can raise this 
  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.7
Committed:	Tue May 1 16:53:38 2012 UTC (13 years ago) by pharris
Content type:	text/plain
Branch:	MAIN
Changes since 1.6:	+7 -7 lines
Log Message:	Fixed Castor Jet Bug
#	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.4	fJetPtMin(0) , //We need to lower this
19			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.4	if(iJet->Pt() < fJetPtMin) return false;
158	pharris	1.7	if(fabs(iJet->Eta()) > 4.99) return false;
159	pharris	1.5
160	pharris	1.4	double lMVA = MVAValue (iJet,iVertex,iVertices,iJetCorrector,iPileupEnergyDensity);
161			double lPt = correctedPt(iJet, iJetCorrector,iPileupEnergyDensity);
162	pharris	1.7	if(lPt > 50) return true; //==> we can raise this
163
164	pharris	1.4	int lPtId = 0;
165	pharris	1.7	if(lPt > 10 && lPt < 20) lPtId = 1;
166			if(lPt > 20 && lPt < 30) lPtId = 2;
167	pharris	1.4	if(lPt > 30 ) lPtId = 3;
168
169			int lEtaId = 0;
170			if(fabs(iJet->Eta()) > 2.5 && fabs(iJet->Eta()) < 2.75) lEtaId = 1;
171			if(fabs(iJet->Eta()) > 2.75 && fabs(iJet->Eta()) < 3.0 ) lEtaId = 2;
172			if(fabs(iJet->Eta()) > 3.0 && fabs(iJet->Eta()) < 5.0 ) lEtaId = 3;
173
174			double lMVACut = fMVACut[lPtId][lEtaId];
175			if(lMVA < lMVACut) return false;
176	pharris	1.1	return true;
177	pharris	1.4	//if( fabs(JetTools::impactParameter(iJet,iVertex,true)) < 0.2) return true;
178			//if(correctedPt(iJet,iJetCorrector,iPileupEnergyDensity) < fJetPtMin && iJet->TrackCountingHighEffBJetTagsDisc() == -100) return false;
179			//This line is a bug in the Met training
180			//if(lMVA < -0.8) return false;
181			//if(lMVA < -0.5 && fabs(iJet->Eta()) > 3.0) return false;
182	pharris	1.1	}
183			//--------------------------------------------------------------------------------------------------
184	pharris	1.4	Bool_t JetIDMVA::pass(const PFJet iJet,const Vertex iVertex,const VertexCol *iVertices) {
185	pharris	1.2	if(!JetTools::passPFLooseId(iJet)) return false;
186	pharris	1.5	if(iJet->Pt() < fJetPtMin) return false;
187	pharris	1.7	if(iJet->Pt() > 50) return true; //==> we can raise this
188			if(fabs(iJet->Eta()) > 4.99) return false;
189	pharris	1.4	double lMVA = MVAValue(iJet,iVertex,iVertices);
190
191			int lPtId = 0;
192			if(iJet->Pt() > 10 && iJet->Pt() < 20) lPtId = 1;
193			if(iJet->Pt() > 20 && iJet->Pt() < 30) lPtId = 2;
194			if(iJet->Pt() > 30 ) lPtId = 3;
195
196			int lEtaId = 0;
197			if(fabs(iJet->Eta()) > 2.5 && fabs(iJet->Eta()) < 2.75) lEtaId = 1;
198			if(fabs(iJet->Eta()) > 2.75 && fabs(iJet->Eta()) < 3.0 ) lEtaId = 2;
199			if(fabs(iJet->Eta()) > 3.0 && fabs(iJet->Eta()) < 5.0 ) lEtaId = 3;
200
201			double lMVACut = fMVACut[lPtId][lEtaId];
202			if(lMVA < lMVACut) return false;
203	pharris	1.2	return true;
204	pharris	1.4	//if(lMVA < -0.8) return false;
205			//if(lMVA < -0.5 && fabs(iJet->Eta()) > 3.0) return false;
206			//if(iJet->Pt() < fJetPtMin && iJet->TrackCountingHighEffBJetTagsDisc() == -100) return false; //This line is a bug in the Met training
207			//if( fabs(JetTools::impactParameter(iJet,iVertex,true)) < 0.2) return true;
208	pharris	1.2	}
209			//--------------------------------------------------------------------------------------------------
210	pharris	1.4	Double_t JetIDMVA::MVAValue(const PFJet iJet,const Vertex iVertex,const VertexCol *iVertices, //Vertex here is the PV
211	pharris	1.1	FactorizedJetCorrector *iJetCorrector,
212			const PileupEnergyDensityCol *iPileupEnergyDensity,
213			Bool_t printDebug) {
214
215			if (!fIsInitialized) {
216			std::cout << "Error: JetIDMVA not properly initialized.\n";
217			return -9999;
218			}
219			if(!JetTools::passPFLooseId(iJet)) return -2.;
220
221			//set all input variables
222	pharris	1.4	fNVtx = iVertices->GetEntries();
223	pharris	1.1	fJPt1 = correctedPt(iJet,iJetCorrector,iPileupEnergyDensity);
224			fJEta1 = iJet->RawMom().Eta();
225			fJPhi1 = iJet->RawMom().Phi();
226	pharris	1.4	fJD01 = JetTools::impactParameter(iJet,iVertex);
227			fJDZ1 = JetTools::impactParameter(iJet,iVertex,true);
228			fBeta = JetTools::Beta(iJet,iVertex,fDZCut);
229			fBetaStar = JetTools::betaStar(iJet,iVertex,iVertices,fDZCut);
230			fNCharged = iJet->ChargedMultiplicity();
231			fNNeutrals = iJet->NeutralMultiplicity();
232
233			fDRMean = JetTools::dRMean(iJet,-1);
234			fFrac01 = JetTools::frac (iJet,0.1,0. ,-1);
235			fFrac02 = JetTools::frac (iJet,0.2,0.1,-1);
236			fFrac03 = JetTools::frac (iJet,0.3,0.2,-1);
237			fFrac04 = JetTools::frac (iJet,0.4,0.3,-1);
238			fFrac05 = JetTools::frac (iJet,0.5,0.4,-1);
239
240			double lMVA = 0;
241			if(fJPt1 < 10) lMVA = fReader->EvaluateMVA( fLowPtMethodName );
242			if(fJPt1 > 10) lMVA = fReader->EvaluateMVA( fHighPtMethodName );
243	pharris	1.1	if (printDebug == kTRUE) {
244			std::cout << "Debug Jet MVA: "
245	pharris	1.4	<< fNVtx << " "
246			<< fJPt1 << " "
247			<< fJEta1 << " "
248			<< fJPhi1 << " "
249			<< fJD01 << " "
250			<< fJDZ1 << " "
251			<< fBeta << " "
252			<< fBetaStar << " "
253			<< fNCharged << " "
254			<< fNNeutrals << " "
255			<< fDRMean << " "
256			<< fFrac01 << " "
257			<< fFrac02 << " "
258			<< fFrac03 << " "
259			<< fFrac04 << " "
260			<< fFrac05
261	pharris	1.1	<< " === : === "
262			<< lMVA << " "
263			<< std::endl;
264			}
265
266			return lMVA;
267			}
268	pharris	1.4	Double_t JetIDMVA::MVAValue(const PFJet iJet,const Vertex iVertex, const VertexCol *iVertices,//Vertex here is the PV
269	pharris	1.2	Bool_t printDebug) {
270
271			if (!fIsInitialized) {
272			std::cout << "Error: JetIDMVA not properly initialized.\n";
273			return -9999;
274			}
275			if(!JetTools::passPFLooseId(iJet)) return -2.;
276
277			//set all input variables
278	pharris	1.4	fNVtx = iVertices->GetEntries();
279	pharris	1.2	fJPt1 = iJet->Pt();
280			fJEta1 = iJet->RawMom().Eta();
281			fJPhi1 = iJet->RawMom().Phi();
282	pharris	1.4	fJD01 = JetTools::impactParameter(iJet,iVertex);
283			fJDZ1 = JetTools::impactParameter(iJet,iVertex,true);
284			fBeta = JetTools::Beta(iJet,iVertex,fDZCut);
285			fBetaStar = JetTools::betaStar(iJet,iVertex,iVertices,fDZCut);
286			fNCharged = iJet->ChargedMultiplicity();
287			fNNeutrals = iJet->NeutralMultiplicity();
288
289			fDRMean = JetTools::dRMean(iJet,-1);
290			fFrac01 = JetTools::frac (iJet,0.1,0. ,-1);
291			fFrac02 = JetTools::frac (iJet,0.2,0.1,-1);
292			fFrac03 = JetTools::frac (iJet,0.3,0.2,-1);
293			fFrac04 = JetTools::frac (iJet,0.4,0.3,-1);
294			fFrac05 = JetTools::frac (iJet,0.5,0.4,-1);
295
296			double lMVA = 0;
297			if(fJPt1 < 10) lMVA = fReader->EvaluateMVA( fLowPtMethodName );
298			if(fJPt1 > 10) lMVA = fReader->EvaluateMVA( fHighPtMethodName );
299	pharris	1.5
300	pharris	1.2	if (printDebug == kTRUE) {
301			std::cout << "Debug Jet MVA: "
302	pharris	1.4	<< fNVtx << " "
303			<< fJPt1 << " "
304			<< fJEta1 << " "
305			<< fJPhi1 << " "
306			<< fJD01 << " "
307			<< fJDZ1 << " "
308			<< fBeta << " "
309			<< fBetaStar << " "
310			<< fNCharged << " "
311			<< fNNeutrals << " "
312			<< fDRMean << " "
313			<< fFrac01 << " "
314			<< fFrac02 << " "
315			<< fFrac03 << " "
316			<< fFrac04 << " "
317			<< fFrac05
318	pharris	1.2	<< " === : === "
319			<< lMVA << " "
320			<< std::endl;
321			}
322
323			return lMVA;
324			}
325	pharris	1.1	Double_t JetIDMVA::correctedPt(const PFJet iJet, FactorizedJetCorrector iJetCorrector,const PileupEnergyDensityCol *iPUEnergyDensity) {
326			const FourVectorM rawMom = iJet->RawMom();
327			iJetCorrector->setJetEta(rawMom.Eta());
328			iJetCorrector->setJetPt (rawMom.Pt());
329			iJetCorrector->setJetPhi(rawMom.Phi());
330			iJetCorrector->setJetE (rawMom.E());
331			iJetCorrector->setRho (iPUEnergyDensity->At(0)->RhoHighEta());
332			iJetCorrector->setJetA (iJet->JetArea());
333			iJetCorrector->setJetEMF(-99.0);
334			Double_t correction = iJetCorrector->getCorrection();
335			return rawMom.Pt()*correction;
336			}