NDFDCVNBDTCutSelector.C
Go to the documentation of this file.
1 #include "CAFAna/Core/LoadFromFile.h"
2 #include "CAFAna/Core/Sample.h"
3 #include "CAFAna/Core/Registry.h"
4 #include "NuXAna/Vars/FitVarsSterile.h"
5 #include "CAFAna/Cuts/QuantileCuts.h"
6 #include "CAFAna/Systs/SystConcat.h"
7 #include "CAFAna/Systs/XSecSystLists.h"
8 #include "CAFAna/Prediction/PredictionInterp.h"
9 #include "CAFAna/Prediction/CovarianceMatrix.h"
10 #include "CAFAna/Experiment/GaussianConstraint.h"
11 #include "CAFAna/Analysis/Exposures.h"
12 
13 #include "3FlavorAna/Cuts/NumuCuts2020.h"
14 #include "3FlavorAna/Cuts/NueCuts2020.h"
15 
16 #include "NuXAna/Cuts/NusCuts18.h"
17 #include "NuXAna/Cuts/NusCuts20.h"
18 #include "NuXAna/Prediction/NDPredictionSterile.h"
19 #include "NuXAna/Prediction/FDPredictionSterile.h"
20 #include "NuXAna/Prediction/PredictionGeneratorNuX.h"
21 #include "NuXAna/Analysis/SystMaker.h"
22 #include "NuXAna/macros/Nus20/Selection/FOMUtilities.h"
23 #include "NuXAna/macros/Nus20/Selection/DecorrelateMatrix.h"
24 
25 #include "TH1.h"
26 #include "TH2.h"
27 #include "TH3.h"
28 #include "TMatrixT.h"
29 
30 enum Types {
31  Signal,
32  MCTotal,
33  Total
34 };
35 
36 // add in statistical uncertainty
37 // ---------------------------------------------------------------------------
38 void AddStatisticalUncertainty(std::vector<double>* totunc, TH1D* allTotE){
39  for (int i = 0; i < totunc->size(); ++i){
40  double systFrac = totunc->at(i)/allTotE->GetBinContent(i+1);
41  double statFrac = std::sqrt(allTotE->GetBinContent(i+1))/allTotE->GetBinContent(i+1);
42  totunc->at(i) = allTotE->GetBinContent(i+1)*std::sqrt(std::pow(systFrac,2) + std::pow(statFrac,2));
43  }
44 
45  return;
46 }
47 
48 // print bin contents of histogram
49 // ---------------------------------------------------------------------------
50 void PrintHist(TH1D* h){
51  std::cout << "hist name: " << h->GetName() << std::endl;
52  for (int i = 0; i < h->GetNbinsX(); ++i)
53  std::cout << h->GetBinContent(i+1) << std::endl;
54 }
55 
56 // apply projection to get selected events for a given cut value
57 //-----------------------------------------------------------------------------
58 TH1D* GetSelectedEvents(TH2D* h, std::string det, int cutBin, Types typ){
59 
60  std::string histName = det+"proj"+"_"+std::to_string(typ)+"_"+std::to_string(cutBin);
61  TH1D* proj = h->ProjectionX(histName.c_str(), cutBin+1, h->GetNbinsY()+1);
62 
63  return proj;
64 }
65 
66 TH1D* GetSelectedEvents(TH3D* h, std::string det, int cutBin, int bdtCutBin, Types typ){
67 
68  std::string histName = det+"proj"+"_"+std::to_string(typ)+"_"+std::to_string(cutBin)+"_"+std::to_string(bdtCutBin);
69  TH1D* proj = h->ProjectionX(histName.c_str(), cutBin+1, h->GetNbinsY()+1, bdtCutBin+1, h->GetNbinsZ()+1);
70 
71  return proj;
72 }
73 
74 // concatenate the nd and fd histograms
75 // ----------------------------------------------------------------------------
76 TH1D* ConcatenateHists(TH1D* nd, TH1D* fd, Types typ){
77 
78  int ndBins = nd->GetNbinsX();
79  int fdBins = fd->GetNbinsX();
80  int totBins = ndBins + fdBins;
81 
82  std::string name = "allMCE" + std::to_string(typ);
83  TH1D* h = new TH1D(name.c_str(), ";Energy Bin;", totBins, 0, totBins);
84 
85  for (int i = 0; i < totBins; ++i){
86  if (i < ndBins)
87  h->SetBinContent(i+1, nd->GetBinContent(i+1));
88  else
89  h->SetBinContent(i+1, fd->GetBinContent(i+1-ndBins));
90  }
91 
92  //PrintHist(nd);
93  //PrintHist(fd);
94  //PrintHist(h);
95 
96  return h;
97 }
98 
99 // covnert fractional covariance matrix to full covariance matrix
100 // including the statistical uncertainty
101 // ----------------------------------------------------------------------------
102 TMatrixD FracToFull(TMatrixD* fracCov, TH1D* spectra){
103 
104  TMatrixD inMat = *fracCov;
105 
106  int nBins = spectra->GetNbinsX();
107  TMatrixD fullCov(nBins, nBins);
108 
109  for (int i =0; i < nBins; ++i){
110  for (int j = 0; j < nBins; ++j){
111 
112  // multiply out systematic error
113  fullCov(i,j) = inMat(i,j) *
114  spectra->GetBinContent(i+1) *
115  spectra->GetBinContent(j+1);
116 
117  // add in variance from statistical uncertainties
118  if (i == j)
119  fullCov(i,j) = fullCov(i,j) + spectra->GetBinContent(i+1);
120 
121  }
122  }
123 
124  //fullCov.Print();
125 
126  return fullCov;
127 
128 }
129 
130 // main function
131 // ----------------------------------------------------------------------------
132 void NDFDCVNBDTCutSelector(std::string file){
133 
134  // load covariance matrix
135  TFile* covmxFile = new TFile("/pnfs/nova/persistent/analysis/nux/nus20/selection/covmx.root");
136 
137  TMatrixD* fraccovmx = (TMatrixD*)covmxFile->Get("fraccovmx");
138 
139  std::string fdTotalMC = "FD/kCaloECVNnc_looseptpNCCosRejBDTG20PreselQualityFiducialTotalMC";
140  std::string fdCosmic = "FD/kCaloECVNnc_looseptpNCCosRejBDTG20PreselQualityFiducialTotalCosmic";
141  std::string fdSignalMC = "FD/kCaloECVNnc_looseptpNCCosRejBDTG20PreselQualityFiducialSignalMC";
142  std::string ndTotalMC = "ND/kCaloEkCVNnc_looseptpPreselQualityFiducialTotalMC";
143  std::string ndSignalMC = "ND/kCaloEkCVNnc_looseptpPreselQualityFiducialSignalMC";
144 
145  TFile* spectraFile = new TFile(file.c_str(), "read");
146 
147  TH3D* hCVNBDTEnergyTotalMCFD
148  = (TH3D*)spectraFile->Get(fdTotalMC.c_str());
149 
150  TH3D* hCVNBDTEnergySignalMCFD
151  = (TH3D*)spectraFile->Get(fdSignalMC.c_str());
152 
153  TH3D* hCVNBDTEnergyCosmicFD
154  = (TH3D*)spectraFile->Get(fdCosmic.c_str());
155 
156  TH2D* hCVNEnergyND
157  = (TH2D*)spectraFile->Get(ndTotalMC.c_str());
158 
159  TH2D* hCVNEnergySignalND
160  = (TH2D*)spectraFile->Get(ndSignalMC.c_str());
161 
162 
163  // now need to scale to some POT
164 
165  // produce histograms for each possible BDT cut value
166  std::vector<TH2D*> fomValHists;
167  for (int i = 0; i < hCVNBDTEnergyTotalMCFD->GetNbinsZ(); ++i){
168  std::string nameString = "fomVal" + std::to_string(hCVNBDTEnergyTotalMCFD->GetZaxis()->GetBinLowEdge(i+1));
169  fomValHists.push_back(new TH2D(nameString.c_str(), ";ND CVN cut value;FD CVN cut value", 100, 0, 1, 100, 0, 1));
170  }
171  TH1D* ndMCE; TH1D* fdMCE; TH1D* fdTotalE; TH1D* allMCE; TH1D* allCosE; TH1D* allTotE;
172  TH1D* ndSigE; TH1D* fdSigE; TH1D* allSigE;
173  TH1D* fdCosE;
174 
175  for (int iND = 0; iND < 100; ++iND){
176  ndMCE = GetSelectedEvents(hCVNEnergyND , "nd" , iND, Types::MCTotal);
177  ndSigE = GetSelectedEvents(hCVNEnergySignalND, "nd" , iND, Types::Signal);
178 
179  std::cout << iND << std::endl;
180 
181  for (int iBDT = 0; iBDT < hCVNBDTEnergyTotalMCFD->GetNbinsZ(); ++iBDT){
182 
183  for (int iFD = 0; iFD < 100; ++iFD){
184  fdMCE = GetSelectedEvents(hCVNBDTEnergyTotalMCFD , "fd", iFD, iBDT, Types::MCTotal);
185  fdSigE = GetSelectedEvents(hCVNBDTEnergySignalMCFD, "fd", iFD, iBDT, Types::Signal);
186  fdCosE = GetSelectedEvents(hCVNBDTEnergyCosmicFD , "fd", iFD, iBDT, Types::Total);
187  fdTotalE = fdMCE;
188  fdTotalE->Add(fdCosE);
189 
190  // concatenate histograms and get the full matrix
191  allMCE = ConcatenateHists(ndMCE, fdMCE , Types::MCTotal);
192  allSigE = ConcatenateHists(ndSigE, fdSigE , Types::Signal);
193  allTotE = ConcatenateHists(ndMCE, fdTotalE, Types::Total);
194 
195  //TCanvas* c1 = new TCanvas("c1", "c1", 500, 500);
196  //c1->SetLogy(1);
197  //c1->SetLogz(0);
198  //allMCE->GetYaxis()->SetRangeUser(1, 10e6);
199  //allMCE->Draw();
200  //std::string allMCEName = "allMCE_ND_" + std::to_string(iND) +
201  // "_FD_" + std::to_string(iFD) +
202  // "_BDT_" + std::to_string(iBDT) + ".png";
203  //if (iND == 10 || iFD == 10)
204  // c1->SaveAs(allMCEName.c_str());
205 
206  TMatrixD fullcovmx = FracToFull(fraccovmx, allMCE);
207  //std::string covmxName = "covmx_ND_" + std::to_string(iND) +
208  // "_FD_" + std::to_string(iFD) +
209  // "_BDT_" + std::to_string(iBDT) + ".png";
210  //c1->SetLogy(0);
211  //c1->SetLogz(1);
212  //fullcovmx.Draw("colz");
213  //if (iND == 10 || iFD == 10)
214  // c1->SaveAs(covmxName.c_str());
215 
216  // get full systematic uncertainty
217  std::vector<double> fulldecorrelatedunc =
218  GetDecorrelatedUncertainty(&fullcovmx,
219  ndMCE->GetNbinsX(),
220  fdMCE->GetNbinsX());
221 
222  AddStatisticalUncertainty(&fulldecorrelatedunc, allTotE);
223 
224  //TH1D* unc = new TH1D("fulldecorr", ";Energy Bin;Decorrelated Uncertainty",
225  // fulldecorrelatedunc.size(), 0, fulldecorrelatedunc.size());
226 
227  //for (int ib = 0; ib < unc->GetNbinsX(); ++ib){
228  // unc->SetBinContent(ib+1, fulldecorrelatedunc.at(ib));
229  //}
230 
231  //c1->SetLogy(1);
232  //c1->SetLogz(0);
233  //unc->GetYaxis()->SetRangeUser(0.1, 10e5);
234  //unc->Draw();
235  //std::string uncName = "uncertainty_ND_" + std::to_string(iND) + "_FD_" + std::to_string(iFD) + ".png";
236  //if (iND == 10 || iFD == 10)
237  // c1->SaveAs(uncName.c_str());
238 
239  double fomVal = calcBinnedSOverSSB(allSigE, fulldecorrelatedunc);
240 
241  fomValHists.at(iBDT)->SetBinContent(iND+1,iFD+1, fomVal);
242 
243  delete allMCE;
244  delete allTotE;
245  delete allSigE;
246  //delete unc;
247  //delete c1;
248 
249  }
250  }
251  }
252 
253  TFile* outfile = new TFile("outfilefom.root", "recreate");
254  outfile->cd();
255 
256  double val = 0;
257  float selBinValI = 0;
258  float selBinValJ = 0;
259  float SelBDTBinVal = 0;
260  int selBinI = 0;
261  int selBinJ = 0;
262  int selBinBDT = 0;
263  for(int iBDT = 0; iBDT < fomValHists.size(); ++iBDT){
264 
265  TH2D* fomVals = fomValHists.at(iBDT);
266  std::cout << "fomVals(50) pre: " << fomVals->GetBinContent(50, 50) << std::endl;
267  fomVals->Write();
268 
269  TCanvas* c1 = new TCanvas("c1", "c1", 500, 500);
270  c1->SetRightMargin(0.18);
271  fomVals->GetZaxis()->SetRangeUser(1.0,5.5);
272  fomVals->Draw("colz");
273  std::string name = (std::string)fomVals->GetName();
274  c1->SaveAs((name+".png").c_str());
275  c1->SaveAs((name+".pdf").c_str());
276 
277  int bdtBin = iBDT+(int)(hCVNBDTEnergyTotalMCFD->GetZaxis()->GetBinLowEdge(1)*100);
278 
279  for (int i = 0; i < fomVals->GetNbinsX(); ++i){
280  std::cout << "i: " << i << std::endl;
281  for (int j = 0; j < fomVals->GetNbinsY(); ++j){
282  if (fomVals->GetBinContent(i+1,j+1) >= val){
283  val = fomVals->GetBinContent(i+1, j+1);
284  SelBDTBinVal = hCVNBDTEnergyTotalMCFD->GetZaxis()->GetBinLowEdge(iBDT+1);
285  selBinBDT = bdtBin+1;
286  selBinValI = fomVals->GetXaxis()->GetBinLowEdge(i+1);
287  selBinI = i+1;
288  selBinValJ = fomVals->GetYaxis()->GetBinLowEdge(j+1);
289  selBinJ = j+1;
290  }
291  }
292  }
293  }
294 
295  std::cout
296  << "FOM is maximised by" << val
297  << " at ND: " << selBinValI << " (" << selBinI << ") "
298  << ", FD: " << selBinValJ << " (" << selBinJ << ") "
299  << ", BDT: " << SelBDTBinVal << " (" << selBinBDT << ") "
300  << std::endl;
301 
302  // now make projected histograms
303  TH1D* bdtProj = new TH1D("bdtProj", ";BDT Score; FOM (projected)" , 100, 0, 1);
304  TH1D* ndProj = new TH1D("ndProj" , ";ND CVN Score; FOM (projected)", 100, 0, 1);
305  TH1D* fdProj = new TH1D("fdProj" , ";FD CVN Score; FOM (projected)", 100, 0, 1);
306  for(int iBDT = 0; iBDT < fomValHists.size(); ++iBDT){
307  TH2D* fomVals = fomValHists.at(iBDT);
308 
309  int bdtBin = iBDT+1+(int)(hCVNBDTEnergyTotalMCFD->GetZaxis()->GetBinLowEdge(1)*100);
310 
311  bdtProj->SetBinContent(bdtBin, fomVals->GetBinContent(selBinI, selBinJ));
312 
313  std::cout << "set bin " << bdtBin << " to " << fomVals->GetBinContent(selBinI, selBinJ) << std::endl;
314 
315  if (bdtBin == selBinBDT){
316  for (int i = 0; i < fomVals->GetNbinsX()+1; ++i){
317  ndProj->SetBinContent(i, fomVals->GetBinContent(i, selBinJ));
318  std::cout << "set ndProj " << i << " " << fomVals->GetBinContent(i, selBinJ) << std::endl;
319  }
320  for (int j = 0; j < fomVals->GetNbinsY()+1; ++j){
321  std::cout << "set fdProj " << j << " " << fomVals->GetBinContent(selBinI, j) << std::endl;
322  fdProj->SetBinContent(j, fomVals->GetBinContent(selBinI, j));
323  }
324  }
325  }
326 
327  std::cout << "projections have a maximum of: (bdt: " << bdtProj->GetMaximum()
328  << ", fd: " << fdProj->GetMaximum()
329  << ", nd: " << ndProj->GetMaximum() << ")"
330  << std::endl;
331 
332  TCanvas *c1 = new TCanvas("c1", "c1", 600, 400);
333  c1->cd();
334  bdtProj->SetLineWidth(2);
335  bdtProj->Draw("hist");
336  c1->SaveAs("bdtProj.png");
337  c1->SaveAs("bdtProj.pdf");
338  ndProj->SetLineWidth(2);
339  ndProj->Draw("hist");
340  c1->SaveAs("ndProj.png");
341  c1->SaveAs("ndProj.pdf");
342  fdProj->SetLineWidth(2);
343  fdProj->Draw("hist");
344  c1->SaveAs("fdProj.png");
345  c1->SaveAs("fdProj.pdf");
346 }
name
const XML_Char * name
Definition: expat.h:151
plotROC.nBins
int nBins
Definition: plotROC.py:16
nd
bool nd
Definition: attenuation_calibration_comparisons.C:31
GetSelectedEvents
TH1D * GetSelectedEvents(TH2D *h, std::string det, int cutBin, Types typ)
Definition: NDFDCVNBDTCutSelector.C:58
AddStatisticalUncertainty
void AddStatisticalUncertainty(std::vector< double > *totunc, TH1D *allTotE)
Definition: NDFDCVNBDTCutSelector.C:38
std::sqrt
T sqrt(T number)
Definition: d0nt_math.hpp:156
cet::pow
constexpr T pow(T x)
Definition: pow.h:75
fd
bool fd
Definition: attenuation_calibration_comparisons.C:32
calcBinnedSOverSSB
double calcBinnedSOverSSB(TH1D *projSig, std::vector< double > totUnc)
Definition: FOMUtilities.h:161
make_syst_table_plots.h
h
Definition: make_syst_table_plots.py:240
ConcatenateHists
TH1D * ConcatenateHists(TH1D *nd, TH1D *fd, Types typ)
Definition: NDFDCVNBDTCutSelector.C:76
calib::j
const double j
Definition: BetheBloch.cxx:29
om::cout
OStream cout
Definition: OStream.cxx:6
PrintHist
void PrintHist(TH1D *h)
Definition: NDFDCVNBDTCutSelector.C:50
proj
Float_t proj
Definition: plot.C:35
file
TFile * file
Definition: cellShifts.C:17
demo5.c1
c1
Definition: demo5.py:24
art::to_string
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
FracToFull
TMatrixD FracToFull(TMatrixD *fracCov, TH1D *spectra)
Definition: NDFDCVNBDTCutSelector.C:102
febshutoff_auto.val
val
Definition: febshutoff_auto.py:277
outfile
FILE * outfile
Definition: dump_event.C:13
GetDecorrelatedUncertainty
std::vector< double > GetDecorrelatedUncertainty(TMatrixD *inmat, int ndbins, int fdbins)
Definition: DecorrelateMatrix.h:116
NDFDCVNBDTCutSelector
void NDFDCVNBDTCutSelector(std::string file)
Definition: NDFDCVNBDTCutSelector.C:132
string
enum BeamMode string
Generated on Mon Nov 30 2020 05:15:01 for NOvA by   doxygen 1.8.11