goodness_of_fit.C

Go to the documentation of this file.

#include "CAFAna/Analysis/Calcs.h"
#include "CAFAna/Analysis/Exposures.h"
#include "CAFAna/Fit/Fit.h"
#include "CAFAna/Analysis/Plots.h"
#include "3FlavorAna/Plotting/NuePlotStyle.h"
#include "CAFAna/Analysis/Style.h"
#include "CAFAna/Experiment/Dmsq32Constraint.h"
#include "CAFAna/Experiment/MultiExperiment.h"
#include "CAFAna/Experiment/ReactorExperiment.h"
#include "CAFAna/Experiment/SingleSampleExperiment.h"
#include "CAFAna/Fit/FrequentistSurface.h"
#include "CAFAna/FC/FCSurface.h"
#include "3FlavorAna/Prediction/PredictionSystJoint2018.h"
#include "CAFAna/Prediction/PredictionCombinePeriods.h"
#include "3FlavorAna/Systs/NumuSysts.h"
#include "CAFAna/Vars/FitVars.h"
#include "OscLib/IOscCalc.h"
#include "3FlavorAna/Ana2019/nue/joint_fit_2019_loader_tools.h"

#include "TCanvas.h"
#include "TBox.h"
#include "TColor.h"
#include "TGraph.h"
#include "TVectorD.h"
#include "TF1.h"
#include "TLegend.h"
#include "TText.h"
#include "TLatex.h"
#include "TPad.h"
#include "TLine.h"
#include "TMarker.h"
#include "TStyle.h"
#include "TSystem.h"
#include "TGaxis.h"

#include <algorithm>
#include <vector>
#include <string>

using namespace ana;

void goodness_of_fit(int num = 1, bool throwexp = true, bool corrSysts = false,
                     TString options="joint_realData_both", TString nametag = "test_fake2019")
{
  bool nueOnly = options.Contains("nueOnly");
  bool numuOnly = options.Contains("numuOnly");
  bool joint = options.Contains("joint");
  assert (nueOnly || numuOnly || joint);

  bool FHCOnly = options.Contains("FHCOnly");
  bool RHCOnly = options.Contains("RHCOnly");
  bool both = options.Contains("both");
  assert (FHCOnly || RHCOnly || both);

  bool fake2019 = nametag.Contains("fake2019");

  TString outdir = "./";

  if(throwexp){ 
    //////////////////////////////////////////////////
    // Load Nue and Numu experiments
    //////////////////////////////////////////////////
    //need numu only for prestage seeds

    struct predictions {
           const string name;
           const IPrediction * pred;
           std::pair <TH1D*, double> cos;
           double pot;
    };

    std::vector <predictions> preds;
    std::vector <predictions> preds_numu_only;
    std::vector <Spectrum * > data;
    std::vector <Spectrum * > data_numu_only;
    std::vector <const IExperiment * > expts;
    std::vector <const IExperiment * > expts_numu_only;

    auto calc_fake = DefaultOscCalc();
    double th23, dmsq, dcp;
    if(fake2019) {th23 = 0.565; dmsq = 2.48e-3; dcp = 1.99;}

    if(!numuOnly) {
      if(FHCOnly || both ) {
         preds.push_back({"nue_fhc", GetNuePrediction2019("combo", DefaultOscCalc(), corrSysts, "fhc", false, true), GetNueCosmics2019("fhc", true), GetPOT()});
      } 
      if(RHCOnly || both ) {
         preds.push_back({"nue_fhc", GetNuePrediction2019("prop", DefaultOscCalc(), corrSysts, "rhc", false, true), GetNueCosmics2019("rhc", true), GetPOT(false)});
      } 
    }
   int nnumu = 4;
   if(!nueOnly) {
      if(FHCOnly || both ) {
        auto numu_preds = GetNumuPredictions2019(nnumu, corrSysts, "fhc", true);
        auto numu_cosmics = GetNumuCosmics2019(nnumu, "fhc", true);
        for (int i = 0; i< nnumu; i++ ){
            preds.push_back({"numu_fhc_"+std::to_string(i+1), numu_preds[i], numu_cosmics[i], GetPOT()});
            preds_numu_only.push_back({"numu_fhc_"+std::to_string(i+1), numu_preds[i], numu_cosmics[i], GetPOT()});
        }
      }
      if(RHCOnly || both ) {
        auto numu_preds = GetNumuPredictions2019(nnumu, corrSysts, "rhc", true);
        auto numu_cosmics = GetNumuCosmics2019(nnumu, "rhc", true);
        for (int i = 0; i< nnumu; i++ ){
            preds.push_back({"numu_rhc_"+std::to_string(i+1), numu_preds[i], numu_cosmics[i], GetPOT(false)});
            preds_numu_only.push_back({"numu_rhc_"+std::to_string(i+1), numu_preds[i], numu_cosmics[i], GetPOT(false)});
        }
       }
    }
     
     TFile* histNH_file = new TFile("/cvmfs/nova.osgstorage.org/analysis/nue/Ana2019/goodness/hist_contours_2019_"+options+"_onlyNH_systs.root","read");
     TFile* histIH_file = new TFile("/cvmfs/nova.osgstorage.org/analysis/nue/Ana2019/goodness/hist_contours_2019_"+options+"_onlyIH_systs.root","read");

     TH2F* histNH_pdf = new TH2F("histNH_pdf","histNH_pdf", 30, 0, 2*M_PI, 30, 0.3, 0.7);
     TH2F* histIH_pdf = new TH2F("histIH_pdf","histIH_pdf", 30, 0, 2*M_PI, 30, 0.3, 0.7);

     auto histNH = (TH2F*)histNH_file->Get("delta_th23_NH");
     auto histIH = (TH2F*)histIH_file->Get("delta_th23_IH");

     for(int i = 1; i<=30; i++){
         for(int j = 1; j<=30; j++){
             double binNH = histNH->GetBinContent(i,j);
             histNH_pdf->SetBinContent(i, j, exp(-binNH/2.0));
             double binIH = histIH->GetBinContent(i,j);
             histIH_pdf->SetBinContent(i, j, exp(-binIH/2.0));
         }
     }

     new TCanvas();
     histNH_pdf->Draw("colz");
     gPad->Print(outdir+"pdf_transform_hist2d_delta_th23_nh_"+options+"_"+nametag+".pdf");

     new TCanvas();
     histIH_pdf->Draw("colz");
     gPad->Print(outdir+"pdf_transform_hist2d_delta_th23_ih_"+options+"_"+nametag+".pdf");

     cout<<"integrals: "<<histNH_pdf->Integral()<<" and ih is "<<histIH_pdf->Integral()<<endl;
     double sep = histNH_pdf->Integral()/(histNH_pdf->Integral()+histIH_pdf->Integral());
    
     cout<<"probablilty to get the NH is "<<histNH_pdf->Integral()/(histNH_pdf->Integral()+histIH_pdf->Integral())<<" and to get the ih is "<<histIH_pdf->Integral()/(histNH_pdf->Integral()+histIH_pdf->Integral())<<endl;

     TFile* dmsq_file = new TFile("/cvmfs/nova.osgstorage.org/analysis/nue/Ana2019/goodness/hist_slices_2019_"+options+"_systs_dmsq_noOct.root","read");

     TH1F* histdmsqNH_pdf = new TH1F("histdmsqNH_pdf", "histdmsqNH_pdf", 60, 2e-3, 3e-3);
     TH1F* histdmsqIH_pdf = new TH1F("histdmsqIH_pdf", "histdmsqIH_pdf", 60, -3e-3, -2e-3);

     auto histdmsqNH = (TH1F*)dmsq_file->Get("slice_dmsq_NH");
     auto histdmsqIH = (TH1F*)dmsq_file->Get("slice_dmsq_IH");

     for(int i = 1; i<=60; i++){
        double binNH = histdmsqNH->GetBinContent(i);
        histdmsqNH_pdf->SetBinContent(i, exp(-binNH/2.0));
        double binIH = histdmsqIH->GetBinContent(i);
        histdmsqIH_pdf->SetBinContent(i, exp(-binIH/2.0));
     }

     new TCanvas();
     histdmsqNH_pdf->Draw("hist");
     gPad->Print(outdir+"pdf_transform_dmsq32_nh_"+options+"_"+nametag+".pdf");

     new TCanvas();
     histdmsqIH_pdf->Draw("hist");
     gPad->Print(outdir+"pdf_transform_dmsq32_ih_"+options+"_"+nametag+".pdf");

     auto filename = outdir+ "throw_exp_"+options + (corrSysts?"_with_syst_":"_just_stat_")+std::to_string(num)+"_exp_"+nametag+".root";

     TFile* file = new TFile(filename,"recreate");

  for (int j = 0; j<num; j++){

    file->cd();

    TDirectory* d = file->mkdir(("ex_"+std::to_string(j)).c_str());
    file->cd(("ex_"+std::to_string(j)).c_str());

    TRandom3 rnd(0);

    double hie = rnd.Uniform(1);
    cout<<"hie is "<<hie<<" hie separator is "<<sep<<endl;
    bool isNH;
    if(hie<sep) isNH = true;
    else isNH = false;
    double thisdcp, thisth23, thisdmsq, thisth13;
    if (isNH) {
        histNH_pdf->GetRandom2(thisdcp, thisth23);
        thisdmsq= histdmsqNH_pdf->GetRandom();
     }
     else{
        histIH_pdf->GetRandom2(thisdcp, thisth23);
        thisdmsq = histdmsqIH_pdf->GetRandom();
     }
     thisth13 = rnd.Gaus(0.082, 0.004);

     calc_fake->SetdCP(thisdcp*M_PI);
     calc_fake->SetTh23(asin(sqrt(thisth23)));
     calc_fake->SetDmsq32(thisdmsq);
     calc_fake->SetTh13(asin(sqrt(thisth13))/2);

     TVectorD v(4);
     v[0] = thisdmsq;
     v[1] = thisth23;
     v[2] = thisdcp;
     v[3] = thisth13;
     v.Write((std::to_string(j)+"_osc_pars").c_str());
     cout<<"osc pars in exp "<<j<<" are "<<thisdmsq<<" "<<thisth23<<" "<<thisdcp<<endl;

     struct systshifts{
        double shift;
        TString systname;
     };
     std::vector <systshifts> shifts;

     auto systs = GetJointFitSystematicList(corrSysts, false, false, true, true, true);

    for (auto &s:systs){
         auto sh = rnd.Gaus(0, 1);
         shifts.push_back({sh, s->ShortName().c_str()});
         cout<<s->ShortName().c_str()<<"("<<sh<<"), ";
         TVectorD vs(1);
         vs[0] = sh;
         vs.Write((std::to_string(j)+"_"+s->ShortName()).c_str());
    }
    cout<<"\n\n";

    for(int i = 0; i < (int) preds.size(); ++i){

      double POT = preds[i].pot;

      SystShifts sshifts;

      if(corrSysts){
       if(both){
         if(i==0) systs = getAllAna2018Systs(kNueAna2018, true, kFHC);
         if(i==1) systs = getAllAna2018Systs(kNueAna2018, true, kRHC);
         if(i>= 2) systs = getAllAna2018Systs(kNumuAna2018); // doesn't play any role
        }
        if(FHCOnly || RHCOnly){
           if(i==0) systs = getAllAna2018Systs(kNueAna2018, true, (FHCOnly?kFHC:kRHC));
           if(i>= 1) systs = getAllAna2018Systs(kNumuAna2018);
        }
         for (auto &s:systs){         
             for(auto &sh:shifts){ if(s->ShortName() == sh.systname) {sshifts.SetShift(s, sh.shift); cout<<s->ShortName()<<" "<<sh.shift<<", ";} }
         }
      }
      cout<<"\n\n";

      auto tempmock = preds[i].pred->PredictSyst(calc_fake, sshifts).MockData(POT);
      auto tempfake = preds[i].pred->PredictSyst(calc_fake, sshifts).FakeData(POT);

      tempmock += Spectrum(preds[i].cos.first, POT,0).MockData(POT);
      tempfake += Spectrum(preds[i].cos.first, POT,0);

      auto thisdatamock = new Spectrum(tempmock);
      thisdatamock->SaveTo(d, (std::to_string(i)+"_check_mockdata").c_str());
      auto thisdatafake = new Spectrum(tempfake);
      thisdatafake->SaveTo(d, (std::to_string(i)+"_check_fakedata").c_str());
      auto thisdata =thisdatamock;
      thisdata->SaveTo(d, (std::to_string(i)+"_mockdata").c_str());
      cout<<i<<" "<<preds[i].name<<" POT "<<POT<<" tot MC "<<preds[i].pred->Predict(calc_fake).Integral(POT)<<" cos "<<preds[i].cos.first->Integral()<<" cos er "<<preds[i].cos.second<<" analyze data "<<thisdata->Integral(POT)<<endl;
      expts.push_back(new SingleSampleExperiment(preds[i].pred, *thisdata, preds[i].cos.first, preds[i].cos.second));
      //auto single_exp = new SingleSampleExperiment(preds[i].pred, *thisdata, preds[i].cos.first, preds[i].cos.second); // this is heavy
      //single_exp->SaveTo(d, ("single_exp_"+std::to_string(i)).c_str()); // this is heavy
    }

    ////////////////////////////////////////////////////////////
    // Add constraints, make experiments
    ////////////////////////////////////////////////////////////
   
    std::cout << "\nCreating multiexperiment\n" << std::endl;

    expts.push_back(WorldReactorConstraint2017()); std::cout << "Adding WorldReactorConstraint2017\n";
    if(nueOnly) {
      std::cout << "Adding Dmsq32ConstraintPDG2017\n";
      expts.push_back(&kDmsq32ConstraintPDG2017);
    }
    std::cout << "Creating Multiexperiment with a total of "
              << expts.size() << " experiments\n\n" << std::endl;
    auto exptThis = new MultiExperiment(expts);
    //file->Close();

    ////////////////////////////////////////////////////////////
    // Systematics
    ////////////////////////////////////////////////////////////
    std::cout << "Systematics for the fit:\n";
    systs = GetJointFitSystematicList(corrSysts, numuOnly, nueOnly, true, true, true);

    std::cout << "\n\nSystematic correlations...\n";
    if(!nueOnly && ! numuOnly && corrSysts){
      if(FHCOnly){
        exptThis->SetSystCorrelations(0, GetCorrelations(true, true));
        auto notfornumu = GetCorrelations(false, true);
        for(int i =0; i < nnumu; ++i) exptThis->SetSystCorrelations(i+1, notfornumu);   
      }
      if(RHCOnly){
        exptThis->SetSystCorrelations(0, GetCorrelations(true, false));
        auto notfornumu = GetCorrelations(false, true);
        for(int i =0; i < nnumu; ++i) exptThis->SetSystCorrelations(i+1, notfornumu);   
      }
      if(both){
        exptThis->SetSystCorrelations(0, GetCorrelations(true, true));
        exptThis->SetSystCorrelations(1, GetCorrelations(true, false));
        auto notfornumu = GetCorrelations(false, true);
        for(int i =0; i < 8; ++i) exptThis->SetSystCorrelations(i+2, notfornumu);       
      } 
    }

    std::cout << "Systematics for the numu only fit:\n";
    auto systs_numu_only = GetJointFitSystematicList(corrSysts, false, true, true, true, true);

    ////////////////////////////////////////////////////////////
    // Fit
    ////////////////////////////////////////////////////////////
    std::cout << "Starting the fit" << std::endl;

    osc::IOscCalcAdjustable* calc = DefaultOscCalc();

    SystShifts auxShifts = SystShifts::Nominal();

    // In case some systs need different seeds
    std::vector <SystShifts> seedShifts = {};

//////////////////////////////////////////////////////////////////////
///////////////////////// Seed controller ////////////////////////////
//////////////////////////////////////////////////////////////////////


   struct th23helper{
     std::vector<double> seeds;
     const IFitVar * var;
     TString label;
   };

   std::vector <th23helper> th23seeds;

   th23helper temp;
   temp.seeds = {0.3};
   temp.var = kFitSinSqTheta23LowerOctant;
   temp.label = "LO";
   bool octantSlice =true;
   //for NH:
  if(octantSlice) {
    th23seeds.push_back( { {0.499, 0.45}, kFitSinSqTheta23LowerOctant, "LO"});
    th23seeds.push_back( { {0.501, 0.55}, kFitSinSqTheta23UpperOctant, "UO"});
   }
   else th23seeds.push_back({ {0.45, 0.5, 0.55}, &kFitSinSqTheta23, ""});

   std::vector<double> delta_seeds = {0, 0.5, 1., 1.5};
   std::vector<double> th23_old_seeds = {0.45, 0.5, 0.55};

/////////////////////////////////////////////////////////////////////////
////////////////////////// Different best fit searches //////////////////
/////////////////////////////////////////////////////////////////////////

    file->cd(("ex_"+std::to_string(j)).c_str());
    // Find the best fit points
    TVectorD lchi(5);
    double minchi23 = 1E20;
    for(int hie:{-1, 1}){
      for (auto & thseed:th23seeds){

        std::cout << "\n\nFinding best fit " << (hie > 0 ? "NH " : "IH ")
                  << thseed.label << ", "
                  << "ssth23 seeds ";
        for (auto const & s:thseed.seeds) std::cout << s << ", ";
        std::cout << std::endl;
        std::vector <const IFitVar*> fitvars = {&kFitDeltaInPiUnits,
                                                thseed.var, 
                                                &kFitDmSq32Scaled, 
                                                &kFitSinSq2Theta13};
        Fitter fit23(exptThis, fitvars, systs);
        auto thisminchi = fit23.Fit(calc, auxShifts,
                                    {
                                      { &kFitDmSq32Scaled,   {hie*2.5} },
                                      { thseed.var,          thseed.seeds },
                                      { &kFitDeltaInPiUnits, delta_seeds },
                                      { &kFitSinSq2Theta13, {0.082} }
                                    }, seedShifts, 
                                    Fitter::kQuiet);
        minchi23= min(minchi23, thisminchi);
        cout<<minchi23<<endl;

        lchi[0] = thisminchi;
        lchi[1] = kFitDeltaInPiUnits.GetValue(calc);
        lchi[2] = thseed.var->GetValue(calc);
        lchi[3] = kFitDmSq32Scaled.GetValue(calc);
        lchi[4] = kFitSinSq2Theta13.GetValue(calc);
        lchi.Write((hie>0?"NH_":"IH_")+thseed.label+"_osc_pars");

        ResetOscCalcToDefault(calc);
        auxShifts.ResetToNominal();
      }
    }
    TVectorD chi(1);
    chi[0] = minchi23;
    chi.Write("best_chi");
    expts.clear();
    
    std::cout << "\nFound overall minchi " << minchi23 << "\n\n";
 }
 file->Close();
}

if(!throwexp){

  int min = 110; int max = 250; int bins = 140;
  if(RHCOnly || FHCOnly) {min = 40; max = 140; bins = 100;}
  TH1F* chisq = new TH1F("chisq", "best chisq in each mock exp", bins, min, max);  
     double LLreal;
   if(both) LLreal = 157.146;
   if(RHCOnly) LLreal = 82.1438;
   if(FHCOnly) LLreal = 73.2594;
   int counter = 0;

   outdir = "/pnfs/nova/scratch/users/lkolupae/2019/goodness/";
   
  for (int j=0; j<num; j++){
      auto filename = outdir+ "throw_exp_"+options + (corrSysts?"_with_syst_":"_just_stat_")+"5_exp_thr"+std::to_string(j)+"_fake2019.root";
      TFile* file = new TFile(filename,"read");
      for (int k =0; k<5; k++){
        auto min =* (TVectorD*)file->Get(("ex_"+std::to_string(k)+"/best_chi").c_str());
        double minchi = min[0];
        if(minchi>LLreal) counter++;
        chisq->Fill(minchi);
        cout<<minchi<<endl; 
      }
   }  

   TLine* line = new TLine(LLreal, 0, LLreal, 20);
   line->SetLineColor(kRed);
    cout<<"chisq intgr "<<chisq->Integral()<<endl;
    cout<<"events with N> data chisq "<<counter;
   cout<<" out of "<<chisq->Integral()<<endl;
   new TCanvas;
   chisq->Draw("hist");
   line->Draw();
   gPad->Print("exp_"+options+".pdf");

}


}