syst_table_fit.C

Go to the documentation of this file.

#include "3FlavorAna/Ana2020/joint_fit_2020_loader_tools.h"
#include "CAFAna/Experiment/AtmConstraint.h"
#include "CAFAna/Experiment/MultiExperiment.h"
#include "CAFAna/Experiment/ReactorExperiment.h"
#include "CAFAna/Experiment/SolarConstraints.h"
#include "3FlavorAna/Systs/3FlavorAna2020Systs.h"
#include "CAFAna/Fit/Fit.h"
#include "3FlavorAna/Systs/NumuSysts.h"
#include "CAFAna/Systs/Systs.h"
#include "3FlavorAna/Systs/NueAcceptSysts.h"
#include "3FlavorAna/Systs/MichelTaggingSyst.h"
#include "CAFAna/Analysis/Plots.h"

#include <fstream>

using namespace ana;

void GetLimits( const TGraph* plot, float& hi, float& low )

{
  const int N = plot->GetN();
  low = plot->GetX()[0];
  hi = plot->GetX()[N-1];
  float min = plot->GetY()[0];
  float bf = plot->GetX()[0];
  int bfp = 0;

  for ( int i = 0; i < N; ++i )
  {
    float y = plot->GetY()[i];
    if ( y < min ) { min = y; bf = plot->GetX()[i]; bfp = i; }
  }

  std::cout << "min = " << min << ", bf = " << bf << ", bfp = " << bfp << std::endl;
  const float chi2th = 1.0; // Chi squared threshold

  for ( int i = bfp; i > 0; --i )
  {
    float y1 = plot->GetY()[i-1];
    float y2 = plot->GetY()[i];
    if ( y1 > chi2th && y2 < chi2th )
    { 
      float x1 = plot->GetX()[i-1];
      float x2 = plot->GetX()[i];
      low = x1 + ( chi2th - y1 ) * ( x2 - x1 ) / ( y2 - y1 );
      std::cout << "Crossed chisq of " << chi2th << ", descending, with x value of: " << low << std::endl;
      break;
    }
  }

  for ( int i = bfp; i < N-1; ++i )
  {
    float y1 = plot->GetY()[i];
    float y2 = plot->GetY()[i+1];
    if ( y1 < chi2th && y2 > chi2th )
    {
      float x1 = plot->GetX()[i];
      float x2 = plot->GetX()[i+1];
      hi = x1 + ( chi2th - y1 ) * ( x2 - x1 ) / ( y2 - y1 );
      std::cout << "Crossed chisq of " << chi2th << ", ascending, with x value of: " << hi << std::endl;
      break;
    }
  }
}


void simulation()
{
  TLatex* prelim = new TLatex(.9, .93, "NOvA Simulation");
  prelim->SetTextColor(kGray+1);
  prelim->SetNDC();
  prelim->SetTextSize(1.5/30.);
  prelim->SetTextAlign(32);
  prelim->Draw();
}

void WriteUncertainty( ofstream& txt, double plus, double minus )
{ 
  txt<< " & +" << plus << " / -" << minus;
  //txt << Form( " & +%.3e / -%.3e", plus, minus );
}

typedef std::vector<std::pair<const ISyst*,const ISyst*> > Correlations;
typedef std::pair<SingleSampleExperiment*,Correlations> ExptCorr;

void AddExptCorrNue( std::vector<ExptCorr>& exptCorrs, osc::IOscCalc* calc, std::string beam, bool ptExtrap, bool fakeNDData )
{
  bool isFHC = beam == "fhc";
  bool corrSysts = true;
  bool isNue = true;
  double POT = isFHC ? kAna2020FHCPOT : kAna2020RHCPOT;
  std::string decomp = isFHC ? "combo" : "prop";
  if ( !ptExtrap ) decomp += "_noPt";

  auto pred = GetNuePrediction2020( decomp, DefaultOscCalc(), corrSysts, beam, fakeNDData );
  auto cosmics = GetNueCosmics2020( beam );
  Spectrum* fakeData = GetFakeData( pred, calc, POT, cosmics.first );
  SingleSampleExperiment* expt = new SingleSampleExperiment( pred, *fakeData, *cosmics.first, cosmics.second, true );
  exptCorrs.emplace_back( expt, GetCorrelations( isNue, isFHC, ptExtrap ) );
}

void AddExptCorrNumu( std::vector<ExptCorr>& exptCorrs, osc::IOscCalc* calc, std::string beam, bool ptExtrap, bool fakeNDData )
{ 
  bool isFHC = beam == "fhc";
  bool useSysts = true;
  bool isNue = false; 
  double POT = isFHC ? kAna2020FHCPOT : kAna2020RHCPOT;
  int nEhadQuants = 4;
  std::string decomp = ptExtrap ? "" : "noPt";

  auto preds = GetNumuPredictions2020( nEhadQuants, decomp, DefaultOscCalc(), useSysts, beam, fakeNDData );
  auto cosmics = GetNumuCosmics2020( nEhadQuants, beam );
  auto corrs = GetCorrelations( isNue, isFHC, ptExtrap );

  for ( int i = 0; i < nEhadQuants; ++i )
  {
    Spectrum* fakeData = GetFakeData( preds[i], calc, POT, cosmics[i].first );
    SingleSampleExperiment* expt = new SingleSampleExperiment( preds[i], *fakeData, *cosmics[i].first, cosmics[i].second, true );
    exptCorrs.emplace_back( expt, corrs );
  }
}

void syst_table_fit_new( const TString option )
{

  const bool full     = option.Contains( "full" );
  const bool ptExtrap = option.Contains( "ptExtrap" );
  const bool fakeNDData = option.Contains( "fakeND" );
  const bool seeds    = option.Contains( "seeds" );

  std::string outStr = "";
  if ( ptExtrap ) outStr += "ptExtrap_";
  if ( fakeNDData ) outStr += "fakeND_";
  if ( seeds )    outStr += "withSeeds_";
  outStr += full ? "full" : "short";
  //std::cout << "outStr = " << outStr << std::endl; 

  std::string cwd = gSystem->GetWorkingDirectory();
  std::string outDir = cwd + "/syst_tables";
  if ( gSystem->AccessPathName( outDir.c_str() ) ) gSystem->mkdir( outDir.c_str() );

  // 2019 Best Fit
  //std::string calcfile = "/nova/ana/nu_e_ana/Ana2019/Results/BestFits/bestfits_joint_realData_both_systs.root";
  //TFile *calcf = new TFile(calcfile.c_str());
  //auto calctemp = LoadFrom<osc::IOscCalc>(calcf, "NH_UO_calc").release();
  //auto calc = (osc::IOscCalcAdjustable*)calctemp->Copy();

  //double dmsq32 = kFitDmSq32.GetValue(calc);
  //double ssth23 = kFitSinSqTheta23.GetValue(calc);
  //double dCP = kFitDeltaInPiUnits.GetValue(calc);
  // 2020 Best Fit
  /*double dmsq32 = 2.406e-3;
  double ssth23 = .568;
  double ss2th13 = .085;
  double dCP = .820;*/
  
  // 2020 Best Fit
  std::string calcfile = "/nova/ana/3flavor/Ana2020/Fits/RealData/WithSeeds/bestfits_joint_realData_both_systs.root";

  TFile *calcf = new TFile(calcfile.c_str());
  auto calctemp = LoadFrom<osc::IOscCalc>(calcf, "NH_UO_calc").release();
  auto calc = (osc::IOscCalcAdjustable*)calctemp->Copy();

  double dmsq32 = kFitDmSq32.GetValue(calc);
  double ssth23 = kFitSinSqTheta23.GetValue(calc);
  double dCP = kFitDeltaInPiUnits.GetValue(calc);
  double ss2th13 = kFitSinSq2Theta13.GetValue(calc);

  std::cout << "Best fit points:" << std::endl;
  std::cout << "dmsq32 = " << dmsq32 << std::endl;
  std::cout << "ssth23 = " << ssth23 << std::endl;
  std::cout << "dCP = " << dCP << std::endl;

  // Experiments and correlations
  std::vector<ExptCorr> exptCorrs;
  AddExptCorrNue( exptCorrs, calc, "fhc", ptExtrap, fakeNDData );
  AddExptCorrNue( exptCorrs, calc, "rhc", ptExtrap, fakeNDData );
  AddExptCorrNumu( exptCorrs, calc, "fhc", ptExtrap, fakeNDData );
  AddExptCorrNumu( exptCorrs, calc, "rhc", ptExtrap, fakeNDData );

  auto exptThis = new MultiExperiment();
  for ( unsigned int i = 0; i < exptCorrs.size(); ++i )
  {
    std::cout << "Adding experiment " << i+1 << std::endl;
    exptThis->Add( exptCorrs[i].first );
    exptThis->SetSystCorrelations( i, exptCorrs[i].second );
  }
  exptThis->Add( WorldReactorConstraint2019() );

  // All Joint Systs
  bool smallgenie = true;
  bool isFit = true;
  std::vector<const ISyst*> allSysts = get3FlavorAna2020AllSysts( EAnaType2020::k3FlavorAna, smallgenie, BeamType2020::kBoth, isFit, ptExtrap );

  // Syst Categories
  std::vector<const ISyst*> xsecSysts;
  Add3FlavorAna2020XSecSysts( xsecSysts );
  std::vector<const ISyst*> fluxSysts;
  Add3FlavorAna2020BeamSysts( fluxSysts );
  std::vector<const ISyst*> calibSysts;
  Add3FlavorAna2020CalibSysts( calibSysts );
  std::vector<const ISyst*> lightSysts;
  Add3FlavorAna2020LightSysts( lightSysts );
  std::vector<const ISyst*> leptonRecoSysts;
  Add3FlavorAna2020MuEnergySysts( leptonRecoSysts );
  Add3FlavorAna2020LeptonAngleSysts( leptonRecoSysts, ptExtrap );
  std::vector<const ISyst*> neutronSysts;
  Add3FlavorAna2020NeutronSysts( neutronSysts );
  std::vector<const ISyst*> nearFarSysts;
  Add3FlavorAna2020NueAcceptSysts( nearFarSysts, EAnaType2020::k3FlavorAna, BeamType2020::kBoth, ptExtrap );
  Add3FlavorAna2020MichelTagSysts( nearFarSysts );
  Add3FlavorAna2020NormSysts( nearFarSysts );
  AddNonLoadable2020Systs( nearFarSysts );

  std::map< std::string, std::vector<const ISyst*> > mapSystVecs;
  mapSystVecs[ "Neutrino Cross Sections" ] = xsecSysts;
  mapSystVecs[ "Beam Flux"               ] = fluxSysts;
  mapSystVecs[ "Detector Calibration"    ] = calibSysts;
  mapSystVecs[ "Detector Response"       ] = lightSysts;
  mapSystVecs[ "Lepton Reconstruction"   ] = leptonRecoSysts;
  mapSystVecs[ "Neutron Uncertainty"     ] = neutronSysts;
  mapSystVecs[ "Near-Far Uncor."         ] = nearFarSysts;

  std::map  <const IFitVar *, std::vector <double> > th23_seeds   = {{&kFitDeltaInPiUnits, {dCP}}, { &kFitDmSq32, {dmsq32}}};
  std::map  <const IFitVar *, std::vector <double> > delta_seeds  = {{&kFitSinSqTheta23, {ssth23}}, { &kFitDmSq32, {dmsq32}}, {&kFitSinSq2Theta13, {.085}}};
  std::map  <const IFitVar *, std::vector <double> > dmsq32_seeds = {{&kFitSinSqTheta23, {ssth23}}, {&kFitDeltaInPiUnits, {dCP}}};
  if(seeds)
  {
    delta_seeds  = { {&kFitSinSqTheta23, {.466, .5, ssth23}}, { &kFitDmSq32, {dmsq32}}, {&kFitSinSq2Theta13, {.085}} };
    th23_seeds   = { {&kFitDeltaInPiUnits, {.0, .5, 1.0, 1.5}}, { &kFitDmSq32, {dmsq32}} };
    dmsq32_seeds = { {&kFitSinSqTheta23, {.466, .5, ssth23}}, {&kFitDeltaInPiUnits, {.0, .5, 1.0, 1.5}}};
  }
    
  // Define what to profile over
  const int steps = 60;
  struct ProfDef{
    const IFitVar * fitvar;
    double minx;
    double maxx;
    std::vector<const IFitVar *> profvars;
    std::map  <const IFitVar *, std::vector <double> >profseeds;
    //const SeedList profseeds;
    double bf;
    int oom;
    std::string latex;
    std::string root;
    std::string suffix;
  };

  std::vector<ProfDef> profiles = {
    {&kFitSinSqTheta23, 0.3, 0.7, {&kFitDmSq32, &kFitDeltaInPiUnits, &kFitSinSq2Theta13}, th23_seeds,
        ssth23, 0, "$\\mathrm{sin}^2\\theta_{23}$", "Uncertainty in sin^{2}#theta_{23}", "th23"},
    {&kFitDeltaInPiUnits, -.5, 2.2, {&kFitDmSq32, &kFitSinSqTheta23, &kFitSinSq2Theta13}, delta_seeds,
        dCP, 0, "$\\delta_{CP}/\\pi$", "Uncertainty in #delta_{CP}/#pi", "dcp"},
    {&kFitDmSq32, 2e-3, 3e-3, {&kFitDeltaInPiUnits, &kFitSinSqTheta23, &kFitSinSq2Theta13}, dmsq32_seeds,
        dmsq32, 3, "$\\vert\\Delta m^2_{32}\\vert$ ($\\times 10^{-3}~\\mathrm{eV}^2$)", "Uncertainty in #Delta m^{2}_{32} (#times10^{-3} eV^{2})", "dmsq"}
  };

  // outputs
  TFile* fOut = new TFile( ( outDir + "/syst_table_profiles_" + outStr + ".root" ).c_str(), "RECREATE" );

  ofstream txt;
  txt.open( ( outDir + "/syst_table_" + outStr + ".txt" ).c_str() );

  txt<<std::setprecision(2);

  // Profiles n things
  SetFakeCalc(calc, ssth23, dmsq32, dCP);
  kFitSinSq2Theta13.SetValue(calc, ss2th13);

  double range;
  float hi, low;
  std::map<std::string, double> statup;
  std::map<std::string, double> statdn;
  std::map<std::string,std::map<std::string,std::pair<double,double>>> bars;

  txt<< "Source of Uncertainty";
  for(auto prof: profiles) txt<< " & " << prof.latex;
  txt<< " \\\\" << std::endl;
  txt<< "\\hline" << std::endl;   

  // Stat only
  for(auto prof: profiles){
    SetFakeCalc(calc, ssth23, dmsq32, dCP);
    kFitSinSq2Theta13.SetValue(calc, ss2th13);
    
    TGraph* plot = (TGraph*)Profile(exptThis, calc,
                                prof.fitvar,
                                steps, prof.minx, prof.maxx, -1,
                                prof.profvars,
                {},
                prof.profseeds);
                
    plot -> Write((prof.suffix+"_stat").c_str());

    GetLimits(plot, hi, low);

    double up = hi-prof.bf;
    double dn = prof.bf-low;

    statup[prof.suffix] = up;
    statdn[prof.suffix] = dn;
  }

  // Systs
  if(!full){
    for(auto const &systVec: mapSystVecs) {
      if ( systVec.second.empty() )
      { 
        std::cout << "Skipping " << systVec.first << " with empty syst vector" << std::endl;
        continue;
      }
      txt<< systVec.first ;
      
      for(auto prof: profiles){
        SetFakeCalc(calc, ssth23, dmsq32, dCP);
    kFitSinSq2Theta13.SetValue(calc, ss2th13);
        
        TGraph* systplot = (TGraph*)Profile(exptThis, calc,
                                        prof.fitvar,
                                        steps, prof.minx, prof.maxx, -1,
                                        prof.profvars,
                                        systVec.second,
                    prof.profseeds);
        systplot -> Write((prof.suffix+"_" + systVec.first).c_str());
        
        GetLimits(systplot, hi, low);

        double plus = sqrt( fabs( (hi-prof.bf)*(hi-prof.bf) - (statup[prof.suffix]*statup[prof.suffix]) ) ) * pow(10,prof.oom);
        double minus = sqrt( fabs( (prof.bf-low)*(prof.bf-low) - (statdn[prof.suffix]*statdn[prof.suffix]) ) ) * pow(10,prof.oom);

        WriteUncertainty( txt, plus, minus );

        bars[prof.suffix][systVec.first] = {minus,plus};
      }
      txt<< " \\\\" <<std::endl;
    }    
  }

  if(full){
    for(const ISyst* syst: allSysts) {
      txt<< syst->LatexName() ;
      for(auto prof: profiles){
            SetFakeCalc(calc, ssth23, dmsq32, dCP);
        kFitSinSq2Theta13.SetValue(calc, ss2th13);
        
            TGraph* systplot = (TGraph*)Profile(exptThis, calc,
                                        prof.fitvar,
                                        steps, prof.minx, prof.maxx, -1,
                                        prof.profvars,
                                        {syst},
                    prof.profseeds);
                    
            systplot -> Write((prof.suffix+"_" + syst->ShortName()).c_str());
        
            GetLimits(systplot, hi, low);

        double plus = sqrt( fabs( (hi-prof.bf)*(hi-prof.bf) - (statup[prof.suffix]*statup[prof.suffix]) ) ) * pow(10,prof.oom);
        double minus = sqrt( fabs( (prof.bf-low)*(prof.bf-low) - (statdn[prof.suffix]*statdn[prof.suffix]) ) ) * pow(10,prof.oom);

        WriteUncertainty( txt, plus, minus );
        
        bars[prof.suffix][syst->ShortName()] = {minus, plus};
      }
      txt<< " \\\\" <<std::endl;
    }
  }

  txt<< "\\hline" <<std::endl;

  // Total uncertainty
  txt<< "Systematic Uncertainty" ;
  for(auto prof: profiles){
    SetFakeCalc(calc, ssth23, dmsq32, dCP);
    
    TGraph* totplot = (TGraph*)Profile(exptThis, calc,
                                   prof.fitvar,
                                   steps, prof.minx, prof.maxx, -1,
                                   prof.profvars,
                                   allSysts,
                   prof.profseeds);
    totplot -> Write((prof.suffix+"_totError").c_str());
    
    GetLimits(totplot, hi, low);

    double plus = sqrt( fabs( (hi-prof.bf)*(hi-prof.bf) - (statup[prof.suffix]*statup[prof.suffix]) ) ) * pow(10,prof.oom);
    double minus = sqrt( fabs( (prof.bf-low)*(prof.bf-low) - (statdn[prof.suffix]*statdn[prof.suffix]) ) ) * pow(10,prof.oom);

    WriteUncertainty( txt, plus, minus );

  }
  txt<< " \\\\" <<std::endl;

  txt<< "Statistical Uncertainty" ;
  for(auto prof: profiles){

    double plus = statup[prof.suffix] * pow(10,prof.oom);
    double minus = statdn[prof.suffix] * pow(10,prof.oom);
    WriteUncertainty( txt, plus, minus );
  }
  std::cout<< std::endl;

  fOut -> Close();
  delete fOut;

  for(auto prof: profiles){
    TCanvas *c;
    ErrorBarChart(bars[prof.suffix],{statdn[prof.suffix]*pow(10,prof.oom),statup[prof.suffix]*pow(10,prof.oom)},prof.root);
    gPad->SetLeftMargin(0.35);
    gPad->SetBottomMargin(0.15);
    simulation();
    gPad->Print( ( outDir + "/syst_table_" + prof.suffix + "_" + outStr + ".pdf" ).c_str() );
    gPad->Print( ( outDir + "/syst_table_" + prof.suffix + "_" + outStr + ".png" ).c_str() );
    gPad->Print( ( outDir + "/syst_table_" + prof.suffix + "_" + outStr + ".C"   ).c_str() );
  }

}