futureSig_reach_singlePOTcombo_syst.C

Go to the documentation of this file.

#include "CAFAna/Analysis/CountingExperimentStatistics.h"
#include "CAFAna/Analysis/Exposures.h"
#include "CAFAna/Fit/MinuitFitter.h"
#include "CAFAna/Analysis/Style.h"
#include "CAFAna/Vars/FitVars.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/LoadFromFile.h"
#include "CAFAna/Core/Utilities.h"
#include "CAFAna/Analysis/Calcs.h"
#include "CAFAna/Experiment/CountingExperiment.h"
#include "CAFAna/Experiment/SingleSampleExperiment.h"
#include "CAFAna/Experiment/MultiExperiment.h"
#include "CAFAna/Experiment/Dmsq32Constraint.h"
#include "CAFAna/Experiment/AtmConstraint.h"
#include "CAFAna/Experiment/ReactorExperiment.h"
#include "CAFAna/Core/Progress.h"
#include "CAFAna/Prediction/IPrediction.h"
#include "CAFAna/Prediction/PredictionInterp.h"
#include "CAFAna/Prediction/PredictionExtrap.h"
#include "CAFAna/Prediction/PredictionCombinePeriods.h"
#include "CAFAna/Systs/SystComponentScale.h"
#include "OscLib/OscCalcPMNSOpt.h"
#include "OscLib/OscCalcDumb.h"
#include "Utilities/func/MathUtil.h"

using namespace ana;

#include "TFile.h"
#include "TH2.h"
#include "TMath.h"
#include "TCanvas.h"
#include "TGraph.h"
#include "TLatex.h"
#include "TLegend.h"
#include "TBox.h"
#include "TLine.h"
#include "TSystem.h"
#include "TColor.h"

#include <iostream>
#include <iomanip>
#include "Utilities/rootlogon.C"

#include "futureSig_systs.h"

osc::IOscCalcAdjustable* MyDefaultOscCalc();
void ResetOscCalcToMyDefault(osc::IOscCalcAdjustable* calc);

std::string MakeSuffix(const std::string PID, const std::string fakeHieDelta,
                       const double ssth23,   const std::string potID, 
                       const std::string systType, const std::string version);

TGraph * GetSignGraph(TH1D* hist, std::vector <int> potYears, 
                      std::vector <std::pair <int, int> > potCombos,
                      Int_t color, Int_t style);
void PaintReachCanvas(TString plottitle_split);


void futureSig_reach_singlePOTcombo_syst(const bool makeFile = false,
                                         const std::string version = "v01",
                                         const std::string fakeHieDelta = "NH3pi2",
                                         const std::string fakeOctant = "UO",
                                         const std::string systType = "_12",
                                         const std::string potBaseID = "comboA",
                                         const std::string potCompID = "none")
{
  bool stats = (systType=="_stats");  
  double ssth23 = ((fakeOctant=="UO")?0.625:0.403);
  if(fakeOctant =="MAX") ssth23 = 0.5;
  
  auto suffix = MakeSuffix("CVN",fakeHieDelta,ssth23,potBaseID,systType,version);

  TString filename = "root_reach/hist_novareach_"+ suffix;
  TString filenameCompare = ("root_reach/hist_novareach_"+MakeSuffix("CVN",fakeHieDelta,ssth23,potCompID,systType,version));
  TString plotfolder = version+"/";
  gSystem->mkdir("root_reach");
  gSystem->mkdir(plotfolder);

  
  std::string hieDelStr = "Normal #delta_{CP}=3#pi/2";
  if(fakeHieDelta == "IH3pi2") hieDelStr = "Inverted #delta_{CP}=3#pi/2";
  
  TString plottitle = TString::Format("%s, sin^{2}#theta_{23}=%.3f, #Deltam^{2}_{32}=2.5#times10^{-3}eV^{2}, sin^{2}#theta_{13}=0.022",hieDelStr.c_str(),ssth23);
  TString plottitle_split = TString::Format("#splitline{%s, sin^{2}#theta_{23}=%.3f}{#Deltam^{2}_{32}=2.5#times10^{-3}eV^{2}, sin^{2}#theta_{13}=0.022}",hieDelStr.c_str(),ssth23);

  std::vector <int> potYears =                     {2016,2017, 2018,  2019,    2020,  2021,   2022,   2023,   2024};
  std::vector <std::pair <int, int> > potComboT = {{6,0},{9,3}}; // for tests
  std::vector <std::pair <int, int> > potComboA = {{6,0},{9,3},{9,9},{12,12},{15,15},{18,18},{21,21},{24,24},{27,27}};
  std::vector <std::pair <int, int> > potComboB = {{6,0},{12,0},{18,0},{24,0},{30,0},{36,0},{42,0},{48,0},{54,0}};
  std::vector <std::pair <int, int> > potComboC = {{6,0},{9,3},{9,9},{15,9},{18,12},{21,15},{24,18},{27,21},{30,24}};
  std::vector <std::pair <int, int> > potComboD = {{6,0},{9,3},{9,9},{13,12},{17,15},{21,18},{25,21},{29,24},{33,27}}; 
  std::vector <std::pair <int, int> > potComboE = {{6,0},{9,3},{9,9},{13,12},{16,16},{20,19},{23,23},{27,26},{30,30}};
  
  std::vector <std::pair <int, int> > potComboBase;
  std::vector <std::pair <int, int> > potComboCompare;

  if(potBaseID == "comboT") potComboBase = potComboT;
  if(potBaseID == "comboA") potComboBase = potComboA;
  if(potBaseID == "comboB") potComboBase = potComboB;
  if(potBaseID == "comboC") potComboBase = potComboC;
  if(potBaseID == "comboD") potComboBase = potComboD;
  if(potBaseID == "comboE") potComboBase = potComboE;

  if(potCompID == "comboA") potComboCompare = potComboA;
  if(potCompID == "comboB") potComboCompare = potComboB;
  if(potCompID == "comboC") potComboCompare = potComboC;
  if(potCompID == "comboD") potComboCompare = potComboD;
  if(potCompID == "comboE") potComboCompare = potComboE;
  
  bool compare = (potCompID!="none");
  
  if(makeFile){
    TFile * file = new TFile (filename+".root","recreate");

    //This format here looks weird, but matches the reach_optimization 
    TH2D * sign_Hie = new TH2D ("hie","hie; Total POT (#times10^{20}); RHC POT (#times10^{20})",
                                67,5.5,72.5,67,-0.5,66.5);
    TH2D * sign_Max = new TH2D ("max","max; Total POT (#times10^{20}); RHC POT (#times10^{20})",
                                67,5.5,72.5,67,-0.5,66.5);  
    TH2D * sign_Oct = new TH2D ("oct","oct; Total POT (#times10^{20}); RHC POT (#times10^{20})",
                                67,5.5,72.5,67,-0.5,66.5);  
    TH2D * sign_Cpv = new TH2D ("cpv","cpv; Total POT (#times10^{20}); RHC POT (#times10^{20})",
                                67,5.5,72.5,67,-0.5,66.5);  

    TString dir = "root_predictions/";
    TFile * file1 = new TFile( dir+"predInterp_nue_FHC1_"+version+".root","read");   
    TFile * file2 = new TFile( dir+"predInterp_nue_FHC2_"+version+".root","read");   
    TFile * file3 = new TFile( dir+"predInterp_nue_RHC_"+version+".root","read");   
    TFile * file4 = new TFile( dir+"predInterp_numu_FHC1_full_"+version+".root","read");   
    TFile * file5 = new TFile( dir+"predInterp_numu_FHC2_full_"+version+".root","read");   
    TFile * file6 = new TFile( dir+"predInterp_numu_RHC_full_"+version+".root","read");   
    
    auto predFHC1 = 
      PredictionInterp::LoadFrom(file1, "nue_fhcPrediction_p123b").release();
    auto predFHC2 = 
      PredictionInterp::LoadFrom(file2, "nue_fhcPrediction_p3c").release();
    auto predRHC  = 
      PredictionInterp::LoadFrom(file3, "nue_rhcPrediction").release();
    
    auto prudFHC1 = 
      PredictionInterp::LoadFrom(file4, "numu_fhcPrediction_p123b_full").release();
    auto prudFHC2 = 
      PredictionInterp::LoadFrom(file5, "numu_fhcPrediction_p3c_full").release();
    auto prudRHC  = 
      PredictionInterp::LoadFrom(file6, "numu_rhcPrediction_full").release();

    delete file1, file2, file3, file4, file5, file6;

    for(auto potCombo:potComboBase){
      
      int fhc = potCombo.first;
      int rhc = potCombo.second;
      int totpot = fhc+rhc;
      
      std::cerr << "\nFinding best fit points for " << fhc << "p" << rhc 
                << "ssth23 " << ssth23 <<"\n\n";
      
      double p1POT = kSecondAnaPeriod1POT + kSecondAnaPeriod2POT +  kSecondAnaEpoch3bPOT;
      double p2POT = fhc*1E20 - 3.41E20;
      double potFHC = p1POT + p2POT;
      double potRHC = rhc*1E20;
      
      auto predFHC = new PredictionCombinePeriods({{predFHC1,p1POT},{predFHC2,p2POT}});
      auto prudFHC = new PredictionCombinePeriods({{prudFHC1,p1POT},{prudFHC2,p2POT}});
        
      osc::IOscCalcAdjustable* calc;
      calc = MyDefaultOscCalc();  
      calc->SetTh23(asin(sqrt(ssth23)));
      if(fakeHieDelta == "IH3pi2") calc->SetDmsq32(-1*fabs(calc->GetDmsq32()));
      
      auto obsFHC = predFHC->Predict(calc).FakeData(potFHC);
      auto obsRHC = predRHC->Predict(calc).FakeData(potRHC);
      auto ubsFHC = prudFHC->Predict(calc).FakeData(potFHC);
      auto ubsRHC = prudRHC->Predict(calc).FakeData(potRHC);

      auto exptFHC = new SingleSampleExperiment(predFHC, obsFHC);
      auto exptRHC = new SingleSampleExperiment(predRHC, obsRHC);

      auto uxptFHC = new SingleSampleExperiment(prudFHC, ubsFHC);
      auto uxptRHC = new SingleSampleExperiment(prudRHC, ubsRHC); 

      MultiExperiment * expt;
      if(rhc > 0 ) expt = new MultiExperiment({exptFHC, uxptFHC, exptRHC, uxptRHC,
            WorldReactorConstraint2015()});
      else expt = new MultiExperiment({exptFHC, uxptFHC, WorldReactorConstraint2015()});
      
      ////////////////////////////////////////////////////////////      
      //Systematics
      ////////////////////////////////////////////////////////////
      auto mySysts = GetSystsList(systType);
      //set up correlations
      //0 nue, 1 numu; 2 nue RHC, 3 numu rhc
      std::vector< std::pair< const ISyst *, const ISyst * >> corr_nue_FHC;
      std::vector< std::pair< const ISyst *, const ISyst * >> corr_numu_FHC;
      std::vector< std::pair< const ISyst *, const ISyst * >> corr_nue_RHC;
      std::vector< std::pair< const ISyst *, const ISyst * >> corr_numu_RHC;

      // nue FHC should not have any numu or wrong sign
      //similar logic for all
      for (auto nosys:systsNumuRHC)     corr_nue_FHC.push_back({nosys,NULL});
      expt->SetSystCorrelations(0,corr_nue_FHC);

      for (auto nosys:systsNueRHC)      corr_numu_FHC.push_back({nosys,NULL});
      expt->SetSystCorrelations(1,corr_numu_FHC);

      if(rhc>0){
        for (auto nosys:systsNumuFHC)   corr_nue_RHC.push_back({nosys,NULL});
        expt->SetSystCorrelations(2,corr_nue_RHC);
        
        for (auto nosys:systsNueFHC)    corr_numu_RHC.push_back({nosys,NULL});
        expt->SetSystCorrelations(3,corr_numu_RHC);
      } 

      if(stats)    {  mySysts.clear(); std::cout << "WARNING STATS ONLY ACTIVATED \n"; }

     ////////////////////////////////////////////////////////////
     // Fit
     ////////////////////////////////////////////////////////////
     
      MinuitFitter fit23(expt, {&kFitDeltaInPiUnits,&kFitSinSqTheta23,
                          &kFitDmSq32,&kFitSinSq2Theta13},mySysts);

      double minchi_NH_UO=1E20;
      double minchi_NH_LO=1E20;
      double minchi_IH_UO=1E20;
      double minchi_IH_LO=1E20;
      
      for(int hie:{-1,+1}){
        for(double seed:{0.3,0.7}){
          double minchi_temp = 1E20;
          calc = MyDefaultOscCalc();  
          kFitSinSqTheta23.SetValue(calc, seed);
          calc->SetDmsq32(hie*fabs(calc->GetDmsq32()));
          minchi_temp = std::min(minchi_temp, fit23.Fit(calc)->EvalMetricVal());
          const double dcp = calc->GetdCP();
          for(int n = 1; n <= 3; ++n){
            calc->SetdCP(dcp+n*M_PI/2);
            minchi_temp = std::min(minchi_temp, fit23.Fit(calc)->EvalMetricVal());
          }
          if(hie>0 && seed<0.5) minchi_NH_LO = minchi_temp;
          if(hie>0 && seed>0.5) minchi_NH_UO = minchi_temp;
          if(hie<0 && seed<0.5) minchi_IH_LO = minchi_temp;
          if(hie<0 && seed>0.5) minchi_IH_UO = minchi_temp;
          
          std::cout << "\n \n Found minchi " << (hie>0?"NH ": "IH ") 
                    << (seed<0.5? "LO ":"UO ") << minchi_temp << "\n\n";
        }
      }
      
      double minchi_NH = std::min(minchi_NH_UO,minchi_NH_LO);
      double minchi_IH = std::min(minchi_IH_UO,minchi_IH_LO);
      double minchi_UO = std::min(minchi_NH_UO,minchi_IH_UO);
      double minchi_LO = std::min(minchi_NH_LO,minchi_IH_LO);
      double minchi_best = std::min(minchi_NH,minchi_IH);
      
      //fixed ssth23
      MinuitFitter fit23Max(expt, {&kFitDeltaInPiUnits,&kFitDmSq32,&kFitSinSq2Theta13},mySysts);
      double minchi_Max=1E20;
      
      for(int hie:{-1,+1}){
        calc = MyDefaultOscCalc();  
        kFitSinSqTheta23.SetValue(calc, 0.5);
        calc->SetDmsq32(hie*fabs(calc->GetDmsq32()));
        minchi_Max = std::min(minchi_Max,fit23Max.Fit(calc)->EvalMetricVal());
        const double dcp2 = calc->GetdCP();
        for(int n = 1; n <= 3; ++n){
          calc->SetdCP(dcp2+n*M_PI/2);
          minchi_Max = std::min(minchi_Max, fit23Max.Fit(calc)->EvalMetricVal());
        }
      }
      // fixed dCP 
      MinuitFitter fit23CP(expt, {&kFitSinSqTheta23,&kFitDmSq32,&kFitSinSq2Theta13},mySysts);
      double minchi_CPV=1E20;
      
      for(int hie:{-1,+1}){
        for(double seed:{0.3,0.7}){
          for(double dCP:{0,1}){
            calc = MyDefaultOscCalc();  
            kFitSinSqTheta23.SetValue(calc, seed);
            kFitDeltaInPiUnits.SetValue(calc, dCP);
            calc->SetDmsq32(hie*fabs(calc->GetDmsq32()));
            minchi_CPV = min(minchi_CPV, fit23CP.Fit(calc)->EvalMetricVal());
          }
        }
      }
      
      std::cout << "\n \n Found minchi NH " << minchi_NH 
                << " IH " << minchi_IH
                << " UO " << minchi_UO
                << " LO " << minchi_LO
                << " Max " << minchi_Max
                << " CPV " << minchi_CPV
                << std::endl;
      
      sign_Hie->SetBinContent(sign_Hie->FindBin(totpot,rhc),fabs(minchi_IH - minchi_NH));
      sign_Oct->SetBinContent(sign_Oct->FindBin(totpot,rhc),fabs(minchi_LO - minchi_UO));
      sign_Max->SetBinContent(sign_Max->FindBin(totpot,rhc),minchi_Max - minchi_best);
      sign_Cpv->SetBinContent(sign_Cpv->FindBin(totpot,rhc),minchi_CPV - minchi_best);
    } //end potCombo

    sign_Hie->SetTitle("#Delta #chi^{2} Wrong hie. rejection - "+ plottitle);
    sign_Oct->SetTitle("#Delta #chi^{2} Octant resolution - "+plottitle);
    sign_Max->SetTitle("#Delta #chi^{2} Maximal mix. rejection - "+plottitle);
    sign_Cpv->SetTitle("#Delta #chi^{2} CP violation - "+plottitle);
    
    sign_Hie->SetOption("colz");
    sign_Oct->SetOption("colz");
    sign_Max->SetOption("colz");
    sign_Cpv->SetOption("colz");
    
    file->cd();
    sign_Hie->Write();
    sign_Oct->Write();
    sign_Max->Write();
    sign_Cpv->Write();

    std::cout << "Saved significances to " << file->GetName() << std::endl;
    return;
  }//end makefile


  else{
    std::vector <TString> sign_label={"max","hie","oct","cpv"};
    std::vector <TString> sign_name={"Max. mixing","Hierarchy","Octant","CPV"};
    
    Int_t potColors[]={
      TColor::GetColor("#6c71c4"), TColor::GetColor("#268bd2"),
      TColor::GetColor("#2aa198"), TColor::GetColor("#859900"),
      TColor::GetColor("#b58900"), TColor::GetColor("#cb4b16"),
      TColor::GetColor("#dc322f"), TColor::GetColor("#d33682")
    };
    Int_t lineStyle[]={1,7,10,9};    
    
    ////////////////////////////////////////////////////////////
    
    PaintReachCanvas(plottitle_split);
    
    TLegend * leg = new TLegend (0.12, 0.52, 0.47, 0.85);
    leg->SetTextSize(0.055);
    leg->SetFillStyle(0);
    leg->SetMargin(1.3*leg->GetMargin());
    leg->SetHeader("NOvA joint #nu_{e}+#nu_{#mu}");
    
    //Get graphs and plot
    TFile * file = new TFile (filename+".root","read");    
    int sigIdx=0;
    for(auto sign:sign_label){
      auto hist= (TH1D*) file->Get(sign);
      auto gr = GetSignGraph(hist, potYears,potComboBase,potColors[sigIdx*2],lineStyle[sigIdx]);
      gr->Draw("L same");
      leg->AddEntry(gr,sign_name[sigIdx],"l");
      sigIdx++;
    }
    //Compare to a different POT combo
    if(compare){
      TFile * fileCompare = new TFile (filenameCompare+".root","read");    
      if(fileCompare->IsZombie()) {std::cout << "can t find " << fileCompare->GetName() 
                                             << "\n\n"; return;}
      sigIdx=0;
      for(TString sign:{"hie","oct","max","cpv"}){
        auto hist= (TH1D*) fileCompare->Get(sign);
        Int_t trans_color = TColor::GetColorTransparent(potColors[sigIdx*2], 0.3);
        auto gr = GetSignGraph(hist, potYears,potComboCompare,trans_color,kSolid);
        gr->Draw("L same");
        sigIdx++;
      }
    }//end compare
    
    leg->Draw();
    if(!compare) gPad->Print(plotfolder+"sign_reach_"+suffix+".pdf");
    else gPad->Print(plotfolder+"sign_reach_"+suffix+"vs"+potCompID+".pdf");      
    
    return;    
  }//end !makefile
  
}//end all

////////////////////////////////////////////////////////////
// Osc. parameters as chosen by conveners
////////////////////////////////////////////////////////////

void ResetOscCalcToMyDefault(osc::IOscCalcAdjustable* calc)
{
  calc->SetL(810);
  calc->SetRho(2.84); // g/cm^3                                                                                                                                                                        
  calc->SetDmsq21(7.4e-5);
  calc->SetTh12(asin(sqrt(.30)));
  calc->SetTh23(M_PI/4);
  calc->SetDmsq32(2.5e-3); 
  calc->SetTh13(asin(sqrt(.022)));
  calc->SetdCP(1.5*M_PI);
}

osc::IOscCalcAdjustable* MyDefaultOscCalc()
{
  osc::IOscCalcAdjustable* ret = new osc::OscCalcPMNSOpt;
  ResetOscCalcToMyDefault(ret);
  return ret;
}

////////////////////////////////////////////////////////////
// Osc. parameters as chosen by conveners
////////////////////////////////////////////////////////////

TGraph * GetSignGraph(TH1D* hist, std::vector <int> potYears, 
                      std::vector <std::pair <int, int> > potCombos,
                      Int_t color, Int_t style){
  TGraph * gr = new TGraph();
  int yearIdx = 0;
  for(auto potCombo:potCombos){
    int fhc = potCombo.first;
    int rhc = potCombo.second;
    int totpot = fhc+rhc;
    double this_chisq = hist->GetBinContent(hist->FindBin(totpot,rhc));
    double this_sign = sqrt(std::max(this_chisq,0.));
    gr->SetPoint(gr->GetN(),potYears[yearIdx], this_sign);
    yearIdx++;  
  }//end years
  gr->SetLineWidth(3);
  gr->SetLineColor(color);
  gr->SetLineStyle(style);                                
  return gr;
}
////////////////////////////////////////////////////////////
void PaintReachCanvas(TString plottitle_split){
  new TCanvas;
  gPad->SetTopMargin(0.14);
  TH2D * axes = new TH2D ("ax",";Year;Significance (#sigma)",18,2015.5,2024.5,100,0,5.9);
  axes->GetXaxis()->CenterTitle();
  axes->GetYaxis()->CenterTitle();
  axes->GetYaxis()->SetTitleOffset(0.6);
  axes->Draw();
  TLatex * ltx = new TLatex();
  ltx->SetTextAlign(13);
  ltx->SetTextSize(axes->GetXaxis()->GetLabelSize());
  ltx->DrawLatexNDC(0.12,0.99,plottitle_split.Data() );
  
  TLatex* ltx2 = new TLatex();
  ltx2->SetTextSize(1.2*axes->GetXaxis()->GetLabelSize());
  ltx2->DrawLatexNDC(0.35,0.16,"#splitline{2016 analysis techniques with projected}{systematic uncertainty improvements}");
  gPad->SetGridx();
  gPad->SetGridy();
  Simulation();
}
     
std::string MakeSuffix(const std::string PID, const std::string fakeHieDelta,
                       const double ssth23,   const std::string potID, 
                       const std::string systType, const std::string version)
{
  std::string suffix= PID;
  suffix+="_"+fakeHieDelta;
  if(ssth23<1)  suffix+=(std::to_string(ssth23)).replace(0,2,"_0p").erase(6,4);
  else suffix+="_0pXX";
  bool stats = (systType=="_stats");
  suffix+=("_"+ potID+(stats?"":"_syst")+systType);
  suffix+="_"+version;
  return suffix;
}