GeniePCASyst.cxx
Go to the documentation of this file.
2 
4 
5 // 2020 cuts
8 // 2018 cuts
11 
12 #include "CAFAna/Cuts/TruthCuts.h"
13 
15 
16 #include "TFile.h"
17 #include "TH1.h"
18 #include "TKey.h"
19 
20 namespace ana
21 {
22  // 2020 versions
23  /*
24  GeniePCASyst* GetGeniePrincipals2020(int PCIdx)
25  {
26  if(PCIdx >= 100){
27  std::cerr << "Error : Component Index too large!" << std::endl;
28  abort();
29  }
30 
31  static GeniePCASyst* Cache[100] = {0};
32  std::string PCIdxStr = TString::Format("%02d", PCIdx).Data();
33  if(!Cache[PCIdx])
34  Cache[PCIdx] = new GeniePCASyst((FindCAFAnaDir()+"/data/xs/genie_all_pc_shifts_fn_2020.root"),
35  "genie_all_pc"+PCIdxStr, "Genie PC "+std::to_string(PCIdx));
36 
37  return Cache[PCIdx];
38  }
39  */
40 
42  {
43  if(PCIdx >= 100){
44  std::cerr << "Error : Component Index too large!" << std::endl;
45  abort();
46  }
47 
48  static GeniePCASyst* Cache[100] = {0};
49  std::string PCIdxStr = TString::Format("%02d", PCIdx).Data();
50  if(!Cache[PCIdx])
51  Cache[PCIdx] = new GeniePCASyst((FindCAFAnaDir()+"/data/xs/genie_small_pc_shifts_fn_2020.root"),
52  "genie_small_pc"+PCIdxStr, "Genie PC "+std::to_string(PCIdx));
53 
54  return Cache[PCIdx];
55  }
56 
57  // 2018 versions
59  {
60  if(PCIdx >= 100){
61  std::cerr << "Error : Component Index too large!" << std::endl;
62  abort();
63  }
64 
65  static GeniePCASyst* Cache[100] = {0};
66  std::string PCIdxStr = TString::Format("%02d", PCIdx).Data();
67  if(!Cache[PCIdx])
68  Cache[PCIdx] = new GeniePCASyst((FindCAFAnaDir()+"/data/xs/genie_all_pc_shifts_fn_2018.root"),
69  "genie_all_pc"+PCIdxStr, "Genie PC "+std::to_string(PCIdx));
70 
71  return Cache[PCIdx];
72  }
73 
75  {
76  if(PCIdx >= 100){
77  std::cerr << "Error : Component Index too large!" << std::endl;
78  abort();
79  }
80 
81  static GeniePCASyst* Cache[100] = {0};
82  std::string PCIdxStr = TString::Format("%02d", PCIdx).Data();
83  if(!Cache[PCIdx])
84  Cache[PCIdx] = new GeniePCASyst((FindCAFAnaDir()+"/data/xs/genie_small_pc_shifts_fn_2018.root"),
85  "genie_small_pc"+PCIdxStr, "Genie PC "+std::to_string(PCIdx));
86 
87  return Cache[PCIdx];
88  }
89 
91  const std::string& shortname,
92  const std::string& latexname)
93  : ISyst(shortname,latexname),
94  fFileName(fname), fHistName(shortname), fHistos()
95  {
96  }
97 
99  {
100  for(int i = 0; i < kNumFluxType; ++i){
101  for(int j = 0; j < kNumDets; ++j){
102  for(int k = 0; k < kNumSels; ++k){
103  for(int l = 0; l < kNumFlavors; ++l){
104  for(int m = 0; m < kNumSigns; ++m){
105  delete fHistos[i][j][k][l][m];
106  fHistos[i][j][k][l][m] = 0;
107  }
108  }
109  }
110  }
111  }
112  }
113 
115  {
116  if(fHistos[0][0][0][0][0]) return;
117 
118  TFile fin(fFileName.c_str(),"READ");
119  if(fin.IsZombie()){
120  std::cerr << "Warning: couldn't open " << fFileName << std::endl;
121  return;
122  }
123 
124  TIter iterHist(gDirectory->GetListOfKeys());
125  TKey* keyHist;
126  int foundHisto=0;
127 
128  while((keyHist = (TKey*)iterHist())) {
129  TString histName = keyHist->GetName(); // Get a histogram name
130  if(histName.Contains(fHistName)){ //is it the desired syst
131  foundHisto++;
132 
133  int flux=0;
134  if(histName.Contains("FHC")) flux = kFHC;
135  if(histName.Contains("RHC")) flux = kRHC;
136 
137  int det = 0;
138  if(histName.Contains("ND")) det = kND;
139  if(histName.Contains("FD")) det = kFD;
140 
141  int sel = 0;
142  if(histName.Contains("Numu")) sel = kNumuSel;
143  if(histName.Contains("Nue" )) sel = kNueSel;
144 
145  int flav = 0;
146  if(histName.Contains("SigQE" )) flav = kSigQE;
147  if(histName.Contains("SigNonQE")) flav = kSigNonQE;
148  if(histName.Contains("Other" )) flav = kBkgFlav;
149  if(histName.Contains("NC" )) flav = kNC;
150 
151  int sign = 0;
152  if(histName.Contains("minus")) sign = kMinus;
153  if(histName.Contains("plus" )) sign = kPlus;
154 
155  //store relevant histograms
156  fHistos[flux][det][sel][flav][sign] = (TH1D*) fin.Get(histName)->Clone();
157  // disassociate it from the file it came from
158  fHistos[flux][det][sel][flav][sign]->SetDirectory(0);
159  }
160  }
161 
162  if (foundHisto==0) {
163  std::cerr << "Genie PC systematic histogram " << fHistName
164  << " is not in file " << fFileName <<"; abort!" << std::endl;
165  abort();
166  }
167 
168  fin.Close();
169  }
170 
172  ESign sign) const
173  {
175 
176  if(sr->mc.nnu == 0) return 1;
177 
178  int flux = sr->spill.isRHC ? kRHC : kFHC;
179 
180  int det = 0;
181  switch (sr->hdr.det){
182  case 1 : det = kND; break;
183  case 2 : det = kFD; break;
184  default: std::cerr <<"Wrong detector; ignore" << std::endl; return 1;
185  };
186 
187  // 2018 remnants
188  //const Cut NumuNDSel = kNumuCutND2018;
189  //const Cut NueNDSel = kNue2018NDCVNSsb;
190  //const Cut NumuFDSel = kNumuCutFD2018;
191  //const Cut NueFDSel = kNue2018FDAllSamples;
192 
193  // 2020 version
194  const Cut NumuNDSel = kNumu2020ND;
195  const Cut NueNDSel = kNue2020ND;
196  const Cut NumuFDSel = kNumu2020FD;
197  const Cut NueFDSel = kNue2020FD;
198 
199  int sel = -1;
200  if(det == kND){
201  if(NumuNDSel(sr)) sel = kNumuSel;
202  if(NueNDSel(sr)) sel = kNueSel;
203  }
204  if(det == kFD){
205  if(NumuFDSel(sr)) sel = kNumuSel;
206  if(NueFDSel(sr)) sel = kNueSel;
207  }
208  if(sel == -1) return 1;
209 
210  int flav = -1;
211  if(sel == kNumuSel){
212  if(kIsNumu(sr) && !kIsNC(sr) && kIsQE(sr)) flav = kSigQE;
213  if(kIsNumu(sr) && !kIsNC(sr) && !kIsQE(sr)) flav = kSigNonQE;
214  if(kIsNC(sr)) flav = kNC;
215  }
216  if(sel == kNueSel){
217  if(kIsNue(sr) && !kIsNC(sr) && kIsQE(sr)) flav = kSigQE;
218  if(kIsNue(sr) && !kIsNC(sr) && !kIsQE(sr)) flav = kSigNonQE;
219  if(kIsNumu(sr) && !kIsNC(sr)) flav = kBkgFlav;
220  if(kIsNC(sr)) flav = kNC;
221  }
222  if(flav == -1) return 1;
223 
224 
225  TH1D* h = fHistos[flux][det][sel][flav][sign];
226  if (h == 0){
227  std::cerr << fHistName+": Can't find desired histogram; ignore"<< std:: endl;
228  return 1;
229  }
230 
231  double energy = sr->mc.nu[0].E;
232 
233  if (energy > h->GetXaxis()->GetXmin() &&
234  energy < h->GetXaxis()->GetXmax() ){
235  return h->Interpolate(energy);
236  }
237 
238  if (energy > h->GetXaxis()->GetXmax())
239  return h->GetBinContent(h->GetNbinsX());
240 
241  return 1;
242  }
243 
245  caf::SRProxy* sr, double& weight) const
246  {
247  if(sigma == 0) return;
248 
249  // Need different histogram for plus or minus sigma
250  ESign sign = kPlus;
251  if(sigma<0) sign = kMinus;
252 
253  weight *=1+(GetWeight(sr, sign)-1.)*TMath::Abs(sigma);
254  weight = std::max(0., weight);
255  }
256 
257  TH1D* GeniePCASyst::getHist(int f,int d,int s,int fl, int p)
258  {
260 
261  return fHistos[f][d][s][fl][p];
262  }
263 }
const Cut kIsQE
Definition: TruthCuts.cxx:104
caf::Proxy< caf::SRSpill > spill
Definition: SRProxy.h:2102
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:2096
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:2085
caf::Proxy< std::vector< caf::SRNeutrino > > nu
Definition: SRProxy.h:574
const Cut kNumu2020ND
Definition: NumuCuts2020.h:57
string fFileName
Definition: nue_pid_effs.C:43
OStream cerr
Definition: OStream.cxx:7
caf::Proxy< short int > nnu
Definition: SRProxy.h:573
std::string FindCAFAnaDir()
Definition: Utilities.cxx:208
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:2097
std::string fFileName
Definition: GeniePCASyst.h:37
GeniePCASyst(const std::string &fname, const std::string &shortname, const std::string &latexname)
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
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:1332
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
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:231
GeniePCASyst * GetGeniePrincipals2020Small(int PCIdx)