runCheatDecomp.C
Go to the documentation of this file.
3 #include "CAFAna/Core/Binning.h"
4 #include "CAFAna/Core/HistAxis.h"
5 #include "CAFAna/Vars/Vars.h"
6 #include "CAFAna/Core/Spectrum.h"
9 #include "CAFAna/Core/ISyst.h"
10 #include "CAFAna/Systs/Systs.h"
11 #include "CAFAna/Core/Utilities.h"
12 
13 #include "nugen/NuReweight/ReweightLabels.h"
14 
15 
16 
17 #include <iostream>
18 #include <sstream>
19 #include <map>
20 #include <utility>
21 
22 #include "TBrowser.h"
23 #include "TFile.h"
25 
26 using namespace ana;
27 
28 
29 void runCheatDecomp(TString fileListTextFileName)
30 {
31 
32 
33  std::map<std::string, rwgt::ReweightLabel_t> parameters;
34  // NCEL tweaking parameters:
35  parameters["ReweightMaNCEL"] = rwgt::fReweightMaNCEL; ///< tweak Ma NCEL, affects dsigma(NCEL)/dQ2 both in shape and normalization
36  parameters["ReweightEtaNCEL"] = rwgt::fReweightEtaNCEL; ///< tweak NCEL strange axial form factor eta, affects dsigma(NCEL)/dQ2 both in shape and normalization
37  // CCQE tweaking parameters:
38  parameters["ReweightNormCCQE"] = rwgt::fReweightNormCCQE; ///< tweak CCQE normalization (energy independent)
39  parameters["ReweightNormCCQEenu"] = rwgt::fReweightNormCCQEenu; ///< tweak CCQE normalization (maintains dependence on neutrino energy)
40  parameters["ReweightMaCCQEshape"] = rwgt::fReweightMaCCQEshape; ///< tweak Ma CCQE, affects dsigma(CCQE)/dQ2 in shape only (normalized to constant integral)
41  parameters["ReweightMaCCQE"] = rwgt::fReweightMaCCQE; ///< tweak Ma CCQE, affects dsigma(CCQE)/dQ2 both in shape and normalization
42  parameters["ReweightVecCCQEshape"] = rwgt::fReweightVecCCQEshape;///< tweak elastic nucleon form factors (BBA/default -> dipole) - shape only effect of dsigma(CCQE)/dQ2
43  // Resonance neutrino-production tweaking parameters:
44  parameters["ReweightNormCCRES"] = rwgt::fReweightNormCCRES; ///< tweak CCRES normalization
45  parameters["ReweightMaCCRESshape"] = rwgt::fReweightMaCCRESshape; ///< tweak Ma CCRES, affects d2sigma(CCRES)/dWdQ2 in shape only (normalized to constant integral)
46  parameters["ReweightMvCCRESshape"] = rwgt::fReweightMvCCRESshape; ///< tweak Mv CCRES, affects d2sigma(CCRES)/dWdQ2 in shape only (normalized to constant integral)
47  parameters["ReweightMaCCRES"] = rwgt::fReweightMaCCRES; ///< tweak Ma CCRES, affects d2sigma(CCRES)/dWdQ2 both in shape and normalization
48  parameters["ReweightMvCCRES"] = rwgt::fReweightMvCCRES; ///< tweak Mv CCRES, affects d2sigma(CCRES)/dWdQ2 both in shape and normalization
49 
50  parameters["ReweightNormNCRES"] = rwgt::fReweightNormNCRES; ///< tweak NCRES normalization
51  parameters["ReweightMaNCRESshape"] = rwgt::fReweightMaNCRESshape; ///< tweak Ma NCRES, affects d2sigma(NCRES)/dWdQ2 in shape only (normalized to constant integral)
52  parameters["ReweightMvNCRESshape"] = rwgt::fReweightMvNCRESshape; ///< tweak Mv NCRES, affects d2sigma(NCRES)/dWdQ2 in shape only (normalized to constant integral)
53  parameters["ReweightMaNCRES"] = rwgt::fReweightMaNCRES; ///< tweak Ma NCRES, affects d2sigma(NCRES)/dWdQ2 both in shape and normalization
54  parameters["ReweightMvNCRES"] = rwgt::fReweightMvNCRES; ///< tweak Mv NCRES, affects d2sigma(NCRES)/dWdQ2 both in shape and normalization
55 
56  // Coherent pion production tweaking parameters:
57  parameters["ReweightMaCOHpi"] = rwgt::fReweightMaCOHpi; ///< tweak Ma for COH pion production
58  parameters["ReweightR0COHpi"] = rwgt::fReweightR0COHpi; ///< tweak R0 for COH pion production
59  // Non-resonance background tweaking parameters:
60  parameters["ReweightRvpCC1pi"] = rwgt::fReweightRvpCC1pi; ///< tweak the 1pi non-RES bkg in the RES region, for v+p CC
61  parameters["ReweightRvpCC2pi"] = rwgt::fReweightRvpCC2pi; ///< tweak the 2pi non-RES bkg in the RES region, for v+p CC
62  parameters["ReweightRvpNC1pi"] = rwgt::fReweightRvpNC1pi; ///< tweak the 1pi non-RES bkg in the RES region, for v+p NC
63  parameters["ReweightRvpNC2pi"] = rwgt::fReweightRvpNC2pi; ///< tweak the 2pi non-RES bkg in the RES region, for v+p NC
64  parameters["ReweightRvnCC1pi"] = rwgt::fReweightRvnCC1pi; ///< tweak the 1pi non-RES bkg in the RES region, for v+n CC
65  parameters["ReweightRvnCC2pi"] = rwgt::fReweightRvnCC2pi; ///< tweak the 2pi non-RES bkg in the RES region, for v+n CC
66  parameters["ReweightRvnNC1pi"] = rwgt::fReweightRvnNC1pi; ///< tweak the 1pi non-RES bkg in the RES region, for v+n NC
67  parameters["ReweightRvnNC2pi"] = rwgt::fReweightRvnNC2pi; ///< tweak the 2pi non-RES bkg in the RES region, for v+n NC
68  parameters["ReweightRvbarpCC1pi"] = rwgt::fReweightRvbarpCC1pi; ///< tweak the 1pi non-RES bkg in the RES region, for vbar+p CC
69  parameters["ReweightRvbarpCC2pi"] = rwgt::fReweightRvbarpCC2pi; ///< tweak the 2pi non-RES bkg in the RES region, for vbar+p CC
70  parameters["ReweightRvbarpNC1pi"] = rwgt::fReweightRvbarpNC1pi; ///< tweak the 1pi non-RES bkg in the RES region, for vbar+p NC
71  parameters["ReweightRvbarpNC2pi"] = rwgt::fReweightRvbarpNC2pi; ///< tweak the 2pi non-RES bkg in the RES region, for vbar+p NC
72  parameters["ReweightRvbarnCC1pi"] = rwgt::fReweightRvbarnCC1pi; ///< tweak the 1pi non-RES bkg in the RES region, for vbar+n CC
73  parameters["ReweightRvbarnCC2pi"] = rwgt::fReweightRvbarnCC2pi; ///< tweak the 2pi non-RES bkg in the RES region, for vbar+n CC
74  parameters["ReweightRvbarnNC1pi"] = rwgt::fReweightRvbarnNC1pi; ///< tweak the 1pi non-RES bkg in the RES region, for vbar+n NC
75  parameters["ReweightRvbarnNC2pi"] = rwgt::fReweightRvbarnNC2pi; ///< tweak the 2pi non-RES bkg in the RES region, for vbar+n NC
76  // DIS tweaking parameters - applied for DIS events with (Q2>Q2o, W>Wo),
77  // typically Q2okReweight =1GeV^2, WokReweight =1.7-2.0GeV
78  parameters["ReweightAhtBY"] = rwgt::fReweightAhtBY; ///< tweak the Bodek-Yang model parameter A_{ht} - incl. both shape and normalization effect
79  parameters["ReweightBhtBY"] = rwgt::fReweightBhtBY; ///< tweak the Bodek-Yang model parameter B_{ht} - incl. both shape and normalization effect
80  parameters["ReweightCV1uBY"] = rwgt::fReweightCV1uBY; ///< tweak the Bodek-Yang model parameter CV1u - incl. both shape and normalization effect
81  parameters["ReweightCV2uBY"] = rwgt::fReweightCV2uBY; ///< tweak the Bodek-Yang model parameter CV2u - incl. both shape and normalization effect
82  parameters["ReweightAhtBYshape"] = rwgt::fReweightAhtBYshape; ///< tweak the Bodek-Yang model parameter A_{ht} - shape only effect to d2sigma(DIS)/dxdy
83  parameters["ReweightBhtBYshape"] = rwgt::fReweightBhtBYshape; ///< tweak the Bodek-Yang model parameter B_{ht} - shape only effect to d2sigma(DIS)/dxdy
84  parameters["ReweightCV1uBYshape"] = rwgt::fReweightCV1uBYshape; ///< tweak the Bodek-Yang model parameter CV1u - shape only effect to d2sigma(DIS)/dxdy
85  parameters["ReweightCV2uBYshape"] = rwgt::fReweightCV2uBYshape; ///< tweak the Bodek-Yang model parameter CV2u - shape only effect to d2sigma(DIS)/dxdy
86  parameters["ReweightNormDISCC"] = rwgt::fReweightNormDISCC; ///< tweak the inclusive DIS CC normalization (not currently working in genie)
87  parameters["ReweightRnubarnuCC"] = rwgt::fReweightRnubarnuCC; ///< tweak the ratio of \sigma(\bar\nu CC) / \sigma(\nu CC) (not currently working in genie)
88  parameters["ReweightDISNuclMod"] = rwgt::fReweightDISNuclMod; ///< tweak DIS nuclear modification (shadowing, anti-shadowing, EMC). Does not appear to be working in GENIE at the moment
89  //
90  parameters["ReweightNC"] = rwgt::fReweightNC; ///<
91 
92 
93  //
94  // Hadronization (free nucleon target)
95  //
96 
97  parameters["ReweightAGKY_xF1pi"] = rwgt::fReweightAGKY_xF1pi; ///< tweak xF distribution for low multiplicity (N + pi) DIS f/s produced by AGKY
98  parameters["ReweightAGKY_pT1pi"] = rwgt::fReweightAGKY_pT1pi; ///< tweak pT distribution for low multiplicity (N + pi) DIS f/s produced by AGKY
99 
100 
101  //
102  // Medium-effects to hadronization
103  //
104 
105  parameters["ReweightFormZone"] = rwgt::fReweightFormZone; ///< tweak formation zone
106 
107 
108  //
109  // Intranuclear rescattering systematics.
110  // There are 2 sets of parameters:
111  // - parameters that control the total rescattering probability, P(total)
112  // - parameters that control the fraction of each process (`fate'), given a total rescat. prob., P(fate|total)
113  // These parameters are considered separately for pions and nucleons.
114  //
115 
116  parameters["ReweightMFP_pi"] = rwgt::fReweightMFP_pi; ///< tweak mean free path for pions
117  parameters["ReweightMFP_N"] = rwgt::fReweightMFP_N; ///< tweak mean free path for nucleons
118  parameters["ReweightFrCEx_pi"] = rwgt::fReweightFrCEx_pi; ///< tweak charge exchange probability for pions, for given total rescattering probability
119  parameters["ReweightFrElas_pi"] = rwgt::fReweightFrElas_pi; ///< tweak elastic probability for pions, for given total rescattering probability
120  parameters["ReweightFrInel_pi"] = rwgt::fReweightFrInel_pi; ///< tweak inelastic probability for pions, for given total rescattering probability
121  parameters["ReweightFrAbs_pi"] = rwgt::fReweightFrAbs_pi; ///< tweak absorption probability for pions, for given total rescattering probability
122  parameters["ReweightFrPiProd_pi"] = rwgt::fReweightFrPiProd_pi; ///< tweak pion production probability for pions, for given total rescattering probability
123  parameters["ReweightFrCEx_N"] = rwgt::fReweightFrCEx_N; ///< tweak charge exchange probability for nucleons, for given total rescattering probability
124  parameters["ReweightFrElas_N"] = rwgt::fReweightFrElas_N; ///< tweak elastic probability for nucleons, for given total rescattering probability
125  parameters["ReweightFrInel_N"] = rwgt::fReweightFrInel_N; ///< tweak inelastic probability for nucleons, for given total rescattering probability
126  parameters["ReweightFrAbs_N"] = rwgt::fReweightFrAbs_N; ///< tweak absorption probability for nucleons, for given total rescattering probability
127  parameters["ReweightFrPiProd_N"] = rwgt::fReweightFrPiProd_N; ///< tweak pion production probability for nucleons, for given total rescattering probability
128 
129  //
130  // Nuclear model
131  //
132 
133  parameters["ReweightCCQEPauliSupViaKF"] = rwgt::fReweightCCQEPauliSupViaKF; ///<
134  parameters["ReweightCCQEMomDistroFGtoSF"] = rwgt::fReweightCCQEMomDistroFGtoSF; ///<
135 
136  //
137  // Resonance decays
138  //
139 
140  parameters["ReweightBR1gamma"] = rwgt::fReweightBR1gamma; ///< tweak Resonance -> X + gamma branching ratio, eg Delta+(1232) -> p gamma
141  parameters["ReweightBR1eta"] = rwgt::fReweightBR1eta; ///< tweak Resonance -> X + eta branching ratio, eg N+(1440) -> p eta
142  parameters["ReweightTheta_Delta2Npi"] = rwgt::fReweightTheta_Delta2Npi; ///< distort pi angular distribution in Delta -> N + pi
143 
144 
145 
146  int nParameters = 0;
147  for(auto parameter:parameters) nParameters++;
148  std::cout << "N Parameters : " << nParameters << std::endl;
149 
150  std::string fileListTextFile = (std::string)fileListTextFileName;
151  assert(fileListTextFile.find(".txt") != std::string::npos && "Failed to find .txt in argument. We're going to to replace it with .root for the output file, so it needs to be there.");
152  std::vector<std::string> fileList = LoadFileList(fileListTextFile);
153  SpectrumLoader loader (fileList);
154 
155 
156 
157  std::map<std::string, std::map<int, SystShifts*> > shifts;
158  std::vector<int> sigmas;
159  sigmas.push_back(-2);
160  sigmas.push_back(-1);
161  sigmas.push_back(1);
162  sigmas.push_back(2);
163 
164  for (const auto & parameter:parameters)
165  {
166  for(const auto & sigma:sigmas)
167  {
168 
169  shifts[parameter.first][sigma] = new SystShifts(GetGenieKnobSyst(parameter.second), sigma);
170 
171  }
172  }
173 
174  const Binning bins = Binning::Simple(50, 0, 5);
175 
176  std::map<std::string, std::pair<Cut*, HistAxis*> > samples;
177 // samples["qe"] = std::pair<Cut*, HistAxis* >(new Cut(kNumuNCRej && kNumuContainND && kNumuQuality && kQE), new HistAxis("QE Energy [GeV]", bins, kQEE));
178 // samples["nonqe"] = std::pair<Cut*, HistAxis* >(new Cut(kNumuNCRej && kNumuContainND && kNumuQuality && !kQE), new HistAxis("Non-QE Energy [GeV]", bins, kCCE));
179  samples["allcc"] = std::pair<Cut*, HistAxis* >(new Cut(kNumuNCRej && kNumuContainND && kNumuQuality), new HistAxis("Reconstructed Neutrino Energy [GeV]", bins, kCCE));
180 
181  std::map<std::string, CheatDecomp*> decomps;
182 
183  std::map<std::string, std::map<std::string, std::map<int, CheatDecomp*> > > shiftedDecomps;
184 
185 
186 
187  for(const auto & sample:samples){
188  std::cout << "Sample: " << sample.first << std::endl;
189  decomps[sample.first] = new CheatDecomp(loader, *sample.second.second, *sample.second.first);
190  for(const auto & parameter:parameters)
191  {
192  for(const auto & shift:shifts[parameter.first])
193  {
194  shiftedDecomps[sample.first][parameter.first][shift.first] = new CheatDecomp(loader, *sample.second.second, *sample.second.first, *shift.second);
195 
196  }
197  }
198  }
199 
200  loader.Go();
201 
202 
203  std::map<std::string, TDirectory*> dirs;
204 
205  TFile outputFile(fileListTextFile.replace(fileListTextFile.find(".txt"), 4, ".root").c_str(), "RECREATE");
206  for(const auto & decomp:decomps){
207 
208 
209  dirs[decomp.first] = outputFile.mkdir(decomp.first.c_str());
210  decomp.second->SaveTo(dirs[decomp.first]);
211 
212  }
213 
214  for(const auto & parameter:shiftedDecomps)
215  {
216  for(const auto & shift:parameter.second)
217  {
218  for(const auto & shiftedDecomp:shift.second)
219  {
220  std::stringstream dirName("");
221  dirName << parameter.first << "_" << shift.first << "_" << shiftedDecomp.first ;
222  dirs[dirName.str()] = outputFile.mkdir(dirName.str().c_str());
223  shiftedDecomp.second->SaveTo(dirs[dirName.str()]);
224  }
225  }
226 
227  }
228 }
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
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
GenericCut< caf::SRProxy > Cut
Representation of a cut (selection) to be applied to a caf::StandardRecord object.
Definition: Cut.h:96
const Cut kNumuContainND([](const caf::SRProxy *sr){return( sr->trk.kalman.ntracks > sr->trk.kalman.idxremid &&sr->slc.ncellsfromedge > 1 &&sr->slc.firstplane > 1 &&sr->slc.lastplane< 212 &&sr->trk.kalman.tracks[0].start.Z()< 1150 &&( sr->trk.kalman.tracks[0].stop.Z()< 1275 ||sr->sel.contain.kalyposattrans< 55) &&( sr->energy.numu.ndhadcalcatE +sr->energy.numu.ndhadcaltranE)< 0.03 &&sr->sel.contain.kalfwdcellnd > 4 &&sr->sel.contain.kalbakcellnd > 8);})
Definition: NumuCuts.h:22
list fileList
Definition: cafExposure.py:30
std::vector< std::string > LoadFileList(const std::string &listfile)
Read list of input files from a text file, one per line.
Definition: Utilities.cxx:435
const Var kCCE
Definition: NumuVars.h:21
void runCheatDecomp(TString fileListTextFileName)
virtual void Go() override
Load all the registered spectra.
const Cut kNumuNCRej([](const caf::SRProxy *sr){return(sr->sel.remid.pid >0.75);})
Definition: NumuCuts.h:24
loader
Definition: demo0.py:10
double sigma(TH1F *hist, double percentile)
OStream cout
Definition: OStream.cxx:6
const Binning bins
Definition: NumuCC_CPiBin.h:8
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
GenericHistAxis< Var > HistAxis
Definition: HistAxis.h:111
assert(nhit_max >=nhit_nbins)
const Cut kNumuQuality
Definition: NumuCuts.h:18
const NOvARwgtSyst * GetGenieKnobSyst(rwgt::ReweightLabel_t knobIdx, std::string altName, std::string altLabel)
Convenience function to get a GENIE knob syst. (Allows using the GENIE knob name & description as the...
Definition: XSecSysts.cxx:119
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:38
Just return the ND truth spectra as the decomposition.
Definition: CheatDecomp.h:10