NueAcceptSysts.cxx
Go to the documentation of this file.
5 
7 
9 
10 #include "TFile.h"
11 #include "TH1.h"
12 #include "TKey.h"
13 #include "TMath.h"
14 
15 #include <iostream>
16 
17 namespace ana
18 {
19  /// FHC
23  const NueAcceptSystSignalKin2020FHC kNueAcceptSystSignalKinPtExtrap2020FHC(ana::kExtrapPt, "accept_signalkin_pTextrap_FHC_2020");
24 
25 
26  /// RHC
30  const NueAcceptSystSignalKin2020RHC kNueAcceptSystSignalKinPtExtrap2020RHC(ana::kExtrapPt, "accept_signalkin_pTextrap_RHC_2020");
31 
32 
33  /// FHC implementation
35  {
36  // Nue acceptance systs of bkg from docdb-27935 ND subcomponents extrapolation
37  // maximal core / peripheral FHC (RHC) difference: +- 1.5% / +-1.2% (+- 4.1% / +- 2.3%)
38 
39  /// FHC values
40  const double core = .015;
41  const double peri = .012;
42 
43  // Check if FHC, otherwise left unaltered
44  if(kIsRHC(sr)) return;
45 
46  // Check if in FD, otherwise left unaltered
47  if(sr->hdr.det != caf::kFARDET) return;
48 
49  // Check if neutrino, but not numu->nue signal, otherwise left unaltered
50  if(sr->mc.nnu == 0) return;
51  if(abs(sr->mc.nu[0].pdg) == 12 && abs(sr->mc.nu[0].pdgorig) == 14 && sr->mc.nu[0].iscc) return;
52  // Left unaltered if tau neutrino
53  if(abs(sr->mc.nu[0].pdg) == 16 && sr->mc.nu[0].iscc) return;
54 
55  if(sigma == 0) return;
56 
57  double kFDBkgCorr = core;
58 
59  // Check whether in peripheral or core sample
60  if(kNue2018FDPeripheral(sr))
61  kFDBkgCorr = peri;
62 
63  weight *= 1+kFDBkgCorr*sigma;
64  }
65 
67  {
68  // Check if FHC, otherwise left unaltered
69  if(kIsRHC(sr)) return;
70 
71  // Check if in FD, otherwise left unaltered
72  if(sr->hdr.det != caf::kFARDET) return;
73 
74  // Check if signal or WS
75  if(sr->mc.nnu == 0) return;
76  if(
77  abs(sr->mc.nu[0].pdg) != 12 ||
78  abs(sr->mc.nu[0].pdgorig) != 14 ||
79  !(sr->mc.nu[0].iscc)
80  ) return;
81 
82  if(!fWeightHistFHC){
83 
84  // FD in file from docdb-27935 ExtrapSyst technote
85  // Can be remade with 3FlavorAna/Ana2018/ExtrapSysts/ scripts
86  // Using only correction of the weight with max impact, i.e. Cos weight for FHC and trueQ2 weight for RHC, see technote
87 
88  /// FHC
89  const std::string kWeightsFnameFHC = FindCAFAnaDir()+"/nue/Ana2018/AcceptSysts/FD_KinematicsCorrection_FHC.root";
90  const std::string kWeightNameFHC = "Cos_FD_KinematicsCorrection_FHC";
91 
92  TFile weightsFileFHC (kWeightsFnameFHC.c_str(),"read");
93 
94  if(weightsFileFHC.IsZombie()){
95  std::cerr << "Warning: couldn't open " << kWeightsFnameFHC << std::endl;
96  abort();
97  }
98 
99  fWeightHistFHC = (TH1*) weightsFileFHC.Get(kWeightNameFHC.c_str())->Clone();
100 
101  // disassociate the histogram from files
102  fWeightHistFHC -> SetDirectory(0);
103 
104  weightsFileFHC.Close();
105  }
106 
107  if (kNue2018FDAllSamples(sr))
108  {
109  double kAnaBin = kNue2018AnaBin(sr);
110  double kFDSignalWeight = 1;
111 
112  if(fWeightHistFHC->GetBinContent(kAnaBin) != 0)
113  kFDSignalWeight = fWeightHistFHC->GetBinContent(kAnaBin);
114 
115  weight *= 1+(kFDSignalWeight-1)*sigma;
116  }
117  }
118 
120  {
121  // Check if FHC, otherwise left unaltered
122  if(kIsRHC(sr)) return;
123 
124  // Check if in FD, otherwise left unaltered
125  if(sr->hdr.det != caf::kFARDET) return;
126 
127  // Check if signal or WS
128  if(sr->mc.nnu == 0) return;
129  if(
130  abs(sr->mc.nu[0].pdg) != 12 ||
131  abs(sr->mc.nu[0].pdgorig) != 14 ||
132  !(sr->mc.nu[0].iscc)
133  ) return;
134 
135  if(!fWeightHistFHC){
136 
137  // FD in file from docdb-43917 ExtrapSyst technote
138  // Can be remade with 3FlavorAna/Ana2020/AcceptSysts/ scripts
139  // Using only correction of the weight with max impact, which is PtP weight for FHC with pT extrap and trueQ2 without pT extrap, see technote
140 
141  /// FHC
142  std::string kWeightsFnameFHC;
143  std::string kWeightNameFHC;
144 
145  if(extrapVar == ana::kExtrapEmpty){
146  kWeightsFnameFHC = FindCAFAnaDir()+"/data/3flavor/FD_KinematicsCorrection_FHC.root";
147  kWeightNameFHC = "trueQ2_FD_KinematicsCorrection_FHC";
148  }
149  else if(extrapVar == ana::kExtrapPt){
150  kWeightsFnameFHC = FindCAFAnaDir()+"/data/3flavor/FD_KinematicsCorrection_pTExtrap_FHC.root";
151  kWeightNameFHC = "PtP_FD_KinematicsCorrection_FHC";
152  }
153  else{
154  std::cerr << "Error: Only pT extrap or no extrap options are available for acceptance syst!" << std::endl;
155  abort();
156  }
157 
158  TFile weightsFileFHC (kWeightsFnameFHC.c_str(),"read");
159 
160  if(weightsFileFHC.IsZombie()){
161  std::cerr << "Warning: couldn't open " << kWeightsFnameFHC << std::endl;
162  abort();
163  }
164 
165  fWeightHistFHC = (TH1*) weightsFileFHC.Get(kWeightNameFHC.c_str())->Clone();
166 
167  // disassociate the histogram from files
168  fWeightHistFHC -> SetDirectory(0);
169 
170  weightsFileFHC.Close();
171  }
172 
173  if (kNue2020FDAllSamples(sr))
174  {
175  double kAnaBin = kNue2020AnaBin(sr) +1 ; //Need to add +1 to get the right bin# for GetBinContent unlike kNue2020Axis where we Fill the histogram
176  double kFDSignalWeight = 1;
177 
178  if(fWeightHistFHC->GetBinContent(kAnaBin) != 0)
179  kFDSignalWeight = fWeightHistFHC->GetBinContent(kAnaBin);
180 
181  weight *= 1+(kFDSignalWeight-1)*sigma;
182  }
183  }
184 
185  /// RHC implementation
187  {
188  // Nue acceptance systs of bkg from docdb-27935 ND subcomponents extrapolation
189  // maximal core / peripheral FHC (RHC) difference: +- 1.5% / +-1.2% (+- 4.1% / +- 2.3%)
190 
191  /// RHC values
192  const double core = .041;
193  const double peri = .023;
194 
195  // Check if RHC, otherwise left unaltered
196  if(!kIsRHC(sr)) return;
197 
198  // Check if in FD, otherwise left unaltered
199  if(sr->hdr.det != caf::kFARDET) return;
200 
201  // Check if neutrino, but not numu->nue signal, otherwise left unaltered
202  if(sr->mc.nnu == 0) return;
203  if(abs(sr->mc.nu[0].pdg) == 12 && abs(sr->mc.nu[0].pdgorig) == 14 && sr->mc.nu[0].iscc) return;
204  // Left unaltered if tau neutrino
205  if(abs(sr->mc.nu[0].pdg) == 16 && sr->mc.nu[0].iscc) return;
206 
207  if(sigma == 0) return;
208 
209  double kFDBkgCorr = core;
210 
211  // Check whether in peripheral or core sample
212  if(kNue2018FDPeripheral(sr))
213  kFDBkgCorr = peri;
214 
215  weight *= 1+kFDBkgCorr*sigma;
216  }
217 
219  {
220  // Check if RHC, otherwise left unaltered
221  if(!kIsRHC(sr)) return;
222 
223  // Check if in FD, otherwise left unaltered
224  if(sr->hdr.det != caf::kFARDET) return;
225 
226  // Check if signal or WS
227  if(sr->mc.nnu == 0) return;
228  if(
229  abs(sr->mc.nu[0].pdg) != 12 ||
230  abs(sr->mc.nu[0].pdgorig) != 14 ||
231  !(sr->mc.nu[0].iscc)
232  ) return;
233 
234  if(!fWeightHistRHC){
235 
236  // FD in file from docdb-27935 ExtrapSyst technote
237  // Can be remade with CAFAna/nue/Ana2018/ExtrapSysts/ scripts
238  // Using only correction of the weight with max impact, i.e. Cos weight for FHC and trueQ2 weight for RHC, see technote
239 
240  /// RHC
241  const std::string kWeightsFnameRHC = FindCAFAnaDir()+"/nue/Ana2018/AcceptSysts/FD_KinematicsCorrection_RHC.root";
242  const std::string kWeightNameRHC = "trueQ2_FD_KinematicsCorrection_RHC";
243 
244  TFile weightsFileRHC (kWeightsFnameRHC.c_str(),"read");
245 
246  if(weightsFileRHC.IsZombie()){
247  std::cerr << "Warning: couldn't open " << kWeightsFnameRHC << std::endl;
248  abort();
249  }
250 
251  fWeightHistRHC = (TH1*) weightsFileRHC.Get(kWeightNameRHC.c_str())->Clone();
252 
253  //disassociate the histogram from files
254  fWeightHistRHC -> SetDirectory(0);
255 
256  weightsFileRHC.Close();
257  }
258 
259  if (kNue2018FDAllSamples(sr))
260  {
261  double kAnaBin = kNue2018AnaBin(sr);
262  double kFDSignalWeight = 1;
263 
264  if(fWeightHistRHC->GetBinContent(kAnaBin) != 0)
265  kFDSignalWeight = fWeightHistRHC->GetBinContent(kAnaBin);
266 
267  weight *= 1+(kFDSignalWeight-1)*sigma;
268  }
269  }
270 
272  {
273  // Check if RHC, otherwise left unaltered
274  if(!kIsRHC(sr)) return;
275 
276  // Check if in FD, otherwise left unaltered
277  if(sr->hdr.det != caf::kFARDET) return;
278 
279  // Check if signal or WS
280  if(sr->mc.nnu == 0) return;
281  if(
282  abs(sr->mc.nu[0].pdg) != 12 ||
283  abs(sr->mc.nu[0].pdgorig) != 14 ||
284  !(sr->mc.nu[0].iscc)
285  ) return;
286 
287  if(!fWeightHistRHC){
288 
289  // FD in file from docdb-43917 ExtrapSyst technote
290  // Can be remade with 3FlavorAna/Ana2020/AcceptSysts/ scripts
291  // Using only correction of the weight with max impact, which is pT/p for RHC without pT extrap and cos for RHC with pT extrap, see technote
292 
293  /// RHC
294  std::string kWeightsFnameRHC;
295  std::string kWeightNameRHC;
296 
297  if(extrapVar == ana::kExtrapEmpty){
298  kWeightsFnameRHC = FindCAFAnaDir()+"/data/3flavor/FD_KinematicsCorrection_RHC.root";
299  kWeightNameRHC = "PtP_FD_KinematicsCorrection_RHC";
300  }
301  else if(extrapVar == ana::kExtrapPt){
302  kWeightsFnameRHC = FindCAFAnaDir()+"/data/3flavor/FD_KinematicsCorrection_pTExtrap_RHC.root";
303  kWeightNameRHC = "Cos_FD_KinematicsCorrection_RHC";
304  }
305  else{
306  std::cerr << "Error: Only pT extrap and no extrap options are available for acceptance syst!" << std::endl;
307  abort();
308  }
309 
310  TFile weightsFileRHC (kWeightsFnameRHC.c_str(),"read");
311 
312  if(weightsFileRHC.IsZombie()){
313  std::cerr << "Warning: couldn't open " << kWeightsFnameRHC << std::endl;
314  abort();
315  }
316 
317  fWeightHistRHC = (TH1*) weightsFileRHC.Get(kWeightNameRHC.c_str())->Clone();
318 
319  //disassociate the histogram from files
320  fWeightHistRHC -> SetDirectory(0);
321 
322  weightsFileRHC.Close();
323  }
324 
325  if (kNue2020FDAllSamples(sr))
326  {
327  double kAnaBin = kNue2020AnaBin(sr) + 1; //Need to add +1 to get the right bin# for GetBinContent unlike kNue2020Axis where we Fill the histogram
328  double kFDSignalWeight = 1;
329 
330  if(fWeightHistRHC->GetBinContent(kAnaBin) != 0)
331  kFDSignalWeight = fWeightHistRHC->GetBinContent(kAnaBin);
332 
333  weight *= 1+(kFDSignalWeight-1)*sigma;
334  }
335  }
336 }
Far Detector at Ash River.
Definition: SREnums.h:11
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
const NueAcceptSystSignalKin2018RHC kNueAcceptSystSignalKin2018RHC
const Var weight
caf::Proxy< caf::SRHeader > hdr
Definition: SRProxy.h:2137
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
caf::Proxy< std::vector< caf::SRNeutrino > > nu
Definition: SRProxy.h:618
void Shift(double sigma, caf::SRProxy *sr, double &weight) const override
FHC implementation.
OStream cerr
Definition: OStream.cxx:7
void abs(TH1 *hist)
caf::Proxy< short int > nnu
Definition: SRProxy.h:617
void Shift(double sigma, caf::SRProxy *sr, double &weight) const override
RHC implementation.
void Shift(double sigma, caf::SRProxy *sr, double &weight) const override
Perform the systematic shift.
const NueAcceptSystSignalKin2020FHC kNueAcceptSystSignalKin2020FHC
std::string FindCAFAnaDir()
Definition: Utilities.cxx:204
const Var kNue2018AnaBin([](const caf::SRProxy *sr){int selBin=kNue2018SelectionBin(sr);float nuE=kNueEnergy2018(sr);int nuEBin=nuE/0.5;assert(nuEBin<=8 &&"An event with nuE > 4.5 should never happen");int anaBin=9 *selBin+nuEBin;return anaBin;})
Use this Analysis Binning for Ana2018, official Binning.
Definition: NueCuts2018.h:171
void Shift(double sigma, caf::SRProxy *sr, double &weight) const override
Perform the systematic shift.
const NueAcceptSystBkg2018RHC kNueAcceptSystBkg2018RHC
RHC.
void Shift(double sigma, caf::SRProxy *sr, double &weight) const override
Perform the systematic shift.
const NueAcceptSystSignalKin2018FHC kNueAcceptSystSignalKin2018FHC
caf::StandardRecord * sr
const NueAcceptSystSignalKin2020FHC kNueAcceptSystSignalKinPtExtrap2020FHC(ana::kExtrapPt,"accept_signalkin_pTextrap_FHC_2020")
const NueAcceptSystBkg2018FHC kNueAcceptSystBkg2018FHC
FHC.
const Cut kNue2018FDPeripheral(kNue2018FDPeripheralFunc)
double sigma(TH1F *hist, double percentile)
caf::Proxy< caf::SRTruthBranch > mc
Definition: SRProxy.h:2138
const Cut kIsRHC([](const caf::SRProxy *sr){return sr->spill.isRHC;})
Definition: Vars.h:16
const Var kNue2020AnaBin([](const caf::SRProxy *sr){int selBin=kNue2020SelectionBin(sr);float nuE=kNueEnergy2020(sr);int nuEBin=nuE/0.5;assert(nuEBin<=8 &&"An event with nuE > 4.5 should never happen");int anaBin=9 *selBin+nuEBin;return anaBin;})
Use this Analysis Binning for Ana2020, official Binning.
Definition: NueCuts2020.h:191
const Cut kNue2018FDAllSamples
Definition: NueCuts2018.h:77
const Cut kNue2020FDAllSamples
Definition: NueCuts2020.h:84
const NueAcceptSystSignalKin2020RHC kNueAcceptSystSignalKinPtExtrap2020RHC(ana::kExtrapPt,"accept_signalkin_pTextrap_RHC_2020")
void Shift(double sigma, caf::SRProxy *sr, double &weight) const override
Perform the systematic shift.
caf::Proxy< caf::Det_t > det
Definition: SRProxy.h:232
enum BeamMode string
const NueAcceptSystSignalKin2020RHC kNueAcceptSystSignalKin2020RHC