MECTuningUtils.h
Go to the documentation of this file.
1 #pragma once
2 
3 #include "CAFAna/Core/Loaders.h"
4 #include <TObjString.h>
7 
8 namespace ana
9 {
10 
11  const Var GetWeightFromShifts( const SystShifts& sys_shifts )
12  {
13  const Var kWeightFromShifts( [ sys_shifts ]( const caf::SRProxy *sr )
14  {
15  caf::SRProxy* face_palm = const_cast<caf::SRProxy*>( sr );
16  double weight = 1.0;
17  sys_shifts.Shift( face_palm, weight );
18  return weight;
19  });
20  return kWeightFromShifts;
21  }
22 
23  // this thing basically needs to be an IPrediction that gives
24  // the correct components of the spectrum.
25  // or at least have a PredictComponent() that can do
26  // (Flavors::kAll, Current::kBoth, Sign::kBoth)
27  // when asked, since that's the one that you need
28  // to get PredictSyst() to work.
30  {
31  public:
33  const HistAxis& axis,
34  const Cut& cut,
35  const SystShifts& shift = kNoShift,
36  const Var& wei = kUnweighted)
37  : fSpec(loaders.GetLoader(caf::kNEARDET, Loaders::kMC),
38  axis, cut, shift, wei)
39  {};
40 
42  : fSpec(spec)
43  {};
44 
45  Spectrum Predict(osc::IOscCalculator* calc) const override;
46 
48  Flavors::Flavors_t flav,
50  Sign::Sign_t sign) const override;
51 
52  void SaveTo(TDirectory* dir) const override;
53  static std::unique_ptr<TrivialPrediction> LoadFrom(TDirectory* dir);
54 
55 
56  private:
58  };
60  {
61  return fSpec;
62  }
64  Flavors::Flavors_t flav,
66  Sign::Sign_t sign) const
67  {
68  // basically the whole prediction is regarded as being one component
69  if (flav == Flavors::kNuMuToNuMu && curr == Current::kCC && sign == Sign::kBoth)
70  {
71  return fSpec;
72  }
73  else
74  {
75  auto spec = Spectrum(fSpec);
76  spec.Clear();
77  return std::move(spec);
78  }
79  }
80  void TrivialPrediction::SaveTo(TDirectory* dir) const
81  {
82  TDirectory* tmp = gDirectory;
83 
84  dir->cd();
85  TObjString("TrivialPrediction").Write("type", TObject::kOverwrite);
86  fSpec.SaveTo(dir->mkdir("spec"));
87 
88  tmp->cd();
89  }
90  std::unique_ptr<TrivialPrediction> TrivialPrediction::LoadFrom(TDirectory* dir)
91  {
92  TDirectory * specDir = dir->GetDirectory("spec");
93  assert(specDir);
94  return std::make_unique<TrivialPrediction>(*Spectrum::LoadFrom(specDir));
95  }
96 
97  //----------------------------------------------------------------------
98 
100  {
101  public:
102  NDPredGenerator(const HistAxis axis, const Cut cut, const Var wei = kUnweighted )
103  : fAxis(axis), fCut(cut), fWei(wei)
104  {};
105 
106  virtual ~NDPredGenerator() {};
107 
108  std::unique_ptr<IPrediction> Generate(
109  Loaders& loaders,
110  const SystShifts& shiftMC = kNoShift ) const override;
111 
112  private:
114  const Cut fCut;
115  const Var fWei;
116  };
117  std::unique_ptr<IPrediction> NDPredGenerator::Generate(
118  Loaders& loaders,
119  const SystShifts& shiftMC ) const
120  {
121  return std::unique_ptr<IPrediction>( new TrivialPrediction(
122  loaders, fAxis, fCut, shiftMC, fWei ) );
123  }
124 
125 
126 // Derived SingleSampleExpt to blow up chisq and let fitter ignore more penalties
128  {
129  public:
131  const Spectrum& data,
132  double scaleFactor = 1)
133  : SingleSampleExperiment(pred, data), fScaleFactor(scaleFactor)
134  {}
135 
137  const SystShifts& syst = SystShifts::Nominal()) const override
138  {
139  SingleSampleExperiment expt(fMC,fData);
140  return fScaleFactor * expt.ChiSq(osc, syst);
141  }
142 
143  private:
144  double fScaleFactor;
145  };
146 
147 //==================================================================
148 // Q0Q3 bins utilities
149 //==================================================================
150  class IRescaledSigmaSyst : public ISyst
151  {
152  public:
153  IRescaledSigmaSyst(const std::string & shortName, const std::string & latexName, double sigmaScale=1.0)
154  : ISyst(shortName, latexName), fSigmaScale(sigmaScale) {};
155  double GetSigmaScale() const { return fSigmaScale; };
156  void SetSigmaScale(double sc) { fSigmaScale = sc; };
157  protected:
158  double fSigmaScale;
159  };
160 
162  {
163  public:
164  CompNormSyst(const Cut & selCut, double sigmaScale=1.0)
165  : fID(sInstanceCount++),
166  IRescaledSigmaSyst("CompNormShift_" + std::to_string(sInstanceCount), "Component Normalization Shift " + std::to_string(sInstanceCount), sigmaScale),
167  fSelCut(selCut)
168  {}
169  void Shift(double sigma, caf::SRProxy* sr, double& weight) const override;
170 
171  private:
172  const Cut fSelCut;
173  unsigned int fID;
174 
175  static unsigned int sInstanceCount;
176  };
177  unsigned int CompNormSyst::sInstanceCount = 0;
178  void CompNormSyst::Shift(double sigma, caf::SRProxy* sr, double& weight) const
179  {
180  if (!fSelCut(sr))
181  return;
182 
183  double wgt = 1 + fSigmaScale * sigma;
184  if (wgt < 0)
185  wgt = 0;
186  weight *= wgt;
187  }
188 
189  const Cut GetCutIsFitMEC( const bool isRHC )
190  {
191  const Cut kIsFitMEC( [ isRHC ]( const caf::SRProxy *sr )
192  {
193  if ( !kIsDytmanMEC( sr ) || !sr->mc.nu[0].iscc ) return false;
194  if ( !isRHC && sr->mc.nu[0].pdg == 14 ) return true;
195  else if ( isRHC && sr->mc.nu[0].pdg == -14 ) return true;
196  else return false;
197  });
198  return kIsFitMEC;
199  }
200 
201  struct Q3Q0Bin
202  {
203  float loQ3;
204  float hiQ3;
205  float loQ0;
206  float hiQ0;
207  };
208 
209  Cut Q3Q0CutFactory( float loQ3, float hiQ3, float loQ0, float hiQ0 )
210  {
211  return std::move( Cut ( [ loQ3, hiQ3, loQ0, hiQ0 ]( const caf::SRProxy * sr )
212  {
213  double q3 = kTrueQ3( sr );
214  double q0 = kTrueQ0( sr );
215  return q3 > loQ3 && q3 < hiQ3 && q0 > loQ0 && q0 < hiQ0;
216  }
217  ));
218  }
219 
220 //===============================================================================================================
221 // Valencia MEC 2D Gaussian Enhancement
222 //===============================================================================================================
224 
225 double CalcMECGaussEnh( const double q0, const double q3, const MECGaussEnhParam shift_param, const double shift_sigma )
226 {
227  double norm = 10.0;
228  double mean_q0 = 0.25;
229  double mean_q3 = 0.50;
230  double sigma_q0 = 0.07;
231  double sigma_q3 = 0.13;
232  double corr = 0.8;
233 
234  switch( shift_param )
235  {
236  case kGauss2DNorm :
237  norm += shift_sigma * 0.5;
238  if ( norm < 0 ) norm = 0.0;
239  break;
240  case kGauss2DMeanQ0 :
241  mean_q0 += shift_sigma * 0.025;
242  break;
243  case kGauss2DMeanQ3 :
244  mean_q3 += shift_sigma * 0.05;
245  break;
246  case kGauss2DSigmaQ0 :
247  sigma_q0 += shift_sigma * 0.015;
248  if ( sigma_q0 < 0.01 ) sigma_q0 = 0.01;
249  break;
250  case kGauss2DSigmaQ3 :
251  sigma_q3 += shift_sigma * 0.02;
252  if ( sigma_q3 < 0.01 ) sigma_q3 = 0.01;
253  break;
254  case kGauss2DCorr :
255  corr += shift_sigma * 0.2;
256  if ( corr < -0.99 ) corr = -0.99;
257  else if ( corr > 0.99 ) corr = 0.99;
258  }
259 // reminder: http://mathworld.wolfram.com/BivariateNormalDistribution.html
260  double z = pow( ( q0 - mean_q0 ) / sigma_q0, 2 ) + pow( ( q3 - mean_q3 ) / sigma_q3, 2 )
261  - 2 * corr * ( q0 - mean_q0 ) * ( q3 - mean_q3 ) / ( sigma_q0 * sigma_q3 );
262 
263  return 1.0 + norm * exp( -0.5 * z / ( 1 - corr * corr ) );
264 }
265 
266 const Var kMECGaussEnh( []( const caf::SRProxy* sr )
267 {
268  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?
269  double q0 = kTrueQ0( sr );
270  double q3 = kTrueQ3( sr );
271  return CalcMECGaussEnh( q0, q3, kGauss2DNorm, 0 ); // Shift = 0 gives nominal suppression
272 });
273 
274 class MECGaussEnhSyst : public ISyst
275 {
276  public :
277  MECGaussEnhSyst( const MECGaussEnhParam& shift_param, const std::string& shift_param_name )
278  : ISyst( "MECGaussEnhSyst" + shift_param_name, "MEC 2D Gauss Syst " + shift_param_name ),
279  fShiftParam( shift_param )
280  {}
281 
282  void Shift( double sigma, caf::SRProxy* sr, double& weight ) const override
283  {
284  //if ( !kIsNumuCC || !kIsDytmanMEC ) return;
285  if ( !kIsNumuCC( sr ) || sr->mc.nu[0].mode != caf::kMEC ) return;
286  // This will work for Empirical MEC. Will it also work for Valencia MEC?
287 
288  double q0 = kTrueQ0( sr );
289  double q3 = kTrueQ3( sr );
290 
291  double wgt_nominal = CalcMECGaussEnh( q0, q3, fShiftParam, 0 );
292  double wgt_shift = CalcMECGaussEnh( q0, q3, fShiftParam, sigma );
293 
294  weight *= wgt_shift / wgt_nominal;
295  }
296 
297  private:
299 };
300 
301 
302 
303 //===============================================================================================================
304 // Valencia MEC double 2D Gaussian Enhancement
305 //===============================================================================================================
309 
310 double CalcMECDoubleGaussEnh( const double q0, const double q3, const MECDoubleGaussEnhParam shift_param, const double shift_sigma)
311 { /*// params found for mp5 files, can be taken as initial seed
312  double norm = 12.78;
313  double mean_q0 = 0.33;
314  double mean_q3 = 0.8;
315  double sigma_q0 = 0.1;
316  double sigma_q3 = 0.29;
317  double corr = 0.91;
318  double norm_2 = 38.8;
319  double mean_q0_2 = 0.036;
320  double mean_q3_2 = 0.46;
321  double sigma_q0_2 = 0.042;
322  double sigma_q3_2 = 0.233;
323  double corr_2 = 0.71;
324  double baseline = 1.0;
325  */
326 // params of prod5 files fit
327  double norm = 14.9925;
328  double mean_q0 = 0.365515;
329  double mean_q3 = 0.888951;
330  double sigma_q0 = 0.0905725;
331  double sigma_q3 = 0.340095;
332  double corr = 0.808604;
333  double norm_2 = 47.4552;
334  double mean_q0_2 = 0.0355624;
335  double mean_q3_2 = 0.444703;
336  double sigma_q0_2 = 0.0587598;
337  double sigma_q3_2 = 0.269792;
338  double corr_2 = 0.706661;
339  double baseline = 0.13100;
340 
341  switch( shift_param )
342  {
343  case kGauss2DNorm_1:
344  norm += shift_sigma * 0.5;
345  if ( norm < 0 ) norm = 0.0;
346  break;
347  case kGauss2DMeanQ0_1:
348  mean_q0 += shift_sigma * 0.025;
349  if ( mean_q0 < 0 ) mean_q0 = 0.0;
350  break;
351  case kGauss2DMeanQ3_1:
352  mean_q3 += shift_sigma * 0.05;
353  if ( mean_q3 < 0 ) mean_q0 = 0.0;
354  break;
355  case kGauss2DSigmaQ0_1:
356  sigma_q0 += shift_sigma * 0.015;
357  if ( sigma_q0 < 0.01 ) sigma_q0 = 0.01;
358  break;
359  case kGauss2DSigmaQ3_1:
360  sigma_q3 += shift_sigma * 0.02;
361  if ( sigma_q3 < 0.01 ) sigma_q3 = 0.01;
362  break;
363  case kGauss2DCorr_1 :
364  corr += shift_sigma * 0.2;
365  if ( corr < -0.99 ) corr = -0.99;
366  else if ( corr > 0.99 ) corr = 0.99;
367  break;
368  case kGauss2DNorm_2 :
369  norm_2 += shift_sigma * 2;
370  if ( norm_2 < 0 ) norm_2 = 0.0;
371  break;
372  case kGauss2DMeanQ0_2 :
373  mean_q0_2 += shift_sigma * 0.025;
374  if ( mean_q0_2 < 0 ) mean_q0_2=0;
375  break;
376  case kGauss2DMeanQ3_2 :
377  mean_q3_2 += shift_sigma * 0.05;
378  if ( mean_q3_2 < 0 ) mean_q3_2 = 0;
379  break;
380  case kGauss2DSigmaQ0_2 :
381  sigma_q0_2 += shift_sigma * 0.015;
382  if ( sigma_q0_2 < 0.01 ) sigma_q0_2 = 0.01;
383  break;
384  case kGauss2DSigmaQ3_2 :
385  sigma_q3_2 += shift_sigma * 0.02;
386  if ( sigma_q3_2 < 0.01 ) sigma_q3_2 = 0.01;
387  break;
388  case kGauss2DCorr_2 :
389  corr_2 += shift_sigma * 0.2;
390  if ( corr_2 < -0.99 ) corr_2 = -0.99;
391  if ( corr_2 > 0.99 ) corr_2 = 0.99;
392  break;
393  case kBaseline :
394  baseline += shift_sigma*0.05;
395  //if (baseline < 0.6) baseline=0.6;
396 
397  }
398 // reminder: http://mathworld.wolfram.com/BivariateNormalDistribution.html
399  double z = pow( ( q0 - mean_q0 ) / sigma_q0, 2 ) + pow( ( q3 - mean_q3 ) / sigma_q3, 2 )
400  - 2 * corr * ( q0 - mean_q0 ) * ( q3 - mean_q3 ) / ( sigma_q0 * sigma_q3 );
401 
402  double z_2 = pow( ( q0 - mean_q0_2 ) / sigma_q0_2, 2 ) + pow( ( q3 - mean_q3_2 ) / sigma_q3_2, 2 )
403  - 2 * corr_2 * ( q0 - mean_q0_2 ) * ( q3 - mean_q3_2 ) / ( sigma_q0_2 * sigma_q3_2 );
404 
405  double weight= baseline + norm * exp( -0.5 * z / ( 1 - corr * corr ) ) + norm_2* exp( -0.5 * z_2 / ( 1 - corr_2 * corr_2 ) );
406  if (weight<0) weight=0;
407  return weight;
408 
409 }
410 
411 const std::map< std::string , MECDoubleGaussEnhParam > DoubleGaussMap{
412 {"MECDoubleGaussEnhSystNorm_1", kGauss2DNorm_1},
413 {"MECDoubleGaussEnhSystMeanQ0_1", kGauss2DMeanQ0_1},
414 {"MECDoubleGaussEnhSystMeanQ3_1", kGauss2DMeanQ3_1},
415 {"MECDoubleGaussEnhSystSigmaQ0_1", kGauss2DSigmaQ0_1},
416 {"MECDoubleGaussEnhSystSigmaQ3_1", kGauss2DSigmaQ3_1 },
417 {"MECDoubleGaussEnhSystCorr_1", kGauss2DCorr_1 },
418 {"MECDoubleGaussEnhSystNorm_2", kGauss2DNorm_2 },
419 {"MECDoubleGaussEnhSystMeanQ0_2", kGauss2DMeanQ0_2},
420 {"MECDoubleGaussEnhSystMeanQ3_2", kGauss2DMeanQ3_2},
421 {"MECDoubleGaussEnhSystSigmaQ0_2", kGauss2DSigmaQ0_2},
422 {"MECDoubleGaussEnhSystSigmaQ3_2", kGauss2DSigmaQ3_2},
423 {"MECDoubleGaussEnhSystCorr_2", kGauss2DCorr_2 },
424 {"MECDoubleGaussEnhSystBaseline", kBaseline}
425 };
426 
427 double CalcMECDoubleGaussEnhShiftedParam( std::string shift_param, double shift_sigma)
428 {
429 
430  MECDoubleGaussEnhParam shifted_param = DoubleGaussMap.find(shift_param)->second;
431  /*// params found for mp5 files, can be taken as initial seed
432  double norm = 12.78;
433  double mean_q0 = 0.33;
434  double mean_q3 = 0.8;
435  double sigma_q0 = 0.1;
436  double sigma_q3 = 0.29;
437  double corr = 0.91;
438  double norm_2 = 38.8;
439  double mean_q0_2 = 0.036;
440  double mean_q3_2 = 0.46;
441  double sigma_q0_2 = 0.042;
442  double sigma_q3_2 = 0.233;
443  double corr_2 = 0.71;
444  double baseline = 1.0;
445  */
446 // params of prod5 files fit
447  double norm = 14.9925;
448  double mean_q0 = 0.365515;
449  double mean_q3 = 0.888951;
450  double sigma_q0 = 0.0905725;
451  double sigma_q3 = 0.340095;
452  double corr = 0.808604;
453  double norm_2 = 47.4552;
454  double mean_q0_2 = 0.0355624;
455  double mean_q3_2 = 0.444703;
456  double sigma_q0_2 = 0.0587598;
457  double sigma_q3_2 = 0.269792;
458  double corr_2 = 0.706661;
459  double baseline = 0.13100;
460 
461  double param = 0.0;
462 
463  switch( shifted_param )
464  {
465  case kGauss2DNorm_1:
466  norm += shift_sigma * 0.5; //.5
467  if ( norm < 0 ) norm = 0.0;
468  param = norm;
469  break;
470  case kGauss2DMeanQ0_1:
471  mean_q0 += shift_sigma * 0.025;
472  if ( mean_q0 < 0 ) mean_q0 = 0.0;
473  param = mean_q0;
474  break;
475  case kGauss2DMeanQ3_1:
476  mean_q3 += shift_sigma * 0.05;
477  if ( mean_q3 < 0 ) mean_q3 = 0.0;
478  param = mean_q3;
479  break;
480  case kGauss2DSigmaQ0_1:
481  sigma_q0 += shift_sigma * 0.015;
482  if ( sigma_q0 < 0.001 ) sigma_q0 = 0.001;
483  param =sigma_q0;
484  break;
485  case kGauss2DSigmaQ3_1:
486  sigma_q3 += shift_sigma * 0.02;
487  if ( sigma_q3 < 0.001 ) sigma_q3 = 0.001;
488  param = sigma_q3;
489  break;
490  case kGauss2DCorr_1 :
491  corr += shift_sigma * 0.2;
492  if ( corr < -0.99 ) corr = -0.99;
493  else if ( corr > 0.99 ) corr = 0.99;
494  param=corr;
495  break;
496  case kGauss2DNorm_2 :
497  norm_2 += shift_sigma * 2;
498  if ( norm_2 < 0 ) norm_2 = 0.0;
499  param = norm_2;
500  break;
501  case kGauss2DMeanQ0_2 :
502  mean_q0_2 += shift_sigma * 0.025;
503  if ( mean_q0_2 < 0 ) mean_q0_2 = 0.0;
504  param = mean_q0_2;
505  break;
506  case kGauss2DMeanQ3_2 :
507  mean_q3_2 += shift_sigma * 0.05;
508  if ( mean_q3_2 < 0 ) mean_q3_2 = 0.0;
509  param = mean_q3_2;
510  break;
511  case kGauss2DSigmaQ0_2 :
512  sigma_q0_2 += shift_sigma * 0.015;
513  if ( sigma_q0_2 < 0.001 ) sigma_q0_2 = 0.001;
514  param = sigma_q0_2;
515  break;
516  case kGauss2DSigmaQ3_2 :
517  sigma_q3_2 += shift_sigma * 0.02;
518  if ( sigma_q3_2 < 0.001 ) sigma_q3_2 = 0.001;
519  param = sigma_q3_2;
520  break;
521  case kGauss2DCorr_2 :
522  corr_2 += shift_sigma * 0.2;
523  if ( corr_2 < -0.99 ) corr_2 = -0.99;
524  if ( corr_2 > 0.99 ) corr_2 = 0.99;
525  param = corr_2;
526  break;
527  case kBaseline :
528  baseline += shift_sigma*0.05;
529  // if (baseline < 0.6) baseline=0.6;
530  param = baseline;
531  }
532  return param;
533 }
534 
535 
536 const Var kMECDoubleGaussEnh( []( const caf::SRProxy* sr )
537 {
538  if ( !kIsNumuCC( sr ) || sr->mc.nu[0].mode != caf::kMEC ) return 1.0;
539  double q0 = kTrueQ0( sr );
540  double q3 = kTrueQ3( sr );
541  return CalcMECDoubleGaussEnh( q0, q3, kGauss2DNorm_1, 0 ); // Shift = 0 gives nominal suppression
542 });
543 
544 double simpleMECDoubleGaussEnh( const double q0, const double q3,
545  double norm, double mean_q0 , double mean_q3, double sigma_q0, double sigma_q3, double corr,
546  double norm_2, double mean_q0_2 , double mean_q3_2, double sigma_q0_2, double sigma_q3_2, double corr_2,
547  double baseline)
548 {
549  double z = pow( ( q0 - mean_q0 ) / sigma_q0, 2 ) + pow( ( q3 - mean_q3 ) / sigma_q3, 2 )
550  - 2 * corr * ( q0 - mean_q0 ) * ( q3 - mean_q3 ) / ( sigma_q0 * sigma_q3 );
551  double z_2 = pow( ( q0 - mean_q0_2 ) / sigma_q0_2, 2 ) + pow( ( q3 - mean_q3_2 ) / sigma_q3_2, 2 )
552  - 2 * corr_2 * ( q0 - mean_q0_2 ) * ( q3 - mean_q3_2 ) / ( sigma_q0_2 * sigma_q3_2 );
553 
554  return baseline + norm * exp( -0.5 * z / ( 1 - corr * corr ) ) + norm_2* exp( -0.5 * z_2 / ( 1 - corr_2 * corr_2 ) );
555 }
556 
557 const Var kMECDoubleGaussEnhSimple(double norm, double mean_q0 , double mean_q3, double sigma_q0, double sigma_q3, double corr,
558  double norm_2, double mean_q0_2 , double mean_q3_2, double sigma_q0_2, double sigma_q3_2, double corr_2,
559  double baseline )
560 {
561  const Var kMECDoubleGaussEnhSimpleVar( [norm, mean_q0 , mean_q3, sigma_q0, sigma_q3, corr,
562  norm_2, mean_q0_2 , mean_q3_2, sigma_q0_2, sigma_q3_2, corr_2, baseline]( const caf::SRProxy* sr )
563  {
564  // no shifts taken in this weights, we calculate a simple
565  if ( !kIsNumuCC( sr ) || sr->mc.nu[0].mode != caf::kMEC ) return 1.0;
566  double q0 = kTrueQ0( sr );
567  double q3 = kTrueQ3( sr );
568  return simpleMECDoubleGaussEnh( q0, q3,
569  norm, mean_q0 , mean_q3, sigma_q0, sigma_q3, corr,
570  norm_2, mean_q0_2 , mean_q3_2, sigma_q0_2, sigma_q3_2, corr_2, baseline );
571  });
572  return kMECDoubleGaussEnhSimpleVar;
573 };
574 
575 
577 {
578  public :
579  MECDoubleGaussEnhSyst( const MECDoubleGaussEnhParam& shift_param, const std::string& shift_param_name )
580  : ISyst( "MECDoubleGaussEnhSyst" + shift_param_name, "MEC 2D Gauss Syst " + shift_param_name ),
581  fShiftParam( shift_param )
582  {}
583 
584  void Shift( double sigma, caf::SRProxy* sr, double& weight ) const override
585  {
586  //if ( !kIsNumuCC || !kIsDytmanMEC ) return;
587  if ( !kIsNumuCC( sr ) || sr->mc.nu[0].mode != caf::kMEC ) return;
588  // This will work for Empirical MEC. Will it also work for Valencia MEC?
589 
590  double q0 = kTrueQ0( sr );
591  double q3 = kTrueQ3( sr );
592 
593  double wgt_nominal = CalcMECDoubleGaussEnh( q0, q3, fShiftParam, 0 );
594  double wgt_shift = CalcMECDoubleGaussEnh( q0, q3, fShiftParam, sigma );
595 
596  weight *= wgt_shift / wgt_nominal;
597  }
598 
599  private:
601 };
602 
603 
604 
605 //===============================================================================================================
606 // MINOS Resonance Suppression
607 //===============================================================================================================
608 
610 
611 double CalcMinosResSupp( const double Q2, const MinosResSuppParam shift_param, const double shift_sigma )
612 {
613  double A = 1.010;
614  double Q0 = 0.156;
615 
616  switch( shift_param )
617  {
618  case kMinosResSuppNorm :
619  A += shift_sigma * 0.1;
620  if ( A < 0 ) A = 0.0;
621  break;
622  case kMinosResSuppQ0 :
623  Q0 += shift_sigma * 0.03;
624  if ( Q0 < 0.01 ) Q0 = 0.01;
625  }
626 
627  double supp = A / ( 1 + exp( 1 - sqrt( Q2 ) / Q0 ) );
628  if ( supp > 1 ) supp = 1.0;
629 
630  return supp;
631 }
632 
633 const Var kMinosResSupp( []( const caf::SRProxy* sr )
634 {
635  if ( !kIsNumuCC( sr ) || sr->mc.nu[0].mode != kIsRes( sr ) || sr->mc.nu[0].tgtA == 1 ) return 1.0; // Exclude hydrogen
636  double Q2 = kTrueQ2( sr );
637  return CalcMinosResSupp( Q2, kMinosResSuppNorm, 0 ); // Shift = 0 gives nominal suppression
638 });
639 
640 class MinosResSuppSyst : public ISyst
641 {
642  public :
643  MinosResSuppSyst( const MinosResSuppParam& shift_param, const std::string& shift_param_name )
644  : ISyst( "MinosResSuppSyst" + shift_param_name, "MINOS RES Supp Syst " + shift_param_name ),
645  fShiftParam( shift_param )
646  {}
647 
648  void Shift( double sigma, caf::SRProxy* sr, double& weight ) const override
649  {
650  if ( !kIsNumuCC( sr ) || sr->mc.nu[0].mode != kIsRes( sr ) || sr->mc.nu[0].tgtA == 1 ) return; // Exclude hydrogen
651 
652  double Q2 = kTrueQ2( sr );
653 
654  double wgt_nominal = CalcMinosResSupp( Q2, fShiftParam, 0 );
655  double wgt_shift = CalcMinosResSupp( Q2, fShiftParam, sigma );
656 
657  weight *= wgt_shift / wgt_nominal;
658  }
659 
660  private:
662 };
663 
664 //=============================================================================
665 // Extra Weight
666 //=============================================================================
667 
668 const Var kWeightPionDeuteriumTune( []( const caf::SRProxy* sr )
669 {
670  double weight = 1.0;
671 
672  if ( kIsNumuCC( sr ) && kIsRes( sr ) )
673  {
674  if ( sr->mc.nu[0].rwgt.genie.size() <= rwgt::fReweightMaCCRES - 1 )
675  {
676  if ( sr->mc.nu[0].isvtxcont )
677  throw std::runtime_error( "kWeightPionDeuteriumTune: Cannot do MA CCRES rescaling without GENIE reweights available." );
678  else
679  return 1.0;
680  }
681  // GENIE single pion production deuterium fit Eur. Phys. J. C (2016) 76: 474
682  // Scale down MA CCRES from GENIE default 1.12 GeV (+/-1 sigma = +/-20%) to Deuterium fit 0.94+/-0.05 GeV
683  const double correctionInSigma = ( 1.12 - 0.94 ) / 0.2;
684  weight *= 1. + correctionInSigma * ( sr->mc.nu[0].rwgt.genie[rwgt::fReweightMaCCRES].minus1sigma - 1. );
685 
686  // Apply 1.15 normalization correction from Deuterium fit
687  weight *= 1.15;
688  }
689  else if ( kIsNumuCC( sr ) && kIsDIS( sr ) && kTrueW( sr ) < 1.7 )
690  {
691  // Apply non-RES 1pi normalization from deuterium fit
692  if ( sr->mc.nu[0].npiminus + sr->mc.nu[0].npiplus + sr->mc.nu[0].npizero == 1 )
693  {
694  weight *= 0.43;
695  }
696  }
697  return weight;
698 });
699 
700 }
701 
const MinosResSuppParam fShiftParam
MECDoubleGaussEnhParam
const Cut kIsNumuCC(CCFlavSel(14, 14))
Select CC .
MinosResSuppSyst(const MinosResSuppParam &shift_param, const std::string &shift_param_name)
double CalcMECDoubleGaussEnh(const double q0, const double q3, const MECDoubleGaussEnhParam shift_param, const double shift_sigma)
TrivialPrediction(const Spectrum &spec)
osc::OscCalculatorDumb calc
Oscillation analysis framework, runs over CAF files outside of ART.
General interface to oscillation calculators.
Definition: FwdDeclare.h:15
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:991
VectorProxy< SRNeutrinoProxy > nu
Definition: SRProxy.h:1890
const Var weight
Spectrum Predict(osc::IOscCalculator *calc) const override
const HistAxis fAxis
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
const Cut kIsRes
Definition: TruthCuts.cxx:111
Collection of SpectrumLoaders for many configurations.
Definition: Loaders.h:17
constexpr T pow(T x)
Definition: pow.h:75
const Var kMECDoubleGaussEnh([](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 CalcMECDoubleGaussEnh(q0, q3, kGauss2DNorm_1, 0);})
virtual double ChiSq(osc::IOscCalculatorAdjustable *osc, const SystShifts &syst=SystShifts::Nominal()) const override
const Cut kIsDIS
Definition: TruthCuts.cxx:118
Proxy for StandardRecord.
Definition: SRProxy.h:2237
double corr(TGraph *g, double thresh)
Float_t tmp
Definition: plot.C:36
MECDoubleGaussEnhSyst(const MECDoubleGaussEnhParam &shift_param, const std::string &shift_param_name)
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);})
static SystShifts Nominal()
Definition: SystShifts.h:34
const Color_t kMC
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
Spectrum PredictComponent(osc::IOscCalculator *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const override
const Var kMinosResSupp([](const caf::SRProxy *sr){if(!kIsNumuCC(sr)||sr->mc.nu[0].mode!=kIsRes(sr)||sr->mc.nu[0].tgtA==1) return 1.0;double Q2=kTrueQ2(sr);return CalcMinosResSupp(Q2, kMinosResSuppNorm, 0);})
const Var kTrueQ2
Definition: TruthVars.h:27
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:37
MECGaussEnhParam
void Shift(double sigma, caf::SRProxy *sr, double &weight) const override
Perform the systematic shift.
const XML_Char const XML_Char * data
Definition: expat.h:268
const Cut kIsDytmanMEC
Definition: TruthCuts.cxx:187
double CalcMECGaussEnh(const double q0, const double q3, const MECGaussEnhParam shift_param, const double shift_sigma)
Charged-current interactions.
Definition: IPrediction.h:39
Definition: Shift.h:6
expt
Definition: demo5.py:34
General interface to any calculator that lets you set the parameters.
const Var kTrueQ3
Definition: TruthVars.h:38
double CalcMinosResSupp(const double Q2, const MinosResSuppParam shift_param, const double shift_sigma)
if(dump)
CompNormSyst(const Cut &selCut, double sigmaScale=1.0)
void Shift(double sigma, caf::SRProxy *sr, double &weight) const override
Perform the systematic shift.
Cut Q3Q0CutFactory(float loQ3, float hiQ3, float loQ0, float hiQ0)
double CalcMECDoubleGaussEnhShiftedParam(std::string shift_param, double shift_sigma)
Near Detector in the NuMI cavern.
unsigned int fID
double simpleMECDoubleGaussEnh(const double q0, const double q3, std::map< std::string, double > DGPars)
Definition: XsecTunes.cxx:61
void SaveTo(TDirectory *dir) const override
const Var kMECDoubleGaussEnhSimple(double norm, double mean_q0, double mean_q3, double sigma_q0, double sigma_q3, double corr, double norm_2, double mean_q0_2, double mean_q3_2, double sigma_q0_2, double sigma_q3_2, double corr_2, double baseline)
std::unique_ptr< IPrediction > Generate(Loaders &loaders, const SystShifts &shiftMC=kNoShift) const override
const ana::Var wgt
static unsigned int sInstanceCount
z
Definition: test.py:28
void SetSigmaScale(double sc)
Oscillation probability calculators.
Definition: Calcs.h:5
MECGaussEnhSyst(const MECGaussEnhParam &shift_param, const std::string &shift_param_name)
def sign(x)
Definition: canMan.py:204
const SystShifts kNoShift
Definition: SystShifts.h:112
(&#39; survival&#39;)
Definition: IPrediction.h:19
NDPredGenerator(const HistAxis axis, const Cut cut, const Var wei=kUnweighted)
static std::unique_ptr< TrivialPrediction > LoadFrom(TDirectory *dir)
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
const Cut cut
Definition: exporter_fd.C:30
static const double A
Definition: Units.h:82
IRescaledSigmaSyst(const std::string &shortName, const std::string &latexName, double sigmaScale=1.0)
const Var kWeightPionDeuteriumTune([](const caf::SRProxy *sr){double weight=1.0;if(kIsNumuCC(sr)&&kIsRes(sr)){if(sr->mc.nu[0].rwgt.genie.size()<=rwgt::fReweightMaCCRES-1){if(sr->mc.nu[0].isvtxcont) throw std::runtime_error("kWeightPionDeuteriumTune: Cannot do MA CCRES rescaling without GENIE reweights available.");else return 1.0;} const double correctionInSigma=(1.12-0.94)/0.2;weight *=1.+correctionInSigma *(sr->mc.nu[0].rwgt.genie[rwgt::fReweightMaCCRES].minus1sigma-1.);weight *=1.15;}else if(kIsNumuCC(sr)&&kIsDIS(sr)&&kTrueW(sr)< 1.7){if(sr->mc.nu[0].npiminus+sr->mc.nu[0].npiplus+sr->mc.nu[0].npizero==1){weight *=0.43;}}return weight;})
Float_t norm
const Cut GetCutIsFitMEC(const bool isRHC)
double ChiSq(osc::IOscCalculatorAdjustable *osc, const SystShifts &syst=SystShifts::Nominal()) const override
TDirectory * dir
Definition: macro.C:5
BigChi2SingleSampleExperiment(const IPrediction *pred, const Spectrum &data, double scaleFactor=1)
const MECGaussEnhParam fShiftParam
std::vector< Loaders * > loaders
Definition: syst_header.h:385
string syst
Definition: plotSysts.py:176
SRTruthBranchProxy mc
Definition: SRProxy.h:2253
const MECDoubleGaussEnhParam fShiftParam
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.
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
void Shift(double sigma, caf::SRProxy *sr, double &weight) const override
Perform the systematic shift.
Given loaders and an MC shift, Generate() generates an IPrediction.
This module creates Common Analysis Files.
Definition: CSVMaker.h:8
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir)
Definition: Spectrum.cxx:1055
double GetSigmaScale() const
const Var kTrueW
Definition: TruthVars.h:22
const std::map< std::string, MECDoubleGaussEnhParam > DoubleGaussMap
const Var kUnweighted
The simplest possible Var, always 1. Used as a default weight.
Definition: Var.h:100
MinosResSuppParam
TrivialPrediction(Loaders &loaders, const HistAxis &axis, const Cut &cut, const SystShifts &shift=kNoShift, const Var &wei=kUnweighted)
static constexpr Double_t sr
Definition: Munits.h:164
Compare a single data spectrum to the MC + cosmics expectation.
void SaveTo(TDirectory *dir) const
Definition: Spectrum.cxx:1029
Template for Var and SpillVar.
Definition: Var.h:16
void Shift(caf::SRProxy *sr, double &weight) const
Definition: SystShifts.cxx:154
const Var GetWeightFromShifts(const SystShifts &sys_shifts)