plot_nuTree.C
Go to the documentation of this file.
4 #include "CAFAna/Core/Ratio.h"
6 #include "CAFAna/Core/Spectrum.h"
8 #include "CAFAna/Vars/Vars.h"
10 #include "CAFAna/Extrap/IExtrap.h"
11 using namespace ana;
12 
13 #include "Utilities/func/MathUtil.h"
14 
15 #include "OscLib/IOscCalc.h"
16 #include "OscLib/OscCalcDumb.h"
17 #include "CAFAna/Analysis/Calcs.h"
19 
20 #include "TCanvas.h"
21 #include "TFile.h"
22 #include "TH1.h"
23 #include "TH2.h"
24 #include "THStack.h"
25 #include <iostream>
26 #include <cstring>
27 #include "TStyle.h"
28 #include "TLegend.h"
29 #include "TPad.h"
30 #include "Rtypes.h"//syu
31 #include "TPaveStats.h"
32 #include "../Utils/MyPlot.h"
33 #include "Utilities/rootlogon.C"
34 
36 {
37  struct HistDef
38  {
40  };
41 
42  int kNumVars =3;
43  HistDef defs[3] = {
44 "trueE",
45 "trueEFull",
46 "recoE",
47  };
48 
49  if(std::strstr(fname.c_str(),"nuTree")){
50  kNumVars =2;
51  }
52 
53  std::vector<int> StackColor = {kRed,kAzure, kGreen+2, kBlack, kOrange+7, kPink+9, kViolet -5,49 };
54 
55 
56  const int kNumComps = 5;
57  const int kNumAnas = 1;
58  std::string compNames[kNumComps] = {"Total","FromPi","FromKa","FromPro","FromOthers"};
59 
60  const std::string anaNames[kNumAnas] = {/*"Total","Signal","BeamNue",*/"NumuCC"/*,"NC"*/};
62 
63  int LineStyles[2]={1,7};
64  Spectrum* FDNu[kNumVars][kNumComps][kNumAnas];
65  Spectrum* NDNu[kNumVars][kNumComps][kNumAnas];
66 
67  Spectrum* FDNuBar[kNumVars][kNumComps][kNumAnas];
68  Spectrum* NDNuBar[kNumVars][kNumComps][kNumAnas];
69 
70  IPrediction* preds[kNumVars][kNumComps];
71  IPrediction* prednue[kNumVars][kNumComps];
72  IPrediction* prednumu[kNumVars][kNumComps];
73  IExtrap* signue[kNumVars][kNumComps];
74  osc::IOscCalc* nocalc = new osc::NoOscillations;
75  string dirname;
76 
77  for(int varIdx = 0; varIdx < kNumVars; ++varIdx){
78  const char* varName = defs[varIdx].name.c_str();
79  for(int compIdx = 0; compIdx<kNumComps; ++compIdx){
80  const char* compname = compNames[compIdx].c_str();
81  if(std::strstr(fname.c_str(),"recTree")){
82  preds[varIdx][compIdx] = LoadFromFile<IPrediction>(fname, TString::Format("%s/%s",varName,compname).Data()).release();
83  prednue[varIdx][compIdx] = LoadFromFile<IPrediction>(fname, TString::Format("%s/ndnue_%s",varName,compname).Data()).release();
84  prednumu[varIdx][compIdx] = LoadFromFile<IPrediction>(fname, TString::Format("%s/ndnumu_%s",varName,compname).Data()).release();
85  signue[varIdx][compIdx] = LoadFromFile<IExtrap>(fname, TString::Format("%s/%s/extrap",varName,compname).Data()).release();
86  }
87 
88  for(int anaIdx = 0; anaIdx < kNumAnas; ++anaIdx){
89  const char* ananame = anaNames[anaIdx].c_str();
90 
91  if(std::strstr(fname.c_str(),"nuTree")){//load nu spec
92  dirname = "nuTree";
93 std::cout<< "loading nuTree..."<<std::endl;
94 
95  if(false)//anaIdx ==0) //load signal
96  {
97  FDNu[varIdx][compIdx][anaIdx] = LoadFromFile<Spectrum>(fname, TString::Format("%s/FDFS_%s_Signal_Nu",varName,compname).Data()).release();
98  FDNuBar[varIdx][compIdx][anaIdx] = LoadFromFile<Spectrum>(fname, TString::Format("%s/FDFS_%s_Signal_NuBar",varName,compname).Data()).release();
99  NDNu[varIdx][compIdx][anaIdx] = LoadFromFile<Spectrum>(fname, TString::Format("%s/ND_%s_NumuCC_Nu",varName,compname).Data()).release();
100  NDNuBar[varIdx][compIdx][anaIdx] = LoadFromFile<Spectrum>(fname, TString::Format("%s/ND_%s_NumuCC_NuBar",varName,compname).Data()).release();
101  }
102  else{
103  FDNu[varIdx][compIdx][anaIdx] = LoadFromFile<Spectrum>(fname, TString::Format("%s/FDNS_%s_%s_Nu",varName,compname,ananame).Data()).release();
104  FDNuBar[varIdx][compIdx][anaIdx] = LoadFromFile<Spectrum>(fname, TString::Format("%s/FDNS_%s_%s_NuBar",varName,compname,ananame).Data()).release();
105 
106  NDNu[varIdx][compIdx][anaIdx] = LoadFromFile<Spectrum>(fname, TString::Format("%s/ND_%s_%s_Nu",varName,compname,ananame).Data()).release();
107  NDNuBar[varIdx][compIdx][anaIdx] = LoadFromFile<Spectrum>(fname, TString::Format("%s/ND_%s_%s_NuBar",varName,compname,ananame).Data()).release();
108  }
109  }
110 
111  else{
112  dirname = "recTree";
113 //hack hack
114  if (anaIdx == 0){
115  FDNu[varIdx][compIdx][anaIdx] = new Spectrum ((signue[varIdx][compIdx]->NueAppComponent()).Unoscillated());
116  FDNuBar[varIdx][compIdx][anaIdx] = new Spectrum ((signue[varIdx][compIdx]->AntiNueAppComponent()).Unoscillated());
117  NDNu[varIdx][compIdx][anaIdx] =new Spectrum (prednumu[varIdx][compIdx]->PredictComponent(nocalc,
118  Flavors::kNuMuToNuMu,Current::kCC,Sign::kNu));
119  NDNuBar[varIdx][compIdx][anaIdx] =new Spectrum (prednumu[varIdx][compIdx]->PredictComponent(nocalc,
120  Flavors::kNuMuToNuMu,Current::kCC,Sign::kAntiNu));
121  }
122  else{
123  FDNu[varIdx][compIdx][anaIdx] = new Spectrum(preds[varIdx][compIdx] ->PredictComponent(nocalc,
124  flavors[anaIdx],
125  Current::kCC,
126  Sign::kNu));
127  FDNuBar[varIdx][compIdx][anaIdx] = new Spectrum(preds[varIdx][compIdx] ->PredictComponent(nocalc,
128  flavors[anaIdx],
129  Current::kCC,
130  Sign::kAntiNu));
131  NDNu[varIdx][compIdx][anaIdx] = new Spectrum(prednue[varIdx][compIdx] ->PredictComponent(nocalc,
132  flavors[anaIdx],
133  Current::kCC,
134  Sign::kNu));
135  NDNuBar[varIdx][compIdx][anaIdx] = new Spectrum( prednue[varIdx][compIdx] ->PredictComponent(nocalc,
136  flavors[anaIdx],
137  Current::kCC,
138  Sign::kAntiNu));
139  }
140  }
141  }
142  }
143  }
144 
145 //---- ana bins
146  for(int varIdx = 0; varIdx < kNumVars; ++varIdx){
147  const std::string varName = defs[varIdx].name;
148 
149  for(int anaIdx=0; anaIdx<kNumAnas; ++anaIdx){//all cc
150  std::string ananame = anaNames[anaIdx];
151 rootlogon();
152  TCanvas* CompNu = new TCanvas();
153  TLegend* LegCompNu = new TLegend(.5,.6,.95,.9);
154  TCanvas* CompNuBar = new TCanvas();
155  TLegend* LegCompNuBar = new TLegend(.5,.6,.95,.9);
156 
157  TCanvas* NDCompNu = new TCanvas();
158  TLegend* NDLegCompNu = new TLegend(.45,.45,.9,.9);
159  TCanvas* NDCompNuBar = new TCanvas();
160  TLegend* NDLegCompNuBar = new TLegend(.45,.45,.9,.9);
161 rootlogon();
162 
163  TCanvas* FNNuRatios = new TCanvas();
164  TLegend* FNNuRatiosLeg = new TLegend(.7,.6,.95,.9);
165  FNNuRatiosLeg->SetFillStyle(0);
166  TH1* hSumNu = new TH1F();
167  TCanvas* FNNuBarRatios = new TCanvas();
168  TLegend* FNNuBarRatiosLeg = new TLegend(.7,.6,.95,.9);
169  FNNuBarRatiosLeg->SetFillStyle(0);
170  TH1* hSumNuBarF = new TH1F();
171  TH1* hSumNuBarN = new TH1F();
172 
173  for(int compIdx =0; compIdx<kNumComps; ++compIdx){// different parents
174  std::string compname= compNames[compIdx];
175 rootlogon();
176  // FD and ND all ancestors in one plots, 4 in total
177  CompNu->cd();
178  TH1* hFDNuComp = OneCompPlot(FDNu[varIdx][compIdx][anaIdx], " FD Nu Decomp by Ancestors", ananame+" from "+compname, StackColor[compIdx],7, LegCompNu);
179  if(compIdx==0) hFDNuComp->Draw("hist");
180  else hFDNuComp->Draw("hist same");
181 
182  CompNuBar->cd();
183  TH1* hFDNuBarComp = OneCompPlot(FDNuBar[varIdx][compIdx][anaIdx], " FD NuBar Decomp by Ancestors", ananame+"Bar from "+compname, StackColor[compIdx],1, LegCompNuBar);
184  if(compIdx==0) hFDNuBarComp->Draw("hist");
185  else hFDNuBarComp->Draw("hist same");
186 
187  NDCompNu->cd();
188 // TH1* hNDNuComp = OneCompPlot(NDNu[varIdx][compIdx][anaIdx], "ND Nu Decomp by Ancestors", ananame+" from "+compname, StackColor[compIdx],7, NDLegCompNu);
189  TH1* hNDNuComp = OneCompPlot(NDNu[varIdx][compIdx][anaIdx], "ND Nu Decomp by Ancestors", ananame+" from "+compname, StackColor[compIdx],1, NDLegCompNu);
190  if(compIdx==0) hNDNuComp->Draw("hist");
191  else hNDNuComp->Draw("hist same");
192 
193  NDCompNuBar->cd();
194  TH1* hNDNuBarComp = OneCompPlot(NDNuBar[varIdx][compIdx][anaIdx], "ND NuBar Decomp by Ancestors", ananame+"Bar from "+compname, StackColor[compIdx],1, NDLegCompNuBar);
195  if(compIdx==0) hNDNuBarComp->Draw("hist");
196  else hNDNuBarComp->Draw("hist same");
197 
198  // FD/ND ratio for each ancestor, comparing Nu and NuBar
199  TCanvas* FN = new TCanvas();
200  TLegend* LegFN = new TLegend(.6, .6, .95 ,.9);
201  TH1* FN_Nu = Ratio(FDNu[varIdx][compIdx][anaIdx], NDNu[varIdx][compIdx][anaIdx],"FD/ND", compname+" Ancestor", ananame+" Nu",StackColor[compIdx],7, LegFN );
202  TH1* FN_NuBar = Ratio(FDNuBar[varIdx][compIdx][anaIdx], NDNuBar[varIdx][compIdx][anaIdx],"FD/ND", compname+" Ancestor", ananame+" NuBar",StackColor[compIdx],1, LegFN );
203  DrawHists(FN_Nu, FN_NuBar);
204  LegFN->Draw();
205  rootlogon();
206  gPad->Update();
207  FN->SaveAs(("plots/"+dirname+"/FN_ratio_"+varName+"_"+compname+"_"+ananame+"_ancestor.pdf").c_str());
208 
209  // FD/ND all ratios in 1 plot
210  FNNuRatios->cd();
211  if(compIdx == 0 ) FN_Nu->Draw("hist");
212  else FN_Nu->Draw("hist same");
213  FNNuRatiosLeg->AddEntry(FN_Nu, compname.c_str(), "l");
214 
215  FNNuBarRatios->cd();
216  if(compIdx == 0 ) FN_NuBar->Draw("hist");
217  else FN_NuBar->Draw("hist same");
218  FNNuBarRatiosLeg->AddEntry(FN_NuBar, compname.c_str(), "l");
219 
220  TCanvas* NNbar = new TCanvas();
221  TLegend* LegNNbar = new TLegend(.6,.6,.95,.9);
222  TH1* NuNbarND = Ratio(NDNu[varIdx][compIdx][anaIdx], NDNuBar[varIdx][compIdx][anaIdx], ananame+" Nu/NuBar", compname+" Ancestor", "ND", StackColor[compIdx],1, LegNNbar );
223  TH1* NuNbarFD = Ratio(FDNu[varIdx][compIdx][anaIdx], FDNuBar[varIdx][compIdx][anaIdx], ananame+" Nu/NuBar", compname+" Ancestor", "FD",StackColor[compIdx],7, LegNNbar );
224  DrawHists(NuNbarND,NuNbarFD);
225  LegNNbar->Draw();
226  rootlogon();
227  gPad->Update();
228  NNbar->SaveAs(("plots/"+dirname+"/NuNuBar_ratio_"+varName+"_"+compname+"_"+ananame+"_ancestor.pdf").c_str());
229 
230  }
231  FNNuRatios->cd();
232  FNNuRatiosLeg->Draw();
233  gPad->Update();
234  FNNuRatios->SaveAs(("plots/All_FN_NuRatios_"+varName+"_"+ananame+".pdf").c_str());
235 
236  FNNuBarRatios->cd();
237  FNNuBarRatiosLeg->Draw();
238  gPad->Update();
239  FNNuBarRatios->SaveAs(("plots/All_FN_NuBarRatios_"+varName+"_"+ananame+".pdf").c_str());
240 
241  // FD and ND ancestors plots
242  CompNu->cd();
243  LegCompNu->Draw();
244  gPad->Update();
245  CompNu->SaveAs(("plots/"+dirname+"/FD_Nu_ancestor_"+varName+"_"+ananame+".pdf").c_str());
246 
247  CompNuBar->cd();
248  LegCompNuBar->Draw();
249  gPad->Update();
250  CompNuBar->SaveAs(("plots/"+dirname+"/FD_NuBar_ancestor_"+varName+"_"+ananame+".pdf").c_str());
251 
252  NDCompNu->cd();
253  NDLegCompNu->Draw();
254  gPad->Update();
255  NDCompNu->SaveAs(("plots/"+dirname+"/ND_Nu_ancestor_"+varName+"_"+ananame+".pdf").c_str());
256 
257  NDCompNuBar->cd();
258  NDLegCompNuBar->Draw();
259  gPad->Update();
260  NDCompNuBar->SaveAs(("plots/"+dirname+"/ND_NuBar_ancestor_"+varName+"_"+ananame+".pdf").c_str());
261 
262 
263 }
264 }
265 }
enum BeamMode kOrange
const int kNumVars
Definition: vars.h:14
enum BeamMode kRed
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
Antineutrinos-only.
Definition: IPrediction.h:50
(&#39; appearance&#39;)
Definition: IPrediction.h:18
(&#39;beam &#39;)
Definition: IPrediction.h:15
General interface to oscillation calculators.
Definition: StanTypedefs.h:22
std::string name
Definition: NuePlotLists.h:12
void plot_nuTree(std::string fname)
Definition: plot_nuTree.C:35
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
_NoOscillations< double > NoOscillations
Definition: IOscCalc.h:67
Charged-current interactions.
Definition: IPrediction.h:39
std::vector< std::string > flavors
const HistDef defs[kNumVars]
Definition: vars.h:15
std::vector< float > Spectrum
Definition: Constants.h:610
Represent the ratio between two spectra.
Definition: Ratio.h:8
Neutrinos-only.
Definition: IPrediction.h:49
OStream cout
Definition: OStream.cxx:6
(&#39; survival&#39;)
Definition: IPrediction.h:19
Interface to extrapolation procedures.
Definition: IExtrap.h:8
enum BeamMode kViolet
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
void rootlogon()
Definition: rootlogon.C:7
enum BeamMode kGreen
enum BeamMode string