run_joint_fit_2020_contours.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 "Utilities/rootlogon.C"
#include "OscLib/IOscCalc.h"

#include "3FlavorAna/Ana2020/joint_fit_2020_loader_tools.h"

#include "TCanvas.h"
#include "TMarker.h"
#include "TBox.h"
#include "TLatex.h"
#include "TColor.h"
#include "TGraph.h"
#include "TVectorD.h"
#include "TF1.h"
#include "TLegend.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 run_joint_fit_2020_contours(bool corrSysts = false,
                                 bool runOnGrid = false,
                                 TString options="joint_fake2019_both_onlyIH",
                                 bool dmsqSurf = false,
                                 bool th13Surf = false)
{
  bool debug = 1;
  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 = options.Contains("fake2019");
  bool realData = options.Contains("realData");

  bool onlyNH = options.Contains("onlyNH");
  bool onlyIH = options.Contains("onlyIH");

  TString tag;
  if(FHCOnly) tag="FHCOnly/";
  if(RHCOnly) tag="RHCOnly/";
  if(both) tag = "RHC_and_FHC/";

  auto suffix = options;
  if(dmsqSurf) suffix+= "_dmsq";
  if(th13Surf) suffix+= "_th13";
  if(corrSysts) suffix+="_systs";
  else suffix+="_stats";

  TString outdir = "/nova/ana/nu_e_ana/Ana2020/prefits/";//"/nova/ana/nu_e_ana/Ana2020/Results/"+tag+"contours/FCInputs/";
  if (runOnGrid) outdir = "./";
  TString outFCfilename (outdir +"contours_FCInput_2020_" + suffix);

  // outdir="/nova/ana/nu_e_ana/Ana2020/Results/"+tag+"contours/";
  /*if(th13Surf) outdir += "th13_delta/";
  else if (dmsqSurf) outdir += "th23_dmsq/";
  else outdir += "delta_th23/";
  outdir += corrSysts?"syst/":"stat/";
  if (runOnGrid) outdir = "./";*/
  TString outfilename   (outdir +  "hist_contours_2020_"    + suffix);

  //////////////////////////////////////////////////
  // Load Nue and Numu experiments
  //////////////////////////////////////////////////
  //need numu only for prestage seeds

  //std::string decomp = "noPt"; // use Pt extrap all the times
  std::string decomp = ""; // use Pt extrap all the times
  std::vector <ana::predictions> preds = LoadPredictions (corrSysts, runOnGrid, decomp, nueOnly || joint, numuOnly || joint, FHCOnly || both, RHCOnly || both);
  std::vector <ana::predictions> preds_numu_only = LoadPredictions (corrSysts, runOnGrid, decomp, false, numuOnly || joint, FHCOnly || both, RHCOnly || both);
  std::vector <Spectrum * > data = LoadRealData (runOnGrid, nueOnly || joint, numuOnly || joint, FHCOnly || both, RHCOnly || both);
  std::vector <Spectrum * > data_numu_only = LoadRealData (runOnGrid, false, numuOnly || joint, FHCOnly || both, RHCOnly || both);
  std::vector <const IExperiment * > expts;
  std::vector <const IExperiment * > expts_numu_only;

  int nnumu = 4;

  auto calc_fake = DefaultOscCalc();
  if(fake2019) SetFakeCalc(calc_fake, 0.565, 2.48e-3, 0);
  else if(!realData) {std::cerr << "need setting for data\n"; exit(1);}

  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, preds[i].livetime);
      if(realData) thisdata = data[i];
      auto hcos = preds[i].cos.first->ToTH1(preds[i].livetime, kBlack, kSolid, kLivetime);
      cout<<i<<" "<<preds[i].name<<" POT "<<POT<<" tot MC "<<preds[i].pred->Predict(calc_fake).Integral(POT)<<
               " cos "<<hcos->Integral()<<" cos er "<<preds[i].cos.second<<" analyze data "<<thisdata->Integral(thisdata->POT())<<endl;
      expts.push_back(new SingleSampleExperiment(preds[i].pred, *thisdata, *preds[i].cos.first, preds[i].cos.second, true));
      if(debug){
         new TCanvas();
         DataMCComparison (*thisdata, preds[i].pred, DefaultOscCalc());
         hcos->SetMarkerStyle(34);
         hcos->SetMarkerColor(kCosmicBackgroundColor);
         hcos->Draw("hist p same");
         gPad->Print(outdir + (corrSysts?"syst/":"stat/")+ "debug_predictions_" + suffix + "_" + std::to_string(i) + ".pdf");
      }//end debug
  }

  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, preds_numu_only[i].livetime);
      if(realData) thisdata = data_numu_only[i];
      auto hcos = preds_numu_only[i].cos.first->ToTH1(preds_numu_only[i].livetime, kBlack, kSolid, kLivetime);
      cout<<i<<" "<<preds_numu_only[i].name<<" POT "<<POT<<" tot MC "<<preds_numu_only[i].pred->Predict(calc_fake).Integral(POT)<<
               " cos "<<hcos->Integral()<<" cos er "<<preds_numu_only[i].cos.second<<" analyze data "<<thisdata->Integral(thisdata->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, true));
  }

  //////////////////////////////////////////////////
  // Add constraints, make experiments
  //////////////////////////////////////////////////
  std::cout << "\nCreating multiexperiment\n" << std::endl;
  if(!th13Surf){
      std::cout << "Adding WorldReactorConstraint2019\n";
      expts.push_back(WorldReactorConstraint2019());
  }
  if(nueOnly) {
     std::cout << "Adding Dmsq32ConstraintPDG2019\n";
     expts.push_back(&kDmsq32ConstraintPDG2019);
  }
  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);

  std::cout<<"Clean verctors"<<std::endl;
  preds.clear();
  preds_numu_only.clear();

  ////////////////////////////////////////////////////////////
  // Systematics
  ////////////////////////////////////////////////////////////

  bool ptExtrap = true;
  if (decomp == "noPt")  ptExtrap = false; 
  auto systs = GetJointFitSystematicList(corrSysts, nueOnly, numuOnly, FHCOnly||both, RHCOnly||both, true, ptExtrap);

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

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

  if (corrSysts && numuOnly && both){
      auto notfornumufhc = GetCorrelations(false, true, ptExtrap);
      for(int i =0; i < 4; ++i) {
          exptThis->SetSystCorrelations(i, notfornumufhc);
          exptThis_numu_only->SetSystCorrelations(i, notfornumufhc);
      }
      auto notfornumurhc = GetCorrelations(false, false, ptExtrap);
      for(int i =0; i < 4; ++i) {
          exptThis->SetSystCorrelations(i+4, notfornumurhc);
          exptThis_numu_only->SetSystCorrelations(i+4, notfornumurhc);
      }
  }

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

  osc::IOscCalcAdjustable* calc = DefaultOscCalc();

  SystShifts auxShifts = SystShifts::Nominal();

  std::vector <SystShifts> seedShifts = {};
  if (corrSysts) {
    for (double systshift:{-2, +2}){
      SystShifts tempShifts (&kAnaLightlevelNDSyst, systshift);  // attempt to fix flip between positive and negative pulls
      seedShifts.push_back(tempShifts);
    }
  }
  //////////////////////////////////////////////////////////////////////
  /////////////////////////// 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";
      MinuitFitter fitnumu_only(exptThis_numu_only, {&kFitSinSqTheta23, &kFitDmSq32Scaled}, {});

      auto thisminchi = fitnumu_only.Fit(calc, auxShifts ,
                                        SeedList({{&kFitDmSq32Scaled,{hie*2.5}},
                                        {&kFitSinSqTheta23, {0.4} }}), {},// seedShifts
                                         IFitter::kVerbose)->EvalMetricVal();

      if(thisminchi < minchi_numu) minchi_numu = thisminchi;
  }
  pre_seed_th23 = kFitSinSqTheta23.GetValue(calc);
  pre_seed_dmsq = kFitDmSq32Scaled.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> th23_all_seeds_NH = {numu_pre_seedLONH, 0.5, numu_pre_seedUONH};

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

  double minchi23 = 1E20;
  for(int hie:{-1,1}){
      for (auto & thseed: th23seeds){
           std::cout << "\n\nFinding best fit " << thseed.label
                  << (hie>0? " NH " : " IH ") <<"\n"<<endl;
           for (auto const & s:thseed.seeds) std::cout << s << ", ";
           std::cout << std::endl;
           cout<<"pre seed from numu is "<<pre_seed_dmsq<<endl;
           std::vector <const IFitVar*> fitvars = {&kFitDeltaInPiUnits,
                                                   thseed.var,
                                                   &kFitDmSq32Scaled,
                                                   &kFitSinSq2Theta13};

          //if (th13Surf) fitvars = {&kFitDeltaInPiUnits,thseed.var,&kFitDmSq32Scaled};

           MinuitFitter fit23(exptThis, fitvars, systs);
           auto thisminchi = fit23.Fit(calc, auxShifts, SeedList({
                                       {&kFitDmSq32Scaled,{hie*pre_seed_dmsq}},
                                       {thseed.var, thseed.seeds},
                                       {&kFitDeltaInPiUnits, delta_seeds}
                                       }),
                                       {}, IFitter::kVerbose)->EvalMetricVal();
           minchi23= min(minchi23, thisminchi);
           ResetOscCalcToDefault(calc);
           auxShifts.ResetToNominal();
      }//end th23
  }//end hie
  std::cout << "\nFound overall minchi " << minchi23 << std::endl;

  cout<<"Check with oldstyle seeds when we didn't split into octant"<<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,
                                              &kFitDmSq32Scaled,
                                              &kFitSinSq2Theta13};

        //if (th13Surf) fitvars = {&kFitDeltaInPiUnits, &kFitSinSqTheta23, &kFitDmSq32Scaled};

      MinuitFitter fit23(exptThis, fitvars, systs);
      auto thisminchi = fit23.Fit(calc,auxShifts , SeedList({
                                  {&kFitDmSq32Scaled, {hie*fabs(kFitDmSq32Scaled.GetValue(calc))}},
                                  {&kFitSinSqTheta23, th23_old_seeds},
                                  {&kFitDeltaInPiUnits, delta_seeds}
                                  }),seedShifts)->EvalMetricVal();
      minchi23test= min(minchi23test, thisminchi);
      ResetOscCalcToDefault(calc);
      auxShifts.ResetToNominal();
  }//end hie
  std::cout << "\nFound overall test minchi " << minchi23test << std::endl;

  TFile * outfile = new TFile(outfilename+".root","recreate");
  TFile * outFCfile = new TFile(outFCfilename+".root","recreate");

  outfile->cd();
  TVectorD v(1);
  v[0] = minchi23;
  cout<<"minchi wrote into the file "<<minchi23<<endl;
  v.Write("overall_minchi");

  int steps = 30;//10;// 5;
  if(dmsqSurf){
     std::cerr << "\n WARNING Using 20 bins for dmsq surface\n\n";
     steps = 20;
  }
  //Now create all the surfaces
  for(int hie: {-1, +1}){

      if((onlyNH && hie<0) || (onlyIH && hie>0)) continue;

      std::cout << "Starting surface " << (hie>0? "NH ": "IH") << "\n\n";
      cout<<"the delta seeds are ";
      for (auto const & s:delta_seeds) std::cout << s << ", ";
      std::string hieStr = hie>0 ? "NH" : "IH";

      if(!th13Surf && ! dmsqSurf){
         ResetOscCalcToDefault(calc);
         calc->SetDmsq32(hie*fabs(calc->GetDmsq32()));
         double low = (nueOnly)?0:0.3;
         double high = (nueOnly)?1:0.7;
         cout<<"low "<<low<<" high "<<high<<endl;
         FrequentistSurface surf23(exptThis, calc,
                                   &kFitDeltaInPiUnits,steps, 0, 2,
                                   &kFitSinSqTheta23,  steps, low, high,
                                   {&kFitDmSq32Scaled,&kFitSinSq2Theta13},systs,
                                   SeedList(), seedShifts
                                  );

         auto surf1=surf23.ToTH2(minchi23);
         outfile->cd();
         surf1->Write((TString)"delta_th23_"+hieStr);
         surf23.SaveTo(outfile, "surf_delta_th23_"+hieStr);

         outFCfile->cd();
         std::vector<TH2*> profhists = surf23.GetProfiledHists();
         //Osc parameters are special case
         profhists[0]->Write((hieStr+"_DmSq32").c_str());
         profhists[1]->Write((hieStr+"_SinSq2Theta13").c_str());
         for (int i = 0; i < int(systs.size()); i++){
              profhists[i+2]->Write((hieStr+"_"+systs[i]->ShortName()).c_str());
         }
      }
      if(!th13Surf && dmsqSurf){
         calc = DefaultOscCalc();
         calc->SetDmsq32(hie*fabs(calc->GetDmsq32()));
         double lowth23 = (RHCOnly)?0.2:0.35;
         double highth23 = (RHCOnly)?1:0.7;
         cout<<"lowth23 "<<lowth23<<" highth23 "<<highth23<<endl;
         double lowdmsq = (hie>0?(RHCOnly?2.0:2.1):(RHCOnly?-3.0:-2.7));
         double highdmsq = (hie>0?(RHCOnly?3.0:2.7):(RHCOnly?-2.0:-2.2));
         cout<<"lowdmsq "<<lowdmsq<<" highdmsq "<<highdmsq<<endl;
         FrequentistSurface surf23m(exptThis, calc,
                                    &kFitSinSqTheta23, steps,lowth23, highth23,
                                    &kFitDmSq32Scaled, steps, lowdmsq, highdmsq,
                                    {&kFitSinSq2Theta13, &kFitDeltaInPiUnits},systs,
                                    SeedList({{&kFitDeltaInPiUnits, delta_seeds}}));
         auto surf6=surf23m.ToTH2(minchi23);
         outfile->cd();
         surf6->Write((TString)"th23_dm32_"+(hie>0?"NH":"IH"));
         surf23m.SaveTo(outfile, std::string("surf_th23_dm32_")+(hie>0?"NH":"IH"));

         outFCfile->cd();
         std::vector<TH2*> profhists = surf23m.GetProfiledHists();
         //Osc parameters are special case
         profhists[0]->Write((hieStr+"_SinSq2Theta13").c_str());
         profhists[1]->Write((hieStr+"_DeltaCPpi").c_str());
         for (int i = 0; i < int(systs.size()); i++){
              profhists[i+2]->Write((hieStr+"_"+systs[i]->ShortName()).c_str());
         }

      }
      if(th13Surf){
         calc = DefaultOscCalc();
         calc->SetDmsq32(hie*fabs(calc->GetDmsq32()));

         FrequentistSurface surf13(exptThis, calc,
                                   &kFitSinSq2Theta13,  steps, 0, (RHCOnly?1:0.35),
                                   &kFitDeltaInPiUnits,steps, 0, 2,
                                   {&kFitSinSqTheta23, &kFitDmSq32Scaled}, systs,
                                   SeedList({{&kFitSinSqTheta23, th23_all_seeds_NH}}), seedShifts);
         auto surf4 = surf13.ToTH2(minchi23);
         outfile->cd();
         surf4->Write((TString)"th13_delta_"+(hie>0?"NH":"IH"));
         surf13.SaveTo(outfile, std::string("surf_th13_delta_")+(hie>0?"NH":"IH"));

         outFCfile->cd();
         std::vector<TH2*> profhists = surf13.GetProfiledHists();
         //Osc parameters are special case
         profhists[0]->Write((hieStr+"_SinSqTheta23").c_str());
         profhists[1]->Write((hieStr+"_DmSq32").c_str());
         for (int i = 0; i < int(systs.size()); i++){
              profhists[i+2]->Write((hieStr+"_"+systs[i]->ShortName()).c_str());
         }

      }
  } //end hie

  delete outfile;
  std::cout << "\nSurfaces saved to " << outfilename << std::endl;
  return;

}