mec_tuning_preds_2020.C
Go to the documentation of this file.
1 /*
2  * Fill the useful spectra and prediction for mec tuning fit
3  * either with 2D gaussian (MINERvA-style),
4  * two 2D gaussians with / without overall normalization
5  * or lots of normalization bin shifts (like 2018/9 analyses )
6  *
7  */
8 
9 #include "TCanvas.h"
10 #include "TGraphErrors.h"
11 #include "TH1D.h"
12 #include "TH2D.h"
13 #include "THStack.h"
14 #include "TLegend.h"
15 
16 #include "CAFAna/Core/Loaders.h"
17 #include "CAFAna/Core/Spectrum.h"
23 #include "CAFAna/Core/Var.h"
24 #include "CAFAna/Cuts/SpillCuts.h"
25 #include "CAFAna/Cuts/TruthCuts.h"
33 
34 #include "CAFAna/Core/ISyst.h"
35 #include "CAFAna/Systs/Systs.h"
36 #include "CAFAna/Systs/XSecSysts.h"
37 #include "CAFAna/Systs/RPASysts.h"
38 #include "CAFAna/Systs/RESSysts.h"
39 
40 #include "OscLib/OscCalcPMNSOpt.h"
42 
44 
47 
48 #include <fstream>
49 
51  {
52  return Form( "/nova/ana/users/%s/mec_tuning_2020_systs", getenv( "USER" ) );
53  }
54 
57  {
58  std::string numu_mec_weights_file_name = "/pnfs/nova/persistent/users/mcasales/prod5/mec_fit_hists_fhc_bins_2020_chis100_ext.root";
59  //miniprod5/mec_fit_hists_fhc_bins_mp5.root" ;
60  std::cout << "Loading FHC MEC weights from file:\n" << numu_mec_weights_file_name << std::endl;
61  TFile *numu_mec_weights_file = TFile::Open( ana::pnfs2xrootd(numu_mec_weights_file_name).c_str());
62  hist_input_mec_weights_numu = (TH2D*)numu_mec_weights_file->Get( "mec_weights_smoothed" );
63  }
64 
65  const ana::Var kWeightTunedNumuMEC( []( const caf::SRProxy * sr )
66  {
67  if ( !ana::kIsDytmanMEC( sr ) || !sr->mc.nu[0].iscc ) return 1.0;
68  if ( sr->mc.nu[0].pdg != 14 ) return 1.0;
69  if ( hist_input_mec_weights_numu == NULL ) throw std::runtime_error( "FHC MEC true (q0,q3) weights histogram is NULL" );
70 
71  double q0 = ana::kTrueQ0( sr );
72  double q3 = ana::kTrueQ3( sr );
73 
74  if ( q0 < hist_input_mec_weights_numu->GetYaxis()->GetXmin() ||
75  q0 > hist_input_mec_weights_numu->GetYaxis()->GetXmax() ||
76  q3 < hist_input_mec_weights_numu->GetXaxis()->GetXmin() ||
77  q3 > hist_input_mec_weights_numu->GetXaxis()->GetXmax() ) return 1.0;
78 
79  double weight = hist_input_mec_weights_numu->Interpolate( q3, q0 );
80  //double weight = hist_input_mec_weights_numu->GetBinContent( hist_input_mec_weights_numu->FindFixBin( q3, q0 ) );
81  if ( weight < 0 ) weight = 0.0;
82  return weight;
83  });
84 
85 
86  const ana::Cut kNumu2020PIDLoosePTP( [](const caf::SRProxy* sr)
87  {
88  return (sr->sel.remid.pid > 0.7 && sr->sel.cvnloosepreselptp.numuid > 0.82);
89  });
90 
91 
92 
94 (
95  const std::string beam = "rhc", // so far only fhc
96  const int test = 1,
97  const std::string def="mcnp",
98  const std::string fit="doublegaussB", // gauss (for 2D gaussian only), gausssupp (2D gaussian + RPA/QE suppression) , bins (2018-style 200-norm bins)
99  const std::string shift_non_mec = "" //up, down
100 )
101 {
102  if ( beam != "fhc" && beam != "rhc" ) throw std::runtime_error( "Beam must be \"fhc\" or \"rhc\"" );
103  if ( fit!="gausssupp" && fit != "gauss" && fit != "bins" && fit != "doublegauss" && fit != "doublegaussB"
104  && fit != "doublegaussB1" && fit!= "doublegaussSupp" && fit!= "doublegaussSupp2" ) throw std::runtime_error( "Fit must be \"gauss\", \"bins\" or \"doublegauss\" " );
105  const bool isRHC = beam == "rhc";
106  const bool MECbins= isRHC && fit=="bins";
107  if ( MECbins ) LoadWeightsTunedNumuMEC();
108  const bool Gauss2 = ( fit == "doublegauss" || fit== "doublegaussSupp" || fit=="doublegaussSupp2" || fit=="doublegaussB" || fit=="doublegaussB1" );
109  const bool Gauss = fit =="gauss";
110 
111  ana::SystShifts kShiftsNonMEC;
112  if ( shift_non_mec != "down" && shift_non_mec != "up" )
113  {
114  kShiftsNonMEC = ana::kNoShift;
115  }
116  else if ( shift_non_mec == "down" )
117  {
118  kShiftsNonMEC.SetShift( ana::GetGenieKnobSyst(rwgt::fReweightZNormCCQE), +1);
119  kShiftsNonMEC.SetShift( ana::GetGenieKnobSyst(rwgt::fReweightZExpA1CCQE), +1);
120  kShiftsNonMEC.SetShift( ana::GetGenieKnobSyst(rwgt::fReweightZExpA2CCQE), +1);
121  kShiftsNonMEC.SetShift( ana::GetGenieKnobSyst(rwgt::fReweightZExpA3CCQE), +1);
122  kShiftsNonMEC.SetShift( ana::GetGenieKnobSyst(rwgt::fReweightZExpA4CCQE), +1);
123  kShiftsNonMEC.SetShift( &ana::kRPACCQESuppSyst2020, +1 );
124  kShiftsNonMEC.SetShift( &ana::kRPACCQEEnhSyst2020, +1 );
125  kShiftsNonMEC.SetShift( ana::GetGenieKnobSyst( rwgt::fReweightMaCCRES ), -1 );
126  kShiftsNonMEC.SetShift( ana::GetGenieKnobSyst( rwgt::fReweightMvCCRES ), -1 );
127  kShiftsNonMEC.SetShift( &ana::kRESLowQ2SuppressionSyst2020, +1 );
128  }
129  else if ( shift_non_mec == "up" )
130  {
131  kShiftsNonMEC.SetShift( ana::GetGenieKnobSyst(rwgt::fReweightZNormCCQE), -1);
132  kShiftsNonMEC.SetShift( ana::GetGenieKnobSyst(rwgt::fReweightZExpA1CCQE), -1);
133  kShiftsNonMEC.SetShift( ana::GetGenieKnobSyst(rwgt::fReweightZExpA2CCQE), -1);
134  kShiftsNonMEC.SetShift( ana::GetGenieKnobSyst(rwgt::fReweightZExpA3CCQE), -1);
135  kShiftsNonMEC.SetShift( ana::GetGenieKnobSyst(rwgt::fReweightZExpA4CCQE), -1);
136  kShiftsNonMEC.SetShift( &ana::kRPACCQESuppSyst2020, -1 );
137  kShiftsNonMEC.SetShift( &ana::kRPACCQEEnhSyst2020, -1 );
138  kShiftsNonMEC.SetShift( ana::GetGenieKnobSyst( rwgt::fReweightMaCCRES ), +1 );
139  kShiftsNonMEC.SetShift( ana::GetGenieKnobSyst( rwgt::fReweightMvCCRES ), +1 );
140  kShiftsNonMEC.SetShift( &ana::kRESLowQ2SuppressionSyst2020, -1 );
141  }
142 
143  else throw std::runtime_error( Form( "Unrecognized shift: \"%s\"", shift_non_mec.c_str() ) );
144 
145  const ana::Var wgt_non_mec = ana::kUnweighted; // * ana::kMinosResSupp;
146  const ana::Var wgt_shift_non_mec = !shift_non_mec.empty() ? ana::GetWeightFromShifts( kShiftsNonMEC ) : ana::kUnweighted;
147  const ana::Var wgt_mec_numu = MECbins ? kWeightTunedNumuMEC : ana::kUnweighted;
148  const ana::Var wgt_mec_2gauss = Gauss2 ? ana::kMECDoubleGaussEnh : ana::kUnweighted;
149  const ana::Var wgt_mec_2gauss_UP = (Gauss2 && shift_non_mec=="up") ? ana::kMECDoubleGaussEnhUP : ana::kUnweighted;
150  const ana::Var wgt_mec_2gauss_DOWN = (Gauss2 && shift_non_mec=="down") ? ana::kMECDoubleGaussEnhDOWN : ana::kUnweighted;
151  const ana::Var wgt_mec_gauss = Gauss ? ana::kMECGaussEnh : ana::kUnweighted;
152  const ana::Var wgt_fit = ana::kPPFXFluxCVWgt * wgt_non_mec * wgt_shift_non_mec * wgt_mec_numu
153  * wgt_mec_gauss * wgt_mec_2gauss * wgt_mec_2gauss_UP * wgt_mec_2gauss_DOWN;
154 
156 // const ana::Cut kSelCuts = ana::kNumu2020ND;
157  const ana::Cut kSelCuts = kNumu2020NDLoosePTP;
158 
159 // --------------------------- could be more adaptable -----------------------------------------
160  std::vector<const ana::ISyst*> systs_ptrs;
161  // need a less awkward way to declare this // extended range
162  std::vector<std::unique_ptr<const ana::CompNormSyst>> systs= ana::systsQ3Q0(24, 0.0, 1.2, 24, 0.0, 1.2,
163  ana::GetCutIsFitMEC( isRHC ) && kSelCuts);
164  if (fit=="bins"){
165  for ( const auto & s : systs ){
166  systs_ptrs.push_back( s.get() );
167  }
168  }
169  if (fit=="gauss"){ systs_ptrs=ana::systs_params_gauss();}
170  if (fit=="gausssupp"){ systs_ptrs=ana::systs_params_gauss_supp(); }
171  if (fit=="doublegauss"){ systs_ptrs=ana::systs_params_double_gauss(); }
172  if (fit=="doublegaussB" || fit=="doublegaussB1"){
173  if (shift_non_mec=="up") systs_ptrs=ana::systs_params_double_gauss_extraUP();
174  else if (shift_non_mec=="down") systs_ptrs=ana::systs_params_double_gauss_extraDOWN();
175  else systs_ptrs=ana::systs_params_double_gauss_extra();
176  }
177  if (fit== "doublegaussSupp") { systs_ptrs= ana::systs_params_double_gauss_supp();}
178  if (fit== "doublegaussSupp2"){ systs_ptrs= ana::systs_params_double_gauss_supp();}
179 
180  std::cout << "Fitting (systs_ptrs.size()) "<< systs_ptrs.size()<< " systematic parameters" << std::endl;
181 // -------------------------------------------------------------------------------------------
182 // datasets
183  int strided=0;
184  int stridem=0;
185  std::string defn_data = "";
186  std::string defn_mc = "";
187 if ( beam == "fhc" )
188  {
189  defn_data = "defname: prod_sumdecaf_R19-11-18-prod5reco.d.f.h.l_nd_numi_fhc_full_v1_goodruns_numu2020";
190  defn_mc = "def_snapshot prod_sumdecaf_R19-11-18-prod5reco.d.h.l_nd_genie_N1810j0211a_nonswap_fhc_nova_v08_full_v1_numu2020";
191  //"def_snapshot prod_caf_R19-11-18-prod5reco.d_nd_genie_N1810j0211a_nonswap_fhc_nova_v08_full_v1";
192  // if(fit != "bins" ) strided=0; stridem=5;
193  if (def=="mcnp") defn_mc = "defname: prod_caf_R19-11-18-prod5reco.muremove-hotfix.a_nd_genie_N1810j0211a_nonswap_genierw_fhc_nova_v08_full_mcnp_v1";
194  // defn_mc = "defname: prod_caf_R19-11-18-prod5reco.muremove-hotfix.a_nd_genie_N1810j0211a_nonswap_fhc_nova_v08_full_mcnp_v1";
195  }
196  else if ( beam == "rhc" )
197  {
198  defn_data = "def_snapshot prod_sumdecaf_R19-11-18-prod5reco.g_nd_numi_rhc_full_v1_goodruns_numu2020";
199  defn_mc = "def_snapshot prod_sumdecaf_R19-11-18-prod5reco.d_nd_genie_N1810j0211a_nonswap_rhc_nova_v08_full_v1_numu2020";
200  if (def=="mcnp") defn_mc = "def_snapshot prod_caf_R19-11-18-prod5reco.muremove-hotfix.a_nd_genie_N1810j0211a_nonswap_genierw_rhc_nova_v08_full_mcnp_v1";
201  //defn_mc = "def_snapshot prod_caf_R19-11-18-prod5reco.d_nd_genie_N1810j0211a_nonswap_rhc_nova_v08_full_v1";
202  //if (fit !="bins" ) strided=0; stridem=5;
203  }
204 
205  defn_data += test ? " with limit 1" : " with stride "+ std::to_string(strided);
206  defn_mc += test ? " with limit 1" : " with stride "+std::to_string(stridem);
207 
208  if ( defn_data.empty() ) throw std::runtime_error( "Data SAM Definition is empty" );
209  if ( defn_mc .empty() ) throw std::runtime_error( "MC SAM Definition is empty" );
210 
212  loaders.SetLoaderPath( defn_mc , caf::kNEARDET, ana::Loaders::kMC );
213  loaders.SetLoaderPath( defn_data, caf::kNEARDET, ana::Loaders::kData );
215 
216  //Define range for data-mc agreement comparison;
217  double q0_max=0;
218  double q3_max=0;
219  double q0_min=0;
220  double q3_min=0;
221  int binsq0=0;
222  int binsq3=0;
223  if (fit=="bins" ){
224  q3_max=1.4; binsq3=28;
225  q0_max=0.4; binsq0=20;
226  }
227  if (fit=="gauss" || fit=="gausssupp" || fit=="doublegauss" || fit=="doublegaussB" || fit=="doublegaussSupp" || fit=="doublegaussSupp2"){
228  q3_max=1.4; binsq3=28;
229  q0_max=0.4; binsq0=20;//-1;
230  // q0_min=0.02;
231  }
232  // removing first bin of ehad which never agrees >:(
233  if (fit=="doublegaussB1"){
234  q3_max=1.4; binsq3=28;
235  q0_max=0.4; binsq0=20-2;
236  q0_min=0.02*2;
237  }
238  ana::HistAxis axis_fit( "Reco |#vec{q}| (GeV)" , ana::Binning::Simple( binsq3, q3_min, q3_max ), ana::kRecoQmag,
239  "Reco E_{had, vis} (GeV)", ana::Binning::Simple( binsq0, q0_min, q0_max ), ana::kNumuHadVisE );
240 
241  ana::Spectrum spec_fit_data
242  (
244  axis_fit,
245  kSelCuts
246  );
247 
248  ana::Spectrum spec_fit_mc_raw
249  (
251  axis_fit,
252  kSelCuts,
254  wgt_fit
255  );
256 
257  ana::Spectrum spec_fit_non_mec_raw
258  (
260  axis_fit,
261  kSelCuts && !ana::kIsDytmanMEC,
262  ana::kNoShift,
263  wgt_fit
264  );
265 
266  ana::Spectrum spec_fit_mec_raw
267  (
269  axis_fit,
270  kSelCuts && ana::kIsDytmanMEC,
271  ana::kNoShift,
272  wgt_fit
273  );
274 
275  ana::Spectrum spec_fit_mc_flux
276  (
278  axis_fit,
279  kSelCuts,
282  );
283 
284 
285  osc::NoOscillations no_osc_calc;
286  ana::NoOscPredictionGenerator pred_gen( loaders.GetLoader(caf::kNEARDET,ana::Loaders::kMC), axis_fit, kSelCuts , wgt_fit);
287  // create a prediction that only has MEC events to see the how MEC only changes
288  ana::NoOscPredictionGenerator pred_genMEC( loaders.GetLoader(caf::kNEARDET,ana::Loaders::kMC), axis_fit, kSelCuts && ana::kIsDytmanMEC , wgt_fit);
289 
290  ana::PredictionInterp pred_interp( systs_ptrs, &no_osc_calc, pred_gen, loaders );
291  ana::PredictionInterp pred_interpMEC( systs_ptrs, &no_osc_calc, pred_genMEC, loaders );
292 
293  loaders.Go();
294 
295  const std::string out_dir = GetMECTuningDirectory();
296  if ( gSystem->AccessPathName( out_dir.c_str() ) ) gSystem->mkdir( out_dir.c_str(), true );
297 
298  std::string file_name = "mec_tuning_predictions_"+beam+"_"+fit+"_"+shift_non_mec+"2020";
299  if (def=="mcnp") file_name = "mec_tuning_predictions_"+beam+"_"+fit+"_"+shift_non_mec+"2020_mcnp_genierw";
300  if(test==1) file_name=out_dir +"/"+ file_name + "_test";
301  file_name=file_name + ".root";
302 
303  TFile file( file_name.c_str(), "recreate" );
304  file.cd();
305  spec_fit_data.SaveTo( & file, "spec_data" ) ;
306  spec_fit_non_mec_raw.SaveTo( & file,"spec_mc_nonmec_raw") ;
307  spec_fit_mec_raw.SaveTo( & file,"spec_mc_mec_raw") ;
308  spec_fit_mc_raw.SaveTo( & file,"spec_mc_raw") ;
309  pred_interp.SaveTo( & file,"pred_interp") ;
310  pred_interpMEC.SaveTo( & file,"pred_interpMEC") ;
311  spec_fit_mc_flux.SaveTo( & file, "spec_mc_flux");
312 
313  // store your parameters
314  TH1D * init_h = new TH1D( "initial pars", "Gaussian parameters",
315  systs_ptrs.size() , 0., systs_ptrs.size());
316 
317 if (fit=="doublegaussB1" || fit=="doublegaussB" || fit=="doublegauss" || fit=="doublegaussSupp" || fit=="doublegaussSupp2" ){
318  int bin=1;
319  std::string initial_pars_file_name = out_dir + "/mec_fit_initial_" + beam + "_" + fit + "_"+shift_non_mec+"2020.txt";//_bigchisq.txt";
320  if (def=="mcnp") initial_pars_file_name = out_dir + "/mec_fit_initial_" + beam + "_" + fit + "_"+shift_non_mec+"2020_mcnp.txt";
321  std::ofstream fit_initial_file( initial_pars_file_name.c_str(), std::ofstream::out );
322  for ( auto & syst : systs_ptrs )
323  {
324  fit_initial_file << syst->ShortName() << "old param = " <<
325  ana::CalcMECDoubleGaussEnhShiftedParam( syst->ShortName(), 0 , shift_non_mec)<< std::endl;
326  init_h->SetBinContent(bin, ana::CalcMECDoubleGaussEnhShiftedParam( syst->ShortName(),0, shift_non_mec));
327  init_h->GetXaxis()->SetBinLabel(bin,(syst->ShortName()).c_str());
328  bin++;
329  }
330  }
331 
332  init_h->Write("initial_gauss");
333 
334  file.Close();
335  std::cout << "Wrote file: "<< file_name << std::endl;
336 
337 
338 }
caf::Proxy< caf::SRCVNResult > cvnloosepreselptp
Definition: SRProxy.h:1254
Near Detector underground.
Definition: SREnums.h:10
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;})
Pass neutrinos through unchanged.
Definition: IOscCalc.h:45
Implements systematic errors by interpolation between shifted templates.
std::vector< const ISyst * > systs_params_double_gauss_extra()
double CalcMECDoubleGaussEnhShiftedParam(std::string shift_param, double shift_sigma, std::string shifted)
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
correl_xv GetYaxis() -> SetDecimals()
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
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);})
const Var weight
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
std::vector< const ISyst * > systs_params_double_gauss()
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
Collection of SpectrumLoaders for many configurations.
Definition: Loaders.h:23
caf::Proxy< std::vector< caf::SRNeutrino > > nu
Definition: SRProxy.h:618
caf::Proxy< float > pid
Definition: SRProxy.h:1136
std::string pnfs2xrootd(std::string loc, bool unauth)
Definition: UtilsExt.cxx:237
const ana::Cut kNumu2020PIDLoosePTP([](const caf::SRProxy *sr){return(sr->sel.remid.pid > 0.7 &&sr->sel.cvnloosepreselptp.numuid > 0.82);})
const Var kTrueQ0
Definition: TruthVars.h:32
std::vector< const ISyst * > systs_params_double_gauss_extraDOWN()
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);})
const Var kNumuHadVisE([](const caf::SRProxy *sr){return kNumuHadCalE(sr)+kNumuHadTrkE(sr);})
Definition: NumuVars.h:124
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
std::string GetMECTuningDirectory()
void Go()
Call Go() on all the loaders.
Definition: Loaders.cxx:162
const Cut kIsDytmanMEC
Definition: TruthCuts.cxx:187
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);})
const XML_Char * s
Definition: expat.h:262
void SetSpillCut(const SpillCut &cut)
Definition: Loaders.cxx:121
virtual void SaveTo(TDirectory *dir, const std::string &name) const override
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))
std::vector< const ISyst * > systs_params_gauss()
void LoadWeightsTunedNumuMEC()
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)
std::vector< const ISyst * > systs_params_double_gauss_supp()
caf::StandardRecord * sr
caf::Proxy< caf::SRRemid > remid
Definition: SRProxy.h:1269
float bin[41]
Definition: plottest35.C:14
const NOvARwgtSyst kRPACCQEEnhSyst2020("RPAShapeenh2020","RPA shape: higher-Q^{2} enhancement (2020)", novarwgt::kRPACCQEEnhSyst2020)
Definition: RPASysts.h:17
const NOvARwgtSyst kRESLowQ2SuppressionSyst2020("LowQ2RESSupp2020","RES low-Q^2 suppression", novarwgt::kRESLowQ2SuppressionSyst2020)
Definition: RESSysts.h:10
std::vector< const ISyst * > systs_params_double_gauss_extraUP()
const SystShifts kNoShift
Definition: SystShifts.cxx:21
OStream cout
Definition: OStream.cxx:6
caf::Proxy< caf::SRTruthBranch > mc
Definition: SRProxy.h:2138
void mec_tuning_preds_2020(const std::string beam="rhc", const int test=1, const std::string def="mcnp", const std::string fit="doublegaussB", const std::string shift_non_mec="")
const Cut GetCutIsFitMEC(const bool isRHC)
const NOvARwgtSyst kRPACCQESuppSyst2020("RPAShapesupp2020","RPA shape: low-Q^{2} suppression (2020)", novarwgt::kRPACCQESuppSyst2020)
Definition: RPASysts.h:18
std::vector< Loaders * > loaders
Definition: syst_header.h:386
const Var kMECDoubleGaussEnh([](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 CalcMECDoubleGaussEnh(q0, q3, Q2, kGauss2DNorm_1, 0);})
TFile * file
Definition: cellShifts.C:17
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
caf::Proxy< caf::SRIDBranch > sel
Definition: SRProxy.h:2141
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
const Cut kNumuQuality
Definition: NumuCuts.h:18
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 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
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
caf::Proxy< float > numuid
Definition: SRProxy.h:907
std::vector< const ISyst * > systs_params_gauss_supp()
TH2D * hist_input_mec_weights_numu
const Var GetWeightFromShifts(const SystShifts &sys_shifts)
enum BeamMode string