diff --git a/CMakeLists.txt b/CMakeLists.txt index 10892da..8f21629 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) @@ -74,7 +70,7 @@ find_package(MaCh3) if(NOT MaCh3_FOUND) CPMFindPackage( NAME MaCh3 - GIT_TAG "develop" + GIT_TAG "picker24/feature/GenericBinningAndReturnKinematicParameterTidy" GITHUB_REPOSITORY mach3-software/MaCh3 ) else() @@ -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) 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..a2df937 100644 --- a/configs/Samples/SamplePDFDune_FHC_numuselec.yaml +++ b/configs/Samples/SamplePDFDune_FHC_numuselec.yaml @@ -10,16 +10,39 @@ 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: "ELepRec" + Uniform: [25,0,5] + Title: "E_{lep.}^{Rec} [GeV]" + - VarStr: "EHadRec" + Uniform: [25,0,5] + Title: "E_{had.}^{Vis} [GeV]" + - VarStr: "RecoNeutrinoEnergy" + Uniform: [10,0.5,5.5] + Title: "E_{#nu}^{Rec} [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 +50,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/samplePDFDUNE/StructsDUNE.h b/samplePDFDUNE/StructsDUNE.h index fc87e48..d406f4a 100644 --- a/samplePDFDUNE/StructsDUNE.h +++ b/samplePDFDUNE/StructsDUNE.h @@ -6,85 +6,91 @@ 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 true_q0; + std::vector true_q3; + + 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; + bool * 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 2919f81..5635f0d 100644 --- a/samplePDFDUNE/samplePDFDUNEAtm.cpp +++ b/samplePDFDUNE/samplePDFDUNEAtm.cpp @@ -66,15 +66,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); @@ -136,49 +136,28 @@ void samplePDFDUNEAtm::setupFDMC(int iSample) { } } -const double* samplePDFDUNEAtm::GetPointerToKinematicParameter(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: - MACH3LOG_ERROR("Unknown KinPar: {}",KinPar); + MACH3LOG_ERROR("Unknown KinematicParameter: {}",KinematicParameter); throw MaCh3Exception(__FILE__, __LINE__); } - - return KinematicValue; -} - -const double* samplePDFDUNEAtm::GetPointerToKinematicParameter(double KinematicVariable, int iSample, int iEvent) { - KinematicTypes KinPar = (KinematicTypes) std::round(KinematicVariable); - return GetPointerToKinematicParameter(KinPar,iSample,iEvent); -} - -const double* samplePDFDUNEAtm::GetPointerToKinematicParameter(std::string KinematicParameter, int iSample, int iEvent) { - KinematicTypes KinPar = static_cast(ReturnKinematicParameterFromString(KinematicParameter)); - return GetPointerToKinematicParameter(KinPar,iSample,iEvent); -} - -double samplePDFDUNEAtm::ReturnKinematicParameter(double KinematicVariable, int iSample, int iEvent) { - return *GetPointerToKinematicParameter(KinematicVariable, iSample, iEvent); } -double samplePDFDUNEAtm::ReturnKinematicParameter(std::string KinematicParameter, int iSample, int iEvent) { - return *GetPointerToKinematicParameter(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 7cd4653..7aeb797 100644 --- a/samplePDFDUNE/samplePDFDUNEAtm.h +++ b/samplePDFDUNE/samplePDFDUNEAtm.h @@ -28,18 +28,15 @@ class samplePDFDUNEAtm : virtual public samplePDFFDBase double CalcXsecWeightFunc(int iSample, int iEvent) {return 1.;} void applyShifts(int iSample, int iEvent) {} - const double* GetPointerToKinematicParameter(KinematicTypes KinPar, int iSample, int iEvent); - const double* GetPointerToKinematicParameter(double KinematicVariable, int iSample, int iEvent); - const double* GetPointerToKinematicParameter(std::string KinematicParameter, int iSample, int iEvent); + double const& ReturnKinematicParameterByReference(int KinematicParameter, int iSample, int iEvent); + double ReturnKinematicParameter(int 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); + std::vector ReturnKinematicParameterBinning(int KinematicParameter); + + int ReturnKinematicParameterFromString(std::string KinematicStr); + std::string ReturnStringFromKinematicParameter(int KinematicVariable); - std::vector dunemcSamples; + std::vector dunemcSamples; bool IsELike; }; diff --git a/samplePDFDUNE/samplePDFDUNEBeamFD.cpp b/samplePDFDUNE/samplePDFDUNEBeamFD.cpp index 339de75..6f79c89 100644 --- a/samplePDFDUNE/samplePDFDUNEBeamFD.cpp +++ b/samplePDFDUNE/samplePDFDUNEBeamFD.cpp @@ -53,7 +53,10 @@ void samplePDFDUNEBeamFD::Init() { em_res_fd_pos = -999; cvn_numu_fd_pos = -999; cvn_nue_fd_pos = -999; - + + // create dunemc storage + dunemcSamples.resize(nSamples); + nFDDetectorSystPointers = funcParsIndex.size(); std::unordered_map FDDetectorSystPointersMap; FDDetectorSystPointers = std::vector(nFDDetectorSystPointers); @@ -165,13 +168,11 @@ void samplePDFDUNEBeamFD::Init() { void samplePDFDUNEBeamFD::SetupSplines() { ///@todo move all of the spline setup into core - if(XsecCov->GetNumParamsFromDetID(SampleDetID, kSpline) > 0){ + 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)); - splinesDUNE* DUNESplines = new splinesDUNE(XsecCov); - splineFile = (splineFDBase*)DUNESplines; + splineFile = new splinesDUNE(XsecCov); InitialiseSplineObject(); - } - else{ + } else { MACH3LOG_INFO("Found {} splines for this sample so I will not load or evaluate splines", XsecCov->GetNumParamsFromDetID(SampleDetID, kSpline)); splineFile = nullptr; } @@ -197,24 +198,25 @@ void samplePDFDUNEBeamFD::SetupWeightPointers() { int samplePDFDUNEBeamFD::setupExperimentMC(int iSample) { - dunemc_base *duneobj = &(dunemcSamples[iSample]); + auto &duneobj = dunemcSamples[iSample]; + int nutype = sample_nutype[iSample]; int oscnutype = sample_oscnutype[iSample]; bool signal = sample_signal[iSample]; MACH3LOG_INFO("-------------------------------------------------------------------"); - MACH3LOG_INFO("input file: {}", mc_files[iSample]); + MACH3LOG_INFO("input file: {}", mc_files[iSample].native()); _sampleFile = new TFile(mc_files[iSample].c_str(), "READ"); _data = (TTree*)_sampleFile->Get("caf"); if(_data){ - MACH3LOG_INFO("Found \"caf\" tree in {}", mc_files[iSample]); + MACH3LOG_INFO("Found \"caf\" tree in {}", mc_files[iSample].native()); MACH3LOG_INFO("With number of entries: {}", _data->GetEntries()); } else{ - MACH3LOG_ERROR("Could not find \"caf\" tree in {}", mc_files[iSample]); - throw MaCh3Exception(__FILE__, __LINE__); + MACH3LOG_ERROR("Could not find \"caf\" tree in {}", mc_files[iSample].native()); + throw MaCh3Exception(__FILE__, __LINE__); } _data->SetBranchStatus("*", 0); @@ -276,7 +278,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){ @@ -285,106 +300,122 @@ 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); - 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; // 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.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] = (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.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] = (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_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_etru[i] = (double)_ev; - duneobj->rw_theta[i] = (double)_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_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) { @@ -470,7 +501,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 @@ -522,112 +553,109 @@ 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]; + 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: - MACH3LOG_ERROR("Did not recognise Kinematic Parameter type"); - MACH3LOG_ERROR("Was given a Kinematic Variable of {}", KinematicVariable); - throw MaCh3Exception(__FILE__, __LINE__); + 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()); + } } - - 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: - MACH3LOG_ERROR("Did not recognise Kinematic Parameter type..."); - throw MaCh3Exception(__FILE__, __LINE__); - } - - return KinematicValue; -} + 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: { -const double* samplePDFDUNEBeamFD::GetPointerToKinematicParameter(std::string KinematicParameter, int iSample, int iEvent) { - KinematicTypes KinPar = static_cast(ReturnKinematicParameterFromString(KinematicParameter)); - double* KinematicValue = nullptr; - - 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: - MACH3LOG_ERROR("Did not recognise Kinematic Parameter type..."); - throw MaCh3Exception(__FILE__, __LINE__); - } - - return KinematicValue; -} + 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){ if (KinematicParameterStr.find("TrueNeutrinoEnergy") != std::string::npos) {return kTrueNeutrinoEnergy;} @@ -638,91 +666,73 @@ 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;} + 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()); } -const double* samplePDFDUNEBeamFD::GetPointerToKinematicParameter(double KinematicVariable, int iSample, int iEvent) { - KinematicTypes KinPar = (KinematicTypes) std::round(KinematicVariable); - double* KinematicValue = nullptr; +std::string samplePDFDUNEBeamFD::ReturnStringFromKinematicParameter( + int KinematicParameter) { - switch(KinPar){ - case kTrueNeutrinoEnergy: - KinematicValue = &dunemcSamples[iSample].rw_etru[iEvent]; - break; + switch (KinematicParameter) { case kRecoNeutrinoEnergy: - KinematicValue = &dunemcSamples[iSample].rw_erec_shifted[iEvent]; - break; + return "RecoNeutrinoEnergy"; + case kTrueNeutrinoEnergy: + return "RecoNeutrinoEnergy"; case kTrueXPos: - KinematicValue = &dunemcSamples[iSample].rw_vtx_x[iEvent]; - break; + return "TrueXPos"; case kTrueYPos: - KinematicValue = &dunemcSamples[iSample].rw_vtx_y[iEvent]; - break; + return "TrueYPos"; case kTrueZPos: - KinematicValue = &dunemcSamples[iSample].rw_vtx_z[iEvent]; - break; + return "TrueZPos"; case kCVNNumu: - KinematicValue = &dunemcSamples[iSample].rw_cvnnumu_shifted[iEvent]; - break; + return "CVNNumu"; case kCVNNue: - KinematicValue = &dunemcSamples[iSample].rw_cvnnue_shifted[iEvent]; - break; - default: + 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: { MACH3LOG_ERROR("Did not recognise Kinematic Parameter type..."); throw MaCh3Exception(__FILE__, __LINE__); } - - return KinematicValue; -} - -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; - default: - break; } - - return KinematicString; } 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]); } } @@ -806,9 +816,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 2dc477b..0f86875 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,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); @@ -41,23 +57,23 @@ 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 ReturnKinematicParameter (double KinematicVariable, int iSample, int iEvent); - double ReturnKinematicParameter(std::string KinematicParameter, int iSample, int iEvent); - const double* GetPointerToKinematicParameter(std::string KinematicParameter, int iSample, int iEvent); - const double* GetPointerToKinematicParameter(double KinematicVariable, int iSample, int iEvent); + double const &ReturnKinematicParameterByReference(int KinematicParameter, + int iSample, int iEvent); + double ReturnKinematicParameter(int KinematicParameter, int iSample, + int iEvent); - std::vector ReturnKinematicParameterBinning(std::string KinematicParameter); - inline int ReturnKinematicParameterFromString(std::string KinematicParameterStr); - inline std::string ReturnStringFromKinematicParameter(int KinematicParameterStr); + std::vector ReturnKinematicParameterBinning(int KinematicParameter); + 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); // dunemc - std::vector dunemcSamples; + std::vector dunemcSamples; double pot; diff --git a/samplePDFDUNE/samplePDFDUNEBeamND.cpp b/samplePDFDUNE/samplePDFDUNEBeamND.cpp index 4ba29e7..b70ad4a 100644 --- a/samplePDFDUNE/samplePDFDUNEBeamND.cpp +++ b/samplePDFDUNE/samplePDFDUNEBeamND.cpp @@ -165,24 +165,23 @@ void samplePDFDUNEBeamND::SetupWeightPointers() { } int samplePDFDUNEBeamND::setupExperimentMC(int iSample) { - dunemc_base *duneobj = &(dunendmcSamples[iSample]); int nutype = sample_nutype[iSample]; int oscnutype = sample_oscnutype[iSample]; bool signal = sample_signal[iSample]; - MACH3LOG_INFO("-------------------------------------------------------------------"); - MACH3LOG_INFO("input file: {}", mc_files.at(iSample)); + std::cout << "-------------------------------------------------------------------" << std::endl; + std::cout << "input file: " << mc_files[iSample] << std::endl; - _sampleFile = new TFile(mc_files.at(iSample).c_str(), "READ"); + _sampleFile = new TFile(mc_files[iSample].c_str(), "READ"); _data = (TTree*)_sampleFile->Get("caf"); if(_data){ - MACH3LOG_INFO("Found \"caf\" tree in {}", mc_files[iSample]); - MACH3LOG_INFO("With number of entries: {}", _data->GetEntries()); + std::cout << "Found mtuple tree is " << mc_files[iSample] << std::endl; + std::cout << "N of entries: " << _data->GetEntries() << std::endl; } else{ - MACH3LOG_ERROR("Could not find \"caf\" tree in {}", mc_files[iSample]); + MACH3LOG_ERROR("Could not find \"caf\" tree in {}", mc_files[iSample].native()); throw MaCh3Exception(__FILE__, __LINE__); } @@ -242,66 +241,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; @@ -316,41 +315,22 @@ int samplePDFDUNEBeamND::setupExperimentMC(int iSample) { return duneobj->nEvents; } - -const double* samplePDFDUNEBeamND::GetPointerToKinematicParameter(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: MACH3LOG_ERROR("Did not recognise Kinematic Parameter type..."); throw MaCh3Exception(__FILE__, __LINE__); } - return KinematicValue; } -const double* samplePDFDUNEBeamND::GetPointerToKinematicParameter(double KinematicVariable, int iSample, int iEvent) { - KinematicTypes KinPar = (KinematicTypes) std::round(KinematicVariable); - return GetPointerToKinematicParameter(KinPar,iSample,iEvent); -} - -const double* samplePDFDUNEBeamND::GetPointerToKinematicParameter(std::string KinematicParameter, int iSample, int iEvent) { - KinematicTypes KinPar = static_cast(ReturnKinematicParameterFromString(KinematicParameter)); - return GetPointerToKinematicParameter(KinPar,iSample,iEvent); -} - -double samplePDFDUNEBeamND::ReturnKinematicParameter(double KinematicVariable, int iSample, int iEvent) { - return *GetPointerToKinematicParameter(KinematicVariable, iSample, iEvent); -} - -double samplePDFDUNEBeamND::ReturnKinematicParameter(std::string KinematicParameter, int iSample, int iEvent) { - return *GetPointerToKinematicParameter(KinematicParameter, iSample, iEvent); +double samplePDFDUNEBeamND::ReturnKinematicParameter(int KinematicParameter, int iSample, int iEvent) { + return ReturnKinematicParameterByReference(KinematicParameter, iSample, iEvent); } void samplePDFDUNEBeamND::setupFDMC(int iSample) { @@ -434,7 +414,8 @@ void samplePDFDUNEBeamND::applyShifts(int iSample, int iEvent) { } -std::vector samplePDFDUNEBeamND::ReturnKinematicParameterBinning(std::string KinematicParameterStr) + +std::vector samplePDFDUNEBeamND::ReturnKinematicParameterBinning(int KinematicParameter) { std::vector binningVector; return binningVector; diff --git a/samplePDFDUNE/samplePDFDUNEBeamND.h b/samplePDFDUNE/samplePDFDUNEBeamND.h index 3699aed..6683145 100644 --- a/samplePDFDUNE/samplePDFDUNEBeamND.h +++ b/samplePDFDUNE/samplePDFDUNEBeamND.h @@ -35,14 +35,11 @@ class samplePDFDUNEBeamND : virtual public samplePDFFDBase void SetupWeightPointers(); void SetupSplines(); - const double* GetPointerToKinematicParameter(KinematicTypes KinPar, int iSample, int iEvent); - const double* GetPointerToKinematicParameter(double KinematicVariable, int iSample, int iEvent); - const double* GetPointerToKinematicParameter(std::string KinematicParameter, int iSample, int iEvent); + double const &ReturnKinematicParameterByReference(int KinematicParameter, int iSample, int iEvent); + double ReturnKinematicParameter(int 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(int KinematicParameter); - std::vector ReturnKinematicParameterBinning(std::string KinematicParameter); int ReturnKinematicParameterFromString(std::string KinematicParameterStr); std::string ReturnStringFromKinematicParameter(int KinematicParameter); @@ -50,7 +47,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 98f1412..b375f65 100644 --- a/samplePDFDUNE/samplePDFDUNEBeamNDGar.cpp +++ b/samplePDFDUNE/samplePDFDUNEBeamNDGar.cpp @@ -37,7 +37,6 @@ void samplePDFDUNEBeamNDGar::SetupSplines() { return; } - void samplePDFDUNEBeamNDGar::SetupWeightPointers() { for (int i = 0; i < (int)dunendgarmcSamples.size(); ++i) { for (int j = 0; j < dunendgarmcSamples[i].nEvents; ++j) { @@ -54,24 +53,23 @@ void samplePDFDUNEBeamNDGar::SetupWeightPointers() { } int samplePDFDUNEBeamNDGar::setupExperimentMC(int iSample) { - dunemc_base *duneobj = &(dunendgarmcSamples[iSample]); int nutype = sample_nutype[iSample]; int oscnutype = sample_oscnutype[iSample]; bool signal = sample_signal[iSample]; MACH3LOG_INFO("-------------------------------------------------------------------"); - MACH3LOG_INFO("Input File: {}", mc_files.at(iSample)); + MACH3LOG_INFO("Input File: {}", mc_files.at(iSample).native()); _sampleFile = new TFile(mc_files.at(iSample).c_str(), "READ"); _data = (TTree*)_sampleFile->Get("cafTree"); if(_data){ - MACH3LOG_INFO("Found \"caf\" tree in {}", mc_files[iSample]); + MACH3LOG_INFO("Found \"caf\" tree in {}", mc_files[iSample].native()); MACH3LOG_INFO("With number of entries: {}", _data->GetEntries()); } else{ - MACH3LOG_ERROR("Could not find \"caf\" tree in {}", mc_files[iSample]); + MACH3LOG_ERROR("Could not find \"caf\" tree in {}", mc_files[iSample].native()); throw MaCh3Exception(__FILE__, __LINE__); } @@ -87,47 +85,47 @@ int samplePDFDUNEBeamNDGar::setupExperimentMC(int iSample) { duneobj->signal = signal; // 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; @@ -249,79 +247,48 @@ int samplePDFDUNEBeamNDGar::setupExperimentMC(int iSample) { return duneobj->nEvents; } -const double* samplePDFDUNEBeamNDGar::GetPointerToKinematicParameter(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: MACH3LOG_ERROR("Did not recognise Kinematic Parameter type..."); throw MaCh3Exception(__FILE__, __LINE__); } - return KinematicValue; -} - -const double* samplePDFDUNEBeamNDGar::GetPointerToKinematicParameter(double KinematicVariable, int iSample, int iEvent) { - KinematicTypes KinPar = (KinematicTypes) std::round(KinematicVariable); - return GetPointerToKinematicParameter(KinPar,iSample,iEvent); -} - -const double* samplePDFDUNEBeamNDGar::GetPointerToKinematicParameter(std::string KinematicParameter, int iSample, int iEvent) { - KinematicTypes KinPar = static_cast(ReturnKinematicParameterFromString(KinematicParameter)); - return GetPointerToKinematicParameter(KinPar,iSample,iEvent); } -double samplePDFDUNEBeamNDGar::ReturnKinematicParameter(double KinematicVariable, int iSample, int iEvent) { - return *GetPointerToKinematicParameter(KinematicVariable, iSample, iEvent); -} - -double samplePDFDUNEBeamNDGar::ReturnKinematicParameter(std::string KinematicParameter, int iSample, int iEvent) { - return *GetPointerToKinematicParameter(KinematicParameter, iSample, iEvent); +double samplePDFDUNEBeamNDGar::ReturnKinematicParameter(int KinematicParameter, int iSample, int iEvent) { + return ReturnKinematicParameterByReference(KinematicParameter,iSample,iEvent); } void samplePDFDUNEBeamNDGar::setupFDMC(int iSample) { @@ -342,14 +309,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 d9914d8..ebefb0c 100644 --- a/samplePDFDUNE/samplePDFDUNEBeamNDGar.h +++ b/samplePDFDUNE/samplePDFDUNEBeamNDGar.h @@ -40,20 +40,16 @@ class samplePDFDUNEBeamNDGar : virtual public samplePDFFDBase double CalcXsecWeightFunc(int iSample, int iEvent) {return 1.;} void applyShifts(int iSample, int iEvent) {} - const double* GetPointerToKinematicParameter(KinematicTypes KinPar, int iSample, int iEvent); - const double* GetPointerToKinematicParameter(double KinematicVariable, int iSample, int iEvent); - const double* GetPointerToKinematicParameter(std::string KinematicParameter, int iSample, int iEvent); + double const& ReturnKinematicParameterByReference(int KinematicParameter, int iSample, int iEvent); + double ReturnKinematicParameter(int 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(int KinematicParameter); - std::vector ReturnKinematicParameterBinning(std::string KinematicParameter); - std::vector ReturnKinematicParameterBinning(KinematicTypes KinematicParameter); int ReturnKinematicParameterFromString(std::string KinematicParameterStr); std::string ReturnStringFromKinematicParameter(int KinematicParameter); // dunendmc - std::vector dunendgarmcSamples; + std::vector dunendgarmcSamples; TFile *_sampleFile; TTree *_data; 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; }; diff --git a/src/EventRates.cpp b/src/EventRates.cpp index 91c64b0..42f03d0 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,55 +31,93 @@ void Write1DHistogramsToFile(std::string OutFileName, std::vector Histogr return; } -void Write1DHistogramsToPdf(std::string OutFileName, std::vector Histograms) { - - // Now write out the saved hsitograms to file +void Write1DHistogramsToPdf(std::string OutFileName, + std::vector Histograms) { + // 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){ - MACH3LOG_ERROR("Usage: bin/EventRatesDUNEBeam config.cfg"); +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; + covarianceXsec *xsec = nullptr; + covarianceOsc *osc = nullptr; + + // #################################################################################### + // Create samplePDFFD objects - //#################################################################################### - //Create samplePDFFD objects - - std::vector DUNEPdfs; - MakeMaCh3DuneInstance(fitMan.get(), DUNEPdfs, xsec, osc); + std::vector DUNEPdfs; + MakeMaCh3DuneInstance(fitMan.get(), DUNEPdfs, xsec, osc); - std::vector DUNEHists; - for(auto Sample : DUNEPdfs){ + 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()) { + + 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->Draw("COLZ"); + gc1->Print("GenericBinTest.pdf"); + + 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"); + } + } + } + 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); }