LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
Functions
CalcPropLikeRatio.C File Reference

Go to the source code of this file.

Functions

void FormatCanvas (TCanvas *canvas)
 
void FormatHistLong (TH1D *hist, int colour, int style)
 
void FormatHist (TH1D *hist, int colour, int style)
 
void CalcPropLikeRatio (const int nProbBins=25)
 

Function Documentation

void CalcPropLikeRatio ( const int  nProbBins = 25)

Definition at line 8 of file CalcPropLikeRatio.C.

References FormatCanvas(), FormatHist(), FormatHistLong(), and util::size().

9 {
10 
11  const int nParticles = 5;
12 
13  TString testFileName =
14  "/data/t2k/phsxvg/DUNEPid/data/4APA/EventSamples/train_muon_realflux_split_2.root";
15 
16  TString trainFileName[nParticles] = {
17  "/data/t2k/phsxvg/DUNEPid/data/4APA/TrainMVA/mvaPlots_muon_all.root",
18  "/data/t2k/phsxvg/DUNEPid/data/4APA/TrainMVA/mvaPlots_electron_all.root",
19  "/data/t2k/phsxvg/DUNEPid/data/4APA/TrainMVA/mvaPlots_proton_all.root",
20  "/data/t2k/phsxvg/DUNEPid/data/4APA/TrainMVA/mvaPlots_pich_all.root",
21  "/data/t2k/phsxvg/DUNEPid/data/4APA/TrainMVA/mvaPlots_photon_all.root"};
22 
23  TString weightFileName[nParticles] = {
24  "/data/t2k/phsxvg/DUNEPid/data/4APA/weights/muon_all_BDT.weights.xml",
25  "/data/t2k/phsxvg/DUNEPid/data/4APA/weights/electron_all_BDT.weights.xml",
26  "/data/t2k/phsxvg/DUNEPid/data/4APA/weights/proton_all_BDT.weights.xml",
27  "/data/t2k/phsxvg/DUNEPid/data/4APA/weights/pich_all_BDT.weights.xml",
28  "/data/t2k/phsxvg/DUNEPid/data/4APA/weights/photon_all_BDT.weights.xml"};
29 
30  TFile* fileForBinning = TFile::Open(trainFileName[0]);
31  fileForBinning->cd("Method_BDT/BDT");
32 
33  int numbins = MVA_BDT_Train_S->GetNbinsX();
34  const int nbins = numbins;
35 
36  fileForBinning->Close();
37 
38  double likeRatioS[nParticles][nbins] = {0};
39 
40  double binWidth[nParticles], MVAMin[nParticles], MVAMax[nParticles];
41 
42  for (int p = 0; p < nParticles; p++) {
43  cout << endl;
44  cout << "Calculating likelihood ratios from file " << trainFileName[p] << endl;
45  cout << endl;
46 
47  TFile* trainFile = TFile::Open(trainFileName[p]);
48  trainFile->cd("Method_BDT/BDT");
49 
50  if (MVA_BDT_Train_S->GetNbinsX() != nbins) {
51  cerr << "MVA_BDT_Train_S->GetNbinsX() not equal to nbins in file " << trainFileName[p]
52  << endl;
53  exit(1);
54  }
55  if (MVA_BDT_Train_B->GetNbinsX() != nbins) {
56  cerr << "MVA_BDT_Train_B->GetNbinsX() not equal to nbins in file " << trainFileName[p]
57  << endl;
58  exit(1);
59  }
60 
61  binWidth[p] = MVA_BDT_Train_S->GetBinWidth(1);
62  MVAMin[p] = MVA_BDT_Train_S->GetBinLowEdge(1);
63  MVAMax[p] = MVA_BDT_Train_S->GetBinLowEdge(nbins) + MVA_BDT_Train_S->GetBinWidth(nbins);
64 
65  if (MVA_BDT_Train_B->GetBinWidth(1) != binWidth[p]) {
66  cerr << "MVA_BDT_Train_B->GetBinWidth(1) not equal to binWidth in file " << trainFileName[p]
67  << endl;
68  exit(1);
69  }
70  if (MVA_BDT_Train_B->GetBinLowEdge(1) != MVAMin[p]) {
71  cerr << "MVA_BDT_Train_B->GetBinLowEdge(1) not equal to MVAMin in file " << trainFileName[p]
72  << endl;
73  exit(1);
74  }
75  if (MVA_BDT_Train_B->GetBinLowEdge(nbins) + MVA_BDT_Train_B->GetBinWidth(nbins) != MVAMax[p]) {
76  cerr << "MVA_BDT_Train_B->GetBinLowEdge(nbins) + MVA_BDT_Train_B->GetBinWidth(nbins) not "
77  "equal to MVAMax in file "
78  << trainFileName[p] << endl;
79  exit(1);
80  }
81 
82  for (int i = 1; i <= nbins; i++) {
83  likeRatioS[p][i] =
84  MVA_BDT_Train_S->GetBinContent(i) / TMath::Max(0.01, MVA_BDT_Train_B->GetBinContent(i));
85  }
86  }
87 
88  TFile* testfile = TFile::Open(testFileName);
89  testfile->cd();
90  TTree* PIDVar = (TTree*)testfile->Get("mvaTree");
91 
92  float evalRatio, coreHaloRatio, concentration, conicalness, dEdxStart, dEdxEnd, dEdxEndRatio;
93 
94  PIDVar->SetBranchAddress("evalRatio", &evalRatio);
95  PIDVar->SetBranchAddress("coreHaloRatio", &coreHaloRatio);
96  PIDVar->SetBranchAddress("concentration", &concentration);
97  PIDVar->SetBranchAddress("conicalness", &conicalness);
98  PIDVar->SetBranchAddress("dEdxStart", &dEdxStart);
99  PIDVar->SetBranchAddress("dEdxEnd", &dEdxEnd);
100  PIDVar->SetBranchAddress("dEdxEndRatio", &dEdxEndRatio);
101 
102  double MVAValue;
103  int MVABin;
104 
105  float evalRatioR, coreHaloRatioR, concentrationR, conicalnessR, dEdxStartR, dEdxEndR,
106  dEdxEndRatioR;
107 
108  vector<double> likeRatio[nParticles];
109  vector<double>::iterator likeIter;
110 
111  TH1D* hMVAValue[nParticles] = {0};
112  for (int p = 0; p < nParticles; p++)
113  hMVAValue[p] = new TH1D(Form("hMVAValue%d", p), "", nbins, -0.5, 0.5);
114 
115  for (int p = 0; p < nParticles; p++) {
116  cout << endl;
117  cout << "Reading weights from file " << weightFileName[p] << endl;
118 
119  TMVA::Reader* reader = new TMVA::Reader();
120 
121  reader->AddVariable("evalRatio", &evalRatioR);
122  reader->AddVariable("coreHaloRatio", &coreHaloRatioR);
123  reader->AddVariable("concentration", &concentrationR);
124  reader->AddVariable("conicalness", &conicalnessR);
125  reader->AddVariable("dEdxStart", &dEdxStartR);
126  reader->AddVariable("dEdxEnd", &dEdxEndR);
127  reader->AddVariable("dEdxEndRatio", &dEdxEndRatioR);
128 
129  reader->BookMVA("BDT", weightFileName[p]);
130 
131  for (int i = 0; i < PIDVar->GetEntries(); ++i) {
132  PIDVar->GetEntry(i);
133 
134  evalRatioR = evalRatio;
135  coreHaloRatioR = coreHaloRatio;
136  concentrationR = concentration;
137  conicalnessR = conicalness;
138  dEdxStartR = dEdxStart;
139  dEdxEndR = dEdxEnd;
140  dEdxEndRatioR = dEdxEndRatio;
141 
142  MVAValue = reader->EvaluateMVA("BDT");
143 
144  MVABin = 1 + int((MVAValue - MVAMin[p]) / binWidth[p]);
145 
146  hMVAValue[p]->Fill(MVAValue);
147 
148  likeRatio[p].push_back(likeRatioS[p][MVABin]);
149  }
150 
151  delete reader;
152  }
153 
154  TCanvas* MVAVal = new TCanvas("MVAVal", "MVAVal", 800, 600);
155  FormatCanvas(MVAVal);
156  hMVAValue[0]->GetXaxis()->SetTitle("BDT variable");
157  hMVAValue[0]->GetYaxis()->SetRangeUser(0, 250);
158  FormatHistLong(hMVAValue[0], 2, 1);
159  hMVAValue[0]->Draw();
160  FormatHist(hMVAValue[1], 4, 1);
161  hMVAValue[1]->Draw("same");
162  FormatHist(hMVAValue[2], 6, 1);
163  hMVAValue[2]->Draw("same");
164  FormatHist(hMVAValue[3], 8, 1);
165  hMVAValue[3]->Draw("same");
166  FormatHist(hMVAValue[4], 1, 1);
167  hMVAValue[4]->Draw("same");
168 
169  TLegend* lM = new TLegend(0.54, 0.57, 0.84, 0.87);
170  lM->AddEntry(hMVAValue[0], "Muon weights", "L");
171  lM->AddEntry(hMVAValue[1], "Electron weights", "L");
172  lM->AddEntry(hMVAValue[2], "Proton weights", "L");
173  lM->AddEntry(hMVAValue[3], "#pi^{#pm} weights", "L");
174  lM->AddEntry(hMVAValue[4], "Photon weights", "L");
175  lM->SetBorderSize(0);
176  lM->Draw();
177 
178  int likeVectorSize = likeRatio[0].size();
179  for (int p = 1; p < nParticles; p++) {
180  if (likeRatio[p].size() != likeVectorSize) {
181  cerr << "Likelihood ratio vector sizes not equal" << endl;
182  exit(1);
183  }
184  }
185 
186  int vectorIndex = 0;
187 
188  double particleLikeRatio[nParticles] = {0};
189  double totalLikeRatio = 0.0;
190 
191  TH1D* hParticlePLR[nParticles] = {0};
192 
193  for (int p = 0; p < nParticles; p++)
194  hParticlePLR[p] = new TH1D(Form("hParticlePLR%d", p), "", nProbBins, 0.0, 1.0);
195 
196  TH1D* hLikeRatio[nParticles] = {0};
197 
198  for (int p = 0; p < nParticles; p++)
199  hLikeRatio[p] = new TH1D(Form("hLikeRatio%d", p), "", 40, -10.0, 10.0);
200 
201  while (vectorIndex < likeRatio[0].size()) {
202  totalLikeRatio = 0.0;
203 
204  for (int p = 0; p < nParticles; p++) {
205  likeIter = likeRatio[p].begin() + vectorIndex;
206  particleLikeRatio[p] = (*likeIter);
207  totalLikeRatio += (*likeIter);
208  }
209 
210  for (int p = 0; p < nParticles; p++) {
211  //if particleLikeRatio[p] / totalLikeRatio = 1.0, this goes to overflow bin of hParticlePLR[p]
212  //for this reason, set it to be no more than 0.9999
213  hParticlePLR[p]->Fill(TMath::Min(0.9999, particleLikeRatio[p] / totalLikeRatio));
214  }
215 
216  ++vectorIndex;
217  }
218 
219  TCanvas* ParticlePLR = new TCanvas("ParticlePLR", "ParticlePLR", 800, 600);
220  FormatCanvas(ParticlePLR);
221  hParticlePLR[0]->GetXaxis()->SetTitle("Proportion of likelihood ratio");
222  hParticlePLR[0]->GetYaxis()->SetRangeUser(0, 800);
223  FormatHistLong(hParticlePLR[0], 2, 1);
224  hParticlePLR[0]->Draw();
225  FormatHist(hParticlePLR[1], 4, 1);
226  hParticlePLR[1]->Draw("same");
227  FormatHist(hParticlePLR[2], 6, 1);
228  hParticlePLR[2]->Draw("same");
229  FormatHist(hParticlePLR[3], 8, 1);
230  hParticlePLR[3]->Draw("same");
231  FormatHist(hParticlePLR[4], 1, 1);
232  hParticlePLR[4]->Draw("same");
233 
234  TLegend* lP = new TLegend(0.6, 0.55, 0.83, 0.85);
235  lP->AddEntry(hParticlePLR[0], "Muon", "L");
236  lP->AddEntry(hParticlePLR[1], "Electron", "L");
237  lP->AddEntry(hParticlePLR[2], "Proton", "L");
238  lP->AddEntry(hParticlePLR[3], "#pi^{#pm}", "L");
239  lP->AddEntry(hParticlePLR[4], "Photon", "L");
240  lP->SetBorderSize(0);
241  lP->Draw();
242 
243 } //end of CalcPropLikeRatio()
iterator
intermediate_table::iterator iterator
Definition: intermediate_table.cc:27
FormatCanvas
void FormatCanvas(TCanvas *canvas)
Definition: CalcPropLikeRatio.C:245
util::size
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:101
FormatHistLong
void FormatHistLong(TH1D *hist, int colour, int style)
Definition: CalcPropLikeRatio.C:254
FormatHist
void FormatHist(TH1D *hist, int colour, int style)
Definition: CalcPropLikeRatio.C:267
void FormatCanvas ( TCanvas *  canvas)

Definition at line 245 of file CalcPropLikeRatio.C.

Referenced by CalcPropLikeRatio().

246 {
247  canvas->Divide(1);
248  canvas->cd(1);
249  gPad->SetLeftMargin(0.18);
250  gPad->SetBottomMargin(0.15);
251  gPad->SetRightMargin(0.15);
252 }
void FormatHist ( TH1D *  hist,
int  colour,
int  style 
)

Definition at line 267 of file CalcPropLikeRatio.C.

Referenced by CalcPropLikeRatio().

268 {
269  hist->SetTitle("");
270  hist->SetLineColor(colour);
271  hist->SetLineWidth(2);
272  hist->SetLineStyle(style);
273  hist->SetStats(0);
274 }
hist
TH2F * hist
Definition: plot.C:134
void FormatHistLong ( TH1D *  hist,
int  colour,
int  style 
)

Definition at line 254 of file CalcPropLikeRatio.C.

Referenced by CalcPropLikeRatio().

255 {
256  hist->SetTitle("");
257  hist->SetLineColor(colour);
258  hist->SetLineWidth(2);
259  hist->SetLineStyle(style);
260  hist->SetStats(0);
261  hist->GetXaxis()->SetTitleOffset(1.2);
262  hist->GetYaxis()->SetTitleOffset(1.4);
263  hist->GetXaxis()->CenterTitle();
264  hist->GetYaxis()->CenterTitle();
265 }
hist
TH2F * hist
Definition: plot.C:134