MECSystsNux.h
Go to the documentation of this file.
1 #pragma once
2 
5 #include <TObjString.h>
7 
8 namespace ana
9 {
10 
11 
12 //===============================================================================================================
13 // Valencia MEC double 2D Gaussian Enhancement
14 //===============================================================================================================
20 
21 double CalcMECDoubleGaussEnhNux( const double q0, const double q3, const MECDoubleGaussEnhParamNux shift_param, const double shift_sigma)
22 {
23 // params of prod5 files fit
24 // double norm = 14.9925;
25 // double mean_q0 = 0.365515;
26 // double mean_q3 = 0.888951;
27 // double sigma_q0 = 0.0905725;
28 // double sigma_q3 = 0.340095;
29 // double corr = 0.808604;
30 // double norm_2 = 47.4552;
31 // double mean_q0_2 = 0.0355624;
32 // double mean_q3_2 = 0.444703;
33 // double sigma_q0_2 = 0.0587598;
34 // double sigma_q3_2 = 0.269792;
35 // double corr_2 = 0.706661;
36 // double baseline = 0.13100;
37 
38 
39  //params of MC to MC fit, meant to be the nominal...
40 
41  double norm = 15.0104;
42  double mean_q0 = 0.365683;
43  double mean_q3 = 0.89007;
44  double sigma_q0 = 0.090725;
45  double sigma_q3 = 0.340702;
46  double corr = 0.808929;
47  double norm_2 = 47.477;
48  double mean_q0_2 = 0.0356445;
49  double mean_q3_2 = 0.444196;
50  double sigma_q0_2 = 0.0588446;
51  double sigma_q3_2 = 0.270499;
52  double corr_2 = 0.709473;
53  double baseline = 0.131061;
54 
55 
56  switch( shift_param )
57  {
58  case kGauss2DNorm_1Nux:
59  norm += shift_sigma * 0.5;
60  if ( norm < 0 ) norm = 0.0;
61  break;
63  mean_q0 += shift_sigma * 0.025;
64  if ( mean_q0 < 0 ) mean_q0 = 0.0;
65  break;
67  mean_q3 += shift_sigma * 0.05;
68  if ( mean_q3 < 0 ) mean_q0 = 0.0;
69  break;
71  sigma_q0 += shift_sigma * 0.015;
72  if ( sigma_q0 < 0.01 ) sigma_q0 = 0.01;
73  break;
75  sigma_q3 += shift_sigma * 0.02;
76  if ( sigma_q3 < 0.01 ) sigma_q3 = 0.01;
77  break;
78  case kGauss2DCorr_1Nux :
79  corr += shift_sigma * 0.2;
80  if ( corr < -0.99 ) corr = -0.99;
81  else if ( corr > 0.99 ) corr = 0.99;
82  break;
83  case kGauss2DNorm_2Nux :
84  norm_2 += shift_sigma * 2;
85  if ( norm_2 < 0 ) norm_2 = 0.0;
86  break;
87  case kGauss2DMeanQ0_2Nux :
88  mean_q0_2 += shift_sigma * 0.025;
89  if ( mean_q0_2 < 0 ) mean_q0_2=0;
90  break;
91  case kGauss2DMeanQ3_2Nux :
92  mean_q3_2 += shift_sigma * 0.05;
93  if ( mean_q3_2 < 0 ) mean_q3_2 = 0;
94  break;
96  sigma_q0_2 += shift_sigma * 0.015;
97  if ( sigma_q0_2 < 0.01 ) sigma_q0_2 = 0.01;
98  break;
100  sigma_q3_2 += shift_sigma * 0.02;
101  if ( sigma_q3_2 < 0.01 ) sigma_q3_2 = 0.01;
102  break;
103  case kGauss2DCorr_2Nux :
104  corr_2 += shift_sigma * 0.2;
105  if ( corr_2 < -0.99 ) corr_2 = -0.99;
106  if ( corr_2 > 0.99 ) corr_2 = 0.99;
107  break;
108  case kBaselineNux :
109  baseline += shift_sigma*0.05;
110  //if (baseline < 0.6) baseline=0.6;
111 
112  }
113 // reminder: http://mathworld.wolfram.com/BivariateNormalDistribution.html
114  double z = pow( ( q0 - mean_q0 ) / sigma_q0, 2 ) + pow( ( q3 - mean_q3 ) / sigma_q3, 2 )
115  - 2 * corr * ( q0 - mean_q0 ) * ( q3 - mean_q3 ) / ( sigma_q0 * sigma_q3 );
116 
117  double z_2 = pow( ( q0 - mean_q0_2 ) / sigma_q0_2, 2 ) + pow( ( q3 - mean_q3_2 ) / sigma_q3_2, 2 )
118  - 2 * corr_2 * ( q0 - mean_q0_2 ) * ( q3 - mean_q3_2 ) / ( sigma_q0_2 * sigma_q3_2 );
119 
120  double weight= baseline + norm * exp( -0.5 * z / ( 1 - corr * corr ) ) + norm_2* exp( -0.5 * z_2 / ( 1 - corr_2 * corr_2 ) );
121  if (weight<0) weight=0;
122 // std::cout<<"The weight for this shift "<<shift_param<<" is "<<weight<<std::endl;
123  return weight;
124 
125 }
126 
127 const std::map< std::string , MECDoubleGaussEnhParamNux > DoubleGaussMapNux{
128 {"MECDoubleGaussEnhSystNorm_1Nux", kGauss2DNorm_1Nux},
129 {"MECDoubleGaussEnhSystMeanQ0_1Nux", kGauss2DMeanQ0_1Nux},
130 {"MECDoubleGaussEnhSystMeanQ3_1Nux", kGauss2DMeanQ3_1Nux},
131 {"MECDoubleGaussEnhSystSigmaQ0_1Nux", kGauss2DSigmaQ0_1Nux},
132 {"MECDoubleGaussEnhSystSigmaQ3_1Nux", kGauss2DSigmaQ3_1Nux },
133 {"MECDoubleGaussEnhSystCorr_1Nux", kGauss2DCorr_1Nux },
134 {"MECDoubleGaussEnhSystNorm_2Nux", kGauss2DNorm_2Nux },
135 {"MECDoubleGaussEnhSystMeanQ0_2Nux", kGauss2DMeanQ0_2Nux},
136 {"MECDoubleGaussEnhSystMeanQ3_2Nux", kGauss2DMeanQ3_2Nux},
137 {"MECDoubleGaussEnhSystSigmaQ0_2Nux", kGauss2DSigmaQ0_2Nux},
138 {"MECDoubleGaussEnhSystSigmaQ3_2Nux", kGauss2DSigmaQ3_2Nux},
139 {"MECDoubleGaussEnhSystCorr_2Nux", kGauss2DCorr_2Nux },
140 {"MECDoubleGaussEnhSystBaselineNux", kBaselineNux}
141 };
142 
143 double CalcMECDoubleGaussEnhShiftedParamNux( std::string shift_param, double shift_sigma)
144 {
145 
146  MECDoubleGaussEnhParamNux shifted_param = DoubleGaussMapNux.find(shift_param)->second;
147 
148 // params of prod5 files fit
149 // double norm = 14.9925;
150 // double mean_q0 = 0.365515;
151 // double mean_q3 = 0.888951;
152 // double sigma_q0 = 0.0905725;
153 // double sigma_q3 = 0.340095;
154 // double corr = 0.808604;
155 // double norm_2 = 47.4552;
156 // double mean_q0_2 = 0.0355624;
157 // double mean_q3_2 = 0.444703;
158 // double sigma_q0_2 = 0.0587598;
159 // double sigma_q3_2 = 0.269792;
160 // double corr_2 = 0.706661;
161 // double baseline = 0.13100;
162 
163 
164 
165  //params of MC to MC fit, meant to be the nominal...
166 
167  double norm = 15.0104;
168  double mean_q0 = 0.365683;
169  double mean_q3 = 0.89007;
170  double sigma_q0 = 0.090725;
171  double sigma_q3 = 0.340702;
172  double corr = 0.808929;
173  double norm_2 = 47.477;
174  double mean_q0_2 = 0.0356445;
175  double mean_q3_2 = 0.444196;
176  double sigma_q0_2 = 0.0588446;
177  double sigma_q3_2 = 0.270499;
178  double corr_2 = 0.709473;
179  double baseline = 0.131061;
180 
181  double param = 0.0;
182 
183  switch( shifted_param )
184  {
185  case kGauss2DNorm_1Nux:
186  norm += shift_sigma * 0.5; //.5
187  if ( norm < 0 ) norm = 0.0;
188  param = norm;
189  break;
190  case kGauss2DMeanQ0_1Nux:
191  mean_q0 += shift_sigma * 0.025;
192  if ( mean_q0 < 0 ) mean_q0 = 0.0;
193  param = mean_q0;
194  break;
195  case kGauss2DMeanQ3_1Nux:
196  mean_q3 += shift_sigma * 0.05;
197  if ( mean_q3 < 0 ) mean_q3 = 0.0;
198  param = mean_q3;
199  break;
201  sigma_q0 += shift_sigma * 0.015;
202  if ( sigma_q0 < 0.001 ) sigma_q0 = 0.001;
203  param =sigma_q0;
204  break;
206  sigma_q3 += shift_sigma * 0.02;
207  if ( sigma_q3 < 0.001 ) sigma_q3 = 0.001;
208  param = sigma_q3;
209  break;
210  case kGauss2DCorr_1Nux :
211  corr += shift_sigma * 0.2;
212  if ( corr < -0.99 ) corr = -0.99;
213  else if ( corr > 0.99 ) corr = 0.99;
214  param=corr;
215  break;
216  case kGauss2DNorm_2Nux :
217  norm_2 += shift_sigma * 2;
218  if ( norm_2 < 0 ) norm_2 = 0.0;
219  param = norm_2;
220  break;
221  case kGauss2DMeanQ0_2Nux :
222  mean_q0_2 += shift_sigma * 0.025;
223  if ( mean_q0_2 < 0 ) mean_q0_2 = 0.0;
224  param = mean_q0_2;
225  break;
226  case kGauss2DMeanQ3_2Nux :
227  mean_q3_2 += shift_sigma * 0.05;
228  if ( mean_q3_2 < 0 ) mean_q3_2 = 0.0;
229  param = mean_q3_2;
230  break;
231  case kGauss2DSigmaQ0_2Nux :
232  sigma_q0_2 += shift_sigma * 0.015;
233  if ( sigma_q0_2 < 0.001 ) sigma_q0_2 = 0.001;
234  param = sigma_q0_2;
235  break;
236  case kGauss2DSigmaQ3_2Nux :
237  sigma_q3_2 += shift_sigma * 0.02;
238  if ( sigma_q3_2 < 0.001 ) sigma_q3_2 = 0.001;
239  param = sigma_q3_2;
240  break;
241  case kGauss2DCorr_2Nux :
242  corr_2 += shift_sigma * 0.2;
243  if ( corr_2 < -0.99 ) corr_2 = -0.99;
244  if ( corr_2 > 0.99 ) corr_2 = 0.99;
245  param = corr_2;
246  break;
247  case kBaselineNux :
248  baseline += shift_sigma*0.05;
249  // if (baseline < 0.6) baseline=0.6;
250  param = baseline;
251  }
252 
253  return param;
254 }
255 
256 
257 const Var kMECDoubleGaussEnhNux( []( const caf::SRProxy* sr )
258 {
259  if ( !kIsNumuCC( sr ) || sr->mc.nu[0].mode != caf::kMEC ) return 1.0;
260  double q0 = kTrueQ0( sr );
261  double q3 = kTrueQ3( sr );
262  return CalcMECDoubleGaussEnhNux( q0, q3, kGauss2DNorm_1Nux, 0 ); // Shift = 0 gives nominal suppression
263 });
264 
265 const Var kMECDoubleGaussEnhNux_One( []( const caf::SRProxy* sr )
266 {
267  if ( !kIsNumuCC( sr ) || sr->mc.nu[0].mode != caf::kMEC ) return 1.0;
268  double q0 = kTrueQ0( sr );
269  double q3 = kTrueQ3( sr );
270  return CalcMECDoubleGaussEnhNux( q0, q3, kGauss2DNorm_1Nux, 1 ); // Shift = 0 gives nominal suppression
271 });
272 
274 {
275  public :
276  MECDoubleGaussEnhSystNux( const MECDoubleGaussEnhParamNux& shift_param, const std::string& shift_param_name )
277  : ISyst( "MECDoubleGaussEnhSystNux" + shift_param_name, "MEC 2D Gauss Syst " + shift_param_name ),
278  fShiftParam( shift_param )
279  {}
280 
281  void Shift( double sigma, caf::SRProxy* sr, double& weight ) const override
282  {
283  //This will need to be expanded when we have NC MEC Events in the new processing.
284  if( (kIsNumuCC( sr ) || kIsNue( sr)) && sr->mc.nu[0].mode == caf::kMEC ){
285  double q0 = kTrueQ0( sr );
286  double q3 = kTrueQ3( sr );
287 
288  double wgt_nominal = CalcMECDoubleGaussEnhNux( q0, q3, fShiftParam, 0 );
289  double wgt_shift = CalcMECDoubleGaussEnhNux( q0, q3, fShiftParam, sigma );
290 
291  weight *= wgt_shift / wgt_nominal;
292  }
293  else return;
294  }
295 
296  private:
298 };
299 
300 
301 }
302 
const Cut kIsNumuCC(CCFlavSel(14, 14))
Select CC .
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
MECDoubleGaussEnhParamNux
Definition: MECSystsNux.h:15
const std::map< std::string, MECDoubleGaussEnhParamNux > DoubleGaussMapNux
Definition: MECSystsNux.h:127
const Var weight
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2125
caf::Proxy< std::vector< caf::SRNeutrino > > nu
Definition: SRProxy.h:617
constexpr T pow(T x)
Definition: pow.h:75
double corr(TGraph *g, double thresh)
const Var kTrueQ0
Definition: TruthVars.h:32
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
double CalcMECDoubleGaussEnhNux(const double q0, const double q3, const MECDoubleGaussEnhParamNux shift_param, const double shift_sigma)
Definition: MECSystsNux.h:21
const Var kTrueQ3
Definition: TruthVars.h:38
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:10
const Cut kIsNue([](const caf::SRProxy *sr){return(sr->mc.nnu==1 &&abs(sr->mc.nu[0].pdg)==12);})
Definition: TruthCuts.h:9
caf::StandardRecord * sr
z
Definition: test.py:28
double CalcMECDoubleGaussEnhShiftedParamNux(std::string shift_param, double shift_sigma)
Definition: MECSystsNux.h:143
double sigma(TH1F *hist, double percentile)
caf::Proxy< caf::SRTruthBranch > mc
Definition: SRProxy.h:2137
const MECDoubleGaussEnhParamNux fShiftParam
Definition: MECSystsNux.h:297
const Var kMECDoubleGaussEnhNux_One([](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 CalcMECDoubleGaussEnhNux(q0, q3, kGauss2DNorm_1Nux, 1);})
Float_t norm
const Var kMECDoubleGaussEnhNux([](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 CalcMECDoubleGaussEnhNux(q0, q3, kGauss2DNorm_1Nux, 0);})
void Shift(double sigma, caf::SRProxy *sr, double &weight) const override
Perform the systematic shift.
Definition: MECSystsNux.h:281
MECDoubleGaussEnhSystNux(const MECDoubleGaussEnhParamNux &shift_param, const std::string &shift_param_name)
Definition: MECSystsNux.h:276
enum BeamMode string