MakeThePlots.C
Go to the documentation of this file.
1 #include "TCanvas.h"
2 #include "TFile.h"
3 #include "TGraph.h"
4 #include "TGraphAsymmErrors.h"
5 #include "TH1D.h"
6 #include "TLine.h"
7 #include "TLatex.h"
8 #include "TLegend.h"
9 #include "TStyle.h"
10 #include "TSystem.h"
11 
13 #include "CAFAna/Analysis/Plots.h"
14 
15 #include "CAFAna/Core/Binning.h"
16 #include "CAFAna/Core/Spectrum.h"
18 
19 #include "OscLib/OscCalcPMNSOpt.h"
20 #include "CAFAna/Analysis/Calcs.h"
22 
23 #include "Utilities/rootlogon.C"
24 
25 #include "TFeldmanCousins.h"
26 
27 #include <iostream>
28 #include <sstream>
29 #include <fstream>
30 #include <vector>
31 #include <string>
32 
33 #include <TROOT.h>
34 #include <TStyle.h>
35 
36 using namespace ana;
37 
38 //......................................................
39 struct LoadedHistograms {
40  // *** Some key ones first.
41  TH1D* hReconstEnergy;
42  TH1D* hMuonEnergy ;
43  TH1D* hMuonEnPerHit ;
44  TH1D* hHadronEnergy ;
45  TH1D* hHadroEnPerHit;
46  TH1D* hHadFracEnergy;
47  TH1D* hHitsPerSlice ;
48  //*
49  TH1D* hSliceTimeFull;
50  TH1D* hSliceTimeZoom;
51  TH1D* hRemIDScore ;
52  TH1D* hCVNCosmicScor;
53  TH1D* hCVNNuMuIDScor;
55  TH1D* hNuMuContPID ;
56  TH1D* hSANuMuContPID;
57  // *** Some important Truth based ones.
58  TH1D* hTrueNuEnergy ;
59  TH1D* hReTrOverTrNuE;
60 
61  // *** Basic positional plots.
62  TH1D* hTrkStartX;
63  TH1D* hTrkStartY;
64  TH1D* hTrkStartZ;
65  TH1D* hTrkEndX ;
66  TH1D* hTrkEndY ;
67  TH1D* hTrkEndZ ;
68  TH1D* hTrkLenXY ;
69 
70  // *** Track based ones...
71  TH1D* hNumKalTracks ;
72  TH1D* hNumKalTrHits ;
73  TH1D* hRatKalHitSlc ;
74  TH1D* hKalTrLength ;
75  TH1D* hKalTrBeamAng ;
80  TH1D* hKalTrDir_Y ;
82 
83  // *** Hit based ones...
84  TH1D* hFirstHitCell ;
85  TH1D* hLastHitCell ;
91 
92  // *** RemID score / inputs
93  TH1D* hRemIDScatLLH ;
94  TH1D* hRemIDdEdxLLH ;
96 
97  // *** Some calibration quantities
98  TH1D* hCor_RecEn[4];
99  TH1D* hCor_MuoEn[4];
100  TH1D* hCor_HadEn[4];
101  TH1D* hCor_HFrEn[4];
102  //*/
103 };
104 
105 //......................................................
106 double POTNom = 1;
107 double LivNom = 1;
108 bool IsFHC = true;
109 bool AllPer = true;
110 std::string sFHC = "fhc" ;
111 std::string sPer = "full";
112 // --- Cuts
113 const int NumCut = 4;
115  "_QualCuts",
116  "_DetCuts" ,
117  "_CosmCuts",
118  "_AllCuts"
119 };
120 // --- Quantiles
121 const int NumQuant = 5;
123  "_AllQuants",
124  "_Quant1",
125  "_Quant2",
126  "_Quant3",
127  "_Quant4"
128 };
129 // --- Some variables that are global, but change.
130 unsigned int NPlot = 5;
131 // --- What colours, and legend names do I want?
132 const int Colours [5] = { kBlack , kViolet-5 , kGreen-6 , kGray+1 , kAzure+1 };
133 std::string LegNames[5] = { "FD data", "Best Fit Pred", "Wrong Sign", "Total Bkg.", "Cosmic Bkg." };
134 
135 // --- Some other things.
136 const bool SingleCut = true;
137 const bool SingleQuant = false;
138 const bool PlotMean = false;
139 const bool PlotTest = false;
140 const bool SaveRatio = false;
141 const bool SaveHists = true;
142 
143 // --- Calibration labels.
144 const std::string CornName[4] = { "TopL", "TopR", "BotL", "BotR" };
145 
146 //......................................................
147 LoadedHistograms LoadFile_GetHists( TFile* InFile, std::string CTier, int FileType );
148 //......................................................
149 TH1D* GetSpectToHist( TFile* MyF, std::string LoadName, std::string Axis, int FileType );
150 //......................................................
151 void ProducePlots( std::vector< LoadedHistograms > AllTheHists, std::string CTier );
152 //......................................................
153 TGraph* graphAsymmError(TH1* hDataScaled, TH1* hData);
154 //......................................................
155 TH1D* MakeRatio( TH1D* num, TH1D* denom, int Col, std::string FType );
156 //......................................................
157 void SetRange( std::vector<TH1D*> Histos, std::vector<TH1D*> Ratios, // Pass the histograms
158  double XLow, double XHigh, // What X values do I want?
159  double H_YLow, double H_YHigh, double R_YLow, double R_YHigh ); // What Y values do I want?
160 //......................................................
161 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 );
162 //......................................................
163 void MakeSplitCans( std::string PlotName, std::string CutTier, std::string CanTitle, std::string PlotDir, TFile* OutFile, std::vector<TH1D*> Hists);
164 //......................................................
165 void PlotHistProp( TH1D* Hist, TLegend* L );
166 //......................................................
167 void CornerLabel(std::string Str);
168 //......................................................
169 void MakeTextFile( std::string Cut, std::string PTitle, std::string Name );
170 //......................................................
171 
172 //................................................................................
173 //................................................................................
174 void MakeThePlots( bool Horn, bool FullDataset = true, bool UsingCAFs = false ) {
175 
176  std::cout << "\n\nSetting IsFHC to " << Horn << "\n\n" << std::endl;
177  if (Horn) {IsFHC = true ; sFHC = "fhc"; LegNames[2] = "Wrong Sign: #bar{#nu_{#mu}} CC"; }
178  else {IsFHC = false; sFHC = "rhc"; LegNames[2] = "Wrong Sign: #nu_{#mu} CC"; }
179 
180  if (FullDataset) { AllPer = true; sPer = "full"; }
181  else { AllPer = false; sPer = "Epochs_7d_8b"; }
182 
183  // ---- First off, lets set the styles...
184  gStyle->SetOptStat(0);
185  gROOT->SetStyle("novaStyle");
186 
187  std::string sCAF_Concat = "res-concat";
188  if (UsingCAFs) sCAF_Concat = "caf";
189 
190  std::string BaseCAFFileDir="/nova/ana/users/karlwarb/Ana19/FDValidPlots/";
191 
192  // --- Declare a vector of my file name structures
193  std::string MontFileLoc = BaseCAFFileDir+"CompPlots_FD_PredNoExtrap_"+sFHC+"_full_"+sCAF_Concat+".root";
194  std::string NuMIFileLoc = BaseCAFFileDir+"CompPlots_FD_NuMI_" +sFHC+"_"+sPer+"_"+sCAF_Concat+".root";
195  std::string CosmFileLoc = BaseCAFFileDir+"CompPlots_FD_Cosmics_" +sFHC+"_"+sPer+"_"+sCAF_Concat+".root";
196 
197  // --- I need to open my files and figure out POT and Livetime for Data...
198  // -- MC
199  TFile *MontFile = TFile::Open( MontFileLoc.c_str() );
200  // -- NuMI
201  TFile *NuMIFile = TFile::Open( NuMIFileLoc.c_str() );
202  std::unique_ptr<Spectrum> TempNuMISpec = Spectrum::LoadFrom(NuMIFile, "sSliceTimeFull_AllCuts_AllQuants" ) ;
203  POTNom = TempNuMISpec->POT();
204  // -- Cosmics
205  TFile *CosmFile = TFile::Open( CosmFileLoc.c_str() );
206  std::unique_ptr<Spectrum> TempCosmSpec = Spectrum::LoadFrom(CosmFile, "sSliceTimeFull_AllCuts_AllQuants" ) ;
207  LivNom = TempCosmSpec->Livetime();
208  std::cout << "\n\nMy Livetime is now " << POTNom << " Exposure is now " << LivNom << "\n\n" << std::endl;
209 
210  // --- Now to load all of my histograms!
211  for (int cut=0; cut<NumCut; ++cut) {
212  // Check if just want to run a single cut?
213  //if (SingleCut && cut != 0) continue;
214  if (SingleCut && cut != 3) continue;
215 
216  for (int quant=0; quant<NumQuant; ++quant) {
217  if (SingleQuant && quant!=0) continue;
218 
219  // --- What is the Cut + Quant combination?
220  std::string CutQuant = CutNames[cut] + QuantNames[quant];
221  std::cout << "\n\nNow to look at " << CutQuant << std::endl;
222 
223  // --- Load the Monte Carlo histograms.
224  LoadedHistograms PredHists = LoadFile_GetHists( MontFile, CutQuant, 0 );
225  // --- Load the Monte Carlo histograms.
226  LoadedHistograms WrSiHists = LoadFile_GetHists( MontFile, CutQuant, 1 );
227  // --- Load the Monte Carlo histograms.
228  LoadedHistograms BeamHists = LoadFile_GetHists( MontFile, CutQuant, 2 );
229  // --- Load the NuMI histograms.
230  LoadedHistograms NuMIHists = LoadFile_GetHists( NuMIFile, CutQuant, 3 );
231  // --- Load the Cosmic histograms.
232  LoadedHistograms CosmHists = LoadFile_GetHists( CosmFile, CutQuant, 4 );
233 
234  // --- Make my plots
235  std::cout << "Now to make a bunch of plots..." << std::endl;
236  std::vector<LoadedHistograms> AllMyHists;
237  AllMyHists.push_back( NuMIHists );
238  AllMyHists.push_back( PredHists );
239  AllMyHists.push_back( WrSiHists );
240  AllMyHists.push_back( BeamHists );
241  AllMyHists.push_back( CosmHists );
242  // --- Add my Horn polarity information to CutQuant
243  CutQuant = "_"+sFHC+CutQuant;
244  // --- And now call ProducePlots!
245  ProducePlots( AllMyHists, CutQuant );
246 
247  } // Loop through cuts.
248  } // Loop through quantiles
249  // --- Finally, return.
250  return;
251 } // DoThePlots
252 
253 //......................................................
254 LoadedHistograms LoadFile_GetHists( TFile* InFile, std::string CutTier, int FileType ) {
255  // --- Declare the vector of LoadedHistograms which I am going to return...
256  LoadedHistograms MyHists;
257 
258  std::string YAxisStr;
259  double AxisPOT = 0;
260  if ( IsFHC ) YAxisStr = "Events / 9.48#times10^{20} POT-equiv";
261  if (!IsFHC ){ YAxisStr = "Events / 12.33#times10^{20} POT-equiv";
262  if (sPer == "Epochs_7d_8b" ) YAxisStr = "Events / 3.55#times10^{20} POT-equiv";
263  }
264 
265 
266  // --- Grab my spectra...
267  // *** Some key ones first.
268  MyHists.hReconstEnergy = GetSpectToHist( InFile, "sReconstEnergy"+CutTier, ";Reconstructed neutrino energy (GeV);" +YAxisStr+" / 0.1 GeV", FileType );
269  MyHists.hMuonEnergy = GetSpectToHist( InFile, "sMuonEnergy" +CutTier, ";Reconstructed muon energy (GeV);" +YAxisStr, FileType );
270  MyHists.hMuonEnPerHit = GetSpectToHist( InFile, "sMuonEnPerHit" +CutTier, ";Energy per hit of the muon (GeV);" +YAxisStr, FileType );
271  MyHists.hHadronEnergy = GetSpectToHist( InFile, "sHadronEnergy" +CutTier, ";Hadronic energy in the slice (GeV);" +YAxisStr, FileType );
272  MyHists.hHadroEnPerHit = GetSpectToHist( InFile, "sHadroEnPerHit"+CutTier, ";Hadronic energy per hit in the slice (GeV);"+YAxisStr, FileType );
273  MyHists.hHadFracEnergy = GetSpectToHist( InFile, "sHadFracEnergy"+CutTier, ";Hadronic energy fraction in the slice;" +YAxisStr, FileType );
274  MyHists.hHitsPerSlice = GetSpectToHist( InFile, "sHitsPerSlice" +CutTier, ";Number of hits in the slice;" +YAxisStr, FileType );
275  //*
276  MyHists.hSliceTimeFull = GetSpectToHist( InFile, "sSliceTimeFull"+CutTier, ";Mean time of the slice over full window (ns);"+YAxisStr, FileType );
277  MyHists.hSliceTimeZoom = GetSpectToHist( InFile, "sSliceTimeZoom"+CutTier, ";Mean time of the slice over beam window (ns);"+YAxisStr, FileType );
278  MyHists.hRemIDScore = GetSpectToHist( InFile, "sRemIDScore" +CutTier, ";RemID score;" +YAxisStr, FileType );
279  MyHists.hCVNCosmicScor = GetSpectToHist( InFile, "sCVNCosmicScor"+CutTier, ";CVN Cosmic score;" +YAxisStr, FileType );
280  MyHists.hCVNNuMuIDScor = GetSpectToHist( InFile, "sCVNNuMuIDScor"+CutTier, ";CVN #nu_{#mu} score;" +YAxisStr, FileType );
281  MyHists.hCVNNuMuID2017 = GetSpectToHist( InFile, "sCVNNuMuID2017"+CutTier, ";2017 CVN #nu_{#mu} score;" +YAxisStr, FileType );
282  MyHists.hNuMuContPID = GetSpectToHist( InFile, "sNuMuContPID" +CutTier, ";Third analysis BDT score;" +YAxisStr, FileType );
283  MyHists.hSANuMuContPID = GetSpectToHist( InFile, "sSANuMuContPID"+CutTier, ";Second analysis BDT score;" +YAxisStr, FileType );
284  // *** Some important Truth based ones.
285  MyHists.hTrueNuEnergy = GetSpectToHist( InFile, "sTrueNuEnergy" +CutTier, ";True neutrino energy in the slice (GeV);"+YAxisStr, FileType );
286  MyHists.hReTrOverTrNuE = GetSpectToHist( InFile, "sReTrOverTrNuE"+CutTier, ";(Reco-True)/True neutrino energy in the slice (GeV);"+YAxisStr, FileType);
287 
288  // *** Basic positional plots.
289  MyHists.hTrkStartX = GetSpectToHist( InFile, "sTrkStartX"+CutTier, ";Initial X position of the muon track (m);"+YAxisStr, FileType );
290  MyHists.hTrkStartY = GetSpectToHist( InFile, "sTrkStartY"+CutTier, ";Initial Y position of the muon track (m);"+YAxisStr, FileType );
291  MyHists.hTrkStartZ = GetSpectToHist( InFile, "sTrkStartZ"+CutTier, ";Initial Z position of the muon track (m);"+YAxisStr, FileType );
292  MyHists.hTrkEndX = GetSpectToHist( InFile, "sTrkEndX" +CutTier, ";Initial X position of the muon track (m);"+YAxisStr, FileType );
293  MyHists.hTrkEndY = GetSpectToHist( InFile, "sTrkEndY" +CutTier, ";Initial Y position of the muon track (m);"+YAxisStr, FileType );
294  MyHists.hTrkEndZ = GetSpectToHist( InFile, "sTrkEndZ" +CutTier, ";Initial Z position of the muon track (m);"+YAxisStr, FileType );
295  MyHists.hTrkLenXY = GetSpectToHist( InFile, "sTrkLenXY" +CutTier, ";Muon track length in the XY plane (m);"+YAxisStr, FileType );
296 
297  // *** Track based ones...
298  MyHists.hNumKalTracks = GetSpectToHist( InFile, "sNumKalTracks" +CutTier, ";Number of kalman tracks;"+YAxisStr, FileType );
299  MyHists.hNumKalTrHits = GetSpectToHist( InFile, "sNumKalTrHits" +CutTier, ";Number of hits in kalman trakcs;"+YAxisStr, FileType );
300  MyHists.hRatKalHitSlc = GetSpectToHist( InFile, "sRatKalHitSlc" +CutTier, ";Ratio of kalman track hits to total hits;"+YAxisStr, FileType );
301  MyHists.hKalTrLength = GetSpectToHist( InFile, "sKalTrLength" +CutTier, ";Kalman track length (m);"+YAxisStr, FileType );
302  MyHists.hKalTrBeamAng = GetSpectToHist( InFile, "sKalTrBeamAng" +CutTier, ";#cos(#theta) between kalman track and the beam;"+YAxisStr, FileType);
303  MyHists.hKalMostFwdCel = GetSpectToHist( InFile, "sKalMostFwdCel"+CutTier, ";Most forward cell in kalman track;"+YAxisStr, FileType );
304  MyHists.hKalMostBakCel = GetSpectToHist( InFile, "sKalMostBakCel"+CutTier, ";Most backward cell in kalman track;"+YAxisStr, FileType );
305  MyHists.hKalTrVer_MaxY = GetSpectToHist( InFile, "sKalTrVer_MaxY"+CutTier, ";Maximal Y vertex position in a slice (m);"+YAxisStr, FileType );
306  MyHists.hKalTrVer_MaxZ = GetSpectToHist( InFile, "sKalTrVer_MaxZ"+CutTier, ";Maximal Z vertex position in a slice (m);"+YAxisStr, FileType );
307  MyHists.hKalTrDir_Y = GetSpectToHist( InFile, "sKalTrDir_Y" +CutTier, ";Kalman track direction in Y; Direction in Y;"+YAxisStr, FileType );
308  MyHists.hScattKalTrLen = GetSpectToHist( InFile, "sScattKalTrLen"+CutTier, ";Sum of angular scattering over track length;"+YAxisStr, FileType );
309 
310  // *** Hit based ones...
311  MyHists.hFirstHitCell = GetSpectToHist( InFile, "sFirstHitCell" +CutTier, ";First cell hit in the slice;"+YAxisStr, FileType );
312  MyHists.hLastHitCell = GetSpectToHist( InFile, "sLastHitCell" +CutTier, ";Last cell hit in the slice;"+YAxisStr, FileType );
313  MyHists.hMaxActivity_Y = GetSpectToHist( InFile, "sMaxActivity_Y"+CutTier, ";Maximum position of activity in Y (m);"+YAxisStr, FileType );
314  MyHists.hMinActivity_Y = GetSpectToHist( InFile, "sMinActivity_Y"+CutTier, ";Minimum position of activity in Y (m);"+YAxisStr, FileType );
315  MyHists.hMaxActivity_Z = GetSpectToHist( InFile, "sMaxActivity_Z"+CutTier, ";Maximum position of activity in Z (m);"+YAxisStr, FileType );
316  MyHists.hMinActivity_Z = GetSpectToHist( InFile, "sMinActivity_Z"+CutTier, ";Minimum position of activity in Z (m);"+YAxisStr, FileType );
317  MyHists.hMinCellToEdge = GetSpectToHist( InFile, "sMinCellToEdge"+CutTier, ";Minimum number of cells hit to detector edge;"+YAxisStr, FileType );
318 
319  // *** RemID score / inputs
320  MyHists.hRemIDScatLLH = GetSpectToHist( InFile, "sRemIDScatLLH" +CutTier, ";Log-likelihood of scattering;"+YAxisStr, FileType );
321  MyHists.hRemIDdEdxLLH = GetSpectToHist( InFile, "sRemIDdEdxLLH" +CutTier, ";Log-likelihood of dEdx;"+YAxisStr, FileType );
322  MyHists.hRemIDMeasFrac = GetSpectToHist( InFile, "sRemIDMeasFrac"+CutTier, ";Log-likelihood of RemID;"+YAxisStr, FileType );
323 
324  // *** Some Calibration plots
325  for (int cor=0; cor<4; ++cor) {
326  MyHists.hCor_RecEn[cor] = GetSpectToHist( InFile, "sCor_RecEn"+CutTier+"_"+CornName[cor], ";Reconstructed En. for corner " +CornName[cor]+" (GeV);"+YAxisStr, FileType );
327  MyHists.hCor_MuoEn[cor] = GetSpectToHist( InFile, "sCor_MuoEn"+CutTier+"_"+CornName[cor], ";Muon En. for corner " +CornName[cor]+" (GeV);"+YAxisStr, FileType );
328  MyHists.hCor_HadEn[cor] = GetSpectToHist( InFile, "sCor_HadEn"+CutTier+"_"+CornName[cor], ";Hadronic En. for corner " +CornName[cor]+" (GeV);"+YAxisStr, FileType );
329  MyHists.hCor_HFrEn[cor] = GetSpectToHist( InFile, "sCor_HFrEn"+CutTier+"_"+CornName[cor], ";Hadronic En. fraction for corner "+CornName[cor]+" (GeV);"+YAxisStr, FileType );
330  }
331  //*/
332 
333  // --- Finally, return MyHists.
334  return MyHists;
335 }
336 
337 //......................................................
338 TH1D* GetSpectToHist( TFile* MyF, std::string LoadName, std::string Axis, int FileType ) {
339  // --- Declare my histogram.
340  TH1D* MyHist;
341 
342  // --- Am I looking at the NoExtrap MC?
343  if ( FileType == 0 || FileType == 1 || FileType == 2 ) {
344  std::unique_ptr<PredictionNoExtrap> TempSpec = PredictionNoExtrap::LoadFrom( MyF, LoadName.c_str() ) ;
345  // --- Need to set my oscillation calculator...
348  calc.SetdCP ( M_PI*2 ); // preliminary Ana2017
349  calc.SetDmsq32( 0.00248 );
350  calc.SetTh23 ( asin(sqrt(0.565)) );
351  /*
352  if (IsFHC) {
353  calc.SetDmsq32(2.45e-3);
354  calc.SetTh23(asin(sqrt(0.51)));
355  } else {
356  calc.SetDmsq32(2.57e-3);
357  calc.SetTh23(asin(sqrt(0.36)));
358  }
359  */
360 
361  // --- Am I looking at The nominal prediction?
362  if ( FileType == 0 ) {
363  MyHist = TempSpec -> Predict(&calc).ToTH1( POTNom );
364  //std::cout << "\t Looking at Best Fit Pred, integral is " << MyHist->Integral() << std::endl;
365  }
366  // --- Am I looking at the wrong sign component?
367  if ( FileType == 1 ) {
369  if (IsFHC) wrongs = Sign::kAntiNu;
370  MyHist = TempSpec->PredictComponent( &calc, Flavors::kAllNuMu, Current::kCC, wrongs ).ToTH1( POTNom );
371  //std::cout << "\t Looking at Wrong Sign component, integral is " << MyHist->Integral() << std::endl;
372  }
373  // --- Am I looking at the Beam background?
374  if ( FileType == 2 ) {
375  // Want everything other than muon neutrinos
376  TH1D* AllNue = TempSpec->PredictComponent( &calc, Flavors::kAllNuE , Current::kCC, Sign::kBoth ).ToTH1( POTNom );
377  TH1D* AllTau = TempSpec->PredictComponent( &calc, Flavors::kAllNuTau , Current::kCC, Sign::kBoth ).ToTH1( POTNom );
378  TH1D* AllNC = TempSpec->PredictComponent( &calc, Flavors::kAll , Current::kNC, Sign::kBoth ).ToTH1( POTNom );
379  TH1D* NueNumu = TempSpec->PredictComponent( &calc, Flavors::kNuEToNuMu, Current::kCC, Sign::kBoth ).ToTH1( POTNom );
380  // Now want to add all of those histograms together...
381  MyHist = (TH1D*)AllNue->Clone();
382  MyHist -> Add( AllTau );
383  MyHist -> Add( AllNC );
384  MyHist -> Add( NueNumu );
385  //std::cout << "\t Looking at Beam background, integral is " << MyHist->Integral() << std::endl;
386  }
387  } else {
388  // --- Not looking at MC, so looking at data.
389  std::unique_ptr<Spectrum> TempSpec = Spectrum::LoadFrom( MyF, LoadName.c_str() ) ;
390  // --- NuMI data.
391  if ( FileType == 3 ) {
392  MyHist = TempSpec -> ToTH1( POTNom );
393  //std::cout << "\t Looking at NuMI data, integral is " << MyHist->Integral() << std::endl;
394  }
395  // --- Cosmics.
396  if ( FileType == 4 ) {
397  MyHist = TempSpec -> ToTH1( LivNom, kLivetime );
398  //std::cout << "\t Looking at Cosmic data, integral is " << MyHist->Integral() << std::endl;
399  }
400  }
401 
402  // --- Make my plots pretty, and return.
403  CenterTitles ( MyHist);
404  MyHist->SetMinimum( 0 );
405  MyHist->SetMaximum( MyHist->GetMaximum()*1.05 );
406  MyHist->SetTitle ( Axis.c_str() );
407  return MyHist;
408 
409 }
410 
411 //......................................................
412 TGraph* graphAsymmError(TH1* hDataScaled, TH1* hData) {
413  for(int bin = 1; bin <= hData->GetNbinsX(); bin++) {
414  if(hData->GetBinContent(bin)==0){
415  hData->SetBinContent(bin,0.001);
416  hDataScaled->SetBinContent(bin,0.001);
417  }
418  }
419  TFeldmanCousins fc(0.6827); // to get correct poisson error bars
420  TGraphAsymmErrors *datapoisson = new TGraphAsymmErrors(hDataScaled);
421  for(int bin = 1; bin <= hData->GetNbinsX(); bin++) {
422  float binWidth = hData->GetBinWidth(bin);
423  float scaleFactor = 0.1/binWidth;
424  float y = hData->GetBinContent(bin);
425  float posup = (fc.CalculateUpperLimit(y,0));
426  float posdn = (fc.CalculateLowerLimit(y,0));
427  float errup0 = (fc.CalculateUpperLimit(y,0)-y);
428  float errdn0 = (y-fc.CalculateLowerLimit(y,0));
429  float errup = scaleFactor * (fc.CalculateUpperLimit(y,0)-y);
430  float errdn = scaleFactor * (y-fc.CalculateLowerLimit(y,0));
431  datapoisson->SetPointEYhigh(bin-1,errup);
432  datapoisson->SetPointEYlow(bin-1,errdn);
433  }
434  datapoisson->SetMarkerStyle(kFullCircle);
435  datapoisson->SetMarkerColor(1);
436  datapoisson->SetMarkerSize(1);
437  datapoisson->SetLineWidth(2);
438  return datapoisson;
439 }
440 
441 //......................................................
442 TH1D* MakeRatio( TH1D* num, TH1D* denom, int Col, std::string PlotName ) {
443  // Declare a new histogram.
444  TH1D* ratHist = new TH1D( TString(PlotName)+TString("_ratio_")+TString(std::to_string(Col)),
445  TString(";")+TString(num->GetXaxis()->GetTitle())+TString(";Ratio"),
446  num->GetNbinsX(),
447  num->GetXaxis()->GetXmin(),
448  num->GetXaxis()->GetXmax()
449  );
450  // Don't want any exponents?
451  ratHist->GetYaxis()->SetNoExponent();
452  ratHist->GetXaxis()->SetNoExponent();
453  // Divide the two input histograms...
454  ratHist->Divide( num , denom );
455  // Set Line and Marker colors / styles.
456  ratHist->SetLineColor ( Col );
457  ratHist->SetMarkerColor( Col );
458  ratHist->SetMarkerStyle( 20 );
459  ratHist->SetMarkerSize ( 0 );
460  // Figure out what the range should be...
461  ratHist->GetYaxis()->SetRangeUser( 0, 2 );
462  // Return my shiny new ratio plot!
463  CenterTitles(ratHist);
464  return ratHist;
465 }
466 
467 //......................................................
468 void SetRange( std::vector<TH1D*> Histos, std::vector<TH1D*> Ratios, // Pass the histograms
469  double XLow, double XHigh, // What X values do I want?
470  double H_YLow, double H_YHigh, double R_YLow, double R_YHigh ) { // What Y values do I want?
471  // --- First off, set the X ranges.
472  if ( (XLow != -1) && (XHigh != -1) ) {
473  for (unsigned int hist=0; hist<NPlot; ++hist) {
474  Histos[hist] -> GetXaxis() -> SetRangeUser( XLow, XHigh );
475  Ratios[hist] -> GetXaxis() -> SetRangeUser( XLow, XHigh );
476  }
477  }
478 
479  // --- And, the Y axes, but check that they aren't -1!
480  // Hists
481  if ( (H_YLow != -1) && (H_YHigh != -1) ) {
482  for (unsigned int hist=0; hist<NPlot; ++hist) {
483  Histos[hist] -> GetYaxis() -> SetRangeUser( H_YLow, H_YHigh );
484  }
485  }
486  // Ratios.
487  if ( (R_YLow != -1) && (R_YHigh != -1) ) {
488  for (unsigned int hist=0; hist<NPlot; ++hist) {
489  Ratios[hist] -> GetYaxis() -> SetRangeUser( R_YLow, R_YHigh );
490  }
491  }
492  // --- Return
493  return;
494 }
495 
496 //......................................................
497 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 ) {
498 
499  XLow = -1; XHigh = -1;
500  H_YLow = 0; H_YHigh = MaxVal;
501  R_YLow = 0.5; R_YHigh = 1.5;
502  SetLogy = false;
503  bool IsQuant = ( (Name.find("AllQuants") != std::string::npos) == true ? false:true);
504 
505  if ( (Name.find("ReconstEnergy") != std::string::npos) ) {
506  if (IsFHC) { if (!IsQuant) { H_YHigh = 12; } else { H_YHigh = 5; } }
507  else { if (!IsQuant) { H_YHigh = 10; } else { H_YHigh = 6; } }
508  }
509 
510  if ( (Name.find("MuonEnergy") != std::string::npos) ) {
511  if (IsFHC) { if (!IsQuant) { H_YHigh = 17; } else { H_YHigh = 10; } }
512  else { if (!IsQuant) { H_YHigh = 15; } else { H_YHigh = 8 ; } }
513  }
514 
515  if ( (Name.find("MuonEnPerHit") != std::string::npos) ) {
516  if (IsFHC) { if (!IsQuant) { H_YHigh = 50; } else { H_YHigh = 20; } }
517  else { if (!IsQuant) { H_YHigh = 50; } else { H_YHigh = 20; } }
518  XLow = 0.01; XHigh = 0.03;
519  }
520 
521  if ( (Name.find("HadronEnergy") != std::string::npos) ) {
522  if (IsFHC) { if (!IsQuant) { H_YHigh = 22; } else { H_YHigh = 18; } }
523  else { if (!IsQuant) { H_YHigh = 22; } else { H_YHigh = 15; } }
524  }
525 
526  if ( (Name.find("HadroEnPerHit") != std::string::npos) ) {
527  if (IsFHC) { if (!IsQuant) { H_YHigh = 20; } else { H_YHigh = 15; } }
528  else { if (!IsQuant) { H_YHigh = 20; } else { H_YHigh = 12; } }
529  XLow = 0.00; XHigh = 0.025;
530  }
531 
532  if ( (Name.find("HadFracEnergy") != std::string::npos) ) {
533  if (IsFHC) { if (!IsQuant) { H_YHigh = 18; } else { H_YHigh = 18; } }
534  else { if (!IsQuant) { H_YHigh = 20; } else { H_YHigh = 20; } }
535  }
536 
537  if ( (Name.find("HitsPerSlice") != std::string::npos) ) {
538  XLow = 0.; XHigh = 300;
539  if (IsFHC) { if (!IsQuant) { H_YHigh = 30; } else { H_YHigh = 15; } }
540  else { if (!IsQuant) { H_YHigh = 25; } else { H_YHigh = 10; } }
541  }
542 
543  if ( (Name.find("NumKalTracks") != std::string::npos) ) { XLow = 0.; XHigh = 5; }
544 
545  if ( (Name.find("NumKalTrHits") != std::string::npos) ) {
546  if (IsFHC) { if (!IsQuant) { H_YHigh = 30; } else { H_YHigh = 15; } }
547  else { if (!IsQuant) { H_YHigh = 25; } else { H_YHigh = 15; } }
548  XLow = 0.; XHigh = 300;
549  }
550 
551  if ( (Name.find("RemIDScore" ) != std::string::npos) ) {
552  if (IsFHC) { if (!IsQuant) { H_YHigh = 30; } else { H_YHigh = 15; } }
553  else { if (!IsQuant) { H_YHigh = 80; } else { H_YHigh = 30; } }
554  XLow = 0.8; XHigh = 1;
555  }
556 
557  if ( (Name.find("CVNNuMuIDScor" ) != std::string::npos) ) {
558  if (IsFHC) { if (!IsQuant) { H_YHigh = 30; } else { H_YHigh = 15; } }
559  else { if (!IsQuant) { H_YHigh = 40; } else { H_YHigh = 15; } }
560  XLow = 0.5; XHigh = 1;
561  }
562 
563  if ( (Name.find("CVNNuMuID2017" ) != std::string::npos) ) {
564  if (IsFHC) { if (!IsQuant) { H_YHigh = 30; } else { H_YHigh = 15; } }
565  else { if (!IsQuant) { H_YHigh = 30; } else { H_YHigh = 10; } }
566  XLow = 0.5; XHigh = 1;
567  }
568 
569  if ( (Name.find("CVNCosmicScor") != std::string::npos) ) {
570  if (IsFHC) { if (!IsQuant) { H_YHigh = 30 ; } else { H_YHigh = 15; } }
571  else { if (!IsQuant) { H_YHigh = 125; } else { H_YHigh = 40; } }
572  XLow = 0 ; XHigh = 0.1;
573  }
574 
575  if ( (Name.find("NuMuContPID" ) != std::string::npos) ) {
576  if (IsFHC) { if (!IsQuant) { H_YHigh = 30; } else { H_YHigh = 15; } }
577  else { if (!IsQuant) { H_YHigh = 30; } else { H_YHigh = 10; } }
578  XLow = 0.4; XHigh = 0.9;
579  }
580 
581  if ( (Name.find("Activity_Y" ) != std::string::npos) ) { H_YLow = 0.1; }
582  if ( (Name.find("Activity_Z" ) != std::string::npos) ) { H_YLow = 0.1; }
583  if ( (Name.find("KalMost" ) != std::string::npos) ) { H_YLow = 0.1; }
584 
585  if ( (Name.find("TrkStart") != std::string::npos) || (Name.find("TrkEnd" ) != std::string::npos) ) {
586  if (IsFHC) { H_YLow = 0; H_YHigh = 15; }
587  else { if (!IsQuant) { H_YHigh = 20; } else { H_YHigh = 10; } }
588  }
589 
590  if ( (Name.find("TrkLenXY") != std::string::npos) ) {
591  if (IsFHC) { H_YLow = 0; H_YHigh = 15; }
592  else { if (!IsQuant) { H_YHigh = 30; } else { H_YHigh = 10; } }
593  XLow = 0; XHigh = 10;
594  }
595 
596  if ( (Name.find("KalTrLength") != std::string::npos) ) {
597  if (IsFHC) { if (!IsQuant) { H_YHigh = 15; } else { H_YHigh = 8; } }
598  else { if (!IsQuant) { H_YHigh = 15; } else { H_YHigh = 8; } }
599  }
600 
601  if ( (Name.find("KalTrAng") != std::string::npos) ) {
602  if (IsFHC) { if (!IsQuant) { H_YHigh = 15; } else { H_YHigh = 8; } }
603  else { if (!IsQuant) { H_YHigh = 60; } else { H_YHigh = 30; } }
604  XLow = 0.5; XHigh = 1;
605  }
606 
607  if ( (Name.find("KalMost") != std::string::npos) ) {
608  if (IsFHC) { if (!IsQuant) { H_YHigh = 20; } else { H_YHigh = 8; } }
609  else { if (!IsQuant) { H_YHigh = 20; } else { H_YHigh = 8; } }
610  }
611 
612  if ( (Name.find("KalTrVer") != std::string::npos) ) {
613  if (IsFHC) { if (!IsQuant) { H_YHigh = 20; } else { H_YHigh = 8; } }
614  else { if (!IsQuant) { H_YHigh = 15; } else { H_YHigh = 5; } }
615  }
616 
617  if ( (Name.find("KalTrDir_Y") != std::string::npos) ) {
618  if (IsFHC) { if (!IsQuant) { H_YHigh = 20; } else { H_YHigh = 8; } }
619  else { if (!IsQuant) { H_YHigh = 25; } else { H_YHigh = 10; } }
620  }
621 
622  if ( (Name.find("ScattKalTrLen") != std::string::npos) ) {
623  if (IsFHC) { if (!IsQuant) { H_YHigh = 20; } else { H_YHigh = 8; } }
624  else { if (!IsQuant) { H_YHigh = 40; } else { H_YHigh = 15; } }
625  }
626 
627  if ( (Name.find("HitCell") != std::string::npos) ) {
628  if (IsFHC) { if (!IsQuant) { H_YHigh = 20; } else { H_YHigh = 8; } }
629  else { if (!IsQuant) { H_YHigh = 60; } else { H_YHigh = 20; } }
630  XLow = 0; XHigh = 400;
631  }
632 
633  if ( (Name.find("Activity_Y") != std::string::npos) ) {
634  if (IsFHC) { if (!IsQuant) { H_YHigh = 20; } else { H_YHigh = 8; } }
635  else { if (!IsQuant) { H_YHigh = 20; } else { H_YHigh = 5; } }
636  }
637 
638  if ( (Name.find("Activity_Z") != std::string::npos) ) {
639  if (IsFHC) { if (!IsQuant) { H_YHigh = 20; } else { H_YHigh = 8; } }
640  else { if (!IsQuant) { H_YHigh = 15; } else { H_YHigh = 8; } }
641  }
642 
643  if ( (Name.find("CellToEdge") != std::string::npos) ) {
644  if (IsFHC) { if (!IsQuant) { H_YHigh = 20; } else { H_YHigh = 8; } }
645  else { if (!IsQuant) { H_YHigh = 15; } else { H_YHigh = 5; } }
646  }
647 
648  if ( (Name.find("RemIDScat") != std::string::npos) ) {
649  if (IsFHC) { if (!IsQuant) { H_YHigh = 20; } else { H_YHigh = 8; } }
650  else { if (!IsQuant) { H_YHigh = 40; } else { H_YHigh = 15; } }
651  XLow = -0.2; XHigh = 0.3;
652  }
653 
654  if ( (Name.find("RemIDdEdx") != std::string::npos) ) {
655  if (IsFHC) { if (!IsQuant) { H_YHigh = 20; } else { H_YHigh = 8; } }
656  else { if (!IsQuant) { H_YHigh = 30; } else { H_YHigh = 15; } }
657  XLow = -0.2; XHigh = 0.5;
658  }
659 
660  if ( (Name.find("RemIDMeas") != std::string::npos) ) {
661  if (IsFHC) { if (!IsQuant) { H_YHigh = 20; } else { H_YHigh = 8; } }
662  else { if (!IsQuant) { H_YHigh = 30; } else { H_YHigh = 20; } }
663  XLow = 0.5; XHigh = 1.0;
664  }
665 
666 
667 
668  /*
669  std::cout << "\t\tParams are; Xlow " << XLow << ", XHigh: " << XHigh << ", H_YLow: " << H_YLow << ", H_YHigh: " << H_YHigh
670  << ", R_YLow: " << R_YLow << ", R_YHigh " << R_YHigh << ", SetLogy " << SetLogy
671  << std::endl;
672  */
673  return;
674 }
675 
676 //......................................................
677 void PlotHistProp( TH1D* Hist, TLegend* L ) {
678  double Mean_Nom = Hist->GetMean();
679  double RMS_Nom = Hist->GetRMS();
680  double Int_Nom = Hist->Integral();
681  L -> AddEntry( Hist, Form("Mean: %.3f", Mean_Nom), "" );
682  L -> AddEntry( Hist, Form("RMS: %.3f" , RMS_Nom ), "" );
683  L -> AddEntry( Hist, Form("Int: %.3f" , Int_Nom ), "" );
684  return;
685 }
686 //......................................................
688  TLatex* CornLab = new TLatex(.1, .93, Str.c_str());
689  CornLab->SetTextColor(kGray+1);
690  CornLab->SetNDC();
691  CornLab->SetTextSize (2/30.);
692  CornLab->SetTextAlign(11);
693  CornLab->Draw();
694 }
695 //..................................
697  // Does this canvas have a ratio on it?
698  std::string WRat = ", with data/mc ratio.";
699  if (Name.find("_hist") != std::string::npos) WRat = ".";
700  // Remove the underscores from my cut tier
701  while (Cut.find("_") != std::string::npos)
702  Cut.replace(Cut.find("_"),1,", "); //CTier = "Cosmic_Rej";
703  // Determine my caption
704  std::string Cap = "Plot showing the distribution of " + PTitle + WRat + " {VarName, Horn, Cut, Quantile} is; " + Cut;
705  std::cout << "\nFile: " << Name << "\n\t Has caption: " << Cap << ".\n" << std::endl;
706  // Make the output text file
707  std::ofstream TxtOut ( Name.c_str(), std::ofstream::out );
708  TxtOut << Cap;
709  TxtOut.close();
710  // return
711  return;
712 }
713 //......................................................
714 void MakeSplitCans( std::string PlotName, std::string CutTier, std::string CanTitle, std::string PlotDir, TFile* OutFile, std::vector<TH1D*> Hists) {
715  PlotName += CutTier;
716  /*
717  Hists[0] == NuMI Data
718  Hists[1] == Best Fit Pred
719  Hists[2] == Wrong Sign Bkg
720  Hists[3] == Total Beam Bkg
721  Hists[4] == Cosmic Bkg
722  */
723  // Make a copy of NuMI data before any changes...
724  TH1D* NuMIClone = (TH1D*)Hists[0]->Clone();
725 
726  // Add cosmic to nominal and wrong sign.
727  Hists[1] -> Add( Hists[4] ); // Add cosmics to Prediction
728  Hists[2] -> Add( Hists[4] ); // Add cosmics to Wrong Sign
729 
730  // --- What do I want to rebin?
731  double ReFac = 4;
732  if ( (PlotName.find("ReconstEnergy") != std::string::npos) ||
733  (PlotName.find("RemIDScore" ) != std::string::npos) ||
734  (PlotName.find("CVNCosmicScor") != std::string::npos) ||
735  (PlotName.find("NumKalTracks" ) != std::string::npos) ||
736  (PlotName.find("RemIDScattLLH") != std::string::npos) ||
737  (PlotName.find("RemIDdEdxLLH" ) != std::string::npos) ||
738  (PlotName.find("RemIDMeasFrac") != std::string::npos)
739  ) {
740  ReFac = 0;
741  }
742  if ( (PlotName.find("HitsPerSlice" ) != std::string::npos) || // 0 to 300
743  (PlotName.find("MuonEnPerHit" ) != std::string::npos) || // 0 to 0.04
744  (PlotName.find("HadroEnPerHit") != std::string::npos) || // 0 to 0.04
745  (PlotName.find("SliceTimeZoom") != std::string::npos) ||
746  (PlotName.find("CVNNuMuIDScor") != std::string::npos) ||
747  (PlotName.find("CVNNuMuID2017") != std::string::npos) ||
748  (PlotName.find("SANuMuContPID") != std::string::npos) ||
749  (PlotName.find("TrkLenXY" ) != std::string::npos) ||
750  (PlotName.find("NumKalTrHits" ) != std::string::npos)
751  ) {
752  ReFac = 2;
753  }
754  if ( (PlotName.find("TrkStart") != std::string::npos) || // 0 to 300
755  (PlotName.find("TrkEnd" ) != std::string::npos) // 0 to 300
756  ) {
757  ReFac = 5;
758  }
759  if (ReFac != 0) {
760  for (unsigned int hist=0; hist<NPlot; ++hist) {
761  Hists[hist] -> Rebin( ReFac );
762  }
763  }
764 
765  // --- Set the colour of the histograms...
766  for (unsigned int hist=0; hist<NPlot; ++hist) {
767  Hists[hist] -> SetLineColor( Colours[hist] );
768  if (hist != 1) Hists[hist] -> SetFillColor( Colours[hist] );
769  }
770 
771  // --- If Reco energy then normalise.
772  if (PlotName.find("ReconstEnergy") != std::string::npos) {
773  for (unsigned int hist=0; hist<NPlot; ++hist) {
774  Hists[hist] -> Scale( 0.1, "width" );
775  }
776  }
777 
778  // --- Make ratios.
779  std::vector<TH1D*> Ratios;
780  for (unsigned int hist=0; hist<NPlot; ++hist) {
781  Ratios.push_back( MakeRatio( Hists[hist], Hists[0], Colours[hist], PlotName ) );
782  }
783 
784  // --- What axis ranges do I want to use?
785  double XL, XH, H_YL, H_YH, R_YL, R_YH;
786  bool SetLogy;
787  FindAxisRange( XL, XH, H_YL, H_YH, R_YL, R_YH, 1.5*Hists[1]->GetMaximum(), SetLogy, PlotName );
788  // --- And set the ranges...
789  SetRange( Hists, Ratios, XL, XH, H_YL, H_YH, R_YL, R_YH);
790 
791  // For DataMC I want to make a graph of the 0th shift. If not MCData then I'll just ignore this...
792  TGraph* grData;
793  if (PlotName.find("ReconstEnergy") != std::string::npos) {
794  grData = graphAsymmError(Hists[0], NuMIClone);
795  } else {
796  grData = (TGraph*)GraphWithPoissonErrors(Hists[0], false, false);
797  }
798  grData -> SetMarkerStyle(kFullCircle);
799  grData -> SetLineWidth(2);
800 
801  // --- Now, let's make a histogram with just the plot.
802  std::string HCanNa = PlotName + "_hist";
803  TCanvas* HCan = new TCanvas( HCanNa.c_str(), HCanNa.c_str() );
804  // Where do I want my legend...
805  double HistLeg_X1 = 0.55, HistLeg_Y1 = 0.05, HistLeg_X2 = 0.85, HistLeg_Y2 = 0.88;
806  // If not plotting mean vals, then want smaller legend.
807  if (!PlotMean) {
808  HistLeg_Y1 = 0.55;
809  }
810  // Do I want the legend on the left instead?
811  if ( //(PlotName.find("ReconstEnergy") != std::string::npos) ||
812  ( PlotName.find("RatKalHitSlc" ) != std::string::npos) ||
813  ( PlotName.find("RemIDScore" ) != std::string::npos) ||
814  ( PlotName.find("CVNNuMuID" ) != std::string::npos) ||
815  ( PlotName.find("NuMuContPID" ) != std::string::npos) ||
816  ((PlotName.find("MuonEnergy" ) != std::string::npos) && (PlotName.find("Quant4") == std::string::npos)) ||
817  ( PlotName.find("TrkStart" ) != std::string::npos) || // 0 to 300
818  ( PlotName.find("TrkEnd" ) != std::string::npos) // 0 to 300
819  ) {
820  HistLeg_X1 = 0.10;
821  HistLeg_X2 = 0.56;
822  }
823  // Make the legend where I've decided...
824  TLegend *HistLeg = new TLegend( HistLeg_X1, HistLeg_Y1, HistLeg_X2, HistLeg_Y2 );
825  HistLeg -> SetFillStyle (0);
826  HistLeg -> SetBorderSize(0);
827  HistLeg -> SetTextSize( 0.8*Hists[0]->GetXaxis()->GetTitleSize() );
828 
829  // --- And now start drawing things!!!
830  Hists[1] -> DrawCopy("hist");
831  for (unsigned int hist=0; hist<NPlot; ++hist) {
832  if (hist==0) {
833  HistLeg -> AddEntry( grData, LegNames[0].c_str(), "lep" );
834  //if (PlotMean) PlotHistProp( Hists[0], HistLeg );
835  } else if (hist==1) {
836  Hists[hist] -> Draw("hist same");
837  HistLeg -> AddEntry( Hists[hist], LegNames[hist].c_str(), "l" );
838  } else {
839  Hists[hist] -> Draw("hist same");
840  HistLeg -> AddEntry( Hists[hist], LegNames[hist].c_str(), "f" );
841  //if (PlotMean) PlotHistProp( Hists[hist], HistLeg );
842  }
843  }
844  // --- Add which quantile
845  if (PlotName.find("Quant1") != std::string::npos) { HistLeg -> AddEntry( grData, "Quartile 1", "" ); }
846  else if (PlotName.find("Quant2") != std::string::npos) { HistLeg -> AddEntry( grData, "Quartile 2", "" ); }
847  else if (PlotName.find("Quant3") != std::string::npos) { HistLeg -> AddEntry( grData, "Quartile 3", "" ); }
848  else if (PlotName.find("Quant4") != std::string::npos) { HistLeg -> AddEntry( grData, "Quartile 4", "" ); }
849  // --- Add whether FHC or RHC
850  if (IsFHC ) CornerLabel( "Neutrino beam" );
851  else CornerLabel( "Antineutrino beam" );
852  // --- Redraw the axes, and the data.
853  grData -> Draw("ep same");
854  // Do I want the canvas to be loagarithmic?
855  if ( SetLogy ) HCan -> SetLogy();
856  // And, draw the legend!
857  HistLeg -> Draw();
858  Preliminary();
859 
860  /*
861  // --- Get ready to make the canvas with ratios.
862  std::string RatCanNa = FType + CTier + DirNa;
863  TCanvas* RatCan = new TCanvas( RatCanNa.c_str(), CanTi.c_str() );
864  // Sort out the canvas.
865  RatCan -> SetBottomMargin( 0. );
866  double Spl = 0.4;
867  TPad* P1 = new TPad( "Temp_1", "", 0.0, Spl, 1.0, 1.0, 0 );
868  TPad* P2 = new TPad( "Temp_2", "", 0.0, 0.0, 1.0, Spl, 0 );
869  P2 -> SetRightMargin (.03);
870  P2 -> SetTopMargin (.00);
871  P2 -> SetBottomMargin(.22);
872  P2 -> SetLeftMargin (.08);
873  P2 -> Draw();
874  P1 -> SetRightMargin (.03);
875  P1 -> SetLeftMargin (.08);
876  P1 -> SetTopMargin (.10);
877  P1 -> SetBottomMargin(.03);
878  P1 -> Draw();
879  // Set some label sizes.
880  double Lb1 = 0.07;
881  double Lb2 = 0.10;
882  // --- First, draw the ratios so cd onto Pad2
883  P2 -> cd();
884  // Set axis ranges etc.
885  Rat[0] -> GetYaxis() -> SetTitleSize( Lb2 );
886  Rat[0] -> GetYaxis() -> SetTitleOffset(0.4);
887  Rat[0] -> GetYaxis() -> SetLabelSize( Lb2 );
888  Rat[0] -> GetXaxis() -> SetTitleSize( Lb2 );
889  Rat[0] -> GetXaxis() -> SetLabelSize( Lb2 );
890  Rat[0] -> Draw();
891  // Add a reference line
892  //TLine *line = new TLine( Nom->GetXaxis()->GetXmin(), 1, Nom->GetXaxis()->GetXmax(), 1 );
893  double x[2] = { Nom->GetXaxis()->GetXmin(), Nom->GetXaxis()->GetXmax() };
894  double y[2] = { 1, 1 };
895  TGraph *line = new TGraph( 2, x, y );
896  line->SetLineColor(kGray);
897  line->SetLineWidth( 3 );
898  line->SetLineStyle( 2 );
899  line->Draw("l same");
900  // Redraw Rat[0]
901  Rat[0] -> Draw("axis same");
902  if (MySys != "DataMC") {
903  for (unsigned int sh=1; sh<NPlot; ++sh) {
904  Rat[sh] -> Draw("same");
905  }
906  }
907  // --- And cd into Pad1
908  P1 -> cd();
909  // Where do I want my legend...
910  double Leg_X1 = 0.62, Leg_Y1 = 0.05, Leg_X2 = 0.94, Leg_Y2 = 0.85;
911  // If not plotting mean vals, then want smaller legend.
912  if (!PlotMean) Leg_Y1 = 0.4;
913  // Do I want the legend on the left instead?
914  if ( (DirNa.find("RatKalHitSlc") != std::string::npos) ||
915  (DirNa.find("RemIDScore") != std::string::npos) ||
916  (DirNa.find("_cCVNNuMuID") != std::string::npos) ||
917  (DirNa.find("_cNuMuContPID")!= std::string::npos)
918  ) {
919  Leg_X1 = 0.1;
920  Leg_X2 = 0.32;
921  }
922 
923  // --- Special Legend positions for DataMC.
924  // For Start/EndPos plots want to split the legend.
925 
926  // Make the legend where I've decided...
927  TLegend *Leg = new TLegend( Leg_X1, Leg_Y1, Leg_X2, Leg_Y2 );
928  Leg -> SetFillStyle (0);
929  Leg -> SetBorderSize(0);
930  Leg -> SetTextSize( 1.25*HistLeg -> GetTextSize() );
931  // --- Draw nominal first!
932  Nom -> GetYaxis() -> SetTitleSize( Lb1 );
933  Nom -> GetYaxis() -> SetLabelSize( Lb1 );
934  Nom -> GetYaxis() -> SetTitleOffset( 0.5 );
935  // Remove the x axis labels
936  Nom -> GetXaxis() -> SetLabelSize (0 );
937  Nom -> GetXaxis() -> SetLabelOffset(99);
938  Nom -> DrawCopy("hist");
939  // --- If DataMC want my Data as first thing on the legendto the legend
940  if (MySys == "DataMC" ) Leg -> AddEntry( grData, LegNames[1].c_str(), "lep" );
941  // --- Add nominal to legend
942  Leg -> AddEntry( Nom, LegNames[0].c_str(), "l" );
943  if (PlotMean) PlotHistProp( Nom, Leg );
944  // --- Now draw the shifts!
945  for (unsigned int sh=0; sh<NPlot; ++sh) {
946  if (MySys == "DataMC" ) {
947  if ( sh != 0 ) {
948  Shi[sh] -> DrawCopy("hist same");
949  Leg -> AddEntry( Shi[sh], LegNames[sh+1].c_str(), "f" );
950  }
951  } else {
952  Shi[sh] -> DrawCopy("hist same");
953  Leg -> AddEntry( Shi[sh], LegNames[sh+1].c_str(), "l" );
954  }
955  if (PlotMean) PlotHistProp( Shi[sh], Leg );
956  }
957  // Want to draw data last
958  if (MySys == "DataMC") {
959  grData -> Draw("ep same");
960  }
961  Nom -> Draw( "axis same" );
962  // Do I want Pad 1 to be loagarithmic?
963  if ( SetLogy ) P1 -> SetLogy();
964  // And, draw the legend!
965  Leg -> Draw();
966  Preliminary();
967  */
968 
969  // --- Now I can save my things...
970  std::vector<std::string> Exten = { ".png", ".pdf", ".txt" };
971  for (size_t ex=0; ex<Exten.size(); ++ex) {
972  // Determine the name of my file.
973  std::string HistoCanSaveLoc = PlotDir + HCan ->GetName() + Exten[ex];
974  // Do something special if I have a .txt file
975  if ( Exten[ex] == ".txt" ) {
976  MakeTextFile( PlotName, Hists[0]->GetXaxis()->GetTitle(), HistoCanSaveLoc );
977  } else {
978  // Save the canvases if not.
979  HCan -> SaveAs( HistoCanSaveLoc.c_str() );
980  }
981  // And write the canvases to file too.
982  OutFile -> cd();
983  HCan -> Write( HCan->GetName() );
984  }
985  // --- Return.
986  return;
987 }
988 
989 //.....................................................
990 void ProducePlots( std::vector< LoadedHistograms > AllTheHists, std::string CTier ) {
991 
992  // --- Make sure that my directory exists.
993  std::string TopDir = "Plots_"+sFHC+"_"+sPer+"/";
994  std::string PloDir = TopDir+CTier+"/";
995  gSystem -> MakeDirectory( TopDir.c_str() );
996  gSystem -> MakeDirectory( PloDir.c_str() );
997 
998  // --- Make my output file.
999  std::string OName = PloDir+"Merged"+CTier+".root";
1000  TFile* OutFile = TFile::Open( OName.c_str(), "RECREATE" );
1001  std::cout << "I want to make a file called " << OName << std::endl;
1002 
1003  // --- Now let's start making all of the plots!!!!
1004  // -----------------------------------------------
1005  std::vector<TH1D*> MyHists;
1006 
1007  // *** Some key ones first.
1008  // The reconstructed energy
1009  //MyHists.clear(); for (unsigned int sh=0; sh<NPlot; ++sh) { MyHists.push_back( AllTheHists[sh].hReconstEnergy ); }
1010  //MakeSplitCans( "ReconstEnergy", CTier, ("Total reconstructed energy for "+CTier), PloDir, OutFile, MyHists );
1011 
1012  // The reconstructed Hadronic energy
1013  MyHists.clear(); for (unsigned int sh=0; sh<NPlot; ++sh) { MyHists.push_back( AllTheHists[sh].hHadronEnergy ); }
1014  MakeSplitCans( "HadronEnergy" , CTier, ("Hadronic energy for "+CTier) , PloDir, OutFile, MyHists );
1015 
1016  // The reconstructed hadronic energy per hit
1017  MyHists.clear(); for (unsigned int sh=0; sh<NPlot; ++sh) { MyHists.push_back( AllTheHists[sh].hHadroEnPerHit ); }
1018  MakeSplitCans( "HadroEnPerHit", CTier, ("Hadronic energy per hit for "+CTier) , PloDir, OutFile, MyHists );
1019 
1020  // The Hadronic energy fraction
1021  MyHists.clear(); for (unsigned int sh=0; sh<NPlot; ++sh) { MyHists.push_back( AllTheHists[sh].hHadFracEnergy ); }
1022  MakeSplitCans( "HadFracEnergy", CTier, ("Hadronic energy fraction for "+CTier) , PloDir, OutFile, MyHists );
1023 
1024  // The reconstructed muon energy
1025  MyHists.clear(); for (unsigned int sh=0; sh<NPlot; ++sh) { MyHists.push_back( AllTheHists[sh].hMuonEnergy ); }
1026  MakeSplitCans( "MuonEnergy" , CTier, ("Muon energy for "+CTier) , PloDir, OutFile, MyHists );
1027  /*
1028  // The reconstructed muon energy per hit
1029  MyHists.clear(); for (unsigned int sh=0; sh<NPlot; ++sh) { MyHists.push_back( AllTheHists[sh].hMuonEnPerHit ); }
1030  MakeSplitCans( "MuonEnPerHit" , CTier, ("Muon energy per hit for "+CTier) , PloDir, OutFile, MyHists );
1031 
1032  // Number of hits per slice
1033  MyHists.clear(); for (unsigned int sh=0; sh<NPlot; ++sh) { MyHists.push_back( AllTheHists[sh].hHitsPerSlice ); }
1034  MakeSplitCans( "HitsPerSlice" , CTier, ("Hits per slice for "+CTier) , PloDir, OutFile, MyHists );
1035 
1036  //*
1037  // Mean time of slice Full
1038  MyHists.clear(); for (unsigned int sh=0; sh<NPlot; ++sh) { MyHists.push_back( AllTheHists[sh].hSliceTimeFull ); }
1039  MakeSplitCans( "SliceTimeFull", CTier , ("Mean time of slice full window "+CTier), PloDir, OutFile, MyHists );
1040  // Mean time of slice Zoom
1041  MyHists.clear(); for (unsigned int sh=0; sh<NPlot; ++sh) { MyHists.push_back( AllTheHists[sh].hSliceTimeZoom ); }
1042  MakeSplitCans( "SliceTimeZoom", CTier , ("Mean time of slice beam window "+CTier), PloDir, OutFile, MyHists );
1043  // The RemID score
1044  MyHists.clear(); for (unsigned int sh=0; sh<NPlot; ++sh) { MyHists.push_back( AllTheHists[sh].hRemIDScore ); }
1045  MakeSplitCans( "RemIDScore" , CTier , ("The RemID score "+CTier) , PloDir, OutFile, MyHists );
1046  // CVN Cosmic score
1047  MyHists.clear(); for (unsigned int sh=0; sh<NPlot; ++sh) { MyHists.push_back( AllTheHists[sh].hCVNCosmicScor ); }
1048  MakeSplitCans( "CVNCosmicScor" , CTier , ("CVN Cosmic score "+CTier) , PloDir, OutFile, MyHists );
1049  // CVN NuMu score
1050  MyHists.clear(); for (unsigned int sh=0; sh<NPlot; ++sh) { MyHists.push_back( AllTheHists[sh].hCVNNuMuIDScor ); }
1051  MakeSplitCans( "CVNNuMuIDScor" , CTier , ("CVN NuMu score "+CTier) , PloDir, OutFile, MyHists );
1052  // CVN NuMu 2017 score
1053  MyHists.clear(); for (unsigned int sh=0; sh<NPlot; ++sh) { MyHists.push_back( AllTheHists[sh].hCVNNuMuID2017 ); }
1054  MakeSplitCans( "CVNNuMuID2017" , CTier , ("2017 CVN NuMu score "+CTier) , PloDir, OutFile, MyHists );
1055  // Third analysis BDT
1056  MyHists.clear(); for (unsigned int sh=0; sh<NPlot; ++sh) { MyHists.push_back( AllTheHists[sh].hNuMuContPID ); }
1057  MakeSplitCans( "NuMuContPID" , CTier , ("Third analysis BDT "+CTier) , PloDir, OutFile, MyHists );
1058  // Second analysis BDT cut
1059  MyHists.clear(); for (unsigned int sh=0; sh<NPlot; ++sh) { MyHists.push_back( AllTheHists[sh].hSANuMuContPID ); }
1060  MakeSplitCans( "SANuMuContPID" , CTier , ("Second analysis BDT "+CTier) , PloDir, OutFile, MyHists );
1061 
1062  // *** Some important Truth based ones.
1063  // The true neutrino energy
1064  MyHists.clear(); for (unsigned int sh=0; sh<NPlot; ++sh) { MyHists.push_back( AllTheHists[sh].hTrueNuEnergy ); }
1065  MakeSplitCans( "TrueNuEnergy" , CTier , ("True neutrino energy for "+CTier) , PloDir, OutFile, MyHists );
1066  // The Reco - True / True neutrino energy
1067  MyHists.clear(); for (unsigned int sh=0; sh<NPlot; ++sh) { MyHists.push_back( AllTheHists[sh].hReTrOverTrNuE ); }
1068  MakeSplitCans( "ReTrOverTrNuE", CTier , ("Re-Tr/Te neutrino energy for "+CTier) , PloDir, OutFile, MyHists );
1069 
1070  // *** Basic positional plots
1071  // Start X
1072  MyHists.clear(); for (unsigned int sh=0; sh<NPlot; ++sh) { MyHists.push_back( AllTheHists[sh].hTrkStartX ); }
1073  MakeSplitCans( "TrkStartX" , CTier , ("The Start X Pos "+CTier), PloDir, OutFile, MyHists );
1074  // Start Y
1075  MyHists.clear(); for (unsigned int sh=0; sh<NPlot; ++sh) { MyHists.push_back( AllTheHists[sh].hTrkStartY ); }
1076  MakeSplitCans( "TrkStartY" , CTier , ("The Start Y Pos "+CTier), PloDir, OutFile, MyHists );
1077  // Start Z
1078  MyHists.clear(); for (unsigned int sh=0; sh<NPlot; ++sh) { MyHists.push_back( AllTheHists[sh].hTrkStartZ ); }
1079  MakeSplitCans( "TrkStartZ" , CTier , ("The Start Z Pos "+CTier), PloDir, OutFile, MyHists );
1080  // End X
1081  MyHists.clear(); for (unsigned int sh=0; sh<NPlot; ++sh) { MyHists.push_back( AllTheHists[sh].hTrkEndX ); }
1082  MakeSplitCans( "TrkEndX" , CTier , ("The End X Pos "+CTier) , PloDir, OutFile, MyHists );
1083  // End Y
1084  MyHists.clear(); for (unsigned int sh=0; sh<NPlot; ++sh) { MyHists.push_back( AllTheHists[sh].hTrkEndY ); }
1085  MakeSplitCans( "TrkEndY" , CTier , ("The End Y Pos "+CTier) , PloDir, OutFile, MyHists );
1086  // End Z
1087  MyHists.clear(); for (unsigned int sh=0; sh<NPlot; ++sh) { MyHists.push_back( AllTheHists[sh].hTrkEndZ ); }
1088  MakeSplitCans( "TrkEndZ" , CTier , ("The End Z Pos "+CTier) , PloDir, OutFile, MyHists );
1089  // Len XY
1090  MyHists.clear(); for (unsigned int sh=0; sh<NPlot; ++sh) { MyHists.push_back( AllTheHists[sh].hTrkLenXY); }
1091  MakeSplitCans( "TrkLenXY" , CTier , ("Length in XY " +CTier) , PloDir, OutFile, MyHists );
1092 
1093  // *** Track based ones...
1094  // Number of Kal trk in slice
1095  MyHists.clear(); for (unsigned int sh=0; sh<NPlot; ++sh) { MyHists.push_back( AllTheHists[sh].hNumKalTracks ); }
1096  MakeSplitCans( "NumKalTracks" , CTier , ("Number of kalman track hits "+CTier) , PloDir, OutFile, MyHists );
1097  // Number of hits in Kal track
1098  MyHists.clear(); for (unsigned int sh=0; sh<NPlot; ++sh) { MyHists.push_back( AllTheHists[sh].hNumKalTrHits ); }
1099  MakeSplitCans( "NumKalTrHits" , CTier , ("Number of hits in kalman track "+CTier), PloDir, OutFile, MyHists );
1100  // Ratio of Kal hits in slice
1101  MyHists.clear(); for (unsigned int sh=0; sh<NPlot; ++sh) { MyHists.push_back( AllTheHists[sh].hRatKalHitSlc ); }
1102  MakeSplitCans( "RatKalHitSlc" , CTier , ("Ratio of kalman tr hit in slc "+CTier), PloDir, OutFile, MyHists );
1103  // Length of Kalman track
1104  MyHists.clear(); for (unsigned int sh=0; sh<NPlot; ++sh) { MyHists.push_back( AllTheHists[sh].hKalTrLength ); }
1105  MakeSplitCans( "KalTrLength" , CTier , ("Length of kalman track "+CTier) , PloDir, OutFile, MyHists );
1106  // Angle Kal track & beam
1107  MyHists.clear(); for (unsigned int sh=0; sh<NPlot; ++sh) { MyHists.push_back( AllTheHists[sh].hKalTrBeamAng ); }
1108  MakeSplitCans( "KalTrBeamAng" , CTier , ("Angle between kal track & beam "+CTier), PloDir, OutFile, MyHists );
1109  // Most forward Kal tr hit
1110  MyHists.clear(); for (unsigned int sh=0; sh<NPlot; ++sh) { MyHists.push_back( AllTheHists[sh].hKalMostFwdCel ); }
1111  MakeSplitCans( "KalMostFwdCel", CTier , ("Most forward kalman track hit "+CTier) , PloDir, OutFile, MyHists );
1112  // Most backward Kal tr hit
1113  MyHists.clear(); for (unsigned int sh=0; sh<NPlot; ++sh) { MyHists.push_back( AllTheHists[sh].hKalMostBakCel ); }
1114  MakeSplitCans( "KalMostBakCel", CTier , ("Most backward kalman track hit "+CTier), PloDir, OutFile, MyHists );
1115  // Highest Y Kalman tr vertex
1116  MyHists.clear(); for (unsigned int sh=0; sh<NPlot; ++sh) { MyHists.push_back( AllTheHists[sh].hKalTrVer_MaxY ); }
1117  MakeSplitCans( "KalTrVer_MaxY" , CTier , ("Highest kalman track vertex Y "+CTier), PloDir, OutFile, MyHists );
1118  // Highest Z Kalman tr vertex
1119  MyHists.clear(); for (unsigned int sh=0; sh<NPlot; ++sh) { MyHists.push_back( AllTheHists[sh].hKalTrVer_MaxZ ); }
1120  MakeSplitCans( "KalTrVer_MaxZ" , CTier , ("Highest kalman track vertex Z "+CTier), PloDir, OutFile, MyHists );
1121  // Y direction of Kal track
1122  MyHists.clear(); for (unsigned int sh=0; sh<NPlot; ++sh) { MyHists.push_back( AllTheHists[sh].hKalTrDir_Y ); }
1123  MakeSplitCans( "KalTrDir_Y" , CTier , ("Y direction of kalman track "+CTier) , PloDir, OutFile, MyHists );
1124  // Scattering over tr length
1125  MyHists.clear(); for (unsigned int sh=0; sh<NPlot; ++sh) { MyHists.push_back( AllTheHists[sh].hScattKalTrLen ); }
1126  MakeSplitCans( "ScattKalTrLen" , CTier , ("Sum of scattering over track "+CTier) , PloDir, OutFile, MyHists );
1127 
1128  // *** Hit based ones...
1129  // First cell that is hit
1130  MyHists.clear(); for (unsigned int sh=0; sh<NPlot; ++sh) { MyHists.push_back( AllTheHists[sh].hFirstHitCell ); }
1131  MakeSplitCans( "FirstHitCell" , CTier , ("First hit in the cell "+CTier) , PloDir, OutFile, MyHists );
1132  // Last cell that is hit
1133  MyHists.clear(); for (unsigned int sh=0; sh<NPlot; ++sh) { MyHists.push_back( AllTheHists[sh].hLastHitCell ); }
1134  MakeSplitCans( "LastHitCell" , CTier , ("Last hit in the cell "+CTier) , PloDir, OutFile, MyHists );
1135  // Maximum Y pos of activty
1136  MyHists.clear(); for (unsigned int sh=0; sh<NPlot; ++sh) { MyHists.push_back( AllTheHists[sh].hMaxActivity_Y ); }
1137  MakeSplitCans( "MaxActivity_Y", CTier , ("Maximum Y position of activity "+CTier) , PloDir, OutFile, MyHists );
1138  // Minimum Y pos of activty
1139  MyHists.clear(); for (unsigned int sh=0; sh<NPlot; ++sh) { MyHists.push_back( AllTheHists[sh].hMinActivity_Y ); }
1140  MakeSplitCans( "MinActivity_Y", CTier , ("Minimum Y position of activity "+CTier) , PloDir, OutFile, MyHists );
1141  // Maximum Z pos of activty
1142  MyHists.clear(); for (unsigned int sh=0; sh<NPlot; ++sh) { MyHists.push_back( AllTheHists[sh].hMaxActivity_Z ); }
1143  MakeSplitCans( "MaxActivity_Z" , CTier , ("Maximum Z position of activity "+CTier), PloDir, OutFile, MyHists );
1144  // Minimum Z pos of activty
1145  MyHists.clear(); for (unsigned int sh=0; sh<NPlot; ++sh) { MyHists.push_back( AllTheHists[sh].hMinActivity_Z ); }
1146  MakeSplitCans( "MinActivity_Z" , CTier , ("Minimum Y position of activity "+CTier), PloDir, OutFile, MyHists );
1147  // Minimum cells to edge
1148  MyHists.clear(); for (unsigned int sh=0; sh<NPlot; ++sh) { MyHists.push_back( AllTheHists[sh].hMinCellToEdge ); }
1149  MakeSplitCans( "MinCellToEdge" , CTier , ("Minmum number of cell to edge "+CTier), PloDir, OutFile, MyHists );
1150 
1151  // *** RemID score / inputs
1152  // The RemID scattering over track length
1153  MyHists.clear(); for (unsigned int sh=0; sh<NPlot; ++sh) { MyHists.push_back( AllTheHists[sh].hRemIDScatLLH ); }
1154  MakeSplitCans( "RemIDScatLLH" , CTier , ("The RemID scattering over tr l "+CTier), PloDir, OutFile, MyHists );
1155  // The RemID dE/dx log-likelihood
1156  MyHists.clear(); for (unsigned int sh=0; sh<NPlot; ++sh) { MyHists.push_back( AllTheHists[sh].hRemIDdEdxLLH ); }
1157  MakeSplitCans( "RemIDdEdxLLH" , CTier , ("The RemID dE/dx "+CTier) , PloDir, OutFile, MyHists );
1158  // The RemID fraction of planes which were measured.
1159  MyHists.clear(); for (unsigned int sh=0; sh<NPlot; ++sh) { MyHists.push_back( AllTheHists[sh].hRemIDMeasFrac ); }
1160  MakeSplitCans( "RemIDMeasFrac", CTier , ("The RemID frac of planes hit "+CTier) , PloDir, OutFile, MyHists );
1161 
1162  // *** Some Calibration plots
1163  for (int cor=0; cor<4; ++cor) {
1164  // Reco Energy
1165  MyHists.clear(); for (unsigned int sh=0; sh<NPlot; ++sh) { MyHists.push_back( AllTheHists[sh].hCor_RecEn[cor] ); }
1166  MakeSplitCans( "Cor_RecEn"+CornName[cor], CTier , ("RecoE_"+CTier+"_"+CornName[cor]), PloDir, OutFile, MyHists );
1167  // Muon Energy
1168  MyHists.clear(); for (unsigned int sh=0; sh<NPlot; ++sh) { MyHists.push_back( AllTheHists[sh].hCor_MuoEn[cor] ); }
1169  MakeSplitCans( "Cor_MuoEn"+CornName[cor], CTier , ("MuonE_"+CTier+"_"+CornName[cor]), PloDir, OutFile, MyHists );
1170  // Hadronic Energy
1171  MyHists.clear(); for (unsigned int sh=0; sh<NPlot; ++sh) { MyHists.push_back( AllTheHists[sh].hCor_HadEn[cor] ); }
1172  MakeSplitCans( "Cor_HadEn"+CornName[cor], CTier , ("HadE_" +CTier+"_"+CornName[cor]), PloDir, OutFile, MyHists );
1173  // Hadronic Energy Fraction
1174  MyHists.clear(); for (unsigned int sh=0; sh<NPlot; ++sh) { MyHists.push_back( AllTheHists[sh].hCor_HFrEn[cor] ); }
1175  MakeSplitCans( "Cor_HFrEn"+CornName[cor], CTier , ("HadEF_"+CTier+"_"+CornName[cor]), PloDir, OutFile, MyHists );
1176  }
1177  //*/
1178  return;
1179 }
1180 
1181 //......................................................
TH1D * hHadroEnPerHit
Definition: DoThePlots.C:41
tree Draw("slc.nhit")
TH1D * hKalTrVer_MaxZ
Definition: MakeThePlots.C:79
TH1D * hKalMostBakCel
Definition: MakeThePlots.C:77
const bool SingleCut
Definition: MakeThePlots.C:136
bin1_2sigma SetFillColor(3)
TH1D * hNumKalTracks
Definition: MakeThePlots.C:71
std::string LegNames[5]
Definition: MakeThePlots.C:133
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
TH1D * GetSpectToHist(TFile *MyF, std::string LoadName, std::string Axis, int FileType)
Definition: MakeThePlots.C:338
const int NumCut
Definition: MakeThePlots.C:113
void ResetOscCalcToDefault(osc::IOscCalcAdjustable *calc)
Reset calculator to default assumptions for all parameters.
Definition: Calcs.cxx:23
Antineutrinos-only.
Definition: IPrediction.h:50
const bool SingleQuant
Definition: MakeThePlots.C:137
correl_xv GetYaxis() -> SetDecimals()
TH1D * hRemIDdEdxLLH
Definition: MakeThePlots.C:94
TH1D * hMuonEnergy
Definition: DoThePlots.C:38
TH1D * hReconstEnergy
Definition: DoThePlots.C:35
TH1D * hMinActivity_Y
Definition: MakeThePlots.C:87
T sqrt(T number)
Definition: d0nt_math.hpp:156
prelim SetTextSize(2/30.)
TH1D * hKalTrVer_MaxY
Definition: MakeThePlots.C:78
TH1D * hHadronEnergy
Definition: DoThePlots.C:40
Definition: Hist.h:29
TH1D * MakeRatio(TH1D *num, TH1D *denom, int Col, std::string FType)
Definition: MakeThePlots.C:442
TH1D * hTrkStartX
Definition: DoThePlots.C:54
TH2 * Rebin(TH2 *hist, std::vector< double > edges)
TH1D * hRemIDMeasFrac
Definition: MakeThePlots.C:95
TH1D * hCor_HadEn[4]
Definition: MakeThePlots.C:100
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)
Definition: MakeThePlots.C:497
void MakeSplitCans(std::string PlotName, std::string CutTier, std::string CanTitle, std::string PlotDir, TFile *OutFile, std::vector< TH1D * > Hists)
Definition: MakeThePlots.C:714
void CenterTitles(TH1 *histo)
Definition: Plots.cxx:1481
TGraphAsymmErrors * GraphWithPoissonErrors(const TH1 *h, bool noErrorsXaxis, bool drawEmptyBins)
Calculate statistical errors appropriate for small Poisson numbers.
Definition: Plots.cxx:910
unsigned int NPlot
Definition: MakeThePlots.C:130
TH1D * hReTrOverTrNuE
Definition: DoThePlots.C:51
TH1D * hTrkStartZ
Definition: DoThePlots.C:56
fVtxDx SetMarkerStyle(20)
double POTNom
Definition: MakeThePlots.C:106
const std::string CutNames[NumCut]
Definition: MakeThePlots.C:114
osc::OscCalcDumb calc
ntuple SetFillStyle(1001)
TH1D * hRemIDScatLLH
Definition: MakeThePlots.C:93
TH1D * hMinCellToEdge
Definition: MakeThePlots.C:90
#define M_PI
Definition: SbMath.h:34
bool AllPer
Definition: MakeThePlots.C:109
double LivNom
Definition: MakeThePlots.C:107
TH1D * hCor_MuoEn[4]
Definition: MakeThePlots.C:99
void PlotHistProp(TH1D *Hist, TLegend *L)
Definition: MakeThePlots.C:677
TH1D * hCVNNuMuID2017
Definition: MakeThePlots.C:54
TH1D * hCor_RecEn[4]
Definition: MakeThePlots.C:98
TH1D * hNumKalTrHits
Definition: MakeThePlots.C:72
void MakeTextFile(std::string Cut, std::string PTitle, std::string Name)
Definition: MakeThePlots.C:696
TH1D * hMaxActivity_Y
Definition: MakeThePlots.C:86
TH1D * hTrueNuEnergy
Definition: DoThePlots.C:50
Charged-current interactions.
Definition: IPrediction.h:39
legend SetBorderSize(0)
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:535
static constexpr double L
correl_xv GetXaxis() -> SetDecimals()
hmean SetLineWidth(2)
TGraph * graphAsymmError(TH1 *hDataScaled, TH1 *hData)
Definition: MakeThePlots.C:412
leg AddEntry(GRdata,"data","p")
hmean SetLineColor(4)
Optimized version of OscCalcPMNS.
Definition: StanTypedefs.h:31
void MakeThePlots(bool Horn, bool FullDataset=true, bool UsingCAFs=false)
Definition: MakeThePlots.C:174
_Cut< caf::SRProxy > Cut
Representation of a cut (selection) to be applied to a caf::StandardRecord object.
Definition: Cut.h:96
std::string PlotDir
void SetRange(std::vector< TH1D * > Histos, std::vector< TH1D * > Ratios, double XLow, double XHigh, double H_YLow, double H_YHigh, double R_YLow, double R_YHigh)
Definition: MakeThePlots.C:468
std::string sPer
Definition: MakeThePlots.C:111
TH1D * hTrkLenXY
Definition: DoThePlots.C:60
std::string sFHC
Definition: MakeThePlots.C:110
const bool PlotMean
Definition: MakeThePlots.C:138
TH1D * hMuonEnPerHit
Definition: DoThePlots.C:39
float bin[41]
Definition: plottest35.C:14
TH1D * hHitsPerSlice
Definition: DoThePlots.C:43
void ProducePlots(std::vector< LoadedHistograms > AllTheHists, std::string CTier)
Definition: MakeThePlots.C:990
TH1D * hCVNCosmicScor
Definition: DoThePlots.C:45
TH1D * hMinActivity_Z
Definition: MakeThePlots.C:89
void Preliminary()
Neutrinos-only.
Definition: IPrediction.h:49
TH1D * hRatKalHitSlc
Definition: MakeThePlots.C:73
OStream cout
Definition: OStream.cxx:6
const Cut cut
Definition: exporter_fd.C:30
TH1D * hCor_HFrEn[4]
Definition: MakeThePlots.C:101
TH1D * hScattKalTrLen
Definition: MakeThePlots.C:81
TH1D * hCVNNuMuIDScor
Definition: DoThePlots.C:46
TH1D * hKalMostFwdCel
Definition: MakeThePlots.C:76
TH1D * hSliceTimeZoom
Definition: DoThePlots.C:37
TH1D * hNuMuContPID
Definition: DoThePlots.C:47
bool IsFHC
Definition: MakeThePlots.C:108
Sign::Sign_t wrongs
Definition: saveFDMCHists.C:27
int num
Definition: f2_nu.C:119
enum BeamMode kViolet
LoadedHistograms LoadFile_GetHists(TFile *InFile, std::string CTier, int FileType)
Definition: MakeThePlots.C:254
TH1D * hFirstHitCell
Definition: MakeThePlots.C:84
const bool SaveHists
Definition: MakeThePlots.C:141
Neutral-current interactions.
Definition: IPrediction.h:40
void CornerLabel(std::string Str)
Definition: MakeThePlots.C:687
TH1D * hMaxActivity_Z
Definition: MakeThePlots.C:88
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
const bool PlotTest
Definition: MakeThePlots.C:139
TH1D * hSANuMuContPID
Definition: DoThePlots.C:48
const bool SaveRatio
Definition: MakeThePlots.C:140
TFile * OutFile
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
All neutrinos, any flavor.
Definition: IPrediction.h:26
(&#39; appearance&#39;)
Definition: IPrediction.h:16
TH1D * hHadFracEnergy
Definition: DoThePlots.C:42
TH1D * hRemIDScore
Definition: DoThePlots.C:44
cosmicTree SaveAs("cosmicTree.root")
static std::unique_ptr< PredictionNoExtrap > LoadFrom(TDirectory *dir, const std::string &name)
const std::string QuantNames[NumQuant]
Definition: MakeThePlots.C:122
enum BeamMode kGreen
simulatedPeak Scale(1/simulationNormalisationFactor)
c cd(1)
const int NumQuant
Definition: MakeThePlots.C:121
const std::string CornName[4]
Definition: MakeThePlots.C:144
TH1D * hSliceTimeFull
Definition: DoThePlots.C:36
gPad SetLogy()
TH1D * hKalTrBeamAng
Definition: DoThePlots.C:47
T asin(T number)
Definition: d0nt_math.hpp:60
gm Write()
TH1D * hTrkStartY
Definition: DoThePlots.C:55
const int Colours[5]
Definition: MakeThePlots.C:132
enum BeamMode string