neutronE_plot.C
Go to the documentation of this file.
1 #include "TCanvas.h"
2 #include "TFile.h"
3 #include "TH1D.h"
4 #include "TLatex.h"
5 #include "TLegend.h"
6 #include "TStyle.h"
7 #include "TSystem.h"
8 #include "TF1.h"
9 #include "TH2.h"
10 
12 #include "CAFAna/Analysis/Plots.h"
13 #include "CAFAna/Core/Spectrum.h"
15 
16 #include <iostream>
17 #include <sstream>
18 #include <fstream>
19 #include <vector>
20 #include <string>
21 #include <TROOT.h>
22 #include <TStyle.h>
23 
24 using namespace ana;
25 
26 // ****************************** Some very broad things ******************************
27 
28 // --- Define the number of GENIE interaction types I'm looking at....
29 const std::vector<std::string> GENIEStr = { "All" }; const unsigned int nGENIE = 1;
30 //const std::vector<std::string> GENIEStr = { "All", "MEC", "DIS", "Resonance", "Coherent", "QuasiElactic" }; const unsigned int nGENIE = 6;
31 
32 // --- Define my broad cuts
33 std::vector<std::string> FD_CutNames = {"FD_2018-PID", "FD_2018-Cont", "FD_2018!PID", "FD_2018!Cont", "FD_2018full"}; const unsigned int nFD_CutNames = 5;
34 std::vector<std::string> ND_CutNames = {"ND_2018-PID", "ND_2018-Cont", "ND_2018!PID", "ND_2018!Cont", "ND_2018full"}; const unsigned int nND_CutNames = 5;
35 unsigned int FCut_FD = nFD_CutNames - 1; unsigned int FCut_ND = nND_CutNames - 1;
36 
37 // --- Am I running quantile boundaries?
38 std::vector<std::string> QuantNames = { "AllQuant", "Quant1", "Quant2", "Quant3", "Quant4" }; const unsigned int nQuantNames = 5;
39 
40 // --- FHC or RHC?
41 std::vector<std::string> HCNames = {"fhc", "rhc"}; const unsigned int nHC = 2;
42 
43 // --- ND or FD? (for labeling)
44 std::vector<std::string> DetNames = {"ND", "FD"};
45 
46 // ****************************** Define a struct ******************************
47 struct MySpectra {
48  // *** Far detector histograms
49  TH1D* NeutrinoEnergy_FD [nGENIE][nFD_CutNames][nQuantNames][nHC];
50  TH1D* TotPrimEnergy_FD [nGENIE][nFD_CutNames][nQuantNames][nHC];
51  TH1D* TotPrimEnergynKEp20_FD[nGENIE][nFD_CutNames][nQuantNames][nHC];
52  TH1D* TotPrimEnergynKEm20_FD[nGENIE][nFD_CutNames][nQuantNames][nHC];
53  TH1D* NeutMult_FD [nGENIE][nFD_CutNames][nQuantNames][nHC];
54  TH1D* NeutEnergy_FD [nGENIE][nFD_CutNames][nQuantNames][nHC];
55  TH1D* NeutEFrac_FD [nGENIE][nFD_CutNames][nQuantNames][nHC];
56  TH1D* NeutPngCalE_FD [nGENIE][nFD_CutNames][nQuantNames][nHC];
57  TH1D* NeutPngCalEpH_FD [nGENIE][nFD_CutNames][nQuantNames][nHC];
58  TH1D* NeutPng2DCalE_FD [nGENIE][nFD_CutNames][nQuantNames][nHC];
59  TH1D* NeutPng2DCalEpH_FD [nGENIE][nFD_CutNames][nQuantNames][nHC];
60  TH1D* NeutPngWtCalE_FD [nGENIE][nFD_CutNames][nQuantNames][nHC];
61  TH1D* NeutPngWtCalEpH_FD [nGENIE][nFD_CutNames][nQuantNames][nHC];
62  TH1D* NeutPng2DWtCalE_FD [nGENIE][nFD_CutNames][nQuantNames][nHC];
63  TH1D* NeutPng2DWtCalEpH_FD [nGENIE][nFD_CutNames][nQuantNames][nHC];
64  TH2* NEFvsCVN_FD [nGENIE][nFD_CutNames][nQuantNames][nHC];
65  TH2* NEFvsReMId_FD [nGENIE][nFD_CutNames][nQuantNames][nHC];
66 
67 
68  // *** Near detector histograms
69  TH1D* NeutrinoEnergy_ND [nGENIE][nND_CutNames][nQuantNames][nHC];
70  TH1D* TotPrimEnergy_ND [nGENIE][nND_CutNames][nQuantNames][nHC];
71  TH1D* TotPrimEnergynKEp20_ND[nGENIE][nND_CutNames][nQuantNames][nHC];
72  TH1D* TotPrimEnergynKEm20_ND[nGENIE][nND_CutNames][nQuantNames][nHC];
73  TH1D* NeutMult_ND [nGENIE][nND_CutNames][nQuantNames][nHC];
74  TH1D* NeutEnergy_ND [nGENIE][nND_CutNames][nQuantNames][nHC];
75  TH1D* NeutEFrac_ND [nGENIE][nND_CutNames][nQuantNames][nHC];
76  TH1D* NeutPngCalE_ND [nGENIE][nND_CutNames][nQuantNames][nHC];
77  TH1D* NeutPngCalEpH_ND [nGENIE][nND_CutNames][nQuantNames][nHC];
78  TH1D* NeutPng2DCalE_ND [nGENIE][nND_CutNames][nQuantNames][nHC];
79  TH1D* NeutPng2DCalEpH_ND [nGENIE][nND_CutNames][nQuantNames][nHC];
80  TH1D* NeutPngWtCalE_ND [nGENIE][nND_CutNames][nQuantNames][nHC];
81  TH1D* NeutPngWtCalEpH_ND [nGENIE][nND_CutNames][nQuantNames][nHC];
82  TH1D* NeutPng2DWtCalE_ND [nGENIE][nND_CutNames][nQuantNames][nHC];
83  TH1D* NeutPng2DWtCalEpH_ND [nGENIE][nND_CutNames][nQuantNames][nHC];
84  TH2* NEFvsCVN_ND [nGENIE][nND_CutNames][nQuantNames][nHC];
85  TH2* NEFvsReMId_ND [nGENIE][nND_CutNames][nQuantNames][nHC];
86 
87 };
88 
89 // ****************************** Define functions ******************************
90 TH1D* SpecToHist( TFile* MyF, std::string SpecName, double POT, std::string AxisTitle );
91 TH2* Spec2DToHist( TFile* MyF, std::string SpecName, double POT, std::string AxisTitle );
92 void FitGauss(TH1D* h);
93 void MakeComparisonPlots( std::vector<TH1D*> TheHists, std::string PlotName, std::vector<std::string> Labels);
94 void Print2D(TH2* hist, std::string PlotName);
95 
96 // ****************************** Main code block ******************************
97 void neutronE_plot() {
98  // ---- First off, lets set the styles...
99  gStyle->SetOptStat(0);
100  gROOT->SetStyle("novaStyle");
101 
102  // --- Where are my base input files?
103 // std::string BaseCAFFileDir="/nova/app/users/karlwarb/Workspace/Ana2018/EnergyContainment/";
104 
105  // --- Make my iteration of MySpectra.
106  MySpectra sp;
107  // --- Loop through my input files...
108  for (unsigned int HC=0; HC<nHC; ++HC) {
109  // --- Open my files.
110  std::string FDName = "neutronKE_FD"+HCNames[HC]+".root"; TFile *FDFile = TFile::Open( FDName.c_str() );
111  std::string NDName = "neutronKE_ND"+HCNames[HC]+".root"; TFile *NDFile = TFile::Open( NDName.c_str() );
112  // --- Determine the POT
113  double POT = kAna2018FHCPOT;
114  if (HCNames[HC] == "rhc") POT = kAna2018RHCPOT;
115  // *** Loop through GENIE interaction types ***
116  for (unsigned int gen=0; gen<nGENIE; ++gen) {
117  // Loop through the quantile cuts.
118  for (unsigned int quant=0; quant<nQuantNames; ++quant) {
119  // --- Loop through the FD cuts first...
120  for(unsigned int det=0; det<nFD_CutNames; ++det) {
121  // Figure out what the end of my spectra names is.
122  std::string Comb = GENIEStr[gen]+"_"+FD_CutNames[det]+"_"+QuantNames[quant];
123  // Load my spectra
124  // 1D
125  // energy
126  sp.NeutrinoEnergy_FD [gen][det][quant][HC] = SpecToHist(FDFile, "NeutrinoE_"+Comb, POT, ";True #nu E (GeV);Events");
127  sp.TotPrimEnergy_FD [gen][det][quant][HC] = SpecToHist(FDFile, "TotPrimE_"+Comb, POT, ";Total prim E (GeV);Events");
128  sp.TotPrimEnergynKEp20_FD[gen][det][quant][HC] = SpecToHist(FDFile, "TotPrimEnKEp20_"+Comb, POT, ";Total prim E (GeV), n KE +20%;Events");
129  sp.TotPrimEnergynKEm20_FD[gen][det][quant][HC] = SpecToHist(FDFile, "TotPrimEnKEm20_"+Comb, POT, ";Total prim E (GeV), n KE -20%;Events");
130  sp.NeutMult_FD [gen][det][quant][HC] = SpecToHist(FDFile, "NeutMult"+Comb, POT, ";Neutron Multiplicity;Events");
131  sp.NeutEnergy_FD [gen][det][quant][HC] = SpecToHist(FDFile, "NeutEnergy"+Comb, POT, ";True Neutron KE (GeV);Events");
132  sp.NeutEFrac_FD [gen][det][quant][HC] = SpecToHist(FDFile, "NeutEFrac"+Comb, POT, ";True Neutron KE/True Neutrino Energy;Events");
133  sp.NeutPngCalE_FD [gen][det][quant][HC] = SpecToHist(FDFile, "NeutPngCalE"+Comb, POT, ";Neutron Prong calE (GeV);Events");
134  sp.NeutPngCalEpH_FD [gen][det][quant][HC] = SpecToHist(FDFile, "NeutPngCalEpH"+Comb, POT, ";Neutron Prong calE per Hit (MeV);Events");
135  sp.NeutPng2DCalE_FD [gen][det][quant][HC] = SpecToHist(FDFile, "NeutPng2DCalE"+Comb, POT, ";Neutron Prong (inc. 2D) calE (GeV);Events");
136  sp.NeutPng2DCalEpH_FD [gen][det][quant][HC] = SpecToHist(FDFile, "NeutPng2DCalEpH"+Comb, POT, ";Neutron Prong (inc.2D) calE per Hit (MeV);Events");
137  sp.NeutPngWtCalE_FD [gen][det][quant][HC] = SpecToHist(FDFile, "NeutPngWtCalE"+Comb, POT, ";Neutron Prong weighted calE;Events");
138  sp.NeutPngWtCalEpH_FD [gen][det][quant][HC] = SpecToHist(FDFile, "NeutPngWtCalEpH"+Comb, POT, ";Neutron Prong weighted calE per hit;Events");
139  sp.NeutPng2DWtCalE_FD [gen][det][quant][HC] = SpecToHist(FDFile, "NeutPng2DWtCalE"+Comb, POT, ";Neutron Prong (inc. 2D) weighted calE;Events");
140  sp.NeutPng2DWtCalEpH_FD [gen][det][quant][HC] = SpecToHist(FDFile, "NeutPng2DWtCalEpH"+Comb, POT, ";Neutron Prong (inc. 2D) weighted calE per hit;Events");
141  sp.NEFvsCVN_FD [gen][det][quant][HC] = Spec2DToHist(FDFile,"NEFvsCVN_"+Comb, POT, ";CVN muon id;True Neutron KE/True Neutrino E");
142  sp.NEFvsReMId_FD [gen][det][quant][HC] = Spec2DToHist(FDFile,"NEFvsReMId_"+Comb, POT, ";ReMId score;True Neutron KE/True Neutrino E");
143 
144  } // nFD_CutNames
145  // --- Now, loop through my ND cuts...
146  for(unsigned int det=0; det<nND_CutNames; ++det) {
147  // Figure out what the end of my spectra names is.
148  std::string Comb = GENIEStr[gen]+"_"+ND_CutNames[det]+"_"+QuantNames[quant];
149  // Load my spectra
150  // 1D
151  // energy
152  sp.NeutrinoEnergy_ND [gen][det][quant][HC] = SpecToHist(NDFile, "NeutrinoE_"+Comb, POT, ";True #nu E (GeV);Events");
153  sp.TotPrimEnergy_ND [gen][det][quant][HC] = SpecToHist(NDFile, "TotPrimE_"+Comb, POT, ";Total prim E (GeV);Events");
154  sp.TotPrimEnergynKEp20_ND[gen][det][quant][HC] = SpecToHist(NDFile, "TotPrimEnKEp20_"+Comb, POT, ";Total prim E (GeV), n KE +20%;Events");
155  sp.TotPrimEnergynKEm20_ND[gen][det][quant][HC] = SpecToHist(NDFile, "TotPrimEnKEm20_"+Comb, POT, ";Total prim E (GeV), n KE -20%;Events");
156  sp.NeutMult_ND [gen][det][quant][HC] = SpecToHist(NDFile, "NeutMult"+Comb, POT, ";Neutron Multiplicity;Events");
157  sp.NeutEnergy_ND [gen][det][quant][HC] = SpecToHist(NDFile, "NeutEnergy"+Comb, POT, ";True Neutron KE (GeV);Events");
158  sp.NeutEFrac_ND [gen][det][quant][HC] = SpecToHist(NDFile, "NeutEFrac"+Comb, POT, ";True Neutron KE/True Neutrino Energy;Events");
159  sp.NeutPngCalE_ND [gen][det][quant][HC] = SpecToHist(NDFile, "NeutPngCalE"+Comb, POT, ";Neutron Prong calE (GeV);Events");
160  sp.NeutPngCalEpH_ND [gen][det][quant][HC] = SpecToHist(NDFile, "NeutPngCalEpH"+Comb, POT, ";Neutron Prong calE per Hit (MeV);Events");
161  sp.NeutPng2DCalE_ND [gen][det][quant][HC] = SpecToHist(NDFile, "NeutPng2DCalE"+Comb, POT, ";Neutron Prong (inc. 2D) calE (GeV);Events");
162  sp.NeutPng2DCalEpH_ND [gen][det][quant][HC] = SpecToHist(NDFile, "NeutPng2DCalEpH"+Comb, POT, ";Neutron Prong (inc.2D) calE per Hit (MeV);Events");
163  sp.NeutPngWtCalE_ND [gen][det][quant][HC] = SpecToHist(NDFile, "NeutPngWtCalE"+Comb, POT, ";Neutron Prong weighted calE;Events");
164  sp.NeutPngWtCalEpH_ND [gen][det][quant][HC] = SpecToHist(NDFile, "NeutPngWtCalEpH"+Comb, POT, ";Neutron Prong weighted calE per hit;Events");
165  sp.NeutPng2DWtCalE_ND [gen][det][quant][HC] = SpecToHist(NDFile, "NeutPng2DWtCalE"+Comb, POT, ";Neutron Prong (inc. 2D) weighted calE;Events");
166  sp.NeutPng2DWtCalEpH_ND [gen][det][quant][HC] = SpecToHist(NDFile, "NeutPng2DWtCalEpH"+Comb, POT, ";Neutron Prong (inc. 2D) weighted calE per hit;Events");
167  sp.NEFvsCVN_ND [gen][det][quant][HC] = Spec2DToHist(NDFile,"NEFvsCVN_"+Comb, POT, ";CVN muon id;True Neutron KE/True Neutrino E");
168  sp.NEFvsReMId_ND [gen][det][quant][HC] = Spec2DToHist(NDFile,"NEFvsReMId_"+Comb, POT, ";ReMId score;True Neutron KE/True Neutrino E");
169 
170  } // nND_CutNames
171  } // nQuantNames
172  } // nGENIE
173  } // files
174 
175  // ...................................................... Make my combined plots ......................................................
176 
177  //final cut plots
178  for (unsigned int HC=0; HC<nHC; ++HC){
179  for (unsigned int gen=0; gen<nGENIE; ++gen) {
180  for (unsigned int quant=0; quant<nQuantNames; ++quant) {
181  // --- Now I can make the combination string...
182  std::string Comb = "_"+GENIEStr[gen]+"_"+QuantNames[quant]+"_"+HCNames[HC];
183  // --- FD First.
184  // Compare true neutrino energy and total energy from primaries
185  MakeComparisonPlots({sp.NeutrinoEnergy_FD[gen][FCut_FD][quant][HC], sp.TotPrimEnergy_FD[gen][FCut_FD][quant][HC]}, "NeutrinoandPrimaries_FD"+Comb, {"True #nu E", "Tot. E from prim."});
186  // put all primary combos on one plot
187  MakeComparisonPlots({sp.TotPrimEnergy_FD[gen][FCut_FD][quant][HC], sp.TotPrimEnergynKEp20_FD[gen][FCut_FD][quant][HC], sp.TotPrimEnergynKEm20_FD[gen][FCut_FD][quant][HC]}, "Primariespm20_FD"+Comb, {"Tot. E from prim.", "Tot. E from prim +20% nKE", "Tot. E from prim -20% nKE"});
188  // --- n/nu frac no pid vs. REMID, CVN
189  Print2D(sp.NEFvsCVN_FD[gen][0][quant][HC], "NEFvsCVN_FD_noPID"+Comb);
190  Print2D(sp.NEFvsReMId_FD[gen][0][quant][HC], "NEFvsReMId_FD_noPID"+Comb);
191 
192  // compare n/nu frac for full sel, !pid
193  MakeComparisonPlots({sp.NeutEFrac_FD[gen][2][quant][HC], sp.NeutEFrac_FD[gen][4][quant][HC]}, "NeutEFrac_PID"+Comb, {FD_CutNames[2], FD_CutNames[4]});
194  // compare n/nu frac for full sel, !cont
195  MakeComparisonPlots({sp.NeutEFrac_FD[gen][3][quant][HC], sp.NeutEFrac_FD[gen][4][quant][HC]}, "NeutEFrac_Cont"+Comb, {FD_CutNames[3], FD_CutNames[4]});
196 
197 
198  // --- ND second
199  // Compare true neutrino energy and total energy from primaries
200  MakeComparisonPlots({sp.NeutrinoEnergy_ND[gen][FCut_ND][quant][HC], sp.TotPrimEnergy_ND[gen][FCut_ND][quant][HC]}, "NeutrinoandPrimaries_ND"+Comb, {"True #nu E", "Tot. E from prim."});
201  // put all primary combos on one plot
202  MakeComparisonPlots({sp.TotPrimEnergy_ND[gen][FCut_ND][quant][HC], sp.TotPrimEnergynKEp20_ND[gen][FCut_ND][quant][HC], sp.TotPrimEnergynKEm20_ND[gen][FCut_ND][quant][HC]}, "Primariespm20_ND"+Comb, {"Tot. E from prim.", "Tot. E from prim +20% nKE", "Tot. E from prim -20% nKE"});
203  // --- n/nu frac no pid vs. REMID, CVN
204  // use -PID cut
205  Print2D(sp.NEFvsCVN_ND[gen][0][quant][HC], "NEFvsCVN_ND_noPID"+Comb);
206  Print2D(sp.NEFvsReMId_ND[gen][0][quant][HC], "NEFvsReMId_ND_noPID"+Comb);
207  // compare n/nu frac for full sel, !pid
208  MakeComparisonPlots({sp.NeutEFrac_ND[gen][2][quant][HC], sp.NeutEFrac_ND[gen][4][quant][HC]}, "NeutEFrac_PID"+Comb, {ND_CutNames[2], ND_CutNames[4]});
209  // compare n/nu frac for full sel, !cont
210  MakeComparisonPlots({sp.NeutEFrac_ND[gen][3][quant][HC], sp.NeutEFrac_ND[gen][4][quant][HC]}, "NeutEFrac_Cont"+Comb, {ND_CutNames[3], ND_CutNames[4]});
211 
212 
213  // --- Compare ND and FD --- need to area normalize
214  MakeComparisonPlots({sp.NeutMult_ND[gen][FCut_ND][quant][HC], sp.NeutMult_FD[gen][FCut_FD][quant][HC]}, "NeutMult_NDFD"+Comb, DetNames);
215  MakeComparisonPlots({sp.NeutEnergy_ND[gen][FCut_ND][quant][HC], sp.NeutEnergy_FD[gen][FCut_FD][quant][HC]}, "NeutEnergy_NDFD"+Comb, DetNames);
216  MakeComparisonPlots({sp.NeutEFrac_ND[gen][FCut_ND][quant][HC], sp.NeutEFrac_FD[gen][FCut_FD][quant][HC]}, "NeutEFrac_NDFD"+Comb, DetNames);
217  MakeComparisonPlots({sp.NeutPngCalE_ND[gen][FCut_ND][quant][HC], sp.NeutPngCalE_FD[gen][FCut_FD][quant][HC]}, "NeutPngCalE_NDFD"+Comb, DetNames);
218  MakeComparisonPlots({sp.NeutPngCalEpH_ND[gen][FCut_ND][quant][HC], sp.NeutPngCalEpH_FD[gen][FCut_FD][quant][HC]}, "NeutPngCalEpH_NDFD"+Comb, DetNames);
219  MakeComparisonPlots({sp.NeutPng2DCalE_ND[gen][FCut_ND][quant][HC], sp.NeutPng2DCalE_FD[gen][FCut_FD][quant][HC]}, "NeutPng2DCalE_NDFD"+Comb, DetNames);
220  MakeComparisonPlots({sp.NeutPng2DCalEpH_ND[gen][FCut_ND][quant][HC], sp.NeutPng2DCalEpH_FD[gen][FCut_FD][quant][HC]}, "NeutPng2DCalEpH_NDFD"+Comb, DetNames);
221  MakeComparisonPlots({sp.NeutPngWtCalE_ND[gen][FCut_ND][quant][HC], sp.NeutPngWtCalE_FD[gen][FCut_FD][quant][HC]}, "NeutPngWtCalE_NDFD"+Comb, DetNames);
222  MakeComparisonPlots({sp.NeutPngWtCalEpH_ND[gen][FCut_ND][quant][HC], sp.NeutPngWtCalEpH_FD[gen][FCut_FD][quant][HC]}, "NeutPngWtCalEpH_NDFD"+Comb, DetNames);
223  MakeComparisonPlots({sp.NeutPng2DWtCalE_ND[gen][FCut_ND][quant][HC], sp.NeutPng2DWtCalE_FD[gen][FCut_FD][quant][HC]}, "NeutPng2DWtCalE_NDFD"+Comb, DetNames);
224  MakeComparisonPlots({sp.NeutPng2DWtCalEpH_ND[gen][FCut_ND][quant][HC], sp.NeutPng2DWtCalEpH_FD[gen][FCut_FD][quant][HC]}, "NeutPng2DWtCalEpH_NDFD"+Comb, DetNames);
225 
226  } // Loop over quants.
227  } // Loop through GENIE interaction types.
228  }// FHC/RHC
229 
230  // compare FHC and RHC
231  for(unsigned int gen=0; gen<nGENIE; ++gen) {
232  for (unsigned int quant=0; quant<nQuantNames; ++quant) {
233  std::string Comb = "_"+GENIEStr[gen]+"_"+QuantNames[quant];
234  // --- ND
235  // --- all neutron characteristics
236  MakeComparisonPlots({sp.NeutMult_ND[gen][FCut_ND][quant][0], sp.NeutMult_ND[gen][FCut_ND][quant][1]}, "NeutMult_ND_FHCRHC"+Comb, HCNames);
237  MakeComparisonPlots({sp.NeutEnergy_ND[gen][FCut_ND][quant][0], sp.NeutEnergy_ND[gen][FCut_ND][quant][1]}, "NeutEnergy_ND_FHCRHC"+Comb, HCNames);
238  MakeComparisonPlots({sp.NeutEFrac_ND[gen][FCut_ND][quant][0], sp.NeutEFrac_ND[gen][FCut_ND][quant][1]}, "NeutEFrac_ND_FHCRHC"+Comb, HCNames);
239  MakeComparisonPlots({sp.NeutPngCalE_ND[gen][FCut_ND][quant][0],sp.NeutPngCalE_ND[gen][FCut_ND][quant][1]}, "NeutPngCalE_ND_FHCRHC"+Comb, HCNames);
240  MakeComparisonPlots({sp.NeutPngCalEpH_ND[gen][FCut_ND][quant][0], sp.NeutPngCalEpH_ND[gen][FCut_ND][quant][1]}, "NeutPngCalEpH_ND_FHCRHC"+Comb, HCNames);
241  MakeComparisonPlots({sp.NeutPng2DCalE_ND[gen][FCut_ND][quant][0], sp.NeutPng2DCalE_ND[gen][FCut_ND][quant][1]}, "NeutPng2DCalE_ND_FHCRHC"+Comb, HCNames);
242  MakeComparisonPlots({sp.NeutPng2DCalEpH_ND[gen][FCut_ND][quant][0], sp.NeutPng2DCalEpH_ND[gen][FCut_ND][quant][1]}, "NeutPng2DCalEpH_ND_FHCRHC"+Comb, HCNames);
243  MakeComparisonPlots({sp.NeutPngWtCalE_ND[gen][FCut_ND][quant][0], sp.NeutPngWtCalE_ND[gen][FCut_ND][quant][1]}, "NeutPngWtCalE_ND_FHCRHC"+Comb, HCNames);
244  MakeComparisonPlots({sp.NeutPngWtCalEpH_ND[gen][FCut_ND][quant][0], sp.NeutPngWtCalEpH_ND[gen][FCut_ND][quant][1]}, "NeutPngWtCalEpH_ND_FHCRHC"+Comb, HCNames);
245  MakeComparisonPlots({sp.NeutPng2DWtCalE_ND[gen][FCut_ND][quant][0], sp.NeutPng2DWtCalE_ND[gen][FCut_ND][quant][1]}, "NeutPng2DWtCalE_ND_FHCRHC"+Comb, HCNames);
246  MakeComparisonPlots({sp.NeutPng2DWtCalEpH_ND[gen][FCut_ND][quant][0], sp.NeutPng2DWtCalEpH_ND[gen][FCut_ND][quant][1]}, "NeutPng2DWtCalEpH_ND_FHCRHC"+Comb, HCNames);
247 
248  // FD
249  // --- all neutron characteristics
250  MakeComparisonPlots({sp.NeutMult_FD[gen][FCut_FD][quant][0], sp.NeutMult_FD[gen][FCut_FD][quant][1]}, "NeutMult_FD_FHCRHC"+Comb, HCNames);
251  MakeComparisonPlots({sp.NeutEnergy_FD[gen][FCut_FD][quant][0], sp.NeutEnergy_FD[gen][FCut_FD][quant][1]}, "NeutEnergy_FD_FHCRHC"+Comb, HCNames);
252  MakeComparisonPlots({sp.NeutEFrac_FD[gen][FCut_FD][quant][0], sp.NeutEFrac_FD[gen][FCut_FD][quant][1]}, "NeutEFrac_FD_FHCRHC"+Comb, HCNames);
253  MakeComparisonPlots({sp.NeutPngCalE_FD[gen][FCut_FD][quant][0],sp.NeutPngCalE_FD[gen][FCut_FD][quant][1]}, "NeutPngCalE_FD_FHCRHC"+Comb, HCNames);
254  MakeComparisonPlots({sp.NeutPngCalEpH_FD[gen][FCut_FD][quant][0], sp.NeutPngCalEpH_FD[gen][FCut_FD][quant][1]}, "NeutPngCalEpH_FD_FHCRHC"+Comb, HCNames);
255  MakeComparisonPlots({sp.NeutPng2DCalE_FD[gen][FCut_FD][quant][0], sp.NeutPng2DCalE_FD[gen][FCut_FD][quant][1]}, "NeutPng2DCalE_FD_FHCRHC"+Comb, HCNames);
256  MakeComparisonPlots({sp.NeutPng2DCalEpH_FD[gen][FCut_FD][quant][0], sp.NeutPng2DCalEpH_FD[gen][FCut_FD][quant][1]}, "NeutPng2DCalEpH_FD_FHCRHC"+Comb, HCNames);
257  MakeComparisonPlots({sp.NeutPngWtCalE_FD[gen][FCut_FD][quant][0], sp.NeutPngWtCalE_FD[gen][FCut_FD][quant][1]}, "NeutPngWtCalE_FD_FHCRHC"+Comb, HCNames);
258  MakeComparisonPlots({sp.NeutPngWtCalEpH_FD[gen][FCut_FD][quant][0], sp.NeutPngWtCalEpH_FD[gen][FCut_FD][quant][1]}, "NeutPngWtCalEpH_FD_FHCRHC"+Comb, HCNames);
259  MakeComparisonPlots({sp.NeutPng2DWtCalE_FD[gen][FCut_FD][quant][0], sp.NeutPng2DWtCalE_FD[gen][FCut_FD][quant][1]}, "NeutPng2DWtCalE_FD_FHCRHC"+Comb, HCNames);
260  MakeComparisonPlots({sp.NeutPng2DWtCalEpH_FD[gen][FCut_FD][quant][0], sp.NeutPng2DWtCalEpH_FD[gen][FCut_FD][quant][1]}, "NeutPng2DWtCalEpH_FD_FHCRHC"+Comb, HCNames);
261 
262  }
263  }
264 
265 
266  return;
267 }
268 
269 //......................................................
270 TH1D* SpecToHist( TFile* MyF, std::string SpecName, double POT, std::string AxisTitle ) {
271  std::unique_ptr<Spectrum> TempSpec = Spectrum::LoadFrom(MyF, SpecName.c_str() ) ;
272  TH1D* MyHist = TempSpec -> ToTH1( POT );
273  ana::CenterTitles ( MyHist);
274  MyHist->SetMinimum( 0 );
275  MyHist->SetMaximum( MyHist->GetMaximum()*1.1 );
276  MyHist->SetTitle ( AxisTitle.c_str() );
277  //area norm
278 // MyHist->Scale(1./MyHist->Integral());
279 
280  // I want to check what my maximum bin is...
281  // this is clever but if we're trying to compare two sets of things and one has the zero bin and one doesn't it becomes difficult
282 // int MaxBin = MyHist->GetMaximumBin();
283 // int MaxVal = MyHist->GetMaximum();
284 // int SecBin = MyHist->GetBinContent(2);
285 // if ( MaxBin == 1 && (MaxVal/10 > SecBin) ) {
286 // double NewMax = 0;
287 // for (int bin=2; bin<MyHist->GetNbinsX(); ++bin) {
288 // if (MyHist->GetBinContent(bin) > NewMax) NewMax = 1.1*MyHist->GetBinContent(bin);
289 // }
290 // MyHist->SetMaximum( NewMax );
291 //
292 // }
293 
294  // comment this out to keep the 0bin
295 // MyHist->GetXaxis()->SetRange(2,MyHist->GetNbinsX());
296 
297  return MyHist;
298 }
299 
300 //......................................................
301 TH2* Spec2DToHist( TFile* MyF, std::string SpecName, double POT, std::string AxisTitle ) {
302  std::unique_ptr<Spectrum> TempSpec = Spectrum::LoadFrom(MyF, SpecName.c_str() ) ;
303  TH2* MyHist = TempSpec -> ToTH2( POT );
304  ana::CenterTitles ( MyHist);
305  MyHist->SetMinimum( 0 );
306  MyHist->SetMaximum( MyHist->GetMaximum()*1.1 );
307  MyHist->SetTitle ( AxisTitle.c_str() );
308 
309 // // borrowed from Michael -- thanks!!
310 // // Loop over the histos to normalize them in vertical strips one bin wide.
311 // unsigned int NXbins = 0.0;
312 // unsigned int NYbins = 0.0;
313 //
314 // // W0, X0
315 // NXbins = MyHist->GetNbinsX();
316 // NYbins = MyHist->GetNbinsY();
317 //
318 // for(unsigned int i = 0; i <= NXbins+1; ++i) {
319 // float integral = 0.0;
320 // for(unsigned int j = 2; j <= NYbins; ++j) // 0 bin swamping the plot... remove the 0bin
321 // integral += MyHist->GetBinContent(i,j);
322 // if(integral == 0.0) continue;
323 // for(unsigned int j = 0; j <= NYbins+1; ++j)
324 // MyHist->SetBinContent(i,j,MyHist->GetBinContent(i,j)/integral);
325 // }
326 //
327 // // the 0 bin is swamping the plot... remove it
328 // MyHist->GetYaxis()->SetRange(2,NYbins);
329 //
330 // // since it's normalized, need to set the z-axis to go from 0-1
331 // MyHist->GetZaxis()->SetRangeUser(0.0, 1.0);
332  return MyHist;
333 }
334 
335 //......................................................
336 void FitGauss(TH1D* h)
337 {
338  double mean, bin1, bin2, fwhm;
339  mean = h->GetMean();
340  bin1 = h->FindFirstBinAbove(h->GetMaximum()/2);
341  bin2 = h->FindLastBinAbove(h->GetMaximum()/2);
342  fwhm = h->GetBinCenter(bin2 ) - h->GetBinCenter(bin1 );
343 
344  int ilo = 1;
345  int ihi = h->GetNbinsX();
346 
347  int imx = h->GetMaximumBin();
348  double mx = h->GetBinContent(imx);
349 
350  int ileft;
351  for (ileft=imx; ileft>=ilo; --ileft){
352  if (h->GetBinContent(ileft)<0.5*mx) {
353  ++ileft;
354  break;
355  }
356  }
357 
358  int iright;
359  for (iright=imx; iright<=ihi; ++iright) {
360  if (h->GetBinContent(iright)<0.5*mx) {
361  break;
362  }
363  }
364 
365  TF1 *myFit = new TF1("myFit","gaus", h->GetBinLowEdge(ileft), h->GetBinLowEdge(iright));
366  h->Fit(myFit,"R q");
367  myFit->SetLineColor(h->GetLineColor());
368 // myFit->Draw("same");
369  std::cout<<"Mean: " << myFit->GetParameter(1) << "+/-" << myFit->GetParError(1) <<std::endl;
370  return;
371 }
372 
373 //......................................................
374 void MakeComparisonPlots( std::vector<TH1D*> TheHists, std::string PlotName, std::vector<std::string> Labels ) {
375 
376  // --- Make my canvas!
377  TCanvas* ThisCan = new TCanvas( PlotName.c_str(), PlotName.c_str() );
378  // --- Legend time!
379  double HistLeg_X1 = 0.57, HistLeg_Y1 = 0.65, HistLeg_X2 = 0.88, HistLeg_Y2 = 0.88;
380  TLegend* Leg = new TLegend(HistLeg_X1, HistLeg_Y1, HistLeg_X2, HistLeg_Y2 );
381 
382  double max = 0.0;
383  // --- Loop through hists and plot them!
384  for (size_t hist=0; hist<TheHists.size(); ++hist) {
385  TheHists[hist] -> SetLineColor( 1+hist );
386 
387  if(TheHists[hist]->GetMaximum() > max) max = TheHists[hist]->GetMaximum();
388  if (hist==0) {
389  TheHists[hist] -> Draw("hist");
390 // std::cout<<"Label: "<<Labels[hist]<<std::endl;
391 // FitGauss(TheHists[hist]);
392  }
393  else{
394  TheHists[hist] -> Draw("hist same");
395 // std::cout<<"Label: "<<Labels[hist]<<std::endl;
396 // FitGauss(TheHists[hist]);
397  }
398  Leg -> AddEntry( TheHists[hist], Labels[hist].c_str() );
399  }
400  Leg -> SetFillStyle(0);
401  Leg -> SetBorderSize(0);
402 
403  Leg -> Draw();
404  // set the new maximum
405  TheHists[0]->SetMaximum(max*1.1);
406 
407  // --- Do I want to have a log axis?
408  if ( PlotName.find("CompareUncontainEnergies") != std::string::npos) {
409  for (size_t hist=0; hist<TheHists.size(); ++hist) TheHists[hist] -> GetYaxis() -> SetRangeUser( 10, TheHists[0]->GetMaximum() );
410  ThisCan -> SetLogy();
411  }
412 
413 
414  // --- Save the canvas.
415  ThisCan -> SaveAs( TString("Plots/")+TString(ThisCan->GetName())+TString(".pdf") );
416  ThisCan -> SaveAs( TString("Plots/")+TString(ThisCan->GetName())+TString(".png") );
417 
418  // return
419  return;
420 }
421 
422 void Print2D(TH2* hist, std::string PlotName)
423 {
424  // --- Make my canvas!
425  TCanvas* ThisCan = new TCanvas( PlotName.c_str(), PlotName.c_str() );
426  hist->Draw("colz");
427  ThisCan -> SaveAs( TString("Plots/")+TString(ThisCan->GetName())+TString(".pdf") );
428  ThisCan -> SaveAs( TString("Plots/")+TString(ThisCan->GetName())+TString(".png") );
429  return;
430 
431 }
432 
TH1D * NeutPng2DWtCalEpH_ND[nGENIE][nND_CutNames][nQuantNames][nHC]
Definition: neutronE_plot.C:83
tree Draw("slc.nhit")
TH1D * TotPrimEnergy_ND[nGENIE][nND_CutNames][nQuantNames][nHC]
Definition: neutronE_plot.C:70
TH1D * TotPrimEnergy_FD[nGENIE][nFD_CutNames][nQuantNames][nHC]
Definition: neutronE_plot.C:50
unsigned int FCut_FD
Definition: neutronE_plot.C:35
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
TH1D * TotPrimEnergynKEm20_ND[nGENIE][nND_CutNames][nQuantNames][nHC]
Definition: neutronE_plot.C:72
std::vector< std::string > DetNames
Definition: neutronE_plot.C:44
TH1D * TotPrimEnergynKEp20_ND[nGENIE][nND_CutNames][nQuantNames][nHC]
Definition: neutronE_plot.C:71
correl_xv GetYaxis() -> SetDecimals()
TH1D * NeutEnergy_ND[nGENIE][nND_CutNames][nQuantNames][nHC]
Definition: neutronE_plot.C:74
TH1D * SpecToHist(TFile *MyF, std::string SpecName, double POT, std::string AxisTitle)
void FitGauss(TH1D *h)
const unsigned int nGENIE
Definition: neutronE_plot.C:29
const unsigned int nHC
Definition: neutronE_plot.C:41
TH1D * NeutMult_FD[nGENIE][nFD_CutNames][nQuantNames][nHC]
Definition: neutronE_plot.C:53
TH1D * NeutPng2DWtCalEpH_FD[nGENIE][nFD_CutNames][nQuantNames][nHC]
Definition: neutronE_plot.C:63
TH1D * NeutPngCalEpH_FD[nGENIE][nFD_CutNames][nQuantNames][nHC]
Definition: neutronE_plot.C:57
TH1D * NeutEFrac_FD[nGENIE][nFD_CutNames][nQuantNames][nHC]
Definition: neutronE_plot.C:55
void CenterTitles(TH1 *histo)
Definition: Plots.cxx:1483
TH1D * NeutPngWtCalEpH_FD[nGENIE][nFD_CutNames][nQuantNames][nHC]
Definition: neutronE_plot.C:61
TH1D * NeutrinoEnergy_FD[nGENIE][nFD_CutNames][nQuantNames][nHC]
ntuple SetFillStyle(1001)
TH2 * Spec2DToHist(TFile *MyF, std::string SpecName, double POT, std::string AxisTitle)
const unsigned int nQuantNames
Definition: neutronE_plot.C:38
TH1D * NeutPng2DCalEpH_FD[nGENIE][nFD_CutNames][nQuantNames][nHC]
Definition: neutronE_plot.C:59
const double kAna2018RHCPOT
Definition: Exposures.h:208
TH1D * NeutPngCalE_ND[nGENIE][nND_CutNames][nQuantNames][nHC]
Definition: neutronE_plot.C:76
TH1D * TotPrimEnergynKEm20_FD[nGENIE][nFD_CutNames][nQuantNames][nHC]
Definition: neutronE_plot.C:52
unsigned int FCut_ND
Definition: neutronE_plot.C:35
legend SetBorderSize(0)
TH1D * NeutPng2DWtCalE_ND[nGENIE][nND_CutNames][nQuantNames][nHC]
Definition: neutronE_plot.C:82
TH1D * NeutPngWtCalE_FD[nGENIE][nFD_CutNames][nQuantNames][nHC]
Definition: neutronE_plot.C:60
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:563
TH1D * NeutPng2DCalEpH_ND[nGENIE][nND_CutNames][nQuantNames][nHC]
Definition: neutronE_plot.C:79
TH1D * NeutrinoEnergy_ND[nGENIE][nND_CutNames][nQuantNames][nHC]
leg AddEntry(GRdata,"data","p")
hmean SetLineColor(4)
TH1D * NeutEFrac_ND[nGENIE][nND_CutNames][nQuantNames][nHC]
Definition: neutronE_plot.C:75
TH1D * NeutPng2DWtCalE_FD[nGENIE][nFD_CutNames][nQuantNames][nHC]
Definition: neutronE_plot.C:62
std::vector< std::string > ND_CutNames
Definition: neutronE_plot.C:34
const unsigned int nFD_CutNames
Definition: neutronE_plot.C:33
TH1D * NeutMult_ND[nGENIE][nND_CutNames][nQuantNames][nHC]
Definition: neutronE_plot.C:73
std::vector< std::string > HCNames
Definition: neutronE_plot.C:41
std::vector< std::string > FD_CutNames
Definition: neutronE_plot.C:33
TH2 * NEFvsReMId_FD[nGENIE][nFD_CutNames][nQuantNames][nHC]
Definition: neutronE_plot.C:65
OStream cout
Definition: OStream.cxx:6
std::vector< double > POT
TH1D * TotPrimEnergynKEp20_FD[nGENIE][nFD_CutNames][nQuantNames][nHC]
Definition: neutronE_plot.C:51
TH1D * NeutPngCalE_FD[nGENIE][nFD_CutNames][nQuantNames][nHC]
Definition: neutronE_plot.C:56
TH2 * NEFvsCVN_FD[nGENIE][nFD_CutNames][nQuantNames][nHC]
Definition: neutronE_plot.C:64
TH1D * NeutPngCalEpH_ND[nGENIE][nND_CutNames][nQuantNames][nHC]
Definition: neutronE_plot.C:77
TH1D * NeutEnergy_FD[nGENIE][nFD_CutNames][nQuantNames][nHC]
Definition: neutronE_plot.C:54
std::vector< std::string > QuantNames
Definition: neutronE_plot.C:38
const unsigned int nND_CutNames
Definition: neutronE_plot.C:34
void Print2D(TH2 *hist, std::string PlotName)
TH2 * NEFvsReMId_ND[nGENIE][nND_CutNames][nQuantNames][nHC]
Definition: neutronE_plot.C:85
const double kAna2018FHCPOT
Definition: Exposures.h:207
void neutronE_plot()
Definition: neutronE_plot.C:97
TH2 * NEFvsCVN_ND[nGENIE][nND_CutNames][nQuantNames][nHC]
Definition: neutronE_plot.C:84
cosmicTree SaveAs("cosmicTree.root")
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:66
void MakeComparisonPlots(std::vector< TH1D * > TheHists, std::string PlotName, std::vector< std::string > Labels)
TH1D * NeutPngWtCalEpH_ND[nGENIE][nND_CutNames][nQuantNames][nHC]
Definition: neutronE_plot.C:81
TH2 * ToTH2(const Spectrum &s, double exposure, ana::EExposureType expotype, const Binning &binsx, const Binning &binsy, ana::EBinType bintype)
For use with Var2D.
Definition: UtilsExt.cxx:115
TH1D * NeutPng2DCalE_FD[nGENIE][nFD_CutNames][nQuantNames][nHC]
Definition: neutronE_plot.C:58
TH1D * NeutPngWtCalE_ND[nGENIE][nND_CutNames][nQuantNames][nHC]
Definition: neutronE_plot.C:80
gPad SetLogy()
TH1D * NeutPng2DCalE_ND[nGENIE][nND_CutNames][nQuantNames][nHC]
Definition: neutronE_plot.C:78
const std::vector< std::string > GENIEStr
Definition: neutronE_plot.C:29
enum BeamMode string