LArSoft  v06_85_00
Liquid Argon Software toolkit - http://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(), and FormatHistLong().

11 {
12 
13  const int nParticles = 5;
14 
15  TString testFileName = "/data/t2k/phsxvg/DUNEPid/data/4APA/EventSamples/train_muon_realflux_split_2.root";
16 
17  TString trainFileName[nParticles] = {"/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] = {"/data/t2k/phsxvg/DUNEPid/data/4APA/weights/muon_all_BDT.weights.xml",
24  "/data/t2k/phsxvg/DUNEPid/data/4APA/weights/electron_all_BDT.weights.xml",
25  "/data/t2k/phsxvg/DUNEPid/data/4APA/weights/proton_all_BDT.weights.xml",
26  "/data/t2k/phsxvg/DUNEPid/data/4APA/weights/pich_all_BDT.weights.xml",
27  "/data/t2k/phsxvg/DUNEPid/data/4APA/weights/photon_all_BDT.weights.xml"};
28 
29  TFile *fileForBinning = TFile::Open(trainFileName[0]);
30  fileForBinning->cd("Method_BDT/BDT");
31 
32  int numbins = MVA_BDT_Train_S->GetNbinsX();
33  const int nbins = numbins;
34 
35  fileForBinning->Close();
36 
37  double likeRatioS[nParticles][nbins] = {0};
38 
39  double binWidth[nParticles], MVAMin[nParticles], MVAMax[nParticles];
40 
41  for(int p=0; p<nParticles; p++)
42  {
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  {
52  cerr<<"MVA_BDT_Train_S->GetNbinsX() not equal to nbins in file "<<trainFileName[p]<<endl;
53  exit(1);
54  }
55  if(MVA_BDT_Train_B->GetNbinsX() != nbins)
56  {
57  cerr<<"MVA_BDT_Train_B->GetNbinsX() not equal to nbins in file "<<trainFileName[p]<<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  {
67  cerr<<"MVA_BDT_Train_B->GetBinWidth(1) not equal to binWidth in file "<<trainFileName[p]<<endl;
68  exit(1);
69  }
70  if(MVA_BDT_Train_B->GetBinLowEdge(1) != MVAMin[p])
71  {
72  cerr<<"MVA_BDT_Train_B->GetBinLowEdge(1) not equal to MVAMin in file "<<trainFileName[p]<<endl;
73  exit(1);
74  }
75  if(MVA_BDT_Train_B->GetBinLowEdge(nbins) + MVA_BDT_Train_B->GetBinWidth(nbins) != MVAMax[p])
76  {
77  cerr<<"MVA_BDT_Train_B->GetBinLowEdge(nbins) + MVA_BDT_Train_B->GetBinWidth(nbins) not equal to MVAMax in file "<<trainFileName[p]<<endl;
78  exit(1);
79  }
80 
81  for(int i=1; i<=nbins; i++)
82  {
83  likeRatioS[p][i] = MVA_BDT_Train_S->GetBinContent(i) / TMath::Max(0.01, MVA_BDT_Train_B->GetBinContent(i));
84  }
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, dEdxEndRatioR;
106 
107  vector<double> likeRatio[nParticles];
108  vector<double>::iterator likeIter;
109 
110  TH1D *hMVAValue[nParticles] = {0};
111  for(int p=0; p<nParticles; p++)
112  hMVAValue[p] = new TH1D(Form("hMVAValue%d", p), "", nbins, -0.5, 0.5);
113 
114  for(int p=0; p<nParticles; p++)
115  {
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  {
133  PIDVar->GetEntry(i);
134 
135  evalRatioR = evalRatio;
136  coreHaloRatioR = coreHaloRatio;
137  concentrationR = concentration;
138  conicalnessR = conicalness;
139  dEdxStartR = dEdxStart;
140  dEdxEndR = dEdxEnd;
141  dEdxEndRatioR = dEdxEndRatio;
142 
143  MVAValue = reader->EvaluateMVA("BDT");
144 
145  MVABin = 1 + int((MVAValue - MVAMin[p]) / binWidth[p]);
146 
147  hMVAValue[p]->Fill(MVAValue);
148 
149  likeRatio[p].push_back(likeRatioS[p][MVABin]);
150 
151  }
152 
153  delete reader;
154 
155  }
156 
157  TCanvas *MVAVal = new TCanvas("MVAVal", "MVAVal", 800, 600);
158  FormatCanvas(MVAVal);
159  hMVAValue[0]->GetXaxis()->SetTitle("BDT variable");
160  hMVAValue[0]->GetYaxis()->SetRangeUser(0, 250);
161  FormatHistLong(hMVAValue[0], 2, 1);
162  hMVAValue[0]->Draw();
163  FormatHist(hMVAValue[1], 4, 1);
164  hMVAValue[1]->Draw("same");
165  FormatHist(hMVAValue[2], 6, 1);
166  hMVAValue[2]->Draw("same");
167  FormatHist(hMVAValue[3], 8, 1);
168  hMVAValue[3]->Draw("same");
169  FormatHist(hMVAValue[4], 1, 1);
170  hMVAValue[4]->Draw("same");
171 
172  TLegend *lM = new TLegend(0.54, 0.57, 0.84, 0.87);
173  lM->AddEntry(hMVAValue[0], "Muon weights", "L");
174  lM->AddEntry(hMVAValue[1], "Electron weights", "L");
175  lM->AddEntry(hMVAValue[2], "Proton weights", "L");
176  lM->AddEntry(hMVAValue[3], "#pi^{#pm} weights", "L");
177  lM->AddEntry(hMVAValue[4], "Photon weights", "L");
178  lM->SetBorderSize(0);
179  lM->Draw();
180 
181  int likeVectorSize = likeRatio[0].size();
182  for(int p=1; p<nParticles; p++)
183  {
184  if(likeRatio[p].size() != likeVectorSize)
185  {
186  cerr<<"Likelihood ratio vector sizes not equal"<<endl;
187  exit(1);
188  }
189  }
190 
191  int vectorIndex = 0;
192 
193  double particleLikeRatio[nParticles] = {0};
194  double totalLikeRatio = 0.0;
195 
196  TH1D *hParticlePLR[nParticles] = {0};
197 
198  for(int p=0; p<nParticles; p++)
199  hParticlePLR[p] = new TH1D(Form("hParticlePLR%d", p), "", nProbBins, 0.0, 1.0);
200 
201  TH1D *hLikeRatio[nParticles] = {0};
202 
203  for(int p=0; p<nParticles; p++)
204  hLikeRatio[p] = new TH1D(Form("hLikeRatio%d", p), "", 40, -10.0, 10.0);
205 
206  while(vectorIndex < likeRatio[0].size())
207  {
208  totalLikeRatio = 0.0;
209 
210  for(int p=0; p<nParticles; p++)
211  {
212  likeIter = likeRatio[p].begin() + vectorIndex;
213  particleLikeRatio[p] = (*likeIter);
214  totalLikeRatio += (*likeIter);
215  }
216 
217  for(int p=0; p<nParticles; p++)
218  {
219  //if particleLikeRatio[p] / totalLikeRatio = 1.0, this goes to overflow bin of hParticlePLR[p]
220  //for this reason, set it to be no more than 0.9999
221  hParticlePLR[p]->Fill(TMath::Min(0.9999, particleLikeRatio[p] / totalLikeRatio));
222  }
223 
224  ++vectorIndex;
225  }
226 
227  TCanvas *ParticlePLR = new TCanvas("ParticlePLR", "ParticlePLR", 800, 600);
228  FormatCanvas(ParticlePLR);
229  hParticlePLR[0]->GetXaxis()->SetTitle("Proportion of likelihood ratio");
230  hParticlePLR[0]->GetYaxis()->SetRangeUser(0, 800);
231  FormatHistLong(hParticlePLR[0], 2, 1);
232  hParticlePLR[0]->Draw();
233  FormatHist(hParticlePLR[1], 4, 1);
234  hParticlePLR[1]->Draw("same");
235  FormatHist(hParticlePLR[2], 6, 1);
236  hParticlePLR[2]->Draw("same");
237  FormatHist(hParticlePLR[3], 8, 1);
238  hParticlePLR[3]->Draw("same");
239  FormatHist(hParticlePLR[4], 1, 1);
240  hParticlePLR[4]->Draw("same");
241 
242  TLegend *lP = new TLegend(0.6, 0.55, 0.83, 0.85);
243  lP->AddEntry(hParticlePLR[0], "Muon", "L");
244  lP->AddEntry(hParticlePLR[1], "Electron", "L");
245  lP->AddEntry(hParticlePLR[2], "Proton", "L");
246  lP->AddEntry(hParticlePLR[3], "#pi^{#pm}", "L");
247  lP->AddEntry(hParticlePLR[4], "Photon", "L");
248  lP->SetBorderSize(0);
249  lP->Draw();
250 
251 }//end of CalcPropLikeRatio()
iterator
intermediate_table::iterator iterator
Definition: intermediate_table.cc:24
FormatCanvas
void FormatCanvas(TCanvas *canvas)
Definition: CalcPropLikeRatio.C:253
FormatHistLong
void FormatHistLong(TH1D *hist, int colour, int style)
Definition: CalcPropLikeRatio.C:262
FormatHist
void FormatHist(TH1D *hist, int colour, int style)
Definition: CalcPropLikeRatio.C:275
void FormatCanvas ( TCanvas *  canvas)

Definition at line 253 of file CalcPropLikeRatio.C.

Referenced by CalcPropLikeRatio().

254 {
255  canvas->Divide(1);
256  canvas->cd(1);
257  gPad->SetLeftMargin(0.18);
258  gPad->SetBottomMargin(0.15);
259  gPad->SetRightMargin(0.15);
260 }
void FormatHist ( TH1D *  hist,
int  colour,
int  style 
)

Definition at line 275 of file CalcPropLikeRatio.C.

Referenced by CalcPropLikeRatio().

276 {
277  hist->SetTitle("");
278  hist->SetLineColor(colour);
279  hist->SetLineWidth(2);
280  hist->SetLineStyle(style);
281  hist->SetStats(0);
282 }
hist
TH2F * hist
Definition: plot.C:136
void FormatHistLong ( TH1D *  hist,
int  colour,
int  style 
)

Definition at line 262 of file CalcPropLikeRatio.C.

Referenced by CalcPropLikeRatio().

263 {
264  hist->SetTitle("");
265  hist->SetLineColor(colour);
266  hist->SetLineWidth(2);
267  hist->SetLineStyle(style);
268  hist->SetStats(0);
269  hist->GetXaxis()->SetTitleOffset(1.2);
270  hist->GetYaxis()->SetTitleOffset(1.4);
271  hist->GetXaxis()->CenterTitle();
272  hist->GetYaxis()->CenterTitle();
273 }
hist
TH2F * hist
Definition: plot.C:136