diff --git a/src/FCN/SampleList.cxx b/src/FCN/SampleList.cxx index ff323b7..e0f3fcf 100644 --- a/src/FCN/SampleList.cxx +++ b/src/FCN/SampleList.cxx @@ -1,1440 +1,1431 @@ #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_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_1Dthetapi_nu.h" +#include "ArgoNeuT_CC1Pi_XSec_1Dthetamupi_antinu.h" #include "ArgoNeuT_CC1Pi_XSec_1Dthetamupi_nu.h" -#include "ArgoNeuT_CC1Pi_XSec_1Dpmu_antinu.h" -#include "ArgoNeuT_CC1Pi_XSec_1Dthetamu_antinu.h" #include "ArgoNeuT_CC1Pi_XSec_1Dthetapi_antinu.h" -#include "ArgoNeuT_CC1Pi_XSec_1Dthetamupi_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_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_1DQ2_antinu.h" #include "MiniBooNE_CCQE_XSec_1DQ2_nu.h" #include "MiniBooNE_CCQE_XSec_2DTcos_antinu.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_2D_nu.h" #include "MINERvA_CC0pi_XSec_1D_2018_nu.h" +#include "MINERvA_CC0pi_XSec_2D_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_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_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.h" +#include "T2K_CC0pi_XSec_2DPcos_nu_I.h" +#include "T2K_CC0pi_XSec_2DPcos_nu_II.h" // T2K CC-inclusive with full acceptance 2018 #include "T2K_CCinc_XSec_2DPcos_nu_nonuniform.h" // T2K STV CC0pi 2018 -#include "T2K_CC0pi_XSec_2DPcos_nu_nonuniform.h" -#include "T2K_CC0pinp_STV_XSec_1Ddpt_nu.h" -#include "T2K_CC0pinp_STV_XSec_1Ddphit_nu.h" -#include "T2K_CC0pinp_STV_XSec_1Ddat_nu.h" #include "T2K_CC0pi1p_XSec_3DPcoscos_nu_nonuniform.h" -#include "T2K_CC0pinp_ifk_XSec_3Dinfp_nu.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_2Dpmucosmu_nu.h" -#include "T2K_CC1pip_CH_XSec_1Dppi_nu.h" -#include "T2K_CC1pip_CH_XSec_1Dthpi_nu.h" -#include "T2K_CC1pip_CH_XSec_1Dthmupi_nu.h" -#include "T2K_CC1pip_CH_XSec_1DQ2_nu.h" #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" //#include "T2K_CC1pip_CH_XSec_1Dthq3pi_nu.h" //#include "T2K_CC1pip_CH_XSec_1DWrec_nu.h" //#include "T2K_CC1pip_CH_XSec_1Dq3_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" - #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_MuPiNoVA_1Dthetamu_nu.h" -#include "SciBooNE_CCCOH_MuPiNoVA_1Dpmu_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" #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 "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."); + << " shared object libraries."); } -DynamicSampleFactory* DynamicSampleFactory::glblDSF = NULL; +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) { +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"); + 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"); + 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 *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..."); + << ent->d_name + << " checking for relevant methods..."); - void* dlobj = + void *dlobj = dlopen((dirpath + ent->d_name).c_str(), RTLD_NOW | RTLD_GLOBAL); - char const* dlerr_cstr = dlerror(); + 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); + << (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); + << (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); + << (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); + << (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."); + << 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); + 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); + 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."); + << 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() { +DynamicSampleFactory &DynamicSampleFactory::Get() { if (!glblDSF) { glblDSF = new DynamicSampleFactory(); } return *glblDSF; } void DynamicSampleFactory::Print() { - std::map > ManifestSamples; + std::map> ManifestSamples; - for (std::map >::iterator smp_it = + 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 = + 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) { +bool DynamicSampleFactory::HasSample(std::string const &name) { return Samples.count(name); } -bool DynamicSampleFactory::HasSample(nuiskey& samplekey) { +bool DynamicSampleFactory::HasSample(nuiskey &samplekey) { return HasSample(samplekey.GetS("name")); } -MeasurementBase* DynamicSampleFactory::CreateSample(nuiskey& samplekey) { +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); + 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, +MeasurementBase *CreateSample(std::string name, std::string file, std::string type, std::string fkdt, - FitWeight* rw) { + 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) { +MeasurementBase *CreateSample(nuiskey samplekey) { #ifdef __USE_DYNSAMPLES__ if (DynamicSampleFactory::Get().HasSample(samplekey)) { NUIS_LOG(SAM, "Instantiating dynamic sample..."); - MeasurementBase* ds = DynamicSampleFactory::Get().CreateSample(samplekey); + 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(); + 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 -*/ + /* + 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")) { + 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_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")) { + /* + 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 - */ + /* + 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_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")) { + } else if (!name.compare("MINERvA_CC0pi_XSec_2Dptpz_nu")) { return (new MINERvA_CC0pi_XSec_2D_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")) { + } 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")) { + } 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") || - !name.compare("T2K_CC0pi_XSec_2DPcos_nu_I") || - !name.compare("T2K_CC0pi_XSec_2DPcos_nu_II")) { - return (new T2K_CC0pi_XSec_2DPcos_nu(samplekey)); + 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_nonuniform")) { - return (new T2K_CC0pi_XSec_2DPcos_nu_nonuniform(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)); /* 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_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_1Dppi_nu")) { - return (new T2K_CC1pip_CH_XSec_1Dppi_nu(samplekey)); - - } else if (!name.compare("T2K_CC1pip_CH_XSec_1Dthpi_nu")) { + } 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")) { + } 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")) { + } 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")) { + } 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")) { + } else if (!name.compare("T2K_CC1pip_CH_XSec_1DCosThAdler_nu")) { return (new T2K_CC1pip_CH_XSec_1DCosThAdler_nu(samplekey)); // Maybe something for the future: were in Raquel's thesis //} else if (!name.compare("T2K_CC1pip_CH_XSec_1Dq3_nu")) { - //return (new T2K_CC1pip_CH_XSec_1Dq3_nu(file, rw, type, fkdt)); + // return (new T2K_CC1pip_CH_XSec_1Dq3_nu(file, rw, type, fkdt)); //} else if (!name.compare("T2K_CC1pip_CH_XSec_1Dthq3pi_nu")) { - //return (new T2K_CC1pip_CH_XSec_1Dthq3pi_nu(file, rw, type, fkdt)); + // return (new T2K_CC1pip_CH_XSec_1Dthq3pi_nu(file, rw, type, fkdt)); //} else if (!name.compare("T2K_CC1pip_CH_XSec_1DWrec_nu")) { - //return (new T2K_CC1pip_CH_XSec_1DWrec_nu(file, rw, type, fkdt)); + // return (new T2K_CC1pip_CH_XSec_1DWrec_nu(file, rw, type, fkdt)); /* 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_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)); /* 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.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 diff --git a/src/FCN/sample_list.xml b/src/FCN/sample_list.xml index d8301e1..d6d9442 100644 --- a/src/FCN/sample_list.xml +++ b/src/FCN/sample_list.xml @@ -1,415 +1,415 @@ - + \ No newline at end of file diff --git a/src/T2K/CMakeLists.txt b/src/T2K/CMakeLists.txt index 339a1f3..e410de2 100644 --- a/src/T2K/CMakeLists.txt +++ b/src/T2K/CMakeLists.txt @@ -1,102 +1,102 @@ # 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 . ################################################################################ set(IMPLFILES -T2K_CC0pi_XSec_2DPcos_nu.cxx -T2K_CC0pi_XSec_2DPcos_nu_nonuniform.cxx +T2K_CC0pi_XSec_2DPcos_nu_I.cxx +T2K_CC0pi_XSec_2DPcos_nu_II.cxx T2K_CCinc_XSec_2DPcos_nu_nonuniform.cxx T2K_CC1pip_CH_XSec_2Dpmucosmu_nu.cxx T2K_CC1pip_CH_XSec_1Dppi_nu.cxx T2K_CC1pip_CH_XSec_1Dthpi_nu.cxx T2K_CC1pip_CH_XSec_1Dthmupi_nu.cxx T2K_CC1pip_CH_XSec_1DQ2_nu.cxx T2K_CC1pip_CH_XSec_1DAdlerPhi_nu.cxx T2K_CC1pip_CH_XSec_1DCosThAdler_nu.cxx T2K_CC1pip_H2O_XSec_1DEnuDelta_nu.cxx T2K_CC1pip_H2O_XSec_1DEnuMB_nu.cxx T2K_CC1pip_H2O_XSec_1Dcosmu_nu.cxx T2K_CC1pip_H2O_XSec_1Dcosmupi_nu.cxx T2K_CC1pip_H2O_XSec_1Dcospi_nu.cxx T2K_CC1pip_H2O_XSec_1Dpmu_nu.cxx T2K_CC1pip_H2O_XSec_1Dppi_nu.cxx T2K_CC0pinp_STV_XSec_1Ddpt_nu.cxx T2K_CC0pinp_STV_XSec_1Ddphit_nu.cxx T2K_CC0pinp_STV_XSec_1Ddat_nu.cxx T2K_CC0pi1p_XSec_3DPcoscos_nu_nonuniform.cxx T2K_CC0pinp_ifk_XSec_3Dinfp_nu.cxx T2K_CC0pinp_ifk_XSec_3Dinfa_nu.cxx T2K_CC0pinp_ifk_XSec_3Dinfip_nu.cxx T2K_SignalDef.cxx ) set(HEADERFILES -T2K_CC0pi_XSec_2DPcos_nu.h -T2K_CC0pi_XSec_2DPcos_nu_nonuniform.cxx -T2K_CCinc_XSec_2DPcos_nu_nonuniform.cxx +T2K_CC0pi_XSec_2DPcos_nu_I.h +T2K_CC0pi_XSec_2DPcos_nu_II.h +T2K_CCinc_XSec_2DPcos_nu_nonuniform.h T2K_CC1pip_CH_XSec_2Dpmucosmu_nu.h T2K_CC1pip_CH_XSec_1Dppi_nu.h T2K_CC1pip_CH_XSec_1Dthpi_nu.h T2K_CC1pip_CH_XSec_1Dthmupi_nu.h T2K_CC1pip_CH_XSec_1DQ2_nu.h T2K_CC1pip_CH_XSec_1DAdlerPhi_nu.h T2K_CC1pip_CH_XSec_1DCosThAdler_nu.h T2K_CC1pip_H2O_XSec_1DEnuDelta_nu.h T2K_CC1pip_H2O_XSec_1DEnuMB_nu.h T2K_CC1pip_H2O_XSec_1Dcosmu_nu.h T2K_CC1pip_H2O_XSec_1Dcosmupi_nu.h T2K_CC1pip_H2O_XSec_1Dcospi_nu.h T2K_CC1pip_H2O_XSec_1Dpmu_nu.h T2K_CC1pip_H2O_XSec_1Dppi_nu.h T2K_CC0pinp_STV_XSec_1Ddpt_nu.h T2K_CC0pinp_STV_XSec_1Ddphit_nu.h T2K_CC0pinp_STV_XSec_1Ddat_nu.h T2K_CC0pi1p_XSec_3DPcoscos_nu_nonuniform.h T2K_CC0pinp_ifk_XSec_3Dinfp_nu.h T2K_CC0pinp_ifk_XSec_3Dinfa_nu.h T2K_CC0pinp_ifk_XSec_3Dinfip_nu.h T2K_SignalDef.h ) set(LIBNAME expT2K) if(CMAKE_BUILD_TYPE MATCHES DEBUG) add_library(${LIBNAME} STATIC ${IMPLFILES}) else(CMAKE_BUILD_TYPE MATCHES RELEASE) add_library(${LIBNAME} SHARED ${IMPLFILES}) endif() include_directories(${MINIMUM_INCLUDE_DIRECTORIES}) set_target_properties(${LIBNAME} PROPERTIES VERSION "${NUISANCE_VERSION_MAJOR}.${NUISANCE_VERSION_MINOR}.${NUISANCE_VERSION_REVISION}") #set_target_properties(${LIBNAME} PROPERTIES LINK_FLAGS ${ROOT_LD_FLAGS}) if(DEFINED PROJECTWIDE_EXTRA_DEPENDENCIES) add_dependencies(${LIBNAME} ${PROJECTWIDE_EXTRA_DEPENDENCIES}) endif() install(TARGETS ${LIBNAME} DESTINATION lib) #Can uncomment this to install the headers... but is it really neccessary? #install(FILES ${HEADERFILES} DESTINATION include) set(MODULETargets ${MODULETargets} ${LIBNAME} PARENT_SCOPE) diff --git a/src/T2K/T2K_CC0pi1p_XSec_3DPcoscos_nu_nonuniform.cxx b/src/T2K/T2K_CC0pi1p_XSec_3DPcoscos_nu_nonuniform.cxx index 94277b2..82ba031 100644 --- a/src/T2K/T2K_CC0pi1p_XSec_3DPcoscos_nu_nonuniform.cxx +++ b/src/T2K/T2K_CC0pi1p_XSec_3DPcoscos_nu_nonuniform.cxx @@ -1,226 +1,252 @@ // 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 . -*******************************************************************************/ + * 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) { //******************************************************************** +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"; + 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.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())); + // 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){ +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){ +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; + TLorentzVector Pp = event->GetHMFSParticle(2212)->fP; - //double pmu = Pmu.Vect().Mag()/1000.; - double pp = Pp.Vect().Mag()/1000.; + // 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(){ +void T2K_CC0pi1p_XSec_3DPcoscos_nu_nonuniform::FillHistograms() { Measurement1D::FillHistograms(); - if (Signal){ - fMCHist_Fine2D->Fill( fXVar, fYVar, Weight ); - FillMCSlice( fXVar, fYVar, fZVar, Weight ); - + 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(){ +void T2K_CC0pi1p_XSec_3DPcoscos_nu_nonuniform::ConvertEventRates() { - for (int i = 0; i < 4; i++){ + 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); - + // 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)); + 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::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(){ +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(); - + 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"); + 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); + 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"); + + 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++){ + 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; + (*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"); + 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 + 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)); + // 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); + // 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) ); + // 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))); + 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(fDataHist_Slices[i], kCMD_Write); SetAutoProcessTH1(fMCHist_Slices[i]); fMCHist_Slices[i]->ClearBinContents(); } return; }; diff --git a/src/T2K/T2K_CC0pi_XSec_2DPcos_nu.cxx b/src/T2K/T2K_CC0pi_XSec_2DPcos_nu.cxx deleted file mode 100644 index 4196e6f..0000000 --- a/src/T2K/T2K_CC0pi_XSec_2DPcos_nu.cxx +++ /dev/null @@ -1,170 +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_CC0pi_XSec_2DPcos_nu.h" - - - -//******************************************************************** -T2K_CC0pi_XSec_2DPcos_nu::T2K_CC0pi_XSec_2DPcos_nu(nuiskey samplekey) { -//******************************************************************** - - // Sample overview --------------------------------------------------- - std::string descrip = "T2K_CC0pi_XSec_2DPcos_nu sample. \n" \ - "Target: CH \n" \ - "Flux: MINERvA Medium Energy FHC numu \n" \ - "Signal: CC-inclusive with theta < 20deg \n"; - - // Setup common settings - fSettings = LoadSampleSettings(samplekey); - fSettings.SetDescription(descrip); - fSettings.SetXTitle("P_{#mu} (GeV)"); - fSettings.SetYTitle("cos#theta_{#mu}"); - fSettings.SetZTitle("d^{2}#sigma/dP_{#mu}dcos#theta_{#mu} (cm^{2}/GeV)"); - fSettings.SetAllowedTypes("DIAG,FULL/FREE,SHAPE,FIX/SYSTCOV/STATCOV","FIX"); - fSettings.SetEnuRange(0.0, 10.0); - fSettings.DefineAllowedTargets("C,H"); - - if (fName == "T2K_CC0pi_XSec_2DPcos_nu_I") fAnalysis = 1; - else fAnalysis = 2; - - - // CCQELike plot information - fSettings.SetTitle("T2K_CC0pi_XSec_2DPcos_nu"); - fSettings.DefineAllowedSpecies("numu"); - - forwardgoing = (fSettings.GetS("type").find("REST") != std::string::npos); - - FinaliseSampleSettings(); - - // Scaling Setup --------------------------------------------------- - // ScaleFactor automatically setup for DiffXSec/cm2/Nucleon - fScaleFactor = ((GetEventHistogram()->Integral("width")/(fNEvents+0.)) * 1E-38 / (TotalIntegratedFlux())); - - // Setup Histograms - SetHistograms(); - - // Final setup --------------------------------------------------- - FinaliseMeasurement(); - -}; - - -bool T2K_CC0pi_XSec_2DPcos_nu::isSignal(FitEvent *event){ - return SignalDef::isT2K_CC0pi(event, EnuMin, EnuMax, forwardgoing); -}; - -void T2K_CC0pi_XSec_2DPcos_nu::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())); - - fXVar = pmu; - fYVar = CosThetaMu; - - return; -}; - -void T2K_CC0pi_XSec_2DPcos_nu::SetHistograms(){ - - fIsSystCov = fSettings.GetS("type").find("SYSTCOV") != std::string::npos; - fIsStatCov = fSettings.GetS("type").find("STATCOV") != std::string::npos; - fIsNormCov = fSettings.GetS("type").find("NORMCOV") != std::string::npos; - fNDataPointsX = 12; - fNDataPointsY = 10; - - // Open file - std::string infile = FitPar::GetDataBase()+"/T2K/CC0pi/T2K_CC0PI_2DPmuCosmu_Data.root"; - TFile* rootfile = new TFile(infile.c_str(), "READ"); - TH2D* tempcov = NULL; - - // ANALYSIS 2 - if (fAnalysis == 2){ - - // Get Data - fDataHist = (TH2D*) rootfile->Get("analysis2_data"); - fDataHist->SetDirectory(0); - fDataHist->SetNameTitle((fName + "_data").c_str(), - (fName + "_data" + fPlotTitles).c_str()); - // For some reason the error on the data in the data release is 1E-20 - // That is wrong, so set it to zero here - for (int i = 0; i < fDataHist->GetXaxis()->GetNbins()+1; ++i) { - for (int j = 0; j < fDataHist->GetYaxis()->GetNbins()+1; ++j) { - fDataHist->SetBinError(i+1, j+1, 0.0); - } - } - - // Get Map - fMapHist = (TH2I*) rootfile->Get("analysis2_map"); - fMapHist->SetDirectory(0); - fMapHist->SetNameTitle((fName + "_map").c_str(), - (fName + "_map" + fPlotTitles).c_str()); - - // Get Syst/Stat Covar - TH2D* tempsyst = (TH2D*) rootfile->Get("analysis2_systcov"); - TH2D* tempstat = (TH2D*) rootfile->Get("analysis2_statcov"); - TH2D* tempnorm = (TH2D*) rootfile->Get("analysis2_normcov"); - - // Create covar [Default is both] - tempcov = (TH2D*) tempsyst->Clone(); - tempcov->Reset(); - - if (fIsSystCov) tempcov->Add(tempsyst); - if (fIsStatCov) tempcov->Add(tempstat); - if (fIsNormCov) tempcov->Add(tempnorm); - - if (!fIsSystCov && !fIsStatCov && !fIsNormCov){ - tempcov->Add(tempsyst); - tempcov->Add(tempstat); - tempcov->Add(tempnorm); - } - } - - - if (!tempcov){ - NUIS_ABORT("TEMPCOV NOT SET"); - } - - // Setup Covar - int nbins = tempcov->GetNbinsX(); - fFullCovar = new TMatrixDSym(nbins); - - for (int i = 0; i < nbins; i++){ - for (int j = 0; j < nbins; j++){ - (*fFullCovar)(i,j) = tempcov->GetBinContent(i+1,j+1); - } - } - covar = StatUtils::GetInvert(fFullCovar); - fDecomp = StatUtils::GetDecomp(covar); - - // Set Data Errors - StatUtils::SetDataErrorFromCov(fDataHist, fFullCovar, fMapHist, 1E-38); - - // Remove root file - rootfile->Close(); - return; -}; diff --git a/src/T2K/T2K_CC0pi_XSec_2DPcos_nu.h b/src/T2K/T2K_CC0pi_XSec_2DPcos_nu.h deleted file mode 100644 index afcbb61..0000000 --- a/src/T2K/T2K_CC0pi_XSec_2DPcos_nu.h +++ /dev/null @@ -1,67 +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 . -*******************************************************************************/ -#ifndef T2K_CC0PI_2DPCOS_NU_H_SEEN -#define T2K_CC0PI_2DPCOS_NU_H_SEEN - -#include "Measurement2D.h" - -class T2K_CC0pi_XSec_2DPcos_nu : public Measurement2D { -public: - - /// Basic Constructor. - /// /brief Parses two different measurements. - /// - /// T2K_CC0pi_XSec_2DPcos_nu -> T2K CC0PI Analysis 2 - /// T2K_CC0pi_XSec_2DPcos_nu_I -> T2K CC0PI Analysis 1 - /// T2K_CC0pi_XSec_2DPcos_nu_II -> T2K CC0PI Analysis 2 - T2K_CC0pi_XSec_2DPcos_nu(nuiskey samplekey); - - /// Virtual Destructor - ~T2K_CC0pi_XSec_2DPcos_nu() {}; - - /// Numu CC0PI Signal Definition - /// - /// /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); - - /// Have to do a weird event scaling for analysis 1 - // void ConvertEventRates(); - - private: - - bool forwardgoing; - bool only_allowed_particles; - bool numu_event; - double numu_energy; - int particle_pdg; - double pmu, CosThetaMu; - int fAnalysis; - - bool fIsSystCov, fIsStatCov, fIsNormCov; - -}; - -#endif diff --git a/src/T2K/T2K_CC0pi_XSec_2DPcos_nu_I.cxx b/src/T2K/T2K_CC0pi_XSec_2DPcos_nu_I.cxx new file mode 100644 index 0000000..9f3e0c8 --- /dev/null +++ b/src/T2K/T2K_CC0pi_XSec_2DPcos_nu_I.cxx @@ -0,0 +1,229 @@ +// 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_CC0pi_XSec_2DPcos_nu_I.h" + +//******************************************************************** +T2K_CC0pi_XSec_2DPcos_nu_I::T2K_CC0pi_XSec_2DPcos_nu_I(nuiskey samplekey) { + //******************************************************************** + + // Sample overview --------------------------------------------------- + std::string descrip = "T2K_CC0pi_XSec_2DPcos_nu_I sample. \n" + "Target: CH \n" + "Flux: T2K 2.5 degree off-axis (ND280) \n" + "Signal: CC0pi\n" + "https://journals.aps.org/prd/abstract/10.1103/" + "PhysRevD.93.112012 Analysis I"; + + // Setup common settings + fSettings = LoadSampleSettings(samplekey); + fSettings.SetDescription(descrip); + fSettings.SetXTitle("P_{#mu} (GeV)"); + fSettings.SetYTitle("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"); + + // CCQELike plot information + fSettings.SetTitle("T2K_CC0pi_XSec_2DPcos_nu_I"); + fSettings.DefineAllowedSpecies("numu"); + + FinaliseSampleSettings(); + + // Scaling Setup --------------------------------------------------- + // ScaleFactor automatically setup for DiffXSec/cm2/Nucleon + fScaleFactor = ((GetEventHistogram()->Integral("width") / (fNEvents + 0.)) * + 1E-38 / (TotalIntegratedFlux())); + + // Setup Histograms + SetHistograms(); + + // Final setup --------------------------------------------------- + FinaliseMeasurement(); +}; + +bool T2K_CC0pi_XSec_2DPcos_nu_I::isSignal(FitEvent *event) { + return SignalDef::isT2K_CC0pi(event, EnuMin, EnuMax, SignalDef::kAnalysis_I); +}; + +void T2K_CC0pi_XSec_2DPcos_nu_I::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())); + + fXVar = pmu; + fYVar = CosThetaMu; + + return; +}; + +void T2K_CC0pi_XSec_2DPcos_nu_I::FillHistograms() { + + Measurement1D::FillHistograms(); + if (Signal) { + fMCHist_Fine2D->Fill(fXVar, fYVar, Weight); + FillMCSlice(fXVar, fYVar, 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_CC0pi_XSec_2DPcos_nu_I::ConvertEventRates() { + + for (int i = 0; i < 9; i++) { + fMCHist_Slices[i]->GetSumw2(); + } + + // Do standard conversion. + Measurement1D::ConvertEventRates(); + + // First scale MC slices also by their width in Y + fMCHist_Slices[0]->Scale(1.0 / 1.00); + fMCHist_Slices[1]->Scale(1.0 / 0.60); + fMCHist_Slices[2]->Scale(1.0 / 0.10); + fMCHist_Slices[3]->Scale(1.0 / 0.10); + fMCHist_Slices[4]->Scale(1.0 / 0.05); + fMCHist_Slices[5]->Scale(1.0 / 0.05); + fMCHist_Slices[6]->Scale(1.0 / 0.04); + fMCHist_Slices[7]->Scale(1.0 / 0.04); + fMCHist_Slices[8]->Scale(1.0 / 0.02); + + // Now Convert into 1D list + fMCHist->Reset(); + int bincount = 0; + for (int i = 0; i < 9; i++) { + for (int j = 0; j < fDataHist_Slices[i]->GetNbinsX(); 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_CC0pi_XSec_2DPcos_nu_I::FillMCSlice(double x, double y, double w) { + + if (y >= -1.0 and y < 0.0) + fMCHist_Slices[0]->Fill(x, w); + else if (y >= 0.0 and y < 0.6) + fMCHist_Slices[1]->Fill(x, w); + else if (y >= 0.6 and y < 0.7) + fMCHist_Slices[2]->Fill(x, w); + else if (y >= 0.7 and y < 0.8) + fMCHist_Slices[3]->Fill(x, w); + else if (y >= 0.8 and y < 0.85) + fMCHist_Slices[4]->Fill(x, w); + else if (y >= 0.85 and y < 0.90) + fMCHist_Slices[5]->Fill(x, w); + else if (y >= 0.90 and y < 0.94) + fMCHist_Slices[6]->Fill(x, w); + else if (y >= 0.94 and y < 0.98) + fMCHist_Slices[7]->Fill(x, w); + else if (y >= 0.98 and y <= 1.00) + fMCHist_Slices[8]->Fill(x, w); +} + +void T2K_CC0pi_XSec_2DPcos_nu_I::SetHistograms() { + + // Read in 1D Data Histograms + fInputFile = new TFile( + (FitPar::GetDataBase() + "/T2K/CC0pi/T2K_CC0PI_2DPmuCosmu_Data.root") + .c_str(), + "READ"); + // fInputFile->ls(); + + // Read in 1D Data + fDataHist = (TH1D *)fInputFile->Get("datahist"); + + fMCHist_Fine2D = new TH2D("T2K_CC0pi_XSec_2DPcos_nu_I_Fine2D", + "T2K_CC0pi_XSec_2DPcos_nu_I_Fine2D", 400, 0.0, 30.0, + 100, -1.0, 1.0); + SetAutoProcessTH1(fMCHist_Fine2D); + + TH2D *tempcov = (TH2D *)fInputFile->Get("analysis1_totcov"); + + 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); + } + } + covar = StatUtils::GetInvert(fFullCovar); + fDecomp = StatUtils::GetDecomp(fFullCovar); + + // Read in 2D Data + fDataPoly = (TH2Poly *)fInputFile->Get("datapoly"); + fDataPoly->SetNameTitle("T2K_CC0pi_XSec_2DPcos_nu_I_datapoly", + "T2K_CC0pi_XSec_2DPcos_nu_I_datapoly"); + SetAutoProcessTH1(fDataPoly, kCMD_Write); + fDataHist->Reset(); + + // Read in 2D Data Slices and Make MC Slices + int bincount = 0; + for (int i = 0; i < 9; i++) { + + // Get Data Histogram + // fInputFile->ls(); + fDataHist_Slices.push_back( + (TH1D *)fInputFile->Get(Form("dataslice_%i", i))->Clone()); + fDataHist_Slices[i]->SetNameTitle( + Form("T2K_CC0pi_XSec_2DPcos_nu_I_data_Slice%i", i), + (Form("T2K_CC0pi_XSec_2DPcos_nu_I_data_Slice%i", i))); + + // Loop over nbins and set errors from covar + for (int j = 0; j < fDataHist_Slices[i]->GetNbinsX(); j++) { + fDataHist_Slices[i]->SetBinError( + j + 1, sqrt((*fFullCovar)(bincount, bincount)) * 1E-38); + + // std::cout << "Setting data hist " << + // fDataHist_Slices[i]->GetBinContent(j+1) << " " << + // fDataHist_Slices[i]->GetBinError(j+1) << std::endl; + fDataHist->SetBinContent(bincount + 1, + fDataHist_Slices[i]->GetBinContent(j + 1)); + fDataHist->SetBinError(bincount + 1, + fDataHist_Slices[i]->GetBinError(j + 1)); + + bincount++; + } + + // Make MC Clones + fMCHist_Slices.push_back((TH1D *)fDataHist_Slices[i]->Clone()); + fMCHist_Slices[i]->SetNameTitle( + Form("T2K_CC0pi_XSec_2DPcos_nu_I_MC_Slice%i", i), + (Form("T2K_CC0pi_XSec_2DPcos_nu_I_MC_Slice%i", i))); + + SetAutoProcessTH1(fDataHist_Slices[i], kCMD_Write); + SetAutoProcessTH1(fMCHist_Slices[i]); + // fMCHist_Slices[i]->Reset(); + } + + return; +}; diff --git a/src/T2K/T2K_CC0pi_XSec_2DPcos_nu_nonuniform.h b/src/T2K/T2K_CC0pi_XSec_2DPcos_nu_I.h similarity index 76% rename from src/T2K/T2K_CC0pi_XSec_2DPcos_nu_nonuniform.h rename to src/T2K/T2K_CC0pi_XSec_2DPcos_nu_I.h index fa0b8c9..269a6de 100644 --- a/src/T2K/T2K_CC0pi_XSec_2DPcos_nu_nonuniform.h +++ b/src/T2K/T2K_CC0pi_XSec_2DPcos_nu_I.h @@ -1,86 +1,72 @@ // 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_CC0PI_2DPCOS_NU_NONUNIFORM_H_SEEN #define T2K_CC0PI_2DPCOS_NU_NONUNIFORM_H_SEEN #include "Measurement1D.h" #include "TH2Poly.h" #include "MeasurementVariableBox2D.h" -class T2K_CC0pi_XSec_2DPcos_nu_nonuniform : public Measurement1D { +class T2K_CC0pi_XSec_2DPcos_nu_I : public Measurement1D { public: - /// Basic Constructor. - /// /brief Parses two different measurements. - /// - /// T2K_CC0pi_XSec_2DPcos_nu_nonuniform -> T2K CC0PI Analysis 2 - /// T2K_CC0pi_XSec_2DPcos_nu_nonuniform_I -> T2K CC0PI Analysis 1 - /// T2K_CC0pi_XSec_2DPcos_nu_nonuniform_II -> T2K CC0PI Analysis 2 - T2K_CC0pi_XSec_2DPcos_nu_nonuniform(nuiskey samplekey); + T2K_CC0pi_XSec_2DPcos_nu_I(nuiskey samplekey); /// Virtual Destructor - ~T2K_CC0pi_XSec_2DPcos_nu_nonuniform() {}; + ~T2K_CC0pi_XSec_2DPcos_nu_I() {}; /// Numu CC0PI Signal Definition /// /// /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(); /// \brief Create Q2 Box to save correction info inline MeasurementVariableBox* CreateBox(){ return new MeasurementVariableBox2D(); }; private: - bool forwardgoing; - bool only_allowed_particles; - bool numu_event; - double numu_energy; - int particle_pdg; - double pmu, CosThetaMu; - int fAnalysis; - bool fIsSystCov, fIsStatCov, fIsNormCov; TH2Poly* fDataPoly; TH2Poly* fMCPoly; TFile* fInputFile; TH2D* fMCHist_Fine2D; std::vector fMCHist_Slices; std::vector fDataHist_Slices; void FillMCSlice(double x, double y, double w); }; #endif diff --git a/src/T2K/T2K_CC0pi_XSec_2DPcos_nu_II.cxx b/src/T2K/T2K_CC0pi_XSec_2DPcos_nu_II.cxx new file mode 100644 index 0000000..3c71549 --- /dev/null +++ b/src/T2K/T2K_CC0pi_XSec_2DPcos_nu_II.cxx @@ -0,0 +1,168 @@ +// 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_CC0pi_XSec_2DPcos_nu_II.h" + +//******************************************************************** +T2K_CC0pi_XSec_2DPcos_nu_II::T2K_CC0pi_XSec_2DPcos_nu_II(nuiskey samplekey) { + //******************************************************************** + + // Sample overview --------------------------------------------------- + std::string descrip = "T2K_CC0pi_XSec_2DPcos_nu_II sample. \n" + "Target: CH \n" + "Flux: T2K 2.5 degree off-axis (ND280) \n" + "Signal: CC0pi with p_mu > 200 MeV\n" + " cth_mu > 0 \n" + "https://journals.aps.org/prd/abstract/10.1103/" + "PhysRevD.93.112012 Analysis II"; + + // Setup common settings + fSettings = LoadSampleSettings(samplekey); + fSettings.SetDescription(descrip); + fSettings.SetXTitle("P_{#mu} (GeV)"); + fSettings.SetYTitle("cos#theta_{#mu}"); + fSettings.SetZTitle("d^{2}#sigma/dP_{#mu}dcos#theta_{#mu} (cm^{2}/GeV)"); + fSettings.SetAllowedTypes("DIAG,FULL/FREE,SHAPE,FIX/SYSTCOV/STATCOV", "FIX"); + fSettings.SetEnuRange(0.0, 10.0); + fSettings.DefineAllowedTargets("C,H"); + + // CCQELike plot information + fSettings.SetTitle("T2K_CC0pi_XSec_2DPcos_nu_II"); + fSettings.DefineAllowedSpecies("numu"); + + FinaliseSampleSettings(); + + // Scaling Setup --------------------------------------------------- + // ScaleFactor automatically setup for DiffXSec/cm2/Nucleon + fScaleFactor = ((GetEventHistogram()->Integral("width") / (fNEvents + 0.)) * + 1E-38 / (TotalIntegratedFlux())); + + // Setup Histograms + SetHistograms(); + + // Final setup --------------------------------------------------- + FinaliseMeasurement(); +}; + +bool T2K_CC0pi_XSec_2DPcos_nu_II::isSignal(FitEvent *event) { + return SignalDef::isT2K_CC0pi(event, EnuMin, EnuMax, SignalDef::kAnalysis_II); +}; + +void T2K_CC0pi_XSec_2DPcos_nu_II::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())); + + fXVar = pmu; + fYVar = CosThetaMu; + + return; +}; + +void T2K_CC0pi_XSec_2DPcos_nu_II::SetHistograms() { + + fIsSystCov = fSettings.GetS("type").find("SYSTCOV") != std::string::npos; + fIsStatCov = fSettings.GetS("type").find("STATCOV") != std::string::npos; + fIsNormCov = fSettings.GetS("type").find("NORMCOV") != std::string::npos; + fNDataPointsX = 12; + fNDataPointsY = 10; + + // Open file + std::string infile = + FitPar::GetDataBase() + "/T2K/CC0pi/T2K_CC0PI_2DPmuCosmu_Data.root"; + TFile *rootfile = new TFile(infile.c_str(), "READ"); + TH2D *tempcov = NULL; + + // Get Data + fDataHist = (TH2D *)rootfile->Get("analysis2_data"); + fDataHist->SetDirectory(0); + fDataHist->SetNameTitle((fName + "_data").c_str(), + (fName + "_data" + fPlotTitles).c_str()); + // For some reason the error on the data in the data release is 1E-20 + // That is wrong, so set it to zero here + for (int i = 0; i < fDataHist->GetXaxis()->GetNbins() + 1; ++i) { + for (int j = 0; j < fDataHist->GetYaxis()->GetNbins() + 1; ++j) { + fDataHist->SetBinError(i + 1, j + 1, 0.0); + } + } + + // Get Map + fMapHist = (TH2I *)rootfile->Get("analysis2_map"); + fMapHist->SetDirectory(0); + fMapHist->SetNameTitle((fName + "_map").c_str(), + (fName + "_map" + fPlotTitles).c_str()); + + // Get Syst/Stat Covar + TH2D *tempsyst = 0; + rootfile->GetObject("analysis2_systcov", tempsyst); + TH2D *tempstat = 0; + rootfile->GetObject("analysis2_statcov", tempstat); + TH2D *tempnorm = 0; + rootfile->GetObject("analysis2_normcov", tempnorm); + + if (!tempsyst) { + NUIS_ABORT("tempsyst not found"); + } + + // Create covar [Default is both] + tempcov = (TH2D *)tempsyst->Clone(); + tempcov->Reset(); + + if (fIsSystCov) { + tempcov->Add(tempsyst); + } + if (fIsStatCov) { + tempcov->Add(tempstat); + } + if (fIsNormCov) { + tempcov->Add(tempnorm); + } + + if (!fIsSystCov && !fIsStatCov && !fIsNormCov) { + tempcov->Add(tempsyst); + tempcov->Add(tempstat); + tempcov->Add(tempnorm); + } + + // Setup Covar + int nbins = tempcov->GetNbinsX(); + fFullCovar = new TMatrixDSym(nbins); + + for (int i = 0; i < nbins; i++) { + for (int j = 0; j < nbins; j++) { + (*fFullCovar)(i, j) = tempcov->GetBinContent(i + 1, j + 1); + } + } + covar = StatUtils::GetInvert(fFullCovar); + fDecomp = StatUtils::GetDecomp(covar); + + // Set Data Errors + StatUtils::SetDataErrorFromCov(fDataHist, fFullCovar, fMapHist, 1E-38); + + // Remove root file + rootfile->Close(); +}; diff --git a/src/T2K/T2K_CC0pi_XSec_2DPcos_nu_II.h b/src/T2K/T2K_CC0pi_XSec_2DPcos_nu_II.h new file mode 100644 index 0000000..a0c4503 --- /dev/null +++ b/src/T2K/T2K_CC0pi_XSec_2DPcos_nu_II.h @@ -0,0 +1,49 @@ +// 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_CC0PI_2DPCOS_NU_H_SEEN +#define T2K_CC0PI_2DPCOS_NU_H_SEEN + +#include "Measurement2D.h" + +class T2K_CC0pi_XSec_2DPcos_nu_II : public Measurement2D { +public: + /// Basic Constructor. + T2K_CC0pi_XSec_2DPcos_nu_II(nuiskey samplekey); + + /// Virtual Destructor + ~T2K_CC0pi_XSec_2DPcos_nu_II(){}; + + /// Numu CC0PI Signal Definition + /// + /// /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); + +private: + bool fIsSystCov, fIsStatCov, fIsNormCov; +}; + +#endif diff --git a/src/T2K/T2K_CC0pi_XSec_2DPcos_nu_nonuniform.cxx b/src/T2K/T2K_CC0pi_XSec_2DPcos_nu_nonuniform.cxx deleted file mode 100644 index 014ed85..0000000 --- a/src/T2K/T2K_CC0pi_XSec_2DPcos_nu_nonuniform.cxx +++ /dev/null @@ -1,215 +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_CC0pi_XSec_2DPcos_nu_nonuniform.h" - - - -//******************************************************************** -T2K_CC0pi_XSec_2DPcos_nu_nonuniform::T2K_CC0pi_XSec_2DPcos_nu_nonuniform(nuiskey samplekey) { -//******************************************************************** - - // Sample overview --------------------------------------------------- - std::string descrip = "T2K_CC0pi_XSec_2DPcos_nu_nonuniform sample. \n" \ - "Target: CH \n" \ - "Flux: MINERvA Medium Energy FHC numu \n" \ - "Signal: CC-inclusive with theta < 20deg \n"; - - // Setup common settings - fSettings = LoadSampleSettings(samplekey); - fSettings.SetDescription(descrip); - fSettings.SetXTitle("P_{#mu} (GeV)"); - fSettings.SetYTitle("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_CC0pi_XSec_2DPcos_nu_nonuniform"); - fSettings.DefineAllowedSpecies("numu"); - - forwardgoing = (fSettings.GetS("type").find("REST") != std::string::npos); - - FinaliseSampleSettings(); - - // Scaling Setup --------------------------------------------------- - // ScaleFactor automatically setup for DiffXSec/cm2/Nucleon - fScaleFactor = ((GetEventHistogram()->Integral("width")/(fNEvents+0.)) * 1E-38 / (TotalIntegratedFlux())); - - // Setup Histograms - SetHistograms(); - - // Final setup --------------------------------------------------- - FinaliseMeasurement(); - -}; - - -bool T2K_CC0pi_XSec_2DPcos_nu_nonuniform::isSignal(FitEvent *event){ - return SignalDef::isT2K_CC0pi(event, EnuMin, EnuMax, forwardgoing); -}; - -void T2K_CC0pi_XSec_2DPcos_nu_nonuniform::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())); - - fXVar = pmu; - fYVar = CosThetaMu; - - return; -}; - -void T2K_CC0pi_XSec_2DPcos_nu_nonuniform::FillHistograms(){ - - Measurement1D::FillHistograms(); - if (Signal){ - fMCHist_Fine2D->Fill( fXVar, fYVar, Weight ); - FillMCSlice( fXVar, fYVar, 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_CC0pi_XSec_2DPcos_nu_nonuniform::ConvertEventRates(){ - - for (int i = 0; i < 9; i++){ - fMCHist_Slices[i]->GetSumw2(); - } - - // Do standard conversion. - Measurement1D::ConvertEventRates(); - - // First scale MC slices also by their width in Y - fMCHist_Slices[0]->Scale(1.0 / 1.00); - fMCHist_Slices[1]->Scale(1.0 / 0.60); - fMCHist_Slices[2]->Scale(1.0 / 0.10); - fMCHist_Slices[3]->Scale(1.0 / 0.10); - fMCHist_Slices[4]->Scale(1.0 / 0.05); - fMCHist_Slices[5]->Scale(1.0 / 0.05); - fMCHist_Slices[6]->Scale(1.0 / 0.04); - fMCHist_Slices[7]->Scale(1.0 / 0.04); - fMCHist_Slices[8]->Scale(1.0 / 0.02); - - - // Now Convert into 1D list - fMCHist->Reset(); - int bincount = 0; - for (int i = 0; i < 9; i++){ - for (int j = 0; j < fDataHist_Slices[i]->GetNbinsX(); 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_CC0pi_XSec_2DPcos_nu_nonuniform::FillMCSlice(double x, double y, double w){ - - if (y >= -1.0 and y < 0.0) fMCHist_Slices[0]->Fill(x,w); - else if (y >= 0.0 and y < 0.6) fMCHist_Slices[1]->Fill(x,w); - else if (y >= 0.6 and y < 0.7) fMCHist_Slices[2]->Fill(x,w); - else if (y >= 0.7 and y < 0.8) fMCHist_Slices[3]->Fill(x,w); - else if (y >= 0.8 and y < 0.85) fMCHist_Slices[4]->Fill(x,w); - else if (y >= 0.85 and y < 0.90) fMCHist_Slices[5]->Fill(x,w); - else if (y >= 0.90 and y < 0.94) fMCHist_Slices[6]->Fill(x,w); - else if (y >= 0.94 and y < 0.98) fMCHist_Slices[7]->Fill(x,w); - else if (y >= 0.98 and y <= 1.00) fMCHist_Slices[8]->Fill(x,w); -} - - -void T2K_CC0pi_XSec_2DPcos_nu_nonuniform::SetHistograms(){ - - // Read in 1D Data Histograms - fInputFile = new TFile( (FitPar::GetDataBase() + "/T2K/CC0pi/T2K_CC0PI_2DPmuCosmu_Data.root").c_str(),"READ"); - //fInputFile->ls(); - - // Read in 1D Data - fDataHist = (TH1D*) fInputFile->Get("datahist"); - - fMCHist_Fine2D = new TH2D("T2K_CC0pi_XSec_2DPcos_nu_nonuniform_Fine2D","T2K_CC0pi_XSec_2DPcos_nu_nonuniform_Fine2D", 400, 0.0,30.0,100,-1.0,1.0); - SetAutoProcessTH1(fMCHist_Fine2D); - - TH2D* tempcov = (TH2D*) fInputFile->Get("analysis1_totcov"); - - 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); - } - } - covar = StatUtils::GetInvert(fFullCovar); - fDecomp = StatUtils::GetDecomp(fFullCovar); - - // Read in 2D Data - fDataPoly = (TH2Poly*) fInputFile->Get("datapoly"); - fDataPoly->SetNameTitle("T2K_CC0pi_XSec_2DPcos_nu_nonuniform_datapoly","T2K_CC0pi_XSec_2DPcos_nu_nonuniform_datapoly"); - SetAutoProcessTH1(fDataPoly, kCMD_Write); - fDataHist->Reset(); - - // Read in 2D Data Slices and Make MC Slices - int bincount = 0; - for (int i = 0; i < 9; i++){ - - // Get Data Histogram - //fInputFile->ls(); - fDataHist_Slices.push_back((TH1D*)fInputFile->Get(Form("dataslice_%i",i))->Clone()); - fDataHist_Slices[i]->SetNameTitle(Form("T2K_CC0pi_XSec_2DPcos_nu_nonuniform_data_Slice%i",i), - (Form("T2K_CC0pi_XSec_2DPcos_nu_nonuniform_data_Slice%i",i))); - - // Loop over nbins and set errors from covar - for (int j = 0; j < fDataHist_Slices[i]->GetNbinsX(); j++){ - fDataHist_Slices[i]->SetBinError(j+1, sqrt((*fFullCovar)(bincount,bincount)) * 1E-38); - - //std::cout << "Setting data hist " << fDataHist_Slices[i]->GetBinContent(j+1) << " " << fDataHist_Slices[i]->GetBinError(j+1) << std::endl; - fDataHist->SetBinContent(bincount+1, fDataHist_Slices[i]->GetBinContent(j+1) ); - fDataHist->SetBinError(bincount+1, fDataHist_Slices[i]->GetBinError(j+1) ); - - bincount++; - } - - // Make MC Clones - fMCHist_Slices.push_back((TH1D*) fDataHist_Slices[i]->Clone()); - fMCHist_Slices[i]->SetNameTitle(Form("T2K_CC0pi_XSec_2DPcos_nu_nonuniform_MC_Slice%i",i), - (Form("T2K_CC0pi_XSec_2DPcos_nu_nonuniform_MC_Slice%i",i))); - - SetAutoProcessTH1(fDataHist_Slices[i],kCMD_Write); - SetAutoProcessTH1(fMCHist_Slices[i]); - // fMCHist_Slices[i]->Reset(); - } - - return; -}; diff --git a/src/T2K/T2K_CC0pinp_STV_XSec_1Ddat_nu.cxx b/src/T2K/T2K_CC0pinp_STV_XSec_1Ddat_nu.cxx index 956401e..c9d1b62 100644 --- a/src/T2K/T2K_CC0pinp_STV_XSec_1Ddat_nu.cxx +++ b/src/T2K/T2K_CC0pinp_STV_XSec_1Ddat_nu.cxx @@ -1,90 +1,90 @@ // 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_CC0pinp_STV_XSec_1Ddat_nu.h" //******************************************************************** T2K_CC0pinp_STV_XSec_1Ddat_nu::T2K_CC0pinp_STV_XSec_1Ddat_nu( nuiskey samplekey) { //******************************************************************** // Sample overview --------------------------------------------------- - std::string descrip = "T2K_CC0pinp_STV_XSec_1Ddat_nu sample. \n" + std::string descrip = "T2K_CC0pinp_STV_XSec_1Ddpt_nu sample. \n" "Target: CH \n" "Flux: T2K 2.5 degree off-axis (ND280) \n" - "Signal: CC0piNp (N>=1) with 450M eV < p_p < 1 GeV \n" + "Signal: CC0piNp (N>=1) with 450 MeV < p_p < 1 GeV \n" " p_mu > 250 MeV \n" - " cth_p > 0.6 \n" - " cth_mu > -0.4 \n" + " cth_p > 0.4 \n" + " cth_mu > -0.6 \n" "https://doi.org/10.1103/PhysRevD.98.032003 \n"; // Setup common settings fSettings = LoadSampleSettings(samplekey); fSettings.SetDescription(descrip); fSettings.SetXTitle("#delta#it{#alpha}_{T} (GeV c^{-1})"); fSettings.SetYTitle( "#frac{d#sigma}{d#delta#it{#alpha}_{T}} (cm^{2} nucleon^{-1} rads^{-1})"); fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/DIAG"); fSettings.SetEnuRange(0.0, 50.0); fSettings.DefineAllowedTargets("C,H"); // CCQELike plot information fSettings.SetTitle("T2K_CC0pinp_STV_XSec_1Ddat_nu"); // fSettings.SetDataInput( GeneralUtils::GetTopLevelDir() + // "/data/T2K/T2K_CC0pinp_STV_XSec_1Ddat_nu.dat"); fSettings.SetDataInput(FitPar::GetDataBase() + "/T2K/CC0pi/STV/datResults.root;Result"); fSettings.SetCovarInput(FitPar::GetDataBase() + "/T2K/CC0pi/STV/datResults.root;Correlation_Matrix"); fSettings.DefineAllowedSpecies("numu"); FinaliseSampleSettings(); // Scaling Setup --------------------------------------------------- // ScaleFactor automatically setup for DiffXSec/cm2/Nucleon fScaleFactor = GetEventHistogram()->Integral("width") * double(1E-38) / (double(fNEvents) * TotalIntegratedFlux("width")); // Plot Setup ------------------------------------------------------- // SetDataFromTextFile( fSettings.GetDataInput() ); // ScaleData(1E-38); // SetCovarFromDiagonal(); SetDataFromRootFile(fSettings.GetDataInput()); SetCorrelationFromRootFile(fSettings.GetCovarInput()); // Final setup --------------------------------------------------- FinaliseMeasurement(); }; void T2K_CC0pinp_STV_XSec_1Ddat_nu::FillEventVariables(FitEvent *event) { fXVar = FitUtils::Get_STV_dalphat(event, 14, true); return; }; //******************************************************************** bool T2K_CC0pinp_STV_XSec_1Ddat_nu::isSignal(FitEvent *event) //******************************************************************** { return SignalDef::isT2K_CC0pi_STV(event, EnuMin, EnuMax); } diff --git a/src/T2K/T2K_CC0pinp_STV_XSec_1Ddphit_nu.cxx b/src/T2K/T2K_CC0pinp_STV_XSec_1Ddphit_nu.cxx index cdc60dd..9dda5ff 100644 --- a/src/T2K/T2K_CC0pinp_STV_XSec_1Ddphit_nu.cxx +++ b/src/T2K/T2K_CC0pinp_STV_XSec_1Ddphit_nu.cxx @@ -1,90 +1,90 @@ // 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_CC0pinp_STV_XSec_1Ddphit_nu.h" //******************************************************************** T2K_CC0pinp_STV_XSec_1Ddphit_nu::T2K_CC0pinp_STV_XSec_1Ddphit_nu( nuiskey samplekey) { //******************************************************************** // Sample overview --------------------------------------------------- - std::string descrip = "T2K_CC0pinp_STV_XSec_1Ddphit_nu sample. \n" + std::string descrip = "T2K_CC0pinp_STV_XSec_1Ddpt_nu sample. \n" "Target: CH \n" "Flux: T2K 2.5 degree off-axis (ND280) \n" - "Signal: CC0piNp (N>=1) with 450M eV < p_p < 1 GeV \n" + "Signal: CC0piNp (N>=1) with 450 MeV < p_p < 1 GeV \n" " p_mu > 250 MeV \n" - " cth_p > 0.6 \n" - " cth_mu > -0.4 \n" + " cth_p > 0.4 \n" + " cth_mu > -0.6 \n" "https://doi.org/10.1103/PhysRevD.98.032003 \n"; // Setup common settings fSettings = LoadSampleSettings(samplekey); fSettings.SetDescription(descrip); fSettings.SetXTitle("#delta#it{#phi}_{T} (GeV c^{-1})"); fSettings.SetYTitle( "#frac{d#sigma}{d#delta#it{#phi}_{T}} (cm^{2} nucleon^{-1} rads^{-1})"); fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/DIAG"); fSettings.SetEnuRange(0.0, 50.0); fSettings.DefineAllowedTargets("C,H"); // CCQELike plot information fSettings.SetTitle("T2K_CC0pinp_STV_XSec_1Ddphit_nu"); // fSettings.SetDataInput( GeneralUtils::GetTopLevelDir() + // "/data/T2K/T2K_CC0pinp_STV_XSec_1Ddphit_nu.dat"); fSettings.SetDataInput(FitPar::GetDataBase() + "/T2K/CC0pi/STV/dphitResults.root;Result"); fSettings.SetCovarInput( FitPar::GetDataBase() + "/T2K/CC0pi/STV/dphitResults.root;Correlation_Matrix"); fSettings.DefineAllowedSpecies("numu"); FinaliseSampleSettings(); // Scaling Setup --------------------------------------------------- // ScaleFactor automatically setup for DiffXSec/cm2/Nucleon fScaleFactor = GetEventHistogram()->Integral("width") * double(1E-38) / (double(fNEvents) * TotalIntegratedFlux("width")); // Plot Setup ------------------------------------------------------- // SetDataFromTextFile( fSettings.GetDataInput() ); // ScaleData(1E-38); // SetCovarFromDiagonal(); SetDataFromRootFile(fSettings.GetDataInput()); SetCorrelationFromRootFile(fSettings.GetCovarInput()); // Final setup --------------------------------------------------- FinaliseMeasurement(); }; void T2K_CC0pinp_STV_XSec_1Ddphit_nu::FillEventVariables(FitEvent *event) { fXVar = FitUtils::Get_STV_dphit(event, 14, true); return; }; //******************************************************************** bool T2K_CC0pinp_STV_XSec_1Ddphit_nu::isSignal(FitEvent *event) //******************************************************************** { return SignalDef::isT2K_CC0pi_STV(event, EnuMin, EnuMax); } diff --git a/src/T2K/T2K_CC0pinp_STV_XSec_1Ddpt_nu.cxx b/src/T2K/T2K_CC0pinp_STV_XSec_1Ddpt_nu.cxx index 809db8a..65d6559 100644 --- a/src/T2K/T2K_CC0pinp_STV_XSec_1Ddpt_nu.cxx +++ b/src/T2K/T2K_CC0pinp_STV_XSec_1Ddpt_nu.cxx @@ -1,90 +1,90 @@ // 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_CC0pinp_STV_XSec_1Ddpt_nu.h" //******************************************************************** T2K_CC0pinp_STV_XSec_1Ddpt_nu::T2K_CC0pinp_STV_XSec_1Ddpt_nu( nuiskey samplekey) { //******************************************************************** // Sample overview --------------------------------------------------- std::string descrip = "T2K_CC0pinp_STV_XSec_1Ddpt_nu sample. \n" "Target: CH \n" "Flux: T2K 2.5 degree off-axis (ND280) \n" - "Signal: CC0piNp (N>=1) with 450M eV < p_p < 1 GeV \n" + "Signal: CC0piNp (N>=1) with 450 MeV < p_p < 1 GeV \n" " p_mu > 250 MeV \n" - " cth_p > 0.6 \n" - " cth_mu > -0.4 \n" + " cth_p > 0.4 \n" + " cth_mu > -0.6 \n" "https://doi.org/10.1103/PhysRevD.98.032003 \n"; // Setup common settings fSettings = LoadSampleSettings(samplekey); fSettings.SetDescription(descrip); fSettings.SetXTitle("#delta#it{p}_{T} (GeV c^{-1})"); fSettings.SetYTitle( "#frac{d#sigma}{d#delta#it{p}_{T}} (cm^{2} nucleon^{-1} GeV^{-1} c)"); fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/DIAG"); fSettings.SetEnuRange(0.0, 50.0); fSettings.DefineAllowedTargets("C,H"); // CCQELike plot information fSettings.SetTitle("T2K_CC0pinp_STV_XSec_1Ddpt_nu"); // fSettings.SetDataInput( GeneralUtils::GetTopLevelDir() + // "/data/T2K/T2K_CC0pinp_STV_XSec_1Ddpt_nu.dat"); fSettings.SetDataInput(FitPar::GetDataBase() + "/T2K/CC0pi/STV/dptResults.root;Result"); fSettings.SetCovarInput(FitPar::GetDataBase() + "/T2K/CC0pi/STV/dptResults.root;Correlation_Matrix"); fSettings.DefineAllowedSpecies("numu"); FinaliseSampleSettings(); // Scaling Setup --------------------------------------------------- // ScaleFactor automatically setup for DiffXSec/cm2/Nucleon fScaleFactor = GetEventHistogram()->Integral("width") * double(1E-38) / (double(fNEvents) * TotalIntegratedFlux("width")); // Plot Setup ------------------------------------------------------- // SetDataFromTextFile( fSettings.GetDataInput() ); // SetCovarFromDiagonal(); // ScaleData(1E-38); SetDataFromRootFile(fSettings.GetDataInput()); SetCorrelationFromRootFile(fSettings.GetCovarInput()); // SetCovarianceFromRootFile(fSettings.GetCovarInput() ); // Final setup --------------------------------------------------- FinaliseMeasurement(); }; void T2K_CC0pinp_STV_XSec_1Ddpt_nu::FillEventVariables(FitEvent *event) { fXVar = FitUtils::Get_STV_dpt(event, 14, true) / 1000.0; return; }; //******************************************************************** bool T2K_CC0pinp_STV_XSec_1Ddpt_nu::isSignal(FitEvent *event) //******************************************************************** { return SignalDef::isT2K_CC0pi_STV(event, EnuMin, EnuMax); } diff --git a/src/T2K/T2K_CC0pinp_ifk_XSec_3Dinfa_nu.cxx b/src/T2K/T2K_CC0pinp_ifk_XSec_3Dinfa_nu.cxx index 151c738..93513fe 100644 --- a/src/T2K/T2K_CC0pinp_ifk_XSec_3Dinfa_nu.cxx +++ b/src/T2K/T2K_CC0pinp_ifk_XSec_3Dinfa_nu.cxx @@ -1,208 +1,229 @@ // 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 . -*******************************************************************************/ + * 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_CC0pinp_ifk_XSec_3Dinfa_nu.h" - - -//******************************************************************** -T2K_CC0pinp_ifk_XSec_3Dinfa_nu::T2K_CC0pinp_ifk_XSec_3Dinfa_nu(nuiskey samplekey) { //******************************************************************** +T2K_CC0pinp_ifk_XSec_3Dinfa_nu::T2K_CC0pinp_ifk_XSec_3Dinfa_nu( + nuiskey samplekey) { + //******************************************************************** // Sample overview --------------------------------------------------- - std::string descrip = "T2K_CC0pinp_ifk_XSec_3Dinfa_nu sample. \n" \ - "Target: CH \n" \ - "Flux: T2K 2.5 degree off-axis (ND280) \n" \ - "Signal: CC0piNp (N>=1) with p_p>450MeV and cthp>0.4 \n"; + std::string descrip = "T2K_CC0pinp_ifk_XSec_3Dinfa_nu sample. \n" + "Target: CH \n" + "Flux: T2K 2.5 degree off-axis (ND280) \n" + "Signal: CC0piNp (N>=1) with p_p>450MeV and cthp>0.4 \n" + "https://doi.org/10.1103/PhysRevD.98.032003 \n"; // Setup common settings fSettings = LoadSampleSettings(samplekey); fSettings.SetDescription(descrip); fSettings.SetXTitle("#Delta #theta"); fSettings.SetYTitle("p_#mu"); 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.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; - outOfBoundsMC = 0.0; + outOfBoundsMC = 0.0; // CCQELike plot information fSettings.SetTitle("T2K_CC0pinp_ifk_XSec_3Dinfa_nu"); fSettings.DefineAllowedSpecies("numu"); FinaliseSampleSettings(); // Scaling Setup --------------------------------------------------- // ScaleFactor automatically setup for DiffXSec/cm2/Nucleon - fScaleFactor = ((GetEventHistogram()->Integral("width")/(fNEvents+0.)) / (TotalIntegratedFlux())); - //fScaleFactor = ((GetEventHistogram()->Integral("width")/(fNEvents+0.)) * 10 / (TotalIntegratedFlux())); + fScaleFactor = ((GetEventHistogram()->Integral("width") / (fNEvents + 0.)) / + (TotalIntegratedFlux())); + // fScaleFactor = ((GetEventHistogram()->Integral("width")/(fNEvents+0.)) * 10 + // / (TotalIntegratedFlux())); - fSettings.SetDataInput( FitPar::GetDataBase() + "/T2K/CC0pi/STV/infkResults_origBin.root;result_a" ); - SetDataFromRootFile( fSettings.GetDataInput() ); + fSettings.SetDataInput(FitPar::GetDataBase() + + "/T2K/CC0pi/STV/infkResults_origBin.root;result_a"); + SetDataFromRootFile(fSettings.GetDataInput()); - fSettings.SetCovarInput( FitPar::GetDataBase() + "/T2K/CC0pi/STV/infkResults_origBin.root;cor_a" ); - SetCorrelationFromRootFile(fSettings.GetCovarInput() ); - //SetCovarFromRootFile(FitPar::GetDataBase() + "/T2K/CC0pi/infkResults_origBin.root", "cov_a" ); - + fSettings.SetCovarInput(FitPar::GetDataBase() + + "/T2K/CC0pi/STV/infkResults_origBin.root;cor_a"); + SetCorrelationFromRootFile(fSettings.GetCovarInput()); + // SetCovarFromRootFile(FitPar::GetDataBase() + + // "/T2K/CC0pi/infkResults_origBin.root", "cov_a" ); // Setup Histograms SetHistograms(); // Final setup --------------------------------------------------- FinaliseMeasurement(); - }; - -bool T2K_CC0pinp_ifk_XSec_3Dinfa_nu::isSignal(FitEvent *event){ +bool T2K_CC0pinp_ifk_XSec_3Dinfa_nu::isSignal(FitEvent *event) { return SignalDef::isT2K_CC0pi_ifk(event, EnuMin, EnuMax); }; -void T2K_CC0pinp_ifk_XSec_3Dinfa_nu::FillEventVariables(FitEvent* event){ +void T2K_CC0pinp_ifk_XSec_3Dinfa_nu::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; + TLorentzVector Pp = event->GetHMFSParticle(2212)->fP; - double pmu = Pmu.Vect().Mag()/1000.; - //double pp = Pp.Vect().Mag()/1000.; + 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())); double infp_CosThetaP = FitUtils::cthpInfK(Pmu, CosThetaMu, 25, true); - double delta_a = acos(CosThetaP)*180/TMath::Pi() - acos(infp_CosThetaP)*180/TMath::Pi(); + double delta_a = acos(CosThetaP) * 180 / TMath::Pi() - + acos(infp_CosThetaP) * 180 / TMath::Pi(); // std::cout << "CosTheta P is: " << CosThetaP << std::endl; - // std::cout << "Inferred CosTheta P is: " << infp_CosThetaP << std::endl << std::endl; + // std::cout << "Inferred CosTheta P is: " << infp_CosThetaP << std::endl << + // std::endl; -/* - TVector3 tp_inf = FitUtils::tppInfK(Pmu, CosThetaMu, 25, true); + /* + TVector3 tp_inf = FitUtils::tppInfK(Pmu, CosThetaMu, 25, true); - //std::cout << "Theta P is: " << (Pp.Vect()).Theta() << std::endl; - //std::cout << "Inferred Theta P is: " << tp_inf.Theta() << std::endl << std::endl; + //std::cout << "Theta P is: " << (Pp.Vect()).Theta() << std::endl; + //std::cout << "Inferred Theta P is: " << tp_inf.Theta() << std::endl << + std::endl; - double delta_a = (Pp.Vect()).Theta()*180/TMath::Pi() - tp_inf.Theta()*180/TMath::Pi(); + double delta_a = (Pp.Vect()).Theta()*180/TMath::Pi() - + tp_inf.Theta()*180/TMath::Pi(); -*/ + */ fXVar = delta_a; fYVar = pmu; fZVar = CosThetaMu; return; }; -void T2K_CC0pinp_ifk_XSec_3Dinfa_nu::FillHistograms(){ +void T2K_CC0pinp_ifk_XSec_3Dinfa_nu::FillHistograms() { Measurement1D::FillHistograms(); - if (Signal){ - FillMCSlice( fXVar, fYVar, fZVar, Weight ); + if (Signal) { + FillMCSlice(fXVar, fYVar, fZVar, Weight); } - } +void T2K_CC0pinp_ifk_XSec_3Dinfa_nu::ConvertEventRates() { -void T2K_CC0pinp_ifk_XSec_3Dinfa_nu::ConvertEventRates(){ - - for (int i = 0; i < 7; i++){ + for (int i = 0; i < 7; 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); - + // 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(); - //The first bin in the histogram in underflow, so set this and start bincount at 1 + // The first bin in the histogram in underflow, so set this and start bincount + // at 1 fMCHist->SetBinContent(1, outOfBoundsMC); - int bincount = 1; - for (int i = 0; i < 7; i++){ - for (int j = 0; j < fDataHist_Slices[i]->GetNbinsX(); j++){ - fMCHist->SetBinContent(bincount+1, fMCHist_Slices[i]->GetBinContent(j+1)); - //fMCHist->SetBinError(bincount+1, fMCHist_Slices[i]->GetBinError(j+1)); + int bincount = 1; + for (int i = 0; i < 7; i++) { + for (int j = 0; j < fDataHist_Slices[i]->GetNbinsX(); 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_CC0pinp_ifk_XSec_3Dinfa_nu::FillMCSlice(double x, double y, double z, double w){ +void T2K_CC0pinp_ifk_XSec_3Dinfa_nu::FillMCSlice(double x, double y, double z, + double w) { // x is delta_a // y is pmu // z is CosThetaMu - if (z <= -0.6) fMCHist_Slices[0]->Fill(x,w); - else if (z >= -0.6 and z < 0.0 and y < 0.25) fMCHist_Slices[1]->Fill(x,w); - else if (z >= -0.6 and z < 0.0 and y > 0.25) fMCHist_Slices[2]->Fill(x,w); - else if (z >= 0.0 and y < 0.25) fMCHist_Slices[3]->Fill(x,w); - else if (z >= 0.0 and z < 0.8 and y >= 0.25) fMCHist_Slices[4]->Fill(x,w); - else if (z >= 0.8 and z < 1.0 and y >= 0.25 and y < 0.75) fMCHist_Slices[5]->Fill(x,w); - else if (z >= 0.8 and z < 1.0 and y >= 0.75) fMCHist_Slices[6]->Fill(x,w); - else outOfBoundsMC += w; - + if (z <= -0.6) + fMCHist_Slices[0]->Fill(x, w); + else if (z >= -0.6 and z < 0.0 and y < 0.25) + fMCHist_Slices[1]->Fill(x, w); + else if (z >= -0.6 and z < 0.0 and y > 0.25) + fMCHist_Slices[2]->Fill(x, w); + else if (z >= 0.0 and y < 0.25) + fMCHist_Slices[3]->Fill(x, w); + else if (z >= 0.0 and z < 0.8 and y >= 0.25) + fMCHist_Slices[4]->Fill(x, w); + else if (z >= 0.8 and z < 1.0 and y >= 0.25 and y < 0.75) + fMCHist_Slices[5]->Fill(x, w); + else if (z >= 0.8 and z < 1.0 and y >= 0.75) + fMCHist_Slices[6]->Fill(x, w); + else + outOfBoundsMC += w; } - -void T2K_CC0pinp_ifk_XSec_3Dinfa_nu::SetHistograms(){ +void T2K_CC0pinp_ifk_XSec_3Dinfa_nu::SetHistograms() { // Read in 1D Data Histograms - fInputFile = new TFile( (FitPar::GetDataBase() + "/T2K/CC0pi/STV/infkResults_origBin.root").c_str(),"READ"); - //fInputFile->ls(); - + fInputFile = new TFile( + (FitPar::GetDataBase() + "/T2K/CC0pi/STV/infkResults_origBin.root") + .c_str(), + "READ"); + // fInputFile->ls(); + // Read in 1D Data - fDataHist = (TH1D*) fInputFile->Get("result_a"); - fDataHist->SetNameTitle("T2K_CC0pinp_ifk_XSec_3Dinfa_nu_data", "T2K_CC0pinp_ifk_XSec_3Dinfa_nu_data"); - SetAutoProcessTH1(fDataHist,kCMD_Write); - + fDataHist = (TH1D *)fInputFile->Get("result_a"); + fDataHist->SetNameTitle("T2K_CC0pinp_ifk_XSec_3Dinfa_nu_data", + "T2K_CC0pinp_ifk_XSec_3Dinfa_nu_data"); + SetAutoProcessTH1(fDataHist, kCMD_Write); + // Read in 2D Data Slices and Make MC Slices - for (int i = 0; i < 7; i++){ //both y and z slices + for (int i = 0; i < 7; i++) { // both y and z slices // Get Data Histogram - //fInputFile->ls(); - fDataHist_Slices.push_back((TH1D*)fInputFile->Get(Form("resultBin%i_a",i))->Clone()); - fDataHist_Slices[i]->SetNameTitle(Form("T2K_CC0pinp_ifk_XSec_3Dinfa_nu_data_Slice%i",i), (Form("T2K_CC0pinp_ifk_XSec_3Dinfa_nu_data_Slice%i",i))); - + // fInputFile->ls(); + fDataHist_Slices.push_back( + (TH1D *)fInputFile->Get(Form("resultBin%i_a", i))->Clone()); + fDataHist_Slices[i]->SetNameTitle( + Form("T2K_CC0pinp_ifk_XSec_3Dinfa_nu_data_Slice%i", i), + (Form("T2K_CC0pinp_ifk_XSec_3Dinfa_nu_data_Slice%i", i))); + // Make MC Clones - fMCHist_Slices.push_back((TH1D*) fDataHist_Slices[i]->Clone()); - fMCHist_Slices[i]->SetNameTitle(Form("T2K_CC0pinp_ifk_XSec_3Dinfa_nu_MC_Slice%i",i), (Form("T2K_CC0pinp_ifk_XSec_3Dinfa_nu_MC_Slice%i",i))); + fMCHist_Slices.push_back((TH1D *)fDataHist_Slices[i]->Clone()); + fMCHist_Slices[i]->SetNameTitle( + Form("T2K_CC0pinp_ifk_XSec_3Dinfa_nu_MC_Slice%i", i), + (Form("T2K_CC0pinp_ifk_XSec_3Dinfa_nu_MC_Slice%i", i))); - SetAutoProcessTH1(fDataHist_Slices[i],kCMD_Write); + SetAutoProcessTH1(fDataHist_Slices[i], kCMD_Write); SetAutoProcessTH1(fMCHist_Slices[i]); fMCHist_Slices[i]->Reset(); } return; }; diff --git a/src/T2K/T2K_CC0pinp_ifk_XSec_3Dinfip_nu.cxx b/src/T2K/T2K_CC0pinp_ifk_XSec_3Dinfip_nu.cxx index 62aed1f..701aa3f 100644 --- a/src/T2K/T2K_CC0pinp_ifk_XSec_3Dinfip_nu.cxx +++ b/src/T2K/T2K_CC0pinp_ifk_XSec_3Dinfip_nu.cxx @@ -1,198 +1,217 @@ // 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 . -*******************************************************************************/ + * 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_CC0pinp_ifk_XSec_3Dinfip_nu.h" - - -//******************************************************************** -T2K_CC0pinp_ifk_XSec_3Dinfip_nu::T2K_CC0pinp_ifk_XSec_3Dinfip_nu(nuiskey samplekey) { //******************************************************************** +T2K_CC0pinp_ifk_XSec_3Dinfip_nu::T2K_CC0pinp_ifk_XSec_3Dinfip_nu( + nuiskey samplekey) { + //******************************************************************** // Sample overview --------------------------------------------------- - std::string descrip = "T2K_CC0pinp_ifk_XSec_3Dinfip_nu sample. \n" \ - "Target: CH \n" \ - "Flux: T2K 2.5 degree off-axis (ND280) \n" \ - "Signal: CC0piNp (N>=1) with p_p>450MeV and cthp>0.4 \n"; + std::string descrip = "T2K_CC0pinp_ifk_XSec_3Dinfip_nu sample. \n" + "Target: CH \n" + "Flux: T2K 2.5 degree off-axis (ND280) \n" + "Signal: CC0piNp (N>=1) with p_p>450MeV and cthp>0.4 \n" + "https://doi.org/10.1103/PhysRevD.98.032003 \n"; // Setup common settings fSettings = LoadSampleSettings(samplekey); fSettings.SetDescription(descrip); fSettings.SetXTitle("|#Delta p|"); fSettings.SetYTitle("p_#mu"); 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.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; - + outOfBoundsMC = 0.0; // CCQELike plot information fSettings.SetTitle("T2K_CC0pinp_ifk_XSec_3Dinfip_nu"); fSettings.DefineAllowedSpecies("numu"); FinaliseSampleSettings(); // Scaling Setup --------------------------------------------------- // ScaleFactor automatically setup for DiffXSec/cm2/Nucleon - fScaleFactor = ((GetEventHistogram()->Integral("width")/(fNEvents+0.)) / (TotalIntegratedFlux())); - //fScaleFactor = ((GetEventHistogram()->Integral("width")/(fNEvents+0.)) * 10 / (TotalIntegratedFlux())); + fScaleFactor = ((GetEventHistogram()->Integral("width") / (fNEvents + 0.)) / + (TotalIntegratedFlux())); + // fScaleFactor = ((GetEventHistogram()->Integral("width")/(fNEvents+0.)) * 10 + // / (TotalIntegratedFlux())); - fSettings.SetDataInput( FitPar::GetDataBase() + "/T2K/CC0pi/STV/infkResults_origBin.root;result_tp" ); - SetDataFromRootFile( fSettings.GetDataInput() ); + fSettings.SetDataInput(FitPar::GetDataBase() + + "/T2K/CC0pi/STV/infkResults_origBin.root;result_tp"); + SetDataFromRootFile(fSettings.GetDataInput()); - fSettings.SetCovarInput( FitPar::GetDataBase() + "/T2K/CC0pi/STV/infkResults_origBin.root;cor_tp" ); - SetCorrelationFromRootFile(fSettings.GetCovarInput() ); - //SetCovarFromRootFile(FitPar::GetDataBase() + "/T2K/CC0pi/infkResults_origBin.root", "cov_tp" ); + fSettings.SetCovarInput(FitPar::GetDataBase() + + "/T2K/CC0pi/STV/infkResults_origBin.root;cor_tp"); + SetCorrelationFromRootFile(fSettings.GetCovarInput()); + // SetCovarFromRootFile(FitPar::GetDataBase() + + // "/T2K/CC0pi/infkResults_origBin.root", "cov_tp" ); // Setup Histograms SetHistograms(); // Final setup --------------------------------------------------- FinaliseMeasurement(); - }; - -bool T2K_CC0pinp_ifk_XSec_3Dinfip_nu::isSignal(FitEvent *event){ +bool T2K_CC0pinp_ifk_XSec_3Dinfip_nu::isSignal(FitEvent *event) { return SignalDef::isT2K_CC0pi_ifk(event, EnuMin, EnuMax); }; -void T2K_CC0pinp_ifk_XSec_3Dinfip_nu::FillEventVariables(FitEvent* event){ +void T2K_CC0pinp_ifk_XSec_3Dinfip_nu::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; + TLorentzVector Pp = event->GetHMFSParticle(2212)->fP; - double pmu = Pmu.Vect().Mag()/1000.; - //double pp = Pp.Vect().Mag()/1000.; + double pmu = Pmu.Vect().Mag() / 1000.; + // double pp = Pp.Vect().Mag()/1000.; double CosThetaMu = cos(Pnu.Vect().Angle(Pmu.Vect())); TVector3 tp_inf = FitUtils::tppInfK(Pmu, CosThetaMu, 25, true); - TVector3 Pp_mev(Pp.X()/1000,Pp.Y()/1000,Pp.Z()/1000); + TVector3 Pp_mev(Pp.X() / 1000, Pp.Y() / 1000, Pp.Z() / 1000); - TVector3 delta_tp = tp_inf-Pp_mev; + TVector3 delta_tp = tp_inf - Pp_mev; - //std::cout << "Proton 3 mom is: " << std::endl; + // std::cout << "Proton 3 mom is: " << std::endl; //(Pp.Vect()).Print("all"); - //std::cout << "Inferred Proton 3 mom is: " << std::endl; - //tp_inf.Print("all"); - //std::cout << " " << std::endl; + // std::cout << "Inferred Proton 3 mom is: " << std::endl; + // tp_inf.Print("all"); + // std::cout << " " << std::endl; fXVar = delta_tp.Mag(); fYVar = pmu; fZVar = CosThetaMu; return; }; -void T2K_CC0pinp_ifk_XSec_3Dinfip_nu::FillHistograms(){ +void T2K_CC0pinp_ifk_XSec_3Dinfip_nu::FillHistograms() { Measurement1D::FillHistograms(); - if (Signal){ - FillMCSlice( fXVar, fYVar, fZVar, Weight ); + if (Signal) { + FillMCSlice(fXVar, fYVar, fZVar, Weight); } - } +void T2K_CC0pinp_ifk_XSec_3Dinfip_nu::ConvertEventRates() { -void T2K_CC0pinp_ifk_XSec_3Dinfip_nu::ConvertEventRates(){ - - for (int i = 0; i < 7; i++){ + for (int i = 0; i < 7; 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); - + // 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(); - //The first bin in the histogram in underflow, so set this and start bincount at 1 + // The first bin in the histogram in underflow, so set this and start bincount + // at 1 fMCHist->SetBinContent(1, outOfBoundsMC); - int bincount = 1; - for (int i = 0; i < 7; i++){ - for (int j = 0; j < fDataHist_Slices[i]->GetNbinsX(); j++){ - fMCHist->SetBinContent(bincount+1, fMCHist_Slices[i]->GetBinContent(j+1)); - //fMCHist->SetBinError(bincount+1, fMCHist_Slices[i]->GetBinError(j+1)); + int bincount = 1; + for (int i = 0; i < 7; i++) { + for (int j = 0; j < fDataHist_Slices[i]->GetNbinsX(); 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_CC0pinp_ifk_XSec_3Dinfip_nu::FillMCSlice(double x, double y, double z, double w){ +void T2K_CC0pinp_ifk_XSec_3Dinfip_nu::FillMCSlice(double x, double y, double z, + double w) { // x is delta_tp // y is pmu // z is CosThetaMu - if (z <= -0.6) fMCHist_Slices[0]->Fill(x,w); - else if (z >= -0.6 and z < 0.0 and y < 0.25) fMCHist_Slices[1]->Fill(x,w); - else if (z >= -0.6 and z < 0.0 and y > 0.25) fMCHist_Slices[2]->Fill(x,w); - else if (z >= 0.0 and y < 0.25) fMCHist_Slices[3]->Fill(x,w); - else if (z >= 0.0 and z < 0.8 and y >= 0.25) fMCHist_Slices[4]->Fill(x,w); - else if (z >= 0.8 and z < 1.0 and y >= 0.25 and y < 0.75) fMCHist_Slices[5]->Fill(x,w); - else if (z >= 0.8 and z < 1.0 and y >= 0.75) fMCHist_Slices[6]->Fill(x,w); - else outOfBoundsMC += w; - + if (z <= -0.6) + fMCHist_Slices[0]->Fill(x, w); + else if (z >= -0.6 and z < 0.0 and y < 0.25) + fMCHist_Slices[1]->Fill(x, w); + else if (z >= -0.6 and z < 0.0 and y > 0.25) + fMCHist_Slices[2]->Fill(x, w); + else if (z >= 0.0 and y < 0.25) + fMCHist_Slices[3]->Fill(x, w); + else if (z >= 0.0 and z < 0.8 and y >= 0.25) + fMCHist_Slices[4]->Fill(x, w); + else if (z >= 0.8 and z < 1.0 and y >= 0.25 and y < 0.75) + fMCHist_Slices[5]->Fill(x, w); + else if (z >= 0.8 and z < 1.0 and y >= 0.75) + fMCHist_Slices[6]->Fill(x, w); + else + outOfBoundsMC += w; } - -void T2K_CC0pinp_ifk_XSec_3Dinfip_nu::SetHistograms(){ +void T2K_CC0pinp_ifk_XSec_3Dinfip_nu::SetHistograms() { // Read in 1D Data Histograms - fInputFile = new TFile( (FitPar::GetDataBase() + "/T2K/CC0pi/STV/infkResults_origBin.root").c_str(),"READ"); - //fInputFile->ls(); - + fInputFile = new TFile( + (FitPar::GetDataBase() + "/T2K/CC0pi/STV/infkResults_origBin.root") + .c_str(), + "READ"); + // fInputFile->ls(); + // Read in 1D Data - fDataHist = (TH1D*) fInputFile->Get("result_tp"); - fDataHist->SetNameTitle("T2K_CC0pinp_ifk_XSec_3Dinfip_nu_data", "T2K_CC0pinp_ifk_XSec_3Dinfip_nu_data"); - SetAutoProcessTH1(fDataHist,kCMD_Write); - + fDataHist = (TH1D *)fInputFile->Get("result_tp"); + fDataHist->SetNameTitle("T2K_CC0pinp_ifk_XSec_3Dinfip_nu_data", + "T2K_CC0pinp_ifk_XSec_3Dinfip_nu_data"); + SetAutoProcessTH1(fDataHist, kCMD_Write); + // Read in 2D Data Slices and Make MC Slices - for (int i = 0; i < 7; i++){ //both y and z slices + for (int i = 0; i < 7; i++) { // both y and z slices // Get Data Histogram - //fInputFile->ls(); - fDataHist_Slices.push_back((TH1D*)fInputFile->Get(Form("resultBin%i_tp",i))->Clone()); - fDataHist_Slices[i]->SetNameTitle(Form("T2K_CC0pinp_ifk_XSec_3Dinfip_nu_data_Slice%i",i), (Form("T2K_CC0pinp_ifk_XSec_3Dinfip_nu_data_Slice%i",i))); + // fInputFile->ls(); + fDataHist_Slices.push_back( + (TH1D *)fInputFile->Get(Form("resultBin%i_tp", i))->Clone()); + fDataHist_Slices[i]->SetNameTitle( + Form("T2K_CC0pinp_ifk_XSec_3Dinfip_nu_data_Slice%i", i), + (Form("T2K_CC0pinp_ifk_XSec_3Dinfip_nu_data_Slice%i", i))); // Make MC Clones - fMCHist_Slices.push_back((TH1D*) fDataHist_Slices[i]->Clone()); - fMCHist_Slices[i]->SetNameTitle(Form("T2K_CC0pinp_ifk_XSec_3Dinfip_nu_MC_Slice%i",i), (Form("T2K_CC0pinp_ifk_XSec_3Dinfip_nu_MC_Slice%i",i))); + fMCHist_Slices.push_back((TH1D *)fDataHist_Slices[i]->Clone()); + fMCHist_Slices[i]->SetNameTitle( + Form("T2K_CC0pinp_ifk_XSec_3Dinfip_nu_MC_Slice%i", i), + (Form("T2K_CC0pinp_ifk_XSec_3Dinfip_nu_MC_Slice%i", i))); - SetAutoProcessTH1(fDataHist_Slices[i],kCMD_Write); + SetAutoProcessTH1(fDataHist_Slices[i], kCMD_Write); SetAutoProcessTH1(fMCHist_Slices[i]); fMCHist_Slices[i]->Reset(); } return; }; diff --git a/src/T2K/T2K_CC0pinp_ifk_XSec_3Dinfp_nu.cxx b/src/T2K/T2K_CC0pinp_ifk_XSec_3Dinfp_nu.cxx index 78a22d5..c5a3fc3 100644 --- a/src/T2K/T2K_CC0pinp_ifk_XSec_3Dinfp_nu.cxx +++ b/src/T2K/T2K_CC0pinp_ifk_XSec_3Dinfp_nu.cxx @@ -1,187 +1,188 @@ // 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_CC0pinp_ifk_XSec_3Dinfp_nu.h" //******************************************************************** T2K_CC0pinp_ifk_XSec_3Dinfp_nu::T2K_CC0pinp_ifk_XSec_3Dinfp_nu(nuiskey samplekey) { //******************************************************************** // Sample overview --------------------------------------------------- std::string descrip = "T2K_CC0pinp_ifk_XSec_3Dinfp_nu sample. \n" \ "Target: CH \n" \ "Flux: T2K 2.5 degree off-axis (ND280) \n" \ - "Signal: CC0piNp (N>=1) with p_p>450MeV and cthp>0.4 \n"; + "Signal: CC0piNp (N>=1) with p_p>450MeV and cthp>0.4 \n" + "https://doi.org/10.1103/PhysRevD.98.032003 \n"; // Setup common settings fSettings = LoadSampleSettings(samplekey); fSettings.SetDescription(descrip); fSettings.SetXTitle("#Delta p"); fSettings.SetYTitle("p_#mu"); 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; outOfBoundsMC = 0.0; // CCQELike plot information fSettings.SetTitle("T2K_CC0pinp_ifk_XSec_3Dinfp_nu"); fSettings.DefineAllowedSpecies("numu"); FinaliseSampleSettings(); // Scaling Setup --------------------------------------------------- // ScaleFactor automatically setup for DiffXSec/cm2/Nucleon fScaleFactor = ((GetEventHistogram()->Integral("width")/(fNEvents+0.)) / (TotalIntegratedFlux())); //fScaleFactor = ((GetEventHistogram()->Integral("width")/(fNEvents+0.)) * 10 / (TotalIntegratedFlux())); fSettings.SetDataInput( FitPar::GetDataBase() + "/T2K/CC0pi/STV/infkResults_origBin.root;result_p" ); SetDataFromRootFile( fSettings.GetDataInput() ); fSettings.SetCovarInput( FitPar::GetDataBase() + "/T2K/CC0pi/STV/infkResults_origBin.root;cor_p" ); SetCorrelationFromRootFile(fSettings.GetCovarInput() ); //SetCovarFromRootFile(FitPar::GetDataBase() + "/T2K/CC0pi/infkResults_origBin.root", "cov_p" ); // Setup Histograms SetHistograms(); // Final setup --------------------------------------------------- FinaliseMeasurement(); }; bool T2K_CC0pinp_ifk_XSec_3Dinfp_nu::isSignal(FitEvent *event){ return SignalDef::isT2K_CC0pi_ifk(event, EnuMin, EnuMax); }; void T2K_CC0pinp_ifk_XSec_3Dinfp_nu::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 CosThetaMu = cos(Pnu.Vect().Angle(Pmu.Vect())); double delta_p = Pp.Vect().Mag()/1000. - FitUtils::ppInfK(Pmu, CosThetaMu, 25, true); fXVar = delta_p; fYVar = pmu; fZVar = CosThetaMu; return; }; void T2K_CC0pinp_ifk_XSec_3Dinfp_nu::FillHistograms(){ Measurement1D::FillHistograms(); if (Signal){ FillMCSlice( fXVar, fYVar, fZVar, Weight ); } } void T2K_CC0pinp_ifk_XSec_3Dinfp_nu::ConvertEventRates(){ for (int i = 0; i < 7; 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(); //The first bin in the histogram in underflow, so set this and start bincount at 1 fMCHist->SetBinContent(1, outOfBoundsMC); int bincount = 1; for (int i = 0; i < 7; i++){ for (int j = 0; j < fDataHist_Slices[i]->GetNbinsX(); 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_CC0pinp_ifk_XSec_3Dinfp_nu::FillMCSlice(double x, double y, double z, double w){ // x is delta_p // y is pmu // z is CosThetaMu if (z <= -0.6) fMCHist_Slices[0]->Fill(x,w); else if (z >= -0.6 and z < 0.0 and y < 0.25) fMCHist_Slices[1]->Fill(x,w); else if (z >= -0.6 and z < 0.0 and y > 0.25) fMCHist_Slices[2]->Fill(x,w); else if (z >= 0.0 and y < 0.25) fMCHist_Slices[3]->Fill(x,w); else if (z >= 0.0 and z < 0.8 and y >= 0.25) fMCHist_Slices[4]->Fill(x,w); else if (z >= 0.8 and z < 1.0 and y >= 0.25 and y < 0.75) fMCHist_Slices[5]->Fill(x,w); else if (z >= 0.8 and z < 1.0 and y >= 0.75) fMCHist_Slices[6]->Fill(x,w); else outOfBoundsMC += w; } void T2K_CC0pinp_ifk_XSec_3Dinfp_nu::SetHistograms(){ // Read in 1D Data Histograms fInputFile = new TFile( (FitPar::GetDataBase() + "/T2K/CC0pi/STV/infkResults_origBin.root").c_str(),"READ"); //fInputFile->ls(); // Read in 1D Data fDataHist = (TH1D*) fInputFile->Get("result_p"); fDataHist->SetNameTitle("T2K_CC0pinp_ifk_XSec_3Dinfp_nu_data", "T2K_CC0pinp_ifk_XSec_3Dinfp_nu_data"); SetAutoProcessTH1(fDataHist,kCMD_Write); // Read in 2D Data Slices and Make MC Slices for (int i = 0; i < 7; i++){ //both y and z slices // Get Data Histogram //fInputFile->ls(); fDataHist_Slices.push_back((TH1D*)fInputFile->Get(Form("resultBin%i_p",i))->Clone()); fDataHist_Slices[i]->SetNameTitle(Form("T2K_CC0pinp_ifk_XSec_3Dinfp_nu_data_Slice%i",i), (Form("T2K_CC0pinp_ifk_XSec_3Dinfp_nu_data_Slice%i",i))); // Make MC Clones fMCHist_Slices.push_back((TH1D*) fDataHist_Slices[i]->Clone()); fMCHist_Slices[i]->SetNameTitle(Form("T2K_CC0pinp_ifk_XSec_3Dinfp_nu_MC_Slice%i",i), (Form("T2K_CC0pinp_ifk_XSec_3Dinfp_nu_MC_Slice%i",i))); SetAutoProcessTH1(fDataHist_Slices[i],kCMD_Write); SetAutoProcessTH1(fMCHist_Slices[i]); fMCHist_Slices[i]->Reset(); } return; }; diff --git a/src/T2K/T2K_SignalDef.cxx b/src/T2K/T2K_SignalDef.cxx index b407cb6..24edfa1 100644 --- a/src/T2K/T2K_SignalDef.cxx +++ b/src/T2K/T2K_SignalDef.cxx @@ -1,343 +1,342 @@ // 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) { +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; } } int npicuts = 0; npicuts += bool(pscuts & kPionFwdHighMom); npicuts += bool(pscuts & kPionHighMom); npicuts += bool(pscuts & kPionHighEff); if (npicuts != 1) { NUIS_ABORT( "isCC1pip_T2K_arxiv1909_03936 signal definition passed incompatible " "pion phase space cuts. Should be either kMuonHighEff, or one of " "kPionFwdHighMom,kPionHighMom, or kPionHighEff"); } TLorentzVector Ppip = event->GetHMFSParticle(211)->fP; double cos_th_pi = cos(FitUtils::th(Pnu, Ppip)); double p_pi = FitUtils::p(Ppip) * 1000; if (pscuts & kPionFwdHighMom) { return ((cos_th_pi > 0) && (p_pi > 200)); } if (pscuts & kPionHighMom) { return (p_pi > 200); } if (pscuts & kMuonHighEff) { return ((cos_th_pi > 0.0) && (p_pi > 200)); } return false; }; -bool isT2K_CC0pi(FitEvent *event, double EnuMin, double EnuMax, - bool forwardgoing) { +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 (forwardgoing) { - if (CosThetaMu < 0.0 || p_mu < 200) + 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 - TN328 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 && CosThetaMu < 0.85) return false; if (Pmu < 670 && CosThetaMu < 0.88) return false; if (Pmu < 800 && CosThetaMu < 0.9) return false; if (Pmu < 1000 && CosThetaMu < 0.91) return false; if (Pmu < 1380 && CosThetaMu < 0.92) return false; if (Pmu < 2010 && CosThetaMu < 0.95) return false; return true; } bool isT2K_CC0piNp(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 < 2) 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 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; } 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 bfc68f5..aec9b64 100644 --- a/src/T2K/T2K_SignalDef.h +++ b/src/T2K/T2K_SignalDef.h @@ -1,50 +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), kMuonHighEff = (1 << 1), kPionFwdHighMom = (1 << 2), kPionHighMom = (1 << 3), kPionHighEff = (1 << 4) }; bool isCC1pip_T2K_arxiv1909_03936(FitEvent *event, double EnuMin, double EnuMax, - int); + int cuts); + +enum PRD93112012_Ana { + kAnalysis_I, + kAnalysis_II, +}; bool isT2K_CC0pi(FitEvent *event, double EnuMin, double EnuMax, - bool forwardgoing); + 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