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

using namespace ana;

void prongcvn_optimization()
{

  std::map<TString, SpectrumLoader *> loaders;    
  loaders["kNominal"]          = new SpectrumLoader(NominalMC_entire_RHC_prod4);
  loaders["kCalibNeg"]         = new SpectrumLoader(DWNCalibFiles_full_RHC_prod4);
  loaders["kCalibPos"]         = new SpectrumLoader(UPCalibFiles_full_RHC_prod4);
  loaders["kCalibShape"]       = new SpectrumLoader(ShapeCalibFiles_full_RHC_prod4);
  loaders["kLightUpCalibDown"] = new SpectrumLoader(UPLightFiles_full_RHC_prod4);
  loaders["kLightDownCalibUp"] = new SpectrumLoader(DWNLightFiles_full_RHC_prod4);
  loaders["kCherenkov"]        = new SpectrumLoader(CherenkovFiles_full_RHC_prod4);

  // define a hist axis to optimize
  const HistAxis * opti_axis = new HistAxis("Max ProngCVN Electron Score",
                                            Binning::Simple(50, 0, 1),
                                            kvProngCVN);
  // Define selection and signal cuts
  const Cut * kSignal = new Cut(CutFromNuTruthCut(kIsNuebarCCST && kTrueDetectorST));

  const Cut * kSelection = new Cut(kcDQ                    &&                              
                                   kcFrontPlanes           &&
                                   kcNHits                 &&
                                   kTrueDetector           &&
                                   kcNueCCIncContainment   &&
                                   kcNueCCIncFiducialLoose &&
                                   kcMuonIDCut);

  // 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
  GenieMultiverseParameters genie_multiverse(100, 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("prongcvn_optimization.root", "recreate");
  optimize->SaveTo(output, "prongcvn_optimize");    
  output->Close();

}