16 #include "CAFAna/Core/Binning.h" 20 #include "3FlavorAna/Core/HistAxes.h" 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 )
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() ) );
58 const bool finer_bins =
true;
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";
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() );
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 );
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";
100 if ( sr->
mc.
nu[0].pdg > 0 )
102 if (
hist_mec_wgt_numu == NULL )
throw std::runtime_error(
"Numu MEC weights histogram is NULL" );
104 if ( q0 < hist_mec_wgt_numu->
GetYaxis()->GetXmin() ||
115 if (
hist_mec_wgt_numubar == NULL )
throw std::runtime_error(
"Numubar MEC weights histogram is NULL" );
117 if ( q0 < hist_mec_wgt_numubar->
GetYaxis()->GetXmin() ||
126 if ( weight < 0 ) weight = 0.0;
129 return kWeightTunedMEC;
135 if ( sr->
mc.
nu[0].tgtA == 1 )
return 1.0;
137 return 1.010 / ( 1 +
exp( 1 -
sqrt(
Q2 ) / 0.156 ) );
142 if ( !
kIsRes( sr ) )
return 1.0;
143 if ( sr->
mc.
nu[0].tgtA == 1 )
return 1.0;
155 if (
Q2 > x3 )
return 1.0;
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 );
185 const double np_frac_empirical = 0.89;
186 const double np_frac_valencia = 2.0 / 3.0;
189 weight *= np_frac_valencia / np_frac_empirical;
191 weight *= ( 1 - np_frac_valencia ) / ( 1 - np_frac_empirical );
198 return sr->
mc.
nnu == 0;
203 return sr->
mc.
nu[0].tgtA == 1;
221 return W > 1.2 && W < 1.5;
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 );
257 const Var kWResolution( [ fractional, true_w_eqn ](
const caf::SRProxy * sr )
260 const double recoW =
kRecoW( sr );
261 if ( trueW > 0 && recoW > 0 )
263 double wres = recoW - trueW;
264 if ( fractional ) wres /= trueW;
290 fSpectra[
"TotMC" ] = std::make_unique<Spectrum>(
loader,
axis,
cut , shiftMC, wgtMC );
295 fSpectra[
"Other" ] = std::make_unique<Spectrum>(
loader,
axis, cut &&
kIsOther , shiftMC, wgtMC );
327 fSpectra[
"MEC_Numu" ] = std::make_unique<Spectrum>(
loader,
axis, cut &&
kIsNumuCCMEC , shiftMC, wgtMC );
349 fSpectra[
"Data" ] = std::make_unique<Spectrum>(
loader,
axis,
cut );
355 for (
const auto &
s : fSpectra )
357 double pot =
s.second.get()->POT();
359 int ndim =
s.second.get()->NDimensions();
362 s.second.get()->ToTH1( pot )->Write( out_name.c_str() );
364 else if ( ndim == 2 )
366 s.second.get()->ToTH2( pot )->Write( out_name.c_str() );
368 else throw std::runtime_error( Form(
"Write operation for %d dimensions not defined", ndim ) );
373 double POT()
const {
return fSpectra.begin()->second.get()->POT(); }
378 std::map<std::string,std::unique_ptr<Spectrum>>
fSpectra;
389 if ( type !=
"data" && type !=
"mc" )
394 const bool isMC = ( type ==
"mc" );
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() );
410 else if ( beam ==
"rhc" )
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() );
416 if ( defn.empty() )
throw std::runtime_error(
"SAM Definition is empty" );
424 std::map< std::string, HistAxis >
axes;
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" );
458 "RecoEhadVis_vs_RecoQ3",
487 std::string title_wres_abs =
"W_{reco} - W_{true} (GeV)";
488 std::string title_wres_frac =
"( W_{reco} - W_{true} ) / W_{true}";
540 std::string InDirQuant =
"/pnfs/nova/persistent/analysis/numu/Ana2018/provisional/quantiles";
541 std::string InFileQuant = InDirQuant +
"/quantiles__" + beam +
"_full__numu2018.root";
544 TH2* FDSpec2D = (TH2*)inFile->FindObjectAny(
"FDSpec2D");
572 std::vector<Distribution> dists;
573 dists.reserve( 100 );
575 for (
const auto &
a : axes )
579 dists.emplace_back(
a.first, isMC, loader,
a.second, kSelCuts, kWeightCV,
kNoShift );
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 );
587 if ( vars_with_quantiles.find(
a.first ) != vars_with_quantiles.end() )
589 for (
unsigned int icut = 0; icut < HadEFracQuantCuts.size(); ++icut )
591 std::string dist_quant_name =
a.first + Form(
"_EhadQuantile%d", icut + 1 );
593 dists.emplace_back( dist_quant_name, isMC, loader,
a.second, kSelCuts && HadEFracQuantCuts[ icut ], kWeightCV,
kNoShift );
605 if (
a.first ==
"RecoQ2" )
607 dists.emplace_back(
a.first +
"_WSelected", isMC, loader,
a.second, kSelCuts &&
kWSelection, kWeightCV,
kNoShift );
652 double pot = dists.front().POT();
654 TTree potTree(
"pot",
"pot" );
655 potTree.Branch( Form(
"pot_%s", type.c_str() ), &pot, Form(
"pot_%s/D", type.c_str() ) );
659 if( gSystem->AccessPathName( out_dir.c_str() ) ) gSystem->mkdir( out_dir.c_str(), true );
661 std::string out_file_name = Form(
"hists_2018_%s_%s_%s", beam.c_str(), type.c_str(),
period.c_str() );
666 TFile outputFile( Form(
"%s.root", out_file_name.c_str() ),
"RECREATE" );
670 for (
const auto &
d : dists ) {
d.Write(); }
const Cut kIsNumuCC(CCFlavSel(14, 14))
Select CC .
const Var kMINERvA_Wgt_MEC
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)
_HistAxis< Var > HistAxis
const NOvARwgtSyst kMECShapeSyst2018RPAFixAntiNu("MECShape2018NuRPAfixAntiNu","MEC 2018 shape (RPA fix), antineutrinos", novarwgt::kMECQ0Q3RespSystNubar2018_RPAfix)
Cuts and Vars for the 2020 FD DiF Study.
double Q2(const Interaction *const i)
correl_xv GetYaxis() -> SetDecimals()
Proxy for caf::SRNeutrino.
Simple record of shifts applied to systematic parameters.
Float_t x1[n_points_granero]
caf::Proxy< caf::SRHeader > hdr
TFile * file_mec_wgt_numu
const Var kRecoQmag([](const caf::SRProxy *sr){return sqrt(kRecoQ2(sr)+util::sqr(kHadE(sr)));})
Reconstructed three-momentum transfer.
const Var kEmpiricalMECtoValenciaMECWgt
const NOvARwgtSyst kRPACCQESuppSyst2019("RPAShapesupp2019","RPA shape: low-Q^{2} suppression (2019)", novarwgt::kRPACCQESuppSyst2019)
Proxy for caf::StandardRecord.
Distribution(std::string name, bool isMC, SpectrumLoaderBase &loader, const HistAxis &axis, const Cut &cut, const Var &wgtMC=kUnweighted, const SystShifts &shiftMC=kNoShift)
caf::Proxy< std::vector< caf::SRNeutrino > > nu
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 ...
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?
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
caf::Proxy< short int > nnu
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)
const Var kTrkLength([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-5.f;return sr->trk.kalman.tracks[0].len/100;})
const Var kNumuHadVisE([](const caf::SRProxy *sr){return kNumuHadCalE(sr)+kNumuHadTrkE(sr);})
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 Cut kIsNumubarCCMEC
const Cut kIsMuonAntiNeutrino([](const caf::SRProxy *sr){return sr->hdr.ismc &&sr->mc.nnu==1 &&sr->mc.nu[0].pdg==-14;})
const ana::Var GetTunedMECWeight(const bool interpolate=false)
const ana::Var GetWResolution(const bool fractional=false, const bool true_w_eqn=false)
const Cut kCutOnHydrogen([](const caf::SRProxy *sr){return sr->mc.nu[0].tgtA==1;})
fvar< T > exp(const fvar< T > &x)
void LoadTunedMECWeights(const std::string shift_non_mec="", const bool extra_stat_mc=false, const bool no_smoothing=false, const bool mec_refit=false)
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;})
const NOvARwgtSyst kRPARESSyst2019("RPAShapeRES2019","RPA shape: resonance events", novarwgt::kRPARESSyst2019)
const HistAxis kNumuCCOptimisedAxis("Reconstructed Neutrino Energy (GeV)", kNumuCCEOptimisedBinning, kCCE)
HistAxis that implements optimised numuCCE from L. Vinton. See docdb 16332. This was close to 'custC'...
const SystShifts kNoShift
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
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)
Var VarFromNuTruthVar(const NuTruthVar &stv, double _default)
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;})
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)
TFile * file_mec_wgt_numubar
const SpillCut kStandardSpillCuts
Apply this unless you're doing something special.
This module creates Common Analysis Files.
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));})
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={})
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 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;})
const Var kXSecCVWgt2018RPAFix_noDIS