make_starplots.C File Reference

#include "3FlavorAna/Ana2018/numu/includes.h"
#include "../numu_tools.h"
#include "../settings.h"
#include "includes_starplots.h"
#include "CAFAna/Core/Progress.h"
#include "TCanvas.h"
#include "TFile.h"
#include "TGraph.h"
#include "TGraphAsymmErrors.h"
#include "TH1.h"
#include "TH2.h"
#include "TLatex.h"
#include "TLegend.h"
#include "TStyle.h"

Go to the source code of this file.

Functions
void	RestartCalculator (osc::IOscCalcAdjustable *calc, std::string anaName, bool nh)
void	make_starplots (const std::string anaName="rhcfhc", const int quantN=0, bool usecosmics=true, bool nh=true)

Variables
std::string	dataName = "fakedata"
std::vector< Style_t >	styleline = {kSolid, kSolid, kSolid, kSolid, kSolid, kSolid, kSolid, kSolid, kSolid}
std::vector< Color_t >	color = {kMagenta, kGreen+1, kOrange+1, kBlue-7, kCyan+1, kRed-4, kBlue+1, kGray+1, kOrange+3}

Function Documentation

void make_starplots	(	const std::string	anaName = "rhcfhc",
		const int	quantN = 0,
		bool	usecosmics = true,
		bool	nh = true
	)

Definition at line 53 of file make_starplots.C.

References AddEntry(), file_size_ana::axes, calc, canvas(), cd(), Close(), color, make_associated_cosmic_defs::cosmics, om::cout, ana::DefaultOscCalc(), DrawLegendBFNull(), allTimeWatchdog::endl, HTMLTools::entry(), FakeData(), file, ana::IFitter::Fit(), MECModelEnuComparisons::g, ana::Gaussian90Percent2D(), ana::getAllAna2018Systs(), FitSinSqTheta23UO::GetValue(), ana::ConstrainedFitVarWithPrior::GetValue(), FitSinSqTheta23LO::GetValue(), horn_one, horn_two, inFile, ana::kFitDmSq32Scaled, ana::kFitSinSqTheta23, kFitSinSqTheta23LO, kFitSinSqTheta23UO, ana::kNuMu, ana::kNumuAna2018, ana::IFitter::kQuiet, ana::kSolarConstraintsPDG2017, MECModelEnuComparisons::leg, livetime, livetime_fhc, livetime_one, livetime_two, maxy, MyCalibSysts, MyFluxSysts, MyGeniePCASysts, MyGenieRwgtFrSysts, MyGenieRwgtSysts, MyMECSysts, MyOtherSysts, MyXSecSysts, confusionMatrixTree::out, period_one, pot, pot_fhc, pot_one, pot_two, cacheDefinitionData::prog, runNovaSAM::release, RestartCalculator(), SetFillColor(), SetFillStyle(), ana::SystShifts::SetShift(), sigma(), string, styleline, systs, totquant, twobeams, and ana::WorldReactorConstraint2017().

{

  float minx      = 0.35, maxx      = 0.65;
  float minx_zoom = 0.58, maxx_zoom = 0.60;
  float miny      = 2.30, maxy      = 2.70;
  float miny_zoom = 2.43, maxy_zoom = 2.50;

  if(anaName=="rhc"){
    twobeams = false;
    totquant = 4;
    minx      = 0.55; maxx      = 0.70;
    minx_zoom = 0.62; maxx_zoom = 0.68;
    miny      = 2.20; maxy      = 3.00;
    miny_zoom = 2.555;maxy_zoom = 2.65;
  }
  if(anaName=="fhc"){
    horn_one = "fhc";
    period_one = "full";
    pot_one = pot_fhc;
    livetime_one = livetime_fhc;
    twobeams = false;
    totquant = 4;
    minx      = 0.35; maxx      = 0.65;
    minx_zoom = 0.48; maxx_zoom = 0.55;
    miny      = 2.20; maxy      = 2.70;
    miny_zoom = 2.40; maxy_zoom = 2.55;
  }


  // calculator
  osc::IOscCalcAdjustable* calc = DefaultOscCalc();
  RestartCalculator(calc, anaName, nh);

  // full systematics
  std::vector< const ISyst* > systs;
  systs = getAllAna2018Systs(kNumuAna2018,true);

  std::map< TString, std::vector<const ISyst*> > mapSystLists;
  mapSystLists["GenieRwgt"] = MyGenieRwgtSysts;
  mapSystLists["GenieRwgtFr"] = MyGenieRwgtFrSysts;
  mapSystLists["MEC"] = MyMECSysts;
  mapSystLists["XSec"] = MyXSecSysts;
  mapSystLists["GeniePCA"] = MyGeniePCASysts;
  mapSystLists["FluxSysts"] = MyFluxSysts ;
  mapSystLists["CalibSysts"] = MyCalibSysts;
  mapSystLists["OtherSysts"] = MyOtherSysts ;

  // get contour                       
  std::string inFile = "../results/surfprof_realdata_cosmics_"+anaName+"__extrap_stats_2018calc_nh__numu2018.root";
  TFile* fSurfSens = new TFile(inFile.c_str());
  FrequentistSurface* surfaceSens = LoadFrom< FrequentistSurface >(fSurfSens,"surface").release();
  fSurfSens -> Close();

  double bestFitDm = surfaceSens->GetBestFitY();
  double bestFitSsq = surfaceSens->GetBestFitX();
  std::cout<< "Best fit from result contour: ssqth23: " << bestFitSsq
           << ", dm^2_32: " << bestFitDm <<std::endl;
  delete fSurfSens;


  //inputs
  TString dout_name      = "../plots_blessing/";
  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 = "";


  //predictions systs framework
  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::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));
  }


  // cosmics
  std::string cosmics = "nocosmics";
  std::vector<Spectrum> s_cosmics;
  if(usecosmics){
    pot = pot_one; livetime = livetime_one;
    cosmics = "cosmics";
    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());
    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, livetime));
    s_cosmics.push_back(Spectrum (cosmics_quant2_one, pot, livetime));
    s_cosmics.push_back(Spectrum (cosmics_quant3_one, pot, livetime));
    s_cosmics.push_back(Spectrum (cosmics_quant4_one, pot, livetime));
    if(twobeams){
      pot = pot_two; livetime = livetime_two;
      std::cout << "" << horn_two << std::endl;
      std::string fcosmics_name_two = dcosm_name + "cosmics_" + horn_two + "__numu2018.root";
      TFile fcosmics_two(fcosmics_name_two.c_str());
      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, livetime));
      s_cosmics.push_back(Spectrum (cosmics_quant2_two, pot, livetime));
      s_cosmics.push_back(Spectrum (cosmics_quant3_two, pot, livetime));
      s_cosmics.push_back(Spectrum (cosmics_quant4_two, pot, livetime));
    }
  }



  // actually do something
  for(auto const& thisSystList:  mapSystLists){

    std::string NameZoom = "starplot_"+anaName+"_"+(std::string)thisSystList.first+"__zoom";
    TCanvas* canvas_zoom = new TCanvas("canvas_zoom", (NameZoom).c_str());
    TString out_zoom = dout_name+NameZoom.c_str();

    std::string Name = "starplot_"+anaName+"_"+(std::string)thisSystList.first;
    TCanvas* canvas = new TCanvas("canvas", (Name).c_str());
    TString out = dout_name+Name.c_str();

    TH2F *axes_zoom = new TH2F("setaxes_zoom","", 30, minx_zoom, maxx_zoom, 50, miny_zoom, maxy_zoom);
    axes_zoom->GetXaxis()->SetTitle("sin^{2}#theta_{23}");
    axes_zoom->GetYaxis()->SetTitle("#Deltam^{2}_{32} (10^{-3} eV^{2})");
    axes_zoom->GetXaxis()->SetTitleOffset(0.95);
    axes_zoom->GetYaxis()->SetTitleOffset(0.85);
    axes_zoom->GetXaxis()->SetTitleSize(0.055);
    axes_zoom->GetYaxis()->SetTitleSize(0.055);
    axes_zoom->GetXaxis()->SetLabelSize(0.04);
    axes_zoom->GetYaxis()->SetLabelSize(0.04);
    axes_zoom->GetXaxis()->CenterTitle();
    axes_zoom->GetYaxis()->CenterTitle();
    axes_zoom->SetTitle("");
    axes_zoom->GetYaxis()->SetDecimals();

    TH2F *axes = new TH2F("setaxes","" ,38 ,minx, maxx, 52, miny, maxy);
    axes->GetXaxis()->SetTitle("sin^{2}#theta_{23}");
    axes->GetYaxis()->SetTitle("#Deltam^{2}_{32} (10^{-3} eV^{2})");
    axes->GetXaxis()->SetTitleOffset(0.95);
    axes->GetYaxis()->SetTitleOffset(0.85);
    axes->GetXaxis()->SetTitleSize(0.055);
    axes->GetYaxis()->SetTitleSize(0.055);
    axes->GetXaxis()->SetLabelSize(0.04);
    axes->GetYaxis()->SetLabelSize(0.04);
    axes->GetXaxis()->CenterTitle();
    axes->GetYaxis()->CenterTitle();
    axes->SetTitle("");
    axes->GetYaxis()->SetDecimals();
    axes->Draw();
    axes->Draw();


    canvas_zoom->cd();
    gPad->SetFillStyle(0);
    gPad->SetTopMargin(0.08);
    gPad->SetRightMargin(0.08);
    gPad->SetBottomMargin(0.12);
    gPad->SetLeftMargin(0.12);
    gStyle->SetTitleOffset(0.95,"x");
    gStyle->SetTitleOffset(0.85,"y");
    axes_zoom->Draw();
    surfaceSens->DrawContour(Gaussian90Percent2D(*surfaceSens), kSolid, kBlack);
    surfaceSens->DrawBestFit(kBlack, kFullStar);
    /*
    TH2* h2d = (TH2*) surfaceSens -> ToTH2();
    for(int binx = 1; binx<h2d->GetNbinsX(); binx++)
      for(int biny = 1; biny<h2d->GetNbinsY(); biny++)
        h2d-> SetBinContent(binx, biny, 0);
    h2d-> GetYaxis()->SetRangeUser(miny_zoom,maxy_zoom);
    h2d-> GetXaxis()->SetRangeUser(minx_zoom,maxx_zoom);
    h2d-> GetYaxis()->CenterTitle();
    h2d-> GetXaxis()->CenterTitle();
    h2d-> Draw();
    */


    canvas->cd();
    gPad->SetFillStyle(0);
    gPad->SetTopMargin(0.08);
    gPad->SetRightMargin(0.08);
    gPad->SetBottomMargin(0.12);
    gPad->SetLeftMargin(0.12);
    gStyle->SetTitleOffset(0.95,"x");
    gStyle->SetTitleOffset(0.85,"y");
    axes->Draw();
    surfaceSens->DrawContour(Gaussian90Percent2D(*surfaceSens), kSolid, kBlack);                                                                                                  
    surfaceSens->DrawBestFit(kBlack, kFullStar);   
    DrawLegendBFNull(surfaceSens->GetBestFitX(), surfaceSens->GetBestFitY(), kBlack, kFullStar);

    TLegend *leg = new TLegend(0.40, 0.65, 0.70, 0.85);
    TLegendEntry *entry = leg->AddEntry("NULL","NOvA NH #nu_{#mu}+#bar{#nu}_{#mu} Systematic","h");
    /*
    entry = legend->AddEntry("zero","90% CL Sensitivity","L");
    entry->SetLineColor(kBlack);
    entry->SetLineStyle(kSolid);
    entry->SetLineWidth(3);
    */
    int graphId = 0;


    std::vector<const ISyst*> systList = thisSystList.second;

    for(const ISyst* s: systList){
      Progress prog("Making star plot for "+s->ShortName());
      TGraph* g = new TGraph;
      TGraph* gLO = new TGraph;
      for(double sigma = -1; sigma < +1.05; sigma += 0.2){

        SystShifts systListShift;
        systListShift.SetShift(s, sigma);
        RestartCalculator(calc, anaName, nh);

        // fake data
        std::cout << "\nfilling fake data" << std::endl;
        std::vector<Spectrum> s_fakedata;
        for(int quant = 0; quant < totquant; quant++){
          std::cout << "quantile " << quant+1 << std::endl;
          if(quant<4){ pot = pot_one; livetime = livetime_one;}
          else {pot = pot_two; livetime = livetime_two;}
          s_fakedata.push_back( predictions[quant]->PredictSyst(calc, systListShift).FakeData(pot) );
        }

        // single sample experiments
        std::vector <const IExperiment*> experiments;
        if(quantN==0){
          for(int quant = 0; quant < totquant; quant++){
            if(usecosmics){
              s_fakedata[quant] += s_cosmics[quant];
              experiments.push_back(new SingleSampleExperiment(predictions[quant], s_fakedata[quant], s_cosmics[quant], 0.1));
            }
            else experiments.push_back(new SingleSampleExperiment(predictions[quant], s_fakedata[quant]));
          }
        }
        else{
          if(usecosmics){
            s_fakedata[quantN-1] += s_cosmics[quantN-1];
            experiments.push_back(new SingleSampleExperiment(predictions[quantN-1], s_fakedata[quantN-1], s_cosmics[quantN-1], 0.1));
          }
          else experiments.push_back(new SingleSampleExperiment(predictions[quantN-1], s_fakedata[quantN-1]));
        }

        experiments.push_back(&kSolarConstraintsPDG2017);
        experiments.push_back(WorldReactorConstraint2017());
        MultiExperiment exptSyst(experiments);

        std::vector <const IFitVar*> fitvarsUO = {&kFitSinSqTheta23UO, &kFitDmSq32Scaled};
        std::vector <const IFitVar*> fitvarsLO = {&kFitSinSqTheta23LO, &kFitDmSq32Scaled};

        RestartCalculator(calc, anaName, nh);
        MinuitFitter fitter(&exptSyst, fitvarsUO);
        fitter.Fit(calc, IFitter::kQuiet);

        g->SetPoint(g->GetN(),
                    kFitSinSqTheta23UO.GetValue(calc),
                    kFitDmSq32Scaled.GetValue(calc));

        std::cout<< "\nsigma: " << sigma <<std::endl;
        std::cout << kFitSinSqTheta23.GetValue(calc) << " " << kFitDmSq32Scaled.GetValue(calc) << std::endl;

        RestartCalculator(calc, anaName, nh);
        MinuitFitter fitterLO(&exptSyst, fitvarsLO);
        fitterLO.Fit(calc, IFitter::kQuiet);

        gLO->SetPoint(gLO->GetN(),
                      kFitSinSqTheta23LO.GetValue(calc),
                      kFitDmSq32Scaled.GetValue(calc));

        std::cout << kFitSinSqTheta23LO.GetValue(calc) << " " << kFitDmSq32Scaled.GetValue(calc) << std::endl;
        prog.SetProgress((sigma+1)/2.);

      }

      g->SetLineColor(color[graphId]);
      g->SetLineStyle(styleline[graphId]); 
      gLO->SetLineColor(color[graphId]);
      gLO->SetLineStyle(styleline[graphId]);
      graphId++;

      g->SetTitle((s->LatexName()).c_str());
      g->SetFillColor(10);
      g->SetMarkerStyle(0);
      g->SetLineWidth(3);
      gLO->SetTitle((s->LatexName()).c_str());
      gLO->SetFillColor(10);
      gLO->SetMarkerStyle(0);
      gLO->SetLineWidth(3);

      leg-> AddEntry(g, (s->LatexName()).c_str(), "l");
      leg-> SetFillColor(0);
      leg-> SetFillStyle(0);

      canvas-> cd();
      g  ->Draw("l");
      gLO->Draw("l");
      surfaceSens->DrawBestFit(kBlack,kFullStar);

      canvas_zoom-> cd();
      g  ->Draw("l");
      gLO->Draw("l");
      surfaceSens->DrawBestFit(kBlack,kFullStar);
    }

    canvas-> cd();
    leg->Draw();
    canvas->SaveAs((out+".root"));
    canvas->SaveAs((out+".pdf"));
    canvas->SaveAs((out+".png"));
    std::ofstream file;
    file.open (out + ".txt");
    file << "90% confidence level sensitivity at NOvA's combined neutrino + antineutrino 2018 stats only best fit. Each point in the coloured lines represent the new best fit when a single systematic is applied to the fit: from start to end, a shift to the systematic is applied to cover a +-1$\sigma$ range.";
    file.close();
    canvas->Close();

    canvas_zoom-> cd();
    leg->Draw();
    canvas_zoom->SaveAs((out_zoom+".root"));
    canvas_zoom->SaveAs((out_zoom+".pdf"));
    canvas_zoom->SaveAs((out_zoom+".png"));
    std::ofstream file_zoom;
    file_zoom.open (out_zoom + ".txt");
    file_zoom << "Zoom in the 90% confidence level sensitivity at NOvA's combined neutrino + antineutrino 2018 stats only best fit. Each point in the coloured lines represent the new best fit when a single systematic is applied to the fit: from start to end, a shift to the systematic is applied to cover a +-1$\sigma$ range.";
    file_zoom.close();
    canvas_zoom->Close();

  }


} // End of function

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

Definition at line 28 of file make_starplots.C.

References std::asin(), delta_fhc_stats, delta_rhc_stats, delta_rhcfhc_stats, ana::ResetOscCalcToDefault(), osc::_IOscCalcAdjustable< T >::SetDmsq32(), osc::_IOscCalcAdjustable< T >::SetTh23(), sin_fhc_stats, sin_rhc_stats, sin_rhcfhc_stats, and std::sqrt().

Referenced by make_starplots().

{
  double sinHere = sin_rhcfhc_stats;
  double deltaHere = delta_rhcfhc_stats;
  // naively changing hierarchy
  // this isnt totally right
  double deltaSign = 1.;
  if(nh==false) deltaSign = -1.;

  if(anaName=="rhc"){
    sinHere = sin_rhc_stats;
    deltaHere = delta_rhc_stats*deltaSign;
  }
  if(anaName=="fhc"){
    sinHere = sin_fhc_stats;
    deltaHere = delta_fhc_stats*deltaSign;
  }
  ana::ResetOscCalcToDefault(calc);
  calc->SetTh23(asin(sqrt(sinHere)));
  calc->SetDmsq32(deltaHere);
}

Variable Documentation

std::vector<Color_t> color = {kMagenta, kGreen+1, kOrange+1, kBlue-7, kCyan+1, kRed-4, kBlue+1, kGray+1, kOrange+3}

Definition at line 25 of file make_starplots.C.

Referenced by make_starplots().

std::string dataName = "fakedata"

Definition at line 22 of file make_starplots.C.

Referenced by doTemplateFit(), getTimePeakPlots(), MakeValidation(), makeXSecPlots1D(), makeXSecPlots2D(), makeXSecPlots_TemplateFit(), and plot_contprof().

std::vector<Style_t> styleline = {kSolid, kSolid, kSolid, kSolid, kSolid, kSolid, kSolid, kSolid, kSolid}

Definition at line 24 of file make_starplots.C.

Referenced by make_starplots(), and plot_contprof().