plot_predictions.C File Reference

#include "includes.h"
#include "numu_tools.h"
#include "3FlavorAna/Plotting/NuePlotStyle.h"

Go to the source code of this file.

Functions
void	HornSettings (int quantile, double &pot, double &livetime, Sign::Sign_t &wrongs, double &maxy)
void	plot_predictions (std::string anaName="fhc", std::string calculator="2018calc", bool nh=true, bool fake=true)

Variables
const double	pot_fhc = kAna2018FHCPOT
const double	pot_rhc = kAna2018RHCPOT
const double	livetime_fhc = kAna2018FHCLivetime
const double	livetime_rhc = kAna2018RHCLivetime
double	pot = 0
double	pot_one = pot_rhc
double	pot_two = pot_fhc
double	livetime = 0
double	livetime_one = livetime_rhc
double	livetime_two = livetime_fhc
double	maxy = 0
double	maxy_one = 1.999
double	maxy_two = 3.999
Sign::Sign_t	wrongs = Sign::kNu
Sign::Sign_t	ws_one = Sign::kNu
Sign::Sign_t	ws_two = Sign::kAntiNu
std::string	extrapName = "extrap"
std::string	horn_one = "rhc"
std::string	horn_two = "fhc"
std::string	period_one = "full"
std::string	period_two = "full"
bool	twobeams = true
int	totquant = 8
bool	DianaIsLazy = true
bool	errorBand = true
Color_t	color_band = kViolet-9
Color_t	color_pred = kViolet+1
Color_t	color_ws = kBlue-7
Color_t	color_nc = kBlue-7
Color_t	color_data = kBlack
Color_t	color_cos = kGreen+1
Color_t	color_bkg = kGray+1
Style_t	line_ws = kSolid
Style_t	line_nc = kSolid
Style_t	line_data = kSolid
Style_t	line_pred = kSolid
Style_t	line_cos = kSolid
Style_t	line_bkg = kSolid
Style_t	marker_ws = kFullCircle
Style_t	marker_nc = kFullSquare
Style_t	marker_data = kFullCircle
Style_t	marker_pred = kFullCircle
Style_t	marker_cos = kFullTriangleUp
Style_t	marker_bkg = kFullCircle

Function Documentation

void HornSettings	(	int	quantile,
		double &	pot,
		double &	livetime,
		Sign::Sign_t &	wrongs,
		double &	maxy
	)

Definition at line 69 of file plot_predictions.C.

References livetime_one, livetime_two, maxy_one, maxy_two, pot_one, pot_two, ws_one, and ws_two.

                                                                                                {
  if(quantile<4){
    pot = pot_one;
    livetime = livetime_one;
    wrongs = ws_one;
    maxy = maxy_one;
  }
  else{
    pot = pot_two;
    livetime = livetime_two;
    wrongs = ws_two;
    maxy = maxy_two;
  }

}// end HornSettings

void plot_predictions	(	std::string	anaName = "fhc",
		std::string	calculator = "2018calc",
		bool	nh = true,
		bool	fake = true
	)

Definition at line 87 of file plot_predictions.C.

References std::asin(), calc, canvas(), ana::CenterTitles(), color_band, color_bkg, color_cos, color_data, color_nc, color_pred, color_ws, osc::_NoOscillations< T >::Copy(), CornerLabel(), make_associated_cosmic_defs::cosmics, om::cout, ana::DefaultOscCalc(), DianaIsLazy, drawQuantLabel(), allTimeWatchdog::endl, HTMLTools::entry(), errorBand, extrapName, FakeData(), ana::FillWithDimColor(), ana::getAllAna2018Systs(), make_mec_shifts_plots::GetMaximum(), horn_one, horn_two, HornSettings(), MECModelEnuComparisons::i, ana::Flavors::kAll, ana::Flavors::kAllNuMu, ana::Sign::kAntiNu, ana::Sign::kBoth, ana::Current::kCC, ana::Current::kNC, ana::kNuMu, ana::kNumuAna2018, make_mec_shifts_plots::legend, line_bkg, line_cos, line_data, line_nc, line_pred, line_ws, livetime, livetime_fhc, livetime_one, livetime_two, ana::MakeHistCanvasReady_Quant(), marker_bkg, marker_cos, marker_data, marker_nc, marker_pred, marker_ws, maxy, maxy_one, runNovaSAM::outName, pad1, pad2, pad3, period_one, ana::PimpHist(), ana::PlotWithSystErrorBand_Quant(), pot, pot_fhc, pot_one, pot_two, preliminary(), osc::_IOscCalcAdjustable< T >::SetDmsq32(), ana::SystShifts::SetShift(), osc::_IOscCalcAdjustable< T >::SetTh23(), ana::SplitCanvasQuant(), std::sqrt(), string, systs, totquant, twobeams, wrongs, and ws_one.

{
  

  bool preliminary = false;

  if(anaName=="rhc"){
    totquant = 4;
    twobeams = false;
  }
  // if fhc, change settings so the first 4 quantiles are fhc instead of rhc
  // dont do anything if rhc as the first 4 quantiles are rhc already
  if(anaName=="fhc"){
    horn_one = "fhc";
    period_one = "full";
    pot_one = pot_fhc;
    livetime_one = livetime_fhc;
    ws_one = Sign::kAntiNu;
    maxy_one = 3.999;
    twobeams = false;
    totquant = 4;
  }

  // calculator                                                                                                                                                                    
  std::string calcName = "2018calc";
  std::string hierarchy = "";
  osc::IOscCalcAdjustable* calc = DefaultOscCalc();
  if(nh){
    hierarchy = "nh";
    calc->SetTh23(asin(sqrt(0.558)));
    calc->SetDmsq32(0.00244);
  }
  if(!nh){
    hierarchy = "ih";
    calc->SetTh23(asin(sqrt(0.558)));
    calc->SetDmsq32(-0.00244);
  }
  osc::IOscCalcAdjustable* calcfake = DefaultOscCalc();
  calcfake->SetTh23(asin(sqrt(0.4)));
  calcfake->SetDmsq32(0.00244);
  osc::NoOscillations osc;
  auto calcNoOsc = osc.Copy();


  // print on screen /////////////////////////////////////////////////////////                                                                                                     
  std::cout << "\n\nmake table with predicted number of events\n" << std::endl;
  std::cout << "ana name: " << anaName << std::endl;
  std::cout << "calculator: " << calcName << ", hierarchy: " << hierarchy << std::endl;


  //inputs                                                                                                                                                                         
  TString dout_name      = "";
  std::string ddata_name = "/nova/ana/nu_mu_ana/Ana2018/Data/";
  std::string dcosm_name = "/nova/ana/nu_mu_ana/Ana2018/Cosmics/";
  std::string dpred_name = "/nova/ana/nu_mu_ana/Ana2018/Predictions/";
  std::string fpred_name = "";
  std::string fpred_suffix = "";


  // cosmics                                                                                                                                                                       
  std::string cosmics = "nocosmics";
  std::vector<Spectrum> s_cosmics;
  std::vector<TH1*> hCosmics;
  std::cout << "\nloading cosmics" << std::endl;
  std::cout << "" << horn_one << std::endl;
  std::string fcosmics_name_one = dcosm_name + "cosmics_" + horn_one + "__numu2018.root";
  TFile fcosmics_one(fcosmics_name_one.c_str());
  hCosmics.push_back((TH1*)fcosmics_one.Get("cosmics_q0")); // remember, quant0 = all                                                                                        
  hCosmics.push_back((TH1*)fcosmics_one.Get("cosmics_q1"));
  hCosmics.push_back((TH1*)fcosmics_one.Get("cosmics_q2"));
  hCosmics.push_back((TH1*)fcosmics_one.Get("cosmics_q3"));
  hCosmics.push_back((TH1*)fcosmics_one.Get("cosmics_q4")); 
  auto *cosmics_quant1_one = (TH1*)fcosmics_one.Get("cosmics_q1");
  auto *cosmics_quant2_one = (TH1*)fcosmics_one.Get("cosmics_q2");
  auto *cosmics_quant3_one = (TH1*)fcosmics_one.Get("cosmics_q3");
  auto *cosmics_quant4_one = (TH1*)fcosmics_one.Get("cosmics_q4");
  s_cosmics.push_back(Spectrum (cosmics_quant1_one, pot_one, livetime_one));
  s_cosmics.push_back(Spectrum (cosmics_quant2_one, pot_one, livetime_one));
  s_cosmics.push_back(Spectrum (cosmics_quant3_one, pot_one, livetime_one));
  s_cosmics.push_back(Spectrum (cosmics_quant4_one, pot_one, livetime_one));
  if(twobeams){
    std::string fcosmics_name_two = dcosm_name + "cosmics_" + horn_two + "__numu2018.root";
    TFile fcosmics_two(fcosmics_name_two.c_str());
    hCosmics.push_back((TH1*)fcosmics_two.Get("cosmics_q0")); // remember, quant0 = all                                                                                     
    hCosmics.push_back((TH1*)fcosmics_two.Get("cosmics_q1"));
    hCosmics.push_back((TH1*)fcosmics_two.Get("cosmics_q2"));
    hCosmics.push_back((TH1*)fcosmics_two.Get("cosmics_q3"));
    hCosmics.push_back((TH1*)fcosmics_two.Get("cosmics_q4"));
    auto *cosmics_quant1_two = (TH1*)fcosmics_two.Get("cosmics_q1");
    auto *cosmics_quant2_two = (TH1*)fcosmics_two.Get("cosmics_q2");
    auto *cosmics_quant3_two = (TH1*)fcosmics_two.Get("cosmics_q3");
    auto *cosmics_quant4_two = (TH1*)fcosmics_two.Get("cosmics_q4");
    s_cosmics.push_back(Spectrum (cosmics_quant1_two, pot_two, livetime_two));
    s_cosmics.push_back(Spectrum (cosmics_quant2_two, pot_two, livetime_two));
    s_cosmics.push_back(Spectrum (cosmics_quant3_two, pot_two, livetime_two));
    s_cosmics.push_back(Spectrum (cosmics_quant4_two, pot_two, livetime_two));
  }
  if(twobeams){
    for(int quantId = 5; quantId < totquant+1; quantId++){
      hCosmics[0]->Add(hCosmics[quantId]);
    }
  }
  



  std::cout << "\nloading predictions" << std::endl;
  const std::string dpred_name_fhc= dpred_name+"pred_numuconcat_fhc__numu2018.root";
  const std::string dpred_name_rhc= dpred_name+"pred_numuconcat_rhc__numu2018.root";
  std::string dpred_name_one = dpred_name_rhc;
  std::string dpred_name_two = dpred_name_fhc;
  if(anaName=="fhc") dpred_name_one = dpred_name_fhc;

  std::vector<PredictionInterp*> predictions;
  auto systs = getAllAna2018Systs(kNumuAna2018);
  int systSize = systs.size();

  //predictions here
  std::cout << "\nloading predictions" << std::endl;
  std::vector<PredictionSystJoint2018*> predictions;
  ENu2018ExtrapType extrap = kNuMu;
  for(int quant = 1; quant < totquant+1; quant++){
    std::cout << "quantile " << quant << std::endl;
    if(quant<5) predictions.push_back(new PredictionSystJoint2018(extrap, calc, horn_one, quant, systs, dpred_name_one));
    else predictions.push_back(new PredictionSystJoint2018(extrap, calc, horn_two, quant-4, systs, dpred_name_two));
  }


  /*
  //predictions old style                                                                                                                                             
  std::cout << "\nloading predictions" << std::endl;
  for(int quant = 1; quant < totquant+1; quant++){
    std::cout << "quantile " << quant << std::endl;
    if(quant<5){
      fpred_suffix = extrapName + "_stats__" + horn_one + "_" + calcName + "_" + hierarchy + "__" + period_one + "__numu2018.root";
      fpred_name = horn_one + "/predictions/" + extrapName + "/" + Form("prediction_quant%d__", quant);
      std::string input = dpred_name + fpred_name + fpred_suffix;
      TFile* file       = new TFile(input.c_str());
      predictions.push_back(LoadFrom<PredictionInterp>(file, "prediction")).release();
      file->Close();
    }
    else{
      fpred_suffix = extrapName + "_stats__" + horn_two + "_" + calcName + "_" + hierarchy + "__" + period_two + "__numu2018.root";
      fpred_name = horn_two + "/predictions/" + extrapName + "/" + Form("prediction_quant%d__", quant-4);
      std::string input = dpred_name + fpred_name + fpred_suffix;
      TFile* file       = new TFile(input.c_str());
      predictions.push_back(LoadFrom<PredictionInterp>(file, "prediction")).release();
      file->Close();
    }
  }//loop quant
  */


  // real or fake fake data                                                                                                                                           
  std::cout << "\nData" << std::endl;
  std::vector<Spectrum> s_predictions;
  std::vector<Spectrum> s_fakedata;

  std::cout << "\nloading fake data" << std::endl;
  for(int quant = 0; quant < totquant; quant++){
    std::cout << "quantile " << quant+1 << std::endl;
    HornSettings(quant, pot, livetime, wrongs, maxy);
    s_predictions.push_back(predictions[quant]->Predict(calcfake));
    s_fakedata.push_back(s_predictions[quant].FakeData(pot));
  }
  for(int quant = 0; quant < totquant; quant++) s_fakedata[quant] += s_cosmics[quant];
  
  
  std::vector<TH1*> hDataUnscaled;
  std::vector<TH1*> hData;
  std::vector<TH1*> hNC;
  std::vector<TH1*> hCC;
  std::vector<TH1*> hNuMu;
  std::vector<TH1*> hNoOsc;
  std::vector<TH1*> hPred;
  std::vector<TH1*> hBkg;
  std::vector<TH1*> hAllBkg;
  std::vector<TH1*> hWS;  


  HornSettings(0, pot, livetime, wrongs, maxy);
  hData.push_back(s_fakedata[0].ToTH1(pot));
  hDataUnscaled.push_back(s_fakedata[0].ToTH1(pot));  
  hNC.push_back(predictions[0]->PredictComponent(calc, Flavors::kAll, Current::kNC, Sign::kBoth).ToTH1(pot));
  hCC.push_back(predictions[0]->PredictComponent(calc, Flavors::kAll, Current::kCC, Sign::kBoth).ToTH1(pot));
  hNuMu.push_back(predictions[0]->PredictComponent(calc, Flavors::kAllNuMu, Current::kCC, Sign::kBoth).ToTH1(pot));
  hBkg.push_back(predictions[0]->Predict(calc).ToTH1(pot));
  hAllBkg.push_back(predictions[0]->Predict(calc).ToTH1(pot));
  hPred.push_back(predictions[0]->Predict(calc).ToTH1(pot));
  hNoOsc.push_back(predictions[0]->Predict(calcNoOsc).ToTH1(pot));
  hWS.push_back(predictions[0]->PredictComponent(calc, Flavors::kAllNuMu, Current::kCC, wrongs).ToTH1(pot));

  for(int quantId = 0; quantId < totquant; quantId++){
    HornSettings(quantId, pot, livetime, wrongs, maxy);
    hData.push_back(s_fakedata[quantId].ToTH1(pot));
    hDataUnscaled.push_back(s_fakedata[quantId].ToTH1(pot));
    hNC.push_back(predictions[quantId]->PredictComponent(calc, Flavors::kAll, Current::kNC, Sign::kBoth).ToTH1(pot));
    hCC.push_back(predictions[quantId]->PredictComponent(calc, Flavors::kAll, Current::kCC, Sign::kBoth).ToTH1(pot));
    hNuMu.push_back(predictions[quantId]->PredictComponent(calc, Flavors::kAllNuMu, Current::kCC, Sign::kBoth).ToTH1(pot));
    hBkg.push_back(predictions[quantId]->Predict(calc).ToTH1(pot));
    hAllBkg.push_back(predictions[quantId]->Predict(calc).ToTH1(pot));
    hPred.push_back(predictions[quantId]->Predict(calc).ToTH1(pot));
    hNoOsc.push_back(predictions[quantId]->Predict(calcNoOsc).ToTH1(pot));
    hWS.push_back(predictions[quantId]->PredictComponent(calc, Flavors::kAllNuMu, Current::kCC, wrongs).ToTH1(pot));
    if(quantId > 0){ // already added the first quantile
      hDataUnscaled[0] ->Add(s_fakedata[quantId].ToTH1(pot));
      hData[0]  ->Add(s_fakedata[quantId].ToTH1(pot));
      hNC[0]    ->Add(predictions[quantId]->PredictComponent(calc, Flavors::kAll, Current::kNC, Sign::kBoth).ToTH1(pot));
      hCC[0]    ->Add(predictions[quantId]->PredictComponent(calc, Flavors::kAll, Current::kCC, Sign::kBoth).ToTH1(pot));
      hNuMu[0]  ->Add(predictions[quantId]->PredictComponent(calc, Flavors::kAllNuMu, Current::kCC, Sign::kBoth).ToTH1(pot));
      hBkg[0]   ->Add(predictions[quantId]->Predict(calc).ToTH1(pot));
      hAllBkg[0]->Add(predictions[quantId]->Predict(calc).ToTH1(pot));
      hPred[0]  ->Add(predictions[quantId]->Predict(calc).ToTH1(pot));
      hNoOsc[0] ->Add(predictions[quantId]->Predict(calcNoOsc).ToTH1(pot));
      hWS[0]    ->Add(predictions[quantId]->PredictComponent(calc, Flavors::kAllNuMu, Current::kCC, wrongs).ToTH1(pot));    
    }
  }
  // fake stack: draw bkg, cosmics, ws in this order
  // add bkg and cosmics to ws
  // add cosmics to bkg
  for(int quantId = 0; quantId < totquant+1; quantId++){
    hBkg[quantId]->Add(hNuMu[quantId], -1);
    hAllBkg[quantId]->Add(hNuMu[quantId], -1);
    hAllBkg[quantId]->Add(hCosmics[quantId]);
    hWS[quantId]->Add(hAllBkg[quantId]);
    hPred[quantId]->Add(hCosmics[quantId]);
  }



  // scale everything, everything!                                                                                                                                                  
  std::cout << "\nScaling histograms to 0.1 GeV" << std::endl;
  for(int quantId = 0; quantId < totquant+1; quantId++){
    hData[quantId]   ->Scale(0.1, "width");
    hCosmics[quantId]->Scale(0.1, "width");
    hNC[quantId]     ->Scale(0.1, "width");
    hCC[quantId]     ->Scale(0.1, "width");
    hBkg[quantId]    ->Scale(0.1, "width");
    hAllBkg[quantId] ->Scale(0.1, "width");
    hNuMu[quantId]   ->Scale(0.1, "width");
    hPred[quantId]   ->Scale(0.1, "width");
    hNoOsc[quantId]  ->Scale(0.1, "width");
    hWS[quantId]     ->Scale(0.1, "width");
  }



  std::vector<Spectrum> sPrediction;                                                                                                                                               
  for(int quantId = 0; quantId < totquant+1; quantId++){
    if(quantId==0){
      HornSettings(0, pot, livetime, wrongs, maxy);
      Spectrum prediction(hPred[quantId], pot, livetime);                                                                                                  
      sPrediction.push_back(prediction);                                                                                                                                           
    }                                                                                                                                                                              
    else{                                                                                                                                                                          
      HornSettings(quantId-1, pot, livetime, wrongs, maxy);
      Spectrum prediction(hPred[quantId], pot, livetime);                                                                                              
      sPrediction.push_back(prediction);                                                                                                                            
    }                                                                                                                                                                              
  }                   


  // save up and dn shifts for each quantile
  // looks messy but really makes it less difficult                                                                                                                                
  std::vector<Spectrum> ups[totquant+1], dns[totquant+1];
  std::cout << "\nFilling ups and dns\n" << std::endl;

  if(errorBand){
  for(int quantId = 0; quantId < totquant+1; quantId++){

    std::cout << "quantile " << quantId << std::endl;
    // spectrums for syst band best fit                                                                                                                                            
    for(const ISyst* syst: systs){
      SystShifts shifts;

      if(quantId==0){
        HornSettings(0, pot, livetime, wrongs, maxy);
        shifts.SetShift(syst, +1);
        TH1 * htempup = predictions[0]->PredictSyst(calc, shifts).ToTH1(pot);
        for(int i = 1; i < totquant; i++){htempup->Add(predictions[i]->PredictSyst(calc, shifts).ToTH1(pot));}
        htempup->Scale(0.1, "width");
        htempup->Add(hCosmics[0]);
        Spectrum specup(htempup, pot, livetime);
        ups[quantId].push_back(specup);

        shifts.SetShift(syst, -1);
        TH1 * htempdn = predictions[0]->PredictSyst(calc, shifts).ToTH1(pot);
        for(int i = 1; i < totquant; i++){htempdn->Add(predictions[i]->PredictSyst(calc, shifts).ToTH1(pot));}
        htempdn->Scale(0.1, "width");
        htempdn->Add(hCosmics[0]);
        Spectrum specdn(htempdn, pot, livetime);
        dns[quantId].push_back(specdn);
      }// all quantiles

      else{
        HornSettings(quantId-1, pot, livetime, wrongs, maxy);
        shifts.SetShift(syst, +1);
        TH1* htempup = predictions[quantId-1]->PredictSyst(calc, shifts).ToTH1(pot);
        htempup->Scale(0.1, "width");
        htempup->Add(hCosmics[quantId]);
        Spectrum specup(htempup, pot, livetime);
        ups[quantId].push_back(specup);
        shifts.SetShift(syst, -1);
        TH1 * htempdn = predictions[quantId-1]->PredictSyst(calc, shifts).ToTH1(pot);
        htempdn->Scale(0.1, "width");
        htempdn->Add(hCosmics[quantId]);
        Spectrum specdn(htempdn, pot, livetime);
        dns[quantId].push_back(specdn);
      }
    } // shifts
  }// for quantiles
  }

  

  // pimping  
  for(int quantId = 0; quantId< 1; quantId++){
    PimpHist(hData[quantId], line_data, color_data, 3, marker_data, color_data, 1);
    PimpHist(hWS[quantId], line_ws, color_ws, 1, marker_ws, color_ws, 1);
    PimpHist(hNC[quantId], line_nc, color_nc, 1, marker_nc, color_nc, 1);
    PimpHist(hPred[quantId], line_pred, color_pred, 3, marker_pred, color_pred, 1);
    PimpHist(hCosmics[quantId], line_cos, color_cos, 1, marker_cos, color_cos, 1);
    PimpHist(hAllBkg[quantId], line_bkg, color_bkg, 1, marker_bkg, color_bkg, 1);
    PimpHist(hBkg[quantId], line_bkg, color_bkg, 1, marker_bkg, color_bkg, 1);
    FillWithDimColor(hWS[quantId], 0);
    FillWithDimColor(hCosmics[quantId], 0);
    FillWithDimColor(hAllBkg[quantId], 0);
  }// for quantId  
  //root is not behaving for the quantiles
  //set different line width and marker size
  for(int quantId = 1; quantId< totquant+1; quantId++){
    PimpHist(hData[quantId], line_data, color_data, 2, marker_data, color_data, 2);
    PimpHist(hWS[quantId], line_ws, color_ws, 2, marker_ws, color_ws, 2);
    PimpHist(hNC[quantId], line_nc, color_nc, 2, marker_nc, color_nc, 2);
    PimpHist(hPred[quantId], line_pred, color_pred, 2, marker_pred, color_pred, 2);
    PimpHist(hCosmics[quantId], line_cos, color_cos, 2, marker_cos, color_cos, 2);
    PimpHist(hAllBkg[quantId], line_bkg, color_bkg, 2, marker_bkg, color_bkg, 2);
    PimpHist(hBkg[quantId], line_bkg, color_bkg, 2, marker_bkg, color_bkg, 2);
    FillWithDimColor(hWS[quantId], 0);
    FillWithDimColor(hCosmics[quantId], 0);
    FillWithDimColor(hAllBkg[quantId], 0);
  }// for quantId


  // for later use in the 4pad plots
  TLatex* xTitle = new TLatex(0.5, 0.03, "Reconstructed Energy (GeV)");
  xTitle->SetNDC();
  xTitle->SetTextAlign(22);
  TLatex* yTitle = new TLatex(0.02, 0.5, "Events / 0.1 GeV");
  yTitle->SetNDC();
  yTitle->SetTextAlign(22);
  yTitle->SetTextAngle(90.);


  const int cx = 600, cy = 500;
  // all quantiles together
  // equals quant 0
  if(DianaIsLazy){
    TString outName = dout_name + "prediction_spectra__" + anaName;
    std::string idName = anaName + "_" + extrapName;
    int quantId = 0;
    TCanvas *canvas = new TCanvas("canvas","canvas", cx, cy);
    canvas->cd();
    ana::CenterTitles(hPred[0]);
    hPred[0]->GetYaxis()->SetTitle("Events / 0.1 GeV");
    hPred[0]->GetYaxis()->SetRangeUser(0.0, 1.3 * hPred[0]->GetMaximum());
    HornSettings(0, pot, livetime, wrongs, maxy);
    if(errorBand) PlotWithSystErrorBand_Quant(quantId, sPrediction[quantId], ups[quantId], dns[quantId], pot, color_pred, color_band, hPred[quantId]->GetMaximum(), true);
    else hPred[quantId]->Draw("hist same");
    hWS[quantId]->Draw("hist same");
    hAllBkg[quantId]->Draw("hist same");
    hCosmics[quantId]->Draw("hist same");
    hPred[quantId]->Draw("hist same");
    //hAllBkg[quantId]->Draw("l p same");
    //hCosmics[quantId]->Draw("hist p same");
    TLegend *legend = new TLegend(0.58,0.60,0.78,0.85,"","NDC");
    legend->AddEntry(hPred[quantId], "Prediction","l");
    TLegendEntry *entry=legend->AddEntry("error","1-#sigma syst. range","bf");
    entry->SetFillStyle(1001);
    entry->SetFillColor(color_band);
    entry->SetLineColor(color_band);
    legend->AddEntry(hAllBkg[quantId], "Total bkgd.","bf");
    legend->AddEntry(hCosmics[quantId], "Cosmic bkgd.","bf");
    legend->AddEntry(hWS[quantId], "WS","bf");
    //legend->AddEntry(hBkg[quantId], "Beam bkg.","lp");
    //legend->AddEntry(hCosmics[quantId], "Cosmic bkg.","lp");
    //legend->AddEntry(hWS[quantId], "WS","l");
    legend->SetTextSize(0.05); //no border for legend                                                                                                                           
    legend->SetBorderSize(0); //no border for legend                                                                                                                             
    legend->SetFillStyle(0);  //fill colour is white                                                                                                                             
    legend->Draw();
    CornerLabel(idName);
    drawQuantLabel(0);
    gPad->RedrawAxis();
    gPad->Update();
    canvas->cd();
    if(errorBand){
      xTitle->Draw();
      yTitle->Draw();
    }
    
    canvas->SaveAs(outName + ".eps");
    canvas->SaveAs(outName + ".png");
    canvas->SaveAs(outName + ".pdf");

  } // plot all quantiles  


  if(DianaIsLazy){

    TCanvas *canvas;
    TPad *pad1, *pad2, *pad3, *pad4;
    SplitCanvasQuant(canvas, pad1, pad2, pad3, pad4);
    TString outName = dout_name + "prediction_spectra__" + anaName;
    std::string idName = anaName + "_" + extrapName;
    for(int quantId = 1; quantId < 5; quantId++){
      if(quantId==1) pad1->cd();
      if(quantId==2) pad2->cd();
      if(quantId==3) pad3->cd();
      if(quantId==4) pad4->cd();
      HornSettings(quantId-1, pot, livetime, wrongs, maxy);
      MakeHistCanvasReady_Quant(quantId, hPred[quantId], maxy);
      if(errorBand) PlotWithSystErrorBand_Quant(quantId, sPrediction[quantId], ups[quantId], dns[quantId], pot, color_pred, color_band, maxy, true);
      else hPred[quantId]->Draw("hist same");
      hWS[quantId]->Draw("hist same");
      hAllBkg[quantId]->Draw("hist same");
      hCosmics[quantId]->Draw("hist same");
      hPred[quantId]->Draw("hist same");
      //hAllBkg[quantId]->Draw("l p same");
      //hCosmics[quantId]->Draw("hist p same");
      //grData[quantId]  ->Draw("p0 same");
      drawQuantLabel(quantId);
      if(quantId==1){
        pad1->cd();
        TLegend *legend = new TLegend(0.09,0.70,0.27,0.88,"","NDC");
        legend->AddEntry(hPred[quantId], "Prediction","l");
        TLegendEntry *entry=legend->AddEntry("error","1-#sigma syst. range","bf");
        entry->SetFillStyle(1001);
        entry->SetFillColor(color_band);
        entry->SetLineColor(color_band);
        legend->AddEntry(hAllBkg[quantId], "Total bkgd.","bf");
        legend->AddEntry(hCosmics[quantId], "Cosmic bkgd.","bf");
        legend->AddEntry(hWS[quantId], "WS","bf");
        //legend->AddEntry(hBkg[quantId], "Beam bkg.","lp");
        //legend->AddEntry(hCosmics[quantId], "Cosmic bkg.","lp");
        //legend->AddEntry(hWS[quantId], "WS","l");
        legend->SetTextSize(0.025); //no border for legend                                                                                                                         
        legend->SetBorderSize(0); //no border for legend                                                                                                                           
        legend->SetFillStyle(0);  //fill colour is white                                                                                                                           
        legend->Draw();
        CornerLabel(idName);
        gPad->Update();
      }
    }// for quantiles    


    for(auto & p:{pad1,pad2,pad3,pad4}){
      p->RedrawAxis();
    }
    canvas->Update();
    canvas->cd();
    xTitle->Draw();
    yTitle->Draw();
    canvas->SaveAs(outName + "_quant14.png");
    canvas->SaveAs(outName + "_quant14.eps");
    canvas->SaveAs(outName + "_quant14.pdf");

  } // plot quantiles 1 to 4



  if(DianaIsLazy && twobeams){
    TCanvas *canvas;
    TPad *pad1, *pad2, *pad3, *pad4;
    SplitCanvasQuant(canvas, pad1, pad2, pad3, pad4);
    TString outName = dout_name + "datapred_bestfit__" + anaName;
    std::string idName = anaName + "_" + extrapName;
    for(int quantId = 5; quantId <totquant+1; quantId++){ // revers prevents pad overlaping    
      if(quantId==5) pad1->cd();
      if(quantId==6) pad2->cd();
      if(quantId==7) pad3->cd();
      if(quantId==8) pad4->cd();
      HornSettings(quantId-1, pot, livetime, wrongs, maxy);
      MakeHistCanvasReady_Quant(quantId, hPred[quantId], maxy);
      if(errorBand) PlotWithSystErrorBand_Quant(quantId-4, sPrediction[quantId], ups[quantId], dns[quantId], pot, color_pred, color_band, maxy, true);
      else hPred[quantId]->Draw("hist same");
      hWS[quantId]->Draw("hist same");
      hAllBkg[quantId]->Draw("hist same");
      hCosmics[quantId]->Draw("hist same");
      hPred[quantId]->Draw("hist same");
      drawQuantLabel(quantId-4);
      if(quantId==1){
        pad1->cd();
        TLegend *legend = new TLegend(0.09,0.70,0.27,0.88,"","NDC");
        legend->AddEntry(hPred[quantId], "Prediction","l");
        TLegendEntry *entry=legend->AddEntry("error","1-#sigma syst. range","bf");
        entry->SetFillStyle(1001);
        entry->SetFillColor(color_band);
        entry->SetLineColor(color_band);
        legend->AddEntry(hAllBkg[quantId], "Total bkgd.","bf");
        legend->AddEntry(hCosmics[quantId], "Cosmic bkgd.","bf");
        legend->AddEntry(hWS[quantId], "WS","bf");
        //legend->AddEntry(hNC[quantId], "NC bkg.","lp");
        //legend->AddEntry(hBkg[quantId], "Beam bkg.","lp");
        //legend->AddEntry(hCosmics[quantId], "Cosmic bkg.","lp");
        //legend->AddEntry(hWS[quantId], "WS","l");
        legend->SetTextSize(0.025); //no border for legend
        legend->SetBorderSize(0); //no border for legend                                                                                                                           
        legend->SetFillStyle(0);  //fill colour is white
        legend->Draw();
        CornerLabel(idName);
        gPad->Update();
      }
    }// for quantiles
  
    for(auto & p:{pad1,pad2,pad3,pad4}){
      p->RedrawAxis();
    }
    canvas->Update();
    canvas->cd();
    xTitle->Draw();
    yTitle->Draw();    
    canvas->SaveAs(outName + "_quant58.png");
    canvas->SaveAs(outName + "_quant58.eps");
    canvas->SaveAs(outName + "_quant58.pdf");
    
  } // plot quantiles 1 to 4


  
} // End of function

Variable Documentation

Color_t color_band = kViolet-9

Definition at line 46 of file plot_predictions.C.

Referenced by plot_predictions().

Color_t color_bkg = kGray+1

Definition at line 52 of file plot_predictions.C.

Referenced by plot_predictions().

Color_t color_cos = kGreen+1

Definition at line 51 of file plot_predictions.C.

Referenced by plot_predictions().

Color_t color_data = kBlack

Definition at line 50 of file plot_predictions.C.

Referenced by plot_predictions().

Color_t color_nc = kBlue-7

Definition at line 49 of file plot_predictions.C.

Referenced by plot_predictions().

Color_t color_pred = kViolet+1

Definition at line 47 of file plot_predictions.C.

Referenced by plot_predictions().

Color_t color_ws = kBlue-7

Definition at line 48 of file plot_predictions.C.

Referenced by plot_predictions().

bool DianaIsLazy = true

Definition at line 41 of file plot_predictions.C.

bool errorBand = true

Definition at line 42 of file plot_predictions.C.

std::string extrapName = "extrap"

Definition at line 32 of file plot_predictions.C.

std::string horn_one = "rhc"

Definition at line 33 of file plot_predictions.C.

Referenced by plot_predictions().

std::string horn_two = "fhc"

Definition at line 34 of file plot_predictions.C.

Referenced by plot_predictions().

Style_t line_bkg = kSolid

Definition at line 58 of file plot_predictions.C.

Referenced by plot_predictions().

Style_t line_cos = kSolid

Definition at line 57 of file plot_predictions.C.

Referenced by plot_predictions().

Style_t line_data = kSolid

Definition at line 55 of file plot_predictions.C.

Referenced by plot_predictions().

Style_t line_nc = kSolid

Definition at line 54 of file plot_predictions.C.

Referenced by plot_predictions().

Style_t line_pred = kSolid

Definition at line 56 of file plot_predictions.C.

Referenced by plot_predictions().

Style_t line_ws = kSolid

Definition at line 53 of file plot_predictions.C.

Referenced by plot_predictions().

double livetime = 0

Definition at line 21 of file plot_predictions.C.

Referenced by plot_predictions().

const double livetime_fhc = kAna2018FHCLivetime

Definition at line 12 of file plot_predictions.C.

Referenced by plot_predictions().

double livetime_one = livetime_rhc

Definition at line 22 of file plot_predictions.C.

Referenced by HornSettings(), and plot_predictions().

const double livetime_rhc = kAna2018RHCLivetime

Definition at line 13 of file plot_predictions.C.

double livetime_two = livetime_fhc

Definition at line 23 of file plot_predictions.C.

Referenced by HornSettings(), and plot_predictions().

Style_t marker_bkg = kFullCircle

Definition at line 64 of file plot_predictions.C.

Referenced by plot_predictions().

Style_t marker_cos = kFullTriangleUp

Definition at line 63 of file plot_predictions.C.

Referenced by plot_predictions().

Style_t marker_data = kFullCircle

Definition at line 61 of file plot_predictions.C.

Referenced by plot_predictions().

Style_t marker_nc = kFullSquare

Definition at line 60 of file plot_predictions.C.

Referenced by plot_predictions().

Style_t marker_pred = kFullCircle

Definition at line 62 of file plot_predictions.C.

Referenced by plot_predictions().

Style_t marker_ws = kFullCircle

Definition at line 59 of file plot_predictions.C.

Referenced by plot_predictions().

double maxy = 0

Definition at line 24 of file plot_predictions.C.

double maxy_one = 1.999

Definition at line 25 of file plot_predictions.C.

double maxy_two = 3.999

Definition at line 26 of file plot_predictions.C.

std::string period_one = "full"

Definition at line 35 of file plot_predictions.C.

Referenced by plot_predictions().

std::string period_two = "full"

Definition at line 36 of file plot_predictions.C.

double pot = 0

Definition at line 18 of file plot_predictions.C.

Referenced by plot_datamcpred(), and plot_predictions().

const double pot_fhc = kAna2018FHCPOT

Definition at line 10 of file plot_predictions.C.

Referenced by plot_predictions().

double pot_one = pot_rhc

Definition at line 19 of file plot_predictions.C.

Referenced by HornSettings(), and plot_predictions().

const double pot_rhc = kAna2018RHCPOT

Definition at line 11 of file plot_predictions.C.

double pot_two = pot_fhc

Definition at line 20 of file plot_predictions.C.

Referenced by HornSettings(), and plot_predictions().

int totquant = 8

Definition at line 40 of file plot_predictions.C.

Referenced by plot_predictions().

bool twobeams = true

Definition at line 39 of file plot_predictions.C.

Referenced by plot_predictions().

Sign::Sign_t wrongs = Sign::kNu

Definition at line 27 of file plot_predictions.C.

Referenced by plot_predictions().

Sign::Sign_t ws_one = Sign::kNu

Definition at line 28 of file plot_predictions.C.

Referenced by HornSettings(), and plot_predictions().

Sign::Sign_t ws_two = Sign::kAntiNu

Definition at line 29 of file plot_predictions.C.

Referenced by HornSettings().