19 #include "TTreeReader.h" 20 #include "TTreeReaderValue.h" 71 void PrintInfo(
int Entry,
int TotEnt,
double IsNuMuCC,
double TrueNuE,
double RecoNuE,
double EventWgt,
double CosVeto,
72 double CVNMuonLooPTP,
double ReMId,
double CVNProngMuon,
double CosRej2020 );
74 void CalcFOM( TH3D* Sig, TH3D*
Bkg, TH3D* Data,
int bintrpid,
int bincvn,
int bincos,
bool IsFD,
75 double &FOM,
double &TotSig,
double &
TotBkg );
79 bool PassCutPoint(
double PIDScore,
double PIDCut,
double CVNScore,
double CVNCut,
80 int Quant,
int CosVeto, std::vector<int> cQuantVec,
double RecoE,
bool isDip );
82 void FillMyPlots (
double CVN,
double ReMId,
double CosRej,
int Quant, TH1D* Hists[9],
double RecoE,
double Wgt );
86 double RecoEBins[
HistBins+1] = { 0, 0.75, 1, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2, 2.25, 2.5, 2.75, 3, 3.5, 4, 5 };
98 gStyle->SetOptStat(0);
102 if (!IsFD) { sIsFD =
"ND", sData =
"Data"; }
105 if (!IsFHC) sIsFHC =
"RHC";
108 if (!IsReMId) sTrPID =
"ProngCVN";
111 if (!IsLoosePTP){sCVN =
"OldPresel";
std::cout <<
"\nHave you reenabled OldPresel CVN?\n" <<
std::endl;
return; }
114 std::vector<int> cQuantVec;
115 if (QuantEnum == 0 ) { sQuant =
"AllQuant"; cQuantVec = {1,2,3,4}; }
116 else if (QuantEnum == 123) { sQuant =
"Quant123"; cQuantVec = {1,2,3 }; }
117 else if (QuantEnum == 4 ) { sQuant =
"Quant4" ; cQuantVec = { 4}; }
118 else {
std::cout <<
"\n\nQuantEnum must be 0, 123 or 4. You gave me " << QuantEnum <<
"\n\n" <<
std::endl;
return; }
120 std::string PIDOutString = sIsFD+
"-"+sIsFHC+
"-"+sTrPID+
"-"+sCVN+
"-"+sQuant;
123 std::string MCFile =
"/nova/ana/users/karlwarb/Analysis2020/PIDOpt/TrimCAF/200204/TrimCAF-2020-"+sIsFD+
"-"+sIsFHC+
"-"+sMC +
".root";
124 std::string DataFile =
"/nova/ana/users/karlwarb/Analysis2020/PIDOpt/TrimCAF/200204/TrimCAF-2020-"+sIsFD+
"-"+sIsFHC+
"-"+sData+
".root";
131 std::cout <<
"\n Is this FD? " << IsFD <<
"\t Is this FHC? " << IsFHC <<
"\t QuantEnum? " << QuantEnum
133 <<
"\n\t MCFile is " << MCFile
134 <<
"\n\t DataFile is " << DataFile
135 <<
"\n\t Am I using ReMId as my Tracking PID? " << IsReMId <<
" --> sTrPID is " << sTrPID
136 <<
"\n\t Am I using the Loose PTP trained CVN? " << IsLoosePTP <<
" --> sCVN is " << sCVN
137 <<
"\n Therefore PIDOutString is " << PIDOutString
142 gSystem -> MakeDirectory( OutDir.c_str() );
145 TFile *MyMCSamp =
new TFile( MCFile .c_str() );
146 TFile *MyDataSamp =
new TFile( DataFile.c_str() );
149 TTreeReader MC_Reader(
"TrimTree", MyMCSamp);
150 const int MC_TotalTree = MC_Reader.GetEntries(
true);
151 TTreeReaderValue<double> MC_Run (MC_Reader,
"Run" );
152 TTreeReaderValue<double> MC_Evt (MC_Reader,
"Evt" );
153 TTreeReaderValue<double> MC_PPFXWgt (MC_Reader,
"PPFXWgt" );
154 TTreeReaderValue<double> MC_XSecWgt (MC_Reader,
"XSecWgt" );
155 TTreeReaderValue<double> MC_TrueNuE (MC_Reader,
"TrueNuE" );
156 TTreeReaderValue<double> MC_IsNuMuCC (MC_Reader,
"IsNuMuCC" );
157 TTreeReaderValue<double> MC_IntMode (MC_Reader,
"IntMode" );
158 TTreeReaderValue<double> MC_DumbOsc (MC_Reader,
"DumbOsc" );
159 TTreeReaderValue<double> MC_BestFit2019 (MC_Reader,
"BestFit2019" );
160 TTreeReaderValue<double> MC_RecoNuE (MC_Reader,
"RecoNuE" );
164 TTreeReaderValue<double> MC_LooPTPQuant (MC_Reader,
"LooPTPQuant" );
165 TTreeReaderValue<double> MC_CosVeto (MC_Reader,
"CosVeto" );
166 TTreeReaderValue<double> MC_CVNMuonLooPTP(MC_Reader,
"CVNMuonLooPTP");
167 TTreeReaderValue<double> MC_CVNCosmLooPTP(MC_Reader,
"CVNCosmLooPTP");
168 TTreeReaderValue<double> MC_ReMId (MC_Reader,
"ReMId" );
169 TTreeReaderValue<double> MC_CVNProngMuon (MC_Reader,
"CVNProngMuon" );
170 TTreeReaderValue<double> MC_CosRej2020 (MC_Reader,
"CosRej2020" );
174 TTreeReader Data_Reader(
"TrimTree", MyDataSamp);
175 const int Data_TotalTree = Data_Reader.GetEntries(
true);
176 TTreeReaderValue<double> Data_Run (Data_Reader,
"Run" );
177 TTreeReaderValue<double> Data_Evt (Data_Reader,
"Evt" );
178 TTreeReaderValue<double> Data_PPFXWgt (Data_Reader,
"PPFXWgt" );
179 TTreeReaderValue<double> Data_XSecWgt (Data_Reader,
"XSecWgt" );
180 TTreeReaderValue<double> Data_TrueNuE (Data_Reader,
"TrueNuE" );
181 TTreeReaderValue<double> Data_IsNuMuCC (Data_Reader,
"IsNuMuCC" );
182 TTreeReaderValue<double> Data_IntMode (Data_Reader,
"IntMode" );
183 TTreeReaderValue<double> Data_DumbOsc (Data_Reader,
"DumbOsc" );
184 TTreeReaderValue<double> Data_BestFit2019 (Data_Reader,
"BestFit2019" );
185 TTreeReaderValue<double> Data_RecoNuE (Data_Reader,
"RecoNuE" );
189 TTreeReaderValue<double> Data_LooPTPQuant (Data_Reader,
"LooPTPQuant" );
190 TTreeReaderValue<double> Data_CosVeto (Data_Reader,
"CosVeto" );
191 TTreeReaderValue<double> Data_CVNMuonLooPTP(Data_Reader,
"CVNMuonLooPTP");
192 TTreeReaderValue<double> Data_CVNCosmLooPTP(Data_Reader,
"CVNCosmLooPTP");
193 TTreeReaderValue<double> Data_ReMId (Data_Reader,
"ReMId" );
194 TTreeReaderValue<double> Data_CVNProngMuon (Data_Reader,
"CVNProngMuon" );
195 TTreeReaderValue<double> Data_CosRej2020 (Data_Reader,
"CosRej2020" );
200 TH3D* hMC_Sig_FullE_Hist =
new TH3D(
"MCSig-FullE-Hist",
"Monte Carlo Signal; TrPID Score; CVN Loose PTP Score; CosRej",
202 TH3D* hMC_Bkg_FullE_Hist =
new TH3D(
"MCBkg-FullE-Hist",
"Monte Carlo Background; TrPID Score; CVN Loose PTP Score; CosRej",
204 TH3D* hData_FullE_Hist =
new TH3D(
"Data-FullE-Hist" ,
"Cosmic Events; TrPID Score; CVN Loose PTP Score; CosRej",
207 TH3D* hMC_Sig_Dip_Hist =
new TH3D(
"MCSig-Dip-Hist",
"Monte Carlo Signal; TrPID Score; CVN Loose PTP Score; CosRej",
209 TH3D* hMC_Bkg_Dip_Hist =
new TH3D(
"MCBkg-Dip-Hist",
"Monte Carlo Background; TrPID Score; CVN Loose PTP Score; CosRej",
211 TH3D* hData_Dip_Hist =
new TH3D(
"Data-Dip-Hist" ,
"Cosmic Events; TrPID Score; CVN Loose PTP Score; CosRej",
215 TH1D* hRecoNuE_Sig[9];
216 TH1D* hRecoNuE_Bkg[9];
217 TH1D* hRecoNuE_Cos[9];
218 for (
int hh=0;
hh<9; ++
hh) {
221 if (
hh==0) { app =
"FOMFull-AllQ"; Col = 1; }
222 else if (
hh==1) { app =
"FOMFull-Q123"; Col = 1; }
223 else if (
hh==2) { app =
"FOMFull-Q4" ; Col = 1; }
224 else if (
hh==3) { app =
"FOMDip-AllQ" ; Col = 2; }
225 else if (
hh==4) { app =
"FOMDip-Q123" ; Col = 2; }
226 else if (
hh==5) { app =
"FOMDip-Q4" ; Col = 2; }
227 else if (
hh==6) { app =
"BadCut-AllQ" ; Col = 4; }
228 else if (
hh==7) { app =
"BadCut-Q123" ; Col = 4; }
229 else if (
hh==8) { app =
"BadCut-Q4" ; Col = 4; }
230 std::string axis =
";Reconstructed Neutrino Enegry (GeV);Events / 0.1 GeV";
231 hRecoNuE_Sig[
hh] =
new TH1D( TString(
"hSig-")+TString(app), axis.c_str(),
HistBins,
RecoEBins );
232 hRecoNuE_Bkg[
hh] =
new TH1D( TString(
"hBkg-")+TString(app), axis.c_str(),
HistBins,
RecoEBins );
233 hRecoNuE_Cos[
hh] =
new TH1D( TString(
"hCos-")+TString(app), axis.c_str(),
HistBins,
RecoEBins );
240 std::map< std::string, MyFOMParams > FOMMap;
246 while (MC_Reader.Next()) {
248 int Entry = MC_Reader.GetCurrentEntry();
251 double EventWgt = *MC_PPFXWgt * *MC_XSecWgt * *MC_BestFit2019;
255 PrintInfo( Entry, MC_TotalTree, *MC_IsNuMuCC, *MC_TrueNuE, *MC_RecoNuE, EventWgt, *MC_CosVeto,
256 *MC_CVNMuonLooPTP, *MC_ReMId, *MC_CVNProngMuon, *MC_CosRej2020 );
260 if (*MC_IsNuMuCC)
FillMyPlots( *MC_CVNMuonLooPTP, *MC_ReMId, *MC_CosRej2020, (
int)*MC_LooPTPQuant, hRecoNuE_Sig, *MC_RecoNuE, EventWgt );
261 else FillMyPlots( *MC_CVNMuonLooPTP, *MC_ReMId, *MC_CosRej2020, (
int)*MC_LooPTPQuant, hRecoNuE_Bkg, *MC_RecoNuE, EventWgt );
264 double MyTrPIDScore = *MC_ReMId;
265 if (!IsReMId) MyTrPIDScore = *MC_CVNProngMuon;
267 double MyCVNScore = *MC_CVNMuonLooPTP;
270 if (
PassCutPoint( MyTrPIDScore,
MinTrPID, MyCVNScore,
MinCVN, *MC_LooPTPQuant, *MC_CosVeto, cQuantVec, *MC_RecoNuE,
false ) ) {
272 if (*MC_IsNuMuCC) hMC_Sig_FullE_Hist->Fill( MyTrPIDScore, MyCVNScore, *MC_CosRej2020, EventWgt );
273 else hMC_Bkg_FullE_Hist->Fill( MyTrPIDScore, MyCVNScore, *MC_CosRej2020, EventWgt );
276 if (
PassCutPoint( MyTrPIDScore,
MinTrPID, MyCVNScore,
MinCVN, *MC_LooPTPQuant, *MC_CosVeto, cQuantVec, *MC_RecoNuE,
true ) ) {
277 if (*MC_IsNuMuCC) hMC_Sig_Dip_Hist->Fill( MyTrPIDScore, MyCVNScore, *MC_CosRej2020, EventWgt );
278 else hMC_Bkg_Dip_Hist->Fill( MyTrPIDScore, MyCVNScore, *MC_CosRej2020, EventWgt );
283 while (Data_Reader.Next()) {
285 int Entry = Data_Reader.GetCurrentEntry();
290 PrintInfo( Entry, Data_TotalTree, *Data_IsNuMuCC, *Data_TrueNuE, *Data_RecoNuE, EventWgt, *Data_CosVeto,
291 *Data_CVNMuonLooPTP, *Data_ReMId, *Data_CVNProngMuon, *Data_CosRej2020 );
295 FillMyPlots( *Data_CVNMuonLooPTP, *Data_ReMId, *Data_CosRej2020, (
int)*Data_LooPTPQuant, hRecoNuE_Cos, *Data_RecoNuE, EventWgt );
298 double MyTrPIDScore = *Data_ReMId;
299 if (!IsReMId) MyTrPIDScore = *Data_CVNProngMuon;
301 double MyCVNScore = *Data_CVNMuonLooPTP;
304 if (
PassCutPoint( MyTrPIDScore,
MinTrPID, MyCVNScore,
MinCVN, *Data_LooPTPQuant, *Data_CosVeto, cQuantVec, *Data_RecoNuE,
false ) ) {
305 hData_FullE_Hist->Fill( MyTrPIDScore, MyCVNScore, *Data_CosRej2020, EventWgt );
308 if (
PassCutPoint( MyTrPIDScore,
MinTrPID, MyCVNScore,
MinCVN, *Data_LooPTPQuant, *Data_CosVeto, cQuantVec, *Data_RecoNuE,
true ) ) {
309 hData_Dip_Hist->Fill ( MyTrPIDScore, MyCVNScore, *Data_CosRej2020, EventWgt );
315 for (
int hh=0;
hh<9; ++
hh) {
321 std::cout << setprecision(3) <<
"\n After scaling, the different cut hists are as follows;" 322 <<
"\n\n\t FOM Full --> ReMId 0.3, CVN 0.85, CosRej 0.45." 323 <<
"\n\t\t AllQuant: Sig " << hRecoNuE_Sig[0]->Integral() <<
", Bkg " << hRecoNuE_Bkg[0]->Integral() <<
", Cos " << hRecoNuE_Cos[0]->Integral()
324 <<
"\n\t\t Quant123: Sig " << hRecoNuE_Sig[1]->Integral() <<
", Bkg " << hRecoNuE_Bkg[1]->Integral() <<
", Cos " << hRecoNuE_Cos[1]->Integral()
325 <<
"\n\t\t Quant4 : Sig " << hRecoNuE_Sig[2]->Integral() <<
", Bkg " << hRecoNuE_Bkg[2]->Integral() <<
", Cos " << hRecoNuE_Cos[2]->Integral()
326 <<
"\n\t\t --> Q4Frac is " << hRecoNuE_Sig[2]->Integral() / hRecoNuE_Sig[0]->Integral()
328 <<
"\n\n\t FOM Dip --> ReMId 0.32, CVN 0.91, CosRej 0.45." 329 <<
"\n\t\t AllQuant: Sig " << hRecoNuE_Sig[3]->Integral() <<
", Bkg " << hRecoNuE_Bkg[3]->Integral() <<
", Cos " << hRecoNuE_Cos[3]->Integral()
330 <<
"\n\t\t Quant123: Sig " << hRecoNuE_Sig[4]->Integral() <<
", Bkg " << hRecoNuE_Bkg[4]->Integral() <<
", Cos " << hRecoNuE_Cos[4]->Integral()
331 <<
"\n\t\t Quant4 : Sig " << hRecoNuE_Sig[5]->Integral() <<
", Bkg " << hRecoNuE_Bkg[5]->Integral() <<
", Cos " << hRecoNuE_Cos[5]->Integral()
332 <<
"\n\t\t --> Q4Frac is " << hRecoNuE_Sig[5]->Integral() / hRecoNuE_Sig[3]->Integral()
334 <<
"\n\n\t Finally, current cut --> ReMId 0.7, CVN 0.82, CosRej 0.47." 335 <<
"\n\t\t AllQuant: Sig " << hRecoNuE_Sig[6]->Integral() <<
", Bkg " << hRecoNuE_Bkg[6]->Integral() <<
", Cos " << hRecoNuE_Cos[6]->Integral()
336 <<
"\n\t\t Quant123: Sig " << hRecoNuE_Sig[7]->Integral() <<
", Bkg " << hRecoNuE_Bkg[7]->Integral() <<
", Cos " << hRecoNuE_Cos[7]->Integral()
337 <<
"\n\t\t Quant4 : Sig " << hRecoNuE_Sig[8]->Integral() <<
", Bkg " << hRecoNuE_Bkg[8]->Integral() <<
", Cos " << hRecoNuE_Cos[8]->Integral()
338 <<
"\n\t\t --> Q4Frac is " << hRecoNuE_Sig[8]->Integral() / hRecoNuE_Sig[6]->Integral()
341 TCanvas *cFOMFull =
new TCanvas(TString(PIDOutString)+TString(
"-FOMFull_QuickRecoE"),
"Reco E for the cut optimised for FOM over All E");
342 hRecoNuE_Sig[0]->Scale(0.1,
"width"); hRecoNuE_Sig[0] ->
Draw();
343 hRecoNuE_Bkg[0]->Scale(0.1,
"width"); hRecoNuE_Bkg[0] ->
Draw(
"same");
344 hRecoNuE_Cos[0]->Scale(0.1,
"width"); hRecoNuE_Cos[0] ->
Draw(
"same");
345 cFOMFull ->
SaveAs( TString(OutDir)+TString(cFOMFull ->
GetName())+TString(
".pdf") );
347 TCanvas *cFOMDip =
new TCanvas(TString(PIDOutString)+TString(
"-FOMDip_QuickRecoE"),
"Reco E for the cut optimised for FOM in Dip");
348 hRecoNuE_Sig[3]->Scale(0.1,
"width"); hRecoNuE_Sig[3] ->
Draw();
349 hRecoNuE_Bkg[3]->Scale(0.1,
"width"); hRecoNuE_Bkg[3] ->
Draw(
"same");
350 hRecoNuE_Cos[3]->Scale(0.1,
"width"); hRecoNuE_Cos[3] ->
Draw(
"same");
351 cFOMDip ->
SaveAs( TString(OutDir)+TString(cFOMDip ->
GetName())+TString(
".pdf") );
353 TCanvas *cFOMCurrCut =
new TCanvas(TString(PIDOutString)+TString(
"-FOMCurrCut_QuickRecoE"),
"Reco E for the current cut");
354 hRecoNuE_Sig[6]->Scale(0.1,
"width"); hRecoNuE_Sig[6] ->
Draw();
355 hRecoNuE_Bkg[6]->Scale(0.1,
"width"); hRecoNuE_Bkg[6] ->
Draw(
"same");
356 hRecoNuE_Cos[6]->Scale(0.1,
"width"); hRecoNuE_Cos[6] ->
Draw(
"same");
357 cFOMCurrCut ->
SaveAs( TString(OutDir)+TString(cFOMCurrCut->GetName())+TString(
".pdf") );
359 TCanvas *cFOMFull_Q123 =
new TCanvas(TString(PIDOutString)+TString(
"-FOMFull_Q123"),
"Reco E for the cut optimised for FOM over All E - Q123");
360 hRecoNuE_Sig[1]->Scale(0.1,
"width"); hRecoNuE_Sig[1] ->
Draw();
361 hRecoNuE_Bkg[1]->Scale(0.1,
"width"); hRecoNuE_Bkg[1] ->
Draw(
"same");
362 hRecoNuE_Cos[1]->Scale(0.1,
"width"); hRecoNuE_Cos[1] ->
Draw(
"same");
363 cFOMFull_Q123->
SaveAs( TString(OutDir)+TString(cFOMFull_Q123->GetName())+TString(
".pdf") );
365 TCanvas *cFOMFull_Q4 =
new TCanvas(TString(PIDOutString)+TString(
"-FOMFull_Q4"),
"Reco E for the cut optimised for FOM over All E - Q4");
366 hRecoNuE_Sig[2]->Scale(0.1,
"width"); hRecoNuE_Sig[2] ->
Draw();
367 hRecoNuE_Bkg[2]->Scale(0.1,
"width"); hRecoNuE_Bkg[2] ->
Draw(
"same");
368 hRecoNuE_Cos[2]->Scale(0.1,
"width"); hRecoNuE_Cos[2] ->
Draw(
"same");
369 cFOMFull_Q4->
SaveAs( TString(OutDir)+TString(cFOMFull_Q4->GetName())+TString(
".pdf") );
379 double FOMMin_f, FOMMax_f, FOMMin_d, FOMMax_d;
380 int nbinsFOM_f, nbinsFOM_d;
381 if ( IsFHC && QuantEnum == 0 ) { FOMMin_f = 6.; FOMMax_f = 12.; nbinsFOM_f = 60; FOMMin_d = 3.; FOMMax_d = 5.5; nbinsFOM_d = 70; }
382 else if ( IsFHC && QuantEnum == 123 ) { FOMMin_f = 6.; FOMMax_f = 10.; nbinsFOM_f = 40; FOMMin_d = 3.; FOMMax_d = 4.2; nbinsFOM_d = 24; }
383 else if ( IsFHC && QuantEnum == 4 ) { FOMMin_f = 3.; FOMMax_f = 6.5; nbinsFOM_f = 35; FOMMin_d = 1.; FOMMax_d = 3.5; nbinsFOM_d = 50; }
384 else if (!IsFHC && QuantEnum == 0 ) { FOMMin_f = 6.; FOMMax_f = 10.; nbinsFOM_f = 40; FOMMin_d = 2.; FOMMax_d = 4.5; nbinsFOM_d = 50; }
385 else if (!IsFHC && QuantEnum == 123 ) { FOMMin_f = 6.; FOMMax_f = 8.5; nbinsFOM_f = 25; FOMMin_d = 2.; FOMMax_d = 4 ; nbinsFOM_d = 40; }
386 else if (!IsFHC && QuantEnum == 4 ) { FOMMin_f = 3.; FOMMax_f = 5.5; nbinsFOM_f = 25; FOMMin_d = 0.; FOMMax_d = 2.5; nbinsFOM_d = 50; }
387 TH2 *hFOMFull_Bkg =
new TH2D(
"FOMFull_Bkg_Hist",
"Full FOM vs Number of Bkg Events; Total Background; FOM", 150, 0, 30, nbinsFOM_f, FOMMin_f, FOMMax_f);
388 TH2 *hFOMDip_Bkg =
new TH2D(
"FOMDip_Bkg_Hist" ,
"Dip FOM vs Number of Bkg Events; Total Background; FOM" , 150, 0, 15, nbinsFOM_d, FOMMin_d, FOMMax_d);
392 TH2 *hGoodFOM_TrPID_CVN =
new TH2D(
"GoodFOM_TrPID_CVN" ,
"Density of TrackPID vs CVN Cuts above FOM Thresholds; Track PID Cut; CVN Cut",
394 TH2 *hGoodFOM_TrPID_CosRej =
new TH2D(
"GoodFOM_TrPID_CosRej",
"Density of TrackPID vs CosRej Cuts above FOM Thresholds; Track PID Cut; CosRej Cut",
396 TH2 *hGoodFOM_CVN_CosRej =
new TH2D(
"GoodFOM_CVN_CosRej" ,
"Density of CosRej vs CVN Cuts above FOM Thresholds; CVN Cut; CosRej Cut",
400 double MaxFOM_Full = 0, MaxFOM_Dip = 0;
402 double FOMInt_Full = 0, FOMInt_Dip = 0;
404 if ( IsFHC && QuantEnum == 0 ) { FOMInt_Full = 10.8; FOMInt_Dip = 4.8; }
405 else if ( IsFHC && QuantEnum == 123 ) { FOMInt_Full = 9.1 ; FOMInt_Dip = 3.9; }
406 else if ( IsFHC && QuantEnum == 4 ) { FOMInt_Full = 5.9 ; FOMInt_Dip = 2.8; }
408 else if (!IsFHC && QuantEnum == 0 ) { FOMInt_Full = 9.5 ; FOMInt_Dip = 3.9; }
409 else if (!IsFHC && QuantEnum == 123 ) { FOMInt_Full = 7.9 ; FOMInt_Dip = 3.4; }
410 else if (!IsFHC && QuantEnum == 4 ) { FOMInt_Full = 5.2 ; FOMInt_Dip = 2.0; }
413 std::cout <<
"\tNow looking at combinations with trpid " << trpid <<
". FOMMap currently has size " << FOMMap.size() <<
std::endl;
417 int bin_trpid = hData_FullE_Hist->GetXaxis()->FindBin( trpid );
418 int bin_cvn = hData_FullE_Hist->GetYaxis()->FindBin(
cvn );
419 int bin_cos = hData_FullE_Hist->GetZaxis()->FindBin(
cos );
434 if ( TempFOM.
FOM_Full > FOMInt_Full || TempFOM.
FOM_Dip > FOMInt_Dip ) {
436 std::ostringstream CutString;
437 CutString << std::setprecision(2) << std::fixed <<
"TrPID-" << trpid <<
"-CVN-" <<
cvn <<
"-CosRej-" <<
cos;
440 || ( bin_trpid==1 && bin_cvn==46 && bin_cos==6 )
441 || ( bin_trpid==1 && bin_cvn==52 && bin_cos==6 )
442 || ( bin_trpid==41 && bin_cvn==43 && bin_cos==8 )
445 <<
"Looking at My TrPID " << trpid <<
", My CVN " <<
cvn <<
", CosRej2020 " << cos <<
" --> " << CutString.str()
446 <<
"\n\t bin_rem is " << bin_trpid <<
", bin_cvn is " << bin_cvn <<
", bin_cos is " << bin_cos
452 FOMMap[ CutString.str() ] = TempFOM;
454 hGoodFOM_TrPID_CVN ->
Fill( trpid,
cvn );
455 hGoodFOM_TrPID_CosRej ->
Fill( trpid, cos );
456 hGoodFOM_CVN_CosRej ->
Fill(
cvn , cos );
464 TCanvas* cFOMFull_Bkg =
new TCanvas(TString(PIDOutString)+TString(
"-FOMFull_Bkg"),
"FOM over All E vs Total Bkg");
465 hFOMFull_Bkg->Draw(
"colz");
466 TCanvas* cFOMDip_Bkg =
new TCanvas(TString(PIDOutString)+TString(
"-FOMDip_Bkg") ,
"FOM in Dip Region vs Total Bkg");
467 hFOMDip_Bkg ->Draw(
"colz");
470 TCanvas *cGoodFOM_TrPID_CVN =
new TCanvas(TString(PIDOutString)+TString(
"-GoodFOM_TrPID_CVN") ,
"Density of TrackPID vs CVN Cuts above FOM Thresholds" );
471 hGoodFOM_TrPID_CVN ->
Draw(
"colz");
472 TCanvas *cGoodFOM_TrPID_CosRej =
new TCanvas(TString(PIDOutString)+TString(
"-GoodFOM_TrPID_CosRej"),
"Density of TrackPID vs CosRej Cuts above FOM Thresholds" );
473 hGoodFOM_TrPID_CosRej ->
Draw(
"colz");
474 TCanvas *cGoodFOM_CVN_CosRej =
new TCanvas(TString(PIDOutString)+TString(
"-GoodFOM_CVN_CosRej") ,
"Density of CVN vs CosRej Cuts above FOM Thresholds" );
475 hGoodFOM_CVN_CosRej ->
Draw(
"colz");
478 cFOMFull_Bkg ->
SaveAs( TString(OutDir)+TString(cFOMFull_Bkg ->
GetName())+TString(
".pdf") );
479 cFOMDip_Bkg ->
SaveAs( TString(OutDir)+TString(cFOMDip_Bkg ->
GetName())+TString(
".pdf") );
480 cGoodFOM_TrPID_CVN ->
SaveAs( TString(OutDir)+TString(cGoodFOM_TrPID_CVN ->
GetName())+TString(
".pdf") );
481 cGoodFOM_TrPID_CosRej->
SaveAs( TString(OutDir)+TString(cGoodFOM_TrPID_CosRej->GetName())+TString(
".pdf") );
482 cGoodFOM_CVN_CosRej ->
SaveAs( TString(OutDir)+TString(cGoodFOM_CVN_CosRej ->
GetName())+TString(
".pdf") );
486 std::cout <<
"\n After all of that. The MaxFOM_Full is " << MaxFOM_Full <<
", the MaxFOM_Dip is " << MaxFOM_Dip
488 <<
" Lets look at what cuts are best for different criteria." 490 const int PrintLim = 10;
int Printed;
496 std::multimap< double, std::string, std::greater<double> > FOMMap_SortFullFOM;
497 for (std::pair< std::string, MyFOMParams > FOMMapIt : FOMMap) {
498 FOMMap_SortFullFOM.insert(
std::make_pair( FOMMapIt.second.FOM_Full, FOMMapIt.first ) );
500 for (std::pair< double, std::string > SortFOMIt : FOMMap_SortFullFOM) {
501 if (Printed > PrintLim)
break;
510 std::multimap< double, std::string, std::greater<double> > FOMMap_SortDipFOM;
511 for (std::pair< std::string, MyFOMParams > FOMMapIt : FOMMap) {
512 FOMMap_SortDipFOM.insert(
std::make_pair( FOMMapIt.second.FOM_Dip, FOMMapIt.first ) );
514 for (std::pair< double, std::string > SortFOMIt : FOMMap_SortDipFOM) {
515 if (Printed > PrintLim)
break;
524 std::multimap< double, std::string, std::greater<double> > FOMMap_SortSig_Dip;
525 for (std::pair< std::string, MyFOMParams > FOMMapIt : FOMMap) {
526 FOMMap_SortSig_Dip.insert(
std::make_pair( FOMMapIt.second.TotalSig_Dip, FOMMapIt.first ) );
528 for (std::pair< double, std::string > SortFOMIt : FOMMap_SortSig_Dip) {
529 if (Printed > PrintLim)
break;
538 std::multimap< double, std::string > FOMMap_SortBkg_Dip;
539 for (std::pair< std::string, MyFOMParams > FOMMapIt : FOMMap) {
540 FOMMap_SortBkg_Dip.insert(
std::make_pair( FOMMapIt.second.TotalBkg_Dip, FOMMapIt.first ) );
542 for (std::pair< double, std::string > SortFOMIt : FOMMap_SortBkg_Dip) {
543 if (Printed > PrintLim)
break;
553 void PrintInfo(
int Entry,
int TotEnt,
double IsNuMuCC,
double TrueNuE,
double RecoNuE,
double EventWgt,
double CosVeto,
554 double CVNMuonLooPTP,
double ReMId,
double CVNProngMuon,
double CosRej2020 ) {
557 <<
"\n Looking at entry " << Entry <<
" of " << TotEnt <<
" ---> " << 100*(double)Entry/(
double)TotEnt <<
"% complete... \n\t" 558 <<
" IsNuMuCC " << IsNuMuCC <<
", TrueNuE " << TrueNuE <<
", RecoNuE " << RecoNuE <<
", EventWgt " << EventWgt
559 <<
", Pass CosVeto? " << CosVeto
560 <<
", CVNMuonLooPTP " << CVNMuonLooPTP
562 <<
", ReMId " << ReMId
563 <<
", CVNProngMuon " << CVNProngMuon
564 <<
", CosRej2020 " << CosRej2020
569 void CalcFOM( TH3D* Sig, TH3D*
Bkg, TH3D* Data,
int bintrpid,
int bincvn,
int bincos,
bool IsFD,
double &FOM,
570 double &TotSig,
double &
TotBkg ) {
578 TotBkg = BkgInt + DataInt;
581 FOM = TotSig /
std::pow( TotSig + TotBkg, 0.5 );
587 std::cout <<
"\t"<<Key << setprecision(4)
588 <<
", FOM_Full " << MainMap[ Key ].FOM_Full
589 <<
", TotalSig_Full: " << MainMap[ Key ].TotalSig_Full
590 <<
", TotalBkg_Full: " << MainMap[ Key ].TotalBkg_Full
591 <<
", FOM_Dip " << MainMap[ Key ].FOM_Dip
592 <<
", TotalSig_Dip: " << MainMap[ Key ].TotalSig_Dip
593 <<
", TotalBkg_Dip: " << MainMap[ Key ].TotalBkg_Dip
598 bool PassCutPoint(
double PIDScore,
double PIDCut,
double CVNScore,
double CVNCut,
599 int Quant,
int CosVeto, std::vector<int> cQuantVec,
double RecoE,
bool isDip ) {
602 if ( PIDScore >= PIDCut && CVNScore >= CVNCut
606 for (
size_t cq=0; cq < cQuantVec.size(); ++cq) {
608 if (Quant == cQuantVec[cq]) {
610 if ( (RecoE > 1 && RecoE < 2) || !isDip )
return true;
618 void FillMyPlots (
double CVN,
double ReMId,
double CosRej,
int Quant, TH1D* Hists[9],
double RecoE,
double Wgt ) {
620 if ( ReMId > 0.3 && CVN > 0.85 && CosRej > 0.45 ) {
621 Hists[0] ->
Fill( RecoE, Wgt );
622 if (Quant != 4) Hists[1] ->
Fill( RecoE, Wgt );
623 else Hists[2] ->
Fill( RecoE, Wgt );
626 if ( ReMId > 0.32 && CVN > 0.91 && CosRej > 0.45 ) {
627 Hists[3] ->
Fill( RecoE, Wgt );
628 if (Quant != 4) Hists[4] ->
Fill( RecoE, Wgt );
629 else Hists[5] ->
Fill( RecoE, Wgt );
632 if ( ReMId > 0.7 && CVN > 0.82 && CosRej > 0.47 ) {
633 Hists[6] ->
Fill( RecoE, Wgt );
634 if (Quant != 4) Hists[7] ->
Fill( RecoE, Wgt );
635 else Hists[8] ->
Fill( RecoE, Wgt );
const double FHC_POT_Scale
double Integral(const Spectrum &s, const double pot, int cvnregion)
const double RHC_Liv_Scale
bool PassCutPoint(double PIDScore, double PIDCut, double CVNScore, double CVNCut, int Quant, int CosVeto, std::vector< int > cQuantVec, double RecoE, bool isDip)
double TotBkg(std::vector< double > integrals)
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
const double FHC_Liv_Scale
Defines an enumeration for prong classification.
void PrintInfo(int Entry, int TotEnt, double IsNuMuCC, double TrueNuE, double RecoNuE, double EventWgt, double CosVeto, double CVNMuonLooPTP, double ReMId, double CVNProngMuon, double CosRej2020)
std::string GetName(int i)
void CalcFOM(TH3D *Sig, TH3D *Bkg, TH3D *Data, int bintrpid, int bincvn, int bincos, bool IsFD, double &FOM, double &TotSig, double &TotBkg)
correl_yv Fill(-(dy[iP-1][iC-1]), hyv->GetBinContent(iP, iC))
double RecoEBins[HistBins+1]
const double RHC_POT_Scale
void PrintMyMap(std::map< std::string, MyFOMParams > MainMap, std::string Key)
cosmicTree SaveAs("cosmicTree.root")
void NuMu2020_MakePIDPlots(bool IsFD, bool IsFHC, int QuantEnum, bool IsReMId=true, bool IsLoosePTP=true)
simulatedPeak Scale(1/simulationNormalisationFactor)
void FillMyPlots(double CVN, double ReMId, double CosRej, int Quant, TH1D *Hists[9], double RecoE, double Wgt)