example_macro.C
Go to the documentation of this file.
1 // This script fills spectra using G4Rwgt weights.
2 // To make plots run example_plot.C over the output of this script.
3 // Author: Cathal Sweeney - csweeney@fnal.gov
4 
5 #include "CAFAna/Core/Binning.h"
6 #include "CAFAna/Core/Spectrum.h"
8 #include "CAFAna/Core/Var.h"
9 
10 #include "TFile.h"
11 
14 
15 //#include "G4Rwgt/plotting_scripts/Geant4WeightVars.h"
17 
18 using namespace ana;
19 
21 {
22  struct VarStruct {
25  Binning bins;
26  const NuTruthVar var;
27 
28  };
29 
30  std::vector<VarStruct> var_vec{
31  {"True neutrino energy", "Enu", Binning::Simple(100,0.,10.), kNuE_NT},
32  {"True pion energy", "Epi", Binning::Simple(100,0.,5.), kPiE_NT}
33  };
34 
35 
36  std::vector<std::pair<std::string, const NuTruthCut>> kMyCut_vec{
37  {"all" , kTrue_NT}
38 /*
39  {"numucc", kIsCC_NT && kIsNuMu_NT},
40  {"numubarcc", kIsCC_NT && kIsNuMuBar_NT},
41  {"nuecc", kIsCC_NT && kIsNuE_NT},
42  {"nuebarcc", kIsCC_NT && kIsNuEBar_NT},
43  {"nc", !kIsCC_NT}
44 */
45  };
46 
47  std::vector<std::string> systName_vec{
48  "tot"//,
49 /*
50  "cex",
51  "dcex",
52  "qe",
53  "abs",
54  "prod",
55 */
56  };
57 
58  std::vector<const ReinteractionSyst*> syst_vec{
59  &kPiplusTotSyst/*,
60  &kPiplusCexSyst,
61  &kPiplusDcexSyst,
62  &kPiplusQeSyst,
63  &kPiplusAbsSyst,
64  &kPiplusProdSyst
65  */
66  };
67 
68 
70 
71 
72  SpectrumLoader* loader = new SpectrumLoader("prod_caf_R20-11-25-prod5.1reco.a_nd_genie_N1810j0211a_nonswap_fhc_nova_v08_full_v1");
73 
74 
75  std::string outName = "outFile.root";
76 
77 
78  std::vector<std::pair<std::string, Spectrum *>> mySpec_vec; // Holds label and unweighted (nominal) spectra
79  std::vector<std::pair<std::string, Spectrum *>> myUniv_vec; // Holds label and weighted spectra
80 
81 
82  for(uint iVar=0; iVar<var_vec.size(); iVar++){
83  std::string var_label = var_vec[iVar].label;
84  std::string shortName = var_vec[iVar].shortName;
85 
86  const NuTruthVar nuTruthVar = var_vec[iVar].var;
87  Var var = VarFromNuTruthVar(nuTruthVar);
88 
89  Binning bins = var_vec[iVar].bins;
90 
91  for(uint iCut=0; iCut<kMyCut_vec.size(); iCut++){
92  std::string cut_label = kMyCut_vec[iCut].first;
93  const NuTruthCut nuTruthCut = kMyCut_vec[iCut].second && kSignal;
94  Cut cut = CutFromNuTruthCut(nuTruthCut);
95 
96 
97 
98  std::string label = shortName + "_" + cut_label;
99 
100 
101  Spectrum * sNom = new Spectrum(var_label, bins,
102  *loader, var,
103  cut);
104 
105 
106  mySpec_vec.push_back( {label + "_nom", sNom} );
107 
108 
109  for(uint iSyst=0; iSyst<syst_vec.size(); iSyst++){
110  std::string syst_label = systName_vec[iSyst];
111  const ReinteractionSyst* syst = syst_vec[iSyst];
112 
113 
114  Spectrum * sPlus1Sigma = new Spectrum(var_label, bins,
115  *loader, var,
116  cut,
117  SystShifts(syst, 1.00)); // Sigma is 20%, so here we set knob to 1.2
118  myUniv_vec.push_back( {label + "_" + syst_label + "_plus1sigma", sPlus1Sigma} );
119 
120 
121  Spectrum * sMinus1Sigma = new Spectrum(var_label, bins,
122  *loader, var,
123  cut,
124  SystShifts(syst, -1.00)); // set knob to 0.8
125  myUniv_vec.push_back( {label + "_" + syst_label + "_minus1sigma", sMinus1Sigma} );
126 
127  } // end for(iSyst)
128  } // end for(iCut)
129  }// end for(iVar)
130 
131  //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
132 
133  loader->Go();
134 
135 
136  TFile * outFile = new TFile( outName.c_str(), "recreate");
137 
138  TDirectory* base_dir = outFile->mkdir("foo");
139 
140  for(uint i=0; i< mySpec_vec.size(); ++i){
141  std::string name = mySpec_vec[i].first;
142  mySpec_vec[i].second->SaveTo( base_dir, name.c_str() );
143  }//end for(i)
144 
145  for(uint i=0; i< myUniv_vec.size(); ++i){
146  std::string name = myUniv_vec[i].first;
147  myUniv_vec[i].second->SaveTo( base_dir, name.c_str() );
148  }//end for(i)
149 
150 
151 
152 }
153 
const XML_Char * name
Definition: expat.h:151
const NuTruthCut kHasPiplus_NT([](const caf::SRNeutrinoProxy *nu){int nPiplus=0;int nPrim=nu->prim.size();for(int i=0;i< nPrim;i++){if(nu->prim[i].pdg==211) nPiplus+=1;}if(nPiplus > 0) return true;else return false;})
Represent the binning of a Spectrum&#39;s x-axis.
Definition: Binning.h:16
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
Base class for GEANT4 reinteraction systs.
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
NuTruthVar kNuE_NT([](const caf::SRNeutrinoProxy *nu){float nuE=-5.;nuE=nu->E;if(nuE< 0) std::cout<< "BAR \n";return nuE;})
const NuTruthCut kTrueFiducial_NT([](const caf::SRNeutrinoProxy *nu){return(nu->vtx.X()< my_vtxmax.X()&& nu->vtx.X() > my_vtxmin.X()&& nu->vtx.Y()< my_vtxmax.Y()&& nu->vtx.Y() > my_vtxmin.Y()&& nu->vtx.Z()< my_vtxmax.Z()&& nu->vtx.Z() > my_vtxmin.Z());})
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
const char * label
const NuTruthCut kTrue_NT([](const caf::SRNeutrinoProxy *nu){return true;})
TFile * outFile
Definition: PlotXSec.C:135
virtual void Go() override
Load all the registered spectra.
loader
Definition: demo0.py:10
std::vector< float > Spectrum
Definition: Constants.h:728
const Binning bins
Definition: NumuCC_CPiBin.h:8
const Cut kSignal
Definition: SINCpi0_Cuts.h:372
const ReinteractionSyst kPiplusTotSyst("PiplusTot","#pi^{+} Total", 211, EG4RwChannel::kTot)
const Cut cut
Definition: exporter_fd.C:30
Var VarFromNuTruthVar(const NuTruthVar &stv, double _default)
Definition: Var.cxx:7
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
void example_macro()
Definition: example_macro.C:20
string shortName
THUMBNAIL BLOCK: We need to make a thumbnail for each.
Template for Cut and SpillCut.
Definition: Cut.h:15
Template for Var and SpillVar.
NuTruthVar kPiE_NT([](const caf::SRNeutrinoProxy *nu){double maxE=-5.0;int nPrim=nu->prim.size();for(int i=0;i< nPrim;i++){auto &prim=nu->prim[i];if(prim.pdg!=211) continue;double energy=prim.p.E;if(energy > maxE) maxE=energy;}if(maxE< 0) std::cout<< "FOO \n";return maxE;})
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:107
Cut CutFromNuTruthCut(const NuTruthCut &stc)
Definition: Cut.cxx:7
enum BeamMode string
unsigned int uint