GeniePCASyst.cxx

Go to the documentation of this file.

#include "3FlavorAna/Systs/GeniePCASyst.h"

#include "CAFAna/Core/Utilities.h"

// 2020 cuts
#include "3FlavorAna/Cuts/NumuCuts2020.h"
#include "3FlavorAna/Cuts/NueCuts2020.h"
// 2018 cuts
#include "3FlavorAna/Cuts/NumuCuts2018.h"
#include "3FlavorAna/Cuts/NueCuts2018.h"

#include "CAFAna/Cuts/TruthCuts.h"

#include "StandardRecord/Proxy/SRProxy.h"

#include "TFile.h"
#include "TH1.h"
#include "TKey.h"

namespace ana
{
  // 2020 versions
  /*
  GeniePCASyst* GetGeniePrincipals2020(int PCIdx)
  {
    if(PCIdx >= 100){
      std::cerr << "Error : Component Index too large!" << std::endl;
      abort();
    }

    static GeniePCASyst* Cache[100] = {0};
    std::string PCIdxStr = TString::Format("%02d", PCIdx).Data();
    if(!Cache[PCIdx])
      Cache[PCIdx] = new GeniePCASyst((FindCAFAnaDir()+"/data/xs/genie_all_pc_shifts_fn_2020.root"),
                                      "genie_all_pc"+PCIdxStr, "Genie PC "+std::to_string(PCIdx));

    return Cache[PCIdx];
  }
  */

  GeniePCASyst* GetGeniePrincipals2020Small(int PCIdx)
  {
    if(PCIdx >= 100){
      std::cerr << "Error : Component Index too large!" << std::endl;
      abort();
    }

    static GeniePCASyst* Cache[100] = {0};
    std::string PCIdxStr = TString::Format("%02d", PCIdx).Data();
    if(!Cache[PCIdx])
      Cache[PCIdx] = new GeniePCASyst((FindCAFAnaDir()+"/data/xs/genie_small_pc_shifts_fn_2020.root"),
                                      "genie_small_pc"+PCIdxStr, "Genie PC "+std::to_string(PCIdx));

    return Cache[PCIdx];
  }

  // 2018 versions
  GeniePCASyst* GetGeniePrincipals2018(int PCIdx)
  {
    if(PCIdx >= 100){
      std::cerr << "Error : Component Index too large!" << std::endl;
      abort();
    }

    static GeniePCASyst* Cache[100] = {0};
    std::string PCIdxStr = TString::Format("%02d", PCIdx).Data();
    if(!Cache[PCIdx])
      Cache[PCIdx] = new GeniePCASyst((FindCAFAnaDir()+"/data/xs/genie_all_pc_shifts_fn_2018.root"),
                                      "genie_all_pc"+PCIdxStr, "Genie PC "+std::to_string(PCIdx));

    return Cache[PCIdx];
  }

  GeniePCASyst* GetGeniePrincipals2018Small(int PCIdx)
  {
    if(PCIdx >= 100){
      std::cerr << "Error : Component Index too large!" << std::endl;
      abort();
    }

    static GeniePCASyst* Cache[100] = {0};
    std::string PCIdxStr = TString::Format("%02d", PCIdx).Data();
    if(!Cache[PCIdx])
      Cache[PCIdx] = new GeniePCASyst((FindCAFAnaDir()+"/data/xs/genie_small_pc_shifts_fn_2018.root"),
                                      "genie_small_pc"+PCIdxStr, "Genie PC "+std::to_string(PCIdx));

    return Cache[PCIdx];
  }

  GeniePCASyst::GeniePCASyst(const std::string& fname,
                     const std::string& shortname,
                     const std::string& latexname)
    : ISyst(shortname,latexname),
      fFileName(fname), fHistName(shortname), fHistos()
  {
  }

  GeniePCASyst::~GeniePCASyst()
  {
    for(int i = 0; i < kNumFluxType; ++i){
      for(int j = 0; j < kNumDets; ++j){
        for(int k = 0; k < kNumSels; ++k){
          for(int l = 0; l < kNumFlavors; ++l){
            for(int m = 0; m < kNumSigns; ++m){
              delete fHistos[i][j][k][l][m];
              fHistos[i][j][k][l][m] = 0;
            }
          }
        }
      }
    }
  }

  void GeniePCASyst::InitializeHistograms() const
  {
    if(fHistos[0][0][0][0][0]) return;

    TFile fin(fFileName.c_str(),"READ");
    if(fin.IsZombie()){
      std::cerr << "Warning: couldn't open " << fFileName << std::endl;
      return;
    }

    TIter iterHist(gDirectory->GetListOfKeys());
    TKey* keyHist;
    int foundHisto=0;

    while((keyHist = (TKey*)iterHist())) {
      TString histName = keyHist->GetName(); // Get a histogram name
      if(histName.Contains(fHistName)){   //is it the desired syst
        foundHisto++;

        int flux=0;
        if(histName.Contains("FHC")) flux = kFHC;
        if(histName.Contains("RHC")) flux = kRHC;

        int det = 0;
        if(histName.Contains("ND")) det = kND;
        if(histName.Contains("FD")) det = kFD;

        int sel = 0;
        if(histName.Contains("Numu")) sel = kNumuSel;
        if(histName.Contains("Nue" )) sel = kNueSel;

        int flav = 0;
        if(histName.Contains("SigQE"   )) flav = kSigQE;
        if(histName.Contains("SigNonQE")) flav = kSigNonQE;
        if(histName.Contains("Other"   )) flav = kBkgFlav;
        if(histName.Contains("NC"      )) flav = kNC;

        int sign = 0;
        if(histName.Contains("minus")) sign = kMinus;
        if(histName.Contains("plus" )) sign = kPlus;

        //store relevant histograms
        fHistos[flux][det][sel][flav][sign] = (TH1D*) fin.Get(histName)->Clone();
        // disassociate it from the file it came from
        fHistos[flux][det][sel][flav][sign]->SetDirectory(0);
      }
    }
    
    if (foundHisto==0) {
      std::cerr << "Genie PC systematic histogram " << fHistName
                << " is not in file " << fFileName <<"; abort!" << std::endl;
      abort();
    }

    fin.Close();
  }

  double GeniePCASyst::GetWeight(const caf::SRProxy* sr,
                             ESign sign) const
  {
    InitializeHistograms();

    if(sr->mc.nnu == 0) return 1;

    int flux = sr->spill.isRHC ? kRHC : kFHC;

    int det = 0;
    switch (sr->hdr.det){
    case 1 : det = kND; break;
    case 2 : det = kFD; break;
    default: std::cerr <<"Wrong detector; ignore" << std::endl; return 1;
    };

    // 2018 remnants
    //const Cut NumuNDSel = kNumuCutND2018;
    //const Cut NueNDSel = kNue2018NDCVNSsb;
    //const Cut NumuFDSel = kNumuCutFD2018;
    //const Cut NueFDSel = kNue2018FDAllSamples;
    
    // 2020 version
    const Cut NumuNDSel = kNumu2020ND;
    const Cut NueNDSel = kNue2020ND;
    const Cut NumuFDSel = kNumu2020FD;
    const Cut NueFDSel = kNue2020FD;

    int sel = -1;
    if(det == kND){
      if(NumuNDSel(sr)) sel = kNumuSel;
      if(NueNDSel(sr))  sel = kNueSel;
    }
    if(det == kFD){
      if(NumuFDSel(sr)) sel = kNumuSel;
      if(NueFDSel(sr))  sel = kNueSel;
    }
    if(sel == -1) return 1;

    int flav = -1;
    if(sel == kNumuSel){
      if(kIsNumu(sr) && !kIsNC(sr) && kIsQE(sr)) flav = kSigQE;
      if(kIsNumu(sr) && !kIsNC(sr) && !kIsQE(sr)) flav = kSigNonQE;
      if(kIsNC(sr)) flav = kNC;
    }
    if(sel == kNueSel){
      if(kIsNue(sr) && !kIsNC(sr) && kIsQE(sr)) flav = kSigQE;
      if(kIsNue(sr) && !kIsNC(sr) && !kIsQE(sr)) flav = kSigNonQE;
      if(kIsNumu(sr) && !kIsNC(sr)) flav = kBkgFlav;
      if(kIsNC(sr)) flav = kNC;
    }
    if(flav == -1) return 1;


    TH1D* h = fHistos[flux][det][sel][flav][sign];
    if (h == 0){
      std::cerr << fHistName+": Can't find desired histogram; ignore"<< std:: endl;
      return 1;
    }

    double energy = sr->mc.nu[0].E;

    if (energy > h->GetXaxis()->GetXmin() &&
        energy < h->GetXaxis()->GetXmax() ){
      return h->Interpolate(energy);
    }

    if (energy > h->GetXaxis()->GetXmax())
      return h->GetBinContent(h->GetNbinsX());

    return 1;
  }

  void GeniePCASyst::Shift(double sigma,
                           caf::SRProxy* sr, double& weight) const
  {
    if(sigma == 0) return;

    // Need different histogram for plus or minus sigma
    ESign sign = kPlus;
    if(sigma<0) sign = kMinus;

    weight *=1+(GetWeight(sr, sign)-1.)*TMath::Abs(sigma);
    weight = std::max(0., weight);
  }

  TH1D* GeniePCASyst::getHist(int f,int d,int s,int fl, int p)
  {
    InitializeHistograms();
    
    return fHistos[f][d][s][fl][p];
  }
}