| 1 |
#include <vector>
|
| 2 |
#include <TFile.h>
|
| 3 |
#include <TRandom3.h>
|
| 4 |
#include "TMVA/Tools.h"
|
| 5 |
#include "TMVA/Reader.h"
|
| 6 |
#include "FWCore/PythonParameterSet/interface/MakeParameterSets.h"
|
| 7 |
#include "FWCore/ParameterSet/interface/ParameterSet.h"
|
| 8 |
#include "MitPhysics/Utils/interface/TauIsoMVA.h"
|
| 9 |
#include "MitCommon/MathTools/interface/MathUtils.h"
|
| 10 |
#include "CondFormats/EgammaObjects/interface/GBRForest.h"
|
| 11 |
#include "Cintex/Cintex.h"
|
| 12 |
#include <utility>
|
| 13 |
|
| 14 |
ClassImp(mithep::TauIsoMVA)
|
| 15 |
|
| 16 |
using namespace mithep;
|
| 17 |
|
| 18 |
//--------------------------------------------------------------------------------------------------
|
| 19 |
TauIsoMVA::TauIsoMVA() {
|
| 20 |
fReader = 0;
|
| 21 |
}
|
| 22 |
//--------------------------------------------------------------------------------------------------
|
| 23 |
TauIsoMVA::~TauIsoMVA() {
|
| 24 |
delete fReader;
|
| 25 |
}
|
| 26 |
//--------------------------------------------------------------------------------------------------
|
| 27 |
void TauIsoMVA::Initialize( TString iWeightFile) {
|
| 28 |
fReader = new TMVA::Reader("!Color:!Silent");
|
| 29 |
fReader->AddVariable("tau_nchiso", (Float_t *)0);
|
| 30 |
fReader->AddVariable("tau_ngiso", (Float_t *)0);
|
| 31 |
fReader->AddVariable("tau_nneuiso", (Float_t *)0);
|
| 32 |
fReader->AddVariable("ring_ch_0*tau_pt", (Float_t *)0);
|
| 33 |
fReader->AddVariable("ring_ch_1*tau_pt", (Float_t *)0);
|
| 34 |
fReader->AddVariable("ring_ch_2*tau_pt", (Float_t *)0);
|
| 35 |
fReader->AddVariable("ring_ch_3*tau_pt", (Float_t *)0);
|
| 36 |
fReader->AddVariable("ring_ch_4*tau_pt", (Float_t *)0);
|
| 37 |
fReader->AddVariable("ring_g_0*tau_pt", (Float_t *)0);
|
| 38 |
fReader->AddVariable("ring_g_1*tau_pt", (Float_t *)0);
|
| 39 |
fReader->AddVariable("ring_g_2*tau_pt", (Float_t *)0);
|
| 40 |
fReader->AddVariable("ring_g_3*tau_pt", (Float_t *)0);
|
| 41 |
fReader->AddVariable("ring_g_4*tau_pt", (Float_t *)0);
|
| 42 |
fReader->AddVariable("ring_neu_0*tau_pt", (Float_t *)0);
|
| 43 |
fReader->AddVariable("ring_neu_1*tau_pt", (Float_t *)0);
|
| 44 |
fReader->AddVariable("ring_neu_2*tau_pt", (Float_t *)0);
|
| 45 |
fReader->AddVariable("ring_neu_3*tau_pt", (Float_t *)0);
|
| 46 |
fReader->AddVariable("ring_neu_4*tau_pt", (Float_t *)0);
|
| 47 |
fReader->AddVariable("shape_ch_eta", (Float_t *)0);
|
| 48 |
fReader->AddVariable("shape_ch_phi", (Float_t *)0);
|
| 49 |
fReader->AddVariable("shape_ch_etaeta", (Float_t *)0);
|
| 50 |
fReader->AddVariable("shape_ch_phiphi", (Float_t *)0);
|
| 51 |
fReader->AddVariable("shape_ch_etaphi", (Float_t *)0);
|
| 52 |
fReader->AddVariable("shape_g_eta", (Float_t *)0);
|
| 53 |
fReader->AddVariable("shape_g_phi", (Float_t *)0);
|
| 54 |
fReader->AddVariable("shape_g_etaeta", (Float_t *)0);
|
| 55 |
fReader->AddVariable("shape_g_phiphi", (Float_t *)0);
|
| 56 |
fReader->AddVariable("shape_g_etaphi", (Float_t *)0);
|
| 57 |
fReader->AddVariable("shape_neu_eta", (Float_t *)0);
|
| 58 |
fReader->AddVariable("shape_neu_phi", (Float_t *)0);
|
| 59 |
fReader->AddVariable("shape_neu_etaeta", (Float_t *)0);
|
| 60 |
fReader->AddVariable("shape_neu_phiphi", (Float_t *)0);
|
| 61 |
fReader->AddVariable("shape_neu_etaphi", (Float_t *)0);
|
| 62 |
fReader->AddVariable("rho", (Float_t *)0);
|
| 63 |
fReader->AddSpectator("tau_pt", (Float_t *)0);
|
| 64 |
fReader->AddSpectator("tau_eta", (Float_t *)0);
|
| 65 |
fReader->AddSpectator("tau_iso", (Float_t *)0);
|
| 66 |
fReader->AddSpectator("gen_pt", (Float_t *)0);
|
| 67 |
fReader->AddSpectator("jet_pt", (Float_t *)0);
|
| 68 |
fReader->AddSpectator("pv", (Float_t *)0);
|
| 69 |
fReader->BookMVA ("BDTG", iWeightFile);
|
| 70 |
fGBR = false;
|
| 71 |
}
|
| 72 |
void TauIsoMVA::InitializeGBR( TString iWeightFile) {
|
| 73 |
ROOT::Cintex::Cintex::Enable();
|
| 74 |
TFile *lForest = new TFile(iWeightFile,"READ");
|
| 75 |
fGBRReader = (GBRForest*)lForest->Get("gbrfTauIso");
|
| 76 |
lForest->Close();
|
| 77 |
fGBR = true;
|
| 78 |
}
|
| 79 |
double TauIsoMVA::MVAValue(const PFTau *iTau,double iRho) {
|
| 80 |
IsoRings lRings = computeIsoRings(iTau);
|
| 81 |
std::vector<float> mvainput = lRings.getVector();
|
| 82 |
mvainput.push_back( iRho);
|
| 83 |
if(!fGBR) mvainput.insert(mvainput.end(), 6, 0);
|
| 84 |
if(!fGBR) return fReader->EvaluateMVA(mvainput, "BDTG");
|
| 85 |
return fGBRReader->GetClassifier(&mvainput[0]);
|
| 86 |
}
|
| 87 |
TauIsoMVA::IsoRings TauIsoMVA::computeIsoRings(const PFTau *iTau) {
|
| 88 |
std::vector<int> niso (3);
|
| 89 |
std::vector<std::vector<float> > rings (3, std::vector<float>(5));
|
| 90 |
std::vector<std::vector<float> > shapes (3, std::vector<float>(5));
|
| 91 |
std::vector<float> isoptsum(3);
|
| 92 |
|
| 93 |
for(UInt_t i0 = 0; i0 < iTau->NIsoPFCandS(); i0++) {
|
| 94 |
const PFCandidate *pCand = iTau->IsoPFCand(i0);
|
| 95 |
float deta = iTau->Eta() - pCand->Eta();
|
| 96 |
float dphi = MathUtils::DeltaPhi((double)iTau ->Phi(), (double)pCand->Phi());
|
| 97 |
float dr = fabs(MathUtils::DeltaR ((double)iTau ->Phi(), (double)iTau ->Eta(),
|
| 98 |
(double)pCand->Phi(), (double)pCand->Eta()));
|
| 99 |
int pftype = 0;
|
| 100 |
|
| 101 |
|
| 102 |
if(pCand->Charge() != 0) pftype = 0;
|
| 103 |
else if(pCand->PFType() == PFCandidate::eGamma) pftype = 1;
|
| 104 |
else pftype = 2;
|
| 105 |
|
| 106 |
// Number of isolation candidates by type
|
| 107 |
|
| 108 |
niso[pftype]++;
|
| 109 |
|
| 110 |
// Isolation Rings
|
| 111 |
|
| 112 |
if(dr < 0.1) rings[pftype][0] += pCand->Pt();
|
| 113 |
else if(dr < 0.2) rings[pftype][1] += pCand->Pt();
|
| 114 |
else if(dr < 0.3) rings[pftype][2] += pCand->Pt();
|
| 115 |
else if(dr < 0.4) rings[pftype][3] += pCand->Pt();
|
| 116 |
else if(dr < 0.5) rings[pftype][4] += pCand->Pt();
|
| 117 |
|
| 118 |
// Angle Shape Variables
|
| 119 |
|
| 120 |
shapes[pftype][0] += pCand->Pt() * deta;
|
| 121 |
shapes[pftype][1] += pCand->Pt() * dphi;
|
| 122 |
shapes[pftype][2] += pCand->Pt() * deta*deta;
|
| 123 |
shapes[pftype][3] += pCand->Pt() * dphi*dphi;
|
| 124 |
shapes[pftype][4] += pCand->Pt() * deta*dphi;
|
| 125 |
isoptsum[pftype] += pCand->Pt();
|
| 126 |
}
|
| 127 |
|
| 128 |
// Mean and variance of angle variables are weighted by pT
|
| 129 |
|
| 130 |
for(uint i = 0; i < shapes.size(); i++)
|
| 131 |
{
|
| 132 |
for(uint j = 0; j < shapes[i].size(); j++)
|
| 133 |
{
|
| 134 |
shapes[i][j] = isoptsum[i] > 0 ? fabs(shapes[i][j]/isoptsum[i]) : 0;
|
| 135 |
}
|
| 136 |
}
|
| 137 |
|
| 138 |
// Fill IsoRing object
|
| 139 |
|
| 140 |
IsoRings isoRings;
|
| 141 |
isoRings.niso = niso;
|
| 142 |
isoRings.rings = rings;
|
| 143 |
isoRings.shapes = shapes;
|
| 144 |
|
| 145 |
return isoRings;
|
| 146 |
}
|
