fiducial_optimization.C

Go to the documentation of this file.

#pragma once

#include "CAFAna/XSec/CutOptimization.h"

// analysis cuts/vars
#include "NDAna/nuebarcc_inc/NuebarCCIncCuts.h"
#include "NDAna/nuebarcc_inc/NuebarCCIncExtra.h"
#include "3FlavorAna/Vars/NueVars.h"
#include "CAFAna/Cuts/SpillCuts.h"

// for data pot
#include "CAFAna/Analysis/Exposures.h"

// Multiverse includes
#include "CAFAna/Systs/XSecSystLists.h"
#include "CAFAna/Vars/XsecTunes.h"
#include "CAFAna/Vars/PPFXWeights.h"
#include "CAFAna/XSec/GenieMultiverseSyst.h"
#include "CAFAna/XSec/FluxMultiverseSyst.h"

// temporary measure to use prod4 trained muonid
#include "NDAna/muonid/NDXSecMuonPID.h"
namespace nuebarccinc {
  const ana::Cut kMuonIDProd4Cut = ana::kMuonIDProd4 < -0.2;
}

void load_libs_muonid()
{
  const int ret = gSystem->Load("libNDAnamuonid");
  if(ret != 0) exit(1);
}

using namespace ana;
void fiducial_optimization(int axis = 0)
{
  load_libs_muonid();
  
  // axis = 0, 1, 2 for x, y, z respectively
  std::map<TString, SpectrumLoader *> loaders;    
  loaders["kNominal"]          = new SpectrumLoader(nuebarccinc::PROD5_MC_RHC_NOMINAL    );
  loaders["kCalibNeg"]         = new SpectrumLoader(nuebarccinc::PROD5_MC_RHC_CALIB_DOWN );
  loaders["kCalibPos"]         = new SpectrumLoader(nuebarccinc::PROD5_MC_RHC_CALIB_UP   );
  loaders["kCalibShape"]       = new SpectrumLoader(nuebarccinc::PROD5_MC_RHC_CALIB_SHAPE);
  loaders["kLightUpCalibDown"] = new SpectrumLoader(nuebarccinc::PROD5_MC_RHC_LIGHT_UP   );
  loaders["kLightDownCalibUp"] = new SpectrumLoader(nuebarccinc::PROD5_MC_RHC_LIGHT_DOWN );
  loaders["kCherenkov"]        = new SpectrumLoader(nuebarccinc::PROD5_MC_RHC_CHERENKOV  );


  std::vector<std::string> opti_labels = {"fiducial_vtx_x", "fiducial_vtx_y", "fiducial_vtx_z"};
  // define a hist axis to optimize
  std::vector<const HistAxis *> opti_axes = {new HistAxis("Vertex in X Direction (cm)",
                                                          Binning::Simple(40, -200, 200),
                                                          kVtxX),
                                             new HistAxis("Vertex in Y Direction (cm)",
                                                          Binning::Simple(40, -200, 200),
                                                          kVtxY),
                                             new HistAxis("Vertex in Z Direction (cm)",
                                                          Binning::Simple(80, 0, 1200),
                                                          kVtxZ)};
  // Define selection and signal cuts
  const Cut * kSignal = new Cut(CutFromNuTruthCut((nuebarccinc::kNuebarCCST || nuebarccinc::kNueCCST ) && 
                                                  nuebarccinc::kTrueDetectorST));
  // This demo is for fiducial cut optimization so the selection
  // should be all preselection minus fiducial constraints
  const Cut * kSelection = new Cut(nuebarccinc::kDQ                          && 
                                   nuebarccinc::kFrontPlanes                 &&
                                   nuebarccinc::kNHits                       &&
                                   nuebarccinc::kMuonIDProd4Cut              &&
                                   nuebarccinc::kNuebarCCIncContainmentLoose && 
                                   nuebarccinc::kNuebarCCIncFiducialLoose    &&
                                   nuebarccinc::kCVNNueIDCut);
  

  // central value weight to be applied to all samples
  const Var * kCVWeight = new Var(kPPFXFluxCVWgt * kXSecCVWgt2020);

  // nominal definition
  SystematicDef * sNominal =
    new SystematicDef("sNominal",
                      loaders["kNominal"],
                      opti_axes[axis],
                      kCVWeight);

  // systematic definitions
  SystematicDef * sLightlevelUpDown = 
    new SystematicDef("sLightlevelUpDown",
                      loaders["kLightUpCalibDown"],
                      loaders["kLightDownCalibUp"],
                      opti_axes[axis],
                      kCVWeight);
  SystematicDef * sCalibPosNeg = 
    new SystematicDef("sCalibPosNeg",
                      loaders["kCalibPos"],
                      loaders["kCalibNeg"],
                      opti_axes[axis],
                      kCVWeight);
  SystematicDef * sCalibShape = 
    new SystematicDef("sCalibShape",
                      loaders["kCalibShape"], 
                      opti_axes[axis],
                      kCVWeight);
  SystematicDef * sCherenkov = 
    new SystematicDef("sCherenkov",
                      loaders["kCherenkov"], 
                      opti_axes[axis],
                      kCVWeight);

  // Setup the genie multiverses   
  // We take flux uncertainty to be constant, so no need for PPFX universe
  GenieMultiverseParameters genie_multiverse(300, getAllXsecSysts_2020());
  std::vector<SystShifts> genie_shifts = genie_multiverse.GetSystShifts();
    
  SystematicDef * sGenieMultiverse = 
    new SystematicDef("sGenieMultiverse",
                      loaders["kNominal"],
                      opti_axes[axis],
                      kCVWeight);
                                                            
  // create the optimizer. 
  CutOptimization * optimizer;
  optimizer = new CutOptimization(opti_axes[axis],
                                  kSignal,
                                  kSelection);
  optimizer->DefineNominal(sNominal);
  optimizer->DefineSystematic(sLightlevelUpDown);
  optimizer->DefineSystematic(sCalibPosNeg     );
  optimizer->DefineSystematic(sCalibShape      );
  optimizer->DefineSystematic(sCherenkov       );
  optimizer->DefineMultiverseSystematic(sGenieMultiverse, genie_shifts);
  // optionally set spill cuts for all loaders 
  optimizer->SetSpillCuts(kStandardSpillCuts);


  // Let the loaders go
  optimizer->Go();

  // save to file
  TFile * output = new TFile(("fiducial_optimization_" + opti_labels[axis] + ".root").c_str(), "recreate");
  optimizer->SaveTo(output, opti_labels[axis].c_str());    
  output->Close();

}