ppfx_smooth_weights_make.C
Go to the documentation of this file.
2 
3 #include "CAFAna/Core/Spectrum.h"
6 #include "CAFAna/Core/HistAxis.h"
7 #include "CAFAna/Vars/Vars.h"
10 #include "CAFAna/Core/Binning.h"
11 #include "CAFAna/Core/Cut.h"
12 #include "CAFAna/Core/SystShifts.h"
13 #include "CAFAna/Cuts/Cuts.h"
14 #include "CAFAna/Cuts/SpillCuts.h"
16 #include "CAFAna/Cuts/TruthCuts.h"
17 #include "CAFAna/Cuts/NueCutsSecondAna.h"
18 #include "CAFAna/Cuts/TimingCuts.h"
19 #include "CAFAna/Systs/Systs.h"
20 
21 #include "TStopwatch.h"
22 #include "CAFAna/Core/Utilities.h"
23 #include "Utilities/func/MathUtil.h"
24 
25 #include "TFile.h"
26 #include "TRandom.h"
27 #include "TCanvas.h"
28 
29 #include <fstream>
30 #include <iostream>
31 #include <cmath>
32 
33 
34 using namespace ana;
35 
37 {
38 
39  std::vector <std::string> comps = {
40  "numu",
41  "nue",
42  "anumu",
43  "anue",
44  };
45 
46  std::vector<std::string> directories;
47  for(int i = 0; i < (int)comps.size(); i++){
48  directories.push_back("unweighted_FHC_ND_"+comps[i]);
49  directories.push_back("unweighted_FHC_FD_"+comps[i]);
50  directories.push_back("ppfx_cv_FHC_ND_"+comps[i]);
51  directories.push_back("ppfx_cv_FHC_FD_"+comps[i]);
52  directories.push_back("unweighted_RHC_ND_"+comps[i]);
53  directories.push_back("unweighted_RHC_FD_"+comps[i]);
54  directories.push_back("ppfx_cv_RHC_ND_"+comps[i]);
55  directories.push_back("ppfx_cv_RHC_FD_"+comps[i]);
56 
57  for(int j = 0; j < 100; j++){
58  int dig0 = j % 10;
59  int dig1 = (j/10) % 10;
60  directories.push_back("ppfx_univ"+std::to_string(dig1)+std::to_string(dig0)+"_FHC_ND_"+comps[i]);
61  directories.push_back("ppfx_univ"+std::to_string(dig1)+std::to_string(dig0)+"_FHC_FD_"+comps[i]);
62  directories.push_back("ppfx_univ"+std::to_string(dig1)+std::to_string(dig0)+"_RHC_ND_"+comps[i]);
63  directories.push_back("ppfx_univ"+std::to_string(dig1)+std::to_string(dig0)+"_RHC_FD_"+comps[i]);
64  }
65  }
66 
67  std::string filename = "ppfx_weights.root";
68  TFile* fread = new TFile(filename.c_str(), "READ");
69 
70  std::map<std::string, TH1*> hists;
71 
72  std::for_each(directories.begin(), directories.end(),
73  [&](std::string dir){
74  auto spect = Spectrum::LoadFrom(fread, dir.c_str());
75  TH1* h = (TH1*)spect->ToTH1(spect->POT());
76  h->Rebin(4);
77  hists.insert(std::make_pair(dir, h));
78 
79  });
80 
81  std::map<std::string, TH1*> hist_ratios_smooth;
82  std::map<std::string, TH1*> hist_ratios_unsmooth;
83  TFile* fout = new TFile("ppfx_smooth_multiverse_weights_2018.root", "recreate");
84  fout->cd();
85  std::for_each(hists.begin(), hists.end(),
86  [&](std::pair<std::string, TH1*> hist){
87  if(hist.first.find("unweighted") != string::npos)
88  return;
89 
91  std::string flav;
93 
94  if(hist.first.find("FHC") != string::npos)
95  horn = "FHC";
96  if(hist.first.find("RHC") != string::npos)
97  horn = "RHC";
98 
99  if(hist.first.find("ND") != string::npos)
100  det = "ND";
101  if(hist.first.find("FD") != string::npos)
102  det = "FD";
103 
104  if(hist.first.find("numu") != string::npos)
105  flav = "numu";
106  if(hist.first.find("nue") != string::npos)
107  flav = "nue";
108  if(hist.first.find("anumu") != string::npos)
109  flav = "anumu";
110  if(hist.first.find("anue") != string::npos)
111  flav = "anue";
112 
113  TH1* h = (TH1*)hist.second->Clone("");
114  std::string denom = "unweighted_"+horn+"_"+det+"_"+flav;
115  h->Divide(hists[denom]);
116  TH1* h_unsmooth = (TH1*) h->Clone("");
117  h_unsmooth->GetYaxis()->SetRangeUser(0.6, 1.2);
118  hist_ratios_unsmooth.insert(std::make_pair(hist.first, h_unsmooth));
119 
120  h->Smooth(1);
121  h->GetYaxis()->SetRangeUser(0.6, 1.2);
122  hist_ratios_smooth.insert(std::make_pair(hist.first, h));
123 
124  std::string plusname = hist.first+"_plussigma";
125  std::string minusname = hist.first+"_minussigma";
126  h->SetName(plusname.c_str());
127  h->Write();
128 
129  TH1* hminus = (TH1*)h->Clone(minusname.c_str());
130  hminus->Write();
131  });
132 
133  // // plotting
134  // // fd
135  // for(int i = 0; i < (int)comps.size(); i++){
136  // TCanvas *c1 = new TCanvas("c1", "c1", 200, 10, 1200, 600);
137  // c1->Divide(2,1);
138  // for(int j = 0; j < 10; j++){
139  // int dig0 = j % 10;
140  // int dig1 = (j/10) % 10;
141  // TH1* h_univ_fd_smooth = (TH1*) hist_ratios_smooth["ppfx_univ"+std::to_string(dig1)+std::to_string(dig0)+"_RHC_FD_"+comps[i]]->Clone("");
142  // TH1* h_univ_fd_unsmooth = (TH1*) hist_ratios_unsmooth["ppfx_univ"+std::to_string(dig1)+std::to_string(dig0)+"_RHC_FD_"+comps[i]]->Clone("");
143  //
144  // c1->cd(1);
145  // h_univ_fd_smooth->Draw("hist same");
146  // h_univ_fd_smooth->SetLineColor(dig0);
147  //
148  // c1->cd(2);
149  // h_univ_fd_unsmooth->Draw("hist same");
150  // h_univ_fd_unsmooth->SetLineColor(dig0);
151  // }
152  // TH1* h_cv_fd_smooth = (TH1*) hist_ratios_smooth["ppfx_cv_RHC_FD_"+comps[i]]->Clone("");
153  // TH1* h_cv_fd_unsmooth = (TH1*) hist_ratios_unsmooth["ppfx_cv_RHC_FD_"+comps[i]]->Clone("");
154  //
155  // c1->cd(1);
156  // h_cv_fd_smooth->Draw("hist same");
157  // h_cv_fd_smooth->SetLineColor(36);
158  //
159  // c1->cd(2);
160  // h_cv_fd_unsmooth->Draw("hist same");
161  // h_cv_fd_unsmooth->SetLineColor(36);
162  // c1->Update();
163  // std::string picname = "smooth_plots/fd_RHC_"+comps[i]+".pdf";
164  // std::string picname2 = "fd_"+comps[i]+".root";
165  // c1->Print(picname.c_str());
166  // // c1->Print(picname2.c_str());
167  // }
168  //
169  // // nd
170  // for(int i = 0; i < (int)comps.size(); i++){
171  // TCanvas *c1 = new TCanvas("c1", "c1", 200, 10, 1200, 600);
172  // c1->Divide(2,1);
173  // for(int j = 0; j < 10; j++){
174  // int dig0 = j % 10;
175  // int dig1 = (j/10) % 10;
176  // TH1* h_univ_nd_smooth = (TH1*) hist_ratios_smooth["ppfx_univ"+std::to_string(dig1)+std::to_string(dig0)+"_RHC_ND_"+comps[i]]->Clone("");
177  // TH1* h_univ_nd_unsmooth = (TH1*) hist_ratios_unsmooth["ppfx_univ"+std::to_string(dig1)+std::to_string(dig0)+"_RHC_ND_"+comps[i]]->Clone("");
178  //
179  // c1->cd(1);
180  // h_univ_nd_smooth->Draw("hist same");
181  // h_univ_nd_smooth->SetLineColor(dig0);
182  //
183  // c1->cd(2);
184  // h_univ_nd_unsmooth->Draw("hist same");
185  // h_univ_nd_unsmooth->SetLineColor(dig0);
186  // }
187  // TH1* h_cv_nd_smooth = (TH1*) hist_ratios_smooth["ppfx_cv_RHC_ND_"+comps[i]]->Clone("");
188  // TH1* h_cv_nd_unsmooth = (TH1*) hist_ratios_unsmooth["ppfx_cv_RHC_ND_"+comps[i]]->Clone("");
189  //
190  // c1->cd(1);
191  // h_cv_nd_smooth->Draw("hist same");
192  // h_cv_nd_smooth->SetLineColor(36);
193  //
194  // c1->cd(2);
195  // h_cv_nd_unsmooth->Draw("hist same");
196  // h_cv_nd_unsmooth->SetLineColor(36);
197  // c1->Update();
198  // std::string picname = "smooth_plots/nd_RHC_"+comps[i]+".pdf";
199  // std::string picname2 = "nd_"+comps[i]+".root";
200  // c1->Print(picname.c_str());
201  // // c1->Print(picname2.c_str());
202  // }
203  //
204  // // fn 10
205  // for(int i = 0; i < (int)comps.size(); i++){
206  // TCanvas *c1 = new TCanvas("c1", "c1", 200, 10, 1200, 600);
207  // c1->Divide(2,1);
208  // TH1* h_cv_nd_smooth = (TH1*) hist_ratios_smooth["ppfx_cv_RHC_ND_"+comps[i]]->Clone("");
209  // TH1* h_cv_fd_smooth = (TH1*) hist_ratios_smooth["ppfx_cv_RHC_FD_"+comps[i]]->Clone("");
210  // h_cv_fd_smooth->Divide(h_cv_nd_smooth);
211  //
212  // TH1* h_cv_nd_unsmooth = (TH1*) hist_ratios_unsmooth["ppfx_cv_RHC_ND_"+comps[i]]->Clone("");
213  // TH1* h_cv_fd_unsmooth = (TH1*) hist_ratios_unsmooth["ppfx_cv_RHC_FD_"+comps[i]]->Clone("");
214  // h_cv_fd_unsmooth->Divide(h_cv_nd_unsmooth);
215  //
216  // for(int j = 0; j < 10; j++){
217  // int dig0 = j % 10;
218  // int dig1 = (j/10) % 10;
219  // TH1* h_univ_nd_smooth = (TH1*) hist_ratios_smooth["ppfx_univ"+std::to_string(dig1)+std::to_string(dig0)+"_RHC_ND_"+comps[i]]->Clone("");
220  // TH1* h_univ_fd_smooth = (TH1*) hist_ratios_smooth["ppfx_univ"+std::to_string(dig1)+std::to_string(dig0)+"_RHC_FD_"+comps[i]]->Clone("");
221  // h_univ_fd_smooth->Divide(h_univ_nd_smooth);
222  // h_univ_fd_smooth->Divide(h_cv_fd_smooth);
223  //
224  // TH1* h_univ_nd_unsmooth = (TH1*) hist_ratios_unsmooth["ppfx_univ"+std::to_string(dig1)+std::to_string(dig0)+"_RHC_ND_"+comps[i]]->Clone("");
225  // TH1* h_univ_fd_unsmooth = (TH1*) hist_ratios_unsmooth["ppfx_univ"+std::to_string(dig1)+std::to_string(dig0)+"_RHC_FD_"+comps[i]]->Clone("");
226  // h_univ_fd_unsmooth->Divide(h_univ_nd_unsmooth);
227  // h_univ_fd_unsmooth->Divide(h_cv_fd_unsmooth);
228  //
229  // c1->cd(1);
230  // h_univ_fd_smooth->GetYaxis()->SetRangeUser(0.8, 1.2);
231  // h_univ_fd_smooth->Draw("hist same");
232  // h_univ_fd_smooth->SetLineColor(dig0);
233  //
234  // c1->cd(2);
235  // h_univ_fd_unsmooth->GetYaxis()->SetRangeUser(0.8, 1.2);
236  // h_univ_fd_unsmooth->Draw("hist same");
237  // h_univ_fd_unsmooth->SetLineColor(dig0);
238  // }
239  // c1->Update();
240  // std::string picname = "smooth_plots/fn_double_ratio_10_"+comps[i]+".pdf";
241  // std::string picname2 = "fn_double_ratio_10_"+comps[i]+".root";
242  // c1->Print(picname.c_str());
243  // // c1->Print(picname2.c_str());
244  //
245  // }
246 
247 }
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
string filename
Definition: shutoffs.py:106
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
Definition: DataMCLoad.C:336
TString hists[nhists]
Definition: bdt_com.C:3
void ppfx_smooth_weights_make()
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:535
const double j
Definition: BetheBloch.cxx:29
std::vector< std::string > comps
TDirectory * dir
Definition: macro.C:5
TH1F * hminus
Definition: Xsec_final.C:76
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
enum BeamMode string