NuMu2020_BasicPIDPlots_Plot.C
Go to the documentation of this file.
1 #include "TCanvas.h"
2 #include "TFile.h"
3 #include "TH1D.h"
4 #include "TFeldmanCousins.h"
5 #include "TGraph.h"
6 #include "TLine.h"
7 #include "TLatex.h"
8 #include "TLegend.h"
9 #include "TStyle.h"
10 #include "TSystem.h"
11 
12 #include "CAFAna/Analysis/Calcs.h"
14 #include "CAFAna/Core/Spectrum.h"
16 
18 
19 #include "OscLib/OscCalcPMNSOpt.h"
20 #include "Utilities/rootlogon.C"
21 
22 #include <iostream>
23 #include <sstream>
24 #include <fstream>
25 #include <vector>
26 #include <string>
27 
28 #include <TROOT.h>
29 #include <TStyle.h>
30 
31 using namespace ana;
32 
33 //......................................................
35  TH1D* hNuRecoE ;
36  /*
37  TH1D* hReMIdScore ;
38  TH1D* hProngScoremuon;
39  TH1D* hCosRejScoreP4 ;
40  TH1D* hCosRejScoreP5 ;
41  TH1D* hOldPreselmuon ;
42  TH1D* hOldPreselcosm ;
43  TH1D* hLoosePTPmuon ;
44  TH1D* hLoosePTPcosm ;
45  */
46 };
47 //......................................................
49  TH1D* hNuRecoE ;
50  /*
51  TH1D* hReMIdScore ;
52  TH1D* hProngScoremuon;
53  TH1D* hCosRejScoreP4 ;
54  TH1D* hCVNmuon ;
55  */
56 };
57 
58 
59 //......................................................
60 // --- Some global vars ---
61 bool bIsFHC = true; std::string sIsFHC = "";
62 bool bIsFD = true; std::string sIsFD = "";
63 // --- What to normalise to?
64 double POTNorm = 1e20;
65 double LivNorm = 1e20;
66 // --- Random Vars
67 size_t NPlot = 1;
68 // --- Output
70 TFile* OutFile;
71 
72 //......................................................
73 LoadedHistograms_2020 LoadFile_GetHists_2020( TFile* InFile, std::string CutTier, bool isData );
74 LoadedHistograms_2019 LoadFile_GetHists_2019( TFile* InFile, std::string CutTier, bool isData );
75 //......................................................
76 TH1D* GetSpectToHist( TFile* MyF, std::string LoadName, std::string Axis, bool isData );
77 //......................................................
78 void ProducePlots_2020( LoadedHistograms_2020 NoTrue_2020, LoadedHistograms_2020 TrueNu_2020,
79  LoadedHistograms_2020 TrNuMu_2020, LoadedHistograms_2020 AnNuMu_2020,
80  LoadedHistograms_2019 NoTrue_2019, LoadedHistograms_2020 Data_2020 , std::string CutTier );
81 //......................................................
82 void ProducePlots_2019( LoadedHistograms_2019 NoTrue_2019, LoadedHistograms_2019 TrueNu_2019,
83  LoadedHistograms_2019 TrNuMu_2019, LoadedHistograms_2019 AnNuMu_2019,
84  LoadedHistograms_2020 NoTrue_2020, LoadedHistograms_2019 Data_2019 , std::string CutTier );
85 //......................................................
86 TH1D* MakeRatio( TH1D* num, TH1D* denom, int Col, std::string FType );
87 //......................................................
88 void SetRange( TH1D* MC, TH1D* Data, // Pass the histograms
89  double XLow, double XHigh, // What X values do I want?
90  double H_YLow, double H_YHigh ); // What Y values do I want?
91 //......................................................
92 void FindAxisRange( double &XLow, double &XHigh, double &H_YLow, double &H_YHigh, double &R_YLow, double &R_YHigh, double MaxVal, bool &SetLogy, std::string Name );
93 //......................................................
94 void MakeSplitCans( TH1D* MCHist, TH1D* DataHist, std::string PlotName );
95 //......................................................
96 void PlotHistProp( TH1D* Hist, TLegend* L );
97 //......................................................
98 
99 //................................................................................
100 //................................................................................
101 void NuMu2020_BasicPIDPlots_Plot( bool isfhc, bool isfd, std::string period = "full" ) {
102 
103  // ---- First off, lets set the styles...
104  gStyle->SetOptStat(0);
105  gROOT->SetStyle("novaStyle");
106 
107  // --- Set my global variables and print them out to the screen...
108  bIsFHC = isfhc ; sIsFHC = "fhc" ; if (!isfhc ) {sIsFHC = "rhc" ; }
109  bIsFD = isfd ; sIsFD = "fd" ; if (!isfd ) {sIsFD = "nd" ;}
110  std::cout << "\n================================= \n"
111  << "\n\t bIsFHC --> " << bIsFHC << " --> " << sIsFHC
112  << "\n\t bIsFD --> " << bIsFD << " --> " << sIsFD
113  << "\n\t Period --> " << period
114  << "\n================================= \n"
115  << std::endl;
116 
117  // --- What do I normalise to? Just take the Ana19 numbers...
119  else if ( !bIsFHC ) { POTNorm = kAna2019RHCPOT; LivNorm = kAna2019RHCLivetime; } // RHC
120 
121  // --- Get my cut vector sorted...
122  std::vector<std::pair<Cut, std::string> > My2020_Mont = BasicPIDPlots_Cuts ( bIsFHC, bIsFD, false );
123  std::vector<std::pair<Cut, std::string> > My2020_Data = BasicPIDPlots_Cuts ( bIsFHC, bIsFD, true );
124  std::vector<std::pair<Cut, std::string> > My2019_Mont = BasicPIDPlots2019_Cuts( bIsFHC, bIsFD, false );
125  std::vector<std::pair<Cut, std::string> > My2019_Data = BasicPIDPlots2019_Cuts( bIsFHC, bIsFD, true );
126  // --- Get my cut vector size sorted...
127  size_t NumCut2020_Mont = My2020_Mont.size();
128  size_t NumCut2020_Data = My2020_Data.size();
129  size_t NumCut2019_Mont = My2019_Mont.size();
130  size_t NumCut2019_Data = My2019_Data.size();
131  // --- Write the size of the cuts
132  std::cout << "\nThe sizes of the Cut vectors are;"
133  << "\n\t NumCut2020_Mont: " << NumCut2020_Mont << "\t NumCut2020_Data: " << NumCut2020_Data
134  << "\n\t NumCut2019_Mont: " << NumCut2019_Mont << "\t NumCut2019_Data: " << NumCut2019_Data
135  << std::endl;
136 
137  // --- Where do my input files live?
138  std::string BaseCAFFileDir="/nova/ana/users/karlwarb/Analysis2020/PIDOpt/BasicPlot/200128/";
139 
140  // --- What are my input files called?
141  std::string InputMont_2020 = BaseCAFFileDir+"BasicPIDPlots_2020_"+sIsFD+"_mc_" +sIsFHC+"_"+period+".root";
142  std::string InputData_2020 = BaseCAFFileDir+"BasicPIDPlots_2020_"+sIsFD+"_cosm_"+sIsFHC+"_"+period+".root";
143  std::string InputMont_2019 = BaseCAFFileDir+"BasicPIDPlots_2019_"+sIsFD+"_mc_" +sIsFHC+"_"+period+".root";
144  std::string InputData_2019 = BaseCAFFileDir+"BasicPIDPlots_2019_"+sIsFD+"_cosm_"+sIsFHC+"_"+period+".root";
145  // --- Open the files.
146  std::cout << "\n Going to load input files: "
147  << "\n\t" << InputMont_2020 << "\n\t" << InputData_2020
148  << "\n\t" << InputMont_2019 << "\n\t" << InputData_2019 << std::endl;
149  TFile *MontFile_2020 = TFile::Open( InputMont_2020.c_str() );
150  TFile *MontFile_2019 = TFile::Open( InputMont_2019.c_str() );
151  TFile *DataFile_2020 = TFile::Open( InputData_2020.c_str() );
152  TFile *DataFile_2019 = TFile::Open( InputData_2019.c_str() );
153 
154  // --- Declare my various Loaded Histograms...
155  // 2020 MC
156  std::vector< LoadedHistograms_2020 > LH_MC_NoTrue_2020;
157  std::vector< LoadedHistograms_2020 > LH_MC_TrueNu_2020;
158  std::vector< LoadedHistograms_2020 > LH_MC_TrNuMu_2020;
159  std::vector< LoadedHistograms_2020 > LH_MC_AnNuMu_2020;
160  // 2020 Data
161  std::vector< LoadedHistograms_2020 > LH_Data_2020;
162  // 2019 MC
163  std::vector< LoadedHistograms_2019 > LH_MC_NoTrue_2019;
164  std::vector< LoadedHistograms_2019 > LH_MC_TrueNu_2019;
165  std::vector< LoadedHistograms_2019 > LH_MC_TrNuMu_2019;
166  std::vector< LoadedHistograms_2019 > LH_MC_AnNuMu_2019;
167  // 2019 Data
168  std::vector< LoadedHistograms_2019 > LH_Data_2019;
169 
170  // --- Now to load all of my histograms!
171  // --- 2020 MC...
172  std::cout << "\nLets load the 2020 Monte Carlo." << std::endl;
173  for (size_t cut=0; cut<NumCut2020_Mont; ++cut) {
174  // --- What is the Cut + Quant combination?
175  std::string CutName = My2020_Mont[cut].second;
176  //std::cout << "\n\nNow to look at " << CutName << std::endl;
177  // --- Load the Monte Carlo histograms.
178  if (CutName.find("AntiNuMu") != std::string::npos) LH_MC_AnNuMu_2020.push_back( LoadFile_GetHists_2020( MontFile_2020, CutName, false ) );
179  else if (CutName.find("TrueNuMu") != std::string::npos) LH_MC_TrNuMu_2020.push_back( LoadFile_GetHists_2020( MontFile_2020, CutName, false ) );
180  else if (CutName.find("TrueNu" ) != std::string::npos) LH_MC_TrueNu_2020.push_back( LoadFile_GetHists_2020( MontFile_2020, CutName, false ) );
181  else if (CutName.find("NoTrue" ) != std::string::npos) LH_MC_NoTrue_2020.push_back( LoadFile_GetHists_2020( MontFile_2020, CutName, false ) );
182  }
183  // --- 2020 Data
184  std::cout << "\nLets load the 2020 Data." << std::endl;
185  for (size_t cut=0; cut<NumCut2020_Data; ++cut) {
186  // --- What is the Cut + Quant combination?
187  std::string CutName = My2020_Data[cut].second;
188  //std::cout << "\n\nNow to look at " << CutName << std::endl;
189  // --- Load the Data histograms.
190  LH_Data_2020.push_back( LoadFile_GetHists_2020( DataFile_2020, CutName, true ) );
191  }
192 
193  // --- 2019 MC...
194  std::cout << "\nLets load the 2019 Monte Carlo." << std::endl;
195  for (size_t cut=0; cut<NumCut2019_Mont; ++cut) {
196  // --- What is the Cut + Quant combination?
197  std::string CutName = My2019_Mont[cut].second;
198  //std::cout << "\n\nNow to look at " << CutName << std::endl;
199  // --- Load the Monte Carlo histograms.
200  if (CutName.find("AntiNuMu") != std::string::npos) LH_MC_AnNuMu_2019.push_back( LoadFile_GetHists_2019( MontFile_2019, CutName, false ) );
201  else if (CutName.find("TrueNuMu") != std::string::npos) LH_MC_TrNuMu_2019.push_back( LoadFile_GetHists_2019( MontFile_2019, CutName, false ) );
202  else if (CutName.find("TrueNu" ) != std::string::npos) LH_MC_TrueNu_2019.push_back( LoadFile_GetHists_2019( MontFile_2019, CutName, false ) );
203  else if (CutName.find("NoTrue" ) != std::string::npos) LH_MC_NoTrue_2019.push_back( LoadFile_GetHists_2019( MontFile_2019, CutName, false ) );
204  }
205  // --- 2019 Data
206  std::cout << "\nLets load the 2019 Data." << std::endl;
207  for (size_t cut=0; cut<NumCut2019_Data; ++cut) {
208  // --- What is the Cut + Quant combination?
209  std::string CutName = My2019_Data[cut].second;
210  //std::cout << "\n\nNow to look at " << CutName << std::endl;
211  // --- Load the Data histograms.
212  LH_Data_2019.push_back( LoadFile_GetHists_2019( DataFile_2019, CutName, true ) );
213  }
214  // --- Check the sizes of all of my LoadedHistogram vectors...
215  std::cout << "\nThe size of my various LH vectors are as follows;"
216  << "\n\t LH_MC_NoTrue_2020 " << LH_MC_NoTrue_2020.size()
217  << "\n\t LH_MC_TrueNu_2020 " << LH_MC_TrueNu_2020.size()
218  << "\n\t LH_MC_TrNuMu_2020 " << LH_MC_TrNuMu_2020.size()
219  << "\n\t LH_MC_AnNuMu_2020 " << LH_MC_AnNuMu_2020.size()
220  << "\n\t LH_MC_NoTrue_2019 " << LH_MC_NoTrue_2019.size()
221  << "\n\t LH_MC_TrueNu_2019 " << LH_MC_TrueNu_2019.size()
222  << "\n\t LH_MC_TrNuMu_2019 " << LH_MC_TrNuMu_2019.size()
223  << "\n\t LH_MC_AnNuMu_2019 " << LH_MC_AnNuMu_2019.size()
224  << "\n\t LH_Data_2020 " << LH_Data_2020 .size()
225  << "\n\t LH_Data_2019 " << LH_Data_2019 .size()
226  << std::endl;
227 
228  // --- Now that I've loaded everything, I need to make my plots!
229  /*
230  std::cout << "\n\n Now to plot things..." << std::endl;
231  for (size_t CutInd=0; CutInd<LH_MC_NoTrue_2020.size(); ++CutInd) {
232 
233  // Set the 2020 Cut Factor
234  unsigned int cut = CutInd*4;
235  // Set the 2019 Cut Factor.
236  unsigned int cut19 = CutInd *4, CutInd19 = CutInd;
237  if (CutInd==2){cut19 = (CutInd-1)*4; CutInd19 = CutInd-1;}
238  // Check that the elements are set correctly.
239 
240  //std::cout << "\t What elements should I plot?"
241  // << "\n\t\t 2020 MC: "<<My2020_Mont[cut ].second<<", "<<My2020_Mont[cut+1 ].second<<", "<<My2020_Mont[cut +2].second<<", "<<My2020_Mont[cut +3].second
242  // << "\n\t\t 2020 Data: "<<My2020_Data[CutInd ].second
243  // << "\n\t\t 2019 MC: "<<My2019_Mont[cut19 ].second<<", "<<My2019_Mont[cut19+1].second<<", "<<My2019_Mont[cut19+2].second<<", "<<My2019_Mont[cut19+3].second
244  // << "\n\t\t 2019 Data: "<<My2019_Data[CutInd19].second
245  // << std::endl;
246 
247  // --- Now to figure out how to make my plots...
248  // Want the "CutInd" elements of each LH vector, and the "cut" element of the CutName vector
249  ProducePlots_2020( LH_MC_NoTrue_2020[ CutInd ], // 2020 - NoTrue
250  LH_MC_TrueNu_2020[ CutInd ], // 2020 - TrueNu
251  LH_MC_TrNuMu_2020[ CutInd ], // 2020 - TrNuMu
252  LH_MC_AnNuMu_2020[ CutInd ], // 2020 - AnNuMu
253  LH_MC_NoTrue_2019[ CutInd19 ], // 2019 - NoTrue
254  LH_Data_2020 [ CutInd ], // 2020 - Data
255  My2020_Mont[cut].second ); // CutTier.
256  if (CutInd!=2) {
257  ProducePlots_2019( LH_MC_NoTrue_2019[ CutInd19 ], // 2019 - NoTrue
258  LH_MC_TrueNu_2019[ CutInd19 ], // 2019 - TrueNu
259  LH_MC_TrNuMu_2019[ CutInd19 ], // 2019 - TrNuMu
260  LH_MC_AnNuMu_2019[ CutInd19 ], // 2019 - AnNuMu
261  LH_MC_NoTrue_2020[ CutInd ], // 2020 - NoTrue
262  LH_Data_2019 [ CutInd19 ], // 2019 - Data
263  My2019_Mont[cut19].second ); // CutTier.
264  }
265  }
266  */
267  // --- Finally, return.
268  return;
269 } // DoThePlots
270 
271 //......................................................
272 LoadedHistograms_2020 LoadFile_GetHists_2020( TFile* InFile, std::string CutTier, bool isData ) {
273 
274  // --- Set my Y axis
275  std::string YAxisStr = "Number of Events";
276 
277  // --- Declare the vector of LoadedHistograms_2020 which I am going to return...
278  LoadedHistograms_2020 MyHists;
279  // --- Grab my spectra...
280  MyHists.hNuRecoE = GetSpectToHist( InFile, "RecoNuE" +CutTier, ";Reco Nu E;" +YAxisStr+"/0.1 GeV", isData);
281  /*
282  MyHists.hReMIdScore = GetSpectToHist( InFile, "ReMIdScore" +CutTier, ";ReMId Score;" +YAxisStr, isData);
283  MyHists.hProngScoremuon= GetSpectToHist( InFile, "ProngCVNMuonScore" +CutTier, ";Prong CVN Muon Score;" +YAxisStr, isData);
284  MyHists.hCosRejScoreP4 = GetSpectToHist( InFile, "Prod4CosRejScore" +CutTier, ";Prod4 CosRej Score;" +YAxisStr, isData);
285  MyHists.hCosRejScoreP5 = GetSpectToHist( InFile, "Prod5CosRejScore" +CutTier, ";Prod5 CosRej Score;" +YAxisStr, isData);
286  MyHists.hOldPreselmuon = GetSpectToHist( InFile, "OldPreselMuonScore"+CutTier, ";Old Presel Muon Score;" +YAxisStr, isData);
287  MyHists.hOldPreselcosm = GetSpectToHist( InFile, "OldPreselCosmScore"+CutTier, ";Old Presel Cosmic Score;"+YAxisStr, isData);
288  MyHists.hLoosePTPmuon = GetSpectToHist( InFile, "LoosePtpMuonScore" +CutTier, ";Loose Pt/p Muon Score;" +YAxisStr, isData);
289  MyHists.hLoosePTPcosm = GetSpectToHist( InFile, "LoosePtpCosmScore" +CutTier, ";Loose Pt/p Cosmic Score;"+YAxisStr, isData);
290  */
291  // --- Finally, return MyHists.
292  return MyHists;
293 }
294 
295 //......................................................
296 LoadedHistograms_2019 LoadFile_GetHists_2019( TFile* InFile, std::string CutTier, bool isData ) {
297 
298  // --- Set my Y axis
299  std::string YAxisStr = "Number of Events";
300 
301  // --- Declare the vector of LoadedHistograms_2020 which I am going to return...
302  LoadedHistograms_2019 MyHists;
303  // --- Grab my spectra...
304  TH1D* hCosRejScoreP4 ;
305  TH1D* hCVNmuon ;
306  MyHists.hNuRecoE = GetSpectToHist( InFile, "RecoNuE" +CutTier, ";Reco Nu E;" +YAxisStr+"/0.1 GeV", isData);
307  /*
308  MyHists.hReMIdScore = GetSpectToHist( InFile, "ReMIdScore" +CutTier, ";ReMId Score;" +YAxisStr, isData);
309  MyHists.hProngScoremuon= GetSpectToHist( InFile, "ProngCVNMuonScore" +CutTier, ";Prong CVN Muon Score;" +YAxisStr, isData);
310  MyHists.hCosRejScoreP4 = GetSpectToHist( InFile, "Prod4CosRejScore" +CutTier, ";Prod4 CosRej Score;" +YAxisStr, isData);
311  MyHists.hCVNmuon = GetSpectToHist( InFile, "CVNMuonScore" +CutTier, ";Prod3Train Muon Score;" +YAxisStr, isData);
312  */
313  // --- Finally, return MyHists.
314  return MyHists;
315 }
316 
317 //......................................................
318 TH1D* GetSpectToHist( TFile* MyF, std::string LoadName, std::string Axis, bool isData ) {
319  MyF->cd();
320  // --- Declare my histogram.
321  TH1D* MyHist;
322  // --- Am I looking at the NoExtrap MC?
323  if ( !isData && bIsFD ) {
324  std::unique_ptr<PredictionNoExtrap> TempSpec = PredictionNoExtrap::LoadFrom( MyF, LoadName.c_str() ) ;
325 
326  // --- Set my Oscillation Cal
329  calc.SetdCP ( 0 );
330  calc.SetDmsq32( 0.00248 );
331  calc.SetTh23 ( asin(sqrt(0.565)) );
332  MyHist = TempSpec -> Predict(&calc).ToTH1( POTNorm );
333  //MyHist = TempSpec -> PredictUnoscillated().ToTH1( POTNorm ); std::cout << "\n\t\t!!!UNOSCILLATED SPECTRUM!!!" << std::endl;
334  }
335  // ---- If just looking at a spectra.
336  else {
337  std::unique_ptr<Spectrum> TempSpec = Spectrum::LoadFrom( MyF, LoadName.c_str() ) ;
338  if (bIsFD) MyHist = TempSpec -> ToTH1( LivNorm, kLivetime );
339  else MyHist = TempSpec -> ToTH1( POTNorm );
340  }
341 
342  // --- Scale the RecoNuE plots.
343  if (LoadName.find("RecoNuE") != std::string::npos) {
344  //MyHist -> Scale( 0.1, "width" );
345  }
346 
347  // --- Make my plots pretty, and return.
348  CenterTitles ( MyHist);
349  MyHist->SetMinimum( 0 );
350  MyHist->SetMaximum( MyHist->GetMaximum()*1.05 );
351  MyHist->SetTitle ( Axis.c_str() );
352 
353  std::cout << "\tLoadName is " << LoadName << ", integral is " << MyHist->Integral() << std::endl;
354 
355  return MyHist;
356 }
357 
358 //......................................................
359 TH1D* MakeRatio( TH1D* num, TH1D* denom, int Col, std::string PlotName ) {
360  // Declare a new histogram.
361  TH1D* ratHist = new TH1D( TString(PlotName)+TString("_ratio_")+TString(std::to_string(Col)),
362  TString(";")+TString(num->GetXaxis()->GetTitle())+TString(";Ratio"),
363  num->GetNbinsX(),
364  num->GetXaxis()->GetXmin(),
365  num->GetXaxis()->GetXmax()
366  );
367  // Don't want any exponents?
368  ratHist->GetYaxis()->SetNoExponent();
369  ratHist->GetXaxis()->SetNoExponent();
370  // Divide the two input histograms...
371  ratHist->Divide( num , denom );
372  // Set Line and Marker colors / styles.
373  ratHist->SetLineColor ( Col );
374  ratHist->SetMarkerColor( Col );
375  ratHist->SetMarkerStyle( 20 );
376  ratHist->SetMarkerSize ( 0 );
377  // Figure out what the range should be...
378  ratHist->GetYaxis()->SetRangeUser( 0.8, 1.2 );
379  // Return my shiny new ratio plot!
380  CenterTitles(ratHist);
381  return ratHist;
382 }
383 
384 //......................................................
385 void SetRange( TH1D* MC, TH1D* Data, // Pass the histograms
386  double XLow , double XHigh, // What X values do I want?
387  double H_YLow, double H_YHigh ) { // What Y values do I want?
388  // --- First off, set the X ranges.
389  if ( (XLow != -1) && (XHigh != -1) ) {
390  MC -> GetXaxis() -> SetRangeUser( XLow, XHigh );
391  if (!bIsFD) Data -> GetXaxis() -> SetRangeUser( XLow, XHigh );
392  }
393 
394  // --- And, the Y axes, but check that they aren't -1!
395  // Hists
396  if ( (H_YLow != -1) && (H_YHigh != -1) ) {
397  MC -> GetYaxis() -> SetRangeUser( H_YLow, H_YHigh );
398  if (!bIsFD) Data -> GetYaxis() -> SetRangeUser( H_YLow, H_YHigh );
399  }
400  // --- Return
401  return;
402 }
403 
404 //......................................................
405 void FindAxisRange( double &XLow, double &XHigh, double &H_YLow, double &H_YHigh, double &R_YLow, double &R_YHigh, double MaxVal, bool &SetLogy, std::string Name ) {
406 
407  XLow = -1; XHigh = -1;
408  H_YLow = 0; H_YHigh = MaxVal;
409  R_YLow = 0.5; R_YHigh = 1.5;
410  SetLogy = false;
411  if (!bIsFD) {
412  SetLogy = true;
413  H_YLow = 10;
414  }
415 
416  bool IsQuant = ( (Name.find("AllQuants") != std::string::npos) == true ? false:true);
417 
418  if ( (Name.find("Muon") != std::string::npos) ) {
419  if ( bIsFD && bIsFHC && Name.find("QualCont") != std::string::npos) { H_YLow = 0; H_YHigh = 200; }
420  if ( bIsFD && !bIsFHC && Name.find("QualCont") != std::string::npos) { H_YLow = 0; H_YHigh = 120; }
421  if ( bIsFD && bIsFHC && Name.find("Full") != std::string::npos) { H_YLow = 0; H_YHigh = 100; }
422  if ( bIsFD && !bIsFHC && Name.find("Full") != std::string::npos) { H_YLow = 0; H_YHigh = 100; }
423  }
424 
425  if ( (Name.find("RecoE") != std::string::npos) ) {
426  if ( bIsFD && bIsFHC && Name.find("QualCont") != std::string::npos) { H_YLow = 0; H_YHigh = 30; }
427  if ( bIsFD && !bIsFHC && Name.find("QualCont") != std::string::npos) { H_YLow = 0; H_YHigh = 20; }
428  if ( bIsFD && bIsFHC && Name.find("Full") != std::string::npos) { H_YLow = 0; H_YHigh = 15; }
429  if ( bIsFD && !bIsFHC && Name.find("Full") != std::string::npos) { H_YLow = 0; H_YHigh = 10; }
430 }
431 
432  if ( (Name.find("ReMId") != std::string::npos) ) {
433  if ( bIsFD && bIsFHC ) { H_YLow = 0; H_YHigh = 100; }
434  if ( bIsFD && !bIsFHC ) { H_YLow = 0; H_YHigh = 80; }
435  }
436 
437  /*
438  std::cout << "\t\tParams are; Xlow " << XLow << ", XHigh: " << XHigh << ", H_YLow: " << H_YLow << ", H_YHigh: " << H_YHigh
439  << ", R_YLow: " << R_YLow << ", R_YHigh " << R_YHigh << ", SetLogy " << SetLogy
440  << std::endl;
441  */
442  return;
443 }
444 
445 //......................................................
446 void PlotHistProp( TH1D* Hist, TLegend* L ) {
447  double Mean_Nom = Hist->GetMean();
448  double RMS_Nom = Hist->GetRMS();
449  double Int_Nom = Hist->Integral();
450  L -> AddEntry( Hist, Form("Mean: %.3f", Mean_Nom), "" );
451  L -> AddEntry( Hist, Form("RMS: %.3f" , RMS_Nom ), "" );
452  L -> AddEntry( Hist, Form("Int: %.3f" , Int_Nom ), "" );
453  return;
454 }
455 
456 //......................................................
457 void MakeSplitCans( TH1D* NoTrue, TH1D* TrueNu, TH1D* TrNuMu, TH1D* AnNuMu, TH1D* NoTrueAlt, TH1D* Data, std::string PlotName ) {
458  // ------ What histograms am I passing to the function ------
459  // NoTrue --> The Monte Carlo histogram with no truth cuts applied
460  // TrueNu --> The Monte Carlo histogram with a true neutrino cut applied
461  // TrNuMu --> The Monte Carlo histogram with a true muon neutrino CC cut applied
462  // AnNuMu --> The Monte Carlo histogram with a true anti-muon neutrino CC cut applied
463  // NoTrueAlt --> The Monte Carlo histogram with no truth cuts applied [ For the other Prod file ]
464  // Data --> The Data histogram [ For FD this is empty ].
465  // ----------------------------------------------------------
466  std::cout << "---> In MakeSplitCans .... " << PlotName << std::endl;
467 
468  // --- Find some integrals, mainly for the RecoNuE plots.
469  double Int_NoTrue, Int_TrueNu, Int_TrNuMu, Int_AnNuMu, Int_NoTrueAlt;
470  double Dip_NoTrue, Dip_TrueNu, Dip_TrNuMu, Dip_AnNuMu, Dip_NoTrueAlt;
471  Int_NoTrue = Int_TrueNu = Int_TrNuMu = Int_AnNuMu = Int_NoTrueAlt = 0;
472  Dip_NoTrue = Dip_TrueNu = Dip_TrNuMu = Dip_AnNuMu = Dip_NoTrueAlt = 0;
473  int MinBin, MaxBin;
474  MinBin = MaxBin = 0;
475  // If RecoNuE do the integrals.
476  if ((PlotName.find("RecoE") != std::string::npos) ) {
477  Int_NoTrue = NoTrue ->Integral();
478  Int_TrueNu = TrueNu ->Integral();
479  Int_TrNuMu = TrNuMu ->Integral();
480  Int_AnNuMu = AnNuMu ->Integral();
481  Int_NoTrueAlt = NoTrueAlt->Integral();
482 
483  MinBin = NoTrue -> GetXaxis() -> FindBin( 1.0 );
484  MaxBin = NoTrue -> GetXaxis() -> FindBin( 2.0 );
485  std::cout << "\t MinBin = " << MinBin << ". \t MaxBin = " << MaxBin << std::endl;
486 
487  Dip_NoTrue = NoTrue ->Integral( MinBin, MaxBin );
488  Dip_TrueNu = TrueNu ->Integral( MinBin, MaxBin );
489  Dip_TrNuMu = TrNuMu ->Integral( MinBin, MaxBin );
490  Dip_AnNuMu = AnNuMu ->Integral( MinBin, MaxBin );
491  Dip_NoTrueAlt = NoTrueAlt->Integral( MinBin, MaxBin );
492 
493  std::cout << "\n\t NoTrue Full: " << Int_NoTrue << ", Dip: " << Dip_NoTrue
494  << "\n\t TrueNu Full: " << Int_TrueNu << ", Dip: " << Dip_TrueNu
495  << "\n\t TrNuMu Full: " << Int_TrNuMu << ", Dip: " << Dip_TrNuMu
496  << "\n\t AnNuMu Full: " << Int_AnNuMu << ", Dip: " << Dip_AnNuMu
497  << "\n\t NoTrueAlt Full: " << Int_NoTrueAlt << ", Dip: " << Dip_NoTrueAlt
498  << std::endl;
499 
500  double FOM_Full = Int_TrNuMu / sqrt( Int_NoTrue );
501  double FOM_Dip = Dip_TrNuMu / sqrt( Dip_NoTrue );
502  if (!bIsFHC) {
503  FOM_Full = Int_AnNuMu / sqrt( Int_NoTrue );
504  FOM_Dip = Dip_AnNuMu / sqrt( Dip_NoTrue );
505  }
506  std::cout << "\t FOM_Full = " << FOM_Full << ". \n\t FOM_Dip = " << FOM_Dip << std::endl;
507  }
508 
509  // Add the wrong sign NuMu to the correct sign one.
510  if (bIsFHC) TrNuMu->Add( AnNuMu );
511  else AnNuMu->Add( TrNuMu );
512 
513  // --- Set my histogram colours...
514  NoTrue -> SetLineColor( 2 ); NoTrue -> SetFillColor( 2 ); NoTrue -> SetFillStyle( 3003 );
515  //TrueNu -> SetLineColor( 2 ); TrueNu -> SetFillColor( 2 ); TrueNu -> SetFillStyle( 1001 );
516  TrNuMu -> SetLineColor( 4 ); TrNuMu -> SetFillColor( 4 ); TrNuMu -> SetFillStyle( 1001 );
517  AnNuMu -> SetLineColor( 6 ); AnNuMu -> SetFillColor( 6 ); AnNuMu -> SetFillStyle( 1001 );
518  NoTrueAlt-> SetLineColor( 1 ); NoTrueAlt -> SetMarkerColor( 1 ); NoTrueAlt -> SetMarkerStyle( 20 );
519  if (!bIsFD) { Data -> SetMarkerColor( 1 ); Data -> SetMarkerStyle( 3 ); }
520 
521  // --- What axis ranges do I want to use?
522  double XL, XH, H_YL, H_YH, R_YL, R_YH;
523  bool SetLogy;
524  FindAxisRange( XL, XH, H_YL, H_YH, R_YL, R_YH, 1.5*NoTrue->GetMaximum(), SetLogy, PlotName );
525  // --- And set the ranges...
526  SetRange( NoTrue, Data, XL, XH, H_YL, H_YH );
527 
528  // --- Now, let's make a histogram with just the plot.
529  std::string HistCanNa = PlotName;
530  TCanvas* HistCan = new TCanvas( HistCanNa.c_str(), HistCanNa.c_str() );
531 
532  // --- And now start drawing things!!!
533  NoTrue -> Draw("hist");
534  //TrueNu -> Draw("same");
535  if (bIsFHC) { TrNuMu -> Draw("same"); AnNuMu -> Draw("same"); }
536  else { AnNuMu -> Draw("same"); TrNuMu -> Draw("same"); }
537  NoTrueAlt -> Draw("same");
538  if (!bIsFD) Data -> Draw( "same" );
539 
540  // --- Add whether FHC or RHC
541  if (bIsFHC ) CornerLabel( "Neutrino beam" );
542  else CornerLabel( "Antineutrino beam" );
543 
544  // Do I want the canvas to be loagarithmic?
545  if ( SetLogy ) {
546  TrNuMu -> GetYaxis() -> SetRangeUser( H_YL, H_YH );
547  AnNuMu -> GetYaxis() -> SetRangeUser( H_YL, H_YH );
548  NoTrueAlt -> GetYaxis() -> SetRangeUser( H_YL, H_YH );
549  HistCan -> SetLogy();
550  }
551 
552  // --- Where do I want my legend...
553  double HistLeg_X1 = 0.45, HistLeg_Y1 = 0.55, HistLeg_X2 = 0.85, HistLeg_Y2 = 0.88;
554  // Do I want the legend on the left instead?
555  if ( (PlotName.find("CosRej") != std::string::npos) ) {
556  HistLeg_X1 = 0.10;
557  HistLeg_X2 = 0.56;
558  }
559  // --- Make the legend where I've decided...
560  TLegend *HistLeg = new TLegend( HistLeg_X1, HistLeg_Y1, HistLeg_X2, HistLeg_Y2 );
561  HistLeg -> SetFillStyle (0);
562  HistLeg -> SetBorderSize(0);
563  HistLeg -> SetTextSize (NoTrue ->GetXaxis()->GetTitleSize() );
564  // And now fill it.
565  HistLeg -> AddEntry( NoTrue , "Monte Carlo");
566  //HistLeg -> AddEntry( TrueNu , "True Nu Cut");
567  if (bIsFHC) { HistLeg -> AddEntry( TrNuMu, "NuMu CC Cut" ); HistLeg -> AddEntry( AnNuMu, "Anti-NuMu CC Cut"); }
568  else { HistLeg -> AddEntry( AnNuMu, "Anti-NuMu CC Cut"); HistLeg -> AddEntry( TrNuMu, "NuMu CC Cut" ); }
569  HistLeg -> AddEntry( NoTrueAlt , "Alt Monte Carlo");
570  if (!bIsFD) { HistLeg -> AddEntry( Data, "Data"); }
571  // And draw the legend
572  HistLeg -> Draw();
573 
574  // --- Now I can save my things...
575  HistCan -> SaveAs( TString(OutDir)+TString(HistCan->GetName())+TString(".pdf") );
576  //HistCan -> SaveAs( TString(OutDir)+TString(HistCan->GetName())+TString(".png") );
577  OutFile -> cd();
578  HistCan -> Write( HistCan->GetName() );
579 
580 
581  // -----------------------------------------------------------
582  // --- Now for the ratio part of the plots....
583  // --- Provided that I'm looking at the ND...
584  // -----------------------------------------------------------
585  /*
586  if (!bIsFD) {
587  // Make my ratio plot
588  TH1D* Ratio = MakeRatio( MCHist, Data, 2, PlotName );
589 
590  std::string RatCanNa = PlotName+"_ratio";
591  TCanvas *RatCan = new TCanvas( RatCanNa.c_str(), RatCanNa.c_str() );
592  // --- If looking at cosmics then I want the plot to be log y
593  Ratio -> GetYaxis() -> SetTitle("Ratio");
594  Ratio -> Draw();
595  // --- Add a reference line
596  double x[2] = { Ratio->GetXaxis()->GetXmin(), Ratio->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  Ratio -> Draw("axis same");
605 
606  // ------------------------------
607  // --- Now merge the two canvases
608  std::string MergeCanNa = PlotName+"_merged";
609  TCanvas *MergeCan = new TCanvas( MergeCanNa.c_str(), MergeCanNa.c_str() );
610  MergeCan->Divide(1,2);
611  HistCan -> SetRightMargin(.03); HistCan -> SetLeftMargin(.07); HistCan -> SetTopMargin(.10); HistCan -> SetBottomMargin(.00);
612  RatCan -> SetRightMargin(.03); RatCan -> SetLeftMargin(.07); RatCan -> SetTopMargin(.00); RatCan -> SetBottomMargin(.12);
613 
614  // And now the axes...
615  double Lb1 = 0.07;
616  MCHist -> GetYaxis() -> SetTitleSize( Lb1 ); MCHist -> GetYaxis() -> SetLabelSize( Lb1 ); MCHist -> GetYaxis() -> SetTitleOffset( 0.45 );
617  Ratio -> GetYaxis() -> SetTitleSize( Lb1 ); Ratio -> GetYaxis() -> SetLabelSize( Lb1 ); Ratio -> GetYaxis() -> SetTitleOffset( 0.45 );
618  Ratio -> GetXaxis() -> SetTitleSize( Lb1 ); Ratio -> GetXaxis() -> SetLabelSize( Lb1 );
619 
620  MergeCan->cd(1);
621  HistCan->DrawClonePad();
622  MergeCan->cd(2);
623  RatCan ->DrawClonePad();
624 
625  // Finally, save the merged canvas.
626  MergeCan -> SaveAs( TString(OutDir)+TString(MergeCanNa)+TString(".pdf" ) );
627  //MergeCan -> SaveAs( TString(OutDir)+TString(MergeCanNa)+TString(".png" ) );
628  OutFile -> cd();
629  MergeCan -> Write ( MergeCanNa.c_str() );
630  }
631  */
632  // --- Return.
633  return;
634 }
635 
636 //.....................................................
638  LoadedHistograms_2020 TrNuMu_2020, LoadedHistograms_2020 AnNuMu_2020,
639  LoadedHistograms_2019 NoTrue_2019, LoadedHistograms_2020 Data_2020 , std::string CutTier ) {
640 
641  // --- Make my output directory.
642  OutDir = "Plots-2020-"+sIsFD+"-"+sIsFHC+"/";
643  std::string OutName = OutDir+"BasicPIDPlots-"+sIsFD+"-"+sIsFHC+".root";
644  gSystem -> MakeDirectory( OutDir.c_str() );
645  OutFile = TFile::Open( OutName.c_str(), "RECREATE" );
646 
647  // --- Make a quick string to combine some important nomenclature
648  std::string EndCan = sIsFD+"_"+sIsFHC+CutTier+"_2020";
649 
650  std::cout << "\n In ProducePlots_2020. I passed the cut " << CutTier << std::endl;
651 
652  // --- Now let's start making all of the plots!!!!
653  // -----------------------------------------------
654  MakeSplitCans( NoTrue_2020.hNuRecoE , TrueNu_2020.hNuRecoE , TrNuMu_2020.hNuRecoE , AnNuMu_2020.hNuRecoE ,
655  NoTrue_2019.hNuRecoE , Data_2020 .hNuRecoE , "hNuRecoE_" +EndCan );
656  /*
657  MakeSplitCans( NoTrue_2020.hReMIdScore , TrueNu_2020.hReMIdScore , TrNuMu_2020.hReMIdScore , AnNuMu_2020.hReMIdScore ,
658  NoTrue_2019.hReMIdScore , Data_2020 .hReMIdScore , "hReMIdScore_" +EndCan );
659 
660  MakeSplitCans( NoTrue_2020.hProngScoremuon, TrueNu_2020.hProngScoremuon, TrNuMu_2020.hProngScoremuon, AnNuMu_2020.hProngScoremuon,
661  NoTrue_2019.hProngScoremuon, Data_2020 .hProngScoremuon, "hProngCVNMuonScore_" +EndCan );
662 
663  MakeSplitCans( NoTrue_2020.hOldPreselmuon , TrueNu_2020.hOldPreselmuon , TrNuMu_2020.hOldPreselmuon , AnNuMu_2020.hOldPreselmuon ,
664  NoTrue_2019.hCVNmuon , Data_2020 .hOldPreselmuon , "hOldPreselMuonScore_"+EndCan );
665 
666  MakeSplitCans( NoTrue_2020.hLoosePTPmuon , TrueNu_2020.hLoosePTPmuon , TrNuMu_2020.hLoosePTPmuon , AnNuMu_2020.hLoosePTPmuon ,
667  NoTrue_2019.hCVNmuon , Data_2020 .hLoosePTPmuon , "hLoosePtpMuonScore_" +EndCan );
668  if (bIsFD) {
669  // CosRej
670  MakeSplitCans( NoTrue_2020.hCosRejScoreP4 , TrueNu_2020.hCosRejScoreP4 , TrNuMu_2020.hCosRejScoreP4 , AnNuMu_2020.hCosRejScoreP4 ,
671  NoTrue_2019.hCosRejScoreP4 , Data_2020 .hCosRejScoreP4 , "hProd4CosRejScore_" +EndCan );
672 
673  MakeSplitCans( NoTrue_2020.hCosRejScoreP5 , TrueNu_2020.hCosRejScoreP5 , TrNuMu_2020.hCosRejScoreP5 , AnNuMu_2020.hCosRejScoreP5 ,
674  NoTrue_2019.hCosRejScoreP4 , Data_2020 .hCosRejScoreP5 , "hProd5CosRejScore_" +EndCan );
675  // Cosmic Score
676  MakeSplitCans( NoTrue_2020.hOldPreselcosm , TrueNu_2020.hOldPreselcosm , TrNuMu_2020.hOldPreselcosm , AnNuMu_2020.hOldPreselcosm ,
677  NoTrue_2019.hCVNmuon , Data_2020 .hOldPreselcosm , "hOldPreselCosmScore_"+EndCan );
678 
679  MakeSplitCans( NoTrue_2020.hLoosePTPcosm , TrueNu_2020.hLoosePTPcosm , TrNuMu_2020.hLoosePTPcosm , AnNuMu_2020.hLoosePTPcosm ,
680  NoTrue_2019.hCVNmuon , Data_2020 .hLoosePTPcosm , "hLoosePtpCosmScore_" +EndCan );
681  }
682  */
683  return;
684 }
685 
686 //.....................................................
688  LoadedHistograms_2019 TrNuMu_2019, LoadedHistograms_2019 AnNuMu_2019,
689  LoadedHistograms_2020 NoTrue_2020, LoadedHistograms_2019 Data_2019 , std::string CutTier ) {
690 
691  // --- Make my output directory.
692  OutDir = "Plots-2019-"+sIsFD+"-"+sIsFHC+"/";
693  std::string OutName = OutDir+"BasicPIDPlots-"+sIsFD+"-"+sIsFHC+".root";
694  gSystem -> MakeDirectory( OutDir.c_str() );
695  OutFile = TFile::Open( OutName.c_str(), "RECREATE" );
696 
697  // --- Make a quick string to combine some important nomenclature
698  std::string EndCan = sIsFD+"_"+sIsFHC+CutTier+"_2019";
699 
700  std::cout << "\n In ProducePlots_2019. I passed the cut " << CutTier << std::endl;
701 
702  // --- Now let's start making all of the plots!!!!
703  // -----------------------------------------------
704  MakeSplitCans( NoTrue_2019.hNuRecoE , TrueNu_2019.hNuRecoE , TrNuMu_2019.hNuRecoE , AnNuMu_2019.hNuRecoE ,
705  NoTrue_2020.hNuRecoE , Data_2019.hNuRecoE , "hNuRecoE_" +EndCan );
706 
707  return;
708 }
709 
710 //......................................................
tree Draw("slc.nhit")
bin1_2sigma SetFillColor(3)
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
std::string sIsFD
correl_xv GetYaxis() -> SetDecimals()
void SetRange(TH1D *MC, TH1D *Data, double XLow, double XHigh, double H_YLow, double H_YHigh)
const double kAna2019RHCLivetime
Definition: Exposures.h:227
TH1D * GetSpectToHist(TFile *MyF, std::string LoadName, std::string Axis, bool isData)
T sqrt(T number)
Definition: d0nt_math.hpp:156
prelim SetTextSize(2/30.)
void FindAxisRange(double &XLow, double &XHigh, double &H_YLow, double &H_YHigh, double &R_YLow, double &R_YHigh, double MaxVal, bool &SetLogy, std::string Name)
void CornerLabel(std::string &s)
Definition: numu_tools.h:145
std::vector< std::pair< Cut, std::string > > BasicPIDPlots_Cuts(bool isFHC, bool isFD, bool isData)
Definition: Hist.h:29
bool isfd
Definition: geo.cxx:18
void CenterTitles(TH1 *histo)
Definition: Plots.cxx:1481
fVtxDx SetMarkerStyle(20)
osc::OscCalcDumb calc
ntuple SetFillStyle(1001)
void NuMu2020_BasicPIDPlots_Plot(bool isfhc, bool isfd, std::string period="full")
double LivNorm
legend SetBorderSize(0)
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:535
static constexpr double L
std::string sIsFHC
correl_xv GetXaxis() -> SetDecimals()
std::string OutDir
LoadedHistograms_2020 LoadFile_GetHists_2020(TFile *InFile, std::string CutTier, bool isData)
leg AddEntry(GRdata,"data","p")
hmean SetLineColor(4)
Optimized version of OscCalcPMNS.
Definition: StanTypedefs.h:31
h1 SetMarkerColor(4)
void ProducePlots_2019(LoadedHistograms_2019 NoTrue_2019, LoadedHistograms_2019 TrueNu_2019, LoadedHistograms_2019 TrNuMu_2019, LoadedHistograms_2019 AnNuMu_2019, LoadedHistograms_2020 NoTrue_2020, LoadedHistograms_2019 Data_2019, std::string CutTier)
double POTNorm
const double kAna2019FHCLivetime
Definition: Exposures.h:226
void PlotHistProp(TH1D *Hist, TLegend *L)
OStream cout
Definition: OStream.cxx:6
const double kAna2019RHCPOT
Definition: Exposures.h:224
const Cut cut
Definition: exporter_fd.C:30
LoadedHistograms_2019 LoadFile_GetHists_2019(TFile *InFile, std::string CutTier, bool isData)
int num
Definition: f2_nu.C:119
const double kAna2019FHCPOT
Definition: Exposures.h:223
TFile * OutFile
TH1D * MakeRatio(TH1D *num, TH1D *denom, int Col, std::string FType)
void MakeSplitCans(TH1D *MCHist, TH1D *DataHist, std::string PlotName)
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
cosmicTree SaveAs("cosmicTree.root")
static std::unique_ptr< PredictionNoExtrap > LoadFrom(TDirectory *dir, const std::string &name)
c cd(1)
void ProducePlots_2020(LoadedHistograms_2020 NoTrue_2020, LoadedHistograms_2020 TrueNu_2020, LoadedHistograms_2020 TrNuMu_2020, LoadedHistograms_2020 AnNuMu_2020, LoadedHistograms_2019 NoTrue_2019, LoadedHistograms_2020 Data_2020, std::string CutTier)
std::vector< std::pair< Cut, std::string > > BasicPIDPlots2019_Cuts(bool isFHC, bool isFD, bool isData)
gPad SetLogy()
T asin(T number)
Definition: d0nt_math.hpp:60
gm Write()
enum BeamMode string