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