multiverse_reweighting.C
Go to the documentation of this file.
1 #ifdef __CINT__
3  int nfiles = INT_MAX,
4  int nuniverses = 10,
5  string klist = "knob_confs/maccqe_not_varied.conf",
6  string flistpn = "/nova/ana/users/slin/sam_definition_file_lists/prod_caf_R16-12-20-prod3recopreview.d_nd_genie_nonswap_fhc_nova_v08_full_v3.txt")
7 {
8  std::cout << "Sorry, you must run in compiled mode" << std::endl;
9 }
10 #else
11 
12 #include <fstream>
13 #include <sstream>
15 #include "CAFAna/Core/Spectrum.h"
16 #include "CAFAna/Core/Utilities.h"
17 #include "CAFAna/Cuts/TruthCuts.h"
18 #include "CAFAna/Systs/Systs.h"
20 #include "TCanvas.h"
21 #include "THStack.h"
22 #include "TROOT.h"
23 
24 using namespace ana;
25 
27 
28 map<string, rwgt::ReweightLabel_t> name2knob {
29  {"ReweightNull",rwgt::kReweightNull},
30  {"ReweightMaNCEL",rwgt::fReweightMaNCEL},
31  {"ReweightEtaNCEL",rwgt::fReweightEtaNCEL},
32  {"ReweightNormCCQE",rwgt::fReweightNormCCQE},
33  {"ReweightNormCCQEenu",rwgt::fReweightNormCCQEenu},
34  {"ReweightMaCCQEshape",rwgt::fReweightMaCCQEshape},
35  {"ReweightMaCCQE",rwgt::fReweightMaCCQE},
36  {"ReweightVecCCQEshape",rwgt::fReweightVecCCQEshape},
37  {"ReweightNormCCRES",rwgt::fReweightNormCCRES},
38  {"ReweightMaCCRESshape",rwgt::fReweightMaCCRESshape},
39  {"ReweightMvCCRESshape",rwgt::fReweightMvCCRESshape},
40  {"ReweightMaCCRES",rwgt::fReweightMaCCRES},
41  {"ReweightMvCCRES",rwgt::fReweightMvCCRES},
42  {"ReweightNormNCRES",rwgt::fReweightNormNCRES},
43  {"ReweightMaNCRESshape",rwgt::fReweightMaNCRESshape},
44  {"ReweightMvNCRESshape",rwgt::fReweightMvNCRESshape},
45  {"ReweightMaNCRES",rwgt::fReweightMaNCRES},
46  {"ReweightMvNCRES",rwgt::fReweightMvNCRES},
47  {"ReweightMaCOHpi",rwgt::fReweightMaCOHpi},
48  {"ReweightR0COHpi",rwgt::fReweightR0COHpi},
49  {"ReweightRvpCC1pi",rwgt::fReweightRvpCC1pi},
50  {"ReweightRvpCC2pi",rwgt::fReweightRvpCC2pi},
51  {"ReweightRvpNC1pi",rwgt::fReweightRvpNC1pi},
52  {"ReweightRvpNC2pi",rwgt::fReweightRvpNC2pi},
53  {"ReweightRvnCC1pi",rwgt::fReweightRvnCC1pi},
54  {"ReweightRvnCC2pi",rwgt::fReweightRvnCC2pi},
55  {"ReweightRvnNC1pi",rwgt::fReweightRvnNC1pi},
56  {"ReweightRvnNC2pi",rwgt::fReweightRvnNC2pi},
57  {"ReweightRvbarpCC1pi",rwgt::fReweightRvbarpCC1pi},
58  {"ReweightRvbarpCC2pi",rwgt::fReweightRvbarpCC2pi},
59  {"ReweightRvbarpNC1pi",rwgt::fReweightRvbarpNC1pi},
60  {"ReweightRvbarpNC2pi",rwgt::fReweightRvbarpNC2pi},
61  {"ReweightRvbarnCC1pi",rwgt::fReweightRvbarnCC1pi},
62  {"ReweightRvbarnCC2pi",rwgt::fReweightRvbarnCC2pi},
63  {"ReweightRvbarnNC1pi",rwgt::fReweightRvbarnNC1pi},
64  {"ReweightRvbarnNC2pi",rwgt::fReweightRvbarnNC2pi},
65  {"ReweightAhtBY",rwgt::fReweightAhtBY},
66  {"ReweightBhtBY",rwgt::fReweightBhtBY},
67  {"ReweightCV1uBY",rwgt::fReweightCV1uBY},
68  {"ReweightCV2uBY",rwgt::fReweightCV2uBY},
69  {"ReweightAhtBYshape",rwgt::fReweightAhtBYshape},
70  {"ReweightBhtBYshape",rwgt::fReweightBhtBYshape},
71  {"ReweightCV1uBYshape",rwgt::fReweightCV1uBYshape},
72  {"ReweightCV2uBYshape",rwgt::fReweightCV2uBYshape},
73  {"ReweightNormDISCC",rwgt::fReweightNormDISCC},
74  {"ReweightRnubarnuCC",rwgt::fReweightRnubarnuCC},
75  {"ReweightDISNuclMod",rwgt::fReweightDISNuclMod},
76  {"ReweightNC",rwgt::fReweightNC},
77  {"ReweightAGKY_xF1pi",rwgt::fReweightAGKY_xF1pi},
78  {"ReweightAGKY_pT1pi",rwgt::fReweightAGKY_pT1pi},
79  {"ReweightFormZone",rwgt::fReweightFormZone},
80  {"ReweightMFP_pi",rwgt::fReweightMFP_pi},
81  {"ReweightMFP_N",rwgt::fReweightMFP_N},
82  {"ReweightFrCEx_pi",rwgt::fReweightFrCEx_pi},
83  {"ReweightFrElas_pi",rwgt::fReweightFrElas_pi},
84  {"ReweightFrInel_pi",rwgt::fReweightFrInel_pi},
85  {"ReweightFrAbs_pi",rwgt::fReweightFrAbs_pi},
86  {"ReweightFrPiProd_pi",rwgt::fReweightFrPiProd_pi},
87  {"ReweightFrCEx_N",rwgt::fReweightFrCEx_N},
88  {"ReweightFrElas_N",rwgt::fReweightFrElas_N},
89  {"ReweightFrInel_N",rwgt::fReweightFrInel_N},
90  {"ReweightFrAbs_N",rwgt::fReweightFrAbs_N},
91  {"ReweightFrPiProd_N",rwgt::fReweightFrPiProd_N},
92  {"ReweightCCQEPauliSupViaKF",rwgt::fReweightCCQEPauliSupViaKF},
93  {"ReweightCCQEMomDistroFGtoSF",rwgt::fReweightCCQEMomDistroFGtoSF},
94  {"ReweightBR1gamma",rwgt::fReweightBR1gamma},
95  {"ReweightBR1eta",rwgt::fReweightBR1eta},
96  {"ReweightTheta_Delta2Npi",rwgt::fReweightTheta_Delta2Npi},
97  {"ReweightZNormCCQE",rwgt::fReweightZNormCCQE},
98  {"ReweightZExpA1CCQE",rwgt::fReweightZExpA1CCQE},
99  {"ReweightZExpA2CCQE",rwgt::fReweightZExpA2CCQE},
100  {"ReweightZExpA3CCQE",rwgt::fReweightZExpA3CCQE},
101  {"ReweightZExpA4CCQE",rwgt::fReweightZExpA4CCQE},
102  {"ReweightAxFFCCQEshape",rwgt::fReweightAxFFCCQEshape}
103 };
104 
105 unsigned int count_enabled_knobs(map<rwgt::ReweightLabel_t, bool>& table)
106 {
107  unsigned int counter = 0;
108  for(auto it: table)
109  if(it.second == true) counter++;
110 
111  return counter;
112 }
113 
114 map<rwgt::ReweightLabel_t, bool> build_knob_switch(string fpn)
115 {
116  /// Check if the table of enabled knobs exists.
117  /// If it does not, default to all enabled.
118  ifstream inf(fpn.c_str());
119  bool fileexists = inf.good();
120  map<rwgt::ReweightLabel_t, bool> knob_switch;
121  for(auto it: name2knob)
122  knob_switch[it.second] = true;
123 
124  if(fileexists)
125  {
126  while (inf)
127  {
128  string line;
129  getline(inf, line);
130  string name;
131  bool value;
132  stringstream sline(line);
133  sline >> name >> value;
134  if(name2knob.find(name) != name2knob.end()) knob_switch[name2knob[name]] = value;
135  }
136  }
137  else cout << "Cannot open the knob_switch.conf file. Default to all enabled." << endl;
138 
139  return knob_switch;
140 }
141 
142 Spectrum* create_universe(SpectrumLoader& loader, TH1F* hrn, map<rwgt::ReweightLabel_t, bool>& kswitch)
143 {
144  TRandom3 normalDist(0);
145  map<const ISyst*, double> parsInCurUni;
146  //~ for(int knobIdx = rwgt::fReweightMaNCEL; knobIdx != rwgt::fReweightAxFFCCQEshape; knobIdx++)
147  for(int knobIdx = rwgt::kReweightNull; knobIdx != rwgt::fReweightAxFFCCQEshape; knobIdx++)
148  {
149  /// The rwgt::fReweightZExpA?CCQE results in a NAN condition. Ignore it for now.
150  if(knobIdx == (int)rwgt::fReweightZExpA1CCQE) continue;
151  if(knobIdx == (int)rwgt::fReweightZExpA2CCQE) continue;
152  if(knobIdx == (int)rwgt::fReweightZExpA3CCQE) continue;
153  if(knobIdx == (int)rwgt::fReweightZExpA4CCQE) continue;
154 
155  rwgt::ReweightLabel_t knoblabel = (rwgt::ReweightLabel_t)knobIdx;
156  if(kswitch[knoblabel] == false) continue;
157  double nsigma = normalDist.Gaus();
158  parsInCurUni[GetGenieRwgtSyst(knoblabel)] = nsigma;
159  hrn->Fill(nsigma);
160  }
161  cout << "Universe created. " << parsInCurUni.size() << " knobs enabled." << endl;
162  return new Spectrum(loader, kRecoEStandardBinning, total_selection, parsInCurUni);
163 }
164 
166  int nfiles = INT_MAX,
167  int nuniverses = 10,
168  string klist = "knob_confs/maccqe_not_varied.conf",
169  string flistpn = "/nova/ana/users/slin/sam_definition_file_lists/prod_caf_R16-12-20-prod3recopreview.d_nd_genie_nonswap_fhc_nova_v08_full_v3.txt")
170 {
171  /// load file names into a vector
172  vector<string> flist = LoadFileList(flistpn);
173  if(!flist.size())
174  {
175  cout << "Not able to load the file list. Quit." << endl;
176  return;
177  }
178  /// only run over the specified number of files
179  flist.resize(min(min(nfiles, INT_MAX), (int)flist.size()));
180 
181  /// string to specify knob configuration
182  string knobstr = klist;
183  knobstr = knobstr.erase(knobstr.find_last_of("."), string::npos);
184  knobstr = knobstr.substr(knobstr.rfind("/")+1, string::npos);
185 
186  /// spectrum loaders for MC
187  SpectrumLoader mcloader(flist);
188 
189  /// containers for results
190  vector<Spectrum*> reco_en_spectra;
191  TH1F* rnBookkeeper = new TH1F("hrn", "random number bookkeeper", 100, -5, 5);
192 
193  /// build a table of knob switches
194  map<rwgt::ReweightLabel_t, bool> knob_switch = build_knob_switch(klist);
195 
196  /// store the nominal universe as the first element
197  for(int i = 0; i < nuniverses; i++)
198  {
199  if(i == 0) reco_en_spectra.push_back(new Spectrum(mcloader, kRecoEStandardBinning, total_selection, kNoShift));
200  else reco_en_spectra.push_back(create_universe(mcloader, rnBookkeeper, knob_switch));
201  //~ cout << i+1 << " universes created. ";
202  //~ cout << count_enabled_knobs(knob_switch) << " knobs enabled." << endl;
203  }
204 
205  /// load spectra
206  mcloader.Go();
207  cout << reco_en_spectra.size() << " universes created" << endl;
208 
209  /// store resulting histograms
210  TFile* toutf = new TFile("multiverse_histograms.root", "recreate");
211  THStack* hs = new THStack("hs", Form("%d universes, %s", reco_en_spectra.size(), knobstr.c_str()));
212  for(int i = reco_en_spectra.size()-1; i >= 0; i--)
213  {
214  TH1* hcur = reco_en_spectra[i]->ToTH1(reco_en_spectra[i]->POT());
215  if(i == 0)
216  {
217  hcur->SetLineColor(kRed);
218  hcur->SetMarkerColor(kRed);
219  }
220  else
221  {
222  hcur->SetLineColor(kGreen);
223  hcur->SetMarkerColor(kGreen);
224  }
225  hcur->SetName(Form("henu%d",i));
226  hcur->Write();
227  hs->Add(hcur);
228  }
229  toutf->Close();
230 
231  hs->Draw("nostack");
232  hs->GetHistogram()->GetXaxis()->SetTitle("Reco neutrino energy (GeV)");
233  gPad->Update();
234  gPad->Modified();
235  gROOT->ProcessLine(".! mkdir -p output");
236  gPad->Print(Form("output/%s_%duniverses.png", knobstr.c_str(), reco_en_spectra.size()));
237  new TCanvas();
238  rnBookkeeper->Draw("hist");
239 }
240 
241 #endif
const XML_Char * name
Definition: expat.h:151
enum BeamMode kRed
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
TFile * inf
Definition: drawXsec.C:1
set< int >::iterator it
unsigned int count_enabled_knobs(map< rwgt::ReweightLabel_t, bool > &table)
Spectrum * create_universe(SpectrumLoader &loader, TH1F *hrn, map< rwgt::ReweightLabel_t, bool > &kswitch)
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
std::vector< std::string > LoadFileList(const std::string &listfile)
Read list of input files from a text file, one per line.
Definition: Utilities.cxx:210
void multiverse_reweighting(int nfiles=INT_MAX, int nuniverses=10, string klist="knob_confs/maccqe_not_varied.conf", string flistpn="/nova/ana/users/slin/sam_definition_file_lists/prod_caf_R16-12-20-prod3recopreview.d_nd_genie_nonswap_fhc_nova_v08_full_v3.txt")
const XML_Char int const XML_Char * value
Definition: expat.h:331
virtual void Go() override
Load all the registered spectra.
loader
Definition: demo0.py:10
std::vector< float > Spectrum
Definition: Constants.h:610
const SystShifts kNoShift
Definition: SystShifts.cxx:21
OStream cout
Definition: OStream.cxx:6
std::vector< double > POT
static float min(const float a, const float b, const float c)
Definition: absgeo.cxx:45
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
map< string, rwgt::ReweightLabel_t > name2knob
const Cut total_selection
const int nuniverses
map< rwgt::ReweightLabel_t, bool > build_knob_switch(string fpn)
enum BeamMode kGreen
const Cut kAllNumuCCCuts