Systs.cxx
Go to the documentation of this file.
1 #include "CAFAna/Systs/Systs.h"
2 
3 #include "CAFAna/Cuts/BeamCuts.h"
5 
7 
8 #include "TFile.h"
9 #include "TH1.h"
10 
11 namespace ana
12 {
13  // Instantiate the systematics from the header
26 
27  const std::string kNueTASystDir = FindCAFAnaDir()+"/data/nd_systs";
28 
29  const std::string kNueTACalibFile = kNueTASystDir+"/test_calibup_fornue.root";
30  const std::string kNumuTACalibFile = kNueTASystDir+"/test_calibup_fornumu.root";
31 
32  const std::string kNueTACalibShapeFile = kNueTASystDir+"/test_calib_shape_fornue.root";
33  const std::string kNumuTACalibShapeFile = kNueTASystDir+"/test_calib_shape_fornumu.root";
34 
35  const std::string kNueTALightFile = kNueTASystDir+"/test_lightup_fornue.root";
36  const std::string kNumuTALightFile = kNueTASystDir+"/test_lightup_fornumu.root";
37 
38  const std::string kNueTACherFile = kNueTASystDir+"/test_cherenkov_fornue.root";
39  const std::string kNumuTACherFile = kNueTASystDir+"/test_cherenkov_fornumu.root";
40 
41  //the tau systematic description is in the Nue executive summary p. 25 of doc-26699-v9
42  const NuTruthSystComponentScale kTauScaleSyst("NuTauScale", "#nu_{#tau} Scale", kIsTau_NT && !kIsNC_NT, 0.6, NuTruthSystComponentScale::kLinear);
43 
44  //----------------------------------------------------------------------
45  void NormSyst::
46  Shift(double sigma, caf::SRProxy* /*sr*/, double& weight) const
47  {
48  weight *= 1+.10*sigma;
49  }
50 
51  //----------------------------------------------------------------------
52  void RelNormSyst::
53  Shift(double sigma, caf::SRProxy* sr, double& weight) const
54  {
55  weight *= 1 + .025*sigma // 5% difference between detectors
56  *( (sr->hdr.det==caf::kFARDET) ? 1 : -1 );
57  }
58 
59  //----------------------------------------------------------------------
60  void NCScaleSyst::
61  Shift(double sigma, caf::SRProxy* sr, double& weight) const
62  {
63  if(sr->mc.nnu == 0) return;
64  if(!sr->mc.nu[0].iscc) weight *= 1+.2*sigma;
65  }
66 
67  //----------------------------------------------------------------------
69  Shift(double sigma, caf::SRProxy* sr, double& weight) const
70  {
71  if(sr->mc.nnu == 0) return;
72  if(sr->mc.nu[0].iscc && abs(sr->mc.nu[0].pdgorig) == 12)
73  weight *= 1+.2*sigma;
74  }
75 
76  //----------------------------------------------------------------------
78  Shift(double sigma, caf::SRProxy* sr, double& weight) const
79  {
80  if(sr->mc.nnu == 0) return;
81  if(sr->mc.nu[0].iscc && abs(sr->mc.nu[0].pdgorig) == 14 && abs(sr->mc.nu[0].pdg) == 14)
82  weight *= 1+.2*sigma;
83  }
84 
85  //----------------------------------------------------------------------
86  void KaonScaleSyst::
87  Shift(double sigma, caf::SRProxy* sr, double& weight) const
88  {
89  if(sr->mc.nnu == 0) return;
90  if(kIsAllFromKaons(sr))
91  weight *= 1+.3*sigma;
92  }
93 
94  //----------------------------------------------------------------------
96  Shift(double sigma, caf::SRProxy* sr, double& weight) const
97  {
98  // 10% shift for 1 sigma
99  double alpha = sigma*0.1;
100  // This is the analytic formual of the ratio of two Poissons
101  // The hardcoded number is from calculate_avg_nslcs.C
102  weight *= pow(1 + alpha, sr->slc.nnonnoise) * exp( -3.91 * alpha);
103  }
104 
105  //----------------------------------------------------------------------
106  void LEMScaleSyst::
107  Shift(double sigma, caf::SRProxy* sr, double& weight) const
108  {
109  const double scale = 1 + 0.02*sigma; //shifts ~5% across boundary
110 
111  sr->sel.lem.pid *= scale;
112  }
113 
114  //----------------------------------------------------------------------
115  void RemIDScaleSyst::
116  Shift(double sigma, caf::SRProxy* sr, double& weight) const
117  {
118  const double scale = 1 + 0.07*sigma; //shifts ~5% across boundary
119 
120  sr->sel.remid.pid *= scale;
121  }
122 
123  //----------------------------------------------------------------------
124  void RockMuonNormSyst::
125  Shift(double sigma, caf::SRProxy* sr, double& weight) const
126  {
127  if (sr->hdr.det != 1) return; // Only reweight rock muons in the ND
128  if (sr->mc.nnu == 0) return;
129 
130  // No easy way to get this info in the CAFs, the proper thing to do is
131  // to add a isVtxContained flag in the truth branch of the CAFs, but I
132  if (fabs(sr->mc.nu[0].vtx.X()) > 200 ||
133  fabs(sr->mc.nu[0].vtx.Y()) > 200 ||
134  sr->mc.nu[0].vtx.Z() < 0 ||
135  sr->mc.nu[0].vtx.Z() > 1300)
136  weight *= 1 + 0.30*sigma;
137  }
138 
139  //----------------------------------------------------------------------
140  void ContainmentSyst::
141  Shift(double sigma, caf::SRProxy* sr, double& weight) const
142  {
143  const double dist = std::min( {sr->sel.contain.nplanestoback*6.,
144  sr->sel.contain.nplanestofront*6.,
145  sr->slc.ncellsfromedge*4., 0.} );
146 
147  weight *= exp(-dist/50) // 50cm characteristic distance
148  * .50 * sigma; // 50% reweighting at edge
149  }
150 
151  //----------------------------------------------------------------------
153  Shift(double sigma, caf::SRProxy* sr, double& weight) const
154  {
155  if(sr->mc.nnu == 0) return;
156  bool isDIS = sr->mc.nu[0].mode == 2;
157  double W2 = sr->mc.nu[0].W2;
158  if(isDIS && W2>4) weight *= 1+.15*sigma;
159  }
160 }
161 
162 
caf::Proxy< float > pid
Definition: SRProxy.h:928
Far Detector at Ash River.
Definition: SREnums.h:11
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
const LEMScaleSyst kLEMScaleSyst
Definition: Systs.cxx:21
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
void Shift(double sigma, caf::SRProxy *sr, double &weight) const override
Perform the systematic shift.
Definition: Systs.cxx:125
const KaonScaleSyst kKaonScaleSyst
Definition: Systs.cxx:19
const Var weight
Relative (uncorrelated FD to ND) normalization systematic.
Definition: Systs.h:36
caf::Proxy< caf::SRHeader > hdr
Definition: SRProxy.h:2119
const std::string kNumuTACalibShapeFile
Definition: Systs.cxx:33
caf::Proxy< caf::SRContain > contain
Definition: SRProxy.h:1230
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2108
void Shift(double sigma, caf::SRProxy *sr, double &weight) const override
Perform the systematic shift.
Definition: Systs.cxx:107
caf::Proxy< std::vector< caf::SRNeutrino > > nu
Definition: SRProxy.h:597
constexpr T pow(T x)
Definition: pow.h:75
const NuTruthSystComponentScale kTauScaleSyst("NuTauScale","#nu_{#tau} Scale", kIsTau_NT &&!kIsNC_NT, 0.6, NuTruthSystComponentScale::kLinear)
100% uncertainty scale on taus
Definition: Systs.h:176
const BeamNueScaleSyst kBeamNueScaleSyst
Definition: Systs.cxx:17
caf::Proxy< float > pid
Definition: SRProxy.h:1115
const DISHighWScaleSyst kDISHighWScaleSyst
Definition: Systs.cxx:25
void Shift(double sigma, caf::SRProxy *sr, double &weight) const override
Perform the systematic shift.
Definition: Systs.cxx:87
void abs(TH1 *hist)
LEM-shift systematic.
Definition: Systs.h:109
caf::Proxy< short int > nnu
Definition: SRProxy.h:596
void Shift(double sigma, caf::SRProxy *sr, double &weight) const override
Perform the systematic shift.
Definition: Systs.cxx:153
const std::string kNumuTALightFile
Definition: Systs.cxx:36
RemID-shift systematic.
Definition: Systs.h:122
const NuTruthCut kIsTau_NT([](const caf::SRNeutrinoProxy *truth){return(abs(truth->pdg)==16);})
Definition: TruthCuts.h:12
std::string FindCAFAnaDir()
Definition: Utilities.cxx:208
const std::string kNueTALightFile
Definition: Systs.cxx:35
Uncertainty in the scale of a single component of the spectrum.
const BeamIntensitySyst kBeamIntensitySyst
Definition: Systs.cxx:20
double dist
Definition: runWimpSim.h:113
void Shift(double sigma, caf::SRProxy *sr, double &weight) const override
Perform the systematic shift.
Definition: Systs.cxx:141
const NormSyst kNormSyst
Definition: Systs.cxx:14
void Shift(double sigma, caf::SRProxy *sr, double &weight) const override
Perform the systematic shift.
Definition: Systs.cxx:78
Double_t scale
Definition: plot.C:25
const ContainmentSyst kContainmentSyst
Definition: Systs.cxx:24
void Shift(double sigma, caf::SRProxy *sr, double &weight) const override
Perform the systematic shift.
Definition: Systs.cxx:53
Rock muon normalization - reweights all events with vtx outside det.
Definition: Systs.h:135
Absolute scale of beam &#39;s systematic.
Definition: Systs.h:60
caf::Proxy< int > nplanestofront
Definition: SRProxy.h:808
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:10
const std::string kNueTASystDir
Definition: Systs.cxx:27
caf::Proxy< unsigned int > ncellsfromedge
Definition: SRProxy.h:1292
void Shift(double sigma, caf::SRProxy *sr, double &weight) const override
Perform the systematic shift.
Definition: Systs.cxx:61
caf::Proxy< unsigned int > nnonnoise
Definition: SRProxy.h:1297
caf::StandardRecord * sr
caf::Proxy< int > nplanestoback
Definition: SRProxy.h:807
caf::Proxy< caf::SRRemid > remid
Definition: SRProxy.h:1248
Absolute Neutral Current scale systematic.
Definition: Systs.h:48
void Shift(double sigma, caf::SRProxy *sr, double &weight) const override
Perform the systematic shift.
Definition: Systs.cxx:116
double sigma(TH1F *hist, double percentile)
caf::Proxy< caf::SRLem > lem
Definition: SRProxy.h:1239
caf::Proxy< caf::SRTruthBranch > mc
Definition: SRProxy.h:2120
const NuTruthCut kIsNC_NT([](const caf::SRNeutrinoProxy *truth){return!truth->iscc;})
Definition: TruthCuts.h:15
const std::string kNueTACalibFile
Definition: Systs.cxx:29
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
const std::string kNueTACalibShapeFile
Definition: Systs.cxx:32
Absolute normalization systematic.
Definition: Systs.h:24
caf::Proxy< caf::SRSlice > slc
Definition: SRProxy.h:2124
void Shift(double sigma, caf::SRProxy *sr, double &weight) const override
Perform the systematic shift.
Definition: Systs.cxx:69
const std::string kNueTACherFile
Definition: Systs.cxx:38
void Shift(double sigma, caf::SRProxy *sr, double &weight) const override
Perform the systematic shift.
Definition: Systs.cxx:96
const NumuCCScaleSyst kNumuCCScaleSyst
Definition: Systs.cxx:18
const RelNormSyst kRelNormSyst
Definition: Systs.cxx:15
const std::string kNumuTACalibFile
Definition: Systs.cxx:30
caf::Proxy< caf::SRIDBranch > sel
Definition: SRProxy.h:2123
const RockMuonNormSyst kRockMuonNormSyst
Definition: Systs.cxx:23
T min(const caf::Proxy< T > &a, T b)
const Cut kIsAllFromKaons([](const caf::SRProxy *sr){if(sr->mc.nnu==0) return false;assert(sr->mc.nnu==1);return( (abs(sr->mc.nu[0].beam.tptype)==310)|| (abs(sr->mc.nu[0].beam.tptype)==130)|| (abs(sr->mc.nu[0].beam.tptype)==311)|| (abs(sr->mc.nu[0].beam.tptype)==321));})
Definition: BeamCuts.h:30
Shift beam intensity - reweights nslcs dist to fake a different intensity.
Definition: Systs.h:96
const NCScaleSyst kNCScaleSyst
Definition: Systs.cxx:16
const std::string kNumuTACherFile
Definition: Systs.cxx:39
void Shift(double sigma, caf::SRProxy *sr, double &weight) const override
Perform the systematic shift.
Definition: Systs.cxx:46
caf::Proxy< caf::Det_t > det
Definition: SRProxy.h:231
const RemIDScaleSyst kRemIDScaleSyst
Definition: Systs.cxx:22