demo_cut_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 "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"

using namespace ana;


void demo_cut_optimization(std::string input_file_name = "")
{
  if(input_file_name == "" ) {
    std::map<TString, SpectrumLoader *> loaders;    
    loaders["kNominal"]          = new SpectrumLoader("prod_caf_R17-11-14-prod4reco.remid-hotfix.b_nd_genie_nonswap_fhc_nova_v08_full_v1");
    loaders["kCalibNeg"]         = new SpectrumLoader("prod_caf_R17-11-14-prod4reco.remid-hotfix.b_nd_genie_nonswap_fhc_nova_v08_full_calib-shift-nd-xyview-neg-offset_v1_good_seventh");
    loaders["kCalibPos"]         = new SpectrumLoader("prod_caf_R17-11-14-prod4reco.remid-hotfix.b_nd_genie_nonswap_fhc_nova_v08_full_calib-shift-nd-xyview-pos-offset_v1_good_seventh");
    loaders["kCalibShape"]       = new SpectrumLoader("prod_caf_R17-11-14-prod4reco.remid-hotfix.b_nd_genie_nonswap_fhc_nova_v08_full_calib-shift-nd-func_v1_good_seventh");
    loaders["kLightUpCalibDown"] = new SpectrumLoader("prod_caf_R17-11-14-prod4reco.remid-hotfix.b_nd_genie_nonswap_fhc_nova_v08_full_lightmodel-lightup-calibdown_v1_good_seventh");
    loaders["kLightDownCalibUp"] = new SpectrumLoader("prod_caf_R17-11-14-prod4reco.remid-hotfix.b_nd_genie_nonswap_fhc_nova_v08_full_lightmodel-lightdown-calibup_v1_good_seventh");
    loaders["kCherenkov"]        = new SpectrumLoader("prod_caf_R17-11-14-prod4reco.remid-hotfix.b_nd_genie_nonswap_fhc_nova_v08_full_ckv-proton-shift-down_v1_good_seventh");

    // define a hist axis to optimize
    const HistAxis * opti_axis = new HistAxis("Vertex in X Direction (cm)",
                                              Binning::Simple(50, -200, 200),
                                              kVtxX);
    // Define selection and signal cuts
    const Cut * kSignal = new Cut(CutFromNuTruthCut(kIsNueCCST && kTrueDetectorST));
    // This demo is for fiducial cut optimization so the selection
    // should be all preselection minus fiducial constraints
    const Cut * kSelection = new Cut(kcDQ && 
                                     kcNueCCIncContainment && 
                                     kcFrontPlanes &&
                                     kcNHits);

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

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

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

    // Setup the genie multiverses   
    // We take flux uncertainty to be constant, so no need for PPFX universe
    // for the sake of example, just do 10 universes
    GenieMultiverseParameters genie_multiverse(10, getAllXsecSysts_2018());
    std::vector<SystShifts> genie_shifts = genie_multiverse.GetSystShifts();
    
    SystematicDef * sGenieMultiverse =
      new SystematicDef("sGenieMultiverse",
                        loaders["kNominal"],
                        opti_axis,
                        kCVWeight);
                                                                
    // create the optimizer. 
    CutOptimization * optimize = new CutOptimization(opti_axis,
                                                     kSignal,
                                                     kSelection);
    optimize->DefineNominal(sNominal);
    optimize->DefineSystematic(sLightlevelUpDown);
    optimize->DefineSystematic(sCalibPosNeg);
    optimize->DefineSystematic(sCalibShape);
    optimize->DefineSystematic(sCherenkov);
    optimize->DefineMultiverseSystematic(sGenieMultiverse, genie_shifts);

    // optionally set spill cuts for all loaders 
    optimize->SetSpillCuts(kStandardSpillCuts);

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

    // save to file
    TFile * output = new TFile("demo_cut_optimization.root", "recreate");
    optimize->SaveTo(output->mkdir("vtx_cut_optimize"));    
    output->Close();
  }
  else {
    // load input file provided in argument
    TFile * input = TFile::Open(input_file_name.c_str(), "read");
    CutOptimization * optimize = CutOptimization::LoadFrom(input->GetDirectory("vtx_cut_optimize")).release();
    input->Close();

    
    // save resulting histograms back to file.
    TFile * results = new TFile("demo_cut_optimization_results.root", "recreate");
    results->mkdir("vtx_cut_optimize/results");
    TDirectory * save_dir = results->GetDirectory("vtx_cut_optimize/results");

    // directory to save plots
    std::string plot_dump = "demo_cut_optimization_plots";

    optimize->Optimize(CutOptimization::kBinByBin,        // OptimizationMethod_t
                       CutOptimization::kdSigmaOverSigma, // FOM_t
                       kAna2018FHCPOT,                    // scale to data pot
                       save_dir,
                       plot_dump,
                       true);
    delete optimize;
  }
  
}