Selection_Opt_Multi.C
Go to the documentation of this file.
1 /* Create histograms for the SINCPi0 selection study
2  Bryan Ramson, Fermilab PostDoc, bjrams87@fnal.gov, May 2020
3  Interactive Use: $ cafe -ss --limit 10 Selection.C
4  Grid Use: submit_cafana.py -n 60 -r R19-10-30-final-prod4.b -t /nova/app/users/btapiaor/ND-run-files -o /pnfs/nova/scratch/users/btapiaor/ND-cafana/1pngE /nova/app/users/btapiaor/ND-run-files/NDAna/SINCpi0/testVar.C
5  Needs a tarball of your current configuration here: /nova/app/users/btapiaor/ND-run-files/NDAna/SINCpi0
6 */
7 #pragma once
8 
10 
11 // analysis cuts/vars
12 //#include "NDAna/nuebarcc_inc/NuebarCCIncCuts.h"
13 //#include "NDAna/numucc_inc/NumuCCIncCuts.h"
14 //#include "NDAna/numucc_inc/NumuCCIncVars.h"
15 //#include "3FlavorAna/Vars/NueVars.h"
16 #include "CAFAna/Cuts/SpillCuts.h"
17 
18 // for data pot
20 #include <iostream>
21 
22 // multiverse includes
24 #include "CAFAna/Vars/XsecTunes.h"
28 
32 
33 using namespace ana;
34 
35 void Selection_Opt_Multi(){//std::string input_file_name = ""){
36  //if(input_file_name == ""){
37  std::map<TString, SpectrumLoader *> loaders;
38  loaders["kNominal"] = new SpectrumLoader("prod_caf_R20-11-25-prod5.1reco.a_nd_genie_N1810j0211a_nonswap_fhc_nova_v08_full_v1_batch1");
39 
40  //define a hist axis to optimize
41  const HistAxis * cvn5labelmuonid_axis = new HistAxis("CVN5LabelMuonIDCut", Binning::Simple(25,0.0,1.0), CVN5LabelMuonID);
42  const HistAxis * cvn2labelpngone_axis = new HistAxis("CVN2LabelPngOneCut", Binning::Simple(25,0.0,1.0), CVN2LabelEMIDPngOne);
43  const HistAxis * cvn2labelpngtwo_axis = new HistAxis("CVN2LabelPngTwoCut", Binning::Simple(25,0.0,1.0), CVN2LabelEMIDPngTwo);
44 
45  const Cut *kOptSignal = new Cut(kTrueNCur && kTruePi0 && kTruePi0100MeV);
46  const Cut *kOptContSel = new Cut(kContCut);
47  const Cut *kOptPreSel = new Cut(kPreSelCut);
48 
49  const Var *kCVWeight = new Var(kXSecCVWgt2020);
50 
51  SystematicDef *sNominal_cvn5labelmuonid = new SystematicDef("sNominal_cvn5labelmuonid", loaders["kNominal"], cvn5labelmuonid_axis, kCVWeight);
52  SystematicDef *sNominal_cvn2labelpngone = new SystematicDef("sNominal_cvn2labelpngone", loaders["kNominal"], cvn2labelpngone_axis, kCVWeight);
53  SystematicDef *sNominal_cvn2labelpngtwo = new SystematicDef("sNominal_cvn2labelpngtwo", loaders["kNominal"], cvn2labelpngtwo_axis, kCVWeight);
54 
55  GenieMultiverseParameters genie_multiverse_cvn5labelmuonid(300, getAllXsecSysts_2020());
56  std::vector<SystShifts> genie_shifts_cvn5labelmuonid = genie_multiverse_cvn5labelmuonid.GetSystShifts();
57  GenieMultiverseParameters genie_multiverse_cvn2labelpngone(300, getAllXsecSysts_2020());
58  std::vector<SystShifts> genie_shifts_cvn2labelpngone = genie_multiverse_cvn2labelpngone.GetSystShifts();
59  GenieMultiverseParameters genie_multiverse_cvn5labelpngtwo(300, getAllXsecSysts_2020());
60  std::vector<SystShifts> genie_shifts_cvn2labelpngtwo = genie_multiverse_cvn5labelpngtwo.GetSystShifts();
61 
62  SystematicDef *sGenieMultiverse_cvn5labelmuonid = new SystematicDef("sGenieMultiverse_cvn5labelmuonid", loaders["kNominal"], cvn5labelmuonid_axis, kCVWeight);
63  SystematicDef *sGenieMultiverse_cvn2labelpngone = new SystematicDef("sGenieMultiverse_cvn2labelpngone", loaders["kNominal"], cvn2labelpngone_axis, kCVWeight);
64  SystematicDef *sGenieMultiverse_cvn2labelpngtwo = new SystematicDef("sGenieMultiverse_cvn2labelpngtwo", loaders["kNominal"], cvn2labelpngtwo_axis, kCVWeight);
65 
66  CutOptimization *optimize_cvn5labelmuonid = new CutOptimization(cvn5labelmuonid_axis, kOptSignal, kOptContSel);
67  CutOptimization *optimize_cvn2labelpngone = new CutOptimization(cvn2labelpngone_axis, kOptSignal, kOptPreSel);
68  CutOptimization *optimize_cvn2labelpngtwo = new CutOptimization(cvn2labelpngtwo_axis, kOptSignal, kOptPreSel);
69 
70  optimize_cvn5labelmuonid->DefineNominal(sNominal_cvn5labelmuonid);
71  optimize_cvn2labelpngone->DefineNominal(sNominal_cvn2labelpngone);
72  optimize_cvn2labelpngtwo->DefineNominal(sNominal_cvn2labelpngtwo);
73 
74  optimize_cvn5labelmuonid->DefineMultiverseSystematic(sGenieMultiverse_cvn5labelmuonid, genie_shifts_cvn5labelmuonid);
75  optimize_cvn2labelpngone->DefineMultiverseSystematic(sGenieMultiverse_cvn2labelpngone, genie_shifts_cvn2labelpngone);
76  optimize_cvn2labelpngtwo->DefineMultiverseSystematic(sGenieMultiverse_cvn2labelpngtwo, genie_shifts_cvn2labelpngtwo);
77 
78  optimize_cvn5labelmuonid->Go();
79  optimize_cvn2labelpngone->Go();
80  optimize_cvn2labelpngtwo->Go();
81 
82  TFile *output = new TFile("test_opt.root","recreate");
83 
84  optimize_cvn5labelmuonid->SaveTo(output, "cvn5labelmuonid_cut_optimize");
85  optimize_cvn2labelpngone->SaveTo(output, "cvn2labelpngone_cut_optimize");
86  optimize_cvn2labelpngtwo->SaveTo(output, "cvn2labelpngtwo_cut_optimize");
87 
88  output->Close();
89  /* }
90  else {
91  //load input file provided in argument
92  TFile *input = TFile::Open(input_file_name.c_str(), "read");
93  CutOptimization *optimize = CutOptimization::LoadFrom(input->GetDirectory("vtx_cut_optimize"),"vtx_cut_optimize").release();
94  input->Close();
95 
96  //save resulting histograms back to file.
97  TFile *results = new TFile("demo_cut_optimization_results.root","recreate");
98  results->mkdir("vtx_cut_optimize/results");
99  TDirectory *save_dir = results->GetDirectory("vtx_cut_optimize");
100 
101  //directory to save plots
102  std::string plot_dump = "demo_cut_optimization_plots";
103 
104  optimize->Optimize(CutOptimization::kBinByBin, // OptimizationMethod_t
105  CutOptimization::kdSigmaOverSigma, // FOM_t
106  kAna2018FHCPOT, // scale to data pot
107  save_dir,
108  plot_dump,
109  true);
110  delete optimize;
111 
112  }*/
113 }
ofstream output
void Selection_Opt_Multi()
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
const Cut kTruePi0100MeV([](const caf::SRProxy *sr){if(sr->vtx.elastic.IsValid==false) return false;if(sr->mc.nnu==0) return false;if(sr->mc.nu.size()< 1) return false;assert(sr->mc.nnu==1);if(sr->mc.nu[0].prim.size()==0) return false;float Pi0Mass=0.134977;float KE=-10.0;float TE=-99.0;int nprim=sr->mc.nu[0].prim.size();int npi0=0;for(int i=0;i< nprim;i++){if(sr->mc.nu[0].prim[i].pdg==111){TE=sr->mc.nu[0].prim[i].p.E;KE=pow(pow(TE, 2.0)-pow(Pi0Mass, 2.0), 0.5);if(KE >=0.1) npi0++;}}if(npi0){return true;}return false;})
const ana::Var kCVWeight
const Cut kContCut
Definition: SINCpi0_Cuts.h:396
void DefineNominal(SystematicDef *nominal)
const Cut kPreSelCut
Definition: SINCpi0_Cuts.h:397
std::vector< const ISyst * > getAllXsecSysts_2020()
void SaveTo(TDirectory *dir, const std::string &name) const
void DefineMultiverseSystematic(SystematicDef *syst, std::vector< Var > mv_weights)
_Var< caf::SRProxy > Var
Definition: Var.h:7
const Var CVN5LabelMuonID([](const caf::SRProxy *sr){if(!(sr->vtx.elastic.IsValid)) return-999.9f;if(sr->vtx.elastic.fuzzyk.npng==0) return-999.9f;float maxscore=-999.9f;for(unsigned int i=0;i< sr->vtx.elastic.fuzzyk.npng;i++) if(maxscore< sr->vtx.elastic.fuzzyk.png[i].spprongcvnpart5label.muonid) maxscore=sr->vtx.elastic.fuzzyk.png[i].spprongcvnpart5label.muonid;return(float(maxscore));})
GenericSystematicDef< Cut, Var > SystematicDef
const Cut kTruePi0([](const caf::SRProxy *sr){if(sr->vtx.elastic.IsValid==false) return false;if(sr->mc.nnu==0) return false;if(sr->mc.nu.size()< 1) return false;assert(sr->mc.nnu==1);if(sr->mc.nu[0].prim.size()==0) return false;float Pi0Mass=0.134977;float KE=-10.0;float TE=-99.0;int nprim=sr->mc.nu[0].prim.size();int npi0=0;for(int i=0;i< nprim;i++){if(sr->mc.nu[0].prim[i].pdg==111){TE=sr->mc.nu[0].prim[i].p.E;KE=pow(pow(TE, 2.0)-pow(Pi0Mass, 2.0), 0.5);if(KE >=0.0) npi0++;}}if(npi0){return true;}return false;})
const Cut kTrueNCur([](const caf::SRProxy *sr){if(sr->vtx.elastic.IsValid==false) return false;if(sr->mc.nnu==0) return false;if(sr->mc.nu.size()< 1) return false;assert(sr->mc.nnu==1);return(sr->mc.nu[0].iscc==0);})
const Var CVN2LabelEMIDPngOne([](const caf::SRProxy *sr){if(!(sr->vtx.elastic.IsValid)) return-999.9f;if(sr->vtx.elastic.fuzzyk.npng==0) return-999.9f;if(sr->vtx.elastic.fuzzyk.npng< 1) return-999.9f;return(float(sr->vtx.elastic.fuzzyk.png[0].spprongcvnpartnumuccemid.emid));})
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
_Cut< caf::SRProxy, caf::SRSpillProxy > Cut
Definition: Cut.h:9
std::vector< Loaders * > loaders
Definition: syst_header.h:386
std::vector< SystShifts > GetSystShifts()
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:107
_HistAxis< Var > HistAxis
Definition: HistAxis.h:10
const Var CVN2LabelEMIDPngTwo([](const caf::SRProxy *sr){if(!(sr->vtx.elastic.IsValid)) return-999.9f;if(sr->vtx.elastic.fuzzyk.npng< 2) return-999.9f;return(float(sr->vtx.elastic.fuzzyk.png[1].spprongcvnpartnumuccemid.emid));})
const Var kXSecCVWgt2020
Definition: XsecTunes.h:105