NuMu2020_MakePIDPlots.C
Go to the documentation of this file.
1 /*
2  Karl W: Written for the 2020 PID NuMu Optimisation.
3  To run:
4  root
5  .x NuMu2020_MakePIDPlots.C( <bool for FD> ,
6  <bool for Horn Current>,
7  <enum for Quantile to use>, { 0 = AllQuants, 123 = Quants 1,2,3 and 4 = Just Quant 4 }
8  <bool for ReMId/ProngCVN>, Defaulted to true
9  <bool for LoosePTP/OldPresel> Defaulted to true )
10  Will run over either the ND or the FD.
11  Will run over either FHC or RHC data.
12  Will fetch either ReMId or Prong CVN Scores
13  Will fetch either LoosePTP or OldPresel trained CVN Scores.
14 */
15 
16 #include "TFile.h"
17 #include "TH1.h"
18 #include "TH3.h"
19 #include "TTreeReader.h"
20 #include "TTreeReaderValue.h"
21 #include "TSystem.h"
22 
23 #include <iostream>
24 #include <iomanip>
25 #include <stdio.h>
26 #include <string>
27 
28 #include <TStyle.h>
29 #include <TCanvas.h>
30 
31 //#include "Utilities/rootlogon.C"
32 
33 using namespace std;
34 
35 // ---------------------------------------------------------------------------------------------------------------------
36 // --- What are the scales to use?
37 // The official POTs and Livetimes are taken from Analysis/Exposures.h, and are the Ana19 values.
38 // The production POTs and Livetimes are taken from a Spectra ran over the entire dataset in NuMu2020_BasicPlots_*.C.
39 const double FHC_POT_Scale = 9.47964e20 / 1.7980696e24;
40 const double FHC_Liv_Scale = 438.1584 / 57385.145;
41 
42 const double RHC_POT_Scale = 12.3319e20 / 4.5758019e24;
43 const double RHC_Liv_Scale = 317.0466 / 43780.253;
44 
45 double MyPOTScale = 0;
46 double MyLivScale = 0;
47 // ---------------------------------------------------------------------------------------------------------------------
48 const double MinTrPID = 0.3, MaxTrPID = 0.95; const int nbinTrPID = 65;
49 const double MinCVN = 0.4, MaxCVN = 0.95; const int nbinCVN = 55;
50 const double MinCosRej = 0.4, MaxCosRej = 0.65; const int nbinCosRej = 25;
51 // ---------------------------------------------------------------------------------------------------------------------
52 const bool bPrintInfo = false; // Dirty, but easy to toggle how much info I want...
53 
54 // ---------------------------------------------------------------------------------------------------------------------
55 // A struct to store key FOM values.
56 struct MyFOMParams {
57  // Numbers for the full dist.
58  double FOM_Full;
59  double TotalSig_Full;
60  double TotalBkg_Full;
61  // Numbers for the dip region.
62  double FOM_Dip;
63  double TotalSig_Dip;
64  double TotalBkg_Dip;
65 };
66 // ---------------------------------------------------------------------------------------------------------------------
67 
68 // ---------------------------------------------------------------------------------------------------------------------
69 // Declare any functions that I'm going to use later
70 // ---------------------------------------------------------------------------------------------------------------------
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 );
73 // ---------------------------------------------------------------------------------------------------------------------
74 void CalcFOM( TH3D* Sig, TH3D* Bkg, TH3D* Data, int bintrpid, int bincvn, int bincos, bool IsFD,
75  double &FOM, double &TotSig, double &TotBkg );
76 // ---------------------------------------------------------------------------------------------------------------------
77 void PrintMyMap( std::map< std::string, MyFOMParams > MainMap, std::string Key );
78 // ---------------------------------------------------------------------------------------------------------------------
79 bool PassCutPoint( double PIDScore, double PIDCut, double CVNScore, double CVNCut,
80  int Quant, int CosVeto, std::vector<int> cQuantVec, double RecoE, bool isDip );
81 // ---------------------------------------------------------------------------------------------------------------------
82 void FillMyPlots ( double CVN, double ReMId, double CosRej, int Quant, TH1D* Hists[9], double RecoE, double Wgt );
83 // ---------------------------------------------------------------------------------------------------------------------
84 // --- Make the axis for my TH1Ds if making Reco Nu E plots.
85 const int HistBins = 19;
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 };
87 
88 // ---------------------------------------------------------------------------------------------------------------------
89 // The main class
90 // ---------------------------------------------------------------------------------------------------------------------
91 void NuMu2020_MakePIDPlots( bool IsFD, bool IsFHC, int QuantEnum, bool IsReMId=true, bool IsLoosePTP=true ) {
92  // --- Quant Enum
93  // 0 - all quants together.
94  // 123 - quantiles 1, 2 and 3 together.
95  // 4 - quantile 4 by itself.
96 
97  // Remove stat boxes.
98  gStyle->SetOptStat(0);
99 
100  // FD or ND? This affects whether the "data" is real data, or cosmic data.
101  std::string sIsFD = "FD", sData = "Cosm", sMC = "MC";
102  if (!IsFD) { sIsFD = "ND", sData = "Data"; }
103  // FHC or RHC?
104  std::string sIsFHC = "FHC";
105  if (!IsFHC) sIsFHC = "RHC";
106  // Am I using ReMId or Prong CVN?
107  std::string sTrPID = "ReMId";
108  if (!IsReMId) sTrPID = "ProngCVN";
109  // Am I using the LoosePTP CVN or the OldPresel CVN?
110  std::string sCVN = "LoosePTP";
111  if (!IsLoosePTP){sCVN = "OldPresel"; std::cout << "\nHave you reenabled OldPresel CVN?\n" << std::endl; return; }
112  // What is my quantile enum?
113  std::string sQuant = "";
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; }
119  // Put it all together to make a string
120  std::string PIDOutString = sIsFD+"-"+sIsFHC+"-"+sTrPID+"-"+sCVN+"-"+sQuant;
121 
122  // Where are the files that I'm loading in?
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";
125 
126  // --> Set my POT and Livetime accordingly to the above..
127  if (IsFHC) { MyPOTScale = FHC_POT_Scale; MyLivScale = FHC_Liv_Scale; }
129 
130  // --- Print out some key info to check things are correct.
131  std::cout << "\n Is this FD? " << IsFD << "\t Is this FHC? " << IsFHC << "\t QuantEnum? " << QuantEnum
132  << "\n\t --> POT Scale is " << MyPOTScale << ", Livetime Scale is " << MyLivScale
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
138  << "\n" << std::endl;
139 
140  // --- Define my OutDir and make sure it exists.
141  std::string OutDir = "Plots-"+PIDOutString+"/";
142  gSystem -> MakeDirectory( OutDir.c_str() );
143 
144  // ---- Define my TFiles.
145  TFile *MyMCSamp = new TFile( MCFile .c_str() );
146  TFile *MyDataSamp = new TFile( DataFile.c_str() );
147 
148  // --- Load the MC TTreeReader, and set its branches
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" );
161  //TTreeReaderValue<double> MC_OldPreQuant (MC_Reader, "OldPreQuant" );
162  //TTreeReaderValue<double> MC_CVNMuonOldPre(MC_Reader, "CVNMuonOldPre");
163  //TTreeReaderValue<double> MC_CVNCosmOldPre(MC_Reader, "CVNCosmOldPre");
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" );
171  //TTreeReaderValue<double> MC_CosRej2018 (MC_Reader, "CosRej2018" );
172 
173  // --- Load the Data TTreeReader, and set its branches
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" );
186  //TTreeReaderValue<double> Data_OldPreQuant (Data_Reader, "OldPreQuant" );
187  //TTreeReaderValue<double> Data_CVNMuonOldPre(Data_Reader, "CVNMuonOldPre");
188  //TTreeReaderValue<double> Data_CVNCosmOldPre(Data_Reader, "CVNCosmOldPre");
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" );
196  //TTreeReaderValue<double> Data_CosRej2018 (Data_Reader, "CosRej2018" );
197 
198  // --- Make my 3D histograms.
199  // First the full 0 to 5 GeV histograms.
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",
206  // Next the dip region hists - 1 to 2 GeV.
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",
213 
214 
215  TH1D* hRecoNuE_Sig[9];
216  TH1D* hRecoNuE_Bkg[9];
217  TH1D* hRecoNuE_Cos[9];
218  for (int hh=0; hh<9; ++hh) {
219  std::string app = "";
220  int Col = 0;
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 );
234  hRecoNuE_Sig[hh] -> SetLineColor( 1 );
235  hRecoNuE_Bkg[hh] -> SetLineColor( 2 );
236  hRecoNuE_Cos[hh] -> SetLineColor( 4 );
237  }
238 
239  // Declare two maps to store the FOMs.
240  std::map< std::string, MyFOMParams > FOMMap;
241 
242  int NFills = 0;
243 
244  // Just access the data as if myPx and myPy were iterators (note the '*' in front of them):
245  // ---- Monte Carlo Reader Loop
246  while (MC_Reader.Next()) {
247  // What entry am I looking at?
248  int Entry = MC_Reader.GetCurrentEntry();
249  // What is my Event Weight?
250  //double EventWgt = *MC_PPFXWgt * *MC_XSecWgt * *MC_DumbOsc;
251  double EventWgt = *MC_PPFXWgt * *MC_XSecWgt * *MC_BestFit2019;
252  //double EventWgt = *MC_PPFXWgt * *MC_XSecWgt;
253  // Do I want to write out some useful information?
254  if ( bPrintInfo && Entry % 50000 == 0 ) {
255  PrintInfo( Entry, MC_TotalTree, *MC_IsNuMuCC, *MC_TrueNuE, *MC_RecoNuE, EventWgt, *MC_CosVeto,
256  *MC_CVNMuonLooPTP, *MC_ReMId, *MC_CVNProngMuon, *MC_CosRej2020 );
257  }
258 
259  // --- Fill my Reco plots...
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 );
262 
263  // Set the Track PID score that I want to look at...
264  double MyTrPIDScore = *MC_ReMId;
265  if (!IsReMId) MyTrPIDScore = *MC_CVNProngMuon;
266  // Set the CVN score that I want to look at...
267  double MyCVNScore = *MC_CVNMuonLooPTP;
268  //if (!IsLoosePTP) MyCVNScore = *MC_CVNMuonOldPre;
269  // ---- Does this event pass the cut points? If so, fill the Full E plots.
270  if ( PassCutPoint( MyTrPIDScore, MinTrPID, MyCVNScore, MinCVN, *MC_LooPTPQuant, *MC_CosVeto, cQuantVec, *MC_RecoNuE, false ) ) {
271  ++NFills; // increment the NFills Var.
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 );
274  }
275  // ---- Does this event pass the cut points? If so, fill the Dip E plots.
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 );
279  } // Only fill if passes a loose presel.
280  } // MC Reader
281 
282  // ---- Data Reader Loop
283  while (Data_Reader.Next()) {
284  // What entry am I looking at?
285  int Entry = Data_Reader.GetCurrentEntry();
286  // What is my Event Weight? --> Always 1 for Data.
287  double EventWgt = 1;
288  // Do I want to write out some useful information?
289  if ( bPrintInfo && Entry % 10000 == 0 ) {
290  PrintInfo( Entry, Data_TotalTree, *Data_IsNuMuCC, *Data_TrueNuE, *Data_RecoNuE, EventWgt, *Data_CosVeto,
291  *Data_CVNMuonLooPTP, *Data_ReMId, *Data_CVNProngMuon, *Data_CosRej2020 );
292  }
293 
294  // --- Fill my Reco plots...
295  FillMyPlots( *Data_CVNMuonLooPTP, *Data_ReMId, *Data_CosRej2020, (int)*Data_LooPTPQuant, hRecoNuE_Cos, *Data_RecoNuE, EventWgt );
296 
297  // Set the Track PID score that I want to look at...
298  double MyTrPIDScore = *Data_ReMId;
299  if (!IsReMId) MyTrPIDScore = *Data_CVNProngMuon;
300  // Set the CVN score that I want to look at...
301  double MyCVNScore = *Data_CVNMuonLooPTP;
302  //if (!IsLoosePTP) MyCVNScore = *Data_CVNMuonOldPre;
303  // ---- Does this event pass the cut points? If so, fill the Full E plots.
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 );
306  }
307  // ---- Does this event pass the cut points? If so, fill the Dip E plots.
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 );
310  }
311  } // Data Reader
312 
313  // ----- Now do my dirty plotting things....
314 
315  for (int hh=0; hh<9; ++hh) {
316  hRecoNuE_Sig[hh] -> Scale( MyPOTScale );
317  hRecoNuE_Bkg[hh] -> Scale( MyPOTScale );
318  hRecoNuE_Cos[hh] -> Scale( MyLivScale );
319  }
320 
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()
327 
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()
333 
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()
339  << "\n\n" << std::endl;
340 
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") );
346 
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") );
352 
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") );
358 
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") );
364 
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") );
370 
371  //return;
372 
373  // ---- AAAnnnd scale the plots...
374  hMC_Sig_FullE_Hist->Scale(MyPOTScale); hMC_Sig_Dip_Hist->Scale(MyPOTScale);
375  hMC_Bkg_FullE_Hist->Scale(MyPOTScale); hMC_Bkg_Dip_Hist->Scale(MyPOTScale);
376  hData_FullE_Hist ->Scale(MyLivScale); hData_Dip_Hist ->Scale(MyLivScale);
377 
378  // --- Declare some 2D FOM vs background plots.
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);
389 
390  // --- Declare some VarX vs Var Y plots to give a feel for the density of good cut points.
391  // First the full E Dist.
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",
398 
399  // ---- Figure out my FOMs.
400  double MaxFOM_Full = 0, MaxFOM_Dip = 0;
401  // ---- What are my FOM thresholds for an "interesting event"
402  double FOMInt_Full = 0, FOMInt_Dip = 0;
403  // FHC
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; }
407  // RHC
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; }
411  // And loop through my combinations.
412  for ( double trpid = MinTrPID ; trpid < MaxTrPID ; trpid=trpid+0.01) {
413  std::cout << "\tNow looking at combinations with trpid " << trpid << ". FOMMap currently has size " << FOMMap.size() << std::endl;
414  for ( double cvn = MinCVN ; cvn < MaxCVN ; cvn =cvn +0.01) {
415  for (double cos = MinCosRej; cos < MaxCosRej; cos =cos +0.01) {
416  // Calculate the bins to integrate over.
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 );
420  // Set the variables to send to CalcFOM.
421  MyFOMParams TempFOM;
422  // Calculate my FOMs.
423  CalcFOM( hMC_Sig_FullE_Hist, hMC_Bkg_FullE_Hist, hData_FullE_Hist, bin_trpid, bin_cvn, bin_cos, false, TempFOM.FOM_Full, TempFOM.TotalSig_Full, TempFOM.TotalBkg_Full );
424 
425  CalcFOM( hMC_Sig_Dip_Hist , hMC_Bkg_Dip_Hist , hData_Dip_Hist , bin_trpid, bin_cvn, bin_cos, true , TempFOM.FOM_Dip , TempFOM.TotalSig_Dip , TempFOM.TotalBkg_Dip );
426  // Is this the largest FOM?
427  if (TempFOM.FOM_Full > MaxFOM_Full) MaxFOM_Full = TempFOM.FOM_Full;
428  if (TempFOM.FOM_Dip > MaxFOM_Dip ) MaxFOM_Dip = TempFOM.FOM_Dip;
429  // Fill my 2D plots.
430  hFOMFull_Bkg -> Fill( TempFOM.TotalBkg_Full, TempFOM.FOM_Full );
431  hFOMDip_Bkg -> Fill( TempFOM.TotalBkg_Dip , TempFOM.FOM_Dip );
432 
433  // --- If the FOM is good, and only if it is good, add it to my Map for closer inspection later...
434  if ( TempFOM.FOM_Full > FOMInt_Full || TempFOM.FOM_Dip > FOMInt_Dip ) {
435  // What is my cut string for this combination?
436  std::ostringstream CutString;
437  CutString << std::setprecision(2) << std::fixed << "TrPID-" << trpid << "-CVN-" << cvn << "-CosRej-" << cos;
438  // Write out the cut poisiton info.
439  if (bPrintInfo
440  || ( bin_trpid==1 && bin_cvn==46 && bin_cos==6 ) // 0.3, 0.85, 0.45
441  || ( bin_trpid==1 && bin_cvn==52 && bin_cos==6 ) // 0.3, 0.91, 0.45
442  || ( bin_trpid==41 && bin_cvn==43 && bin_cos==8 ) // 0.7, 0.82, 0.47
443  ) {
444  std::cout << setprecision(4)
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
447  << "\n\t TotSig_Full is " << TempFOM.TotalSig_Full << ", TotBkg_Full is " << TempFOM.TotalBkg_Full << " ---> FOM_Full is " << TempFOM.FOM_Full
448  << "\n\t TotSig_Dip is " << TempFOM.TotalSig_Dip << ", TotBkg_Dip is " << TempFOM.TotalBkg_Dip << " ---> FOM_Dip is " << TempFOM.FOM_Dip
449  << std::endl;
450  }
451  // Fill my maps.
452  FOMMap[ CutString.str() ] = TempFOM;
453  // Fill the 2D VarX vs VarY plots.
454  hGoodFOM_TrPID_CVN -> Fill( trpid, cvn );
455  hGoodFOM_TrPID_CosRej -> Fill( trpid, cos );
456  hGoodFOM_CVN_CosRej -> Fill( cvn , cos );
457  }
458 
459  } // loop over cosrej
460  } // loop over cvn
461  } // loop over remid
462 
463  // --- Canvas for the FOM vs Bkg plots.
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");
468 
469  // --- Canvas for the VarX vs VarY Good FOM density plots.
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");
476 
477  // --- Save the canvases.
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") );
483 
484 
485  // ---- After all of that how many combinations were considered "good"?
486  std::cout << "\n After all of that. The MaxFOM_Full is " << MaxFOM_Full << ", the MaxFOM_Dip is " << MaxFOM_Dip
487  << "\n I considered " << FOMMap.size() << " of " << nbinTrPID*nbinCVN*nbinCosRej << " cut combinations as good."
488  << " Lets look at what cuts are best for different criteria."
489  << std::endl;
490  const int PrintLim = 10; int Printed;
491 
492  // ------------------------------------------------
493  // --- Invert the map to be <FOM_Full, CutString>.
494  std::cout << "\nFirst, Lets sort the Map by FOM_Full..." << std::endl;
495  Printed = 0;
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 ) );
499  }
500  for (std::pair< double, std::string > SortFOMIt : FOMMap_SortFullFOM) {
501  if (Printed > PrintLim) break;
502  PrintMyMap( FOMMap, SortFOMIt.second );
503  ++Printed;
504  }
505 
506  // ------------------------------------------------
507  // --- Invert the map to be <FOM_Dip, CutString>.
508  std::cout << "\nNext, Lets sort the Map by FOM_Dip..." << std::endl;
509  Printed = 0;
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 ) );
513  }
514  for (std::pair< double, std::string > SortFOMIt : FOMMap_SortDipFOM) {
515  if (Printed > PrintLim) break;
516  PrintMyMap( FOMMap, SortFOMIt.second );
517  ++Printed;
518  }
519 
520  // ------------------------------------------------
521  // --- Invert the map to be <TotalSig_Dip, CutString>.
522  std::cout << "\nNext, Lets sort the Map by TotalSig_Dip..." << std::endl;
523  Printed = 0;
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 ) );
527  }
528  for (std::pair< double, std::string > SortFOMIt : FOMMap_SortSig_Dip) {
529  if (Printed > PrintLim) break;
530  PrintMyMap( FOMMap, SortFOMIt.second );
531  ++Printed;
532  }
533 
534  // ------------------------------------------------
535  // --- Invert the map to be <TotalBkg_Dip, CutString>.
536  std::cout << "\nNext, Lets sort the Map by TotalBkg_Dip..." << std::endl;
537  Printed = 0;
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 ) );
541  }
542  for (std::pair< double, std::string > SortFOMIt : FOMMap_SortBkg_Dip) {
543  if (Printed > PrintLim) break;
544  PrintMyMap( FOMMap, SortFOMIt.second );
545  ++Printed;
546  }
547 
548  return;
549 } // Main Class.
550 
551 // ---------------------------------------------------------------------------------------------------------------------
552 // ---------------------------------------------------------------------------------------------------------------------
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 ) {
555 
556  std::cout << setprecision(4)
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
561  //<< ", CVNMuonOldPre " << CVNMuonOldPre
562  << ", ReMId " << ReMId
563  << ", CVNProngMuon " << CVNProngMuon
564  << ", CosRej2020 " << CosRej2020
565  << std::endl;
566 } // PrintInfo
567 // ---------------------------------------------------------------------------------------------------------------------
568 // ---------------------------------------------------------------------------------------------------------------------
569 void CalcFOM( TH3D* Sig, TH3D* Bkg, TH3D* Data, int bintrpid, int bincvn, int bincos, bool IsFD, double &FOM,
570  double &TotSig, double &TotBkg ) {
571  // Calc the integrals for each histogram.
572  double SigInt = Sig -> Integral( bintrpid, nbinTrPID+1, bincvn, nbinCVN+1, bincos, nbinCosRej+1 );
573  double BkgInt = Bkg -> Integral( bintrpid, nbinTrPID+1, bincvn, nbinCVN+1, bincos, nbinCosRej+1 );
574  double DataInt = Data -> Integral( bintrpid, nbinTrPID+1, bincvn, nbinCVN+1, bincos, nbinCosRej+1 );
575 
576  // Set the inputs.
577  TotSig = SigInt;
578  TotBkg = BkgInt + DataInt;
579 
580  // Calc the FOM
581  FOM = TotSig / std::pow( TotSig + TotBkg, 0.5 );
582 
583 } // CalcFOM
584 // ---------------------------------------------------------------------------------------------------------------------
585 // ---------------------------------------------------------------------------------------------------------------------
586 void PrintMyMap( std::map< std::string, MyFOMParams > MainMap, std::string Key ) {
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
594  << std::endl;
595 } // PrintMyMap
596 // ---------------------------------------------------------------------------------------------------------------------
597 // ---------------------------------------------------------------------------------------------------------------------
598 bool PassCutPoint( double PIDScore, double PIDCut, double CVNScore, double CVNCut,
599  int Quant, int CosVeto, std::vector<int> cQuantVec, double RecoE, bool isDip ) {
600 
601  // Easiest check is on the TrPID and CVN score.
602  if ( PIDScore >= PIDCut && CVNScore >= CVNCut
603  && CosVeto == 1
604  ) {
605  // Loop through cQuantVec
606  for (size_t cq=0; cq < cQuantVec.size(); ++cq) {
607  // Check that this quant is the one I'm interested in.
608  if (Quant == cQuantVec[cq]) {
609  // Finally, if in the dip check the energy range, if not, return true.
610  if ( (RecoE > 1 && RecoE < 2) || !isDip ) return true;
611  }
612  }
613  }
614  return false;
615 } // PassCutPoint
616 // ---------------------------------------------------------------------------------------------------------------------
617 // ---------------------------------------------------------------------------------------------------------------------
618 void FillMyPlots ( double CVN, double ReMId, double CosRej, int Quant, TH1D* Hists[9], double RecoE, double Wgt ) {
619  // -- First check FOM Full
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 );
624  }
625  // -- Then check FOM Dip
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 );
630  }
631  // -- Finally, check the currently implemented cut
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 );
636  }
637  return;
638 }
639 // ---------------------------------------------------------------------------------------------------------------------
640 // ---------------------------------------------------------------------------------------------------------------------
const double FHC_POT_Scale
tree Draw("slc.nhit")
std::string sIsFD
const int HistBins
double Integral(const Spectrum &s, const double pot, int cvnregion)
bool IsFHC
const double RHC_Liv_Scale
const double MinCosRej
constexpr T pow(T x)
Definition: pow.h:75
double MyLivScale
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)
Definition: DataMCLoad.C:336
const double FHC_Liv_Scale
Defines an enumeration for prong classification.
const double MinTrPID
void PrintInfo(int Entry, int TotEnt, double IsNuMuCC, double TrueNuE, double RecoNuE, double EventWgt, double CosVeto, double CVNMuonLooPTP, double ReMId, double CVNProngMuon, double CosRej2020)
const int nbinTrPID
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)
double MyPOTScale
const double MaxCosRej
std::string sIsFHC
correl_yv Fill(-(dy[iP-1][iC-1]), hyv->GetBinContent(iP, iC))
double RecoEBins[HistBins+1]
const double MinCVN
std::string OutDir
const double RHC_POT_Scale
hmean SetLineColor(4)
const Cut Bkg
OStream cout
Definition: OStream.cxx:6
void PrintMyMap(std::map< std::string, MyFOMParams > MainMap, std::string Key)
const int nbinCVN
app
Definition: demo.py:189
T cos(T number)
Definition: d0nt_math.hpp:78
def CVN(num_classes)
Definition: models.py:37
const double MaxCVN
cosmicTree SaveAs("cosmicTree.root")
const int nbinCosRej
void ReMId()
Definition: ReMId.C:47
void NuMu2020_MakePIDPlots(bool IsFD, bool IsFHC, int QuantEnum, bool IsReMId=true, bool IsLoosePTP=true)
simulatedPeak Scale(1/simulationNormalisationFactor)
const double MaxTrPID
hh[ixs]
Definition: PlotSingle.C:6
void FillMyPlots(double CVN, double ReMId, double CosRej, int Quant, TH1D *Hists[9], double RecoE, double Wgt)
const bool bPrintInfo
enum BeamMode string