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,
                      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);
  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,
                               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.10
Committed:	Sat May 12 07:26:56 2012 UTC (12 years, 11 months ago) by ceballos
Content type:	text/plain
Branch:	MAIN
CVS Tags:	Mit_028, Mit_027, Mit_027a
Changes since 1.9:	+30 -9 lines
Log Message:	adding rho functionality
#	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	ceballos	1.10	const PileupEnergyDensityCol *iPileupEnergyDensity,
155			RhoUtilities::RhoType type) {
156	pharris	1.4
157	pharris	1.1	if(!JetTools::passPFLooseId(iJet)) 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	ceballos	1.10	double lPt = correctedPt(iJet, iJetCorrector,iPileupEnergyDensity,type);
162	pharris	1.8	if(lPt < fJetPtMin) return false;
163	pharris	1.7
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			}
178			//--------------------------------------------------------------------------------------------------
179	pharris	1.4	Bool_t JetIDMVA::pass(const PFJet iJet,const Vertex iVertex,const VertexCol *iVertices) {
180	pharris	1.2	if(!JetTools::passPFLooseId(iJet)) return false;
181	pharris	1.5	if(iJet->Pt() < fJetPtMin) return false;
182	pharris	1.7	if(fabs(iJet->Eta()) > 4.99) return false;
183	pharris	1.4	double lMVA = MVAValue(iJet,iVertex,iVertices);
184
185			int lPtId = 0;
186			if(iJet->Pt() > 10 && iJet->Pt() < 20) lPtId = 1;
187			if(iJet->Pt() > 20 && iJet->Pt() < 30) lPtId = 2;
188			if(iJet->Pt() > 30 ) lPtId = 3;
189
190			int lEtaId = 0;
191			if(fabs(iJet->Eta()) > 2.5 && fabs(iJet->Eta()) < 2.75) lEtaId = 1;
192			if(fabs(iJet->Eta()) > 2.75 && fabs(iJet->Eta()) < 3.0 ) lEtaId = 2;
193			if(fabs(iJet->Eta()) > 3.0 && fabs(iJet->Eta()) < 5.0 ) lEtaId = 3;
194
195			double lMVACut = fMVACut[lPtId][lEtaId];
196			if(lMVA < lMVACut) return false;
197	pharris	1.2	return true;
198	pharris	1.4	//if(lMVA < -0.8) return false;
199			//if(lMVA < -0.5 && fabs(iJet->Eta()) > 3.0) return false;
200			//if(iJet->Pt() < fJetPtMin && iJet->TrackCountingHighEffBJetTagsDisc() == -100) return false; //This line is a bug in the Met training
201			//if( fabs(JetTools::impactParameter(iJet,iVertex,true)) < 0.2) return true;
202	pharris	1.2	}
203			//--------------------------------------------------------------------------------------------------
204	pharris	1.4	Double_t JetIDMVA::MVAValue(const PFJet iJet,const Vertex iVertex,const VertexCol *iVertices, //Vertex here is the PV
205	pharris	1.1	FactorizedJetCorrector *iJetCorrector,
206			const PileupEnergyDensityCol *iPileupEnergyDensity,
207			Bool_t printDebug) {
208
209			if (!fIsInitialized) {
210			std::cout << "Error: JetIDMVA not properly initialized.\n";
211			return -9999;
212			}
213			if(!JetTools::passPFLooseId(iJet)) return -2.;
214
215			//set all input variables
216	pharris	1.4	fNVtx = iVertices->GetEntries();
217	pharris	1.1	fJPt1 = correctedPt(iJet,iJetCorrector,iPileupEnergyDensity);
218			fJEta1 = iJet->RawMom().Eta();
219			fJPhi1 = iJet->RawMom().Phi();
220	pharris	1.4	fJD01 = JetTools::impactParameter(iJet,iVertex);
221			fJDZ1 = JetTools::impactParameter(iJet,iVertex,true);
222			fBeta = JetTools::Beta(iJet,iVertex,fDZCut);
223			fBetaStar = JetTools::betaStar(iJet,iVertex,iVertices,fDZCut);
224			fNCharged = iJet->ChargedMultiplicity();
225			fNNeutrals = iJet->NeutralMultiplicity();
226
227			fDRMean = JetTools::dRMean(iJet,-1);
228			fFrac01 = JetTools::frac (iJet,0.1,0. ,-1);
229			fFrac02 = JetTools::frac (iJet,0.2,0.1,-1);
230			fFrac03 = JetTools::frac (iJet,0.3,0.2,-1);
231			fFrac04 = JetTools::frac (iJet,0.4,0.3,-1);
232			fFrac05 = JetTools::frac (iJet,0.5,0.4,-1);
233
234			double lMVA = 0;
235			if(fJPt1 < 10) lMVA = fReader->EvaluateMVA( fLowPtMethodName );
236			if(fJPt1 > 10) lMVA = fReader->EvaluateMVA( fHighPtMethodName );
237	pharris	1.1	if (printDebug == kTRUE) {
238			std::cout << "Debug Jet MVA: "
239	pharris	1.4	<< fNVtx << " "
240			<< fJPt1 << " "
241			<< fJEta1 << " "
242			<< fJPhi1 << " "
243			<< fJD01 << " "
244			<< fJDZ1 << " "
245			<< fBeta << " "
246			<< fBetaStar << " "
247			<< fNCharged << " "
248			<< fNNeutrals << " "
249			<< fDRMean << " "
250			<< fFrac01 << " "
251			<< fFrac02 << " "
252			<< fFrac03 << " "
253			<< fFrac04 << " "
254			<< fFrac05
255	pharris	1.1	<< " === : === "
256			<< lMVA << " "
257			<< std::endl;
258			}
259
260			return lMVA;
261			}
262	pharris	1.4	Double_t JetIDMVA::MVAValue(const PFJet iJet,const Vertex iVertex, const VertexCol *iVertices,//Vertex here is the PV
263	pharris	1.2	Bool_t printDebug) {
264
265			if (!fIsInitialized) {
266			std::cout << "Error: JetIDMVA not properly initialized.\n";
267			return -9999;
268			}
269			if(!JetTools::passPFLooseId(iJet)) return -2.;
270
271			//set all input variables
272	pharris	1.4	fNVtx = iVertices->GetEntries();
273	pharris	1.2	fJPt1 = iJet->Pt();
274			fJEta1 = iJet->RawMom().Eta();
275			fJPhi1 = iJet->RawMom().Phi();
276	pharris	1.4	fJD01 = JetTools::impactParameter(iJet,iVertex);
277			fJDZ1 = JetTools::impactParameter(iJet,iVertex,true);
278			fBeta = JetTools::Beta(iJet,iVertex,fDZCut);
279			fBetaStar = JetTools::betaStar(iJet,iVertex,iVertices,fDZCut);
280			fNCharged = iJet->ChargedMultiplicity();
281			fNNeutrals = iJet->NeutralMultiplicity();
282
283			fDRMean = JetTools::dRMean(iJet,-1);
284			fFrac01 = JetTools::frac (iJet,0.1,0. ,-1);
285			fFrac02 = JetTools::frac (iJet,0.2,0.1,-1);
286			fFrac03 = JetTools::frac (iJet,0.3,0.2,-1);
287			fFrac04 = JetTools::frac (iJet,0.4,0.3,-1);
288			fFrac05 = JetTools::frac (iJet,0.5,0.4,-1);
289
290			double lMVA = 0;
291			if(fJPt1 < 10) lMVA = fReader->EvaluateMVA( fLowPtMethodName );
292			if(fJPt1 > 10) lMVA = fReader->EvaluateMVA( fHighPtMethodName );
293	pharris	1.5
294	pharris	1.2	if (printDebug == kTRUE) {
295			std::cout << "Debug Jet MVA: "
296	pharris	1.4	<< fNVtx << " "
297			<< fJPt1 << " "
298			<< fJEta1 << " "
299			<< fJPhi1 << " "
300			<< fJD01 << " "
301			<< fJDZ1 << " "
302			<< fBeta << " "
303			<< fBetaStar << " "
304			<< fNCharged << " "
305			<< fNNeutrals << " "
306			<< fDRMean << " "
307			<< fFrac01 << " "
308			<< fFrac02 << " "
309			<< fFrac03 << " "
310			<< fFrac04 << " "
311			<< fFrac05
312	pharris	1.2	<< " === : === "
313			<< lMVA << " "
314			<< std::endl;
315			}
316
317			return lMVA;
318			}
319	ceballos	1.10	Double_t JetIDMVA::correctedPt(const PFJet iJet, FactorizedJetCorrector iJetCorrector,
320			const PileupEnergyDensityCol *iPUEnergyDensity,
321			RhoUtilities::RhoType type) {
322			Double_t Rho = 0.0;
323			switch(type) {
324			case RhoUtilities::MIT_RHO_VORONOI_HIGH_ETA:
325			Rho = iPUEnergyDensity->At(0)->Rho();
326			break;
327			case RhoUtilities::MIT_RHO_VORONOI_LOW_ETA:
328			Rho = iPUEnergyDensity->At(0)->RhoLowEta();
329			break;
330			case RhoUtilities::MIT_RHO_RANDOM_HIGH_ETA:
331			Rho = iPUEnergyDensity->At(0)->RhoRandom();
332			break;
333			case RhoUtilities::MIT_RHO_RANDOM_LOW_ETA:
334			Rho = iPUEnergyDensity->At(0)->RhoRandomLowEta();
335			break;
336			case RhoUtilities::CMS_RHO_RHOKT6PFJETS:
337			Rho = iPUEnergyDensity->At(0)->RhoKt6PFJets();
338			break;
339			default:
340			// use the old default
341			Rho = iPUEnergyDensity->At(0)->Rho();
342			break;
343			}
344
345	pharris	1.1	const FourVectorM rawMom = iJet->RawMom();
346			iJetCorrector->setJetEta(rawMom.Eta());
347			iJetCorrector->setJetPt (rawMom.Pt());
348			iJetCorrector->setJetPhi(rawMom.Phi());
349			iJetCorrector->setJetE (rawMom.E());
350	ceballos	1.10	iJetCorrector->setRho (Rho);
351	pharris	1.1	iJetCorrector->setJetA (iJet->JetArea());
352			iJetCorrector->setJetEMF(-99.0);
353			Double_t correction = iJetCorrector->getCorrection();
354			return rawMom.Pt()*correction;
355			}