GeniePCASyst.cxx
Go to the documentation of this file.
2 
4 
5 // The following is to allow the CMake build
6 // to find the necessary files
7 #ifdef NOVACMAKE
8  #include "cetlib/search_path.h"
9  #include "cetlib_except/exception.h"
10 #endif
11 
12 // 2020 cuts
15 // 2018 cuts
18 
19 #include "CAFAna/Cuts/TruthCuts.h"
20 
22 
23 #include "TFile.h"
24 #include "TH1.h"
25 #include "TKey.h"
26 
27 #include <iostream>
28 
29 namespace ana
30 {
31  // 2020 versions
32  /*
33  GeniePCASyst* GetGeniePrincipals2020(int PCIdx)
34  {
35  if(PCIdx >= 100){
36  std::cerr << "Error : Component Index too large!" << std::endl;
37  abort();
38  }
39 
40  static GeniePCASyst* Cache[100] = {0};
41  std::string PCIdxStr = TString::Format("%02d", PCIdx).Data();
42  if(!Cache[PCIdx])
43  Cache[PCIdx] = new GeniePCASyst((FindCAFAnaDir()+"/data/xs/genie_all_pc_shifts_fn_2020.root"),
44  "genie_all_pc"+PCIdxStr, "Genie PC "+std::to_string(PCIdx));
45 
46  return Cache[PCIdx];
47  }
48  */
49 
51  {
52  if(PCIdx >= 100){
53  std::cerr << "Error : Component Index too large!" << std::endl;
54  abort();
55  }
56 
57  static GeniePCASyst* Cache[100] = {0};
58  std::string PCIdxStr = TString::Format("%02d", PCIdx).Data();
59  if(!Cache[PCIdx]){
60 #ifdef NOVACMAKE
61  std::string filePath;
62 
63  cet::search_path sp("FW_SEARCH_PATH");
64  if( !sp.find_file("CAFAna/data/xs/genie_small_pc_shifts_fn_2020.root", filePath) )
65  throw cet::exception("CAFToEventList")
66  << "Cannot find pca file ";
67  Cache[PCIdx] = new GeniePCASyst(filePath,
68  "genie_small_pc" + PCIdxStr, "Genie PC " + std::to_string(PCIdx));
69 #else
70  Cache[PCIdx] = new GeniePCASyst((FindCAFAnaDir() + "/data/xs/genie_small_pc_shifts_fn_2020.root"),
71  "genie_small_pc" + PCIdxStr, "Genie PC " + std::to_string(PCIdx));
72 #endif
73 
74  }
75  return Cache[PCIdx];
76  }
77 
78  // 2018 versions
80  {
81  if(PCIdx >= 100){
82  std::cerr << "Error : Component Index too large!" << std::endl;
83  abort();
84  }
85 
86  static GeniePCASyst* Cache[100] = {0};
87  std::string PCIdxStr = TString::Format("%02d", PCIdx).Data();
88  if(!Cache[PCIdx])
89  Cache[PCIdx] = new GeniePCASyst((FindCAFAnaDir()+"/data/xs/genie_all_pc_shifts_fn_2018.root"),
90  "genie_all_pc"+PCIdxStr, "Genie PC "+std::to_string(PCIdx));
91 
92  return Cache[PCIdx];
93  }
94 
96  {
97  if(PCIdx >= 100){
98  std::cerr << "Error : Component Index too large!" << std::endl;
99  abort();
100  }
101 
102  static GeniePCASyst* Cache[100] = {0};
103  std::string PCIdxStr = TString::Format("%02d", PCIdx).Data();
104  if(!Cache[PCIdx])
105  Cache[PCIdx] = new GeniePCASyst((FindCAFAnaDir()+"/data/xs/genie_small_pc_shifts_fn_2018.root"),
106  "genie_small_pc"+PCIdxStr, "Genie PC "+std::to_string(PCIdx));
107 
108  return Cache[PCIdx];
109  }
110 
112  const std::string& shortname,
113  const std::string& latexname)
114  : ISyst(shortname,latexname),
115  fFileName(fname), fHistName(shortname), fHistos()
116  {
117  }
118 
120  {
121  for(int i = 0; i < kNumFluxType; ++i){
122  for(int j = 0; j < kNumDets; ++j){
123  for(int k = 0; k < kNumSels; ++k){
124  for(int l = 0; l < kNumFlavors; ++l){
125  for(int m = 0; m < kNumSigns; ++m){
126  delete fHistos[i][j][k][l][m];
127  fHistos[i][j][k][l][m] = 0;
128  }
129  }
130  }
131  }
132  }
133  }
134 
136  {
137  if(fHistos[0][0][0][0][0]) return;
138 
139  TFile fin(fFileName.c_str(),"READ");
140  if(fin.IsZombie()){
141  std::cerr << "Warning: couldn't open " << fFileName << std::endl;
142  return;
143  }
144 
145  TIter iterHist(gDirectory->GetListOfKeys());
146  TKey* keyHist;
147  int foundHisto=0;
148 
149  while((keyHist = (TKey*)iterHist())) {
150  TString histName = keyHist->GetName(); // Get a histogram name
151  if(histName.Contains(fHistName)){ //is it the desired syst
152  foundHisto++;
153 
154  int flux=0;
155  if(histName.Contains("FHC")) flux = kFHC;
156  if(histName.Contains("RHC")) flux = kRHC;
157 
158  int det = 0;
159  if(histName.Contains("ND")) det = kND;
160  if(histName.Contains("FD")) det = kFD;
161 
162  int sel = 0;
163  if(histName.Contains("Numu")) sel = kNumuSel;
164  if(histName.Contains("Nue" )) sel = kNueSel;
165 
166  int flav = 0;
167  if(histName.Contains("SigQE" )) flav = kSigQE;
168  if(histName.Contains("SigNonQE")) flav = kSigNonQE;
169  if(histName.Contains("Other" )) flav = kBkgFlav;
170  if(histName.Contains("NC" )) flav = kNC;
171 
172  int sign = 0;
173  if(histName.Contains("minus")) sign = kMinus;
174  if(histName.Contains("plus" )) sign = kPlus;
175 
176  //store relevant histograms
177  fHistos[flux][det][sel][flav][sign] = (TH1D*) fin.Get(histName)->Clone();
178  // disassociate it from the file it came from
179  fHistos[flux][det][sel][flav][sign]->SetDirectory(0);
180  }
181  }
182 
183  if (foundHisto==0) {
184  std::cerr << "Genie PC systematic histogram " << fHistName
185  << " is not in file " << fFileName <<"; abort!" << std::endl;
186  abort();
187  }
188 
189  fin.Close();
190  }
191 
193  ESign sign) const
194  {
196 
197  if(sr->mc.nnu == 0) return 1;
198 
199  int flux = sr->spill.isRHC ? kRHC : kFHC;
200 
201  int det = 0;
202  switch (sr->hdr.det){
203  case 1 : det = kND; break;
204  case 2 : det = kFD; break;
205  default: std::cerr <<"Wrong detector; ignore" << std::endl; return 1;
206  };
207 
208  // 2018 remnants
209  //const Cut NumuNDSel = kNumuCutND2018;
210  //const Cut NueNDSel = kNue2018NDCVNSsb;
211  //const Cut NumuFDSel = kNumuCutFD2018;
212  //const Cut NueFDSel = kNue2018FDAllSamples;
213 
214  // 2020 version
215  const Cut NumuNDSel = kNumu2020ND;
216  const Cut NueNDSel = kNue2020ND;
217  const Cut NumuFDSel = kNumu2020FD;
218  const Cut NueFDSel = kNue2020FD;
219 
220  int sel = -1;
221  if(det == kND){
222  if(NumuNDSel(sr)) sel = kNumuSel;
223  if(NueNDSel(sr)) sel = kNueSel;
224  }
225  if(det == kFD){
226  if(NumuFDSel(sr)) sel = kNumuSel;
227  if(NueFDSel(sr)) sel = kNueSel;
228  }
229  if(sel == -1) return 1;
230 
231  int flav = -1;
232  if(sel == kNumuSel){
233  if(kIsNumu(sr) && !kIsNC(sr) && kIsQE(sr)) flav = kSigQE;
234  if(kIsNumu(sr) && !kIsNC(sr) && !kIsQE(sr)) flav = kSigNonQE;
235  if(kIsNC(sr)) flav = kNC;
236  }
237  if(sel == kNueSel){
238  if(kIsNue(sr) && !kIsNC(sr) && kIsQE(sr)) flav = kSigQE;
239  if(kIsNue(sr) && !kIsNC(sr) && !kIsQE(sr)) flav = kSigNonQE;
240  if(kIsNumu(sr) && !kIsNC(sr)) flav = kBkgFlav;
241  if(kIsNC(sr)) flav = kNC;
242  }
243  if(flav == -1) return 1;
244 
245 
246  TH1D* h = fHistos[flux][det][sel][flav][sign];
247  if (h == 0){
248  std::cerr << fHistName+": Can't find desired histogram; ignore"<< std:: endl;
249  return 1;
250  }
251 
252  double energy = sr->mc.nu[0].E;
253 
254  if (energy > h->GetXaxis()->GetXmin() &&
255  energy < h->GetXaxis()->GetXmax() ){
256  return h->Interpolate(energy);
257  }
258 
259  if (energy > h->GetXaxis()->GetXmax())
260  return h->GetBinContent(h->GetNbinsX());
261 
262  return 1;
263  }
264 
266  caf::SRProxy* sr, double& weight) const
267  {
268  if(sigma == 0) return;
269 
270  // Need different histogram for plus or minus sigma
271  ESign sign = kPlus;
272  if(sigma<0) sign = kMinus;
273 
274  weight *=1+(GetWeight(sr, sign)-1.)*TMath::Abs(sigma);
275  weight = std::max(0., weight);
276  }
277 
278  TH1D* GeniePCASyst::getHist(int f,int d,int s,int fl, int p)
279  {
281 
282  return fHistos[f][d][s][fl][p];
283  }
284 }
const Cut kIsQE
Definition: TruthCuts.cxx:104
caf::Proxy< caf::SRSpill > spill
Definition: SRProxy.h:2143
void InitializeHistograms() const
T max(const caf::Proxy< T > &a, T b)
TString fin
Definition: Style.C:24
std::string fHistName
Definition: GeniePCASyst.h:38
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
const Var weight
GeniePCASyst * GetGeniePrincipals2018(int PCIdx)
caf::Proxy< caf::SRHeader > hdr
Definition: SRProxy.h:2137
const Cut kIsNumu([](const caf::SRProxy *sr){return(sr->mc.nnu==1 &&abs(sr->mc.nu[0].pdg)==14);})
Definition: TruthCuts.h:10
const char * p
Definition: xmltok.h:285
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
caf::Proxy< std::vector< caf::SRNeutrino > > nu
Definition: SRProxy.h:618
const Cut kNumu2020ND
Definition: NumuCuts2020.h:57
string fFileName
Definition: nue_pid_effs.C:43
OStream cerr
Definition: OStream.cxx:7
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
caf::Proxy< short int > nnu
Definition: SRProxy.h:617
std::string find_file(std::string const &filename) const
std::string FindCAFAnaDir()
Definition: Utilities.cxx:203
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
const Cut kNue2020ND
Definition: NueCuts2020.h:178
Loaders::FluxType flux
const XML_Char * s
Definition: expat.h:262
TH1D * fHistos[kNumFluxType][kNumDets][kNumSels][kNumFlavors][kNumSigns]
Definition: GeniePCASyst.h:40
const Cut kIsNC([](const caf::SRProxy *sr){if(sr->mc.nnu==0) return false;assert(sr->mc.nnu==1);return!sr->mc.nu[0].iscc;})
Is this a Neutral Current event?
Definition: TruthCuts.h:8
double energy
Definition: plottest35.C:25
const Cut kIsNue([](const caf::SRProxy *sr){return(sr->mc.nnu==1 &&abs(sr->mc.nu[0].pdg)==12);})
Definition: TruthCuts.h:9
Float_t d
Definition: plot.C:236
caf::StandardRecord * sr
const double j
Definition: BetheBloch.cxx:29
const Cut kNumu2020FD
Definition: NumuCuts2020.h:59
double GetWeight(const caf::SRProxy *sr, ESign sign) const
const Cut kNue2020FD
Definition: NueCuts2020.h:65
double sigma(TH1F *hist, double percentile)
GeniePCASyst * GetGeniePrincipals2018Small(int PCIdx)
caf::Proxy< caf::SRTruthBranch > mc
Definition: SRProxy.h:2138
std::string fFileName
Definition: GeniePCASyst.h:37
GeniePCASyst(const std::string &fname, const std::string &shortname, const std::string &latexname)
void Shift(double sigma, caf::SRProxy *sr, double &weight) const override
Perform the systematic shift.
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
caf::Proxy< bool > isRHC
Definition: SRProxy.h:1376
std::string to_string(ModuleType const mt)
Definition: ModuleType.h:34
def sign(x)
Definition: canMan.py:197
TH1D * getHist(int f, int d, int s, int fl, int p)
caf::Proxy< caf::Det_t > det
Definition: SRProxy.h:232
GeniePCASyst * GetGeniePrincipals2020Small(int PCIdx)
enum BeamMode string