Plot_Evals_BDTMLP_Algs.C
Go to the documentation of this file.
1 //......................................................
5 
6 #include "CAFAna/Core/Loaders.h"
7 #include "CAFAna/Core/Spectrum.h"
10 
13 
15 
16 #include "Utilities/rootlogon.C"
17 
18 #include "OscLib/OscCalcPMNSOpt.h"
19 
20 #include "TCanvas.h"
21 #include "TColor.h"
22 #include "TFile.h"
23 #include "TGraph.h"
24 #include "TGraphAsymmErrors.h"
25 #include "TH1.h"
26 #include "TH2.h"
27 #include "TLatex.h"
28 #include "TLegend.h"
29 #include "TStyle.h"
30 #include "TSystem.h"
31 
32 #include <fstream>
33 #include <iostream>
34 #include <iomanip>
35 #include <string>
36 
37 bool isFHC = false; // Global variable for whether I am running over FHC or RHC?
38 std::string sFHC = "RHC"; // Global variable for string related to FHC or RHC.
39 
40 bool isBPF = false;
41 std::string sBPF = "Kal";
42 //......................................................
43 using namespace ana;
44 
45 //......................................................
46 TH1D* GetHistVectors( TFile* InF, std::string InDir, bool Cosm );
47 //......................................................
48 void WriteTable( double CosRejCosmInt, double CosRejMontInt, std::pair<TH1D*,TH1D*> Data, std::pair<TH1D*,TH1D*> Mont, std::string Name );
49 //......................................................
50 double FOMCalc( TH1D* Cosm, TH1D* Beam, bool DipFOM );
51 //......................................................
52 void SetBin( TH1D* hist, size_t ind, double Val );
53 //......................................................
54 void SetHistProp( TH1D* hist, int Colour, double Width, double Offset );
55 //......................................................
56 TH1D* MakeRatio( TH1D* num, TH1D* denom, int Col, int Sty, std::string FType );
57 //......................................................
58 void MakeOverlayCanvas( TH1D* BaseSpect, TH1D* DefCosRej, TH1D* BDTOpt, TH1D* BDTMatch, TH1D* MLPOpt, TH1D* MLPMatch, bool isCosm, std::string PlotTag );
59 //......................................................
60 void MakeVennDiagram ( TH1D* DefCosRej, TH1D* MyCut, TH1D* Overlap, bool isData, std::string CutName );
61 //......................................................
62 void MakeTextFile( std::string OutName, bool BarChart );
63 //......................................................
64 void CornerLabel(std::string Str);
65 //......................................................
66 // --- Define the names of my cuts.
67 std::vector<std::string> CutNames;
68 int NCuts = 0;
69 
70 // ---- Define my output directory
72 
73 double POTNorm = 0;
74 double LivNorm = 0;
75 
78 
79 //......................................................
80 void Plot_Evals_BDTMLP_Algs( std::string sIdentifier, bool fhc, bool BPF, bool Quant ) {
81  // Set my global variables.
82  isFHC = fhc;
83  sFHC = (isFHC == true ? "fhc":"rhc" );
84 
85  isBPF = BPF;
86  sBPF = (isBPF == true ? "BPF":"Kal" );
87 
88  if ( fhc && sIdentifier == "Period1" ) {
91  }
92  else if ( fhc && sIdentifier == "Period2" ) {
95  }
96  else if ( fhc && sIdentifier == "HighGain" ) {
99  }
100  else if (!fhc && sIdentifier == "HighGain" ) {
103  }
104  else {
105  std::cout << "\n\n I don't recognise sIdentifier = " << sIdentifier << " quitting.\n\n" << std::endl;
106  return;
107  }
108 
109  gStyle->SetHistMinimumZero();
110 
111  // Where are my input files?
112  //*
113 
114  std::string InBase = "/nova/ana/users/karlwarb/PostAna18/CosRej/VerifyBDT-190422/";
115 
116  //std::string MontNa = "BDTSpectra-"+sFHC+"-"+sBPF+"-"+sIdentifier+"-BeamMC-Evaluate.root";
117  //std::string DataNa = "BDTSpectra-"+sFHC+"-"+sBPF+"-"+sIdentifier+"-Cosmics-Evaluate.root";
118  std::string MontNa = "BDTSpectra-"+sFHC+"-"+sIdentifier+"-BeamMC-Evaluate.root";
119  std::string DataNa = "BDTSpectra-"+sFHC+"-"+sIdentifier+"-Cosmics-Evaluate.root";
120 
121  std::string MontLoc = InBase + MontNa;
122  std::string DataLoc = InBase + DataNa;
123 
124  std::cerr << "Loading;"
125  << "\n\t Monte Carlo - " << MontLoc
126  << "\n\t Data Files - " << DataLoc
127  << std::endl;
128 
129  TFile *MontFile = TFile::Open( MontLoc.c_str() );
130  TFile *DataFile = TFile::Open( DataLoc.c_str() );
131 
132  // --- Get my output file sorted out...
133  OutDir = "Plots-"+sFHC+"-"+sIdentifier+"/";
134  gSystem -> MakeDirectory( OutDir.c_str() );
135  std::string OutName = OutDir+"VerifyPlots_"+sFHC+".root";
136  TFile *OutFile = new TFile(OutName.c_str(), "RECREATE");
137 
138  // -----------------------------------------------------------
139  // --- What cuts do I want to use?
140  // -----------------------------------------------------------
141  Cut BaseCut = kNoCut;
142  std::vector<Cut> Cuts;
143  MakeCutVec( sIdentifier, isFHC, isBPF, Quant, BaseCut, Cuts, CutNames ); // The MC/Data flag means nothing here...
144  NCuts = (int)Cuts.size();
145  std::cout << "Out of my function. CutNames has size " << NCuts << std::endl;
146 
147  // -----------------------------------------------------------
148  // --- Now I want to extract my histograms...
149  // -----------------------------------------------------------
150  // --- First though I want to pull out the `MyBaseCut` and `DefCosRej` plots.
151  std::vector<TH1D*> RecoNuE_Data_MyBaseCut; std::vector<TH1D*> RecoNuE_Mont_MyBaseCut;
152  std::vector<TH1D*> RecoNuE_Data_DefCosRej; std::vector<TH1D*> RecoNuE_Mont_DefCosRej;
153  // Data first.
154  RecoNuE_Data_MyBaseCut.push_back( GetHistVectors( DataFile, "RecoNuE_MyBaseCut", true ) );
155  RecoNuE_Data_DefCosRej.push_back( GetHistVectors( DataFile, "RecoNuE_DefCosRej", true ) );
156  // Now Monte Carlo
157  RecoNuE_Mont_MyBaseCut.push_back( GetHistVectors( MontFile, "RecoNuE_MyBaseCut", false ) );
158  RecoNuE_Mont_DefCosRej.push_back( GetHistVectors( MontFile, "RecoNuE_DefCosRej", false ) );
159  // What were they?
160  std::cerr << "\n Got my base level numbers..."
161  << "\n\t RecoNuE_Data_MyBaseCut -> " << RecoNuE_Data_MyBaseCut[0]->Integral() << " && RecoNuE_Mont_MyBaseCut -> " << RecoNuE_Mont_MyBaseCut[0]->Integral()
162  << "\n\t RecoNuE_Data_DefCosRej -> " << RecoNuE_Data_DefCosRej[0]->Integral() << " && RecoNuE_Mont_DefCosRej -> " << RecoNuE_Mont_DefCosRej[0]->Integral()
163  << std::endl;
164  // --- Now load the quantiles....
165  if ( Quant ) {
166  for ( size_t cc=0; cc < 4; ++cc ) {
167  std::string QuantStr = "_Quant"+std::to_string(cc+1);
168  // Data first.
169  RecoNuE_Data_MyBaseCut.push_back( GetHistVectors( DataFile, "RecoNuE_MyBaseCut"+QuantStr, true ) );
170  RecoNuE_Data_DefCosRej.push_back( GetHistVectors( DataFile, "RecoNuE_DefCosRej"+QuantStr, true ) );
171  // Now Monte Carlo
172  RecoNuE_Mont_MyBaseCut.push_back( GetHistVectors( MontFile, "RecoNuE_MyBaseCut"+QuantStr, false ) );
173  RecoNuE_Mont_DefCosRej.push_back( GetHistVectors( MontFile, "RecoNuE_DefCosRej"+QuantStr, false ) );
174  // What were they?
175  std::cerr << "\n Got my base level numbers for " << QuantStr << "..."
176  << "\n\t RecoNuE_Data_MyBaseCut -> " << RecoNuE_Data_MyBaseCut[cc+1]->Integral() << " && RecoNuE_Mont_MyBaseCut -> " << RecoNuE_Mont_MyBaseCut[cc+1]->Integral()
177  << "\n\t RecoNuE_Data_DefCosRej -> " << RecoNuE_Data_DefCosRej[cc+1]->Integral() << " && RecoNuE_Mont_DefCosRej -> " << RecoNuE_Mont_DefCosRej[cc+1]->Integral()
178  << std::endl;
179  }
180  }
181 
182 
183  // --- Define my histograms.
184  std::vector< std::string > CutPointNames;
185  /* Data */
186  std::vector< std::pair<TH1D*,TH1D*> > RecoNuE_BDTOpt_Data;
187  std::vector< std::pair<TH1D*,TH1D*> > RecoNuE_BDTMat_Data;
188  std::vector< std::pair<TH1D*,TH1D*> > RecoNuE_MLPOpt_Data;
189  std::vector< std::pair<TH1D*,TH1D*> > RecoNuE_MLPMat_Data;
190  /* MC */
191  std::vector< std::pair<TH1D*,TH1D*> > RecoNuE_BDTOpt_Mont;
192  std::vector< std::pair<TH1D*,TH1D*> > RecoNuE_BDTMat_Mont;
193  std::vector< std::pair<TH1D*,TH1D*> > RecoNuE_MLPOpt_Mont;
194  std::vector< std::pair<TH1D*,TH1D*> > RecoNuE_MLPMat_Mont;
195 
196  // --- Go through my cut loop vector....
197  int Offset = 2;
198  if (Quant) Offset = 10;
199  for (int cc=Offset; cc<NCuts; cc=cc+2) {
200  // --- Which cut am I looking at now?
201  std::cout << "\nLooking at cut level " << CutNames[cc] << " and " << CutNames[cc+1] << std::endl;
202  if ( (!Quant && ( cc==Offset || cc==Offset+2 || cc==Offset+4 || cc==Offset+6 ) ) ||
203  ( Quant && ( cc==Offset || cc==Offset+10 || cc==Offset+20 || cc==Offset+30 ) ) ) {
204  CutPointNames.push_back( CutNames[cc] );
205  }
206 
207  // Extract the spectra for data.
208  TH1D* CosRej_Data = GetHistVectors( DataFile, "RecoNuE_"+CutNames[cc ], true );
209  TH1D* Overlap_Data = GetHistVectors( DataFile, "RecoNuE_"+CutNames[cc+1], true );
210 
211  // Extract the spectra, and push back into vector of doubles -- Monte Carlo.
212  TH1D* CosRej_Mont = GetHistVectors( MontFile, "RecoNuE_"+CutNames[cc ], false );
213  TH1D* Overlap_Mont = GetHistVectors( MontFile, "RecoNuE_"+CutNames[cc+1], false );
214 
215  std::cout << "\t CosRej_Data -> " << CosRej_Data->Integral() << " && Overlap_Data -> " << Overlap_Data->Integral() << "\n"
216  << "\t CosRej_Mont -> " << CosRej_Mont->Integral() << " && Overlap_Mont -> " << Overlap_Mont->Integral()
217  << std::endl;
218 
219  // -- BDT Optimised cut.
220  if ( (!Quant && cc == Offset ) || ( Quant && cc >= Offset && cc < Offset+10 ) ) {
221  RecoNuE_BDTOpt_Data.push_back( std::make_pair( CosRej_Data, Overlap_Data ) );
222  RecoNuE_BDTOpt_Mont.push_back( std::make_pair( CosRej_Mont, Overlap_Mont ) );
223  }
224  // -- MLP Optimised cut
225  else if ( (!Quant && cc == Offset+2) || ( Quant && cc >= Offset+10 && cc < Offset+20 ) ) {
226  RecoNuE_MLPOpt_Data.push_back( std::make_pair( CosRej_Data, Overlap_Data ) );
227  RecoNuE_MLPOpt_Mont.push_back( std::make_pair( CosRej_Mont, Overlap_Mont ) );
228  }
229  // -- BDT Matched cut
230  else if ( (!Quant && cc == Offset+4) || ( Quant && cc >= Offset+20 && cc < Offset+30 ) ) {
231  RecoNuE_BDTMat_Data.push_back( std::make_pair( CosRej_Data, Overlap_Data ) );
232  RecoNuE_BDTMat_Mont.push_back( std::make_pair( CosRej_Mont, Overlap_Mont ) );
233  }
234  // -- MLP Matched cut
235  else if ( (!Quant && cc == Offset+6) || ( Quant && cc >= Offset+30 && cc < Offset+40 ) ) {
236  RecoNuE_MLPMat_Data.push_back( std::make_pair( CosRej_Data, Overlap_Data ) );
237  RecoNuE_MLPMat_Mont.push_back( std::make_pair( CosRej_Mont, Overlap_Mont ) );
238  } else { std::cout << "\t\t Beyond my range, something has gone wrong..." << std::endl; return; }
239 
240  //OutFile->cd();
241  //MakeVennDiagram( RecoNuE_Data_DefCosRej, RecoNuE_Data.back().first, RecoNuE_Data.back().second, true , CutNames[cc] );
242  //MakeVennDiagram( RecoNuE_Mont_DefCosRej, RecoNuE_Mont.back().first, RecoNuE_Mont.back().second, false, CutNames[cc] );
243  }
244 
245  std::cout << "\n\nLet's check my vector sizes... Starting with CutPointNames --> " << CutPointNames.size()
246  << "\n\t RecoNuE_Data_MyBaseCut = " << RecoNuE_Data_MyBaseCut.size() << ", RecoNuE_Mont_MyBaseCut = " << RecoNuE_Mont_MyBaseCut.size()
247  << "\n\t RecoNuE_Data_DefCosRej = " << RecoNuE_Data_DefCosRej.size() << ", RecoNuE_Mont_DefCosRej = " << RecoNuE_Mont_DefCosRej.size()
248  << "\n\t RecoNuE_BDTOpt_Data = " << RecoNuE_BDTOpt_Data.size() << ", RecoNuE_BDTOpt_Mont = " << RecoNuE_BDTOpt_Mont.size()
249  << "\n\t RecoNuE_MLPOpt_Data = " << RecoNuE_MLPOpt_Data.size() << ", RecoNuE_MLPOpt_Mont = " << RecoNuE_MLPOpt_Mont.size()
250  << "\n\t RecoNuE_BDTMat_Data = " << RecoNuE_BDTMat_Data.size() << ", RecoNuE_BDTMat_Mont = " << RecoNuE_BDTMat_Mont.size()
251  << "\n\t RecoNuE_MLPMat_Data = " << RecoNuE_MLPMat_Data.size() << ", RecoNuE_MLPMat_Mont = " << RecoNuE_MLPMat_Mont.size()
252  << std::endl;
253 
254  // --- Now lets start drawing things!!!
255  std::cerr << "\n\nNow lets start plotting things..." << std::endl;
256  OutFile -> cd();
257 
258  // ===========================================================
259  // I want this next block in a loop through my vectors...
260  // ===========================================================
261 
262  for (size_t VecSize=0; VecSize<RecoNuE_Data_MyBaseCut.size(); ++VecSize) {
263  std::string sQuant = "Full";
264  if (VecSize==1) sQuant = "Quant1";
265  else if (VecSize==2) sQuant = "Quant2";
266  else if (VecSize==3) sQuant = "Quant3";
267  else if (VecSize==4) sQuant = "Quant4";
268 
269  // -----------------------------------------------------------
270  // --- What is my reference FOM?
271  // -----------------------------------------------------------
272  double FOM_Full_Def = FOMCalc( RecoNuE_Data_DefCosRej[VecSize], RecoNuE_Mont_DefCosRej[VecSize], false );
273  double FOM_Dip_Def = FOMCalc( RecoNuE_Data_DefCosRej[VecSize], RecoNuE_Mont_DefCosRej[VecSize], true );
274 
275  // -----------------------------------------------------------
276  // --- First of all lets look at the BDT values.
277  // -----------------------------------------------------------
278  double Int_DataMyBaseCut = RecoNuE_Data_MyBaseCut[VecSize]->Integral();
279  double Int_MontMyBaseCut = RecoNuE_Mont_MyBaseCut[VecSize]->Integral();
280  double Int_DataDefCosRej = RecoNuE_Data_DefCosRej[VecSize]->Integral();
281  double Int_MontDefCosRej = RecoNuE_Mont_DefCosRej[VecSize]->Integral();
282 
283  std::cout << std::setprecision(3) << std::fixed << "\n My input directory is " << InBase << "\n"
284  << "\n\t Cut Name \t\t Cosmic Data \t Overlap \t(frac) \t\t Monte Carlo \t Overlap \t (frac) \t FOM \t FOM-Dip"
285  << "\n\n"
286  // Events which pass the full selection, other than Cosmic Rejection.
287  << "\t Pass all PID \t"
288  << "\t " << Int_DataMyBaseCut << "\t\t " << "\t\t (N/A) \t"
289  << "\t " << std::setprecision(2) << std::fixed << Int_MontMyBaseCut << "\t " << "\t (N/A) \t \n"
290  // Events which pass the default cosmic rejection.
291  << "\t DefaultCosRej \t"
292  << "\t " << std::setprecision(3) << std::fixed << Int_DataDefCosRej << "\t\t " << "\t\t (N/A) \t"
293  << "\t " << std::setprecision(2) << std::fixed << Int_MontDefCosRej << "\t " << "\t (N/A) \t\t\t"
294  << "\t " << FOM_Full_Def << "\t " << FOM_Dip_Def << "\n"
295  << std::endl;
296 
297  // The BDT cut for the "optimal" signal / bkg ratio.
298  WriteTable( Int_DataDefCosRej, Int_MontDefCosRej, RecoNuE_BDTOpt_Data[VecSize], RecoNuE_BDTOpt_Mont[VecSize], CutPointNames[0] );
299  // The BDT cut for the matched signal efficiency.
300  WriteTable( Int_DataDefCosRej, Int_MontDefCosRej, RecoNuE_BDTMat_Data[VecSize], RecoNuE_BDTMat_Mont[VecSize], CutPointNames[2] );
301  std::cout << "" << std::endl;
302  // The MLP cut for the "optimal" signal / bkg ratio.
303  WriteTable( Int_DataDefCosRej, Int_MontDefCosRej, RecoNuE_MLPOpt_Data[VecSize], RecoNuE_MLPOpt_Mont[VecSize], CutPointNames[1] );
304  // The MLP cut for the matched signal efficiency.
305  WriteTable( Int_DataDefCosRej, Int_MontDefCosRej, RecoNuE_MLPMat_Data[VecSize], RecoNuE_MLPMat_Mont[VecSize], CutPointNames[3] );
306  std::cout << "\n\n" << std::endl;
307 
308  // --- Make a canvas showing RecoNuE for each cut tier.
309  MakeOverlayCanvas( RecoNuE_Data_MyBaseCut[VecSize] , RecoNuE_Data_DefCosRej[VecSize],
310  RecoNuE_BDTOpt_Data [VecSize].first, RecoNuE_BDTMat_Data [VecSize].first,
311  RecoNuE_MLPOpt_Data [VecSize].first, RecoNuE_MLPMat_Data [VecSize].first,
312  true , sFHC+"-"+sBPF+"-"+sIdentifier+"-"+sQuant+"-Cosm" );
313 
314  // --- Make a canvas with all of the MC plots.
315  MakeOverlayCanvas( RecoNuE_Mont_MyBaseCut[VecSize] , RecoNuE_Mont_DefCosRej[VecSize],
316  RecoNuE_BDTOpt_Mont [VecSize].first, RecoNuE_BDTMat_Mont [VecSize].first,
317  RecoNuE_MLPOpt_Mont [VecSize].first, RecoNuE_MLPMat_Mont [VecSize].first,
318  false, sFHC+"-"+sBPF+"-"+sIdentifier+"-"+sQuant+"-Mont" );
319  }
320 
321  return;
322 }
323 
324 //......................................................
325 TH1D* GetHistVectors( TFile* InF, std::string InDir, bool Cosm ) {
326  InF -> cd();
327  TH1D* TempHist;
328 
329  //std::cout << "In GetHistVectors: InDir is " << InDir << ", Cosm is " << Cosm << std::endl;
330 
331  // --- Convert my object to a TH1D.
332  if ( Cosm ) {
333  // --- If Cosmics then just have a spectrum
334  std::unique_ptr<Spectrum> Temp = Spectrum::LoadFrom( InF -> GetDirectory( InDir.c_str() ) );
335  TempHist = Temp -> ToTH1( NormToLiv, kLivetime );
336  } else {
337  // --- If MC have to convert from PredNoExtrap...
338  std::unique_ptr<PredictionNoExtrap> TempSpec = PredictionNoExtrap::LoadFrom( InF -> GetDirectory( InDir.c_str() ) );
339  // --- Need to set my oscillation calculator...
342  calc.SetdCP ( M_PI*0.915869 ); // preliminary Ana2017
343  calc.SetDmsq32( 0.00248809 );
344  calc.SetTh23 ( asin(sqrt(0.587685)) );
345  // --- Now I can convert to a histogram
346  TempHist = TempSpec -> Predict(&calc).ToTH1( NormToPOT );
347  }
348 
349  // --- If I have Reco NuE then want to bin normalise.
350  if (InDir.find("RecoNuE") != std::string::npos) {
351  //TempHist -> Scale( 0.1, "width" );
352  TempHist -> GetYaxis() -> SetTitle("Events / 0.1 Gev");
353  TempHist -> GetYaxis() -> CenterTitle();
354  TempHist -> GetXaxis() -> CenterTitle();
355  }
356 
357  // --- And return!
358  return TempHist;
359 }
360 
361 //......................................................
362 void WriteTable( double CosRejCosmInt, double CosRejMontInt, std::pair<TH1D*,TH1D*> Data, std::pair<TH1D*,TH1D*> Mont, std::string Name ) {
363 
364  double DataInt = Data.first ->Integral();
365  double DataOver = Data.second->Integral();
366  double MontInt = Mont.first ->Integral();
367  double MontOver = Mont.second->Integral();
368 
369  double Data_Denom = ( ( CosRejCosmInt - DataOver ) + DataOver + ( DataInt - DataOver ) );
370  double Mont_Denom = ( ( CosRejMontInt - MontOver ) + MontOver + ( MontInt - MontOver ) );
371 
372  // FOM.
373  double FOM_Full = FOMCalc( Data.first, Mont.first, false );
374  double FOM_Dip = FOMCalc( Data.first, Mont.first, true );
375 
376  // Text output.
377  std::cout << "\t "<< Name
378  << std::setprecision(3) << std::fixed << "\t "<< DataInt << "\t\t "<< DataOver << "\t\t ("<< DataOver/Data_Denom<<" ) "
379  << std::setprecision(2) << std::fixed << "\t "<< MontInt << "\t\t "<< MontOver << "\t\t ("<< MontOver/Mont_Denom<<" ) "
380  << std::setprecision(3) << std::fixed << "\t "<< FOM_Full << "\t " << FOM_Dip
381  << std::endl;
382 }
383 
384 //......................................................
385 double FOMCalc( TH1D* Cosm, TH1D* Beam, bool DipFOM ) {
386  double Cosm_Int = Cosm -> Integral();
387  double Beam_Int = Beam -> Integral();
388 
389  // -- If working out for the dip then need to reduce my integral range.
390  if (DipFOM) {
391  int bmin = Cosm -> GetXaxis() -> FindBin(1.0);
392  int bmax = Cosm -> GetXaxis() -> FindBin(2.0);
393 
394  Cosm_Int = Cosm -> Integral( bmin, bmax );
395  Beam_Int = Beam -> Integral( bmin, bmax );
396  }
397 
398  double FOM = Beam_Int / std::pow( Cosm_Int + Beam_Int, 0.5 );
399  return FOM;
400 }
401 
402 //......................................................
403 void SetBin( TH1D* hist, size_t ind, double Val ) {
404  hist -> SetBinContent( ind, Val );
405  std::string BinText = CutNames[ NCuts-ind ];
406  while (BinText.find("_") != std::string::npos) BinText.replace(BinText.find("_"),1," ");
407  hist -> GetXaxis() -> SetBinLabel( ind, BinText.c_str() );
408  return;
409 }
410 //......................................................
411 void SetHistProp( TH1D* hist, int Colour, double Width, double Offset ) {
412  hist -> SetFillColor( Colour );
413  hist -> SetBarWidth ( Width );
414  hist -> SetBarOffset( Offset );
415  hist -> SetMinimum ( 0.1 );
416  hist -> SetStats(0);
417  for (int bi=0; bi<hist->GetNbinsY(); ++bi) hist->SetBinContent( bi, 0 );
418  return;
419 }
420 //......................................................
421 TH1D* MakeRatio( TH1D* num, TH1D* denom, int Col, int Sty, std::string FType ) {
422  // Declare a new histogram.
423  /*
424  TH1D* ratHist = new TH1D( TString(FType)+TString("_ratio")+TString(std::to_string(Col))+TString(num->GetName()),
425  TString(";")+TString(num->GetXaxis()->GetTitle())+TString(";Ratio"),
426  num->GetNbinsX(),
427  num->GetXaxis()->GetXmin(),
428  num->GetXaxis()->GetXmax()
429  );
430  */
431  std::string ratName = num->GetName()+FType+"_ratio";
432  TH1D* ratHist = (TH1D*)num->Clone ( ratName.c_str() );
433  // Don't want any exponents?
434  ratHist->GetYaxis()->SetNoExponent();
435  ratHist->GetXaxis()->SetNoExponent();
436  // Divide the two input histograms...
437  ratHist->Divide( num , denom );
438  // Set Line and Marker colors / styles.
439  ratHist->SetLineColor ( Col );
440  ratHist->SetMarkerColor( Col );
441  ratHist->SetMarkerStyle( Sty );
442  //ratHist->SetMarkerSize ( 3 );
443  // Figure out what the range should be...
444  ratHist->GetYaxis()->SetRangeUser( 0.7, 1.3 );
445  // Return my shiny new ratio plot!
446  return ratHist;
447 }
448 //..................................
449 void MakeTextFile( std::string OutName, bool BarChart ) {
450  // Determine output name, and caption.
451  std::string FNa = "Plots/"+OutName+".txt";
452  while (OutName.find("_") != std::string::npos) OutName.replace(OutName.find("_"),1," ");
453  std::string Cap = "Plot showing the number of events passing cuts for " + OutName;
454  if (BarChart) {
455  std::string DaCap = "for expected data rates, taken from MC and Cosmic sidebands.";
456  if (OutName.find("Data") != std::string::npos) DaCap = "with NuMi data shown, along with estimates from MC and Cosmic sidebands.";
457  Cap = "A bar chart showing the number of events changes as subsequent cuts are applied, " + DaCap;
458  }
459  Cap += " Plot is for "+sFHC+".";
460  std::cout << "\nFile name: " << FNa << "\n\tCaption: " << Cap << std::endl;
461  // Write to file.
462  std::ofstream TxtOut ( FNa.c_str(), std::ofstream::out );
463  TxtOut << Cap;
464  TxtOut.close();
465  // Done.
466  return;
467 }
468 //......................................................
470  TLatex* CornLab = new TLatex(.15, .93, Str.c_str());
471  CornLab->SetTextColor(kGray+1);
472  CornLab->SetNDC();
473  CornLab->SetTextSize (2/30.);
474  CornLab->SetTextAlign(11);
475  CornLab->Draw();
476 }
477 //......................................................
478 void MakeOverlayCanvas( TH1D* BaseSpect, TH1D* DefCosRej, TH1D* BDTOpt, TH1D* BDTMatch, TH1D* MLPOpt, TH1D* MLPMatch, bool isCosm, std::string PlotTag ) {
479 
480  // --- Declare a canvas.
481  std::string NameForm = "NeutrinoE-"+PlotTag;
482  TCanvas *HistCan = new TCanvas( NameForm.c_str(), TString("Reco NuE numbers")+TString(NameForm) );
483  HistCan -> SetLeftMargin (0.15);
484  HistCan -> SetRightMargin(0.05);
485 
486  // --- Set the colours
487  std::vector<int> VecCols = { 1, 2, 4, 6, 8 };
488 
489  BaseSpect -> SetLineColor( kGray );
490  DefCosRej -> SetLineColor( VecCols[0] );
491  BDTOpt -> SetLineColor( VecCols[1] );
492  BDTMatch -> SetLineColor( VecCols[2] );
493  MLPOpt -> SetLineColor( VecCols[3] );
494  MLPMatch -> SetLineColor( VecCols[4] );
495 
496  // --- Make the ratio histograms
497  TH1D* ratBDTOpt = MakeRatio( BDTOpt , DefCosRej, VecCols[1], 20, PlotTag );
498  TH1D* ratBDTMatch = MakeRatio( BDTMatch, DefCosRej, VecCols[2], 21, PlotTag );
499  TH1D* ratMLPOpt = MakeRatio( MLPOpt , DefCosRej, VecCols[3], 22, PlotTag );
500  TH1D* ratMLPMatch = MakeRatio( MLPMatch, DefCosRej, VecCols[4], 23, PlotTag );
501 
502  // --- Get the integrals.
503  double mBaseSpect = BaseSpect -> Integral();
504  double mDefCosRej = DefCosRej -> Integral();
505  double mBDTOpt = BDTOpt -> Integral();
506  double mBDTMatch = BDTMatch -> Integral();
507  double mMLPOpt = MLPOpt -> Integral();
508  double mMLPMatch = MLPMatch -> Integral();
509 
510  // Ratio of surviving events.
511  double rBDTOpt = mBDTOpt / mDefCosRej;
512  double rBDTMatch = mBDTMatch / mDefCosRej;
513  double rMLPOpt = mMLPOpt / mDefCosRej;
514  double rMLPMatch = mMLPMatch / mDefCosRej;
515 
516  // --- Bin width normalise my plots!
517  BaseSpect -> Scale( 0.1, "width" );
518  DefCosRej -> Scale( 0.1, "width" );
519  BDTOpt -> Scale( 0.1, "width" );
520  BDTMatch -> Scale( 0.1, "width" );
521  MLPOpt -> Scale( 0.1, "width" );
522  MLPMatch -> Scale( 0.1, "width" );
523 
524  // --- If cosmics then set the y axis limit to 2.
525  if ( isCosm && PlotTag.find("fhc") != std::string::npos && PlotTag.find("HighGain") != std::string::npos ) { BaseSpect -> GetYaxis() -> SetRangeUser( 0, 0.2 ); }
526  else if ( isCosm && PlotTag.find("fhc") != std::string::npos && PlotTag.find("Period1" ) != std::string::npos ) { BaseSpect -> GetYaxis() -> SetRangeUser( 0, 0.1 ); }
527  else if ( isCosm && PlotTag.find("fhc") != std::string::npos && PlotTag.find("Period2" ) != std::string::npos ) { BaseSpect -> GetYaxis() -> SetRangeUser( 0, 0.2 ); }
528  else if ( isCosm && PlotTag.find("rhc") != std::string::npos && PlotTag.find("HighGain") != std::string::npos ) { BaseSpect -> GetYaxis() -> SetRangeUser( 0, 0.1 ); }
529 
530  if (!isCosm && PlotTag.find("fhc") != std::string::npos && PlotTag.find("HighGain") != std::string::npos ) {
531  if ( PlotTag.find("Full") != std::string::npos ) { BaseSpect -> GetYaxis() -> SetRangeUser( 0, 6.50 ); }
532  else { BaseSpect -> GetYaxis() -> SetRangeUser( 0, 1.70 ); }
533  }
534  if (!isCosm && PlotTag.find("fhc") != std::string::npos && PlotTag.find("Period1" ) != std::string::npos ) {
535  if ( PlotTag.find("Full") != std::string::npos ) { BaseSpect -> GetYaxis() -> SetRangeUser( 0, 0.60 ); }
536  else { BaseSpect -> GetYaxis() -> SetRangeUser( 0, 0.15 ); }
537  }
538  if (!isCosm && PlotTag.find("fhc") != std::string::npos && PlotTag.find("Period2" ) != std::string::npos ) {
539  if ( PlotTag.find("Full") != std::string::npos ) { BaseSpect -> GetYaxis() -> SetRangeUser( 0, 3.00 ); }
540  else { BaseSpect -> GetYaxis() -> SetRangeUser( 0, 0.75 ); }
541  }
542  else if (!isCosm && PlotTag.find("rhc") != std::string::npos && PlotTag.find("HighGain") != std::string::npos ) {
543  if ( PlotTag.find("Full") != std::string::npos ) { BaseSpect -> GetYaxis() -> SetRangeUser( 0, 3.60 ); }
544  else { BaseSpect -> GetYaxis() -> SetRangeUser( 0, 0.90 ); }
545  }
546 
547 
548  // Draw the plots!
549  BaseSpect -> Draw("hist");
550  DefCosRej -> Draw("hist same");
551 
552  BDTOpt -> Draw("same p");
553  BDTMatch -> Draw("same p");
554  MLPOpt -> Draw("same p");
555  MLPMatch -> Draw("same p");
556 
557  // Add a legend to the plot.
558  TLegend *l = new TLegend(0.65,0.03,0.93,0.88);
559  l->AddEntry(BaseSpect,"Pass PID Cuts","l");
560  l->AddEntry(BaseSpect, Form("Count: %.3f",mBaseSpect), "" );
561 
562  l->AddEntry(DefCosRej,"Def Cos Rej","l");
563  l->AddEntry(DefCosRej, Form("Count: %.3f",mDefCosRej), "" );
564 
565  l->AddEntry(BDTOpt ,"BDT Opt Val","l");
566  l->AddEntry(BDTOpt , Form("Count (ratio): %.3f, (%.3f)",mBDTOpt , rBDTOpt ), "" );
567 
568  l->AddEntry(BDTMatch ,"BDT Match" ,"l");
569  l->AddEntry(BDTMatch , Form("Count (ratio): %.3f, (%.3f)",mBDTMatch , rBDTMatch ), "" );
570 
571  l->AddEntry(MLPOpt ,"MLP Opt Val","l");
572  l->AddEntry(MLPOpt , Form("Count (ratio): %.3f, (%.3f)",mMLPOpt , rMLPOpt ), "" );
573 
574  l->AddEntry(MLPMatch ,"MLP Match" ,"l");
575  l->AddEntry(MLPMatch , Form("Count (ratio): %.3f, (%.3f)",mMLPMatch , rMLPMatch ), "" );
576 
577  l->Draw();
578  // Add some hangings on
579  if (isFHC) CornerLabel("Neutrino beam");
580  else CornerLabel("Antineutrino beam");
581 
582  // -----------------------------------------------------------
583  // --- Now for the ratio part of the plots....
584  // -----------------------------------------------------------
585  std::string NameForm_rat = NameForm+"_ratio";
586  TCanvas *RatCan = new TCanvas( NameForm_rat.c_str(), TString("Reco NuE numbers")+TString(NameForm_rat) );
587  // --- If looking at cosmics then I want the plot to be log y
588  ratBDTOpt -> GetYaxis() -> SetTitle("Ratio");
589  if (isCosm) {
590  ratBDTOpt -> GetYaxis() -> SetRangeUser( 5e-2, 20 );
591  ratBDTOpt -> GetYaxis() -> SetTitle("Ratio (Log axis)");
592  RatCan -> SetLogy();
593  }
594  ratBDTOpt -> Draw("hist p");
595  // --- Add a reference line
596  double x[2] = { ratBDTOpt->GetXaxis()->GetXmin(), ratBDTOpt->GetXaxis()->GetXmax() };
597  double y[2] = { 1, 1 };
598  TGraph *line = new TGraph( 2, x, y );
599  line->SetLineColor(kGray);
600  line->SetLineWidth( 3 );
601  line->SetLineStyle( 2 );
602  line->Draw("l same");
603  // --- Finally, draw!
604  ratBDTOpt -> Draw("axis hist p same");
605  ratBDTMatch -> Draw("hist p same");
606  ratMLPOpt -> Draw("hist p same");
607  ratMLPMatch -> Draw("hist p same");
608 
609  // -----------------------------------------------------------
610  // --- Now merge the two canvases....
611  // -----------------------------------------------------------
612  std::string NameForm_merge = NameForm+"_merged";
613  TCanvas *MergeCan = new TCanvas( NameForm_merge.c_str(), TString("Reco NuE numbers")+TString(NameForm_merge) );
614  MergeCan->Divide(1,2);
615  HistCan -> SetRightMargin(.03); HistCan -> SetLeftMargin(.07); HistCan -> SetTopMargin(.10); HistCan -> SetBottomMargin(.00);
616  RatCan -> SetRightMargin(.03); RatCan -> SetLeftMargin(.07); RatCan -> SetTopMargin(.00); RatCan -> SetBottomMargin(.12);
617 
618  // And now the axes...
619  double Lb1 = 0.07;
620  BaseSpect -> GetYaxis() -> SetTitleSize( Lb1 ); BaseSpect -> GetYaxis() -> SetLabelSize( Lb1 ); BaseSpect -> GetYaxis() -> SetTitleOffset( 0.45 );
621  ratBDTOpt -> GetYaxis() -> SetTitleSize( Lb1 ); ratBDTOpt -> GetYaxis() -> SetLabelSize( Lb1 ); ratBDTOpt -> GetYaxis() -> SetTitleOffset( 0.45 );
622  ratBDTOpt -> GetXaxis() -> SetTitleSize( Lb1 ); ratBDTOpt -> GetXaxis() -> SetLabelSize( Lb1 );
623 
624  MergeCan->cd(1);
625  HistCan->DrawClonePad();
626  MergeCan->cd(2);
627  RatCan ->DrawClonePad();
628 
629  // Finally, save the merged canvas.
630  MergeCan -> SaveAs( TString(OutDir)+TString(NameForm_merge)+TString(".pdf" ) );
631  MergeCan -> SaveAs( TString(OutDir)+TString(NameForm_merge)+TString(".png" ) );
632  MergeCan -> Write ( NameForm_merge.c_str() );
633 
634  return;
635 }
636 //......................................................
637 void MakeVennDiagram( TH1D* DefCosRej, TH1D* MyCut, TH1D* Overlap, bool isData, std::string CutName ) {
638 
639  // --- Am I looking at Data or MC? Also make my name and title...
640  std::string sData = "Data";
641  if (!isData) sData = "Mont";
642  std::string ThisName = sData+"_Overlap_"+CutName;
643  std::string TheTitle = "Overlap in "+sData+" of DefCosRej and "+CutName;
644 
645  std::cout << "\n\n Now looking at CutName " << CutName
646  << "\n\t DefCosRej is " << DefCosRej->Integral() << ", MyCut is " << MyCut->Integral() << ", Overlap is " << Overlap->Integral()
647  << std::endl;
648 
649  double nDefOnly = DefCosRej->Integral() - Overlap->Integral();
650  double nOverlap = Overlap ->Integral();
651  double nNewOnly = MyCut ->Integral() - Overlap->Integral();
652 
653  // --- Make the Canvas.
654  TCanvas *cOverlap = new TCanvas( TString(ThisName), TString(TheTitle) );
655  if (!isData) cOverlap->SetLogy();
656  // --- Make the Histogram.
657  TH1D* hOverlap = new TH1D( TString("h")+TString(ThisName), TString(TheTitle), 3, 0, 2 );
658  hOverlap->GetXaxis()->SetBinLabel( 1, "DefCosRej" ); hOverlap->SetBinContent( 1, nDefOnly );
659  hOverlap->GetXaxis()->SetBinLabel( 2, "Overlap" ); hOverlap->SetBinContent( 2, nOverlap );
660  hOverlap->GetXaxis()->SetBinLabel( 3, TString(CutName) ); hOverlap->SetBinContent( 3, nNewOnly );
661  hOverlap->Draw();
662  // --- Make the Legend.
663  TLatex text;
664  text.SetTextAlign(22);
665  text.SetTextSize(0.04);
666  //TLegend* place = AutoPlaceLegend(0.0001, 0.0001, 0.54);
667  std::string HornText = "#bf{Neutrino mode}";
668  if (!isFHC) HornText = "#bf{Antineutrino mode}";
669  std::string DataText = "Default Cos Rej: " + std::to_string( nDefOnly );
670  std::string OverText = "Overlap of Algs: " + std::to_string( nOverlap );
671  std::string NewCText = CutName +": " + std::to_string( nNewOnly );
672  std::string Words = "#splitline{"+HornText+"}{#splitline{"+DataText+"}{#splitline{"+OverText+"}{"+NewCText+"}}}";
673  text.DrawLatexNDC( 0.3, 0.78, Words.c_str() );
674  // --- Save the canvas!
675  cOverlap -> SaveAs( TString(OutDir)+TString(ThisName)+TString(".pdf" ) );
676  cOverlap -> SaveAs( TString(OutDir)+TString(ThisName)+TString(".png" ) );
677  cOverlap -> Write ( TString(ThisName) );
678  //cOverlap -> Close();
679 }
680 //......................................................
681 
682 
683 /*
684  //......................................................
685  // --- Make some area normalised plots.
686  //......................................................
687  // Clone the histos.
688  TH1D* QualCont_Norm = (TH1D*)QualCont ->Clone(); QualCont_Norm -> Scale(1/mQualCont );
689  TH1D* Ana18CosRej_Norm = (TH1D*)Ana18CosRej ->Clone(); Ana18CosRej_Norm -> Scale(1/mAna18CosRej );
690  TH1D* MLPDefCosRej_Norm = (TH1D*)MLPDefCosRej->Clone(); MLPDefCosRej_Norm -> Scale(1/mMLPDefCosRej);
691  TH1D* MLPCorCosRej_Norm = (TH1D*)MLPCorCosRej->Clone(); MLPCorCosRej_Norm -> Scale(1/mMLPCorCosRej);
692  TH1D* BDTDefCosRej_Norm = (TH1D*)BDTDefCosRej->Clone(); BDTDefCosRej_Norm -> Scale(1/mBDTDefCosRej);
693  TH1D* BDTCorCosRej_Norm = (TH1D*)BDTCorCosRej->Clone(); BDTCorCosRej_Norm -> Scale(1/mBDTCorCosRej);
694  double HistMax = 1.2 * std::max( QualCont_Norm -> GetMaximum(),
695  std::max( Ana18CosRej_Norm -> GetMaximum(),
696  std::max( MLPDefCosRej_Norm -> GetMaximum(),
697  std::max( MLPCorCosRej_Norm -> GetMaximum(),
698  std::max( BDTDefCosRej_Norm -> GetMaximum(), BDTCorCosRej_Norm -> GetMaximum() ) ) ) ) );
699  Ana18CosRej_Norm -> GetYaxis() -> SetRangeUser( 0, HistMax );
700  // Make canvas, and draw plots and legend
701  TCanvas *AreaCan = new TCanvas( TString(NameForm)+TString("_Area"), "Reco NuE numbers Area Norm" );
702  AreaCan -> SetLeftMargin (0.15);
703  AreaCan -> SetRightMargin(0.05);
704  QualCont_Norm -> Draw("hist");
705  Ana18CosRej_Norm -> Draw("hist same");
706  MLPDefCosRej_Norm -> Draw("hist same");
707  MLPCorCosRej_Norm -> Draw("hist same");
708  BDTDefCosRej_Norm -> Draw("hist same");
709  BDTCorCosRej_Norm -> Draw("hist same");
710  l->Draw();
711  // Add some hangings on
712  if (isFHC) CornerLabel("Neutrino beam");
713  else CornerLabel("Antineutrino beam");
714  // Finally, save the canvas.
715  AreaCan -> SaveAs( TString(OutDir)+TString(NameForm)+TString("_Area")+TString(".pdf" ) );
716  AreaCan -> SaveAs( TString(OutDir)+TString(NameForm)+TString("_Area")+TString(".png" ) );
717  AreaCan -> Write ( TString(NameForm)+TString("_Area") );
718 
719 */
tree Draw("slc.nhit")
bin1_2sigma SetFillColor(3)
std::string sFHC
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
void ResetOscCalcToDefault(osc::IOscCalcAdjustable *calc)
Reset calculator to default assumptions for all parameters.
Definition: Calcs.cxx:23
TH1D * MakeRatio(TH1D *num, TH1D *denom, int Col, int Sty, std::string FType)
std::string sBPF
const double kSecondAnaPeriod2POT
Definition: Exposures.h:74
correl_xv GetYaxis() -> SetDecimals()
double Integral(const Spectrum &s, const double pot, int cvnregion)
T sqrt(T number)
Definition: d0nt_math.hpp:156
void SetHistProp(TH1D *hist, int Colour, double Width, double Offset)
constexpr T pow(T x)
Definition: pow.h:75
const double kAna2017Period5POT
Definition: Exposures.h:171
def BarChart(node, total, depth=0, startpos=0)
OStream cerr
Definition: OStream.cxx:7
void WriteTable(double CosRejCosmInt, double CosRejMontInt, std::pair< TH1D *, TH1D * > Data, std::pair< TH1D *, TH1D * > Mont, std::string Name)
double Width(Resonance_t res)
resonance width (GeV)
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 kSecondAnaPeriod1POT
Definition: Exposures.h:73
gargamelle SetTitle("Gargamelle #nu_{e} CC data")
double FOMCalc(TH1D *Cosm, TH1D *Beam, bool DipFOM)
osc::OscCalcDumb calc
const double kAna2017Period3Livetime
Definition: Exposures.h:199
#define M_PI
Definition: SbMath.h:34
const double kAna2018RHCPOT
Definition: Exposures.h:208
void CornerLabel(std::string Str)
double POTNorm
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:535
bool isFHC
void MakeCutVec(std::string sIdentifier, bool isFHC, bool isBPF, bool Quants, Cut &kBaseCut, std::vector< Cut > &Cuts, std::vector< std::string > &CutNames)
correl_xv GetXaxis() -> SetDecimals()
hmean SetLineColor(4)
Optimized version of OscCalcPMNS.
Definition: StanTypedefs.h:31
void MakeVennDiagram(TH1D *DefCosRej, TH1D *MyCut, TH1D *Overlap, bool isData, std::string CutName)
void Plot_Evals_BDTMLP_Algs(std::string sIdentifier, bool fhc, bool BPF, bool Quant)
void SetTh23(const T &th23) override
std::vector< std::string > CutNames
const double kAna2017Period5Livetime
Definition: Exposures.h:198
OStream cout
Definition: OStream.cxx:6
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 SetDmsq32(const T &dmsq32) override
void SetBin(TH1D *hist, size_t ind, double Val)
const double kSecondAnaPeriod1Livetime
Definition: Exposures.h:99
TH1D * GetHistVectors(TFile *InF, std::string InDir, bool Cosm)
void cc()
Definition: test_ana.C:28
int num
Definition: f2_nu.C:119
void SetdCP(const T &dCP) override
double LivNorm
bool isBPF
void MakeOverlayCanvas(TH1D *BaseSpect, TH1D *DefCosRej, TH1D *BDTOpt, TH1D *BDTMatch, TH1D *MLPOpt, TH1D *MLPMatch, bool isCosm, std::string PlotTag)
double NormToPOT
const double kSecondAnaPeriod2Livetime
Definition: Exposures.h:100
TFile * OutFile
histo SetMinimum(1.0e-9)
const double kAna2018FHCPOT
Definition: Exposures.h:207
correl_xv SetStats(0)
void Beam(bool isRHC)
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
cosmicTree SaveAs("cosmicTree.root")
Float_t e
Definition: plot.C:35
static std::unique_ptr< PredictionNoExtrap > LoadFrom(TDirectory *dir, const std::string &name)
void MakeTextFile(std::string OutName, bool BarChart)
std::string OutDir
simulatedPeak Scale(1/simulationNormalisationFactor)
const double kAna2017Period3POT
Definition: Exposures.h:172
c cd(1)
const Cut kNoCut
The simplest possible cut: pass everything, used as a default.
Definition: Cut.h:109
double NormToLiv
const double kAna2018FHCLivetime
Definition: Exposures.h:213
gPad SetLogy()
const double kAna2018RHCLivetime
Definition: Exposures.h:214
T asin(T number)
Definition: d0nt_math.hpp:60
gm Write()
enum BeamMode string