ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Utils/src/TauIsoMVA.cc
Revision: 1.3
Committed: Thu Mar 28 13:44:17 2013 UTC (12 years, 1 month ago) by arapyan
Content type: text/plain
Branch: MAIN
CVS Tags: Mit_029c, Mit_029b, Mit_029a, HEAD
Changes since 1.2: +5 -1 lines
Log Message:
proper gbrforest initialization

File Contents

# User Rev Content
1 pharris 1.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 arapyan 1.3 #include "CondFormats/EgammaObjects/interface/GBRForest.h"
11     #include "Cintex/Cintex.h"
12     #include <utility>
13 pharris 1.1
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 pharris 1.2 fGBR = false;
71     }
72 arapyan 1.3 void TauIsoMVA::InitializeGBR( TString iWeightFile) {
73     ROOT::Cintex::Cintex::Enable();
74 pharris 1.2 TFile *lForest = new TFile(iWeightFile,"READ");
75     fGBRReader = (GBRForest*)lForest->Get("gbrfTauIso");
76     lForest->Close();
77     fGBR = true;
78 pharris 1.1 }
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 pharris 1.2 if(!fGBR) mvainput.insert(mvainput.end(), 6, 0);
84     if(!fGBR) return fReader->EvaluateMVA(mvainput, "BDTG");
85     return fGBRReader->GetClassifier(&mvainput[0]);
86 pharris 1.1 }
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     }