DoThePlots.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 "Utilities/rootlogon.C"
20 
21 #include <iostream>
22 #include <sstream>
23 #include <fstream>
24 #include <vector>
25 #include <string>
26 
27 #include <TROOT.h>
28 #include <TStyle.h>
29 
30 using namespace ana;
31 
32 //......................................................
34  // *** Some key ones first.
38  TH1D* hMuonEnergy ;
39  TH1D* hMuonEnPerHit ;
40  TH1D* hHadronEnergy ;
43  TH1D* hHitsPerSlice ;
44  TH1D* hRemIDScore ;
47  TH1D* hNuMuContPID ;
49  // *** Some important Truth based ones.
50  TH1D* hTrueNuEnergy ;
52 
53  // *** Basic positional plots.
54  TH1D* hTrkStartX;
55  TH1D* hTrkStartY;
56  TH1D* hTrkStartZ;
57  TH1D* hTrkEndX ;
58  TH1D* hTrkEndY ;
59  TH1D* hTrkEndZ ;
60  TH1D* hTrkLenXY ;
61  /*
62  // *** Track based ones...
63  TH1D* hNumKalTracks ;
64  TH1D* hNumKalTrHits ;
65  TH1D* hRatKalHitSlc ;
66  TH1D* hKalTrLength ;
67  TH1D* hKalTrBeamAng ;
68  TH1D* hKalMostFwdCel;
69  TH1D* hKalMostBakCel;
70  TH1D* hKalTrVer_MaxY;
71  TH1D* hKalTrVer_MaxZ;
72  TH1D* hKalTrDir_Y ;
73  TH1D* hScattKalTrLen;
74 
75  // *** Hit based ones...
76  TH1D* hFirstHitCell ;
77  TH1D* hLastHitCell ;
78  TH1D* hMaxActivity_Y;
79  TH1D* hMinActivity_Y;
80  TH1D* hMaxActivity_Z;
81  TH1D* hMinActivity_Z;
82  TH1D* hMinCellToEdge;
83 
84  // *** Visible energy quanta.
85  // Slice
86  TH1D* hVisEInSlice ;
87  TH1D* hNHitInSlice ;
88  TH1D* hVisEPerHitSlc;
89  // Track
90  TH1D* hVisEInTrack ;
91  TH1D* hNHitInTrack ;
92  TH1D* hVisEPerHitTrk;
93  // Hadronic
94  TH1D* hVisEInHadron ;
95  TH1D* hNHitInHadron ;
96  TH1D* hVisEPerHitHad;
97 
98  // *** RemID score / inputs
99  TH1D* hRemIDScatLLH ;
100  TH1D* hRemIDdEdxLLH ;
101  TH1D* hRemIDMeasFrac;
102 
103  // *** Some calibration quantities
104  TH1D* hCor_RecEn[4];
105  TH1D* hCor_MuoEn[4];
106  TH1D* hCor_HadEn[4];
107  TH1D* hCor_HFrEn[4];
108  TH1D* hCor_HFrEn[4];
109  */
110 };
111 
112 //......................................................
114  // *** Some key ones first.
115  TCanvas* cReconstEnergy;
116  TCanvas* cSliceTimeFull;
117  TCanvas* cSliceTimeZoom;
118  TCanvas* cMuonEnergy ;
119  TCanvas* cMuonEnPerHit ;
120  TCanvas* cHadronEnergy ;
121  TCanvas* cHadroEnPerHit;
122  TCanvas* cHadFracEnergy;
123  TCanvas* cHitsPerSlice ;
124  TCanvas* cRemIDScore ;
125  TCanvas* cCVNCosmicScor;
126  TCanvas* cCVNNuMuIDScor;
127  TCanvas* cNuMuContPID ;
128  TCanvas* cSANuMuContPID;
129  // *** Some important Truth based ones.
130  TCanvas* cTrueNuEnergy ;
131  TCanvas* cReTrOverTrNuE;
132 
133  // *** Basic positional plots.
134  TCanvas* cTrkStartX;
135  TCanvas* cTrkStartY;
136  TCanvas* cTrkStartZ;
137  TCanvas* cTrkEndX ;
138  TCanvas* cTrkEndY ;
139  TCanvas* cTrkEndZ ;
140  TCanvas* cTrkLenXY ;
141  /*
142  // *** Track based ones...
143  TCanvas* cNumKalTracks ;
144  TCanvas* cNumKalTrHits ;
145  TCanvas* cRatKalHitSlc ;
146  TCanvas* cKalTrLength ;
147  TCanvas* cKalTrBeamAng ;
148  TCanvas* cKalMostFwdCel;
149  TCanvas* cKalMostBakCel;
150  TCanvas* cKalTrVer_MaxY;
151  TCanvas* cKalTrVer_MaxZ;
152  TCanvas* cKalTrDir_Y ;
153  TCanvas* cScattKalTrLen;
154 
155  // *** Hit based ones...
156  TCanvas* cFirstHitCell ;
157  TCanvas* cLastHitCell ;
158  TCanvas* cMaxActivity_Y;
159  TCanvas* cMinActivity_Y;
160  TCanvas* cMaxActivity_Z;
161  TCanvas* cMinActivity_Z;
162  TCanvas* cMinCellToEdge;
163 
164  // *** Visible energy quanta.
165  // Slice
166  TCanvas* cVisEInSlice ;
167  TCanvas* cNHitInSlice ;
168  TCanvas* cVisEPerHitSlc;
169  // Track
170  TCanvas* cVisEInTrack ;
171  TCanvas* cNHitInTrack ;
172  TCanvas* cVisEPerHitTrk;
173  // Hadronic
174  TCanvas* cVisEInHadron ;
175  TCanvas* cNHitInHadron ;
176  TCanvas* cVisEPerHitHad;
177 
178  // *** RemID score / inputs
179  TCanvas* cRemIDScatLLH ;
180  TCanvas* cRemIDdEdxLLH ;
181  TCanvas* cRemIDMeasFrac;
182 
183  // *** Some calibration quantities
184  TCanvas* cCor_RecEn[4];
185  TCanvas* cCor_MuoEn[4];
186  TCanvas* cCor_HadEn[4];
187  TCanvas* cCor_HFrEn[4];
188  TCanvas* cCor_HFrEn[4];
189  */
190 };
191 
192 //......................................................
193 const int Colours[4] = { 2, 1, 4, 6 };
194 std::string LegNames[4] = { "Nom", "Sh1", "Sh2", "Sh3"};
197 // --- Cuts
198 const int NumCut = 4;
200  "_QualCuts",
201  "_DetCuts" ,
202  "_CosmCuts",
203  "_NeuCCuts"
204 };
205 // --- Quantiles
206 const int NumQuant = 5;
208  "_AllQuants",
209  "_Quant1",
210  "_Quant2",
211  "_Quant3",
212  "_Quant4"
213 };
214 // --- How many systs do I have?
215 unsigned int NumSysts = 0;
216 // --- Some other things.
218 const bool SingleCut = true;
219 const bool SingleQuant = true;
220 const bool PlotMean = false;
221 const bool PlotTest = false;
222 // --- Use Nominal sample for light?
223 const bool UseLightNominal = false;
224 
225 // --- Calibration labels.
226 const std::string CornName[4] = { "TopL", "TopR", "BotL", "BotR" };
227 
228 //......................................................
229 std::vector<LoadedHistograms> LoadFile_GetHists( std::vector<std::string> FileLocs, std::string CTier );
230 //......................................................
231 TH1D* GetSpectToHist( TFile* MyF, std::string HNam, std::string CTier, std::string Axis );
232 //......................................................
233 TH1D* MakeRatio( TH1D* num, TH1D* denom, int Col, std::string FType );
234 //......................................................
235 void SetRange( TH1D* H1, std::vector<TH1D*> Shift, std::vector<TH1D*> Ratio, // Pass the histograms
236  double XLow, double XHigh, // What X values do I want?
237  double H_YLow, double H_YHigh, double R_YLow, double R_YHigh ); // What Y values do I want?
238 //......................................................
239 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 );
240 void FindAxisRange_DataMC( double &XLow, double &XHigh, double &H_YLow, double &H_YHigh, double &R_YLow, double &R_YHigh, double MaxVal, bool &SetLogy, std::string Name );
241 //......................................................
242 void MakeSplitCans( TCanvas* Can, std::string DirNa, std::string CanTi, std::string FType, std::string CTier, TFile* OutF, TH1D* Nom, std::vector<TH1D*> Shi);
243 //......................................................
244 void ProducePlots( LoadedHistograms NomHists, std::vector<LoadedHistograms> SubSysts, std::string FType, std::string CTier );
245 //......................................................
246 void PlotHistProp( TH1D* Hist, TLegend* L );
247 //......................................................
248 void Preliminary();
249 //......................................................
250 void MakeTextFile( std::string Cut, std::string PTitle, std::string Name );
251 //......................................................
252 void DoThePlots() {
253  std::cout << "\n\nPlease specify a systematic, eg DataMC, LightShift, CalibShift\n\n" << std::endl;
254  return;
255 }
256 
257 //......................................................
258 void DoThePlots( std::string WhichSyst ) {
259  // ---- First off, lets set the styles...
260  gStyle->SetOptStat(0);
261  gROOT->SetStyle("novaStyle");
262 
263 
264  // --- And set MySys
265  MySys = WhichSyst;
266 
267  std::string BaseCAFFileDir="/nova/ana/users/karlwarb/Systematic_2017/";
268  if (MySys == "DataMC")
269  BaseCAFFileDir="/nova/ana/users/karlwarb/NuMu_BoxOpening2017/Comparison/";
270 
271  // --- Declare a vector of my file name structures
272  std::vector<std::string> NomNames;
273 
274  if (MySys == "LightShift") {
275  // Push back NomNames.
276  NomNames.push_back( "FD_NonSwap_AP" );
277  NomNames.push_back( "ND_NonSwap_AP" );
278  // Set NumSysts
279  NumSysts = 3;
280  // Legend Names
281  LegNames[0] = "Nominal"; if (!UseLightNominal) LegNames[0] = "No Shift";
282  LegNames[1] = "Light Up, Calib Down";
283  LegNames[2] = "Light Down, Calib Up";
284  LegNames[3] = "No shift" ; if (!UseLightNominal) LegNames[3] = "No Cherenkov variation";
285  }
286  else if (MySys == "CalibShift") {
287  // Push back NomNames.
288  NomNames.push_back( "FD_NonSwap_AP" );
289  NomNames.push_back( "ND_NonSwap_AP" );
290  // Set NumSysts
291  NumSysts = 3;
292  // Legend Names
293  LegNames[0] = "Nominal";
294  LegNames[1] = "Negative offset";
295  LegNames[2] = "Positive offset";
296  LegNames[3] = "Shape variation";
297  }
298  else if (MySys == "DataMC") {
299  // Push back NomNames.
300  NomNames.push_back( "Short_AP" );
301  //NomNames.push_back( "P1" );
302  //NomNames.push_back( "P2" );
303  //NomNames.push_back( "P3" );
304  //NomNames.push_back( "P5" );
305  // Set NumSysts
306  NumSysts = 2;
307  // Legend Names
308  LegNames[0] = "Prediction";
309  LegNames[1] = "FD data";
310  LegNames[2] = "Cosmic Bkg";
311  }
312 
313  // --- I need to figure out POT and Livetime for Data...
314  if (MySys == "DataMC") {
315  std::string DataFileLoc = BaseCAFFileDir+"Concat_Data_"+NomNames[0]+".root";
316  TFile *TempDataFile = TFile::Open( DataFileLoc.c_str() );
317  std::unique_ptr<Spectrum> TempDataSpec = Spectrum::LoadFrom(TempDataFile, "sReconstEnergy_NeuCCuts_AllQuants" ) ;
318  POTNom = TempDataSpec->POT();
319  std::string SideFileLoc = BaseCAFFileDir+"Concat_Cosmics_"+NomNames[0]+".root";
320  TFile *TempSideFile = TFile::Open( SideFileLoc.c_str() );
321  std::unique_ptr<Spectrum> TempSideSpec = Spectrum::LoadFrom(TempSideFile, "sReconstEnergy_NeuCCuts_AllQuants" ) ;
322  LivNom = TempSideSpec->Livetime();
323  }
324 
325  std::cout << "\n\nMy Livetime is now " << POTNom << " ("<<kAna2017POT<<") Exposure is now " << LivNom << " ("<<kAna2017Livetime<<")\n\n" << std::endl;
326 
327  // --- Declare a vector of nominal file names, and fill them from NomNames
328  std::vector<std::string> NominalFiles;
329  for (size_t ff=0; ff<NomNames.size(); ++ff) {
330  if (MySys == "LightShift" && !UseLightNominal ) {
331  NominalFiles.push_back( BaseCAFFileDir+"Output_"+NomNames[ff]+"_"+MySys+"_NoShift.root" );
332  }
333  else if (MySys == "CalibShift") {
334  NominalFiles.push_back( BaseCAFFileDir+"Output_"+NomNames[ff]+"_"+MySys+"_Nominal.root" );
335  }
336  else if (MySys == "DataMC") {
337  NominalFiles.push_back( BaseCAFFileDir+"Concat_MC_"+NomNames[ff]+".root" );
338  }
339  }
340  // --- Declare a vector of shifted files for each nominal file...
341  std::vector<std::string> SystematFiles;
342  for (size_t ff=0; ff<NomNames.size(); ++ff) {
343  // Light levels shifts.
344  if (MySys == "LightShift" ) {
345  SystematFiles.push_back( BaseCAFFileDir+"Output_"+NomNames[ff]+"_"+MySys+"_LightUp_CalibDown.root" );
346  SystematFiles.push_back( BaseCAFFileDir+"Output_"+NomNames[ff]+"_"+MySys+"_LightDown_CalibUp.root" );
347  if ( UseLightNominal ) {
348  SystematFiles.push_back( BaseCAFFileDir+"Output_"+NomNames[ff]+"_"+MySys+"_NoShift.root" );
349  } else {
350  SystematFiles.push_back( BaseCAFFileDir+"Output_"+NomNames[ff]+"_"+MySys+"_CherenkovVariance.root" );
351  }
352  } else if (MySys == "CalibShift") {
353  // Calibration shifts.
354  SystematFiles.push_back( BaseCAFFileDir+"Output_"+NomNames[ff]+"_"+MySys+"_XYView_NegShift.root" );
355  SystematFiles.push_back( BaseCAFFileDir+"Output_"+NomNames[ff]+"_"+MySys+"_XYView_PosShift.root" );
356  SystematFiles.push_back( BaseCAFFileDir+"Output_"+NomNames[ff]+"_"+MySys+"_ShapeDifference.root" );
357  }
358  else if (MySys == "DataMC") {
359  // DataMC comparisons.
360  //SystematFiles.push_back( BaseCAFFileDir+"Output_FD_NuMi_"+NomNames[ff]+"_Data.root" );
361  SystematFiles.push_back( BaseCAFFileDir+"Concat_Data_"+NomNames[ff]+".root" );
362  SystematFiles.push_back( BaseCAFFileDir+"Concat_Cosmics_"+NomNames[ff]+".root" );
363  }
364  }
365 
366  for (size_t qq=0; qq<NominalFiles.size(); ++qq) {
367  std::cout << "\nFor qq = " << qq << " I have;" << "\n\tnominal file; " << NominalFiles[qq];
368  for (unsigned int sh=0; sh<NumSysts; ++sh) {
369  std::cout << "\n\tFile " << sh << " is" << SystematFiles[qq*sh];
370  }
371  std::cout << "." << std::endl;
372  }
373 
374  for (int quant=0; quant<NumQuant; ++quant) {
375  if (SingleQuant && quant!=0) continue;
376  for (int cut=0; cut<NumCut; ++cut) {
377  // Check if just want to run a single cut?
378  //if (SingleCut && cut != 0) continue;
379  if (SingleCut && cut != 3) continue;
380 
381  // --- What is the Cut + Quant combination?
382  std::string CutQuant = CutNames[cut] + QuantNames[quant];
383 
384  // Now loop through the files and load the histograms.
385  std::cout << "\nMade my vector of Nominal files. I have a total of " << NominalFiles.size() << " files." << std::endl;
386  std::vector<LoadedHistograms> NominalHists = LoadFile_GetHists( NominalFiles , CutQuant );
387 
388  // Now loop through the files and load the histograms.
389  std::cout << "\nMade my vector of Systematic Files files. I have a total of " << SystematFiles.size() << " files." << std::endl;
390  std::vector<LoadedHistograms> SystematHists = LoadFile_GetHists( SystematFiles, CutQuant );
391 
392  for (size_t NomL=0; NomL<NominalHists.size(); ++NomL) {
393  // --- What do I want my Name to be?
394  std::string LegName = "";
395  if (MySys == "LightShift" ) { LegName = NomNames[NomL]+"_Light"; }
396  else if (MySys == "CalibShift" ) { LegName = NomNames[NomL]+"_Shape"; }
397  else if (MySys == "DataMC" ) { LegName = NomNames[NomL]+"_DataMC";}
398 
399  std::cout << "\nI will be loading things from NomL " << NomL << std::endl;
400  std::vector<LoadedHistograms> SubSysts;
401  for (unsigned int sh=0; sh<NumSysts; ++sh) {
402  std::cout << "\t Will load element " << (NumSysts*NomL)+sh << std::endl;
403  SubSysts.push_back(SystematHists[ (NumSysts*NomL)+sh ]);
404  }
405  std::cout << "My subvector has size " << SubSysts.size() << "\n\n" << std::endl;
406 
407  // --- Make my plots
408  ProducePlots( NominalHists[NomL], SubSysts, LegName, CutQuant );
409  } // Loop through nominal vector.
410  } // Loop through cuts.
411  } // Loop through quantiles
412  // --- Finally, return.
413  return;
414 } // DoThePlots
415 
416 //......................................................
417 std::vector<LoadedHistograms> LoadFile_GetHists( std::vector<std::string> FileLocs, std::string CutTier ) {
418  // --- Declare the vector of LoadedHistograms which I am going to return...
419  std::vector<LoadedHistograms> MyHists;
420  // --- Loop through the file list...
421  for (size_t FileL=0; FileL < FileLocs.size(); ++FileL) {
422  std::cerr << "\tLooking at entry " << FileL << " -- " << FileLocs[FileL] << std::endl;
423  // --- Make a new instance of LoadedHistograms.
424  LoadedHistograms TempHists;
425  // --- Load the current file.
426  TFile *File1 = TFile::Open( FileLocs[FileL].c_str() );
427  // --- Grab my spectra...
428  // *** Some key ones first.
429  TempHists.hReconstEnergy = GetSpectToHist( File1, "sReconstEnergy", CutTier, ";Reconstructed #nu_{#mu} energy in the slice (GeV); Number of events" );
430  TempHists.hSliceTimeFull = GetSpectToHist( File1, "sSliceTimeFull", CutTier, ";Mean time of the slice over full window (ns); Number of events" );
431  TempHists.hSliceTimeZoom = GetSpectToHist( File1, "sSliceTimeZoom", CutTier, ";Mean time of the slice over beam window (ns); Number of events" );
432  TempHists.hMuonEnergy = GetSpectToHist( File1, "sMuonEnergy" , CutTier, ";Reconstructed muon energy (GeV); Number of events" );
433  TempHists.hMuonEnPerHit = GetSpectToHist( File1, "sMuonEnPerHit" , CutTier, ";Energy per hit of the muon (GeV); Number of events" );
434  TempHists.hHadronEnergy = GetSpectToHist( File1, "sHadronEnergy" , CutTier, ";Hadronic energy in the slice (GeV); Number of events" );
435  TempHists.hHadroEnPerHit = GetSpectToHist( File1, "sHadroEnPerHit", CutTier, ";Hadronic energy per hit in the slice (GeV); Number of events" );
436  TempHists.hHadFracEnergy = GetSpectToHist( File1, "sHadFracEnergy", CutTier, ";Hadronic energy fraction in the slice; Number of events" );
437  TempHists.hHitsPerSlice = GetSpectToHist( File1, "sHitsPerSlice" , CutTier, ";Mumber of hits in the slice; Number of events" );
438  TempHists.hRemIDScore = GetSpectToHist( File1, "sRemIDScore" , CutTier, ";RemID score; Number of events" );
439  TempHists.hCVNCosmicScor = GetSpectToHist( File1, "sCVNCosmicScor", CutTier, ";CVN Cosmic score; Number of events" );
440  TempHists.hCVNNuMuIDScor = GetSpectToHist( File1, "sCVNNuMuIDScor", CutTier, ";CVN #nu_{#mu} score; Number of events" );
441  TempHists.hNuMuContPID = GetSpectToHist( File1, "sNuMuContPID" , CutTier, ";Third analysis BDT score; Number of events" );
442  TempHists.hSANuMuContPID = GetSpectToHist( File1, "sSANuMuContPID", CutTier, ";Second analysis BDT score; Number of events" );
443  // *** Some important Truth based ones.
444  TempHists.hTrueNuEnergy = GetSpectToHist( File1, "sTrueNuEnergy" , CutTier, ";True neutrino energy in the slice (GeV); Number of events" );
445  TempHists.hReTrOverTrNuE = GetSpectToHist( File1, "sReTrOverTrNuE", CutTier, ";(Reco-True)/True neutrino energy in the slice (GeV); Number of events");
446 
447  // *** Basic positional plots.
448  TempHists.hTrkStartX = GetSpectToHist( File1, "sTrkStartX", CutTier, ";Initial X position of the muon track (m); Number of events" );
449  TempHists.hTrkStartY = GetSpectToHist( File1, "sTrkStartY", CutTier, ";Initial Y position of the muon track (m); Number of events" );
450  TempHists.hTrkStartZ = GetSpectToHist( File1, "sTrkStartZ", CutTier, ";Initial Z position of the muon track (m); Number of events" );
451  TempHists.hTrkEndX = GetSpectToHist( File1, "sTrkEndX" , CutTier, ";Initial X position of the muon track (m); Number of events" );
452  TempHists.hTrkEndY = GetSpectToHist( File1, "sTrkEndY" , CutTier, ";Initial Y position of the muon track (m); Number of events" );
453  TempHists.hTrkEndZ = GetSpectToHist( File1, "sTrkEndZ" , CutTier, ";Initial Z position of the muon track (m); Number of events" );
454  TempHists.hTrkLenXY = GetSpectToHist( File1, "sTrkLenXY" , CutTier, ";Muon track length in the XY plane (m); Number of events" );
455  /*
456  // *** Track based ones...
457  TempHists.hNumKalTracks = GetSpectToHist( File1, "sNumKalTracks" , CutTier, ";Number of kalman tracks; Number of events" );
458  TempHists.hNumKalTrHits = GetSpectToHist( File1, "sNumKalTrHits" , CutTier, ";Number of hits in kalman trakcs; Number of events" );
459  TempHists.hRatKalHitSlc = GetSpectToHist( File1, "sRatKalHitSlc" , CutTier, ";Ratio of kalman track hits to total hits; Number of events" );
460  TempHists.hKalTrLength = GetSpectToHist( File1, "sKalTrLength" , CutTier, ";Kalman track length (m); Number of events" );
461  TempHists.hKalTrBeamAng = GetSpectToHist( File1, "sKalTrBeamAng" , CutTier, ";#cos(#theta) between kalman track and the beam; Number of events");
462  TempHists.hKalMostFwdCel = GetSpectToHist( File1, "sKalMostFwdCel", CutTier, ";Most forward cell in kalman track; Number of events" );
463  TempHists.hKalMostBakCel = GetSpectToHist( File1, "sKalMostBakCel", CutTier, ";Most backward cell in kalman track; Number of events" );
464  TempHists.hKalTrVer_MaxY = GetSpectToHist( File1, "sKalTrVer_MaxY", CutTier, ";Maximal Y vertex position in a slice (m); Number of events" );
465  TempHists.hKalTrVer_MaxZ = GetSpectToHist( File1, "sKalTrVer_MaxZ", CutTier, ";Maximal Z vertex position in a slice (m); Number of events" );
466  TempHists.hKalTrDir_Y = GetSpectToHist( File1, "sKalTrDir_Y" , CutTier, ";Kalman track direction in Y; Direction in Y; Number of events" );
467  TempHists.hScattKalTrLen = GetSpectToHist( File1, "sScattKalTrLen", CutTier, ";Sum of angular scattering over track length; Number of events" );
468 
469  // *** Hit based ones...
470  TempHists.hFirstHitCell = GetSpectToHist( File1, "sFirstHitCell" , CutTier, ";First cell hit in the slice; Number of events" );
471  TempHists.hLastHitCell = GetSpectToHist( File1, "sLastHitCell" , CutTier, ";Last cell hit in the slice; Number of events" );
472  TempHists.hMaxActivity_Y = GetSpectToHist( File1, "sMaxActivity_Y", CutTier, ";Maximum position of activity in Y (m); Number of events" );
473  TempHists.hMinActivity_Y = GetSpectToHist( File1, "sMinActivity_Y", CutTier, ";Minimum position of activity in Y (m); Number of events" );
474  TempHists.hMaxActivity_Z = GetSpectToHist( File1, "sMaxActivity_Z", CutTier, ";Maximum position of activity in Z (m); Number of events" );
475  TempHists.hMinActivity_Z = GetSpectToHist( File1, "sMinActivity_Z", CutTier, ";Minimum position of activity in Z (m); Number of events" );
476  TempHists.hMinCellToEdge = GetSpectToHist( File1, "sMinCellToEdge", CutTier, ";Minimum number of cells hit to detector edge; Number of events" );
477 
478  // *** Visible energy quanta.
479  // Slice
480  TempHists.hVisEInSlice = GetSpectToHist( File1, "sVisEInSlice" , CutTier, ";Visible energy in slice (GeV); Number of events" );
481  TempHists.hNHitInSlice = GetSpectToHist( File1, "sNHitInSlice" , CutTier, ";Number of hits in slice; Number of events" );
482  TempHists.hVisEPerHitSlc = GetSpectToHist( File1, "sVisEPerHitSlc", CutTier, ";Visible energy per hit in slice (GeV); Number of events" );
483  // Track
484  TempHists.hVisEInTrack = GetSpectToHist( File1, "sVisEInTrack" , CutTier, ";Visible energy in track (GeV); Number of events" );
485  TempHists.hNHitInTrack = GetSpectToHist( File1, "sNHitInTrack" , CutTier, ";Number of hits in track; Number of events" );
486  TempHists.hVisEPerHitTrk = GetSpectToHist( File1, "sVisEPerHitTrk", CutTier, ";Visibile energy per hit in track (GeV); Number of events" );
487  // Hadronic
488  TempHists.hVisEInHadron = GetSpectToHist( File1, "sVisEInHadron" , CutTier, ";Visible hadronic energy (GeV); Number of events" );
489  TempHists.hNHitInHadron = GetSpectToHist( File1, "sNHitInHadron" , CutTier, ";Number of hadronic hits; Number of events" );
490  TempHists.hVisEPerHitHad = GetSpectToHist( File1, "sVisEPerHitHad", CutTier, ";Visible hadronic energy per hit (GeV); Number of events" );
491 
492  // *** RemID score / inputs
493  TempHists.hRemIDScatLLH = GetSpectToHist( File1, "sRemIDScatLLH" , CutTier, ";Log-likelihood of scattering; Number of events" );
494  TempHists.hRemIDdEdxLLH = GetSpectToHist( File1, "sRemIDdEdxLLH" , CutTier, ";Log-likelihood of dEdx; Number of events" );
495  TempHists.hRemIDMeasFrac = GetSpectToHist( File1, "sRemIDMeasFrac", CutTier, ";Log-likelihood of RemID; Number of events" );
496 
497  // *** Some Calibration plots
498  for (int cor=0; cor<4; ++cor) {
499  TempHists.hCor_RecEn[cor] = GetSpectToHist( File1, "sCor_RecEn", CutTier+"_"+CornName[cor], ";Reconstructed En. for corner "+CornName[cor]+" (GeV); Number of events");
500  TempHists.hCor_MuoEn[cor] = GetSpectToHist( File1, "sCor_MuoEn", CutTier+"_"+CornName[cor], ";Muon En. for corner "+CornName[cor]+" (GeV); Number of events" );
501  TempHists.hCor_HadEn[cor] = GetSpectToHist( File1, "sCor_HadEn", CutTier+"_"+CornName[cor], ";Hadronic En. for corner "+CornName[cor]+" (GeV); Number of events" );
502  TempHists.hCor_HFrEn[cor] = GetSpectToHist( File1, "sCor_HFrEn", CutTier+"_"+CornName[cor], ";Hadronic En. fraction for corner "+CornName[cor]+" (GeV); Number of events");
503  }
504  */
505  // --- Finally, push back TempHists into my vector.
506  MyHists.push_back(TempHists);
507  }
508  return MyHists;
509 }
510 
511 //......................................................
512 TH1D* GetSpectToHist( TFile* MyF, std::string HNam, std::string CTier, std::string Axis ) {
513  std::string LoadName = HNam + CTier;
514  std::unique_ptr<Spectrum> TempSpec = Spectrum::LoadFrom(MyF, LoadName.c_str() ) ;
515  if (TempSpec->POT() == 0) {
516  TH1D* MyHist = TempSpec -> ToTH1( LivNom, kLivetime );
517  CenterTitles ( MyHist);
518  MyHist->SetMinimum( 0 );
519  MyHist->SetMaximum( MyHist->GetMaximum()*1.05 );
520  MyHist->SetTitle ( Axis.c_str() );
521  return MyHist;
522  } else {
523  TH1D* MyHist = TempSpec -> ToTH1( POTNom );
524  CenterTitles ( MyHist);
525  MyHist->SetMinimum( 0 );
526  MyHist->SetMaximum( MyHist->GetMaximum()*1.05 );
527  MyHist->SetTitle ( Axis.c_str() );
528  return MyHist;
529  }
530 }
531 
532 //......................................................
533 TH1D* MakeRatio( TH1D* num, TH1D* denom, int Col, std::string FType ) {
534  // Declare a new histogram.
535  TH1D* ratHist = new TH1D( TString(FType)+TString("_ratio")+TString(std::to_string(Col))+TString(num->GetName()),
536  TString(";")+TString(num->GetXaxis()->GetTitle())+TString(";Ratio"),
537  num->GetNbinsX(),
538  num->GetXaxis()->GetXmin(),
539  num->GetXaxis()->GetXmax()
540  );
541  // Don't want any exponents?
542  ratHist->GetYaxis()->SetNoExponent();
543  ratHist->GetXaxis()->SetNoExponent();
544  // Divide the two input histograms...
545  ratHist->Divide( num , denom );
546  // Set Line and Marker colors / styles.
547  ratHist->SetLineColor ( Col );
548  ratHist->SetMarkerColor( Col );
549  ratHist->SetMarkerStyle( 20 );
550  ratHist->SetMarkerSize ( 0 );
551  // Figure out what the range should be...
552  ratHist->GetYaxis()->SetRangeUser( 0, 2 );
553  // Return my shiny new ratio plot!
554  return ratHist;
555 }
556 
557 //......................................................
558 void SetRange( TH1D* H1, std::vector<TH1D*> Shift, std::vector<TH1D*> Ratio, // Pass the histograms
559  double XLow, double XHigh, // What X values do I want?
560  double H_YLow, double H_YHigh, double R_YLow, double R_YHigh ) { // What Y values do I want?
561  // --- First off, set the X ranges.
562  if ( (XLow != -1) && (XHigh != -1) ) {
563  // Hists
564  H1 -> GetXaxis() -> SetRangeUser( XLow, XHigh );
565  for (unsigned int sh=0; sh<NumSysts; ++sh) {
566  Shift[sh] -> GetXaxis() -> SetRangeUser( XLow, XHigh );
567  Ratio[sh] -> GetXaxis() -> SetRangeUser( XLow, XHigh );
568  }
569  }
570 
571  // --- And, the Y axes, but check that they aren't -1!
572  // Hists
573  if ( (H_YLow != -1) && (H_YHigh != -1) ) {
574  H1 -> GetYaxis() -> SetRangeUser( H_YLow, H_YHigh );
575  for (unsigned int sh=0; sh<NumSysts; ++sh) {
576  Shift[sh] -> GetYaxis() -> SetRangeUser( H_YLow, H_YHigh );
577  }
578  }
579  // Ratios.
580  if ( (R_YLow != -1) && (R_YHigh != -1) ) {
581  for (unsigned int sh=0; sh<NumSysts; ++sh) {
582  Ratio[sh] -> GetYaxis() -> SetRangeUser( R_YLow, R_YHigh );
583  }
584  }
585  // --- Return
586  return;
587 }
588 
589 //......................................................
590 void FindAxisRange_DataMC( double &XLow, double &XHigh, double &H_YLow, double &H_YHigh, double &R_YLow, double &R_YHigh, double MaxVal, bool &SetLogy, std::string Name ) {
591 
592  XLow = -1;
593  XHigh = -1;
594  H_YLow = 0;
595  H_YHigh = MaxVal;
596  R_YLow = 0.5;
597  R_YHigh = 1.5;
598  SetLogy = false;
599 
600  if ( (Name.find("ReconstEnergy") != std::string::npos) ) { H_YHigh = 25; }
601  if ( (Name.find("HadFracEnergy") != std::string::npos) ) { H_YHigh = 20; }
602 
603  if ( (Name.find("HitsPerSlice") != std::string::npos) ) { XLow = 0.; XHigh = 300; }
604  if ( (Name.find("NumKalTracks") != std::string::npos) ) { XLow = 0.; XHigh = 5; }
605  if ( (Name.find("NumKalTrHits") != std::string::npos) ) { XLow = 0.; XHigh = 300; }
606  if ( (Name.find("MuonEnPerHit") != std::string::npos) ) { XLow = 0.01; XHigh = 0.03; H_YHigh = 50; }
607  if ( (Name.find("HadroEnPerHit") != std::string::npos) ) { XLow = 0.; XHigh = 0.04; H_YHigh = 50; }
608 
609  if ( (Name.find("CVNNuMuIDScor") != std::string::npos) ) { XLow = 0.5; XHigh = 1 ; H_YHigh = 80; }
610  if ( (Name.find("CVNCosmicScor") != std::string::npos) ) { XLow = 0 ; XHigh = 0.5; H_YHigh = 150;}
611  if ( (Name.find("RemIDScore" ) != std::string::npos) ) { XLow = 0.5; XHigh = 1 ; H_YHigh = 110;}
612  if ( (Name.find("NuMuContPID" ) != std::string::npos) ) { XLow = 0.4; XHigh = 0.9; H_YHigh = 60; }
613 
614  if ( (Name.find("Activity_Y" ) != std::string::npos) ) { H_YLow = 0.1; }
615  if ( (Name.find("Activity_Z" ) != std::string::npos) ) { H_YLow = 0.1; }
616  if ( (Name.find("KalMost" ) != std::string::npos) ) { H_YLow = 0.1; }
617 
618  if ( (Name.find("TrkStart" ) != std::string::npos) ) { H_YLow = 0; H_YHigh = 20; }
619  if ( (Name.find("TrkStartZ") != std::string::npos) ) { H_YLow = 0; H_YHigh = 23; }
620  if ( (Name.find("TrkEnd" ) != std::string::npos) ) { H_YLow = 0; H_YHigh = 20; }
621  /*
622  std::cout << "\t\tParams are; Xlow " << XLow << ", XHigh: " << XHigh << ", H_YLow: " << H_YLow << ", H_YHigh: " << H_YHigh
623  << ", R_YLow: " << R_YLow << ", R_YHigh " << R_YHigh << ", SetLogy " << SetLogy
624  << std::endl;
625  */
626  return;
627 }
628 
629 //......................................................
630 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 ) {
631  XLow = -1;
632  XHigh = -1;
633  H_YLow = 1e-2;
634  H_YHigh = MaxVal;
635  R_YLow = 0.5;
636  R_YHigh = 1.5;
637  SetLogy = false;
638 
639  if ( (Name.find("CVNNuMuIDScor") != std::string::npos) ) { XLow = 0.5; XHigh = 1 ; H_YLow = 0.1; SetLogy = true; }
640  if ( (Name.find("CVNCosmicScor") != std::string::npos) ) { XLow = 0 ; XHigh = 0.5; H_YLow = 0.1; SetLogy = true; }
641  if ( (Name.find("RemIDScore" ) != std::string::npos) ) { XLow = 0.5; XHigh = 1 ; H_YLow = 0.1; SetLogy = true; R_YLow = 0.8; R_YHigh = 1.2; }
642  if ( (Name.find("NuMuContPID" ) != std::string::npos) ) { XLow = 0.5; XHigh = 1 ; H_YLow = 0.1; SetLogy = true; }
643 
644  if ( (Name.find("Activity_Y" ) != std::string::npos) ) { H_YLow = 0.1; SetLogy = true; }
645  if ( (Name.find("Activity_Z" ) != std::string::npos) ) { H_YLow = 0.1; SetLogy = true; }
646  if ( (Name.find("KalMost" ) != std::string::npos) ) { H_YLow = 0.1; SetLogy = true; }
647 
648  if ( Name.find("VisEIn") != std::string::npos ) { XLow = -1; XHigh = -1; R_YLow = 0.8; R_YHigh = 1.2; }
649  if ( Name.find("NHitIn") != std::string::npos ) { XLow = 0; XHigh = 300; R_YLow = 0.8; R_YHigh = 1.2; }
650  if ( Name.find("VisEPerHit") != std::string::npos ) { XLow = 0.008; XHigh = 0.014; R_YLow = 0.8; R_YHigh = 1.2; }
651 
652  if ( Name.find("MuonEnPerHit" ) != std::string::npos ) { XLow = 0.01 ; XHigh = 0.03; R_YLow = 0.8; R_YHigh = 1.5; }
653  if ( Name.find("HadroEnPerHit") != std::string::npos ) { XLow = 0.001; XHigh = 0.03; R_YLow = 0.8; R_YHigh = 1.2; }
654  if ( Name.find("MuonEnergy" ) != std::string::npos ) { XLow = 0.01 ; XHigh = 3.0 ; R_YLow = 0.8; R_YHigh = 1.2; }
655 
656  if ( Name.find("ReTrOverTrNuE") != std::string::npos ) { XLow = -0.2 ; XHigh = 0.2 ; R_YLow = 0.7; R_YHigh = 1.3; }
657  if ( Name.find("HitsPerSlice") != std::string::npos ) { XLow = 0 ; XHigh = 200 ; R_YLow = 0.8; R_YHigh = 1.2; }
658  if ( Name.find("NumKalTrHits") != std::string::npos ) { XLow = 0 ; XHigh = 200 ; R_YLow = 0.8; R_YHigh = 1.2; }
659  /*
660  std::cout << "\t\tParams are; Xlow " << XLow << ", XHigh: " << XHigh << ", H_YLow: " << H_YLow << ", H_YHigh: " << H_YHigh
661  << ", R_YLow: " << R_YLow << ", R_YHigh " << R_YHigh << ", SetLogy " << SetLogy
662  << std::endl;
663  */
664  return;
665 }
666 //......................................................
667 void PlotHistProp( TH1D* Hist, TLegend* L ) {
668  double Mean_Nom = Hist->GetMean();
669  double RMS_Nom = Hist->GetRMS();
670  double Int_Nom = Hist->Integral();
671  L -> AddEntry( Hist, Form("Mean: %.3f", Mean_Nom), "" );
672  L -> AddEntry( Hist, Form("RMS: %.3f" , RMS_Nom ), "" );
673  L -> AddEntry( Hist, Form("Int: %.3f" , Int_Nom ), "" );
674  return;
675 }
676 //......................................................
677 // Put a "NOvA Preliminary" tag in the corner
678 void Preliminary() {
679  TLatex* prelim = new TLatex(.91, .95, "NOvA Preliminary");
680  prelim->SetTextColor(kBlue);
681  prelim->SetNDC();
682  prelim->SetTextSize(2/30.);
683  prelim->SetTextAlign(32);
684  prelim->Draw();
685 }
686 //..................................
688  // Does this canvas have a ratio on it?
689  std::string WRat = ", with data/mc ratio";
690  if (Name.find("_hist") != std::string::npos) WRat = ".";
691  // Remove the underscores from my cut tier
692  while (Cut.find("_") != std::string::npos)
693  Cut.replace(Cut.find("_"),1," "); //CTier = "Cosmic_Rej";
694  // Determine my caption
695  std::string Cap = "Plot showing the distribution of " + PTitle + " for Cut level and Quantile:" + Cut + WRat;
696  std::cout << "\nFile: " << Name << "\n\t Has caption: " << Cap << ".\n" << std::endl;
697  // Make the output text file
698  std::ofstream TxtOut ( Name.c_str(), std::ofstream::out );
699  TxtOut << Cap;
700  TxtOut.close();
701  // return
702  return;
703 }
704 //......................................................
705 void MakeSplitCans( TCanvas* Can, std::string DirNa, std::string CanTi, std::string FType, std::string CTier, TFile* OutF, TH1D* Nom, std::vector<TH1D*> Shi) {
706  // If DataMC, add cosmic to nominal.
707  if (MySys == "DataMC") {
708  Nom -> Add( Shi[1] );
709  }
710 
711  // --- What do I want to rebin?
712  double ReFac = 0;
713  if ( MySys == "DataMC" ) {
714  ReFac = 4;
715  if ( (DirNa.find("HitsPerSlice") != std::string::npos) || // 0 to 300
716  (DirNa.find("MuonEnPerHit") != std::string::npos) || // 0 to 0.04
717  (DirNa.find("NumKalTracks") != std::string::npos) // 0 to 5
718  ) {
719  ReFac = 0;
720  }
721 
722  if ( (DirNa.find("TrkStart") != std::string::npos) || // 0 to 300
723  (DirNa.find("TrkEnd" ) != std::string::npos) // 0 to 300
724  ) {
725  ReFac = 5;
726  }
727  }
728  if (ReFac != 0) {
729  Nom -> Rebin( ReFac );
730  for (unsigned int sh=0; sh<NumSysts; ++sh) {
731  Shi[sh] -> Rebin( ReFac );
732  }
733  }
734 
735  // --- Set the Colours of my histograms, and center titles.
736  // First Nominal
737  Nom -> SetLineColor( Colours[0] );
738  CenterTitles( Nom );
739  double MaxHVal = Nom->GetMaximum();
740  // Then the shifts.
741  std::vector<TH1D*> Rat;
742  for (unsigned int sh=0; sh<NumSysts; ++sh) {
743  if ( MaxHVal < Shi[sh]->GetMaximum() ) MaxHVal = Shi[sh]->GetMaximum(); // Find the Max value.
744  Shi[sh] -> SetLineColor( Colours[sh+1] ); // Set the colours properly.
745  Rat.push_back( MakeRatio( Shi[sh], Nom, Colours[sh+1], FType+CTier ) ); // Make the ratio.
746  CenterTitles( Rat.back() );
747  }
748 
749  // --- What axis ranges do I want to use? Find the ranges, and then set them.
750  double XL, XH, H_YL, H_YH, R_YL, R_YH;
751  bool SetLogy = false;
752  if (MySys == "DataMC") {
753  FindAxisRange_DataMC( XL, XH, H_YL, H_YH, R_YL, R_YH, 1.5*MaxHVal, SetLogy, DirNa );
754  } else {
755  FindAxisRange( XL, XH, H_YL, H_YH, R_YL, R_YH, 1.1*MaxHVal, SetLogy, DirNa );
756  }
757  SetRange( Nom, Shi, Rat, XL, XH, H_YL, H_YH, R_YL, R_YH);
758 
759  // For DataMC I want to make a graph of the 0th shift. If not MCData then I'll just ignore this...
760  TGraphAsymmErrors* grData = GraphWithPoissonErrors(Shi[0], false, false);
761  grData -> SetMarkerStyle(kFullCircle); grData -> SetLineWidth(2);
762  Shi[0] -> SetMarkerStyle(kFullCircle); Shi[0] -> SetLineWidth(2);
763 
764  // --- Now, let's make a histogram with just the plot.
765  std::string HCanNa = FType + CTier + DirNa + "_hist";
766  std::string HCanTi = HCanNa + "_hist";
767  TCanvas* HCan = new TCanvas( HCanNa.c_str(), HCanTi.c_str() );
768  // Where do I want my legend...
769  double HistLeg_X1 = 0.57, HistLeg_Y1 = 0.05, HistLeg_X2 = 0.92, HistLeg_Y2 = 0.88;
770  // If not plotting mean vals, then want smaller legend.
771  if (!PlotMean) HistLeg_Y1 = 0.65;
772  // Do I want the legend on the left instead?
773  if ( (DirNa.find("RatKalHitSlc") != std::string::npos) ||
774  (DirNa.find("RemIDScore") != std::string::npos) ||
775  (DirNa.find("_cCVNNuMuID") != std::string::npos) ||
776  (DirNa.find("_cNuMuContPID")!= std::string::npos)
777  ) {
778  HistLeg_X1 = 0.1;
779  HistLeg_X2 = 0.55;
780  }
781  // Make the legend where I've decided...
782  TLegend *HistLeg = new TLegend( HistLeg_X1, HistLeg_Y1, HistLeg_X2, HistLeg_Y2 );
783  HistLeg -> SetFillStyle (0);
784  HistLeg -> SetBorderSize(0);
785  HistLeg -> SetTextSize( Nom->GetXaxis()->GetTitleSize() );
786  // If DataMC want my Data as first thing on the legendto the legend
787  if (MySys == "DataMC" ) HistLeg -> AddEntry( grData, LegNames[1].c_str(), "lep" );
788  // Draw Nominal
789  Nom -> DrawCopy("hist");
790  HistLeg -> AddEntry( Nom, LegNames[0].c_str(), "l" );
791  if (PlotMean) PlotHistProp( Nom, HistLeg );
792  // Now draw the shifts!
793  for (unsigned int sh=0; sh<NumSysts; ++sh) {
794  if (MySys == "DataMC" && sh ==0 ) continue;
795  else {
796  Shi[sh] -> DrawCopy("hist same");
797  HistLeg -> AddEntry( Shi[sh], LegNames[sh+1].c_str(), "l" );
798  }
799  if (PlotMean) PlotHistProp( Shi[sh], HistLeg );
800  }
801  // If DataMC, then want to draw data last, and add POT
802  if (MySys == "DataMC") {
803  grData -> Draw("ep same");
804  HistLeg -> AddEntry( grData, Form("%.2f#times10^{20} POT", POTNom/1e20), "" );
805  }
806  // Do I want the canvas to be loagarithmic?
807  if ( SetLogy ) HCan -> SetLogy();
808  // And, draw the legend!
809  HistLeg -> Draw();
810  Preliminary();
811 
812  // --- Get ready to make the canvas with ratios.
813  std::string CanNa = FType + CTier + DirNa;
814  Can = new TCanvas( CanNa.c_str(), CanTi.c_str() );
815  // Sort out the canvas.
816  Can -> SetBottomMargin( 0. );
817  double Spl = 0.4;
818  TPad* P1 = new TPad( "Temp_1", "", 0.0, Spl, 1.0, 1.0, 0 );
819  TPad* P2 = new TPad( "Temp_2", "", 0.0, 0.0, 1.0, Spl, 0 );
820  P2 -> SetRightMargin (.03);
821  P2 -> SetTopMargin (.00);
822  P2 -> SetBottomMargin(.22);
823  P2 -> SetLeftMargin (.08);
824  P2 -> Draw();
825  P1 -> SetRightMargin (.03);
826  P1 -> SetLeftMargin (.08);
827  P1 -> SetTopMargin (.10);
828  P1 -> SetBottomMargin(.03);
829  P1 -> Draw();
830  // Set some label sizes.
831  double Lb1 = 0.07;
832  double Lb2 = 0.10;
833  // --- First, draw the ratios so cd onto Pad2
834  P2 -> cd();
835  // Set axis ranges etc.
836  Rat[0] -> GetYaxis() -> SetTitleSize( Lb2 );
837  Rat[0] -> GetYaxis() -> SetTitleOffset(0.4);
838  Rat[0] -> GetYaxis() -> SetLabelSize( Lb2 );
839  Rat[0] -> GetXaxis() -> SetTitleSize( Lb2 );
840  Rat[0] -> GetXaxis() -> SetLabelSize( Lb2 );
841  Rat[0] -> Draw();
842  // Add a reference line
843  //TLine *line = new TLine( Nom->GetXaxis()->GetXmin(), 1, Nom->GetXaxis()->GetXmax(), 1 );
844  double x[2] = { Nom->GetXaxis()->GetXmin(), Nom->GetXaxis()->GetXmax() };
845  double y[2] = { 1, 1 };
846  TGraph *line = new TGraph( 2, x, y );
847  line->SetLineColor(kGray);
848  line->SetLineWidth( 3 );
849  line->SetLineStyle( 2 );
850  line->Draw("l same");
851  // Redraw Rat[0]
852  Rat[0] -> Draw("axis same");
853  if (MySys != "DataMC") {
854  for (unsigned int sh=1; sh<NumSysts; ++sh) {
855  Rat[sh] -> Draw("same");
856  }
857  }
858  // --- And cd into Pad1
859  P1 -> cd();
860  // Where do I want my legend...
861  double Leg_X1 = 0.62, Leg_Y1 = 0.05, Leg_X2 = 0.94, Leg_Y2 = 0.85;
862  // If not plotting mean vals, then want smaller legend.
863  if (!PlotMean) Leg_Y1 = 0.4;
864  // Do I want the legend on the left instead?
865  if ( (DirNa.find("RatKalHitSlc") != std::string::npos) ||
866  (DirNa.find("RemIDScore") != std::string::npos) ||
867  (DirNa.find("_cCVNNuMuID") != std::string::npos) ||
868  (DirNa.find("_cNuMuContPID")!= std::string::npos)
869  ) {
870  Leg_X1 = 0.1;
871  Leg_X2 = 0.32;
872  }
873 
874  // --- Special Legend positions for DataMC.
875  // For Start/EndPos plots want to split the legend.
876 
877  // Make the legend where I've decided...
878  TLegend *Leg = new TLegend( Leg_X1, Leg_Y1, Leg_X2, Leg_Y2 );
879  Leg -> SetFillStyle (0);
880  Leg -> SetBorderSize(0);
881  Leg -> SetTextSize( 1.5*HistLeg -> GetTextSize() );
882  // --- Draw nominal first!
883  Nom -> GetYaxis() -> SetTitleSize( Lb1 );
884  Nom -> GetYaxis() -> SetLabelSize( Lb1 );
885  Nom -> GetYaxis() -> SetTitleOffset( 0.5 );
886  // Remove the x axis labels
887  Nom -> GetXaxis() -> SetLabelSize (0 );
888  Nom -> GetXaxis() -> SetLabelOffset(99);
889  Nom -> DrawCopy("hist");
890  // --- If DataMC want my Data as first thing on the legendto the legend
891  if (MySys == "DataMC" ) Leg -> AddEntry( grData, LegNames[1].c_str(), "lep" );
892  // --- Add nominal to legend
893  Leg -> AddEntry( Nom, LegNames[0].c_str(), "l" );
894  // --- Now draw the shifts!
895  for (unsigned int sh=0; sh<NumSysts; ++sh) {
896  if (MySys == "DataMC" && sh ==0) continue;
897  else {
898  Shi[sh] -> DrawCopy("hist same");
899  Leg -> AddEntry( Shi[sh], LegNames[sh+1].c_str(), "l" );
900  }
901  if (PlotMean) PlotHistProp( Shi[sh], Leg );
902  }
903  // Want to draw data last
904  if (MySys == "DataMC") {
905  grData -> Draw("ep same");
906  Leg -> AddEntry( grData, Form("%.2f#times10^{20} POT", POTNom/1e20), "" );
907  }
908  // Do I want Pad 1 to be loagarithmic?
909  if ( SetLogy ) P1 -> SetLogy();
910  // And, draw the legend!
911  Leg -> Draw();
912  Preliminary();
913 
914  // I need to make sure that my directory structure exists...
915  std::string SystDir = "Plots/"+ MySys +"/";
916  std::string CutsDir = SystDir + CTier +"/";
917  std::string CansDir = CutsDir + DirNa +"/";
918  gSystem -> MakeDirectory( SystDir.c_str() );
919  gSystem -> MakeDirectory( CutsDir.c_str() );
920  gSystem -> MakeDirectory( CansDir.c_str() );
921  // Now I can save my things...
922  std::vector<std::string> Exten = { ".png", ".pdf", ".txt" };
923  for (size_t ex=0; ex<Exten.size(); ++ex) {
924  // Determine the name of my file.
925  std::string HCanSaveLoc = CansDir + HCan->GetName() + Exten[ex];
926  std::string CanSaveLoc = CansDir + Can->GetName() + Exten[ex];
927  // Do something special if I have a .txt file
928  if ( Exten[ex] == ".txt" ) {
929  MakeTextFile( CTier, Nom->GetXaxis()->GetTitle(), HCanSaveLoc );
930  MakeTextFile( CTier, Nom->GetXaxis()->GetTitle(), CanSaveLoc );
931  } else {
932  // Save the canvases if not.
933  HCan -> SaveAs( HCanSaveLoc.c_str() );
934  Can -> SaveAs( CanSaveLoc .c_str() );
935  }
936  // And write the canvases to file too.
937  OutF -> cd();
938  HCan -> Write( HCan->GetName() );
939  Can -> Write( Can ->GetName() );
940  }
941  // --- Return.
942  return;
943 }
944 
945 //.....................................................
946 void ProducePlots( LoadedHistograms NomHists, std::vector<LoadedHistograms> SubSysts, std::string FType, std::string CTier ) {
947 
948  // Make the name a bit easier to read.
949  if (CTier.find("NeuCCuts") != std::string::npos)
950  CTier.replace(CTier.find("NeuCCuts"),8,"CosmicRej"); //CTier = "Cosmic_Rej";
951 
952  std::string OName = "Plots/Merged_"+FType+CTier+".root";
953  std::cout << "I want to make a file called " << OName << std::endl;
954  TFile* OutFile = TFile::Open( OName.c_str(), "RECREATE" );
955 
956  LoadedCanvases TheseCans;
957  std::vector<TH1D*> MyHists;
958  // *** Some key ones first.
959  // The reconstructed energy
960  MyHists.clear(); for (unsigned int sh=0; sh<NumSysts; ++sh) { MyHists.push_back( SubSysts[sh].hReconstEnergy ); }
961  MakeSplitCans( TheseCans.cReconstEnergy, "_cReconstEnergy", ("Total reconstructed energy for "+FType+" "+CTier), FType, CTier, OutFile, NomHists.hReconstEnergy , MyHists );
962 
963  // --- Do I only want to test things out plotting wise?
964  if (PlotTest) return;
965 
966  // Mean time of slice Full
967  MyHists.clear(); for (unsigned int sh=0; sh<NumSysts; ++sh) { MyHists.push_back( SubSysts[sh].hSliceTimeFull ); }
968  MakeSplitCans( TheseCans.cSliceTimeFull, "_cSliceTimeFull", ("Mean time of slice full window "+FType+" "+CTier), FType, CTier, OutFile, NomHists.hSliceTimeFull , MyHists );
969  // Mean time of slice Zoom
970  MyHists.clear(); for (unsigned int sh=0; sh<NumSysts; ++sh) { MyHists.push_back( SubSysts[sh].hSliceTimeZoom ); }
971  MakeSplitCans( TheseCans.cSliceTimeZoom, "_cSliceTimeZoom", ("Mean time of slice beam window "+FType+" "+CTier), FType, CTier, OutFile, NomHists.hSliceTimeZoom , MyHists );
972  // The reconstructed muon energy
973  MyHists.clear(); for (unsigned int sh=0; sh<NumSysts; ++sh) { MyHists.push_back( SubSysts[sh].hMuonEnergy ); }
974  MakeSplitCans( TheseCans.cMuonEnergy , "_cMuonEnergy" , ("Muon energy for "+FType+" "+CTier) , FType, CTier, OutFile, NomHists.hMuonEnergy , MyHists );
975  // The reconstructed muon energy per hit
976  MyHists.clear(); for (unsigned int sh=0; sh<NumSysts; ++sh) { MyHists.push_back( SubSysts[sh].hMuonEnPerHit ); }
977  MakeSplitCans( TheseCans.cMuonEnPerHit , "_cMuonEnPerHit" , ("Muon energy per hit for "+FType+" "+CTier) , FType, CTier, OutFile, NomHists.hMuonEnPerHit , MyHists );
978  // The reconstructed Hadronic energy
979  MyHists.clear(); for (unsigned int sh=0; sh<NumSysts; ++sh) { MyHists.push_back( SubSysts[sh].hHadronEnergy ); }
980  MakeSplitCans( TheseCans.cHadronEnergy , "_cHadronEnergy" , ("Hadronic energy for "+FType+" "+CTier) , FType, CTier, OutFile, NomHists.hHadronEnergy , MyHists );
981  // The reconstructed hadronic energy per hit
982  MyHists.clear(); for (unsigned int sh=0; sh<NumSysts; ++sh) { MyHists.push_back( SubSysts[sh].hHadroEnPerHit ); }
983  MakeSplitCans( TheseCans.cHadroEnPerHit, "_cHadroEnPerHit", ("Hadronic energy per hit for "+FType+" "+CTier) , FType, CTier, OutFile, NomHists.hHadroEnPerHit , MyHists );
984  // The Hadronic energy fraction
985  MyHists.clear(); for (unsigned int sh=0; sh<NumSysts; ++sh) { MyHists.push_back( SubSysts[sh].hHadFracEnergy ); }
986  MakeSplitCans( TheseCans.cHadFracEnergy, "_cHadFracEnergy", ("Hadronic energy fraction for "+FType+" "+CTier) , FType, CTier, OutFile, NomHists.hHadFracEnergy , MyHists );
987  // Number of hits per slice
988  MyHists.clear(); for (unsigned int sh=0; sh<NumSysts; ++sh) { MyHists.push_back( SubSysts[sh].hHitsPerSlice ); }
989  MakeSplitCans( TheseCans.cHitsPerSlice , "_cHitsPerSlice" , ("Hits per slice for "+FType+" "+CTier) , FType, CTier, OutFile, NomHists.hHitsPerSlice , MyHists );
990  // The RemID score
991  MyHists.clear(); for (unsigned int sh=0; sh<NumSysts; ++sh) { MyHists.push_back( SubSysts[sh].hRemIDScore ); }
992  MakeSplitCans( TheseCans.cRemIDScore , "_cRemIDScore" , ("The RemID score "+FType+" "+CTier) , FType, CTier, OutFile, NomHists.hRemIDScore , MyHists );
993  // CVN Cosmic score
994  MyHists.clear(); for (unsigned int sh=0; sh<NumSysts; ++sh) { MyHists.push_back( SubSysts[sh].hCVNCosmicScor ); }
995  MakeSplitCans( TheseCans.cCVNCosmicScor, "_cCVNCosmicScor" , ("CVN Cosmic score "+FType+" "+CTier) , FType, CTier, OutFile, NomHists.hCVNCosmicScor , MyHists );
996  // CVN NuMu score
997  MyHists.clear(); for (unsigned int sh=0; sh<NumSysts; ++sh) { MyHists.push_back( SubSysts[sh].hCVNNuMuIDScor ); }
998  MakeSplitCans( TheseCans.cCVNNuMuIDScor, "_cCVNNuMuIDScor" , ("CVN NuMu score "+FType+" "+CTier) , FType, CTier, OutFile, NomHists.hCVNNuMuIDScor , MyHists );
999  // Third analysis BDT
1000  MyHists.clear(); for (unsigned int sh=0; sh<NumSysts; ++sh) { MyHists.push_back( SubSysts[sh].hNuMuContPID ); }
1001  MakeSplitCans( TheseCans.cNuMuContPID , "_cNuMuContPID" , ("Third analysis BDT "+FType+" "+CTier) , FType, CTier, OutFile, NomHists.hNuMuContPID , MyHists );
1002  // Second analysis BDT cut
1003  MyHists.clear(); for (unsigned int sh=0; sh<NumSysts; ++sh) { MyHists.push_back( SubSysts[sh].hSANuMuContPID ); }
1004  MakeSplitCans( TheseCans.cSANuMuContPID, "_cSANuMuContPID" , ("Second analysis BDT "+FType+" "+CTier) , FType, CTier, OutFile, NomHists.hSANuMuContPID, MyHists );
1005 
1006  // *** Some important Truth based ones.
1007  // The true neutrino energy
1008  MyHists.clear(); for (unsigned int sh=0; sh<NumSysts; ++sh) { MyHists.push_back( SubSysts[sh].hTrueNuEnergy ); }
1009  MakeSplitCans( TheseCans.cTrueNuEnergy , "_cTrueNuEnergy" , ("True neutrino energy for "+FType+" "+CTier) , FType, CTier, OutFile, NomHists.hTrueNuEnergy , MyHists );
1010  // The Reco - True / True neutrino energy
1011  MyHists.clear(); for (unsigned int sh=0; sh<NumSysts; ++sh) { MyHists.push_back( SubSysts[sh].hReTrOverTrNuE ); }
1012  MakeSplitCans( TheseCans.cReTrOverTrNuE, "_cReTrOverTrNuE", ("Re-Tr/Te neutrino energy for "+FType+" "+CTier) , FType, CTier, OutFile, NomHists.hReTrOverTrNuE , MyHists );
1013 
1014  // *** Basic positional plots
1015  // Start X
1016  MyHists.clear(); for (unsigned int sh=0; sh<NumSysts; ++sh) { MyHists.push_back( SubSysts[sh].hTrkStartX ); }
1017  MakeSplitCans( TheseCans.cTrkStartX, "_cTrkStartX" , ("The Start X Pos "+FType+" "+CTier), FType, CTier, OutFile, NomHists.hTrkStartX , MyHists );
1018  // Start Y
1019  MyHists.clear(); for (unsigned int sh=0; sh<NumSysts; ++sh) { MyHists.push_back( SubSysts[sh].hTrkStartY ); }
1020  MakeSplitCans( TheseCans.cTrkStartY, "_cTrkStartY" , ("The Start Y Pos "+FType+" "+CTier), FType, CTier, OutFile, NomHists.hTrkStartY , MyHists );
1021  // Start Z
1022  MyHists.clear(); for (unsigned int sh=0; sh<NumSysts; ++sh) { MyHists.push_back( SubSysts[sh].hTrkStartZ ); }
1023  MakeSplitCans( TheseCans.cTrkStartZ, "_cTrkStartZ" , ("The Start Z Pos "+FType+" "+CTier), FType, CTier, OutFile, NomHists.hTrkStartZ , MyHists );
1024  // End X
1025  MyHists.clear(); for (unsigned int sh=0; sh<NumSysts; ++sh) { MyHists.push_back( SubSysts[sh].hTrkEndX ); }
1026  MakeSplitCans( TheseCans.cTrkEndX , "_cTrkEndX" , ("The End X Pos "+FType+" "+CTier) , FType, CTier, OutFile, NomHists.hTrkEndX , MyHists );
1027  // End Y
1028  MyHists.clear(); for (unsigned int sh=0; sh<NumSysts; ++sh) { MyHists.push_back( SubSysts[sh].hTrkEndY ); }
1029  MakeSplitCans( TheseCans.cTrkEndY , "_cTrkEndY" , ("The End Y Pos "+FType+" "+CTier) , FType, CTier, OutFile, NomHists.hTrkEndY , MyHists );
1030  // End Z
1031  MyHists.clear(); for (unsigned int sh=0; sh<NumSysts; ++sh) { MyHists.push_back( SubSysts[sh].hTrkEndZ ); }
1032  MakeSplitCans( TheseCans.cTrkEndZ , "_cTrkEndZ" , ("The End Z Pos "+FType+" "+CTier) , FType, CTier, OutFile, NomHists.hTrkEndZ , MyHists );
1033  // Len XY
1034  MyHists.clear(); for (unsigned int sh=0; sh<NumSysts; ++sh) { MyHists.push_back( SubSysts[sh].hTrkLenXY); }
1035  MakeSplitCans( TheseCans.cTrkLenXY , "_cTrkLenXY" , ("Length in XY " +FType+" "+CTier) , FType, CTier, OutFile, NomHists.hTrkLenXY , MyHists );
1036 
1037 
1038  /*
1039  // *** Track based ones...
1040  // Number of Kal trk in slice
1041  MyHists.clear(); for (unsigned int sh=0; sh<NumSysts; ++sh) { MyHists.push_back( SubSysts[sh].hNumKalTracks ); }
1042  MakeSplitCans( TheseCans.cNumKalTracks , "_cNumKalTracks" , ("Number of kalman track hits "+FType+" "+CTier) , FType, CTier, OutFile, NomHists.hNumKalTracks , MyHists );
1043  // Number of hits in Kal track
1044  MyHists.clear(); for (unsigned int sh=0; sh<NumSysts; ++sh) { MyHists.push_back( SubSysts[sh].hNumKalTrHits ); }
1045  MakeSplitCans( TheseCans.cNumKalTrHits , "_cNumKalTrHits" , ("Number of hits in kalman track "+FType+" "+CTier), FType, CTier, OutFile, NomHists.hNumKalTrHits , MyHists );
1046  // Ratio of Kal hits in slice
1047  MyHists.clear(); for (unsigned int sh=0; sh<NumSysts; ++sh) { MyHists.push_back( SubSysts[sh].hRatKalHitSlc ); }
1048  MakeSplitCans( TheseCans.cRatKalHitSlc , "_cRatKalHitSlc" , ("Ratio of kalman tr hit in slc "+FType+" "+CTier), FType, CTier, OutFile, NomHists.hRatKalHitSlc , MyHists );
1049  // Length of Kalman track
1050  MyHists.clear(); for (unsigned int sh=0; sh<NumSysts; ++sh) { MyHists.push_back( SubSysts[sh].hKalTrLength ); }
1051  MakeSplitCans( TheseCans.cKalTrLength , "_cKalTrLength" , ("Length of kalman track "+FType+" "+CTier) , FType, CTier, OutFile, NomHists.hKalTrLength , MyHists );
1052  // Angle Kal track & beam
1053  MyHists.clear(); for (unsigned int sh=0; sh<NumSysts; ++sh) { MyHists.push_back( SubSysts[sh].hKalTrBeamAng ); }
1054  MakeSplitCans( TheseCans.cKalTrBeamAng , "_cKalTrBeamAng" , ("Angle between kal track & beam "+FType+" "+CTier), FType, CTier, OutFile, NomHists.hKalTrBeamAng , MyHists );
1055  // Most forward Kal tr hit
1056  MyHists.clear(); for (unsigned int sh=0; sh<NumSysts; ++sh) { MyHists.push_back( SubSysts[sh].hKalMostFwdCel ); }
1057  MakeSplitCans( TheseCans.cKalMostFwdCel, "_cKalMostFwdCel", ("Most forward kalman track hit "+FType+" "+CTier) , FType, CTier, OutFile, NomHists.hKalMostFwdCel , MyHists );
1058  // Most backward Kal tr hit
1059  MyHists.clear(); for (unsigned int sh=0; sh<NumSysts; ++sh) { MyHists.push_back( SubSysts[sh].hKalMostBakCel ); }
1060  MakeSplitCans( TheseCans.cKalMostBakCel, "_cKalMostBakCel", ("Most backward kalman track hit "+FType+" "+CTier), FType, CTier, OutFile, NomHists.hKalMostBakCel , MyHists );
1061  // Highest Y Kalman tr vertex
1062  MyHists.clear(); for (unsigned int sh=0; sh<NumSysts; ++sh) { MyHists.push_back( SubSysts[sh].hKalTrVer_MaxY ); }
1063  MakeSplitCans( TheseCans.cKalTrVer_MaxY, "_cKalTrVer_MaxY" , ("Highest kalman track vertex Y "+FType+" "+CTier), FType, CTier, OutFile, NomHists.hKalTrVer_MaxY , MyHists );
1064  // Highest Z Kalman tr vertex
1065  MyHists.clear(); for (unsigned int sh=0; sh<NumSysts; ++sh) { MyHists.push_back( SubSysts[sh].hKalTrVer_MaxZ ); }
1066  MakeSplitCans( TheseCans.cKalTrVer_MaxZ, "_cKalTrVer_MaxZ" , ("Highest kalman track vertex Z "+FType+" "+CTier), FType, CTier, OutFile, NomHists.hKalTrVer_MaxZ, MyHists );
1067  // Y direction of Kal track
1068  MyHists.clear(); for (unsigned int sh=0; sh<NumSysts; ++sh) { MyHists.push_back( SubSysts[sh].hKalTrDir_Y ); }
1069  MakeSplitCans( TheseCans.cKalTrDir_Y , "_cKalTrDir_Y" , ("Y direction of kalman track "+FType+" "+CTier) , FType, CTier, OutFile, NomHists.hKalTrDir_Y , MyHists );
1070  // Scattering over tr length
1071  MyHists.clear(); for (unsigned int sh=0; sh<NumSysts; ++sh) { MyHists.push_back( SubSysts[sh].hScattKalTrLen ); }
1072  MakeSplitCans( TheseCans.cScattKalTrLen, "_cScattKalTrLen" , ("Sum of scattering over track "+FType+" "+CTier) , FType, CTier, OutFile, NomHists.hScattKalTrLen , MyHists );
1073 
1074  // *** Hit based ones...
1075  // First cell that is hit
1076  MyHists.clear(); for (unsigned int sh=0; sh<NumSysts; ++sh) { MyHists.push_back( SubSysts[sh].hFirstHitCell ); }
1077  MakeSplitCans( TheseCans.cFirstHitCell , "_cFirstHitCell" , ("First hit in the cell "+FType+" "+CTier) , FType, CTier, OutFile, NomHists.hFirstHitCell , MyHists );
1078  // Last cell that is hit
1079  MyHists.clear(); for (unsigned int sh=0; sh<NumSysts; ++sh) { MyHists.push_back( SubSysts[sh].hLastHitCell ); }
1080  MakeSplitCans( TheseCans.cLastHitCell , "_cLastHitCell" , ("Last hit in the cell "+FType+" "+CTier) , FType, CTier, OutFile, NomHists.hLastHitCell , MyHists );
1081  // Maximum Y pos of activty
1082  MyHists.clear(); for (unsigned int sh=0; sh<NumSysts; ++sh) { MyHists.push_back( SubSysts[sh].hMaxActivity_Y ); }
1083  MakeSplitCans( TheseCans.cMaxActivity_Y, "_cMaxActivity_Y", ("Maximum Y position of activity "+FType+" "+CTier), FType, CTier, OutFile, NomHists.hMaxActivity_Y , MyHists );
1084  // Minimum Y pos of activty
1085  MyHists.clear(); for (unsigned int sh=0; sh<NumSysts; ++sh) { MyHists.push_back( SubSysts[sh].hMinActivity_Y ); }
1086  MakeSplitCans( TheseCans.cMinActivity_Y, "_cMinActivity_Y", ("Minimum Y position of activity "+FType+" "+CTier), FType, CTier, OutFile, NomHists.hMinActivity_Y , MyHists );
1087  // Maximum Z pos of activty
1088  MyHists.clear(); for (unsigned int sh=0; sh<NumSysts; ++sh) { MyHists.push_back( SubSysts[sh].hMaxActivity_Z ); }
1089  MakeSplitCans( TheseCans.cMaxActivity_Z, "_cMaxActivity_Z" , ("Maximum Z position of activity "+FType+" "+CTier), FType, CTier, OutFile, NomHists.hMaxActivity_Z, MyHists );
1090  // Minimum Z pos of activty
1091  MyHists.clear(); for (unsigned int sh=0; sh<NumSysts; ++sh) { MyHists.push_back( SubSysts[sh].hMinActivity_Z ); }
1092  MakeSplitCans( TheseCans.cMinActivity_Z, "_cMinActivity_Z" , ("Minimum Y position of activity "+FType+" "+CTier), FType, CTier, OutFile, NomHists.hMinActivity_Z, MyHists );
1093  // Minimum cells to edge
1094  MyHists.clear(); for (unsigned int sh=0; sh<NumSysts; ++sh) { MyHists.push_back( SubSysts[sh].hMinCellToEdge ); }
1095  MakeSplitCans( TheseCans.cMinCellToEdge, "_cMinCellToEdge" , ("Minmum number of cell to edge "+FType+" "+CTier), FType, CTier, OutFile, NomHists.hMinCellToEdge , MyHists );
1096 
1097  // *** Visible energy quanta.
1098  // VisE In Slice
1099  MyHists.clear(); for (unsigned int sh=0; sh<NumSysts; ++sh) { MyHists.push_back( SubSysts[sh].hVisEInSlice ); }
1100  MakeSplitCans( TheseCans.cVisEInSlice , "_cVisEInSlice" , ("Visible energy in slice "+FType+" "+CTier) , FType, CTier, OutFile, NomHists.hVisEInSlice , MyHists );
1101  // NHit In Slice
1102  MyHists.clear(); for (unsigned int sh=0; sh<NumSysts; ++sh) { MyHists.push_back( SubSysts[sh].hNHitInSlice ); }
1103  MakeSplitCans( TheseCans.cNHitInSlice , "_cNHitInSlice" , ("Hits in slice "+FType+" "+CTier) , FType, CTier, OutFile, NomHists.hNHitInSlice , MyHists );
1104  // VisE Per Hit in Slice
1105  MyHists.clear(); for (unsigned int sh=0; sh<NumSysts; ++sh) { MyHists.push_back( SubSysts[sh].hVisEPerHitSlc ); }
1106  MakeSplitCans( TheseCans.cVisEPerHitSlc, "_cVisEPerHitSlc" , ("Visible energy per hit in slice "+FType+" "+CTier), FType, CTier, OutFile, NomHists.hVisEPerHitSlc, MyHists );
1107  // Vis E In Track
1108  MyHists.clear(); for (unsigned int sh=0; sh<NumSysts; ++sh) { MyHists.push_back( SubSysts[sh].hVisEInTrack ); }
1109  MakeSplitCans( TheseCans.cVisEInTrack , "_cVisEInTrack" , ("Visible track energy "+FType+" "+CTier) , FType, CTier, OutFile, NomHists.hVisEInTrack , MyHists );
1110  // NHit In Track
1111  MyHists.clear(); for (unsigned int sh=0; sh<NumSysts; ++sh) { MyHists.push_back( SubSysts[sh].hNHitInTrack ); }
1112  MakeSplitCans( TheseCans.cNHitInTrack , "_cNHitInTrack" , ("Number of track hits "+FType+" "+CTier) , FType, CTier, OutFile, NomHists.hNHitInTrack , MyHists );
1113  // VisE Per hit In Track
1114  MyHists.clear(); for (unsigned int sh=0; sh<NumSysts; ++sh) { MyHists.push_back( SubSysts[sh].hVisEPerHitTrk ); }
1115  MakeSplitCans( TheseCans.cVisEPerHitTrk, "_cVisEPerHitTrk" , ("Visible energy per track hit "+FType+" "+CTier) , FType, CTier, OutFile, NomHists.hVisEPerHitTrk, MyHists );
1116  // Hadronic VisE
1117  MyHists.clear(); for (unsigned int sh=0; sh<NumSysts; ++sh) { MyHists.push_back( SubSysts[sh].hVisEInHadron ); }
1118  MakeSplitCans( TheseCans.cVisEInHadron , "_cVisEInHadron" , ("Visible hadronic energy "+FType+" "+CTier) , FType, CTier, OutFile, NomHists.hVisEInHadron , MyHists );
1119  // Hadronic NHit
1120  MyHists.clear(); for (unsigned int sh=0; sh<NumSysts; ++sh) { MyHists.push_back( SubSysts[sh].hNHitInHadron ); }
1121  MakeSplitCans( TheseCans.cNHitInHadron , "_cNHitInHadron" , ("Number of hadronic hits "+FType+" "+CTier) , FType, CTier, OutFile, NomHists.hNHitInHadron , MyHists );
1122  // Hadronic VisE Per Hit
1123  MyHists.clear(); for (unsigned int sh=0; sh<NumSysts; ++sh) { MyHists.push_back( SubSysts[sh].hVisEPerHitHad ); }
1124  MakeSplitCans( TheseCans.cVisEPerHitHad, "_cVisEPerHitHad" , ("Visible hadronic energy per hit "+FType+" "+CTier), FType, CTier, OutFile, NomHists.hVisEPerHitHad, MyHists );
1125 
1126  // *** RemID score / inputs
1127  // The RemID scattering over track length
1128  MyHists.clear(); for (unsigned int sh=0; sh<NumSysts; ++sh) { MyHists.push_back( SubSysts[sh].hRemIDScatLLH ); }
1129  MakeSplitCans( TheseCans.cRemIDScatLLH , "_cRemIDScatLLH" , ("The RemID scattering over tr l "+FType+" "+CTier), FType, CTier, OutFile, NomHists.hRemIDScatLLH , MyHists );
1130  // The RemID dE/dx log-likelihood
1131  MyHists.clear(); for (unsigned int sh=0; sh<NumSysts; ++sh) { MyHists.push_back( SubSysts[sh].hRemIDdEdxLLH ); }
1132  MakeSplitCans( TheseCans.cRemIDdEdxLLH , "_cRemIDdEdxLLH" , ("The RemID dE/dx "+FType+" "+CTier) , FType, CTier, OutFile, NomHists.hRemIDdEdxLLH , MyHists );
1133  // The RemID fraction of planes which were measured.
1134  MyHists.clear(); for (unsigned int sh=0; sh<NumSysts; ++sh) { MyHists.push_back( SubSysts[sh].hRemIDMeasFrac ); }
1135  MakeSplitCans( TheseCans.cRemIDMeasFrac , "_cRemIDMeasFrac", ("The RemID frac of planes hit "+FType+" "+CTier) , FType, CTier, OutFile, NomHists.hRemIDMeasFrac , MyHists );
1136 
1137  // *** Some Calibration plots
1138  for (int cor=0; cor<4; ++cor) {
1139  // Reco Energy
1140  MyHists.clear(); for (unsigned int sh=0; sh<NumSysts; ++sh) { MyHists.push_back( SubSysts[sh].hCor_RecEn[cor] ); }
1141  MakeSplitCans( TheseCans.cCor_RecEn[cor], "_cCor_RecEn"+CornName[cor], ("RecoE_"+FType+"_"+CTier+"_"+CornName[cor]), FType, CTier, OutFile, NomHists.hCor_RecEn[cor] , MyHists );
1142  // Muon Energy
1143  MyHists.clear(); for (unsigned int sh=0; sh<NumSysts; ++sh) { MyHists.push_back( SubSysts[sh].hCor_MuoEn[cor] ); }
1144  MakeSplitCans( TheseCans.cCor_MuoEn[cor], "_cCor_MuoEn"+CornName[cor], ("MuonE_"+FType+"_"+CTier+"_"+CornName[cor]), FType, CTier, OutFile, NomHists.hCor_MuoEn[cor] , MyHists );
1145  // Hadronic Energy
1146  MyHists.clear(); for (unsigned int sh=0; sh<NumSysts; ++sh) { MyHists.push_back( SubSysts[sh].hCor_HadEn[cor] ); }
1147  MakeSplitCans( TheseCans.cCor_HadEn[cor], "_cCor_HadEn"+CornName[cor], ("HadE_" +FType+"_"+CTier+"_"+CornName[cor]), FType, CTier, OutFile, NomHists.hCor_HadEn[cor] , MyHists );
1148  // Hadronic Energy Fraction
1149  MyHists.clear(); for (unsigned int sh=0; sh<NumSysts; ++sh) { MyHists.push_back( SubSysts[sh].hCor_HFrEn[cor] ); }
1150  MakeSplitCans( TheseCans.cCor_HFrEn[cor], "_cCor_HFrEn"+CornName[cor], ("HadEF_"+FType+"_"+CTier+"_"+CornName[cor]), FType, CTier, OutFile, NomHists.hCor_HFrEn[cor] , MyHists );
1151  }
1152  */
1153  return;
1154 }
1155 
1156 //......................................................
1157 
1158 //......................................................
void MakeTextFile(std::string Cut, std::string PTitle, std::string Name)
Definition: DoThePlots.C:687
std::string LegNames[4]
Definition: DoThePlots.C:194
void ProducePlots(LoadedHistograms NomHists, std::vector< LoadedHistograms > SubSysts, std::string FType, std::string CTier)
Definition: DoThePlots.C:946
TH1D * hHadroEnPerHit
Definition: DoThePlots.C:41
tree Draw("slc.nhit")
void SetRange(TH1D *H1, std::vector< TH1D * > Shift, std::vector< TH1D * > Ratio, double XLow, double XHigh, double H_YLow, double H_YHigh, double R_YLow, double R_YHigh)
Definition: DoThePlots.C:558
TCanvas * cTrueNuEnergy
Definition: DoThePlots.C:130
const bool PlotTest
Definition: DoThePlots.C:221
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
#define H1(NAME, LABEL, N, X1, X2)
TCanvas * cTrkStartZ
Definition: DoThePlots.C:136
correl_xv GetYaxis() -> SetDecimals()
TH1D * hMuonEnergy
Definition: DoThePlots.C:38
TH1D * hReconstEnergy
Definition: DoThePlots.C:35
TCanvas * cTrkStartX
Definition: DoThePlots.C:134
void PlotHistProp(TH1D *Hist, TLegend *L)
Definition: DoThePlots.C:667
const int Colours[4]
Definition: DoThePlots.C:193
TCanvas * cSliceTimeFull
Definition: DoThePlots.C:116
prelim SetTextSize(2/30.)
TCanvas * cTrkStartY
Definition: DoThePlots.C:135
TH1D * hHadronEnergy
Definition: DoThePlots.C:40
Definition: Hist.h:29
std::string MySys
Definition: DoThePlots.C:217
const bool UseLightNominal
Definition: DoThePlots.C:223
TH1D * hTrkStartX
Definition: DoThePlots.C:54
unsigned int NumSysts
Definition: DoThePlots.C:215
TH2 * Rebin(TH2 *hist, std::vector< double > edges)
OStream cerr
Definition: OStream.cxx:7
const double kAna2017POT
Definition: Exposures.h:174
void MakeSplitCans(TCanvas *Can, std::string DirNa, std::string CanTi, std::string FType, std::string CTier, TFile *OutF, TH1D *Nom, std::vector< TH1D * > Shi)
Definition: DoThePlots.C:705
TCanvas * cTrkEndY
Definition: DoThePlots.C:138
std::vector< LoadedHistograms > LoadFile_GetHists(std::vector< std::string > FileLocs, std::string CTier)
Definition: DoThePlots.C:417
TCanvas * cMuonEnPerHit
Definition: DoThePlots.C:119
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
TH1D * hReTrOverTrNuE
Definition: DoThePlots.C:51
TH1D * hTrkStartZ
Definition: DoThePlots.C:56
fVtxDx SetMarkerStyle(20)
ntuple SetFillStyle(1001)
TCanvas * cHadFracEnergy
Definition: DoThePlots.C:122
void FindAxisRange_DataMC(double &XLow, double &XHigh, double &H_YLow, double &H_YHigh, double &R_YLow, double &R_YHigh, double MaxVal, bool &SetLogy, std::string Name)
Definition: DoThePlots.C:590
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: DoThePlots.C:630
TH1D * hTrueNuEnergy
Definition: DoThePlots.C:50
TH1D * MakeRatio(TH1D *num, TH1D *denom, int Col, std::string FType)
Definition: DoThePlots.C:533
Definition: Shift.h:6
std::string GetName(int i)
const double kAna2017Livetime
Definition: Exposures.h:200
legend SetBorderSize(0)
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:535
TCanvas * cReTrOverTrNuE
Definition: DoThePlots.C:131
static constexpr double L
TCanvas * cMuonEnergy
Definition: DoThePlots.C:118
correl_xv GetXaxis() -> SetDecimals()
hmean SetLineWidth(2)
leg AddEntry(GRdata,"data","p")
hmean SetLineColor(4)
TCanvas * cTrkEndZ
Definition: DoThePlots.C:139
TH1D * hTrkLenXY
Definition: DoThePlots.C:60
const bool PlotMean
Definition: DoThePlots.C:220
const int NumQuant
Definition: DoThePlots.C:206
TCanvas * cHitsPerSlice
Definition: DoThePlots.C:123
TH1D * hMuonEnPerHit
Definition: DoThePlots.C:39
const std::string QuantNames[NumQuant]
Definition: DoThePlots.C:207
TH1D * hHitsPerSlice
Definition: DoThePlots.C:43
TLatex * prelim
Definition: Xsec_final.C:133
TCanvas * cSANuMuContPID
Definition: DoThePlots.C:128
TH1D * hCVNCosmicScor
Definition: DoThePlots.C:45
Represent the ratio between two spectra.
Definition: Ratio.h:8
const bool SingleQuant
Definition: DoThePlots.C:219
TCanvas * cTrkLenXY
Definition: DoThePlots.C:140
OStream cout
Definition: OStream.cxx:6
TCanvas * cTrkEndX
Definition: DoThePlots.C:137
TCanvas * cNuMuContPID
Definition: DoThePlots.C:127
const std::string CornName[4]
Definition: DoThePlots.C:226
const int NumCut
Definition: DoThePlots.C:198
const bool SingleCut
Definition: DoThePlots.C:218
const Cut cut
Definition: exporter_fd.C:30
TH1D * hCVNNuMuIDScor
Definition: DoThePlots.C:46
TH1D * hSliceTimeZoom
Definition: DoThePlots.C:37
void DoThePlots()
Definition: DoThePlots.C:252
TH1D * hNuMuContPID
Definition: DoThePlots.C:47
TCanvas * cCVNCosmicScor
Definition: DoThePlots.C:125
TCanvas * cReconstEnergy
Definition: DoThePlots.C:115
int num
Definition: f2_nu.C:119
TCanvas * cRemIDScore
Definition: DoThePlots.C:124
void Preliminary()
Definition: DoThePlots.C:678
TH1D * hSANuMuContPID
Definition: DoThePlots.C:48
TFile ff[ntarg]
Definition: Style.C:26
TFile * OutFile
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
TH1D * hHadFracEnergy
Definition: DoThePlots.C:42
TH1D * hRemIDScore
Definition: DoThePlots.C:44
cosmicTree SaveAs("cosmicTree.root")
TCanvas * cSliceTimeZoom
Definition: DoThePlots.C:117
Float_t e
Definition: plot.C:35
TCanvas * cCVNNuMuIDScor
Definition: DoThePlots.C:126
enum BeamMode kBlue
TH1D * GetSpectToHist(TFile *MyF, std::string HNam, std::string CTier, std::string Axis)
Definition: DoThePlots.C:512
c cd(1)
TCanvas * cHadronEnergy
Definition: DoThePlots.C:120
double LivNom
Definition: DoThePlots.C:196
TCanvas * cHadroEnPerHit
Definition: DoThePlots.C:121
TH1D * hSliceTimeFull
Definition: DoThePlots.C:36
gPad SetLogy()
const std::string CutNames[NumCut]
Definition: DoThePlots.C:199
gm Write()
TH1D * hTrkStartY
Definition: DoThePlots.C:55
double POTNom
Definition: DoThePlots.C:195
enum BeamMode string