MECTuningUtils2020.h
Go to the documentation of this file.
1 #pragma once
2 
3 #include "CAFAna/Core/Loaders.h"
4 #include <TObjString.h>
5 
6 namespace jw
7 {
8  class IRescaledSigmaSyst : public ana::ISyst
9  {
10  public:
11  IRescaledSigmaSyst(const std::string & shortName, const std::string & latexName, double sigmaScale=1.0)
12  : ana::ISyst(shortName, latexName), fSigmaScale(sigmaScale) {};
13  double GetSigmaScale() const { return fSigmaScale; };
14  void SetSigmaScale(double sc) { fSigmaScale = sc; };
15  protected:
16  double fSigmaScale;
17  };
18 
19  class CompNormSyst : public IRescaledSigmaSyst
20  {
21  public:
22  CompNormSyst(const ana::Cut & selCut, double sigmaScale=1.0)
23  : fID(sInstanceCount++),
24  IRescaledSigmaSyst("CompNormShift_" + std::to_string(sInstanceCount), "Component Normalization Shift " + std::to_string(sInstanceCount), sigmaScale),
25  fSelCut(selCut)
26  {}
27  void Shift(double sigma, caf::SRProxy* sr, double& weight) const override;
28 
29  private:
30  const ana::Cut fSelCut;
31  unsigned int fID;
32 
33  static unsigned int sInstanceCount;
34  };
35  unsigned int CompNormSyst::sInstanceCount = 0;
36  void CompNormSyst::Shift(double sigma, caf::SRProxy* sr, double& weight) const
37  {
38  if (!fSelCut(sr))
39  return;
40 
41  double wgt = 1 + fSigmaScale * sigma;
42  if (wgt < 0)
43  wgt = 0;
44  weight *= wgt;
45  }
46 
47 }
48 
49 //===============================================================================================================
50 // Valencia MEC 2D Gaussian Enhancement
51 //===============================================================================================================
52 
53 namespace ana
54 {
55 
57 
58 double CalcMECGaussEnh( const double q0, const double q3, const MECGaussEnhParam shift_param, const double shift_sigma )
59 {
60 
61  double norm = 10.0;
62  double mean_q0 = 0.25;
63  double mean_q3 = 0.50;
64  double sigma_q0 = 0.07;
65  double sigma_q3 = 0.13;
66  double corr = 0.0;
67 
68  switch( shift_param )
69  {
70  case kGauss2DNorm :
71  norm += shift_sigma * 2.0;
72  if ( norm < 0 ) norm = 0.0;
73  break;
74  case kGauss2DMeanQ0 :
75  mean_q0 += shift_sigma * 0.05;
76  break;
77  case kGauss2DMeanQ3 :
78  mean_q3 += shift_sigma * 0.1;
79  break;
80  case kGauss2DSigmaQ0 :
81  sigma_q0 += shift_sigma * 0.015;
82  if ( sigma_q0 < 0.01 ) sigma_q0 = 0.01;
83  break;
84  case kGauss2DSigmaQ3 :
85  sigma_q3 += shift_sigma * 0.02;
86  if ( sigma_q3 < 0.01 ) sigma_q3 = 0.01;
87  break;
88  case kGauss2DCorr :
89  corr += shift_sigma * 0.2;
90  if ( corr < -0.99 ) corr = -0.99;
91  else if ( corr > 0.99 ) corr = 0.99;
92  }
93 // reminder: http://mathworld.wolfram.com/BivariateNormalDistribution.html
94  double z = pow( ( q0 - mean_q0 ) / sigma_q0, 2 ) + pow( ( q3 - mean_q3 ) / sigma_q3, 2 )
95  - 2 * corr * ( q0 - mean_q0 ) * ( q3 - mean_q3 ) / ( sigma_q0 * sigma_q3 );
96 
97  return 1.0 + norm * exp( -0.5 * z / ( 1 - corr * corr ) );
98 }
99 
100 class MECGaussEnhSyst : public ISyst
101 {
102  public :
103  MECGaussEnhSyst( const MECGaussEnhParam& shift_param, const std::string& shift_param_name )
104  : ISyst( "MECGaussEnhSyst" + shift_param_name, "MEC 2D Gauss Syst " + shift_param_name ),
105  fShiftParam( shift_param )
106  {}
107 
108  void Shift( double sigma, caf::SRProxy* sr, double& weight ) const override
109  {
110  //if ( !kIsNumuCC || !kIsDytmanMEC ) return;
111  if ( !kIsNumuCC( sr ) || sr->mc.nu[0].mode != caf::kMEC ) return; // This will work for Empirical MEC. Will it also work for Valencia MEC?
112 
113  double q0 = kTrueQ0( sr );
114  double q3 = kTrueQ3( sr );
115 
116  double wgt_nominal = CalcMECGaussEnh( q0, q3, fShiftParam, 0 );
117  double wgt_shift = CalcMECGaussEnh( q0, q3, fShiftParam, sigma );
118 
119  weight *= wgt_shift / wgt_nominal;
120  }
121 
122  private:
123  const MECGaussEnhParam fShiftParam;
124 };
125 
126 const Var kMECGaussEnh( []( const caf::SRProxy* sr )
127 {
128  if ( !kIsNumuCC( sr ) || sr->mc.nu[0].mode != caf::kMEC ) return 1.0; // This will work for Empirical MEC. Will it also work for Valencia MEC?
129  double q0 = kTrueQ0( sr );
130  double q3 = kTrueQ3( sr );
131  return CalcMECGaussEnh( q0, q3, kGauss2DNorm, 0 ); // Shift = 0 gives nominal suppression
132 });
133 
134 //===============================================================================================================
135 // MINOS Resonance Suppression
136 //===============================================================================================================
137 
139 
140 double CalcMinosResSupp( const double Q2, const MinosResSuppParam shift_param, const double shift_sigma )
141 {
142  double A = 1.010;
143  double Q0 = 0.156;
144 
145  switch( shift_param )
146  {
147  case kMinosResSuppNorm :
148  A += shift_sigma * 0.1;
149  if ( A < 0 ) A = 0.0;
150  break;
151  case kMinosResSuppQ0 :
152  Q0 += shift_sigma * 0.03;
153  if ( Q0 < 0.01 ) Q0 = 0.01;
154  }
155 
156  double supp = A / ( 1 + exp( 1 - sqrt( Q2 ) / Q0 ) );
157  if ( supp > 1 ) supp = 1.0;
158 
159  return supp;
160 }
161 
162 class MinosResSuppSyst : public ISyst
163 {
164  public :
165  MinosResSuppSyst( const MinosResSuppParam& shift_param, const std::string& shift_param_name )
166  : ISyst( "MinosResSuppSyst" + shift_param_name, "MINOS RES Supp Syst " + shift_param_name ),
167  fShiftParam( shift_param )
168  {}
169 
170  void Shift( double sigma, caf::SRProxy* sr, double& weight ) const override
171  {
172  if ( !kIsNumuCC( sr ) || sr->mc.nu[0].mode != kIsRes( sr ) || sr->mc.nu[0].tgtA == 1 ) return; // Exclude hydrogen
173 
174  double Q2 = kTrueQ2( sr );
175 
176  double wgt_nominal = CalcMinosResSupp( Q2, fShiftParam, 0 );
177  double wgt_shift = CalcMinosResSupp( Q2, fShiftParam, sigma );
178 
179  weight *= wgt_shift / wgt_nominal;
180  }
181 
182  private:
183  const MinosResSuppParam fShiftParam;
184 };
185 
186 }
MinosResSuppSyst(const MinosResSuppParam &shift_param, const std::string &shift_param_name)
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:991
const Var weight
CompNormSyst(const ana::Cut &selCut, double sigmaScale=1.0)
const Cut kIsRes
Definition: TruthCuts.cxx:111
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2038
T sqrt(T number)
Definition: d0nt_math.hpp:156
caf::Proxy< std::vector< caf::SRNeutrino > > nu
Definition: SRProxy.h:570
constexpr T pow(T x)
Definition: pow.h:75
virtual void Shift(double sigma, caf::SRProxy *sr, double &weight) const
Perform the systematic shift.
Definition: ISyst.cxx:26
double corr(TGraph *g, double thresh)
const Var kTrueQ0
Definition: TruthVars.h:32
const Var kMECGaussEnh([](const caf::SRProxy *sr){if(!kIsNumuCC(sr)||sr->mc.nu[0].mode!=caf::kMEC) return 1.0;double q0=kTrueQ0(sr);double q3=kTrueQ3(sr);return CalcMECGaussEnh(q0, q3, kGauss2DNorm, 0);})
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
const Var kTrueQ2
Definition: TruthVars.h:27
MECGaussEnhParam
double CalcMECGaussEnh(const double q0, const double q3, const MECGaussEnhParam shift_param, const double shift_sigma)
void SetSigmaScale(double sc)
const Var kTrueQ3
Definition: TruthVars.h:38
double CalcMinosResSupp(const double Q2, const MinosResSuppParam shift_param, const double shift_sigma)
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:10
const ana::Var wgt
z
Definition: test.py:28
double sigma(TH1F *hist, double percentile)
MECGaussEnhSyst(const MECGaussEnhParam &shift_param, const std::string &shift_param_name)
caf::Proxy< caf::SRTruthBranch > mc
Definition: SRProxy.h:2050
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
static const double A
Definition: Units.h:82
void Shift(double sigma, caf::SRProxy *sr, double &weight) const override
Perform the systematic shift.
Definition: mec_tuning.C:207
Float_t norm
void Shift(double sigma, caf::SRProxy *sr, double &weight) const override
Perform the systematic shift.
string shortName
THUMBNAIL BLOCK: We need to make a thumbnail for each.
double GetSigmaScale() const
void Shift(double sigma, caf::SRProxy *sr, double &weight) const override
Perform the systematic shift.
ISyst(const std::string &shortName, const std::string &latexName)
Definition: ISyst.cxx:10
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
static unsigned int sInstanceCount
Definition: mec_tuning.C:204
MinosResSuppParam
IRescaledSigmaSyst(const std::string &shortName, const std::string &latexName, double sigmaScale=1.0)
static constexpr Double_t sr
Definition: Munits.h:164