5 #include "CAFAna/Cuts/QuantileCuts.h" 6 #include "CAFAna/Systs/SystConcat.h" 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);
52 for (
int i = 0;
i < h->GetNbinsX(); ++
i)
61 TH1D*
proj = h->ProjectionX(histName.c_str(), cutBin+1, h->GetNbinsY()+1);
69 TH1D*
proj = h->ProjectionX(histName.c_str(), cutBin+1, h->GetNbinsY()+1, bdtCutBin+1, h->GetNbinsZ()+1);
78 int ndBins = nd->GetNbinsX();
79 int fdBins = fd->GetNbinsX();
80 int totBins = ndBins + fdBins;
83 TH1D*
h =
new TH1D(name.c_str(),
";Energy Bin;", totBins, 0, totBins);
85 for (
int i = 0;
i < totBins; ++
i){
87 h->SetBinContent(
i+1, nd->GetBinContent(
i+1));
89 h->SetBinContent(
i+1, fd->GetBinContent(
i+1-ndBins));
106 int nBins = spectra->GetNbinsX();
113 fullCov(
i,
j) = inMat(
i,
j) *
114 spectra->GetBinContent(
i+1) *
115 spectra->GetBinContent(
j+1);
119 fullCov(
i,
j) = fullCov(
i,
j) + spectra->GetBinContent(
i+1);
135 TFile* covmxFile =
new TFile(
"/pnfs/nova/persistent/analysis/nux/nus20/selection/covmx.root");
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";
145 TFile* spectraFile =
new TFile(file.c_str(),
"read");
147 TH3D* hCVNBDTEnergyTotalMCFD
148 = (TH3D*)spectraFile->Get(fdTotalMC.c_str());
150 TH3D* hCVNBDTEnergySignalMCFD
151 = (TH3D*)spectraFile->Get(fdSignalMC.c_str());
153 TH3D* hCVNBDTEnergyCosmicFD
154 = (TH3D*)spectraFile->Get(fdCosmic.c_str());
157 = (TH2D*)spectraFile->Get(ndTotalMC.c_str());
159 TH2D* hCVNEnergySignalND
160 = (TH2D*)spectraFile->Get(ndSignalMC.c_str());
166 std::vector<TH2D*> fomValHists;
167 for (
int i = 0;
i < hCVNBDTEnergyTotalMCFD->GetNbinsZ(); ++
i){
169 fomValHists.push_back(
new TH2D(nameString.c_str(),
";ND CVN cut value;FD CVN cut value", 100, 0, 1, 100, 0, 1));
171 TH1D* ndMCE; TH1D* fdMCE; TH1D* fdTotalE; TH1D* allMCE; TH1D* allCosE; TH1D* allTotE;
172 TH1D* ndSigE; TH1D* fdSigE; TH1D* allSigE;
175 for (
int iND = 0; iND < 100; ++iND){
181 for (
int iBDT = 0; iBDT < hCVNBDTEnergyTotalMCFD->GetNbinsZ(); ++iBDT){
183 for (
int iFD = 0; iFD < 100; ++iFD){
188 fdTotalE->Add(fdCosE);
217 std::vector<double> fulldecorrelatedunc =
241 fomValHists.at(iBDT)->SetBinContent(iND+1,iFD+1, fomVal);
253 TFile*
outfile =
new TFile(
"outfilefom.root",
"recreate");
257 float selBinValI = 0;
258 float selBinValJ = 0;
259 float SelBDTBinVal = 0;
263 for(
int iBDT = 0; iBDT < fomValHists.size(); ++iBDT){
265 TH2D* fomVals = fomValHists.at(iBDT);
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");
274 c1->SaveAs((name+
".png").c_str());
275 c1->SaveAs((name+
".pdf").c_str());
277 int bdtBin = iBDT+(
int)(hCVNBDTEnergyTotalMCFD->GetZaxis()->GetBinLowEdge(1)*100);
279 for (
int i = 0;
i < fomVals->GetNbinsX(); ++
i){
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);
288 selBinValJ = fomVals->GetYaxis()->GetBinLowEdge(
j+1);
296 <<
"FOM is maximised by" << val
297 <<
" at ND: " << selBinValI <<
" (" << selBinI <<
") " 298 <<
", FD: " << selBinValJ <<
" (" << selBinJ <<
") " 299 <<
", BDT: " << SelBDTBinVal <<
" (" << selBinBDT <<
") " 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);
309 int bdtBin = iBDT+1+(
int)(hCVNBDTEnergyTotalMCFD->GetZaxis()->GetBinLowEdge(1)*100);
311 bdtProj->SetBinContent(bdtBin, fomVals->GetBinContent(selBinI, selBinJ));
313 std::cout <<
"set bin " << bdtBin <<
" to " << fomVals->GetBinContent(selBinI, selBinJ) <<
std::endl;
315 if (bdtBin == selBinBDT){
316 for (
int i = 0;
i < fomVals->GetNbinsX()+1; ++
i){
317 ndProj->SetBinContent(
i, fomVals->GetBinContent(
i, selBinJ));
320 for (
int j = 0;
j < fomVals->GetNbinsY()+1; ++
j){
322 fdProj->SetBinContent(
j, fomVals->GetBinContent(selBinI,
j));
327 std::cout <<
"projections have a maximum of: (bdt: " << bdtProj->GetMaximum()
328 <<
", fd: " << fdProj->GetMaximum()
329 <<
", nd: " << ndProj->GetMaximum() <<
")" 332 TCanvas *
c1 =
new TCanvas(
"c1",
"c1", 600, 400);
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");
TH1D * GetSelectedEvents(TH2D *h, std::string det, int cutBin, Types typ)
void AddStatisticalUncertainty(std::vector< double > *totunc, TH1D *allTotE)
double calcBinnedSOverSSB(TH1D *projSig, std::vector< double > totUnc)
TH1D * ConcatenateHists(TH1D *nd, TH1D *fd, Types typ)
TMatrixD FracToFull(TMatrixD *fracCov, TH1D *spectra)
std::string to_string(ModuleType const mt)
std::vector< double > GetDecorrelatedUncertainty(TMatrixD *inmat, int ndbins, int fdbins)
void NDFDCVNBDTCutSelector(std::string file)