make_xsec_wgts_2018_hists.C
Go to the documentation of this file.
1 #ifdef __CINT__
3 (
4  const std::string beam,
5  const std::string type,
6  const std::string period = "full"
7  //const std::string shift_non_mec = "",
8 )
9 {
10  std::cout << "Sorry, you must run in compiled mode" << std::endl;
11 }
12 #else
13 
14 //#include "StandardRecord/Proxy/SRProxy.h"
15 
16 #include "CAFAna/Core/Binning.h"
17 #include "CAFAna/Core/Spectrum.h"
19 #include "CAFAna/Core/Utilities.h"
20 #include "3FlavorAna/Core/HistAxes.h"
21 #include "CAFAna/Core/Loaders.h"
22 #include "CAFAna/Cuts/Cuts.h"
23 #include "CAFAna/Cuts/SpillCuts.h"
28 #include "CAFAna/Vars/Vars.h"
30 #include "CAFAna/Vars/TruthVars.h"
33 #include "CAFAna/Vars/HistAxes.h"
34 //#include "CAFAna/Vars/XsecTunes.h"
36 #include "CAFAna/Systs/XSecSysts.h"
37 //#include "MCReweight/func/HistWrapper.h"
38 
39 #include "TTree.h"
40 #include "TH2.h"
41 
42 #include <fstream>
43 
44 using namespace ana;
45 using namespace caf;
46 
47 
48 TFile* file_mec_wgt_numu = NULL;
49 TFile* file_mec_wgt_numubar = NULL;
50 TH2D* hist_mec_wgt_numu = NULL;
51 TH2D* hist_mec_wgt_numubar = NULL;
52 
53 void LoadTunedMECWeights( const std::string shift_non_mec = "", const bool extra_stat_mc = false, const bool no_smoothing = false, const bool mec_refit = false )
54 {
55  // "up" and "down" refer to shifting the total MC in the (QEL,RES) region up and down in (q0,q3)
56  if ( !shift_non_mec.empty() && shift_non_mec != "Down" && shift_non_mec != "Up" ) throw std::runtime_error( Form( "Unrecognized shift: \"%s\"", shift_non_mec.c_str() ) );
57 
58  const bool finer_bins = true;
59  std::string fit_tag = "";
60  if ( !shift_non_mec.empty() ) fit_tag += "_ShiftNonMEC" + shift_non_mec;
61  if ( finer_bins ) fit_tag += "_FinerBins";
62  if ( extra_stat_mc ) fit_tag += "_ExtraStatMC";
63 
64  std::string base_indir = mec_refit ? "/nova/app/users/mislivec/mec_weights" : "/nova/app/users/mislivec/MEC_rwgt/mec_fit_q0q3/test";
65  std::string weights_file_name_numu = Form( "%s/FHC_MECFit_2018%s/tuned_mec_weights.root", base_indir.c_str(), fit_tag.c_str() );
66  std::string weights_file_name_numubar = Form( "%s/RHC_MECFit_2018%s/tuned_mec_weights.root", base_indir.c_str(), fit_tag.c_str() );
67 
68  if ( gSystem->AccessPathName( weights_file_name_numu.c_str() ) ) throw std::runtime_error( "Numu MEC weights file not found:\n" + weights_file_name_numu );
69  if ( gSystem->AccessPathName( weights_file_name_numubar.c_str() ) ) throw std::runtime_error( "Numu MEC weights file not found:\n" + weights_file_name_numubar );
70 
71  file_mec_wgt_numu = new TFile( weights_file_name_numu.c_str() );
72  file_mec_wgt_numubar = new TFile( weights_file_name_numubar.c_str() );
73 
74  std::cout << "MEC weight files:" << std::endl;
75  std::cout << file_mec_wgt_numu->GetName() << std::endl;
76  std::cout << file_mec_wgt_numubar->GetName() << std::endl;
77 
78  std::string hist_name_numu = no_smoothing ? "numu_mec_weights" : "numu_mec_weights_smoothed";
79  std::string hist_name_numubar = no_smoothing ? "numubar_mec_weights" : "numubar_mec_weights_smoothed";
80 
81  hist_mec_wgt_numu = (TH2D*)file_mec_wgt_numu ->Get( hist_name_numu.c_str() );
82  hist_mec_wgt_numubar = (TH2D*)file_mec_wgt_numubar->Get( hist_name_numubar.c_str() );
83 
84  //hist_mec_wgt_numu = (TH2D*)file_mec_wgt_numu ->Get( "numu_mec_weights_smoothed" );
85  //hist_mec_wgt_numubar = (TH2D*)file_mec_wgt_numubar->Get( "numubar_mec_weights_smoothed" );
86 
87 }
88 
89 const ana::Var GetTunedMECWeight( const bool interpolate = false )
90 {
91  const Var kWeightTunedMEC( [ interpolate ]( const caf::SRProxy * sr )
92  {
93  if ( !ana::kIsDytmanMEC( sr ) ) return 1.0;
94 
95  double q0 = kTrueQ0( sr );
96  double q3 = kTrueQ3( sr );
97 
98  double weight = 1.0;
99 
100  if ( sr->mc.nu[0].pdg > 0 )
101  {
102  if ( hist_mec_wgt_numu == NULL ) throw std::runtime_error( "Numu MEC weights histogram is NULL" );
103 
104  if ( q0 < hist_mec_wgt_numu->GetYaxis()->GetXmin() ||
105  q0 > hist_mec_wgt_numu->GetYaxis()->GetXmax() ||
106  q3 < hist_mec_wgt_numu->GetXaxis()->GetXmin() ||
107  q3 > hist_mec_wgt_numu->GetXaxis()->GetXmax() ) return 1.0;
108 
109  weight = interpolate ?
110  hist_mec_wgt_numu->Interpolate( q3, q0 ) :
111  hist_mec_wgt_numu->GetBinContent( hist_mec_wgt_numu->FindFixBin( q3, q0 ) );
112  }
113  else
114  {
115  if ( hist_mec_wgt_numubar == NULL ) throw std::runtime_error( "Numubar MEC weights histogram is NULL" );
116 
117  if ( q0 < hist_mec_wgt_numubar->GetYaxis()->GetXmin() ||
118  q0 > hist_mec_wgt_numubar->GetYaxis()->GetXmax() ||
119  q3 < hist_mec_wgt_numubar->GetXaxis()->GetXmin() ||
120  q3 > hist_mec_wgt_numubar->GetXaxis()->GetXmax() ) return 1.0;
121 
122  weight = interpolate ?
123  hist_mec_wgt_numubar->Interpolate( q3, q0 ) :
124  hist_mec_wgt_numubar->GetBinContent( hist_mec_wgt_numubar->FindFixBin( q3, q0 ) );
125  }
126  if ( weight < 0 ) weight = 0.0;
127  return weight;
128  });
129  return kWeightTunedMEC;
130 }
131 
132 const Var kMinosRESSupp( []( const caf::SRProxy * sr )
133 {
134  if ( !kIsRes( sr ) || !kIsNumuCC( sr ) ) return 1.0;
135  if ( sr->mc.nu[0].tgtA == 1 ) return 1.0; // Exclude hydrogen
136  double Q2 = kTrueQ2( sr );
137  return 1.010 / ( 1 + exp( 1 - sqrt( Q2 ) / 0.156 ) ); // From PhysRevD.91.012005
138 });
139 
140 const Var kMINERvALowQ2RESSupp( []( const caf::SRProxy * sr )
141 {
142  if ( !kIsRes( sr ) ) return 1.0;
143  if ( sr->mc.nu[0].tgtA == 1 ) return 1.0; // Exclude hydrogen
144 
145  // arXiv:1903.01558
146  // Joint FrAbs Tune
147  double R1 = 0.32;
148  double R2 = 0.50;
149 
150  double x1 = 0.0;
151  double x2 = 0.35;
152  double x3 = 0.7;
153 
154  double Q2 = kTrueQ2( sr );
155  if ( Q2 > x3 ) return 1.0;
156 
157  double R = R2 * ( Q2 - x1 ) * ( Q2 - x3 ) / ( ( x2 - x1 ) * ( x2 - x3 ) ) + ( Q2 - x1 ) * ( Q2 - x2 ) / ( ( x3 - x1 ) * ( x3 - x2 ) );
158  double weight = 1 - ( 1 - R1 ) * pow( 1 - R, 2 );
159  if ( weight < 0 )
160  {
161  std::cout << "kMINERvALowQ2RESSupp = " << weight << std::endl;
162  std::cout << "Q2 = " << Q2 << std::endl;
163  weight = 0;
164  }
165  return weight;
166 });
167 
168 /*
169 const Var kHitNucBin( []( const caf::SRProxy * sr )
170 {
171  if ( !kIsDytmanMEC( sr ) ) return -1.0;
172  else if ( kHitNuc( sr ) == 2000000200 ) return 0.0;
173  else if ( kHitNuc( sr ) == 2000000201 ) return 1.0;
174  else if ( kHitNuc( sr ) == 2000000202 ) return 2.0;
175  else return 3.0;
176 });
177 */
178 
180 {
181  if ( !kIsDytmanMEC( sr ) ) return 1.0;
182 
183  double weight = kEmpiricalMECtoValenciaMECWgt( sr );
184 
185  const double np_frac_empirical = 0.89;
186  const double np_frac_valencia = 2.0 / 3.0;
187 
188  if ( ana::kHitNuc( sr ) == 2000000201 ) // np
189  weight *= np_frac_valencia / np_frac_empirical;
190  else
191  weight *= ( 1 - np_frac_valencia ) / ( 1 - np_frac_empirical );
192 
193  return weight;
194 });
195 
196 const Cut kCutDISAnomaly( []( const caf::SRProxy* sr )
197 {
198  return sr->mc.nnu == 0;
199 });
200 
201 const Cut kCutOnHydrogen( []( const caf::SRProxy* sr )
202 {
203  return sr->mc.nu[0].tgtA == 1;
204 });
205 
206 /*
207 const Var kTrueY( []( const caf::SRProxy* sr )
208 {
209  return ( sr->mc.nnu == 0 ) ? -1.0 : float( sr->mc.nu[0].y );
210 });
211 */
212 
214 
217 
218 const Cut kWSelection( []( const caf::SRProxy* sr )
219 {
220  double W = kRecoW( sr );
221  return W > 1.2 && W < 1.5;
222 });
223 
224 const Cut kIsMuonNeutrino( []( const caf::SRProxy* sr )
225 {
226  return sr->hdr.ismc && sr->mc.nnu == 1 && sr->mc.nu[0].pdg == 14;
227 });
228 
229 const Cut kIsMuonAntiNeutrino( []( const caf::SRProxy* sr )
230 {
231  return sr->hdr.ismc && sr->mc.nnu == 1 && sr->mc.nu[0].pdg == -14;
232 });
233 
234 /*
235 SystShifts kShiftsMECShapeUp;
236 kShiftsMECShapeUp.SetShift( &kMECShapeSyst2018RPAFixNu, +1 );
237 kShiftsMECShapeUp.SetShift( &kMECShapeSyst2018RPAFixAntiNu, +1 );
238 
239 SystShifts kShiftsMECShapeDown;
240 kShiftsMECShapeDown.SetShift( &kMECShapeSyst2018RPAFixNu, -1 );
241 kShiftsMECShapeDown.SetShift( &kMECShapeSyst2018RPAFixAntiNu, -1 );
242 */
243 
244 //const Var kWeightTunedMECInterp = GetTunedMECWeight( true );
245 
246 const NuTruthVar kTrueW_Eqn_NT( []( const caf::SRNeutrinoProxy * nu )
247 {
248  const double M_p = 0.938272;
249  double W2 = M_p * M_p + 2 * M_p * ( nu->y * nu->E ) - nu->q2;
250  if ( W2 > 0 ) return sqrt( W2 );
251  else return -5.;
252 });
254 
255 const ana::Var GetWResolution( const bool fractional = false, const bool true_w_eqn = false )
256 {
257  const Var kWResolution( [ fractional, true_w_eqn ]( const caf::SRProxy * sr )
258  {
259  const double trueW = true_w_eqn ? kTrueW_Eqn( sr ) : kTrueW( sr );
260  const double recoW = kRecoW( sr );
261  if ( trueW > 0 && recoW > 0 )
262  {
263  double wres = recoW - trueW;
264  if ( fractional ) wres /= trueW;
265  return wres;
266  }
267  return -10.;
268  });
269  return kWResolution;
270 }
271 
273 {
274 
275  public :
276 
279  bool isMC,
281  const HistAxis& axis,
282  const Cut& cut,
283  const Var& wgtMC = kUnweighted,
284  const SystShifts& shiftMC = kNoShift )
285  {
286  fName = name;
287 
288  if ( isMC )
289  {
290  fSpectra[ "TotMC" ] = std::make_unique<Spectrum>( loader, axis, cut , shiftMC, wgtMC );
291  fSpectra[ "QEL" ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumuCC && kIsQE , shiftMC, wgtMC );
292  fSpectra[ "MEC" ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumuCC && kIsDytmanMEC, shiftMC, wgtMC );
293  fSpectra[ "RES" ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumuCC && kIsRes , shiftMC, wgtMC );
294  fSpectra[ "DIS" ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumuCC && kIsDIS , shiftMC, wgtMC );
295  fSpectra[ "Other" ] = std::make_unique<Spectrum>( loader, axis, cut && kIsOther , shiftMC, wgtMC );
296 
297  //fSpectra[ "Numu" ] = std::make_unique<Spectrum>( loader, axis, cut && kIsMuonNeutrino , shiftMC, wgtMC );
298  //fSpectra[ "Numubar" ] = std::make_unique<Spectrum>( loader, axis, cut && kIsMuonAntiNeutrino , shiftMC, wgtMC );
299 
300  /*
301  fSpectra[ "QEL_WSelected" ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumuCC && kIsQE && kWSelection , shiftMC, wgtMC );
302  fSpectra[ "MEC_WSelected" ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumuCC && kIsDytmanMEC && kWSelection, shiftMC, wgtMC );
303  fSpectra[ "RES_WSelected" ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumuCC && kIsRes && kWSelection, shiftMC, wgtMC );
304  fSpectra[ "DIS_WSelected" ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumuCC && kIsDIS && kWSelection, shiftMC, wgtMC );
305  fSpectra[ "Other_WSelected" ] = std::make_unique<Spectrum>( loader, axis, cut && kIsOther && kWSelection, shiftMC, wgtMC );
306  */
307 
308  //fSpectra[ "TotMC_MECShapeUp" ] = std::make_unique<Spectrum>( loader, axis, cut, kShiftsMECShapeUp , wgtMC );
309  //fSpectra[ "TotMC_MECShapeDown" ] = std::make_unique<Spectrum>( loader, axis, cut, kShiftsMECShapeDown, wgtMC );
310 
311  fSpectra[ "TotMC_Nominal" ] = std::make_unique<Spectrum>( loader, axis, cut , shiftMC, kPPFXFluxCVWgt );
312  fSpectra[ "QEL_Nominal" ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumuCC && kIsQE , shiftMC, kPPFXFluxCVWgt );
313  fSpectra[ "MEC_Nominal" ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumuCC && kIsDytmanMEC, shiftMC, kPPFXFluxCVWgt );
314  fSpectra[ "RES_Nominal" ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumuCC && kIsRes , shiftMC, kPPFXFluxCVWgt );
315  fSpectra[ "DIS_Nominal" ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumuCC && kIsDIS , shiftMC, kPPFXFluxCVWgt );
316  fSpectra[ "Other_Nominal" ] = std::make_unique<Spectrum>( loader, axis, cut && kIsOther , shiftMC, kPPFXFluxCVWgt );
317 
318  /*
319  fSpectra[ "Numu_Nominal" ] = std::make_unique<Spectrum>( loader, axis, cut && kIsMuonNeutrino , shiftMC, kPPFXFluxCVWgt );
320  fSpectra[ "Numubar_Nominal" ] = std::make_unique<Spectrum>( loader, axis, cut && kIsMuonAntiNeutrino , shiftMC, kPPFXFluxCVWgt );
321 
322  fSpectra[ "NonMEC_Nominal" ] = std::make_unique<Spectrum>( loader, axis, cut && !kIsDytmanMEC , shiftMC, kPPFXFluxCVWgt );
323  fSpectra[ "NonMEC_Nominal_Numu" ] = std::make_unique<Spectrum>( loader, axis, cut && !kIsDytmanMEC && kIsMuonNeutrino , shiftMC, kPPFXFluxCVWgt );
324  fSpectra[ "NonMEC_Nominal_Numubar" ] = std::make_unique<Spectrum>( loader, axis, cut && !kIsDytmanMEC && kIsMuonAntiNeutrino, shiftMC, kPPFXFluxCVWgt );
325  */
326 
327  fSpectra[ "MEC_Numu" ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumuCCMEC , shiftMC, wgtMC );
328  fSpectra[ "MEC_Numubar" ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumubarCCMEC , shiftMC, wgtMC );
329  fSpectra[ "MEC_Nominal_Numu" ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumuCCMEC , shiftMC, kPPFXFluxCVWgt );
330  fSpectra[ "MEC_Nominal_Numubar" ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumubarCCMEC , shiftMC, kPPFXFluxCVWgt );
331 
332  fSpectra[ "MEC_Valencia" ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumuCC && kIsDytmanMEC, shiftMC, kPPFXFluxCVWgt * kEmpiricalMECtoValenciaMECWgt );
333  fSpectra[ "MEC_Valencia_Numu" ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumuCCMEC , shiftMC, kPPFXFluxCVWgt * kEmpiricalMECtoValenciaMECWgt );
334  fSpectra[ "MEC_Valencia_Numubar" ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumubarCCMEC , shiftMC, kPPFXFluxCVWgt * kEmpiricalMECtoValenciaMECWgt );
335 
336  fSpectra[ "MEC_Valencia_FixNP" ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumuCC && kIsDytmanMEC, shiftMC, kPPFXFluxCVWgt * kEmpiricalMECtoValenciaMECWgt_FixNP );
337  fSpectra[ "MEC_Valencia_FixNP_Numu" ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumuCCMEC , shiftMC, kPPFXFluxCVWgt * kEmpiricalMECtoValenciaMECWgt_FixNP );
338  fSpectra[ "MEC_Valencia_FixNP_Numubar" ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumubarCCMEC , shiftMC, kPPFXFluxCVWgt * kEmpiricalMECtoValenciaMECWgt_FixNP );
339 
340  fSpectra[ "MEC_MINERvA" ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumuCC && kIsDytmanMEC, shiftMC, kPPFXFluxCVWgt * kMINERvA_Wgt_MEC );
341  fSpectra[ "MEC_MINERvA_Numu" ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumuCCMEC , shiftMC, kPPFXFluxCVWgt * kMINERvA_Wgt_MEC );
342  fSpectra[ "MEC_MINERvA_Numubar" ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumubarCCMEC , shiftMC, kPPFXFluxCVWgt * kMINERvA_Wgt_MEC );
343 
344  fSpectra[ "RES_MINOS" ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumuCC && kIsRes, shiftMC, kPPFXFluxCVWgt * kMinosRESSupp );
345  fSpectra[ "RES_MINERvA" ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumuCC && kIsRes, shiftMC, kPPFXFluxCVWgt * kMINERvALowQ2RESSupp );
346  }
347  else
348  {
349  fSpectra[ "Data" ] = std::make_unique<Spectrum>( loader, axis, cut );
350  }
351  }
352 
353  void Write() const
354  {
355  for ( const auto & s : fSpectra )
356  {
357  double pot = s.second.get()->POT();
358  std::string out_name = Form( "%s_%s", fName.c_str(), s.first.c_str() );
359  int ndim = s.second.get()->NDimensions();
360  if ( ndim == 1 )
361  {
362  s.second.get()->ToTH1( pot )->Write( out_name.c_str() );
363  }
364  else if ( ndim == 2 )
365  {
366  s.second.get()->ToTH2( pot )->Write( out_name.c_str() );
367  }
368  else throw std::runtime_error( Form( "Write operation for %d dimensions not defined", ndim ) );
369  //s.second.get()->ToTH1( s.second.get()->POT() )->Write( Form( "%s_%s", fName.c_str(), s.first.c_str() ) );
370  }
371  }
372 
373  double POT() const { return fSpectra.begin()->second.get()->POT(); }
374 
375  private :
376 
378  std::map<std::string,std::unique_ptr<Spectrum>> fSpectra;
379 };
380 
382 (
383  const std::string beam,
384  const std::string type,
385  const std::string period = "full"
386  //const std::string shift_non_mec = "",
387 )
388 {
389  if ( type != "data" && type != "mc" )
390  {
391  std::cout << "type must be \"data\" or \"mc\"" << std::endl;
392  return;
393  }
394  const bool isMC = ( type == "mc" );
395 
396  /*
397  const bool extra_stat_mc = false;
398  const bool no_smoothing = false;
399  const bool mec_refit = false;
400  LoadTunedMECWeights( shift_non_mec, extra_stat_mc, no_smoothing, mec_refit );
401  */
402 
403  std::string defn = "";
404  if ( beam == "fhc" )
405  {
406  defn = isMC ?
407  Form( "prod_sumdecaf_R17-11-14-prod4reco.d_nd_genie_nonswap_fhc_nova_v08_%s_v1_numu2018", period.c_str() ) :
408  Form( "prod_sumdecaf_R17-09-05-prod4recopreview.f_nd_numi_fhc_%s_v1_addShortSimpleCVN_goodruns_numu2018", period.c_str() );
409  }
410  else if ( beam == "rhc" )
411  {
412  defn = isMC ?
413  Form( "prod_sumdecaf_R17-11-14-prod4reco.e_nd_genie_nonswap_rhc_nova_v08_%s_v1_numu2018", period.c_str() ) :
414  Form( "prod_sumdecaf_R17-09-05-prod4recopreview.f_nd_numi_rhc_%s_v1_addShortSimpleCVN_goodruns_numu2018", period.c_str() );
415  }
416  if ( defn.empty() ) throw std::runtime_error( "SAM Definition is empty" );
417 
418  //defn = Form( "defname: %s with limit 1", defn.c_str() );
419  std::cout << "defn = " << defn << std::endl;
420 
421  SpectrumLoader loader( defn );
423 
424  std::map< std::string, HistAxis > axes;
425  axes.emplace( "RecoEhad" , HistAxis( "Reco E_{had} (GeV)", Binning::Simple( 50, 0., 2. ), kHadE ) );
426  axes.emplace( "RecoEhadVis" , HistAxis( "Visible E_{had} (GeV)", Binning::Simple( 50, 0., 1. ), kNumuHadVisE ) );
427  axes.emplace( "RecoEmu" , HistAxis( "Reco E_{#mu} (GeV)", Binning::Simple( 60, 0., 3. ), kMuE ) );
428  axes.emplace( "RecoEnu" , HistAxis( "Reco E_{#nu} (GeV)", Binning::Simple( 40, 0., 4. ), kCCE ) );
429  axes.emplace( "RecoQ2" , HistAxis( "Reco Q^{2} (GeV^{2})", Binning::Simple( 50, 0., 1. ), kRecoQ2 ) );
430  axes.emplace( "RecoQ3" , HistAxis( "Reco |#vec{q}| (GeV)", Binning::Simple( 40, 0., 2. ), kRecoQmag ) );
431  axes.emplace( "RecoW" , HistAxis( "Reco W (GeV)", Binning::Simple( 60, 0., 3. ), kRecoW ) );
432  axes.emplace( "RecoEhadFrac" , HistAxis( "Reco E_{had} Fraction", Binning::Simple( 50, 0., 1. ), kHadEFrac ) );
433  axes.emplace( "RecoCosThetaMu" , HistAxis( "Reco cos(#theta_{#mu})", Binning::Simple( 50, 0., 1. ), kCosNumi ) );
434  axes.emplace( "PrimTrkLen" , HistAxis( "Primary Track Length (m)", Binning::Simple( 60, 0., 15.), kTrkLength ) );
435 
436  std::set< std::string > vars_with_quantiles;
437  vars_with_quantiles.insert( "RecoEhad" );
438  vars_with_quantiles.insert( "RecoEhadVis" );
439  vars_with_quantiles.insert( "RecoEmu" );
440  vars_with_quantiles.insert( "RecoEnu" );
441  vars_with_quantiles.insert( "RecoQ2" );
442  vars_with_quantiles.insert( "RecoQ3" );
443  vars_with_quantiles.insert( "RecoW" );
444  vars_with_quantiles.insert( "RecoCosThetaMu" );
445 
446  if ( isMC )
447  {
448  //axes.emplace( "NucleonPair", HistAxis( "Nucleon Pair" , Binning::Simple( 5, -1, 4 ), kHitNucBin ) );
449  axes.emplace( "TrueQ2" , HistAxis( "True Q^{2} (GeV^{2})", Binning::Simple( 50, 0., 2. ), kTrueQ2 ) );
450  //axes.emplace( "TrueQ0" , HistAxis( "True q_{0} (GeV)" , Binning::Simple( 200, 0., 2. ), kTrueQ0 ) );
451  //axes.emplace( "TrueQ3" , HistAxis( "True |#vec{q}| (GeV)", Binning::Simple( 200, 0., 2. ), kTrueQ3 ) );
452  //axes.emplace( "TrueW" , HistAxis( "True W (GeV)" , Binning::Simple( 60, 0., 3. ), kTrueW ) );
453  //axes.emplace( "TrueY" , HistAxis( "True y" , Binning::Simple( 50, 0., 1. ), kTrueY ) );
454  }
455 
456  axes.emplace
457  (
458  "RecoEhadVis_vs_RecoQ3",
459  HistAxis( "Reco |#vec{q}| (GeV)", Binning::Simple( 40, 0., 2. ), kRecoQmag, "Visible E_{had} (GeV)", Binning::Simple( 50, 0., 1. ), kNumuHadVisE )
460  );
461 
462  /*
463  axes.emplace
464  (
465  "Fine_RecoEhadVis_vs_RecoQ3",
466  HistAxis( "Reco |#vec{q}| (GeV)", Binning::Simple( 200, 0., 2. ), kRecoQmag, "Visible E_{had} (GeV)", Binning::Simple( 100, 0., 1. ), kNumuHadVisE )
467  );
468  */
469 
470  if ( isMC )
471  {
472 
473  axes.emplace
474  (
475  "TrueQ0_vs_TrueQ3",
476  HistAxis( "True |#vec{q}| (GeV)", Binning::Simple( 80, 0., 2. ), kTrueQ3, "True q_{0} (GeV)", Binning::Simple( 80,0, 2. ), kTrueQ0 )
477  );
478 
479  /*
480  axes.emplace
481  (
482  "Fine_TrueQ0_vs_TrueQ3",
483  HistAxis( "True |#vec{q}| (GeV)", Binning::Simple( 160, 0., 2. ), kTrueQ3, "True q_{0} (GeV)", Binning::Simple( 160, 0., 2. ), kTrueQ0 )
484  );
485  */
486 
487  std::string title_wres_abs = "W_{reco} - W_{true} (GeV)";
488  std::string title_wres_frac = "( W_{reco} - W_{true} ) / W_{true}";
489 
490  /*
491  axes.emplace
492  (
493  "WRes_Abs",
494  HistAxis( "True W (GeV)", Binning::Simple( 100, 0., 5. ), kTrueW , title_wres_abs, Binning::Simple( 1000, -10., 10. ), GetWResolution( false, false ) )
495  );
496 
497  axes.emplace
498  (
499  "WRes_Abs_Eqn",
500  HistAxis( "True W (GeV)", Binning::Simple( 100, 0., 5. ), kTrueW_Eqn, title_wres_abs, Binning::Simple( 1000, -10., 10. ), GetWResolution( false, true ) )
501  );
502 
503  axes.emplace
504  (
505  "WRes_Frac",
506  HistAxis( "True W (GeV)", Binning::Simple( 100, 0., 5. ), kTrueW , title_wres_frac, Binning::Simple( 1000, -10., 10. ), GetWResolution( true, false ) )
507  );
508 
509  axes.emplace
510  (
511  "WRes_Frac_Eqn",
512  HistAxis( "True W (GeV)", Binning::Simple( 100, 0., 5. ), kTrueW_Eqn, title_wres_frac, Binning::Simple( 1000, -10., 10. ), GetWResolution( true, true ) )
513  );
514  */
515 
516  /*
517  axes.emplace
518  (
519  "RecoEhadVis_vs_TrueQ0",
520  HistAxis( "True q_{0} (GeV)", Binning::Simple( 40,0, 1. ), kTrueQ0, "Visible E_{had} (GeV)", Binning::Simple( 50, 0., 1. ), kNumuHadVisE )
521  );
522 
523  axes.emplace
524  (
525  "RecoQ3_vs_TrueQ3",
526  HistAxis( "True |#vec{q}| (GeV)", Binning::Simple( 40, 0., 2. ), kTrueQ3, "Reco |#vec{q}| (GeV)", Binning::Simple( 40, 0., 2. ), kRecoQmag )
527  );
528 
529  axes.emplace
530  (
531  "RecoEhad_vs_TrueQ0",
532  HistAxis( "True q_{0} (GeV)", Binning::Simple( 40, 0, 2. ), kTrueQ0, "Reco E_{had} (GeV)", Binning::Simple( 40, 0., 2. ), kHadE )
533  );
534  */
535 
536  }
537 
538  // Hadronic energy fraction quantiles
539 
540  std::string InDirQuant = "/pnfs/nova/persistent/analysis/numu/Ana2018/provisional/quantiles";
541  std::string InFileQuant = InDirQuant + "/quantiles__" + beam + "_full__numu2018.root";
542  std::cout << "\n\n --- get quantile cuts from file ---" << std::endl;
543  TFile* inFile = TFile::Open( pnfs2xrootd(InFileQuant).c_str() );
544  TH2* FDSpec2D = (TH2*)inFile->FindObjectAny("FDSpec2D");
545  std::vector<Cut> HadEFracQuantCuts = QuantileCutsFromTH2(FDSpec2D, kNumuCCOptimisedAxis, kHadEFracAxis, 4);
546  inFile->Close();
547  delete inFile;
548 
549  const Cut kSelCuts = kNumuCutND2018;
550  const Var kWeightCV = kPPFXFluxCVWgt * kXSecCVWgt2018RPAFix_noDIS;
551 
552  SystShifts kShiftsMECShapeUp;
553  kShiftsMECShapeUp.SetShift( &kMECShapeSyst2018RPAFixNu, +1 );
554  kShiftsMECShapeUp.SetShift( &kMECShapeSyst2018RPAFixAntiNu, +1 );
555 
556  SystShifts kShiftsMECShapeDown;
557  kShiftsMECShapeDown.SetShift( &kMECShapeSyst2018RPAFixNu, -1 );
558  kShiftsMECShapeDown.SetShift( &kMECShapeSyst2018RPAFixAntiNu, -1 );
559 
560  SystShifts kShiftsRPACCQESuppUp;
561  kShiftsRPACCQESuppUp.SetShift( &kRPACCQESuppSyst2019, +1 );
562 
563  SystShifts kShiftsRPACCQESuppDown;
564  kShiftsRPACCQESuppDown.SetShift( &kRPACCQESuppSyst2019, -1 );
565 
566  SystShifts kShiftsRPARESUp;
567  kShiftsRPARESUp.SetShift( &kRPARESSyst2019, +1 );
568 
569  SystShifts kShiftsRPARESDown;
570  kShiftsRPARESDown.SetShift( &kRPARESSyst2019, -1 );
571 
572  std::vector<Distribution> dists;
573  dists.reserve( 100 );
574 
575  for ( const auto & a : axes )
576  {
577  //if ( a.first != "RecoEhadVis" && a.first != "TrueQ0_vs_TrueQ3" ) continue;
578 
579  dists.emplace_back( a.first, isMC, loader, a.second, kSelCuts, kWeightCV, kNoShift );
580 
581  if ( isMC )
582  {
583  dists.emplace_back( a.first + "_MECShapeUp" , isMC, loader, a.second, kSelCuts, kWeightCV, kShiftsMECShapeUp );
584  dists.emplace_back( a.first + "_MECShapeDown", isMC, loader, a.second, kSelCuts, kWeightCV, kShiftsMECShapeDown );
585  }
586  //if ( a.first == "RecoQ2" || a.first == "RecoEnu" )
587  if ( vars_with_quantiles.find( a.first ) != vars_with_quantiles.end() )
588  {
589  for ( unsigned int icut = 0; icut < HadEFracQuantCuts.size(); ++icut )
590  {
591  std::string dist_quant_name = a.first + Form( "_EhadQuantile%d", icut + 1 );
592 
593  dists.emplace_back( dist_quant_name, isMC, loader, a.second, kSelCuts && HadEFracQuantCuts[ icut ], kWeightCV, kNoShift );
594  /*
595  dists.emplace_back( dist_quant_name + "_MECShapeUp" , isMC, loader, a.second, kSelCuts && HadEFracQuantCuts[ icut ], kWeightCV, kShiftsMECShapeUp );
596  dists.emplace_back( dist_quant_name + "_MECShapeDown" , isMC, loader, a.second, kSelCuts && HadEFracQuantCuts[ icut ], kWeightCV, kShiftsMECShapeDown );
597  dists.emplace_back( dist_quant_name + "_RPACCQESuppUp" , isMC, loader, a.second, kSelCuts && HadEFracQuantCuts[ icut ], kWeightCV, kShiftsRPACCQESuppUp );
598  dists.emplace_back( dist_quant_name + "_RPACCQESuppDown", isMC, loader, a.second, kSelCuts && HadEFracQuantCuts[ icut ], kWeightCV, kShiftsRPACCQESuppDown );
599  dists.emplace_back( dist_quant_name + "_RPARESUp" , isMC, loader, a.second, kSelCuts && HadEFracQuantCuts[ icut ], kWeightCV, kShiftsRPARESUp );
600  dists.emplace_back( dist_quant_name + "_RPARESDown" , isMC, loader, a.second, kSelCuts && HadEFracQuantCuts[ icut ], kWeightCV, kShiftsRPARESDown );
601  */
602  }
603  }
604 
605  if ( a.first == "RecoQ2" )
606  {
607  dists.emplace_back( a.first + "_WSelected", isMC, loader, a.second, kSelCuts && kWSelection, kWeightCV, kNoShift );
608  }
609  }
610 
611  /*
612  SystShifts kShiftsQELike;
613  kShiftsQELike.SetShift( &kMAQEGenieReducedSyst2018, +1 );
614  kShiftsQELike.SetShift( &kRPACCQESuppSyst2019, +1 );
615  kShiftsQELike.SetShift( &kRPACCQEEnhSyst2019, +1 );
616  kShiftsQELike.SetShift( GetGenieRwgtSyst( rwgt::fReweightCCQEPauliSupViaKF ), -1 );
617  kShiftsQELike.SetShift( GetGenieRwgtSyst( rwgt::fReweightMaCCRES ), -1 );
618  kShiftsQELike.SetShift( GetGenieRwgtSyst( rwgt::fReweightMvCCRES ), -1 );
619  //kShiftsQELike.SetShift( &kDirectRelHadEScaleSyst2017, -1 );
620  kShiftsQELike.SetShift( &kMECShapeSyst2018RPAFixNu, +1 );
621  kShiftsQELike.SetShift( &kMECShapeSyst2018RPAFixAntiNu, +1 );
622 
623  SystShifts kShiftsRESLike;
624  kShiftsRESLike.SetShift( &kMAQEGenieReducedSyst2018, -1 );
625  kShiftsRESLike.SetShift( &kRPACCQESuppSyst2019, -1 );
626  kShiftsRESLike.SetShift( &kRPACCQEEnhSyst2019, -1 );
627  kShiftsRESLike.SetShift( GetGenieRwgtSyst( rwgt::fReweightCCQEPauliSupViaKF ), +1 );
628  kShiftsRESLike.SetShift( GetGenieRwgtSyst( rwgt::fReweightMaCCRES ), +1 );
629  kShiftsRESLike.SetShift( GetGenieRwgtSyst( rwgt::fReweightMvCCRES ), +1 );
630  //kShiftsRESLike.SetShift( &kDirectRelHadEScaleSyst2017, +1 );
631  kShiftsRESLike.SetShift( &kRPARESSyst2019, +1 ); // Turns off RES RPA
632  kShiftsRESLike.SetShift( &kMECShapeSyst2018RPAFixNu, -1 );
633  kShiftsRESLike.SetShift( &kMECShapeSyst2018RPAFixAntiNu, -1 );
634 
635  std::vector< std::string > shift_vars;
636  shift_vars.push_back( "RecoEhadVis" );
637  shift_vars.push_back( "RecoQ3" );
638  shift_vars.push_back( "TrueQ0_vs_TrueQ3" );
639 
640  if ( isMC )
641  {
642  for ( const auto & v : shift_vars )
643  {
644  dists.emplace_back( v + "_ShiftsQELike" , isMC, loader, axes.at( v ), kSelCuts, kWeightCV, kShiftsQELike );
645  dists.emplace_back( v + "_ShiftsRESLike", isMC, loader, axes.at( v ), kSelCuts, kWeightCV, kShiftsRESLike );
646  }
647  }
648  */
649 
650  loader.Go();
651 
652  double pot = dists.front().POT();
653 
654  TTree potTree( "pot", "pot" );
655  potTree.Branch( Form( "pot_%s", type.c_str() ), &pot, Form( "pot_%s/D", type.c_str() ) );
656  potTree.Fill();
657 
658  std::string out_dir = "hist_files";
659  if( gSystem->AccessPathName( out_dir.c_str() ) ) gSystem->mkdir( out_dir.c_str(), true );
660 
661  std::string out_file_name = Form( "hists_2018_%s_%s_%s", beam.c_str(), type.c_str(), period.c_str() );
662  //if ( isMC && !shift_non_mec.empty() )
663  // out_file_name += "_ShiftNonMEC" + shift_non_mec;
664 
665  //TFile outputFile( Form( "%s/%s.root", out_dir.c_str(), out_file_name.c_str() ), "RECREATE" );
666  TFile outputFile( Form( "%s.root", out_file_name.c_str() ), "RECREATE" );
667 
668  potTree.Write();
669 
670  for ( const auto & d : dists ) { d.Write(); }
671 
672  outputFile.Close();
673 
674 }
675 
676 #endif
const Var kHadE
Definition: NumuVars.h:23
const Cut kIsQE
Definition: TruthCuts.cxx:104
const XML_Char * name
Definition: expat.h:151
const Cut kIsNumuCCMEC
const Cut kIsNumuCC(CCFlavSel(14, 14))
Select CC .
const Var kMINERvA_Wgt_MEC
Definition: GenieWeights.h:232
const ana::Var kRecoQ2([](const caf::SRProxy *sr){const double M_mu_sqrd=util::sqr(0.1056);double E_mu=kMuE(sr);double p_mu=sqrt(util::sqr(E_mu)-M_mu_sqrd);return 2 *kCCE(sr)*(E_mu-p_mu *kCosNumi(sr))-M_mu_sqrd;})
Reconstructed four-momentum transfer invariant (Q^2)
Definition: NumuVars.h:146
_HistAxis< Var > HistAxis
Definition: HistAxis.h:103
const NOvARwgtSyst kMECShapeSyst2018RPAFixAntiNu("MECShape2018NuRPAfixAntiNu","MEC 2018 shape (RPA fix), antineutrinos", novarwgt::kMECQ0Q3RespSystNubar2018_RPAfix)
Definition: MECSysts.h:29
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:991
correl_xv GetYaxis() -> SetDecimals()
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
const Var weight
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
Float_t x1[n_points_granero]
Definition: compare.C:5
caf::Proxy< caf::SRHeader > hdr
Definition: SRProxy.h:2096
TFile * file_mec_wgt_numu
const Var kRecoQmag([](const caf::SRProxy *sr){return sqrt(kRecoQ2(sr)+util::sqr(kHadE(sr)));})
Reconstructed three-momentum transfer.
Definition: NumuVars.h:149
const Var kEmpiricalMECtoValenciaMECWgt
Definition: GenieWeights.h:190
const Cut kIsRes
Definition: TruthCuts.cxx:111
const Cut kNumuCutND2018
Definition: NumuCuts2018.h:41
const NOvARwgtSyst kRPACCQESuppSyst2019("RPAShapesupp2019","RPA shape: low-Q^{2} suppression (2019)", novarwgt::kRPACCQESuppSyst2019)
Definition: RPASysts.h:16
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2085
Distribution(std::string name, bool isMC, SpectrumLoaderBase &loader, const HistAxis &axis, const Cut &cut, const Var &wgtMC=kUnweighted, const SystShifts &shiftMC=kNoShift)
T sqrt(T number)
Definition: d0nt_math.hpp:156
caf::Proxy< std::vector< caf::SRNeutrino > > nu
Definition: SRProxy.h:574
constexpr T pow(T x)
Definition: pow.h:75
const HistAxis kHadEFracAxis("E_{had.} / E_{#nu}", Binning::Simple(200, 0, 1), kHadEFrac)
HistAxis that implements Hadronic Energy fraction binning used by L.Vinton to derive Hadronic Energy ...
Definition: HistAxes.h:30
std::map< std::string, std::unique_ptr< Spectrum > > fSpectra
const Cut kIsAntiNu([](const caf::SRProxy *sr){if(sr->mc.nnu==0) return false;assert(sr->mc.nnu==1);return sr->mc.nu[0].pdg< 0;})
Is this truly an antineutrino?
Definition: TruthCuts.h:53
const Cut kIsOther
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
Definition: AddMC.C:8
caf::Proxy< short int > nnu
Definition: SRProxy.h:573
const Cut kIsDIS
Definition: TruthCuts.cxx:118
const Var kMINERvALowQ2RESSupp([](const caf::SRProxy *sr){if(!kIsRes(sr)) return 1.0;if(sr->mc.nu[0].tgtA==1) return 1.0; double R1=0.32;double R2=0.50;double x1=0.0;double x2=0.35;double x3=0.7;double Q2=kTrueQ2(sr);if(Q2 > x3) return 1.0;double R=R2 *(Q2-x1)*(Q2-x3)/((x2-x1)*(x2-x3))+(Q2-x1)*(Q2-x2)/((x3-x1)*(x3-x2));double weight=1-(1-R1)*pow(1-R, 2);if(weight< 0){std::cout<< "kMINERvALowQ2RESSupp = "<< weight<< std::endl;std::cout<< "Q2 = "<< Q2<< std::endl;weight=0;}return weight;})
const Cut kWSelection([](const caf::SRProxy *sr){double W=kRecoW(sr);return W > 1.2 &&W< 1.5;})
void SetSpillCut(const SpillCut &cut)
std::string pnfs2xrootd(std::string loc, bool unauth)
Definition: UtilsExt.cxx:237
const Var kTrueQ0
Definition: TruthVars.h:32
const Var kTrkLength([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-5.f;return sr->trk.kalman.tracks[0].len/100;})
Definition: NumuVars.h:65
const Var kNumuHadVisE([](const caf::SRProxy *sr){return kNumuHadCalE(sr)+kNumuHadTrkE(sr);})
Definition: NumuVars.h:124
const NuTruthVar kTrueW_Eqn_NT([](const caf::SRNeutrinoProxy *nu){const double M_p=0.938272;double W2=M_p *M_p+2 *M_p *(nu->y *nu->E)-nu->q2;if(W2 > 0) return sqrt(W2);else return-5.;})
const Var kTrueQ2
Definition: TruthVars.h:27
const Var kHitNuc
Definition: TruthVars.h:19
const Cut kIsNumubarCCMEC
const Cut kIsMuonAntiNeutrino([](const caf::SRProxy *sr){return sr->hdr.ismc &&sr->mc.nnu==1 &&sr->mc.nu[0].pdg==-14;})
const Cut kIsDytmanMEC
Definition: TruthCuts.cxx:187
ifstream inFile
Definition: AnaPlotMaker.h:34
const XML_Char * s
Definition: expat.h:262
const ana::Var GetTunedMECWeight(const bool interpolate=false)
const ana::Var GetWResolution(const bool fractional=false, const bool true_w_eqn=false)
const Var kTrueQ3
Definition: TruthVars.h:38
TH2D * hist_mec_wgt_numu
if(dump)
const Cut kCutOnHydrogen([](const caf::SRProxy *sr){return sr->mc.nu[0].tgtA==1;})
const double a
const Var kHadEFrac
Definition: NumuVars.h:24
const Var kTrueW_Eqn
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:10
const Var kCCE
Definition: NumuVars.h:21
#define pot
Float_t d
Definition: plot.C:236
void LoadTunedMECWeights(const std::string shift_non_mec="", const bool extra_stat_mc=false, const bool no_smoothing=false, const bool mec_refit=false)
#define R(x)
caf::StandardRecord * sr
virtual void Go() override
Load all the registered spectra.
const Var kEmpiricalMECtoValenciaMECWgt_FixNP([](const caf::SRProxy *sr){if(!kIsDytmanMEC(sr)) return 1.0;double weight=kEmpiricalMECtoValenciaMECWgt(sr);const double np_frac_empirical=0.89;const double np_frac_valencia=2.0/3.0;if(ana::kHitNuc(sr)==2000000201) weight *=np_frac_valencia/np_frac_empirical;else weight *=(1-np_frac_valencia)/(1-np_frac_empirical);return weight;})
loader
Definition: demo0.py:10
const NOvARwgtSyst kRPARESSyst2019("RPAShapeRES2019","RPA shape: resonance events", novarwgt::kRPARESSyst2019)
Definition: RPASysts.h:24
const HistAxis kNumuCCOptimisedAxis("Reconstructed Neutrino Energy (GeV)", kNumuCCEOptimisedBinning, kCCE)
HistAxis that implements optimised numuCCE from L. Vinton. See docdb 16332. This was close to &#39;custC&#39;...
Definition: HistAxes.h:24
const SystShifts kNoShift
Definition: SystShifts.cxx:21
OStream cout
Definition: OStream.cxx:6
void interpolate(double *buf, double *t, long int ncm, long int na, double *position, double *velocity)
Base class for the various types of spectrum loader.
const Cut kCutDISAnomaly([](const caf::SRProxy *sr){return sr->mc.nnu==0;})
caf::Proxy< caf::SRTruthBranch > mc
Definition: SRProxy.h:2097
const Var kRecoW([](const caf::SRProxy *sr){const double M_p=0.938272;double WSq=M_p *M_p+2 *M_p *(kCCE(sr)-kMuE(sr))-kRecoQ2(sr);if(WSq > 0) return sqrt(WSq);else return-5.;})
Reconstructed invariant mass (W)
Definition: NumuVars.h:152
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
const Cut cut
Definition: exporter_fd.C:30
Var VarFromNuTruthVar(const NuTruthVar &stv, double _default)
Definition: Var.cxx:7
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
TH2D * hist_mec_wgt_numubar
const Cut kIsMuonNeutrino([](const caf::SRProxy *sr){return sr->hdr.ismc &&sr->mc.nnu==1 &&sr->mc.nu[0].pdg==14;})
caf::Proxy< bool > ismc
Definition: SRProxy.h:241
std::vector< Cut > QuantileCutsFromTH2(TH2 *quantileHist, const HistAxis &independentAxis, const HistAxis &quantileAxis, const unsigned int &numQuantiles, const bool verbose)
: Do the same as the QuantileCuts function but taking in the TH2 instead of making it...
const NOvARwgtSyst kMECShapeSyst2018RPAFixNu("MECShape2018NuRPAfixNu","MEC 2018 shape (RPA fix), neutrinos", novarwgt::kMECQ0Q3RespSystNu2018_RPAfix)
Definition: MECSysts.h:28
TFile * file_mec_wgt_numubar
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
This module creates Common Analysis Files.
Definition: FileReducer.h:10
const Var kMinosRESSupp([](const caf::SRProxy *sr){if(!kIsRes(sr)||!kIsNumuCC(sr)) return 1.0;if(sr->mc.nu[0].tgtA==1) return 1.0;double Q2=kTrueQ2(sr);return 1.010/(1+exp(1-sqrt(Q2)/0.156));})
const Var kMuE
Definition: NumuVars.h:22
void make_xsec_wgts_2018_hists(const std::string beam, const std::string type, const std::string period="full")
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:114
const Var kTrueW
Definition: TruthVars.h:22
#define W(x)
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 Var kCosNumi([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks > 0 &&sr->trk.kalman.idxremid!=999){if(sr->hdr.det==1){return sr->trk.kalman.tracks[0].dir.Dot(beamDirND);}if(sr->hdr.det==2){return sr->trk.kalman.tracks[0].dir.Dot(beamDirFD);}}return-5.f;})
Definition: NumuVars.h:43
const Var kXSecCVWgt2018RPAFix_noDIS
Definition: XsecTunes.h:57