NuWROSysts.cxx
Go to the documentation of this file.
2 
3 #include "CAFAna/Core/Progress.h"
5 
7 
8 #include <cassert>
9 #include <iostream>
10 
11 #include "TFile.h"
12 #include "TH1.h"
13 #include "TTree.h"
14 
15 namespace ana
16 {
17  std::map<NuWROSyst::Key_t, TH1*> NuWROSyst::fgGenieHists;
18  std::map<NuWROSyst::Key_t, TH1*> NuWROSyst::fgNuWroHists;
20 
26 
27  //----------------------------------------------------------------------
29  : ISyst(ShortNameHelper(rwVar), LatexNameHelper(rwVar)),
30  fReweightVar(rwVar)
31  {
32  ++fgNInstances;
33  }
34 
35  //----------------------------------------------------------------------
37  {
38  --fgNInstances;
39  // Crashes at shutdown. Maybe because load_libs means everything winds up
40  // loaded twice? Leak isn't that big, just leave it.
41  /*
42  if(fgNInstances == 0){
43  // Only clean up when there are no instances left
44  for(auto it: fgGenieHists) delete it.second;
45  for(auto it: fgNuWroHists) delete it.second;
46  // In case a new instance is created
47  fgGenieHists.clear();
48  fgNuWroHists.clear();
49  }
50  if(fgNInstances < 0){
51  std::cerr << "NuWROSysts::fgNInstances is " << fgNInstances
52  << ". This should never happen" << std::endl;
53  }
54  */
55  }
56 
57  //----------------------------------------------------------------------
59  {
60  // Already called
61  if(!fgGenieHists.empty()) return;
62 
63  Progress prog("Loading NuWRO histograms");
64 
65  for(int pdg: {12, 14}){
66 
67  TFile f(TString::Format("$NOVA_ANA/caf/genie_vs_nuwro_%s.root", (pdg == 12) ? "nue" : "numu"));
68 
69  assert(!f.IsZombie());
70 
71  DontAddDirectory guard;
72 
73  for(bool genie: {false, true}){
74  TTree* tr = (TTree*)f.Get(genie ? "genie" : "nuwro");
75  assert(tr);
76  std::map<Key_t, TH1*>& hists = genie ? fgGenieHists : fgNuWroHists;
77 
78  for(int mode = 0; mode < 4; ++mode){
79  for(bool cc: {false, true}){
80  hists[{pdg, mode, cc, EReweightVar::kEnu }] = new TH1F("", "", 100, 0, 10);
81  hists[{pdg, mode, cc, EReweightVar::kX }] = new TH1F("", "", 100, 0, 5);
82  hists[{pdg, mode, cc, EReweightVar::kY }] = new TH1F("", "", 100, 0, 1);
83  hists[{pdg, mode, cc, EReweightVar::kQ2 }] = new TH1F("", "", 100, 0, 5);
84  hists[{pdg, mode, cc, EReweightVar::kMode}] = new TH1F("", "", 1, 0, 1);
85  }
86  }
87 
88  bool qel, res, dis, cc;
89  tr->SetBranchAddress("qel", &qel);
90  tr->SetBranchAddress("res", &res);
91  tr->SetBranchAddress("dis", &dis);
92  tr->SetBranchAddress("cc", &cc);
93 
94  double Ev, x, y, Q2;
95  tr->SetBranchAddress("Ev", &Ev);
96  tr->SetBranchAddress("x", &x);
97  tr->SetBranchAddress("y", &y);
98  tr->SetBranchAddress("Q2", &Q2);
99 
100  const unsigned int I = tr->GetEntries();
101  for(unsigned int i = 0; i < I; ++i){
102  tr->GetEntry(i);
103 
104  int mode = -1;
105  if(qel) mode = caf::kQE;
106  if(res) mode = caf::kRes;
107  if(dis) mode = caf::kDIS;
108  if(mode == -1) continue; // coherent?
109 
110  hists[{pdg, mode, cc, EReweightVar::kEnu }]->Fill(Ev);
111  hists[{pdg, mode, cc, EReweightVar::kX }]->Fill(x );
112  hists[{pdg, mode, cc, EReweightVar::kY }]->Fill(y );
113  hists[{pdg, mode, cc, EReweightVar::kQ2 }]->Fill(Q2);
114  hists[{pdg, mode, cc, EReweightVar::kMode}]->Fill(.5);
115 
116  if(i%10000 == 0) prog.SetProgress(i/(4.*I) +
117  (genie ? .25 : 0) +
118  ((pdg == 14) ? .5 : 0));
119  } // end for i
120  } // end for genie
121  } // end for pdg
122 
123  prog.Done();
124  }
125 
126  //----------------------------------------------------------------------
128  {
129  switch(v){
130  case EReweightVar::kEnu : return "nuwroEnu";
131  case EReweightVar::kX : return "nuwrox";
132  case EReweightVar::kY : return "nuwroy";
133  case EReweightVar::kQ2 : return "nuwroq2";
134  case EReweightVar::kMode: return "nuwroqmode";
135  }
136  assert(0 && "Unknown reweight");
137  }
138 
139  //----------------------------------------------------------------------
141  {
142  switch(v){
143  case EReweightVar::kEnu : return "NuWRO reweight on neutrino energy";
144  case EReweightVar::kX : return "NuWRO reweight on hadronic x";
145  case EReweightVar::kY : return "NuWRO reweight on hadronic y";
146  case EReweightVar::kQ2 : return "NuWRO reweight on Q^{2}";
147  case EReweightVar::kMode: return "NuWRO reweight on interaction mode";
148  }
149  assert(0 && "Unknown reweight");
150  }
151 
152  //----------------------------------------------------------------------
153  void NuWROSyst::Shift(double sigma,
154  caf::SRProxy* sr, double& weight) const
155  {
157 
158  if(sr->mc.nnu == 0) return;
159  //Leave weight unaltered
160  if(sigma == 0) return;
161 
162  const caf::SRNeutrinoProxy& nu = sr->mc.nu[0];
163 
164  if(nu.mode < caf::kQE || nu.mode > caf::kDIS) return; // probably COH
165 
166  if(abs(nu.pdg) != 12 && abs(nu.pdg) != 14) return; // don't reweight taus
167 
168  const Key_t key = {abs(nu.pdg), nu.mode, nu.iscc, fReweightVar};
169 
170  TH1* hgen = fgGenieHists[key];
171  TH1* hnuw = fgNuWroHists[key];
172 
173  assert(hgen && hnuw);
174 
175  double nuvar = -1;
176  if(fReweightVar == EReweightVar::kEnu ) nuvar = nu.E;
177  if(fReweightVar == EReweightVar::kX ) nuvar = nu.x;
178  if(fReweightVar == EReweightVar::kY ) nuvar = nu.y;
179  if(fReweightVar == EReweightVar::kQ2 ) nuvar = nu.q2;
180  if(fReweightVar == EReweightVar::kMode) nuvar = .5;
181  assert(nuvar != -1);
182 
183  // Off our axis range. Don't reweight
184  if(nuvar < hgen->GetXaxis()->GetXmin()) return;
185  if(nuvar > hgen->GetXaxis()->GetXmax()) return;
186 
187  const double wGenie = hgen->GetBinContent(hgen->FindBin(nuvar));
188  const double wNuWro = hnuw->GetBinContent(hnuw->FindBin(nuvar));
189 
190  // Either division by zero or 0/0. Don't reweight
191  if(wGenie == 0) return;
192 
193  // Do the actual reweight
194  weight *= 1 + (wNuWro/wGenie - 1)*sigma;
195 
196  if(weight < 0) weight = 0;
197  }
198 }
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
caf::Proxy< float > x
Definition: SRProxy.h:575
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:991
static void InitializeHistograms()
Definition: NuWROSysts.cxx:58
const Var weight
Proxy for caf::SRNeutrino.
Definition: SRProxy.h:510
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
static std::map< Key_t, TH1 * > fgNuWroHists
Definition: NuWROSysts.h:56
caf::Proxy< std::vector< caf::SRNeutrino > > nu
Definition: SRProxy.h:618
static int fgNInstances
Definition: NuWROSysts.h:57
void abs(TH1 *hist)
const NuWROSyst kNuWROReweightMode(NuWROSyst::EReweightVar::kMode)
Definition: NuWROSysts.h:64
caf::Proxy< int > mode
Definition: SRProxy.h:543
caf::Proxy< short int > nnu
Definition: SRProxy.h:617
caf::Proxy< float > E
Definition: SRProxy.h:520
NuWROSyst(EReweightVar rwVar)
Definition: NuWROSysts.cxx:28
TString hists[nhists]
Definition: bdt_com.C:3
void Shift(double sigma, caf::SRProxy *sr, double &weight) const override
+1 means match nuWRO, -1 moves in the opposite direction
Definition: NuWROSysts.cxx:153
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
const NuWROSyst kNuWROReweightQ2(NuWROSyst::EReweightVar::kQ2)
Definition: NuWROSysts.h:63
EReweightVar fReweightVar
Definition: NuWROSysts.h:38
std::string ShortNameHelper(EReweightVar v) const
Definition: NuWROSysts.cxx:127
const NuWROSyst kNuWROReweightY(NuWROSyst::EReweightVar::kY)
Definition: NuWROSysts.h:62
EReweightVar
The variable to reweight as a function of.
Definition: NuWROSysts.h:16
correl_xv GetXaxis() -> SetDecimals()
caf::Proxy< bool > iscc
Definition: SRProxy.h:538
caf::StandardRecord * sr
caf::Proxy< float > y
Definition: SRProxy.h:577
double sigma(TH1F *hist, double percentile)
caf::Proxy< caf::SRTruthBranch > mc
Definition: SRProxy.h:2138
void SetProgress(double frac)
Update the progress fraction between zero and one.
Definition: Progress.cxx:41
void cc()
Definition: test_ana.C:28
assert(nhit_max >=nhit_nbins)
A simple ascii-art progress bar.
Definition: Progress.h:9
caf::Proxy< float > q2
Definition: SRProxy.h:557
caf::Proxy< short int > pdg
Definition: SRProxy.h:552
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
const NuWROSyst kNuWROReweightEnu(NuWROSyst::EReweightVar::kEnu)
Definition: NuWROSysts.h:60
Prevent histograms being added to the current directory.
Definition: UtilsExt.h:46
static std::map< Key_t, TH1 * > fgGenieHists
Definition: NuWROSysts.h:56
void Done()
Call this when action is completed.
Definition: Progress.cxx:90
const NuWROSyst kNuWROReweightX(NuWROSyst::EReweightVar::kX)
Definition: NuWROSysts.h:61
std::string LatexNameHelper(EReweightVar v) const
Definition: NuWROSysts.cxx:140
enum BeamMode string