MECTuningUtils.h
Go to the documentation of this file.
1 #pragma once
2 
3 #include "CAFAna/Core/Loaders.h"
4 #include <TObjString.h>
5 #include "3FlavorAna/Cuts/NumuCuts2020.h"
6 #include "CAFAna/Vars/TruthVars.h"
7 #include "CAFAna/Experiment/SingleSampleExperiment.h"
8 #include "StandardRecord/Proxy/SRProxy.h"
9 
10 #include "DoubleGaussParams.h"
11 
12 
13 namespace ana
14 {
15 //-------------------------------------------------------------------------------------//
16 // cut used for official tune, not the same as official kNumu2020ND
17 // the fit was made earlier.
18 
19  const Cut kNumu2020PIDLoosePTPxsec(
20  [](const caf::SRProxy* sr)
21  {
22  return (sr->sel.remid.pid > 0.7 && sr->sel.cvnloosepreselptp.numuid > 0.82);
23  }
24  );
25  const Cut kNumu2020NDxsec = kNumuQuality && kNumuContainND2020 && kNumu2020PIDLoosePTPxsec;
26 //-----------------------------------------------------------------------------------------//
27 
28  const Var GetWeightFromShifts( const SystShifts& sys_shifts )
29  {
30  const Var kWeightFromShifts( [ sys_shifts ]( const caf::SRProxy *sr )
31  {
32  caf::SRProxy* face_palm = const_cast<caf::SRProxy*>( sr );
33  double weight = 1.0;
34  sys_shifts.Shift( face_palm, weight );
35  return weight;
36  });
37  return kWeightFromShifts;
38  }
39 
40  // this thing basically needs to be an IPrediction that gives
41  // the correct components of the spectrum.
42  // or at least have a PredictComponent() that can do
43  // (Flavors::kAll, Current::kBoth, Sign::kBoth)
44  // when asked, since that's the one that you need
45  // to get PredictSyst() to work.
46  class TrivialPrediction : public IPrediction
47  {
48  public:
49  TrivialPrediction(Loaders& loaders,
50  const HistAxis& axis,
51  const Cut& cut,
52  const SystShifts& shift = kNoShift,
53  const Var& wei = kUnweighted)
54  : fSpec(loaders.GetLoader(caf::kNEARDET, Loaders::kMC),
55  axis, cut, shift, wei)
56  {};
57 
58  TrivialPrediction(const Spectrum& spec)
59  : fSpec(spec)
60  {};
61 
62  Spectrum Predict(osc::IOscCalc* calc) const override;
63 
64  Spectrum PredictComponent(osc::IOscCalc* calc,
65  Flavors::Flavors_t flav,
66  Current::Current_t curr,
67  Sign::Sign_t sign) const override;
68 
69  void SaveTo(TDirectory* dir, const std::string& name) const override;
70  static std::unique_ptr<TrivialPrediction> LoadFrom(TDirectory* dir, const std::string& name);
71 
72 
73  private:
74  Spectrum fSpec;
75  };
76  Spectrum TrivialPrediction::Predict(osc::IOscCalc* calc) const
77  {
78  return fSpec;
79  }
80  Spectrum TrivialPrediction::PredictComponent(osc::IOscCalc* calc,
81  Flavors::Flavors_t flav,
82  Current::Current_t curr,
83  Sign::Sign_t sign) const
84  {
85  // basically the whole prediction is regarded as being one component
86  if (flav == Flavors::kNuMuToNuMu && curr == Current::kCC && sign == Sign::kBoth)
87  {
88  return fSpec;
89  }
90  else
91  {
92  auto spec = Spectrum(fSpec);
93  spec.Clear();
94  return std::move(spec);
95  }
96  }
97  void TrivialPrediction::SaveTo(TDirectory* dir, const std::string& name) const
98  {
99  TDirectory* tmp = gDirectory;
100 
101  dir = dir->mkdir(name.c_str()); // switch to subdir
102  dir->cd();
103 
104  TObjString("TrivialPrediction").Write("type", TObject::kOverwrite);
105  fSpec.SaveTo(dir, "spec");
106 
107  dir->Write();
108  delete dir;
109 
110  tmp->cd();
111  }
112  std::unique_ptr<TrivialPrediction> TrivialPrediction::LoadFrom(TDirectory* dir, const std::string& name)
113  {
114  dir = dir->GetDirectory(name.c_str()); // switch to subdir
115  assert(dir);
116 
117  return std::make_unique<TrivialPrediction>(*Spectrum::LoadFrom(dir, "spec"));
118  }
119 
120  //----------------------------------------------------------------------
121 
122  class NDPredGenerator: public IPredictionGenerator
123  {
124  public:
125  NDPredGenerator(const HistAxis axis, const Cut cut, const Var wei = kUnweighted )
126  : fAxis(axis), fCut(cut), fWei(wei)
127  {};
128 
129  virtual ~NDPredGenerator() {};
130 
131  std::unique_ptr<IPrediction> Generate(
132  Loaders& loaders,
133  const SystShifts& shiftMC = kNoShift ) const override;
134 
135  private:
136  const HistAxis fAxis;
137  const Cut fCut;
138  const Var fWei;
139  };
140  std::unique_ptr<IPrediction> NDPredGenerator::Generate(
141  Loaders& loaders,
142  const SystShifts& shiftMC ) const
143  {
144  return std::unique_ptr<IPrediction>( new TrivialPrediction(
145  loaders, fAxis, fCut, shiftMC, fWei ) );
146  }
147 
148 
149 // Derived SingleSampleExpt to blow up chisq and let fitter ignore more penalties
150 class BigChi2SingleSampleExperiment : public SingleSampleExperiment
151  {
152  public:
153  BigChi2SingleSampleExperiment(const IPrediction* pred,
154  const Spectrum& data,
155  double scaleFactor = 1)
156  : SingleSampleExperiment(pred, data), fScaleFactor(scaleFactor)
157  {}
158 
159  double ChiSq(osc::IOscCalcAdjustable* osc,
160  const SystShifts& syst = SystShifts::Nominal()) const override
161  {
162  SingleSampleExperiment expt(fMC,fData);
163  return fScaleFactor * expt.ChiSq(osc, syst);
164  }
165 
166  private:
167  double fScaleFactor;
168  };
169 
170 //================================================================================
171 // Q0Q3 bins utilities
172 //================================================================================
173  class IRescaledSigmaSyst : public ISyst
174  {
175  public:
176  IRescaledSigmaSyst(const std::string & shortName, const std::string & latexName, double sigmaScale=1.0)
177  : ISyst(shortName, latexName), fSigmaScale(sigmaScale) {};
178  double GetSigmaScale() const { return fSigmaScale; };
179  void SetSigmaScale(double sc) { fSigmaScale = sc; };
180  protected:
181  double fSigmaScale;
182  };
183 
184  class CompNormSyst : public IRescaledSigmaSyst
185  {
186  public:
187  CompNormSyst(const Cut & selCut, double sigmaScale=1.0)
188  : fID(sInstanceCount++),
189  IRescaledSigmaSyst("CompNormShift_" + std::to_string(sInstanceCount), "Component Normalization Shift " + std::to_string(sInstanceCount), sigmaScale),
190  fSelCut(selCut)
191  {}
192  void Shift(double sigma, caf::SRProxy* sr, double& weight) const override;
193 
194  private:
195  const Cut fSelCut;
196  unsigned int fID;
197 
198  static unsigned int sInstanceCount;
199  };
200  unsigned int CompNormSyst::sInstanceCount = 0;
201  void CompNormSyst::Shift(double sigma, caf::SRProxy* sr, double& weight) const
202  {
203  if (!fSelCut(sr))
204  return;
205 
206  double wgt = 1 + fSigmaScale * sigma;
207  if (wgt < 0)
208  wgt = 0;
209  weight *= wgt;
210  }
211 
212 
213 // Cut for binned fit
214  const Cut GetCutIsFitMEC( const bool isRHC )
215  {
216  const Cut kIsFitMEC( [ isRHC ]( const caf::SRProxy *sr )
217  {
218  if ( !kIsDytmanMEC( sr ) || !sr->mc.nu[0].iscc ) return false;
219  if ( !isRHC && sr->mc.nu[0].pdg == 14 ) return true;
220  else if ( isRHC && sr->mc.nu[0].pdg == -14 ) return true;
221  else return false;
222  });
223  return kIsFitMEC;
224  }
225 
226  struct Q3Q0Bin
227  {
228  float loQ3;
229  float hiQ3;
230  float loQ0;
231  float hiQ0;
232  };
233 
234  Cut Q3Q0CutFactory( float loQ3, float hiQ3, float loQ0, float hiQ0 )
235  {
236  return std::move( Cut ( [ loQ3, hiQ3, loQ0, hiQ0 ]( const caf::SRProxy * sr )
237  {
238  double q3 = kTrueQ3( sr );
239  double q0 = kTrueQ0( sr );
240  return q3 > loQ3 && q3 < hiQ3 && q0 > loQ0 && q0 < hiQ0;
241  }
242  ));
243  }
244 
245 //===============================================================================================================
246 // Valencia MEC 2D Gaussian Enhancement
247 //===============================================================================================================
248 enum MECGaussEnhParam{ kGauss2DNorm, kGauss2DMeanQ0, kGauss2DMeanQ3, kGauss2DSigmaQ0, kGauss2DSigmaQ3, kGauss2DCorr };
249 
250 double CalcMECGaussEnh( const double q0, const double q3, const MECGaussEnhParam shift_param, const double shift_sigma )
251 {
252  double norm = 10.0;
253  double mean_q0 = 0.25;
254  double mean_q3 = 0.50;
255  double sigma_q0 = 0.07;
256  double sigma_q3 = 0.13;
257  double corr = 0.8;
258 
259  switch( shift_param )
260  {
261  case kGauss2DNorm :
262  norm += shift_sigma * 0.5;
263  if ( norm < 0 ) norm = 0.0;
264  break;
265  case kGauss2DMeanQ0 :
266  mean_q0 += shift_sigma * 0.025;
267  break;
268  case kGauss2DMeanQ3 :
269  mean_q3 += shift_sigma * 0.05;
270  break;
271  case kGauss2DSigmaQ0 :
272  sigma_q0 += shift_sigma * 0.015;
273  if ( sigma_q0 < 0.01 ) sigma_q0 = 0.01;
274  break;
275  case kGauss2DSigmaQ3 :
276  sigma_q3 += shift_sigma * 0.02;
277  if ( sigma_q3 < 0.01 ) sigma_q3 = 0.01;
278  break;
279  case kGauss2DCorr :
280  corr += shift_sigma * 0.2;
281  if ( corr < -0.99 ) corr = -0.99;
282  else if ( corr > 0.99 ) corr = 0.99;
283  }
284 // reminder: http://mathworld.wolfram.com/BivariateNormalDistribution.html
285  double z = pow( ( q0 - mean_q0 ) / sigma_q0, 2 ) + pow( ( q3 - mean_q3 ) / sigma_q3, 2 )
286  - 2 * corr * ( q0 - mean_q0 ) * ( q3 - mean_q3 ) / ( sigma_q0 * sigma_q3 );
287 
288  return 1.0 + norm * exp( -0.5 * z / ( 1 - corr * corr ) );
289 }
290 
291 const Var kMECGaussEnh( []( const caf::SRProxy* sr )
292 {
293  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?
294  double q0 = kTrueQ0( sr );
295  double q3 = kTrueQ3( sr );
296  return CalcMECGaussEnh( q0, q3, kGauss2DNorm, 0 ); // Shift = 0 gives nominal suppression
297 });
298 
299 class MECGaussEnhSyst : public ISyst
300 {
301  public :
302  MECGaussEnhSyst( const MECGaussEnhParam& shift_param, const std::string& shift_param_name )
303  : ISyst( "MECGaussEnhSyst" + shift_param_name, "MEC 2D Gauss Syst " + shift_param_name ),
304  fShiftParam( shift_param )
305  {}
306 
307  void Shift( double sigma, caf::SRProxy* sr, double& weight ) const override
308  {
309  //if ( !kIsNumuCC || !kIsDytmanMEC ) return;
310  if ( !kIsNumuCC( sr ) || sr->mc.nu[0].mode != caf::kMEC ) return;
311  // This will work for Empirical MEC. Will it also work for Valencia MEC?
312 
313  double q0 = kTrueQ0( sr );
314  double q3 = kTrueQ3( sr );
315 
316  double wgt_nominal = CalcMECGaussEnh( q0, q3, fShiftParam, 0 );
317  double wgt_shift = CalcMECGaussEnh( q0, q3, fShiftParam, sigma );
318 
319  weight *= wgt_shift / wgt_nominal;
320  }
321 
322  private:
323  const MECGaussEnhParam fShiftParam;
324 };
325 
326 //===============================================================================================================
327 // Valencia MEC double 2D Gaussian Enhancement
328 //===============================================================================================================
329 enum MECDoubleGaussEnhParam{ kGauss2DNorm_1, kGauss2DMeanQ0_1, kGauss2DMeanQ3_1, kGauss2DSigmaQ0_1, kGauss2DSigmaQ3_1, kGauss2DCorr_1,
330  kGauss2DNorm_2, kGauss2DMeanQ0_2, kGauss2DMeanQ3_2, kGauss2DSigmaQ0_2, kGauss2DSigmaQ3_2, kGauss2DCorr_2,
331  kBaseline, kExpA, kExpAlpha, kLogisticC, kLogisticBeta};
332 
333 double CalcMECDoubleGaussEnh( const double q0, const double q3, const double Q2,
334  const MECDoubleGaussEnhParam shift_param, const double shift_sigma)
335 { // official parameters of prod5 commented out
336  // using MCNP first iteration now...
337 
338 // double norm= 10.0;;//14.85;
339 // double mean_q0=1.05;//0.39 ;// 0.36;//0.161587;//0.272638;//0.383689;//0.49474;//0.605791;//0.716842;//0.827893;//0.938944;//
340 // double mean_q3=0.4;//0.95 ;//0.86;
341 // double sigma_q0= 0.1;//0.12 ;// 0.13;
342 // double sigma_q3= 0.04;//0.37;// 0.35;
343 // double corr= 0.9;// 0.89;
344 // double norm_2=40.0;//42.0;
345 // double mean_q0_2= 0.19;// 0.034;
346 // double mean_q3_2= 0.55;//0.45;
347 // double sigma_q0_2= 0.1;// 0.044;
348 // double sigma_q3_2= 0.1;//0.31;
349 // double corr_2= 0.9;// 0.75;
350 // double baseline= 0.;//-0.12;// -0.08;
351 
352 
353  /* Double_t fparams[13] = {0.05, //mean q0
354  0.1, //sigma q0
355  0.4, //mean q3
356  0.04, //sigma q3
357  0.9, //corr
358  0.19, // mean q0 2
359  0.1, // sigma q0 2
360  0.55, // mean q3 2
361  0.01, //sigma q3 2
362  0.9, // corr 2
363  10.0, //baseline
364  14.0, //norm 1
365  40.0};//norm 2
366  */
367 
368 
369 
370  /* double norm= 15.3;;//14.85;
371  double mean_q0=0.39 ;// 0.36;
372  double mean_q3=0.95 ;//0.86;
373  double sigma_q0= 0.12 ;// 0.13;
374  double sigma_q3= 0.37;// 0.35;
375  double corr= 0.87;// 0.89;
376  double norm_2= 48.3;//42.0;
377  double mean_q0_2= 0.033;// 0.034;
378  double mean_q3_2= 0.44;//0.45;
379  double sigma_q0_2= 0.062;// 0.044;
380  double sigma_q3_2= 0.27;//0.31;
381  double corr_2= 0.77;// 0.75;
382  double baseline= -0.08;
383  */
384 
385  double norm = dg::norm;
386  double mean_q0 = dg::mean_q0;
387  double mean_q3 = dg::mean_q3;
388  double sigma_q0 = dg::sigma_q0;
389  double sigma_q3= dg::sigma_q3;
390  double corr= dg::corr;
391  double norm_2= dg::norm_2;
392  double mean_q0_2= dg::mean_q0_2;
393  double mean_q3_2= dg::mean_q3_2;
394  double sigma_q0_2= dg::sigma_q0_2;
395  double sigma_q3_2= dg::sigma_q3_2;
396  double corr_2= dg::corr_2;
397  double baseline= dg::baseline;
398 
399 
400  double expA = 1.0;
401  double expAlpha = 0.3 ;
402  double expLogC = 0.25;
403  double expLogB = 0.01;
404  // std::cout<<"The case here is "<<shift_param<<std::endl;
405  switch( shift_param )
406  {
407  case kGauss2DNorm_1:
408  norm += shift_sigma * dg::norm_shift;//0.5;
409  if ( norm < 0 ) norm = 0.0;
410  break;
411  case kGauss2DMeanQ0_1:
412  mean_q0 += shift_sigma * dg::mean_q0_shift;
413  if ( mean_q0 < 0 ) mean_q0 = 0.0;
414  break;
415  case kGauss2DMeanQ3_1:
416  mean_q3 += shift_sigma * dg::mean_q3_shift;
417  if ( mean_q3 < 0 ) mean_q0 = 0.0;
418  break;
419  case kGauss2DSigmaQ0_1:
420  sigma_q0 += shift_sigma * dg::sigma_q0_shift;
421  if ( sigma_q0 < 0.01 ) sigma_q0 = 0.01;
422  break;
423  case kGauss2DSigmaQ3_1:
424  sigma_q3 += shift_sigma * dg::sigma_q3_shift;
425  if ( sigma_q3 < 0.01 ) sigma_q3 = 0.01;
426  break;
427  case kGauss2DCorr_1 :
428  corr += shift_sigma * dg::corr_shift;
429  if ( corr < -0.99 ) corr = -0.99;
430  else if ( corr > 0.99 ) corr = 0.99;
431  break;
432  case kGauss2DNorm_2 :
433  norm_2 += shift_sigma * dg::norm_2_shift;
434  if ( norm_2 < 0 ) norm_2 = 0.0;
435  break;
436  case kGauss2DMeanQ0_2 :
437  mean_q0_2 += shift_sigma * dg::mean_q0_2_shift;
438  if ( mean_q0_2 < 0 ) mean_q0_2=0;
439  break;
440  case kGauss2DMeanQ3_2 :
441  mean_q3_2 += shift_sigma * dg::mean_q3_2_shift;
442  if ( mean_q3_2 < 0 ) mean_q3_2 = 0;
443  break;
444  case kGauss2DSigmaQ0_2 :
445  sigma_q0_2 += shift_sigma * dg::sigma_q0_2_shift;
446  if ( sigma_q0_2 < 0.01 ) sigma_q0_2 = 0.01;
447  break;
448  case kGauss2DSigmaQ3_2 :
449  sigma_q3_2 += shift_sigma * dg::sigma_q3_2_shift;
450  if ( sigma_q3_2 < 0.01 ) sigma_q3_2 = 0.01;
451  break;
452  case kGauss2DCorr_2 :
453  corr_2 += shift_sigma * dg::corr_2_shift;
454  if ( corr_2 < -0.99 ) corr_2 = -0.99;
455  if ( corr_2 > 0.99 ) corr_2 = 0.99;
456  break;
457  case kBaseline :
458  baseline += shift_sigma * dg::baseline_shift;
459  //if (baseline < 0.6) baseline=0.6;
460  break;
461  case kExpA:
462  expA += shift_sigma * 0.05;
463  if (expA <= 0 ) expA =0.01;
464  if (expA > 1) expA = 1.;
465  break;
466  case kExpAlpha:
467  expAlpha += shift_sigma * 0.025;
468  if (expAlpha < 0.1) expAlpha = 0.1;
469  break;
470  case kLogisticC:
471  expLogC += shift_sigma * 0.01;
472  break;
473  case kLogisticBeta:
474  expLogB += shift_sigma * 0.0025;
475  }
476 // reminder: http://mathworld.wolfram.com/BivariateNormalDistribution.html
477  double z = pow( ( q0 - mean_q0 ) / sigma_q0, 2 ) + pow( ( q3 - mean_q3 ) / sigma_q3, 2 )
478  - 2 * corr * ( q0 - mean_q0 ) * ( q3 - mean_q3 ) / ( sigma_q0 * sigma_q3 );
479 
480  double z_2 = pow( ( q0 - mean_q0_2 ) / sigma_q0_2, 2 ) + pow( ( q3 - mean_q3_2 ) / sigma_q3_2, 2 )
481  - 2 * corr_2 * ( q0 - mean_q0_2 ) * ( q3 - mean_q3_2 ) / ( sigma_q0_2 * sigma_q3_2 );
482 
483  // optional: Aaron's parametrization of suppression
484  //double supp = (1 - expA * exp(- Q2 / expAlpha)) / (1 + exp(-(q0 - expLogC)/expLogB));
485  double supp=1;
486  double weight= supp * ( baseline
487  + norm * exp( -0.5 * z / ( 1 - corr * corr ) )
488  + norm_2* exp( -0.5 * z_2 / ( 1 - corr_2 * corr_2 ) ) ) ;
489 
490  if (weight<0) { //std::cout<< "negative weight!!"<< weight<< std::endl;
491  weight=0;}
492  // std::cout<<"I'm adding a weight of "<<weight<< std::endl;
493  return weight;
494 
495 }
496 
497 /// cluncky way of having a separate calculation of up / down shifted versions
498 double CalcMECDoubleGaussEnhUP( const double q0, const double q3, const double Q2,
499  const MECDoubleGaussEnhParam shift_param, const double shift_sigma)
500 { // official parameters of prod5
501 
502  double norm=17.23;
503  double mean_q0= 0.31;
504  double mean_q3=0.85;
505  double sigma_q0= 0.087;
506  double sigma_q3= 0.45;
507  double corr= 0.95;
508  double norm_2=85.6;
509  double mean_q0_2= 0.039;
510  double mean_q3_2=0.46;
511  double sigma_q0_2= 0.047;
512  double sigma_q3_2=0.26;
513  double corr_2= 0.77;
514  double baseline= -0.19;
515 
516  double expA = 1.0;
517  double expAlpha = 0.3 ;
518  double expLogC = 0.25;
519  double expLogB = 0.01;
520 
521  switch( shift_param )
522  {
523  case kGauss2DNorm_1:
524  norm += shift_sigma * 0.5;
525  if ( norm < 0 ) norm = 0.0;
526  break;
527  case kGauss2DMeanQ0_1:
528  mean_q0 += shift_sigma * 0.025;
529  if ( mean_q0 < 0 ) mean_q0 = 0.0;
530  break;
531  case kGauss2DMeanQ3_1:
532  mean_q3 += shift_sigma * 0.05;
533  if ( mean_q3 < 0 ) mean_q0 = 0.0;
534  break;
535  case kGauss2DSigmaQ0_1:
536  sigma_q0 += shift_sigma * 0.015;
537  if ( sigma_q0 < 0.01 ) sigma_q0 = 0.01;
538  break;
539  case kGauss2DSigmaQ3_1:
540  sigma_q3 += shift_sigma * 0.02;
541  if ( sigma_q3 < 0.01 ) sigma_q3 = 0.01;
542  break;
543  case kGauss2DCorr_1 :
544  corr += shift_sigma * 0.2;
545  if ( corr < -0.99 ) corr = -0.99;
546  else if ( corr > 0.99 ) corr = 0.99;
547  break;
548  case kGauss2DNorm_2 :
549  norm_2 += shift_sigma * 2;
550  if ( norm_2 < 0 ) norm_2 = 0.0;
551  break;
552  case kGauss2DMeanQ0_2 :
553  mean_q0_2 += shift_sigma * 0.025;
554  if ( mean_q0_2 < 0 ) mean_q0_2=0;
555  break;
556  case kGauss2DMeanQ3_2 :
557  mean_q3_2 += shift_sigma * 0.05;
558  if ( mean_q3_2 < 0 ) mean_q3_2 = 0;
559  break;
560  case kGauss2DSigmaQ0_2 :
561  sigma_q0_2 += shift_sigma * 0.015;
562  if ( sigma_q0_2 < 0.01 ) sigma_q0_2 = 0.01;
563  break;
564  case kGauss2DSigmaQ3_2 :
565  sigma_q3_2 += shift_sigma * 0.02;
566  if ( sigma_q3_2 < 0.01 ) sigma_q3_2 = 0.01;
567  break;
568  case kGauss2DCorr_2 :
569  corr_2 += shift_sigma * 0.2;
570  if ( corr_2 < -0.99 ) corr_2 = -0.99;
571  if ( corr_2 > 0.99 ) corr_2 = 0.99;
572  break;
573  case kBaseline :
574  baseline += shift_sigma*0.05;
575  //if (baseline < 0.6) baseline=0.6;
576  break;
577  case kExpA:
578  expA += shift_sigma * 0.05;
579  if (expA <= 0 ) expA =0.01;
580  if (expA > 1) expA = 1.;
581  break;
582  case kExpAlpha:
583  expAlpha += shift_sigma * 0.025;
584  if (expAlpha < 0.1) expAlpha = 0.1;
585  break;
586  case kLogisticC:
587  expLogC += shift_sigma * 0.01;
588  break;
589  case kLogisticBeta:
590  expLogB += shift_sigma * 0.0025;
591  }
592 // reminder: http://mathworld.wolfram.com/BivariateNormalDistribution.html
593  double z = pow( ( q0 - mean_q0 ) / sigma_q0, 2 ) + pow( ( q3 - mean_q3 ) / sigma_q3, 2 )
594  - 2 * corr * ( q0 - mean_q0 ) * ( q3 - mean_q3 ) / ( sigma_q0 * sigma_q3 );
595 
596  double z_2 = pow( ( q0 - mean_q0_2 ) / sigma_q0_2, 2 ) + pow( ( q3 - mean_q3_2 ) / sigma_q3_2, 2 )
597  - 2 * corr_2 * ( q0 - mean_q0_2 ) * ( q3 - mean_q3_2 ) / ( sigma_q0_2 * sigma_q3_2 );
598 
599  // optional: Aaron's parametrization of suppression
600  //double supp = (1 - expA * exp(- Q2 / expAlpha)) / (1 + exp(-(q0 - expLogC)/expLogB));
601  double supp=1;
602  double weight= supp * ( baseline
603  + norm * exp( -0.5 * z / ( 1 - corr * corr ) )
604  + norm_2* exp( -0.5 * z_2 / ( 1 - corr_2 * corr_2 ) ) ) ;
605 
606  if (weight<0) { //std::cout<< "negative weight!!"<< weight<< std::endl;
607  weight=0;}
608  return weight;
609 
610 }
611 
612 double CalcMECDoubleGaussEnhDOWN( const double q0, const double q3, const double Q2,
613  const MECDoubleGaussEnhParam shift_param, const double shift_sigma)
614 { // official parameters of prod5
615 
616  double norm=25.4;
617  double mean_q0= 0.48;
618  double mean_q3=1.00;
619  double sigma_q0= 0.14;
620  double sigma_q3= 0.36;
621  double corr= 0.99;
622  double norm_2=109.7;
623  double mean_q0_2= 0.03;
624  double mean_q3_2=0.44;
625  double sigma_q0_2= 0.04;
626  double sigma_q3_2=0.20;
627  double corr_2= 0.78;
628  double baseline= 0.48;
629 
630  double expA = 1.0;
631  double expAlpha = 0.3 ;
632  double expLogC = 0.25;
633  double expLogB = 0.01;
634 
635  switch( shift_param )
636  {
637  case kGauss2DNorm_1:
638  norm += shift_sigma * 0.5;
639  if ( norm < 0 ) norm = 0.0;
640  break;
641  case kGauss2DMeanQ0_1:
642  mean_q0 += shift_sigma * 0.025;
643  if ( mean_q0 < 0 ) mean_q0 = 0.0;
644  break;
645  case kGauss2DMeanQ3_1:
646  mean_q3 += shift_sigma * 0.05;
647  if ( mean_q3 < 0 ) mean_q0 = 0.0;
648  break;
649  case kGauss2DSigmaQ0_1:
650  sigma_q0 += shift_sigma * 0.015;
651  if ( sigma_q0 < 0.01 ) sigma_q0 = 0.01;
652  break;
653  case kGauss2DSigmaQ3_1:
654  sigma_q3 += shift_sigma * 0.02;
655  if ( sigma_q3 < 0.01 ) sigma_q3 = 0.01;
656  break;
657  case kGauss2DCorr_1 :
658  corr += shift_sigma * 0.2;
659  if ( corr < -0.99 ) corr = -0.99;
660  else if ( corr > 0.99 ) corr = 0.99;
661  break;
662  case kGauss2DNorm_2 :
663  norm_2 += shift_sigma * 2;
664  if ( norm_2 < 0 ) norm_2 = 0.0;
665  break;
666  case kGauss2DMeanQ0_2 :
667  mean_q0_2 += shift_sigma * 0.025;
668  if ( mean_q0_2 < 0 ) mean_q0_2=0;
669  break;
670  case kGauss2DMeanQ3_2 :
671  mean_q3_2 += shift_sigma * 0.05;
672  if ( mean_q3_2 < 0 ) mean_q3_2 = 0;
673  break;
674  case kGauss2DSigmaQ0_2 :
675  sigma_q0_2 += shift_sigma * 0.015;
676  if ( sigma_q0_2 < 0.01 ) sigma_q0_2 = 0.01;
677  break;
678  case kGauss2DSigmaQ3_2 :
679  sigma_q3_2 += shift_sigma * 0.02;
680  if ( sigma_q3_2 < 0.01 ) sigma_q3_2 = 0.01;
681  break;
682  case kGauss2DCorr_2 :
683  corr_2 += shift_sigma * 0.2;
684  if ( corr_2 < -0.99 ) corr_2 = -0.99;
685  if ( corr_2 > 0.99 ) corr_2 = 0.99;
686  break;
687  case kBaseline :
688  baseline += shift_sigma*0.05;
689  //if (baseline < 0.6) baseline=0.6;
690  break;
691  case kExpA:
692  expA += shift_sigma * 0.05;
693  if (expA <= 0 ) expA =0.01;
694  if (expA > 1) expA = 1.;
695  break;
696  case kExpAlpha:
697  expAlpha += shift_sigma * 0.025;
698  if (expAlpha < 0.1) expAlpha = 0.1;
699  break;
700  case kLogisticC:
701  expLogC += shift_sigma * 0.01;
702  break;
703  case kLogisticBeta:
704  expLogB += shift_sigma * 0.0025;
705  }
706 // reminder: http://mathworld.wolfram.com/BivariateNormalDistribution.html
707  double z = pow( ( q0 - mean_q0 ) / sigma_q0, 2 ) + pow( ( q3 - mean_q3 ) / sigma_q3, 2 )
708  - 2 * corr * ( q0 - mean_q0 ) * ( q3 - mean_q3 ) / ( sigma_q0 * sigma_q3 );
709 
710  double z_2 = pow( ( q0 - mean_q0_2 ) / sigma_q0_2, 2 ) + pow( ( q3 - mean_q3_2 ) / sigma_q3_2, 2 )
711  - 2 * corr_2 * ( q0 - mean_q0_2 ) * ( q3 - mean_q3_2 ) / ( sigma_q0_2 * sigma_q3_2 );
712 
713  // optional: Aaron's parametrization of suppression
714  //double supp = (1 - expA * exp(- Q2 / expAlpha)) / (1 + exp(-(q0 - expLogC)/expLogB));
715  double supp=1;
716  double weight= supp * ( baseline
717  + norm * exp( -0.5 * z / ( 1 - corr * corr ) )
718  + norm_2* exp( -0.5 * z_2 / ( 1 - corr_2 * corr_2 ) ) ) ;
719 
720  if (weight<0) { //std::cout<< "negative weight!!"<< weight<< std::endl;
721  weight=0;}
722  return weight;
723 
724 }
725 
726 
727 const std::map< std::string , MECDoubleGaussEnhParam > DoubleGaussMap{
728 {"MECDoubleGaussEnhSystNorm_1", kGauss2DNorm_1},
729 {"MECDoubleGaussEnhSystMeanQ0_1", kGauss2DMeanQ0_1},
730 {"MECDoubleGaussEnhSystMeanQ3_1", kGauss2DMeanQ3_1},
731 {"MECDoubleGaussEnhSystSigmaQ0_1", kGauss2DSigmaQ0_1},
732 {"MECDoubleGaussEnhSystSigmaQ3_1", kGauss2DSigmaQ3_1 },
733 {"MECDoubleGaussEnhSystCorr_1", kGauss2DCorr_1 },
734 {"MECDoubleGaussEnhSystNorm_2", kGauss2DNorm_2 },
735 {"MECDoubleGaussEnhSystMeanQ0_2", kGauss2DMeanQ0_2},
736 {"MECDoubleGaussEnhSystMeanQ3_2", kGauss2DMeanQ3_2},
737 {"MECDoubleGaussEnhSystSigmaQ0_2", kGauss2DSigmaQ0_2},
738 {"MECDoubleGaussEnhSystSigmaQ3_2", kGauss2DSigmaQ3_2},
739 {"MECDoubleGaussEnhSystCorr_2", kGauss2DCorr_2 },
740 {"MECDoubleGaussEnhSystBaseline", kBaseline},
741 
742 {"MECDoubleGaussEnhSystUPNorm_1", kGauss2DNorm_1},
743 {"MECDoubleGaussEnhSystUPMeanQ0_1", kGauss2DMeanQ0_1},
744 {"MECDoubleGaussEnhSystUPMeanQ3_1", kGauss2DMeanQ3_1},
745 {"MECDoubleGaussEnhSystUPSigmaQ0_1", kGauss2DSigmaQ0_1},
746 {"MECDoubleGaussEnhSystUPSigmaQ3_1", kGauss2DSigmaQ3_1 },
747 {"MECDoubleGaussEnhSystUPCorr_1", kGauss2DCorr_1 },
748 {"MECDoubleGaussEnhSystUPNorm_2", kGauss2DNorm_2 },
749 {"MECDoubleGaussEnhSystUPMeanQ0_2", kGauss2DMeanQ0_2},
750 {"MECDoubleGaussEnhSystUPMeanQ3_2", kGauss2DMeanQ3_2},
751 {"MECDoubleGaussEnhSystUPSigmaQ0_2", kGauss2DSigmaQ0_2},
752 {"MECDoubleGaussEnhSystUPSigmaQ3_2", kGauss2DSigmaQ3_2},
753 {"MECDoubleGaussEnhSystUPCorr_2", kGauss2DCorr_2 },
754 {"MECDoubleGaussEnhSystUPBaseline", kBaseline},
755 
756 {"MECDoubleGaussEnhSystDOWNNorm_1", kGauss2DNorm_1},
757 {"MECDoubleGaussEnhSystDOWNMeanQ0_1", kGauss2DMeanQ0_1},
758 {"MECDoubleGaussEnhSystDOWNMeanQ3_1", kGauss2DMeanQ3_1},
759 {"MECDoubleGaussEnhSystDOWNSigmaQ0_1", kGauss2DSigmaQ0_1},
760 {"MECDoubleGaussEnhSystDOWNSigmaQ3_1", kGauss2DSigmaQ3_1 },
761 {"MECDoubleGaussEnhSystDOWNCorr_1", kGauss2DCorr_1 },
762 {"MECDoubleGaussEnhSystDOWNNorm_2", kGauss2DNorm_2 },
763 {"MECDoubleGaussEnhSystDOWNMeanQ0_2", kGauss2DMeanQ0_2},
764 {"MECDoubleGaussEnhSystDOWNMeanQ3_2", kGauss2DMeanQ3_2},
765 {"MECDoubleGaussEnhSystDOWNSigmaQ0_2", kGauss2DSigmaQ0_2},
766 {"MECDoubleGaussEnhSystDOWNSigmaQ3_2", kGauss2DSigmaQ3_2},
767 {"MECDoubleGaussEnhSystDOWNCorr_2", kGauss2DCorr_2 },
768 {"MECDoubleGaussEnhSystDOWNBaseline", kBaseline},
769 
770 {"MECDoubleGaussEnhSystExpA",kExpA},
771 {"MECDoubleGaussEnhSystExpAlpha",kExpAlpha},
772 {"MECDoubleGaussEnhSystLogisticC",kLogisticC},
773 {"MECDoubleGaussEnhSystLogisticBeta",kLogisticBeta},
774 
775 
776 };
777 
778 double CalcMECDoubleGaussEnhShiftedParam( std::string shift_param, double shift_sigma, std::string shifted)
779 {
780 
781  MECDoubleGaussEnhParam shifted_param = DoubleGaussMap.find(shift_param)->second;
782 // double norm, mean_q0, mean_q3, sigma_q0, sigma_q3, corr;
783 // double norm_2, mean_q0_2, mean_q3_2, sigma_q0_2, sigma_q3_2, corr_2, baseline;
784 //
785 // if (shifted=="down"){
786 // //Down
787 // norm=25.4;
788 // mean_q0= 0.48;
789 // mean_q3=1.00;
790 // sigma_q0= 0.14;
791 // sigma_q3= 0.36;
792 // corr= 0.99;
793 // norm_2=109.7;
794 // mean_q0_2= 0.03;
795 // mean_q3_2=0.44;
796 // sigma_q0_2= 0.04;
797 // sigma_q3_2=0.20;
798 // corr_2= 0.78;
799 // baseline= 0.48;
800 //} else if(shifted=="up"){
801 // // up
802 //norm=17.23;
803 //mean_q0= 0.31;
804 //mean_q3=0.85;
805 //sigma_q0= 0.087;
806 //sigma_q3= 0.45;
807 //corr= 0.95;
808 //norm_2=85.6;
809 //mean_q0_2= 0.039;
810 //mean_q3_2=0.46;
811 //sigma_q0_2= 0.047;
812 //sigma_q3_2=0.26;
813 //corr_2= 0.77;
814 //baseline= -0.19;
815 //} else {
816  // official parameters of prod5
817  /* norm= 15.3;;//14.85;
818  mean_q0=0.39 ;// 0.36;
819  mean_q3=0.95 ;//0.86;
820  sigma_q0= 0.12 ;// 0.13;
821  sigma_q3= 0.37;// 0.35;
822  corr= 0.87;// 0.89;
823  norm_2= 48.3;//42.0;
824  mean_q0_2= 0.033;// 0.034;
825  mean_q3_2= 0.44;//0.45;
826  sigma_q0_2= 0.062;// 0.044;
827  sigma_q3_2= 0.27;//0.31;
828  corr_2= 0.77;// 0.75;
829  baseline= -0.08;
830  */
831 
832 // norm= 10.0;//14.85;
833 // mean_q0=1.05;//0.39 ;// 0.36; //0.161587;//0.272638;//0.383689;//0.49474;//0.605791;//0.716842;//0.827893;//0.938944 ;
834 // mean_q3=0.4;//0.95 ;//0.86;
835 // sigma_q0= 0.1;//0.12 ;// 0.13;
836 // sigma_q3= 0.04;//0.37;// 0.35;
837 // corr= 0.9;// 0.89;
838 // norm_2=40.0;//42.0;
839 // mean_q0_2= 0.19;// 0.034;
840 // mean_q3_2= 0.55;//0.45;
841 // sigma_q0_2= 0.1;// 0.044;
842 // sigma_q3_2= 0.1;//0.31;
843 // corr_2= 0.9;// 0.75;
844 // baseline= 0.;//-0.12;// -0.08;
845 
846 
847 
848  double norm = dg::norm;
849  double mean_q0 = dg::mean_q0;
850  double mean_q3 = dg::mean_q3;
851  double sigma_q0 = dg::sigma_q0;
852  double sigma_q3= dg::sigma_q3;
853  double corr= dg::corr;
854  double norm_2= dg::norm_2;
855  double mean_q0_2= dg::mean_q0_2;
856  double mean_q3_2= dg::mean_q3_2;
857  double sigma_q0_2= dg::sigma_q0_2;
858  double sigma_q3_2= dg::sigma_q3_2;
859  double corr_2= dg::corr_2;
860  double baseline= dg::baseline;
861 
862 
863  double expA = 1.0;
864  double expAlpha = 0.3 ;
865  double expLogC = 0.25;
866  double expLogB = 0.01;
867 
868  double param = 0.0;
869 
870  switch( shifted_param )
871  {
872  case kGauss2DNorm_1:
873  norm += shift_sigma * dg::norm_shift;//0.5; //.5
874  if ( norm < 0 ) norm = 0.0;
875  param = norm;
876  break;
877  case kGauss2DMeanQ0_1:
878  mean_q0 += shift_sigma * dg::mean_q0_shift;
879  if ( mean_q0 < 0 ) mean_q0 = 0.0;
880  param = mean_q0;
881  break;
882  case kGauss2DMeanQ3_1:
883  mean_q3 += shift_sigma * dg::mean_q3_shift;
884  if ( mean_q3 < 0 ) mean_q3 = 0.0;
885  param = mean_q3;
886  break;
887  case kGauss2DSigmaQ0_1:
888  sigma_q0 += shift_sigma * dg::sigma_q0_shift;
889  if ( sigma_q0 < 0.001 ) sigma_q0 = 0.001;
890  param = sigma_q0;
891  break;
892  case kGauss2DSigmaQ3_1:
893  sigma_q3 += shift_sigma * dg::sigma_q3_shift;
894  if ( sigma_q3 < 0.001 ) sigma_q3 = 0.001;
895  param = sigma_q3;
896  break;
897  case kGauss2DCorr_1 :
898  corr += shift_sigma * dg::corr_shift;
899  if ( corr < -0.99 ) corr = -0.99;
900  else if ( corr > 0.99 ) corr = 0.99;
901  param = corr;
902  break;
903  case kGauss2DNorm_2 :
904  norm_2 += shift_sigma * dg::norm_2_shift;
905  if ( norm_2 < 0 ) norm_2 = 0.0;
906  param = norm_2;
907  break;
908  case kGauss2DMeanQ0_2 :
909  mean_q0_2 += shift_sigma * dg::mean_q0_2_shift;
910  if ( mean_q0_2 < 0 ) mean_q0_2 = 0.0;
911  param = mean_q0_2;
912  break;
913  case kGauss2DMeanQ3_2 :
914  mean_q3_2 += shift_sigma * dg::mean_q3_2_shift;
915  if ( mean_q3_2 < 0 ) mean_q3_2 = 0.0;
916  param = mean_q3_2;
917  break;
918  case kGauss2DSigmaQ0_2 :
919  sigma_q0_2 += shift_sigma * dg::sigma_q0_2_shift;
920  if ( sigma_q0_2 < 0.001 ) sigma_q0_2 = 0.001;
921  param = sigma_q0_2;
922  break;
923  case kGauss2DSigmaQ3_2 :
924  sigma_q3_2 += shift_sigma * dg::sigma_q3_2_shift;
925  if ( sigma_q3_2 < 0.001 ) sigma_q3_2 = 0.001;
926  param = sigma_q3_2;
927  break;
928  case kGauss2DCorr_2 :
929  corr_2 += shift_sigma * dg::corr_2_shift;
930  if ( corr_2 < -0.99 ) corr_2 = -0.99;
931  if ( corr_2 > 0.99 ) corr_2 = 0.99;
932  param = corr_2;
933  break;
934  case kBaseline :
935  baseline += shift_sigma * dg::baseline_shift;
936  // if (baseline < 0.6) baseline=0.6;
937  param = baseline;
938  break;
939  case kExpA:
940  expA += shift_sigma * 0.05;
941  if (expA <= 0 ) expA = 0.01;
942  if (expA > 1) expA = 1.;
943  param = expA;
944  break;
945  case kExpAlpha:
946  expAlpha += shift_sigma * 0.025;
947  if (expAlpha < 0.001) expAlpha = 0.001;
948  param = expAlpha;
949  break;
950  case kLogisticC:
951  expLogC += shift_sigma * 0.01;
952  param = expLogC;
953  break;
954  case kLogisticBeta:
955  expLogB += shift_sigma * 0.0025;
956  param = expLogB;
957  }
958  return param;
959 }
960 
961 
962 const Var kMECDoubleGaussEnh( []( const caf::SRProxy* sr )
963 {
964 // if ( !kIsNumuCC( sr ) || sr->mc.nu[0].mode != caf::kMEC ) return 1.0;
965  if ( sr->mc.nnu < 1) return 1.0;
966  if ( !(sr->mc.nu[0].iscc) || sr->mc.nu[0].mode != caf::kMEC ) return 1.0;
967  double q0 = kTrueQ0( sr );
968  double q3 = kTrueQ3( sr );
969  double Q2 = kTrueQ2( sr );
970  return CalcMECDoubleGaussEnh( q0, q3, Q2, kGauss2DNorm_1, 0 ); // Shift = 0 gives nominal suppression
971 });
972 
973 const Var kMECDoubleGaussEnhUP( []( const caf::SRProxy* sr )
974 {
975 // if ( !kIsNumuCC( sr ) || sr->mc.nu[0].mode != caf::kMEC ) return 1.0;
976  if ( sr->mc.nnu < 1) return 1.0;
977  if ( !(sr->mc.nu[0].iscc) || sr->mc.nu[0].mode != caf::kMEC ) return 1.0;
978  double q0 = kTrueQ0( sr );
979  double q3 = kTrueQ3( sr );
980  double Q2 = kTrueQ2( sr );
981  return CalcMECDoubleGaussEnhUP( q0, q3, Q2, kGauss2DNorm_1, 0 ); // Shift = 0 gives nominal suppression
982 });
983 
984 const Var kMECDoubleGaussEnhDOWN( []( const caf::SRProxy* sr )
985 {
986 // if ( !kIsNumuCC( sr ) || sr->mc.nu[0].mode != caf::kMEC ) return 1.0;
987  if ( sr->mc.nnu < 1) return 1.0;
988  if ( !(sr->mc.nu[0].iscc) || sr->mc.nu[0].mode != caf::kMEC ) return 1.0;
989  double q0 = kTrueQ0( sr );
990  double q3 = kTrueQ3( sr );
991  double Q2 = kTrueQ2( sr );
992  return CalcMECDoubleGaussEnhDOWN( q0, q3, Q2, kGauss2DNorm_1, 0 ); // Shift = 0 gives nominal suppression
993 });
994 
995 
996 class MECDoubleGaussEnhSyst : public ISyst
997 {
998  public :
999  MECDoubleGaussEnhSyst( const MECDoubleGaussEnhParam& shift_param, const std::string& shift_param_name )
1000  : ISyst( "MECDoubleGaussEnhSyst" + shift_param_name, "MEC 2D Gauss Syst " + shift_param_name ),
1001  fShiftParam( shift_param )
1002  {}
1003  void Shift( double sigma, caf::SRProxy* sr, double& weight ) const override
1004  {
1005 // if ( !kIsNumuCC( sr ) || sr->mc.nu[0].mode != caf::kMEC ) return;
1006  if ( sr->mc.nnu < 1) return;
1007  if ( !(sr->mc.nu[0].iscc) || sr->mc.nu[0].mode != caf::kMEC ) return;
1008  // This will work for Empirical MEC. Will it also work for Valencia MEC?
1009  double q0 = kTrueQ0( sr );
1010  double q3 = kTrueQ3( sr );
1011  double Q2 = kTrueQ2( sr );
1012  double wgt_nominal = CalcMECDoubleGaussEnh( q0, q3, Q2, fShiftParam, 0 );
1013  double wgt_shift = CalcMECDoubleGaussEnh( q0, q3, Q2, fShiftParam, sigma );
1014  weight *= wgt_shift / wgt_nominal;
1015  if (wgt_nominal <= 0 ) weight = 0;
1016  std::cout<<"The weight here is "weight<<std::endl;
1017 
1018  }
1019  private:
1020  const MECDoubleGaussEnhParam fShiftParam;
1021 };
1022 
1023 // Clunky again, UP/DOWN shifted
1024 
1025 class MECDoubleGaussEnhSystUP : public ISyst
1026 {
1027  public :
1028  MECDoubleGaussEnhSystUP( const MECDoubleGaussEnhParam& shift_param, const std::string& shift_param_name )
1029  : ISyst( "MECDoubleGaussEnhSystUP" + shift_param_name, "MEC 2D Gauss Syst (up)" + shift_param_name ),
1030  fShiftParam( shift_param )
1031  {}
1032  void Shift( double sigma, caf::SRProxy* sr, double& weight ) const override
1033  {
1034 // if ( !kIsNumuCC( sr ) || sr->mc.nu[0].mode != caf::kMEC ) return;
1035  if ( sr->mc.nnu < 1) return;
1036  if ( !(sr->mc.nu[0].iscc) || sr->mc.nu[0].mode != caf::kMEC ) return;
1037  // This will work for Empirical MEC. Will it also work for Valencia MEC?
1038  double q0 = kTrueQ0( sr );
1039  double q3 = kTrueQ3( sr );
1040  double Q2 = kTrueQ2( sr );
1041  double wgt_nominal = CalcMECDoubleGaussEnhUP( q0, q3, Q2, fShiftParam, 0 );
1042  double wgt_shift = CalcMECDoubleGaussEnhUP( q0, q3, Q2, fShiftParam, sigma );
1043  weight *= wgt_shift / wgt_nominal;
1044  if (wgt_nominal <= 0 ) weight = 0;
1045 
1046  }
1047  private:
1048  const MECDoubleGaussEnhParam fShiftParam;
1049 };
1050 
1051 
1052 class MECDoubleGaussEnhSystDOWN : public ISyst
1053 {
1054  public :
1055  MECDoubleGaussEnhSystDOWN( const MECDoubleGaussEnhParam& shift_param, const std::string& shift_param_name )
1056  : ISyst( "MECDoubleGaussEnhSystDOWN" + shift_param_name, "MEC 2D Gauss Syst (down)" + shift_param_name ),
1057  fShiftParam( shift_param )
1058  {}
1059  void Shift( double sigma, caf::SRProxy* sr, double& weight ) const override
1060  {
1061 // if ( !kIsNumuCC( sr ) || sr->mc.nu[0].mode != caf::kMEC ) return;
1062  if ( sr->mc.nnu < 1) return;
1063  if ( !(sr->mc.nu[0].iscc) || sr->mc.nu[0].mode != caf::kMEC ) return;
1064  // This will work for Empirical MEC. Will it also work for Valencia MEC?
1065  double q0 = kTrueQ0( sr );
1066  double q3 = kTrueQ3( sr );
1067  double Q2 = kTrueQ2( sr );
1068  double wgt_nominal = CalcMECDoubleGaussEnhDOWN( q0, q3, Q2, fShiftParam, 0 );
1069  double wgt_shift = CalcMECDoubleGaussEnhDOWN( q0, q3, Q2, fShiftParam, sigma );
1070  weight *= wgt_shift / wgt_nominal;
1071  if (wgt_nominal <= 0 ) weight = 0;
1072 
1073  }
1074  private:
1075  const MECDoubleGaussEnhParam fShiftParam;
1076 };
1077 
1078 
1079 //===============================================================================================================
1080 // MINOS Resonance Suppression
1081 //===============================================================================================================
1082 
1083 enum MinosResSuppParam{ kMinosResSuppNorm, kMinosResSuppQ0 };
1084 
1085 double CalcMinosResSupp( const double Q2, const MinosResSuppParam shift_param, const double shift_sigma )
1086 {
1087  double A = 1.010;
1088  double Q0 = 0.156;
1089 
1090  switch( shift_param )
1091  {
1092  case kMinosResSuppNorm :
1093  A += shift_sigma * 0.1;
1094  if ( A < 0 ) A = 0.0;
1095  break;
1096  case kMinosResSuppQ0 :
1097  Q0 += shift_sigma * 0.03;
1098  if ( Q0 < 0.01 ) Q0 = 0.01;
1099  }
1100 
1101  double supp = A / ( 1 + exp( 1 - sqrt( Q2 ) / Q0 ) );
1102  if ( supp > 1 ) supp = 1.0;
1103 
1104  return supp;
1105 }
1106 
1107 const Var kMinosResSupp( []( const caf::SRProxy* sr )
1108 {
1109  if ( !kIsNumuCC( sr ) || sr->mc.nu[0].mode != kIsRes( sr ) || sr->mc.nu[0].tgtA == 1 ) return 1.0; // Exclude hydrogen
1110  double Q2 = kTrueQ2( sr );
1111  return CalcMinosResSupp( Q2, kMinosResSuppNorm, 0 ); // Shift = 0 gives nominal suppression
1112 });
1113 
1114 class MinosResSuppSyst : public ISyst
1115 {
1116  public :
1117  MinosResSuppSyst( const MinosResSuppParam& shift_param, const std::string& shift_param_name )
1118  : ISyst( "MinosResSuppSyst" + shift_param_name, "MINOS RES Supp Syst " + shift_param_name ),
1119  fShiftParam( shift_param )
1120  {}
1121 
1122  void Shift( double sigma, caf::SRProxy* sr, double& weight ) const override
1123  {
1124  if ( !kIsNumuCC( sr ) || sr->mc.nu[0].mode != kIsRes( sr ) || sr->mc.nu[0].tgtA == 1 ) return; // Exclude hydrogen
1125 
1126  double Q2 = kTrueQ2( sr );
1127 
1128  double wgt_nominal = CalcMinosResSupp( Q2, fShiftParam, 0 );
1129  double wgt_shift = CalcMinosResSupp( Q2, fShiftParam, sigma );
1130 
1131  weight *= wgt_shift / wgt_nominal;
1132  }
1133 
1134  private:
1135  const MinosResSuppParam fShiftParam;
1136 };
1137 
1138 //=============================================================================
1139 // Extra Weight
1140 //=============================================================================
1141 
1142 const Var kWeightPionDeuteriumTune( []( const caf::SRProxy* sr )
1143 {
1144  double weight = 1.0;
1145 
1146  if ( kIsNumuCC( sr ) && kIsRes( sr ) )
1147  {
1148  if ( sr->mc.nu[0].rwgt.genie.size() <= rwgt::fReweightMaCCRES - 1 )
1149  {
1150  if ( sr->mc.nu[0].isvtxcont )
1151  throw std::runtime_error( "kWeightPionDeuteriumTune: Cannot do MA CCRES rescaling without GENIE reweights available." );
1152  else
1153  return 1.0;
1154  }
1155  // GENIE single pion production deuterium fit Eur. Phys. J. C (2016) 76: 474
1156  // Scale down MA CCRES from GENIE default 1.12 GeV (+/-1 sigma = +/-20%) to Deuterium fit 0.94+/-0.05 GeV
1157  const double correctionInSigma = ( 1.12 - 0.94 ) / 0.2;
1158  weight *= 1. + correctionInSigma * ( sr->mc.nu[0].rwgt.genie[rwgt::fReweightMaCCRES].minus1sigma - 1. );
1159 
1160  // Apply 1.15 normalization correction from Deuterium fit
1161  weight *= 1.15;
1162  }
1163  else if ( kIsNumuCC( sr ) && kIsDIS( sr ) && kTrueW( sr ) < 1.7 )
1164  {
1165  // Apply non-RES 1pi normalization from deuterium fit
1166  if ( sr->mc.nu[0].npiminus + sr->mc.nu[0].npiplus + sr->mc.nu[0].npizero == 1 )
1167  {
1168  weight *= 0.43;
1169  }
1170  }
1171  return weight;
1172 });
1173 
1174 }
1175 
caf::Proxy< caf::SRIDBranch >::cvnloosepreselptp
caf::Proxy< caf::SRCVNResult > cvnloosepreselptp
Definition: SRProxy.h:1254
ana::MECDoubleGaussEnhParam
MECDoubleGaussEnhParam
Definition: MECTuningSysts.h:14
name
const XML_Char * name
Definition: expat.h:151
plot_validation_datamc.pred
pred
Definition: plot_validation_datamc.py:175
ana::kIsNumuCC
const Cut kIsNumuCC(CCFlavSel(14, 14))
Select CC .
ana::MinosResSuppSyst::MinosResSuppSyst
MinosResSuppSyst(const MinosResSuppParam &shift_param, const std::string &shift_param_name)
Definition: MECTuningUtils.h:1117
ana::TrivialPrediction::TrivialPrediction
TrivialPrediction(const Spectrum &spec)
Definition: MECTuningUtils.h:58
dg::corr
double corr
Definition: DoubleGaussParams.h:11
ana
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
ana::MECDoubleGaussEnhSystDOWN::fShiftParam
const MECDoubleGaussEnhParam fShiftParam
Definition: MECTuningUtils.h:1075
ana::MECDoubleGaussEnhSystDOWN::MECDoubleGaussEnhSystDOWN
MECDoubleGaussEnhSystDOWN(const MECDoubleGaussEnhParam &shift_param, const std::string &shift_param_name)
Definition: MECTuningUtils.h:1055
dg::mean_q0_2_shift
double mean_q0_2_shift
Definition: DoubleGaussParams.h:30
genie::utils::kinematics::Q2
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:991
ana::kMECDoubleGaussEnhDOWN
const Var kMECDoubleGaussEnhDOWN([](const caf::SRProxy *sr){if(sr->mc.nnu< 1) return 1.0;if(!(sr->mc.nu[0].iscc)||sr->mc.nu[0].mode!=caf::kMEC) return 1.0;double q0=kTrueQ0(sr);double q3=kTrueQ3(sr);double Q2=kTrueQ2(sr);return CalcMECDoubleGaussEnhDOWN(q0, q3, Q2, kGauss2DNorm_1, 0);})
ana::SingleSampleExperiment::ChiSq
virtual double ChiSq(osc::IOscCalcAdjustable *osc, const SystShifts &syst=SystShifts::Nominal()) const override
Definition: SingleSampleExperiment.cxx:111
ana::weight
const Var weight
Definition: PIDCutTuningDefs.h:96
make_mec_shifts_plots.isRHC
string isRHC
Definition: make_mec_shifts_plots.py:100
dg::sigma_q0_shift
double sigma_q0_shift
Definition: DoubleGaussParams.h:26
ana::SystShifts
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
ana::BigChi2SingleSampleExperiment::ChiSq
double ChiSq(osc::IOscCalcAdjustable *osc, const SystShifts &syst=SystShifts::Nominal()) const override
Definition: MECTuningUtils.h:159
ana::kIsRes
const Cut kIsRes
Definition: TruthCuts.cxx:111
dg::mean_q3_2_shift
double mean_q3_2_shift
Definition: DoubleGaussParams.h:31
caf::Proxy< caf::StandardRecord >
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
ana::Loaders
Collection of SpectrumLoaders for many configurations.
Definition: Loaders.h:23
std::sqrt
T sqrt(T number)
Definition: d0nt_math.hpp:156
osc::_IOscCalc
General interface to oscillation calculators.
Definition: StanTypedefs.h:22
caf::Proxy< caf::SRTruthBranch >::nu
caf::Proxy< std::vector< caf::SRNeutrino > > nu
Definition: SRProxy.h:618
cet::pow
constexpr T pow(T x)
Definition: pow.h:75
kMC
const Color_t kMC
Definition: make_extrap_figure.C:44
dg::sigma_q3
double sigma_q3
Definition: DoubleGaussParams.h:10
caf::Proxy< caf::SRRemid >::pid
caf::Proxy< float > pid
Definition: SRProxy.h:1136
dg::corr_2_shift
double corr_2_shift
Definition: DoubleGaussParams.h:34
cmf::Spectrum
std::vector< double > Spectrum
Definition: Constants.h:741
dg::mean_q3_shift
double mean_q3_shift
Definition: DoubleGaussParams.h:25
ana::CalcMECDoubleGaussEnhShiftedParam
double CalcMECDoubleGaussEnhShiftedParam(std::string shift_param, double shift_sigma, const GaussParams initialParams)
Definition: MECTuningSysts.cxx:162
caf::Proxy< caf::SRTruthBranch >::nnu
caf::Proxy< short int > nnu
Definition: SRProxy.h:617
ana::kIsDIS
const Cut kIsDIS
Definition: TruthCuts.cxx:118
tmp
Float_t tmp
Definition: plot.C:36
ana::MECDoubleGaussEnhSyst::MECDoubleGaussEnhSyst
MECDoubleGaussEnhSyst(const MECDoubleGaussEnhParam &shift_param, const std::string &shift_param_name)
Definition: MECTuningUtils.h:999
ana::kTrueQ0
const Var kTrueQ0
Definition: TruthVars.h:32
ana::kMECGaussEnh
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);})
ana::SystShifts::Nominal
static SystShifts Nominal()
Definition: SystShifts.h:34
dg::baseline
double baseline
Definition: DoubleGaussParams.h:18
calc
osc::OscCalcDumb calc
Definition: make_nueFDprediction_kinematics.C:70
ana::ISyst
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
dg::sigma_q0_2
double sigma_q0_2
Definition: DoubleGaussParams.h:15
ana::kMinosResSupp
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);})
ana::kTrueQ2
const Var kTrueQ2
Definition: TruthVars.h:27
ana::Var
_Var< caf::SRProxy > Var
Representation of a variable to be retrieved from a caf::StandardRecord object.
Definition: Var.h:74
ana::Spectrum
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
ana::MECGaussEnhParam
MECGaussEnhParam
Definition: MECTuningUtils.h:245
ana::MECDoubleGaussEnhSyst::Shift
void Shift(double sigma, caf::SRProxy *sr, double &weight) const override
Perform the systematic shift.
Definition: MECTuningUtils.h:1003
ana::TrivialPrediction::PredictComponent
Spectrum PredictComponent(osc::IOscCalc *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const override
Definition: MECTuningUtils.h:77
data
const XML_Char const XML_Char * data
Definition: expat.h:268
ana::kIsDytmanMEC
const Cut kIsDytmanMEC
Definition: TruthCuts.cxx:187
ana::kMECDoubleGaussEnhUP
const Var kMECDoubleGaussEnhUP([](const caf::SRProxy *sr){if(sr->mc.nnu< 1) return 1.0;if(!(sr->mc.nu[0].iscc)||sr->mc.nu[0].mode!=caf::kMEC) return 1.0;double q0=kTrueQ0(sr);double q3=kTrueQ3(sr);double Q2=kTrueQ2(sr);return CalcMECDoubleGaussEnhUP(q0, q3, Q2, kGauss2DNorm_1, 0);})
ana::CalcMECGaussEnh
double CalcMECGaussEnh(const double q0, const double q3, const MECGaussEnhParam shift_param, const double shift_sigma)
Definition: MECTuningUtils.h:247
dg::norm_2
double norm_2
Definition: DoubleGaussParams.h:12
dg::mean_q3_2
double mean_q3_2
Definition: DoubleGaussParams.h:14
ana::MECDoubleGaussEnhSystUP::Shift
void Shift(double sigma, caf::SRProxy *sr, double &weight) const override
Perform the systematic shift.
Definition: MECTuningUtils.h:1032
ana::TrivialPrediction::SaveTo
void SaveTo(TDirectory *dir, const std::string &name) const override
Definition: MECTuningUtils.h:94
ana::Current::kCC
Charged-current interactions.
Definition: IPrediction.h:39
Shift
Definition: Shift.h:6
demo5.expt
expt
Definition: demo5.py:34
dg::sigma_q3_shift
double sigma_q3_shift
Definition: DoubleGaussParams.h:27
ana::kTrueQ3
const Var kTrueQ3
Definition: TruthVars.h:38
ana::CalcMinosResSupp
double CalcMinosResSupp(const double Q2, const MinosResSuppParam shift_param, const double shift_sigma)
Definition: MECTuningUtils.h:333
dg::baseline_shift
double baseline_shift
Definition: DoubleGaussParams.h:35
ana::MECDoubleGaussEnhSystDOWN::Shift
void Shift(double sigma, caf::SRProxy *sr, double &weight) const override
Perform the systematic shift.
Definition: MECTuningUtils.h:1059
if
if(dump)
Definition: checkMeasured_Target.C:110
dg::norm_2_shift
double norm_2_shift
Definition: DoubleGaussParams.h:29
dg::sigma_q0
double sigma_q0
Definition: DoubleGaussParams.h:9
dg::corr_shift
double corr_shift
Definition: DoubleGaussParams.h:28
ana::Spectrum::LoadFrom
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:535
ana::CompNormSyst::CompNormSyst
CompNormSyst(const Cut &selCut, double sigmaScale=1.0)
Definition: MECTuningUtils.h:187
ana::kNumu2020NDxsec
const Cut kNumu2020NDxsec
Definition: MECTuningUtils.h:22
ana::CompNormSyst::Shift
void Shift(double sigma, caf::SRProxy *sr, double &weight) const override
Perform the systematic shift.
Definition: MECTuningUtils.h:198
dg::norm_shift
double norm_shift
Definition: DoubleGaussParams.h:23
ana::Q3Q0CutFactory
Cut Q3Q0CutFactory(float loQ3, float hiQ3, float loQ0, float hiQ0)
Definition: MECTuningUtils.h:231
dg::sigma_q3_2_shift
double sigma_q3_2_shift
Definition: DoubleGaussParams.h:33
dg::mean_q0
double mean_q0
Definition: DoubleGaussParams.h:7
ana::MECDoubleGaussEnhSystUP::MECDoubleGaussEnhSystUP
MECDoubleGaussEnhSystUP(const MECDoubleGaussEnhParam &shift_param, const std::string &shift_param_name)
Definition: MECTuningUtils.h:1028
stan::math::exp
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:10
ana::CalcMECDoubleGaussEnh
double CalcMECDoubleGaussEnh(const double q0, const double q3, const double Q2, const MECDoubleGaussEnhParam shift_param, const double shift_sigma, const GaussParams initialParams)
Definition: MECTuningSysts.cxx:25
ana::Cut
_Cut< caf::SRProxy > Cut
Representation of a cut (selection) to be applied to a caf::StandardRecord object.
Definition: Cut.h:96
ana::TrivialPrediction::Predict
Spectrum Predict(osc::IOscCalc *calc) const override
Definition: MECTuningUtils.h:73
sr
caf::StandardRecord * sr
Definition: FillTrainingTrees.h:18
ana::Spectrum::SaveTo
void SaveTo(TDirectory *dir, const std::string &name) const
Definition: Spectrum.cxx:506
caf::Proxy< caf::SRIDBranch >::remid
caf::Proxy< caf::SRRemid > remid
Definition: SRProxy.h:1269
dg::mean_q0_2
double mean_q0_2
Definition: DoubleGaussParams.h:13
ana::NDPredGenerator::Generate
std::unique_ptr< IPrediction > Generate(Loaders &loaders, const SystShifts &shiftMC=kNoShift) const override
Definition: MECTuningUtils.h:137
dg::mean_q3
double mean_q3
Definition: DoubleGaussParams.h:8
wgt
const ana::Var wgt
Definition: specprod_numuccinc.C:37
ana::CompNormSyst::sInstanceCount
static unsigned int sInstanceCount
Definition: MECTuningUtils.h:195
test.z
z
Definition: test.py:28
ana::IRescaledSigmaSyst::SetSigmaScale
void SetSigmaScale(double sc)
Definition: MECTuningUtils.h:179
ana::CalcMECDoubleGaussEnhUP
double CalcMECDoubleGaussEnhUP(const double q0, const double q3, const double Q2, const MECDoubleGaussEnhParam shift_param, const double shift_sigma)
cluncky way of having a separate calculation of up / down shifted versions
Definition: MECTuningUtils.h:498
osc
Oscillation probability calculators.
Definition: Calcs.h:5
sigma
double sigma(TH1F *hist, double percentile)
Definition: median_contours.C:46
ana::DoubleGaussMap
const std::map< std::string, MECDoubleGaussEnhParam > DoubleGaussMap
Definition: MECTuningSysts.h:92
ana::MECGaussEnhSyst::MECGaussEnhSyst
MECGaussEnhSyst(const MECGaussEnhParam &shift_param, const std::string &shift_param_name)
Definition: MECTuningUtils.h:302
ana::kNoShift
const SystShifts kNoShift
Definition: SystShifts.cxx:22
om::cout
OStream cout
Definition: OStream.cxx:6
ana::Flavors::kNuMuToNuMu
(&#39; survival&#39;)
Definition: IPrediction.h:19
ana::NDPredGenerator::NDPredGenerator
NDPredGenerator(const HistAxis axis, const Cut cut, const Var wei=kUnweighted)
Definition: MECTuningUtils.h:125
caf::Proxy< caf::StandardRecord >::mc
caf::Proxy< caf::SRTruthBranch > mc
Definition: SRProxy.h:2138
ana::kMECDoubleGaussEnh
const Var kMECDoubleGaussEnh(const GaussParams initial_params)
Definition: MECTuningSysts.cxx:130
cut
const Cut cut
Definition: exporter_fd.C:30
dg::sigma_q0_2_shift
double sigma_q0_2_shift
Definition: DoubleGaussParams.h:32
genie::units::A
static const double A
Definition: Units.h:82
ana::IRescaledSigmaSyst::IRescaledSigmaSyst
IRescaledSigmaSyst(const std::string &shortName, const std::string &latexName, double sigmaScale=1.0)
Definition: MECTuningUtils.h:176
ana::kWeightPionDeuteriumTune
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;})
norm
Float_t norm
Definition: comparison_ascii.C:71
ana::GetCutIsFitMEC
const Cut GetCutIsFitMEC(const bool isRHC)
Definition: MECTuningUtils.h:211
ana::_HistAxis< Var >
dir
TDirectory * dir
Definition: macro.C:5
dg::norm
double norm
Definition: DoubleGaussParams.h:6
ana::BigChi2SingleSampleExperiment::BigChi2SingleSampleExperiment
BigChi2SingleSampleExperiment(const IPrediction *pred, const Spectrum &data, double scaleFactor=1)
Definition: MECTuningUtils.h:153
loaders
std::vector< Loaders * > loaders
Definition: syst_header.h:386
osc::_IOscCalcAdjustable< double >
ana::assert
assert(nhit_max >=nhit_nbins)
ana::MinosResSuppSyst::Shift
void Shift(double sigma, caf::SRProxy *sr, double &weight) const override
Perform the systematic shift.
Definition: MECTuningUtils.h:1122
galleryMaker.shortName
string shortName
THUMBNAIL BLOCK: We need to make a thumbnail for each.
Definition: galleryMaker.py:469
ana::Sign::kBoth
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
dg::corr_2
double corr_2
Definition: DoubleGaussParams.h:17
ana::IPrediction
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
ana::MECGaussEnhSyst::Shift
void Shift(double sigma, caf::SRProxy *sr, double &weight) const override
Perform the systematic shift.
Definition: MECTuningUtils.h:307
ana::kNumu2020PIDLoosePTPxsec
const Cut kNumu2020PIDLoosePTPxsec([](const caf::SRProxy *sr){return(sr->sel.remid.pid > 0.7 &&sr->sel.cvnloosepreselptp.numuid > 0.82);})
dg::mean_q0_shift
double mean_q0_shift
Definition: DoubleGaussParams.h:24
ana::IPredictionGenerator
Given loaders and an MC shift, Generate() generates an IPrediction.
Definition: PredictionGenerator.h:22
caf::Proxy< caf::StandardRecord >::sel
caf::Proxy< caf::SRIDBranch > sel
Definition: SRProxy.h:2141
caf
This module creates Common Analysis Files.
Definition: FileReducer.h:10
art::to_string
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
ana::CalcMECDoubleGaussEnhDOWN
double CalcMECDoubleGaussEnhDOWN(const double q0, const double q3, const double Q2, const MECDoubleGaussEnhParam shift_param, const double shift_sigma)
Definition: MECTuningUtils.h:612
ana::_Cut< caf::SRProxy >
ana::kNumuQuality
const Cut kNumuQuality
Definition: NumuCuts.h:18
ana::IRescaledSigmaSyst::GetSigmaScale
double GetSigmaScale() const
Definition: MECTuningUtils.h:178
ana::_Var< caf::SRProxy >
ana::MECDoubleGaussEnhSystUP::fShiftParam
const MECDoubleGaussEnhParam fShiftParam
Definition: MECTuningUtils.h:1048
ana::TrivialPrediction::LoadFrom
static std::unique_ptr< TrivialPrediction > LoadFrom(TDirectory *dir, const std::string &name)
Definition: MECTuningUtils.h:109
febshutoff_auto.curr
curr
Definition: febshutoff_auto.py:327
ana::kTrueW
const Var kTrueW
Definition: TruthVars.h:22
ana::kNumuContainND2020
const Cut kNumuContainND2020([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return false;for(unsigned int i=0;i< sr->vtx.elastic.fuzzyk.nshwlid;++i){const caf::SRVector3DProxy &start=sr->vtx.elastic.fuzzyk.png[i].shwlid.start;const caf::SRVector3DProxy &stop=sr->vtx.elastic.fuzzyk.png[i].shwlid.stop;if(std::min(start.X(), stop.X())< -180.0) return false;if(std::max(start.X(), stop.X()) > 180.0) return false;if(std::min(start.Y(), stop.Y())< -180.0) return false;if(std::max(start.Y(), stop.Y()) > 180.0) return false;if(std::min(start.Z(), stop.Z())< 40.0) return false;if(std::max(start.Z(), stop.Z()) > 1525.0) return false;}if(sr->trk.kalman.ntracks< 1) return false;for(unsigned int i=0;i< sr->trk.kalman.ntracks;++i){if(i==sr->trk.kalman.idxremid) continue;else if(sr->trk.kalman.tracks[i].start.Z() > 1275||sr->trk.kalman.tracks[i].stop.Z() > 1275) return false;}return(sr->trk.kalman.ntracks > sr->trk.kalman.idxremid &&sr->slc.firstplane > 1 &&sr->slc.lastplane< 212 &&sr->trk.kalman.tracks[0].start.Z()< 1100 &&(sr->trk.kalman.tracks[0].stop.Z()< 1275 ||sr->sel.contain.kalyposattrans< 55) &&sr->sel.contain.kalfwdcellnd > 5 &&sr->sel.contain.kalbakcellnd > 10);})
Definition: NumuCuts2020.h:31
ana::kUnweighted
const Var kUnweighted
The simplest possible Var, always 1. Used as a default weight.
Definition: Var.h:96
ana::MinosResSuppParam
MinosResSuppParam
Definition: MECTuningUtils.h:331
ana::TrivialPrediction::TrivialPrediction
TrivialPrediction(Loaders &loaders, const HistAxis &axis, const Cut &cut, const SystShifts &shift=kNoShift, const Var &wei=kUnweighted)
Definition: MECTuningUtils.h:49
caf::Proxy< caf::SRCVNResult >::numuid
caf::Proxy< float > numuid
Definition: SRProxy.h:907
dg::sigma_q3_2
double sigma_q3_2
Definition: DoubleGaussParams.h:16
canMan.sign
def sign(x)
Definition: canMan.py:197
ana::SingleSampleExperiment
Compare a single data spectrum to the MC + cosmics expectation.
Definition: SingleSampleExperiment.h:21
ana::GetWeightFromShifts
const Var GetWeightFromShifts(const SystShifts &sys_shifts)
Definition: MECTuningUtils.h:25
string
enum BeamMode string
Generated on Sat Apr 10 2021 14:11:04 for NOvA by   doxygen 1.8.11