From b97067c05e11759733d97d1a22175ca358d82b17 Mon Sep 17 00:00:00 2001 From: Luke Pickering Date: Wed, 16 Oct 2024 11:40:44 +0100 Subject: [PATCH 1/9] replace most pointers with vectors - anyone who said that manual memory management of individual arrays of types was better than using std::vectors was unfortunately lying to you. - Look ma' no leaks - Leaves the bool* as a bare array because the standards committee made a mistake and made std::vector a spicy specialization - removes pointless c-style casts in assignment - removes pointless struct keyword --- samplePDFDUNE/StructsDUNE.h | 135 ++++++++++++----------- samplePDFDUNE/samplePDFDUNEAtm.cpp | 16 +-- samplePDFDUNE/samplePDFDUNEAtm.h | 2 +- samplePDFDUNE/samplePDFDUNEBeamFD.cpp | 132 +++++++++++----------- samplePDFDUNE/samplePDFDUNEBeamFD.h | 2 +- samplePDFDUNE/samplePDFDUNEBeamND.cpp | 96 ++++++++-------- samplePDFDUNE/samplePDFDUNEBeamND.h | 2 +- samplePDFDUNE/samplePDFDUNEBeamNDGar.cpp | 72 ++++++------ samplePDFDUNE/samplePDFDUNEBeamNDGar.h | 2 +- 9 files changed, 233 insertions(+), 226 deletions(-) diff --git a/samplePDFDUNE/StructsDUNE.h b/samplePDFDUNE/StructsDUNE.h index fc87e48..af3d020 100644 --- a/samplePDFDUNE/StructsDUNE.h +++ b/samplePDFDUNE/StructsDUNE.h @@ -6,85 +6,88 @@ struct dunemc_base { int oscnutype; bool signal; // true if signue int nEvents; // how many MC events are there - int *Target; //Target the interaction was on + std::vector Target; //Target the interaction was on - double *rw_erec; - double *rw_erec_shifted; - double *rw_erec_had; - double *rw_erec_lep; - double *rw_yrec; - double *rw_eRecoP; - double *rw_eRecoPip; - double *rw_eRecoPim; - double *rw_eRecoPi0; - double *rw_eRecoN; - - double *rw_LepE; - double *rw_eP; - double *rw_ePip; - double *rw_ePim; - double *rw_ePi0; - double *rw_eN; - - double *rw_etru; - double *rw_mom; - double *rw_theta; - double *rw_Q2; - - double *rw_cvnnumu; - double *rw_cvnnue; - double *rw_cvnnumu_shifted; - double *rw_cvnnue_shifted; - int *rw_reco_nue; - int *rw_reco_numu; - double *rw_berpaacvwgt; - int *rw_isCC; - int *rw_nuPDGunosc; - int *rw_nuPDG; - int *rw_run; - bool *rw_isFHC; - double *rw_vtx_x; - double *rw_vtx_y; - double *rw_vtx_z; + std::vector rw_erec; + std::vector rw_erec_shifted; + std::vector rw_erec_had; + std::vector rw_erec_lep; + std::vector rw_yrec; + std::vector rw_eRecoP; + std::vector rw_eRecoPip; + std::vector rw_eRecoPim; + std::vector rw_eRecoPi0; + std::vector rw_eRecoN; + + std::vector rw_LepE; + std::vector rw_eP; + std::vector rw_ePip; + std::vector rw_ePim; + std::vector rw_ePi0; + std::vector rw_eN; + + std::vector rw_etru; + std::vector rw_mom; + std::vector rw_theta; + std::vector rw_Q2; + + std::vector rw_cvnnumu; + std::vector rw_cvnnue; + std::vector rw_cvnnumu_shifted; + std::vector rw_cvnnue_shifted; + std::vector rw_reco_nue; + std::vector rw_reco_numu; + std::vector rw_berpaacvwgt; + std::vector rw_isCC; + std::vector rw_nuPDGunosc; + std::vector rw_nuPDG; + std::vector rw_run; + std::vector rw_isFHC; + std::vector rw_vtx_x; + std::vector rw_vtx_y; + std::vector rw_vtx_z; double dummy_y; - double *rw_reco_q; - double *reco_numu; + std::vector rw_reco_q; + std::vector reco_numu; double pot_s; double norm_s; - double *beam_w; - double *flux_w; + std::vector beam_w; + std::vector flux_w; - int *mode; - int *isbound; + std::vector mode; + std::vector isbound; - double *rw_truecz; + std::vector rw_truecz; - int *nproton; ///< number of (post-FSI) primary protons - int *nneutron; ///< number of (post-FSI) primary neutrons - int *npip; ///< number of (post-FSI) primary pi+ - int *npim; ///< number of (post-FSI) primary pi- - int *npi0; ///< number of (post-FSI) primary pi0 + std::vector nproton; ///< number of (post-FSI) primary protons + std::vector nneutron; ///< number of (post-FSI) primary neutrons + std::vector npip; ///< number of (post-FSI) primary pi+ + std::vector npim; ///< number of (post-FSI) primary pi- + std::vector npi0; ///< number of (post-FSI) primary pi0 - int *ntruemuon; //number of true muons - int *ntruemuonprim; //number of true primary muons - int *nrecomuon; //number of reconstructed muons - double *nmuonsratio; //number of reco muons divided by number of true muons + std::vector ntruemuon; //number of true muons + std::vector ntruemuonprim; //number of true primary muons + std::vector nrecomuon; //number of reconstructed muons + std::vector nmuonsratio; //number of reco muons divided by number of true muons - double *rw_lep_pT; //transverse lepton momentum - double *rw_lep_pZ; //parallel lepton momentum - double *rw_reco_vtx_x; - double *rw_reco_vtx_y; - double *rw_reco_vtx_z; - double *rw_reco_rad; - double *rw_rad; + std::vector rw_lep_pT; //transverse lepton momentum + std::vector rw_lep_pZ; //parallel lepton momentum + std::vector rw_reco_vtx_x; + std::vector rw_reco_vtx_y; + std::vector rw_reco_vtx_z; + std::vector rw_reco_rad; + std::vector rw_rad; - double *rw_elep_reco; - double *rw_elep_true; + std::vector rw_elep_reco; + std::vector rw_elep_true; - int *nrecoparticles; - bool *in_fdv; + std::vector nrecoparticles; + //unfortunately std::vector is a nightmare, leave this as a pointer + bool * in_fdv; + + std::vector global_bin_number; //for using generic binning, this value holds the global bin number for each event }; // ******************************** diff --git a/samplePDFDUNE/samplePDFDUNEAtm.cpp b/samplePDFDUNE/samplePDFDUNEAtm.cpp index b53bcc6..b76c3eb 100644 --- a/samplePDFDUNE/samplePDFDUNEAtm.cpp +++ b/samplePDFDUNE/samplePDFDUNEAtm.cpp @@ -69,15 +69,15 @@ int samplePDFDUNEAtm::setupExperimentMC(int iSample) { duneobj->oscnutype = sample_oscnutype[iSample]; duneobj->signal = sample_signal[iSample]; - duneobj->mode = new int[duneobj->nEvents]; - duneobj->rw_isCC = new int[duneobj->nEvents]; - duneobj->Target = new int[duneobj->nEvents]; + duneobj->mode.resize(duneobj->nEvents); + duneobj->rw_isCC.resize(duneobj->nEvents); + duneobj->Target.resize(duneobj->nEvents); - duneobj->rw_etru = new double[duneobj->nEvents]; - duneobj->rw_truecz = new double[duneobj->nEvents]; - duneobj->flux_w = new double[duneobj->nEvents]; - duneobj->rw_erec = new double[duneobj->nEvents]; - duneobj->rw_theta = new double[duneobj->nEvents]; + duneobj->rw_etru.resize(duneobj->nEvents); + duneobj->rw_truecz.resize(duneobj->nEvents); + duneobj->flux_w.resize(duneobj->nEvents); + duneobj->rw_erec.resize(duneobj->nEvents); + duneobj->rw_theta.resize(duneobj->nEvents); for (int iEvent=0;iEventnEvents;iEvent++) { Tree->GetEntry(iEvent); diff --git a/samplePDFDUNE/samplePDFDUNEAtm.h b/samplePDFDUNE/samplePDFDUNEAtm.h index 06982f1..d4472ad 100644 --- a/samplePDFDUNE/samplePDFDUNEAtm.h +++ b/samplePDFDUNE/samplePDFDUNEAtm.h @@ -37,7 +37,7 @@ class samplePDFDUNEAtm : virtual public samplePDFFDBase inline int ReturnKinematicParameterFromString(std::string KinematicStr); inline std::string ReturnStringFromKinematicParameter(int KinematicVariable); - std::vector dunemcSamples; + std::vector dunemcSamples; bool IsELike; }; diff --git a/samplePDFDUNE/samplePDFDUNEBeamFD.cpp b/samplePDFDUNE/samplePDFDUNEBeamFD.cpp index c11e53e..70aed7c 100644 --- a/samplePDFDUNE/samplePDFDUNEBeamFD.cpp +++ b/samplePDFDUNE/samplePDFDUNEBeamFD.cpp @@ -305,85 +305,89 @@ int samplePDFDUNEBeamFD::setupExperimentMC(int iSample) { std::cout << "nevents: " << duneobj->nEvents << std::endl; // allocate memory for dunemc variables - duneobj->rw_cvnnumu = new double[duneobj->nEvents]; - duneobj->rw_cvnnue = new double[duneobj->nEvents]; - duneobj->rw_cvnnumu_shifted = new double[duneobj->nEvents]; - duneobj->rw_cvnnue_shifted = new double[duneobj->nEvents]; - duneobj->rw_etru = new double[duneobj->nEvents]; - duneobj->rw_erec = new double[duneobj->nEvents]; - duneobj->rw_erec_shifted = new double[duneobj->nEvents]; - duneobj->rw_erec_had = new double[duneobj->nEvents]; - duneobj->rw_erec_lep = new double[duneobj->nEvents]; - - duneobj->rw_eRecoP = new double[duneobj->nEvents]; - duneobj->rw_eRecoPip = new double[duneobj->nEvents]; - duneobj->rw_eRecoPim = new double[duneobj->nEvents]; - duneobj->rw_eRecoPi0 = new double[duneobj->nEvents]; - duneobj->rw_eRecoN = new double[duneobj->nEvents]; - - duneobj->rw_LepE = new double[duneobj->nEvents]; - duneobj->rw_eP = new double[duneobj->nEvents]; - duneobj->rw_ePip = new double[duneobj->nEvents]; - duneobj->rw_ePim = new double[duneobj->nEvents]; - duneobj->rw_ePi0 = new double[duneobj->nEvents]; - duneobj->rw_eN = new double[duneobj->nEvents]; - - duneobj->rw_theta = new double[duneobj->nEvents]; - duneobj->flux_w = new double[duneobj->nEvents]; - duneobj->rw_isCC = new int[duneobj->nEvents]; - duneobj->rw_nuPDGunosc = new int[duneobj->nEvents]; - duneobj->rw_nuPDG = new int[duneobj->nEvents]; - duneobj->rw_berpaacvwgt = new double[duneobj->nEvents]; - duneobj->rw_vtx_x = new double[duneobj->nEvents]; - duneobj->rw_vtx_y = new double[duneobj->nEvents]; - duneobj->rw_vtx_z = new double[duneobj->nEvents]; - - duneobj->mode = new int[duneobj->nEvents]; - duneobj->Target = new int[duneobj->nEvents]; + duneobj->rw_cvnnumu.resize(duneobj->nEvents); + duneobj->rw_cvnnue.resize(duneobj->nEvents); + duneobj->rw_cvnnumu_shifted.resize(duneobj->nEvents); + duneobj->rw_cvnnue_shifted.resize(duneobj->nEvents); + duneobj->rw_etru.resize(duneobj->nEvents); + duneobj->rw_erec.resize(duneobj->nEvents); + duneobj->rw_erec_shifted.resize(duneobj->nEvents); + duneobj->rw_erec_had.resize(duneobj->nEvents); + duneobj->rw_erec_lep.resize(duneobj->nEvents); + + duneobj->rw_eRecoP.resize(duneobj->nEvents); + duneobj->rw_eRecoPip.resize(duneobj->nEvents); + duneobj->rw_eRecoPim.resize(duneobj->nEvents); + duneobj->rw_eRecoPi0.resize(duneobj->nEvents); + duneobj->rw_eRecoN.resize(duneobj->nEvents); + + duneobj->rw_LepE.resize(duneobj->nEvents); + duneobj->rw_eP.resize(duneobj->nEvents); + duneobj->rw_ePip.resize(duneobj->nEvents); + duneobj->rw_ePim.resize(duneobj->nEvents); + duneobj->rw_ePi0.resize(duneobj->nEvents); + duneobj->rw_eN.resize(duneobj->nEvents); + + duneobj->rw_theta.resize(duneobj->nEvents); + duneobj->flux_w.resize(duneobj->nEvents); + duneobj->rw_isCC.resize(duneobj->nEvents); + duneobj->rw_nuPDGunosc.resize(duneobj->nEvents); + duneobj->rw_nuPDG.resize(duneobj->nEvents); + duneobj->rw_berpaacvwgt.resize(duneobj->nEvents); + duneobj->rw_vtx_x.resize(duneobj->nEvents); + duneobj->rw_vtx_y.resize(duneobj->nEvents); + duneobj->rw_vtx_z.resize(duneobj->nEvents); + + duneobj->global_bin_number.resize(duneobj->nEvents); + + duneobj->mode.resize(duneobj->nEvents); + duneobj->Target.resize(duneobj->nEvents); _data->GetEntry(0); //FILL DUNE STRUCT for (int i = 0; i < duneobj->nEvents; ++i) { // Loop through tree _data->GetEntry(i); - duneobj->rw_cvnnumu[i] = (double)_cvnnumu; - duneobj->rw_cvnnue[i] = (double)_cvnnue; - duneobj->rw_cvnnumu_shifted[i] = (double)_cvnnumu; - duneobj->rw_cvnnue_shifted[i] = (double)_cvnnue; + duneobj->rw_cvnnumu[i] = _cvnnumu; + duneobj->rw_cvnnue[i] = _cvnnue; + duneobj->rw_cvnnumu_shifted[i] = _cvnnumu; + duneobj->rw_cvnnue_shifted[i] = _cvnnue; if (iselike) { - duneobj->rw_erec[i] = (double)_erec_nue; - duneobj->rw_erec_shifted[i] = (double)_erec_nue; - duneobj->rw_erec_had[i] = (double)_erec_had_nue; - duneobj->rw_erec_lep[i] = (double)_erec_lep_nue; + duneobj->rw_erec[i] = _erec_nue; + duneobj->rw_erec_shifted[i] = _erec_nue; + duneobj->rw_erec_had[i] = _erec_had_nue; + duneobj->rw_erec_lep[i] = _erec_lep_nue; } else { - duneobj->rw_erec[i] = (double)_erec; - duneobj->rw_erec_shifted[i] = (double)_erec; - duneobj->rw_erec_had[i] = (double)_erec_had; - duneobj->rw_erec_lep[i] = (double)_erec_lep; + duneobj->rw_erec[i] = _erec; + duneobj->rw_erec_shifted[i] = _erec; + duneobj->rw_erec_had[i] = _erec_had; + duneobj->rw_erec_lep[i] = _erec_lep; } - duneobj->rw_eRecoP[i] = (double)_eRecoP; - duneobj->rw_eRecoPip[i] = (double)_eRecoPip; - duneobj->rw_eRecoPim[i] = (double)_eRecoPim; - duneobj->rw_eRecoPi0[i] = (double)_eRecoPi0; - duneobj->rw_eRecoN[i] = (double)_eRecoN; + duneobj->rw_eRecoP[i] = _eRecoP; + duneobj->rw_eRecoPip[i] = _eRecoPip; + duneobj->rw_eRecoPim[i] = _eRecoPim; + duneobj->rw_eRecoPi0[i] = _eRecoPi0; + duneobj->rw_eRecoN[i] = _eRecoN; - duneobj->rw_LepE[i] = (double)_LepE; - duneobj->rw_eP[i] = (double)_eP; - duneobj->rw_ePip[i] = (double)_ePip; - duneobj->rw_ePim[i] = (double)_ePim; - duneobj->rw_ePi0[i] = (double)_ePi0; - duneobj->rw_eN[i] = (double)_eN; + duneobj->rw_LepE[i] = _LepE; + duneobj->rw_eP[i] = _eP; + duneobj->rw_ePip[i] = _ePip; + duneobj->rw_ePim[i] = _ePim; + duneobj->rw_ePi0[i] = _ePi0; + duneobj->rw_eN[i] = _eN; - duneobj->rw_etru[i] = (double)_ev; - duneobj->rw_theta[i] = (double)_LepNuAngle; + duneobj->rw_etru[i] = _ev; + duneobj->rw_theta[i] = _LepNuAngle; duneobj->rw_isCC[i] = _isCC; duneobj->rw_nuPDGunosc[i] = _nuPDGunosc; duneobj->rw_nuPDG[i] = _nuPDG; - duneobj->rw_berpaacvwgt[i] = (double)_BeRPA_cvwgt; - duneobj->rw_vtx_x[i] = (double)_vtx_x; - duneobj->rw_vtx_y[i] = (double)_vtx_y; - duneobj->rw_vtx_z[i] = (double)_vtx_z; + duneobj->rw_berpaacvwgt[i] = _BeRPA_cvwgt; + duneobj->rw_vtx_x[i] = _vtx_x; + duneobj->rw_vtx_y[i] = _vtx_y; + duneobj->rw_vtx_z[i] = _vtx_z; + + duneobj->global_bin_number[i] = GetGenericBinningGlobalBinNumber(iSample, i); //Assume everything is on Argon for now.... duneobj->Target[i] = 40; diff --git a/samplePDFDUNE/samplePDFDUNEBeamFD.h b/samplePDFDUNE/samplePDFDUNEBeamFD.h index aaa5631..b8f84ae 100644 --- a/samplePDFDUNE/samplePDFDUNEBeamFD.h +++ b/samplePDFDUNE/samplePDFDUNEBeamFD.h @@ -55,7 +55,7 @@ class samplePDFDUNEBeamFD : virtual public samplePDFFDBase void applyShifts(int iSample, int iEvent); // dunemc - std::vector dunemcSamples; + std::vector dunemcSamples; double pot; diff --git a/samplePDFDUNE/samplePDFDUNEBeamND.cpp b/samplePDFDUNE/samplePDFDUNEBeamND.cpp index c9acd94..e80cc34 100644 --- a/samplePDFDUNE/samplePDFDUNEBeamND.cpp +++ b/samplePDFDUNE/samplePDFDUNEBeamND.cpp @@ -242,66 +242,66 @@ int samplePDFDUNEBeamND::setupExperimentMC(int iSample) { duneobj->oscnutype = oscnutype; duneobj->signal = signal; - duneobj->rw_yrec = new double[duneobj->nEvents]; - duneobj->rw_erec_lep = new double[duneobj->nEvents]; - duneobj->rw_erec_had = new double[duneobj->nEvents]; - duneobj->rw_etru = new double[duneobj->nEvents]; - duneobj->rw_erec = new double[duneobj->nEvents]; - duneobj->rw_erec_shifted = new double[duneobj->nEvents]; - duneobj->rw_theta = new double[duneobj->nEvents]; - duneobj->flux_w = new double[duneobj->nEvents]; - duneobj->rw_isCC = new int[duneobj->nEvents]; - duneobj->rw_reco_q = new double[duneobj->nEvents]; - duneobj->rw_nuPDGunosc = new int[duneobj->nEvents]; - duneobj->rw_nuPDG = new int[duneobj->nEvents]; - duneobj->rw_berpaacvwgt = new double[duneobj->nEvents]; - - duneobj->rw_eRecoP = new double[duneobj->nEvents]; - duneobj->rw_eRecoPip = new double[duneobj->nEvents]; - duneobj->rw_eRecoPim = new double[duneobj->nEvents]; - duneobj->rw_eRecoPi0 = new double[duneobj->nEvents]; - duneobj->rw_eRecoN = new double[duneobj->nEvents]; - - duneobj->rw_LepE = new double[duneobj->nEvents]; - duneobj->rw_eP = new double[duneobj->nEvents]; - duneobj->rw_ePip = new double[duneobj->nEvents]; - duneobj->rw_ePim = new double[duneobj->nEvents]; - duneobj->rw_ePi0 = new double[duneobj->nEvents]; - duneobj->rw_eN = new double[duneobj->nEvents]; - - duneobj->mode = new int[duneobj->nEvents]; - duneobj->Target = new int[duneobj->nEvents]; + duneobj->rw_yrec.resize(duneobj->nEvents); + duneobj->rw_erec_lep.resize(duneobj->nEvents); + duneobj->rw_erec_had.resize(duneobj->nEvents); + duneobj->rw_etru.resize(duneobj->nEvents); + duneobj->rw_erec.resize(duneobj->nEvents); + duneobj->rw_erec_shifted.resize(duneobj->nEvents); + duneobj->rw_theta.resize(duneobj->nEvents); + duneobj->flux_w.resize(duneobj->nEvents); + duneobj->rw_isCC.resize(duneobj->nEvents); + duneobj->rw_reco_q.resize(duneobj->nEvents); + duneobj->rw_nuPDGunosc.resize(duneobj->nEvents); + duneobj->rw_nuPDG.resize(duneobj->nEvents); + duneobj->rw_berpaacvwgt.resize(duneobj->nEvents); + + duneobj->rw_eRecoP.resize(duneobj->nEvents); + duneobj->rw_eRecoPip.resize(duneobj->nEvents); + duneobj->rw_eRecoPim.resize(duneobj->nEvents); + duneobj->rw_eRecoPi0.resize(duneobj->nEvents); + duneobj->rw_eRecoN.resize(duneobj->nEvents); + + duneobj->rw_LepE.resize(duneobj->nEvents); + duneobj->rw_eP.resize(duneobj->nEvents); + duneobj->rw_ePip.resize(duneobj->nEvents); + duneobj->rw_ePim.resize(duneobj->nEvents); + duneobj->rw_ePi0.resize(duneobj->nEvents); + duneobj->rw_eN.resize(duneobj->nEvents); + + duneobj->mode.resize(duneobj->nEvents); + duneobj->Target.resize(duneobj->nEvents); _data->GetEntry(0); //FILL DUNE STRUCT for (int i = 0; i < duneobj->nEvents; ++i) { // Loop through tree _data->GetEntry(i); - duneobj->rw_erec[i] = (double)_erec; - duneobj->rw_erec_shifted[i] = (double)_erec; - duneobj->rw_erec_lep[i] = (double)_erec_lep; - duneobj->rw_erec_had[i] = (double)(_erec - _erec_lep); - duneobj->rw_yrec[i] = (double)((_erec - _erec_lep)/_erec); - duneobj->rw_etru[i] = (double)_ev; // in GeV - duneobj->rw_theta[i] = (double)_LepNuAngle; + duneobj->rw_erec[i] = _erec; + duneobj->rw_erec_shifted[i] = _erec; + duneobj->rw_erec_lep[i] = _erec_lep; + duneobj->rw_erec_had[i] = (_erec - _erec_lep); + duneobj->rw_yrec[i] = ((_erec - _erec_lep)/_erec); + duneobj->rw_etru[i] = _ev; // in GeV + duneobj->rw_theta[i] = _LepNuAngle; duneobj->rw_isCC[i] = _isCC; duneobj->rw_reco_q[i] = _reco_q; duneobj->rw_nuPDGunosc[i] = _nuPDGunosc; duneobj->rw_nuPDG[i] = _nuPDG; - duneobj->rw_berpaacvwgt[i] = (double)_BeRPA_cvwgt; + duneobj->rw_berpaacvwgt[i] = _BeRPA_cvwgt; - duneobj->rw_eRecoP[i] = (double)_eRecoP; - duneobj->rw_eRecoPip[i] = (double)_eRecoPip; - duneobj->rw_eRecoPim[i] = (double)_eRecoPim; - duneobj->rw_eRecoPi0[i] = (double)_eRecoPi0; - duneobj->rw_eRecoN[i] = (double)_eRecoN; + duneobj->rw_eRecoP[i] = _eRecoP; + duneobj->rw_eRecoPip[i] = _eRecoPip; + duneobj->rw_eRecoPim[i] = _eRecoPim; + duneobj->rw_eRecoPi0[i] = _eRecoPi0; + duneobj->rw_eRecoN[i] = _eRecoN; - duneobj->rw_LepE[i] = (double)_LepE; - duneobj->rw_eP[i] = (double)_eP; - duneobj->rw_ePip[i] = (double)_ePip; - duneobj->rw_ePim[i] = (double)_ePim; - duneobj->rw_ePi0[i] = (double)_ePi0; - duneobj->rw_eN[i] = (double)_eN; + duneobj->rw_LepE[i] = _LepE; + duneobj->rw_eP[i] = _eP; + duneobj->rw_ePip[i] = _ePip; + duneobj->rw_ePim[i] = _ePim; + duneobj->rw_ePi0[i] = _ePi0; + duneobj->rw_eN[i] = _eN; //Assume everything is on Argon for now.... duneobj->Target[i] = 40; diff --git a/samplePDFDUNE/samplePDFDUNEBeamND.h b/samplePDFDUNE/samplePDFDUNEBeamND.h index fb18b14..d6c1562 100644 --- a/samplePDFDUNE/samplePDFDUNEBeamND.h +++ b/samplePDFDUNE/samplePDFDUNEBeamND.h @@ -48,7 +48,7 @@ class samplePDFDUNEBeamND : virtual public samplePDFFDBase double CalcXsecWeightFunc(int iSample, int iEvent) {return 1.;} void applyShifts(int iSample, int iEvent); - std::vector dunendmcSamples; + std::vector dunendmcSamples; TFile *_sampleFile; TTree *_data; diff --git a/samplePDFDUNE/samplePDFDUNEBeamNDGar.cpp b/samplePDFDUNE/samplePDFDUNEBeamNDGar.cpp index 1ec4ff6..b21cfe3 100644 --- a/samplePDFDUNE/samplePDFDUNEBeamNDGar.cpp +++ b/samplePDFDUNE/samplePDFDUNEBeamNDGar.cpp @@ -90,47 +90,47 @@ int samplePDFDUNEBeamNDGar::setupExperimentMC(int iSample) { std::cout << "nevents: " << duneobj->nEvents << std::endl; // allocate memory for dunendgarmc variables - duneobj->rw_yrec = new double[duneobj->nEvents]; - duneobj->rw_elep_reco = new double[duneobj->nEvents]; - duneobj->rw_etru = new double[duneobj->nEvents]; - duneobj->rw_erec = new double[duneobj->nEvents]; - duneobj->flux_w = new double[duneobj->nEvents]; - duneobj->rw_isCC = new int[duneobj->nEvents]; - duneobj->rw_nuPDGunosc = new int[duneobj->nEvents]; - duneobj->rw_nuPDG = new int[duneobj->nEvents]; - duneobj->rw_berpaacvwgt = new double[duneobj->nEvents]; - - duneobj->mode = new int[duneobj->nEvents]; - - duneobj->nproton = new int[duneobj->nEvents]; - duneobj->nneutron = new int[duneobj->nEvents]; - duneobj->npip = new int[duneobj->nEvents]; - duneobj->npim = new int[duneobj->nEvents]; - duneobj->npi0 = new int[duneobj->nEvents]; - - duneobj->nrecomuon = new int[duneobj->nEvents]; - duneobj->ntruemuon = new int[duneobj->nEvents]; - duneobj->nmuonsratio = new double[duneobj->nEvents]; - duneobj->ntruemuonprim = new int[duneobj->nEvents]; - - duneobj->nrecoparticles = new int[duneobj->nEvents]; + duneobj->rw_yrec.resize(duneobj->nEvents); + duneobj->rw_elep_reco.resize(duneobj->nEvents); + duneobj->rw_etru.resize(duneobj->nEvents); + duneobj->rw_erec.resize(duneobj->nEvents); + duneobj->flux_w.resize(duneobj->nEvents); + duneobj->rw_isCC.resize(duneobj->nEvents); + duneobj->rw_nuPDGunosc.resize(duneobj->nEvents); + duneobj->rw_nuPDG.resize(duneobj->nEvents); + duneobj->rw_berpaacvwgt.resize(duneobj->nEvents); + + duneobj->mode.resize(duneobj->nEvents); + + duneobj->nproton.resize(duneobj->nEvents); + duneobj->nneutron.resize(duneobj->nEvents); + duneobj->npip.resize(duneobj->nEvents); + duneobj->npim.resize(duneobj->nEvents); + duneobj->npi0.resize(duneobj->nEvents); + + duneobj->nrecomuon.resize(duneobj->nEvents); + duneobj->ntruemuon.resize(duneobj->nEvents); + duneobj->nmuonsratio.resize(duneobj->nEvents); + duneobj->ntruemuonprim.resize(duneobj->nEvents); + + duneobj->nrecoparticles.resize(duneobj->nEvents); duneobj->in_fdv = new bool[duneobj->nEvents]; - duneobj->rw_elep_true = new double[duneobj->nEvents]; + duneobj->rw_elep_true.resize(duneobj->nEvents); - duneobj->rw_vtx_x = new double[duneobj->nEvents]; - duneobj->rw_vtx_y = new double[duneobj->nEvents]; - duneobj->rw_vtx_z = new double[duneobj->nEvents]; - duneobj->rw_rad = new double[duneobj->nEvents]; + duneobj->rw_vtx_x.resize(duneobj->nEvents); + duneobj->rw_vtx_y.resize(duneobj->nEvents); + duneobj->rw_vtx_z.resize(duneobj->nEvents); + duneobj->rw_rad.resize(duneobj->nEvents); - duneobj->rw_lep_pT = new double[duneobj->nEvents]; - duneobj->rw_lep_pZ = new double[duneobj->nEvents]; + duneobj->rw_lep_pT.resize(duneobj->nEvents); + duneobj->rw_lep_pZ.resize(duneobj->nEvents); - duneobj->rw_reco_vtx_x = new double[duneobj->nEvents]; - duneobj->rw_reco_vtx_y = new double[duneobj->nEvents]; - duneobj->rw_reco_vtx_z = new double[duneobj->nEvents]; - duneobj->rw_reco_rad = new double[duneobj->nEvents]; + duneobj->rw_reco_vtx_x.resize(duneobj->nEvents); + duneobj->rw_reco_vtx_y.resize(duneobj->nEvents); + duneobj->rw_reco_vtx_z.resize(duneobj->nEvents); + duneobj->rw_reco_rad.resize(duneobj->nEvents); - duneobj->Target = new int[duneobj->nEvents]; + duneobj->Target.resize(duneobj->nEvents); int num_no_ixns =0; int num_no_recparticles = 0; diff --git a/samplePDFDUNE/samplePDFDUNEBeamNDGar.h b/samplePDFDUNE/samplePDFDUNEBeamNDGar.h index 26f35a4..20f3041 100644 --- a/samplePDFDUNE/samplePDFDUNEBeamNDGar.h +++ b/samplePDFDUNE/samplePDFDUNEBeamNDGar.h @@ -51,7 +51,7 @@ class samplePDFDUNEBeamNDGar : virtual public samplePDFFDBase std::string ReturnStringFromKinematicParameter(int KinematicParameter); // dunendmc - std::vector dunendgarmcSamples; + std::vector dunendgarmcSamples; TFile *_sampleFile; TTree *_data; From d20fc36f75da83aaacc8e9fcb56734e355eea541 Mon Sep 17 00:00:00 2001 From: Luke Pickering Date: Wed, 16 Oct 2024 11:41:39 +0100 Subject: [PATCH 2/9] implements simpler ReturnKinematicParameter interface. Adds generic binning - Only FD code gets generic binning currently --- samplePDFDUNE/samplePDFDUNEAtm.cpp | 41 ++---- samplePDFDUNE/samplePDFDUNEAtm.h | 15 +-- samplePDFDUNE/samplePDFDUNEBeamFD.cpp | 156 ++++------------------- samplePDFDUNE/samplePDFDUNEBeamFD.h | 19 +-- samplePDFDUNE/samplePDFDUNEBeamND.cpp | 33 ++--- samplePDFDUNE/samplePDFDUNEBeamND.h | 11 +- samplePDFDUNE/samplePDFDUNEBeamNDGar.cpp | 76 +++-------- samplePDFDUNE/samplePDFDUNEBeamNDGar.h | 12 +- 8 files changed, 89 insertions(+), 274 deletions(-) diff --git a/samplePDFDUNE/samplePDFDUNEAtm.cpp b/samplePDFDUNE/samplePDFDUNEAtm.cpp index b76c3eb..05ce059 100644 --- a/samplePDFDUNE/samplePDFDUNEAtm.cpp +++ b/samplePDFDUNE/samplePDFDUNEAtm.cpp @@ -140,49 +140,28 @@ void samplePDFDUNEAtm::setupFDMC(int iSample) { } } -double* samplePDFDUNEAtm::ReturnKinematicParameterByReference(KinematicTypes KinPar, int iSample, int iEvent) { - double* KinematicValue; +double const &samplePDFDUNEAtm::ReturnKinematicParameterByReference(int KinematicParameter, int iSample, int iEvent) { - switch (KinPar) { + switch (KinematicParameter) { case kTrueNeutrinoEnergy: - KinematicValue = &(dunemcSamples[iSample].rw_etru[iEvent]); - break; + return (dunemcSamples[iSample].rw_etru[iEvent]); case kRecoNeutrinoEnergy: - KinematicValue = &(dunemcSamples[iSample].rw_erec[iEvent]); - break; + return (dunemcSamples[iSample].rw_erec[iEvent]); case kTrueCosZ: - KinematicValue = &(dunemcSamples[iSample].rw_truecz[iEvent]); - break; + return (dunemcSamples[iSample].rw_truecz[iEvent]); case kRecoCosZ: - KinematicValue = &(dunemcSamples[iSample].rw_theta[iEvent]); - break; + return (dunemcSamples[iSample].rw_theta[iEvent]); default: - std::cerr << "Unknown KinPar:" << KinPar << std::endl; + std::cerr << "Unknown KinematicParameter:" << KinematicParameter << std::endl; throw; } - - return KinematicValue; -} - -double* samplePDFDUNEAtm::ReturnKinematicParameterByReference(double KinematicVariable, int iSample, int iEvent) { - KinematicTypes KinPar = (KinematicTypes) std::round(KinematicVariable); - return ReturnKinematicParameterByReference(KinPar,iSample,iEvent); -} - -double* samplePDFDUNEAtm::ReturnKinematicParameterByReference(std::string KinematicParameter, int iSample, int iEvent) { - KinematicTypes KinPar = static_cast(ReturnKinematicParameterFromString(KinematicParameter)); - return ReturnKinematicParameterByReference(KinPar,iSample,iEvent); -} - -double samplePDFDUNEAtm::ReturnKinematicParameter(double KinematicVariable, int iSample, int iEvent) { - return *ReturnKinematicParameterByReference(KinematicVariable, iSample, iEvent); } -double samplePDFDUNEAtm::ReturnKinematicParameter(std::string KinematicParameter, int iSample, int iEvent) { - return *ReturnKinematicParameterByReference(KinematicParameter, iSample, iEvent); +double samplePDFDUNEAtm::ReturnKinematicParameter(int KinematicParameter, int iSample, int iEvent) { + return ReturnKinematicParameterByReference(KinematicParameter, iSample, iEvent); } -std::vector samplePDFDUNEAtm::ReturnKinematicParameterBinning(std::string KinematicParameterStr) { +std::vector samplePDFDUNEAtm::ReturnKinematicParameterBinning(int KinematicParameterStr) { } int samplePDFDUNEAtm::ReturnKinematicParameterFromString(std::string KinematicParameterStr) { diff --git a/samplePDFDUNE/samplePDFDUNEAtm.h b/samplePDFDUNE/samplePDFDUNEAtm.h index d4472ad..7aeb797 100644 --- a/samplePDFDUNE/samplePDFDUNEAtm.h +++ b/samplePDFDUNE/samplePDFDUNEAtm.h @@ -28,14 +28,13 @@ class samplePDFDUNEAtm : virtual public samplePDFFDBase double CalcXsecWeightFunc(int iSample, int iEvent) {return 1.;} void applyShifts(int iSample, int iEvent) {} - double* ReturnKinematicParameterByReference(KinematicTypes KinPar, int iSample, int iEvent); - double* ReturnKinematicParameterByReference(double KinematicVariable, int iSample, int iEvent); - double* ReturnKinematicParameterByReference(std::string KinematicParameter, int iSample, int iEvent); - double ReturnKinematicParameter(double KinematicVariable, int iSample, int iEvent); - double ReturnKinematicParameter(std::string KinematicParameter, int iSample, int iEvent); - std::vector ReturnKinematicParameterBinning(std::string KinematicParameter); - inline int ReturnKinematicParameterFromString(std::string KinematicStr); - inline std::string ReturnStringFromKinematicParameter(int KinematicVariable); + double const& ReturnKinematicParameterByReference(int KinematicParameter, int iSample, int iEvent); + double ReturnKinematicParameter(int KinematicParameter, int iSample, int iEvent); + + std::vector ReturnKinematicParameterBinning(int KinematicParameter); + + int ReturnKinematicParameterFromString(std::string KinematicStr); + std::string ReturnStringFromKinematicParameter(int KinematicVariable); std::vector dunemcSamples; bool IsELike; diff --git a/samplePDFDUNE/samplePDFDUNEBeamFD.cpp b/samplePDFDUNE/samplePDFDUNEBeamFD.cpp index 70aed7c..8466fa2 100644 --- a/samplePDFDUNE/samplePDFDUNEBeamFD.cpp +++ b/samplePDFDUNE/samplePDFDUNEBeamFD.cpp @@ -488,7 +488,7 @@ TH1D* samplePDFDUNEBeamFD::get1DVarHist(KinematicTypes Var1,std::vector< std::ve } TH1D* _h1DVar; - std::vector xBinEdges = ReturnKinematicParameterBinning(ReturnStringFromKinematicParameter(Var1)); + std::vector xBinEdges = ReturnKinematicParameterBinning(Var1); _h1DVar = new TH1D("", "", xBinEdges.size()-1, xBinEdges.data()); //This should be the same as FillArray in core basically, except that @@ -540,111 +540,33 @@ TH1D* samplePDFDUNEBeamFD::get1DVarHist(KinematicTypes Var1,std::vector< std::ve return _h1DVar; } -double samplePDFDUNEBeamFD::ReturnKinematicParameter(double KinematicVariable, int iSample, int iEvent) { - KinematicTypes KinPar = (KinematicTypes) std::round(KinematicVariable); - double KinematicValue = -999; +double const& samplePDFDUNEBeamFD::ReturnKinematicParameterByReference(int KinematicParameter, int iSample, int iEvent) { - switch(KinPar){ + switch(KinematicParameter){ case kTrueNeutrinoEnergy: - KinematicValue = dunemcSamples[iSample].rw_etru[iEvent]; - break; + return dunemcSamples[iSample].rw_etru[iEvent]; case kRecoNeutrinoEnergy: - KinematicValue = dunemcSamples[iSample].rw_erec_shifted[iEvent]; - break; + return dunemcSamples[iSample].rw_erec_shifted[iEvent]; case kTrueXPos: - KinematicValue = dunemcSamples[iSample].rw_vtx_x[iEvent]; - break; + return dunemcSamples[iSample].rw_vtx_x[iEvent]; case kTrueYPos: - KinematicValue = dunemcSamples[iSample].rw_vtx_y[iEvent]; - break; + return dunemcSamples[iSample].rw_vtx_y[iEvent]; case kTrueZPos: - KinematicValue = dunemcSamples[iSample].rw_vtx_z[iEvent]; - break; + return dunemcSamples[iSample].rw_vtx_z[iEvent]; case kCVNNumu: - KinematicValue = dunemcSamples[iSample].rw_cvnnumu_shifted[iEvent]; - break; + return dunemcSamples[iSample].rw_cvnnumu_shifted[iEvent]; case kCVNNue: - KinematicValue = dunemcSamples[iSample].rw_cvnnue_shifted[iEvent]; - break; - case kM3Mode: - KinematicValue = dunemcSamples[iSample].mode[iEvent]; - break; + return dunemcSamples[iSample].rw_cvnnue_shifted[iEvent]; + case kGlobalBinNumber: + return dunemcSamples[iSample].global_bin_number[iEvent]; default: - MACH3LOG_ERROR("Did not recognise Kinematic Parameter type"); - MACH3LOG_ERROR("Was given a Kinematic Variable of {}", KinematicVariable); - throw MaCh3Exception(__FILE__, __LINE__); + std::cout << "[ERROR]: " << __FILE__ << ":" << __LINE__ << " Did not recognise Kinematic Parameter type..." << std::endl; + throw; + } } - - return KinematicValue; -} - -double samplePDFDUNEBeamFD::ReturnKinematicParameter(std::string KinematicParameter, int iSample, int iEvent) { - KinematicTypes KinPar = static_cast(ReturnKinematicParameterFromString(KinematicParameter)); - double KinematicValue = -999; - - switch(KinPar){ - case kTrueNeutrinoEnergy: - KinematicValue = dunemcSamples[iSample].rw_etru[iEvent]; - break; - case kRecoNeutrinoEnergy: - KinematicValue = dunemcSamples[iSample].rw_erec_shifted[iEvent]; - break; - case kTrueXPos: - KinematicValue = dunemcSamples[iSample].rw_vtx_x[iEvent]; - break; - case kTrueYPos: - KinematicValue = dunemcSamples[iSample].rw_vtx_y[iEvent]; - break; - case kTrueZPos: - KinematicValue = dunemcSamples[iSample].rw_vtx_z[iEvent]; - break; - case kCVNNumu: - KinematicValue = dunemcSamples[iSample].rw_cvnnumu_shifted[iEvent]; - break; - case kCVNNue: - KinematicValue = dunemcSamples[iSample].rw_cvnnue_shifted[iEvent]; - break; - default: - std::cout << "[ERROR]: " << __FILE__ << ":" << __LINE__ << " Did not recognise Kinematic Parameter type..." << std::endl; - throw; - } - - return KinematicValue; -} - -const double* samplePDFDUNEBeamFD::ReturnKinematicParameterByReference(std::string KinematicParameter, int iSample, int iEvent) { - KinematicTypes KinPar = static_cast(ReturnKinematicParameterFromString(KinematicParameter)); - double* KinematicValue; - - switch(KinPar){ - case kTrueNeutrinoEnergy: - KinematicValue = &dunemcSamples[iSample].rw_etru[iEvent]; - break; - case kRecoNeutrinoEnergy: - KinematicValue = &dunemcSamples[iSample].rw_erec_shifted[iEvent]; - break; - case kTrueXPos: - KinematicValue = &dunemcSamples[iSample].rw_vtx_x[iEvent]; - break; - case kTrueYPos: - KinematicValue = &dunemcSamples[iSample].rw_vtx_y[iEvent]; - break; - case kTrueZPos: - KinematicValue = &dunemcSamples[iSample].rw_vtx_z[iEvent]; - break; - case kCVNNumu: - KinematicValue = &dunemcSamples[iSample].rw_cvnnumu_shifted[iEvent]; - break; - case kCVNNue: - KinematicValue = &dunemcSamples[iSample].rw_cvnnue_shifted[iEvent]; - break; - default: - std::cout << "[ERROR]: " << __FILE__ << ":" << __LINE__ << " Did not recognise Kinematic Parameter type..." << std::endl; - throw; - } - - return KinematicValue; +double samplePDFDUNEBeamFD::ReturnKinematicParameter(int KinematicParameter, int iSample, int iEvent) { + return ReturnKinematicParameterByReference(KinematicParameter, iSample, iEvent); } int samplePDFDUNEBeamFD::ReturnKinematicParameterFromString(std::string KinematicParameterStr){ @@ -657,41 +579,9 @@ int samplePDFDUNEBeamFD::ReturnKinematicParameterFromString(std::string Kinemati if (KinematicParameterStr.find("CVNNumu") != std::string::npos) {return kCVNNumu;} if (KinematicParameterStr.find("CVNNue") != std::string::npos) {return kCVNNue;} if (KinematicParameterStr.find("M3Mode") != std::string::npos) {return kM3Mode;} + if (KinematicParameterStr.find("global_bin_number") != std::string::npos) {return kGlobalBinNumber;} -} - -const double* samplePDFDUNEBeamFD::ReturnKinematicParameterByReference(double KinematicVariable, int iSample, int iEvent) { - KinematicTypes KinPar = (KinematicTypes) std::round(KinematicVariable); - double* KinematicValue; - switch(KinPar){ - case kTrueNeutrinoEnergy: - KinematicValue = &dunemcSamples[iSample].rw_etru[iEvent]; - break; - case kRecoNeutrinoEnergy: - KinematicValue = &dunemcSamples[iSample].rw_erec_shifted[iEvent]; - break; - case kTrueXPos: - KinematicValue = &dunemcSamples[iSample].rw_vtx_x[iEvent]; - break; - case kTrueYPos: - KinematicValue = &dunemcSamples[iSample].rw_vtx_y[iEvent]; - break; - case kTrueZPos: - KinematicValue = &dunemcSamples[iSample].rw_vtx_z[iEvent]; - break; - case kCVNNumu: - KinematicValue = &dunemcSamples[iSample].rw_cvnnumu_shifted[iEvent]; - break; - case kCVNNue: - KinematicValue = &dunemcSamples[iSample].rw_cvnnue_shifted[iEvent]; - break; - default: - std::cout << "[ERROR]: " << __FILE__ << ":" << __LINE__ << " Did not recognise Kinematic Parameter type..." << std::endl; - throw; - } - - return KinematicValue; } inline std::string samplePDFDUNEBeamFD::ReturnStringFromKinematicParameter(int KinematicParameter) { @@ -720,8 +610,11 @@ inline std::string samplePDFDUNEBeamFD::ReturnStringFromKinematicParameter(int K KinematicString = "CVNNue"; break; case kM3Mode: - KinematicString = "M3Mode"; - break; + KinematicString = "M3Mode"; + break; + case kGlobalBinNumber: + KinematicString = "global_bin_number"; + break; default: break; } @@ -826,9 +719,8 @@ void samplePDFDUNEBeamFD::applyShifts(int iSample, int iEvent) { */ } -std::vector samplePDFDUNEBeamFD::ReturnKinematicParameterBinning(std::string KinematicParameterStr) { +std::vector samplePDFDUNEBeamFD::ReturnKinematicParameterBinning(int KinematicParameter) { std::vector binningVector; - KinematicTypes KinematicParameter = static_cast(ReturnKinematicParameterFromString(KinematicParameterStr)); int nBins = 0; double bin_width = 0; diff --git a/samplePDFDUNE/samplePDFDUNEBeamFD.h b/samplePDFDUNE/samplePDFDUNEBeamFD.h index b8f84ae..d59de6a 100644 --- a/samplePDFDUNE/samplePDFDUNEBeamFD.h +++ b/samplePDFDUNE/samplePDFDUNEBeamFD.h @@ -25,7 +25,7 @@ class samplePDFDUNEBeamFD : virtual public samplePDFFDBase samplePDFDUNEBeamFD(std::string mc_version, covarianceXsec* xsec_cov); ~samplePDFDUNEBeamFD(); - enum KinematicTypes {kTrueNeutrinoEnergy,kRecoNeutrinoEnergy,kTrueXPos,kTrueYPos,kTrueZPos,kCVNNumu,kCVNNue,kM3Mode,kOscChannel}; + enum KinematicTypes {kTrueNeutrinoEnergy,kRecoNeutrinoEnergy,kTrueXPos,kTrueYPos,kTrueZPos,kCVNNumu,kCVNNue,kM3Mode,kGlobalBinNumber,kOscChannel}; //More robust getters to make plots in different variables, mode, osc channel, systematic weighting and with bin range TH1D* get1DVarHist(KinematicTypes Var1, int fModeToFill=-1, int fSampleToFill=-1, int WeightStyle=0, TAxis* Axis=0); @@ -41,15 +41,16 @@ class samplePDFDUNEBeamFD : virtual public samplePDFFDBase /// @todo extract this completely to core ///@brief Setup our spline file, this calls InitialseSplineObject() under the hood void SetupSplines(); + + double const &ReturnKinematicParameterByReference(int KinematicParameter, + int iSample, int iEvent); + double ReturnKinematicParameter(int KinematicParameter, int iSample, + int iEvent); + + std::vector ReturnKinematicParameterBinning(int KinematicParameter); - double ReturnKinematicParameter (double KinematicVariable, int iSample, int iEvent); - double ReturnKinematicParameter(std::string KinematicParameter, int iSample, int iEvent); - std::vector ReturnKinematicParameterBinning(std::string KinematicParameter); - const double* ReturnKinematicParameterByReference(std::string KinematicParameter, int iSample, int iEvent); - const double* ReturnKinematicParameterByReference(double KinematicVariable, int iSample, int iEvent); - - inline int ReturnKinematicParameterFromString(std::string KinematicParameterStr); - inline std::string ReturnStringFromKinematicParameter(int KinematicParameterStr); + int ReturnKinematicParameterFromString(std::string KinematicParameterStr); + std::string ReturnStringFromKinematicParameter(int KinematicParameterStr); //DB functions which could be initialised to do something which is non-trivial double CalcXsecWeightFunc(int iSample, int iEvent) {return 1.;} void applyShifts(int iSample, int iEvent); diff --git a/samplePDFDUNE/samplePDFDUNEBeamND.cpp b/samplePDFDUNE/samplePDFDUNEBeamND.cpp index e80cc34..b37418b 100644 --- a/samplePDFDUNE/samplePDFDUNEBeamND.cpp +++ b/samplePDFDUNE/samplePDFDUNEBeamND.cpp @@ -316,40 +316,22 @@ int samplePDFDUNEBeamND::setupExperimentMC(int iSample) { return duneobj->nEvents; } -double* samplePDFDUNEBeamND::ReturnKinematicParameterByReference(KinematicTypes KinPar, int iSample, int iEvent) { - double* KinematicValue; +double const& samplePDFDUNEBeamND::ReturnKinematicParameterByReference(int KinematicParameter, int iSample, int iEvent) { - switch(KinPar){ + switch(KinematicParameter){ case kTrueNeutrinoEnergy: - KinematicValue = &dunendmcSamples[iSample].rw_etru[iEvent]; - break; + return dunendmcSamples[iSample].rw_etru[iEvent]; case kRecoQ: - KinematicValue = &dunendmcSamples[iSample].rw_reco_q[iEvent]; - break; + return dunendmcSamples[iSample].rw_reco_q[iEvent]; default: std::cout << "[ERROR]: " << __FILE__ << ":" << __LINE__ << " Did not recognise Kinematic Parameter type..." << std::endl; throw; } - return KinematicValue; } -double* samplePDFDUNEBeamND::ReturnKinematicParameterByReference(double KinematicVariable, int iSample, int iEvent) { - KinematicTypes KinPar = (KinematicTypes) std::round(KinematicVariable); - return ReturnKinematicParameterByReference(KinPar,iSample,iEvent); -} - -double* samplePDFDUNEBeamND::ReturnKinematicParameterByReference(std::string KinematicParameter, int iSample, int iEvent) { - KinematicTypes KinPar = static_cast(ReturnKinematicParameterFromString(KinematicParameter)); - return ReturnKinematicParameterByReference(KinPar,iSample,iEvent); -} - -double samplePDFDUNEBeamND::ReturnKinematicParameter(double KinematicVariable, int iSample, int iEvent) { - return *ReturnKinematicParameterByReference(KinematicVariable, iSample, iEvent); -} - -double samplePDFDUNEBeamND::ReturnKinematicParameter(std::string KinematicParameter, int iSample, int iEvent) { - return *ReturnKinematicParameterByReference(KinematicParameter, iSample, iEvent); +double samplePDFDUNEBeamND::ReturnKinematicParameter(int KinematicParameter, int iSample, int iEvent) { + return ReturnKinematicParameterByReference(KinematicParameter, iSample, iEvent); } void samplePDFDUNEBeamND::setupFDMC(int iSample) { @@ -457,7 +439,8 @@ void samplePDFDUNEBeamND::applyShifts(int iSample, int iEvent) { } -std::vector samplePDFDUNEBeamND::ReturnKinematicParameterBinning(std::string KinematicParameterStr) + +std::vector samplePDFDUNEBeamND::ReturnKinematicParameterBinning(int KinematicParameter) { std::cout << "ReturnKinematicVarBinning" << std::endl; std::vector binningVector; diff --git a/samplePDFDUNE/samplePDFDUNEBeamND.h b/samplePDFDUNE/samplePDFDUNEBeamND.h index d6c1562..815c731 100644 --- a/samplePDFDUNE/samplePDFDUNEBeamND.h +++ b/samplePDFDUNE/samplePDFDUNEBeamND.h @@ -35,12 +35,11 @@ class samplePDFDUNEBeamND : virtual public samplePDFFDBase void SetupWeightPointers(); void SetupSplines(); - double* ReturnKinematicParameterByReference(KinematicTypes KinPar, int iSample, int iEvent); - double* ReturnKinematicParameterByReference(double KinematicVariable, int iSample, int iEvent); - double* ReturnKinematicParameterByReference(std::string KinematicParameter, int iSample, int iEvent); - double ReturnKinematicParameter(double KinematicVariable, int iSample, int iEvent); - double ReturnKinematicParameter(std::string KinematicParameter, int iSample, int iEvent); - std::vector ReturnKinematicParameterBinning(std::string KinematicParameter); + double const &ReturnKinematicParameterByReference(int KinematicParameter, int iSample, int iEvent); + double ReturnKinematicParameter(int KinematicParameter, int iSample, int iEvent); + + std::vector ReturnKinematicParameterBinning(int KinematicParameter); + int ReturnKinematicParameterFromString(std::string KinematicParameterStr); std::string ReturnStringFromKinematicParameter(int KinematicParameter); diff --git a/samplePDFDUNE/samplePDFDUNEBeamNDGar.cpp b/samplePDFDUNE/samplePDFDUNEBeamNDGar.cpp index b21cfe3..ae8342d 100644 --- a/samplePDFDUNE/samplePDFDUNEBeamNDGar.cpp +++ b/samplePDFDUNE/samplePDFDUNEBeamNDGar.cpp @@ -254,79 +254,48 @@ int samplePDFDUNEBeamNDGar::setupExperimentMC(int iSample) { return duneobj->nEvents; } -double* samplePDFDUNEBeamNDGar::ReturnKinematicParameterByReference(KinematicTypes KinematicParameter, int iSample, int iEvent) { - double* KinematicValue; +double const& samplePDFDUNEBeamNDGar::ReturnKinematicParameterByReference(int KinematicParameter, int iSample, int iEvent) { switch(KinematicParameter) { case kTrueNeutrinoEnergy: - KinematicValue = &dunendgarmcSamples[iSample].rw_etru[iEvent]; - break; + return dunendgarmcSamples[iSample].rw_etru[iEvent]; case kRecoNeutrinoEnergy: - KinematicValue = &dunendgarmcSamples[iSample].rw_erec[iEvent]; - break; + return dunendgarmcSamples[iSample].rw_erec[iEvent]; case kTrueXPos: - KinematicValue = &dunendgarmcSamples[iSample].rw_vtx_x[iEvent]; - break; + return dunendgarmcSamples[iSample].rw_vtx_x[iEvent]; case kTrueYPos: - KinematicValue = &dunendgarmcSamples[iSample].rw_vtx_y[iEvent]; - break; + return dunendgarmcSamples[iSample].rw_vtx_y[iEvent]; case kTrueZPos: - KinematicValue = &dunendgarmcSamples[iSample].rw_vtx_z[iEvent]; - break; + return dunendgarmcSamples[iSample].rw_vtx_z[iEvent]; case kTrueRad: - KinematicValue = &dunendgarmcSamples[iSample].rw_rad[iEvent]; - break; + return dunendgarmcSamples[iSample].rw_rad[iEvent]; case kNMuonsRecoOverTruth: - KinematicValue = &dunendgarmcSamples[iSample].nmuonsratio[iEvent]; - break; + return dunendgarmcSamples[iSample].nmuonsratio[iEvent]; case kRecoLepEnergy: - KinematicValue = &dunendgarmcSamples[iSample].rw_elep_reco[iEvent]; - break; + return dunendgarmcSamples[iSample].rw_elep_reco[iEvent]; case kTrueLepEnergy: - KinematicValue = &dunendgarmcSamples[iSample].rw_elep_true[iEvent]; - break; + return dunendgarmcSamples[iSample].rw_elep_true[iEvent]; case kRecoXPos: - KinematicValue = &dunendgarmcSamples[iSample].rw_reco_vtx_x[iEvent]; - break; + return dunendgarmcSamples[iSample].rw_reco_vtx_x[iEvent]; case kRecoYPos: - KinematicValue = &dunendgarmcSamples[iSample].rw_reco_vtx_y[iEvent]; - break; + return dunendgarmcSamples[iSample].rw_reco_vtx_y[iEvent]; case kRecoZPos: - KinematicValue = &dunendgarmcSamples[iSample].rw_reco_vtx_z[iEvent]; - break; + return dunendgarmcSamples[iSample].rw_reco_vtx_z[iEvent]; case kRecoRad: - KinematicValue = &dunendgarmcSamples[iSample].rw_reco_rad[iEvent]; - break; + return dunendgarmcSamples[iSample].rw_reco_rad[iEvent]; case kLepPT: - KinematicValue = &dunendgarmcSamples[iSample].rw_lep_pT[iEvent]; - break; + return dunendgarmcSamples[iSample].rw_lep_pT[iEvent]; case kLepPZ: - KinematicValue = &dunendgarmcSamples[iSample].rw_lep_pZ[iEvent]; - break; + return dunendgarmcSamples[iSample].rw_lep_pZ[iEvent]; default: std::cout << "[ERROR]: " << __FILE__ << ":" << __LINE__ << " Did not recognise Kinematic Parameter type..." << std::endl; throw; } - return KinematicValue; } -double* samplePDFDUNEBeamNDGar::ReturnKinematicParameterByReference(double KinematicVariable, int iSample, int iEvent) { - KinematicTypes KinPar = (KinematicTypes) std::round(KinematicVariable); - return ReturnKinematicParameterByReference(KinPar,iSample,iEvent); -} - -double* samplePDFDUNEBeamNDGar::ReturnKinematicParameterByReference(std::string KinematicParameter, int iSample, int iEvent) { - KinematicTypes KinPar = static_cast(ReturnKinematicParameterFromString(KinematicParameter)); - return ReturnKinematicParameterByReference(KinPar,iSample,iEvent); -} - -double samplePDFDUNEBeamNDGar::ReturnKinematicParameter(double KinematicVariable, int iSample, int iEvent) { - return *ReturnKinematicParameterByReference(KinematicVariable, iSample, iEvent); -} - -double samplePDFDUNEBeamNDGar::ReturnKinematicParameter(std::string KinematicParameter, int iSample, int iEvent) { - return *ReturnKinematicParameterByReference(KinematicParameter, iSample, iEvent); +double samplePDFDUNEBeamNDGar::ReturnKinematicParameter(int KinematicParameter, int iSample, int iEvent) { + return ReturnKinematicParameterByReference(KinematicParameter,iSample,iEvent); } void samplePDFDUNEBeamNDGar::setupFDMC(int iSample) { @@ -370,14 +339,9 @@ void samplePDFDUNEBeamNDGar::setupFDMC(int iSample) { } -std::vector samplePDFDUNEBeamNDGar::ReturnKinematicParameterBinning(std::string KinematicParameterStr) { - KinematicTypes KinPar = static_cast(ReturnKinematicParameterFromString(KinematicParameterStr)); - return ReturnKinematicParameterBinning(KinPar); -} - -std::vector samplePDFDUNEBeamNDGar::ReturnKinematicParameterBinning(KinematicTypes KinPar) { +std::vector samplePDFDUNEBeamNDGar::ReturnKinematicParameterBinning(int KinematicParameter) { std::vector binningVector; - switch(KinPar){ + switch(KinematicParameter){ case kTrueNeutrinoEnergy: for(double ibins =0; ibins<10*10; ibins++){ double binval = ibins/10; diff --git a/samplePDFDUNE/samplePDFDUNEBeamNDGar.h b/samplePDFDUNE/samplePDFDUNEBeamNDGar.h index 20f3041..ebefb0c 100644 --- a/samplePDFDUNE/samplePDFDUNEBeamNDGar.h +++ b/samplePDFDUNE/samplePDFDUNEBeamNDGar.h @@ -40,13 +40,11 @@ class samplePDFDUNEBeamNDGar : virtual public samplePDFFDBase double CalcXsecWeightFunc(int iSample, int iEvent) {return 1.;} void applyShifts(int iSample, int iEvent) {} - double* ReturnKinematicParameterByReference(KinematicTypes KinPar, int iSample, int iEvent); - double* ReturnKinematicParameterByReference(double KinematicVariable, int iSample, int iEvent); - double* ReturnKinematicParameterByReference(std::string KinematicParameter, int iSample, int iEvent); - double ReturnKinematicParameter(double KinematicVariable, int iSample, int iEvent); - double ReturnKinematicParameter(std::string KinematicParameter, int iSample, int iEvent); - std::vector ReturnKinematicParameterBinning(std::string KinematicParameter); - std::vector ReturnKinematicParameterBinning(KinematicTypes KinematicParameter); + double const& ReturnKinematicParameterByReference(int KinematicParameter, int iSample, int iEvent); + double ReturnKinematicParameter(int KinematicParameter, int iSample, int iEvent); + + std::vector ReturnKinematicParameterBinning(int KinematicParameter); + int ReturnKinematicParameterFromString(std::string KinematicParameterStr); std::string ReturnStringFromKinematicParameter(int KinematicParameter); From a3cbee59f5f06b579af20c2cc09058471387d369 Mon Sep 17 00:00:00 2001 From: Luke Pickering Date: Wed, 16 Oct 2024 22:47:04 +0100 Subject: [PATCH 3/9] remove double 'check' of ROOT CXX standard and update default to 17 --- CMakeLists.txt | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index fb8ee43..1d7f8d5 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -45,10 +45,6 @@ if(NOT TARGET ROOT::ROOT) cmessage(FATAL_ERROR "MaCh3 Expected dependency target: ROOT::ROOT") endif() -if(DEFINED ROOT_CXX_STANDARD AND ROOT_CXX_STANDARD GREATER CMAKE_CXX_STANDARD) - set(CMAKE_CXX_STANDARD ${ROOT_CXX_STANDARD}) -endif() - ############################ DUNEAnaObj #################################### find_package(duneanaobj) @@ -89,7 +85,7 @@ endif() ############################ C++ Compiler #################################### if (NOT DEFINED CMAKE_CXX_STANDARD OR "${CMAKE_CXX_STANDARD} " STREQUAL " ") - SET(CMAKE_CXX_STANDARD 11) + SET(CMAKE_CXX_STANDARD 17) endif() if(DEFINED ROOT_CXX_STANDARD AND ROOT_CXX_STANDARD GREATER CMAKE_CXX_STANDARD) From 32e1633121f62d3b6b739f0116206c6e4efd9bdb Mon Sep 17 00:00:00 2001 From: Luke Pickering Date: Wed, 16 Oct 2024 22:47:34 +0100 Subject: [PATCH 4/9] adds q0/q3 to structs dune --- samplePDFDUNE/StructsDUNE.h | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/samplePDFDUNE/StructsDUNE.h b/samplePDFDUNE/StructsDUNE.h index af3d020..d406f4a 100644 --- a/samplePDFDUNE/StructsDUNE.h +++ b/samplePDFDUNE/StructsDUNE.h @@ -19,6 +19,9 @@ struct dunemc_base { std::vector rw_eRecoPi0; std::vector rw_eRecoN; + std::vector true_q0; + std::vector true_q3; + std::vector rw_LepE; std::vector rw_eP; std::vector rw_ePip; @@ -42,7 +45,7 @@ struct dunemc_base { std::vector rw_nuPDGunosc; std::vector rw_nuPDG; std::vector rw_run; - std::vector rw_isFHC; + bool * rw_isFHC; std::vector rw_vtx_x; std::vector rw_vtx_y; std::vector rw_vtx_z; From bbad9d441b824702d8f5033e38c92815507b4752 Mon Sep 17 00:00:00 2001 From: Luke Pickering Date: Wed, 16 Oct 2024 22:49:58 +0100 Subject: [PATCH 5/9] use std::filesystem::path for simple/mtuple locations --- samplePDFDUNE/samplePDFDUNEBeamFD.cpp | 27 +++++++++++------------- samplePDFDUNE/samplePDFDUNEBeamND.cpp | 7 +++--- samplePDFDUNE/samplePDFDUNEBeamNDGar.cpp | 7 +++--- splines/splinesDUNE.cpp | 2 +- splines/splinesDUNE.h | 2 +- 5 files changed, 20 insertions(+), 25 deletions(-) diff --git a/samplePDFDUNE/samplePDFDUNEBeamFD.cpp b/samplePDFDUNE/samplePDFDUNEBeamFD.cpp index 8466fa2..8ee3859 100644 --- a/samplePDFDUNE/samplePDFDUNEBeamFD.cpp +++ b/samplePDFDUNE/samplePDFDUNEBeamFD.cpp @@ -168,15 +168,13 @@ void samplePDFDUNEBeamFD::Init() { void samplePDFDUNEBeamFD::SetupSplines() { ///@todo move all of the spline setup into core - if(XsecCov->GetNumParamsFromDetID(SampleDetID, kSpline) > 0){ - MACH3LOG_INFO("Found {} splines for this sample so I will create a spline object", XsecCov->GetNumParamsFromDetID(SampleDetID, kSpline)); - splinesDUNE* DUNESplines = new splinesDUNE(XsecCov); - splineFile = (splineFDBase*)DUNESplines; - InitialiseSplineObject(); - } - else{ - MACH3LOG_INFO("Found {} splines for this sample so I will not load or evaluate splines", XsecCov->GetNumParamsFromDetID(SampleDetID, kSpline)); - splineFile = nullptr; + if(spline_files.size() && (XsecCov->GetNumParamsFromDetID(SampleDetID, kSpline) > 0)){ + MACH3LOG_INFO("Found {} splines for this sample so I will create a spline object", XsecCov->GetNumParamsFromDetID(SampleDetID, kSpline)); + splineFile = new splinesDUNE(XsecCov); + InitialiseSplineObject(); + } else { + MACH3LOG_INFO("Found {} splines for this sample so I will not load or evaluate splines", XsecCov->GetNumParamsFromDetID(SampleDetID, kSpline)); + splineFile = nullptr; } return; @@ -199,24 +197,23 @@ void samplePDFDUNEBeamFD::SetupWeightPointers() { int samplePDFDUNEBeamFD::setupExperimentMC(int iSample) { - const char *sampleFile = (mtuple_files[iSample]).c_str(); - dunemc_base *duneobj = &(dunemcSamples[iSample]); + auto &duneobj = dunemcSamples[iSample]; int nutype = sample_nutype[iSample]; int oscnutype = sample_oscnutype[iSample]; bool signal = sample_signal[iSample]; std::cout << "-------------------------------------------------------------------" << std::endl; - MACH3LOG_INFO("input file: {}", sampleFile); + MACH3LOG_INFO("input file: {}", mtuple_files[iSample].native()); - _sampleFile = new TFile(sampleFile, "READ"); + _sampleFile = new TFile(mtuple_files[iSample].c_str(), "READ"); _data = (TTree*)_sampleFile->Get("caf"); if(_data){ - MACH3LOG_INFO("Found \"caf\" tree in {}", sampleFile); + MACH3LOG_INFO("Found \"caf\" tree in {}", mtuple_files[iSample].native()); MACH3LOG_INFO("With number of entries: {}", _data->GetEntries()); } else{ - MACH3LOG_ERROR("Could not find \"caf\" tree in {}", sampleFile); + MACH3LOG_ERROR("Could not find \"caf\" tree in {}", mtuple_files[iSample].native()); throw MaCh3Exception(__FILE__, __LINE__); } diff --git a/samplePDFDUNE/samplePDFDUNEBeamND.cpp b/samplePDFDUNE/samplePDFDUNEBeamND.cpp index b37418b..afb1739 100644 --- a/samplePDFDUNE/samplePDFDUNEBeamND.cpp +++ b/samplePDFDUNE/samplePDFDUNEBeamND.cpp @@ -169,20 +169,19 @@ void samplePDFDUNEBeamND::SetupWeightPointers() { } int samplePDFDUNEBeamND::setupExperimentMC(int iSample) { - const char *sampleFile = (mtupleprefix+mtuple_files[iSample]+mtuplesuffix).c_str(); dunemc_base *duneobj = &(dunendmcSamples[iSample]); int nutype = sample_nutype[iSample]; int oscnutype = sample_oscnutype[iSample]; bool signal = sample_signal[iSample]; std::cout << "-------------------------------------------------------------------" << std::endl; - std::cout << "input file: " << sampleFile << std::endl; + std::cout << "input file: " << mtuple_files[iSample] << std::endl; - _sampleFile = new TFile(sampleFile, "READ"); + _sampleFile = new TFile(mtuple_files[iSample].c_str(), "READ"); _data = (TTree*)_sampleFile->Get("caf"); if(_data){ - std::cout << "Found mtuple tree is " << sampleFile << std::endl; + std::cout << "Found mtuple tree is " << mtuple_files[iSample] << std::endl; std::cout << "N of entries: " << _data->GetEntries() << std::endl; } diff --git a/samplePDFDUNE/samplePDFDUNEBeamNDGar.cpp b/samplePDFDUNE/samplePDFDUNEBeamNDGar.cpp index ae8342d..0a6c725 100644 --- a/samplePDFDUNE/samplePDFDUNEBeamNDGar.cpp +++ b/samplePDFDUNE/samplePDFDUNEBeamNDGar.cpp @@ -59,19 +59,18 @@ void samplePDFDUNEBeamNDGar::SetupWeightPointers() { } int samplePDFDUNEBeamNDGar::setupExperimentMC(int iSample) { - const char *sampleFile = (mtupleprefix+mtuple_files[iSample]+mtuplesuffix).c_str(); dunemc_base *duneobj = &(dunendgarmcSamples[iSample]); int nutype = sample_nutype[iSample]; int oscnutype = sample_oscnutype[iSample]; bool signal = sample_signal[iSample]; std::cout << "-------------------------------------------------------------------" << std::endl; - std::cout << "input file: " << sampleFile << std::endl; + std::cout << "input file: " << mtuple_files[iSample] << std::endl; - _sampleFile = new TFile(sampleFile, "READ"); + _sampleFile = new TFile(mtuple_files[iSample].c_str(), "READ"); _data = (TTree*)_sampleFile->Get("cafTree"); if(_data){ - std::cout << "Found mtuple tree is " << sampleFile << std::endl; + std::cout << "Found mtuple tree is " << mtuple_files[iSample] << std::endl; std::cout << "N of entries: " << _data->GetEntries() << std::endl; } diff --git a/splines/splinesDUNE.cpp b/splines/splinesDUNE.cpp index 2421ab3..82f8dd0 100644 --- a/splines/splinesDUNE.cpp +++ b/splines/splinesDUNE.cpp @@ -13,7 +13,7 @@ splinesDUNE::~splinesDUNE(){ } //**************************************** -void splinesDUNE::FillSampleArray(std::string SampleName, std::vector OscChanFileNames) +void splinesDUNE::FillSampleArray(std::string SampleName, std::vector OscChanFileNames) // Performs two jobs // 1. Fills indexing/each sample // 2. Creates the big spline vector diff --git a/splines/splinesDUNE.h b/splines/splinesDUNE.h index 45e569b..2f3790e 100644 --- a/splines/splinesDUNE.h +++ b/splines/splinesDUNE.h @@ -9,7 +9,7 @@ class splinesDUNE : virtual public splineFDBase splinesDUNE(covarianceXsec* xsec_cov); virtual ~splinesDUNE(); void SetupSplines(int BinningOpt); - virtual void FillSampleArray(std::string SampleName, std::vector OscChanFileNames) override; + virtual void FillSampleArray(std::string SampleName, std::vector OscChanFileNames) override; virtual std::vector< std::vector > StripDuplicatedModes(std::vector< std::vector > InputVector) override; virtual std::vector< std::vector > GetEventSplines(std::string SampleName, int iOscChan, int EventMode, double Var1Val, double Var2Val, double Var3Val) override; }; From b65f30c8f678c1bd552c0d8d2ea68538d9e6e309 Mon Sep 17 00:00:00 2001 From: Luke Pickering Date: Wed, 16 Oct 2024 22:50:41 +0100 Subject: [PATCH 6/9] adds some additional kinematicparameters to FD samplePDF --- samplePDFDUNE/samplePDFDUNEBeamFD.cpp | 396 ++++++++++++++++---------- samplePDFDUNE/samplePDFDUNEBeamFD.h | 18 +- 2 files changed, 270 insertions(+), 144 deletions(-) diff --git a/samplePDFDUNE/samplePDFDUNEBeamFD.cpp b/samplePDFDUNE/samplePDFDUNEBeamFD.cpp index 8ee3859..3f1d3b8 100644 --- a/samplePDFDUNE/samplePDFDUNEBeamFD.cpp +++ b/samplePDFDUNE/samplePDFDUNEBeamFD.cpp @@ -53,10 +53,7 @@ void samplePDFDUNEBeamFD::Init() { cvn_nue_fd_pos = -999; // create dunemc storage - for (int i=0;i FDDetectorSystPointersMap; @@ -276,7 +273,20 @@ int samplePDFDUNEBeamFD::setupExperimentMC(int iSample) { _data->SetBranchStatus("vtx_y", 1); _data->SetBranchAddress("vtx_y", &_vtx_y); _data->SetBranchStatus("vtx_z", 1); - _data->SetBranchAddress("vtx_z", &_vtx_z); + _data->SetBranchAddress("vtx_z", &_vtx_z); + + _data->SetBranchStatus("NuMomX", 1); + _data->SetBranchAddress("NuMomX", &_NuMomX); + _data->SetBranchStatus("NuMomY", 1); + _data->SetBranchAddress("NuMomY", &_NuMomY); + _data->SetBranchStatus("NuMomZ", 1); + _data->SetBranchAddress("NuMomZ", &_NuMomZ); + _data->SetBranchStatus("LepMomX", 1); + _data->SetBranchAddress("LepMomX", &_LepMomX); + _data->SetBranchStatus("LepMomY", 1); + _data->SetBranchAddress("LepMomY", &_LepMomY); + _data->SetBranchStatus("LepMomZ", 1); + _data->SetBranchAddress("LepMomZ", &_LepMomZ); TH1D* norm = (TH1D*)_sampleFile->Get("norm"); if(!norm){ @@ -287,116 +297,128 @@ int samplePDFDUNEBeamFD::setupExperimentMC(int iSample) { } // now fill the actual variables - duneobj->norm_s = norm->GetBinContent(1); - duneobj->pot_s = pot/norm->GetBinContent(2); + duneobj.norm_s = norm->GetBinContent(1); + duneobj.pot_s = pot/norm->GetBinContent(2); - std::cout<< "pot_s = " << duneobj->pot_s << std::endl; - std::cout<< "norm_s = " << duneobj->norm_s << std::endl; + std::cout<< "pot_s = " << duneobj.pot_s << std::endl; + std::cout<< "norm_s = " << duneobj.norm_s << std::endl; - duneobj->nEvents = _data->GetEntries(); - duneobj->nutype = nutype; - duneobj->oscnutype = oscnutype; - duneobj->signal = signal; + duneobj.nEvents = _data->GetEntries(); + duneobj.nutype = nutype; + duneobj.oscnutype = oscnutype; + duneobj.signal = signal; - std::cout << "signal: " << duneobj->signal << std::endl; - std::cout << "nevents: " << duneobj->nEvents << std::endl; + std::cout << "signal: " << duneobj.signal << std::endl; + std::cout << "nevents: " << duneobj.nEvents << std::endl; // allocate memory for dunemc variables - duneobj->rw_cvnnumu.resize(duneobj->nEvents); - duneobj->rw_cvnnue.resize(duneobj->nEvents); - duneobj->rw_cvnnumu_shifted.resize(duneobj->nEvents); - duneobj->rw_cvnnue_shifted.resize(duneobj->nEvents); - duneobj->rw_etru.resize(duneobj->nEvents); - duneobj->rw_erec.resize(duneobj->nEvents); - duneobj->rw_erec_shifted.resize(duneobj->nEvents); - duneobj->rw_erec_had.resize(duneobj->nEvents); - duneobj->rw_erec_lep.resize(duneobj->nEvents); - - duneobj->rw_eRecoP.resize(duneobj->nEvents); - duneobj->rw_eRecoPip.resize(duneobj->nEvents); - duneobj->rw_eRecoPim.resize(duneobj->nEvents); - duneobj->rw_eRecoPi0.resize(duneobj->nEvents); - duneobj->rw_eRecoN.resize(duneobj->nEvents); - - duneobj->rw_LepE.resize(duneobj->nEvents); - duneobj->rw_eP.resize(duneobj->nEvents); - duneobj->rw_ePip.resize(duneobj->nEvents); - duneobj->rw_ePim.resize(duneobj->nEvents); - duneobj->rw_ePi0.resize(duneobj->nEvents); - duneobj->rw_eN.resize(duneobj->nEvents); - - duneobj->rw_theta.resize(duneobj->nEvents); - duneobj->flux_w.resize(duneobj->nEvents); - duneobj->rw_isCC.resize(duneobj->nEvents); - duneobj->rw_nuPDGunosc.resize(duneobj->nEvents); - duneobj->rw_nuPDG.resize(duneobj->nEvents); - duneobj->rw_berpaacvwgt.resize(duneobj->nEvents); - duneobj->rw_vtx_x.resize(duneobj->nEvents); - duneobj->rw_vtx_y.resize(duneobj->nEvents); - duneobj->rw_vtx_z.resize(duneobj->nEvents); - - duneobj->global_bin_number.resize(duneobj->nEvents); - - duneobj->mode.resize(duneobj->nEvents); - duneobj->Target.resize(duneobj->nEvents); + duneobj.rw_cvnnumu.resize(duneobj.nEvents); + duneobj.rw_cvnnue.resize(duneobj.nEvents); + duneobj.rw_cvnnumu_shifted.resize(duneobj.nEvents); + duneobj.rw_cvnnue_shifted.resize(duneobj.nEvents); + duneobj.rw_etru.resize(duneobj.nEvents); + duneobj.rw_erec.resize(duneobj.nEvents); + duneobj.rw_erec_shifted.resize(duneobj.nEvents); + duneobj.rw_erec_had.resize(duneobj.nEvents); + duneobj.rw_erec_lep.resize(duneobj.nEvents); + + duneobj.true_q0.resize(duneobj.nEvents); + duneobj.true_q3.resize(duneobj.nEvents); + + duneobj.rw_eRecoP.resize(duneobj.nEvents); + duneobj.rw_eRecoPip.resize(duneobj.nEvents); + duneobj.rw_eRecoPim.resize(duneobj.nEvents); + duneobj.rw_eRecoPi0.resize(duneobj.nEvents); + duneobj.rw_eRecoN.resize(duneobj.nEvents); + + duneobj.rw_LepE.resize(duneobj.nEvents); + duneobj.rw_eP.resize(duneobj.nEvents); + duneobj.rw_ePip.resize(duneobj.nEvents); + duneobj.rw_ePim.resize(duneobj.nEvents); + duneobj.rw_ePi0.resize(duneobj.nEvents); + duneobj.rw_eN.resize(duneobj.nEvents); + + duneobj.rw_theta.resize(duneobj.nEvents); + duneobj.flux_w.resize(duneobj.nEvents); + duneobj.rw_isCC.resize(duneobj.nEvents); + duneobj.rw_nuPDGunosc.resize(duneobj.nEvents); + duneobj.rw_nuPDG.resize(duneobj.nEvents); + duneobj.rw_berpaacvwgt.resize(duneobj.nEvents); + duneobj.rw_vtx_x.resize(duneobj.nEvents); + duneobj.rw_vtx_y.resize(duneobj.nEvents); + duneobj.rw_vtx_z.resize(duneobj.nEvents); + + duneobj.global_bin_number.resize(duneobj.nEvents); + + duneobj.mode.resize(duneobj.nEvents); + duneobj.Target.resize(duneobj.nEvents); _data->GetEntry(0); + bool need_global_bin_numbers = (XVarStr == "global_bin_number"); + //FILL DUNE STRUCT - for (int i = 0; i < duneobj->nEvents; ++i) { // Loop through tree + for (int i = 0; i < duneobj.nEvents; ++i) { // Loop through tree _data->GetEntry(i); - duneobj->rw_cvnnumu[i] = _cvnnumu; - duneobj->rw_cvnnue[i] = _cvnnue; - duneobj->rw_cvnnumu_shifted[i] = _cvnnumu; - duneobj->rw_cvnnue_shifted[i] = _cvnnue; + duneobj.rw_cvnnumu[i] = _cvnnumu; + duneobj.rw_cvnnue[i] = _cvnnue; + duneobj.rw_cvnnumu_shifted[i] = _cvnnumu; + duneobj.rw_cvnnue_shifted[i] = _cvnnue; + if (iselike) { - duneobj->rw_erec[i] = _erec_nue; - duneobj->rw_erec_shifted[i] = _erec_nue; - duneobj->rw_erec_had[i] = _erec_had_nue; - duneobj->rw_erec_lep[i] = _erec_lep_nue; + duneobj.rw_erec[i] = _erec_nue; + duneobj.rw_erec_shifted[i] = _erec_nue; + duneobj.rw_erec_had[i] = _erec_had_nue; + duneobj.rw_erec_lep[i] = _erec_lep_nue; } else { - duneobj->rw_erec[i] = _erec; - duneobj->rw_erec_shifted[i] = _erec; - duneobj->rw_erec_had[i] = _erec_had; - duneobj->rw_erec_lep[i] = _erec_lep; + duneobj.rw_erec[i] = _erec; + duneobj.rw_erec_shifted[i] = _erec; + duneobj.rw_erec_had[i] = _erec_had; + duneobj.rw_erec_lep[i] = _erec_lep; } + + duneobj.true_q0[i] = _ev - _LepE; + duneobj.true_q3[i] = (TVector3{_NuMomX, _NuMomY, _NuMomZ} - + TVector3{_LepMomX, _LepMomY, _LepMomZ}) + .Mag(); + + duneobj.rw_eRecoP[i] = _eRecoP; + duneobj.rw_eRecoPip[i] = _eRecoPip; + duneobj.rw_eRecoPim[i] = _eRecoPim; + duneobj.rw_eRecoPi0[i] = _eRecoPi0; + duneobj.rw_eRecoN[i] = _eRecoN; - duneobj->rw_eRecoP[i] = _eRecoP; - duneobj->rw_eRecoPip[i] = _eRecoPip; - duneobj->rw_eRecoPim[i] = _eRecoPim; - duneobj->rw_eRecoPi0[i] = _eRecoPi0; - duneobj->rw_eRecoN[i] = _eRecoN; - - duneobj->rw_LepE[i] = _LepE; - duneobj->rw_eP[i] = _eP; - duneobj->rw_ePip[i] = _ePip; - duneobj->rw_ePim[i] = _ePim; - duneobj->rw_ePi0[i] = _ePi0; - duneobj->rw_eN[i] = _eN; - - duneobj->rw_etru[i] = _ev; - duneobj->rw_theta[i] = _LepNuAngle; - duneobj->rw_isCC[i] = _isCC; - duneobj->rw_nuPDGunosc[i] = _nuPDGunosc; - duneobj->rw_nuPDG[i] = _nuPDG; - duneobj->rw_berpaacvwgt[i] = _BeRPA_cvwgt; - duneobj->rw_vtx_x[i] = _vtx_x; - duneobj->rw_vtx_y[i] = _vtx_y; - duneobj->rw_vtx_z[i] = _vtx_z; - - duneobj->global_bin_number[i] = GetGenericBinningGlobalBinNumber(iSample, i); + duneobj.rw_LepE[i] = _LepE; + duneobj.rw_eP[i] = _eP; + duneobj.rw_ePip[i] = _ePip; + duneobj.rw_ePim[i] = _ePim; + duneobj.rw_ePi0[i] = _ePi0; + duneobj.rw_eN[i] = _eN; + duneobj.rw_etru[i] = _ev; + duneobj.rw_theta[i] = _LepNuAngle; + duneobj.rw_isCC[i] = _isCC; + duneobj.rw_nuPDGunosc[i] = _nuPDGunosc; + duneobj.rw_nuPDG[i] = _nuPDG; + duneobj.rw_berpaacvwgt[i] = _BeRPA_cvwgt; + duneobj.rw_vtx_x[i] = _vtx_x; + duneobj.rw_vtx_y[i] = _vtx_y; + duneobj.rw_vtx_z[i] = _vtx_z; + + if(need_global_bin_numbers){ + duneobj.global_bin_number[i] = GetGenericBinningGlobalBinNumber(iSample, i); + } //Assume everything is on Argon for now.... - duneobj->Target[i] = 40; + duneobj.Target[i] = 40; int mode= TMath::Abs(_mode); - duneobj->mode[i]=SIMBMode_ToMaCh3Mode(mode, _isCC); + duneobj.mode[i]=SIMBMode_ToMaCh3Mode(mode, _isCC); - duneobj->flux_w[i] = 1.0; + duneobj.flux_w[i] = 1.0; } _sampleFile->Close(); - return duneobj->nEvents; + return duneobj.nEvents; } TH1D* samplePDFDUNEBeamFD::get1DVarHist(KinematicTypes Var1, int kModeToFill, int kChannelToFill, int WeightStyle, TAxis* Axis) { @@ -556,15 +578,91 @@ double const& samplePDFDUNEBeamFD::ReturnKinematicParameterByReference(int Kinem return dunemcSamples[iSample].rw_cvnnue_shifted[iEvent]; case kGlobalBinNumber: return dunemcSamples[iSample].global_bin_number[iEvent]; + case kELepRec: { + return dunemcSamples[iSample].rw_erec_lep[iEvent]; + } + case kq0: + return dunemcSamples[iSample].true_q0[iEvent]; + case kq3: + return dunemcSamples[iSample].true_q3[iEvent]; default: - std::cout << "[ERROR]: " << __FILE__ << ":" << __LINE__ << " Did not recognise Kinematic Parameter type..." << std::endl; - throw; + std::stringstream ss; + ss << "[ERROR]: " << __FILE__ << ":" << __LINE__ + << " ReturnKinematicParameterByReference Did not recognise " + "Kinematic Parameter type:" + << KinematicParameter + << ". Is it possibly only available via ReturnKinematicParameter " + "(not by reference?) if you need it here, give it storage in " + "dunemc_base and move it to here." + ; + throw std::runtime_error(ss.str()); } } -double samplePDFDUNEBeamFD::ReturnKinematicParameter(int KinematicParameter, int iSample, int iEvent) { - return ReturnKinematicParameterByReference(KinematicParameter, iSample, iEvent); -} + double samplePDFDUNEBeamFD::ReturnKinematicParameter(int KinematicParameter, + int iSample, + int iEvent) { + switch (KinematicParameter) { + case kERecQE: { + + constexpr double V = 0; // 0 binding energy for now + constexpr double mn = 939.565; // neutron mass + constexpr double mp = 938.272; // proton mass + + double mN_eff = mn - V; + double mN_oth = mp; + + if (dunemcSamples[iSample].rw_nuPDGunosc[iEvent] < + 0) { // if anti-neutrino, swap target/out masses + mN_eff = mp - V; + mN_oth = mn; + } + + double el = dunemcSamples[iSample].rw_erec_lep[iEvent]; + + // this is funky, but don't be scared, it defines an annonymous function + // in place that grabs the lepton mass in MeV when given the neutrino PDG + // and whether the interaction was CC or NC and then immediately calls it. + // It's basically a generalisation of the ternary operator. + double ml = + [](int nupdg, bool isCC) { + switch (std::abs(nupdg)) { + case 12: { + return isCC ? 0.511 : 0; + } + case 14: { + return isCC ? 105.66 : 0; + } + case 16: { + return isCC ? 1777.0 : 0; + } + } + }(dunemcSamples[iSample].rw_nuPDGunosc[iEvent], + dunemcSamples[iSample].rw_isCC[iEvent]); + + double pl = std::sqrt(el*el - ml*ml); // momentum of lepton + + double rEnu = + (2 * mN_eff * el - ml * ml + mN_oth * mN_oth - mN_eff * mN_eff) / + (2 * (mN_eff - el + + pl * std::cos(dunemcSamples[iSample].rw_theta[iEvent]))); + + return rEnu; + } + case kEHadRec: { + + return dunemcSamples[iSample].rw_eRecoP[iEvent] + + dunemcSamples[iSample].rw_eRecoPip[iEvent] + + dunemcSamples[iSample].rw_eRecoPim[iEvent] + + dunemcSamples[iSample].rw_eRecoPi0[iEvent] + + dunemcSamples[iSample].rw_eRecoN[iEvent]; + } + default: { + return ReturnKinematicParameterByReference(KinematicParameter, iSample, + iEvent); + } + } + } int samplePDFDUNEBeamFD::ReturnKinematicParameterFromString(std::string KinematicParameterStr){ @@ -577,62 +675,74 @@ int samplePDFDUNEBeamFD::ReturnKinematicParameterFromString(std::string Kinemati if (KinematicParameterStr.find("CVNNue") != std::string::npos) {return kCVNNue;} if (KinematicParameterStr.find("M3Mode") != std::string::npos) {return kM3Mode;} if (KinematicParameterStr.find("global_bin_number") != std::string::npos) {return kGlobalBinNumber;} - - + if (KinematicParameterStr.find("q0") != std::string::npos) {return kq0;} + if (KinematicParameterStr.find("q3") != std::string::npos) {return kq3;} + if (KinematicParameterStr.find("ERecQE") != std::string::npos) {return kERecQE;} + if (KinematicParameterStr.find("ELepRec") != std::string::npos) {return kELepRec;} + if (KinematicParameterStr.find("EHadRec") != std::string::npos) {return kEHadRec;} + + std::stringstream ss; + ss << "[ERROR]: " << __FILE__ << ":" << __LINE__ + << "failed to translate kinematic parameter string " + << KinematicParameterStr << " to parameter id."; + throw std::runtime_error(ss.str()); } -inline std::string samplePDFDUNEBeamFD::ReturnStringFromKinematicParameter(int KinematicParameter) { - std::string KinematicString = ""; - - switch(KinematicParameter){ - case kRecoNeutrinoEnergy: - KinematicString = "RecoNeutrinoEnergy"; - break; - case kTrueNeutrinoEnergy: - KinematicString = "RecoNeutrinoEnergy"; - break; - case kTrueXPos: - KinematicString= "TrueXPos"; - break; - case kTrueYPos: - KinematicString= "TrueYPos"; - break; - case kTrueZPos: - KinematicString= "TrueZPos"; - break; - case kCVNNumu: - KinematicString = "CVNNumu"; - break; - case kCVNNue: - KinematicString = "CVNNue"; - break; - case kM3Mode: - KinematicString = "M3Mode"; - break; - case kGlobalBinNumber: - KinematicString = "global_bin_number"; - break; - default: - break; - } +std::string samplePDFDUNEBeamFD::ReturnStringFromKinematicParameter( + int KinematicParameter) { - return KinematicString; + switch (KinematicParameter) { + case kRecoNeutrinoEnergy: + return "RecoNeutrinoEnergy"; + case kTrueNeutrinoEnergy: + return "RecoNeutrinoEnergy"; + case kTrueXPos: + return "TrueXPos"; + case kTrueYPos: + return "TrueYPos"; + case kTrueZPos: + return "TrueZPos"; + case kCVNNumu: + return "CVNNumu"; + case kCVNNue: + return "CVNNue"; + case kM3Mode: + return "M3Mode"; + case kGlobalBinNumber: + return "global_bin_number"; + case kq0: + return "q0"; + case kq3: + return "q3"; + case kERecQE: + return "ERecQE"; + case kELepRec: + return "ELepRec"; + case kEHadRec: + return "EHadRec"; + default: { + std::stringstream ss; + ss << "[ERROR]: " << __FILE__ << ":" << __LINE__ + << "failed to get string from parameter id: " << KinematicParameter; + throw std::runtime_error(ss.str()); + } + } } void samplePDFDUNEBeamFD::setupFDMC(int iSample) { - dunemc_base *duneobj = &(dunemcSamples[iSample]); + auto &duneobj = dunemcSamples[iSample]; fdmc_base *fdobj = &(MCSamples[iSample]); - fdobj->nutype = duneobj->nutype; - fdobj->oscnutype = duneobj->oscnutype; - fdobj->signal = duneobj->signal; + fdobj->nutype = duneobj.nutype; + fdobj->oscnutype = duneobj.oscnutype; + fdobj->signal = duneobj.signal; fdobj->SampleDetID = SampleDetID; for(int iEvent = 0 ;iEvent < fdobj->nEvents ; ++iEvent) { - fdobj->rw_etru[iEvent] = &(duneobj->rw_etru[iEvent]); - fdobj->mode[iEvent] = &(duneobj->mode[iEvent]); - fdobj->Target[iEvent] = &(duneobj->Target[iEvent]); - fdobj->isNC[iEvent] = !(duneobj->rw_isCC[iEvent]); + fdobj->rw_etru[iEvent] = &(duneobj.rw_etru[iEvent]); + fdobj->mode[iEvent] = &(duneobj.mode[iEvent]); + fdobj->Target[iEvent] = &(duneobj.Target[iEvent]); + fdobj->isNC[iEvent] = !(duneobj.rw_isCC[iEvent]); } } diff --git a/samplePDFDUNE/samplePDFDUNEBeamFD.h b/samplePDFDUNE/samplePDFDUNEBeamFD.h index d59de6a..4a55fd4 100644 --- a/samplePDFDUNE/samplePDFDUNEBeamFD.h +++ b/samplePDFDUNE/samplePDFDUNEBeamFD.h @@ -25,7 +25,23 @@ class samplePDFDUNEBeamFD : virtual public samplePDFFDBase samplePDFDUNEBeamFD(std::string mc_version, covarianceXsec* xsec_cov); ~samplePDFDUNEBeamFD(); - enum KinematicTypes {kTrueNeutrinoEnergy,kRecoNeutrinoEnergy,kTrueXPos,kTrueYPos,kTrueZPos,kCVNNumu,kCVNNue,kM3Mode,kGlobalBinNumber,kOscChannel}; + enum KinematicTypes { + kTrueNeutrinoEnergy, + kRecoNeutrinoEnergy, + kTrueXPos, + kTrueYPos, + kTrueZPos, + kCVNNumu, + kCVNNue, + kM3Mode, + kGlobalBinNumber, + kOscChannel, + kq0, + kq3, + kERecQE, + kELepRec, + kEHadRec + }; //More robust getters to make plots in different variables, mode, osc channel, systematic weighting and with bin range TH1D* get1DVarHist(KinematicTypes Var1, int fModeToFill=-1, int fSampleToFill=-1, int WeightStyle=0, TAxis* Axis=0); From 9580f7045467c16967c1109e2698b621c1086b2f Mon Sep 17 00:00:00 2001 From: Luke Pickering Date: Thu, 17 Oct 2024 14:26:51 +0100 Subject: [PATCH 7/9] adds example usage of generic binning --- configs/EventRates_Beam.yaml | 9 +- .../Samples/SamplePDFDune_FHC_numuselec.yaml | 52 ++++++--- src/EventRates.cpp | 103 +++++++++++------- 3 files changed, 110 insertions(+), 54 deletions(-) diff --git a/configs/EventRates_Beam.yaml b/configs/EventRates_Beam.yaml index 15c538e..524c52d 100644 --- a/configs/EventRates_Beam.yaml +++ b/configs/EventRates_Beam.yaml @@ -1,11 +1,17 @@ General: OutputFile: "DuneEventRates.root" - DUNESamples: ["configs/Samples/SamplePDFDune_FHC_numuselec.yaml", "configs/Samples/SamplePDFDune_FHC_nueselec.yaml", "configs/Samples/SamplePDFDune_RHC_numuselec.yaml", "configs/Samples/SamplePDFDune_RHC_nueselec.yaml"] + DUNESamples: + - "configs/Samples/SamplePDFDune_FHC_numuselec.yaml" + # - "configs/Samples/SamplePDFDune_FHC_nueselec.yaml" + # - "configs/Samples/SamplePDFDune_RHC_numuselec.yaml" + # - "configs/Samples/SamplePDFDune_RHC_nueselec.yaml" + #Nu-FIT #OscillationParameters: [0.310, 0.582, 0.224, 7.39E-5, 2.5254E-3, -2.498] # T2K-like best-fit OscillationParameters: [0.307, 0.528, 0.0218, 7.53e-5, 2.509e-3, -1.601, 1284.9, 2.848] OscillatorConfigName: "configs/OscillatorObj.yaml" + Systematics: XsecCovFile: ["configs/CovObjs/xsec_covariance_DUNE_systs_2022a_FD_v3_xsec.yaml"] XsecCovName: "xsec_cov" @@ -22,6 +28,7 @@ General: Output: FileName: "TestEventRates.root" OUTPUTNAME: "TestLLH.root" + ProcessMCMC: No Seed: 0 Debug: No diff --git a/configs/Samples/SamplePDFDune_FHC_numuselec.yaml b/configs/Samples/SamplePDFDune_FHC_numuselec.yaml index cdd67d4..bd1c263 100644 --- a/configs/Samples/SamplePDFDune_FHC_numuselec.yaml +++ b/configs/Samples/SamplePDFDune_FHC_numuselec.yaml @@ -10,16 +10,36 @@ SelectionCuts: Bounds: [ 50.0, 1244.0 ] - KinematicStr: "CVNNumu" Bounds: [ 0.5, 999 ] + Binning: - XVarStr : "RecoNeutrinoEnergy" - XVarBins: [0., 0.5, 1., 1.25, 1.5, 1.75, 2., 2.25, 2.5, 2.75, 3., 3.25, 3.5, 3.75, 4., 5., 6., 10.] + Axes: +# Appropriate values for VarStr +# TrueNeutrinoEnergy +# RecoNeutrinoEnergy +# CVNNumu +# CVNNue +# q0 +# q3 +# ERecQE +# ELepRec +# EHadRec + - VarStr: "RecoNeutrinoEnergy" + VarBins: [0,0.5,1,1.25,1.5,1.75,2,2.25,2.5,2.75,3,3.25,3.5,4,4.5,5,6,7,10] + Title: "E_{#nu}^{Rec} [GeV]" + - VarStr: "TrueNeutrinoEnergy" + Uniform: [100,0,10] + Title: "E_{#nu}^{True} [GeV]" + # if you want to use the generic binning with a single + # axis, useful for calculated VarStr that cannot be referenced + ForceGeneric: True + DUNESampleBools: iselike: no InputFiles: mtupleprefix: "inputs/DUNE_CAF_files/FD_FHC_ger_" mtuplesuffix: "_numuselec.root" - splineprefix: "inputs/DUNE_spline_files/FD_FHC_ger_" - splinesuffix: "_numuselec_splines.root" + # splineprefix: "inputs/DUNE_spline_files/FD_FHC_ger_" + # splinesuffix: "_numuselec_splines.root" DetID: 24 NuOsc: NuOscConfigFile: "configs/NuOsc/CUDAProb3Linear.yaml" @@ -27,84 +47,84 @@ NSubSamples: 12 SubSamples: - name: "FHC_numu_x_numu" mtuplefile: "numu_x_numu" - splinefile: "numu_x_numu" + # splinefile: "numu_x_numu" samplevecno: 0 nutype: 14 oscnutype: 14 signal: false - name: "FHC_nue_x_nue" mtuplefile: "nue_x_nue" - splinefile: "nue_x_nue" + # splinefile: "nue_x_nue" samplevecno: 1 nutype: 12 oscnutype: 12 signal: false - name: "FHC_numubar_x_numubar" mtuplefile: "numubar_x_numubar" - splinefile: "numubar_x_numubar" + # splinefile: "numubar_x_numubar" samplevecno: 2 nutype: -14 oscnutype: -14 signal: false - name: "FHC_nuebar_x_nuebar" mtuplefile: "nuebar_x_nuebar" - splinefile: "nuebar_x_nuebar" + # splinefile: "nuebar_x_nuebar" samplevecno: 3 nutype: -12 oscnutype: -12 signal: false - name: "FHC_numu_x_nue" mtuplefile: "numu_x_nue" - splinefile: "numu_x_nue" + # splinefile: "numu_x_nue" samplevecno: 4 nutype: 14 oscnutype: 12 signal: true - name: "FHC_numubar_x_nuebar" mtuplefile: "numubar_x_nuebar" - splinefile: "numubar_x_nuebar" + # splinefile: "numubar_x_nuebar" samplevecno: 5 nutype: -14 oscnutype: -12 signal: true - name: "FHC_nue_x_numu" mtuplefile: "nue_x_numu" - splinefile: "nue_x_numu" + # splinefile: "nue_x_numu" samplevecno: 6 nutype: 12 oscnutype: 14 signal: true - name: "FHC_nuebar_x_numubar" mtuplefile: "nuebar_x_numubar" - splinefile: "nuebar_x_numubar" + # splinefile: "nuebar_x_numubar" samplevecno: 7 nutype: -12 oscnutype: -14 signal: true - name: "FHC_numu_x_nutau" mtuplefile: "numu_x_nutau" - splinefile: "numu_x_nutau" + # splinefile: "numu_x_nutau" samplevecno: 8 nutype: 14 oscnutype: 16 signal: true - name: "FHC_nue_x_nutau" mtuplefile: "nue_x_nutau" - splinefile: "nue_x_nutau" + # splinefile: "nue_x_nutau" samplevecno: 9 nutype: 12 oscnutype: 16 signal: true - name: "FHC_numubar_x_nutaubar" mtuplefile: "numubar_x_nutaubar" - splinefile: "numubar_x_nutaubar" + # splinefile: "numubar_x_nutaubar" samplevecno: 10 nutype: -14 oscnutype: -16 signal: true - name: "FHC_nuebar_x_nutaubar" mtuplefile: "nuebar_x_nutaubar" - splinefile: "nuebar_x_nutaubar" + # splinefile: "nuebar_x_nutaubar" samplevecno: 11 nutype: -12 oscnutype: -16 diff --git a/src/EventRates.cpp b/src/EventRates.cpp index 7e70735..0d422f9 100644 --- a/src/EventRates.cpp +++ b/src/EventRates.cpp @@ -1,25 +1,29 @@ -#include #include #include +#include #include +#include +#include #include #include -#include -#include -#include #include -#include #include +#include +#include + +#include "samplePDF/GenericBinningTools.h" #include "samplePDFDUNE/MaCh3DUNEFactory.h" -void Write1DHistogramsToFile(std::string OutFileName, std::vector Histograms) { +void Write1DHistogramsToFile(std::string OutFileName, + std::vector Histograms) { - // Now write out the saved hsitograms to file - auto OutputFile = std::unique_ptr(new TFile(OutFileName.c_str(), "RECREATE")); + // Now write out the saved hsitograms to file + auto OutputFile = + std::unique_ptr(new TFile(OutFileName.c_str(), "RECREATE")); OutputFile->cd(); - for(auto Hist : Histograms){ + for (auto Hist : Histograms) { Hist->Write(); } OutputFile->Close(); @@ -27,57 +31,82 @@ void Write1DHistogramsToFile(std::string OutFileName, std::vector Histogr return; } -void Write1DHistogramsToPdf(std::string OutFileName, std::vector Histograms) { +void Write1DHistogramsToPdf(std::string OutFileName, + std::vector Histograms) { - // Now write out the saved hsitograms to file + // Now write out the saved hsitograms to file - - //Remove root from end of file + // Remove root from end of file OutFileName.erase(OutFileName.find('.')); - OutFileName+=".pdf"; + OutFileName += ".pdf"; auto c1 = std::unique_ptr(new TCanvas("c1", "c1", 800, 600)); c1->cd(); - c1->Print(std::string(OutFileName+"[").c_str()); - for(auto Hist : Histograms){ + c1->Print(std::string(OutFileName + "[").c_str()); + for (auto Hist : Histograms) { Hist->Draw("HIST"); - c1->Print(OutFileName.c_str()); + c1->Print(OutFileName.c_str()); } - c1->Print(std::string(OutFileName+"]").c_str()); + c1->Print(std::string(OutFileName + "]").c_str()); return; } -int main(int argc, char * argv[]) { - if(argc == 1){ +int main(int argc, char *argv[]) { + if (argc == 1) { std::cout << "Usage: bin/EventRatesDUNEBeam config.cfg" << std::endl; return 1; } auto fitMan = std::unique_ptr(new manager(argv[1])); - covarianceXsec* xsec = nullptr; - covarianceOsc* osc = nullptr; - - //#################################################################################### - //Create samplePDFFD objects - - std::vector DUNEPdfs; - MakeMaCh3DuneInstance(fitMan.get(), DUNEPdfs, xsec, osc); + covarianceXsec *xsec = nullptr; + covarianceOsc *osc = nullptr; + + // #################################################################################### + // Create samplePDFFD objects + + std::vector DUNEPdfs; + MakeMaCh3DuneInstance(fitMan.get(), DUNEPdfs, xsec, osc); + + auto gc1 = std::unique_ptr(new TCanvas("gc1", "gc1", 800, 600)); + gc1->Print("GenericBinTest.pdf["); + + std::vector DUNEHists; + for (auto Sample : DUNEPdfs) { + Sample->reweight(); + + if(Sample->generic_binning.GetNDimensions()){ + auto myhist = GetGenericBinningTH1(*Sample, "myhist"); + myhist->Scale(1,"WIDTH"); + myhist->Draw(); + gc1->Print("GenericBinTest.pdf"); + if(Sample->generic_binning.GetNDimensions() == 2){ + auto myhist2 = GetGenericBinningTH2(*Sample, "myhist2"); + myhist2->Scale(1,"WIDTH"); + myhist2->Draw("COLZ"); + gc1->Print("GenericBinTest.pdf"); + for(auto &slice : GetGenericBinningTH1Slices(*Sample, 0, "myslicehist")){ + slice->Draw("COLZ"); + gc1->Print("GenericBinTest.pdf"); + } + } + } - - std::vector DUNEHists; - for(auto Sample : DUNEPdfs){ - Sample->reweight(); DUNEHists.push_back(Sample->get1DHist()); - std::string EventRateString = fmt::format("{:.2f}", Sample->get1DHist()->Integral()); - MACH3LOG_INFO("Event rate for {} : {:<5}", Sample->GetName(), EventRateString); - + std::string EventRateString = + fmt::format("{:.2f}", Sample->get1DHist()->Integral()); + MACH3LOG_INFO("Event rate for {} : {:<5}", Sample->GetName(), + EventRateString); } - std::string OutFileName = GetFromManager(fitMan->raw()["General"]["OutputFile"], "EventRatesOutput.root"); + gc1->Print("GenericBinTest.pdf]"); + + + std::string OutFileName = GetFromManager( + fitMan->raw()["General"]["OutputFile"], "EventRatesOutput.root"); - Write1DHistogramsToFile(OutFileName, DUNEHists); + Write1DHistogramsToFile(OutFileName, DUNEHists); Write1DHistogramsToPdf(OutFileName, DUNEHists); } From 2057400db7dd601b7182d39c8ab5c86c3d5aeb1e Mon Sep 17 00:00:00 2001 From: Luke Pickering Date: Thu, 17 Oct 2024 15:23:55 +0100 Subject: [PATCH 8/9] point the MaCH3_Core branch at my branch for the time being --- CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 1d7f8d5..71e9bbf 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -70,7 +70,7 @@ find_package(MaCh3) if(NOT MaCh3_FOUND) CPMFindPackage( NAME MaCh3 - GIT_TAG "feature_OscClassMagic" + GIT_TAG "picker24/feature/GenericBinningAndReturnKinematicParameterTidy" GITHUB_REPOSITORY mach3-software/MaCh3 ) else() From 13d7fd0cb21ceb024bc337c9bbf3a32232972dbd Mon Sep 17 00:00:00 2001 From: Luke Pickering Date: Fri, 18 Oct 2024 13:53:01 +0100 Subject: [PATCH 9/9] adds 2d slice example to eventrates --- .../Samples/SamplePDFDune_FHC_numuselec.yaml | 11 +++++---- src/EventRates.cpp | 23 ++++++++++++++----- 2 files changed, 24 insertions(+), 10 deletions(-) diff --git a/configs/Samples/SamplePDFDune_FHC_numuselec.yaml b/configs/Samples/SamplePDFDune_FHC_numuselec.yaml index bd1c263..a2df937 100644 --- a/configs/Samples/SamplePDFDune_FHC_numuselec.yaml +++ b/configs/Samples/SamplePDFDune_FHC_numuselec.yaml @@ -23,12 +23,15 @@ Binning: # ERecQE # ELepRec # EHadRec + - VarStr: "ELepRec" + Uniform: [25,0,5] + Title: "E_{lep.}^{Rec} [GeV]" + - VarStr: "EHadRec" + Uniform: [25,0,5] + Title: "E_{had.}^{Vis} [GeV]" - VarStr: "RecoNeutrinoEnergy" - VarBins: [0,0.5,1,1.25,1.5,1.75,2,2.25,2.5,2.75,3,3.25,3.5,4,4.5,5,6,7,10] + Uniform: [10,0.5,5.5] Title: "E_{#nu}^{Rec} [GeV]" - - VarStr: "TrueNeutrinoEnergy" - Uniform: [100,0,10] - Title: "E_{#nu}^{True} [GeV]" # if you want to use the generic binning with a single # axis, useful for calculated VarStr that cannot be referenced ForceGeneric: True diff --git a/src/EventRates.cpp b/src/EventRates.cpp index 0d422f9..42f03d0 100644 --- a/src/EventRates.cpp +++ b/src/EventRates.cpp @@ -70,24 +70,35 @@ int main(int argc, char *argv[]) { MakeMaCh3DuneInstance(fitMan.get(), DUNEPdfs, xsec, osc); auto gc1 = std::unique_ptr(new TCanvas("gc1", "gc1", 800, 600)); + gStyle->SetOptStat(false); gc1->Print("GenericBinTest.pdf["); std::vector DUNEHists; for (auto Sample : DUNEPdfs) { Sample->reweight(); - if(Sample->generic_binning.GetNDimensions()){ + if (Sample->generic_binning.GetNDimensions()) { + auto myhist = GetGenericBinningTH1(*Sample, "myhist"); - myhist->Scale(1,"WIDTH"); + myhist->Scale(1, "WIDTH"); myhist->Draw(); gc1->Print("GenericBinTest.pdf"); - if(Sample->generic_binning.GetNDimensions() == 2){ + + if (Sample->generic_binning.GetNDimensions() == 2) { auto myhist2 = GetGenericBinningTH2(*Sample, "myhist2"); - myhist2->Scale(1,"WIDTH"); myhist2->Draw("COLZ"); gc1->Print("GenericBinTest.pdf"); - for(auto &slice : GetGenericBinningTH1Slices(*Sample, 0, "myslicehist")){ - slice->Draw("COLZ"); + + for (auto &slice : + GetGenericBinningTH1Slices(*Sample, 0, "myslicehist")) { + slice->Draw(); + gc1->Print("GenericBinTest.pdf"); + } + } + if (Sample->generic_binning.GetNDimensions() == 3) { + for (auto &slice : + GetGenericBinningTH2Slices(*Sample, {0, 1}, "myslicehist")) { + slice->Draw(); gc1->Print("GenericBinTest.pdf"); } }