joint_fit_future_contour_univ.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/Fit/FrequentistSurface.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 "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 "DDTUtilities/func/MathUtil.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 joint_fit_future_contour_univ(int universes=2,
                                   TString uni_tag = "",
                                   double deltaCP = 0.18,
                                   double dmsq32 = 2.51e-3,
                                   double th23 = 0.58,
                                   string sign = "plus",
                                   double FHCexp = kAna2020FHCPOT,
                                   double RHCexp = kAna2020RHCPOT,
                                   double FHClive = kAna2020FHCLivetime,
                                   double RHClive = kAna2020RHCLivetime,
                                   bool corrSysts = false,
                                   TString options="joint_mockData_both_onlyIH",
                                   bool dmsqSurf = false,
                                   bool th13Surf = false)
{

  cout<<"\nLOAD PREDICTIONS\n\n"<<endl;
  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 onlyNH = options.Contains("onlyNH");
  bool onlyIH = options.Contains("onlyIH");

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

  string syst = (corrSysts?"syst":"stat");

  TString outfilename   (syst + "_hist_contours_"    + suffix + "_" + uni_tag);
  TFile * outfile = new TFile(outfilename+".root","recreate");

  //////////////////////////////////////////////////
  // Load Nue and Numu experiments
  //////////////////////////////////////////////////
  // need numu only for prestage seeds
  // make map instead of struct in future

  std::string decomp = ""; // use Pt extrap all the times
  std::vector <ana::predictions> preds = LoadPredictions (corrSysts, true, decomp, nueOnly || joint, numuOnly || joint, FHCOnly || both, RHCOnly || both);
  std::vector <ana::predictions> preds_numu_only = LoadPredictions (corrSysts, true, decomp, false, numuOnly || joint, FHCOnly || both, RHCOnly || both);
  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 (sign =="minus") dmsq32 = -dmsq32;
  SetFakeCalc(calc_fake, th23, dmsq32, deltaCP);

  outfile->cd();
  TVectorD osc(3);
  osc[0] = th23;
  osc[1] = dmsq32;
  osc[2] = deltaCP;
  osc.Write("osc_params");

  for (int j = 0; j< universes; j++){
    cout<<"universe no "<<j<<endl;
    cout<<"clear expts and data vectors"<<endl;
    expts.clear();
    expts_numu_only.clear();
    data.clear();
    data_numu_only.clear();        
    // isn't set for individual fits (only joint nue+numu rhc+fhc)
    if(!numuOnly){
         if(FHCOnly || both ) data.push_back(GenerateFutureData(calc_fake, preds[0].pred, preds[0].cos.first, FHCexp, FHClive));
         if(RHCOnly || both ) data.push_back(GenerateFutureData(calc_fake, preds[1].pred, preds[1].cos.first, RHCexp, RHClive));
     }
    if(!nueOnly){
        if(FHCOnly || both ){
          std::vector <Spectrum> numu_data;
          numu_data.push_back(GenerateFutureData(calc_fake, preds[2].pred, preds[2].cos.first, FHCexp, FHClive));
          numu_data.push_back(GenerateFutureData(calc_fake, preds[3].pred, preds[3].cos.first, FHCexp, FHClive));
          numu_data.push_back(GenerateFutureData(calc_fake, preds[4].pred, preds[4].cos.first, FHCexp, FHClive));
          numu_data.push_back(GenerateFutureData(calc_fake, preds[5].pred, preds[5].cos.first, FHCexp, FHClive));
          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 ){
          std::vector <Spectrum> numu_data;
          numu_data.push_back(GenerateFutureData(calc_fake, preds[6].pred, preds[6].cos.first, RHCexp, RHClive));
          numu_data.push_back(GenerateFutureData(calc_fake, preds[7].pred, preds[7].cos.first, RHCexp, RHClive));
          numu_data.push_back(GenerateFutureData(calc_fake, preds[8].pred, preds[8].cos.first, RHCexp, RHClive));
          numu_data.push_back(GenerateFutureData(calc_fake, preds[9].pred, preds[9].cos.first, RHCexp, RHClive));
          data.insert(data.end(),numu_data.begin(), numu_data.end());
          data_numu_only.insert(data_numu_only.end(),numu_data.begin(), numu_data.end());
        }
    }

    cout<<"check all vector sizes: "<<preds.size()<<" data "<<data.size()<<endl;

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

    for(int i = 0; i < int(preds.size()); ++i){
      auto thisdata = data[i];
      double POT = thisdata.POT();
      double LT = thisdata.Livetime();
      thisdata.SaveTo(d, ("data_exp_"+std::to_string(i)).c_str());
      auto hcos = preds[i].cos.first->ToTH1(LT, 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(POT)<<endl;
      expts.push_back(new SingleSampleExperiment(preds[i].pred, thisdata, *preds[i].cos.first, preds[i].cos.second, true));
    }

    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){
      auto thisdata = data_numu_only[i];
      double POT = thisdata.POT();
      double LT = thisdata.Livetime();
      auto hcos = preds_numu_only[i].cos.first->ToTH1(LT, 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(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 WorldReactorConstraint2017\n";
      expts.push_back(WorldReactorConstraint2019());
    }
    if(nueOnly) {
      std::cout << "Adding Dmsq32ConstraintPDG2017\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);


    ////////////////////////////////////////////////////////////
    // Systematics
    ////////////////////////////////////////////////////////////
    std::cout << "Systematics for the fit:\n";
   
    bool ptExtrap = true;
    auto systs = GetJointFitSystematicList(corrSysts, nueOnly, numuOnly, FHCOnly||both, RHCOnly||both, true, ptExtrap);

    int nnumu = 4;
    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();
    // 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";
      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);

     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; double thisdmsq = -1; double thisth23 = -1; double thisdcp = -1; double thisth13 = -1;
   double curdmsq, curth23, curdcp, curth13;

   for(int hie:{-1,1}){
       for (auto & thseed: th23seeds){
        std::cout << "\n\nFinding best fit "
                  << thseed.label
                  << (hie>0? " NH " : " IH ")
                  << "\n\n";
        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};

        MinuitFitter fit23(exptThis, fitvars, systs);
        auto thisminchi = fit23.Fit(calc, auxShifts, SeedList({
                                   {&kFitDmSq32Scaled,{hie*pre_seed_dmsq}},
                                   {thseed.var, thseed.seeds},
                                   {&kFitDeltaInPiUnits, delta_seeds},
                                   {&kFitSinSq2Theta13, {0.082}} }),
                                   {}, IFitter::kVerbose
                                   );
        curdmsq = kFitDmSq32Scaled.GetValue(calc);
        curth13 = kFitSinSq2Theta13.GetValue(calc);
        curth23 = thseed.var->GetValue(calc);
        curdcp = kFitDeltaInPiUnits.GetValue(calc);
        TVectorD v(5);
        v[0] = thisminchi;
        v[1] = curdmsq;
        v[2] = curth23;
        v[3] = curdcp;
        v[4] = curth13;
        string hietag = (hie>0? "NH": "IH");
        v.Write("chi_"+std::to_string(j)+"_"+hietag+"_"+thseed.label);
        if(thisminchi<minchi23){
           thisdmsq = kFitDmSq32Scaled.GetValue(calc);
           thisth13 = kFitSinSq2Theta13.GetValue(calc);
           thisth23 = thseed.var->GetValue(calc);
           thisdcp = kFitDeltaInPiUnits.GetValue(calc);
           cout <<"this values (dmsq, th13, th23, dcp): "<<thisdmsq<<" "<<thisth13<<" "<<thisth23<<" "<<thisdcp<<endl;
        }
        minchi23= min(minchi23, thisminchi);
        ResetOscCalcToDefault(calc);
        auxShifts.ResetToNominal();
      }//end th23
   }//end hie
   std::cout << "\nFound overall minchi " << minchi23 << std::endl;

   TVectorD v(5);
   v[0] = minchi23;
   v[1] = thisdmsq;
   v[2] = thisth23;
   v[3] = thisdcp;
   v[4] = thisth13;
   cout<<"minchi wrote into the file "<<minchi23<<endl;
   v.Write(("bf_"+std::to_string(j)).c_str());

   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//,
                                  //{},seedShifts
                                  );

        auto surf1=surf23.ToTH2(minchi23);
        surf1->Write((TString)"delta_th23_univ"+std::to_string(j)+"_"+hieStr);
        surf23.SaveTo(d, (TString)"surf_delta_th23_univ"+std::to_string(j)+"_"+hieStr);
      }
      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)); 
        FrequentistSurface surf23m(exptThis, calc,
                                   &kFitSinSqTheta23, steps,lowth23, highth23,
                                   &kFitDmSq32Scaled, steps, lowdmsq, highdmsq,
                                   {&kFitSinSq2Theta13, &kFitDeltaInPiUnits},systs,
                                   SeedList({{&kFitDeltaInPiUnits, delta_seeds}}));
        auto surf6=surf23m.ToTH2(minchi23);
        surf6->Write((TString)"th23_dm32_univ"+std::to_string(j)+"_"+(hie>0?"NH":"IH"));
        surf23m.SaveTo(d, (TString)"surf_th23_dm32_univ"+std::to_string(j)+"_"+(hie>0?"NH":"IH"));

      }
      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_univ"+std::to_string(j)+"_"+(hie>0?"NH":"IH"));
        surf13.SaveTo(d, (TString)"surf_th13_delta_univ"+std::to_string(j)+"_"+(hie>0?"NH":"IH"));
      }
    } //end hie
   }

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

}