ND_DataMC_Comp_Systs.C

Go to the documentation of this file.

// Fills spectra for systematic error band.  Makes both pot norm with full syst band and area norm with shape only syst band
// have to manually make nd_syst_plots directory before running

bool IsFHC=false; // Global variable for whether I am running over FHC or RHC?
bool Bless=true;  // Global variable for whether I am running blessing plots?

#include "3FlavorAna/Ana2018/numu/ND_DataMC_Comp_Systs_header.h"

#include "CAFAna/Cuts/Cuts.h"
#include "3FlavorAna/Cuts/NumuCuts.h"
#include "3FlavorAna/Cuts/NumuCuts2017.h"
#include "3FlavorAna/Cuts/NumuCuts2018.h"
#include "3FlavorAna/Cuts/QuantileCuts.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "CAFAna/Core/Var.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Systs/Systs.h"
#include "3FlavorAna/Systs/NumuSysts.h"
#include "CAFAna/Systs/XSecSystLists.h"
#include "CAFAna/Vars/GenieWeights.h"
#include "3FlavorAna/Vars/HistAxes.h"
#include "3FlavorAna/Vars/NumuVars.h"
#include "CAFAna/Vars/PPFXWeights.h"

#include "TFile.h"
#include "TSystem.h"
#include "3FlavorAna/Vars/NumuVars.h"

using namespace ana;

////////////////////////////////////////
// Making a text file for blessing:
void MakeTextFile( std::string OutName, std::string OutDir ) {
  // Determine output name, and caption.
  std::string FNa = OutDir + OutName + ".txt";

  // Replace underscores in the canvas name
  while (OutName.find("_") != std::string::npos) OutName.replace(OutName.find("_"),1," ");
  // Add a string about it being a cut level.
  if (OutName.find("Qual") != std::string::npos) OutName.insert(OutName.find(" Qual"), ", Cut level:");
  // Add a string about it being a cut level.
  if (OutName.find("pot") != std::string::npos) OutName.replace(OutName.find(" pot"),4,". With POT normalisation.");
  // Add a string about it being a cut level.
  if (OutName.find("area") != std::string::npos) OutName.replace(OutName.find(" area"),5,". With Area normalisation.");
  // Add a string about the Horn current.
  if (OutName.find("fhc")  != std::string::npos) OutName.replace(OutName.find(" fhc"),4,", in Neutrino mode");
  if (OutName.find("rhc")  != std::string::npos) OutName.replace(OutName.find(" rhc"),4,", in Antineutrino mode");

  std::string Cap = "Plot showing the number of ND NuMi events in #nu_{#mu} Ana18 for variable: " + OutName;
  std::cout << "\nFile name: " << FNa << "\n\tCaption: " << Cap << std::endl;
  // Write to file.
  std::ofstream TxtOut ( FNa.c_str(), std::ofstream::out );
  TxtOut << Cap;
  TxtOut.close();
  // Done.
  return;
}

////////////////////////////////////////
// Main class function:
void ND_DataMC_Comp_Systs( bool fhc = true, bool Set1Vars = true, bool IsTest=false, bool IsOnGrid=true, bool Weight18 = true, bool Prod4Cut = true, bool FullCuts = false ) {
  
  // Set my global Horn Current variable...
  IsFHC = fhc;

  std::string sFHC_RHC = (IsFHC    == true ? "fhc"     :"rhc"     );
  std::string sVarSet  = (Set1Vars == true ? "Set1"    : "Set2"   );
  std::string sWeight  = (Weight18 == true ? "wgt2018" :"wgt2017" );
  std::string sProd4   = (Prod4Cut == true ? "CutProd4":"CutProd3");
  std::string sAllCut  = (FullCuts == true ? "_AllCuts": ""       );

  std::cout << "\n\nLots of input parameters, so I'm going to write them out now....;"
            << "\n\t IsFHC      " << (IsFHC     ==true ? "true  --> Run in FHC "          :"false --> Run in RHC"             ) << " ==> sFHC_RHC is " << sFHC_RHC
            << "\n\t Set1Vars   " << (Set1Vars  ==true ? "true  --> Making first 9 plots" :"false --> Making second 9 plots"  ) << " ==> sVarSet is " << sVarSet
            << "\n\t IsTest     " << (IsTest    ==true ? "true  --> Short run use CAFs"   :"false --> Long run use Concat"    ) 
            << "\n\t Weight18   " << (Weight18  ==true ? "true  --> Run using 2018 XSecs" :"false --> Run using 2017 XSecs"   ) << " ==> sWeight is "  << sWeight
            << "\n\t Prod4Cut   " << (Prod4Cut  ==true ? "true  --> Using Prod4 Cuts"     :"false --> Using Prod3 Cuts"       ) << " ==> sProd4 is "   << sProd4
            << "\n\t FullCuts   " << (FullCuts  ==true ? "true  --> Running all cuts"     :"false --> Only running full cuts" ) << " ==> sAllCut is "  << sAllCut
            << "\n\n" << std::endl;

  ////////////////////////////////////////
  // Get stuff for the hadEFrac Cut:
  ////////////////////////////////////////
  // Which file I want to load?
  // Diana's original files
  std::string fdspecfile = "/pnfs/nova/persistent/analysis/numu/Ana2018/official/Quantiles/quantiles__"+sFHC_RHC+"_full__numu2018.root";
    
  std::cerr << "\nI want to load file " << fdspecfile << "\n" << std::endl;

  TFile* inFile = TFile::Open( pnfs2xrootd(fdspecfile).c_str() );
  TH2 *FDSpec2D = (TH2*)inFile->FindObjectAny( "FDSpec2D" ); //TString("FD_full_") + TString(sFHC_RHC) );
  // How many quantiles?
  const int NHadEFracQuantiles = 4; // defines how many divisions of hadEFrac are used 
  // Make my cuts.
  //std::vector<Cut> HadEFracQuantCuts = QuantileAndPIDCutsFromTH2(FDSpec2D, kNumuCCOptimisedAxis, kHadEFracAxis, NHadEFracQuantiles, IsFHC);
  std::vector<Cut> HadEFracQuantCuts = QuantileCutsFromTH2(FDSpec2D, kNumuCCOptimisedAxis, kHadEFracAxis, NHadEFracQuantiles);
  ////////////////////////////////////////
  // Define Axes and cuts.
  ////////////////////////////////////////

  int NumBins = 25;

  const HistAxis EnergyAxis  ( "Reconstructed Muon Energy (GeV)"       , Binning::Simple( NumBins , 0   , 5   ) , kMuE        );
  const HistAxis MuonEPerHit ( "Muon energy per hit in track"          , Binning::Simple( NumBins , 0.01, 0.03) , kTrkEPerNHit);
  const HistAxis LengthAxis  ( "Length of Primary Track (m)"           , Binning::Simple( NumBins , 0   , 16  ) , kTrkLength  );
  const HistAxis HadEAxis    ( "Hadronic Energy (GeV)"                 , Binning::Simple( NumBins , 0   , 3   ) , kHadE       );
  const HistAxis HadEFracAxis( "Hadronic Energy Fraction"              , Binning::Simple( NumBins , 0   , 1.01) , kHadEFrac   );
  const HistAxis HadEFracHRes( "Hadronic Energy Fraction"              , Binning::Simple( 250     , 0   , 1.01) , kHadEFrac   );
  const HistAxis HadEPerHit  ( "Hadronic energy per hit in slice"      , Binning::Simple( NumBins , 0   , 0.04) , kHadEPerNHit);
  const HistAxis NHitsInSlc  ( "Number of hits in the slice"           , Binning::Simple( NumBins , 0   , 300 ) , kNHit       );
  const HistAxis NHitsKalTr  ( "Number of hits in Kalman track"        , Binning::Simple( NumBins , 0   , 300 ) , kTrkNhits   );
  const HistAxis CosNuMiAxis ( "Kalman Track Cos #theta_{NuMI}"        , Binning::Simple( NumBins , 0   , 1.01) , kCosNumi    );
  const HistAxis pTpAxis     ( "p_{T}/p"                               , Binning::Simple( NumBins , 0   , 1.01) , kNumuMuonPtP);
  const HistAxis CVN2017Axis ( "CVN NuMu 2017 score"                   , Binning::Simple( NumBins , 0   , 1.01) , kCVNm2017   );
  const HistAxis CVN2018Axis ( "CVN NuMu score"                        , Binning::Simple( NumBins , 0   , 1.01) , kCVNm       );
  const HistAxis ReMIDAxis   ( "ReMId score"                           , kRemidBinning                          , kRemID      );
  const HistAxis ReMIDFrac   ( "RemID fraction of planes"              , kRemidBinning                          , kReMIdMeasFrac);
  const HistAxis ReMIDScat   ( "ReMId Scattering angle"                , Binning::Simple( NumBins ,-0.25, 0.25) , kReMIdScatLLH );
  const HistAxis ReMIDdEdx   ( "ReMId dE/dx"                           , Binning::Simple( NumBins ,-0.5 , 0.5 ) ,kReMIdDEDxLLH );
  
  std::vector<Cut> TieredCuts;
  std::vector<std::string> CutNames;
  if (Prod4Cut) {
    CutNames.emplace_back("Qual_Cont_PID"); TieredCuts.emplace_back(kNumuQuality && kNumuContainND2017 && kNumuPID2018);
  } else {
    CutNames.emplace_back("Qual_Cont_PID"); TieredCuts.emplace_back(kNumuQuality && kNumuContainND2017 && kNumuPID2017);
  }
  if (FullCuts) {
    CutNames.emplace_back("Qual_Cont"); TieredCuts.emplace_back(kNumuQuality && kNumuContainND2017);
    CutNames.emplace_back("Qual");      TieredCuts.emplace_back(kNumuQuality);
    CutNames.emplace_back("NoCut");     TieredCuts.emplace_back(kNoCut);
  }
  ////////////////////////////////////////
  // Now make the plots....
  ////////////////////////////////////////
  gStyle->SetMarkerStyle(kFullCircle);
  TGaxis::SetMaxDigits(3);

  // --- Configure the definition to use.
  // By default set loaders empty so that it breaks.
  std::string fnameNDMC   = ""; // Dataset name for MC
  std::string fnameNDData = ""; // Dataset name for Data
  // If test.
  if (IsTest) {
    if (IsFHC) {
      fnameNDMC   = "karlwarb_ND_FHC_MC";  //"prod_caf_R17-11-14-prod4reco.e_nd_genie_nonswap_rhc_nova_v08_full_v1";
      fnameNDData = "karlwarb_ND_FHC_Data";//"prod_caf_R17-09-05-prod4recopreview.f_nd_numi_rhc_full_v1_addShortSimpleCVN_goodruns";
    } else {
      fnameNDMC   = "karlwarb_ND_RHC_MC";
      fnameNDData = "karlwarb_ND_RHC_Data";
    }
  } else {
    // If for running for realsies.
    if (IsFHC) {
      fnameNDMC   = "prod_sumdecaf_R17-11-14-prod4reco.d_nd_genie_nonswap_fhc_nova_v08_full_v1_numu2018_stride22";
      //fnameNDMC   = "prod_sumdecaf_R17-11-14-prod4reco.d_nd_genie_nonswap_fhc_nova_v08_full_v1_numu2018_lim1";
      fnameNDData = "prod_sumdecaf_R17-09-05-prod4recopreview.f_nd_numi_fhc_full_v1_addShortSimpleCVN_goodruns_numu2018";
    } else {
      fnameNDMC   = "prod_sumdecaf_R17-11-14-prod4reco.e_nd_genie_nonswap_rhc_nova_v08_full_v1_numu2018_stride14";
      //fnameNDMC   = "prod_sumdecaf_R17-11-14-prod4reco.e_nd_genie_nonswap_rhc_nova_v08_full_v1_numu2018_lim1";
      fnameNDData = "prod_sumdecaf_R17-09-05-prod4recopreview.f_nd_numi_rhc_full_v1_addShortSimpleCVN_goodruns_numu2018";
    }
  }
  // --- Set the output directory.
  std::string fnameODir   = "NDComp_"+sFHC_RHC+"_"+sVarSet+"_"+sWeight+"_"+sProd4+sAllCut; // Name of the output directory
  if (IsTest) fnameODir += "_Test";
  
  std::cout << "\n\n ODir is " << fnameODir << ", I loaded cuts from " << fdspecfile
            << "\n\t fnameNDMC   is " << fnameNDMC
            << "\n\t fnameNDData is " << fnameNDData
            << "\n\n" << std::endl;

  // --- Set the Loaders, along with spill cuts.
  SpectrumLoader loaderMC(fnameNDMC);
  SpectrumLoader loaderData(fnameNDData);
  loaderMC.SetSpillCut(kStandardSpillCuts);
  loaderData.SetSpillCut(kStandardSpillCuts);
  
  // --- Make sure that my directory exists...
  fnameODir += "/";
  if (IsOnGrid) fnameODir = "";
  else          gSystem  -> MakeDirectory( fnameODir.c_str() );
  
  // --- Load the systematics
  std::vector<const ISyst*> systs = getAllDirectNumuSysts2018();

  std::cout << "\n\nsysts has size " << systs.size() << " compared to 2017 systs which has size " << getAllDirectNumuSysts2017().size() << "\n\n" << std::endl;

  // --- Make a vector of "Selections" which are a combination of Quantile && Tiered cuts.
  std::vector<Selection> selections;
  for (size_t cuts=0; cuts < TieredCuts.size(); ++cuts) {
    for (size_t quant=0; quant < HadEFracQuantCuts.size(); ++quant) {
      // --- Want no quantile cut selections too.
      if (quant == 0) {
        selections.emplace_back( TieredCuts[cuts], sFHC_RHC+"_"+CutNames[cuts]+"_QuantAll",
                                 "Selected events pass "+CutNames[cuts]+" all quantiles");
      }
      selections.emplace_back( TieredCuts[cuts] && HadEFracQuantCuts[quant], sFHC_RHC+"_"+CutNames[cuts]+"_Quant"+std::to_string(quant+1),
                               "Selected events pass "+CutNames[cuts]+" quantile "+std::to_string(quant+1)+".");
    }
  }
  
  // --- Make a vector of "TangibleAxis" which have the axis and variables for which I want to make spectra.
  std::cerr << "About to fill variables" << std::endl;
  std::vector<TangibleAxis> variables;
  variables.emplace_back( kNumuCCOptimisedAxis, "NuMuE", "CC energy estimator.");
  if (Set1Vars) {
    variables.emplace_back( EnergyAxis  , "MuonEnergy" , "Reconstructed Muon Energy (GeV)" );
    variables.emplace_back( MuonEPerHit , "MuEPerHit"  , "Muon energy per hit in track"    );
    variables.emplace_back( HadEAxis    , "HadroEnergy", "Hadronic Energy (GeV)"           );
    variables.emplace_back( HadEFracAxis, "HadEFrac"   , "Hadronic Energy Fraction"        );
    variables.emplace_back( HadEPerHit  , "HadEPerHit" , "Hadronic energy per hit in slice");
    variables.emplace_back( LengthAxis  , "KalTrackLen", "Length of Primary Track (m)"     );
    variables.emplace_back( NHitsInSlc  , "NHitsInSlc" , "Number of hits in the slice"     );
    variables.emplace_back( CosNuMiAxis , "TrkCosNuMi" , "Kalman Track Cos #theta_{NuMI}"  );
    variables.emplace_back( pTpAxis     , "pTpNuMi"    , "p_{T}/p"                         );
  } else {
    variables.emplace_back( HadEFracHRes, "HadEFacHRes", "Hadronic Energy Fraction"        );
    variables.emplace_back( NHitsKalTr  , "NHitsKalTr" , "Number of hits in Kalman track"  );
    variables.emplace_back( CosNuMiAxis , "TrkCosNuMi" , "Kalman Track Cos #theta_{NuMI}"  );
    variables.emplace_back( CVN2017Axis , "CVNNuMu2017", "CVN NuMu 2017 score"             );
    variables.emplace_back( CVN2018Axis , "CVNNuMu2018", "CVN NuMu score"                  );
    variables.emplace_back( ReMIDAxis   , "ReMIdScore" , "ReMId score"                     );
    variables.emplace_back( ReMIDFrac   , "ReMIdFrac"  , "RemID fraction of planes"        );
    variables.emplace_back( ReMIDScat   , "ReMIdScatt" , "ReMId Scattering angle"          );
    variables.emplace_back( ReMIDdEdx   , "ReMIddEdx"  , "ReMId dE/dx"                     );
  }
  std::cerr << "Variables has size " << variables.size() << std::endl;

  // --- What is my weight going to be for the spectra?
  //Var MyWeight            = kXSecCVWgt2018 * kPPFXFluxCVWgt;
  Var MyWeight            = kXSecCVWgt2018RPAFix_noDIS * kPPFXFluxCVWgt;
  if (!Weight18) MyWeight = kXSecCVWgt2017 * kPPFXFluxCVWgt;
  
  std::vector<DataMCPair> pairs;
  
  Cut      WrongSign = kIsNu;
  if (fhc) WrongSign = kIsAntiNu;

  pairs.reserve(selections.size() * variables.size());
  for(const auto& sel:selections){
    for(const auto& variable:variables){
      pairs.emplace_back(sel, variable, loaderData, loaderMC, systs, MyWeight, !kIsNumuCC, WrongSign );
    }
  }
  
  std::cerr << "Filling my things..." << std::endl;
  
  loaderMC.Go();
  loaderData.Go();
  
  std::vector<TCanvas*> cans;

  std::cerr << "Filled my things" << std::endl;
  
  
  for(const auto& pair:pairs) {
    
    std::string Name_POT  = pair.ShortName()+"_allSysts_pot";     
    std::string Name_Area = pair.ShortName()+"_allSysts_area";
    std::cout << "\nNow looking pairs, " << pair.ShortName() << " - POT " << Name_POT << ", Area " << Name_Area << ", it has purity " << pair.Purity() << std::endl;
    
    // --- First make the POT Normalised plots.
    cans.push_back( new TCanvas( Name_POT.c_str(),Name_POT.c_str() ) );
    pair.OverlayDataMCSyst();
    WriteBlurb(pair, fnameODir,
               std::string(cans.back()->GetName())+".txt",
               "Full systematic band shown.  ",
               "Near detector data/MC comparison for numu Ana2018." );
    
    // ---- Next make the Area Normalised plots.
    cans.push_back( new TCanvas( Name_Area.c_str(), Name_Area.c_str() ) );
    pair.OverlayDataMCSystNorm();
    WriteBlurb(pair, fnameODir,
               std::string(cans.back()->GetName())+".txt",
               std::string("Each systematically shifted histogram (both up and down) ")
               +std::string("is normalized to the area of the MC distribution, ")
               +std::string("then summed in quadrature." ),
               "Near detector data/MC comparison for numu Ana2018." );
    
  }
  
  // Open a file to save all of the canvases to.
  std::string OutFileName = fnameODir + "NearDet_DataMC_Comp_Systs.root";
  TFile *OutFile = new TFile( OutFileName.c_str(), "RECREATE" );
  OutFile->cd();
  
  for(const auto& can:cans) {
    can->cd();
    Preliminary();
    // -- Set the corner label
    std::string Str = "Neutrino beam";
    if (!IsFHC) Str = "Antineutrino beam";
    TLatex* CornLab = new TLatex(.15, .93, Str.c_str());
    CornLab->SetTextColor(kGray+1);
    CornLab->SetNDC();
    CornLab->SetTextSize (2/30.);
    CornLab->SetTextAlign(11);
    CornLab->Draw();
    //can -> SetLeftMargin(0.13);
    std::string CanNa = fnameODir + can->GetName() + ".pdf";
    can->SaveAs ( CanNa.c_str() );
    can->Write  ( can->GetName() );
    MakeTextFile( can->GetName(), fnameODir );
  }
}