CAF_makeCAFSensitivities_for_FNEX.C

Go to the documentation of this file.

#include "CAFAna/Core/Loaders.h"
#include "CAFAna/Core/LoadFromFile.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/Utilities.h"
#include "CAFAna/Cuts/Cuts.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "CAFAna/Cuts/TimingCuts.h"

#include "CAFAna/Analysis/Calcs.h"
#include "CAFAna/Analysis/Exposures.h"
#include "CAFAna/Fit/MinuitFitter.h"
#include "3FlavorAna/Plotting/NuePlotStyle.h"
#include "CAFAna/Analysis/Style.h"
#include "CAFAna/Fit/FrequentistSurface.h"
#include "CAFAna/Extrap/ModularExtrap.h"
#include "CAFAna/Experiment/AtmConstraint.h"
#include "CAFAna/Experiment/Dmsq32Constraint.h"
#include "CAFAna/Experiment/GaussianConstraint.h"
#include "CAFAna/Experiment/MultiExperiment.h"
#include "CAFAna/Experiment/ReactorExperiment.h"
#include "CAFAna/Experiment/Ssth23Constraint.h"
#include "CAFAna/Experiment/SingleSampleExperiment.h"
#include "CAFAna/Experiment/CountingExperiment.h"
#include "CAFAna/FC/FCSurface.h"
#include "CAFAna/Prediction/PredictionCombinePeriods.h"
#include "CAFAna/Prediction/PredictionExtrap.h"
#include "CAFAna/Vars/FitVars.h"
#include "3FlavorAna/Ana2017/nue/joint_fit_2017_tools.h"

#include "CAFAna/Analysis/Plots.h"

using namespace ana;

#include "Utilities/rootlogon.C"
#include "OscLib/IOscCalc.h"

#include "TCanvas.h"
#include "TVectorD.h"
#include "TBox.h"
#include "TH2.h"
#include "TColor.h"
#include "TGraph.h"
#include "TLegend.h"
#include "TLine.h"
#include "TStyle.h"
#include "TFile.h"


void CAF_makeCAFSensitivities_for_FNEX()
{

  // Set the analysis
  std::string decomp = "prop";
  bool corrSysts     = false;
  bool debug         = true;
  bool th13Surf      = false;
  bool dmsqSurf      = false;
  bool numuOnly      = false;
  bool nueOnly       = false;
  bool joint         = true;
  bool onlyNH        = true;
  bool onlyIH        = false;
  bool fakedata      = false;
  bool realdata      = true;

  // Set the name of the output files
  TString savefile;
  if(fakedata && joint    && onlyNH) savefile = "sensitivity_numu_nue_joint_NH";
  if(fakedata && joint    && onlyIH) savefile = "sensitivity_numu_nue_joint_IH";
  if(fakedata && nueOnly  && onlyNH) savefile = "sensitivity_nueOnly_NH";
  if(fakedata && nueOnly  && onlyIH) savefile = "sensitivity_nueOnly_IH";
  if(fakedata && numuOnly && onlyNH) savefile = "sensitivity_numuOnly_NH";
  if(fakedata && numuOnly && onlyIH) savefile = "sensitivity_numuOnly_IH";
  if(realdata && joint    && onlyNH) savefile = "contours_numu_nue_joint_NH";
  if(realdata && joint    && onlyIH) savefile = "contours_numu_nue_joint_IH";
  if(realdata && nueOnly  && onlyNH) savefile = "contours_nueOnly_NH";
  if(realdata && nueOnly  && onlyIH) savefile = "contours_nueOnly_IH";
  if(realdata && numuOnly && onlyNH) savefile = "contours_numuOnly_NH";
  if(realdata && numuOnly && onlyIH) savefile = "contours_numuOnly_IH";
  //else                            savefile = "";

  // The output directory
  TString outdir     = "/nova/app/users/prabhjot/ReadCAFMakeEventList/";
  TString outfilename   (outdir + "hist_contours_2017_"    + decomp + "_" + savefile + "_2017");

  //We'll save to this file, make pretty plots later
  TFile * outfile = new TFile(outfilename+".root", "recreate");

  //////////////////////////////////////////////////
  //  Load Nue experiments
  //////////////////////////////////////////////////

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

       if(fakedata && onlyNH) SetFakeCalc(calc_fake, 0.477,    2.43e-3,   1.432);
  else if(fakedata && onlyIH) SetFakeCalc(calc_fake, 0.564422, -2.498e-3, 1.45913);

  if(!numuOnly){
    preds   .push_back(GetNuePrediction2017(decomp, calc_fake, corrSysts));
    cosmics .push_back(GetNueCosmics2017());
    if(realdata){
      auto nue_data = GetNueData2017();
      data.insert(data.end(), nue_data);
    }//end of condition that it is a real data
  }//end of condition that we have nue or joint fit

  int nnumu = 4;

  if(!nueOnly){
    auto numu_preds   = GetNumuPredictions2017(nnumu, corrSysts);
    preds.insert(preds.end(), numu_preds.begin(), numu_preds.end());
    auto numu_cosmics = GetNumuCosmics2017(nnumu);
    cosmics.insert(cosmics.end(), numu_cosmics.begin(), numu_cosmics.end());
    if(realdata){
      auto numu_data = GetNumuData2017(nnumu);
      data.insert(data.end(), numu_data.begin(), numu_data.end()); 
    }//end of condition that it is a real data 
  }//end of condition that it is not nue only case


  //get the data
  for(int i = 0; i < int(preds.size()); ++i){
    Spectrum* thisdata;
    if(fakedata) thisdata = GetFakeData (preds[i], calc_fake, kAna2017POT, cosmics[i].first);
    if(realdata) thisdata = data[i];
    expts.push_back(new SingleSampleExperiment(preds[i],
                                               *thisdata,
                                               cosmics[i].first,
                                               cosmics[i].second));

    if(debug){
      new TCanvas();
      DataMCComparison (*thisdata, preds[i], calc_fake);
      cosmics[i].first->SetMarkerStyle(34);
      cosmics[i].first->SetMarkerColor(kCosmicBackgroundColor);
      cosmics[i].first->Draw("hist p same");
      gPad->Print(outdir + "debug_predictions.pdf");
    }//end of debugging
  }//end of loop over predictions

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

  std::cout << "\nCreating multiexperiment\n" << std::endl;
  if(!th13Surf){
    std::cout << "Adding WorldReactorConstraint2017\n";
    expts.push_back(WorldReactorConstraint2017());
  }//end of condition to add th13 world constrainet

  if(nueOnly){
    std::cout << "Adding Dmsq32ConstraintPDG2017\n";
    expts.push_back(&kDmsq32ConstraintPDG2017);
  }//end of condition that it is a nue only analysis

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

  ////////////////////////////////////////////////////////////
  // Systematics - STATS only case
  ////////////////////////////////////////////////////////////
  std::vector<const ISyst*> systs     = {};                      //No systematics

  SystShifts auxShifts                = SystShifts::Nominal();   //Empty systematics - http://nusoft.fnal.gov/nova/novasoft/doxygen/html/classana_1_1SystShifts.html#a23fd2bdae34aaadf646f422d66513403

  std::vector <SystShifts> seedShifts = {};                      //Empty SystShifts seeds

  ////////////////////////////////////////////////////////////
  // Fit
  ////////////////////////////////////////////////////////////

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

  osc::IOscCalcAdjustable* calc = DefaultOscCalc();

  //Find the best fit
  MinuitFitter fit23(exptThis, {&kFitDeltaInPiUnits, &kFitSinSqTheta23, &kFitDmSq32, &kFitSinSq2Theta13}, systs);

  double minchi23 = 1E20;

  for(double thseed : {0.3, 0.7}){

    for(int hie:{-1, 1}){

      std::cout << "\n\nFinding best fit " << (thseed < 0.5 ? "LO " : "UO ") << (hie > 0 ? "NH " : "IH ") << "\n\n";

      //find minimum chisq
      auto thisminchi = fit23.Fit(calc, auxShifts ,{
                                                   {&kFitDmSq32, {hie*fabs(calc->GetDmsq32())} },
                                                   {&kFitSinSqTheta23, {thseed}},
                                                   {&kFitDeltaInPiUnits,{0.33, 1., 1.66}}
                                                   },
                                  seedShifts);

      minchi23 = min(minchi23, thisminchi);

      ResetOscCalcToDefault(calc);
      auxShifts.ResetToNominal();

    }//end of loop over NH and IH hierarchies
  }//end of loop over sinsqth23 values as seed

  std::cout << "\nFound overall minchi " << minchi23 << std::endl;

  outfile->cd();

  TVectorD v(1);
  v[0] = minchi23;
  v.Write("overall_minchi");

  TVectorD bests(4);
  bests[0]= calc_fake->GetdCP();
  bests[1]= calc_fake->GetTh13();
  bests[2]= calc_fake->GetTh23();
  bests[3]= calc_fake->GetDmsq32();
  bests.Write("best_dcp_th13_th23_dmsq32");

  int steps = 30;

  std::cout << "Using 30 bins for dmsq surface\n\n";

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

    //Draw surfaces only if you are not making surfaces for th13 and Dmsq32
    if(!th13Surf && ! dmsqSurf){

      ResetOscCalcToDefault(calc);
      calc->SetDmsq32(hie*fabs(calc->GetDmsq32()));

      FrequentistSurface surf23(exptThis, calc,
                     &kFitDeltaInPiUnits, steps, 0, 2,
                     &kFitSinSqTheta23,  steps, (nueOnly ? 0 : 0.3), (nueOnly ? 1 : 0.7),
                     {&kFitDmSq32, &kFitSinSq2Theta13}, systs,
                     {}, seedShifts);

       auto surf1 = surf23.ToTH2(minchi23);

       std::string hieStr = hie > 0 ? "NH" : "IH";

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

       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());
    }//end of condition that we don't want surfaces for th13 and Dmsq32

    //Draw surfaces if you are making surfaces for dmsq but not for th13
    if(!th13Surf && dmsqSurf){

      calc = DefaultOscCalc();
      calc->SetDmsq32(hie * fabs(calc->GetDmsq32()));

      FrequentistSurface surf23m(exptThis, calc,
                      &kFitSinSqTheta23, steps, 0.3, 0.7,
                      &kFitDmSq32, steps, (hie > 0 ? 2e-3 : -3.0e-3 ), (hie > 0 ? 3.0e-3 : -2e-3),
                      {&kFitSinSq2Theta13, &kFitDeltaInPiUnits}, systs,
                      {{&kFitDeltaInPiUnits, {0.5,1.5}}}, seedShifts);

      auto surf6=surf23m.ToTH2(minchi23);

      outfile->cd();

      surf6->Write((TString)"th23_dm32_"+(hie > 0 ? "NH" : "IH"));

      surf23m.SaveTo(outfile, (TString)"surf_th23_dm32_"+(hie > 0 ? "NH" : "IH"));
    }//end end of condition that we are making surfaces for dmsq but not for th13

    //Draw surfaces for th13 (or draw contours for th13)
    if(th13Surf){

      calc = DefaultOscCalc();
      calc->SetDmsq32(hie*fabs(calc->GetDmsq32()));

      FrequentistSurface surf13(exptThis, calc,
                     &kFitSinSq2Theta13,  steps, 0, 0.35,
                     &kFitDeltaInPiUnits, steps, 0, 2,
                     {&kFitSinSqTheta23, &kFitDmSq32}, systs,
                     {{&kFitSinSqTheta23, {0.3, 0.7}}}, seedShifts);

      auto surf4 = surf13.ToTH2(-1);

      outfile->cd();

      surf4->Write((TString)"th13_delta_"+(hie > 0 ? "NH": "IH"));
      surf13.SaveTo(outfile, (TString)"surf_th13_delta_"+(hie > 0 ? "NH" : "IH"));
    }//end of condition to make surfaces for th13

  }//end of loop over NH and IH hierarchies

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

  //Close file
  outfile->Close();
  std::cout << "Name of the outfile is " << outfilename + ".root" << std::endl;
}//end of main function CAF_makeCAFSensitivities_for_FNEX