MuonCatcherComp_ProdPlots.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 // --- Number of Time Bins and Quartiles.
39 const int NumTimeBins = 6;
40 const int NumQuartile = 5;
41 //......................................................
42 struct LoadedHistograms {
43  // *** Order the plots by [Quartiles][Time Bins]
44  TH1D* MuonRecoE [NumQuartile][NumTimeBins];
45  TH1D* MuonRecoE_MC [NumQuartile][NumTimeBins];
46 
47  TH1D* MuonNHit [NumQuartile][NumTimeBins];
48  TH1D* MuonNHit_MC [NumQuartile][NumTimeBins];
49 
50  TH1D* MuonCalE [NumQuartile][NumTimeBins];
51  TH1D* MuonCalE_MC [NumQuartile][NumTimeBins];
52 
53  TH1D* MuonTrkLen [NumQuartile][NumTimeBins];
54  TH1D* MuonTrkLen_MC[NumQuartile][NumTimeBins];
55 };
56 
57 //......................................................
58 // --- Some important things which need to be configured at run time.
60 bool IsFHC = true;
61 std::string sFHC = "FHC";
62 
63 // --- Some variables that are global, but change.
65 const bool ExtremeBins = true; // Only plot time bins 1 & 6?
66 const bool SingleQuant = false; // Only plot combined quartiles?
67 const bool PlotMean = true; // Plot the Mean, Int, RMS?
68 const bool PlotTest = false; // Only plot one variable?
69 const bool SaveRatio = false; // Produce and save the ratio plots?
70 const bool SaveHists = true; // Produce and save the non-ratio plots?
71 
72 //......................................................
74 //......................................................
75 TH1D* GetSpectToHist( TFile* MyF, double POT, std::string SpecName, std::string Axis );
76 //......................................................
77 void ProducePlots( std::vector<LoadedHistograms> HistVec );
78 //......................................................
79 TH1D* MakeRatio( TH1D* num, TH1D* denom, std::string RatName );
80 //......................................................
81 void SetRange( std::vector<TH1D*> hData, std::vector<TH1D*> hMont, std::vector<TH1D*> hRatio, // Pass the histograms
82  double XLow, double XHigh, double H_YLow, double H_YHigh, double R_YLow, double R_YHigh ); // What Y values do I want?
83 //......................................................
84 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 );
85 //......................................................
86 void MakeSplitCans( std::string PlotNa, std::string QntTxt, TFile* OutF,
87  std::vector<TH1D*> DataHists, std::vector<TH1D*> MontHists, std::vector<std::string> Legends );
88 //......................................................
89 void PlotHistProp( TH1D* Hist, TLegend* L );
90 //......................................................
91 void MakeTextFile( std::string Cut, std::string PTitle, std::string Name );
92 //......................................................
93 void CornerLabel(std::string Str);
94 
95 //......................................................
96 void MuonCatcherComp_ProdPlots( bool Horn ) {
97 
98  IsFHC = Horn;
99  if ( IsFHC) { sFHC = "FHC"; POTNom = kAna2018FHCPOT; }
100  if (!IsFHC) { sFHC = "RHC"; POTNom = kAna2018RHCPOT; }
101 
102  // ---- First off, lets set the styles...
103  gStyle->SetOptStat(0);
104  gROOT->SetStyle("novaStyle");
105 
106  // --- Declare a vector of my file name structures
107  std::string BaseInputDir="/nova/ana/users/mbaird42/prod4_neutron_studies/neutron_shine_ND_studies/booster_batch_slicer_differences/hadded_muon_results_v2/";
108  std::vector<std::string> InputFiles;
109  if ( IsFHC ) {
110  InputFiles.push_back( BaseInputDir+"hadd_neutron_shine_muons_v2_FHC_data.root" );
111  InputFiles.push_back( BaseInputDir+"hadd_neutron_shine_muons_v2_FHC_MC.root" );
112  } else {
113  InputFiles.push_back( BaseInputDir+"hadd_neutron_shine_muons_v2_RHC_data.root" );
114  InputFiles.push_back( BaseInputDir+"hadd_neutron_shine_muons_v2_RHC_MC.root" );
115  }
116 
117  // --- What POT am I normalising to?
118  TFile *TempDataFile = TFile::Open( InputFiles[0].c_str() );
119  std::unique_ptr<Spectrum> TempDataSpec = Spectrum::LoadFrom(TempDataFile, "dir_MuonE_Q1_TB1" ) ;
120  POTNom = TempDataSpec->POT();
121  std::cout << "\n sFHC is " << sFHC << ", and my POT is " << POTNom << "\n" << std::endl;
122 
123  // --- Load all of the histograms that I want.
124  std::vector<LoadedHistograms> RawHistograms;
125  for (size_t InFile=0; InFile<InputFiles.size(); ++InFile) {
126  std::cout << "\n Loading the histograms for file " << InputFiles[InFile] << ".\n" << std::endl;
127  RawHistograms.push_back( LoadFile_GetHists( InputFiles[InFile], POTNom ) );
128  }
129 
130  // --- And now I can make the plots!!!
131  ProducePlots( RawHistograms );
132  // --- Finally, return.
133  return;
134 } // DoThePlots
135 
136 //......................................................
138  // --- Declare the vector of LoadedHistograms which I am going to return...
139  LoadedHistograms MyHists;
140 
141  // --- Load the current file.
142  TFile *InFile = TFile::Open( InputFile.c_str() );
143  // --- Grab my spectra...
144  for ( int Qnt=0; Qnt<NumQuartile; ++Qnt ) {
145  for ( int TB=0; TB<NumTimeBins; ++TB ) {
146  std::string QntTB = "_Q"+std::to_string(Qnt+1)+"_TB"+std::to_string(TB+1);
147 
148  MyHists.MuonRecoE [Qnt][TB] = GetSpectToHist( InFile, POT, "dir_MuonE" +QntTB, ";Reconstructed muon energy (GeV);" );
149  MyHists.MuonRecoE_MC [Qnt][TB] = GetSpectToHist( InFile, POT, "dir_MuonE_MC" +QntTB, ";Reconstructed muon energy (GeV), Muon catcher;" );
150 
151  MyHists.MuonNHit [Qnt][TB] = GetSpectToHist( InFile, POT, "dir_Muon_NHit" +QntTB, ";Number of muon hits;" );
152  MyHists.MuonNHit_MC [Qnt][TB] = GetSpectToHist( InFile, POT, "dir_Muon_NHit_MC" +QntTB, ";Number of muon hits, Muon Catcher;" );
153 
154  MyHists.MuonCalE [Qnt][TB] = GetSpectToHist( InFile, POT, "dir_Muon_calE" +QntTB, ";Calorimetric muon energy (GeV);" );
155  MyHists.MuonCalE_MC [Qnt][TB] = GetSpectToHist( InFile, POT, "dir_Muon_calE_MC" +QntTB, ";Calorimetric muon energy (GeV), Muon catcher;" );
156 
157  MyHists.MuonTrkLen [Qnt][TB] = GetSpectToHist( InFile, POT, "dir_Muon_trk_len" +QntTB, ";Muon Track Length (m);" );
158  MyHists.MuonTrkLen_MC[Qnt][TB] = GetSpectToHist( InFile, POT, "dir_Muon_trk_len_MC"+QntTB, ";Muon Track Length (m), Muon catcher;" );
159  }
160  }
161  return MyHists;
162 }
163 
164 //......................................................
165 TH1D* GetSpectToHist( TFile* MyF, double POT, std::string SpecName, std::string Axis ) {
166 
167  //std::cout << "Loading Spectrum from " << SpecName << std::endl;
168 
169  // --- Load the Spectrum, and convert to a TH1
170  std::unique_ptr<Spectrum> TempSpec = Spectrum::LoadFrom( MyF, SpecName.c_str() ) ;
171  TH1D* MyHist = TempSpec -> ToTH1( POTNom );
172 
173  // --- Make the initial pass of making my histogram pretty...
174  CenterTitles ( MyHist);
175  MyHist->SetMinimum( 0 );
176  MyHist->SetMaximum( MyHist->GetMaximum()*1.05 );
177  MyHist->SetTitle ( Axis.c_str() );
178 
179  // --- Rebin everything...
180  double ReFac = 0;
181  if ( (SpecName.find("MuonE" ) != std::string::npos) ||
182  (SpecName.find("NHit" ) != std::string::npos) ||
183  (SpecName.find("calE" ) != std::string::npos) ||
184  (SpecName.find("trk_len") != std::string::npos)
185  ) {
186  ReFac = 10;
187  }
188  if ( ReFac!=0 )
189  MyHist -> Rebin( ReFac );
190 
191  // --- I want to area normalise everything...
192  MyHist -> Scale( 1/MyHist->Integral() );
193 
194  // -- Return my histogram.
195  return MyHist;
196 }
197 
198 //......................................................
199 TH1D* MakeRatio( TH1D* num, TH1D* denom, std::string RatName ) {
200  // Declare a new histogram.
201  TH1D* ratHist = new TH1D( TString(RatName)+TString("_ratio"),
202  TString(RatName)+TString("_ratio")+TString(";")+TString(num->GetXaxis()->GetTitle())+TString(";Ratio"),
203  num->GetNbinsX(),
204  num->GetXaxis()->GetXmin(),
205  num->GetXaxis()->GetXmax()
206  );
207  // Don't want any exponents?
208  ratHist->GetYaxis()->SetNoExponent();
209  ratHist->GetXaxis()->SetNoExponent();
210  // Divide the two input histograms...
211  ratHist->Divide( num , denom );
212  ratHist->SetMarkerStyle( 20 );
213  ratHist->SetMarkerSize ( 0 );
214  // Figure out what the range should be...
215  ratHist->GetYaxis()->SetRangeUser( 0, 2 );
216  // Return my shiny new ratio plot!
217  return ratHist;
218 }
219 
220 //......................................................
221 void SetRange( std::vector<TH1D*> hData, std::vector<TH1D*> hMont, std::vector<TH1D*> hRatio, // Pass the histograms
222  double XLow, double XHigh, double H_YLow, double H_YHigh, double R_YLow, double R_YHigh ){ // What Y values do I want?
223 
224  for (size_t hh=0; hh<hData.size(); ++hh) {
225  // --- First off, set the X ranges.
226  if ( (XLow != -1) && (XHigh != -1) ) {
227  // Hists
228  hData [hh] -> GetXaxis() -> SetRangeUser( XLow, XHigh );
229  hMont [hh] -> GetXaxis() -> SetRangeUser( XLow, XHigh );
230  hRatio[hh] -> GetXaxis() -> SetRangeUser( XLow, XHigh );
231  }
232 
233  // --- And, the Y axes, but check that they aren't -1!
234  // Hists
235  if ( (H_YLow != -1) && (H_YHigh != -1) ) {
236  hData[hh] -> GetYaxis() -> SetRangeUser( H_YLow, H_YHigh );
237  hMont[hh] -> GetYaxis() -> SetRangeUser( H_YLow, H_YHigh );
238  }
239  // Ratios.
240  if ( (R_YLow != -1) && (R_YHigh != -1) ) {
241  hRatio[hh] -> GetYaxis() -> SetRangeUser( R_YLow, R_YHigh );
242  }
243  }
244  // --- Return
245  return;
246 }
247 
248 //......................................................
249 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 ) {
250  XLow = -1;
251  XHigh = -1;
252  H_YLow = 0;
253  H_YHigh = MaxVal;
254  R_YLow = 0.5;
255  R_YHigh = 1.5;
256  SetLogy = false;
257 
258  //if ( (Name.find("MuonE" ) != std::string::npos) ) { XLow = 0; XHigh = 150; }
259  if ( (Name.find("NHit" ) != std::string::npos) ) { XLow = 0; XHigh = 300; }
260  //if ( (Name.find("CalE" ) != std::string::npos) ) { XLow = 0; XHigh = 150; }
261  //if ( (Name.find("TrkLen") != std::string::npos) ) { XLow = 0; XHigh = 150; }
262 
263  return;
264 }
265 //......................................................
266 void PlotHistProp( TH1D* Hist, TLegend* L ) {
267  double Mean_Nom = Hist->GetMean();
268  double Mean_Err = Hist->GetMeanError();
269  double RMS_Nom = Hist->GetRMS();
270  double Int_Nom = Hist->Integral();
271  L -> AddEntry( Hist, Form("Mean: %.3f #pm %.3f", Mean_Nom, Mean_Err), "" );
272  //L -> AddEntry( Hist, Form("RMS: %.3f" , RMS_Nom ), "" );
273  //L -> AddEntry( Hist, Form("Int: %.3f" , Int_Nom ), "" );
274  return;
275 }
276 //......................................................
278  // Does this canvas have a ratio on it?
279  std::string WRat = ", with data/mc ratio";
280  if (Name.find("_hist") != std::string::npos) WRat = ".";
281  // Remove the underscores from my cut tier
282  while (Cut.find("_") != std::string::npos)
283  Cut.replace(Cut.find("_"),1," "); //CTier = "Cosmic_Rej";
284  // Determine my caption
285  std::string Cap = "Plot showing the distribution of " + PTitle + " for Cut level and Quartile:" + Cut + WRat;
286  if (MySys == "DataMC" && PTitle.find("neutrino energy") == std::string::npos) {
287  Cap += " The simulated entries are simply taken from MC CAF files, as such they are not individually extrapolated / weighted from RecoE.";
288  }
289  std::cout << "\nFile: " << Name << "\n\t Has caption: " << Cap << ".\n" << std::endl;
290  // Make the output text file
291  std::ofstream TxtOut ( Name.c_str(), std::ofstream::out );
292  TxtOut << Cap;
293  TxtOut.close();
294  // return
295  return;
296 }
297 //......................................................
299  TLatex* CornLab = new TLatex(.1, .93, Str.c_str());
300  CornLab->SetTextColor(kGray+1);
301  CornLab->SetNDC();
302  CornLab->SetTextSize (2/30.);
303  CornLab->SetTextAlign(11);
304  CornLab->Draw();
305 }
306 //......................................................
307 void MakeSplitCans( std::string PlotNa, std::string QntTxt, TFile* OutF,
308  std::vector<TH1D*> DataHists, std::vector<TH1D*> MontHists, std::vector<std::string> Legends
309  ) {
310 
311  //std::cout << "Don't want to make Split Cans, returning." << std::endl;
312  //return;
313 
314  std::cout << "At the start of MakeSplitCans, " << PlotNa << ", " << QntTxt << ", DataHists.size() is " << DataHists.size()
315  << ", MontHists.size() is " << MontHists.size() << ", and Legends.size() is " << Legends.size()
316  << std::endl;
317 
318  // -- Make the ratio plots.
319  std::vector<TH1D*> Ratios;
320  for (size_t hh=0; hh<DataHists.size(); ++hh) {
321  Ratios.push_back( MakeRatio( DataHists[hh], MontHists[hh], (PlotNa+QntTxt+Legends[hh]) ) );
322  }
323 
324  // --- Set the Colours of my histograms, and center titles.
325  double MaxHVal = 0;
326  for (size_t hh=0; hh<DataHists.size(); ++hh) {
327  // Center titles.
328  CenterTitles( DataHists[hh] );
329  CenterTitles( MontHists[hh] );
330  CenterTitles( Ratios [hh] );
331  // Set the colour.
332  if (DataHists.size()==2 && hh==0) { // If only two plots, and looking at TB1
333  DataHists[hh] -> SetLineColor( 1 ); DataHists[hh] -> SetMarkerColor( 1 );
334  MontHists[hh] -> SetLineColor( 1 ); Ratios [hh] -> SetLineColor ( 1 );
335  } else if (DataHists.size()==2 && hh==1) { // If only two plots, and looking at TB6
336  DataHists[hh] -> SetLineColor( 2 ); DataHists[hh] -> SetMarkerColor( 2 );
337  MontHists[hh] -> SetLineColor( 2 ); Ratios [hh] -> SetLineColor ( 2 );
338  } else { // If plotting all time bins --> All the colours!
339  DataHists[hh] -> SetLineColor( hh ); DataHists[hh] -> SetMarkerColor( hh );
340  MontHists[hh] -> SetLineColor( hh ); Ratios [hh] -> SetLineColor( hh );
341  }
342  MontHists[hh] -> SetLineStyle( 2 );
343  // Figure out the maximum value.
344  if (DataHists[hh]->GetMaximum() > MaxHVal) MaxHVal = DataHists[hh]->GetMaximum();
345  if (MontHists[hh]->GetMaximum() > MaxHVal) MaxHVal = MontHists[hh]->GetMaximum();
346  }
347 
348  // --- What axis ranges do I want to use? Find the ranges, and then set them.
349  double XL, XH, H_YL, H_YH, R_YL, R_YH;
350  bool SetLogy = false;
351  FindAxisRange( XL, XH, H_YL, H_YH, R_YL, R_YH, 1.1*MaxHVal, SetLogy, PlotNa );
352  SetRange ( DataHists, MontHists, Ratios, XL, XH, H_YL, H_YH, R_YL, R_YH );
353 
354  //......................................................//......................................................
355  // ------------------------------ Get ready to make the canvas without ratios. ------------------------------
356  //......................................................//......................................................
357  std::string HCanNa = PlotNa+"_"+QntTxt+"_"+sFHC+"_hist";
358  TCanvas* HCan = new TCanvas( HCanNa.c_str(), HCanNa.c_str() );
359  // Where do I want my legend...
360  double HistLeg_X1 = 0.55, HistLeg_Y1 = 0.2, HistLeg_X2 = 0.85, HistLeg_Y2 = 0.88;
361  // If not plotting mean vals, then want smaller legend.
362  if (!PlotMean) {
363  HistLeg_Y1 = 0.55;
364  }
365  // Do I want the legend on the left instead?
366  if ( (PlotNa.find("MaxPos") != std::string::npos) // 0 to 300
367  ) {
368  HistLeg_X1 = 0.10;
369  HistLeg_X2 = 0.56;
370  }
371  // Make the legend where I've decided...
372  TLegend *HistLeg = new TLegend( HistLeg_X1, HistLeg_Y1, HistLeg_X2, HistLeg_Y2 );
373  HistLeg -> SetFillStyle (0);
374  HistLeg -> SetBorderSize(0);
375  HistLeg -> SetTextSize( 0.8*DataHists[0]->GetXaxis()->GetTitleSize() );
376 
377  // --- Now plot stuff!
378  for (size_t hh=0; hh<DataHists.size(); ++hh) {
379  // First data
380  if (hh==0) DataHists[hh] -> Draw("lep");
381  else DataHists[hh] -> Draw("lep same");
382  HistLeg -> AddEntry( DataHists[hh], TString("Data ")+TString(Legends[hh].c_str()), "lep" );
383  if (PlotMean) PlotHistProp( DataHists[hh], HistLeg );
384  // Then MC
385  MontHists[hh] -> Draw("hist same");
386  HistLeg -> AddEntry( MontHists[hh], TString("MC ") +TString(Legends[hh].c_str()), "l" );
387  if (PlotMean) PlotHistProp( MontHists[hh], HistLeg );
388  }
389  HistLeg -> AddEntry( DataHists[0], QntTxt.c_str(), "" );
390  DataHists[0] -> Draw( "axis same" );
391  // Do I want the canvas to be loagarithmic?
392  if ( SetLogy ) HCan -> SetLogy();
393  // Finally, draw the legend and whether it is FHC or RHC etc!
394  HistLeg -> Draw();
395  Preliminary();
396  if (IsFHC) CornerLabel( "Neutrino beam" );
397  else CornerLabel( "Antineutrino beam" );
398 
399  //......................................................//......................................................
400  // ------------------------------ Get ready to make the canvas with ratios. ------------------------------
401  //......................................................//......................................................
402  /*
403  std::string RatCanNa = FType + CTier + DirNa;
404  TCanvas* RatCan = new TCanvas( RatCanNa.c_str(), CanTi.c_str() );
405  // Sort out the canvas.
406  RatCan -> SetBottomMargin( 0. );
407  double Spl = 0.4;
408  TPad* P1 = new TPad( "Temp_1", "", 0.0, Spl, 1.0, 1.0, 0 );
409  TPad* P2 = new TPad( "Temp_2", "", 0.0, 0.0, 1.0, Spl, 0 );
410  P2 -> SetRightMargin (.03);
411  P2 -> SetTopMargin (.00);
412  P2 -> SetBottomMargin(.22);
413  P2 -> SetLeftMargin (.08);
414  P2 -> Draw();
415  P1 -> SetRightMargin (.03);
416  P1 -> SetLeftMargin (.08);
417  P1 -> SetTopMargin (.10);
418  P1 -> SetBottomMargin(.03);
419  P1 -> Draw();
420  // Set some label sizes.
421  double Lb1 = 0.07;
422  double Lb2 = 0.10;
423  // --- First, draw the ratios so cd onto Pad2
424  P2 -> cd();
425  // Set axis ranges etc.
426  Rat[0] -> GetYaxis() -> SetTitleSize( Lb2 );
427  Rat[0] -> GetYaxis() -> SetTitleOffset(0.4);
428  Rat[0] -> GetYaxis() -> SetLabelSize( Lb2 );
429  Rat[0] -> GetXaxis() -> SetTitleSize( Lb2 );
430  Rat[0] -> GetXaxis() -> SetLabelSize( Lb2 );
431  Rat[0] -> Draw();
432  // Add a reference line
433  //TLine *line = new TLine( Nom->GetXaxis()->GetXmin(), 1, Nom->GetXaxis()->GetXmax(), 1 );
434  double x[2] = { Nom->GetXaxis()->GetXmin(), Nom->GetXaxis()->GetXmax() };
435  double y[2] = { 1, 1 };
436  TGraph *line = new TGraph( 2, x, y );
437  line->SetLineColor(kGray);
438  line->SetLineWidth( 3 );
439  line->SetLineStyle( 2 );
440  line->Draw("l same");
441  // Redraw Rat[0]
442  Rat[0] -> Draw("axis same");
443  if (MySys != "DataMC") {
444  for (unsigned int sh=1; sh<NumSysts; ++sh) {
445  Rat[sh] -> Draw("same");
446  }
447  }
448  // --- And cd into Pad1
449  P1 -> cd();
450  // Where do I want my legend...
451  double Leg_X1 = 0.62, Leg_Y1 = 0.05, Leg_X2 = 0.94, Leg_Y2 = 0.85;
452  // If not plotting mean vals, then want smaller legend.
453  if (!PlotMean) Leg_Y1 = 0.4;
454  // Do I want the legend on the left instead?
455  if ( (DirNa.find("RatKalHitSlc") != std::string::npos) ||
456  (DirNa.find("RemIDScore") != std::string::npos) ||
457  (DirNa.find("_cCVNNuMuID") != std::string::npos) ||
458  (DirNa.find("_cNuMuContPID")!= std::string::npos)
459  ) {
460  Leg_X1 = 0.1;
461  Leg_X2 = 0.32;
462  }
463 
464  // --- Special Legend positions for DataMC.
465  // For Start/EndPos plots want to split the legend.
466 
467  // Make the legend where I've decided...
468  TLegend *Leg = new TLegend( Leg_X1, Leg_Y1, Leg_X2, Leg_Y2 );
469  Leg -> SetFillStyle (0);
470  Leg -> SetBorderSize(0);
471  Leg -> SetTextSize( 1.25*HistLeg -> GetTextSize() );
472  // --- Draw nominal first!
473  Nom -> GetYaxis() -> SetTitleSize( Lb1 );
474  Nom -> GetYaxis() -> SetLabelSize( Lb1 );
475  Nom -> GetYaxis() -> SetTitleOffset( 0.5 );
476  // Remove the x axis labels
477  Nom -> GetXaxis() -> SetLabelSize (0 );
478  Nom -> GetXaxis() -> SetLabelOffset(99);
479  Nom -> DrawCopy("hist");
480  // --- If DataMC want my Data as first thing on the legendto the legend
481  if (MySys == "DataMC" ) Leg -> AddEntry( grData, LegNames[1].c_str(), "lep" );
482  // --- Add nominal to legend
483  Leg -> AddEntry( Nom, LegNames[0].c_str(), "l" );
484  if (PlotMean) PlotHistProp( Nom, Leg );
485  // --- Now draw the shifts!
486  for (unsigned int sh=0; sh<NumSysts; ++sh) {
487  if (MySys == "DataMC" ) {
488  if ( sh != 0 ) {
489  Shi[sh] -> DrawCopy("hist same");
490  Leg -> AddEntry( Shi[sh], LegNames[sh+1].c_str(), "f" );
491  }
492  } else {
493  Shi[sh] -> DrawCopy("hist same");
494  Leg -> AddEntry( Shi[sh], LegNames[sh+1].c_str(), "l" );
495  }
496  if (PlotMean) PlotHistProp( Shi[sh], Leg );
497  }
498  // Want to draw data last
499  if (MySys == "DataMC") {
500  grData -> Draw("ep same");
501  }
502  Nom -> Draw( "axis same" );
503  // Do I want Pad 1 to be loagarithmic?
504  if ( SetLogy ) P1 -> SetLogy();
505  // And, draw the legend!
506  Leg -> Draw();
507  Preliminary();
508  if (IsFHC) CornerLabel( "Neutrino beam" );
509  else CornerLabel( "Antineutrino beam" );
510  */
511 
512  //......................................................//......................................................
513  // ------------------------------ Now lets all both of those canvases. ------------------------------
514  //......................................................//......................................................
515  std::string PDir = "Plots_MuonCatcher/";
516  std::string HornDir = PDir + sFHC +"/";
517  std::string CansDir = HornDir + PlotNa +"/";
518  gSystem -> MakeDirectory( HornDir.c_str() );
519  gSystem -> MakeDirectory( CansDir.c_str() );
520  // Now I can save my things...
521  if (SaveHists) {
522  std::string HCanSaveLoc = CansDir + HCan->GetName();
523  HCan -> SaveAs( TString(HCanSaveLoc)+TString(".pdf") );
524  HCan -> SaveAs( TString(HCanSaveLoc)+TString(".png") );
525  // And write the canvases to file too.
526  OutF -> cd();
527  HCan -> Write( HCan->GetName() );
528  }
529  /*
530  if (SaveRatio) {
531  for (size_t ex=0; ex<Exten.size(); ++ex) {
532  // Determine the name of my file.
533  std::string RatCanSaveLoc = CansDir + RatCan->GetName() + Exten[ex];
534  // Do something special if I have a .txt file
535  if ( Exten[ex] == ".txt" ) {
536  MakeTextFile( CTier, Nom->GetXaxis()->GetTitle(), RatCanSaveLoc );
537  } else {
538  // Save the canvases if not.
539  RatCan -> SaveAs( RatCanSaveLoc .c_str() );
540  }
541  // And write the canvases to file too.
542  OutF -> cd();
543  RatCan -> Write( RatCan ->GetName() );
544  }
545  */
546  // --- Return.
547  return;
548 }
549 
550 //.....................................................
551 void ProducePlots( std::vector<LoadedHistograms> HistVec ) {
552 
553  // --- I want to make a directory to hold my plots.
554  std::string PDir = "Plots_MuonCatcher/";
555  gSystem -> MakeDirectory( PDir.c_str() );
556  // --- I want to make a file to contain all of my plots.
557  std::string OName = PDir+"Quantile_TimeBins_"+sFHC+"_Plots.root";
558  TFile* OutFile = TFile::Open( OName.c_str(), "RECREATE" );
559  std::cout << "I want to make a file called " << OName << std::endl;
560 
561  // --- I need to loop through all of my histograms and construct the vectors that I want to use to plot things...
562  for ( int Qnt=0; Qnt<NumQuartile; ++ Qnt ) {
563  std::string QuartTxt = "Quartile_"+std::to_string(Qnt+1);
564 
565  std::vector<TH1D* > MyDataHists;
566  std::vector<TH1D* > MyMontHists;
567  std::vector<std::string> MyLegends;
568 
569  // --- Muon Reconstructed Energy
570  MyDataHists.clear(); MyMontHists.clear(); MyLegends.clear();
571  for ( int TB=0; TB<NumTimeBins; ++TB ) {
572  // If plotting extreme bins, only want to fill if TB == 1 or 6
573  if ( ExtremeBins && !(TB+1==1 || TB+1==6) ) continue;
574  // Fill the vectors...
575  MyDataHists.push_back( HistVec[0].MuonRecoE[ Qnt ][ TB ] );
576  MyMontHists.push_back( HistVec[1].MuonRecoE[ Qnt ][ TB ] );
577  MyLegends .push_back( "Time Bin "+std::to_string(TB+1) );
578  }
579  MakeSplitCans( "MuonRecoE", QuartTxt, OutFile, MyDataHists, MyMontHists, MyLegends );
580 
581  // ++++ Want to make my profiles histograms...
582  std::string CanTi = "AvMuonRecoE_Quartile"+std::to_string(Qnt+1)+"_"+sFHC;
583  TCanvas* cMuonEn = new TCanvas( CanTi.c_str(), CanTi.c_str() );
584  TH1D* MuonEn_Data = new TH1D( TString(CanTi)+TString("_Data"), "Average Energy; TimeBin; Average Energy", 6, 0.5, 6.5 );
585  TH1D* MuonEn_Mont = new TH1D( TString(CanTi)+TString("_Mont"), "Average Energy; TimeBin; Average Energy", 6, 0.5, 6.5 );
586  double HistMax = 0, HistMin = 5, HistAv = 0;
587  MuonEn_Mont -> SetLineColor(2);
588  for ( int TB=0; TB<NumTimeBins; ++TB ) {
589  MuonEn_Data -> SetBinContent( TB+1, HistVec[0].MuonRecoE[ Qnt ][ TB ]->GetMean() );
590  MuonEn_Data -> SetBinError ( TB+1, HistVec[0].MuonRecoE[ Qnt ][ TB ]->GetMeanError() );
591  MuonEn_Mont -> SetBinContent( TB+1, HistVec[1].MuonRecoE[ Qnt ][ TB ]->GetMean() );
592  MuonEn_Mont -> SetBinError ( TB+1, HistVec[1].MuonRecoE[ Qnt ][ TB ]->GetMeanError() );
593  if (HistVec[0].MuonRecoE[ Qnt ][ TB ]->GetMean() > HistMax) HistMax = HistVec[0].MuonRecoE[ Qnt ][ TB ]->GetMean();
594  if (HistVec[0].MuonRecoE[ Qnt ][ TB ]->GetMean() < HistMin) HistMin = HistVec[0].MuonRecoE[ Qnt ][ TB ]->GetMean();
595  }
596  HistAv = (HistMin+HistMax)/2;
597  MuonEn_Data->GetYaxis()->SetRangeUser( HistAv-0.02, HistAv+0.02 );
598  MuonEn_Data->Draw();
599  MuonEn_Mont->Draw("same");
600  std::string SaveLoc = "Plots_MuonCatcher/"+sFHC+"/MuonRecoE/"+CanTi;
601  cMuonEn -> SaveAs( TString(SaveLoc)+TString(".pdf") );
602  cMuonEn -> SaveAs( TString(SaveLoc)+TString(".png") );
603 
604  // --- Muon Reconstructed Energy, Muon Catcher
605  MyDataHists.clear(); MyMontHists.clear(); MyLegends.clear();
606  for ( int TB=0; TB<NumTimeBins; ++TB ) {
607  // If plotting extreme bins, only want to fill if TB == 1 or 6
608  if ( ExtremeBins && !(TB+1==1 || TB+1==6) ) continue;
609  // Fill the vectors...
610  MyDataHists.push_back( HistVec[0].MuonRecoE_MC[ Qnt ][ TB ] );
611  MyMontHists.push_back( HistVec[1].MuonRecoE_MC[ Qnt ][ TB ] );
612  MyLegends .push_back( "Time Bin "+std::to_string(TB+1) );
613  }
614  MakeSplitCans( "MuonRecoE_MC", QuartTxt, OutFile, MyDataHists, MyMontHists, MyLegends );
615 
616  // ++++ Want to make my profiles histograms...
617  CanTi = "AvMuonRecoE_MC_Quartile"+std::to_string(Qnt+1)+"_"+sFHC;
618  TCanvas* cMuonEn_MC = new TCanvas( CanTi.c_str(), CanTi.c_str() );
619  TH1D* MuonEn_MC_Data = new TH1D( TString(CanTi)+TString("_Data"), "Average Energy; TimeBin; Average Energy", 6, 0.5, 6.5 );
620  TH1D* MontEn_MC_Mont = new TH1D( TString(CanTi)+TString("_Mont"), "Average Energy; TimeBin; Average Energy", 6, 0.5, 6.5 );
621  HistMax = 0; HistMin = 5; HistAv = 0;
622  MontEn_MC_Mont -> SetLineColor(2);
623  for ( int TB=0; TB<NumTimeBins; ++TB ) {
624  MuonEn_MC_Data -> SetBinContent( TB+1, HistVec[0].MuonRecoE_MC[ Qnt ][ TB ]->GetMean() );
625  MuonEn_MC_Data -> SetBinError ( TB+1, HistVec[0].MuonRecoE_MC[ Qnt ][ TB ]->GetMeanError() );
626  MontEn_MC_Mont -> SetBinContent( TB+1, HistVec[1].MuonRecoE_MC[ Qnt ][ TB ]->GetMean() );
627  MontEn_MC_Mont -> SetBinError ( TB+1, HistVec[1].MuonRecoE_MC[ Qnt ][ TB ]->GetMeanError() );
628  if (HistVec[0].MuonRecoE_MC[ Qnt ][ TB ]->GetMean() > HistMax) HistMax = HistVec[0].MuonRecoE_MC[ Qnt ][ TB ]->GetMean();
629  if (HistVec[0].MuonRecoE_MC[ Qnt ][ TB ]->GetMean() < HistMin) HistMin = HistVec[0].MuonRecoE_MC[ Qnt ][ TB ]->GetMean();
630  }
631  HistAv = (HistMin+HistMax)/2;
632  MuonEn_MC_Data -> GetYaxis()->SetRangeUser( HistAv-0.02, HistAv+0.02 );
633  MuonEn_MC_Data -> Draw();
634  MontEn_MC_Mont -> Draw("same");
635  SaveLoc = "Plots_MuonCatcher/"+sFHC+"/MuonRecoE_MC/"+CanTi;
636  cMuonEn_MC -> SaveAs( TString(SaveLoc)+TString(".pdf") );
637  cMuonEn_MC -> SaveAs( TString(SaveLoc)+TString(".png") );
638 
639  // --- Muon Number of Hits
640  MyDataHists.clear(); MyMontHists.clear(); MyLegends.clear();
641  for ( int TB=0; TB<NumTimeBins; ++TB ) {
642  // If plotting extreme bins, only want to fill if TB == 1 or 6
643  if ( ExtremeBins && !(TB+1==1 || TB+1==6) ) continue;
644  // Fill the vectors...
645  MyDataHists.push_back( HistVec[0].MuonNHit[ Qnt ][ TB ] );
646  MyMontHists.push_back( HistVec[1].MuonNHit[ Qnt ][ TB ] );
647  MyLegends .push_back( "Time Bin "+std::to_string(TB+1) );
648  }
649  MakeSplitCans( "MuonNHit", QuartTxt, OutFile, MyDataHists, MyMontHists, MyLegends );
650 
651  // --- Muon Number of Hits, Muon Catcher
652  MyDataHists.clear(); MyMontHists.clear(); MyLegends.clear();
653  for ( int TB=0; TB<NumTimeBins; ++TB ) {
654  // If plotting extreme bins, only want to fill if TB == 1 or 6
655  if ( ExtremeBins && !(TB+1==1 || TB+1==6) ) continue;
656  // Fill the vectors...
657  MyDataHists.push_back( HistVec[0].MuonNHit_MC[ Qnt ][ TB ] );
658  MyMontHists.push_back( HistVec[1].MuonNHit_MC[ Qnt ][ TB ] );
659  MyLegends .push_back( "Time Bin "+std::to_string(TB+1) );
660  }
661  MakeSplitCans( "MuonNHit_MC", QuartTxt, OutFile, MyDataHists, MyMontHists, MyLegends );
662 
663  // --- Muon Calorimetric Energy
664  MyDataHists.clear(); MyMontHists.clear(); MyLegends.clear();
665  for ( int TB=0; TB<NumTimeBins; ++TB ) {
666  // If plotting extreme bins, only want to fill if TB == 1 or 6
667  if ( ExtremeBins && !(TB+1==1 || TB+1==6) ) continue;
668  // Fill the vectors...
669  MyDataHists.push_back( HistVec[0].MuonCalE[ Qnt ][ TB ] );
670  MyMontHists.push_back( HistVec[1].MuonCalE[ Qnt ][ TB ] );
671  MyLegends .push_back( "Time Bin "+std::to_string(TB+1) );
672  }
673  MakeSplitCans( "MuonCalE", QuartTxt, OutFile, MyDataHists, MyMontHists, MyLegends );
674 
675  // --- Muon Calorimetric Energy, Muon Catcher
676  MyDataHists.clear(); MyMontHists.clear(); MyLegends.clear();
677  for ( int TB=0; TB<NumTimeBins; ++TB ) {
678  // If plotting extreme bins, only want to fill if TB == 1 or 6
679  if ( ExtremeBins && !(TB+1==1 || TB+1==6) ) continue;
680  // Fill the vectors...
681  MyDataHists.push_back( HistVec[0].MuonCalE_MC[ Qnt ][ TB ] );
682  MyMontHists.push_back( HistVec[1].MuonCalE_MC[ Qnt ][ TB ] );
683  MyLegends .push_back( "Time Bin "+std::to_string(TB+1) );
684  }
685  MakeSplitCans( "MuonCalE_MC", QuartTxt, OutFile, MyDataHists, MyMontHists, MyLegends );
686 
687  // --- Muon Track Length
688  MyDataHists.clear(); MyMontHists.clear(); MyLegends.clear();
689  for ( int TB=0; TB<NumTimeBins; ++TB ) {
690  // If plotting extreme bins, only want to fill if TB == 1 or 6
691  if ( ExtremeBins && !(TB+1==1 || TB+1==6) ) continue;
692  // Fill the vectors...
693  MyDataHists.push_back( HistVec[0].MuonTrkLen[ Qnt ][ TB ] );
694  MyMontHists.push_back( HistVec[1].MuonTrkLen[ Qnt ][ TB ] );
695  MyLegends .push_back( "Time Bin "+std::to_string(TB+1) );
696  }
697  MakeSplitCans( "MuonTrkLen", QuartTxt, OutFile, MyDataHists, MyMontHists, MyLegends );
698 
699  // --- Muon Track Length, Muon Catcher
700  MyDataHists.clear(); MyMontHists.clear(); MyLegends.clear();
701  for ( int TB=0; TB<NumTimeBins; ++TB ) {
702  // If plotting extreme bins, only want to fill if TB == 1 or 6
703  if ( ExtremeBins && !(TB+1==1 || TB+1==6) ) continue;
704  // Fill the vectors...
705  MyDataHists.push_back( HistVec[0].MuonTrkLen_MC[ Qnt ][ TB ] );
706  MyMontHists.push_back( HistVec[1].MuonTrkLen_MC[ Qnt ][ TB ] );
707  MyLegends .push_back( "Time Bin "+std::to_string(TB+1) );
708  }
709  MakeSplitCans( "MuonTrkLen_MC", QuartTxt, OutFile, MyDataHists, MyMontHists, MyLegends );
710 
711  }
712  return;
713 }
714 
715 //......................................................
716 
void SetRange(std::vector< TH1D * > hData, std::vector< TH1D * > hMont, std::vector< TH1D * > hRatio, double XLow, double XHigh, double H_YLow, double H_YHigh, double R_YLow, double R_YHigh)
tree Draw("slc.nhit")
void MakeSplitCans(std::string PlotNa, std::string QntTxt, TFile *OutF, std::vector< TH1D * > DataHists, std::vector< TH1D * > MontHists, std::vector< std::string > Legends)
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
void MakeTextFile(std::string Cut, std::string PTitle, std::string Name)
correl_xv GetYaxis() -> SetDecimals()
const bool SaveHists
TH1D * MuonNHit_MC[NumQuartile][NumTimeBins]
hmean SetLineStyle(2)
prelim SetTextSize(2/30.)
double POTNom
Definition: Hist.h:29
LoadedHistograms LoadFile_GetHists(std::string InputFile, double POT)
TH2 * Rebin(TH2 *hist, std::vector< double > edges)
TH1D * MuonRecoE[NumQuartile][NumTimeBins]
void CenterTitles(TH1 *histo)
Definition: Plots.cxx:1481
ntuple SetFillStyle(1001)
TH1D * MuonRecoE_MC[NumQuartile][NumTimeBins]
const double kAna2018RHCPOT
Definition: Exposures.h:208
const bool SaveRatio
const int NumQuartile
legend SetBorderSize(0)
std::string MySys
const bool SingleQuant
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:535
void CornerLabel(std::string Str)
const bool ExtremeBins
static constexpr double L
correl_xv GetXaxis() -> SetDecimals()
leg AddEntry(GRdata,"data","p")
hmean SetLineColor(4)
h1 SetMarkerColor(4)
void ProducePlots(std::vector< LoadedHistograms > HistVec)
TH1D * MakeRatio(TH1D *num, TH1D *denom, std::string RatName)
TH1D * MuonCalE[NumQuartile][NumTimeBins]
void Preliminary()
OStream cout
Definition: OStream.cxx:6
std::vector< double > POT
cout<< "--"<< endl;for(Int_t iP=1;iP<=hyz->GetNbinsX();iP++){for(Int_t iC=1;iC<=hyz->GetNbinsY();iC++){if(hyv->GetBinContent(iP, iC)>-999){goal_hyv-> SetBinContent(iP, iC,-(dy[iP-1][iC-1]))
TH1D * MuonTrkLen[NumQuartile][NumTimeBins]
TH1D * MuonNHit[NumQuartile][NumTimeBins]
const int NumTimeBins
std::string sFHC
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)
TH1D * MuonTrkLen_MC[NumQuartile][NumTimeBins]
const bool PlotTest
int num
Definition: f2_nu.C:119
void PlotHistProp(TH1D *Hist, TLegend *L)
TFile * OutFile
const double kAna2018FHCPOT
Definition: Exposures.h:207
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
TH1D * MuonCalE_MC[NumQuartile][NumTimeBins]
cosmicTree SaveAs("cosmicTree.root")
void MuonCatcherComp_ProdPlots(bool Horn)
simulatedPeak Scale(1/simulationNormalisationFactor)
c cd(1)
TH1D * GetSpectToHist(TFile *MyF, double POT, std::string SpecName, std::string Axis)
hh[ixs]
Definition: PlotSingle.C:6
const bool PlotMean
gPad SetLogy()
gm Write()
enum BeamMode string