mec_tuning.C
Go to the documentation of this file.
1 /*
2  * mec_tuning.C
3  *
4  * Construct numu inclusive plots with Empirical MEC weighted to match Valencia MEC in true q_0
5  *
6  * Created on: Mar 29, 2017
7  * Author: J. Wolcott <jwolcott@fnal.gov>
8  */
9 
10 
11 
12 #include "TCanvas.h"
13 #include "TGraphErrors.h"
14 #include "TH1D.h"
15 #include "TH2D.h"
16 #include "THStack.h"
17 #include "TLegend.h"
18 
19 //#include "CAFAna/Analysis/Ana2017Loaders.h"
21 #include "CAFAna/Core/Loaders.h"
22 #include "CAFAna/Core/Spectrum.h"
27 #include "CAFAna/Core/Var.h"
28 #include "CAFAna/Cuts/SpillCuts.h"
29 #include "CAFAna/Cuts/TruthCuts.h"
36 
37 #include "CAFAna/Core/ISyst.h"
38 #include "CAFAna/Systs/Systs.h"
39 #include "CAFAna/Systs/XSecSysts.h"
40 
41 #include "OscLib/OscCalcPMNSOpt.h"
43 
44 //#include "Weights.h"
45 //#include "DummyExtrapStuff.h"
46 #include "MECTuningUtils.h"
47 
48 #include <fstream>
49 
50 using ana::operator+;
51 
52 namespace jw
53 {
54 
55  const ana::Var GetWeightFromShifts( const ana::SystShifts& sys_shifts )
56  {
57  const ana::Var kWeightFromShifts( [ sys_shifts ]( const caf::SRProxy *sr )
58  {
59  caf::SRProxy* face_palm = const_cast<caf::SRProxy*>( sr );
60  double weight = 1.0;
61  sys_shifts.Shift( face_palm, weight );
62  return weight;
63  });
64  return kWeightFromShifts;
65  }
66 
67  const ana::Cut GetCutIsFitMEC( const bool isRHC )
68  {
69  const ana::Cut kIsFitMEC( [ isRHC ]( const caf::SRProxy * sr )
70  {
71  if ( !ana::kIsDytmanMEC( sr ) || !sr->mc.nu[0].iscc ) return false;
72  if ( !isRHC && sr->mc.nu[0].pdg == 14 ) return true;
73  else if ( isRHC && sr->mc.nu[0].pdg == -14 ) return true;
74  else return false;
75  });
76  return kIsFitMEC;
77  }
78 
79  //===============================================================================
80  // I/O
81  //===============================================================================
82 
84  {
85  return Form( "/nova/ana/users/%s/mec_tuning", getenv( "USER" ) );
86  }
87 
89  {
90  return isRHC ? "numubar" : "numu";
91  }
92 
93  std::string GetFitTag( const bool isRHC, const std::string shift_non_mec )
94  {
95  std::string fit_tag = GetNeutrinoTag( isRHC );
96  if ( !shift_non_mec.empty() ) fit_tag += "_shift_non_mec_" + shift_non_mec;
97  return fit_tag;
98  }
99 
100  std::string GetMECWeightsFileName( const bool isRHC, const std::string shift_non_mec )
101  {
102  return "mec_weights_" + GetFitTag( isRHC, shift_non_mec ) + ".root";
103  }
104 
105  //===============================================================================
106  // Load numu MEC weights for RHC fit
107  //===============================================================================
108 
110  void LoadWeightsTunedNumuMEC( const std::string shift_non_mec )
111  {
112  std::string numu_mec_weights_file_name = GetMECTuningDirectory() + "/" + GetMECWeightsFileName( false, shift_non_mec );
113  if ( gSystem->AccessPathName( numu_mec_weights_file_name.c_str() ) ) throw std::runtime_error( "FHC MEC weights file not found:\n" + numu_mec_weights_file_name );
114  std::cout << "Loading FHC MEC weights from file:\n" << numu_mec_weights_file_name << std::endl;
115  TFile* numu_mec_weights_file = new TFile( numu_mec_weights_file_name.c_str() );
116  hist_input_mec_weights_numu = (TH2D*)numu_mec_weights_file->Get( "numu_mec_weights_smoothed" );
117  }
118 
119  const ana::Var kWeightTunedNumuMEC( []( const caf::SRProxy * sr )
120  {
121  if ( !ana::kIsDytmanMEC( sr ) || !sr->mc.nu[0].iscc ) return 1.0;
122  if ( sr->mc.nu[0].pdg != 14 ) return 1.0;
123  if ( hist_input_mec_weights_numu == NULL ) throw std::runtime_error( "FHC MEC true (q0,q3) weights histogram is NULL" );
124 
125  double q0 = ana::kTrueQ0( sr );
126  double q3 = ana::kTrueQ3( sr );
127 
128  if ( q0 < hist_input_mec_weights_numu->GetYaxis()->GetXmin() ||
129  q0 > hist_input_mec_weights_numu->GetYaxis()->GetXmax() ||
130  q3 < hist_input_mec_weights_numu->GetXaxis()->GetXmin() ||
131  q3 > hist_input_mec_weights_numu->GetXaxis()->GetXmax() ) return 1.0;
132 
133  double weight = hist_input_mec_weights_numu->Interpolate( q3, q0 );
134  //double weight = hist_input_mec_weights_numu->GetBinContent( hist_input_mec_weights_numu->FindFixBin( q3, q0 ) );
135  if ( weight < 0 ) weight = 0.0;
136  return weight;
137  });
138 
139  //===============================================================================
140  // q0 & q3 bins
141  //===============================================================================
142 
143  struct Q3Q0Bin
144  {
145  float loQ3;
146  float hiQ3;
147  float loQ0;
148  float hiQ0;
149  };
150 
151  ana::Cut Q3Q0CutFactory( float loQ3, float hiQ3, float loQ0, float hiQ0 )
152  {
153  return std::move( ana::Cut ( [ loQ3, hiQ3, loQ0, hiQ0 ]( const caf::SRProxy * sr )
154  {
155  double q3 = ana::kTrueQ3( sr );
156  double q0 = ana::kTrueQ0( sr );
157  return q3 > loQ3 && q3 < hiQ3 && q0 > loQ0 && q0 < hiQ0;
158  }
159  ));
160  }
161 
162  /*
163  // we don't want penalties used in fitting
164  class PenaltyFreeSystShifts : public ana::SystShifts
165  {
166  public:
167  virtual ~PenaltyFreeSystShifts() {};
168 
169  //double Penalty() const override;
170  double Penalty() const;
171  };
172  double PenaltyFreeSystShifts::Penalty() const
173  {
174  std::cout << "Penalty is 0!" << std::endl;
175  return 0;
176  }
177  */
178 
180  {
181  public:
182  IRescaledSigmaSyst(const std::string & shortName, const std::string & latexName, double sigmaScale=1.0)
183  : ana::ISyst(shortName, latexName), fSigmaScale(sigmaScale) {};
184  double GetSigmaScale() const { return fSigmaScale; };
185  void SetSigmaScale(double sc) { fSigmaScale = sc; };
186  protected:
187  double fSigmaScale;
188  };
189 
191  {
192  public:
193  CompNormSyst(const ana::Cut & selCut, double sigmaScale=1.0)
194  : fID(sInstanceCount++),
195  IRescaledSigmaSyst("CompNormShift_" + std::to_string(sInstanceCount), "Component Normalization Shift " + std::to_string(sInstanceCount), sigmaScale),
196  fSelCut(selCut)
197  {}
198  void Shift(double sigma, caf::SRProxy* sr, double& weight) const override;
199 
200  private:
202  unsigned int fID;
203 
204  static unsigned int sInstanceCount;
205  };
206  unsigned int CompNormSyst::sInstanceCount = 0;
207  void CompNormSyst::Shift(double sigma, caf::SRProxy* sr, double& weight) const
208  {
209  if (!fSelCut(sr))
210  return;
211 
212  double wgt = 1 + fSigmaScale * sigma;
213  if (wgt < 0)
214  wgt = 0;
215  weight *= wgt;
216  }
217 
219  {
220  public:
221  MECInitStateNPFracShift( double sigmaScale = 1.0 )
222  : IRescaledSigmaSyst( "MECInitStateCompShift", "MEC initial state np fraction", sigmaScale )
223  {};
224 
225  // this guy is a bit hacksy since it doesn't make sense
226  // for the fraction of np pairs to be negative or greater than 1
227  void Shift( double sigma, caf::SRProxy* sr, double& weight ) const override
228  {
229  if ( !ana::kIsDytmanMEC( sr ) ) return;
230 
231  // Correct nominal np fraction for GENIE Empirical MEC
232  const double nominalNPfrac = 0.89;
233  // 10% nominal uncertainty
234  double newNPfrac = nominalNPfrac + 0.1 * fSigmaScale * sigma;
235  // pin to being between 0 and 1
236  if ( newNPfrac > 1 ) newNPfrac = 1;
237  else if ( newNPfrac < 0 ) newNPfrac = 0;
238 
239  if ( ana::kHitNuc( sr ) == 2000000201 ) // np
240  weight *= newNPfrac / nominalNPfrac;
241  else
242  weight *= ( 1 - newNPfrac ) / ( 1 - nominalNPfrac );
243  } // MECInitStateNPFracShift::Shift()
244  }; // class MECInitStateNPFracShift
245 
246 } // namespace jw
247 
248 std::vector<std::string> GetFiles( const std::string playlist, const int n_files_max = -1 )
249 {
250  std::ifstream input( playlist.c_str() );
251  std::vector<std::string> files;
252  std::string s;
253  int count = 0;
254  while( getline( input, s ) )
255  {
256  count++;
257  if ( n_files_max > -1 && count > n_files_max ) break;
258  files.push_back( s );
259  }
260  return files;
261 }
262 
263 ////////////////////////
264 
265 void mec_tuning
266 (
267  const std::string beam,
268  const std::string shift_non_mec = ""
269 )
270 {
271 
272  if ( beam != "fhc" && beam != "rhc" ) throw std::runtime_error( "Beam must be \"fhc\" or \"rhc\"" );
273  const bool isRHC = beam == "rhc";
274 
275  if ( isRHC ) jw::LoadWeightsTunedNumuMEC( shift_non_mec );
276 
277  ana::SystShifts kShiftsNonMEC;
278  if ( shift_non_mec.empty() )
279  {
280  kShiftsNonMEC = ana::kNoShift;
281  }
282  else if ( shift_non_mec == "down" )
283  {
284  kShiftsNonMEC.SetShift( &ana::kMAQEGenieReducedSyst2018, +1 );
285  kShiftsNonMEC.SetShift( &ana::kRPACCQESuppSyst2018, +1 );
286  kShiftsNonMEC.SetShift( &ana::kRPACCQEEnhSyst2018, +1 );
287  kShiftsNonMEC.SetShift( ana::GetGenieKnobSyst( rwgt::fReweightCCQEPauliSupViaKF ), -1 );
288  kShiftsNonMEC.SetShift( ana::GetGenieKnobSyst( rwgt::fReweightMaCCRES ), -1 );
289  kShiftsNonMEC.SetShift( ana::GetGenieKnobSyst( rwgt::fReweightMvCCRES ), -1 );
290  kShiftsNonMEC.SetShift( &ana::kRPARESSyst2018, -1 );
291  //kShiftsNonMEC.SetShift( &ana::kDirectRelHadEScaleSyst2017, -1 );
292  }
293  else if ( shift_non_mec == "up" )
294  {
295  kShiftsNonMEC.SetShift( &ana::kMAQEGenieReducedSyst2018, -1 );
296  kShiftsNonMEC.SetShift( &ana::kRPACCQESuppSyst2018, -1 );
297  kShiftsNonMEC.SetShift( &ana::kRPACCQEEnhSyst2018, -1 );
298  kShiftsNonMEC.SetShift( ana::GetGenieKnobSyst( rwgt::fReweightCCQEPauliSupViaKF ), +1 );
299  kShiftsNonMEC.SetShift( ana::GetGenieKnobSyst( rwgt::fReweightMaCCRES ), +1 );
300  kShiftsNonMEC.SetShift( ana::GetGenieKnobSyst( rwgt::fReweightMvCCRES ), +1 );
301  kShiftsNonMEC.SetShift( &ana::kRPARESSyst2018, +1 );
302  //kShiftsNonMEC.SetShift( &ana::kDirectRelHadEScaleSyst2017, +1 );
303  }
304  else throw std::runtime_error( Form( "Unrecognized shift: \"%s\"", shift_non_mec.c_str() ) );
305 
307  const ana::Var wgt_shift_non_mec = !shift_non_mec.empty() ? jw::GetWeightFromShifts( kShiftsNonMEC ) : ana::kUnweighted;
308  const ana::Var wgt_mec_numu = isRHC ? jw::kWeightTunedNumuMEC : ana::kUnweighted;
309  const ana::Var wgt_fit = ana::kPPFXFluxCVWgt * wgt_non_mec * wgt_shift_non_mec * wgt_mec_numu;
310 
311  const ana::Cut kSelCuts = ana::kNumuCutND2018;
312 
313  const int nbinsQ3 = 20;
314  const float fitMinQ3 = 0.0;
315  const float fitMaxQ3 = 1.0;
316  const float binWidthQ3 = ( fitMaxQ3 - fitMinQ3 ) / nbinsQ3;
317 
318  const int nbinsQ0 = 16;
319  const float fitMinQ0 = 0.0;
320  const float fitMaxQ0 = 0.8;
321  const float binWidthQ0 = ( fitMaxQ0 - fitMinQ0 ) / nbinsQ0;
322 
323  std::vector<jw::Q3Q0Bin> binsQ3Q0;
324  //std::vector<std::unique_ptr<const ana::ISyst>> systsQ3Q0;
325  std::vector<std::unique_ptr<const jw::CompNormSyst>> systsQ3Q0;
326 
327  for ( int i = 0; i < nbinsQ3; ++i )
328  {
329  for ( int j = 0; j < nbinsQ0; ++j )
330  {
331  double q3 = fitMinQ3 + binWidthQ3 * ( i + 0.5 );
332  double q0 = fitMinQ0 + binWidthQ0 * ( j + 0.5 );
333  if ( q3 * q3 < q0 * q0 ) continue;
334  jw::Q3Q0Bin b;
335  b.loQ3 = q3 - binWidthQ3 / 2;
336  b.hiQ3 = q3 + binWidthQ3 / 2;
337  b.loQ0 = q0 - binWidthQ0 / 2;
338  b.hiQ0 = q0 + binWidthQ0 / 2;
339  binsQ3Q0.push_back( b );
340 
341  double sig_scale = 1.0;
342  //double sig_scale = ( q3 < 0.5 && q0 < 0.2 ) ? 2.0 : 1.0;
343 
344  systsQ3Q0.emplace_back( std::make_unique<jw::CompNormSyst>( jw::Q3Q0CutFactory( b.loQ3, b.hiQ3, b.loQ0, b.hiQ0 ) && jw::GetCutIsFitMEC( isRHC ) && kSelCuts, sig_scale ) );
345  }
346  }
347 
348  std::cout << "systsQ3Q0.size() = " << systsQ3Q0.size() << std::endl;
349 
350  std::vector<const ana::ISyst*> systs_ptrs;
351  for ( const auto & s : systsQ3Q0 )
352  {
353  systs_ptrs.push_back( s.get() );
354  }
355 
356  std::string defn_data = "";
357  std::string defn_mc = "";
358  if ( beam == "fhc" )
359  {
360  defn_data = "defname: prod_sumdecaf_R17-09-05-prod4recopreview.f_nd_numi_fhc_period2_v1_addShortSimpleCVN_goodruns_numu2018 with limit 5";
361  defn_mc = "defname: prod_sumdecaf_R17-11-14-prod4reco.d_nd_genie_nonswap_fhc_nova_v08_period2_v1_numu2018 with limit 10";
362  }
363  else if ( beam == "rhc" )
364  {
365  defn_data = "prod_sumdecaf_R17-09-05-prod4recopreview.f_nd_numi_rhc_epoch6a_v1_addShortSimpleCVN_goodruns_numu2018 with limit 5";
366  defn_mc = "prod_sumdecaf_R17-11-14-prod4reco.e_nd_genie_nonswap_rhc_nova_v08_epoch6a_v1_numu2018 with limit 10";
367  }
368  if ( defn_data.empty() ) throw std::runtime_error( "Data SAM Definition is empty" );
369  if ( defn_mc .empty() ) throw std::runtime_error( "MC SAM Definition is empty" );
370 
372  loaders.SetLoaderPath( defn_data, caf::kNEARDET, ana::Loaders::kData );
373  loaders.SetLoaderPath( defn_mc , caf::kNEARDET, ana::Loaders::kMC );
375 
376  ana::HistAxis axis_fit( "Reco |#vec{q}| (GeV)" , ana::Binning::Simple( 20, 0.0, 1.0 ), ana::kRecoQmag,
377  "Reco E_{had, vis} (GeV)", ana::Binning::Simple( 20, 0.0, 0.4 ), ana::kNumuHadVisE );
378 
379  ana::Spectrum spec_fit_data
380  (
382  axis_fit,
383  kSelCuts
384  );
385 
386  ana::Spectrum spec_fit_mc_raw
387  (
389  axis_fit,
390  kSelCuts,
392  wgt_fit
393  );
394 
395  ana::Spectrum spec_fit_non_mec_raw
396  (
398  axis_fit,
399  kSelCuts && !ana::kIsDytmanMEC,
400  ana::kNoShift,
401  wgt_fit
402  );
403 
404  osc::NoOscillations no_osc_calc;
405  jw::NDPredGenerator pred_gen( axis_fit, kSelCuts, wgt_fit );
406  ana::PredictionInterp pred_interp( systs_ptrs, &no_osc_calc, pred_gen, loaders );
407 
408  loaders.Go();
409 
410  ana::SingleSampleExperiment expt( &pred_interp, spec_fit_data );
412  ana::MinuitFitter fitter( &expt, {}, systs_ptrs );
413 
414  ana::SystShifts shifts;
415  //jw::PenaltyFreeSystShifts shifts;
416  for ( auto & syst : systs_ptrs )
417  {
418  shifts.SetShift( syst, 0 );
419  }
420 
421  fitter.Fit( &calc, shifts, ana::IFitter::kVerbose );
422  // the results of the fitting operation are now in the shifts object.
423 
424  //==================================================
425  // Output files
426  //==================================================
427 
428  const std::string out_dir = jw::GetMECTuningDirectory();
429  if ( gSystem->AccessPathName( out_dir.c_str() ) ) gSystem->mkdir( out_dir.c_str(), true );
430 
431  const std::string fit_tag = jw::GetFitTag( isRHC, shift_non_mec );
432 
433  std::string fit_results_file_name = out_dir + "/mec_fit_results_" + fit_tag + ".txt";
434  std::ofstream fit_results_file( fit_results_file_name.c_str(), std::ofstream::out );
435 
436  for ( auto & syst : systs_ptrs )
437  {
438  fit_results_file << syst->ShortName() << " fitted shift = " << shifts.GetShift( syst ) << std::endl;
439  }
440 
441  TH2D* hist_tuned_mec_weights = new TH2D( "tuned_mec_weights", ";True |#vec{q}| (GeV);True q_{0} (GeV);Weight",
442  nbinsQ3, fitMinQ3, fitMaxQ3, nbinsQ0, fitMinQ0, fitMaxQ0 );
443 
444  for ( std::size_t compIdx = 0; compIdx < binsQ3Q0.size(); compIdx++ )
445  {
446  auto b = binsQ3Q0[ compIdx ];
447  float q3 = ( b.loQ3 + b.hiQ3 ) / 2;
448  float q0 = ( b.loQ0 + b.hiQ0 ) / 2;
449  const auto & syst = systsQ3Q0[ compIdx ];
450  double sigma_scale = syst.get()->GetSigmaScale();
451  double scale_factor = 1 + sigma_scale * shifts.GetShift( syst.get() );
452  hist_tuned_mec_weights->SetBinContent( hist_tuned_mec_weights->FindFixBin( q3, q0 ), scale_factor );
453  std::cout << Form( "q3 = %.2f, q0 = %.2f, sigma_scale = %.1f, scale_factor = %.2f", q3, q0, sigma_scale, scale_factor ) << std::endl;
454  }
455 
456  TH2D* hist_tuned_mec_weights_smoothed = (TH2D*)hist_tuned_mec_weights->Clone( "tuned_mec_weights_smoothed" );
457  hist_tuned_mec_weights_smoothed->Smooth( 1, "k3a" );
458 
459  std::string weights_file_name = out_dir + "/" + jw::GetMECWeightsFileName( isRHC, shift_non_mec );
460  TFile weights_file( weights_file_name.c_str(), "recreate" );
461 
462  const std::string neutrino_tag = jw::GetNeutrinoTag( isRHC );
463  hist_tuned_mec_weights ->Write( Form( "%s_mec_weights", neutrino_tag.c_str() ) );
464  hist_tuned_mec_weights_smoothed->Write( Form( "%s_mec_weights_smoothed", neutrino_tag.c_str() ) );
465  if ( isRHC && jw::hist_input_mec_weights_numu != NULL ) jw::hist_input_mec_weights_numu->Write();
466  weights_file.Close();
467 
468  double pot = spec_fit_data.POT();
469  auto h_data = spec_fit_data.ToTH2( pot );
470  auto h_mc_raw = spec_fit_mc_raw.ToTH2( pot );
471  auto h_non_mec_raw = spec_fit_non_mec_raw.ToTH2( pot );
472  ana::Spectrum spec_fit_mc_tuned( pred_interp.PredictSyst( &calc, shifts ) );
473  auto h_mc_tuned = spec_fit_mc_tuned.ToTH2( pot );
474 
475  std::string fit_hists_file_name = out_dir + "/mec_fit_hists_" + fit_tag + ".root";
476  TFile fit_hists_file( fit_hists_file_name.c_str(), "recreate" );
477  fit_hists_file.cd();
478 
479  spec_fit_data.SaveTo(& fit_hists_file, "spec_data" ) ;
480  h_data->Write( "data" );
481  h_mc_raw->Write( "mc_raw" );
482  h_non_mec_raw->Write( "non_mec_raw" );
483  h_mc_tuned->Write( "mc_tuned" );
484 
485  fit_hists_file.Close();
486 
487 }
Near Detector underground.
Definition: SREnums.h:10
Pass neutrinos through unchanged.
Definition: IOscCalc.h:45
std::string GetFitTag(const bool isRHC, const std::string shift_non_mec)
Definition: mec_tuning.C:93
Implements systematic errors by interpolation between shifted templates.
TH2 * ToTH2(double exposure, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Spectrum must be 2D to obtain TH2.
Definition: Spectrum.cxx:165
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
const NOvARwgtSyst kRPACCQEEnhSyst2018("RPAShapeenh2018","RPA shape: higher-Q^{2} enhancement (2018)", novarwgt::kRPACCQEEnhSyst2018)
Definition: RPASysts.h:13
correl_xv GetYaxis() -> SetDecimals()
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
const Var weight
CompNormSyst(const ana::Cut &selCut, double sigmaScale=1.0)
Definition: mec_tuning.C:193
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
const Var kRecoQmag([](const caf::SRProxy *sr){return sqrt(kRecoQ2(sr)+util::sqr(kHadE(sr)));})
Reconstructed three-momentum transfer.
Definition: NumuVars.h:149
const Cut kNumuCutND2018
Definition: NumuCuts2018.h:41
const NOvARwgtSyst kMAQEGenieReducedSyst2018(genie::rew::GSyst::AsString(genie::rew::EGSyst(rwgt::fReweightMaCCQE))++"_reduced_2018", genie::rew::GSyst::AsString(genie::rew::EGSyst(rwgt::fReweightMaCCQE))++"_reduced_2018", novarwgt::kMAQEGenieReducedSyst2018)
2018 &#39;reduced&#39; M_A^QE shift. See documentation in NOvARwgt
Definition: XSecSysts.h:51
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2125
Collection of SpectrumLoaders for many configurations.
Definition: Loaders.h:23
caf::Proxy< std::vector< caf::SRNeutrino > > nu
Definition: SRProxy.h:617
Spectrum PredictSyst(osc::IOscCalc *calc, const SystShifts &shift) const override
const Var kRPAWeightRES2017
Definition: GenieWeights.h:115
MECInitStateNPFracShift(double sigmaScale=1.0)
Definition: mec_tuning.C:221
const Var kTrueQ0
Definition: TruthVars.h:32
unsigned int fID
Definition: mec_tuning.C:202
osc::OscCalcDumb calc
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
const Var kNumuHadVisE([](const caf::SRProxy *sr){return kNumuHadCalE(sr)+kNumuHadTrkE(sr);})
Definition: NumuVars.h:124
float hiQ3
Definition: mec_tuning.C:146
const Var kHitNuc
Definition: TruthVars.h:19
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
const ana::Var kWeightTunedNumuMEC([](const caf::SRProxy *sr){if(!ana::kIsDytmanMEC(sr)||!sr->mc.nu[0].iscc) return 1.0;if(sr->mc.nu[0].pdg!=14) return 1.0;if(hist_input_mec_weights_numu==NULL) throw std::runtime_error("FHC MEC true (q0,q3) weights histogram is NULL");double q0=ana::kTrueQ0(sr);double q3=ana::kTrueQ3(sr);if(q0< hist_input_mec_weights_numu->GetYaxis() ->GetXmin()|| q0 > hist_input_mec_weights_numu->GetYaxis() ->GetXmax()|| q3< hist_input_mec_weights_numu->GetXaxis() ->GetXmin()|| q3 > hist_input_mec_weights_numu->GetXaxis() ->GetXmax()) return 1.0;double weight=hist_input_mec_weights_numu->Interpolate(q3, q0);if(weight< 0) weight=0.0;return weight;})
void Go()
Call Go() on all the loaders.
Definition: Loaders.cxx:162
float hiQ0
Definition: mec_tuning.C:148
const Var kRPAWeightCCQE2018
Definition: GenieWeights.h:96
const Cut kIsDytmanMEC
Definition: TruthCuts.cxx:187
std::string GetMECWeightsFileName(const bool isRHC, const std::string shift_non_mec)
Definition: mec_tuning.C:100
void SetSigmaScale(double sc)
Definition: mec_tuning.C:185
const XML_Char * s
Definition: expat.h:262
void SetSpillCut(const SpillCut &cut)
Definition: Loaders.cxx:121
const ana::Cut fSelCut
Definition: mec_tuning.C:201
const Var kRescaleMAQE
Definition: GenieWeights.h:61
std::vector< std::string > GetFiles(const std::string playlist, const int n_files_max=-1)
Definition: mec_tuning.C:248
Definition: Shift.h:6
expt
Definition: demo5.py:34
const Var kTrueQ3
Definition: TruthVars.h:38
if(dump)
std::vector< std::unique_ptr< const CompNormSyst > > systsQ3Q0(int nbinsQ3=20, float fitMinQ3=0.0, float fitMaxQ3=1.0, int nbinsQ0=16, float fitMinQ0=0.0, float fitMaxQ0=0.8, Cut cut=GetCutIsFitMEC(false))
SpectrumLoaderBase & GetLoader(caf::Det_t det, DataMC datamc, DataSource src=kBeam, SwappingConfig swap=kNonSwap)
Retrieve a specific loader.
Definition: Loaders.cxx:129
std::string getenv(std::string const &name)
float loQ0
Definition: mec_tuning.C:147
Optimized version of OscCalcPMNS.
Definition: StanTypedefs.h:31
#define pot
caf::StandardRecord * sr
const double j
Definition: BetheBloch.cxx:29
const ana::Var wgt
double sigma(TH1F *hist, double percentile)
const SystShifts kNoShift
Definition: SystShifts.cxx:21
OStream cout
Definition: OStream.cxx:6
caf::Proxy< caf::SRTruthBranch > mc
Definition: SRProxy.h:2137
const ana::Var GetWeightFromShifts(const ana::SystShifts &sys_shifts)
Definition: mec_tuning.C:55
void Shift(double sigma, caf::SRProxy *sr, double &weight) const override
Perform the systematic shift.
Definition: mec_tuning.C:227
void Shift(double sigma, caf::SRProxy *sr, double &weight) const override
Perform the systematic shift.
Definition: mec_tuning.C:207
TH2D * hist_input_mec_weights_numu
Definition: mec_tuning.C:109
std::vector< Loaders * > loaders
Definition: syst_header.h:386
std::string GetMECTuningDirectory()
Definition: mec_tuning.C:83
const hit & b
Definition: hits.cxx:21
string shortName
THUMBNAIL BLOCK: We need to make a thumbnail for each.
float loQ3
Definition: mec_tuning.C:145
double GetSigmaScale() const
Definition: mec_tuning.C:184
const NOvARwgtSyst kRPACCQESuppSyst2018("RPAShapesupp2018","RPA shape: low-Q^{2} suppression (2018)", novarwgt::kRPACCQESuppSyst2018)
Definition: RPASysts.h:14
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
void mec_tuning(const std::string beam, const std::string shift_non_mec="")
Definition: mec_tuning.C:266
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
static unsigned int sInstanceCount
Definition: mec_tuning.C:204
const NOvARwgtSyst * GetGenieKnobSyst(rwgt::ReweightLabel_t knobIdx, std::string altName, std::string altLabel)
Convenience function to get a GENIE knob syst. (Allows using the GENIE knob name & description as the...
Definition: XSecSysts.cxx:119
void SetLoaderPath(const std::string &path, caf::Det_t det, DataMC datamc, DataSource src=kBeam, SwappingConfig swap=kNonSwap)
Configure loader via wildcard path.
Definition: Loaders.cxx:25
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:107
const Var kRescaleHighWDIS
Definition: GenieWeights.h:316
const Var kUnweighted
The simplest possible Var, always 1. Used as a default weight.
Definition: Var.h:96
void SetShift(const ISyst *syst, double shift, bool force=false)
Definition: SystShifts.cxx:80
const ana::Cut GetCutIsFitMEC(const bool isRHC)
Definition: mec_tuning.C:67
const Var kFixNonres1Pi
Definition: GenieWeights.h:127
IRescaledSigmaSyst(const std::string &shortName, const std::string &latexName, double sigmaScale=1.0)
Definition: mec_tuning.C:182
ana::Cut Q3Q0CutFactory(float loQ3, float hiQ3, float loQ0, float hiQ0)
Definition: mec_tuning.C:151
void LoadWeightsTunedNumuMEC(const std::string shift_non_mec)
Definition: mec_tuning.C:110
Compare a single data spectrum to the MC + cosmics expectation.
std::string GetNeutrinoTag(const bool isRHC)
Definition: mec_tuning.C:88
void Shift(caf::SRProxy *sr, double &weight) const
Definition: SystShifts.cxx:164
float h_data[xbins][ybins]
const NOvARwgtSyst kRPARESSyst2018("RPAShapeRES2018","RPA shape: resonance events (2018)", novarwgt::kRPARESSyst2018)
Definition: RPASysts.h:21
Perform MINUIT fits in one or two dimensions.
Definition: MinuitFitter.h:17
enum BeamMode string