nuebar_signif.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 "OscLib/IOscCalc.h"
#include "3FlavorAna/Ana2018/nue/joint_fit_2018_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 nuebar_signif(int num = 1, bool throwexp = true, bool corrSysts = false,
                  TString options="fake2017_part1")
{

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

  bool part1 = options.Contains("part1");
  bool part2 = options.Contains("part2");

  auto filename = "throw_exp_nuesignif_"+options + (corrSysts?"_with_syst_":"_just_stat_")+std::to_string(num)+".root";  
 
    //////////////////////////////////////////////////
    // Load Nue and Numu experiments
    //////////////////////////////////////////////////
    //need numu only for prestage seeds
    std::vector <const IPrediction * > preds;
    std::vector <std::pair <TH1D*, double > > cosmics;
    std::vector <Spectrum * > data;
    std::vector <const IExperiment * > expts;

    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);
    //else if(!realData) {std::cerr << "need setting for data\n"; exit(1);}                                   //don't forget about it later

    preds.push_back(GetNuePrediction2018("prop", DefaultOscCalc(), corrSysts, "rhc", false));    //make decomp choosable?
    cosmics.push_back(GetNueCosmics2018("rhc"));        
    data.push_back(GetNueData2018("rhc"));

/// Start with part 1 to get the norm coef.
 if(part1){ 
   if(throwexp) {

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

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

     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");
 
     double POT = kAna2018RHCPOT;
     int i = 0;

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

          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;
 
          cout<<"Set syst. shifts"<<endl;
          auto systs = getAllAna2018Systs(kNueAna2018, true, kRHC);

          SystShifts shifts;
          for (auto &s:systs){

              auto sh = rnd.Gaus(0, 1);
              cout<<s->ShortName()<<"("<<sh<<"), ";
              shifts.SetShift(s, sh);
              TVectorD vs(1);
              vs[0] = sh;
              vs.Write((std::to_string(j)+"_"+s->ShortName()).c_str());

           }
          cout<<"\n\n";

         auto tempnue = preds[i]->PredictComponentSyst(calc_fake, shifts, Flavors::kNuMuToNuE, Current::kCC, Sign::kNu).ToTH1(POT);
         auto tempnuebar =  preds[i]->PredictComponentSyst(calc_fake, shifts, Flavors::kNuMuToNuE, Current::kCC, Sign::kAntiNu).ToTH1(POT);
         auto temptau = preds[i]->PredictComponentSyst(calc_fake, shifts, Flavors::kAllNuTau, Current::kCC, Sign::kBoth).ToTH1(POT);
         auto tempbeam = preds[i]->PredictComponentSyst(calc_fake, shifts, Flavors::kNuEToNuE, Current::kCC, Sign::kBoth).ToTH1(POT);
         auto tempnc = preds[i]->PredictComponentSyst(calc_fake, shifts, Flavors::kAll, Current::kNC, Sign::kBoth).ToTH1(POT);
         auto tempcc = preds[i]->PredictComponentSyst(calc_fake, shifts, Flavors::kAllNuMu, Current::kCC, Sign::kBoth).ToTH1(POT);

        cout<<"mock components: nue "<<  tempnue->Integral()<<" nue bar "<<tempnuebar->Integral()<<" tau "<<temptau->Integral()<<" beam "<<tempbeam->Integral()<<" nc "<<tempnc->Integral()<<" cc "<<tempcc->Integral()<<endl;
        tempnue->Write((std::to_string(i)+"_nue").c_str());
        tempnuebar->Write((std::to_string(i)+"_nuebar").c_str());
        temptau->Write((std::to_string(i)+"_tau").c_str());
        tempbeam->Write((std::to_string(i)+"_beam").c_str());
        tempnc->Write((std::to_string(i)+"_nc").c_str());
        tempcc->Write((std::to_string(i)+"_cc").c_str());

   }
  }


  if(!throwexp){ 

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

    calc_fake->SetdCP(1.20786*M_PI);
    calc_fake->SetTh23(asin(sqrt(0.513445)));
    calc_fake->SetDmsq32(2.46343e-3);
    calc_fake->SetTh13(asin(sqrt(0.0820218))/2);

    double POT = kAna2018RHCPOT;
    auto ws2018fhc =  preds[i]->PredictComponent(calc_fake, Flavors::kNuMuToNuE, Current::kCC, Sign::kNu).ToTH1(POT);
   
    auto tot2018bkg = preds[i]->PredictComponent(calc_fake, Flavors::kAllNuTau, Current::kCC, Sign::kBoth).ToTH1(POT);
    auto beam2018 = preds[i]->PredictComponent(calc_fake, Flavors::kNuEToNuE, Current::kCC, Sign::kBoth).ToTH1(POT);
    tot2018bkg->Add(beam2018);
    auto nc2018 = preds[i]->PredictComponent(calc_fake, Flavors::kAll, Current::kNC, Sign::kBoth).ToTH1(POT);
    tot2018bkg->Add(nc2018);
    auto cc2018 = preds[i]->PredictComponent(calc_fake, Flavors::kAllNuMu, Current::kCC, Sign::kBoth).ToTH1(POT);

    auto Nws2018fhc = ws2018fhc->Integral();
    auto Ntot2018bkg = tot2018bkg->Integral();

    cout<<"2018 ws is "<<ws2018fhc->Integral()<<" and tot bkg is "<<tot2018bkg->Integral()<<endl;

    TH1F* ws = new TH1F("ws", "ws", 100, 0, 10);
    TH1F* bkg = new TH1F("bkg", "bkg", 100, 0, 10);

    TH1F* ws_norm = new TH1F("wsnorm", "ws norm", 100, 0, 10);
    TH1F* bkg_norm = new TH1F("bkgnorm", "bkg norm", 100, 0, 10);

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

      auto thisws = (TH1F*)file->Get(TString::Format("ex_%s/%s_nue", std::to_string(j).c_str(), std::to_string(i).c_str())); 
    
      auto thisbkg = (TH1F*)file->Get(TString::Format("ex_%s/%s_tau", std::to_string(j).c_str(), std::to_string(i).c_str()));

      auto thiscc = (TH1F*)file->Get(TString::Format("ex_%s/%s_cc", std::to_string(j).c_str(), std::to_string(i).c_str())); 
      thisbkg->Add(thiscc);
      auto thisnc = (TH1F*)file->Get(TString::Format("ex_%s/%s_nc", std::to_string(j).c_str(), std::to_string(i).c_str()));
      thisbkg->Add(thisnc);
      auto thisbeam = (TH1F*)file->Get(TString::Format("ex_%s/%s_beam", std::to_string(j).c_str(), std::to_string(i).c_str()));
      thisbkg->Add(thisbeam);

      auto Nws = thisws->Integral();
      auto Nbkg = thisbkg->Integral();

      ws->Fill(Nws);
      bkg->Fill(Nbkg);
    
      ws_norm->Fill(Nws/Nws2018fhc);
      bkg_norm->Fill(Nbkg/Ntot2018bkg);
            
   }

     new TCanvas;
     ws->Draw("hist");
     gPad->Print(("plots/events_ws_"+std::to_string(num)+".pdf").c_str());

     new TCanvas;
     bkg->Draw("hist");
     gPad->Print(("plots/events_bkg_"+std::to_string(num)+".pdf").c_str());

     new TCanvas;
     ws_norm->Draw("hist");
     gPad->Print(("plots/events_ws_norm_"+std::to_string(num)+".pdf").c_str());

     new TCanvas;
     bkg_norm->Draw("hist");
     gPad->Print(("plots/events_bkg_norm_"+std::to_string(num)+".pdf").c_str());
  
  }//end !throw

}//end the part 1

//start with part2 to get the signif. itself 
if(part2){
   
   if(throwexp){

    TFile* file = new TFile("throw_exp_nuesignif_fake2017_part1_with_syst_10000.root","read");
    int i = 0;

    calc_fake->SetdCP(1.20786*M_PI);
    calc_fake->SetTh23(asin(sqrt(0.513445)));
    calc_fake->SetDmsq32(2.46343e-3);
    calc_fake->SetTh13(asin(sqrt(0.0820218))/2);

    double POT = kAna2018RHCPOT;
    auto ws2018fhc =  preds[i]->PredictComponent(calc_fake, Flavors::kNuMuToNuE, Current::kCC, Sign::kNu).ToTH1(POT);

    auto tot2018bkg = preds[i]->PredictComponent(calc_fake, Flavors::kAllNuTau, Current::kCC, Sign::kBoth).ToTH1(POT);
    auto beam2018 = preds[i]->PredictComponent(calc_fake, Flavors::kNuEToNuE, Current::kCC, Sign::kBoth).ToTH1(POT);
    tot2018bkg->Add(beam2018);
    auto nc2018 = preds[i]->PredictComponent(calc_fake, Flavors::kAll, Current::kNC, Sign::kBoth).ToTH1(POT);
    tot2018bkg->Add(nc2018);
    auto cc2018 = preds[i]->PredictComponent(calc_fake, Flavors::kAllNuMu, Current::kCC, Sign::kBoth).ToTH1(POT);
    tot2018bkg->Add(cc2018);

    auto Nws2018fhc = ws2018fhc->Integral();
    auto Ntot2018bkg = tot2018bkg->Integral();

    cout<<"2018 ws is "<<ws2018fhc->Integral()<<" and tot bkg is "<<tot2018bkg->Integral()<<endl;

    TH1F* ws_norm = new TH1F("wsnorm", "ws norm", 50, 0, 10);
    TH1F* bkg_norm = new TH1F("bkgnorm", "bkg norm", 100, 0, 10);

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

      auto thisws = (TH1F*)file->Get(TString::Format("ex_%s/%s_nue", std::to_string(j).c_str(), std::to_string(i).c_str()));

      auto thisbkg = (TH1F*)file->Get(TString::Format("ex_%s/%s_tau", std::to_string(j).c_str(), std::to_string(i).c_str()));
      auto thiscc = (TH1F*)file->Get(TString::Format("ex_%s/%s_cc", std::to_string(j).c_str(), std::to_string(i).c_str()));
      thisbkg->Add(thiscc);
      auto thisnc = (TH1F*)file->Get(TString::Format("ex_%s/%s_nc", std::to_string(j).c_str(), std::to_string(i).c_str()));
      thisbkg->Add(thisnc);
      auto thisbeam = (TH1F*)file->Get(TString::Format("ex_%s/%s_beam", std::to_string(j).c_str(), std::to_string(i).c_str()));
      thisbkg->Add(thisbeam);

      auto Nws = thisws->Integral();
      auto Nbkg = thisbkg->Integral();

      ws_norm->Fill(Nws/Nws2018fhc);
      bkg_norm->Fill(Nbkg/Ntot2018bkg);

   }

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

   auto tempnue = preds[i]->PredictComponent(calc_fake, Flavors::kNuMuToNuE, Current::kCC, Sign::kNu).ToTH1(POT);

   auto tempbkg = preds[i]->PredictComponent(calc_fake, Flavors::kAllNuTau, Current::kCC, Sign::kBoth).ToTH1(POT);
   auto tempbeam = preds[i]->PredictComponent(calc_fake, Flavors::kNuEToNuE, Current::kCC, Sign::kBoth).ToTH1(POT);
   tempbkg->Add(tempbeam);
   auto tempnc = preds[i]->PredictComponent(calc_fake, Flavors::kAll, Current::kNC, Sign::kBoth).ToTH1(POT);
   tempbkg->Add(tempnc);
   auto tempcc = preds[i]->PredictComponent(calc_fake, Flavors::kAllNuMu, Current::kCC, Sign::kBoth).ToTH1(POT);
   tempbkg->Add(tempcc);

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

         auto thistempnue = (TH1*)tempnue->Clone(UniqueName().c_str());
         auto thistempbkg = (TH1*)tempbkg->Clone(UniqueName().c_str());

         auto Nws =  ws_norm->GetRandom();         
         auto Nbkg = bkg_norm->GetRandom();

         cout<<" Nws "<<Nws<<" before scale "<<thistempnue->Integral();
         thistempnue->Scale(Nws);       
         cout<<" after scale "<<thistempnue->Integral()<<endl;
 
         cout<<" Nbkg "<<Nbkg<<" before scale "<<thistempbkg->Integral();
         thistempbkg->Scale(Nbkg);
         cout<<" after scale "<<thistempbkg->Integral()<<endl;

         cout<<"cosmic "<<cosmics[i].first->Integral()<<endl;

         auto total = Spectrum(cosmics[i].first, POT,0);
         total+=Spectrum(thistempnue, POT,0);
         total+=Spectrum(thistempbkg, POT,0);
         
         auto mockexp = total.MockData(POT);       
         cout<<"mock exp "<<mockexp.Integral(POT)<<endl;
         mockexp.SaveTo(writefile, (std::to_string(k)+"_mockexp").c_str());      

    }
  } // end throw fake2

  if(!throwexp){

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

    int i = 0;

    calc_fake->SetdCP(1.20786*M_PI);
    calc_fake->SetTh23(asin(sqrt(0.513445)));
    calc_fake->SetDmsq32(2.46343e-3);
    calc_fake->SetTh13(asin(sqrt(0.0820218))/2);

    double POT = kAna2018RHCPOT;
    auto ws2018fhc =  preds[i]->PredictComponent(calc_fake, Flavors::kNuMuToNuE, Current::kCC, Sign::kNu).ToTH1(POT);

    auto tot2018bkg = preds[i]->PredictComponent(calc_fake, Flavors::kAllNuTau, Current::kCC, Sign::kBoth).ToTH1(POT);
    auto beam2018 = preds[i]->PredictComponent(calc_fake, Flavors::kNuEToNuE, Current::kCC, Sign::kBoth).ToTH1(POT);
    tot2018bkg->Add(beam2018);
    auto nc2018 = preds[i]->PredictComponent(calc_fake, Flavors::kAll, Current::kNC, Sign::kBoth).ToTH1(POT);
    tot2018bkg->Add(nc2018);
    auto cc2018 = preds[i]->PredictComponent(calc_fake, Flavors::kAllNuMu, Current::kCC, Sign::kBoth).ToTH1(POT);
    tot2018bkg->Add(cc2018);
    tot2018bkg->Add(cosmics[i].first);
    tot2018bkg->Add(ws2018fhc);
    cout<<"total 2018 bkg is "<<tot2018bkg->Integral();

    new TCanvas;
   
    auto realdata = data[i]->ToTH1(POT);

    realdata->SetLineColor(kGreen+2);
    realdata->Draw("hist");
    tot2018bkg->Draw("hist same");
    gPad->Print("plots/check_the_2018data_bkg.pdf");

    cout<<" total 2018 rhc data is "<<realdata->Integral();
    auto LLreal = LogLikelihood(tot2018bkg, realdata);
    cout<<" likelihood btw real data and fhc predicted bkg is "<<LLreal<<endl;

   TH1F* chi2 = new TH1F("chi2", "chi2", 40, 0, 40);
   TH1F* integral = new TH1F("integrals", "integrals", 40, 0, 40);

   for (int j = 0; j<num; j++){
       auto thisdata = LoadFromFile<Spectrum>(file, TString::Format("%s_mockexp", std::to_string(j).c_str()).Data()).release()->ToTH1(POT);
       auto thisLL = LogLikelihood(tot2018bkg, thisdata);
       auto integr = thisdata->Integral();
       integral->Fill(integr);
       chi2->Fill(thisLL);
   }

   TFile* writefile = new TFile("finalhisti_"+options+".root","recreate");   
   chi2->Write("histchi");
   writefile->Close();

   cout<<"integral for chiN > LLreal "<<chi2->Integral((int)LLreal,40);

   new TCanvas;
   integral->Draw("hist");
   gPad->Print("plots/integrals_"+options+"_"+std::to_string(num)+"_exp.pdf");

   new TCanvas;
   chi2->Draw("hist"); 
   auto line = new TLine(LLreal, 0, LLreal, 5000);
   line->SetLineColor(kRed);
   line->Draw();
   gPad->Print("plots/chisq_with_"+options+"_"+std::to_string(num)+"_exp.pdf");
  }
 }
}