fake_future_data.C

Go to the documentation of this file.

#include "CAFAna/Analysis/Calcs.h"
#include "CAFAna/Analysis/Exposures.h"
#include "CAFAna/Fit/MinuitFitter.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 "OscLib/IOscCalc.h"
#include "./sensitivity_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 fake_future_data(int num = 1, bool throwexp = true, bool corrSysts = true,
                      TString options="fake2018")
{
  // arguments:
  // num - number of fake experiments you'd like to study
  // throwexp - generation part, if false - read from the file and analyse 
  // corrSysts - stat or syst
  // options - which osc.params to use for fake data

    bool fake2017 = options.Contains("fake2017");
    bool fake2018 = options.Contains("fake2018");

    auto filename = "throw_exp_"+options + (corrSysts?"_with_syst_":"_just_stat_")+std::to_string(num)+".root";  
 
    //////////////////////////////////////////////////
    // Load Nue and Numu experiments
    //////////////////////////////////////////////////


    std::vector <const IPrediction * > preds;
    std::vector <std::pair <TH1D*, double > > cosmics;
    std::vector <const IExperiment * > expts;
    double POT, Livetime;

    auto calc_fake = DefaultOscCalc();
    double th23, dmsq, dcp;
    if(fake2017) {th23 = 0.51; dmsq = 2.44e-3; dcp = 1.21;}
    //SetFakeCalc(calc_fake, 0.51, 2.44e-3, 1.21);
    if(fake2018) {th23 = 0.58; dmsq = 2.51e-3; dcp = 0.17;}
    //SetFakeCalc(calc_fake, 0.58, 2.51e-3, 0.17);

    //nue fhc
    preds.push_back(GetNuePrediction2018("combo", DefaultOscCalc(), corrSysts, "fhc", false));
    cosmics.push_back(GetNueCosmics2018("fhc")); 
    //nue rhc
    preds.push_back(GetNuePrediction2018("prop", DefaultOscCalc(), corrSysts, "rhc", false));    //make decomp choosable?
    cosmics.push_back(GetNueCosmics2018("rhc"));        
    //fhc numu
    int nnumu = 4;
    auto numu_preds = GetNumuPredictions2018(nnumu, corrSysts, "fhc");
    preds.insert(preds.end(),numu_preds.begin(), numu_preds.end());
    auto numu_cosmics = GetNumuCosmics2018(nnumu, "fhc");                     //leave numu cosmic as it is for 2018 ana, probably need to make own cosmic for 1600s of livetime
    cosmics.insert(cosmics.end(),numu_cosmics.begin(), numu_cosmics.end());
    //numu rhc
    numu_preds = GetNumuPredictions2018(nnumu, corrSysts, "rhc");
    preds.insert(preds.end(),numu_preds.begin(), numu_preds.end());
    numu_cosmics = GetNumuCosmics2018(nnumu, "rhc");
    cosmics.insert(cosmics.end(),numu_cosmics.begin(), numu_cosmics.end());     
    
   // make fake data with random osc parameters based on the 2018 fit and random syst.shift
   // then apply stat. fluctuations to the MC and cosmic

   if(throwexp) {

    // Randomize the oscillation parameters, this process will base on the result contous and slices. Use them as inputs and choose random points on the surface/slice.
    // contours - th23 vs delta (take degeneracies into account)
    // slice - dmsq23
   
     TFile* file = new TFile(filename,"recreate"); 

     TFile* histFHCNH = new TFile("/nova/ana/nu_e_ana/Ana2018/Results/FHCOnly/contours/delta_th23/syst/hist_contours_2018_joint_realData_FHCOnly_onlyNHcombo_systs.root","read");
     TFile* histFHCIH = new TFile("/nova/ana/nu_e_ana/Ana2018/Results/FHCOnly/contours/delta_th23/syst/hist_contours_2018_joint_realData_FHCOnly_onlyIHcombo_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*)histFHCNH->Get("delta_th23_NH");
     auto histIH = (TH2F*)histFHCIH->Get("delta_th23_IH");

     //make "pdf" from the chi2 surface
     //
     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("plots/pdf_transform_hist2d_delta_th23_nh.pdf"); 
    
     new TCanvas();
     histIH_pdf->Draw("colz");
     gPad->Print("plots/pdf_transform_hist2d_delta_th23_ih.pdf");
    

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

     // should count the probability to get the ih and nh, also bases on the "pdf". And then throw the hierarchy as uniform distribution

     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* dmsqFHC = new TFile("/nova/ana/nu_e_ana/Ana2018/Results/FHCOnly/slices/syst/hist_slices_2017_joint_realData_FHCOnlycombo_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*)dmsqFHC->Get("slice_dmsq_NH");
     auto histdmsqIH = (TH1F*)dmsqFHC->Get("slice_dmsq_IH");

     // count "pdf" for slice too

     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("plots/pdf_transform_dmsq32_nh.pdf");

     new TCanvas();
     histdmsqIH_pdf->Draw("hist");
     gPad->Print("plots/pdf_transform_dmsq32_ih.pdf");
     
     //start doing experiments, j will count the experiment number

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

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

          TRandom3 rnd(0);

          //throw hierarchy as x in (0,1), compare to the probabilities we got from the surface's pdf

          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;
          // get random osc. parameters from the surface and slice
          if (isNH) {
              histNH_pdf->GetRandom2(thisdcp, thisth23);
              thisdmsq= histdmsqNH_pdf->GetRandom();
          }
          else{
              histIH_pdf->GetRandom2(thisdcp, thisth23);
              thisdmsq = histdmsqIH_pdf->GetRandom();
          }
          // get random th13 for the pdg constraint
          thisth13 = rnd.Gaus(0.082, 0.004);

          // use fake osc. parameters for this experiment  

          calc_fake->SetdCP(thisdcp*M_PI);
          calc_fake->SetTh23(asin(sqrt(thisth23)));
          calc_fake->SetDmsq32(thisdmsq);
          calc_fake->SetTh13(asin(sqrt(thisth13))/2);
          // save to the .root file just in case
          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;

          // randomize the systematic shiftes for this experiment
 
          struct systshifts{
                  double shift;
                  TString systname;
          };
 
          std::vector <systshifts> shifts;

          auto systs = GetJointFitSystematicList(corrSysts, false, false, true, true, true);
          // save to the file and make a lable in order to get the correct shift later
          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";

          // enter the loop for all spectra we analyze, nue fhc, nue rhc, 4 numu fhc, 4 numu rhc = 10 spectra we need to randomize

         for(int i=0; i < 10; i++){
             

             if(i==0) {POT = kFutureFHCPOT; Livetime = kFutureFHCLivetime;}
             if(i==1) {POT = kFutureRHCPOT; Livetime = kFutureRHCLivetime;}
             if(i >= 2 && i < 6) {POT = kFutureFHCPOT; Livetime = kFutureFHCLivetime;}
             if(i >= 6 && i < 10) {POT = kFutureRHCPOT; Livetime = kFutureRHCLivetime;}

             auto check_wo_syst = preds[i]->Predict(calc_fake).ToTH1(POT);
             cout<<" before random syst shifts the integral was "<<check_wo_syst->Integral()<<endl;

             // get systematics for this particular spectra (remember we have some systematics unique for nue/numu or fhc/rhc)
             // and assign them syst shifts we got earlier for this experiment j
             // if systematics are the same for some nue/numu/rhc/fhc spectra the shifts should be the same for them

             SystShifts sshifts;

             if(corrSysts){
                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
              }
             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";
             // make prediction based on the random parameters and random syst shiftes
             auto temphist = preds[i]->PredictSyst(calc_fake, sshifts).ToTH1(POT); // get full spectra for i experiment (make by components and then sum them up?)
             // we have an issue with numu cosmic prediction which was done for the 2018 exposure, so we need to scale it to the future exposure
             auto tempcosm = (TH1D*)cosmics[i].first->Clone(UniqueName().c_str());
             if(i>1) {
                double denom;
                if(i >= 2 && i < 6)  denom = kAna2018FHCLivetime;
                if(i >= 6 && i < 10) denom = kAna2018RHCLivetime;
                cout<<"numu file was made with "<<denom<<" livetime, need to scale to the "<<Livetime<<endl; 
                cout<<"before scaling "<<tempcosm->Integral()<<endl;
                tempcosm->Scale(Livetime/denom);
                cout<<"after scaling "<<tempcosm->Integral()<<endl;
             }
             //sum cosmic and mc prediction
             auto total = Spectrum(tempcosm, POT, 0);
             total+=Spectrum(temphist, POT,0);

             cout<<"mock experiment integral is "<<temphist->Integral()<<" and cosmic is "<<tempcosm->Integral()<<endl;
             //apply poissonian fluctuation
             auto mockexp = total.MockData(POT); 
             cout<<"mock experiment after Pois. fluct "<<mockexp.Integral(POT)<<endl;
             // save to the file; we have fake data now
             mockexp.SaveTo(d, (std::to_string(i)+"_mockexp").c_str());
        }
     }
  }


  if(!throwexp){ 

    // the part for the analyzing the fake data

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

   // loop over all experiments in file

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

       for(int i = 0; i < 10; ++i){
           // read fake data
           auto thisdata = LoadFromFile<Spectrum>(file, TString::Format("ex_%s/%s_mockexp", std::to_string(j).c_str(), std::to_string(i).c_str()).Data()).release();
           if(i==0) {POT = kFutureFHCPOT; Livetime = kFutureFHCLivetime;}
           if(i==1) {POT = kFutureRHCPOT; Livetime = kFutureRHCLivetime;}
           if(i >= 2 && i < 6) {POT = kFutureFHCPOT; Livetime = kFutureFHCLivetime;}
           if(i >= 6 && i < 10) {POT = kFutureRHCPOT; Livetime = kFutureRHCLivetime;}
           // scale the cosmic again, for making the prediction for fit
           auto tempcosm = (TH1D*)cosmics[i].first->Clone(UniqueName().c_str());
           if(i>1) {
                double denom;
                if(i >= 2 && i < 6)  denom = kAna2018FHCLivetime;
                if(i >= 6 && i < 10) denom = kAna2018RHCLivetime;
                tempcosm->Scale(Livetime/denom);
             }
           cout<<i<<" POT "<<POT<<" tot MC "<<preds[i]->Predict(calc_fake).Integral(POT)<<" cos "<<tempcosm->Integral()<<" cos er "<<1/sqrt(tempcosm->Integral())<<" analyze data "<<thisdata->Integral(POT)<<endl;
          // push mc prediction, fake data, scaled cosmic and its error for the fit
           expts.push_back(new SingleSampleExperiment(preds[i],*thisdata, tempcosm,1/sqrt(tempcosm->Integral())));
       }

       ////////////////////////////////////////////////////////////
       // Add constraints, make experiments
       ////////////////////////////////////////////////////////////


        std::cout << "\nCreating multiexperiment\n" << std::endl;

        expts.push_back(WorldReactorConstraint2017()); std::cout << "Adding WorldReactorConstraint2017\n";
        std::cout << "Creating Multiexperiment with a total of "
              << expts.size() << " experiments\n\n" << std::endl;
        auto exptThis = new MultiExperiment(expts);

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

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

         std::cout << "\n\nSystematic correlations...\n";

         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 << "Starting the fit" << std::endl;

         osc::IOscCalcAdjustable* calc = DefaultOscCalc();

         SystShifts auxShifts = SystShifts::Nominal();
         std::vector <SystShifts> seedShifts = {};
   
         // make seeds for the fit

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

          std::vector <th23helper> th23seeds;

          bool octantSlice =true;

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

          // start searching the bf by looping the NH/IH, upper and lower octants. The lowest chi2 will be the best fit

          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};
                   MinuitFitter fit23(exptThis, fitvars, systs);
                   auto thisminchi = fit23.Fit(calc, auxShifts,
                                    {
                                      { &kFitDmSq32Scaled,   {hie*2.5} },
                                      { thseed.var,          thseed.seeds },
                                      { &kFitDeltaInPiUnits, delta_seeds }
                                    }, seedShifts,
                                    IFitter::kQuiet)->EvalMetricVal();
                     minchi23= min(minchi23, thisminchi);

                     // plot the systamtic pulls for this octant and hierarchy
                     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->Update();
                        TLine *l=new TLine(gPad->GetUxmin(),0,gPad->GetUxmax(),0);
                        l->Draw("same"); 
                        gPad->Print("syst_shifts_for_exp_" +std::to_string(j)+
                                  (hie>0? "_NH": "_IH") + thseed.label + ".pdf");
                      }

                     ResetOscCalcToDefault(calc);
                     auxShifts.ResetToNominal();
              }
          }
          std::cout << "\nFound overall minchi " << minchi23 << "\n\n";
          expts.clear();
   }
  }
}