MakeCutFlowPlots.C
Go to the documentation of this file.
1 //......................................................
6 
7 #include "CAFAna/Core/Loaders.h"
8 #include "CAFAna/Core/Spectrum.h"
11 
13 
15 
16 #include "OscLib/OscCalcPMNSOpt.h"
17 #include "Utilities/rootlogon.C"
18 
19 #include "TCanvas.h"
20 #include "TColor.h"
21 #include "TFile.h"
22 #include "TGraph.h"
23 #include "TGraphAsymmErrors.h"
24 #include "TH1.h"
25 #include "TH2.h"
26 #include "TLatex.h"
27 #include "TLegend.h"
28 #include "TStyle.h"
29 #include "TSystem.h"
30 
31 #include <fstream>
32 #include <iostream>
33 #include <string>
34 
35 //......................................................
36 using namespace ana;
37 
38 //......................................................
39 double GetHistInt( TH1D* hist );
40 TH1D* GetHist ( TFile* InF, std::string InDir, int Key );
41 //......................................................
42 void SetBin( TH1D* hist, std::string lab, size_t ind, double Val );
43 //......................................................
44 void SetHistProp( TH1D* hist, int Colour, double Width, double Offset );
45 //......................................................
46 //void CornerLabel(std::string Str);
47 //......................................................
48 void MakeTextFile( std::string OutName, bool BarChart );
49 //......................................................
50 void MakeCanvasForData( std::vector<double> Data,
51  std::vector<double> Sig ,
52  std::vector<double> Wrong,
53  std::vector<double> Beam ,
54  std::vector<double> Cosm ,
55  std::vector<std::string> CutNames, size_t NCuts, bool NuMu );
56 //......................................................
57 void MakeCanvasForSim ( std::vector<double> Sig ,
58  std::vector<double> Wrong,
59  std::vector<double> Beam ,
60  std::vector<double> Cosm ,
61  std::vector<std::string> CutNames, size_t NCuts, bool NuMu );
62 //......................................................
63 void MakeNueTable( std::vector<double> Sig , std::vector<double> WrSi, std::vector<double> Beam,
64  std::vector<double> NuMu, std::vector<double> NuTau, std::vector<double> NC,
65  std::vector<double> Data, std::vector<double> Cosm, std::vector<std::string> CutNames, bool ShowData );
66 //......................................................
67 void MakeNuMuTable( std::vector<double> Sig , std::vector<double> WrSi,
68  std::vector<double> Nue , std::vector<double> NuTau, std::vector<double> NC,
69  std::vector<double> Data, std::vector<double> Cosm, std::vector<std::string> CutNames, bool ShowData );
70 //......................................................
72  // Some large scale changes.
73  if ( Str.find("Nue_Mask") != std::string::npos ) Str = "Ana. Mask";
74  if ( Str.find("Nue_DQ" ) != std::string::npos ) Str = "Data Qual.";
75  if ( Str.find("Nue_Cont") != std::string::npos ) Str = "Containment";
76  if ( Str.find("Nue_En" ) != std::string::npos ) Str = "Energy Cut";
77  if ( Str.find("Nue_CRej") != std::string::npos ) Str = "Cosmic Rej.";
78  //if ( Str.find("Nue_NSlc") != std::string::npos ) Str = "Nearest Slc";
79  if ( Str.find("Nue_NSlc") != std::string::npos ) Str = "Full PID";
80 
81  if ( Str.find("Periph_Rev" ) != std::string::npos ) Str = "Reverse Core";
82  if ( Str.find("Periph_Base") != std::string::npos ) Str = "Basic Cuts";
83  //if ( Str.find("Periph_NSlc") != std::string::npos ) Str = "Nearest Slc";
84  if ( Str.find("Periph_NSlc") != std::string::npos ) Str = "Full PID";
85 
86  if ( Str.find("NuMu_CosmicRej") != std::string::npos ) Str = "Cosmic Rej.";
87  if ( Str.find("NuMu_FullPID" ) != std::string::npos ) Str = "Full PID";
88 
89  // Some general changes.
90  if ( Str.find("NuMu_" ) != std::string::npos ) Str.erase( Str.find("NuMu_" ), 5 );
91  else if ( Str.find("Nue_" ) != std::string::npos ) Str.erase( Str.find("Nue_" ), 4 );
92  else if ( Str.find("Periph_") != std::string::npos ) Str.erase( Str.find("Periph_"), 7 );
93  while (Str.find("_") != std::string::npos) Str.replace(Str.find("_"),1," ");
94  return Str;
95 }
96 //......................................................
97 size_t TrimVector( std::vector<double> &Sig , std::vector<double> &WS , std::vector<double> &Bkg,
98  std::vector<double> &Cosm, std::vector<double> &Data, std::vector<std::string> &Names);
99 //......................................................
100 
101 double POTNorm = 0;
102 double LivNorm = 0;
103 
104 bool isFHC = true;
106 
108 //......................................................
109 // Samp: "full" "periods1235" "periods910"
110 // BFPo: "bf2019" "IHUO2019" "NHLO2019"
111 void MakeCutFlowPlots( bool isfhc, std::string Samp, std::string BFPoint="bf2020" ) {
112 
113  // Set my FHC vars.
114  isFHC = isfhc;
115  sFHC = isFHC? "fhc":"rhc";
116  sBFPoint = BFPoint;
117 
118  // Set my POT and Livetime.
119  if (isFHC) {
122  } else {
125  }
126 
127  // Where are my input files?
128  std::string InBase = "/nova/ana/users/karlwarb/Analysis2020/CutFlow/200615/";
129  std::string MontLoc = InBase + "MonteCarloCutFlow_"+sFHC+".root";
130  std::string CosmLoc = InBase + "CosmicCutFlow_" +sFHC+".root";
131  std::string DataLoc = InBase + "NuMIDataCutFlow_" +sFHC+"_"+Samp+".root";
132 
133  std::cerr << "Loading;"
134  << "\n\t Monte Carlo - " << MontLoc
135  << "\n\t Cosmic File - " << CosmLoc
136  << "\n\t Data Files - " << DataLoc
137  << std::endl;
138 
139  TFile *MontFile = TFile::Open( MontLoc.c_str() );
140  TFile *CosmFile = TFile::Open( CosmLoc.c_str() );
141  TFile *DataFile = TFile::Open( DataLoc.c_str() );
142 
143  // --- Make my cut Vectors.
144  std::vector<std::string> NuMu_CutNames; std::vector<Cut> NuMu_TieredCuts;
145  std::vector<std::string> Nue_CutNames; std::vector<Cut> Nue_TieredCuts;
146  MakeCutFlowVecs( NuMu_CutNames, NuMu_TieredCuts, Nue_CutNames, Nue_TieredCuts );
147  size_t NuMuCuts = NuMu_TieredCuts.size();
148  size_t NueCuts = Nue_TieredCuts .size();
149 
150  // --- Check if I need to rescale my POT and Livetime as I'm looking at a subsample.
151  std::cout << "\n\t POTNorm: " << POTNorm << ", LivNorm: " << LivNorm << std::endl;
152  if (isfhc && Samp != "full" ) {
153  if (Samp == "periods1235") {
156  } else if (Samp == "periods910") {
159  } else {
160  std::cout << "\n You passed Samp: " << Samp << ". I don't recognise that. Exiting." << std::endl;
161  return;
162  }
163  }
164  std::cout << "\n\t After fix, POTNorm: " << POTNorm << ", LivNorm: " << LivNorm << std::endl;
165 
166  // --- Define my histograms.
167  // NuMu double vec
168  std::vector<TH1D* > hBeam_NuMu, hCosm_NuMu, hNuMuCC, hNuMuWS, hNuMuNue, hNuMuNuTau, hNuMuNC, hNotNuMuCC;
169  std::vector<double> nBeam_NuMu, nCosm_NuMu, nNuMuCC, nNuMuWS, nNuMuNue, nNuMuNuTau, nNuMuNC, nNotNuMuCC;
170  // Nue double vec
171  std::vector<TH1D* > hBeam_Nue , hCosm_Nue , hNueCC , hNueWS , hNueNuMu, hNueNuTau , hNueNC , hNotNueCC , hNueBeam;
172  std::vector<double> nBeam_Nue , nCosm_Nue , nNueCC , nNueWS , nNueNuMu, nNueNuTau , nNueNC , nNotNueCC , nNueBeam;
173 
174  // --- Loop through my vector of cut names, and load the histograms.
175  // 0 - Data -> Norm by POT
176  // 1 - Cosm -> Norm by Livetime
177 
178  // 2 - NuMu Signal
179  // 3 - NuMu Wrong Sign
180  // 4 - Nue backgrounds
181  // 5 - NuTau backgrounds
182  // 6 - NC backgrounds
183 
184  // 7 - Nue Signal
185  // 8 - Nue Wrong Sign
186  // 9 - Beam Nues
187  // 10 - NuMu backgrounds
188  // 11 - NuTau backgrounds
189  // 12- NC backgrounds.
190 
191  for (size_t cc=0; cc<NuMuCuts; ++cc) {
192  // --- Data strings.
193  const std::string kThisName = "NuMuCCEnergy_"+NuMu_CutNames[cc];
194  // --- Load the Data
195  hBeam_NuMu.push_back( GetHist( DataFile, kThisName, 0 ) ); nBeam_NuMu.push_back( GetHistInt( hBeam_NuMu.back() ) );
196  hCosm_NuMu.push_back( GetHist( CosmFile, kThisName, 1 ) ); nCosm_NuMu.push_back( GetHistInt( hCosm_NuMu.back() ) );
197  // --- Load the Monte Carlo
198  hNuMuCC .push_back( GetHist( MontFile, kThisName, 2 ) ); nNuMuCC .push_back( GetHistInt( hNuMuCC .back() ) );
199  hNuMuWS .push_back( GetHist( MontFile, kThisName, 3 ) ); nNuMuWS .push_back( GetHistInt( hNuMuWS .back() ) );
200  hNuMuNue .push_back( GetHist( MontFile, kThisName, 4 ) ); nNuMuNue .push_back( GetHistInt( hNuMuNue .back() ) );
201  hNuMuNuTau.push_back( GetHist( MontFile, kThisName, 5 ) ); nNuMuNuTau.push_back( GetHistInt( hNuMuNuTau.back() ) );
202  hNuMuNC .push_back( GetHist( MontFile, kThisName, 6 ) ); nNuMuNC .push_back( GetHistInt( hNuMuNC .back() ) );
203  // --- Make the NotNutNuMuCC hist.
204  TH1D* Temp = (TH1D*)hNuMuNue.back()->Clone(); Temp->Add( hNuMuNuTau.back() ); Temp->Add( hNuMuNC.back() );
205  hNotNuMuCC.push_back( Temp ); nNotNuMuCC.push_back( GetHistInt( hNotNuMuCC.back() ) );
206  }
207  for (size_t cc=0; cc<NueCuts; ++cc) {
208  // --- Data strings.
209  const std::string kThisName = "NueCCEnergy_" +Nue_CutNames [cc];
210  // --- Load the Data
211  hBeam_Nue.push_back( GetHist( DataFile, kThisName, 0 ) ); nBeam_Nue.push_back( GetHistInt( hBeam_Nue.back() ) );
212  hCosm_Nue.push_back( GetHist( CosmFile, kThisName, 1 ) ); nCosm_Nue.push_back( GetHistInt( hCosm_Nue.back() ) );
213  // --- Load the Monte Carlo
214  hNueCC .push_back( GetHist( MontFile, kThisName, 7 ) ); nNueCC .push_back( GetHistInt( hNueCC .back() ) );
215  hNueWS .push_back( GetHist( MontFile, kThisName, 8 ) ); nNueWS .push_back( GetHistInt( hNueWS .back() ) );
216  hNueBeam .push_back( GetHist( MontFile, kThisName, 9 ) ); nNueBeam .push_back( GetHistInt( hNueBeam .back() ) );
217  hNueNuMu .push_back( GetHist( MontFile, kThisName, 10) ); nNueNuMu .push_back( GetHistInt( hNueNuMu .back() ) );
218  hNueNuTau.push_back( GetHist( MontFile, kThisName, 11) ); nNueNuTau.push_back( GetHistInt( hNueNuTau.back() ) );
219  hNueNC .push_back( GetHist( MontFile, kThisName, 12) ); nNueNC .push_back( GetHistInt( hNueNC .back() ) );
220  // --- Make the NotNutNuMuCC hist.
221  TH1D* Temp = (TH1D*)hNueBeam.back()->Clone(); Temp->Add( hNueNuMu.back() ); Temp->Add( hNueNuTau.back() ); Temp->Add( hNueNC.back() );
222  hNotNueCC.push_back( Temp ); nNotNueCC.push_back( GetHistInt( hNotNueCC.back() ) );
223  }
224 
225  // --- Now lets start drawing things!!!
226  std::cerr << "\n\nNow lets start plotting things..." << std::endl;
227 
228  gSystem -> MakeDirectory( "Plots" );
229  std::string OutName = "Plots/CutFlowPlots_"+sFHC+"_Ana2020.root";
230  TFile *OutFile = new TFile(OutName.c_str(), "RECREATE");
231  OutFile -> cd();
232 
233  // --- First, trim the names.
234  for (size_t cc=0; cc<NuMuCuts; ++cc) { CutStr( NuMu_CutNames[cc] ); }
235  for (size_t cc=0; cc<NueCuts ; ++cc) { CutStr( Nue_CutNames [cc] ); }
236 
237  // Make Nue Tables.
238  MakeNueTable( nNueCC, nNueWS, nNueBeam, nNueNuMu, nNueNuTau, nNueNC, nBeam_Nue, nCosm_Nue, Nue_CutNames, true );
239  MakeNueTable( nNueCC, nNueWS, nNueBeam, nNueNuMu, nNueNuTau, nNueNC, nBeam_Nue, nCosm_Nue, Nue_CutNames, false );
240 
241  // Make Numu Tables.
242  MakeNuMuTable( nNuMuCC, nNuMuWS, nNuMuNue, nNuMuNuTau, nNuMuNC, nBeam_NuMu, nCosm_NuMu, NuMu_CutNames, true );
243  MakeNuMuTable( nNuMuCC, nNuMuWS, nNuMuNue, nNuMuNuTau, nNuMuNC, nBeam_NuMu, nCosm_NuMu, NuMu_CutNames, false );
244 
245  // --- Make a canvas for FD Data, TotalPred, BeamBkg, Cosmics
246  //MakeCanvasForData( nBeam_NuMu, nNuMuCC, nNuMuWS, nNotNuMuCC, nCosm_NuMu, NuMu_CutNames, NuMuCuts, true );
247  //MakeCanvasForData( nBeam_Nue , nNueCC , nNueWS , nNotNueCC , nCosm_Nue , Nue_CutNames , NueCuts , false );
248 
249  // --- Make a canvas for event rate without data
250  //MakeCanvasForSim( nNuMuCC, nNuMuWS, nNotNuMuCC, nCosm_NuMu, NuMu_CutNames, NuMuCuts, true );
251  //MakeCanvasForSim( nNueCC , nNueWS , nNotNueCC , nCosm_Nue , Nue_CutNames , NueCuts , false );
252 
253  return;
254 }
255 //......................................................
256 double GetHistInt( TH1D* hist ) {
257  return hist->Integral();
258 }
259 //......................................................
260 TH1D* GetHist( TFile* InF, std::string InDir, int Key ) {
261  // 0 - Data -> Norm by POT
262  // 1 - Cosm -> Norm by Livetime
263 
264  // 2 - NuMu Signal
265  // 3 - NuMu Wrong Sign
266  // 4 - Nue backgrounds
267  // 5 - NuTau backgrounds
268  // 6 - NC backgrounds
269 
270  // 7 - Nue Signal
271  // 8 - Nue Wrong Sign
272  // 9 - Beam Nues
273  // 10 - NuMu backgrounds
274  // 11 - NuTau backgrounds
275  // 12- NC backgrounds.
276 
277  // Loading a data hist, so Spectrum.
278  if ( Key == 0 || Key == 1 ) {
279  std::unique_ptr<Spectrum> Temp = Spectrum::LoadFrom( InF, InDir );
280  if ( Key == 0 ) return Temp -> ToTH1( POTNorm );
281  else return Temp -> ToTH1( LivNorm, kLivetime );
282  }
283  else {
284  // Loading Monte Carlo, so PredNoExtrap.
285  std::unique_ptr<PredictionNoExtrap> TempSpec = PredictionNoExtrap::LoadFrom( InF, InDir ) ;
286 
287  // --- Set the oscillation point
289  if (sBFPoint == "2019" ) { // 2019 Best Fit
290  calc -> SetdCP ( 2*M_PI );
291  calc -> SetTh23 ( asin(sqrt(0.565)) );
292  calc -> SetDmsq32( 2.48e-3 );
293  } else {
294  calc -> SetdCP ( 0.82*M_PI );
295  calc -> SetTh23 ( asin(sqrt(0.568)) );
296  calc -> SetDmsq32( 2.41e-3 );
297  }
298 
299  // --- What neutrino sign to use
300  Sign::Sign_t MySign = Sign::kNu;
301  Sign::Sign_t MyWrongSign = Sign::kAntiNu;
302  if (!isFHC) {
303  MySign = Sign::kAntiNu;
304  MyWrongSign = Sign::kNu;
305  }
306 
307  // --- Now NuMu's.
308  if ( Key == 2 ) { // NuMu Signal
309  return TempSpec->PredictComponent(calc, Flavors::kNuMuToNuMu , Current::kCC, MySign ).ToTH1( POTNorm );
310  } else if ( Key == 3 ) { // NuMu Wrong Sign
311  return TempSpec->PredictComponent(calc, Flavors::kNuMuToNuMu , Current::kCC, MyWrongSign ).ToTH1( POTNorm );
312  } else if ( Key == 4 ) { // NuMu Nue Backgrounds
313  return TempSpec->PredictComponent(calc, Flavors::kAllNuE , Current::kCC, Sign::kBoth ).ToTH1( POTNorm );
314  } else if ( Key == 5 ) { // NuMu NuTau Backgrounds
315  return TempSpec->PredictComponent(calc, Flavors::kNuMuToNuTau, Current::kCC, Sign::kBoth ).ToTH1( POTNorm );
316  } else if ( Key == 6 ) { // NuMu NC Backgrounds
317  return TempSpec->PredictComponent(calc, Flavors::kAll , Current::kNC, Sign::kBoth ).ToTH1( POTNorm );
318  }
319  // ---- Now Nue's.
320  else if ( Key == 7 ) { // Nue Signal
321  return TempSpec->PredictComponent(calc, Flavors::kNuMuToNuE , Current::kCC, MySign ).ToTH1( POTNorm );
322  } else if ( Key == 8 ) { // Nue Wrong Sign
323  return TempSpec->PredictComponent(calc, Flavors::kNuMuToNuE , Current::kCC, MyWrongSign ).ToTH1( POTNorm );
324  } else if ( Key == 9 ) { // Beam Nues
325  return TempSpec->PredictComponent(calc, Flavors::kNuEToNuE , Current::kCC, Sign::kBoth ).ToTH1( POTNorm );
326  } else if ( Key == 10 ) { // Nue NuMu Backgrounds
327  return TempSpec->PredictComponent(calc, Flavors::kAllNuMu , Current::kCC, Sign::kBoth ).ToTH1( POTNorm );
328  } else if ( Key == 11 ) { // Nue NuTau Backgrounds
329  return TempSpec->PredictComponent(calc, Flavors::kNuMuToNuTau, Current::kCC, Sign::kBoth ).ToTH1( POTNorm );
330  } else if ( Key == 12 ) { // Nue NC Backgrounds
331  return TempSpec->PredictComponent(calc, Flavors::kAll , Current::kNC, Sign::kBoth ).ToTH1( POTNorm );
332  }
333 
334  // Shouldn't reach this far down...
335  std::cout << "You passed GetHistVectors:Key " << Key << ", which doesn't fit. Aborting." << std::endl;
336  abort();
337  }
338 }
339 //......................................................
340 void MakeTextFile( std::string OutName, bool BarChart ) {
341  // Determine output name, and caption.
342  std::string FNa = "Plots/"+OutName+".txt";
343  while (OutName.find("_") != std::string::npos) OutName.replace(OutName.find("_"),1," ");
344 
345  std::string Cap = "A bar chart showing the number of events changes as subsequent cuts are applied, ";
346  if (OutName.find("Data") != std::string::npos) {
347  Cap += "with NuMi data shown, along with estimates from MC and Cosmic sidebands.";
348  } else {
349  Cap += "for expected data rates, taken from MC and Cosmic sidebands.";
350  }
351 
352  std::cout << "\nFile name: " << FNa << "\n\tCaption: " << Cap << std::endl;
353 
354  // Write to file.
355  std::ofstream TxtOut ( FNa.c_str(), std::ofstream::out );
356  TxtOut << Cap;
357  TxtOut.close();
358  // Done.
359  return;
360 }
361 //......................................................
362 void MakeNueTable( std::vector<double> Sig , std::vector<double> WrSi, std::vector<double> Beam,
363  std::vector<double> NuMu, std::vector<double> NuTau, std::vector<double> NC,
364  std::vector<double> Data, std::vector<double> Cosm, std::vector<std::string> CutNames, bool ShowData ) {
365  std::string sSig = "$\\nu_e$" , sWrSi = "$\\overline{\\nu_e}$";
366  if (!isFHC) { sSig = "$\\overline{\\nu_e}$"; sWrSi = "$\\nu_e$"; }
367  if (ShowData) {
368  std::cout << "\n..........................................................\n\t Table for Nue, including the Data.\n\n"
369  << "\\begin{table} \n \\centering \n "
370  << "\\caption{Table showing the "+sSig+" analysis selection broken into sequential cuts, and interaction channels, including the selected NuMI data.} "
371  //<< "\n \\scriptsize \n \\begin{tabular}{@{}c|c|cccc cccc|cc|c@{}} \n \\hline \n"
372  //<< " Cut Level & Data & "<<sSig<<" & Tot Bkg & "<<sWrSi<<" & Beam $\\nu_e$ & $\\nu_\\mu$ CC & $\\nu_\\tau$ CC & NC & Cosmics & Effic & Purity & FOM \\\\ \\hline"
373  << "\n \\scriptsize \n \\begin{tabular}{@{}c|c|ccccc|cc@{}} \n \\hline \n"
374  << " Cut Level & Data & "<<sSig<<" & "<<sWrSi<<" & Total Beam Bkg & NC & Cosmics & Effic & Purity \\\\ \\hline"
375  << std::endl;
376  } else {
377  std::cout << "\n..........................................................\n\t Table for Nue, without the Data.\n\n"
378  << "\\begin{table} \n \\centering \n "
379  << "\\caption{Table showing the "+sSig+" analysis selection broken into sequential cuts, and interaction channels.} "
380  //<< "\n \\scriptsize \n \\begin{tabular}{@{}c|cccc cccc|cc|c@{}} \n \\hline \n"
381  //<< " Cut Level & "<<sSig<<" & Tot Bkg & "<<sWrSi<<" & Beam $\\nu_e$ & $\\nu_\\mu$ CC & $\\nu_\\tau$ CC & NC & Cosmics & Effic & Purity & FOM \\\\ \\hline"
382  << "\n \\scriptsize \n \\begin{tabular}{@{}c|ccccc|cc@{}} \n \\hline \n"
383  << " Cut Level & "<<sSig<<" & "<<sWrSi<<" & Total Beam Bkg & NC & Cosmics & Effic & Purity \\\\ \\hline"
384  << std::endl;
385  }
386 
387  for (size_t cc=0; cc<Sig.size(); ++cc) {
388  double TotBkg = WrSi[cc] + Beam[cc] + NuMu[cc] + NuTau[cc] + NC[cc] + Cosm[cc];
389 
390  std::cout.precision(3);
391  std::cout << "\n " << CutNames[cc];
392  if (ShowData) std::cout << " & " << Data[cc];
393  std::cout << " & " << Sig [cc]
394  //<< " & " << TotBkg
395  << " & " << WrSi [cc]
396  //<< " & " << Beam [cc] << " & " << NuMu [cc] << " & " << NuTau[cc]
397  << " & " << Beam [cc] + NuMu [cc] + NuTau[cc]
398  << " & " << NC [cc]
399  << " & " << Cosm [cc];
400  //if (cc==0) { std::cout<< " & & & \\\\"; }
401  if (cc==0) { std::cout<< " & & \\\\"; }
402  else {
403  std::cout << " & " << (Sig[cc]) / (Sig[0])
404  << " & " << (Sig[cc]) / (Sig[cc] + TotBkg)
405  //<< " & " << (Sig[cc]) / pow( Sig[cc]+TotBkg, 0.5 )
406  << " \\\\ ";
407  }
408  if (CutNames[cc] == "Full PID" && cc!=Sig.size()-1) std::cout << "\\hline";
409  }
410  std::cout << "\n\n \\hline \n \\end{tabular} \n \\label{} \n\\end{table}"
411  << "\n\t The analysis efficiency is: " << (Sig[8]+Sig[12]) / Sig[0 ]
412  << "\n\t The analysis purity is : " << (Sig[8]+Sig[12]) / ( Sig[8 ]+WrSi[8 ]+Beam[8 ]+NuMu[8 ]+NuTau[8 ]+NC[8 ]+Cosm[8 ]+
413  Sig[12]+WrSi[12]+Beam[12]+NuMu[12]+NuTau[12]+NC[12]+Cosm[12])
414  << "\n.........................................................." << std::endl;
415 }
416 //......................................................
417 void MakeNuMuTable( std::vector<double> Sig , std::vector<double> WrSi,
418  std::vector<double> Nue , std::vector<double> NuTau, std::vector<double> NC,
419  std::vector<double> Data, std::vector<double> Cosm, std::vector<std::string> CutNames, bool ShowData ) {
420  std::string sSig = "$\\nu_\\mu$" , sWrSi = "$\\overline{\\nu_\\mu}$";
421  if (!isFHC) { sSig = "$\\overline{\\nu_\\mu}$"; sWrSi = "$\\nu_\\mu$"; }
422  if (ShowData) {
423  std::cout << "\n..........................................................\n\t Table for Numu, including the Data.\n\n"
424  << "\\begin{table} \n \\centering \n "
425  << "\\caption{Table showing the "+sSig+" analysis selection broken into sequential cuts, and interaction channels, including the selected NuMI data.} "
426  //<< "\n \\scriptsize \n \\begin{tabular}{@{}c|c|ccc cccc|cc|c@{}} \n \\hline \n"
427  //<< " Cut Level & Data & "<<sSig<<" & "<<sWrSi<<" & Tot Bkg & $\\nu_e$ CC & $\\nu_\\tau$ CC & NC & Cosmics & Effic & Purity & FOM \\\\ \\hline"
428  << "\n \\scriptsize \n \\begin{tabular}{@{}c|c|ccccc|cc@{}} \n \\hline \n"
429  << " Cut Level & Data & "<<sSig<<" & "<<sWrSi<<" & Total Beam Bkg & NC & Cosmics & Effic & Purity \\\\ \\hline"
430  << std::endl;
431  } else {
432  std::cout << "\n..........................................................\n\t Table for Numu, without the Data.\n\n"
433  << "\\begin{table} \n \\centering \n "
434  << "\\caption{Table showing the "+sSig+" analysis selection broken into sequential cuts, and interaction channels.} "
435  //<< "\n \\scriptsize \n \\begin{tabular}{@{}c|ccc cccc|cc|c@{}} \n \\hline \n"
436  //<< " Cut Level & "<<sSig<<" & "<<sWrSi<<" & Tot Bkg & $\\nu_e$ CC & $\\nu_\\tau$ CC & NC & Cosmics & Effic & Purity & FOM \\\\ \\hline"
437  << "\n \\scriptsize \n \\begin{tabular}{@{}c|ccccc|cc@{}} \n \\hline \n"
438  << " Cut Level & "<<sSig<<" & "<<sWrSi<<" & Total Beam Bkg & NC & Cosmics & Effic & Purity \\\\ \\hline"
439  << std::endl;
440  }
441 
442  for (size_t cc=0; cc<Sig.size(); ++cc) {
443  double TotBkg = Nue[cc] + NuTau[cc] + NC[cc] + Cosm[cc];
444 
445  std::cout.precision(3);
446  std::cout << "\n " << CutNames[cc];
447  if (ShowData) std::cout << " & " << Data[cc];
448  std::cout << " & " << Sig [cc]
449  << " & " << WrSi [cc]
450  //<< " & " << TotBkg << " & " << Nue [cc] << " & " << NuTau[cc]
451  << " & " << Nue [cc] + NuTau[cc]
452  << " & " << NC [cc]
453  << " & " << Cosm [cc];
454  //if (cc==0) { std::cout<< " & & & \\\\"; }
455  if (cc==0) { std::cout<< " & & \\\\"; }
456  else {
457  std::cout << " & " << (Sig[cc]) / (Sig[0])
458  << " & " << (Sig[cc]) / (Sig[cc] + TotBkg)
459  //<< " & " << (Sig[cc]) / pow( Sig[cc]+TotBkg, 0.5 )
460  << " \\\\ ";
461  }
462  }
463  std::cout << "\n\n \\hline \n \\end{tabular} \n \\label{} \n\\end{table} \n.........................................................." << std::endl;
464 }
465 //......................................................
466 void SetBin( TH1D* hist, std::string lab, size_t ind, double Val ) {
467  hist -> SetBinContent( ind, Val );
468  while (lab.find("_") != std::string::npos) lab.replace(lab.find("_"),1," ");
469  hist -> GetXaxis() -> SetBinLabel( ind, lab.c_str() );
470  return;
471 }
472 //......................................................
473 void SetHistProp( TH1D* hist, int Colour, double Width, double Offset ) {
474  hist -> SetMinimum ( 0.1 );
475  hist -> SetFillColor( Colour );
476  hist -> SetBarWidth ( Width );
477  hist -> SetBarOffset( Offset );
478  hist -> SetFillStyle( 1001 );
479  hist -> SetStats(0);
480  for (int bi=0; bi<hist->GetNbinsY(); ++bi) hist->SetBinContent( bi, 0.1 );
481  return;
482 }
483 //......................................................
484 void MakeCanvasForData( std::vector<double> Data,
485  std::vector<double> Sig ,
486  std::vector<double> Wrong,
487  std::vector<double> Beam ,
488  std::vector<double> Cosm ,
489  std::vector<std::string> CutNames, size_t NCuts, bool NuMu ) {
490  //gStyle->SetHistMinimumZero();
491  std::string sNuMu = "NuMu";
492  if (!NuMu) sNuMu = "Nue";
493 
494  std::string CanNa = "DataCutFlow_"+sNuMu+"_"+sFHC;
495  std::string PredNa = "DataPred_" +sNuMu;
496  std::string SigNa = "DataSig_" +sNuMu;
497  std::string WrSNa = "DataWrS_" +sNuMu;
498  std::string BeamNa = "DataBeam_"+sNuMu;
499  std::string CosmNa = "DataCosm_"+sNuMu;
500  std::string DataNa = "DataNuMI_"+sNuMu;
501 
502  // --- Declare a canvas.
503  TCanvas *DataCan = new TCanvas( CanNa.c_str(), "Cut Efficiency Canvas for Data vs MC" );
504  DataCan -> SetLogx();
505  DataCan -> SetGridx();
506  DataCan -> SetLeftMargin(0.15);
507  DataCan -> SetRightMargin(0.05);
508 
509  double HistWidth = 0.14;
510 
511  // First Total Predicted flux.
512  TH1D *hTot = new TH1D( PredNa.c_str(), "", NCuts, 0, NCuts );
513  if (isFHC) hTot -> GetYaxis() -> SetRangeUser( 0.1, 1e7 );
514  else hTot -> GetYaxis() -> SetRangeUser( 0.1, 1e7 );
515  SetHistProp( hTot, kRed, HistWidth, 0.78 );
516  for(size_t i=1; i<=NCuts; i++) {
517  double TotBin = Sig[ NCuts-i ] + Cosm[ NCuts-i ];
518  SetBin( hTot, CutNames[NCuts-i], i, TotBin );
519  }
520  // Draw.
521  hTot -> Draw("hist");
522  // Set axis range etc.
523  hTot -> GetXaxis() -> SetLabelSize( 0.055 );
524  hTot -> GetYaxis() -> SetTitle("Events");
525  hTot -> GetYaxis() -> CenterTitle();
526 
527  // Next, Signal events.
528  TH1D *hSig = new TH1D( SigNa.c_str(), "", NCuts, 0, NCuts );
529  SetHistProp( hSig, kViolet+1, HistWidth, 0.64 );
530  for(size_t i=1; i<=NCuts; i++) {
531  SetBin( hSig, CutNames[NCuts-i], i, Sig[ NCuts-i ] );
532  }
533  hSig -> Draw("hist same");
534 
535  // Now Wrong sign
536  TH1D *hWrS = new TH1D(WrSNa.c_str(),"",NCuts, 0, NCuts);
537  SetHistProp( hWrS, kGreen+2, HistWidth, 0.50 );
538  for(size_t i=1; i<=NCuts; i++) {
539  SetBin( hWrS, CutNames[NCuts-i], i, Wrong[ NCuts-i ]);
540  }
541  hWrS -> Draw("hist same");
542 
543  // Now Beam background
544  TH1D *hBea = new TH1D(BeamNa.c_str(),"",NCuts, 0, NCuts);
545  SetHistProp( hBea, kGray+1, HistWidth, 0.36 );
546  for(size_t i=1; i<=NCuts; i++) {
547  SetBin( hBea, CutNames[NCuts-i], i, Beam[ NCuts-i ]);
548  }
549  hBea -> Draw("hist same");
550 
551  // Now Cosmic background
552  TH1D *hCos = new TH1D(CosmNa.c_str(),"",NCuts, 0, NCuts);
553  SetHistProp( hCos, kAzure+1, HistWidth, 0.22 );
554  for(size_t i=1; i<=NCuts; i++) {
555  SetBin( hCos, CutNames[NCuts-i], i, Cosm[ NCuts-i ]);
556  }
557  hCos -> Draw("hist same");
558 
559  // Last Data
560  TH1D *hDat = new TH1D(DataNa.c_str(),"",NCuts, 0, NCuts);
561  SetHistProp( hDat, kBlack, HistWidth, 0.08 );
562  for(size_t i=1; i<=NCuts; i++) {
563  SetBin( hDat, CutNames[NCuts-i], i, Data[ NCuts-i ]);
564  }
565  hDat -> Draw("hist same");
566 
567  // Add a legend to the plot.
568  double ymin = 0.15, ymax = 0.39;
569  if (!NuMu){ymin = 0.32; ymax = 0.56; }
570  TLegend *l = new TLegend(0.6, ymin, 0.9, ymax);
571  l->AddEntry(hTot,"Predicted Num. Events","f");
572  if (NuMu) {
573  if (isFHC) l->AddEntry(hSig,"#nu_{#mu} Signal" ,"f");
574  else l->AddEntry(hSig,"#bar{#nu_{#mu}} Signal" ,"f");
575  if (isFHC) l->AddEntry(hWrS,"Wrong Sign: #bar{#nu_{#mu}} CC","f");
576  else l->AddEntry(hWrS,"Wrong Sign: #nu_{#mu} CC" ,"f");
577  } else {
578  if (isFHC) l->AddEntry(hSig,"#nu_{e} Signal" ,"f");
579  else l->AddEntry(hSig,"#bar{#nu_{e}} Signal" ,"f");
580  if (isFHC) l->AddEntry(hWrS,"Wrong Sign: #bar{#nu_{e}} CC" ,"f");
581  else l->AddEntry(hWrS,"Wrong Sign: #nu_{e} CC" ,"f");
582  }
583  l->AddEntry(hBea,"Beam Background" ,"f");
584  l->AddEntry(hCos,"Cosmic Background" ,"f");
585  l->AddEntry(hDat,"Recorded data" ,"f");
586  l->Draw();
587 
588  // Draw Preliminary()
589  Preliminary();
590  if (isFHC) CornerLabel("Neutrino beam");
591  else CornerLabel("Antineutrino beam");
592 
593  // And redraw axis
594  hTot -> Draw("axis same");
595  DataCan -> Update();
596 
597  // Finally, save the canvas.
598  DataCan -> SaveAs( TString("Plots/")+TString(CanNa)+TString(".pdf" ) );
599  DataCan -> Write ( CanNa.c_str() );
600  MakeTextFile ( CanNa , true );
601 
602  return;
603 }
604 //......................................................
605 void MakeCanvasForSim ( std::vector<double> Sig ,
606  std::vector<double> Wrong,
607  std::vector<double> Beam ,
608  std::vector<double> Cosm ,
609  std::vector<std::string> CutNames, size_t NCuts, bool NuMu ) {
610 
611  std::string sNuMu = "NuMu";
612  if (!NuMu) sNuMu = "Nue";
613 
614  std::string CanNa = "SimCutFlow_"+sNuMu+"_"+sFHC;
615  std::string SigNa = "Sig_" +sNuMu;
616  std::string WrSNa = "WrS_" +sNuMu;
617  std::string BeamNa = "Beam_"+sNuMu;
618  std::string CosmNa = "Cosm_"+sNuMu;
619 
620  // --- Declare a canvas.
621  TCanvas *SimCan = new TCanvas( CanNa.c_str(), "Cut Efficiency Canvas for MC" );
622  SimCan -> SetLogx();
623  SimCan -> SetGridx();
624  SimCan -> SetLeftMargin(0.15);
625  SimCan -> SetRightMargin(0.05);
626 
627  double HistWidth = 0.2;
628 
629  // Signal
630  TH1D *hSig = new TH1D( SigNa.c_str(), "", NCuts, 0, NCuts );
631  SetHistProp( hSig, kViolet+1, HistWidth, 0.7 );
632  for(size_t i=1; i<=NCuts; i++) {
633  SetBin( hSig, CutNames[NCuts-i], i, Sig[ NCuts-i ]);
634  }
635  // Draw.
636  hSig -> Draw("hist");
637  // Set axis range etc.
638  hSig -> GetXaxis() -> SetLabelSize( 0.055 );
639  hSig -> GetYaxis() -> SetTitle("Events");
640  hSig -> GetYaxis() -> CenterTitle();
641  if (isFHC) hSig -> GetYaxis() -> SetRangeUser( 0.1, 1e7 );
642  else hSig -> GetYaxis() -> SetRangeUser( 0.1, 1e7 );
643 
644  // Now Wrong sign
645  TH1D *hWrS = new TH1D( WrSNa.c_str(),"",NCuts, 0, NCuts);
646  SetHistProp( hWrS, kGreen+2, HistWidth, 0.5 );
647  for(size_t i=1; i<=NCuts; i++) {
648  SetBin( hWrS, CutNames[NCuts-i], i, Wrong[ NCuts-i ]);
649  }
650  hWrS -> Draw("hist same");
651 
652  // Now Beam background
653  TH1D *hBeam = new TH1D( BeamNa.c_str(), "", NCuts, 0, NCuts );
654  SetHistProp( hBeam, kGray+1, HistWidth, 0.3 );
655  for(size_t i=1; i<=NCuts; i++) {
656  SetBin( hBeam, CutNames[NCuts-i], i, Beam[ NCuts-i ]);
657  }
658  hBeam -> Draw("hist same");
659 
660  // Now Cosmic background
661  TH1D *hCosm = new TH1D( CosmNa.c_str(), "", NCuts, 0, NCuts );
662  SetHistProp( hCosm, kAzure+1, HistWidth, 0.1 );
663  for(size_t i=1; i<=NCuts; i++) {
664  SetBin( hCosm, CutNames[NCuts-i], i, Cosm[ NCuts-i ]);
665  }
666  hCosm -> Draw("hist same");
667 
668  // Add a legend to the plot.
669  double ymin = 0.15, ymax = 0.35;
670  if (!NuMu){ymin = 0.30; ymax = 0.50; }
671  TLegend *l = new TLegend(0.6, ymin, 0.9, ymax);
672  if (NuMu) {
673  if (isFHC) l->AddEntry(hSig,"#nu_{#mu} Signal" ,"f");
674  else l->AddEntry(hSig,"#bar{#nu_{#mu}} Signal" ,"f");
675  if (isFHC) l->AddEntry(hWrS,"Wrong Sign: #bar{#nu_{#mu}} CC","f");
676  else l->AddEntry(hWrS,"Wrong Sign: #nu_{#mu} CC" ,"f");
677  } else {
678  if (isFHC) l->AddEntry(hSig,"#nu_{e} Signal" ,"f");
679  else l->AddEntry(hSig,"#bar{#nu_{e}} Signal" ,"f");
680  if (isFHC) l->AddEntry(hWrS,"Wrong Sign: #bar{#nu_{e}} CC" ,"f");
681  else l->AddEntry(hWrS,"Wrong Sign: #nu_{e} CC" ,"f");
682  }
683  l->AddEntry(hBeam,"Beam Background" ,"f");
684  l->AddEntry(hCosm,"Cosmic Background","f");
685  l->Draw();
686 
687  // Draw Preliminary()
688  Preliminary();
689  if (isFHC) CornerLabel("Neutrino beam");
690  else CornerLabel("Antineutrino beam");
691 
692  // And redraw axis
693  hSig -> Draw("axis same");
694  SimCan -> Update();
695 
696  // Finally, save the canvas.
697  SimCan -> SaveAs( TString("Plots/")+TString(CanNa)+TString(".pdf") );
698  SimCan -> Write ( CanNa.c_str() );
699  MakeTextFile ( CanNa , true );
700  return;
701 }
702 //......................................................
703 
tree Draw("slc.nhit")
size_t NCuts
Definition: MakeCutFlow.C:50
bin1_2sigma SetFillColor(3)
enum BeamMode kRed
c1 SetGridx()
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
Antineutrinos-only.
Definition: IPrediction.h:50
c1 Update()
(&#39; appearance&#39;)
Definition: IPrediction.h:18
correl_xv GetYaxis() -> SetDecimals()
void SetHistProp(TH1D *hist, int Colour, double Width, double Offset)
(&#39;beam &#39;)
Definition: IPrediction.h:15
T sqrt(T number)
Definition: d0nt_math.hpp:156
void CornerLabel(std::string &s)
Definition: numu_tools.h:145
bool isFHC
void MakeCutFlowVecs(std::vector< std::string > &NuMu_CutNames, std::vector< Cut > &NuMu_TieredCuts, std::vector< std::string > &Nue_CutNames, std::vector< Cut > &Nue_TieredCuts)
Definition: CutFlow_header.h:6
void MakeCanvasForSim(std::vector< double > Sig, std::vector< double > Wrong, std::vector< double > Beam, std::vector< double > Cosm, std::vector< std::string > CutNames, size_t NCuts, bool NuMu)
def BarChart(node, total, depth=0, startpos=0)
OStream cerr
Definition: OStream.cxx:7
double TotBkg(std::vector< double > integrals)
double Width(Resonance_t res)
resonance width (GeV)
std::string sFHC
gargamelle SetTitle("Gargamelle #nu_{e} CC data")
osc::OscCalcDumb calc
double GetHistInt(TH1D *hist)
osc::IOscCalcAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
Definition: Calcs.cxx:49
void MakeCanvasForData(std::vector< double > Data, std::vector< double > Sig, std::vector< double > Wrong, std::vector< double > Beam, std::vector< double > Cosm, std::vector< std::string > CutNames, size_t NCuts, bool NuMu)
ntuple SetFillStyle(1001)
#define M_PI
Definition: SbMath.h:34
Double_t ymax
Definition: plot.C:25
Charged-current interactions.
Definition: IPrediction.h:39
size_t TrimVector(std::vector< double > &Sig, std::vector< double > &WS, std::vector< double > &Bkg, std::vector< double > &Cosm, std::vector< double > &Data, std::vector< std::string > &Names)
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:535
std::string sBFPoint
correl_xv GetXaxis() -> SetDecimals()
const double kAna2020FHCLivetime
Definition: Exposures.h:237
const Cut Bkg
c1 SetLogx()
const double kAna2019FHCLivetime
Definition: Exposures.h:226
void MakeNueTable(std::vector< double > Sig, std::vector< double > WrSi, std::vector< double > Beam, std::vector< double > NuMu, std::vector< double > NuTau, std::vector< double > NC, std::vector< double > Data, std::vector< double > Cosm, std::vector< std::string > CutNames, bool ShowData)
double POTNorm
const double kAna2020FHCPOT
Definition: Exposures.h:233
const double kAna2020RHCPOT
Definition: Exposures.h:235
void Preliminary()
double LivNorm
Neutrinos-only.
Definition: IPrediction.h:49
OStream cout
Definition: OStream.cxx:6
(&#39; survival&#39;)
Definition: IPrediction.h:19
cout<< "--"<< endl;for(Int_t iP=1;iP<=hyz->GetNbinsX();iP++){for(Int_t iC=1;iC<=hyz->GetNbinsY();iC++){if(hyv->GetBinContent(iP, iC)>-999){goal_hyv-> SetBinContent(iP, iC,-(dy[iP-1][iC-1]))
void SetBin(TH1D *hist, std::string lab, size_t ind, double Val)
void MakeNuMuTable(std::vector< double > Sig, std::vector< double > WrSi, std::vector< double > Nue, std::vector< double > NuTau, std::vector< double > NC, std::vector< double > Data, std::vector< double > Cosm, std::vector< std::string > CutNames, bool ShowData)
void cc()
Definition: test_ana.C:28
const double kAna2020RHCLivetime
Definition: Exposures.h:238
enum BeamMode kViolet
Neutral-current interactions.
Definition: IPrediction.h:40
const double kAna2019FHCPOT
Definition: Exposures.h:223
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
Double_t ymin
Definition: plot.C:24
void MakeCutFlowPlots(bool isfhc, std::string Samp, std::string BFPoint="bf2020")
TFile * OutFile
histo SetMinimum(1.0e-9)
correl_xv SetStats(0)
void Beam(bool isRHC)
All neutrinos, any flavor.
Definition: IPrediction.h:26
cosmicTree SaveAs("cosmicTree.root")
Float_t e
Definition: plot.C:35
static std::unique_ptr< PredictionNoExtrap > LoadFrom(TDirectory *dir, const std::string &name)
enum BeamMode kGreen
TH1D * GetHist(TFile *InF, std::string InDir, int Key)
c cd(1)
void MakeTextFile(std::string OutName, bool BarChart)
std::string CutStr(std::string &Str)
T asin(T number)
Definition: d0nt_math.hpp:60
gm Write()
std::vector< std::string > CutNames
Definition: MakeCutFlow.C:49
enum BeamMode string