plot_datapredictions.C File Reference

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

Go to the source code of this file.

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

Variables
std::string	extrapName = "extrap"
bool	DianaIsLazy = true
bool	errorBand = true
double	maxy = 0
double	maxy_one = 4.999
double	maxy_two = 5.999
double	paperscale = 0.8/0.5
double	maxy_pad0 = 0
double	maxy_rhc = 9
double	maxy_fhc = 13

Function Documentation

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

Definition at line 57 of file plot_datapredictions.C.

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

Referenced by plot_datapredictions().

                                                                                                {
  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_datapredictions	(	std::string	anaName = "fhc",
		std::string	calculator = "2018calc",
		bool	nh = true
	)

Definition at line 75 of file plot_datapredictions.C.

References calc, make_mec_shifts_plots::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, file, ana::FillWithDimColor(), ana::getAllAna2018Systs(), GetXaxis(), GetYaxis(), ana::graphAsymmErrorScaled(), 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_fhc, maxy_one, maxy_pad0, maxy_rhc, runNovaSAM::outName, pad1, pad2, pad3, paperscale, period_one, ana::PimpHist(), ana::PlotWithSystErrorBand_Quant(), POT, pot, pot_fhc, pot_one, pot_two, preliminary, PreliminaryBoxOpening(), runNovaSAM::release, RestartCalculator(), ana::SystShifts::SetShift(), ana::SplitCanvasQuant(), systs, totquant, twobeams, wrongs, and ws_one.

{
  
  
  bool preliminary = false;
  std::string anaName2 ="";
  if(anaName=="rhc"){
    anaName2="Antineutrino beam";
    totquant = 4;
    twobeams = false;
    maxy_pad0 = maxy_rhc;
  }
  // 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"){
    anaName2="Neutrino beam";
    horn_one = "fhc";
    period_one = "full";
    pot_one = pot_fhc;
    livetime_one = livetime_fhc;
    ws_one = Sign::kAntiNu;
    maxy_one = 5.999;
    twobeams = false;
    totquant = 4;
    maxy_pad0 = maxy_fhc;
  }

  // calculator
  std::string calcName = "2018calc"; 
  std::string hierarchy = "";
  osc::IOscCalcAdjustable* calc = DefaultOscCalc();
  RestartCalculator(calc, anaName, nh);
  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 = "";

  TFile* fdata_one = new TFile((ddata_name + "fd_dataspectrum_" + horn_one + "_full__numu2018.root").c_str());
  TFile* fdata_two = new TFile((ddata_name + "fd_dataspectrum_" + horn_two + "_full__numu2018.root").c_str());
  std::vector<Spectrum*> s_realdata;
  std::cout << "\nloading data files and spectra" << std::endl;
  for(int quant=1; quant<totquant+1; quant++){
    if(quant<5) s_realdata.push_back(LoadFrom<Spectrum>(fdata_one, Form("fd_data_q%d", quant) ).release());
    else s_realdata.push_back(LoadFrom<Spectrum>(fdata_two, Form("fd_data_q%d", quant-4) ).release());
  }


  double pot1here = s_realdata[0]->POT();
  std::cout << "\n\n pot for this sample = " << pot1here << "\n\n" << std::endl;

  // 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_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< const ISyst* > systs;
  systs = getAllAna2018Systs(kNumuAna2018,true);
  int systSize = systs.size();

  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));
  }


  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_realdata[0]->ToTH1(s_realdata[0]->POT()));
  hDataUnscaled.push_back(s_realdata[0]->ToTH1(s_realdata[0]->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_realdata[quantId]->ToTH1(s_realdata[0]->POT()));
    hDataUnscaled.push_back(s_realdata[quantId]->ToTH1(s_realdata[0]->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_realdata[quantId]->ToTH1(s_realdata[0]->POT()));
      hData[0]  ->Add(s_realdata[quantId]->ToTH1(s_realdata[0]->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]);
  }


  std::cout << "\n\nIntegrated events = " << hData[0]->Integral() << "\n\n" << std::endl;

  // 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");
  }

  // data graphs
  std::cout << "\nDefine data graphs\n" << std::endl;
  std::vector<TGraph*> grData;
  for(int quantId = 0; quantId < totquant+1; quantId++){
    grData.push_back(graphAsymmErrorScaled(hData[quantId], hDataUnscaled[quantId]));
  }



  // predictions again?
  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);
     grData[quantId]->SetMarkerSize(1.);
   }// 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, 1.);
     PimpHist(hWS[quantId], line_ws, color_ws, 2, marker_ws, color_ws, 3);
     PimpHist(hNC[quantId], line_nc, color_nc, 2, marker_nc, color_nc, 3);
     PimpHist(hPred[quantId], line_pred, color_pred, 2, marker_pred, color_pred, 3);
     PimpHist(hCosmics[quantId], line_cos, color_cos, 2, marker_cos, color_cos, 3);
     PimpHist(hAllBkg[quantId], line_bkg, color_bkg, 2, marker_bkg, color_bkg, 3);
     PimpHist(hBkg[quantId], line_bkg, color_bkg, 2, marker_bkg, color_bkg, 3);
     FillWithDimColor(hWS[quantId], 0);
     FillWithDimColor(hCosmics[quantId], 0);
     FillWithDimColor(hAllBkg[quantId], 0);
     grData[quantId]->SetMarkerSize(1.);
   }// 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);
   xTitle->SetTextSize(0.055 /paperscale);
   TLatex* yTitle = new TLatex(0.02, 0.5, "Events / 0.1 GeV");
   yTitle->SetNDC();
   yTitle->SetTextAlign(22);
   yTitle->SetTextAngle(90.);
   yTitle->SetTextSize(0.055 / paperscale);
   
   const int cx = 960, cy = 800;
   // all quantiles together
   // equals quant 0
   if(DianaIsLazy){
     std::ofstream file;
     TString outName = dout_name + "fd_dataprediction_spectra__" + anaName;
     std::string idName = anaName2;
     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, maxy_pad0);
     hPred[0]->GetXaxis()->SetLabelSize(hPred[0]->GetXaxis()->GetLabelSize());
     hPred[0]->GetYaxis()->SetLabelSize(hPred[0]->GetYaxis()->GetLabelSize());
     hPred[0]->Draw("hist");
     HornSettings(0, pot, livetime, wrongs, maxy);
     if(errorBand) PlotWithSystErrorBand_Quant(quantId, sPrediction[quantId], ups[quantId], dns[quantId], pot, color_pred, color_band, maxy_pad0, false);
     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");
     grData[quantId]->Draw("p0 same");
     grData[quantId]->Draw("p0 same");
     TLegend *legend = new TLegend(0.58,0.54,0.89,0.89);
     legend->AddEntry(hData[quantId], "FD Data","lep");
     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);
     if(anaName=="fhc")legend->AddEntry(hWS[quantId], "Wrong Sign:#bar{#nu}_{#mu}CC","bf");
     else legend->AddEntry(hWS[quantId], "Wrong Sign:#nu_{#mu}CC","bf");
     legend->AddEntry(hAllBkg[quantId], "Total bkg.","bf");
     legend->AddEntry(hCosmics[quantId], "Cosmic bkg.","bf");
     legend->SetTextSize(0.04);
     legend->SetBorderSize(0);
     legend->SetFillStyle(0);
     legend->Draw();
     CornerLabel(idName);
     PreliminaryBoxOpening();
     drawQuantLabel(0);
     gPad->RedrawAxis();
     gPad->Update();
     canvas->cd();
     
     canvas->SaveAs(outName + "__"+calcName+".png");
     canvas->SaveAs(outName + "__"+calcName+".pdf");
     canvas->SaveAs(outName + "__"+calcName+".root");
     canvas->SaveAs(outName + "__"+calcName+".C");
     file.open (outName + "__"+calcName+".txt");
     file<< "FD spectra showing data, prediction at best fit with systematics band about the best fit with systematics, and beam and cosmic backgrounds for the NuMI beam in " << anaName << " mode.";
     file.close();

   } // plot all quantiles  
   
   
   if(DianaIsLazy){
     
     std::ofstream file;
     TCanvas *canvas;
     TPad *pad1, *pad2, *pad3, *pad4;
     SplitCanvasQuant(canvas, pad1, pad2, pad3, pad4);
     gStyle->SetLineScalePS(3./paperscale);
     canvas->SetCanvasSize(paperscale * canvas->GetWw(), paperscale * canvas->GetWh());
     TString outName = dout_name + "fd_dataprediction_spectra__" + anaName;
     std::string idName = anaName2;
     for(int quantId:{3,4,1,2}){
       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);
       // hack to avoid 0-5 overlapping in the x-axis
       hPred[quantId]->GetXaxis()->SetLabelSize(hPred[quantId]->GetXaxis()->GetLabelSize()/paperscale);
       hPred[quantId]->GetYaxis()->SetLabelSize(hPred[quantId]->GetYaxis()->GetLabelSize()/paperscale);
       hPred[quantId]->Draw("hist");
       if(quantId==3){
         TGaxis* gax1 = new TGaxis(0., 0., 4., 0., 0.,4.,4, "");
         gax1->SetLabelSize(hPred[quantId]->GetXaxis()->GetLabelSize());
         gax1->SetLabelOffset(hPred[quantId]->GetXaxis()->GetLabelOffset());
         gax1->SetLabelFont(hPred[quantId]->GetXaxis()->GetLabelFont());
         hPred[quantId]->GetXaxis()->SetLabelSize(0.);
         hPred[quantId]->Draw("hist");
         gax1->Draw("same");
       }

       if(errorBand) PlotWithSystErrorBand_Quant(quantId, sPrediction[quantId], ups[quantId], dns[quantId], pot, color_pred, color_band, maxy, false);
       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");
       grData[quantId]->Draw("p0 same");
       drawQuantLabel(quantId);
       if(quantId==1){
         pad1->cd();
         TLegend *legend = new TLegend(0.09,0.73,0.27,0.89,"","NDC");
         legend->AddEntry(hData[quantId], "FD Data","lep");
         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);
         if(anaName=="fhc")legend->AddEntry(hWS[quantId], "Wrong Sign:#bar{#nu}_{#mu}CC","bf");
         else legend->AddEntry(hWS[quantId], "Wrong Sign:#nu_{#mu}CC","bf");
         legend->AddEntry(hAllBkg[quantId], "Total bkg.","bf");
         legend->AddEntry(hCosmics[quantId], "Cosmic bkg.","bf");
         legend->SetTextSize(0.025);
         legend->SetBorderSize(0);
         legend->SetFillStyle(0);
         legend->Draw();
         CornerLabel(idName);
         PreliminaryBoxOpening();
         gPad->Update();
       }
     }// for quantiles    
     
     
     for(auto & p:{pad3,pad4,pad1,pad2}){
       p->RedrawAxis();
     }

     canvas->Update();
     canvas->cd();
     xTitle->Draw();
     yTitle->Draw();
     canvas->SaveAs(outName + "_quant14__"+calcName+".png");
     canvas->SaveAs(outName + "_quant14__"+calcName+".pdf");
     canvas->SaveAs(outName + "_quant14__"+calcName+".root");
     canvas->SaveAs(outName + "_quant14__"+calcName+".C");
     file.open (outName + "_quant14__"+calcName+".txt");
     file<< "FD spectra showing data, prediction at best fit with systematics band about the best fit with systematics, and beam and cosmic backgrounds for the NuMI beam in " << anaName << " mode, separated per quartiles.";
     file.close();

   } // plot quantiles 1 to 4
   
      
} // End of function

void RestartCalculator	(	osc::IOscCalcAdjustable *	calc,
		std::string	anaName,
		bool	nh
	)

Definition at line 24 of file plot_datapredictions.C.

References std::asin(), delta_rhcfhc_joint, ana::ResetOscCalcToDefault(), osc::_IOscCalcAdjustable< T >::SetDmsq32(), osc::_IOscCalcAdjustable< T >::SetTh23(), sin_rhcfhc_joint, and std::sqrt().

Referenced by plot_datapredictions().

{

  // take different values from ../settings.h
  // either numu only or joint BF
  // we use the joint BF for neutrino
  double sinHere = sin_rhcfhc_joint;
  double deltaHere = delta_rhcfhc_joint;
  // naively changing hierarchy
  // this isnt totally right
  /*
  double deltaSign = 1.;
  if(nh==false) deltaSign = -1.;

  if(anaName=="rhc"){
    sinHere = sin_rhcfhc;
    deltaHere = delta_rhcfhc*deltaSign;
  }
  if(anaName=="fhc"){
    sinHere = sin_rhcfhc;
    deltaHere = delta_rhcfhc*deltaSign;
  }
  */
  ana::ResetOscCalcToDefault(calc);
  calc->SetTh23(asin(sqrt(sinHere)));
  calc->SetDmsq32(deltaHere);

}

Variable Documentation

bool DianaIsLazy = true

Definition at line 12 of file plot_datapredictions.C.

Referenced by plot_datapredictions().

bool errorBand = true

Definition at line 13 of file plot_datapredictions.C.

Referenced by plot_datapredictions().

std::string extrapName = "extrap"

Definition at line 10 of file plot_datapredictions.C.

Referenced by make_prediction_extrap(), make_prediction_noextrap(), make_surfprof(), make_surfprof_sensitivity(), and plot_contprof().

double maxy = 0

Definition at line 15 of file plot_datapredictions.C.

Referenced by Arrows(), Beam(), BiCountCanvas(), BiProbCanvas(), Compare(), containment_optimization_plots(), ana::DataMCAreaNormalizedRatio(), ana::DataMCRatio(), demoFitSlices(), demoStarPlots(), draw_percentiles(), DrawContourCanvas(), fiducial_optimization_plots(), fiducial_optimization_plots_2d(), comi::NumiFiltering::filter(), hv::HoughVertexAlg::FindVertex(), joint_fit_2017_slices(), joint_fit_2018_slices(), joint_fit_2019_slices(), make_dst_cosrejbdttrain(), make_starplots(), muonid_optimization_plots(), ana::Nue2018ThreeBinAxis(), ana::CutOptimization::OptimizedSigmaOverSigma(), Plot1D(), plot_datapredictions(), plot_joint_fit_2020_slices(), ana::PlotAll(), PlotEffs(), PlotMultiverse(), PlotSpectra(), print_tables(), calib::BetheBlochFit::ProfileMPV(), and Tutorial2019FitSlices().

double maxy_fhc = 13

Definition at line 22 of file plot_datapredictions.C.

Referenced by plot_datapredictions().

double maxy_one = 4.999

Definition at line 16 of file plot_datapredictions.C.

Referenced by HornSettings(), and plot_datapredictions().

double maxy_pad0 = 0

Definition at line 20 of file plot_datapredictions.C.

Referenced by plot_datapredictions().

double maxy_rhc = 9

Definition at line 21 of file plot_datapredictions.C.

Referenced by plot_datapredictions().

double maxy_two = 5.999

Definition at line 17 of file plot_datapredictions.C.

Referenced by HornSettings().

double paperscale = 0.8/0.5

Definition at line 18 of file plot_datapredictions.C.

Referenced by plot_datapredictions().