diff --git a/data/T2K/CC0pi/STV/multidif_binMap.txt b/data/T2K/CC0pi/STV/multidif_binMap.txt new file mode 100644 index 0000000..59671e0 --- /dev/null +++ b/data/T2K/CC0pi/STV/multidif_binMap.txt @@ -0,0 +1,132 @@ +trueNp +0 + trueCThmu + -1 -0.3 : 0 + -0.3 0.3 + truePmu + 0 0.3 : 1 + 0.3 0.4 : 2 + 0.4 30 : 3 + 0.3 0.6 + truePmu + 0 0.3 : 4 + 0.3 0.4 : 5 + 0.4 0.5 : 6 + 0.5 0.6 : 7 + 0.6 30 : 8 + 0.6 0.7 + truePmu + 0 0.3 : 9 + 0.3 0.4 : 10 + 0.4 0.5 : 11 + 0.5 0.6 : 12 + 0.6 30 : 13 + 0.7 0.8 + truePmu + 0 0.3 : 14 + 0.3 0.4 : 15 + 0.4 0.5 : 16 + 0.5 0.6 : 17 + 0.6 0.7 : 18 + 0.7 0.8 : 19 + 0.8 30 : 20 + 0.8 0.85 + truePmu + 0 0.4 : 21 + 0.4 0.5 : 22 + 0.5 0.6 : 23 + 0.6 0.7 : 24 + 0.7 0.8 : 25 + 0.8 30 : 26 + 0.85 0.9 + truePmu + 0 0.3 : 27 + 0.3 0.4 : 28 + 0.4 0.5 : 29 + 0.5 0.6 : 30 + 0.6 0.7 : 31 + 0.7 0.8 : 32 + 0.8 1 : 33 + 1 30 : 34 + 0.9 0.94 + truePmu + 0 0.4 : 35 + 0.4 0.5 : 36 + 0.5 0.6 : 37 + 0.6 0.7 : 38 + 0.7 0.8 : 39 + 0.8 1.25 : 40 + 1.25 30 : 41 + 0.94 0.98 + truePmu + 0 0.4 : 42 + 0.4 0.5 : 43 + 0.5 0.6 : 44 + 0.6 0.7 : 45 + 0.7 0.8 : 46 + 0.8 1 : 47 + 1 1.25 : 48 + 1.25 1.5 : 49 + 1.5 2 : 50 + 2 30 : 51 + 0.98 1 + truePmu + 0 0.5 : 52 + 0.5 0.65 : 53 + 0.65 0.8 : 54 + 0.8 1.25 : 55 + 1.25 2 : 56 + 2 3 : 57 + 3 5 : 58 + 5 30 : 59 +1 + trueCThmu + -1 -0.3 + trueCThp + -1 0.87 : 60 + 0.87 0.94 : 61 + 0.94 0.97 : 62 + 0.97 1 : 63 + -0.3 0.3 + trueCThp + -1 0.75 : 64 + 0.75 0.85 : 65 + 0.85 0.94 + truePp + 0.5 0.68 : 66 + 0.68 0.78 : 67 + 0.78 0.9 : 68 + 0.9 30 : 69 + 0.94 1 : 70 + 0.3 0.8 + trueCThp + -1 0.3 : 71 + 0.3 0.5 : 72 + 0.5 0.8 + truePp + 0.5 0.6 : 73 + 0.6 0.7 : 74 + 0.7 0.8 : 75 + 0.8 0.9 : 76 + 0.9 30 : 77 + 0.8 1 + truePp + 0.5 0.6 : 78 + 0.6 0.7 : 79 + 0.7 0.8 : 80 + 0.8 1 : 81 + 1 30 : 82 + 0.8 1 + trueCThp + -1 0 : 83 + 0 0.3 : 84 + 0.3 0.8 + truePp + 0.5 0.6 : 85 + 0.6 0.7 : 86 + 0.7 0.8 : 87 + 0.8 0.9 : 88 + 0.9 1.1 : 89 + 1.1 30 : 90 + 0.8 1 : 91 +2 : 92 \ No newline at end of file diff --git a/data/T2K/CC0pi/STV/multidif_results.root b/data/T2K/CC0pi/STV/multidif_results.root new file mode 100644 index 0000000..79eed7c Binary files /dev/null and b/data/T2K/CC0pi/STV/multidif_results.root differ diff --git a/src/FCN/SampleList.cxx b/src/FCN/SampleList.cxx index 76b2584..46b6313 100644 --- a/src/FCN/SampleList.cxx +++ b/src/FCN/SampleList.cxx @@ -1,1475 +1,1477 @@ #include "SampleList.h" #ifndef __NO_ANL__ #include "ANL_CCQE_Evt_1DQ2_nu.h" #include "ANL_CCQE_XSec_1DEnu_nu.h" // ANL CC1ppip #include "ANL_CC1ppip_Evt_1DQ2_nu.h" #include "ANL_CC1ppip_Evt_1DcosmuStar_nu.h" #include "ANL_CC1ppip_Evt_1DcosthAdler_nu.h" #include "ANL_CC1ppip_Evt_1Dphi_nu.h" #include "ANL_CC1ppip_Evt_1Dppi_nu.h" #include "ANL_CC1ppip_Evt_1Dthpr_nu.h" #include "ANL_CC1ppip_XSec_1DEnu_nu.h" #include "ANL_CC1ppip_XSec_1DQ2_nu.h" // ANL CC1npip #include "ANL_CC1npip_Evt_1DQ2_nu.h" #include "ANL_CC1npip_Evt_1DcosmuStar_nu.h" #include "ANL_CC1npip_Evt_1Dppi_nu.h" #include "ANL_CC1npip_XSec_1DEnu_nu.h" // ANL CC1pi0 #include "ANL_CC1pi0_Evt_1DQ2_nu.h" #include "ANL_CC1pi0_Evt_1DcosmuStar_nu.h" #include "ANL_CC1pi0_XSec_1DEnu_nu.h" // ANL NC1npip (mm, exotic!) #include "ANL_NC1npip_Evt_1Dppi_nu.h" // ANL NC1ppim (mm, exotic!) #include "ANL_NC1ppim_Evt_1DcosmuStar_nu.h" #include "ANL_NC1ppim_XSec_1DEnu_nu.h" // ANL CC2pi 1pim1pip (mm, even more exotic!) #include "ANL_CC2pi_1pim1pip_Evt_1Dpmu_nu.h" #include "ANL_CC2pi_1pim1pip_Evt_1Dppim_nu.h" #include "ANL_CC2pi_1pim1pip_Evt_1Dppip_nu.h" #include "ANL_CC2pi_1pim1pip_Evt_1Dpprot_nu.h" #include "ANL_CC2pi_1pim1pip_XSec_1DEnu_nu.h" // ANL CC2pi 1pip1pip (mm, even more exotic!) #include "ANL_CC2pi_1pip1pip_Evt_1Dpmu_nu.h" #include "ANL_CC2pi_1pip1pip_Evt_1Dpneut_nu.h" #include "ANL_CC2pi_1pip1pip_Evt_1DppipHigh_nu.h" #include "ANL_CC2pi_1pip1pip_Evt_1DppipLow_nu.h" #include "ANL_CC2pi_1pip1pip_XSec_1DEnu_nu.h" // ANL CC2pi 1pip1pi0 (mm, even more exotic!) #include "ANL_CC2pi_1pip1pi0_Evt_1Dpmu_nu.h" #include "ANL_CC2pi_1pip1pi0_Evt_1Dppi0_nu.h" #include "ANL_CC2pi_1pip1pi0_Evt_1Dppip_nu.h" #include "ANL_CC2pi_1pip1pi0_Evt_1Dpprot_nu.h" #include "ANL_CC2pi_1pip1pi0_XSec_1DEnu_nu.h" #endif #ifndef __NO_ArgoNeuT__ // ArgoNeuT CC1Pi #include "ArgoNeuT_CC1Pi_XSec_1Dpmu_antinu.h" #include "ArgoNeuT_CC1Pi_XSec_1Dpmu_nu.h" #include "ArgoNeuT_CC1Pi_XSec_1Dthetamu_antinu.h" #include "ArgoNeuT_CC1Pi_XSec_1Dthetamu_nu.h" #include "ArgoNeuT_CC1Pi_XSec_1Dthetamupi_antinu.h" #include "ArgoNeuT_CC1Pi_XSec_1Dthetamupi_nu.h" #include "ArgoNeuT_CC1Pi_XSec_1Dthetapi_antinu.h" #include "ArgoNeuT_CC1Pi_XSec_1Dthetapi_nu.h" // ArgoNeuT CC-inclusive #include "ArgoNeuT_CCInc_XSec_1Dpmu_antinu.h" #include "ArgoNeuT_CCInc_XSec_1Dpmu_nu.h" #include "ArgoNeuT_CCInc_XSec_1Dthetamu_antinu.h" #include "ArgoNeuT_CCInc_XSec_1Dthetamu_nu.h" #endif #ifndef __NO_BNL__ // BNL CCQE #include "BNL_CCQE_Evt_1DQ2_nu.h" #include "BNL_CCQE_XSec_1DEnu_nu.h" // BNL CC1ppip #include "BNL_CC1ppip_Evt_1DQ2_nu.h" #include "BNL_CC1ppip_Evt_1DcosthAdler_nu.h" #include "BNL_CC1ppip_Evt_1Dphi_nu.h" #include "BNL_CC1ppip_XSec_1DEnu_nu.h" // BNL CC1npip #include "BNL_CC1npip_Evt_1DQ2_nu.h" #include "BNL_CC1npip_XSec_1DEnu_nu.h" // BNL CC1pi0 #include "BNL_CC1pi0_Evt_1DQ2_nu.h" #include "BNL_CC1pi0_XSec_1DEnu_nu.h" #endif #ifndef __NO_FNAL__ // FNAL CCQE #include "FNAL_CCQE_Evt_1DQ2_nu.h" // FNAL CC1ppip #include "FNAL_CC1ppip_Evt_1DQ2_nu.h" #include "FNAL_CC1ppip_XSec_1DEnu_nu.h" #include "FNAL_CC1ppip_XSec_1DQ2_nu.h" // FNAL CC1ppim #include "FNAL_CC1ppim_XSec_1DEnu_antinu.h" #endif #ifndef __NO_BEBC__ // BEBC CCQE #include "BEBC_CCQE_XSec_1DQ2_nu.h" // BEBC CC1ppip #include "BEBC_CC1ppip_XSec_1DEnu_nu.h" #include "BEBC_CC1ppip_XSec_1DQ2_nu.h" // BEBC CC1npip #include "BEBC_CC1npip_XSec_1DEnu_nu.h" #include "BEBC_CC1npip_XSec_1DQ2_nu.h" // BEBC CC1pi0 #include "BEBC_CC1pi0_XSec_1DEnu_nu.h" #include "BEBC_CC1pi0_XSec_1DQ2_nu.h" // BEBC CC1npim #include "BEBC_CC1npim_XSec_1DEnu_antinu.h" #include "BEBC_CC1npim_XSec_1DQ2_antinu.h" // BEBC CC1ppim #include "BEBC_CC1ppim_XSec_1DEnu_antinu.h" #include "BEBC_CC1ppim_XSec_1DQ2_antinu.h" #endif #ifndef __NO_GGM__ // GGM CC1ppip #include "GGM_CC1ppip_Evt_1DQ2_nu.h" #include "GGM_CC1ppip_XSec_1DEnu_nu.h" #endif #ifndef __NO_MiniBooNE__ // MiniBooNE CCQE #include "MiniBooNE_CCQE_XSec_1DEnu_nu.h" #include "MiniBooNE_CCQE_XSec_1DQ2_antinu.h" #include "MiniBooNE_CCQE_XSec_1DQ2_nu.h" #include "MiniBooNE_CCQE_XSec_2DTcos_antinu.h" #include "MiniBooNE_CCQE_XSec_2DTcos_nu.h" // MiniBooNE CC1pi+ 1D #include "MiniBooNE_CC1pip_XSec_1DEnu_nu.h" #include "MiniBooNE_CC1pip_XSec_1DQ2_nu.h" #include "MiniBooNE_CC1pip_XSec_1DTpi_nu.h" #include "MiniBooNE_CC1pip_XSec_1DTu_nu.h" // MiniBooNE CC1pi+ 2D #include "MiniBooNE_CC1pip_XSec_2DQ2Enu_nu.h" #include "MiniBooNE_CC1pip_XSec_2DTpiCospi_nu.h" #include "MiniBooNE_CC1pip_XSec_2DTpiEnu_nu.h" #include "MiniBooNE_CC1pip_XSec_2DTuCosmu_nu.h" #include "MiniBooNE_CC1pip_XSec_2DTuEnu_nu.h" // MiniBooNE CC1pi0 #include "MiniBooNE_CC1pi0_XSec_1DEnu_nu.h" #include "MiniBooNE_CC1pi0_XSec_1DQ2_nu.h" #include "MiniBooNE_CC1pi0_XSec_1DTu_nu.h" #include "MiniBooNE_CC1pi0_XSec_1Dcosmu_nu.h" #include "MiniBooNE_CC1pi0_XSec_1Dcospi0_nu.h" #include "MiniBooNE_CC1pi0_XSec_1Dppi0_nu.h" #include "MiniBooNE_NC1pi0_XSec_1Dcospi0_antinu.h" #include "MiniBooNE_NC1pi0_XSec_1Dcospi0_nu.h" #include "MiniBooNE_NC1pi0_XSec_1Dppi0_antinu.h" #include "MiniBooNE_NC1pi0_XSec_1Dppi0_nu.h" // MiniBooNE NC1pi0 //#include "MiniBooNE_NCpi0_XSec_1Dppi0_nu.h" // MiniBooNE NCEL #include "MiniBooNE_NCEL_XSec_Treco_nu.h" #endif #ifndef __NO_MicroBooNE__ #include "MicroBooNE_CCInc_XSec_2DPcos_nu.h" #endif #ifndef __NO_MINERvA__ // MINERvA CCQE #include "MINERvA_CCQE_XSec_1DQ2_antinu.h" #include "MINERvA_CCQE_XSec_1DQ2_joint.h" #include "MINERvA_CCQE_XSec_1DQ2_nu.h" // MINERvA CC0pi #include "MINERvA_CC0pi_XSec_1DEe_nue.h" #include "MINERvA_CC0pi_XSec_1DQ2_nu_proton.h" #include "MINERvA_CC0pi_XSec_1DQ2_nue.h" #include "MINERvA_CC0pi_XSec_1DThetae_nue.h" // 2018 MINERvA CC0pi STV #include "MINERvA_CC0pinp_STV_XSec_1D_nu.h" // 2018 MINERvA CC0pi 2D #include "MINERvA_CC0pi_XSec_1D_2018_nu.h" #include "MINERvA_CC0pi_XSec_2D_nu.h" // #include "MINERvA_CC0pi_XSec_3DptpzTp_nu.h" // 2018 MINERvA CC0pi 2D antinu #include "MINERvA_CC0pi_XSec_2D_antinu.h" // MINERvA CC1pi+ #include "MINERvA_CC1pip_XSec_1DTpi_20deg_nu.h" #include "MINERvA_CC1pip_XSec_1DTpi_nu.h" #include "MINERvA_CC1pip_XSec_1Dth_20deg_nu.h" #include "MINERvA_CC1pip_XSec_1Dth_nu.h" // 2017 data update #include "MINERvA_CC1pip_XSec_1D_2017Update.h" // MINERvA CCNpi+ #include "MINERvA_CCNpip_XSec_1DEnu_nu.h" #include "MINERvA_CCNpip_XSec_1DQ2_nu.h" #include "MINERvA_CCNpip_XSec_1DTpi_nu.h" #include "MINERvA_CCNpip_XSec_1Dpmu_nu.h" #include "MINERvA_CCNpip_XSec_1Dth_nu.h" #include "MINERvA_CCNpip_XSec_1Dthmu_nu.h" // MINERvA CC1pi0 #include "MINERvA_CC1pi0_XSec_1DEnu_antinu.h" #include "MINERvA_CC1pi0_XSec_1DQ2_antinu.h" #include "MINERvA_CC1pi0_XSec_1DTpi0_antinu.h" #include "MINERvA_CC1pi0_XSec_1Dpmu_antinu.h" #include "MINERvA_CC1pi0_XSec_1Dppi0_antinu.h" #include "MINERvA_CC1pi0_XSec_1Dth_antinu.h" #include "MINERvA_CC1pi0_XSec_1Dthmu_antinu.h" // MINERvA CC1pi0 neutrino #include "MINERvA_CC1pi0_XSec_1D_nu.h" // MINERvA CCINC #include "MINERvA_CCinc_XSec_1DEnu_ratio.h" #include "MINERvA_CCinc_XSec_1Dx_ratio.h" #include "MINERvA_CCinc_XSec_2DEavq3_nu.h" // MINERvA CCDIS #include "MINERvA_CCDIS_XSec_1DEnu_ratio.h" #include "MINERvA_CCDIS_XSec_1Dx_ratio.h" // MINERvA CCCOH pion #include "MINERvA_CCCOHPI_XSec_1DEnu_antinu.h" #include "MINERvA_CCCOHPI_XSec_1DEpi_antinu.h" #include "MINERvA_CCCOHPI_XSec_1DQ2_antinu.h" #include "MINERvA_CCCOHPI_XSec_1DEpi_nu.h" #include "MINERvA_CCCOHPI_XSec_1DQ2_nu.h" #include "MINERvA_CCCOHPI_XSec_1Dth_nu.h" #include "MINERvA_CCCOHPI_XSec_joint.h" #include "MINERvA_CC0pi_XSec_1DQ2_TgtRatio_nu.h" #include "MINERvA_CC0pi_XSec_1DQ2_Tgt_nu.h" #endif #ifndef __NO_T2K__ // T2K CC0pi 2016 #include "T2K_CC0pi_XSec_2DPcos_nu_I.h" #include "T2K_CC0pi_XSec_2DPcos_nu_II.h" // T2K CC0pi 2020 arXiv:1908.10249 #include "T2K_CC0pi_XSec_H2O_2DPcos_anu.h" // T2K CC0pi 2020 arXiv:2004.05434 #include "T2K_NuMu_CC0pi_OC_XSec_2DPcos.h" #include "T2K_NuMu_CC0pi_OC_XSec_2DPcos_joint.h" // T2K CC0pi 2020 arXiv:2002.09323 #include "T2K_NuMuAntiNuMu_CC0pi_CH_XSec_2DPcos.h" #include "T2K_NuMuAntiNuMu_CC0pi_CH_XSec_2DPcos_joint.h" // T2K CC-inclusive with full acceptance 2018 #include "T2K_CCinc_XSec_2DPcos_nu_nonuniform.h" // T2K nue CC-inclusive 2019 #include "T2K_nueCCinc_XSec_1Dpe.h" #include "T2K_nueCCinc_XSec_1Dthe.h" #include "T2K_nueCCinc_XSec_1Dpe_joint.h" #include "T2K_nueCCinc_XSec_1Dthe_joint.h" #include "T2K_nueCCinc_XSec_joint.h" // T2K STV CC0pi 2018 -#include "T2K_CC0pi1p_XSec_3DPcoscos_nu_nonuniform.h" +#include "T2K_CC0piWithProtons_XSec_2018_multidif_0p_1p_Np.h" #include "T2K_CC0pinp_STV_XSec_1Ddat_nu.h" #include "T2K_CC0pinp_STV_XSec_1Ddphit_nu.h" #include "T2K_CC0pinp_STV_XSec_1Ddpt_nu.h" #include "T2K_CC0pinp_ifk_XSec_3Dinfa_nu.h" #include "T2K_CC0pinp_ifk_XSec_3Dinfip_nu.h" #include "T2K_CC0pinp_ifk_XSec_3Dinfp_nu.h" // T2K CC1pi+ on CH #include "T2K_CC1pip_CH_XSec_1DAdlerPhi_nu.h" #include "T2K_CC1pip_CH_XSec_1DCosThAdler_nu.h" #include "T2K_CC1pip_CH_XSec_1DQ2_nu.h" #include "T2K_CC1pip_CH_XSec_1Dppi_nu.h" #include "T2K_CC1pip_CH_XSec_1Dthmupi_nu.h" #include "T2K_CC1pip_CH_XSec_1Dthpi_nu.h" #include "T2K_CC1pip_CH_XSec_2Dpmucosmu_nu.h" // T2K CC1pi+ on H2O #include "T2K_CC1pip_H2O_XSec_1DEnuDelta_nu.h" #include "T2K_CC1pip_H2O_XSec_1DEnuMB_nu.h" #include "T2K_CC1pip_H2O_XSec_1Dcosmu_nu.h" #include "T2K_CC1pip_H2O_XSec_1Dcosmupi_nu.h" #include "T2K_CC1pip_H2O_XSec_1Dcospi_nu.h" #include "T2K_CC1pip_H2O_XSec_1Dpmu_nu.h" #include "T2K_CC1pip_H2O_XSec_1Dppi_nu.h" +// add header here + #endif #ifndef __NO_SciBooNE__ // SciBooNE COH studies #include "SciBooNE_CCCOH_1TRK_1DQ2_nu.h" #include "SciBooNE_CCCOH_1TRK_1Dpmu_nu.h" #include "SciBooNE_CCCOH_1TRK_1Dthetamu_nu.h" #include "SciBooNE_CCCOH_MuPiNoVA_1DQ2_nu.h" #include "SciBooNE_CCCOH_MuPiNoVA_1Dpmu_nu.h" #include "SciBooNE_CCCOH_MuPiNoVA_1Dthetamu_nu.h" #include "SciBooNE_CCCOH_MuPiNoVA_1Dthetapi_nu.h" #include "SciBooNE_CCCOH_MuPiNoVA_1Dthetapr_nu.h" #include "SciBooNE_CCCOH_MuPiVA_1DQ2_nu.h" #include "SciBooNE_CCCOH_MuPiVA_1Dpmu_nu.h" #include "SciBooNE_CCCOH_MuPiVA_1Dthetamu_nu.h" #include "SciBooNE_CCCOH_MuPr_1DQ2_nu.h" #include "SciBooNE_CCCOH_MuPr_1Dpmu_nu.h" #include "SciBooNE_CCCOH_MuPr_1Dthetamu_nu.h" #include "SciBooNE_CCCOH_STOPFINAL_1DQ2_nu.h" #include "SciBooNE_CCCOH_STOP_NTrks_nu.h" #include "SciBooNE_CCInc_XSec_1DEnu_nu.h" #endif #ifndef __NO_K2K__ // K2K NC1pi0 #include "K2K_NC1pi0_Evt_1Dppi0_nu.h" #endif // MC Studies #include "ExpMultDist_CCQE_XSec_1DVar_FakeStudy.h" #include "ExpMultDist_CCQE_XSec_2DVar_FakeStudy.h" #include "MCStudy_CCQEHistograms.h" #include "GenericFlux_Tester.h" #include "GenericFlux_Vectors.h" #include "ElectronFlux_FlatTree.h" #include "ElectronScattering_DurhamData.h" #include "MCStudy_KaonPreSelection.h" #include "MCStudy_MuonValidation.h" #include "OfficialNIWGPlots.h" #include "T2K2017_FakeData.h" #include "SigmaEnuHists.h" #include "Simple_Osc.h" #include "Smear_SVDUnfold_Propagation_Osc.h" #include "FitWeight.h" #include "NuisConfig.h" #include "NuisKey.h" #ifdef __USE_DYNSAMPLES__ #include "TRegexp.h" #include // linux #include DynamicSampleFactory::DynamicSampleFactory() : NSamples(0), NManifests(0) { LoadPlugins(); NUIS_LOG(FIT, "Loaded " << NSamples << " from " << NManifests << " shared object libraries."); } DynamicSampleFactory *DynamicSampleFactory::glblDSF = NULL; DynamicSampleFactory::PluginManifest::~PluginManifest() { for (size_t i_it = 0; i_it < Instances.size(); ++i_it) { (*(DSF_DestroySample))(Instances[i_it]); } } std::string EnsureTrailingSlash(std::string const &inp) { if (!inp.length()) { return "/"; } if (inp[inp.length() - 1] == '/') { return inp; } return inp + "/"; } void DynamicSampleFactory::LoadPlugins() { std::vector SearchDirectories; if (Config::HasPar("dynamic_sample.path")) { SearchDirectories = GeneralUtils::ParseToStr(Config::GetParS("dynamic_sample.path"), ":"); } char const *envPath = getenv("NUISANCE_DS_PATH"); if (envPath) { std::vector envPaths = GeneralUtils::ParseToStr(envPath, ":"); for (size_t ep_it = 0; ep_it < envPaths.size(); ++ep_it) { SearchDirectories.push_back(envPaths[ep_it]); } } if (!SearchDirectories.size()) { char const *pwdPath = getenv("PWD"); if (pwdPath) { SearchDirectories.push_back(pwdPath); } } for (size_t sp_it = 0; sp_it < SearchDirectories.size(); ++sp_it) { std::string dirpath = EnsureTrailingSlash(SearchDirectories[sp_it]); NUIS_LOG(FIT, "Searching for dynamic sample manifests in: " << dirpath); Ssiz_t len = 0; DIR *dir; struct dirent *ent; dir = opendir(dirpath.c_str()); if (dir != NULL) { TRegexp matchExp("*.so", true); while ((ent = readdir(dir)) != NULL) { if (matchExp.Index(TString(ent->d_name), &len) != Ssiz_t(-1)) { NUIS_LOG(FIT, "\tFound shared object: " << ent->d_name << " checking for relevant methods..."); void *dlobj = dlopen((dirpath + ent->d_name).c_str(), RTLD_NOW | RTLD_GLOBAL); char const *dlerr_cstr = dlerror(); std::string dlerr; if (dlerr_cstr) { dlerr = dlerr_cstr; } if (dlerr.length()) { NUIS_ERR(WRN, "\tDL Load Error: " << dlerr); continue; } PluginManifest plgManif; plgManif.dllib = dlobj; plgManif.soloc = (dirpath + ent->d_name); plgManif.DSF_NSamples = reinterpret_cast(dlsym(dlobj, "DSF_NSamples")); dlerr = ""; dlerr_cstr = dlerror(); if (dlerr_cstr) { dlerr = dlerr_cstr; } if (dlerr.length()) { NUIS_ERR(WRN, "\tFailed to load symbol \"DSF_NSamples\" from " << (dirpath + ent->d_name) << ": " << dlerr); dlclose(dlobj); continue; } plgManif.DSF_GetSampleName = reinterpret_cast( dlsym(dlobj, "DSF_GetSampleName")); dlerr = ""; dlerr_cstr = dlerror(); if (dlerr_cstr) { dlerr = dlerr_cstr; } if (dlerr.length()) { NUIS_ERR(WRN, "\tFailed to load symbol \"DSF_GetSampleName\" from " << (dirpath + ent->d_name) << ": " << dlerr); dlclose(dlobj); continue; } plgManif.DSF_GetSample = reinterpret_cast( dlsym(dlobj, "DSF_GetSample")); dlerr = ""; dlerr_cstr = dlerror(); if (dlerr_cstr) { dlerr = dlerr_cstr; } if (dlerr.length()) { NUIS_ERR(WRN, "\tFailed to load symbol \"DSF_GetSample\" from " << (dirpath + ent->d_name) << ": " << dlerr); dlclose(dlobj); continue; } plgManif.DSF_DestroySample = reinterpret_cast( dlsym(dlobj, "DSF_DestroySample")); dlerr = ""; dlerr_cstr = dlerror(); if (dlerr_cstr) { dlerr = dlerr_cstr; } if (dlerr.length()) { NUIS_ERR(WRN, "Failed to load symbol \"DSF_DestroySample\" from " << (dirpath + ent->d_name) << ": " << dlerr); dlclose(dlobj); continue; } plgManif.NSamples = (*(plgManif.DSF_NSamples))(); NUIS_LOG(FIT, "\tSuccessfully loaded dynamic sample manifest: " << plgManif.soloc << ". Contains " << plgManif.NSamples << " samples."); for (size_t smp_it = 0; smp_it < plgManif.NSamples; ++smp_it) { char const *smp_name = (*(plgManif.DSF_GetSampleName))(smp_it); if (!smp_name) { NUIS_ABORT("Could not load sample " << smp_it << " / " << plgManif.NSamples << " from " << plgManif.soloc); } if (Samples.count(smp_name)) { NUIS_ERR(WRN, "Already loaded a sample named: \"" << smp_name << "\". cannot load duplciates. This " "sample will be skipped."); continue; } plgManif.SamplesProvided.push_back(smp_name); Samples[smp_name] = std::make_pair(plgManif.soloc, smp_it); NUIS_LOG(FIT, "\t\t" << smp_name); } if (plgManif.SamplesProvided.size()) { Manifests[plgManif.soloc] = plgManif; NSamples += plgManif.SamplesProvided.size(); NManifests++; } else { dlclose(dlobj); } } } closedir(dir); } else { NUIS_ERR(WRN, "Tried to open non-existant directory."); } } } DynamicSampleFactory &DynamicSampleFactory::Get() { if (!glblDSF) { glblDSF = new DynamicSampleFactory(); } return *glblDSF; } void DynamicSampleFactory::Print() { std::map > ManifestSamples; for (std::map >::iterator smp_it = Samples.begin(); smp_it != Samples.end(); ++smp_it) { if (!ManifestSamples.count(smp_it->second.first)) { ManifestSamples[smp_it->second.first] = std::vector(); } ManifestSamples[smp_it->second.first].push_back(smp_it->first); } NUIS_LOG(FIT, "Dynamic sample manifest: "); for (std::map >::iterator m_it = ManifestSamples.begin(); m_it != ManifestSamples.end(); ++m_it) { NUIS_LOG(FIT, "\tLibrary " << m_it->first << " contains: "); for (size_t s_it = 0; s_it < m_it->second.size(); ++s_it) { NUIS_LOG(FIT, "\t\t" << m_it->second[s_it]); } } } bool DynamicSampleFactory::HasSample(std::string const &name) { return Samples.count(name); } bool DynamicSampleFactory::HasSample(nuiskey &samplekey) { return HasSample(samplekey.GetS("name")); } MeasurementBase *DynamicSampleFactory::CreateSample(nuiskey &samplekey) { if (!HasSample(samplekey)) { NUIS_ERR(WRN, "Asked to load unknown sample: \"" << samplekey.GetS("name") << "\"."); return NULL; } std::pair sample = Samples[samplekey.GetS("name")]; NUIS_LOG(SAM, "\tLoading sample " << sample.second << " from " << sample.first); return (*(Manifests[sample.first].DSF_GetSample))(sample.second, &samplekey); } DynamicSampleFactory::~DynamicSampleFactory() { Manifests.clear(); } #endif //! Functions to make it easier for samples to be created and handled. namespace SampleUtils { //! Create a given sample given its name, file, type, fakdata(fkdt) file and the //! current rw engine and push it back into the list fChain. MeasurementBase *CreateSample(std::string name, std::string file, std::string type, std::string fkdt, FitWeight *rw) { nuiskey samplekey = Config::CreateKey("sample"); samplekey.Set("name", name); samplekey.Set("input", file); samplekey.Set("type", type); return CreateSample(samplekey); } MeasurementBase *CreateSample(nuiskey samplekey) { #ifdef __USE_DYNSAMPLES__ if (DynamicSampleFactory::Get().HasSample(samplekey)) { NUIS_LOG(SAM, "Instantiating dynamic sample..."); MeasurementBase *ds = DynamicSampleFactory::Get().CreateSample(samplekey); if (ds) { NUIS_LOG(SAM, "Done."); return ds; } NUIS_ABORT("Failed to instantiate dynamic sample."); } #endif FitWeight *rw = FitBase::GetRW(); std::string name = samplekey.GetS("name"); std::string file = samplekey.GetS("input"); std::string type = samplekey.GetS("type"); std::string fkdt = ""; /* ANL CCQE Samples */ #ifndef __NO_ANL__ if (!name.compare("ANL_CCQE_XSec_1DEnu_nu") || !name.compare("ANL_CCQE_XSec_1DEnu_nu_PRD26") || !name.compare("ANL_CCQE_XSec_1DEnu_nu_PRL31") || !name.compare("ANL_CCQE_XSec_1DEnu_nu_PRD16")) { return (new ANL_CCQE_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("ANL_CCQE_Evt_1DQ2_nu") || !name.compare("ANL_CCQE_Evt_1DQ2_nu_PRL31") || !name.compare("ANL_CCQE_Evt_1DQ2_nu_PRD26") || !name.compare("ANL_CCQE_Evt_1DQ2_nu_PRD16")) { return (new ANL_CCQE_Evt_1DQ2_nu(samplekey)); /* ANL CC1ppip samples */ } else if (!name.compare("ANL_CC1ppip_XSec_1DEnu_nu") || !name.compare("ANL_CC1ppip_XSec_1DEnu_nu_W14Cut") || !name.compare("ANL_CC1ppip_XSec_1DEnu_nu_Uncorr") || !name.compare("ANL_CC1ppip_XSec_1DEnu_nu_W14Cut_Uncorr") || !name.compare("ANL_CC1ppip_XSec_1DEnu_nu_W16Cut_Uncorr")) { return (new ANL_CC1ppip_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("ANL_CC1ppip_XSec_1DQ2_nu")) { return (new ANL_CC1ppip_XSec_1DQ2_nu(samplekey)); } else if (!name.compare("ANL_CC1ppip_Evt_1DQ2_nu") || !name.compare("ANL_CC1ppip_Evt_1DQ2_nu_W14Cut")) { return (new ANL_CC1ppip_Evt_1DQ2_nu(samplekey)); } else if (!name.compare("ANL_CC1ppip_Evt_1Dppi_nu")) { return (new ANL_CC1ppip_Evt_1Dppi_nu(samplekey)); } else if (!name.compare("ANL_CC1ppip_Evt_1Dthpr_nu")) { return (new ANL_CC1ppip_Evt_1Dthpr_nu(samplekey)); } else if (!name.compare("ANL_CC1ppip_Evt_1DcosmuStar_nu")) { return (new ANL_CC1ppip_Evt_1DcosmuStar_nu(samplekey)); } else if (!name.compare("ANL_CC1ppip_Evt_1DcosthAdler_nu")) { return (new ANL_CC1ppip_Evt_1DcosthAdler_nu(samplekey)); } else if (!name.compare("ANL_CC1ppip_Evt_1Dphi_nu")) { return (new ANL_CC1ppip_Evt_1Dphi_nu(samplekey)); /* ANL CC1npip sample */ } else if (!name.compare("ANL_CC1npip_XSec_1DEnu_nu") || !name.compare("ANL_CC1npip_XSec_1DEnu_nu_W14Cut") || !name.compare("ANL_CC1npip_XSec_1DEnu_nu_Uncorr") || !name.compare("ANL_CC1npip_XSec_1DEnu_nu_W14Cut_Uncorr") || !name.compare("ANL_CC1npip_XSec_1DEnu_nu_W16Cut_Uncorr")) { return (new ANL_CC1npip_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("ANL_CC1npip_Evt_1DQ2_nu") || !name.compare("ANL_CC1npip_Evt_1DQ2_nu_W14Cut")) { return (new ANL_CC1npip_Evt_1DQ2_nu(samplekey)); } else if (!name.compare("ANL_CC1npip_Evt_1Dppi_nu")) { return (new ANL_CC1npip_Evt_1Dppi_nu(samplekey)); } else if (!name.compare("ANL_CC1npip_Evt_1DcosmuStar_nu")) { return (new ANL_CC1npip_Evt_1DcosmuStar_nu(samplekey)); /* ANL CC1pi0 sample */ } else if (!name.compare("ANL_CC1pi0_XSec_1DEnu_nu") || !name.compare("ANL_CC1pi0_XSec_1DEnu_nu_W14Cut") || !name.compare("ANL_CC1pi0_XSec_1DEnu_nu_Uncorr") || !name.compare("ANL_CC1pi0_XSec_1DEnu_nu_W14Cut_Uncorr") || !name.compare("ANL_CC1pi0_XSec_1DEnu_nu_W16Cut_Uncorr")) { return (new ANL_CC1pi0_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("ANL_CC1pi0_Evt_1DQ2_nu") || !name.compare("ANL_CC1pi0_Evt_1DQ2_nu_W14Cut")) { return (new ANL_CC1pi0_Evt_1DQ2_nu(samplekey)); } else if (!name.compare("ANL_CC1pi0_Evt_1DcosmuStar_nu")) { return (new ANL_CC1pi0_Evt_1DcosmuStar_nu(samplekey)); /* ANL NC1npip sample */ } else if (!name.compare("ANL_NC1npip_Evt_1Dppi_nu")) { return (new ANL_NC1npip_Evt_1Dppi_nu(samplekey)); /* ANL NC1ppim sample */ } else if (!name.compare("ANL_NC1ppim_XSec_1DEnu_nu")) { return (new ANL_NC1ppim_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("ANL_NC1ppim_Evt_1DcosmuStar_nu")) { return (new ANL_NC1ppim_Evt_1DcosmuStar_nu(samplekey)); /* ANL CC2pi sample */ } else if (!name.compare("ANL_CC2pi_1pim1pip_XSec_1DEnu_nu")) { return (new ANL_CC2pi_1pim1pip_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("ANL_CC2pi_1pim1pip_Evt_1Dpmu_nu")) { return (new ANL_CC2pi_1pim1pip_Evt_1Dpmu_nu(samplekey)); } else if (!name.compare("ANL_CC2pi_1pim1pip_Evt_1Dppip_nu")) { return (new ANL_CC2pi_1pim1pip_Evt_1Dppip_nu(samplekey)); } else if (!name.compare("ANL_CC2pi_1pim1pip_Evt_1Dppim_nu")) { return (new ANL_CC2pi_1pim1pip_Evt_1Dppim_nu(samplekey)); } else if (!name.compare("ANL_CC2pi_1pim1pip_Evt_1Dpprot_nu")) { return (new ANL_CC2pi_1pim1pip_Evt_1Dpprot_nu(samplekey)); } else if (!name.compare("ANL_CC2pi_1pip1pip_XSec_1DEnu_nu")) { return (new ANL_CC2pi_1pip1pip_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("ANL_CC2pi_1pip1pip_Evt_1Dpmu_nu")) { return (new ANL_CC2pi_1pip1pip_Evt_1Dpmu_nu(samplekey)); } else if (!name.compare("ANL_CC2pi_1pip1pip_Evt_1Dpneut_nu")) { return (new ANL_CC2pi_1pip1pip_Evt_1Dpneut_nu(samplekey)); } else if (!name.compare("ANL_CC2pi_1pip1pip_Evt_1DppipHigh_nu")) { return (new ANL_CC2pi_1pip1pip_Evt_1DppipHigh_nu(samplekey)); } else if (!name.compare("ANL_CC2pi_1pip1pip_Evt_1DppipLow_nu")) { return (new ANL_CC2pi_1pip1pip_Evt_1DppipLow_nu(samplekey)); } else if (!name.compare("ANL_CC2pi_1pip1pi0_XSec_1DEnu_nu")) { return (new ANL_CC2pi_1pip1pi0_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("ANL_CC2pi_1pip1pi0_Evt_1Dpmu_nu")) { return (new ANL_CC2pi_1pip1pi0_Evt_1Dpmu_nu(samplekey)); } else if (!name.compare("ANL_CC2pi_1pip1pi0_Evt_1Dppip_nu")) { return (new ANL_CC2pi_1pip1pi0_Evt_1Dppip_nu(samplekey)); } else if (!name.compare("ANL_CC2pi_1pip1pi0_Evt_1Dppi0_nu")) { return (new ANL_CC2pi_1pip1pi0_Evt_1Dppi0_nu(samplekey)); } else if (!name.compare("ANL_CC2pi_1pip1pi0_Evt_1Dpprot_nu")) { return (new ANL_CC2pi_1pip1pi0_Evt_1Dpprot_nu(samplekey)); /* ArgoNeut Samples */ } else #endif #ifndef __NO_ArgoNeuT__ if (!name.compare("ArgoNeuT_CCInc_XSec_1Dpmu_antinu")) { return (new ArgoNeuT_CCInc_XSec_1Dpmu_antinu(samplekey)); } else if (!name.compare("ArgoNeuT_CCInc_XSec_1Dpmu_nu")) { return (new ArgoNeuT_CCInc_XSec_1Dpmu_nu(samplekey)); } else if (!name.compare("ArgoNeuT_CCInc_XSec_1Dthetamu_antinu")) { return (new ArgoNeuT_CCInc_XSec_1Dthetamu_antinu(samplekey)); } else if (!name.compare("ArgoNeuT_CCInc_XSec_1Dthetamu_nu")) { return (new ArgoNeuT_CCInc_XSec_1Dthetamu_nu(samplekey)); } else if (!name.compare("ArgoNeuT_CC1Pi_XSec_1Dpmu_nu")) { return (new ArgoNeuT_CC1Pi_XSec_1Dpmu_nu(samplekey)); } else if (!name.compare("ArgoNeuT_CC1Pi_XSec_1Dthetamu_nu")) { return (new ArgoNeuT_CC1Pi_XSec_1Dthetamu_nu(samplekey)); } else if (!name.compare("ArgoNeuT_CC1Pi_XSec_1Dthetapi_nu")) { return (new ArgoNeuT_CC1Pi_XSec_1Dthetapi_nu(samplekey)); } else if (!name.compare("ArgoNeuT_CC1Pi_XSec_1Dthetamupi_nu")) { return (new ArgoNeuT_CC1Pi_XSec_1Dthetamupi_nu(samplekey)); } else if (!name.compare("ArgoNeuT_CC1Pi_XSec_1Dpmu_antinu")) { return (new ArgoNeuT_CC1Pi_XSec_1Dpmu_antinu(samplekey)); } else if (!name.compare("ArgoNeuT_CC1Pi_XSec_1Dthetamu_antinu")) { return (new ArgoNeuT_CC1Pi_XSec_1Dthetamu_antinu(samplekey)); } else if (!name.compare("ArgoNeuT_CC1Pi_XSec_1Dthetapi_antinu")) { return (new ArgoNeuT_CC1Pi_XSec_1Dthetapi_antinu(samplekey)); } else if (!name.compare("ArgoNeuT_CC1Pi_XSec_1Dthetamupi_antinu")) { return (new ArgoNeuT_CC1Pi_XSec_1Dthetamupi_antinu(samplekey)); /* BNL Samples */ } else #endif #ifndef __NO_BNL__ if (!name.compare("BNL_CCQE_XSec_1DEnu_nu")) { return (new BNL_CCQE_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("BNL_CCQE_Evt_1DQ2_nu")) { return (new BNL_CCQE_Evt_1DQ2_nu(samplekey)); /* BNL CC1ppip samples */ } else if (!name.compare("BNL_CC1ppip_XSec_1DEnu_nu") || !name.compare("BNL_CC1ppip_XSec_1DEnu_nu_Uncorr") || !name.compare("BNL_CC1ppip_XSec_1DEnu_nu_W14Cut") || !name.compare("BNL_CC1ppip_XSec_1DEnu_nu_W14Cut_Uncorr")) { return (new BNL_CC1ppip_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("BNL_CC1ppip_Evt_1DQ2_nu") || !name.compare("BNL_CC1ppip_Evt_1DQ2_nu_W14Cut")) { return (new BNL_CC1ppip_Evt_1DQ2_nu(samplekey)); } else if (!name.compare("BNL_CC1ppip_Evt_1DcosthAdler_nu")) { return (new BNL_CC1ppip_Evt_1DcosthAdler_nu(samplekey)); } else if (!name.compare("BNL_CC1ppip_Evt_1Dphi_nu")) { return (new BNL_CC1ppip_Evt_1Dphi_nu(samplekey)); /* BNL CC1npip samples */ } else if (!name.compare("BNL_CC1npip_XSec_1DEnu_nu") || !name.compare("BNL_CC1npip_XSec_1DEnu_nu_Uncorr")) { return (new BNL_CC1npip_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("BNL_CC1npip_Evt_1DQ2_nu")) { return (new BNL_CC1npip_Evt_1DQ2_nu(samplekey)); /* BNL CC1pi0 samples */ } else if (!name.compare("BNL_CC1pi0_XSec_1DEnu_nu")) { return (new BNL_CC1pi0_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("BNL_CC1pi0_Evt_1DQ2_nu")) { return (new BNL_CC1pi0_Evt_1DQ2_nu(samplekey)); /* FNAL Samples */ } else #endif #ifndef __NO_FNAL__ if (!name.compare("FNAL_CCQE_Evt_1DQ2_nu")) { return (new FNAL_CCQE_Evt_1DQ2_nu(samplekey)); /* FNAL CC1ppip */ } else if (!name.compare("FNAL_CC1ppip_XSec_1DEnu_nu")) { return (new FNAL_CC1ppip_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("FNAL_CC1ppip_XSec_1DQ2_nu")) { return (new FNAL_CC1ppip_XSec_1DQ2_nu(samplekey)); } else if (!name.compare("FNAL_CC1ppip_Evt_1DQ2_nu")) { return (new FNAL_CC1ppip_Evt_1DQ2_nu(samplekey)); /* FNAL CC1ppim */ } else if (!name.compare("FNAL_CC1ppim_XSec_1DEnu_antinu")) { return (new FNAL_CC1ppim_XSec_1DEnu_antinu(samplekey)); /* BEBC Samples */ } else #endif #ifndef __NO_BEBC__ if (!name.compare("BEBC_CCQE_XSec_1DQ2_nu")) { return (new BEBC_CCQE_XSec_1DQ2_nu(samplekey)); /* BEBC CC1ppip samples */ } else if (!name.compare("BEBC_CC1ppip_XSec_1DEnu_nu")) { return (new BEBC_CC1ppip_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("BEBC_CC1ppip_XSec_1DQ2_nu")) { return (new BEBC_CC1ppip_XSec_1DQ2_nu(samplekey)); /* BEBC CC1npip samples */ } else if (!name.compare("BEBC_CC1npip_XSec_1DEnu_nu")) { return (new BEBC_CC1npip_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("BEBC_CC1npip_XSec_1DQ2_nu")) { return (new BEBC_CC1npip_XSec_1DQ2_nu(samplekey)); /* BEBC CC1pi0 samples */ } else if (!name.compare("BEBC_CC1pi0_XSec_1DEnu_nu")) { return (new BEBC_CC1pi0_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("BEBC_CC1pi0_XSec_1DQ2_nu")) { return (new BEBC_CC1pi0_XSec_1DQ2_nu(samplekey)); /* BEBC CC1npim samples */ } else if (!name.compare("BEBC_CC1npim_XSec_1DEnu_antinu")) { return (new BEBC_CC1npim_XSec_1DEnu_antinu(samplekey)); } else if (!name.compare("BEBC_CC1npim_XSec_1DQ2_antinu")) { return (new BEBC_CC1npim_XSec_1DQ2_antinu(samplekey)); /* BEBC CC1ppim samples */ } else if (!name.compare("BEBC_CC1ppim_XSec_1DEnu_antinu")) { return (new BEBC_CC1ppim_XSec_1DEnu_antinu(samplekey)); } else if (!name.compare("BEBC_CC1ppim_XSec_1DQ2_antinu")) { return (new BEBC_CC1ppim_XSec_1DQ2_antinu(samplekey)); /* GGM CC1ppip samples */ } else #endif #ifndef __NO_GGM__ if (!name.compare("GGM_CC1ppip_XSec_1DEnu_nu")) { return (new GGM_CC1ppip_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("GGM_CC1ppip_Evt_1DQ2_nu")) { return (new GGM_CC1ppip_Evt_1DQ2_nu(samplekey)); /* MiniBooNE Samples */ /* CCQE */ } else #endif #ifndef __NO_MiniBooNE__ if (!name.compare("MiniBooNE_CCQE_XSec_1DQ2_nu") || !name.compare("MiniBooNE_CCQELike_XSec_1DQ2_nu")) { return (new MiniBooNE_CCQE_XSec_1DQ2_nu(samplekey)); } else if (!name.compare("MiniBooNE_CCQE_XSec_1DEnu_nu") || !name.compare("MiniBooNE_CCQELike_XSec_1DEnu_nu")) { return (new MiniBooNE_CCQE_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("MiniBooNE_CCQE_XSec_1DQ2_antinu") || !name.compare("MiniBooNE_CCQELike_XSec_1DQ2_antinu") || !name.compare("MiniBooNE_CCQE_CTarg_XSec_1DQ2_antinu")) { return (new MiniBooNE_CCQE_XSec_1DQ2_antinu(samplekey)); } else if (!name.compare("MiniBooNE_CCQE_XSec_2DTcos_nu") || !name.compare("MiniBooNE_CCQELike_XSec_2DTcos_nu")) { return (new MiniBooNE_CCQE_XSec_2DTcos_nu(samplekey)); } else if (!name.compare("MiniBooNE_CCQE_XSec_2DTcos_antinu") || !name.compare("MiniBooNE_CCQELike_XSec_2DTcos_antinu")) { return (new MiniBooNE_CCQE_XSec_2DTcos_antinu(samplekey)); /* MiniBooNE CC1pi+ */ // 1D } else if (!name.compare("MiniBooNE_CC1pip_XSec_1DEnu_nu")) { return (new MiniBooNE_CC1pip_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("MiniBooNE_CC1pip_XSec_1DQ2_nu")) { return (new MiniBooNE_CC1pip_XSec_1DQ2_nu(samplekey)); } else if (!name.compare("MiniBooNE_CC1pip_XSec_1DTpi_nu")) { return (new MiniBooNE_CC1pip_XSec_1DTpi_nu(samplekey)); } else if (!name.compare("MiniBooNE_CC1pip_XSec_1DTu_nu")) { return (new MiniBooNE_CC1pip_XSec_1DTu_nu(samplekey)); // 2D } else if (!name.compare("MiniBooNE_CC1pip_XSec_2DQ2Enu_nu")) { return (new MiniBooNE_CC1pip_XSec_2DQ2Enu_nu(samplekey)); } else if (!name.compare("MiniBooNE_CC1pip_XSec_2DTpiCospi_nu")) { return (new MiniBooNE_CC1pip_XSec_2DTpiCospi_nu(samplekey)); } else if (!name.compare("MiniBooNE_CC1pip_XSec_2DTpiEnu_nu")) { return (new MiniBooNE_CC1pip_XSec_2DTpiEnu_nu(samplekey)); } else if (!name.compare("MiniBooNE_CC1pip_XSec_2DTuCosmu_nu")) { return (new MiniBooNE_CC1pip_XSec_2DTuCosmu_nu(samplekey)); } else if (!name.compare("MiniBooNE_CC1pip_XSec_2DTuEnu_nu")) { return (new MiniBooNE_CC1pip_XSec_2DTuEnu_nu(samplekey)); /* MiniBooNE CC1pi0 */ } else if (!name.compare("MiniBooNE_CC1pi0_XSec_1DEnu_nu")) { return (new MiniBooNE_CC1pi0_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("MiniBooNE_CC1pi0_XSec_1DQ2_nu")) { return (new MiniBooNE_CC1pi0_XSec_1DQ2_nu(samplekey)); } else if (!name.compare("MiniBooNE_CC1pi0_XSec_1DTu_nu")) { return (new MiniBooNE_CC1pi0_XSec_1DTu_nu(samplekey)); } else if (!name.compare("MiniBooNE_CC1pi0_XSec_1Dcosmu_nu")) { return (new MiniBooNE_CC1pi0_XSec_1Dcosmu_nu(samplekey)); } else if (!name.compare("MiniBooNE_CC1pi0_XSec_1Dcospi0_nu")) { return (new MiniBooNE_CC1pi0_XSec_1Dcospi0_nu(samplekey)); } else if (!name.compare("MiniBooNE_CC1pi0_XSec_1Dppi0_nu")) { return (new MiniBooNE_CC1pi0_XSec_1Dppi0_nu(samplekey)); } else if (!name.compare("MiniBooNE_NC1pi0_XSec_1Dcospi0_antinu") || !name.compare("MiniBooNE_NC1pi0_XSec_1Dcospi0_rhc")) { return (new MiniBooNE_NC1pi0_XSec_1Dcospi0_antinu(samplekey)); } else if (!name.compare("MiniBooNE_NC1pi0_XSec_1Dcospi0_nu") || !name.compare("MiniBooNE_NC1pi0_XSec_1Dcospi0_fhc")) { return (new MiniBooNE_NC1pi0_XSec_1Dcospi0_nu(samplekey)); } else if (!name.compare("MiniBooNE_NC1pi0_XSec_1Dppi0_antinu") || !name.compare("MiniBooNE_NC1pi0_XSec_1Dppi0_rhc")) { return (new MiniBooNE_NC1pi0_XSec_1Dppi0_antinu(samplekey)); } else if (!name.compare("MiniBooNE_NC1pi0_XSec_1Dppi0_nu") || !name.compare("MiniBooNE_NC1pi0_XSec_1Dppi0_fhc")) { return (new MiniBooNE_NC1pi0_XSec_1Dppi0_nu(samplekey)); /* MiniBooNE NCEL */ } else if (!name.compare("MiniBooNE_NCEL_XSec_Treco_nu")) { return (new MiniBooNE_NCEL_XSec_Treco_nu(samplekey)); } else #endif #ifndef __NO_MicroBooNE__ /* MicroBooNE Samples */ /* MicroBooNE CCinclusive */ if (!name.compare("MicroBooNE_CCInc_XSec_2DPcos_nu")) { return (new MicroBooNE_CCInc_XSec_2DPcos_nu(samplekey)); } else #endif #ifndef __NO_MINERvA__ /* MINERvA Samples */ if (!name.compare("MINERvA_CCQE_XSec_1DQ2_nu") || !name.compare("MINERvA_CCQE_XSec_1DQ2_nu_20deg") || !name.compare("MINERvA_CCQE_XSec_1DQ2_nu_oldflux") || !name.compare("MINERvA_CCQE_XSec_1DQ2_nu_20deg_oldflux")) { return (new MINERvA_CCQE_XSec_1DQ2_nu(samplekey)); } else if (!name.compare("MINERvA_CCQE_XSec_1DQ2_antinu") || !name.compare("MINERvA_CCQE_XSec_1DQ2_antinu_20deg") || !name.compare("MINERvA_CCQE_XSec_1DQ2_antinu_oldflux") || !name.compare("MINERvA_CCQE_XSec_1DQ2_antinu_20deg_oldflux")) { return (new MINERvA_CCQE_XSec_1DQ2_antinu(samplekey)); } else if (!name.compare("MINERvA_CCQE_XSec_1DQ2_joint_oldflux") || !name.compare("MINERvA_CCQE_XSec_1DQ2_joint_20deg_oldflux") || !name.compare("MINERvA_CCQE_XSec_1DQ2_joint") || !name.compare("MINERvA_CCQE_XSec_1DQ2_joint_20deg")) { return (new MINERvA_CCQE_XSec_1DQ2_joint(samplekey)); } else if (!name.compare("MINERvA_CC0pi_XSec_1DEe_nue")) { return (new MINERvA_CC0pi_XSec_1DEe_nue(samplekey)); } else if (!name.compare("MINERvA_CC0pi_XSec_1DQ2_nue")) { return (new MINERvA_CC0pi_XSec_1DQ2_nue(samplekey)); } else if (!name.compare("MINERvA_CC0pi_XSec_1DThetae_nue")) { return (new MINERvA_CC0pi_XSec_1DThetae_nue(samplekey)); } else if (!name.compare("MINERvA_CC0pinp_STV_XSec_1Dpmu_nu") || !name.compare("MINERvA_CC0pinp_STV_XSec_1Dthmu_nu") || !name.compare("MINERvA_CC0pinp_STV_XSec_1Dpprot_nu") || !name.compare("MINERvA_CC0pinp_STV_XSec_1Dthprot_nu") || !name.compare("MINERvA_CC0pinp_STV_XSec_1Dpnreco_nu") || !name.compare("MINERvA_CC0pinp_STV_XSec_1Ddalphat_nu") || !name.compare("MINERvA_CC0pinp_STV_XSec_1Ddpt_nu") || !name.compare("MINERvA_CC0pinp_STV_XSec_1Ddphit_nu")) { return (new MINERvA_CC0pinp_STV_XSec_1D_nu(samplekey)); } else if (!name.compare("MINERvA_CC0pi_XSec_1DQ2_nu_proton")) { return (new MINERvA_CC0pi_XSec_1DQ2_nu_proton(samplekey)); } else if (!name.compare("MINERvA_CC0pi_XSec_1DQ2_TgtC_nu") || !name.compare("MINERvA_CC0pi_XSec_1DQ2_TgtCH_nu") || !name.compare("MINERvA_CC0pi_XSec_1DQ2_TgtFe_nu") || !name.compare("MINERvA_CC0pi_XSec_1DQ2_TgtPb_nu")) { return (new MINERvA_CC0pi_XSec_1DQ2_Tgt_nu(samplekey)); } else if (!name.compare("MINERvA_CC0pi_XSec_1DQ2_TgtRatioC_nu") || !name.compare("MINERvA_CC0pi_XSec_1DQ2_TgtRatioFe_nu") || !name.compare("MINERvA_CC0pi_XSec_1DQ2_TgtRatioPb_nu")) { return (new MINERvA_CC0pi_XSec_1DQ2_TgtRatio_nu(samplekey)); // Dan Ruterbories measurements of late 2018 } else if (!name.compare("MINERvA_CC0pi_XSec_2Dptpz_nu")) { return (new MINERvA_CC0pi_XSec_2D_nu(samplekey)); // } else if (!name.compare("MINERvA_CC0pi_XSec_3DptpzTp_nu")) { // return (new MINERvA_CC0pi_XSec_3DptpzTp_nu(samplekey)); } else if (!name.compare("MINERvA_CC0pi_XSec_1Dpt_nu") || !name.compare("MINERvA_CC0pi_XSec_1Dpz_nu") || !name.compare("MINERvA_CC0pi_XSec_1DQ2QE_nu") || !name.compare("MINERvA_CC0pi_XSec_1DEnuQE_nu")) { return (new MINERvA_CC0pi_XSec_1D_2018_nu(samplekey)); // C. Patrick's early 2018 measurements } else if (!name.compare("MINERvA_CC0pi_XSec_2Dptpz_antinu") || !name.compare("MINERvA_CC0pi_XSec_2DQ2QEEnuQE_antinu") || !name.compare("MINERvA_CC0pi_XSec_2DQ2QEEnuTrue_antinu")) { return (new MINERvA_CC0pi_XSec_2D_antinu(samplekey)); /* CC1pi+ */ // DONE } else if (!name.compare("MINERvA_CC1pip_XSec_1DTpi_nu") || !name.compare("MINERvA_CC1pip_XSec_1DTpi_nu_20deg") || !name.compare("MINERvA_CC1pip_XSec_1DTpi_nu_fluxcorr") || !name.compare("MINERvA_CC1pip_XSec_1DTpi_nu_20deg_fluxcorr")) { return (new MINERvA_CC1pip_XSec_1DTpi_nu(samplekey)); // DONE } else if (!name.compare("MINERvA_CC1pip_XSec_1Dth_nu") || !name.compare("MINERvA_CC1pip_XSec_1Dth_nu_20deg") || !name.compare("MINERvA_CC1pip_XSec_1Dth_nu_fluxcorr") || !name.compare("MINERvA_CC1pip_XSec_1Dth_nu_20deg_fluxcorr")) { return (new MINERvA_CC1pip_XSec_1Dth_nu(samplekey)); } else if (!name.compare("MINERvA_CC1pip_XSec_1DTpi_nu_2017") || !name.compare("MINERvA_CC1pip_XSec_1Dth_nu_2017") || !name.compare("MINERvA_CC1pip_XSec_1Dpmu_nu_2017") || !name.compare("MINERvA_CC1pip_XSec_1Dthmu_nu_2017") || !name.compare("MINERvA_CC1pip_XSec_1DQ2_nu_2017") || !name.compare("MINERvA_CC1pip_XSec_1DEnu_nu_2017")) { return (new MINERvA_CC1pip_XSec_1D_2017Update(samplekey)); /* CCNpi+ */ } else if (!name.compare("MINERvA_CCNpip_XSec_1Dth_nu") || !name.compare("MINERvA_CCNpip_XSec_1Dth_nu_2015") || !name.compare("MINERvA_CCNpip_XSec_1Dth_nu_2016") || !name.compare("MINERvA_CCNpip_XSec_1Dth_nu_2015_20deg") || !name.compare("MINERvA_CCNpip_XSec_1Dth_nu_2015_fluxcorr") || !name.compare("MINERvA_CCNpip_XSec_1Dth_nu_2015_20deg_fluxcorr")) { return (new MINERvA_CCNpip_XSec_1Dth_nu(samplekey)); } else if (!name.compare("MINERvA_CCNpip_XSec_1DTpi_nu") || !name.compare("MINERvA_CCNpip_XSec_1DTpi_nu_2015") || !name.compare("MINERvA_CCNpip_XSec_1DTpi_nu_2016") || !name.compare("MINERvA_CCNpip_XSec_1DTpi_nu_2015_20deg") || !name.compare("MINERvA_CCNpip_XSec_1DTpi_nu_2015_fluxcorr") || !name.compare( "MINERvA_CCNpip_XSec_1DTpi_nu_2015_20deg_fluxcorr")) { return (new MINERvA_CCNpip_XSec_1DTpi_nu(samplekey)); } else if (!name.compare("MINERvA_CCNpip_XSec_1Dthmu_nu")) { return (new MINERvA_CCNpip_XSec_1Dthmu_nu(samplekey)); } else if (!name.compare("MINERvA_CCNpip_XSec_1Dpmu_nu")) { return (new MINERvA_CCNpip_XSec_1Dpmu_nu(samplekey)); } else if (!name.compare("MINERvA_CCNpip_XSec_1DQ2_nu")) { return (new MINERvA_CCNpip_XSec_1DQ2_nu(samplekey)); } else if (!name.compare("MINERvA_CCNpip_XSec_1DEnu_nu")) { return (new MINERvA_CCNpip_XSec_1DEnu_nu(samplekey)); /* MINERvA CC1pi0 anti-nu */ // Done } else if (!name.compare("MINERvA_CC1pi0_XSec_1Dth_antinu") || !name.compare("MINERvA_CC1pi0_XSec_1Dth_antinu_2015") || !name.compare("MINERvA_CC1pi0_XSec_1Dth_antinu_2016") || !name.compare("MINERvA_CC1pi0_XSec_1Dth_antinu_fluxcorr") || !name.compare("MINERvA_CC1pi0_XSec_1Dth_antinu_2015_fluxcorr") || !name.compare("MINERvA_CC1pi0_XSec_1Dth_antinu_2016_fluxcorr")) { return (new MINERvA_CC1pi0_XSec_1Dth_antinu(samplekey)); } else if (!name.compare("MINERvA_CC1pi0_XSec_1Dppi0_antinu") || !name.compare("MINERvA_CC1pi0_XSec_1Dppi0_antinu_fluxcorr")) { return (new MINERvA_CC1pi0_XSec_1Dppi0_antinu(samplekey)); } else if (!name.compare("MINERvA_CC1pi0_XSec_1DTpi0_antinu")) { return (new MINERvA_CC1pi0_XSec_1DTpi0_antinu(samplekey)); // Done } else if (!name.compare("MINERvA_CC1pi0_XSec_1DQ2_antinu")) { return (new MINERvA_CC1pi0_XSec_1DQ2_antinu(samplekey)); // Done } else if (!name.compare("MINERvA_CC1pi0_XSec_1Dthmu_antinu")) { return (new MINERvA_CC1pi0_XSec_1Dthmu_antinu(samplekey)); // Done } else if (!name.compare("MINERvA_CC1pi0_XSec_1Dpmu_antinu")) { return (new MINERvA_CC1pi0_XSec_1Dpmu_antinu(samplekey)); // Done } else if (!name.compare("MINERvA_CC1pi0_XSec_1DEnu_antinu")) { return (new MINERvA_CC1pi0_XSec_1DEnu_antinu(samplekey)); // MINERvA CC1pi0 nu } else if (!name.compare("MINERvA_CC1pi0_XSec_1DTpi_nu") || !name.compare("MINERvA_CC1pi0_XSec_1Dth_nu") || !name.compare("MINERvA_CC1pi0_XSec_1Dpmu_nu") || !name.compare("MINERvA_CC1pi0_XSec_1Dthmu_nu") || !name.compare("MINERvA_CC1pi0_XSec_1DQ2_nu") || !name.compare("MINERvA_CC1pi0_XSec_1DEnu_nu") || !name.compare("MINERvA_CC1pi0_XSec_1DWexp_nu") || !name.compare("MINERvA_CC1pi0_XSec_1DPPi0Mass_nu") || !name.compare("MINERvA_CC1pi0_XSec_1DPPi0MassDelta_nu") || !name.compare("MINERvA_CC1pi0_XSec_1DCosAdler_nu") || !name.compare("MINERvA_CC1pi0_XSec_1DPhiAdler_nu")) { return (new MINERvA_CC1pi0_XSec_1D_nu(samplekey)); /* CCINC */ } else if (!name.compare("MINERvA_CCinc_XSec_2DEavq3_nu")) { return (new MINERvA_CCinc_XSec_2DEavq3_nu(samplekey)); } else if (!name.compare("MINERvA_CCinc_XSec_1Dx_ratio_C12_CH") || !name.compare("MINERvA_CCinc_XSec_1Dx_ratio_Fe56_CH") || !name.compare("MINERvA_CCinc_XSec_1Dx_ratio_Pb208_CH")) { return (new MINERvA_CCinc_XSec_1Dx_ratio(samplekey)); } else if (!name.compare("MINERvA_CCinc_XSec_1DEnu_ratio_C12_CH") || !name.compare("MINERvA_CCinc_XSec_1DEnu_ratio_Fe56_CH") || !name.compare("MINERvA_CCinc_XSec_1DEnu_ratio_Pb208_CH")) { return (new MINERvA_CCinc_XSec_1DEnu_ratio(samplekey)); /* CCDIS */ } else if (!name.compare("MINERvA_CCDIS_XSec_1Dx_ratio_C12_CH") || !name.compare("MINERvA_CCDIS_XSec_1Dx_ratio_Fe56_CH") || !name.compare("MINERvA_CCDIS_XSec_1Dx_ratio_Pb208_CH")) { return (new MINERvA_CCDIS_XSec_1Dx_ratio(samplekey)); } else if (!name.compare("MINERvA_CCDIS_XSec_1DEnu_ratio_C12_CH") || !name.compare("MINERvA_CCDIS_XSec_1DEnu_ratio_Fe56_CH") || !name.compare("MINERvA_CCDIS_XSec_1DEnu_ratio_Pb208_CH")) { return (new MINERvA_CCDIS_XSec_1DEnu_ratio(samplekey)); /* CC-COH */ } else if (!name.compare("MINERvA_CCCOHPI_XSec_1DEnu_nu")) { return (new MINERvA_CCCOHPI_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("MINERvA_CCCOHPI_XSec_1DEpi_nu")) { return (new MINERvA_CCCOHPI_XSec_1DEpi_nu(samplekey)); } else if (!name.compare("MINERvA_CCCOHPI_XSec_1Dth_nu")) { return (new MINERvA_CCCOHPI_XSec_1Dth_nu(samplekey)); } else if (!name.compare("MINERvA_CCCOHPI_XSec_1DQ2_nu")) { return (new MINERvA_CCCOHPI_XSec_1DQ2_nu(samplekey)); } else if (!name.compare("MINERvA_CCCOHPI_XSec_1DEnu_antinu")) { return (new MINERvA_CCCOHPI_XSec_1DEnu_antinu(samplekey)); } else if (!name.compare("MINERvA_CCCOHPI_XSec_1DEpi_antinu")) { return (new MINERvA_CCCOHPI_XSec_1DEpi_antinu(samplekey)); } else if (!name.compare("MINERvA_CCCOHPI_XSec_1Dth_antinu")) { return (new MINERvA_CCCOHPI_XSec_1Dth_antinu(samplekey)); } else if (!name.compare("MINERvA_CCCOHPI_XSec_1DQ2_antinu")) { return (new MINERvA_CCCOHPI_XSec_1DQ2_antinu(samplekey)); } else if (!name.compare("MINERvA_CCCOHPI_XSec_1DEnu_joint")) { return (new MINERvA_CCCOHPI_XSec_joint(samplekey)); } else if (!name.compare("MINERvA_CCCOHPI_XSec_1DEpi_joint")) { return (new MINERvA_CCCOHPI_XSec_joint(samplekey)); } else if (!name.compare("MINERvA_CCCOHPI_XSec_1Dth_joint")) { return (new MINERvA_CCCOHPI_XSec_joint(samplekey)); } else if (!name.compare("MINERvA_CCCOHPI_XSec_1DQ2_joint")) { return (new MINERvA_CCCOHPI_XSec_joint(samplekey)); /* T2K Samples */ } else #endif #ifndef __NO_T2K__ if (!name.compare("T2K_CC0pi_XSec_2DPcos_nu_I")) { return (new T2K_CC0pi_XSec_2DPcos_nu_I(samplekey)); } else if (!name.compare("T2K_CC0pi_XSec_2DPcos_nu_II")) { return (new T2K_CC0pi_XSec_2DPcos_nu_II(samplekey)); } else if (!name.compare("T2K_CCinc_XSec_2DPcos_nu_nonuniform")) { return (new T2K_CCinc_XSec_2DPcos_nu_nonuniform(samplekey)); } else if (!name.compare("T2K_CC0pi_XSec_H2O_2DPcos_anu")) { return (new T2K_CC0pi_XSec_H2O_2DPcos_anu(samplekey)); } else if (!name.compare("T2K_NuMu_CC0pi_O_XSec_2DPcos") || !name.compare("T2K_NuMu_CC0pi_C_XSec_2DPcos")) { return (new T2K_NuMu_CC0pi_OC_XSec_2DPcos(samplekey)); } else if (!name.compare("T2K_NuMu_CC0pi_OC_XSec_2DPcos_joint")) { return (new T2K_NuMu_CC0pi_OC_XSec_2DPcos_joint(samplekey)); - + } else if (!name.compare("T2K_NuMu_CC0pi_CH_XSec_2DPcos") || !name.compare("T2K_AntiNuMu_CC0pi_CH_XSec_2DPcos")) { return (new T2K_NuMuAntiNuMu_CC0pi_CH_XSec_2DPcos(samplekey)); } else if (!name.compare("T2K_NuMuAntiNuMu_CC0pi_CH_XSec_2DPcos_joint")) { return (new T2K_NuMuAntiNuMu_CC0pi_CH_XSec_2DPcos_joint(samplekey)); - + } else if (!name.compare("T2K_nueCCinc_XSec_1Dpe_FHC") || !name.compare("T2K_nueCCinc_XSec_1Dpe_RHC") || !name.compare("T2K_nuebarCCinc_XSec_1Dpe_RHC")) { return (new T2K_nueCCinc_XSec_1Dpe(samplekey)); } else if (!name.compare("T2K_nueCCinc_XSec_1Dthe_FHC") || !name.compare("T2K_nueCCinc_XSec_1Dthe_RHC") || !name.compare("T2K_nuebarCCinc_XSec_1Dthe_RHC")) { return (new T2K_nueCCinc_XSec_1Dthe(samplekey)); - + } else if (!name.compare("T2K_nueCCinc_XSec_1Dpe_joint")) { return (new T2K_nueCCinc_XSec_1Dpe_joint(samplekey)); } else if (!name.compare("T2K_nueCCinc_XSec_1Dthe_joint")) { return (new T2K_nueCCinc_XSec_1Dthe_joint(samplekey)); } else if (!name.compare("T2K_nueCCinc_XSec_joint")) { return (new T2K_nueCCinc_XSec_joint(samplekey)); /* T2K CC1pi+ CH samples */ // Comment these out for now because we don't have the proper data } else if (!name.compare("T2K_CC1pip_CH_XSec_2Dpmucosmu_nu")) { return (new T2K_CC1pip_CH_XSec_2Dpmucosmu_nu(samplekey)); } else if (!name.compare("T2K_CC1pip_CH_XSec_1Dppi_nu")) { return (new T2K_CC1pip_CH_XSec_1Dppi_nu(samplekey)); } else if (!name.compare("T2K_CC1pip_CH_XSec_1Dthpi_nu")) { return (new T2K_CC1pip_CH_XSec_1Dthpi_nu(samplekey)); } else if (!name.compare("T2K_CC1pip_CH_XSec_1Dthmupi_nu")) { return (new T2K_CC1pip_CH_XSec_1Dthmupi_nu(samplekey)); } else if (!name.compare("T2K_CC1pip_CH_XSec_1DQ2_nu")) { return (new T2K_CC1pip_CH_XSec_1DQ2_nu(samplekey)); } else if (!name.compare("T2K_CC1pip_CH_XSec_1DAdlerPhi_nu")) { return (new T2K_CC1pip_CH_XSec_1DAdlerPhi_nu(samplekey)); } else if (!name.compare("T2K_CC1pip_CH_XSec_1DCosThAdler_nu")) { return (new T2K_CC1pip_CH_XSec_1DCosThAdler_nu(samplekey)); /* T2K CC1pi+ H2O samples */ } else if (!name.compare("T2K_CC1pip_H2O_XSec_1DEnuDelta_nu")) { return (new T2K_CC1pip_H2O_XSec_1DEnuDelta_nu(samplekey)); } else if (!name.compare("T2K_CC1pip_H2O_XSec_1DEnuMB_nu")) { return (new T2K_CC1pip_H2O_XSec_1DEnuMB_nu(samplekey)); } else if (!name.compare("T2K_CC1pip_H2O_XSec_1Dcosmu_nu")) { return (new T2K_CC1pip_H2O_XSec_1Dcosmu_nu(samplekey)); } else if (!name.compare("T2K_CC1pip_H2O_XSec_1Dcosmupi_nu")) { return (new T2K_CC1pip_H2O_XSec_1Dcosmupi_nu(samplekey)); } else if (!name.compare("T2K_CC1pip_H2O_XSec_1Dcospi_nu")) { return (new T2K_CC1pip_H2O_XSec_1Dcospi_nu(samplekey)); } else if (!name.compare("T2K_CC1pip_H2O_XSec_1Dpmu_nu")) { return (new T2K_CC1pip_H2O_XSec_1Dpmu_nu(samplekey)); } else if (!name.compare("T2K_CC1pip_H2O_XSec_1Dppi_nu")) { return (new T2K_CC1pip_H2O_XSec_1Dppi_nu(samplekey)); /* T2K CC0pi + np CH samples */ } else if (!name.compare("T2K_CC0pinp_STV_XSec_1Ddpt_nu")) { return (new T2K_CC0pinp_STV_XSec_1Ddpt_nu(samplekey)); } else if (!name.compare("T2K_CC0pinp_STV_XSec_1Ddphit_nu")) { return (new T2K_CC0pinp_STV_XSec_1Ddphit_nu(samplekey)); } else if (!name.compare("T2K_CC0pinp_STV_XSec_1Ddat_nu")) { return (new T2K_CC0pinp_STV_XSec_1Ddat_nu(samplekey)); - } else if (!name.compare("T2K_CC0pi1p_XSec_3DPcoscos_nu_nonuniform")) { - return (new T2K_CC0pi1p_XSec_3DPcoscos_nu_nonuniform(samplekey)); + } else if (!name.compare("T2K_CC0piWithProtons_2018_multidif_0p_1p_Np") || !name.compare("T2K_CC0piWithProtons_2018_multidif_0p_1p") || !name.compare("T2K_CC0piWithProtons_2018_multidif_0p") || !name.compare("T2K_CC0piWithProtons_2018_multidif_1p")) { + return (new T2K_CC0T2K_CC0piWithProtons_2018_multidif_0p_1p_Np(samplekey)); } else if (!name.compare("T2K_CC0pinp_ifk_XSec_3Dinfp_nu")) { return (new T2K_CC0pinp_ifk_XSec_3Dinfp_nu(samplekey)); } else if (!name.compare("T2K_CC0pinp_ifk_XSec_3Dinfa_nu")) { return (new T2K_CC0pinp_ifk_XSec_3Dinfa_nu(samplekey)); } else if (!name.compare("T2K_CC0pinp_ifk_XSec_3Dinfip_nu")) { return (new T2K_CC0pinp_ifk_XSec_3Dinfip_nu(samplekey)); // SciBooNE COH studies } else #endif #ifndef __NO_SciBooNE__ if (!name.compare("SciBooNE_CCCOH_STOP_NTrks_nu")) { return (new SciBooNE_CCCOH_STOP_NTrks_nu(samplekey)); } else if (!name.compare("SciBooNE_CCCOH_1TRK_1DQ2_nu")) { return (new SciBooNE_CCCOH_1TRK_1DQ2_nu(samplekey)); } else if (!name.compare("SciBooNE_CCCOH_1TRK_1Dpmu_nu")) { return (new SciBooNE_CCCOH_1TRK_1Dpmu_nu(samplekey)); } else if (!name.compare("SciBooNE_CCCOH_1TRK_1Dthetamu_nu")) { return (new SciBooNE_CCCOH_1TRK_1Dthetamu_nu(samplekey)); } else if (!name.compare("SciBooNE_CCCOH_MuPr_1DQ2_nu")) { return (new SciBooNE_CCCOH_MuPr_1DQ2_nu(samplekey)); } else if (!name.compare("SciBooNE_CCCOH_MuPr_1Dpmu_nu")) { return (new SciBooNE_CCCOH_MuPr_1Dpmu_nu(samplekey)); } else if (!name.compare("SciBooNE_CCCOH_MuPr_1Dthetamu_nu")) { return (new SciBooNE_CCCOH_MuPr_1Dthetamu_nu(samplekey)); } else if (!name.compare("SciBooNE_CCCOH_MuPiVA_1DQ2_nu")) { return (new SciBooNE_CCCOH_MuPiVA_1DQ2_nu(samplekey)); } else if (!name.compare("SciBooNE_CCCOH_MuPiVA_1Dpmu_nu")) { return (new SciBooNE_CCCOH_MuPiVA_1Dpmu_nu(samplekey)); } else if (!name.compare("SciBooNE_CCCOH_MuPiVA_1Dthetamu_nu")) { return (new SciBooNE_CCCOH_MuPiVA_1Dthetamu_nu(samplekey)); } else if (!name.compare("SciBooNE_CCCOH_MuPiNoVA_1DQ2_nu")) { return (new SciBooNE_CCCOH_MuPiNoVA_1DQ2_nu(samplekey)); } else if (!name.compare("SciBooNE_CCCOH_MuPiNoVA_1Dthetapr_nu")) { return (new SciBooNE_CCCOH_MuPiNoVA_1Dthetapr_nu(samplekey)); } else if (!name.compare("SciBooNE_CCCOH_MuPiNoVA_1Dthetapi_nu")) { return (new SciBooNE_CCCOH_MuPiNoVA_1Dthetapi_nu(samplekey)); } else if (!name.compare("SciBooNE_CCCOH_MuPiNoVA_1Dthetamu_nu")) { return (new SciBooNE_CCCOH_MuPiNoVA_1Dthetamu_nu(samplekey)); } else if (!name.compare("SciBooNE_CCCOH_MuPiNoVA_1Dpmu_nu")) { return (new SciBooNE_CCCOH_MuPiNoVA_1Dpmu_nu(samplekey)); } else if (!name.compare("SciBooNE_CCCOH_STOPFINAL_1DQ2_nu")) { return (new SciBooNE_CCCOH_STOPFINAL_1DQ2_nu(samplekey)); } else if (!name.compare("SciBooNE_CCInc_XSec_1DEnu_nu") || !name.compare("SciBooNE_CCInc_XSec_1DEnu_nu_NEUT") || !name.compare("SciBooNE_CCInc_XSec_1DEnu_nu_NUANCE")) { return (new SciBooNE_CCInc_XSec_1DEnu_nu(samplekey)); /* K2K Samples */ /* NC1pi0 */ } else #endif #ifndef __NO_K2K__ if (!name.compare("K2K_NC1pi0_Evt_1Dppi0_nu")) { return (new K2K_NC1pi0_Evt_1Dppi0_nu(samplekey)); /* Fake Studies */ } else #endif if (name.find("ExpMultDist_CCQE_XSec_1D") != std::string::npos && name.find("_FakeStudy") != std::string::npos) { return ( new ExpMultDist_CCQE_XSec_1DVar_FakeStudy(name, file, rw, type, fkdt)); } else if (name.find("ExpMultDist_CCQE_XSec_2D") != std::string::npos && name.find("_FakeStudy") != std::string::npos) { return ( new ExpMultDist_CCQE_XSec_2DVar_FakeStudy(name, file, rw, type, fkdt)); } else if (name.find("GenericFlux_") != std::string::npos) { return (new GenericFlux_Tester(name, file, rw, type, fkdt)); } else if (name.find("GenericVectors_") != std::string::npos) { return (new GenericFlux_Vectors(name, file, rw, type, fkdt)); } else if (!name.compare("T2K2017_FakeData")) { return (new T2K2017_FakeData(samplekey)); } else if (!name.compare("MCStudy_CCQE")) { return (new MCStudy_CCQEHistograms(name, file, rw, type, fkdt)); } else if (!name.compare("ElectronFlux_FlatTree")) { return (new ElectronFlux_FlatTree(name, file, rw, type, fkdt)); } else if (name.find("ElectronData_") != std::string::npos) { return new ElectronScattering_DurhamData(samplekey); } else if (name.find("MuonValidation_") != std::string::npos) { return (new MCStudy_MuonValidation(name, file, rw, type, fkdt)); } else if (!name.compare("NIWGOfficialPlots")) { return (new OfficialNIWGPlots(samplekey)); } else if ((name.find("SigmaEnuHists") != std::string::npos) || (name.find("SigmaEnuPerEHists") != std::string::npos)) { return (new SigmaEnuHists(samplekey)); } else if (!name.compare("Simple_Osc")) { return (new Simple_Osc(samplekey)); } else if (!name.compare("Smear_SVDUnfold_Propagation_Osc")) { return (new Smear_SVDUnfold_Propagation_Osc(samplekey)); } else { NUIS_ABORT("Error: No such sample: " << name << std::endl); } // Return NULL if no sample loaded. return NULL; } // namespace SampleUtils } // namespace SampleUtils diff --git a/src/T2K/T2K_CC0pi1p_XSec_3DPcoscos_nu_nonuniform.cxx b/src/T2K/T2K_CC0pi1p_XSec_3DPcoscos_nu_nonuniform.cxx deleted file mode 100644 index 82ba031..0000000 --- a/src/T2K/T2K_CC0pi1p_XSec_3DPcoscos_nu_nonuniform.cxx +++ /dev/null @@ -1,252 +0,0 @@ -// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret - -/******************************************************************************* - * This file is part of NUISANCE. - * - * NUISANCE is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. - * - * NUISANCE is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with NUISANCE. If not, see . - *******************************************************************************/ - -#include "T2K_SignalDef.h" - -#include "T2K_CC0pi1p_XSec_3DPcoscos_nu_nonuniform.h" - -//******************************************************************** -T2K_CC0pi1p_XSec_3DPcoscos_nu_nonuniform:: - T2K_CC0pi1p_XSec_3DPcoscos_nu_nonuniform(nuiskey samplekey) { - //******************************************************************** - - // Sample overview --------------------------------------------------- - std::string descrip = "T2K_CC0pi1p_XSec_3DPcoscos_nu_nonuniform sample. \n" - "Target: CH \n" - "Flux: T2K 2.5 degree off-axis (ND280) \n" - "Signal: CC0piNp (N>=1) with p_p>500MeV \n" - "https://doi.org/10.1103/PhysRevD.98.032003 \n"; - - // Setup common settings - fSettings = LoadSampleSettings(samplekey); - fSettings.SetDescription(descrip); - fSettings.SetXTitle("P_{p} (GeV)"); - fSettings.SetYTitle("cos#theta_{p}"); - fSettings.SetZTitle("cos#theta_{#mu}"); - // fSettings.SetZTitle("d^{2}#sigma/dP_{#mu}dcos#theta_{#mu} (cm^{2}/GeV)"); - fSettings.SetAllowedTypes("FULL,DIAG/FREE,SHAPE,FIX/SYSTCOV/STATCOV", - "FIX/FULL"); - fSettings.SetEnuRange(0.0, 10.0); - fSettings.DefineAllowedTargets("C,H"); - - fAnalysis = 1; - - // CCQELike plot information - fSettings.SetTitle("T2K_CC0pi1p_XSec_3DPcoscos_nu_nonuniform"); - fSettings.DefineAllowedSpecies("numu"); - - FinaliseSampleSettings(); - - // Scaling Setup --------------------------------------------------- - // ScaleFactor automatically setup for DiffXSec/cm2/Nucleon - // fScaleFactor = ((GetEventHistogram()->Integral("width")/(fNEvents+0.)) * - // 1E-38 / (TotalIntegratedFlux())); - fScaleFactor = ((GetEventHistogram()->Integral("width") / (fNEvents + 0.)) * - 10 / (TotalIntegratedFlux())); - - // Setup Histograms - SetHistograms(); - - // Final setup --------------------------------------------------- - FinaliseMeasurement(); -}; - -bool T2K_CC0pi1p_XSec_3DPcoscos_nu_nonuniform::isSignal(FitEvent *event) { - return SignalDef::isT2K_CC0pi1p(event, EnuMin, EnuMax); -}; - -void T2K_CC0pi1p_XSec_3DPcoscos_nu_nonuniform::FillEventVariables( - FitEvent *event) { - - if (event->NumFSParticle(13) == 0 || event->NumFSParticle(2212) == 0) - return; - - TLorentzVector Pnu = event->GetNeutrinoIn()->fP; - TLorentzVector Pmu = event->GetHMFSParticle(13)->fP; - TLorentzVector Pp = event->GetHMFSParticle(2212)->fP; - - // double pmu = Pmu.Vect().Mag()/1000.; - double pp = Pp.Vect().Mag() / 1000.; - double CosThetaMu = cos(Pnu.Vect().Angle(Pmu.Vect())); - double CosThetaP = cos(Pnu.Vect().Angle(Pp.Vect())); - - fXVar = pp; - fYVar = CosThetaP; - fZVar = CosThetaMu; - - return; -}; - -void T2K_CC0pi1p_XSec_3DPcoscos_nu_nonuniform::FillHistograms() { - - Measurement1D::FillHistograms(); - if (Signal) { - fMCHist_Fine2D->Fill(fXVar, fYVar, Weight); - FillMCSlice(fXVar, fYVar, fZVar, Weight); - } -} - -// Modification is needed after the full reconfigure to move bins around -// Otherwise this would need to be replaced by a TH2Poly which is too awkward. -void T2K_CC0pi1p_XSec_3DPcoscos_nu_nonuniform::ConvertEventRates() { - - for (int i = 0; i < 4; i++) { - fMCHist_Slices[i]->GetSumw2(); - } - - // Do standard conversion. - Measurement1D::ConvertEventRates(); - - // First scale MC slices also by their width in Y and Z - // MCHist_Slices[0]->Scale(1.0 / 1.00); - // MCHist_Slices[1]->Scale(1.0 / 0.60); - // MCHist_Slices[2]->Scale(1.0 / 0.10); - // MCHist_Slices[3]->Scale(1.0 / 0.10); - - // Now Convert into 1D list - fMCHist->Reset(); - int bincount = 0; - for (int i = 0; i < 4; i++) { - for (int j = 0; j < fDataHist_Slices[i]->GetNumberOfBins(); j++) { - fMCHist->SetBinContent(bincount + 1, - fMCHist_Slices[i]->GetBinContent(j + 1)); - // fMCHist->SetBinError(bincount+1, fMCHist_Slices[i]->GetBinError(j+1)); - bincount++; - } - } - - return; -} - -void T2K_CC0pi1p_XSec_3DPcoscos_nu_nonuniform::FillMCSlice(double x, double y, - double z, double w) { - - if (z >= -1.0 and z < -0.3) - fMCHist_Slices[0]->Fill(y, x, w); - else if (z >= -0.3 and z < 0.3) - fMCHist_Slices[1]->Fill(y, x, w); - else if (z >= 0.3 and z < 0.8) - fMCHist_Slices[2]->Fill(y, x, w); - else if (z >= 0.8 and z < 1.0) - fMCHist_Slices[3]->Fill(y, x, w); -} - -void T2K_CC0pi1p_XSec_3DPcoscos_nu_nonuniform::SetHistograms() { - - // Read in 1D Data Histograms - fInputFile = new TFile( - (FitPar::GetDataBase() + "/T2K/CC0pi/STV/multidif_3D_pcoscos.root") - .c_str(), - "READ"); - // fInputFile->ls(); - - // Read in 1D Data - fDataHist = (TH1D *)fInputFile->Get("LinResult"); - - fMCHist_Fine2D = new TH2D("T2K_CC0pi1p_XSec_3DPcoscos_nu_nonuniform_Fine2D", - "T2K_CC0pi1p_XSec_3DPcoscos_nu_nonuniform_Fine2D", - 400, 0.0, 30.0, 100, -1.0, 1.0); - SetAutoProcessTH1(fMCHist_Fine2D); - - TH2D *tempcov = (TH2D *)fInputFile->Get("CovMatrix"); - - fFullCovar = new TMatrixDSym(fDataHist->GetNbinsX()); - for (int i = 0; i < fDataHist->GetNbinsX(); i++) { - for (int j = 0; j < fDataHist->GetNbinsX(); j++) { - //(*fFullCovar)(i,j) = tempcov->GetBinContent(i+1, j+1); - (*fFullCovar)(i, j) = tempcov->GetBinContent(i + 1, j + 1) * - fDataHist->GetBinContent(i + 1) * - fDataHist->GetBinContent(j + 1); - if (i == j) - fDataHist->SetBinError(i + 1, sqrt((*fFullCovar)(i, j))); - // if(i==j) std::cout << "For bin " << i+1 << ", relative covariance was " - // << tempcov->GetBinContent(i+1,j+1); if(i==j) std::cout << ". Absolute - // covariance is now " << (*fFullCovar)(i,j) << ", linear xsec is: " << - // fDataHist->GetBinContent(i+1) << std::endl; - } - } - covar = StatUtils::GetInvert(fFullCovar); - fDecomp = StatUtils::GetDecomp(fFullCovar); - - TH1D *linearResult = new TH1D(*fDataHist); - linearResult->SetNameTitle("T2K_CC0pi1p_XSec_3DPcoscos_nu_nonuniform_data", - "T2K_CC0pi1p_XSec_3DPcoscos_nu_nonuniform_data"); - SetAutoProcessTH1(linearResult, kCMD_Write); - - // Read in 3D Data - // for (int c = 0; c < 4; c++){ - // fDataPoly[c] = (TH2Poly*) fInputFile->Get(Form("dataslice_%i",c)); - // fDataPoly[c]->SetNameTitle(Form("T2K_CC0pi1p_XSec_3DPcoscos_nu_nonuniform_datapoly_%i",c),Form("T2K_CC0pi1p_XSec_3DPcoscos_nu_nonuniform_datapoly_%i",c)); - // SetAutoProcessTH1(fDataPoly[c], kCMD_Write); - // fDataHist->Reset(); - // } - - // Read in 2D Data Slices and Make MC Slices - fDataHist->Reset(); - int bincount = 0; - for (int i = 0; i < 4; i++) { // both y and z slices - - // Get Data Histogram - // fInputFile->ls(); - fDataHist_Slices.push_back( - (TH2Poly *)fInputFile->Get(Form("dataslice_%i", i))->Clone()); - fDataHist_Slices[i]->SetNameTitle( - Form("T2K_CC0pi1p_XSec_3DPcoscos_nu_nonuniform_data_Slice%i", i), - (Form("T2K_CC0pi1p_XSec_3DPcoscos_nu_nonuniform_data_Slice%i", i))); - - for (int j = 0; j < fDataHist_Slices[i]->GetNumberOfBins(); j++) { - fDataHist_Slices[i]->SetBinError( - j, sqrt(tempcov->GetBinContent(bincount, bincount))); - fDataHist->SetBinContent(bincount + 1, - fDataHist_Slices[i]->GetBinContent(j + 1)); - fDataHist->SetBinError(bincount + 1, - fDataHist_Slices[i]->GetBinError(j + 1)); - bincount++; - } - - // Loop over nbins and set errors from covar - // for (int j = 0; j < fDataHist_Slices[i]->GetNbinsX(); j++){ - // for (int k = 0; k < fDataHist_Slices[i]->GetNbinsY(); k++){ - // fDataHist_Slices[i]->SetBinError(j+1, k+1, - // sqrt((*fFullCovar)(bincount,bincount)) * 1E-38); - // - // std::cout << "Setting data hist " << - // fDataHist_Slices[i]->GetBinContent(j+1,k+1) << " " << - // fDataHist_Slices[i]->GetBinError(j+1,k+1) << std::endl; - // fDataHist->SetBinContent(bincount+1, - // fDataHist_Slices[i]->GetBinContent(j+1,k+1) ); - // fDataHist->SetBinError(bincount+1, - // fDataHist_Slices[i]->GetBinError(j+1,k+1) ); - // - // bincount++; - // } - // } - - // Make MC Clones - fMCHist_Slices.push_back((TH2Poly *)fDataHist_Slices[i]->Clone()); - fMCHist_Slices[i]->SetNameTitle( - Form("T2K_CC0pi1p_XSec_3DPcoscos_nu_nonuniform_MC_Slice%i", i), - (Form("T2K_CC0pi1p_XSec_3DPcoscos_nu_nonuniform_MC_Slice%i", i))); - - SetAutoProcessTH1(fDataHist_Slices[i], kCMD_Write); - SetAutoProcessTH1(fMCHist_Slices[i]); - fMCHist_Slices[i]->ClearBinContents(); - } - return; -}; diff --git a/src/T2K/T2K_CC0piWithProtons_XSec_2018_multidif_0p_1p_Np.cxx b/src/T2K/T2K_CC0piWithProtons_XSec_2018_multidif_0p_1p_Np.cxx new file mode 100644 index 0000000..f8e015e --- /dev/null +++ b/src/T2K/T2K_CC0piWithProtons_XSec_2018_multidif_0p_1p_Np.cxx @@ -0,0 +1,641 @@ +// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret + +/******************************************************************************* + * This file is part of NUISANCE. + * + * NUISANCE is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * NUISANCE is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with NUISANCE. If not, see . + *******************************************************************************/ + +#include "T2K_SignalDef.h" + +#include "T2K_CC0piWithProtons_2018_multidif_0p_1p_Np.h" + +//******************************************************************** +T2K_CC0piWithProtons_2018_multidif_0p_1p_Np:: + T2K_CC0piWithProtons_2018_multidif_0p_1p_Np(nuiskey samplekey) { + //******************************************************************** + + // Sample overview --------------------------------------------------- + std::string descrip = "T2K_CC0piWithProtons_2018_multidif_0p_1p_Np sample. \n" + "Target: CH \n" + "Flux: T2K 2.5 degree off-axis (ND280) \n" + "Signal: CC0piNp (N>=0) with p_p>500MeV \n" + "https://doi.org/10.1103/PhysRevD.98.032003 \n" + "Data release: https://t2k-experiment.org/results/2018-transverse-cc0pi/ \n"; + + // Setup common settings + fSettings = LoadSampleSettings(samplekey); + fSettings.SetDescription(descrip); + fSettings.SetXTitle("P_{p} (GeV)"); + fSettings.SetYTitle("cos#theta_{p}"); + fSettings.SetZTitle("cos#theta_{#mu}"); + // fSettings.SetZTitle("d^{2}#sigma/dP_{#mu}dcos#theta_{#mu} (cm^{2}/GeV)"); + fSettings.SetAllowedTypes("FULL,DIAG/FREE,SHAPE,FIX/SYSTCOV/STATCOV", + "FIX/FULL"); + fSettings.SetEnuRange(0.0, 10.0); + fSettings.DefineAllowedTargets("C,H"); + + fAnalysis = 1; + + // CCQELike plot information + fSettings.SetTitle("T2K_CC0piWithProtons_2018_multidif_0p_1p_Np"); + fSettings.DefineAllowedSpecies("numu"); + + FinaliseSampleSettings(); + + //TODO: set useCC0pi0p, useCC0pi1p, and useCC0piNp *before this point* + if (fName == "T2K_CC0piWithProtons_2018_multidif_0p_1p_Np"){ + useCC0pi0p = true; + useCC0pi1p = true; + useCC0piNp = true; + } + else if (fName == "T2K_CC0piWithProtons_2018_multidif_0p_1p"){ + useCC0pi0p = true; + useCC0pi1p = true; + } + else if (fName == "T2K_CC0piWithProtons_2018_multidif_0p"){ + useCC0pi0p = true; + } + else if (fName == "T2K_CC0piWithProtons_2018_multidif_1p"){ + useCC0pi1p = true; + } + + + // Scaling Setup --------------------------------------------------- + // ScaleFactor automatically setup for DiffXSec/cm2/Nucleon + // fScaleFactor = ((GetEventHistogram()->Integral("width")/(fNEvents+0.)) * + // 1E-38 / (TotalIntegratedFlux())); + fScaleFactor = ((GetEventHistogram()->Integral("width") / (fNEvents + 0.)) * + 10 / (TotalIntegratedFlux())); + + // Setup Histograms + SetHistograms(); + + // Final setup --------------------------------------------------- + FinaliseMeasurement(); +}; + +bool T2K_CC0piWithProtons_2018_multidif_0p_1p_Np::isSignal(FitEvent *event) { + // If looking at all subsamples, only requirement is that this is a CC0pi event + if (useCC0pi0p && useCC0pi1p && useCC0piNp) + return SignalDef::isCC0pi(event, 14, EnuMin, EnuMax); + + // Otherwise, evaluate each relevant signal definition separately + if (useCC0pi0p && SignalDef::isT2K_CC0pi0p(event, EnuMin, EnuMax)) + return true; + if (useCC0pi1p && SignalDef::isT2K_CC0pi1p(event, EnuMin, EnuMax)) + return true; + if (useCC0piNp && SignalDef::isT2K_CC0piNp(event, EnuMin, EnuMax)) + return true; + + // If you get here, the event has failed and is not signal + return false; +}; + +void T2K_CC0piWithProtons_2018_multidif_0p_1p_Np::FillEventVariables( + FitEvent *event) { + + if (event->NumFSParticle(13) == 0) + return; + + TLorentzVector Pnu = event->GetNeutrinoIn()->fP; + TLorentzVector Pmu = event->GetHMFSParticle(13)->fP; + + double pmu = Pmu.Vect().Mag()/1000.; + double CosThetaMu = cos(Pnu.Vect().Angle(Pmu.Vect())); + + // Dummy proton variables if we don't have a proton + double pp = -999; + double CosThetaP = -999; + // Check if we do have a proton and fill variables + if (event->NumFSParticle(2212) > 0){ + TLorentzVector Pp = event->GetHMFSParticle(2212)->fP; + pp = Pp.Vect().Mag() / 1000.; + CosThetaP = cos(Pnu.Vect().Angle(Pp.Vect())); + } + + // How many protons above threshold? + std::vector protons = event->GetAllFSProton(); + int nProtonsAboveThresh = 0; + for (size_t i = 0; i < protons.size(); i++) { + if (protons[i]->p() > 500) + nProtonsAboveThresh++; + } + + // Get bin number in total 1D histogram + int binnumber = Get1DBin(nProtonsAboveThresh, pmu, CosThetaMu, pp, CosThetaP); + + // I'm hacking this to fit in the Measurement1D framework, but it's going to be super ugly - apologies... + // The 1D histogram handled by NUISANCE is defined in terms of bin number, so that has to be fXVar + fXVar = binnumber; + + // Then fill variables so I can use them to fill the slice histograms in FillHistograms() + fPP = pp; + fPMu = pmu; + fCosThetaP = CosThetaP; + fCosThetaMu = CosThetaMu; + fNp = nProtonsAboveThresh; + + return; +}; + +void T2K_CC0piWithProtons_2018_multidif_0p_1p_Np::FillHistograms() { + + Measurement1D::FillHistograms(); + if (Signal) { + // fMCHist_Fine2D->Fill(fXVar, fYVar, Weight); + if (useCC0pi0p && fNp == 0){ + fMCHist_CC0pi0pCosTheta->Fill(fCosThetaMu, Weight); + } + else if (useCC0pi1p && fNp == 1){ + fMCHist_CC0pi1pCosTheta->Fill(fCosThetaMu, Weight); + } + FillMCSlice(fNp, fPMu, fCosThetaMu, fPP, fCosThetaP, Weight); + } +} + +// Don't implement this for now - copied from T2K_CC0pi1p_XSec_3DPcoscos_nu +// // Modification is needed after the full reconfigure to move bins around +// // Otherwise this would need to be replaced by a TH2Poly which is too awkward. +// void T2K_CC0piWithProtons_2018_multidif_0p_1p_Np::ConvertEventRates() { +// +// for (int i = 0; i < 4; i++) { +// fMCHist_Slices[i]->GetSumw2(); +// } +// +// // Do standard conversion. +// Measurement1D::ConvertEventRates(); +// +// // First scale MC slices also by their width in Y and Z +// // MCHist_Slices[0]->Scale(1.0 / 1.00); +// // MCHist_Slices[1]->Scale(1.0 / 0.60); +// // MCHist_Slices[2]->Scale(1.0 / 0.10); +// // MCHist_Slices[3]->Scale(1.0 / 0.10); +// +// // Now Convert into 1D list +// fMCHist->Reset(); +// int bincount = 0; +// for (int i = 0; i < 4; i++) { +// for (int j = 0; j < fDataHist_Slices[i]->GetNumberOfBins(); j++) { +// fMCHist->SetBinContent(bincount + 1, +// fMCHist_Slices[i]->GetBinContent(j + 1)); +// // fMCHist->SetBinError(bincount+1, fMCHist_Slices[i]->GetBinError(j+1)); +// bincount++; +// } +// } +// +// return; +// } + +void T2K_CC0piWithProtons_2018_multidif_0p_1p_Np::FillMCSlice(int nProtonsAboveThresh, double pmu, double CosThetaMu, double pp, double CosThetaP, double w) { +// Get slice number for 1D CosThetaMu slice + int CosThetaMuSliceNo = GetCosThetaMuSlice(nProtonsAboveThresh, CosThetaMu); + // If sliceno is valid (not negative), fill the relevant slice + if (CosThetaMuSliceNo < 0) return; + // CC0pi0p slices: fill with pmu + if (useCC0pi0p && nProtonsAboveThresh == 0 && CosThetaMuSliceNo < 10){ + fMCHist_Slices[CosThetaMuSliceNo]->Fill(pmu, w); + } + // CC0pi1p slices: fill with CosThetaP + if (useCC0pi1p && nProtonsAboveThresh == 1){ + fMCHist_Slices[CosThetaMuSliceNo]->Fill(CosThetaP, w); + + // If we're looking at CC0pi1p, also fill the CosThetaMu-CosThetaP slices with PP + int CC0pi1p2DSliceNo = GetCC0pi1p2DSlice(nProtonsAboveThresh, CosThetaMu, CosThetaP); + if (CC0pi1p2DSliceNo < 0) return; + fMCHist_Slices[CC0pi1p2DSliceNo]->Fill(pp, w); + } +} + +void T2K_CC0piWithProtons_2018_multidif_0p_1p_Np::SetHistograms() { + + // Read in 1D Data Histograms + fInputFile = new TFile( + (FitPar::GetDataBase() + "/T2K/CC0pi/STV/multidif_results.root") + .c_str(), + "READ"); + // fInputFile->ls(); + + // Read in 1D Data + TH1D *tempDataHist = (TH1D *)fInputFile->Get("Result"); + + // Read in covariance matrix + TH2D *tempcov = (TH2D *)fInputFile->Get("CovarianceMatrix"); + + // The input data and covariance matrix include bins for CC0pi0p, CC0pi1p, and CC0piNp. We may not want them all if we only want to look at one or two of the sub-samples, so go through and only keep the bins we want + // CC0pi0p: bins 1-60 -> 60 bins + // CC0pi1p: bins 61-92 -> 32 bins + // CC0piNp: bin 93 -> 1 bin + int n_binskeep = 0; + if (useCC0pi0p) n_binskeep += 60; + if (useCC0pi1p) n_binskeep += 32; + if (useCC0piNp) n_binskeep += 1; + fDataHist = new TH1D("DataHist", tempDataHist->GetTitle().c_str(),n_binskeep,0,n_binskeep); + fFullCovar = new TMatrixDSym(n_binskeep); + + int i_binskeep = 1; + for (int i_allbins=1; i_allbinsGetNbinsX()+1; i_allbins++){ + if ((i_allbins >=1 && i_allbins <=60) && !useCC0pi0p) continue; + if ((i_allbins >= 61 && i_allbins <=92) && !useCC0pi1p) continue; + if ((i_allbins == 93) && !useCC0piNp) continue; + + fDataHist->SetBinContent(i_binskeep,tempDataHist->GetBinContent(i_allbins)); + fDataHist->SetBinError(i_binskeep,tempDataHist->GetBinError(i_allbins)); + + int j_binskeep = 1; + for (int j_allbins = 1; j_allbins < tempcov->GetNbinsX()+1; j_allbins++){ + if ((j_allbins >=1 && j_allbins <=60) && !useCC0pi0p) continue; + if ((j_allbins >= 61 && j_allbins <=92) && !useCC0pi1p) continue; + if ((j_allbins == 93) && !useCC0piNp) continue; + + (*fFullCovar)(i_binskeep-1,j_binskeep-1) = tempcov->GetBinContent(i_allbins, j_allbins); + } // end loop over j_allbins + + i_binskeep++; + } + + covar = StatUtils::GetInvert(fFullCovar); + fDecomp = StatUtils::GetDecomp(fFullCovar); + + // Make 1D MC histogram + TH1D *linearResult = new TH1D(*fDataHist); + linearResult->SetName("T2K_CC0piWithProtons_2018_multidif_0p_1p_Np_data"); + SetAutoProcessTH1(linearResult, kCMD_Write); + + + // Fine histograms - don't implement for now (this is copied from T2K_CC0pi1p_XSec_3DPcoscos_nu) + // fMCHist_Fine2D = new TH2D("T2K_CC0piWithProtons_2018_multidif_0p_1p_Np_Fine2D", + // "T2K_CC0piWithProtons_2018_multidif_0p_1p_Np_Fine2D", + // 400, 0.0, 30.0, 100, -1.0, 1.0); + // SetAutoProcessTH1(fMCHist_Fine2D); + + + // The code below converts a relative covariance matrix to an absolute one. I think the input is absolute so we don't need it, but come back to this if the results look weird + // for (int i = 0; i < fDataHist->GetNbinsX(); i++) { + // for (int j = 0; j < fDataHist->GetNbinsX(); j++) { + // //(*fFullCovar)(i,j) = tempcov->GetBinContent(i+1, j+1); + // (*fFullCovar)(i, j) = tempcov->GetBinContent(i + 1, j + 1) * + // fDataHist->GetBinContent(i + 1) * + // fDataHist->GetBinContent(j + 1); + // if (i == j) + // fDataHist->SetBinError(i + 1, sqrt((*fFullCovar)(i, j))); + // // if(i==j) std::cout << "For bin " << i+1 << ", relative covariance was " + // // << tempcov->GetBinContent(i+1,j+1); if(i==j) std::cout << ". Absolute + // // covariance is now " << (*fFullCovar)(i,j) << ", linear xsec is: " << + // // fDataHist->GetBinContent(i+1) << std::endl; + // } + // } + + + // Read in data slice histograms and make MC slices + // Slices are stored in slightly different ways for 0p and 1p samples + // CC0pi0p: folder NoProtonsAbove500MeV + // -> TH1D ResultInMuonCosTheta + // -> TH1D MuonCosThetaSlice_i where i = 0 to 9 + // CC0pi1p: folder OneProtonAbove500MeV + // -> TH1D ResultInMuonCosTheta + // -> TH1D MuonCosThetaSlice_1D_i where i = 0 to 3 + // -> TH1D MuCThSlice_1_PCthSlice_0 + // -> TH1D MuCThSlice_2_PCthSlice_0 + // -> TH1D MuCThSlice_2_PCthSlice_1 + // -> TH1D MuCThSlice_3_PCthSlice_0 + // We also have proton multiplicity, in folder ProtonMultiplicity + // -> TH1D Result + // -> TH2D CovarianceMatrix + // TODO add this as a separate sample? Don't implement here + + // CC0pi0p slices + if (useCC0pi0p){ + fDataHist_CC0pi0pCosTheta = (TH1D*)fInputFile->Get("NoProtonsAbove500MeV/ResultInMuonCosTheta")->Clone(); + fMCHist_CC0pi0pCosTheta = fDataHist_CC0pi0pCosTheta->Clone("T2K_2018_CC0pi0p_ResultInMuonCosTheta"); + fMCHist_CC0pi0pCosTheta->Reset(); + SetAutoProcessTH1(fDataHist_CC0pi0pCosTheta, kCMD_Write); + SetAutoProcessTH1(fMCHist_CC0pi0pCosTheta, kCMD_Reset); + + for (int i=0; i<=9; i++){ + fDataHist_Slices.push_back((TH1D*)fInputFile->Get(Form("NoProtonsAbove500MeV/MuonCosThetaSlice_%i", i))->Clone(Form("T2K_2018_CC0pi0p_Data_Slice%i", i))); + fMCHist_Slices.push_back((TH1D*)fDataHist_Slices[i]->Clone(Form("T2K_2018_CC0pi0p_MC_Slice%i", i))); + } // end loop over i + } + + + // CC0pi1p slices + if (useCC0pi1p){ + fDataHist_CC0pi1pCosTheta = (TH1D*)fInputFile->Get("OneProtonAbove500MeV/ResultInMuonCosTheta")->Clone(); + fMCHist_CC0pi1pCosTheta = fDataHist_CC0pi1pCosTheta->Clone("T2K_2018_CC0pi1p_ResultInMuonCosTheta"); + fMCHist_CC0pi1pCosTheta->Reset(); + SetAutoProcessTH1(fDataHist_CC0pi1pCosTheta, kCMD_Write); + SetAutoProcessTH1(fMCHist_CC0pi1pCosTheta, kCMD_Reset); + + for (int i=0; i<=3; i++){ + fDataHist_Slices.push_back((TH1D*)fInputFile->Get(Form("OneProtonAbove500MeV/MuonCosThetaSlice_1D_%i", i))->Clone(Form("T2K_2018_CC0pi1p_Data_MuonCosTh1DSlice%i", i))); + fMCHist_Slices.push_back((TH1D*)fDataHist_Slices[i]->Clone(Form("T2K_2018_CC0pi1p_MC_MuonCosTh1DSlice%i", i))); + } + // Add in the muon costh-p costh slices (which aren't as nicely numbered) + // -> TH1D MuCThSlice_1_PCthSlice_0 + fDataHist_Slices.push_back((TH1D*)fInputFile->Get("MuCThSlice_1_PCthSlice_0")->Clone("T2K_2018_CC0pi1p_Data_MuCThSlice_1_PCthSlice_0")); + fMCHist_Slices.push_back((TH1D*)fDataHist_Slices[i]->Clone("T2K_2018_CC0pi1p_MC_MuCThSlice_1_PCthSlice_0")); + // -> TH1D MuCThSlice_2_PCthSlice_0 + fDataHist_Slices.push_back((TH1D*)fInputFile->Get("MuCThSlice_2_PCthSlice_0")->Clone("T2K_2018_CC0pi1p_Data_MuCThSlice_2_PCthSlice_0")); + fMCHist_Slices.push_back((TH1D*)fDataHist_Slices[i]->Clone("T2K_2018_CC0pi1p_MC_MuCThSlice_2_PCthSlice_0")); + // -> TH1D MuCThSlice_2_PCthSlice_1 + fDataHist_Slices.push_back((TH1D*)fInputFile->Get("MuCThSlice_2_PCthSlice_1")->Clone("T2K_2018_CC0pi1p_Data_MuCThSlice_2_PCthSlice_1")); + fMCHist_Slices.push_back((TH1D*)fDataHist_Slices[i]->Clone("T2K_2018_CC0pi1p_MC_MuCThSlice_2_PCthSlice_1")); + // -> TH1D MuCThSlice_3_PCthSlice_0 + fDataHist_Slices.push_back((TH1D*)fInputFile->Get("MuCThSlice_3_PCthSlice_0")->Clone("T2K_2018_CC0pi1p_Data_MuCThSlice_3_PCthSlice_0")); + fMCHist_Slices.push_back((TH1D*)fDataHist_Slices[i]->Clone("T2K_2018_CC0pi1p_MC_MuCThSlice_3_PCthSlice_0")); + } + + + // Set all slice histograms to auto-process and reset MC histograms + for (int i=0; iClearBinContents(); + SetAutoProcessTH1(fDataHist_Slices[i], kCMD_Write); + SetAutoProcessTH1(fMCHist_Slices[i], kCMD_Reset); + } + return; +}; + +// Yay hardcoding +// Taken from multidif_binMap.txt in data release +int T2K_CC0piWithProtons_2018_multidif_0p_1p_Np::Get1DBin(int nProtonsAboveThresh, double pmu, double CosThetaMu, double pp, double CosThetaP) { + + int binnumber = -999; + + // Calculate bin number (check that we want to use this sample before looking for the correct bin to save computation) + if (nProtonsAboveThresh == 0 && useCC0pi0p){ //CC0pi0p: 2D binning in CosThetaMu--pmu + if (CosThetaMu >= -1 && CosThetaMu <= -0.3) {binnumber = 1; break;} + else if (CosThetaMu > -0.3 && CosThetaMu <= 0.3){ + // Now find bin in pmu + if (pmu >= 0 && pmu <= 0.3) {binnumber = 2; break;} + else if (pmu > 0.3 && pmu <= 0.4){binnumber = 3; break;} + else if (pmu > 0.4 && pmu <= 30){binnumber = 4; break;} + } // end if (CosThetaMu > -0.3 && CosThetaMu <= 0.3) + else if (CosThetaMu > 0.3 && CosThetaMu <= 0.6){ + // Now find bin in pmu + if (pmu >= 0 && pmu <= 0.3){binnumber = 5; break;} + else if (pmu > 0.3 && pmu <= 0.4){binnumber = 6; break;} + else if (pmu > 0.4 && pmu <= 0.5){binnumber = 7; break;} + else if (pmu > 0.5 && pmu <= 0.6){binnumber = 8; break;} + else if (pmu > 0.6 && pmu <= 30){binnumber = 9; break;} + } // end if (CosThetaMu > 0.3 && CosThetaMu <= 0.6) + else if (CosThetaMu > 0.6 && CosThetaMu <= 0.7){ + // Now find bin in pmu + if (pmu >= 0 && pmu <= 0.3){binnumber = 10; break;} + else if (pmu > 0.3 && pmu <= 0.4){binnumber = 11; break;} + else if (pmu > 0.4 && pmu <= 0.5){binnumber = 12; break;} + else if (pmu > 0.5 && pmu <= 0.6){binnumber = 13; break;} + else if (pmu > 0.6 && pmu <= 30){binnumber = 14; break;} + } // end if (CosThetaMu > 0.6 && CosThetaMu <= 0.7) + else if (CosThetaMu > 0.7 && CosThetaMu <= 0.8){ + // Now find bin in pmu + if (pmu >= 0 && pmu <= 0.3){binnumber = 15; break;} + else if (pmu > 0.3 && pmu <= 0.4){binnumber = 16; break;} + else if (pmu > 0.4 && pmu <= 0.5){binnumber = 17; break;} + else if (pmu > 0.5 && pmu <= 0.6){binnumber = 18; break;} + else if (pmu > 0.6 && pmu <= 0.7){binnumber = 19; break;} + else if (pmu > 0.7 && pmu <= 0.8){binnumber = 20; break;} + else if (pmu > 0.8 && pmu <= 30){binnumber = 21; break;} + } // end if (CosThetaMu > 0.7 && CosThetaMu <= 0.8) + else if (CosThetaMu > 0.8 && CosThetaMu <= 0.85){ + // Now find bin in pmu + if (pmu >= 0 && pmu <= 0.4){binnumber = 22; break;} + else if (pmu > 0.4 && pmu <= 0.5){binnumber = 23; break;} + else if (pmu > 0.5 && pmu <= 0.6){binnumber = 24; break;} + else if (pmu > 0.6 && pmu <= 0.7){binnumber = 25; break;} + else if (pmu > 0.7 && pmu <= 0.8){binnumber = 26; break;} + else if (pmu > 0.8 && pmu <= 30){binnumber = 27; break;} + } // end if (CosThetaMu > 0.8 && CosThetaMu <= 0.85) + else if (CosThetaMu > 0.85 && CosThetaMu <= 0.9){ + // Now find bin in pmu + if (pmu >= 0 && pmu <= 0.3){binnumber = 28; break;} + else if (pmu > 0.3 && pmu <= 0.4){binnumber = 29; break;} + else if (pmu > 0.4 && pmu <= 0.5){binnumber = 30; break;} + else if (pmu > 0.5 && pmu <= 0.6){binnumber = 31; break;} + else if (pmu > 0.6 && pmu <= 0.7){binnumber = 32; break;} + else if (pmu > 0.7 && pmu <= 0.8){binnumber = 33; break;} + else if (pmu > 0.8 && pmu <= 1){binnumber = 34; break;} + else if (pmu > 1 && pmu <= 30){binnumber = 35; break;} + } // end if (CosThetaMu > 0.85 && CosThetaMu <= 0.9) + else if (CosThetaMu > 0.9 && CosThetaMu <= 0.94){ + // Now find bin in pmu + if (pmu >= 0 && pmu <= 0.4){binnumber = 36; break;} + else if (pmu > 0.4 && pmu <= 0.5){binnumber = 37; break;} + else if (pmu > 0.5 && pmu <= 0.6){binnumber = 38; break;} + else if (pmu > 0.6 && pmu <= 0.7){binnumber = 39; break;} + else if (pmu > 0.7 && pmu <= 0.8){binnumber = 40; break;} + else if (pmu > 0.8 && pmu <= 1.25){binnumber = 41; break;} + else if (pmu > 1.25 && pmu <= 30){binnumber = 42; break;} + } // end if (CosThetaMu > 0.9 && CosThetaMu <= 0.94) + else if (CosThetaMu > 0.94 && CosThetaMu <= 0.98){ + // Now find bin in pmu + if (pmu >= 0 && pmu <= 0.4){binnumber = 43; break;} + else if (pmu > 0.4 && pmu <= 0.5){binnumber = 44; break;} + else if (pmu > 0.5 && pmu <= 0.6){binnumber = 45; break;} + else if (pmu > 0.6 && pmu <= 0.7){binnumber = 46; break;} + else if (pmu > 0.7 && pmu <= 0.8){binnumber = 47; break;} + else if (pmu > 0.8 && pmu <= 1){binnumber = 48; break;} + else if (pmu > 1 && pmu <= 1.25){binnumber = 49; break;} + else if (pmu > 1.25 && pmu <= 1.5){binnumber = 50; break;} + else if (pmu > 1.5 && pmu <= 2){binnumber = 51; break;} + else if (pmu > 2 && pmu <= 30){binnumber = 52; break;} + } // end if (CosThetaMu > 0.94 && CosThetaMu <= 0.98) + else if (CosThetaMu > 0.98 && CosThetaMu <= 1){ + // Now find bin in pmu + if (pmu >= 0 && pmu <= 0.5){binnumber = 53; break;} + else if (pmu > 0.5 && pmu <= 0.65){binnumber = 54; break;} + else if (pmu > 0.65 && pmu <= 0.8){binnumber = 55; break;} + else if (pmu > 0.8 && pmu <= 1.25){binnumber = 56; break;} + else if (pmu > 1.25 && pmu <= 2){binnumber = 57; break;} + else if (pmu > 2 && pmu <= 3){binnumber = 58; break;} + else if (pmu > 3 && pmu <= 5){binnumber = 59; break;} + else if (pmu > 5 && pmu <= 30){binnumber = 60; break;} + } // end if (CosThetaMu > 0.98 && CosThetaMu <= 1) + } // end (nProtonsAboveThresh == 0) + else if (nProtonsAboveThresh == 1 && useCC0pi1p){ + if (CosThetaMu >= -1 && CosThetaMu <= -0.3){ + // Find bin in CosThetaP + if (CosThetaP >= -1 && CosThetaP <= 0.87){binnumber = 61; break;} + else if (CosThetaP > 0.87 && CosThetaP <= 0.94){binnumber = 62; break;} + else if (CosThetaP > 0.94 && CosThetaP <= 0.97){binnumber = 63; break;} + else if (CosThetaP > 0.97 && CosThetaP <= 1){binnumber = 64; break;} + } // end (CosThetaMu >= -1 && CosThetaMu <= -0.3) + else if (CosThetaMu > -0.3 && CosThetaMu <= 0.3){ + // Find bin in CosThetaP + if (CosThetaP >= -1 && CosThetaP <= 0.75){binnumber = 65; break;} + else if (CosThetaP > 0.75 && CosThetaP <= 0.85){binnumber = 66; break;} + else if (CosThetaP > 0.85 && CosThetaP <= 0.94){ + // Find bin in pp + if (pp >= 0.5 && pp <= 0.68){binnumber = 67; break;} + else if (pp > 0.68 && pp <= 0.78){binnumber = 68; break;} + else if (pp > 0.78 && pp <= 0.9){binnumber = 69; break;} + else if (pp > 0.9 && pp <= 30){binnumber = 70; break;} + } // end if (CosThetaP > 0.85 && CosThetaP <= 0.94) + else if (CosThetaP > 0.94 && CosThetaP <= 1){binnumber = 71; break;} + } // end if (CosThetaMu > -0.3 && CosThetaMu <= 0.3) + else if (CosThetaMu > 0.3 && CosThetaMu <= 0.8){ + // Find bin in CosThetaP + if (CosThetaP >= -1 && CosThetaP <= 0.3){binnumber = 72; break;} + else if (CosThetaP > 0.3 && CosThetaP <= 0.5){binnumber = 73; break;} + else if (CosThetaP > 0.5 && CosThetaP <= 0.8){ + // find bin in pp + if (pp >= 0.5 && pp <= 0.6){binnumber = 74; break;} + else if (pp > 0.6 && pp <= 0.7){binnumber = 75; break;} + else if (pp > 0.7 && pp <= 0.8){binnumber = 76; break;} + else if (pp > 0.8 && pp <= 0.9){binnumber = 77; break;} + else if (pp > 0.9 && pp <= 30){binnumber = 78; break;} + } // end if (CosThetaP > 0.5 && CosThetaP <= 0.8) + else if (CosThetaP > 0.8 && CosThetaP <= 1){ + // find bin in pp + if (pp >= 0.5 && pp <= 0.6){return 79;} + else if (pp > 0.6 && pp <= 0.7){binnumber = 80; break;} + else if (pp > 0.7 && pp <= 0.8){binnumber = 81; break;} + else if (pp > 0.8 && pp <= 1){binnumber = 82; break;} + else if (pp > 1 && pp <= 30){binnumber = 83; break;} + } // end if (CosThetaP > 0.8 && CosThetaP <= 1) + } // end if (CosThetaMu > 0.3 && CosThetaMu <= 0.8) + else if (CosThetaMu > 0.8 && CosThetaMu <= 1){ + // Find bin in CosThetaP + if (CosThetaP >= -1 && CosThetaP <= 0){binnumber = 84; break;} + else if (CosThetaP > 0 && CosThetaP <= 0.3){binnumber = 85; break;} + else if (CosThetaP > 0.3 && CosThetaP <= 0.8){ + // find bin in pp + if (pp >= 0.5 && pp <= 0.6){binnumber = 86; break;} + else if (pp > 0.6 && pp <= 0.7){binnumber = 87; break;} + else if (pp > 0.7 && pp <= 0.8){binnumber = 88; break;} + else if (pp > 0.8 && pp <= 0.9){binnumber = 89; break;} + else if (pp > 0.9 && pp <= 1.1){binnumber = 90; break;} + else if (pp > 1.1 && pp <= 30){binnumber = 91; break;} + } // end if (CosThetaP > 0.3 && CosThetaP <= 0.8) + else if (CosThetaP > 0.8 && CosThetaP <= 1){binnumber = 92; break} + } // end if (CosThetaMu > 0.8 && CosThetaMu <= 1) + } // end if (nProtonsAboveThresh == 1) + else if (nProtonsAboveThresh > 1 && useCC0piNp){ + binnumber = 93; + } + + // If binnumber is still -999, something has gone wrong + if (binnumber == -999){ + std::cout << "ERROR did not find correct 1D bin for an event with nProtonsAboveThresh = " << nProtonsAboveThresh << ", pmu = " << pmu << ", CosThetaMu = " << CosThetaMu << ", pp = " << pp << ", CosThetaP = " << CosThetaP << std::endl; + return -999; + } + + // Now need to work out adjustments for if we don't want to use one or more of the samples + // If all samples are wanted, no adjustments required - the binning stands + // If only CC0pi0p sample is wanted, no adjustments required - CC0pi0p binning is correct and the other bins won't be assigned + // If CC0pi0p and CC0pi1p samples are wanted, no adjustments required + + // If only CC0pi1p sample is wanted or CC0pi1p and CC0piNp samples are wanted, need to subtract off total number of CC0pi0p bins + if (useCC0pi1p && !useCC0pi0p){ + binnumber -= 60; + } + // If only CC0piNp sample is wanted, need to subtract off total number of CC0pi0p and CC0pi1p bins + if (useCC0piNp && !useCC0pi0p && !useCC0pi1p){ + binnumber =- 92; + } + // If CC0pi0p and CC0piNp samples are wanted, need to subtract off CC0pi1p bins from CC0piNp bin number + if (useCC0piNp && useCC0pi0p && !useCC0pi1p && binnumber = 93){ + binnumber = 61; + } + + return binnumber; +}; + +// Yay hardcoding again! +// Taken from multidif_binMap.txt in data release +int T2K_CC0piWithProtons_2018_multidif_0p_1p_Np::GetCosThetaMuSlice(int nProtonsAboveThresh, double CosThetaMu) { + + int slicenumber = -999; + + // if (useCC0pi0p), the first 10 slices will be CC0pi0p slices + // -> TH1D MuonCosThetaSlice_i where i = 0 to 9 + // -> slices 0 to 9 + // if (useCC0pi1p), the next 8 slices will be CC0pi1p slices + // -> TH1D MuonCosThetaSlice_1D_i where i = 0 to 3 + // -> TH1D MuCThSlice_1_PCthSlice_0 + // -> TH1D MuCThSlice_2_PCthSlice_0 + // -> TH1D MuCThSlice_2_PCthSlice_1 + // -> TH1D MuCThSlice_3_PCthSlice_0 + // -> slices 10 to 17 if (useCC0pi0p) + // -> slices 0 to 7 if (!useCC0pi0p) + // In this code, we just find the costhetamu slices + + // Calculate bin number (check that we want to use this sample before looking for the correct bin to save computation) + if (nProtonsAboveThresh == 0 && useCC0pi0p){ //CC0pi0p: 10 slices in CosThetaMu + if (CosThetaMu >= -1 && CosThetaMu <= -0.3) {slicenumber = 0; break;} + else if (CosThetaMu > -0.3 && CosThetaMu <= 0.3){slicenumber = 1; break;} + else if (CosThetaMu > 0.3 && CosThetaMu <= 0.6){slicenumber = 2; break;} + else if (CosThetaMu > 0.6 && CosThetaMu <= 0.7){slicenumber = 3; break;} + else if (CosThetaMu > 0.7 && CosThetaMu <= 0.8){slicenumber = 4; break;} + else if (CosThetaMu > 0.8 && CosThetaMu <= 0.85){slicenumber = 5; break;} + else if (CosThetaMu > 0.85 && CosThetaMu <= 0.9){slicenumber = 6; break;} + else if (CosThetaMu > 0.9 && CosThetaMu <= 0.94){slicenumber = 7; break;} + else if (CosThetaMu > 0.94 && CosThetaMu <= 0.98){slicenumber = 8; break;} + else if (CosThetaMu > 0.98 && CosThetaMu <= 1){slicenumber = 9; break;} + } // end (nProtonsAboveThresh == 0) + else if (nProtonsAboveThresh == 1 && useCC0pi1p){ // CC0pi1p: 8 slices, either in CosThetaMu or CosThetaMu-CosThetaP, depending on the slice + if (CosThetaMu >= -1 && CosThetaMu <= -0.3){slicenumber = 10; break;} + else if (CosThetaMu > -0.3 && CosThetaMu <= 0.3){slicenumber = 11; break;} + else if (CosThetaMu > 0.3 && CosThetaMu <= 0.8){slicenumber = 12; break;} + else if (CosThetaMu > 0.8 && CosThetaMu <= 1){slicenumber = 13; break;} + } // end if (nProtonsAboveThresh == 1) + + // slicenumber may be -999 if nProtonsAboveThresh > 1 + // otherwise, something has gone wrong + if (slicenumber == -999 && nProtonsAboveThresh <= 1){ + std::cout << "ERROR did not find correct 1D CosThetaMu slice for an event with nProtonsAboveThresh = " << nProtonsAboveThresh << ", pmu = " << pmu << ", CosThetaMu = " << CosThetaMu << std::endl; + return -999; + } + + // If useCC0pi0p is false, adjust slice numbers for CC0pi1p + if (useCC0pi1p && !useCC0pi0p){ + slicenumber =- 10; + } + + return slicenumber; +}; + +// Hardcoding one more time +// Taken from multidif_binMap.txt in data release +int T2K_CC0piWithProtons_2018_multidif_0p_1p_Np::GetCC0pi1p2DSlice(int nProtonsAboveThresh, double CosThetaMu, double CosThetaP) { + + int slicenumber = -999; + + // Calculate slice number (note: only 2D CC0pi1p slices in CosThetaMu-CosThetaP are handled here) + // These slices exist for: + // CosThetaMu -0.3 - 0.3, CosThetaP 0.85 - 0.94 + // CosThetaMu 0.3 - 0.8, CosThetaP 0.5 - 0.8 + // CosThetaMu 0.3 - 0.8, CosThetaP 0.8 - 1.0 + // CosThetaMu 0.8 - 1.0, CosThetaP 0.3 - 0.8 + // If useCC0pi0p is true, these will be slices 14-17 + // If useCC0pi0p is false, these will be slices 4-7 + if (nProtonsAboveThresh == 1 && useCC0pi1p){ + if (CosThetaMu > -0.3 && CosThetaMu <= 0.3 && CosThetaP > 0.85 && CosThetaP <= 0.94){slicenumber = 14; break;} + else if (CosThetaMu > 0.3 && CosThetaMu <= 0.8){ + // Find bin in CosThetaP + if (CosThetaP > 0.5 && CosThetaP <= 0.8){slicenumber = 15; break;} + else if (CosThetaP > 0.8 && CosThetaP <= 1){slicenumber = 16; break;} + } // end if (CosThetaMu > 0.3 && CosThetaMu <= 0.8) + else if (CosThetaMu > 0.8 && CosThetaMu <= 1 && CosThetaP > 0.3 && CosThetaP <= 0.8){slicenumber = 17; break;} + } // end if (nProtonsAboveThresh == 1) + + // No check on binnumber = -999 here, because many events won't fall into one of these slices + + // If useCC0pi0p is false, adjust slice numbers for CC0pi1p + if (useCC0pi1p && !useCC0pi0p){ + slicenumber =- 10; + } + + return binnumber; +}; diff --git a/src/T2K/T2K_CC0pi1p_XSec_3DPcoscos_nu_nonuniform.h b/src/T2K/T2K_CC0piWithProtons_XSec_2018_multidif_0p_1p_Np.h similarity index 62% rename from src/T2K/T2K_CC0pi1p_XSec_3DPcoscos_nu_nonuniform.h rename to src/T2K/T2K_CC0piWithProtons_XSec_2018_multidif_0p_1p_Np.h index 3124225..d9bdebd 100644 --- a/src/T2K/T2K_CC0pi1p_XSec_3DPcoscos_nu_nonuniform.h +++ b/src/T2K/T2K_CC0piWithProtons_XSec_2018_multidif_0p_1p_Np.h @@ -1,76 +1,83 @@ // Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret /******************************************************************************* * This file is part of NUISANCE. * * NUISANCE is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * NUISANCE is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with NUISANCE. If not, see . *******************************************************************************/ -#ifndef T2K_CC0PI1P_3DPCOSCOS_NU_NONUNIFORM_H_SEEN -#define T2K_CC0PI1P_3DPCOSCOS_NU_NONUNIFORM_H_SEEN +#ifndef T2K_CC0PIWITHPROTONS_2018_MULTIDIF_0P_1P_NP_H_SEEN +#define T2K_CC0PIWITHPROTONS_2018_MULTIDIF_0P_1P_NP_H_SEEN #include "Measurement1D.h" -#include "TH2Poly.h" -class T2K_CC0pi1p_XSec_3DPcoscos_nu_nonuniform : public Measurement1D { +class T2K_CC0piWithProtons_2018_multidif_0p_1p_Np : public Measurement1D { public: /// Basic Constructor. - T2K_CC0pi1p_XSec_3DPcoscos_nu_nonuniform(nuiskey samplekey); + T2K_CC0piWithProtons_2018_multidif_0p_1p_Np(nuiskey samplekey); /// Virtual Destructor - ~T2K_CC0pi1p_XSec_3DPcoscos_nu_nonuniform() {}; + ~T2K_CC0piWithProtons_2018_multidif_0p_1p_Np() {}; /// Numu CC0PI Signal Definition /// - /// /item + /// /item bool isSignal(FitEvent *nvect); /// Read histograms in a special way because format is different. /// Read from FitPar::GetDataBase()+"/T2K/CC0pi/T2K_CC0PI_2DPmuCosmu_Data.root" void SetHistograms(); /// Bin Tmu CosThetaMu void FillEventVariables(FitEvent* customEvent); // Fill Histograms void FillHistograms(); /// Have to do a weird event scaling for analysis 1 void ConvertEventRates(); private: + bool useCC0pi0p = false; + bool useCC0pi1p = false; + bool useCC0piNp = false; + bool only_allowed_particles; bool numu_event; double numu_energy; int particle_pdg; - double pmu, CosThetaMu; int fAnalysis; + double fPP, fCosThetaP, fCosThetaMu; bool fIsSystCov, fIsStatCov, fIsNormCov; - std::vector fDataPoly; - std::vector fMCPoly; - TFile* fInputFile; TH2D* fMCHist_Fine2D; - std::vector fMCHist_Slices; - std::vector fDataHist_Slices; + TH1D fMCHist_CC0pi0pCosTheta; + TH1D fDataHist_CC0pi0pCosTheta; + TH1D fMCHist_CC0pi1pCosTheta; + TH1D fDataHist_CC0pi1pCosTheta; + std::vector fMCHist_Slices; + std::vector fDataHist_Slices; + + void FillMCSlice(int nProtonsAboveThresh, double pmu, double CosThetaMu, double pp, double CosThetaP, double w); + int Get1DBin(int nProtonsAboveThresh, double pmu, double CosThetaMu, double pp, double CosThetaP); + int GetCosThetaMuSlice(int nProtonsAboveThresh, double CosThetaMu); + int GetCC0pi1p2DSlice(int nProtonsAboveThresh, double CosThetaMu, double CosThetaP); - void FillMCSlice(double x, double y, double z, double w); - }; - + #endif diff --git a/src/T2K/T2K_SignalDef.cxx b/src/T2K/T2K_SignalDef.cxx index 0df7aaf..1fc01a6 100644 --- a/src/T2K/T2K_SignalDef.cxx +++ b/src/T2K/T2K_SignalDef.cxx @@ -1,325 +1,328 @@ // Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret /******************************************************************************* * This file is part of NUISANCE. * * NUISANCE is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * NUISANCE is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with NUISANCE. If not, see . *******************************************************************************/ #include "T2K_SignalDef.h" #include "FitUtils.h" namespace SignalDef { // T2K H2O signal definition // https://doi.org/10.1103/PhysRevD.97.012001 bool isCC1pip_T2K_PRD97_012001(FitEvent *event, double EnuMin, double EnuMax) { if (!isCC1pi(event, 14, 211, EnuMin, EnuMax)) return false; TLorentzVector Pnu = event->GetHMISParticle(14)->fP; TLorentzVector Pmu = event->GetHMFSParticle(13)->fP; TLorentzVector Ppip = event->GetHMFSParticle(211)->fP; double p_mu = FitUtils::p(Pmu) * 1000; double p_pi = FitUtils::p(Ppip) * 1000; double cos_th_mu = cos(FitUtils::th(Pnu, Pmu)); double cos_th_pi = cos(FitUtils::th(Pnu, Ppip)); if (p_mu <= 200 || p_pi <= 200 || cos_th_mu <= 0.3 || cos_th_pi <= 0.3) { return false; } return true; }; // ****************************************************** // T2K CC1pi+ CH analysis (Raquel's thesis) // Has different phase space cuts depending on if using Michel tag or not // // Essentially consists of two samples: one sample which has Michel e (which we // can't get pion direction from); this covers backwards angles quite well. // Measurements including this sample does not have include pion kinematics cuts // one sample which has PID in FGD and TPC // and no Michel e. These are mostly // forward-going so require a pion // kinematics cut // // Essentially, cuts are: // 1 negative muon // 1 positive pion // Any number of nucleons // No other particles in the final state // // https://arxiv.org/abs/1909.03936 bool isCC1pip_T2K_arxiv1909_03936(FitEvent *event, double EnuMin, double EnuMax, int pscuts) { // ****************************************************** if (!isCC1pi(event, 14, 211, EnuMin, EnuMax)) { return false; } TLorentzVector Pnu = event->GetHMISParticle(14)->fP; TLorentzVector Pmu = event->GetHMFSParticle(13)->fP; double cos_th_mu = cos(FitUtils::th(Pnu, Pmu)); if (pscuts == kMuonFwd) { return (cos_th_mu > 0); } double p_mu = FitUtils::p(Pmu) * 1000; if (pscuts & kMuonHighEff) { if ((cos_th_mu <= 0.2) || (p_mu <= 200)) { return false; } } TLorentzVector Ppip = event->GetHMFSParticle(211)->fP; double cos_th_pi = cos(FitUtils::th(Pnu, Ppip)); double p_pi = FitUtils::p(Ppip) * 1000; if ((pscuts & kPionFwd) && (cos_th_pi <= 0)) { return false; } if ((pscuts & kPionVFwd) && (cos_th_pi <= 0.2)) { return false; } if ((pscuts & kPionHighMom) && (p_pi <= 200)) { return false; } return true; }; bool isT2K_CC0pi(FitEvent *event, double EnuMin, double EnuMax, int ana) { // Require a numu CC0pi event if (!isCC0pi(event, 14, EnuMin, EnuMax)) return false; TLorentzVector Pnu = event->GetHMISParticle(14)->fP; TLorentzVector Pmu = event->GetHMFSParticle(13)->fP; double CosThetaMu = cos(Pnu.Vect().Angle(Pmu.Vect())); double p_mu = Pmu.Vect().Mag(); // If we're doing a restricted phase space, Analysis II asks for: // Cos(theta_mu) > 0.0 and p_mu > 200 MeV if (ana == kAnalysis_II) { if ((CosThetaMu < 0.0) || (p_mu < 200)) { return false; } } return true; } bool isT2K_CC0pi1p(FitEvent *event, double EnuMin, double EnuMax) { // Require a numu CC0pi event if (!isCC0pi(event, 14, EnuMin, EnuMax)) return false; // Require at least one FS proton if (event->NumFSParticle(2212) == 0) return false; TLorentzVector pnu = event->GetHMISParticle(14)->fP; TLorentzVector pmu = event->GetHMFSParticle(13)->fP; TLorentzVector pp = event->GetHMFSParticle(2212)->fP; // Proton phase space if (pp.Vect().Mag() < 500) { return false; } // Need exactly one proton with 500 MeV or more momentum std::vector protons = event->GetAllFSProton(); int nProtonsAboveThresh = 0; for (size_t i = 0; i < protons.size(); i++) { if (protons[i]->p() > 500) nProtonsAboveThresh++; } if (nProtonsAboveThresh != 1) return false; return true; } -// CC0pi antinu in the P0D https://arxiv.org/abs/1908.10249 -bool isT2K_CC0piAnuP0D(FitEvent *event, double EnuMin, double EnuMax) { - - // Require a anumu CC0pi event - if (!isCC0pi(event, -14, EnuMin, EnuMax)) - return false; - - TLorentzVector pnu = event->GetHMISParticle(-14)->fP; - TLorentzVector pmu = event->GetHMFSParticle(-13)->fP; - double Pmu = pmu.Vect().Mag(); - double CosThetaMu = cos(pnu.Vect().Angle(pmu.Vect())); - // Muon phase space - if (Pmu < 400 || Pmu > 3410) return false; - if (Pmu < 530 && Pmu>=400 && CosThetaMu<0.84) return false; - if (Pmu < 670 && Pmu>=530 && CosThetaMu<0.85) return false; - if (Pmu < 800 && Pmu>=670 && CosThetaMu<0.88) return false; - if (Pmu < 1000 &&Pmu>=800 && CosThetaMu<0.9) return false; - if (Pmu < 1380 && Pmu>=1000 && CosThetaMu<0.91) return false; - if (Pmu < 2010 && Pmu>=1380 && CosThetaMu<0.92) return false; - if (Pmu < 3410 && Pmu>=2010 && CosThetaMu<0.95) return false; - - return true; -} - bool isT2K_CC0piNp(FitEvent *event, double EnuMin, double EnuMax) { + // In this case, we specifically mean N>2, as defined in the T2K CC0pi 2018 paper // Require a numu CC0pi event if (!isCC0pi(event, 14, EnuMin, EnuMax)) return false; // Require at least one FS proton if (event->NumFSParticle(2212) == 0) return false; TLorentzVector pnu = event->GetHMISParticle(14)->fP; TLorentzVector pmu = event->GetHMFSParticle(13)->fP; TLorentzVector pp = event->GetHMFSParticle(2212)->fP; // Proton phase space if (pp.Vect().Mag() < 500) { return false; } // Need exactly one proton with 500 MeV or more momentum std::vector protons = event->GetAllFSProton(); int nProtonsAboveThresh = 0; for (size_t i = 0; i < protons.size(); i++) { if (protons[i]->p() > 500) nProtonsAboveThresh++; } - if (nProtonsAboveThresh < 2) + if (nProtonsAboveThresh <= 1) return false; return true; } bool isT2K_CC0pi0p(FitEvent *event, double EnuMin, double EnuMax) { // Require a numu CC0pi event if (!isCC0pi(event, 14, EnuMin, EnuMax)) return false; - // Require at least one FS proton + // Return true if no FS proton if (event->NumFSParticle(2212) == 0) - return false; + return true; + // Otherwise, check momentum of highest-momentum proton, and require it is below threshold TLorentzVector pnu = event->GetHMISParticle(14)->fP; TLorentzVector pmu = event->GetHMFSParticle(13)->fP; TLorentzVector pp = event->GetHMFSParticle(2212)->fP; // Proton phase space - if (pp.Vect().Mag() > 500) { - return false; + if (pp.Vect().Mag() < 500) { + return true; } + // If there are FS protons and the highest-momentum one is above threshold, return false + return false; +} + +// CC0pi antinu in the P0D https://arxiv.org/abs/1908.10249 +bool isT2K_CC0piAnuP0D(FitEvent *event, double EnuMin, double EnuMax) { + + // Require a anumu CC0pi event + if (!isCC0pi(event, -14, EnuMin, EnuMax)) + return false; + + TLorentzVector pnu = event->GetHMISParticle(-14)->fP; + TLorentzVector pmu = event->GetHMFSParticle(-13)->fP; + double Pmu = pmu.Vect().Mag(); + double CosThetaMu = cos(pnu.Vect().Angle(pmu.Vect())); + // Muon phase space + if (Pmu < 400 || Pmu > 3410) return false; + if (Pmu < 530 && Pmu>=400 && CosThetaMu<0.84) return false; + if (Pmu < 670 && Pmu>=530 && CosThetaMu<0.85) return false; + if (Pmu < 800 && Pmu>=670 && CosThetaMu<0.88) return false; + if (Pmu < 1000 &&Pmu>=800 && CosThetaMu<0.9) return false; + if (Pmu < 1380 && Pmu>=1000 && CosThetaMu<0.91) return false; + if (Pmu < 2010 && Pmu>=1380 && CosThetaMu<0.92) return false; + if (Pmu < 3410 && Pmu>=2010 && CosThetaMu<0.95) return false; + return true; } bool isT2K_CC0pi_STV(FitEvent *event, double EnuMin, double EnuMax) { // Require a numu CC0pi event if (!isCC0pi(event, 14, EnuMin, EnuMax)) return false; // Require at least one FS proton if (event->NumFSParticle(2212) == 0) return false; TLorentzVector pnu = event->GetHMISParticle(14)->fP; TLorentzVector pmu = event->GetHMFSParticle(13)->fP; TLorentzVector pp = event->GetHMFSParticle(2212)->fP; // Muon phase space // Pmu > 250 MeV, cos(theta_mu) > -0.6 (Sweet phase space!) if ((pmu.Vect().Mag() < 250) || cos(pnu.Vect().Angle(pmu.Vect())) < -0.6) { return false; } // Proton phase space // Pprot > 450 MeV, cos(theta_proton) > 0.4 if ((pp.Vect().Mag() < 450) || (pp.Vect().Mag() > 1E3) || (cos(pnu.Vect().Angle(pp.Vect())) < 0.4)) { return false; } return true; } bool isT2K_CC0pi_ifk(FitEvent *event, double EnuMin, double EnuMax) { // Require a numu CC0pi event if (!isCC0pi(event, 14, EnuMin, EnuMax)) return false; // Require at least one FS proton if (event->NumFSParticle(2212) == 0) return false; TLorentzVector pnu = event->GetHMISParticle(14)->fP; TLorentzVector pmu = event->GetHMFSParticle(13)->fP; TLorentzVector pp = event->GetHMFSParticle(2212)->fP; // Proton phase space // Pprot > 450 MeV, cos(theta_proton) > 0.4 if ((pp.Vect().Mag() < 450) || (cos(pnu.Vect().Angle(pp.Vect())) < 0.4)) { return false; } return true; } bool isT2K_CC0pi_1bin(FitEvent *event, double EnuMin, double EnuMax) { // Require a numu CC0pi event if (!isCC0pi(event, 14, EnuMin, EnuMax)) return false; // Require at least one FS proton if (event->NumFSParticle(2212) == 0) return false; TLorentzVector pnu = event->GetHMISParticle(14)->fP; TLorentzVector pmu = event->GetHMFSParticle(13)->fP; TLorentzVector pp = event->GetHMFSParticle(2212)->fP; // Muon phase space // if ((pmu.Vect().Mag() < 250) || cos(pnu.Vect().Angle(pmu.Vect())) < -0.6) { // return false; //} // Proton phase space if ((pp.Vect().Mag() < 450) || (cos(pnu.Vect().Angle(pp.Vect())) < 0.4)) { return false; } return true; } } // namespace SignalDef diff --git a/src/T2K/T2K_SignalDef.h b/src/T2K/T2K_SignalDef.h index cbc90b2..0d0f965 100644 --- a/src/T2K/T2K_SignalDef.h +++ b/src/T2K/T2K_SignalDef.h @@ -1,56 +1,56 @@ // Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret /******************************************************************************* * This file is part of NUISANCE. * * NUISANCE is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * NUISANCE is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with NUISANCE. If not, see . *******************************************************************************/ #ifndef T2K_SIGNALDEF_H_SEEN #define T2K_SIGNALDEF_H_SEEN #include "SignalDef.h" namespace SignalDef { bool isCC1pip_T2K_PRD97_012001(FitEvent *event, double EnuMin, double EnuMax); enum arxiv1909_03936_PScuts { kMuonFwd = (1 << 0), // cos(th_mu) > 0 kMuonHighEff = (1 << 1), // cos(th_mu) > 0.2, pmu > 200 kPionFwd = (1 << 2), // cos(th_pi) > 0 kPionVFwd = (1 << 3), // cos(th_pi) > 0.2 kPionHighMom = (1 << 4) // ppi > 200 }; bool isCC1pip_T2K_arxiv1909_03936(FitEvent *event, double EnuMin, double EnuMax, int cuts); enum PRD93112012_Ana { kAnalysis_I, kAnalysis_II, }; bool isT2K_CC0pi(FitEvent *event, double EnuMin, double EnuMax, int analysis); - + bool isT2K_CC0piNp(FitEvent *event, double EnuMin, double EnuMax); bool isT2K_CC0pi1p(FitEvent *event, double EnuMin, double EnuMax); bool isT2K_CC0pi0p(FitEvent *event, double EnuMin, double EnuMax); bool isT2K_CC0pi_STV(FitEvent *event, double EnuMin, double EnuMax); bool isT2K_CC0pi_1bin(FitEvent *event, double EnuMin, double EnuMax); bool isT2K_CC0pi_ifk(FitEvent *event, double EnuMin, double EnuMax); bool isT2K_CC0piAnuP0D(FitEvent *event, double EnuMin, double EnuMax); // TN328 } // namespace SignalDef #endif