MakeCAFSensitivities_for_FNEX.C

Go to the documentation of this file.

#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Analysis/Calcs.h"
#include "CAFAna/Analysis/Exposures.h"
#include "CAFAna/Fit/MinuitFitter.h"
#include "CAFAna/Fit/FrequentistSurface.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/Prediction/PredictionExtrap.h"
#include "CAFAna/Vars/FitVars.h"
#include "3FlavorAna/Ana2018/nue/joint_fit_2018_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 "TH2.h"
#include "TStyle.h"
#include "TFile.h"

void MakeCAFSensitivities_for_FNEX()
{

  // Set the conditions we want sensitivites for
  std::string decomp = "prop";    // Extrapolation tye
  bool corrSysts     = false;     // Systematics
  bool debug         = true;      // Data/MC debug plots
  bool th13Surf      = false;     // Th13 surface
  bool dmsqSurf      = false;     // Dmsq32 surface
  bool nueOnly       = true ;     // NuE only analysis
  bool numuOnly      = false;     // NuMu only analysis
  bool joint         = false;     // NuMu-NuE joint analyses
  bool onlyNH        = true;      // Normal hierarchy
  bool fakedata      = true ;     // Fake FD Data
  bool realdata      = false;     // Real FD Data
  bool FHCOnly       = false;     // Neutrino mode only
  bool RHCOnly       = false;     // Anti-neutrino mode only
  bool FHCRHC        = true ;     // Both neutrino and anti-neutrino modes

  bool onlyIH        = false;     // Inverted hierarchy

  if(!onlyNH)  onlyIH    = 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 = "NO_ANALYSIS_MENTIONED";

  // The output directory
  TString outdir     = "/nova/app/users/prabhjot/ReadCAFMakeEventList/development/results/sensitivities/my_macro/";
  TString outfilename   (outdir + "hist_contours_2018_"    + decomp + "_" + savefile + "_2018");

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

  // Set the values of oscillation paraeters for fake FD Data
       if(fakedata && onlyNH) SetFakeCalc(calc_fake, 0.594,  2.49134e-3, 0.805  );
  else if(fakedata && onlyIH) SetFakeCalc(calc_fake, 0.588, -2.51776e-3, 1.507  );

  // Predictions for the NuE only or joint analyses
  if(!numuOnly){
    if(FHCOnly || FHCRHC){
      preds   .push_back(GetNuePrediction2018(decomp, calc_fake, corrSysts, "fhc", false));
      cosmics .push_back(GetNueCosmics2018("fhc"));
    }
    if(RHCOnly || FHCRHC){
      preds   .push_back(GetNuePrediction2018(decomp, calc_fake, corrSysts, "rhc", false));
      cosmics .push_back(GetNueCosmics2018("rhc"));
    }

    if(realdata){
      if(FHCOnly || FHCRHC){
        auto nue_data = GetNueData2018("fhc");
        data.insert(data.end(), nue_data);
      }
      if(RHCOnly || FHCRHC){
        auto nue_data = GetNueData2018("rhc");
        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;

  // Predictions for the NuMu only or joint analyses
  if(!nueOnly){
    if(FHCOnly || FHCRHC){
      auto numu_preds   = GetNumuPredictions2018(nnumu, corrSysts, "fhc");
      preds.insert(preds.end(), numu_preds.begin(), numu_preds.end());
      auto numu_cosmics = GetNumuCosmics2018(nnumu, "fhc");
      cosmics.insert(cosmics.end(), numu_cosmics.begin(), numu_cosmics.end());
    }
    if(RHCOnly || FHCRHC){
      auto numu_preds   = GetNumuPredictions2018(nnumu, corrSysts, "rhc");
      preds.insert(preds.end(), numu_preds.begin(), numu_preds.end());
      auto numu_cosmics = GetNumuCosmics2018(nnumu, "rhc");
      cosmics.insert(cosmics.end(), numu_cosmics.begin(), numu_cosmics.end());
    }
    if(realdata){
      if(FHCOnly || FHCRHC){
        auto numu_data = GetNumuData2018(nnumu, "fhc");
        data.insert(data.end(), numu_data.begin(), numu_data.end());
      }
      if(RHCOnly || FHCRHC){
        auto numu_data = GetNumuData2018(nnumu, "rhc");
        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){

    // Get the POTs
    double           kPOT = 1.0;
         if(FHCOnly) kPOT = kAna2018FHCPOT;
    else if(RHCOnly) kPOT = kAna2018RHCPOT;
    else if(FHCRHC)  kPOT = kAna2018FHCPOT + kAna2018RHCPOT;

    Spectrum* thisdata;
    if(fakedata) thisdata = GetFakeData(preds[i], calc_fake, kPOT, 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
  std::vector <const IFitVar *> fitvars = {&kFitDeltaInPiUnits, &kFitSinSqTheta23, &kFitDmSq32, &kFitSinSq2Theta13  };
  if(th13Surf)                  fitvars = {&kFitDeltaInPiUnits, &kFitSinSqTheta23, &kFitDmSq32                      };

  MinuitFitter fit23(exptThis, fitvars, 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, 0.5, 1, 1.5}}
                                                   },
                                  seedShifts)->EvalMetricVal();

      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 MakeCAFSensitivities_for_FNEX