7 #include "RooUnfoldResponse.h" 8 #include "RooUnfoldBayes.h" 26 TH2F* hDD = (TH2F*)sig->Project3D(
"yx");
28 for(
int i = 1;
i <= sig->GetXaxis()->GetNbins();
i++){
29 for(
int j = 1;
j <= sig->GetYaxis()->GetNbins();
j++){
35 hBkgd1D->Add(hSig1D,-1);
36 int sigbin = hSig1D->FindBin(0.85);
38 float tot_sig = hSig1D->Integral(sigbin,-1);
39 float tot_bkgd = hBkgd1D->Integral(sigbin,-1);
41 if(tot_bkgd <= 0) hDD->SetBinContent(i,j,0);
42 else hDD->SetBinContent(i,j,tot_sig/tot_bkgd);
47 hDD->GetZaxis()->SetTitle(
"Signal/Background (CVNe > 0.85)");
48 hDD->GetXaxis()->SetTitle(
"cos #theta_{e}");
49 hDD->GetYaxis()->SetTitle(
"Electron Energy, E_{e} (GeV)");
50 hDD->GetXaxis()->SetRangeUser(0.75,1);
51 hDD->GetYaxis()->SetRangeUser(1.0,10.0);
52 hDD->GetXaxis()->CenterTitle();
53 hDD->GetYaxis()->CenterTitle();
54 hDD->GetZaxis()->CenterTitle();
55 TCanvas* c88 =
new TCanvas(
"c88",
"c88");
56 c88->SetRightMargin(0.15);
57 c88->SetLeftMargin(0.15);
58 hDD->SetMarkerColor(kWhite);
59 gStyle->SetPaintTextFormat(
"4.2f");
60 hDD->Draw(
"COLZ TEXT");
62 c88->SaveAs(out.c_str());
72 bool shouldSave =
true)
76 TFile* fAnalysis =
new TFile((
sDir+
sAnalysis).c_str(),
"read");
85 std::vector<std::unique_ptr<ana::Spectrum>> ppfx_spects;
86 for(
int iuniv = 0; iuniv < 100; iuniv++){
88 sprintf(name,
"%s_%i_%s_%s",
"ppfx", iuniv,
89 pidName.c_str(), varName.c_str());
90 std::cout <<
"From " << fTemplates->GetName() <<
" loading " << name <<
"....";
119 std::vector<std::unique_ptr<ana::Spectrum>>genie_spects;
120 for(
int iuniv = 0; iuniv < 100; iuniv++){
122 sprintf(name,
"%s_%i_%s_%s",
"genie", iuniv,
123 pidName.c_str(), varName.c_str());
148 std::vector<TH3F*> systs_hists;
149 std::vector<TH3F*> systs_hists_up;
150 std::vector<TH3F*> systs_hists_down;
153 sprintf(name,
"calibpos_%s_%s_%s", pidName.c_str(),
154 varName.c_str(),
chns[0].
name.c_str());
156 (fTemplates->GetDirectory(name)))
159 sprintf(name,
"lightup_%s_%s_%s", pidName.c_str(),
160 varName.c_str(),
chns[0].
name.c_str());
162 (fTemplates->GetDirectory(name)))
164 systs_hists_up.push_back(hGenieUpper);
165 systs_hists_up.push_back(hPPFXUpper);
167 sprintf(name,
"calibneg_%s_%s_%s", pidName.c_str(),
168 varName.c_str(),
chns[0].
name.c_str());
170 (fTemplates->GetDirectory(name)))
173 sprintf(name,
"lightdown_%s_%s_%s", pidName.c_str(),
174 varName.c_str(),
chns[0].
name.c_str());
176 (fTemplates->GetDirectory(name)))
178 systs_hists_down.push_back(hGenieLower);
179 systs_hists_down.push_back(hPPFXLower);
181 sprintf(name,
"calibshape_%s_%s_%s", pidName.c_str(),
182 varName.c_str(),
chns[0].
name.c_str());
184 (fTemplates->GetDirectory(name)))->ToTH3(
pot));
185 sprintf(name,
"ckv_%s_%s_%s", pidName.c_str(),
186 varName.c_str(),
chns[0].
name.c_str());
188 (fTemplates->GetDirectory(name)))->ToTH3(
pot));
190 std::vector<TH3F*> systs_calibshape = {systs_hists[0]};
191 std::vector<TH3F*> systs_ckv = {systs_hists[1]};
192 std::vector<TH3F*> systs_genie_up = {hGenieUpper};
193 std::vector<TH3F*> systs_genie_down = {hGenieLower};
194 std::vector<TH3F*> systs_ppfx_up = {hPPFXUpper};
195 std::vector<TH3F*> systs_ppfx_down = {hPPFXLower};
196 std::vector<TH3F*> systs_light_up = {systs_hists_up[1]};
197 std::vector<TH3F*> systs_light_down = {systs_hists_down[1]};
198 std::vector<TH3F*> systs_calib_up = {systs_hists_up[0]};
199 std::vector<TH3F*> systs_calib_down = {systs_hists_down[0]};
201 std::vector<TH3F*> systs_onesided;
202 std::vector<TH3F*> systs_upper;
203 std::vector<TH3F*> systs_lower;
205 if(systName ==
"tot"){
206 systs_onesided = {systs_hists[0],systs_hists[1]};
207 systs_upper = {systs_hists_up[0], systs_hists_up[1],
208 hGenieUpper,hPPFXUpper};
209 systs_lower = {systs_hists_down[0], systs_hists_down[1],
210 hGenieLower,hPPFXLower};
212 else if(systName ==
"calibshape"){
213 systs_onesided = {systs_hists[0]};
217 else if(systName ==
"ckv"){
218 systs_onesided = {systs_hists[1]};
222 else if(systName ==
"genie"){
224 systs_upper = {hGenieUpper};
225 systs_lower = {hGenieLower};
227 else if(systName ==
"ppfx"){
229 systs_upper = {hPPFXUpper};
230 systs_lower = {hPPFXLower};
232 else if(systName ==
"calib"){
234 systs_upper = {systs_hists_down[0]};
235 systs_lower = {systs_hists_down[0]};
237 else if(systName ==
"light"){
239 systs_upper = {systs_hists_up[1]};
240 systs_lower = {systs_hists_down[1]};
250 std::vector<TH3F*> nominal_hists;
252 sprintf(name,
"nominal_%s_%s_%s", pidName.c_str(),
256 (fTemplates->GetDirectory(name)))
261 TH2F* hPlot2D = (TH2F*)nominal_hists[0]->Project3D(
"yx");
263 hPlot2D->SetTitle(
"");
264 hPlot2D->GetYaxis()->SetTitle(
"Electron Energy, E_{e} (GeV)");
265 hPlot2D->GetXaxis()->SetTitle(
"cos #theta_{e}");
266 hPlot2D->GetYaxis()->CenterTitle();
267 hPlot2D->GetYaxis()->SetTitleOffset(0.45);
268 hPlot2D->GetZaxis()->SetMaxDigits(3);
269 hPlot2D->GetXaxis()->CenterTitle();
270 hPlot2D->GetZaxis()->SetTitle(
"Events/8.09 #times 10^{20} POT");
271 hPlot2D->GetZaxis()->CenterTitle();
272 hPlot2D->GetXaxis()->SetRangeUser(0.75,1.00);
273 hPlot2D->GetYaxis()->SetRangeUser(1.0,6.0);
274 TH2F* hPlot2D_sig = (TH2F*)nominal_hists[2]->Project3D(
"yx");
275 TH2F* hPlot2D_sig_2 = (TH2F*)nominal_hists[7]->Project3D(
"yx");
277 hPlot2D_sig->GetXaxis()->SetRangeUser(0.75,1);
278 hPlot2D_sig->GetYaxis()->SetRangeUser(1.0,10.0);
280 hPlot2D_sig->SetTitle(
"");
281 hPlot2D_sig->GetYaxis()->SetTitle(
"Electron Energy, E_{e} (GeV)");
282 hPlot2D_sig->GetXaxis()->SetTitle(
"cos #theta_{e}");
283 hPlot2D_sig->GetYaxis()->CenterTitle();
284 hPlot2D_sig->GetYaxis()->SetTitleOffset(0.45);
285 hPlot2D_sig->GetZaxis()->SetMaxDigits(3);
286 hPlot2D_sig->GetXaxis()->CenterTitle();
287 hPlot2D_sig->GetZaxis()->SetTitle(
"Signal Events/8.09 #times 10^{20} POT");
288 hPlot2D_sig->GetZaxis()->CenterTitle();
293 TCanvas* c29 =
new TCanvas(
"c28",
"c29");
294 c29->SetRightMargin(0.15);
295 hPlot2D_sig->SetMarkerColor(kWhite);
296 gStyle->SetPaintTextFormat(
"4.0f");
297 hPlot2D_sig->Draw(
"COLZ TEXT");
299 c29->SaveAs(
"signal_mc_2d.png");
303 TH3F* hNominalSignalLike =
305 hNominalSignalLike->Add(nominal_hists[1]);
306 hNominalSignalLike->Add(nominal_hists[7]);
309 "s_b_ratio_signal_region.png");
312 TCanvas *c88 =
new TCanvas(
"c88",
"c88");
314 c88->SetRightMargin(0.15);
316 hPlot2D->Draw(
"COLZ");
318 c88->SaveAs(
"total_mc_2d.png");
322 TH2F* hPlot2D_a = (TH2F*)nominal_hists[0]->Project3D(
"yz");
324 hPlot2D_a->SetTitle(
"");
325 hPlot2D_a->GetYaxis()->SetTitle(
"Electron Energy, E_{e} (GeV)");
326 hPlot2D_a->GetXaxis()->SetTitle(
"Event CVNe");
327 hPlot2D_a->GetYaxis()->CenterTitle();
328 hPlot2D_a->GetYaxis()->SetTitleOffset(0.45);
329 hPlot2D_a->GetZaxis()->SetMaxDigits(3);
330 hPlot2D_a->GetXaxis()->CenterTitle();
331 hPlot2D_a->GetZaxis()->SetTitle(
"Events/8.09 #times 10^{20} POT");
332 hPlot2D_a->GetZaxis()->CenterTitle();
333 hPlot2D_a->GetXaxis()->SetRangeUser(0.,1.00);
334 hPlot2D_a->GetYaxis()->SetRangeUser(1.0,6.0);
336 TCanvas *c89 =
new TCanvas(
"c89",
"c89");
338 c89->SetRightMargin(0.15);
340 hPlot2D_a->Draw(
"COLZ");
342 c89->SaveAs(
"total_mc_2d_electron_pid.png");
346 std::vector<int> intvect = {};
354 std::vector<TH3F*> data_hists;
355 if (
dataName.find(
"ppfx") != std::string::npos) {
356 sprintf(name,
"mc_%s_%s_%s",
dataName.c_str(),
359 else if (
dataName.find(
"genie") != std::string::npos) {
360 sprintf(name,
"mc_%s_%s_%s",
dataName.c_str(),
364 sprintf(name,
"mc_%s_%s_%s",
dataName.c_str(),
368 GetDirectory(name)))->
372 std::vector<TH3F*> nominal_hists_3d;
374 sprintf(name,
"%s_%s_%s_%s",
"mc",
"nominal",
375 "analysis",
chns[
i].name.c_str());
387 std::vector<TH3F*> weighted_hists_3d =
389 systs_upper,systs_lower,nominal_hists_3d,
390 pidName,varName,out);
398 if(weighted_hists_3d.size() < 4){
407 sprintf(name,
"mc_%s_%s_signal_postfit",
dataName.c_str(), systName.c_str());
408 TH3F* hSignalOut = (TH3F*)weighted_hists_3d[2]->
Clone(name);
409 hSignalOut->SetName(name);
413 const int nx = hSignalOut->GetXaxis()->GetNbins();
414 const int ny = hSignalOut->GetYaxis()->GetNbins();
415 const int nz = hSignalOut->GetZaxis()->GetNbins();
418 TH1F* hHolder =
new TH1F(
"hHolder",
"",nx*ny*nz,0,nx*ny*nz);
420 for(
int i = 0;
i < hHolder->GetNbinsX(); ++
i){
421 const int nynz = ny*nz;
422 const int nmodnynz =
i%nynz;
423 const int ix =
i/nynz;
424 const int iy = nmodnynz/nz;
427 const double val = hSignalOut->GetBinContent(ix+1, iy+1, iz+1);
428 const double err = hSignalOut->GetBinError (ix+1, iy+1, iz+1);
430 hHolder->SetBinContent(
i+1,val);
431 hHolder->SetBinError(
i+1,err);
434 std::vector<TH1F*> vHolder;
435 for(
int j = 0;
j < (
int)weighted_hists_3d.size();
j++){
436 TH1F* hHold =
new TH1F(
ana::UniqueName().c_str(),
"",nx*ny*nz,0,nx*ny*nz);
438 for(
int i = 0;
i < hHold->GetNbinsX(); ++
i){
439 const int nynz = ny*nz;
440 const int nmodnynz =
i%nynz;
441 const int ix =
i/nynz;
442 const int iy = nmodnynz/nz;
445 const double val = weighted_hists_3d[
j]->GetBinContent(ix+1, iy+1, iz+1);
446 const double err = weighted_hists_3d[
j]->GetBinError (ix+1, iy+1, iz+1);
448 hHold->SetBinContent(
i+1,val);
449 hHold->SetBinError(
i+1,err);
451 vHolder.push_back(hHold);
452 sprintf(name,
"mc_%s_%s_%s_postfit",
dataName.c_str(),
454 vHolder[
j]->SetName(name);
457 if(shouldSave ==
true){
463 "CrossSection_TemplateFitResults_fakedata_development_stat.root";
464 TFile *
outf =
new TFile((outfile).c_str(),
"update");
466 sprintf(name,
"mc_%s_%s_signal_postfit",
dataName.c_str(), systName.c_str());
467 sSignalOut.
SaveTo(outf, name);
469 std::vector<Spectrum*> sHolder;
470 for(
int j = 0;
j < (
int)vHolder.size();
j++){
471 sprintf(name,
"mc_%s_%s_%s_postfit",
dataName.c_str(),
474 sHolder[
j]->SaveTo(outf, name);
const std::string sTemplates
void PlotSignalBackground(TH3F *sig, TH3F *tot, std::string out)
Cuts and Vars for the 2020 FD DiF Study.
const ana::Binning eelecbins
void makeXSecPlots_TemplateFit(std::string pidName="prongCVN", std::string varName="double", std::string dataName="nominal", std::string systName="tot", bool shouldSave=true)
const SelDef chns[knumchns]
const ana::Binning costhetabins
Representation of a spectrum in any variable, with associated POT.
const ana::Binning pidbins
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
void SaveTo(TDirectory *dir, const std::string &name) const
std::vector< float > Spectrum
std::vector< TH3F * > BinByBinTemplateFit(TH3F *data, std::vector< TH3F * > templates, std::vector< TH3F * > systs_hists, std::vector< TH3F * > systs_hists_up, std::vector< TH3F * > systs_hists_down, std::vector< TH3F * > analysis_templates, std::string pidName, std::string varName, std::string dataName)
const std::string sAnalysis
const Spectrum * UpperSigma(BandOptions opt=kBandFromNominal) const
const Spectrum * LowerSigma(BandOptions opt=kBandFromNominal) const
const double kAna2018FHCPOT
TH3 * ToTH3(const Spectrum &s, double exposure, ana::EExposureType expotype, const Binning &binsx, const Binning &binsy, const Binning &binsz, ana::EBinType bintype)
Same as ToTH2, but with 3 dimensions.
std::string UniqueName()
Return a different string each time, for creating histograms.