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";
  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 true; //==> Castor
  //if(iJet->Pt() > 50)                                return true; //==> we can raise this 
  
  double lMVA = MVAValue   (iJet,iVertex,iVertices,iJetCorrector,iPileupEnergyDensity);
  double lPt  = correctedPt(iJet,                  iJetCorrector,iPileupEnergyDensity);

  int lPtId = 0; 
  if(lPt > 10 && iJet->Pt() < 20) lPtId = 1;
  if(lPt > 20 && iJet->Pt() < 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 true; //==> Castor
  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.5
Committed:	Fri Apr 13 14:24:00 2012 UTC (13 years ago) by pharris
Content type:	text/plain
Branch:	MAIN
Changes since 1.4:	+10 -9 lines
Log Message:	* empty log message *
#	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			std::vector<double> lPt010 = lConfig.getParameter<std::vector<double> >(("Pt010_" +lCutType).c_str());
90			std::vector<double> lPt1020 = lConfig.getParameter<std::vector<double> >(("Pt1020_"+lCutType).c_str());
91			std::vector<double> lPt2030 = lConfig.getParameter<std::vector<double> >(("Pt2030_"+lCutType).c_str());
92			std::vector<double> lPt3050 = lConfig.getParameter<std::vector<double> >(("Pt3050_"+lCutType).c_str());
93			for(int i0 = 0; i0 < 4; i0++) fMVACut[0][i0] = lPt010 [i0];
94			for(int i0 = 0; i0 < 4; i0++) fMVACut[1][i0] = lPt1020[i0];
95			for(int i0 = 0; i0 < 4; i0++) fMVACut[2][i0] = lPt2030[i0];
96			for(int i0 = 0; i0 < 4; i0++) fMVACut[3][i0] = lPt3050[i0];
97	pharris	1.5	//std::cout << " Working Points : << " << std::endl;
98			//for(int i0 = 0; i0 < 4; i0++) for(int i1 = 0; i1 < 4; i1++)
99			// std::cout << " ==> " << i0 << " -- " << i1 << " -- " << fMVACut[i0][i1] << std::endl;
100	pharris	1.1	}
101
102			//--------------------------------------------------------------------------------------------------
103			Double_t JetIDMVA::MVAValue(
104			Float_t iNPV ,
105			Float_t iJPt1 ,
106			Float_t iJEta1 ,
107			Float_t iJPhi1 ,
108			Float_t iJD01 ,
109			Float_t iJDZ1 ,
110	pharris	1.4	Float_t iBeta ,
111			Float_t iBetaStar,
112			Float_t iNCharged,
113			Float_t iNNeutrals,
114			Float_t iDRMean ,
115			Float_t iFrac01 ,
116			Float_t iFrac02 ,
117			Float_t iFrac03 ,
118			Float_t iFrac04 ,
119			Float_t iFrac05
120	pharris	1.1	){
121
122	pharris	1.4	if(!fIsInitialized) {
123	pharris	1.1	std::cout << "Error: JetIDMVA not properly initialized.\n";
124			return -9999;
125			}
126
127	pharris	1.4	fNVtx = iNPV;
128			fJPt1 = iJPt1;
129			fJEta1 = iJEta1;
130			fJPhi1 = fJPhi1;
131			fJD01 = iJD01;
132			fJDZ1 = iJDZ1;
133			fBeta = iBeta;
134			fBetaStar = iBetaStar;
135			fNCharged = iNCharged;
136			fNNeutrals = iNNeutrals;
137			fDRMean = iDRMean;
138			fFrac01 = iFrac01;
139			fFrac02 = iFrac02;
140			fFrac03 = iFrac03;
141			fFrac04 = iFrac04;
142			fFrac05 = iFrac05;
143	pharris	1.1
144			Double_t lMVA = -9999;
145	pharris	1.4	if(iJPt1 < 10) lMVA = fReader->EvaluateMVA( fLowPtMethodName );
146			if(iJPt1 > 10) lMVA = fReader->EvaluateMVA( fHighPtMethodName );
147	pharris	1.1
148			return lMVA;
149			}
150			//--------------------------------------------------------------------------------------------------
151	pharris	1.4	Bool_t JetIDMVA::pass(const PFJet iJet,const Vertex iVertex,const VertexCol *iVertices,
152	pharris	1.1	FactorizedJetCorrector *iJetCorrector,
153			const PileupEnergyDensityCol *iPileupEnergyDensity) {
154	pharris	1.4
155	pharris	1.1	if(!JetTools::passPFLooseId(iJet)) return false;
156	pharris	1.4	if(iJet->Pt() < fJetPtMin) return false;
157	pharris	1.5	if(fabs(iJet->Eta()) > 4.99) return true; //==> Castor
158			//if(iJet->Pt() > 50) return true; //==> we can raise this
159
160	pharris	1.4	double lMVA = MVAValue (iJet,iVertex,iVertices,iJetCorrector,iPileupEnergyDensity);
161			double lPt = correctedPt(iJet, iJetCorrector,iPileupEnergyDensity);
162
163			int lPtId = 0;
164			if(lPt > 10 && iJet->Pt() < 20) lPtId = 1;
165			if(lPt > 20 && iJet->Pt() < 30) lPtId = 2;
166			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			//if(iJet->Pt() > 50) return true; //==> we can raise this
187			if(fabs(iJet->Eta()) > 4.99) return true; //==> Castor
188	pharris	1.4	double lMVA = MVAValue(iJet,iVertex,iVertices);
189
190			int lPtId = 0;
191			if(iJet->Pt() > 10 && iJet->Pt() < 20) lPtId = 1;
192			if(iJet->Pt() > 20 && iJet->Pt() < 30) lPtId = 2;
193			if(iJet->Pt() > 30 ) lPtId = 3;
194
195			int lEtaId = 0;
196			if(fabs(iJet->Eta()) > 2.5 && fabs(iJet->Eta()) < 2.75) lEtaId = 1;
197			if(fabs(iJet->Eta()) > 2.75 && fabs(iJet->Eta()) < 3.0 ) lEtaId = 2;
198			if(fabs(iJet->Eta()) > 3.0 && fabs(iJet->Eta()) < 5.0 ) lEtaId = 3;
199
200			double lMVACut = fMVACut[lPtId][lEtaId];
201			if(lMVA < lMVACut) return false;
202	pharris	1.2	return true;
203	pharris	1.4	//if(lMVA < -0.8) return false;
204			//if(lMVA < -0.5 && fabs(iJet->Eta()) > 3.0) return false;
205			//if(iJet->Pt() < fJetPtMin && iJet->TrackCountingHighEffBJetTagsDisc() == -100) return false; //This line is a bug in the Met training
206			//if( fabs(JetTools::impactParameter(iJet,iVertex,true)) < 0.2) return true;
207	pharris	1.2	}
208			//--------------------------------------------------------------------------------------------------
209	pharris	1.4	Double_t JetIDMVA::MVAValue(const PFJet iJet,const Vertex iVertex,const VertexCol *iVertices, //Vertex here is the PV
210	pharris	1.1	FactorizedJetCorrector *iJetCorrector,
211			const PileupEnergyDensityCol *iPileupEnergyDensity,
212			Bool_t printDebug) {
213
214			if (!fIsInitialized) {
215			std::cout << "Error: JetIDMVA not properly initialized.\n";
216			return -9999;
217			}
218			if(!JetTools::passPFLooseId(iJet)) return -2.;
219
220			//set all input variables
221	pharris	1.4	fNVtx = iVertices->GetEntries();
222	pharris	1.1	fJPt1 = correctedPt(iJet,iJetCorrector,iPileupEnergyDensity);
223			fJEta1 = iJet->RawMom().Eta();
224			fJPhi1 = iJet->RawMom().Phi();
225	pharris	1.4	fJD01 = JetTools::impactParameter(iJet,iVertex);
226			fJDZ1 = JetTools::impactParameter(iJet,iVertex,true);
227			fBeta = JetTools::Beta(iJet,iVertex,fDZCut);
228			fBetaStar = JetTools::betaStar(iJet,iVertex,iVertices,fDZCut);
229			fNCharged = iJet->ChargedMultiplicity();
230			fNNeutrals = iJet->NeutralMultiplicity();
231
232			fDRMean = JetTools::dRMean(iJet,-1);
233			fFrac01 = JetTools::frac (iJet,0.1,0. ,-1);
234			fFrac02 = JetTools::frac (iJet,0.2,0.1,-1);
235			fFrac03 = JetTools::frac (iJet,0.3,0.2,-1);
236			fFrac04 = JetTools::frac (iJet,0.4,0.3,-1);
237			fFrac05 = JetTools::frac (iJet,0.5,0.4,-1);
238
239			double lMVA = 0;
240			if(fJPt1 < 10) lMVA = fReader->EvaluateMVA( fLowPtMethodName );
241			if(fJPt1 > 10) lMVA = fReader->EvaluateMVA( fHighPtMethodName );
242	pharris	1.1	if (printDebug == kTRUE) {
243			std::cout << "Debug Jet MVA: "
244	pharris	1.4	<< fNVtx << " "
245			<< fJPt1 << " "
246			<< fJEta1 << " "
247			<< fJPhi1 << " "
248			<< fJD01 << " "
249			<< fJDZ1 << " "
250			<< fBeta << " "
251			<< fBetaStar << " "
252			<< fNCharged << " "
253			<< fNNeutrals << " "
254			<< fDRMean << " "
255			<< fFrac01 << " "
256			<< fFrac02 << " "
257			<< fFrac03 << " "
258			<< fFrac04 << " "
259			<< fFrac05
260	pharris	1.1	<< " === : === "
261			<< lMVA << " "
262			<< std::endl;
263			}
264
265			return lMVA;
266			}
267	pharris	1.4	Double_t JetIDMVA::MVAValue(const PFJet iJet,const Vertex iVertex, const VertexCol *iVertices,//Vertex here is the PV
268	pharris	1.2	Bool_t printDebug) {
269
270			if (!fIsInitialized) {
271			std::cout << "Error: JetIDMVA not properly initialized.\n";
272			return -9999;
273			}
274			if(!JetTools::passPFLooseId(iJet)) return -2.;
275
276			//set all input variables
277	pharris	1.4	fNVtx = iVertices->GetEntries();
278	pharris	1.2	fJPt1 = iJet->Pt();
279			fJEta1 = iJet->RawMom().Eta();
280			fJPhi1 = iJet->RawMom().Phi();
281	pharris	1.4	fJD01 = JetTools::impactParameter(iJet,iVertex);
282			fJDZ1 = JetTools::impactParameter(iJet,iVertex,true);
283			fBeta = JetTools::Beta(iJet,iVertex,fDZCut);
284			fBetaStar = JetTools::betaStar(iJet,iVertex,iVertices,fDZCut);
285			fNCharged = iJet->ChargedMultiplicity();
286			fNNeutrals = iJet->NeutralMultiplicity();
287
288			fDRMean = JetTools::dRMean(iJet,-1);
289			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			if(fJPt1 < 10) lMVA = fReader->EvaluateMVA( fLowPtMethodName );
297			if(fJPt1 > 10) lMVA = fReader->EvaluateMVA( fHighPtMethodName );
298	pharris	1.5
299	pharris	1.2	if (printDebug == kTRUE) {
300			std::cout << "Debug Jet MVA: "
301	pharris	1.4	<< fNVtx << " "
302			<< fJPt1 << " "
303			<< fJEta1 << " "
304			<< fJPhi1 << " "
305			<< fJD01 << " "
306			<< fJDZ1 << " "
307			<< fBeta << " "
308			<< fBetaStar << " "
309			<< fNCharged << " "
310			<< fNNeutrals << " "
311			<< fDRMean << " "
312			<< fFrac01 << " "
313			<< fFrac02 << " "
314			<< fFrac03 << " "
315			<< fFrac04 << " "
316			<< fFrac05
317	pharris	1.2	<< " === : === "
318			<< lMVA << " "
319			<< std::endl;
320			}
321
322			return lMVA;
323			}
324	pharris	1.1	Double_t JetIDMVA::correctedPt(const PFJet iJet, FactorizedJetCorrector iJetCorrector,const PileupEnergyDensityCol *iPUEnergyDensity) {
325			const FourVectorM rawMom = iJet->RawMom();
326			iJetCorrector->setJetEta(rawMom.Eta());
327			iJetCorrector->setJetPt (rawMom.Pt());
328			iJetCorrector->setJetPhi(rawMom.Phi());
329			iJetCorrector->setJetE (rawMom.E());
330			iJetCorrector->setRho (iPUEnergyDensity->At(0)->RhoHighEta());
331			iJetCorrector->setJetA (iJet->JetArea());
332			iJetCorrector->setJetEMF(-99.0);
333			Double_t correction = iJetCorrector->getCorrection();
334			return rawMom.Pt()*correction;
335			}