13 #include "TGraphErrors.h" 27 #include "CAFAna/Core/Var.h" 61 sys_shifts.
Shift( face_palm, weight );
64 return kWeightFromShifts;
72 if ( !isRHC && sr->
mc.
nu[0].pdg == 14 )
return true;
73 else if ( isRHC && sr->
mc.
nu[0].pdg == -14 )
return true;
85 return Form(
"/nova/ana/users/%s/mec_tuning",
getenv(
"USER" ) );
90 return isRHC ?
"numubar" :
"numu";
96 if ( !shift_non_mec.empty() ) fit_tag +=
"_shift_non_mec_" + shift_non_mec;
102 return "mec_weights_" +
GetFitTag( isRHC, shift_non_mec ) +
".root";
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" );
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" );
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;
133 double weight = hist_input_mec_weights_numu->Interpolate( q3, q0 );
135 if ( weight < 0 ) weight = 0.0;
157 return q3 > loQ3 && q3 < hiQ3 && q0 > loQ0 && q0 <
hiQ0;
183 :
ana::ISyst(shortName, latexName), fSigmaScale(sigmaScale) {};
194 : fID(sInstanceCount++),
212 double wgt = 1 + fSigmaScale *
sigma;
222 :
IRescaledSigmaSyst(
"MECInitStateCompShift",
"MEC initial state np fraction", sigmaScale )
232 const double nominalNPfrac = 0.89;
234 double newNPfrac = nominalNPfrac + 0.1 * fSigmaScale *
sigma;
236 if ( newNPfrac > 1 ) newNPfrac = 1;
237 else if ( newNPfrac < 0 ) newNPfrac = 0;
240 weight *= newNPfrac / nominalNPfrac;
242 weight *= ( 1 - newNPfrac ) / ( 1 - nominalNPfrac );
250 std::ifstream input( playlist.c_str() );
251 std::vector<std::string>
files;
254 while( getline( input, s ) )
257 if ( n_files_max > -1 && count > n_files_max )
break;
258 files.push_back( s );
272 if ( beam !=
"fhc" && beam !=
"rhc" )
throw std::runtime_error(
"Beam must be \"fhc\" or \"rhc\"" );
273 const bool isRHC = beam ==
"rhc";
278 if ( shift_non_mec.empty() )
282 else if ( shift_non_mec ==
"down" )
293 else if ( shift_non_mec ==
"up" )
304 else throw std::runtime_error( Form(
"Unrecognized shift: \"%s\"", shift_non_mec.c_str() ) );
313 const int nbinsQ3 = 20;
314 const float fitMinQ3 = 0.0;
315 const float fitMaxQ3 = 1.0;
316 const float binWidthQ3 = ( fitMaxQ3 - fitMinQ3 ) / nbinsQ3;
318 const int nbinsQ0 = 16;
319 const float fitMinQ0 = 0.0;
320 const float fitMaxQ0 = 0.8;
321 const float binWidthQ0 = ( fitMaxQ0 - fitMinQ0 ) / nbinsQ0;
323 std::vector<jw::Q3Q0Bin> binsQ3Q0;
325 std::vector<std::unique_ptr<const jw::CompNormSyst>>
systsQ3Q0;
327 for (
int i = 0;
i < nbinsQ3; ++
i )
329 for (
int j = 0;
j < nbinsQ0; ++
j )
331 double q3 = fitMinQ3 + binWidthQ3 * (
i + 0.5 );
332 double q0 = fitMinQ0 + binWidthQ0 * (
j + 0.5 );
333 if ( q3 * q3 < q0 * q0 )
continue;
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 );
341 double sig_scale = 1.0;
350 std::vector<const ana::ISyst*> systs_ptrs;
351 for (
const auto &
s : systsQ3Q0 )
353 systs_ptrs.push_back(
s.get() );
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";
363 else if ( beam ==
"rhc" )
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";
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" );
416 for (
auto & syst : systs_ptrs )
429 if ( gSystem->AccessPathName( out_dir.c_str() ) ) gSystem->mkdir( out_dir.c_str(), true );
433 std::string fit_results_file_name = out_dir +
"/mec_fit_results_" + fit_tag +
".txt";
436 for (
auto & syst : systs_ptrs )
438 fit_results_file << syst->ShortName() <<
" fitted shift = " << shifts.GetShift( syst ) <<
std::endl;
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 );
444 for ( std::size_t compIdx = 0; compIdx < binsQ3Q0.size(); compIdx++ )
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;
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" );
460 TFile
weights_file( weights_file_name.c_str(),
"recreate" );
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() ) );
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 );
473 auto h_mc_tuned = spec_fit_mc_tuned.
ToTH2( pot );
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" );
479 spec_fit_data.SaveTo(& fit_hists_file,
"spec_data" ) ;
481 h_mc_raw->Write(
"mc_raw" );
482 h_non_mec_raw->Write(
"non_mec_raw" );
483 h_mc_tuned->Write(
"mc_tuned" );
485 fit_hists_file.Close();
Near Detector underground.
Pass neutrinos through unchanged.
std::string GetFitTag(const bool isRHC, const std::string shift_non_mec)
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.
Cuts and Vars for the 2020 FD DiF Study.
const NOvARwgtSyst kRPACCQEEnhSyst2018("RPAShapeenh2018","RPA shape: higher-Q^{2} enhancement (2018)", novarwgt::kRPACCQEEnhSyst2018)
correl_xv GetYaxis() -> SetDecimals()
CompNormSyst(const ana::Cut &selCut, double sigmaScale=1.0)
Simple record of shifts applied to systematic parameters.
const Var kRecoQmag([](const caf::SRProxy *sr){return sqrt(kRecoQ2(sr)+util::sqr(kHadE(sr)));})
Reconstructed three-momentum transfer.
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 'reduced' M_A^QE shift. See documentation in NOvARwgt
Proxy for caf::StandardRecord.
Collection of SpectrumLoaders for many configurations.
caf::Proxy< std::vector< caf::SRNeutrino > > nu
Spectrum PredictSyst(osc::IOscCalc *calc, const SystShifts &shift) const override
const Var kRPAWeightRES2017
MECInitStateNPFracShift(double sigmaScale=1.0)
Encapsulate code to systematically shift a caf::SRProxy.
const Var kNumuHadVisE([](const caf::SRProxy *sr){return kNumuHadCalE(sr)+kNumuHadTrkE(sr);})
Representation of a spectrum in any variable, with associated POT.
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.
const Var kRPAWeightCCQE2018
std::string GetMECWeightsFileName(const bool isRHC, const std::string shift_non_mec)
void SetSigmaScale(double sc)
void SetSpillCut(const SpillCut &cut)
std::vector< std::string > GetFiles(const std::string playlist, const int n_files_max=-1)
SpectrumLoaderBase & GetLoader(caf::Det_t det, DataMC datamc, DataSource src=kBeam, SwappingConfig swap=kNonSwap)
Retrieve a specific loader.
std::string getenv(std::string const &name)
Optimized version of OscCalcPMNS.
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))
const SystShifts kNoShift
caf::Proxy< caf::SRTruthBranch > mc
const ana::Var GetWeightFromShifts(const ana::SystShifts &sys_shifts)
void Shift(double sigma, caf::SRProxy *sr, double &weight) const override
Perform the systematic shift.
void Shift(double sigma, caf::SRProxy *sr, double &weight) const override
Perform the systematic shift.
TH2D * hist_input_mec_weights_numu
std::string GetMECTuningDirectory()
string shortName
THUMBNAIL BLOCK: We need to make a thumbnail for each.
double GetSigmaScale() const
const NOvARwgtSyst kRPACCQESuppSyst2018("RPAShapesupp2018","RPA shape: low-Q^{2} suppression (2018)", novarwgt::kRPACCQESuppSyst2018)
const SpillCut kStandardSpillCuts
Apply this unless you're doing something special.
void mec_tuning(const std::string beam, const std::string shift_non_mec="")
std::string to_string(ModuleType mt)
static unsigned int sInstanceCount
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...
void SetLoaderPath(const std::string &path, caf::Det_t det, DataMC datamc, DataSource src=kBeam, SwappingConfig swap=kNonSwap)
Configure loader via wildcard path.
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
const Var kRescaleHighWDIS
const Var kUnweighted
The simplest possible Var, always 1. Used as a default weight.
void SetShift(const ISyst *syst, double shift, bool force=false)
const ana::Cut GetCutIsFitMEC(const bool isRHC)
IRescaledSigmaSyst(const std::string &shortName, const std::string &latexName, double sigmaScale=1.0)
ana::Cut Q3Q0CutFactory(float loQ3, float hiQ3, float loQ0, float hiQ0)
void LoadWeightsTunedNumuMEC(const std::string shift_non_mec)
Compare a single data spectrum to the MC + cosmics expectation.
std::string GetNeutrinoTag(const bool isRHC)
void Shift(caf::SRProxy *sr, double &weight) const
float h_data[xbins][ybins]
const NOvARwgtSyst kRPARESSyst2018("RPAShapeRES2018","RPA shape: resonance events (2018)", novarwgt::kRPARESSyst2018)
Perform MINUIT fits in one or two dimensions.