xsec_extrap_plots.C

Go to the documentation of this file.

#include "CAFAna/Analysis/Calcs.h"
#include "CAFAna/Core/SystShifts.h"
#include "CAFAna/Core/Loaders.h"
#include "CAFAna/Cuts/Cuts.h"
#include "CAFAna/Cuts/NueCutsSecondAna.h"
#include "3FlavorAna/Cuts/NumuCuts.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "CAFAna/Extrap/ModularExtrap.h"
#include "3FlavorAna/Prediction/PredictionGenerator3Flavor.h"
#include "CAFAna/Prediction/PredictionExtrap.h"
#include "CAFAna/Prediction/PredictionNoExtrap.h"
#include "CAFAna/Systs/Systs.h"
#include "CAFAna/Vars/GenieWeights.h"
#include "CAFAna/Vars/HistAxes.h"
#include "CAFAna/Vars/PPFXWeights.h"
#include "temp_draw_plots_util.h"

#include <TFile.h>
#include <TString.h>
#include <vector>
#include <iostream>

using namespace ana;

void SetLoaderPaths2017(Loaders& loaders, const bool usedecafs = false)
{
  std::string ldrNDData  = "prod_sumdecaf_R17-03-01-prod3reco.d_nd_numi_fhc_full_v1_goodruns_nue2017";
  std::string ldrNDMC    = "prod_sumdecaf_R17-03-01-prod3reco.d_nd_genie_nonswap_fhc_nova_v08_full_v1_nue2017";
  std::string ldrNonSwap = "prod_sumdecaf_R17-03-01-prod3reco.g_fd_genie_nonswap_fhc_nova_v08_full_v1_nue2017";
  std::string ldrFlxSwap = "prod_sumdecaf_R17-03-01-prod3reco.g_fd_genie_fluxswap_fhc_nova_v08_full_v1_nue2017";
  std::string ldrTauSwap = "prod_sumdecaf_R17-03-01-prod3reco.g_fd_genie_tau_fhc_nova_v08_full_v1_nue2017";
  
  if (usedecafs) {
    ldrNDData  = "prod_decaf_R17-03-01-prod3reco.d_nd_numi_fhc_full_nue_or_numu_or_nus_contain_v1_goodruns";
    ldrNDMC    = "prod_decaf_R17-03-01-prod3reco.d_nd_genie_nonswap_fhc_nova_v08_full_nue_or_numu_or_nus_contain_v1";
    ldrNonSwap = "prod_decaf_R17-03-01-prod3reco.g_fd_genie_nonswap_fhc_nova_v08_full_nue_or_numu_or_nus_contain_v1";
    ldrFlxSwap = "prod_decaf_R17-03-01-prod3reco.g_fd_genie_fluxswap_fhc_nova_v08_full_nue_or_numu_or_nus_contain_v1";
    ldrTauSwap = "prod_decaf_R17-03-01-prod3reco.g_fd_genie_tau_fhc_nova_v08_full_nue_or_numu_or_nus_contain_v1";
  }
  
  
  loaders.SetLoaderPath(ldrNDData, caf::kNEARDET, Loaders::kData, 
                        DataSource::kBeam);
  loaders.SetLoaderPath(ldrNDMC,   caf::kNEARDET, Loaders::kMC, 
                        DataSource::kBeam);
  loaders.SetLoaderPath(ldrNonSwap, caf::kFARDET, Loaders::kMC, 
                        DataSource::kBeam, Loaders::kNonSwap);
  loaders.SetLoaderPath(ldrFlxSwap, caf::kFARDET, Loaders::kMC, 
                        DataSource::kBeam, Loaders::kFluxSwap);
  loaders.SetLoaderPath(ldrTauSwap, caf::kFARDET, Loaders::kMC, 
                        DataSource::kBeam, Loaders::kTauSwap);
  
}


void xsec_extrap_plots(bool createFile = false, 
                       TString suffix = "test")
{
  TString filename = "pred_xsec_fhc_" + suffix + ".root";
  
  struct WeiDef {
    TString name;
    Var weight;
  }; 
  struct ShiDef {
    TString name;
    SystShifts shift;
  };

  std::vector <WeiDef> weights;
  weights.push_back({"ppfx_xsec", kPPFXFluxCVWgt * kXSecCVWgt2017});
  weights.push_back({"ppfx_only", kPPFXFluxCVWgt});
  
  std::vector <ShiDef> shifts;
  shifts.push_back({"no_shift", kNoShift });
  
  for (auto sigma:{+1, -1}) {
    TString sig_str = (sigma == 1 ? "plus" : "minus");

    shifts.push_back({"MEC_enushape_" + sig_str, 
                     SystShifts(&kMECEnuShapeSyst2017, sigma)});
    shifts.push_back({"MEC_q0shape_" + sig_str, 
                     SystShifts(&kMECq0ShapeSyst2017,  sigma)});
    shifts.push_back({"MEC_NPfrac_"  + sig_str, 
                     SystShifts(&kMECInitStateNPFracSyst2017, sigma)});
    shifts.push_back({"MEC_scaleSA_" + sig_str, 
                     SystShifts(&kMECScaleSystSA, sigma)});
    shifts.push_back({"MEC_all3_" + sig_str,  
                     SystShifts({   
                       { &kMECEnuShapeSyst2017,        sigma},
                       { &kMECq0ShapeSyst2017 ,        sigma},
                       { &kMECInitStateNPFracSyst2017, sigma}
                       })
                     });
  } //end sigma

  if (createFile) {
      
    Loaders loaders;
    bool use_decafs = (suffix.Contains("test") || 
                       suffix.Contains("decaf"));
    SetLoaderPaths2017(loaders, use_decafs);
    loaders.SetSpillCut(kStandardSpillCuts);
   
    const HistAxis axisNumu = kNumuNonQEAxisFirstAna; 
    const HistAxis axisNue  = kNueAxisSecondAnaCVN;
   
    std::vector <IPrediction*> v_pred_prop;
    std::vector <IPrediction*> v_pred_nxp;
    std::vector <TString> names;

    // Loop over weights and shifts
    // Don't do all variations: e.g. only nominal weight for shifts   

    for (auto weight:weights) {  
      
      auto generator_nxp = NoExtrapGenerator(axisNue, 
                                             kNueSAFDCVNSsb, 
                                             weight.weight);
      
      for (auto shift:shifts) {

        if (shift.name == "no_shift" || weight.name == "ppfx_xsec") {  

          auto generator_prop = NuePropExtrapGenerator(axisNue, 
                                                       axisNumu,
                                                       kNueSAFDCVNSsb, 
                                                       kNueNDCVNSsb, 
                                                       kNumuND,
                                                       kNoShift, 
                                                       weight.weight);
          
          auto pred_prop = generator_prop.Generate(loaders, 
                                                   shift.shift);
          auto pred_nxp  = generator_nxp .Generate(loaders, 
                                                   shift.shift);
          
          v_pred_prop.emplace_back(pred_prop.release());
          v_pred_nxp. emplace_back(pred_nxp. release());
          
          names.push_back (shift.name + "_" + weight.name);
        }//end if 
      }//end shifts
    }//end weights
    
    loaders.Go();
    
    TFile* file = new TFile(filename, "recreate");
    
    for (int i = 0 ; i < (int) names.size() ; ++i) {
      v_pred_nxp [i]->SaveTo(file, "pred_nxp_"     + names[i]);
      v_pred_prop[i]->SaveTo(file, "pred_xp_prop_" + names[i]);
    }

    delete file;
    std::cout << "Saved predictions to " << filename << std::endl;
  
  }  //end create file
  else {
    // Make some plots
    TFile* file = new TFile (filename, "READ");
    TFile* outfile_xpplots = new TFile("plots/xp/modExtrapNuePlots.root",
                                       "RECREATE") ;
    bool printtable_numu = 1;
    bool printtable_nue = 1;
    bool pidaxes = true;
    bool printtable = 1;

    // Weights are on/off. Compare 1 to nominal

    TDirectory* dir_pred_no, * dir_pred_sh;
    TDirectory* dir_decomp_no, * dir_decomp_sh;
    TString dirname_decomp, plot_title, plot_name;

    for (auto weight:weights) {
      if (weight.name != "ppfx_xsec")  // Skip nominal vs nominal
        for (TString xp: {"nxp", "xp_prop"}) {
          
          dir_pred_no = file->GetDirectory ("pred_" + xp +
                                            "_no_shift_ppfx_xsec");
          dir_pred_sh = file->GetDirectory ("pred_" + xp + 
                                            "_no_shift_" +
                                            weight.name); 
          //ND plots
          if (xp == "xp_prop") {
            //Numu
            dirname_decomp = "extrap/MEextrap/Decomp";
            dir_decomp_no = dir_pred_no->GetDirectory (dirname_decomp);
            dir_decomp_sh = dir_pred_sh->GetDirectory (dirname_decomp);

            plot_title = weight.name + " (NearDet MC);" +
              " Non-Quasielastic Energy Estimator (GeV); Events / Bin";
            plot_name = "compare_numu_ND_" + weight.name;
            
            CompareNDDataMC (dir_decomp_no, dir_decomp_sh, 
                             plot_title, plot_name, "Numu",
                             printtable_numu);

            //nue
            dirname_decomp = "extrap/EEextrap/Decomp/";
            plot_title = weight.name + " (NearDet MC)"+
              "; Neutrino Energy Estimator (GeV) X CVN; Events / Bin";
            plot_name = "compare_nue_ND_" + weight.name;
            
            dir_decomp_no = dir_pred_no->GetDirectory (dirname_decomp);
            dir_decomp_sh = dir_pred_sh->GetDirectory (dirname_decomp);
            CompareNDDataMC (dir_decomp_no, dir_decomp_sh,
                             plot_title, plot_name, "CVN",
                             printtable_nue, pidaxes); 
          }  //end nd plots
          
          //FD Plots
          plot_title = weight.name +
            (xp != "nxp" ? " FarDet Extrap":" (FarDet MC)") +
            ";Neutrino Energy Estimator (GeV) X CVN;"+ "Events / Bin";
          plot_name = "compare_" + xp + "_" + weight.name;
          
          auto pred_no = ana::LoadFrom<IPrediction> (dir_pred_no);
          auto pred_sh = ana::LoadFrom<IPrediction> (dir_pred_sh);
          ComparePredictions (pred_no.get(), pred_sh.get(), 
                              plot_title, plot_name, "CVN",
                              printtable, pidaxes);
          
          //Additional XP plots
          if (xp == "xp_prop") {
            auto xp = (ModularExtrap*)
              ((PredictionExtrap*)pred_sh.release())->GetExtrap();
            xp->SavePlotsNue (outfile_xpplots->mkdir (
                                "plots_xp_" + weight.name), 9E20);
            
            if ("weight.name" == "no_weight") {
              xp = (ModularExtrap*) 
                ((PredictionExtrap*)pred_no.release())->GetExtrap();
              xp->SavePlotsNue (outfile_xpplots->mkdir (
                                  "plots_xp_ppfx_xsec"), 9E20);
            }  //end no weight
          }  //end XP plots
          
        }  //end xp
    }  //end weights 

    TDirectory* dir_pred_pl,   * dir_pred_mi;
    TDirectory* dir_decomp_pl, * dir_decomp_mi;
   
    // Compare plus and minus shifts to nominal at once
    for (TString shiftname: {"MEC_enushape", "MEC_q0shape",
                             "MEC_NPfrac",   "MEC_scaleSA",
                             "MEC_all3"}) {

      for (TString xp: {"nxp","xp_prop"}) {
        
        auto dir_pred_no = file->GetDirectory ("pred_" + xp + 
                                               "_no_shift_ppfx_xsec");
        auto dir_pred_pl = file->GetDirectory ("pred_" + xp + 
                                               "_" + shiftname + 
                                               "_plus_ppfx_xsec"); 
        auto dir_pred_mi = file->GetDirectory ("pred_" + xp +
                                               "_" + shiftname + 
                                               "_minus_ppfx_xsec"); 
        
        // ND plots
        if (xp == "xp_prop") {
          dirname_decomp = "extrap/MEextrap/Decomp/";
          plot_title = shiftname + " (NearDet MC)" +
            "; Non-Quasielastic Energy Estimator (GeV); Events / Bin";
          plot_name = "compare_numu_ND_" + shiftname;

          dir_decomp_no = dir_pred_no->GetDirectory (dirname_decomp);
          dir_decomp_pl = dir_pred_pl->GetDirectory (dirname_decomp);
          dir_decomp_mi = dir_pred_mi->GetDirectory (dirname_decomp);

          CompareNDDataMC (dir_decomp_no, dir_decomp_pl, dir_decomp_mi,
                           plot_title, plot_name, "Numu",
                           printtable_numu);


          dirname_decomp = "extrap/EEextrap/Decomp/";
          plot_title = shiftname + " (NearDet MC)" +
            "; Neutrino Energy Estimator (GeV) X CVN; Events / Bin";
          plot_name = "compare_nue_ND_" + shiftname;

          dir_decomp_no = dir_pred_no->GetDirectory (dirname_decomp);
          dir_decomp_pl = dir_pred_pl->GetDirectory (dirname_decomp);
          dir_decomp_mi = dir_pred_mi->GetDirectory (dirname_decomp);

          CompareNDDataMC (dir_decomp_no, dir_decomp_pl, dir_decomp_mi,
                           plot_title, plot_name, "CVN",
                           printtable_nue, pidaxes); 
        }
        
        // FD Plots
        plot_title = shiftname + 
          (xp != "nxp" ? " FarDet Extrap" : " (FarDet MC)") + 
          ";Neutrino Energy Estimator (GeV) X CVN;"+ "Events / Bin";
        plot_name = "compare_" + xp + "_" + shiftname;

        auto pred_no = ana::LoadFrom <IPrediction> (dir_pred_no);
        auto pred_pl = ana::LoadFrom <IPrediction> (dir_pred_pl);
        auto pred_mi = ana::LoadFrom <IPrediction> (dir_pred_mi);
        ComparePredictions (pred_no.get(), pred_pl.get(), pred_mi.get(),
                            plot_title, plot_name, "CVN",
                            printtable, pidaxes);
        
      }  // end xp
    }  // end weights
  }  // end plots
}