Classes | Functions | Variables
plot_uncertainty.C File Reference
#include <stdio.h>
#include <stdlib.h>
#include <string>
#include "Utilities/rootlogon.C"
#include "CAFAna/Cuts/SpillCuts.h"
#include "CAFAna/Core/Utilities.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Cuts/TimingCuts.h"
#include "CAFAna/Analysis/Plots.h"
#include "CAFAna/Analysis/Style.h"
#include "CAFAna/Core/LoadFromFile.h"
#include "CAFAna/Core/Ratio.h"
#include "CAFAna/XSec/GenieMultiverseSyst.h"
#include "TH1D.h"
#include "TCanvas.h"
#include "TGraph.h"
#include "TMultiGraph.h"
#include "TLine.h"
#include "TLegend.h"
#include "TPad.h"
#include "TGaxis.h"

Go to the source code of this file.

Classes

struct  FileUp
 
struct  DateUp
 

Functions

void NDSimulation ()
 
void PlotWithSystErrorBandWidth (TH1 *&nom, std::vector< TH1 * > &ups, std::vector< TH1 * > &dns, int col, int errCol, float headroom, bool newaxis)
 
void PlotDeltaSigmaWithSigma (TH1 *&nom, TH1 *&nomsel, std::vector< TH1 * > &ups, std::vector< TH1 * > &dns, std::vector< TH1 * > &upssel, std::vector< TH1 * > &dnssel, TH1 *&heff, std::vector< TH1 * > &effups, std::vector< TH1 * > &effdns, TString outname1, TString outname2, TString outname3, TString outname4, TString outname5, TString outname6, TString outname7, TString outname8, TString outname9, Bool_t printme=false)
 
void GetFluxError (const Spectrum *nom, std::vector< ana::Spectrum * > univs, std::vector< ana::Spectrum * > &supdown)
 
void plot_uncertainty ()
 

Variables

const double pot = 8.09e20
 
TString outDir2 = "/nova/app/users/acedeno/tag_releasesS18-06-14/NDAna/ncpi0_semi_inc_png_cvn/NewDataSetSyst"
 
const int File =13
 
const FileUp CurrentFile [File]
 
const int Date = 1
 
const DateUp CurrentDate [Date]
 

Function Documentation

void GetFluxError ( const Spectrum *  nom,
std::vector< ana::Spectrum * >  univs,
std::vector< ana::Spectrum * > &  supdown 
)

Definition at line 428 of file plot_uncertainty.C.

References make_syst_table_plots::ibin, ana::Spectrum::Livetime(), extractScale::mean, pot, cet::pow(), std::sqrt(), and ana::Spectrum::ToTH1().

Referenced by plot_uncertainty().

430 {
431  TH1* hnom = nom->ToTH1(pot);
432  TH1* hup = (TH1*)hnom->Clone("up");
433  TH1* hdown= (TH1*)hnom->Clone("down");
434 
435  for(int ibin=1; ibin <= hnom->GetNbinsX(); ibin++){
436 
437  //mean:
438  double mean = hnom->GetBinContent(ibin);
439  if(mean<=0)
440  continue;
441  //error:
442  double err = 0.0;
443  for(unsigned int iuniv=0; iuniv < univs.size(); iuniv++){
444  err += pow(univs[iuniv]->ToTH1(pot)->GetBinContent(ibin) - mean, 2);
445  }
446  err /= double(univs.size());
447  err = sqrt(err);
448 
449  hup->SetBinContent(ibin, mean+err);
450  hdown->SetBinContent(ibin, mean-err);
451  }
452  supdown.push_back(new Spectrum(hup, pot, nom->Livetime()) );
453  supdown.push_back(new Spectrum(hdown, pot, nom->Livetime()) );
454 }
ana::Spectrum::ToTH1
TH1D * ToTH1(double exposure, Color_t col=kBlack, Style_t style=kSolid, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Histogram made from this Spectrum, scaled to some exposure.
Definition: Spectrum.cxx:148
std::sqrt
T sqrt(T number)
Definition: d0nt_math.hpp:156
cet::pow
constexpr T pow(T x)
Definition: pow.h:75
cmf::Spectrum
std::vector< double > Spectrum
Definition: Constants.h:743
make_syst_table_plots.ibin
ibin
Definition: make_syst_table_plots.py:246
pot
const double pot
Definition: plot_uncertainty.C:37
ana::Spectrum::Livetime
double Livetime() const
Seconds. For informational purposes only. No calculations use this.
Definition: Spectrum.h:230
void NDSimulation ( )

Definition at line 41 of file plot_uncertainty.C.

References prelim.

Referenced by plot_uncertainty(), and PlotDeltaSigmaWithSigma().

42 {
43  TLatex* prelim = new TLatex(.9, .95, "NOvA Simulation");
44  prelim->SetTextColor(kGray+1);
45  prelim->SetNDC();
46  prelim->SetTextSize(2/30.);
47  prelim->SetTextAlign(32);
48  TGaxis::SetMaxDigits(3);
49  prelim->Draw();
50 }
prelim
TLatex * prelim
Definition: Xsec_final.C:133
void plot_uncertainty ( )

***********************************************************************************************************************************************************************************************////

comment this if just caliip or down is needed for plots

comment this if just caliip or down is needed for plots

comment this if just caliip or down is needed for plots

comment this if just caliip or down is needed for plots

***********************************************************************************************************************************************************************************************////

comment this if just caliip or down is needed for plots

comment this if just caliip or down is needed for plots

comment this if just caliip or down is needed for plots

comment this if just caliip or down is needed for plots

***********************************************************************************************************************************************************************************************////

comment this if just caliip or down is needed for plots

comment this if just caliip or down is needed for plots

comment this if just caliip or down is needed for plots

comment this if just caliip or down is needed for plots

***********************************************************************************************************************************************************************************************////

comment this if just caliip or down is needed for plots

comment this if just caliip or down is needed for plots

comment this if just caliip or down is needed for plots

comment this if just caliip or down is needed for plots

TH1* truesignal_ckv = ckv[2]->ToTH1(pot); ///comment this if just caliip or down is needed for plots

***********************************************************************************************************************************************************************************************////

comment this if just caliip or down is needed for plots

comment this if just caliip or down is needed for plots

comment this if just caliip or down is needed for plots

comment this if just caliip or down is needed for plots

comment this if just caliip or down is needed for plots

comment this if just caliip or down is needed for plots

*************************************************************************************************************************************************************************************///// /Total PreSel : kVtxInFiducial && kContainment && k2Prongs && kPi0RemIDND && kNpCVN1gammaID;

comment this if just caliip or down is needed for plots

comment this if just caliip or down is needed for plots

comment this if just caliip or down is needed for plots

comment this if just caliip or down is needed for plots

comment this if just caliip or down is needed for plots

comment this if just caliip or down is needed for plots

comment lightup[3]->ToTH1(pot), ckv[3]->ToTH1(pot), calib_shape[3]->ToTH1(pot), fluxup[3]->ToTH1(pot), xsecup[3]->ToTH1(pot) if plotting just calib.

comment lightup[3]->ToTH1(pot), ckv[3]->ToTH1(pot), calib_shape[3]->ToTH1(pot), fluxup[3]->ToTH1(pot), xsecup[3]->ToTH1(pot) if plotting just calib.

comment this if just caliip or down is needed for plots

comment this if just caliip or down is needed for plots

comment this if just caliip or down is needed for plots

comment this if just caliip or down is needed for plots

comment this if just caliip or down is needed for plots

comment this if just caliip or down is needed for plots

comment this if just caliip or down is needed for plots

comment this if just caliip or down is needed for plots

comment this if just caliip or down is needed for plots

comment this if just caliip or down is needed for plots

comment this if just caliip or down is needed for plots

comment this if just caliip or down is needed for plots

comment this if just caliip or down is needed for plots

comment this if just caliip or down is needed for plots

comment this if just caliip or down is needed for plots

comment this if just caliip or down is needed for plots

comment this if just caliip or down is needed for plots

comment this if just caliip or down is needed for plots

Definition at line 457 of file plot_uncertainty.C.

References om::cout, cut_labels, DateUp::date, allTimeWatchdog::endl, MakeMiniprodValidationCuts::f, fin, plot_validation_datamc::fname, GetFluxError(), ana::GenieMultiverseSpectra::kBandFromNominal, kRed, ana::Spectrum::LoadFrom(), ana::GenieMultiverseSpectra::LowerSigma(), getGoodRuns4SAM::n, NDSimulation(), outDir2, PlotDeltaSigmaWithSigma(), PlotWithSystErrorBandWidth(), pot, PandAna.Demos.pi0_spectra::pur, runNovaSAM::release, std::sqrt(), string, art::to_string(), and ana::GenieMultiverseSpectra::UpperSigma().

458 {
459 
460  //Give input file here:
461  //Fulloutput
462  std::string fname = "/nova/app/users/acedeno/tag_releasesS18-06-14/NDAna/ncpi0_semi_inc_png_cvn/PNG1OPTISyst_newpngcvn_Nominal.root";
463 
464  TFile fin(fname.c_str());
465 
466  const int nvar = 1;
467  //Define variable and different cuts to be used
468  std::vector<std::string> taxes_labels = {"BDTG"};
469  std::vector<std::string> cut_labels = {"totmc","NSel_Cuts","Nsig_fid_only","Nbkg_Cuts"};
470 
471  // loop over all the vars
472  for(int ivar = 0; ivar < nvar; ivar++){
473 
474  // vectors for different cut levels
475  std::vector<ana::Spectrum*> calibup, calibdown;
476  std::vector<ana::Spectrum*> lightup, lightdown;
477  std::vector<ana::Spectrum*> fluxup, fluxdown;
478  std::vector<ana::Spectrum*> ckv, calib_shape;
479  std::vector<ana::Spectrum*> nominal;
480  std::vector<const ana::Spectrum*> xsecup, xsecdown;
481 
482  //for(int icut = 0; icut < 4; icut++){
483  for(int icut = 0; icut < 4; icut++){
484 
485  nominal.push_back(LoadFromFile<Spectrum>(fname, "nominal/"+cut_labels[icut]+"_"+taxes_labels[ivar]).release());
486 
487  //calibup.push_back(LoadFromFile<Spectrum>(fname, "calibpos/"+cut_labels[icut]+"_"+taxes_labels[ivar]).release());
488  //calibdown.push_back(LoadFromFile<Spectrum>(fname, "calibneg/"+cut_labels[icut]+"_"+taxes_labels[ivar]).release());
489 
490  //lightup.push_back(LoadFromFile<Spectrum>(fname, "lightup/"+cut_labels[icut]+"_"+taxes_labels[ivar]).release()); ///comment this if just caliip or down is needed for plots
491  //lightdown.push_back(LoadFromFile<Spectrum>(fname, "lightdown/"+cut_labels[icut]+"_"+taxes_labels[ivar]).release()); ///comment this if just caliip or down is needed for plots
492 
493  //ckv.push_back(LoadFromFile<Spectrum>(fname, "ckv/"+cut_labels[icut]+"_"+taxes_labels[ivar]).release()); ///comment this if just caliip or down is needed for plots
494  //calib_shape.push_back(LoadFromFile<Spectrum>(fname, "calibshape/"+cut_labels[icut]+"_"+taxes_labels[ivar]).release()); ///comment this if just caliip or down is needed for plots
495 
496  // get the univ spectr
497  std::vector<ana::Spectrum*> flux_univs;
498  std::vector<std::unique_ptr<ana::Spectrum> > xsec_univs;
499  for(int iuniv = 0; iuniv < 100; iuniv++){
500 
501  xsec_univs.push_back(Spectrum::LoadFrom((TDirectory*)fin.Get(("xsec_"+cut_labels[icut]+
502  "_"+taxes_labels[ivar]+
503  "/spec"+ std::to_string(iuniv)).c_str() )));
504 
505  flux_univs.push_back(LoadFromFile<Spectrum>(fname, "flux_"+cut_labels[icut]+"_"+
506  taxes_labels[ivar]+"/spec"+ std::to_string(iuniv)).release());
507  }// end loop over universes
508 
509  // for flux error
510  std::vector<ana::Spectrum*> fupdown;
511  GetFluxError( nominal[icut], flux_univs, fupdown);
512  fluxup.push_back(fupdown[0]);
513  fluxdown.push_back(fupdown[1]);
514 
515  // for xsec errors
516  ana::MultiverseSpectra verses(xsec_univs);
517  xsecup.push_back(verses.UpperSigma(MultiverseSpectra::kBandFromNominal));
518  xsecdown.push_back(verses.LowerSigma(MultiverseSpectra::kBandFromNominal));
519 
520  }// end loop over cuts
521 
522  ///***********************************************************************************************************************************************************************************************////
523 
524  //For background kContainment && k2Prongs && kPi0RemIDND && kNpCVN1gammaID && !kNCurrent && !kPi02; (Not Signal Background) (NBkg)
525 
526  TH1* Nbkg_Cuts_xsecdown = xsecdown[3]->ToTH1(pot); ///comment this if just caliip or down is needed for plots
527  TH1* Nbkg_Cuts_xsecup = xsecup[3]->ToTH1(pot); ///comment this if just caliip or down is needed for plots
528 
529  TH1* Nbkg_Cuts_fluxdown = fluxdown[3]->ToTH1(pot); ///comment this if just caliip or down is needed for plots
530  TH1* Nbkg_Cuts_fluxup = fluxup[3]->ToTH1(pot); ///comment this if just caliip or down is needed for plots
531 
532  // TH1* Nbkg_Cuts_calib_shape = calib_shape[3]->ToTH1(pot); ///comment this if just caliip or down is needed for plots
533 
534  // TH1* Nbkg_Cuts_ckv = ckv[3]->ToTH1(pot); ///comment this if just caliip or down is needed for plots
535 
536  // TH1* Nbkg_Cuts_lightup = lightup[3]->ToTH1(pot); ///comment this if just caliip or down is needed for plots
537  // TH1* Nbkg_Cuts_lightdown = lightdown[3]->ToTH1(pot); ///comment this if just caliip or down is needed for plots
538 
539  // TH1* Nbkg_Cuts_calibdown = calibdown[3]->ToTH1(pot);
540  // TH1* Nbkg_Cuts_calibup = calibup[3]->ToTH1(pot);
541 
542 
543  // bkg signal distributions
544  std::vector<TH1*> vector_up_Nbkg_Cuts = {xsecup[3]->ToTH1(pot),fluxup[3]->ToTH1(pot)}; // {fluxup[3]->ToTH1(pot), calibup[3]->ToTH1(pot), lightup[3]->ToTH1(pot), ckv[3]->ToTH1(pot), calib_shape[3]->ToTH1(pot)}; ///comment lightup[3]->ToTH1(pot), ckv[3]->ToTH1(pot), calib_shape[3]->ToTH1(pot), fluxup[3]->ToTH1(pot), xsecup[3]->ToTH1(pot) if plotting just calib.
545  std::vector<TH1*> vector_dwn_Nbkg_Cuts = {xsecdown[3]->ToTH1(pot),fluxdown[3]->ToTH1(pot)};// {fluxdown[3]->ToTH1(pot), xsecdown[3]->ToTH1(pot), calibdown[3]->ToTH1(pot), lightdown[3]->ToTH1(pot),ckv[3]->ToTH1(pot), calib_shape[3]->ToTH1(pot)}; ///comment lightup[3]->ToTH1(pot), ckv[3]->ToTH1(pot), calib_shape[3]->ToTH1(pot), fluxup[3]->ToTH1(pot), xsecup[3]->ToTH1(pot) if plotting just calib.
546 
547 
548  TH1* hnom_Nbkg_Cuts = nominal[3]->ToTH1(pot);
549 
550  std::cout << "*************Bins nominal*************: " << hnom_Nbkg_Cuts->GetXaxis()->GetNbins() << std::endl;
551 
552  hnom_Nbkg_Cuts->SetTitle(";BDTG; Background Events/8.09 #times 10^{20}POT");
553  // hnom_bkg_nopresel->SetTitleSize(";0.05; 0.04");
554 
555  gROOT->ProcessLine(".! mkdir -p plots");
556 
557  TCanvas* Prong_CVN_Gamma_ID_Nbkg_Cuts = new TCanvas("Prong_CVN_Gamma_ID_Nbkg_Cuts", "Prong_CVN_Gamma_ID_Nbkg_Cuts");
558  Prong_CVN_Gamma_ID_Nbkg_Cuts->SetLeftMargin(0.15);
559 
560  PlotWithSystErrorBandWidth(hnom_Nbkg_Cuts,
561  vector_up_Nbkg_Cuts,
562  vector_dwn_Nbkg_Cuts, -1, -1, 1.3f, true);
563 
564  NDSimulation();
565 
566  Prong_CVN_Gamma_ID_Nbkg_Cuts->SaveAs(outDir2+CurrentFile[0].name.c_str()+CurrentDate[0].date.c_str()+".pdf");
567  //Prong_CVN_Gamma_ID_Nbkg_Cuts->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/Updatedfid_containment_all_daisysknobs/Prong_CVN_Gamma_ID_Nbkg_Cuts_0_3_high_Updatedfid_containment_all_daisysknobs.pdf");
568  //Prong_CVN_Gamma_ID_Nbkg_Cuts->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/Prong_CVN_Gamma_ID_Nbkg_Cuts_0_3_high_Combined_OneKnob1_out_3_MaNCRES.pdf");
569  //Prong_CVN_Gamma_ID_Nbkg_Cuts->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/Prong_CVN_Gamma_ID_Nbkg_Cuts_0_3_high.pdf");
570  //Prong_CVN_Gamma_ID_Nbkg_Cuts->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/Calibration/Prong_CVN_Gamma_ID_Nbkg_Cuts_calib_0_3_high.pdf");
571  // Prong_CVN_Gamma_ID_Nbkg_Cuts->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/Calibration_shape/Prong_CVN_Gamma_ID_Nbkg_Cuts__calib_shape_0_3_high.pdf");
572  //Prong_CVN_Gamma_ID_Nbkg_Cuts->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/Cherenkov/Prong_CVN_Gamma_ID_Nbkg_Cuts_Cherenkov_0_3_high.pdf");
573  // Prong_CVN_Gamma_ID_Nbkg_Cuts->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/GENIE/Prong_CVN_Gamma_ID_Nbkg_Cuts_GENIE_0_3_high.pdf");
574  // Prong_CVN_Gamma_ID_Nbkg_Cuts->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/light_level/Prong_CVN_Gamma_ID_Nbkg_Cuts_light_level_0_3_high.pdf");
575  // Prong_CVN_Gamma_ID_Nbkg_Cuts->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/ppfx/Prong_CVN_Gamma_ID_Nbkg_Cuts_ppfx_0_3_high.pdf");
576  ///***********************************************************************************************************************************************************************************************////
577 
578  //For Total Preselection cut: kVtxInFiducial && kContainment && k2Prongs && kPi0RemIDND && kNpCVN1gammaID (totmc)
579 
580  TH1* NSel_Cuts_xsecdown = xsecdown[0]->ToTH1(pot); ///comment this if just caliip or down is needed for plots
581  TH1* NSel_Cuts_xsecup = xsecup[0]->ToTH1(pot); ///comment this if just caliip or down is needed for plots
582 
583  TH1* NSel_Cuts_fluxdown = fluxdown[0]->ToTH1(pot); ///comment this if just caliip or down is needed for plots
584  TH1* NSel_Cuts_fluxup = fluxup[0]->ToTH1(pot); ///comment this if just caliip or down is needed for plots
585 
586  // TH1* NSel_Cuts_calib_shape = calib_shape[0]->ToTH1(pot); ///comment this if just caliip or down is needed for plots
587 
588  // TH1* NSel_Cuts_ckv = ckv[0]->ToTH1(pot); ///comment this if just caliip or down is needed for plots
589 
590  // TH1* NSel_Cuts_lightup = lightup[0]->ToTH1(pot); ///comment this if just caliip or down is needed for plots
591  // TH1* NSel_Cuts_lightdown = lightdown[0]->ToTH1(pot); ///comment this if just caliip or down is needed for plots
592 
593  // TH1* NSel_Cuts_calibdown = calibdown[0]->ToTH1(pot);
594  // TH1* NSel_Cuts_calibup = calibup[0]->ToTH1(pot);
595 
596 
597  // sig signal distributions
598  std::vector<TH1*> vector_up_NSel_Cuts = {xsecup[0]->ToTH1(pot),fluxup[0]->ToTH1(pot)};//{fluxup[0]->ToTH1(pot), xsecup[0]->ToTH1(pot),calibup[0]->ToTH1(pot), lightup[0]->ToTH1(pot), ckv[0]->ToTH1(pot), calib_shape[0]->ToTH1(pot)}; ///comment lightup[3]->ToTH1(pot), ckv[3]->ToTH1(pot), calib_shape[3]->ToTH1(pot), fluxup[3]->ToTH1(pot), xsecup[3]->ToTH1(pot) if plotting just calib.
599  std::vector<TH1*> vector_dwn_NSel_Cuts = {xsecdown[0]->ToTH1(pot),fluxdown[0]->ToTH1(pot)};// {fluxdown[0]->ToTH1(pot), xsecdown[0]->ToTH1(pot),calibdown[0]->ToTH1(pot), lightdown[0]->ToTH1(pot),ckv[0]->ToTH1(pot), calib_shape[0]->ToTH1(pot)}; ///comment lightup[3]->ToTH1(pot), ckv[3]->ToTH1(pot), calib_shape[3]->ToTH1(pot), fluxup[3]->ToTH1(pot), xsecup[3]->ToTH1(pot) if plotting just calib.
600 
601 
602  TH1* hnom_NSel_Cuts = nominal[0]->ToTH1(pot);
603  hnom_NSel_Cuts->SetTitle(";BDTG; Selected Events/8.09 #times 10^{20}POT");
604  // hnom_sig_nopresel->SetTitleSize(";0.05; 0.04");
605  std::cout << "*************Bins nominal*************: " << hnom_NSel_Cuts->GetXaxis()->GetNbins() << std::endl;
606  gROOT->ProcessLine(".! mkdir -p plots");
607 
608  TCanvas* Prong_CVN_Gamma_ID_NSel_Cuts = new TCanvas("Prong_CVN_Gamma_ID_NSel_Cuts", "Prong_CVN_Gamma_ID_NSel_Cuts");
609  Prong_CVN_Gamma_ID_NSel_Cuts ->SetLeftMargin(0.15);
610 
611  PlotWithSystErrorBandWidth(hnom_NSel_Cuts,
612  vector_up_NSel_Cuts,
613  vector_dwn_NSel_Cuts, -1, -1, 1.3f, true);
614 
615  NDSimulation();
616 
617  Prong_CVN_Gamma_ID_NSel_Cuts->SaveAs(outDir2+CurrentFile[1].name.c_str()+CurrentDate[0].date.c_str()+".pdf");
618  // Prong_CVN_Gamma_ID_NSel_Cuts->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/Prong_CVN_Gamma_ID_NSel_Cuts_0_3_high_Combined_OneKnob1_out_3_MaNCRES.pdf");
619  // CvtxX_sig_nopresel->SaveAs("plots/vtxX_sig_nopresel.png");
620  // Prong_CVN_Gamma_ID_NSel_Cuts->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/Prong_CVN_Gamma_ID_NSel_Cuts_0_3_high.pdf");
621  // Prong_CVN_Gamma_ID_NSel_Cuts->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/Calibration/Prong_CVN_Gamma_ID_NSel_Cuts_calib_0_3_high.pdf");
622  //Prong_CVN_Gamma_ID_NSel_Cuts->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/Calibration_shape/Prong_CVN_Gamma_ID_NSel_Cuts__calib_shape_0_3_high.pdf");
623  //Prong_CVN_Gamma_ID_NSel_Cuts->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/Cherenkov/Prong_CVN_Gamma_ID_NSel_Cuts_Cherenkov_0_3_high.pdf");
624  // Prong_CVN_Gamma_ID_NSel_Cuts->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/GENIE/Prong_CVN_Gamma_ID_NSel_Cuts_GENIE_0_3_high.pdf");
625  // Prong_CVN_Gamma_ID_NSel_Cuts->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/light_level/Prong_CVN_Gamma_ID_NSel_Cuts_light_level_0_3_high.pdf");
626  // Prong_CVN_Gamma_ID_NSel_Cuts->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/ppfx/Prong_CVN_Gamma_ID_NSel_Cuts_ppfx_0_3_high.pdf");
627 
628  ///***********************************************************************************************************************************************************************************************////
629 
630  //For Signal cut: kVtxInFiducial && kContainment && k2Prongs && kPi0RemIDND && kNCurrent && kPi0 (0.3geV) ; (Signal PreSel Cuts)
631 
632  TH1* signal_xsecdown = xsecdown[1]->ToTH1(pot); ///comment this if just caliip or down is needed for plots
633  TH1* signal_xsecup = xsecup[1]->ToTH1(pot); ///comment this if just caliip or down is needed for plots
634 
635  TH1* signal_fluxdown = fluxdown[1]->ToTH1(pot); ///comment this if just caliip or down is needed for plots
636  TH1* signal_fluxup = fluxup[1]->ToTH1(pot); ///comment this if just caliip or down is needed for plots
637 
638  // TH1* signal_calib_shape = calib_shape[1]->ToTH1(pot); ///comment this if just caliip or down is needed for plots
639 
640  // TH1* signal_ckv = ckv[1]->ToTH1(pot); ///comment this if just caliip or down is needed for plots
641 
642  // TH1* signal_lightup = lightup[1]->ToTH1(pot); ///comment this if just caliip or down is needed for plots
643  // TH1* signal_lightdown = lightdown[1]->ToTH1(pot); ///comment this if just caliip or down is needed for plots
644 
645  // TH1* signal_calibdown = calibdown[1]->ToTH1(pot);
646  // TH1* signal_calibup = calibup[1]->ToTH1(pot);
647 
648 
649  // sig signal distributions
650  std::vector<TH1*> vector_up_signal = {xsecup[1]->ToTH1(pot),fluxup[1]->ToTH1(pot)};//{fluxup[1]->ToTH1(pot), xsecup[1]->ToTH1(pot), calibup[1]->ToTH1(pot), lightup[1]->ToTH1(pot), ckv[1]->ToTH1(pot), calib_shape[1]->ToTH1(pot)}; ///comment lightup[3]->ToTH1(pot), ckv[3]->ToTH1(pot), calib_shape[3]->ToTH1(pot), fluxup[3]->ToTH1(pot), xsecup[3]->ToTH1(pot) if plotting just calib.
651  std::vector<TH1*> vector_dwn_signal = {xsecdown[1]->ToTH1(pot),fluxdown[1]->ToTH1(pot)};// {fluxdown[1]->ToTH1(pot), xsecdown[1]->ToTH1(pot), calibdown[1]->ToTH1(pot), lightdown[1]->ToTH1(pot), ckv[1]->ToTH1(pot), calib_shape[1]->ToTH1(pot)}; ///comment lightup[3]->ToTH1(pot), ckv[3]->ToTH1(pot), calib_shape[3]->ToTH1(pot), fluxup[3]->ToTH1(pot), xsecup[3]->ToTH1(pot) if plotting just calib.
652 
653  TH1* hnom_signal = nominal[1]->ToTH1(pot);
654  hnom_signal->SetTitle(";BDTG; Signal events/8.09 #times 10^{20}POT");
655  // hnom_signal->SetTitleSize(";0.05; 0.04");
656  std::cout << "*************Bins nominal*************: " << hnom_signal->GetXaxis()->GetNbins() << std::endl;
657  gROOT->ProcessLine(".! mkdir -p plots");
658 
659  TCanvas* Prong_CVN_Gamma_ID_NSig_Cuts = new TCanvas("Prong_CVN_Gamma_ID_NSig_Cuts", "Prong_CVN_Gamma_ID_NSig_Cuts");
660  Prong_CVN_Gamma_ID_NSig_Cuts->SetLeftMargin(0.15);
661 
662  PlotWithSystErrorBandWidth(hnom_signal,
663  vector_up_signal,
664  vector_dwn_signal, -1, -1, 1.3f, true);
665 
666  NDSimulation();
667 
668  Prong_CVN_Gamma_ID_NSig_Cuts->SaveAs(outDir2+CurrentFile[2].name.c_str()+CurrentDate[0].date.c_str()+".pdf");
669  //Prong_CVN_Gamma_ID_NSig_Cuts->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/Prong_CVN_Gamma_ID_NSig_Cuts_0_3_high_Combined_OneKnob1_out_3_MaNCRES.pdf");
670  //Prong_CVN_Gamma_ID_NSig_Cuts->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/Prong_CVN_Gamma_ID_NSig_Cuts_0_3_high.pdf");
671  // Prong_CVN_Gamma_ID_NSig_Cuts->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/Calibration/Prong_CVN_Gamma_ID_NSig_Cuts_calib_0_3_high.pdf");
672  // Prong_CVN_Gamma_ID_NSig_Cuts->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/Calibration_shape/Prong_CVN_Gamma_ID_NSig_Cuts__calib_shape_0_3_high.pdf");
673  // Prong_CVN_Gamma_ID_NSig_Cuts->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/Cherenkov/Prong_CVN_Gamma_ID_NSig_Cuts_Cherenkov_0_3_high.pdf");
674  // Prong_CVN_Gamma_ID_NSig_Cuts->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/GENIE/Prong_CVN_Gamma_ID_NSig_Cuts_GENIE_0_3_high.pdf");
675  // Prong_CVN_Gamma_ID_NSig_Cuts->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/light_level/Prong_CVN_Gamma_ID_NSig_Cuts_light_level_0_3_high.pdf");
676  // Prong_CVN_Gamma_ID_NSig_Cuts->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/ppfx/Prong_CVN_Gamma_ID_NSig_Cuts_ppfx_0_3_high.pdf");
677  ///***********************************************************************************************************************************************************************************************////
678 
679  // for true signal
680 
681 
682  //true signal fiducial (Nsig,fid only): k2Prongs && kVtxInFiducial && kNCurrent && kPi0 (0.3geV);
683 
684  TH1* truesignal_xsecdown = xsecdown[2]->ToTH1(pot); ///comment this if just caliip or down is needed for plots
685  TH1* truesignal_xsecup = xsecup[2]->ToTH1(pot); ///comment this if just caliip or down is needed for plots
686 
687  TH1* truesignal_fluxdown = fluxdown[2]->ToTH1(pot); ///comment this if just caliip or down is needed for plots
688  TH1* truesignal_fluxup = fluxup[2]->ToTH1(pot); ///comment this if just caliip or down is needed for plots
689 
690  //TH1* truesignal_calib_shape = calib_shape[2]->ToTH1(pot); ///comment this if just caliip or down is needed for plots
691 
692  ///TH1* truesignal_ckv = ckv[2]->ToTH1(pot); ///comment this if just caliip or down is needed for plots
693 
694  // TH1* truesignal_lightup = lightup[2]->ToTH1(pot); ///comment this if just caliip or down is needed for plots
695  // TH1* truesignal_lightdown = lightdown[2]->ToTH1(pot); ///comment this if just caliip or down is needed for plots
696 
697  // TH1* truesignal_calibdown = calibdown[2]->ToTH1(pot);
698  // TH1* truesignal_calibup = calibup[2]->ToTH1(pot);
699 
700 
701  // sig signal distributions
702  std::vector<TH1*> vector_up_truesignal = {xsecup[2]->ToTH1(pot),fluxup[2]->ToTH1(pot)};// {fluxup[2]->ToTH1(pot), xsecup[2]->ToTH1(pot), calibup[2]->ToTH1(pot), lightup[2]->ToTH1(pot), ckv[2]->ToTH1(pot), calib_shape[2]->ToTH1(pot)}; ///comment lightup[3]->ToTH1(pot), ckv[3]->ToTH1(pot), calib_shape[3]->ToTH1(pot), fluxup[3]->ToTH1(pot), xsecup[3]->ToTH1(pot) if plotting just calib.
703  std::vector<TH1*> vector_dwn_truesignal = {xsecdown[2]->ToTH1(pot),fluxdown[2]->ToTH1(pot)};// {fluxdown[2]->ToTH1(pot), xsecdown[2]->ToTH1(pot), calibdown[2]->ToTH1(pot), lightdown[2]->ToTH1(pot), ckv[2]->ToTH1(pot), calib_shape[2]->ToTH1(pot)}; ///comment lightup[3]->ToTH1(pot), ckv[3]->ToTH1(pot), calib_shape[3]->ToTH1(pot), fluxup[3]->ToTH1(pot), xsecup[3]->ToTH1(pot) if plotting just calib.
704 
705 
706 
707  TH1* hnom_truesignal = nominal[2]->ToTH1(pot);
708  hnom_truesignal->SetTitle(";BDTG; True signal events/8.09 #times 10^{20}POT");
709  // hnom_truesignal->SetTitleSize(";0.05; 0.04");
710 
711  gROOT->ProcessLine(".! mkdir -p plots");
712 
713  TCanvas* Prong_CVN_Gamma_ID_NtrueSig_Cut= new TCanvas("Prong_CVN_Gamma_ID_NtrueSig_Cut", "Prong_CVN_Gamma_ID_NtrueSig_Cut");
714  Prong_CVN_Gamma_ID_NtrueSig_Cut->SetLeftMargin(0.15);
715 
716  PlotWithSystErrorBandWidth(hnom_truesignal,
717  vector_up_truesignal,
718  vector_dwn_truesignal, -1, -1, 1.3f, true);
719 
720  NDSimulation();
721 
722  Prong_CVN_Gamma_ID_NtrueSig_Cut->SaveAs(outDir2+CurrentFile[3].name.c_str()+CurrentDate[0].date.c_str()+".pdf");
723  //Prong_CVN_Gamma_ID_NtrueSig_Cut->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/Prong_CVN_Gamma_ID_NtrueSig_Cuts_0_3_high_Combined_OneKnob2_out_3_MaNCRES_MaCOHPi.pdf");
724  //Prong_CVN_Gamma_ID_NtrueSig_Cut->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/Prong_CVN_Gamma_ID_NtrueSig_Cuts_0_3_high.pdf");
725  //Prong_CVN_Gamma_ID_NtrueSig_Cut->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/Calibration/Prong_CVN_Gamma_ID_NtrueSig_Cuts_calib_0_3_high.pdf");
726  // Prong_CVN_Gamma_ID_NtrueSig_Cut->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/Calibration_shape/Prong_CVN_Gamma_ID_NtrueSig_Cuts__calib_shape_0_3_high.pdf");
727  // Prong_CVN_Gamma_ID_NtrueSig_Cut->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/Cherenkov/Prong_CVN_Gamma_ID_NtrueSig_Cuts_Cherenkov_0_3_high.pdf");
728  // Prong_CVN_Gamma_ID_NtrueSig_Cut->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/GENIE/Prong_CVN_Gamma_ID_NtrueSig_Cuts_GENIE_0_3_high.pdf");
729  // Prong_CVN_Gamma_ID_NtrueSig_Cut->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/light_level/Prong_CVN_Gamma_ID_NtrueSig_Cuts_light_level_0_3_high.pdf");
730  // Prong_CVN_Gamma_ID_NtrueSig_Cut->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/ppfx/Prong_CVN_Gamma_ID_NtrueSig_Cuts_ppfx_0_3_high.pdf");
731  ///***********************************************************************************************************************************************************************************************////
732 
733 
734 
735  TH1* hNbkg_Cuts = new TH1F("Nbkg_Cuts","Nbkg_Cuts", 100,-1,1);
736 
737 
738  hNbkg_Cuts->GetXaxis()->CenterTitle();
739  hNbkg_Cuts->GetYaxis()->CenterTitle();
740 
741  TH1* hNbkg_Cuts_xsecdown = new TH1F("xsecdown","xsecdown",100,-1,1); ///comment this if just caliip or down is needed for plots
742  TH1* hNbkg_Cuts_xsecup = new TH1F("xsecup","xsecup",100,-1,1); ///comment this if just caliip or down is needed for plots
743  TH1* hNbkg_Cuts_fluxdown = new TH1F("fluxdown","fluxdown",100,-1,1); ///comment this if just caliip or down is needed for plots
744  TH1* hNbkg_Cuts_fluxup = new TH1F("fluxup","fluxup",100,-1,1); ///comment this if just caliip or down is needed for plots
745 
746  // TH1* hNbkg_Cuts_calib_shape = new TH1F("calibshape","calibshape",100,-1,1); ///comment this if just caliip or down is needed for plots
747 
748  // TH1* hNbkg_Cuts_ckv = new TH1F("ckv","ckv",100,-1,1); ///comment this if just caliip or down is needed for plots
749 
750  // TH1* hNbkg_Cuts_lightup = new TH1F("lightup","lightup",100,-1,1); ///comment this if just caliip or down is needed for plots
751  // TH1* hNbkg_Cuts_lightdown = new TH1F("lightdown","lightdown",100,-1,1); ///comment this if just caliip or down is needed for plots
752 
753  // TH1* hNbkg_Cuts_calibdown = new TH1F("calibdown","calibdown",100,-1,1);
754  // TH1* hNbkg_Cuts_calibup = new TH1F("calibup","calibup",100,-1,1) ;
755 
756  TH1* hPur=new TH1F("hPur","",100,-1,1);
757  float pur=0.0;
758  TH1* hFom=new TH1F("hPur","",100,-1,1);
759  float fom=0.0;
760 
761  for(int n = 1; n < hnom_Nbkg_Cuts->GetNbinsX()+1; ++n){
762 
763 
764  const float nNSel_Cuts = hnom_NSel_Cuts->Integral(n, 100);
765  const float nNbkg_Cuts = hnom_Nbkg_Cuts->Integral(n, 100);
766 
767  ////Purity and FOM plots////////////////////////////////////////////////////////////////////////////////////////////////////////
768  float Purity_den = nNSel_Cuts;
769  float Purity_num = nNSel_Cuts-nNbkg_Cuts;
770 
771  float Fom_den = nNSel_Cuts;
772  float Fom_num = nNSel_Cuts-nNbkg_Cuts;
773  if (Purity_num!=0 ){
774  pur=Purity_num/nNSel_Cuts;
775 
776  cout<<"Total Selected:"<<nNSel_Cuts<<endl;
777  cout<<"Background:"<<nNbkg_Cuts<<endl;
778  cout<<"Purity_num:"<<nNSel_Cuts-nNbkg_Cuts<<endl;
779  cout<<"Purity:"<<Purity_num/nNSel_Cuts<<endl;
780  }
781 
782  if (Fom_num!=0 ){
783  fom= Fom_num/sqrt(nNSel_Cuts);
784 
785  cout<<"Total Selected:"<<nNSel_Cuts<<endl;
786  cout<<"Background:"<<nNbkg_Cuts<<endl;
787  cout<<"Fom_num:"<<nNSel_Cuts-nNbkg_Cuts<<endl;
788  cout<<"FOM:"<<Fom_num/sqrt(nNSel_Cuts)<<endl;
789  }
790 
791  hPur->SetBinContent(n,pur);
792  hFom->SetBinContent(n,fom);
793  hNbkg_Cuts->SetBinContent(n,nNbkg_Cuts);
794 
795  // hNbkg_Cuts_calibup->SetBinContent(n, Nbkg_Cuts_calibup->Integral(n,100));
796  // hNbkg_Cuts_calibdown->SetBinContent(n, Nbkg_Cuts_calibdown->Integral(n, 100));
797 
798  hNbkg_Cuts_xsecup->SetBinContent(n, Nbkg_Cuts_xsecup->Integral(n,100)); // comment this if just caliip or down is needed for plots
799  hNbkg_Cuts_xsecdown->SetBinContent(n, Nbkg_Cuts_xsecdown->Integral(n,100)); // comment this if just caliip or down is needed for plots
800 
801  // hNbkg_Cuts_lightup->SetBinContent(n, Nbkg_Cuts_lightup->Integral(n,100)); ///comment this if just caliip or down is needed for plots
802  // hNbkg_Cuts_lightdown->SetBinContent(n, Nbkg_Cuts_lightdown->Integral(n,100)); ///comment this if just caliip or down is needed for plots
803 
804  hNbkg_Cuts_fluxup->SetBinContent(n, Nbkg_Cuts_fluxup->Integral(n,100)); ///comment this if just caliip or down is needed for plots
805  hNbkg_Cuts_fluxdown->SetBinContent(n, Nbkg_Cuts_fluxdown->Integral(n,100)); // comment this if just caliip or down is needed for plots
806 
807  // hNbkg_Cuts_calib_shape->SetBinContent(n,Nbkg_Cuts_calib_shape->Integral(n,100)); ///comment this if just caliip or down is needed for plots
808  // hNbkg_Cuts_ckv->SetBinContent(n, Nbkg_Cuts_ckv->Integral(n,100));
809  ///comment this if just caliip or down is needed for plots
810 
811  }
812 
813 
814  std::vector<TH1*> vector_up_Nbkg_Cuts_binbybin = {hNbkg_Cuts_xsecup,hNbkg_Cuts_fluxup};//{hNbkg_Cuts_fluxup, hNbkg_Cuts_xsecup, hNbkg_Cuts_calibup, hNbkg_Cuts_lightup, hNbkg_Cuts_ckv, hNbkg_Cuts_calib_shape}; ///comment lightup[3]->ToTH1(pot), ckv[3]->ToTH1(pot), calib_shape[3]->ToTH1(pot), fluxup[3]->ToTH1(pot), xsecup[3]->ToTH1(pot) if plotting just calib.
815 
816  std::vector<TH1*> vector_dwn_Nbkg_Cuts_binbybin = {hNbkg_Cuts_xsecdown,hNbkg_Cuts_fluxdown};// {hNbkg_Cuts_fluxdown, hNbkg_Cuts_xsecdown, hNbkg_Cuts_calibdown, hNbkg_Cuts_lightdown, hNbkg_Cuts_ckv, hNbkg_Cuts_calib_shape}; ///comment lightup[3]->ToTH1(pot), ckv[3]->ToTH1(pot), calib_shape[3]->ToTH1(pot), fluxup[3]->ToTH1(pot), xsecup[3]->ToTH1(pot) if plotting just calib.
817 
818 
819  /*-----------------------------------------
820  Plotting Nbkg with systematics
821  -----------------------------------------*/
822  TCanvas* cNbkg_Cuts = new TCanvas("cNbkg_Cuts","cNbkg_Cuts");
823  cNbkg_Cuts->SetLeftMargin(0.15);
824 
825  hNbkg_Cuts->SetTitle(";BDTG Cut value; Background Events/8.09 #times 10^{20} POT");
826  // hbkg_nopresel->SetTitleSize("; 0.05; 0.04")
827 
828 
829  PlotWithSystErrorBandWidth(hNbkg_Cuts,
830  vector_up_Nbkg_Cuts_binbybin,
831  vector_dwn_Nbkg_Cuts_binbybin, -1, -1, 1.3f, true);
832 
833  NDSimulation();
834  cNbkg_Cuts->SaveAs(outDir2+CurrentFile[4].name.c_str()+CurrentDate[0].date.c_str()+".pdf");
835  //cNbkg_Cuts->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/Prong_CVN_Gamma_ID_Cut_WithErrorBand_Nbkg_Cuts_0_3_high_Combined_OneKnob1_out_3_MaNCRES.pdf");
836  //cNbkg_Cuts->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/Prong_CVN_Gamma_ID_Cut_WithErrorBand_Nbkg_Cuts_0_3_high.pdf");
837  //cNbkg_Cuts->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/Calibration/Prong_CVN_Gamma_ID_Nbkg_Cuts_calib_0_3_high.pdf");
838  // cNbkg_Cuts->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/Calibration_shape/Prong_CVN_Gamma_ID_Nbkg_Cuts__calib_shape_0_3_high.pdf");
839  // cNbkg_Cuts->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/Cherenkov/Prong_CVN_Gamma_ID_Nbkg_Cuts_Cherenkov_0_3_high.pdf");
840  // cNbkg_Cuts->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/GENIE/Prong_CVN_Gamma_ID_Nbkg_Cuts_GENIE_0_3_high.pdf");
841  // cNbkg_Cuts->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/light_level/Prong_CVN_Gamma_ID_Nbkg_Cuts_light_level_0_3_high.pdf");
842  // cNbkg_Cuts->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/ppfx/Prong_CVN_Gamma_ID_Nbkg_Cuts_ppfx_0_3_high.pdf");
843 
844  /*-----------------------------------------
845  Plotting Purity
846  -----------------------------------------*/
847  TCanvas* cpur = new TCanvas("cpur","cpur");
848  cpur->SetLeftMargin(0.15);
849  hPur->SetLineColor(kRed);
850  hPur->GetYaxis()->SetTitle("Purity");
851  hPur->GetXaxis()->SetTitle("BDTG Cut Value");
852  hPur->Draw();
853  NDSimulation();
854  cpur->SaveAs(outDir2+CurrentFile[9].name.c_str()+CurrentDate[0].date.c_str()+".pdf");
855  // // cpur->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/Calibration/fom_calib_0_635_high.pdf");
856  // // cpur->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/Calibration_shape/fom_calib_shape_0_635_high.pdf");
857  // // cpur->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/Cherenkov/fom_Cherenkov_0_635_high.pdf");
858  // // cpur->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/GENIE/fom_GENIE_0_635_high.pdf");
859  // // cpur->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/light_level/fom_light_level_0_635_high.pdf");
860  // // cpur->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/ppfx/fom_ppfx_0_635_high.pdf");
861  /*-----------------------------------------
862  Plotting FOM
863  -----------------------------------------*/
864  TCanvas* cfom = new TCanvas("cfom","cfom");
865  cfom->SetLeftMargin(0.15);
866  hFom->SetLineColor(kRed);
867  hFom->GetYaxis()->SetTitle("FOM");
868  hFom->GetXaxis()->SetTitle("BDTG Cut Value");
869  hFom->Draw();
870  NDSimulation();
871  cfom->SaveAs(outDir2+CurrentFile[10].name.c_str()+CurrentDate[0].date.c_str()+".pdf");
872 
873 
874 
875  ///*************************************************************************************************************************************************************************************/////
876  ////Total PreSel : kVtxInFiducial && kContainment && k2Prongs && kPi0RemIDND && kNpCVN1gammaID;
877 
878  //Selected Bin by Bin
879 
880 
881  TH1* hNSel_Cuts = new TH1F("NSel_Cuts","NSel_Cuts", 100,-1,1 );
882 
883  hNSel_Cuts->GetXaxis()->CenterTitle();
884  hNSel_Cuts->GetYaxis()->CenterTitle();
885 
886 
887  TH1* hNSel_Cuts_xsecdown = new TH1F("xsecdown","xsecdown",100,-1,1); ///comment this if just caliip or down is needed for plots
888  TH1* hNSel_Cuts_xsecup = new TH1F("xsecup","xsecup",100,-1,1); ///comment this if just caliip or down is needed for plots
889 
890  TH1* hNSel_Cuts_fluxdown = new TH1F("fluxdown","fluxdown",100,-1,1); ///comment this if just caliip or down is needed for plots
891  TH1* hNSel_Cuts_fluxup = new TH1F("fluxup","fluxup",100,-1,1); ///comment this if just caliip or down is needed for plots
892 
893  // TH1* hNSel_Cuts_calib_shape = new TH1F("calibshape","calibshape",100,-1,1); ///comment this if just caliip or down is needed for plots
894 
895  // TH1* hNSel_Cuts_ckv = new TH1F("ckv","ckv", 100,-1,1 ); ///comment this if just caliip or down is needed for plots
896 
897  // TH1* hNSel_Cuts_lightup = new TH1F("lightup","lightup",100,-1,1); ///comment this if just caliip or down is needed for plots
898  // TH1* hNSel_Cuts_lightdown = new TH1F("lightdown","lightdown", 100,-1,1); ///comment this if just caliip or down is needed for plots
899 
900  // TH1* hNSel_Cuts_calibdown = new TH1F("calibdown","calibdown",100,-1,1 );
901  // TH1* hNSel_Cuts_calibup = new TH1F("calibup","calibup" ,100,-1,1);
902 
903 
904 
905  for(int n = 1; n < hnom_NSel_Cuts->GetNbinsX() + 1; ++n){
906 
907  const float nsel_cuts = hnom_NSel_Cuts->Integral(n,100);
908 
909  hNSel_Cuts->SetBinContent(n, nsel_cuts);
910 
911  // hNSel_Cuts_calibup->SetBinContent(n,NSel_Cuts_calibup->Integral(n,100));
912  // hNSel_Cuts_calibdown->SetBinContent(n, NSel_Cuts_calibdown->Integral(n,100));
913 
914  hNSel_Cuts_xsecup->SetBinContent(n, NSel_Cuts_xsecup->Integral(n,100)); // comment this if just caliip or down is needed for plots
915  hNSel_Cuts_xsecdown->SetBinContent(n, NSel_Cuts_xsecdown->Integral(n,100)); //comment this if just caliip or down is needed for plots
916 
917  // hNSel_Cuts_lightup->SetBinContent(n, NSel_Cuts_lightup->Integral(n,100)); ///comment this if just caliip or down is needed for plots
918  // hNSel_Cuts_lightdown->SetBinContent(n, NSel_Cuts_lightdown->Integral(n,100)); ///comment this if just caliip or down is needed for plots
919 
920  hNSel_Cuts_fluxup->SetBinContent(n, NSel_Cuts_fluxup->Integral(n,100)); ///comment this if just caliip or down is needed for plots
921  hNSel_Cuts_fluxdown->SetBinContent(n, NSel_Cuts_fluxdown->Integral(n,100));
922 
923  // hNSel_Cuts_calib_shape->SetBinContent(n, NSel_Cuts_calib_shape->Integral(n,100)); ///comment this if just caliip or down is needed for plots
924  // hNSel_Cuts_ckv->SetBinContent(n, NSel_Cuts_ckv->Integral(n,100));
925  ///comment this if just caliip or down is needed for plots
926 
927  }
928 
929 
930  std::vector<TH1*> vector_up_NSel_Cuts_binbybin = {hNSel_Cuts_xsecup,hNSel_Cuts_fluxup};// {hNSel_Cuts_fluxup, hNSel_Cuts_xsecup, hNSel_Cuts_calibup, hNSel_Cuts_lightup, hNSel_Cuts_ckv, hNSel_Cuts_calib_shape};
931 ///comment lightup[3]->ToTH1(pot), ckv[3]->ToTH1(pot), calib_shape[3]->ToTH1(pot), fluxup[3]->ToTH1(pot), xsecup[3]->ToTH1(pot) if plotting just calib.
932 
933 
934  std::vector<TH1*> vector_dwn_NSel_Cuts_binbybin = {hNSel_Cuts_xsecdown,hNSel_Cuts_fluxdown};//{hNSel_Cuts_fluxdown, hNSel_Cuts_xsecdown, hNSel_Cuts_calibdown, hNSel_Cuts_lightdown, hNSel_Cuts_ckv, hNSel_Cuts_calib_shape};
935 ///comment lightup[3]->ToTH1(pot), ckv[3]->ToTH1(pot), calib_shape[3]->ToTH1(pot), fluxup[3]->ToTH1(pot), xsecup[3]->ToTH1(pot) if plotting just calib.
936 
937 
938 
939  /*-----------------------------------------
940  Plotting NSel with systematics
941  -----------------------------------------*/
942  TCanvas* cNSel_Cuts = new TCanvas("cNSel_Cuts","cNSel_Cuts");
943  cNSel_Cuts->SetLeftMargin(0.15);
944 
945  hNSel_Cuts->SetTitle(";BDTG Cut Value;Selected Events/8.09 #times 10^{20} POT");
946  // hsig_nopresel->SetTitleSize("; 0.05; 0.04")
947 
948 
949  PlotWithSystErrorBandWidth(hNSel_Cuts,
950  vector_up_NSel_Cuts_binbybin,
951  vector_dwn_NSel_Cuts_binbybin, -1, -1, 1.3f, true);
952 
953  NDSimulation();
954 
955  cNSel_Cuts->SaveAs(outDir2+CurrentFile[5].name.c_str()+CurrentDate[0].date.c_str()+".pdf");
956  //cNSel_Cuts->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/Prong_CVN_Gamma_IDcut_WithErrorBand_NSel_Cuts_0_3_high_Combined_OneKnob1_out_3_MaNCRES.pdf");
957  //cNSel_Cuts->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/Prong_CVN_Gamma_IDcut_WithErrorBand_NSel_Cuts_0_3_high.pdf");
958  // cNSel_Cuts->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/Calibration/Prong_CVN_Gamma_IDcut_WithErrorBand_NSel_Cuts_calib_0_3_high.pdf");
959  // cNSel_Cuts->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/Calibration_shape/Prong_CVN_Gamma_IDcut_WithErrorBand_NSel_Cuts_calib_shape_0_3_high.pdf");
960  // cNSel_Cuts->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/Cherenkov/Prong_CVN_Gamma_IDcut_WithErrorBand_NSel_Cuts_Cherenkov_0_3_high.pdf");
961  // cNSel_Cuts->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/GENIE/Prong_CVN_Gamma_IDcut_WithErrorBand_NSel_Cuts_GENIE_0_3_high.pdf");
962  // cNSel_Cuts->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/light_level/Prong_CVN_Gamma_IDcut_WithErrorBand_NSel_Cuts_light_level_0_3_high.pdf");
963  // cNSel_Cuts->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/ppfx/Prong_CVN_Gamma_IDcut_WithErrorBand_NSel_Cuts_ppfx_0_3_high.pdf");
964 
965 //****************************************************************************************************************************************************************************************////
966 
967 
968  //For Signal: ///Signal Presel : kVtxInFiducial && kContainment && k2Prongs && kPi0RemIDND && kNCurrent && kPi0 (0.3geV) ; (Signal PreSel Cuts)
969 
970  TH1* hsignal = new TH1F("signal","signal", 100,-1,1);
971 
972 
973  hsignal->GetXaxis()->CenterTitle();
974  hsignal->GetYaxis()->CenterTitle();
975 
976 
977  TH1* hsignal_xsecdown = new TH1F("xsecdown","xsecdown", 100,-1,1);
978  TH1* hsignal_xsecup = new TH1F("xsecup","xsecup", 100,-1,1);
979 
980  TH1* hsignal_fluxdown = new TH1F("fluxdown","fluxdown", 100,-1,1);
981  TH1* hsignal_fluxup = new TH1F("fluxup","fluxup", 100,-1,1);
982 
983  // TH1* hsignal_calib_shape = new TH1F("calibshape","calibshape", 100,-1,1);
984 
985  // TH1* hsignal_ckv = new TH1F("ckv","ckv", 100,-1,1);
986 
987  // TH1* hsignal_lightup = new TH1F("lightup","lightup", 100,-1,1);
988  // TH1* hsignal_lightdown = new TH1F("lightdown","lightdown", 100,-1,1);
989 
990  // TH1* hsignal_calibdown = new TH1F("calibdown","calibdown", 100,-1,1);
991  // TH1* hsignal_calibup = new TH1F("calibup","calibup", 100,-1,1);
992 
993 
994 
995 
996  for(int n = 1; n < hnom_signal->GetNbinsX()+1; ++n){
997 
998  const float nsig = hnom_signal->Integral(n,100);
999 
1000  hsignal->SetBinContent(n, nsig);
1001 
1002  // hsignal_calibup->SetBinContent(n,signal_calibup->Integral(n,100));
1003  // hsignal_calibdown->SetBinContent(n, signal_calibdown->Integral(n,100));
1004 
1005  hsignal_xsecup->SetBinContent(n, signal_xsecup->Integral(n,100));
1006  hsignal_xsecdown->SetBinContent(n, signal_xsecdown->Integral(n,100));
1007 
1008  // hsignal_lightup->SetBinContent(n, signal_lightup->Integral(n,100));
1009  // hsignal_lightdown->SetBinContent(n, signal_lightdown->Integral(n,100));
1010 
1011  hsignal_fluxup->SetBinContent(n, signal_fluxup->Integral(n,100));
1012  hsignal_fluxdown->SetBinContent(n, signal_fluxdown->Integral(n,100));
1013 
1014  // hsignal_calib_shape->SetBinContent(n, signal_calib_shape->Integral(n,100));
1015  // hsignal_ckv->SetBinContent(n, signal_ckv->Integral(n,100));
1016 
1017  }
1018 
1019 
1020  std::vector<TH1*> vector_up_signal_binbybin = {hsignal_xsecup,hsignal_fluxup};// {hsignal_fluxup, hsignal_xsecup, hsignal_calibup, hsignal_lightup, hsignal_ckv, hsignal_calib_shape};
1021 
1022  std::vector<TH1*> vector_dwn_signal_binbybin = {hsignal_xsecdown,hsignal_fluxdown};// {hsignal_fluxdown, hsignal_xsecdown,hsignal_calibdown, hsignal_lightdown, hsignal_ckv, hsignal_calib_shape};
1023 
1024 
1025 
1026  /*-----------------------------------------
1027  Plotting signal with systematics
1028  -----------------------------------------*/
1029  TCanvas* csignal = new TCanvas("csignal","csignal");
1030  csignal->SetLeftMargin(0.15);
1031 
1032  hsignal->SetTitle(";BDTG Cut Value;Signal Events/8.09 #times 10^{20} POT");
1033 
1034  PlotWithSystErrorBandWidth(hsignal,
1035  vector_up_signal_binbybin,
1036  vector_dwn_signal_binbybin, -1, -1, 1.3f, true);
1037 
1038  NDSimulation();
1039 
1040  csignal->SaveAs(outDir2+CurrentFile[6].name.c_str()+CurrentDate[0].date.c_str()+".pdf");
1041  //csignal->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/Prong_CVN_Gamma_IDcut_WithErrorBand_NSig_Cuts_0_3_high_Combined_OneKnob1_out_3_MaNCRES.pdf");
1042  //csignal->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/Prong_CVN_Gamma_IDcut_WithErrorBand_NSig_Cuts_0_3_high.pdf");
1043  // csignal->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/Calibration/Prong_CVN_Gamma_IDcut_WithErrorBand_NSig_Cuts_calib_0_3_high.pdf");
1044  // csignal->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/Calibration_shape/Prong_CVN_Gamma_IDcut_WithErrorBand_NSig_Cuts_calib_shape_0_3_high.pdf");
1045  // csignal->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/Cherenkov/Prong_CVN_Gamma_IDcut_WithErrorBand_NSig_Cuts_Cherenkov_0_3_high.pdf");
1046  // csignal->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/GENIE/Prong_CVN_Gamma_IDcut_WithErrorBand_NSig_Cuts_GENIE_0_3_high.pdf");
1047  // csignal->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/light_level/Prong_CVN_Gamma_IDcut_WithErrorBand_NSig_Cuts_light_level_0_3_high.pdf");
1048  // csignal->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/ppfx/Prong_CVN_Gamma_IDcut_WithErrorBand_NSig_Cuts_ppfx_0_3_high.pdf");
1049 
1050 //****************************************************************************************************************************************************************************************////
1051 
1052  //For true Signal
1053 
1054  TH1* htruesignal = new TH1F("","", 100,-1,1);
1055 
1056 
1057  htruesignal->GetXaxis()->CenterTitle();
1058  htruesignal->GetYaxis()->CenterTitle();
1059 
1060 
1061  TH1* htruesignal_xsecdown = new TH1F("xsecdown","xsecdown", 100,-1,1);
1062  TH1* htruesignal_xsecup = new TH1F("xsecup","xsecup", 100,-1,1);
1063 
1064  TH1* htruesignal_fluxdown = new TH1F("fluxdown","fluxdown", 100,-1,1);
1065  TH1* htruesignal_fluxup = new TH1F("fluxup","fluxup", 100,-1,1);
1066 
1067  // TH1* htruesignal_calib_shape = new TH1F("calibshape","calibshape", 100,-1,1);
1068 
1069  // TH1* htruesignal_ckv = new TH1F("ckv","ckv", 100,-1,1);
1070 
1071  // TH1* htruesignal_lightup = new TH1F("lightup","lightup", 100,-1,1);
1072  // TH1* htruesignal_lightdown = new TH1F("lightdown","lightdown", 100,-1,1);
1073 
1074  // TH1* htruesignal_calibdown = new TH1F("calibdown","calibdown", 100,-1,1);
1075  // TH1* htruesignal_calibup = new TH1F("calibup","calibup", 100,-1,1);
1076 
1077 
1078 
1079 
1080  for(int n = 1; n < hnom_truesignal->GetNbinsX()+1; ++n){
1081 
1082  const float nsig = hnom_truesignal->Integral(n,100);
1083 
1084  htruesignal->SetBinContent(n, nsig);
1085 
1086  // htruesignal_calibup->SetBinContent(n,truesignal_calibup->Integral(n,100));
1087  // htruesignal_calibdown->SetBinContent(n, truesignal_calibdown->Integral(n,100));
1088 
1089  htruesignal_xsecup->SetBinContent(n, truesignal_xsecup->Integral(n,100));
1090  htruesignal_xsecdown->SetBinContent(n, truesignal_xsecdown->Integral(n,100));
1091 
1092  // htruesignal_lightup->SetBinContent(n, truesignal_lightup->Integral(n,100));
1093  // htruesignal_lightdown->SetBinContent(n, truesignal_lightdown->Integral(n,100));
1094 
1095  htruesignal_fluxup->SetBinContent(n, truesignal_fluxup->Integral(n,100));
1096  htruesignal_fluxdown->SetBinContent(n, truesignal_fluxdown->Integral(n,100));
1097 
1098  // htruesignal_calib_shape->SetBinContent(n, truesignal_calib_shape->Integral(n,100));
1099  // htruesignal_ckv->SetBinContent(n, truesignal_ckv->Integral(n,100));
1100 
1101  }
1102 
1103 
1104  std::vector<TH1*> vector_up_truesignal_binbybin = {htruesignal_xsecup,htruesignal_fluxup};//{htruesignal_fluxup, htruesignal_xsecup, htruesignal_calibup, htruesignal_lightup, htruesignal_ckv, htruesignal_calib_shape};
1105 
1106  std::vector<TH1*> vector_dwn_truesignal_binbybin = {htruesignal_xsecdown,htruesignal_fluxdown};//{htruesignal_fluxdown, htruesignal_xsecdown, htruesignal_calibdown, htruesignal_lightdown, htruesignal_ckv, htruesignal_calib_shape};
1107 
1108 
1109 
1110 
1111  /*-----------------------------------------
1112  Plotting true signal with systematics
1113  -----------------------------------------*/
1114  TCanvas* ctruesignal = new TCanvas("ctruesignal","ctruesignal");
1115  //ctruesignal->SetLeftMargin(0.15);
1116 
1117  //htruesignal->SetTitle("Prong CVN Gamma ID Cut Value;True signal Events/8.09 #times 10^{20} POT");
1118  htruesignal->GetYaxis()->SetTitle("True signal Events/8.09 #times 10^{20} POT");
1119  htruesignal->SetTitleOffset(0.75,"Y");
1120 
1121  htruesignal->GetXaxis()->SetTitle("BDTG Cut Value");
1122  PlotWithSystErrorBandWidth(htruesignal,
1123  vector_up_truesignal_binbybin,
1124  vector_dwn_truesignal_binbybin, -1, -1, 1.3f, true);
1125 
1126  NDSimulation();
1127 
1128  ctruesignal->SaveAs(outDir2+CurrentFile[7].name.c_str()+CurrentDate[0].date.c_str()+".pdf");
1129  //ctruesignal->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/Prong_CVN_Gamma_IDcut_WithErrorBand_NtrueSig_Cuts_0_3_high_Combined_OneKnob1_out_3_MaNCRES.pdf");
1130  //ctruesignal->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/Prong_CVN_Gamma_IDcut_WithErrorBand_NtrueSig_Cuts_0_3_high.pdf");
1131  // ctruesignal->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/Calibration/Prong_CVN_Gamma_IDcut_WithErrorBand_NtrueSig_Cuts_calib_0_3_high.pdf");
1132  // ctruesignal->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/Calibration_shape/Prong_CVN_Gamma_IDcut_WithErrorBand_NtrueSig_Cuts_calib_shape_0_3_high.pdf");
1133  // ctruesignal->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/Cherenkov/Prong_CVN_Gamma_IDcut_WithErrorBand_NtrueSig_Cuts_Cherenkov_0_3_high.pdf");
1134  // ctruesignal->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/GENIE/Prong_CVN_Gamma_IDcut_WithErrorBand_NtrueSig_Cuts_GENIE_0_3_high.pdf");
1135  // ctruesignal->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/light_level/Prong_CVN_Gamma_IDcut_WithErrorBand_NtrueSig_Cuts_light_level_0_3_high.pdf");
1136  // ctruesignal->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/ppfx/Prong_CVN_Gamma_IDcut_WithErrorBand_NtrueSig_Cuts_ppfx_0_3_high.pdf");
1137  //****************************************************************************************************************************************************************************************////
1138 
1139  //denominator [3]---> 1vtx , denomnator [3]----> 1vtx + containment
1140 
1141  TH1* num_xsecdown = xsecdown[1]->ToTH1(pot); ///comment this if just caliip or down is needed for plots
1142  TH1* num_xsecup = xsecup[1]->ToTH1(pot); ///comment this if just caliip or down is needed for plots
1143 
1144  TH1* denom_xsecdown = xsecdown[2]->ToTH1(pot); ///comment this if just caliip or down is needed for plots
1145  TH1* denom_xsecup = xsecup[2]->ToTH1(pot); ///comment this if just caliip or down is needed for plots
1146 
1147  TH1* num_fluxdown = fluxdown[1]->ToTH1(pot); ///comment this if just caliip or down is needed for plots
1148  TH1* num_fluxup = fluxup[1]->ToTH1(pot); ///comment this if just caliip or down is needed for plots
1149 
1150  TH1* denom_fluxdown = fluxdown[2]->ToTH1(pot); ///comment this if just caliip or down is needed for plots
1151  TH1* denom_fluxup = fluxup[2]->ToTH1(pot); ///comment this if just caliip or down is needed for plots
1152 
1153  // TH1* num_calib_shape = calib_shape[1]->ToTH1(pot); ///comment this if just caliip or down is needed for plots
1154  // TH1* denom_calib_shape = calib_shape[2]->ToTH1(pot); ///comment this if just caliip or down is needed for plots
1155 
1156  // TH1* num_ckv = ckv[1]->ToTH1(pot); ///comment this if just caliip or down is needed for plots
1157  // TH1* denom_ckv = ckv[2]->ToTH1(pot); ///comment this if just caliip or down is needed for plots
1158 
1159  // TH1* num_lightup = lightup[1]->ToTH1(pot); ///comment this if just caliip or down is needed for plots
1160  // TH1* num_lightdown = lightdown[1]->ToTH1(pot); ///comment this if just caliip or down is needed for plots
1161 
1162  // TH1* denom_lightup = lightup[2]->ToTH1(pot); ///comment this if just caliip or down is needed for plots
1163  // TH1* denom_lightdown = lightdown[2]->ToTH1(pot); ///comment this if just caliip or down is needed for plots
1164 
1165  // TH1* num_calibdown = calibdown[1]->ToTH1(pot);
1166  // TH1* num_calibup = calibup[1]->ToTH1(pot);
1167 
1168  // TH1* denom_calibdown = calibdown[2]->ToTH1(pot);
1169  // TH1* denom_calibup = calibup[2]->ToTH1(pot);
1170 
1171 
1172 
1173  TH1* sig = nominal[1]->ToTH1(pot); // Signal Presel : kVtxInFiducial && kContainment && k2Prongs && kPi0RemIDND && kNCurrent && kPi0 && kNpCVN1gammaID (Nominal)
1174  TH1* denome = nominal[2]->ToTH1(pot); // ///true signal fiducial (Nsig,fid only): k2Prongs && kVtxInFiducial && kNCurrent && kPi0 (0.3geV);(Nominal)
1175 
1176  float totsig = denome->Integral(1, 100);
1177  float vdenom_xsecdown = denom_xsecdown->Integral(1, 100); ///comment this if just caliip or down is needed for plots
1178  float vdenom_xsecup = denom_xsecup->Integral(1, 100); ///comment this if just caliip or down is needed for plots
1179  float vdenom_fluxdown = denom_fluxdown->Integral(1, 100); ///comment this if just caliip or down is needed for plots
1180  float vdenom_fluxup = denom_fluxup->Integral(1, 100); // comment this if just caliip or down is needed for plots
1181  // float vdenom_calibdown = denom_calibdown->Integral(1, 100);
1182  // float vdenom_calibup = denom_calibup->Integral(1, 100);
1183  // float vdenom_lightdown = denom_lightdown->Integral(1, 100); ///comment this if just caliip or down is needed for plots
1184  // float vdenom_lightup = denom_lightup->Integral(1, 100); ///comment this if just caliip or down is needed for plots
1185  // float vdenom_calib_shape = denom_calib_shape->Integral(1, 100); ///comment this if just caliip or down is needed for plots
1186  // float vdenom_ckv = denom_ckv->Integral(1, 100); ///comment this if just caliip or down is needed for plots
1187 
1188 std::cout << "******************T=Sig: "<< sig <<std::endl; ////////////////////////////////////////////////////////////////////////
1189  std::cout << "******************T=denome: "<< denome <<std::endl; ////////////////////////////////////////////////////////////////////////
1190 
1191  gPad->RedrawAxis();
1192 
1193 //std::cout << totsig << std::endl;
1194 
1195 TH1* nomeff = new TH1F("eff_",
1196  "; BDTG Cut;Selection Efficiency",
1197  100,-1,1);
1198 
1199 TH1* eff_calibdown = new TH1F("eff_calibdown",
1200  "Efficiency vs Cut Value; BDTG ;Efficiency",
1201  100,-1,1);
1202 
1203 TH1* eff_calibup = new TH1F("eff_calibup",
1204  "Efficiency vs Cut Value;BDTG ;Efficiency",
1205  100,-1,1);
1206 
1207 TH1* eff_lightdown = new TH1F("eff_lightdown",
1208  "Efficiency vs Cut Value; BDTG ;Efficiency",
1209  100,-1,1);
1210 
1211 TH1* eff_lightup = new TH1F("eff_lightup",
1212  "Efficiency vs Cut Value; BDTG ;Efficiency",
1213  100,-1,1 );
1214 
1215 TH1* eff_ckv = new TH1F("eff_ckv",
1216  "Efficiency vs Cut Value;BDTG ;Efficiency",
1217  100,-1,1);
1218 
1219 TH1* eff_calib_shape = new TH1F("eff_calibshape",
1220  "Efficiency vs Cut Value; BDTG ;Efficiency",
1221  100,-1,1);
1222 
1223 TH1* eff_fluxdown = new TH1F("eff_fluxdown",
1224  "Efficiency vs Cut Value; BDTG ;Efficiency",
1225  100,-1,1);
1226 
1227 TH1* eff_fluxup = new TH1F("eff_fluxup",
1228  "Efficiency vs Cut Value; BDTG ;Efficiency",
1229  100,-1,1);
1230 
1231 TH1* eff_xsecdown = new TH1F("eff_xsecdown",
1232  "Efficiency vs Cut Value; BDTG ;Efficiency",
1233  100,-1,1);
1234 
1235 TH1* eff_xsecup = new TH1F("eff_xsecup",
1236  "Efficiency vs Cut Value; BDTG ;Efficiency",
1237  100,-1,1);
1238 
1239 
1240 
1241 
1242 
1243 
1244 
1245  std::cout << "******************Total selected events (sta+sys): "<< totsig <<std::endl; ////////////////////////////////////////////////////////////////////////
1246 
1247 for(int n = 1; n < hnom_Nbkg_Cuts->GetNbinsX() +1; ++n){
1248 
1249  const float n1sig = sig->Integral(n, 100);
1250  const float nxsecdown = num_xsecdown->Integral(n, 100)/vdenom_xsecdown; ///comment this if just caliip or down is needed for plots
1251  const float nxsecup = num_xsecup->Integral(n, 100)/vdenom_xsecup; ///comment this if just caliip or down is needed for plots
1252  const float nfluxdown = num_fluxdown->Integral(n, 100)/vdenom_fluxdown;
1253  const float nfluxup = num_fluxup->Integral(n, 100)/vdenom_fluxup; ///comment this if just caliip or down is needed for plots
1254  // const float ncalibdown = num_calibdown->Integral(n, 100)/vdenom_calibdown;
1255  // const float ncalibup = num_calibup->Integral(n, 100)/vdenom_calibup;
1256  // const float nlightdown = num_lightdown->Integral(n, 100)/vdenom_lightdown; ///comment this if just caliip or down is needed for plots
1257  // const float nlightup = num_lightup->Integral(n, 100)/vdenom_lightup; ///comment this if just caliip or down is needed for plots
1258  // const float ncalib_shape = num_calib_shape->Integral(n, 100)/vdenom_calib_shape; ///comment this if just caliip or down is needed for plots
1259  // const float nckv = num_ckv->Integral(n, 100)/vdenom_ckv;
1260  ///comment this if just caliip or down is needed for plots
1261  std::cout << "******************n1sig: "<< n1sig <<std::endl;
1262  std::cout << "******************Total selected events (sta+sys): "<< totsig <<std::endl;
1263 
1264 
1265  float neff = n1sig/(totsig); //Signal Presel/True vtx fiducial only signal
1266  std::cout << "******************neff: "<< neff <<std::endl; /////////////////////////////////////////////////////////////////////////
1267 
1268 
1269  nomeff->SetBinContent(n, neff);
1270 
1271  std::cout << "******************nomeff: "<< nomeff <<std::endl; //////////////////////////////////////////////////////////////////////
1272  eff_xsecdown->SetBinContent(n, nxsecdown);
1273 
1274  //std::cout << "******************eff_xsecdown: "<< eff_xsecdown<<std::endl; //////////////////////////////////////////////////////////////////////
1275 
1276  eff_xsecup->SetBinContent(n, nxsecup); ///comment this if just caliip or down is needed for plots
1277  eff_fluxdown->SetBinContent(n, nfluxdown); ///comment this if just caliip or down is needed for plots
1278  eff_fluxup->SetBinContent(n, nfluxup); // comment this if just caliip or down is needed for plots
1279  // eff_calibdown->SetBinContent(n, ncalibdown);
1280  // eff_calibup->SetBinContent(n, ncalibup);
1281  // eff_lightdown->SetBinContent(n, nlightdown); ///comment this if just caliip or down is needed for plots
1282  // eff_lightup->SetBinContent(n, nlightup); ///comment this if just caliip or down is needed for plots
1283  // eff_calib_shape->SetBinContent(n, ncalib_shape); ///comment this if just caliip or down is needed for plots
1284  // eff_ckv->SetBinContent(n, nckv);
1285  ///comment this if just caliip or down is needed for plots
1286  }
1287 
1288  std::vector<TH1*> rdn = {eff_xsecdown,eff_fluxdown};// {eff_calibdown, eff_lightdown, eff_ckv, eff_calib_shape, eff_fluxdown, eff_xsecdown}; ///comment lightup[3]->ToTH1(pot), ckv[3]->ToTH1(pot), calib_shape[3]->ToTH1(pot), fluxup[3]->ToTH1(pot), xsecup[3]->ToTH1(pot) if plotting just calib
1289 
1290  std::vector<TH1*> rup = {eff_xsecup,eff_fluxup};//{eff_calibup, eff_lightup, eff_ckv, eff_calib_shape, eff_fluxup, eff_xsecup}; ///comment lightup[3]->ToTH1(pot), ckv[3]->ToTH1(pot), calib_shape[3]->ToTH1(pot), fluxup[3]->ToTH1(pot), xsecup[3]->ToTH1(pot) if plotting just calib
1291 
1292 
1293 
1294  PlotDeltaSigmaWithSigma(hNbkg_Cuts,hNSel_Cuts,
1295  vector_up_Nbkg_Cuts_binbybin,
1296  vector_dwn_Nbkg_Cuts_binbybin,
1297  vector_up_NSel_Cuts_binbybin,
1298  vector_dwn_NSel_Cuts_binbybin,
1299  nomeff, rup, rdn,
1300  "BDTG_sigma_Systematics_optimization_Myfid_cont","BDTG_Nbkg_Cuts_syst_Systematics_optimization_Myfid_cont","BDTG_Nbkg_Cuts_stat__Systematics_optimization_Myfid_cont", "BDTG_NSel_Cuts_syst__Systematics_optimization_Myfid_cont","BDTG_NSel_Cuts_stat_Systematics_optimization_Myfid_cont","BDTG_delta_eff__Systematics_optimization_Myfid_cont","c7__Systematics_BDTG_optimization_Myfid_cont","c8__Systematics_BDTG_optimization_Myfid_cont","c9__Systematics_BDTG_optimization_Myfid_cont", true);
1301 
1302 
1303  TCanvas* ceff = new TCanvas("ceff", "ceff");
1304  ceff->SetLeftMargin(0.15);
1305  PlotWithSystErrorBandWidth(nomeff, rup, rdn, -1, -1, 1.3f, true);
1306 
1307  NDSimulation();
1308 
1309  ceff->SaveAs(outDir2+CurrentFile[8].name.c_str()+CurrentDate[0].date.c_str()+".pdf");
1310  //ceff->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/ceff03GEV_Combined_OneKnob1_out_3_MaNCRES.pdf");
1311  //ceff->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/ceff03GEV.pdf");
1312  // ceff->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/Calibration/Prong_CVN_Gamma_ID_eff_cutvalue_calib_0_3s_high.pdf");
1313  // ceff->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/Calibration_shape/Prong_CVN_Gamma_ID_eff_cutvalue_calib_shape_0_3_high.pdf");
1314  // ceff->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/Cherenkov/Prong_CVN_Gamma_ID_eff_cutvalue_Cherenkov_0_3_high.pdf");
1315  // ceff->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/GENIE/Prong_CVN_Gamma_ID_eff_cutvalue_GENIE_0_3_high.pdf");
1316  // ceff->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/light_level/Prong_CVN_Gamma_ID_eff_cutvalue_light_level_0_3_high.pdf");
1317  // ceff->SaveAs("/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/plots/2prong03high/ppfx/Prong_CVN_Gamma_ID_eff_cutvalue_ppfx_0_3_high.pdf");
1318 
1319  }// end loop over vars
1320 
1321 }
fin
TString fin
Definition: Style.C:24
name
const XML_Char * name
Definition: expat.h:151
plot_validation_datamc.fname
fname
Definition: plot_validation_datamc.py:180
GetFluxError
void GetFluxError(const Spectrum *nom, std::vector< ana::Spectrum * > univs, std::vector< ana::Spectrum * > &supdown)
Definition: plot_uncertainty.C:428
kRed
enum BeamMode kRed
NDSimulation
void NDSimulation()
Definition: plot_uncertainty.C:41
std::sqrt
T sqrt(T number)
Definition: d0nt_math.hpp:156
CurrentDate
const DateUp CurrentDate[Date]
Definition: plot_uncertainty.C:77
PlotWithSystErrorBandWidth
void PlotWithSystErrorBandWidth(TH1 *&nom, std::vector< TH1 * > &ups, std::vector< TH1 * > &dns, int col, int errCol, float headroom, bool newaxis)
Definition: plot_uncertainty.C:82
ana::LoadFrom
std::unique_ptr< T > LoadFrom(TDirectory *dir, const std::string &label)
Definition: LoadFromFile.h:17
cut_labels
std::vector< std::string > cut_labels
Definition: preselection_cutflow.C:22
CurrentFile
const FileUp CurrentFile[File]
Definition: plot_uncertainty.C:62
PandAna.Demos.pi0_spectra.pur
tuple pur
Definition: pi0_spectra.py:135
PlotDeltaSigmaWithSigma
void PlotDeltaSigmaWithSigma(TH1 *&nom, TH1 *&nomsel, std::vector< TH1 * > &ups, std::vector< TH1 * > &dns, std::vector< TH1 * > &upssel, std::vector< TH1 * > &dnssel, TH1 *&heff, std::vector< TH1 * > &effups, std::vector< TH1 * > &effdns, TString outname1, TString outname2, TString outname3, TString outname4, TString outname5, TString outname6, TString outname7, TString outname8, TString outname9, Bool_t printme=false)
Definition: plot_uncertainty.C:140
DateUp::date
std::string date
Definition: plot_uncertainty.C:58
om::cout
OStream cout
Definition: OStream.cxx:6
outDir2
TString outDir2
Definition: plot_uncertainty.C:53
pot
const double pot
Definition: plot_uncertainty.C:37
art::to_string
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
string
enum BeamMode string
void PlotDeltaSigmaWithSigma ( TH1 *&  nom,
TH1 *&  nomsel,
std::vector< TH1 * > &  ups,
std::vector< TH1 * > &  dns,
std::vector< TH1 * > &  upssel,
std::vector< TH1 * > &  dnssel,
TH1 *&  heff,
std::vector< TH1 * > &  effups,
std::vector< TH1 * > &  effdns,
TString  outname1,
TString  outname2,
TString  outname3,
TString  outname4,
TString  outname5,
TString  outname6,
TString  outname7,
TString  outname8,
TString  outname9,
Bool_t  printme = false 
)

uncertainty on Selected events

sigma

bkg syst

bkg Nsel

Numerator bkg Nsel

Denominator bkg Nsel

Sig Nsel

Definition at line 140 of file plot_uncertainty.C.

References demo5::c1, demo5::c2, chisquared::c3, chisquared::c4, om::cout, allTimeWatchdog::endl, kRed, NDSimulation(), outDir, std::sqrt(), and submit_syst::y.

Referenced by plot_uncertainty().

153 {
154 
155  double yMax = nom->GetBinContent(nom->GetMaximumBin()); // location of bin with max value
156 
157 
158 
159 
160 //****************************Total Background **********************************************************///
161 
162  TH1F* hnbkg_nsel = new TH1F("hnbkg_nsel","; BDTG Cut; #frac{#sqrt{N_{bkg}}}{N_{sel}-N_{bkg}}", ////change variable name hnbkg_stat
163  nom->GetNbinsX(), nom->GetXaxis()->GetBinLowEdge(1),
164  nom->GetXaxis()->GetBinUpEdge(nom->GetNbinsX()));
165  TH1F* hnbkg_nsel_num = new TH1F("hnbkg_nsel_num","; BDTG Cut value; #sqrt{N_{bkg}}",
166  nom->GetNbinsX(), nom->GetXaxis()->GetBinLowEdge(1),
167  nom->GetXaxis()->GetBinUpEdge(nom->GetNbinsX()));
168 
169  TH1F* hnbkg_nsel_den = new TH1F("hnbkg_nsel_den",";BDTG Cut value; N_{sel}-N_{bkg}",
170  nom->GetNbinsX(), nom->GetXaxis()->GetBinLowEdge(1),
171  nom->GetXaxis()->GetBinUpEdge(nom->GetNbinsX()));
172 
173 
174  TH1F* hnbkg_syst = new TH1F("hnbkg_syst",";BDTG Cut; #frac{#deltaN^{syst}_{bkg}}{N_{sel}-N_{bkg}}",
175  nom->GetNbinsX(), nom->GetXaxis()->GetBinLowEdge(1),
176  nom->GetXaxis()->GetBinUpEdge(nom->GetNbinsX()));
177  //****************************Total Preselection **********************************************************///
178 
179  TH1F* hnsel_nsel = new TH1F("hnsel_nsel",";BDTG Cut; #frac{#sqrt{N_{sel}}}{N_{sel}-N_{bkg}}", //PreSel
180  nomsel->GetNbinsX(), nomsel->GetXaxis()->GetBinLowEdge(1),
181  nomsel->GetXaxis()->GetBinUpEdge(nom->GetNbinsX()));
182 
183  TH1F* hnsel_nsel_num = new TH1F("hnsel_nsel_num",";BDTG Cut value; #sqrt{N_{sel}}",
184  nomsel->GetNbinsX(), nomsel->GetXaxis()->GetBinLowEdge(1),
185  nomsel->GetXaxis()->GetBinUpEdge(nomsel->GetNbinsX()));
186 
187  TH1F* hnsel_syst = new TH1F("hnsel_syst",";BDTG Cut; #frac{#deltaN^{syst}_{sig}}{N_{sel}-N_{bkg}}",
188  nomsel->GetNbinsX(), nomsel->GetXaxis()->GetBinLowEdge(1),
189  nomsel->GetXaxis()->GetBinUpEdge(nom->GetNbinsX()));
190 
191  TH1F* hepsilon = new TH1F("hepsilon",";BDTG Cut; #frac{#delta#varepsilon}{#varepsilon}",
192  nom->GetNbinsX(), nom->GetXaxis()->GetBinLowEdge(1),
193  nom->GetXaxis()->GetBinUpEdge(nom->GetNbinsX()));
194 
195  TH1F* hsigma = new TH1F("hsigma",";BDTG Cut; #frac{#delta#sigma}{#sigma}",
196  nom->GetNbinsX(), nom->GetXaxis()->GetBinLowEdge(1),
197  nom->GetXaxis()->GetBinUpEdge(nom->GetNbinsX()));
198 
199 
200 
201  // for(int binIdx = 0; binIdx < 51; ++binIdx){
202  //for(int binIdx = 0; binIdx < nom->GetNbinsX()<51; ++binIdx){
203 
204 for(int binIdx = 1; binIdx < nom->GetNbinsX()+1; ++binIdx){
205 
206  const double y = nom ->GetBinContent(binIdx); //Nbkg events
207  const double ysel = nomsel ->GetBinContent(binIdx); //Nsel events
208  const double Neff = heff->GetBinContent(binIdx); // efficiency
209 
210 
211 
212  double deltaNbkg2 = 0;
213  double deltaNsig2 = 0;
214  double deltaeff2 = 0;
215  double errUp = 0;
216  double errDn = 0;
217  double errUpsig = 0;
218  double errDnsig = 0;
219  double delta_eff = 0;
220  double totalbkg = 0;
221  double totalpresel = 0;
222 
223 
224 
225 
226 
227  for(unsigned int systIdx = 0; systIdx < ups.size(); ++systIdx){
228  double deltaNbkg = (ups[systIdx]->GetBinContent(binIdx)-dns[systIdx]->GetBinContent(binIdx))/2;
229  double deltaNsel = (upssel[systIdx]->GetBinContent(binIdx)-dnssel[systIdx]->GetBinContent(binIdx))/2;
230  double deltaeff = (effups[systIdx]->GetBinContent(binIdx)-effdns[systIdx]->GetBinContent(binIdx))/2;
231 
232 
233  deltaeff2 += deltaeff*deltaeff; //uncertainty on efficiency
234 
235  deltaNbkg2 += deltaNbkg*deltaNbkg; //uncertainty on the background
236 
237  deltaNsig2 += deltaNsel*deltaNsel; //uncertainty on the Selection
238  delta_eff+=deltaeff;
239 
240  } // end for systIdx
241  const double s = 0.23;
242  totalbkg = y*s + deltaNbkg2; //Nbkg+(systematic uncertainty in the background)^2
243  totalpresel = ysel; //Nsel
244  double sel_bkg = ysel - y; //Nsel-Nbkg
245  double sel_bkg2 = sel_bkg*sel_bkg; //(Nsel-Nbkg)^2
246  double uncbkg = deltaNbkg2/sel_bkg2;//(Nbkg+(syst.unc.bkg)^2)/(Nsel-Nbkg)^2
247  // double uncbkg = totalbkg/sel_bkg2;
248  double uncsel = totalpresel/sel_bkg2;//Nsel/(Nsel-Nbkg)^2
249  double Neff2 = Neff*Neff; //efficiency^2
250  double delta_eff_new = delta_eff;
251  double unceff = deltaeff2/Neff2; //uncertainty on selection efficiency
252  double uncsigma2 = uncbkg + unceff + uncsel ;
253 
254  if (binIdx==1){
255  cout << "#########Bin number value for 0.90: "<< nom->GetXaxis()->FindBin(0.90);
256  cout << "#########Bin number value for 0.91: "<< nom->GetXaxis()->FindBin(0.91);
257  cout << "#########Bin number value for 0.92: "<< nom->GetXaxis()->FindBin(0.92);
258  cout << "#########Bin number value for 0.93: "<< nom->GetXaxis()->FindBin(0.93);
259  cout << "#########Bin number value for 0.94: "<< nom->GetXaxis()->FindBin(0.94);
260  cout << "#########Bin number value for 0.95: "<< nom->GetXaxis()->FindBin(0.95);
261  cout << "#########Bin number value for 0.96: "<< nom->GetXaxis()->FindBin(0.96);
262  cout << "#########Bin number value for 0.97: "<< nom->GetXaxis()->FindBin(0.97);
263  cout << "#########Bin number value for 0.98: "<< nom->GetXaxis()->FindBin(0.98);
264  cout << "#########Bin number value for 0.99: "<< nom->GetXaxis()->FindBin(0.99);
265  cout << "#########Bin number value for 1: "<< nom->GetXaxis()->FindBin(1);
266 
267  cout << "*********************** bin number**********************: "<< binIdx <<endl;
268  cout << "Bkg number (from bin content)"<< y <<endl;
269  cout << "Total bkg (sta+sys): "<< totalbkg <<endl;
270  cout << "Selected Sig(denominator): "<< sel_bkg <<endl;
271  cout << "Selected Sig square(denominator): "<< sel_bkg2 <<endl;
272  cout << "Stat unc in bkg "<< sqrt(y)/sqrt(sel_bkg2) <<endl;
273  cout << "Sys unc in bkg "<< sqrt(deltaNbkg2)/sqrt(sel_bkg2) <<endl;
274  cout << "Sys unc in bkg - the numerator" << sqrt(deltaNbkg2)<<endl;
275  cout << "Total unc in bkg "<< uncbkg <<endl;
276 
277  cout << "selected number (from bin content): "<< ysel <<endl;
278  cout << "Total selected events (sta+sys): "<< totalpresel <<endl;
279  cout << "Selected Sig(denominator): "<< sel_bkg <<endl;
280  cout << "Selected Sig square(denominator): "<< sel_bkg2 <<endl;
281  cout << "Stat unc in selected "<< sqrt(ysel)/sqrt(sel_bkg2) <<endl;
282  cout << "Sys unc in selected "<< sqrt(deltaNsig2)/sqrt(sel_bkg2) <<endl;
283  cout << "Total unc in selected "<< uncsel <<endl;
284  cout << "Efficiency bin content (from bin content): "<< Neff <<endl;
285  cout << "uncertainty efficiency bin content (from bin content): "<< delta_eff_new<<endl;
286 
287  cout << "Square of Efficiency (from bin content): "<< Neff2 <<endl;
288  cout << "unc Efficiency : "<< unceff <<endl;
289  cout << "sigma uncertainty" << (uncsigma2)<<endl;
290  cout << "square root sigma uncertainty" << sqrt(uncsigma2)<<endl;
291  }
292  hnbkg_syst->Fill(nom->GetXaxis()->GetBinCenter(binIdx), sqrt(deltaNbkg2)/sqrt(sel_bkg2)); //background systematics
293  hnbkg_nsel->Fill(nom->GetXaxis()->GetBinCenter(binIdx), sqrt(y)/sqrt(sel_bkg2)); //uncertainty on the background events
294  hnbkg_nsel_num->Fill(nom->GetXaxis()->GetBinCenter(binIdx), sqrt(y)); //(numerator) Nbkg
295  hnbkg_nsel_den->Fill(nom->GetXaxis()->GetBinCenter(binIdx), sel_bkg);//(denominator) Nsel-Nbkg
296 
297  hnsel_syst->Fill(nomsel->GetXaxis()->GetBinCenter(binIdx), sqrt(deltaNsig2)/sqrt(sel_bkg2));//selected events systematics
298  hnsel_nsel_num->Fill(nom->GetXaxis()->GetBinCenter(binIdx), sqrt(ysel));//numerator Nsel
299  hnsel_nsel->Fill(nom->GetXaxis()->GetBinCenter(binIdx), sqrt(ysel)/sqrt(sel_bkg2));///uncertainty on Selected events
300 
301 
302  hepsilon->Fill(nom->GetXaxis()->GetBinCenter(binIdx), sqrt(deltaeff2)/sqrt(Neff2));
303  hsigma->Fill(nom->GetXaxis()->GetBinCenter(binIdx), sqrt(uncsigma2));
304 
305  } // end for i
306 
307 
308  hnbkg_syst->SetLineColor(kRed);
309  hnbkg_nsel->SetLineColor(kRed);
310  hnbkg_nsel_num->SetLineColor(kRed);
311  hnbkg_nsel_den->SetLineColor(kRed);
312 
313  hnsel_syst->SetLineColor(kRed);
314  hnsel_nsel->SetLineColor(kRed);
315  hnsel_nsel_num->SetLineColor(kRed);
316 
317  hepsilon->SetLineColor(kRed);
318  hsigma->SetLineColor(kRed);
319 
320  hnbkg_syst->GetXaxis()->CenterTitle();
321  hnbkg_syst->GetYaxis()->CenterTitle();
322  hnbkg_nsel->GetXaxis()->CenterTitle();
323  hnbkg_nsel->GetYaxis()->CenterTitle();
324  hnbkg_nsel_num->GetXaxis()->CenterTitle();
325  hnbkg_nsel_num->GetYaxis()->CenterTitle();
326  hnbkg_nsel_den->GetXaxis()->CenterTitle();
327  hnbkg_nsel_den->GetYaxis()->CenterTitle();
328 
329  hnsel_syst->GetXaxis()->CenterTitle();
330  hnsel_syst->GetYaxis()->CenterTitle();
331  hnsel_nsel->GetXaxis()->CenterTitle();
332  hnsel_nsel->GetYaxis()->CenterTitle();
333  hnsel_nsel_num->GetXaxis()->CenterTitle();
334  hnsel_nsel_num->GetYaxis()->CenterTitle();
335 
336  hepsilon->GetXaxis()->CenterTitle();
337  hepsilon->GetYaxis()->CenterTitle();
338  hsigma->GetXaxis()->CenterTitle();
339  hsigma->GetYaxis()->CenterTitle();
340  gPad->RedrawAxis();
341  //nom->Draw("hist ][ same");
342  TCanvas *c1 = new TCanvas("c1","c1"); ///sigma
343  c1->SetLeftMargin(0.15);
344  hsigma->Draw("hist ][");
345 
346  NDSimulation();
347 
348  TCanvas *c2 = new TCanvas("c2","c2"); ///bkg syst
349  c2->SetLeftMargin(0.16);
350  hnbkg_syst->Draw("hist ][");
351  hnbkg_syst->GetYaxis()->SetTitleOffset(1.2);
352 
353  NDSimulation();
354 
355  TCanvas *c2a = new TCanvas("c2a","c2a"); ///bkg Nsel
356  c2a->SetLeftMargin(0.15);
357  hnbkg_nsel->Draw("hist ][");
358  hnbkg_nsel->GetYaxis()->SetTitleOffset(1.1);
359 
360  NDSimulation();
361 
362 
363  TCanvas *c2b = new TCanvas("c2b","c2b"); ///Numerator bkg Nsel
364  c2b->SetLeftMargin(0.15);
365  hnbkg_nsel_num->Draw("hist ][");
366  hnbkg_nsel_num->GetYaxis()->SetTitleOffset(1.1);
367 
368  NDSimulation();
369 
370  TCanvas *c2c = new TCanvas("c2c","c2c"); ///Denominator bkg Nsel
371  c2c->SetLeftMargin(0.15);
372  hnbkg_nsel_den->Draw("hist ][");
373  hnbkg_nsel_den->GetYaxis()->SetTitleOffset(1.1);
374 
375  NDSimulation();
376 
377  TCanvas *c3 = new TCanvas("c3","c3"); ////Sig syst
378  c3->SetLeftMargin(0.16);
379  hnsel_syst->Draw("hist ][");
380  hnsel_syst->GetYaxis()->SetTitleOffset(1.2);
381  NDSimulation();
382 
383  TCanvas *c3a = new TCanvas("c3a","c3a"); ///Sig Nsel
384  c3a->SetLeftMargin(0.15);
385  hnsel_nsel->Draw("hist ][");
386  hnsel_nsel->GetYaxis()->SetTitleOffset(1.1);
387 
388  NDSimulation();
389 
390  TCanvas *c3b = new TCanvas("c3b","c3b");
391  c3b->SetLeftMargin(0.15);
392  hnsel_nsel_num->Draw("hist ][");
393  hnsel_nsel_num->GetYaxis()->SetTitleOffset(1.1);
394 
395  NDSimulation();
396 
397 
398  TCanvas *c4 = new TCanvas("c4", "c4");
399  c4->SetLeftMargin(0.15);
400  hepsilon->Draw("hist ][");
401 
402  NDSimulation();
403 
404 
405  if (printme){
406  TString outDir = "/nova/app/users/acedeno/tag_releasesS18-06-14/NDAna/ncpi0_semi_inc_png_cvn/NewDataSetSyst";
407  c1->SaveAs(outDir+"/"+outname1+CurrentDate[0].date.c_str()+".pdf");
408  c2->SaveAs(outDir+"/"+outname2+CurrentDate[0].date.c_str()+".pdf");
409  c2a->SaveAs(outDir+"/"+outname3+CurrentDate[0].date.c_str()+".pdf");
410  c2b->SaveAs(outDir+"/"+outname4+CurrentDate[0].date.c_str()+".pdf");
411 
412  c2c->SaveAs(outDir+"/"+outname5+CurrentDate[0].date.c_str()+".pdf");
413  c3->SaveAs(outDir+"/"+outname6+CurrentDate[0].date.c_str()+".pdf");
414  c3a->SaveAs(outDir+"/"+outname7+CurrentDate[0].date.c_str()+".pdf");
415  c3b->SaveAs(outDir+"/"+outname8+CurrentDate[0].date.c_str()+".pdf");
416  c4->SaveAs(outDir+"/"+outname9+CurrentDate[0].date.c_str()+".pdf");
417  // c6->SaveAs(outDir+"/"+outname6+".png");
418  // c7->SaveAs(outDir+"/"+outname7+".png");
419  // c8->SaveAs(outDir+"/"+outname8+".png");
420 
421  }
422 
423 
424 }//end of plotwith error band
xml_schema::date
::xsd::cxx::tree::date< char, simple_type > date
Definition: Database.h:186
kRed
enum BeamMode kRed
NDSimulation
void NDSimulation()
Definition: plot_uncertainty.C:41
std::sqrt
T sqrt(T number)
Definition: d0nt_math.hpp:156
CurrentDate
const DateUp CurrentDate[Date]
Definition: plot_uncertainty.C:77
outDir
std::string outDir
Definition: plot_ND_DataMC.C:35
submit_syst.y
y
Definition: submit_syst.py:249
demo5.c2
c2
Definition: demo5.py:33
s
const XML_Char * s
Definition: expat.h:262
om::cout
OStream cout
Definition: OStream.cxx:6
demo5.c1
c1
Definition: demo5.py:24
void PlotWithSystErrorBandWidth ( TH1 *&  nom,
std::vector< TH1 * > &  ups,
std::vector< TH1 * > &  dns,
int  col,
int  errCol,
float  headroom,
bool  newaxis 
)

Definition at line 82 of file plot_uncertainty.C.

References om::cout, allTimeWatchdog::endl, MECModelEnuComparisons::g, hi(), ana::kTotalMCColor, ana::kTotalMCErrorBandColor, lo(), std::sqrt(), std::swap(), w, and submit_syst::y.

Referenced by plot_uncertainty().

87 {
88 
89  if(col == -1){
90  col = kTotalMCColor;
91  errCol = kTotalMCErrorBandColor;
92  }
93  else if(errCol == -1) errCol = col-7; // hopefully a lighter version
94 
95  nom->SetLineColor(col);
96  nom->GetXaxis()->CenterTitle();
97  nom->GetYaxis()->CenterTitle();
98 
99  if(newaxis) nom->Draw("hist ]["); // Set the axes up
100 
101  double yMax = nom->GetBinContent(nom->GetMaximumBin());
102  cout << nom->GetMaximumBin() << "\t" << yMax << endl;
103  TGraphAsymmErrors* g = new TGraphAsymmErrors;
104 
105  for(int binIdx = 0; binIdx < nom->GetNbinsX()+2; ++binIdx){
106  const double y = nom->GetBinContent(binIdx);
107  g->SetPoint(binIdx, nom->GetXaxis()->GetBinCenter(binIdx), y);
108 
109  const double w = nom->GetXaxis()->GetBinWidth(binIdx);
110 
111  double errUp = 0;
112  double errDn = 0;
113  //std::cout << "systIdx \t" << " up \t"<< "y \t" <<"up -y\t"<<" up-y/y\t"<< "dn \t" <<"y-dn\t" << "y-dn/y\t" << "errup \t"<<"errdn \t" << endl;
114  for(unsigned int systIdx = 0; systIdx < ups.size(); ++systIdx){
115  double hii = (ups[systIdx]->GetBinContent(binIdx)-y)/y;
116  double loo = (y-dns[systIdx]->GetBinContent(binIdx))/y;
117  double hi = (ups[systIdx]->GetBinContent(binIdx)-y);
118  double lo = (y-dns[systIdx]->GetBinContent(binIdx));
119 
120  if(lo <= 0 && hi <= 0) std::swap(lo, hi);
121 
122  errUp += hi*hi;
123  errDn += lo*lo;
124  } // end for systIdx
125 
126  g->SetPointError(binIdx, w/2, w/2, sqrt(errDn), sqrt(errUp));
127  std::cout << "error up " << sqrt(errUp) << " and down " << sqrt(errDn) << endl;
128  } // end for i
129  g->SetFillColor(errCol);
130  g->Draw("e2 same");
131  g->GetYaxis()->SetRangeUser(0.00001, headroom*yMax);
132  nom->GetYaxis()->SetRangeUser(0.0001, headroom*yMax);
133  gPad->RedrawAxis();
134  nom->Draw("hist ][ same");
135 
136 
137 }//end of plotwith error band
lo
TSpline3 lo("lo", xlo, ylo, 12,"0")
std::sqrt
T sqrt(T number)
Definition: d0nt_math.hpp:156
ana::kTotalMCErrorBandColor
const Color_t kTotalMCErrorBandColor
Definition: Style.h:17
submit_syst.y
y
Definition: submit_syst.py:249
std::swap
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
Definition: HLTGlobalStatus.h:199
hi
TSpline3 hi("hi", xhi, yhi, 18,"0")
col
Int_t col[ntarg]
Definition: Style.C:29
om::cout
OStream cout
Definition: OStream.cxx:6
ana::kTotalMCColor
const Color_t kTotalMCColor
Definition: Style.h:16
w
Float_t w
Definition: plot.C:20

Variable Documentation

const DateUp CurrentDate[Date]
Initial value:
= {
{"_07-11-18_newpngcvn_png1_BDTGoptimization_NominalwFlux"},
}

Definition at line 77 of file plot_uncertainty.C.

Referenced by bdtstudyplotter(), plotHist_SpectrumCVNID(), resolutionplotter(), and twodvtxcontplotter().

const FileUp CurrentFile[File]
Initial value:
= {
{"/BDTG_Nbkg_Cuts_0_3_high_optimization_Myfid_cont"},
{"/BDTG_NSel_Cuts_0_3_high_optimization_Myfid_cont"},
{"/BDTG_NSig_Cuts_0_3_high_optimization_Myfid_cont"},
{"/BDTG_NtrueSig_Cuts_0_3_high_optimization_Myfid_cont"},
{"/BDTG_Cut_WithErrorBand_Nbkg_Cuts_0_3_high_optimization_Myfid_cont"},
{"/BDTGcut_WithErrorBand_NSel_Cuts_0_3_high_optimization_Myfid_cont"},
{"/BDTGcut_WithErrorBand_NSig_Cuts_0_3_high_optimization_Myfid_cont"},
{"/BDTGcut_WithErrorBand_NtrueSig_Cuts_0_3_high_optimization_Myfid_cont"},
{"/ceff03GEV_BDTG_optimization_Myfid_cont"},
{"/Purity_BDTG_Myfid_cont"},
{"/FOM_BDTG_Myfid_cont"},
}

Definition at line 62 of file plot_uncertainty.C.

Referenced by bdtstudyplotter(), plotHist_SpectrumCVNID(), resolutionplotter(), and twodvtxcontplotter().

const int Date = 1

Definition at line 76 of file plot_uncertainty.C.

Referenced by bdtstudyplotter(), plotHist_SpectrumCVNID(), resolutionplotter(), and twodvtxcontplotter().

const int File =13

Definition at line 61 of file plot_uncertainty.C.

Referenced by bdtstudyplotter(), FEBLoop(), hdf5::HDF5Maker::InitializeOutfile(), plotHist_SpectrumCVNID(), resolutionplotter(), and twodvtxcontplotter().

TString outDir2 = "/nova/app/users/acedeno/tag_releasesS18-06-14/NDAna/ncpi0_semi_inc_png_cvn/NewDataSetSyst"

Definition at line 53 of file plot_uncertainty.C.

Referenced by bdtstudyplotter(), plot_uncertainty(), plotHist_SpectrumCVNID(), resolutionplotter(), and twodvtxcontplotter().

const double pot = 8.09e20

Definition at line 37 of file plot_uncertainty.C.

Referenced by GetFluxError(), and plot_uncertainty().

Generated on Tue Apr 13 2021 14:15:56 for NOvA by   doxygen 1.8.11