joint_fit_2019_bestfit.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/FC/FCSurface.h"
#include "CAFAna/Fit/FrequentistSurface.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 "Utilities/rootlogon.C"

#include "3FlavorAna/Ana2019/nue/joint_fit_2019_loader_tools.h"
//#include "./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>
#include <fstream>

using namespace ana;

void joint_fit_2019_bestfit(bool corrSysts = true,
                            TString options="joint_realData_both")

{
  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 fake2018 = options.Contains("fake2018");
  bool realData = options.Contains("realData"); 

  auto suffix = options;
  if(corrSysts) suffix+="_systs";
  else suffix+="_stats";

  TString outdir="./";//"/nova/ana/nu_e_ana/Ana2019/Results/BestFits/";

  TString outfilename   (outdir +  "bestfits_"    + suffix);

   //////////////////////////////////////////////////
   // 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();
    if(fake2018) SetFakeCalc(calc_fake, 0.58, 2.51e-3, 0.17);
    else if(!realData) {std::cerr << "need setting for data\n"; exit(1);}

    if(!numuOnly) {
      if(FHCOnly || both ) {
         preds.push_back({"nue_fhc", GetNuePrediction2019("combo", DefaultOscCalc(), corrSysts, "fhc", false), GetNueCosmics2019("fhc"), GetPOT()});
      } 
      if(RHCOnly || both ) {
         preds.push_back({"nue_fhc", GetNuePrediction2019("prop", DefaultOscCalc(), corrSysts, "rhc", false), GetNueCosmics2019("rhc"), GetPOT(false)});
      } 
    }

   int nnumu = 4;
   if(!nueOnly) {
      if(FHCOnly || both ) {
        auto numu_preds = GetNumuPredictions2019(nnumu, corrSysts, "fhc");
        auto numu_cosmics = GetNumuCosmics2019(nnumu, "fhc");
        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");
        auto numu_cosmics = GetNumuCosmics2019(nnumu, "rhc");
        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)});
        }
       }
    }


   if(realData){
      if(!numuOnly){
        if(FHCOnly || both ) data.push_back(GetNueData2019("fhc"));
        if(RHCOnly || both ) data.push_back(GetNueData2019("rhc"));
      }
      if(!nueOnly){
        if(FHCOnly || both ){
          auto numu_data = GetNumuData2019(nnumu, "fhc");
          data.insert(data.end(),numu_data.begin(), numu_data.end());
          data_numu_only.insert(data_numu_only.end(),numu_data.begin(), numu_data.end());}
        if(RHCOnly || both ){
          auto numu_data = GetNumuData2019(nnumu, "rhc");
          data.insert(data.end(),numu_data.begin(), numu_data.end());
          data_numu_only.insert(data_numu_only.end(),numu_data.begin(), numu_data.end());}
      }
    }

    for(int i = 0; i < int(preds.size()); ++i){
        double POT = preds[i].pot;
        auto thisdata = GetFakeData (preds[i].pred, calc_fake, POT, preds[i].cos.first);
        if(realData) thisdata = data[i];
        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));
        new TCanvas();
        DataMCComparison (*thisdata, preds[i].pred, DefaultOscCalc());
        preds[i].cos.first->SetMarkerStyle(34);
        preds[i].cos.first->SetMarkerColor(kCosmicBackgroundColor);
        preds[i].cos.first->Draw("hist p same");
        gPad->Print(outdir + "debug_predictions_" + suffix + "_" + std::to_string(i) + ".pdf");
    }

    cout<<"Make numu only experiment to get the seeds"<<endl;
    //make numu only experiment for seeds:

    for(int i = 0; i < (int) preds_numu_only.size(); ++i){
        double POT = preds_numu_only[i].pot;
        auto thisdata = GetFakeData (preds_numu_only[i].pred, calc_fake, POT, preds_numu_only[i].cos.first);
        if(realData) thisdata = data_numu_only[i];
        cout<<i<<" "<<preds_numu_only[i].name<<" POT "<<POT<<" tot MC "<<preds_numu_only[i].pred->Predict(calc_fake).Integral(POT)<<" cos "<<preds_numu_only[i].cos.first->Integral()<<" cos er "<<preds_numu_only[i].cos.second<<" analyze data "<<thisdata->Integral(POT)<<endl;
        expts_numu_only.push_back(new SingleSampleExperiment(preds_numu_only[i].pred, *thisdata, preds_numu_only[i].cos.first, preds_numu_only[i].cos.second));
    }
    ////////////////////////////////////////////////////////////
    // 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);

   std::cout << "Creating Multiexperiment of numu only SimpleExp with a total of "
              << expts_numu_only.size() << " experiments\n\n" << std::endl;
    auto exptThis_numu_only = new MultiExperiment(expts_numu_only);

    ////////////////////////////////////////////////////////////
    // Systematics
    ////////////////////////////////////////////////////////////
    std::cout << "Systematics for the fit:\n";
    auto systs = GetJointFitSystematicList(corrSysts, nueOnly, numuOnly, FHCOnly||both, RHCOnly||both, 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);       
      } 
    }
    if (nueOnly && corrSysts){
       if (FHCOnly) exptThis->SetSystCorrelations(0, GetCorrelations(true, true));
       if (RHCOnly) exptThis->SetSystCorrelations(0, GetCorrelations(true, false));
       if (both) {
          exptThis->SetSystCorrelations(0, GetCorrelations(true, true));
          exptThis->SetSystCorrelations(1, GetCorrelations(true, false));
       }
    }

    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 ////////////////////////////
  //////////////////////////////////////////////////////////////////////

   cout<<"------------------- Start prestage seeds --------------------------"<<endl;

    double minchi_numu = 1E20;
    double pre_seed_th23;
    double pre_seed_dmsq;
    ResetOscCalcToDefault(calc);
    auxShifts.ResetToNominal();

   double maxmix = 0.514;  // from the numu results
   double numu_pre_seedLONH, numu_pre_seedUONH, numu_pre_seedLOIH, numu_pre_seedUOIH, dmsq_numu_pre_seedNH, dmsq_numu_pre_seedIH;

   for(int hie:{1}){
      std::cout << "\n\nFinding test best fit " << (hie>0? "NH " : "IH ") << "\n\n";
      Fitter fitnumu_only(exptThis_numu_only, {&kFitSinSqTheta23, &kFitDmSq32}, {});

      auto thisminchi = fitnumu_only.Fit(calc, auxShifts ,
                                        {{&kFitDmSq32,{hie*2.5e-3}},
                                        {&kFitSinSqTheta23, {0.4} }}, {}, //seedShifts,
                                        Fitter::kVerbose);

     if(thisminchi < minchi_numu) minchi_numu = thisminchi;

     pre_seed_th23 = kFitSinSqTheta23.GetValue(calc);
     pre_seed_dmsq = kFitDmSq32.GetValue(calc);
   }

  numu_pre_seedLONH = ((pre_seed_th23>maxmix)?(2*maxmix-pre_seed_th23):pre_seed_th23);
  numu_pre_seedUONH = ((pre_seed_th23>maxmix)?pre_seed_th23:(2*maxmix-pre_seed_th23));

  ResetOscCalcToDefault(calc);
  auxShifts.ResetToNominal();

  cout<<"------------------- End prestage seeds --------------------------"<<endl;


   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";

   //for NH:
   th23seeds.push_back( { {0.499, numu_pre_seedLONH}, kFitSinSqTheta23LowerOctant, "LO"});
   th23seeds.push_back( { {0.501, numu_pre_seedUONH}, kFitSinSqTheta23UpperOctant, "UO"});

   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 //////////////////
   /////////////////////////////////////////////////////////////////////////

    // Find the best fit points
    TFile * outfile = new TFile(outfilename+".root","recreate");
    outfile->cd();

    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, 
                                                &kFitDmSq32, 
                                                &kFitSinSq2Theta13};
        
        Fitter fit23(exptThis, fitvars, systs);
        cout<<" pre dmsq seed is "<<pre_seed_dmsq<<endl;
        auto thisminchi = fit23.Fit(calc, auxShifts,
                                    {
                                      { &kFitDmSq32,   {hie*pre_seed_dmsq} },
                                      { thseed.var,          thseed.seeds },
                                      { &kFitDeltaInPiUnits, delta_seeds },
                                      { &kFitSinSq2Theta13, {0.082}}
                                    }, seedShifts, 
                                    Fitter::kVerbose);
        minchi23= min(minchi23, thisminchi);

        if(corrSysts){
          // Plot the systematic pulls and label with BFV
          auto shifts = PlotSystShifts(auxShifts);
          TString str = "Best fit " ;
          for (auto &v:fitvars){
            str += TString::Format(" %s=%.3f ",v->LatexName().c_str(),v->GetValue(calc));
          }
          str+= TString::Format(" LL=%.6f", thisminchi);
          shifts->SetTitle(str);
          gPad->Update();
          TLine *l=new TLine(gPad->GetUxmin(),0,gPad->GetUxmax(),0);
          l->Draw("same"); 
          gPad->Print(outdir + "debug_slice_shifts_bestfit_" + suffix +
                      (hie>0? "NH": "IH") + thseed.label + ".pdf");
        }
        TString dirname = (hie > 0 ? "NH_" : "IH_") + thseed.label;
        SaveTo(* (osc::IOscCalc*) calc, outfile, dirname+"_calc");
        auxShifts.SaveTo(outfile, dirname+"_systs");
        ResetOscCalcToDefault(calc);
        auxShifts.ResetToNominal();
      }
    }
    std::cout << "\nFound overall minchi " << minchi23 << "\n\n";
   
    cout<<"Check with oldstyle seeds "<<endl;
    double minchi23test = 1E20;
      for(int hie:{-1,1}){
        std::cout << "\n\nFinding test best fit "
                  << (hie>0? "NH " : "IH ")
                  << "\n\n";
        std::vector <const IFitVar*> fitvars = {&kFitDeltaInPiUnits,
                                                &kFitSinSqTheta23,
                                                &kFitDmSq32,
                                                &kFitSinSq2Theta13};

        Fitter fit23(exptThis, fitvars, systs);
        auto thisminchi = fit23.Fit(calc,auxShifts ,{
            {&kFitDmSq32,{hie*fabs(kFitDmSq32.GetValue(calc))}},
            {&kFitSinSqTheta23, th23_old_seeds},
            {&kFitDeltaInPiUnits, delta_seeds}
          },  seedShifts
          );
        minchi23test= min(minchi23test, thisminchi);
        if(corrSysts){
        auto shifts = PlotSystShifts(auxShifts);
          TString str = "Best fit " ;
          for (auto &v:fitvars){
            str += TString::Format(" %s=%.3f ",v->LatexName().c_str(),v->GetValue(calc));
          }
          str+= TString::Format(" LL=%.6f", thisminchi);
          shifts->SetTitle(str);
          gPad->Print(outdir  + "debug_slice_shifts_bestfit_" + suffix +
                      (hie>0? "NH": "IH") + "test_chisq_with_no_octant_differentiation.pdf");
        }
        
        ResetOscCalcToDefault(calc);
        auxShifts.ResetToNominal();
      }//end hie
    std::cout << "\nFound overall test minchi " << minchi23test << std::endl;


    cout<<"<<<<<<<<<<<<<<<<<<<<<<<<<< Test with only systematics free"<<endl;
    minchi23test = 1E20;
    ResetOscCalcToDefault(calc);
    auxShifts.ResetToNominal();

    if(both){
      calc->SetdCP(1.99969*M_PI);
      calc->SetTh23(asin(sqrt(0.564956)));
      calc->SetDmsq32(0.00248304);
      calc->SetTh13(asin(sqrt(0.0818912))/2);
    }
    if(FHCOnly){
      calc->SetdCP(1.21*M_PI);
      calc->SetTh23(asin(sqrt(0.513)));
      calc->SetDmsq32(0.00246);
      calc->SetTh13(asin(sqrt(0.082))/2);
    }
    if(RHCOnly){
      calc->SetdCP(1.93*M_PI);
      calc->SetTh23(asin(sqrt(0.598)));
      calc->SetDmsq32(0.002515);
      calc->SetTh13(asin(sqrt(0.082))/2);
    }

    Fitter fit23(exptThis, {}, systs);
    auto thisminchi = fit23.Fit(calc,auxShifts ,{}, seedShifts,
          Fitter::kVerbose);
   ofstream txtfile(outdir  + "debug_slice_shifts_bestfit_" + suffix +
                      + "test_chisq_with_osc_par_fixed.txt");
   for (auto &v:systs){
       txtfile<<v->ShortName().c_str()<<"="<<auxShifts.GetShift(v)<<endl;
    }
    txtfile.close();
    minchi23test= min(minchi23test, thisminchi);
    auto shifts = PlotSystShifts(auxShifts);
          TString str = "Best fit " ;
          str+= TString::Format(" LL=%.6f", thisminchi);
          shifts->SetTitle(str);
          gPad->Print(outdir +  "debug_slice_shifts_bestfit_" + suffix +
                      + "test_chisq_with_osc_par_fixed.pdf");
   ResetOscCalcToDefault(calc);
   auxShifts.ResetToNominal();
   cout<<"<<<<<<<<<<<< This is not true BF : "<< minchi23test<<endl;

}