event_reweighting_all_knobs.C
Go to the documentation of this file.
1 #ifdef __CINT__
2 void event_reweighting_all_knobs(int nfiles = INT_MAX, int knobStart = 0, int knobEnd = 1, 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")
3 {
4  std::cout << "Sorry, you must run in compiled mode" << std::endl;
5 }
6 #else
7 
9 #include "CAFAna/Core/Spectrum.h"
10 #include "CAFAna/Core/Utilities.h"
11 #include "CAFAna/Cuts/TruthCuts.h"
12 #include "CAFAna/Systs/Systs.h"
14 #include "TCanvas.h"
15 #include "THStack.h"
16 
17 using namespace ana;
18 
19 map<int, Spectrum*> ShiftOneKnob(int kIdx, SpectrumLoader& loader)
20 {
21  map<int, Spectrum*> rwgtSpectrum;
22 
23  rwgtSpectrum[0] = new Spectrum(loader, kRecoE_Custom_Axis, kIsNumuCC);
24  rwgtSpectrum[1] = new Spectrum(loader, kRecoE_Custom_Axis, kIsNumuCC, SystShifts(GetGenieRwgtSyst((rwgt::ReweightLabel_t)kIdx), +1));
25  rwgtSpectrum[-1] = new Spectrum(loader, kRecoE_Custom_Axis, kIsNumuCC, SystShifts(GetGenieRwgtSyst((rwgt::ReweightLabel_t)kIdx), -1));
26 
27  return rwgtSpectrum;
28 }
29 
30 void event_reweighting_all_knobs(int nfiles = INT_MAX, int knobStart = 0, int knobEnd = 1, 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")
31 {
32  //////////////////////////////////////////////////////////////////////////////
33  // GENIE knob to name map
34  map<rwgt::ReweightLabel_t, string> knobs2names;
35  knobs2names[rwgt::kReweightNull] = "kReweightNull";
36  knobs2names[rwgt::fReweightMaNCEL] = "fReweightMaNCEL";
37  knobs2names[rwgt::fReweightEtaNCEL] = "fReweightEtaNCEL";
38  knobs2names[rwgt::fReweightNormCCQE] = "fReweightNormCCQE";
39  knobs2names[rwgt::fReweightNormCCQEenu] = "fReweightNormCCQEenu";
40  knobs2names[rwgt::fReweightMaCCQEshape] = "fReweightMaCCQEshape";
41  knobs2names[rwgt::fReweightMaCCQE] = "fReweightMaCCQE";
42  knobs2names[rwgt::fReweightVecCCQEshape] = "fReweightVecCCQEshape";
43  knobs2names[rwgt::fReweightNormCCRES] = "fReweightNormCCRES";
44  knobs2names[rwgt::fReweightMaCCRESshape] = "fReweightMaCCRESshape";
45  knobs2names[rwgt::fReweightMvCCRESshape] = "fReweightMvCCRESshape";
46  knobs2names[rwgt::fReweightMaCCRES] = "fReweightMaCCRES";
47  knobs2names[rwgt::fReweightMvCCRES] = "fReweightMvCCRES";
48  knobs2names[rwgt::fReweightNormNCRES] = "fReweightNormNCRES";
49  knobs2names[rwgt::fReweightMaNCRESshape] = "fReweightMaNCRESshape";
50  knobs2names[rwgt::fReweightMvNCRESshape] = "fReweightMvNCRESshape";
51  knobs2names[rwgt::fReweightMaNCRES] = "fReweightMaNCRES";
52  knobs2names[rwgt::fReweightMvNCRES] = "fReweightMvNCRES";
53  knobs2names[rwgt::fReweightMaCOHpi] = "fReweightMaCOHpi";
54  knobs2names[rwgt::fReweightR0COHpi] = "fReweightR0COHpi";
55  knobs2names[rwgt::fReweightRvpCC1pi] = "fReweightRvpCC1pi";
56  knobs2names[rwgt::fReweightRvpCC2pi] = "fReweightRvpCC2pi";
57  knobs2names[rwgt::fReweightRvpNC1pi] = "fReweightRvpNC1pi";
58  knobs2names[rwgt::fReweightRvpNC2pi] = "fReweightRvpNC2pi";
59  knobs2names[rwgt::fReweightRvnCC1pi] = "fReweightRvnCC1pi";
60  knobs2names[rwgt::fReweightRvnCC2pi] = "fReweightRvnCC2pi";
61  knobs2names[rwgt::fReweightRvnNC1pi] = "fReweightRvnNC1pi";
62  knobs2names[rwgt::fReweightRvnNC2pi] = "fReweightRvnNC2pi";
63  knobs2names[rwgt::fReweightRvbarpCC1pi] = "fReweightRvbarpCC1pi";
64  knobs2names[rwgt::fReweightRvbarpCC2pi] = "fReweightRvbarpCC2pi";
65  knobs2names[rwgt::fReweightRvbarpNC1pi] = "fReweightRvbarpNC1pi";
66  knobs2names[rwgt::fReweightRvbarpNC2pi] = "fReweightRvbarpNC2pi";
67  knobs2names[rwgt::fReweightRvbarnCC1pi] = "fReweightRvbarnCC1pi";
68  knobs2names[rwgt::fReweightRvbarnCC2pi] = "fReweightRvbarnCC2pi";
69  knobs2names[rwgt::fReweightRvbarnNC1pi] = "fReweightRvbarnNC1pi";
70  knobs2names[rwgt::fReweightRvbarnNC2pi] = "fReweightRvbarnNC2pi";
71  knobs2names[rwgt::fReweightAhtBY] = "fReweightAhtBY";
72  knobs2names[rwgt::fReweightBhtBY] = "fReweightBhtBY";
73  knobs2names[rwgt::fReweightCV1uBY] = "fReweightCV1uBY";
74  knobs2names[rwgt::fReweightCV2uBY] = "fReweightCV2uBY";
75  knobs2names[rwgt::fReweightAhtBYshape] = "fReweightAhtBYshape";
76  knobs2names[rwgt::fReweightBhtBYshape] = "fReweightBhtBYshape";
77  knobs2names[rwgt::fReweightCV1uBYshape] = "fReweightCV1uBYshape";
78  knobs2names[rwgt::fReweightCV2uBYshape] = "fReweightCV2uBYshape";
79  knobs2names[rwgt::fReweightNormDISCC] = "fReweightNormDISCC";
80  knobs2names[rwgt::fReweightRnubarnuCC] = "fReweightRnubarnuCC";
81  knobs2names[rwgt::fReweightDISNuclMod] = "fReweightDISNuclMod";
82  knobs2names[rwgt::fReweightNC] = "fReweightNC";
83  knobs2names[rwgt::fReweightAGKY_xF1pi] = "fReweightAGKY_xF1pi";
84  knobs2names[rwgt::fReweightAGKY_pT1pi] = "fReweightAGKY_pT1pi";
85  knobs2names[rwgt::fReweightFormZone] = "fReweightFormZone";
86  knobs2names[rwgt::fReweightMFP_pi] = "fReweightMFP_pi";
87  knobs2names[rwgt::fReweightMFP_N] = "fReweightMFP_N";
88  knobs2names[rwgt::fReweightFrCEx_pi] = "fReweightFrCEx_pi";
89  knobs2names[rwgt::fReweightFrElas_pi] = "fReweightFrElas_pi";
90  knobs2names[rwgt::fReweightFrInel_pi] = "fReweightFrInel_pi";
91  knobs2names[rwgt::fReweightFrAbs_pi] = "fReweightFrAbs_pi";
92  knobs2names[rwgt::fReweightFrPiProd_pi] = "fReweightFrPiProd_pi";
93  knobs2names[rwgt::fReweightFrCEx_N] = "fReweightFrCEx_N";
94  knobs2names[rwgt::fReweightFrElas_N] = "fReweightFrElas_N";
95  knobs2names[rwgt::fReweightFrInel_N] = "fReweightFrInel_N";
96  knobs2names[rwgt::fReweightFrAbs_N] = "fReweightFrAbs_N";
97  knobs2names[rwgt::fReweightFrPiProd_N] = "fReweightFrPiProd_N";
98  knobs2names[rwgt::fReweightCCQEPauliSupViaKF] = "fReweightCCQEPauliSupViaKF";
99  knobs2names[rwgt::fReweightCCQEMomDistroFGtoSF] = "fReweightCCQEMomDistroFGtoSF";
100  knobs2names[rwgt::fReweightBR1gamma] = "fReweightBR1gamma";
101  knobs2names[rwgt::fReweightBR1eta] = "fReweightBR1eta";
102  knobs2names[rwgt::fReweightTheta_Delta2Npi] = "fReweightTheta_Delta2Npi";
103  knobs2names[rwgt::fReweightZNormCCQE] = "fReweightZNormCCQE";
104  knobs2names[rwgt::fReweightZExpA1CCQE] = "fReweightZExpA1CCQE";
105  knobs2names[rwgt::fReweightZExpA2CCQE] = "fReweightZExpA2CCQE";
106  knobs2names[rwgt::fReweightZExpA3CCQE] = "fReweightZExpA3CCQE";
107  knobs2names[rwgt::fReweightZExpA4CCQE] = "fReweightZExpA4CCQE";
108  knobs2names[rwgt::fReweightAxFFCCQEshape] = "fReweightAxFFCCQEshape";
109  //////////////////////////////////////////////////////////////////////////////
110 
111  /// check if knob index is in range
112  //~ if(knobStart < rwgt::kReweightNull || knobEnd < rwgt::kReweightNull ||
113  //~ knobStart > rwgt::fReweightAxFFCCQEshape || knobEnd > rwgt::fReweightAxFFCCQEshape)
114  //~ {
115  //~ cout << "GENIE knob index out of range." << endl;
116  //~ return;
117  //~ }
118  knobStart = max(knobStart, (int)rwgt::kReweightNull);
119  knobEnd = min(knobEnd, (int)rwgt::fReweightAxFFCCQEshape);
120 
121  /// load file names into a vector
122  vector<string> flist = LoadFileList(flistpn);
123  if(!flist.size())
124  {
125  cout << "Not able to load the file list. Quit." << endl;
126  return;
127  }
128  /// only run over the specified number of files
129  flist.resize(min(min(nfiles, INT_MAX), (int)flist.size()));
130 
131  /// spectrum loaders for MC
132  SpectrumLoader mcloader(flist);
133 
134  /// containers for results
135  map<int, map<int, Spectrum*> > rwgtSpectra;
136  map<int, THStack*> rwgtHists;
137 
138  /// call for all reweighting, one at a time
139  for(int knobIdx = knobStart; knobIdx <= knobEnd; knobIdx++)
140  rwgtSpectra[knobIdx] = ShiftOneKnob(knobIdx, mcloader);
141 
142  /// load spectra
143  mcloader.Go();
144 
145  /// iterate through the reweighted plots
146  TFile* toutf = new TFile("reweighted_event_spectra.root","recreate");
147  for (auto x: rwgtSpectra)
148  {
149  string cur_rwgt_name = knobs2names[(rwgt::ReweightLabel_t)x.first];
150  new TCanvas(Form("c%d",x.first),Form("c%d",x.first));
151  THStack* hs = new THStack(Form("hs%d",x.first), cur_rwgt_name.c_str());
152  TH1* nom = x.second[0]->ToTH1(x.second[0]->POT());
153  TH1* p1 = x.second[1]->ToTH1(x.second[1]->POT());
154  TH1* m1 = x.second[-1]->ToTH1(x.second[-1]->POT());
155  nom->SetLineColor(kBlack);
156  p1->SetLineColor(kBlue);
157  m1->SetLineColor(kRed);
158  nom->SetName("nom");
159  p1->SetName("p1");
160  m1->SetName("m1");
161  hs->Add(nom);
162  hs->Add(p1);
163  hs->Add(m1);
164  hs->Draw("nostack");
165  /// write to file
166  toutf->cd("/");
167  toutf->mkdir(cur_rwgt_name.c_str());
168  toutf->cd(cur_rwgt_name.c_str());
169  nom->Write();
170  p1->Write();
171  m1->Write();
172  }
173  toutf->Close();
174 }
175 
176 #endif
const Cut kIsNumuCC(CCFlavSel(14, 14))
Select CC .
map< int, Spectrum * > ShiftOneKnob(int kIdx, SpectrumLoader &loader)
enum BeamMode kRed
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
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 event_reweighting_all_knobs(int nfiles=INT_MAX, int knobStart=0, int knobEnd=1, 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")
virtual void Go() override
Load all the registered spectra.
loader
Definition: demo0.py:10
std::vector< float > Spectrum
Definition: Constants.h:610
OStream cout
Definition: OStream.cxx:6
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.
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
enum BeamMode kBlue