diff --git a/src/FCN/SampleList.cxx b/src/FCN/SampleList.cxx index cab1e5c..e5df4b5 100644 --- a/src/FCN/SampleList.cxx +++ b/src/FCN/SampleList.cxx @@ -1,1374 +1,1378 @@ #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 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_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 +// 2018 MINERvA CC0pi STV #include "MINERvA_CC0pinp_STV_XSec_1D_nu.h" +// 2018 MINERvA CC0pi 2D +#include "MINERvA_CC0pi_XSec_2D_nu.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" -#include "MINERvA_CC0pi_XSec_2Dptpx_antinu.h" -#include "MINERvA_CC0pi_XSec_2Dptpx_nu.h" +//#include "MINERvA_CC0pi_XSec_2Dptpx_antinu.h" +//#include "MINERvA_CC0pi_XSec_2Dptpx_nu.h" #endif #ifndef __NO_T2K__ // T2K CC0pi #include "T2K_CC0pi_XSec_2DPcos_nu.h" // T2K STV CC0pi #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_ifk_XSec_3Dinfa_nu.h" #include "T2K_CC0pinp_ifk_XSec_3Dinfip_nu.h" // T2K CC1pi+ on CH #include "T2K_CC1pip_CH_XSec_1DQ2_nu.h" #include "T2K_CC1pip_CH_XSec_1DWrec_nu.h" #include "T2K_CC1pip_CH_XSec_1Dpmu_nu.h" #include "T2K_CC1pip_CH_XSec_1Dppi_nu.h" #include "T2K_CC1pip_CH_XSec_1Dq3_nu.h" #include "T2K_CC1pip_CH_XSec_1Dthmupi_nu.h" #include "T2K_CC1pip_CH_XSec_1Dthpi_nu.h" #include "T2K_CC1pip_CH_XSec_1Dthq3pi_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_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(); QLOG(FIT, "Loaded " << NSamples << " from " << NManifests << " shared object libraries."); } DynamicSampleFactory* DynamicSampleFactory::glblDSF = NULL; DynamicSampleFactory::PluginManifest::~PluginManifest() { for (size_t i_it = 0; i_it < Instances.size(); ++i_it) { (*(DSF_DestroySample))(Instances[i_it]); } } std::string EnsureTrailingSlash(std::string const& inp) { if (!inp.length()) { return "/"; } if (inp[inp.length() - 1] == '/') { return inp; } return inp + "/"; } void DynamicSampleFactory::LoadPlugins() { std::vector SearchDirectories; if (Config::HasPar("dynamic_sample.path")) { SearchDirectories = GeneralUtils::ParseToStr(Config::GetParS("dynamic_sample.path"), ":"); } char const* envPath = getenv("NUISANCE_DS_PATH"); if (envPath) { std::vector envPaths = GeneralUtils::ParseToStr(envPath, ":"); for (size_t ep_it = 0; ep_it < envPaths.size(); ++ep_it) { SearchDirectories.push_back(envPaths[ep_it]); } } if (!SearchDirectories.size()) { char const* pwdPath = getenv("PWD"); if (pwdPath) { SearchDirectories.push_back(pwdPath); } } for (size_t sp_it = 0; sp_it < SearchDirectories.size(); ++sp_it) { std::string dirpath = EnsureTrailingSlash(SearchDirectories[sp_it]); QLOG(FIT, "Searching for dynamic sample manifests in: " << dirpath); Ssiz_t len = 0; DIR* dir; struct dirent* ent; dir = opendir(dirpath.c_str()); if (dir != NULL) { TRegexp matchExp("*.so", true); while ((ent = readdir(dir)) != NULL) { if (matchExp.Index(TString(ent->d_name), &len) != Ssiz_t(-1)) { QLOG(FIT, "\tFound shared object: " << ent->d_name << " checking for relevant methods..."); void* dlobj = dlopen((dirpath + ent->d_name).c_str(), RTLD_NOW | RTLD_GLOBAL); char const* dlerr_cstr = dlerror(); std::string dlerr; if (dlerr_cstr) { dlerr = dlerr_cstr; } if (dlerr.length()) { ERROR(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()) { ERROR(WRN, "\tFailed to load symbol \"DSF_NSamples\" from " << (dirpath + ent->d_name) << ": " << dlerr); dlclose(dlobj); continue; } plgManif.DSF_GetSampleName = reinterpret_cast( dlsym(dlobj, "DSF_GetSampleName")); dlerr = ""; dlerr_cstr = dlerror(); if (dlerr_cstr) { dlerr = dlerr_cstr; } if (dlerr.length()) { ERROR(WRN, "\tFailed to load symbol \"DSF_GetSampleName\" from " << (dirpath + ent->d_name) << ": " << dlerr); dlclose(dlobj); continue; } plgManif.DSF_GetSample = reinterpret_cast( dlsym(dlobj, "DSF_GetSample")); dlerr = ""; dlerr_cstr = dlerror(); if (dlerr_cstr) { dlerr = dlerr_cstr; } if (dlerr.length()) { ERROR(WRN, "\tFailed to load symbol \"DSF_GetSample\" from " << (dirpath + ent->d_name) << ": " << dlerr); dlclose(dlobj); continue; } plgManif.DSF_DestroySample = reinterpret_cast( dlsym(dlobj, "DSF_DestroySample")); dlerr = ""; dlerr_cstr = dlerror(); if (dlerr_cstr) { dlerr = dlerr_cstr; } if (dlerr.length()) { ERROR(WRN, "Failed to load symbol \"DSF_DestroySample\" from " << (dirpath + ent->d_name) << ": " << dlerr); dlclose(dlobj); continue; } plgManif.NSamples = (*(plgManif.DSF_NSamples))(); QLOG(FIT, "\tSuccessfully loaded dynamic sample manifest: " << plgManif.soloc << ". Contains " << plgManif.NSamples << " samples."); for (size_t smp_it = 0; smp_it < plgManif.NSamples; ++smp_it) { char const* smp_name = (*(plgManif.DSF_GetSampleName))(smp_it); if (!smp_name) { THROW("Could not load sample " << smp_it << " / " << plgManif.NSamples << " from " << plgManif.soloc); } if (Samples.count(smp_name)) { ERROR(WRN, "Already loaded a sample named: \"" << smp_name << "\". cannot load duplciates. This " "sample will be skipped."); continue; } plgManif.SamplesProvided.push_back(smp_name); Samples[smp_name] = std::make_pair(plgManif.soloc, smp_it); QLOG(FIT, "\t\t" << smp_name); } if (plgManif.SamplesProvided.size()) { Manifests[plgManif.soloc] = plgManif; NSamples += plgManif.SamplesProvided.size(); NManifests++; } else { dlclose(dlobj); } } } closedir(dir); } else { ERROR(WRN, "Tried to open non-existant directory."); } } } DynamicSampleFactory& DynamicSampleFactory::Get() { if (!glblDSF) { glblDSF = new DynamicSampleFactory(); } return *glblDSF; } void DynamicSampleFactory::Print() { std::map > ManifestSamples; for (std::map >::iterator smp_it = Samples.begin(); smp_it != Samples.end(); ++smp_it) { if (!ManifestSamples.count(smp_it->second.first)) { ManifestSamples[smp_it->second.first] = std::vector(); } ManifestSamples[smp_it->second.first].push_back(smp_it->first); } QLOG(FIT, "Dynamic sample manifest: "); for (std::map >::iterator m_it = ManifestSamples.begin(); m_it != ManifestSamples.end(); ++m_it) { QLOG(FIT, "\tLibrary " << m_it->first << " contains: "); for (size_t s_it = 0; s_it < m_it->second.size(); ++s_it) { QLOG(FIT, "\t\t" << m_it->second[s_it]); } } } bool DynamicSampleFactory::HasSample(std::string const& name) { return Samples.count(name); } bool DynamicSampleFactory::HasSample(nuiskey& samplekey) { return HasSample(samplekey.GetS("name")); } MeasurementBase* DynamicSampleFactory::CreateSample(nuiskey& samplekey) { if (!HasSample(samplekey)) { ERROR(WRN, "Asked to load unknown sample: \"" << samplekey.GetS("name") << "\"."); return NULL; } std::pair sample = Samples[samplekey.GetS("name")]; QLOG(SAM, "\tLoading sample " << sample.second << " from " << sample.first); return (*(Manifests[sample.first].DSF_GetSample))(sample.second, &samplekey); } DynamicSampleFactory::~DynamicSampleFactory() { Manifests.clear(); } #endif //! Functions to make it easier for samples to be created and handled. namespace SampleUtils { //! Create a given sample given its name, file, type, fakdata(fkdt) file and the //! current rw engine and push it back into the list fChain. MeasurementBase* CreateSample(std::string name, std::string file, std::string type, std::string fkdt, FitWeight* rw) { nuiskey samplekey = Config::CreateKey("sample"); samplekey.Set("name", name); samplekey.Set("input", file); samplekey.Set("type", type); return CreateSample(samplekey); } MeasurementBase* CreateSample(nuiskey samplekey) { #ifdef __USE_DYNSAMPLES__ if (DynamicSampleFactory::Get().HasSample(samplekey)) { QLOG(SAM, "Instantiating dynamic sample..."); MeasurementBase* ds = DynamicSampleFactory::Get().CreateSample(samplekey); if (ds) { QLOG(SAM, "Done."); return ds; } THROW("Failed to instantiate dynamic sample."); } #endif FitWeight* rw = FitBase::GetRW(); std::string name = samplekey.GetS("name"); std::string file = samplekey.GetS("input"); std::string type = samplekey.GetS("type"); std::string fkdt = ""; /* ANL CCQE Samples */ #ifndef __NO_ANL__ if (!name.compare("ANL_CCQE_XSec_1DEnu_nu") || !name.compare("ANL_CCQE_XSec_1DEnu_nu_PRD26") || !name.compare("ANL_CCQE_XSec_1DEnu_nu_PRL31") || !name.compare("ANL_CCQE_XSec_1DEnu_nu_PRD16")) { return (new ANL_CCQE_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("ANL_CCQE_Evt_1DQ2_nu") || !name.compare("ANL_CCQE_Evt_1DQ2_nu_PRL31") || !name.compare("ANL_CCQE_Evt_1DQ2_nu_PRD26") || !name.compare("ANL_CCQE_Evt_1DQ2_nu_PRD16")) { return (new ANL_CCQE_Evt_1DQ2_nu(samplekey)); /* ANL CC1ppip samples */ } else if (!name.compare("ANL_CC1ppip_XSec_1DEnu_nu") || !name.compare("ANL_CC1ppip_XSec_1DEnu_nu_W14Cut") || !name.compare("ANL_CC1ppip_XSec_1DEnu_nu_Uncorr") || !name.compare("ANL_CC1ppip_XSec_1DEnu_nu_W14Cut_Uncorr") || !name.compare("ANL_CC1ppip_XSec_1DEnu_nu_W16Cut_Uncorr")) { return (new ANL_CC1ppip_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("ANL_CC1ppip_XSec_1DQ2_nu")) { return (new ANL_CC1ppip_XSec_1DQ2_nu(samplekey)); } else if (!name.compare("ANL_CC1ppip_Evt_1DQ2_nu") || !name.compare("ANL_CC1ppip_Evt_1DQ2_nu_W14Cut")) { return (new ANL_CC1ppip_Evt_1DQ2_nu(samplekey)); } else if (!name.compare("ANL_CC1ppip_Evt_1Dppi_nu")) { return (new ANL_CC1ppip_Evt_1Dppi_nu(samplekey)); } else if (!name.compare("ANL_CC1ppip_Evt_1Dthpr_nu")) { return (new ANL_CC1ppip_Evt_1Dthpr_nu(samplekey)); } else if (!name.compare("ANL_CC1ppip_Evt_1DcosmuStar_nu")) { return (new ANL_CC1ppip_Evt_1DcosmuStar_nu(samplekey)); } else if (!name.compare("ANL_CC1ppip_Evt_1DcosthAdler_nu")) { return (new ANL_CC1ppip_Evt_1DcosthAdler_nu(samplekey)); } else if (!name.compare("ANL_CC1ppip_Evt_1Dphi_nu")) { return (new ANL_CC1ppip_Evt_1Dphi_nu(samplekey)); /* ANL CC1npip sample */ } else if (!name.compare("ANL_CC1npip_XSec_1DEnu_nu") || !name.compare("ANL_CC1npip_XSec_1DEnu_nu_W14Cut") || !name.compare("ANL_CC1npip_XSec_1DEnu_nu_Uncorr") || !name.compare("ANL_CC1npip_XSec_1DEnu_nu_W14Cut_Uncorr") || !name.compare("ANL_CC1npip_XSec_1DEnu_nu_W16Cut_Uncorr")) { return (new ANL_CC1npip_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("ANL_CC1npip_Evt_1DQ2_nu") || !name.compare("ANL_CC1npip_Evt_1DQ2_nu_W14Cut")) { return (new ANL_CC1npip_Evt_1DQ2_nu(samplekey)); } else if (!name.compare("ANL_CC1npip_Evt_1Dppi_nu")) { return (new ANL_CC1npip_Evt_1Dppi_nu(samplekey)); } else if (!name.compare("ANL_CC1npip_Evt_1DcosmuStar_nu")) { return (new ANL_CC1npip_Evt_1DcosmuStar_nu(samplekey)); /* ANL CC1pi0 sample */ } else if (!name.compare("ANL_CC1pi0_XSec_1DEnu_nu") || !name.compare("ANL_CC1pi0_XSec_1DEnu_nu_W14Cut") || !name.compare("ANL_CC1pi0_XSec_1DEnu_nu_Uncorr") || !name.compare("ANL_CC1pi0_XSec_1DEnu_nu_W14Cut_Uncorr") || !name.compare("ANL_CC1pi0_XSec_1DEnu_nu_W16Cut_Uncorr")) { return (new ANL_CC1pi0_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("ANL_CC1pi0_Evt_1DQ2_nu") || !name.compare("ANL_CC1pi0_Evt_1DQ2_nu_W14Cut")) { return (new ANL_CC1pi0_Evt_1DQ2_nu(samplekey)); } else if (!name.compare("ANL_CC1pi0_Evt_1DcosmuStar_nu")) { return (new ANL_CC1pi0_Evt_1DcosmuStar_nu(samplekey)); /* ANL NC1npip sample */ } else if (!name.compare("ANL_NC1npip_Evt_1Dppi_nu")) { return (new ANL_NC1npip_Evt_1Dppi_nu(samplekey)); /* ANL NC1ppim sample */ } else if (!name.compare("ANL_NC1ppim_XSec_1DEnu_nu")) { return (new ANL_NC1ppim_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("ANL_NC1ppim_Evt_1DcosmuStar_nu")) { return (new ANL_NC1ppim_Evt_1DcosmuStar_nu(samplekey)); /* ANL CC2pi sample */ } else if (!name.compare("ANL_CC2pi_1pim1pip_XSec_1DEnu_nu")) { return (new ANL_CC2pi_1pim1pip_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("ANL_CC2pi_1pim1pip_Evt_1Dpmu_nu")) { return (new ANL_CC2pi_1pim1pip_Evt_1Dpmu_nu(samplekey)); } else if (!name.compare("ANL_CC2pi_1pim1pip_Evt_1Dppip_nu")) { return (new ANL_CC2pi_1pim1pip_Evt_1Dppip_nu(samplekey)); } else if (!name.compare("ANL_CC2pi_1pim1pip_Evt_1Dppim_nu")) { return (new ANL_CC2pi_1pim1pip_Evt_1Dppim_nu(samplekey)); } else if (!name.compare("ANL_CC2pi_1pim1pip_Evt_1Dpprot_nu")) { return (new ANL_CC2pi_1pim1pip_Evt_1Dpprot_nu(samplekey)); } else if (!name.compare("ANL_CC2pi_1pip1pip_XSec_1DEnu_nu")) { return (new ANL_CC2pi_1pip1pip_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("ANL_CC2pi_1pip1pip_Evt_1Dpmu_nu")) { return (new ANL_CC2pi_1pip1pip_Evt_1Dpmu_nu(samplekey)); } else if (!name.compare("ANL_CC2pi_1pip1pip_Evt_1Dpneut_nu")) { return (new ANL_CC2pi_1pip1pip_Evt_1Dpneut_nu(samplekey)); } else if (!name.compare("ANL_CC2pi_1pip1pip_Evt_1DppipHigh_nu")) { return (new ANL_CC2pi_1pip1pip_Evt_1DppipHigh_nu(samplekey)); } else if (!name.compare("ANL_CC2pi_1pip1pip_Evt_1DppipLow_nu")) { return (new ANL_CC2pi_1pip1pip_Evt_1DppipLow_nu(samplekey)); } else if (!name.compare("ANL_CC2pi_1pip1pi0_XSec_1DEnu_nu")) { return (new ANL_CC2pi_1pip1pi0_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("ANL_CC2pi_1pip1pi0_Evt_1Dpmu_nu")) { return (new ANL_CC2pi_1pip1pi0_Evt_1Dpmu_nu(samplekey)); } else if (!name.compare("ANL_CC2pi_1pip1pi0_Evt_1Dppip_nu")) { return (new ANL_CC2pi_1pip1pi0_Evt_1Dppip_nu(samplekey)); } else if (!name.compare("ANL_CC2pi_1pip1pi0_Evt_1Dppi0_nu")) { return (new ANL_CC2pi_1pip1pi0_Evt_1Dppi0_nu(samplekey)); } else if (!name.compare("ANL_CC2pi_1pip1pi0_Evt_1Dpprot_nu")) { return (new ANL_CC2pi_1pip1pi0_Evt_1Dpprot_nu(samplekey)); /* ArgoNeut Samples */ } else #endif #ifndef __NO_ArgoNeuT__ if (!name.compare("ArgoNeuT_CCInc_XSec_1Dpmu_antinu")) { return (new ArgoNeuT_CCInc_XSec_1Dpmu_antinu(samplekey)); } else if (!name.compare("ArgoNeuT_CCInc_XSec_1Dpmu_nu")) { return (new ArgoNeuT_CCInc_XSec_1Dpmu_nu(samplekey)); } else if (!name.compare("ArgoNeuT_CCInc_XSec_1Dthetamu_antinu")) { return (new ArgoNeuT_CCInc_XSec_1Dthetamu_antinu(samplekey)); } else if (!name.compare("ArgoNeuT_CCInc_XSec_1Dthetamu_nu")) { return (new ArgoNeuT_CCInc_XSec_1Dthetamu_nu(samplekey)); /* 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)); /* MINERvA Samples */ } else #endif #ifndef __NO_MINERvA__ if (!name.compare("MINERvA_CCQE_XSec_1DQ2_nu") || !name.compare("MINERvA_CCQE_XSec_1DQ2_nu_20deg") || !name.compare("MINERvA_CCQE_XSec_1DQ2_nu_oldflux") || !name.compare("MINERvA_CCQE_XSec_1DQ2_nu_20deg_oldflux")) { return (new MINERvA_CCQE_XSec_1DQ2_nu(samplekey)); } else if (!name.compare("MINERvA_CCQE_XSec_1DQ2_antinu") || !name.compare("MINERvA_CCQE_XSec_1DQ2_antinu_20deg") || !name.compare("MINERvA_CCQE_XSec_1DQ2_antinu_oldflux") || !name.compare("MINERvA_CCQE_XSec_1DQ2_antinu_20deg_oldflux")) { return (new MINERvA_CCQE_XSec_1DQ2_antinu(samplekey)); } else if (!name.compare("MINERvA_CCQE_XSec_1DQ2_joint_oldflux") || !name.compare("MINERvA_CCQE_XSec_1DQ2_joint_20deg_oldflux") || !name.compare("MINERvA_CCQE_XSec_1DQ2_joint") || !name.compare("MINERvA_CCQE_XSec_1DQ2_joint_20deg")) { return (new MINERvA_CCQE_XSec_1DQ2_joint(samplekey)); } else if (!name.compare("MINERvA_CC0pi_XSec_1DEe_nue")) { return (new MINERvA_CC0pi_XSec_1DEe_nue(samplekey)); } else if (!name.compare("MINERvA_CC0pi_XSec_1DQ2_nue")) { return (new MINERvA_CC0pi_XSec_1DQ2_nue(samplekey)); } else if (!name.compare("MINERvA_CC0pi_XSec_1DThetae_nue")) { return (new MINERvA_CC0pi_XSec_1DThetae_nue(samplekey)); } else if (!name.compare("MINERvA_CC0pinp_STV_XSec_1Dpmu_nu") || !name.compare("MINERvA_CC0pinp_STV_XSec_1Dthmu_nu") || !name.compare("MINERvA_CC0pinp_STV_XSec_1Dpprot_nu") || !name.compare("MINERvA_CC0pinp_STV_XSec_1Dthprot_nu") || !name.compare("MINERvA_CC0pinp_STV_XSec_1Dpnreco_nu") || !name.compare("MINERvA_CC0pinp_STV_XSec_1Ddalphat_nu") || !name.compare("MINERvA_CC0pinp_STV_XSec_1Ddpt_nu") || !name.compare("MINERvA_CC0pinp_STV_XSec_1Ddphit_nu")) { return (new MINERvA_CC0pinp_STV_XSec_1D_nu(samplekey)); } else if (!name.compare("MINERvA_CC0pi_XSec_1DQ2_nu_proton")) { return (new MINERvA_CC0pi_XSec_1DQ2_nu_proton(samplekey)); } else if (!name.compare("MINERvA_CC0pi_XSec_1DQ2_TgtC_nu") || !name.compare("MINERvA_CC0pi_XSec_1DQ2_TgtCH_nu") || !name.compare("MINERvA_CC0pi_XSec_1DQ2_TgtFe_nu") || !name.compare("MINERvA_CC0pi_XSec_1DQ2_TgtPb_nu")) { return (new MINERvA_CC0pi_XSec_1DQ2_Tgt_nu(samplekey)); } else if (!name.compare("MINERvA_CC0pi_XSec_1DQ2_TgtRatioC_nu") || !name.compare("MINERvA_CC0pi_XSec_1DQ2_TgtRatioFe_nu") || !name.compare("MINERvA_CC0pi_XSec_1DQ2_TgtRatioPb_nu")) { return (new MINERvA_CC0pi_XSec_1DQ2_TgtRatio_nu(samplekey)); - } else if (!name.compare("MINERvA_CC0pi_XSec_2Dptpx_nu")) { - return (new MINERvA_CC0pi_XSec_2Dptpx_nu(samplekey)); + } else if ( !name.compare("MINERvA_CC0pi_XSec_2Dptpz_nu") || + !name.compare("MINERvA_CC0pi_XSec_2DptQ2_nu")) { + return (new MINERvA_CC0pi_XSec_2D_nu(samplekey)); - } else if (!name.compare("MINERvA_CC0pi_XSec_2Dptpx_antinu")) { - return (new MINERvA_CC0pi_XSec_2Dptpx_antinu(samplekey)); + //} else if (!name.compare("MINERvA_CC0pi_XSec_2Dptpx_antinu")) { + //return (new MINERvA_CC0pi_XSec_2Dptpx_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)); // Done } else if (!name.compare("MINERvA_CCNpip_XSec_1Dthmu_nu")) { return (new MINERvA_CCNpip_XSec_1Dthmu_nu(samplekey)); // Done } else if (!name.compare("MINERvA_CCNpip_XSec_1Dpmu_nu")) { return (new MINERvA_CCNpip_XSec_1Dpmu_nu(samplekey)); // Done } else if (!name.compare("MINERvA_CCNpip_XSec_1DQ2_nu")) { return (new MINERvA_CCNpip_XSec_1DQ2_nu(samplekey)); // Done } 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)); } else if (!name.compare("T2K_CC0pi_XSec_2DPcos_nu_nonuniform")) { return (new T2K_CC0pi_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_1Dpmu_nu")) { return (new T2K_CC1pip_CH_XSec_1Dpmu_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_1DQ2_nu")) { return (new T2K_CC1pip_CH_XSec_1DQ2_nu(file, rw, type, fkdt)); } else if (!name.compare("T2K_CC1pip_CH_XSec_1Dq3_nu")) { return (new T2K_CC1pip_CH_XSec_1Dq3_nu(file, rw, type, fkdt)); } else if (!name.compare("T2K_CC1pip_CH_XSec_1Dthmupi_nu")) { return (new T2K_CC1pip_CH_XSec_1Dthmupi_nu(file, rw, type, fkdt)); } else if (!name.compare("T2K_CC1pip_CH_XSec_1Dthpi_nu")) { return (new T2K_CC1pip_CH_XSec_1Dthpi_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)); } else if (!name.compare("T2K_CC1pip_CH_XSec_1DWrec_nu")) { 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 { ERR(FTL) << "Error: No such sample: " << name << std::endl; exit(-1); return NULL; } // Return NULL if no sample loaded. return NULL; } } diff --git a/src/MINERvA/CMakeLists.txt b/src/MINERvA/CMakeLists.txt index 95fabad..5bd5049 100644 --- a/src/MINERvA/CMakeLists.txt +++ b/src/MINERvA/CMakeLists.txt @@ -1,178 +1,176 @@ # 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 MINERvA_CCQE_XSec_1DQ2_antinu.cxx MINERvA_CCQE_XSec_1DQ2_joint.cxx MINERvA_CCQE_XSec_1DQ2_nu.cxx MINERvA_CC0pi_XSec_1DEe_nue.cxx MINERvA_CC0pi_XSec_1DQ2_nue.cxx MINERvA_CC0pi_XSec_1DQ2_nu_proton.cxx MINERvA_CC0pi_XSec_1DThetae_nue.cxx MINERvA_CC0pinp_STV_XSec_1D_nu.cxx MINERvA_CC1pi0_XSec_1DEnu_antinu.cxx MINERvA_CC1pi0_XSec_1DQ2_antinu.cxx MINERvA_CC1pi0_XSec_1Dpmu_antinu.cxx MINERvA_CC1pi0_XSec_1Dppi0_antinu.cxx MINERvA_CC1pi0_XSec_1DTpi0_antinu.cxx MINERvA_CC1pi0_XSec_1Dth_antinu.cxx MINERvA_CC1pi0_XSec_1Dthmu_antinu.cxx MINERvA_CC1pi0_XSec_1D_nu.cxx MINERvA_CC1pip_XSec_1DTpi_20deg_nu.cxx MINERvA_CC1pip_XSec_1DTpi_nu.cxx MINERvA_CC1pip_XSec_1Dth_20deg_nu.cxx MINERvA_CC1pip_XSec_1Dth_nu.cxx MINERvA_CC1pip_XSec_1D_2017Update.cxx MINERvA_CCNpip_XSec_1DEnu_nu.cxx MINERvA_CCNpip_XSec_1DQ2_nu.cxx MINERvA_CCNpip_XSec_1DTpi_nu.cxx MINERvA_CCNpip_XSec_1Dpmu_nu.cxx MINERvA_CCNpip_XSec_1Dth_nu.cxx MINERvA_CCNpip_XSec_1Dthmu_nu.cxx MINERvA_CCinc_XSec_2DEavq3_nu.cxx MINERvA_CCinc_XSec_1Dx_ratio.cxx MINERvA_CCinc_XSec_1DEnu_ratio.cxx MINERvA_CCinc_XSec_1Dx_nu.cxx MINERvA_CCinc_XSec_1DEnu_nu.cxx MINERvA_CCDIS_XSec_1Dx_ratio.cxx MINERvA_CCDIS_XSec_1DEnu_ratio.cxx MINERvA_CCDIS_XSec_1Dx_nu.cxx MINERvA_CCDIS_XSec_1DEnu_nu.cxx MINERvA_CC0pi_XSec_1DQ2_Tgt_nu.cxx MINERvA_CC0pi_XSec_1DQ2_TgtRatio_nu.cxx -MINERvA_CC0pi_XSec_2Dptpx_nu.cxx -MINERvA_CC0pi_XSec_2Dptpx_antinu.cxx +MINERvA_CC0pi_XSec_2D_nu.cxx MINERvA_CCCOHPI_XSec_1DEnu_nu.cxx MINERvA_CCCOHPI_XSec_1DEpi_nu.cxx MINERvA_CCCOHPI_XSec_1Dth_nu.cxx MINERvA_CCCOHPI_XSec_1DQ2_nu.cxx MINERvA_CCCOHPI_XSec_1DEnu_antinu.cxx MINERvA_CCCOHPI_XSec_1DEpi_antinu.cxx MINERvA_CCCOHPI_XSec_1Dth_antinu.cxx MINERvA_CCCOHPI_XSec_1DQ2_antinu.cxx MINERvA_CCCOHPI_XSec_joint.cxx MINERvAUtils.cxx MINERvA_SignalDef.cxx ) set(HEADERFILES MINERvA_CCQE_XSec_1DQ2_antinu.h MINERvA_CCQE_XSec_1DQ2_joint.h MINERvA_CCQE_XSec_1DQ2_nu.h MINERvA_CC0pi_XSec_1DEe_nue.h MINERvA_CC0pi_XSec_1DQ2_nue.h MINERvA_CC0pi_XSec_1DQ2_nu_proton.h MINERvA_CC0pi_XSec_1DThetae_nue.h MINERvA_CC0pinp_STV_XSec_1D_nu.h MINERvA_CC1pi0_XSec_1DEnu_antinu.h MINERvA_CC1pi0_XSec_1DQ2_antinu.h MINERvA_CC1pi0_XSec_1Dpmu_antinu.h MINERvA_CC1pi0_XSec_1Dppi0_antinu.h MINERvA_CC1pi0_XSec_1DTpi0_antinu.h MINERvA_CC1pi0_XSec_1Dth_antinu.h MINERvA_CC1pi0_XSec_1Dthmu_antinu.h MINERvA_CC1pip_XSec_1DTpi_20deg_nu.h MINERvA_CC1pip_XSec_1DTpi_nu.h MINERvA_CC1pip_XSec_1Dth_20deg_nu.h MINERvA_CC1pip_XSec_1Dth_nu.h MINERvA_CCNpip_XSec_1DEnu_nu.h MINERvA_CCNpip_XSec_1DQ2_nu.h MINERvA_CCNpip_XSec_1DTpi_nu.h MINERvA_CCNpip_XSec_1Dpmu_nu.h MINERvA_CCNpip_XSec_1Dth_nu.h MINERvA_CCNpip_XSec_1Dthmu_nu.h MINERvA_CCinc_XSec_2DEavq3_nu.h MINERvA_CCinc_XSec_1Dx_ratio.h MINERvA_CCinc_XSec_1DEnu_ratio.h MINERvA_CCinc_XSec_1Dx_nu.h MINERvA_CCinc_XSec_1DEnu_nu.h MINERvA_CCDIS_XSec_1Dx_ratio.h MINERvA_CCDIS_XSec_1DEnu_ratio.h MINERvA_CCDIS_XSec_1Dx_nu.h MINERvA_CCDIS_XSec_1DEnu_nu.h MINERvA_CC0pi_XSec_1DQ2_Tgt_nu.h MINERvA_CC0pi_XSec_1DQ2_TgtRatio_nu.h -MINERvA_CC0pi_XSec_2Dptpx_nu.h -MINERvA_CC0pi_XSec_2Dptpx_antinu.h +MINERvA_CC0pi_XSec_2D_nu.h MINERvA_CC1pip_XSec_1D_2017Update.h MINERvA_CCCOHPI_XSec_1DEnu_nu.h MINERvA_CCCOHPI_XSec_1DEpi_nu.h MINERvA_CCCOHPI_XSec_1Dth_nu.h MINERvA_CCCOHPI_XSec_1DQ2_nu.h MINERvA_CCCOHPI_XSec_1DEnu_antinu.h MINERvA_CCCOHPI_XSec_1DEpi_antinu.h MINERvA_CCCOHPI_XSec_1Dth_antinu.h MINERvA_CCCOHPI_XSec_1DQ2_antinu.h MINERvA_CCCOHPI_XSec_joint.h MINERvAUtils.h MINERvA_SignalDef.h MINERvAVariableBoxes.h ) set(LIBNAME expMINERvA) 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/MINERvA/MINERvA_CC0pi_XSec_2D_nu.cxx b/src/MINERvA/MINERvA_CC0pi_XSec_2D_nu.cxx new file mode 100755 index 0000000..f568571 --- /dev/null +++ b/src/MINERvA/MINERvA_CC0pi_XSec_2D_nu.cxx @@ -0,0 +1,464 @@ +//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 . +*******************************************************************************/ + +/* + Authors: Adrian Orea (v1 2017) + Clarence Wret (v2 2018) +*/ + +#include "MINERvA_SignalDef.h" +#include "MINERvA_CC0pi_XSec_2D_nu.h" + +//******************************************************************** +void MINERvA_CC0pi_XSec_2D_nu::SetupDataSettings() { +//******************************************************************** + // Set Distribution + // See header file for enum and some descriptions + std::string name = fSettings.GetS("name"); + // Have two distributions as of summer 2018 + if (!name.compare("MINERvA_CC0pi_XSec_2Dptpz_nu")) fDist = kPtPz; + else if (!name.compare("MINERvA_CC0pi_XSec_2DptQ2_nu")) fDist = kPtQ2; + + // Define what files to use from the dist + std::string datafile = ""; + std::string corrfile = ""; + std::string titles = ""; + std::string distdescript = ""; + std::string histname = ""; + + switch (fDist) { + case (kPtPz): + datafile = "MINERvA/CC0pi_2D/MINERvA_LE_CCQELike_pzmu.root"; + corrfile = "MINERvA/CC0pi_2D/MINERvA_LE_CCQELike_pzmu.root"; + titles = "MINERvA CC0#pi #nu_{#mu} p_{t} p_{z};p_{z} (GeV);p_{t} (GeV);d^{2}#sigma/dP_{t}dP_{z} (cm^{2}/GeV^{2}/nucleon)"; + distdescript = "MINERvA_CC0pi_XSec_2Dptpz_nu sample"; + histname = "h_pzmu_ptmu_data_nobck_unfold_effcor_cross_section_CV_WithErr"; + break; + case (kPtQ2): + datafile = "MINERvA/CC0pi_2D/pt_q2_data.root"; + corrfile = "MINERvA/CC0pi_2D/pt_q2_cov.root"; + titles = "MINERvA CC0#pi #nu_{#mu} p_{t} Q^{2}_{QE};p_t{z} (GeV);Q^{2}_{QE} (GeV^{2});d^{2}#sigma/dP_{t}dQ^{2}_{QE} (cm^{2}/GeV^{3}/nucleon)"; + distdescript = "MINERvA_CC0pi_XSec_2DptQ2_nu sample"; + histname = "h_q2_ptmu_data_nobck_unfold_effcor_cross_section_CV_WithErr"; + break; + default: + THROW("Unknown Analysis Distribution : " << fDist); + } + + fSettings.SetTitle( GeneralUtils::ParseToStr(titles,";")[0] ); + fSettings.SetXTitle( GeneralUtils::ParseToStr(titles,";")[1] ); + fSettings.SetYTitle( GeneralUtils::ParseToStr(titles,";")[2] ); + fSettings.SetYTitle( GeneralUtils::ParseToStr(titles,";")[3] ); + + // Sample overview --------------------------------------------------- + std::string descrip = distdescript + "\n"\ + "Target: CH \n" \ + "Flux: MINERvA Low Energy FHC numu \n" \ + "Signal: CC-0pi \n"; + fSettings.SetDescription(descrip); + + // The input ROOT file + fSettings.SetDataInput( FitPar::GetDataBase() + datafile); + // The covariance matrix ROOT file + //if (fDist == kPtPz) { + fSettings.SetCovarInput( FitPar::GetDataBase() + corrfile); + //} else { + //ERR(WRN) << " no covariance matrix available for PtQ2" << std::endl; + //ERR(WRN) << " ask Dan Ruterbories nicely and he may provide one" << std::endl; + //} + + // Set the data file + SetDataValues(fSettings.GetDataInput(), histname); +} +//******************************************************************** +MINERvA_CC0pi_XSec_2D_nu::MINERvA_CC0pi_XSec_2D_nu(nuiskey samplekey) { +//******************************************************************** + + // Setup common settings + fSettings = LoadSampleSettings(samplekey); + fSettings.SetAllowedTypes("FIX,FREE,SHAPE/FULL,DIAG/MASK", "FIX/FULL"); + fSettings.SetEnuRange(0.0, 100.0); + fSettings.DefineAllowedTargets("C,H"); + fSettings.DefineAllowedSpecies("numu"); + + SetupDataSettings(); + FinaliseSampleSettings(); + + fScaleFactor = (GetEventHistogram()->Integral("width") * 1E-38 / (fNEvents + 0.)) / this->TotalIntegratedFlux(); + + // Data is __NOT__ bin width normalised, so override the default + fIsNoWidth = true; + + // Set the mapping values and the covariance matrix files + //SetMapValuesFromText( fSettings.GetMapInput() ); + // Also have to make our own covariance matrix to exclude the under and overflow + //if (fDist == kPtPz) { + //TMatrixDSym* tempmat = StatUtils::GetCovarFromRootFile(fSettings.GetCovarInput(), "TMatrixDBase"); + TMatrixDSym* tempmat = StatUtils::GetCovarFromRootFile(fSettings.GetCovarInput(), "CovMatrix"); + fFullCovar = tempmat; + // Decomposition is stable for entries that aren't E-xxx + double ScalingFactor = 1E38; + (*fFullCovar) *= ScalingFactor; + + //} else { + //SetCovarFromDiagonal(); + //} + + /* + // Now we cut out every first and last to exclude under and overflow + fFullCovar = new TMatrixDSym(fDataHist->GetXaxis()->GetNbins()*fDataHist->GetYaxis()->GetNbins()); + + // Count the real covariance matrix x and y + int xcounter = 0; + + int xbins = fDataHist->GetXaxis()->GetNbins(); + int ybins = fDataHist->GetYaxis()->GetNbins(); + // Loop over the x bins (underflow adds one, overflow adds one) + for (int i = 0; i < (xbins+2)*(ybins+2); ++i) { + // Skip the under and overflow + if (i < (ybins+2) || i % (ybins+1) == 0 || ((i+1)%(ybins+1) == 0) || i > (ybins+2)*(xbins+1)) { + //std::cout << "Skipping ibin " << i << std::endl; + continue; + } + + // The ycounter + int ycounter = 0; + // For one bin of pT we have pZ^2 bins + for (int j = 0; j < (xbins+2)*(ybins+2); ++j) { + + // Skip the under and overflow + if (j < (ybins+2) || j % (ybins+1) == 0 || ((j+1)%(ybins+1) == 0) || j > (ybins+2)*(xbins+1)) { + //std::cout << "Skipping jbin " << j << std::endl; + continue; + } + + (*fFullCovar)(xcounter, ycounter) = (*tempmat)(i, j); + //std::cout << xcounter << ", " << ycounter << " === " << i << ", " << j << std::endl; + ycounter++; + } + // Check dimensions + if (ycounter != xbins*ybins) { + std::cerr << "Counted " << ycounter << " y bins in cov matrix" << std::endl; + std::cerr << "Whereas there should be " << xbins*ybins << std::endl; + } + xcounter++; + } + // Check dimensions + if (xcounter != xbins*ybins) { + std::cerr << "Counted " << xcounter << " x bins in cov matrix" << std::endl; + std::cerr << "Whereas there should be " << xbins*ybins << std::endl; + } + // Delete the temporary + delete tempmat; + */ + + // Just check that the data error and covariance are the same + for (int i = 0; i < fFullCovar->GetNrows(); ++i) { + for (int j = 0; j < fFullCovar->GetNcols(); ++j) { + // Get the global bin + int xbin1, ybin1, zbin1; + fDataHist->GetBinXYZ(i, xbin1, ybin1, zbin1); + double xlo1 = fDataHist->GetXaxis()->GetBinLowEdge(xbin1); + double xhi1 = fDataHist->GetXaxis()->GetBinLowEdge(xbin1+1); + double ylo1 = fDataHist->GetYaxis()->GetBinLowEdge(ybin1); + double yhi1 = fDataHist->GetYaxis()->GetBinLowEdge(ybin1+1); + if (xlo1 < fDataHist->GetXaxis()->GetBinLowEdge(1) || + ylo1 < fDataHist->GetYaxis()->GetBinLowEdge(1) || + xhi1 > fDataHist->GetXaxis()->GetBinLowEdge(fDataHist->GetXaxis()->GetNbins()+1) || + yhi1 > fDataHist->GetYaxis()->GetBinLowEdge(fDataHist->GetYaxis()->GetNbins()+1)) continue; + double data_error = fDataHist->GetBinError(xbin1, ybin1); + double cov_error = sqrt((*fFullCovar)(i,i)/ScalingFactor); + if (fabs(data_error - cov_error) > 1E-5) { + std::cerr << "Error on data is different to that of covariance" << std::endl; + ERR(FTL) << "Data error: " << data_error << std::endl; + ERR(FTL) << "Cov error: " << cov_error << std::endl; + ERR(FTL) << "Data/Cov: " << data_error/cov_error << std::endl; + ERR(FTL) << "Data-Cov: " << data_error-cov_error << std::endl; + ERR(FTL) << "For x: " << xlo1 << "-" << xhi1 << std::endl; + ERR(FTL) << "For y: " << ylo1 << "-" << yhi1 << std::endl; + throw; + } + } + } + + // Now can make the inverted covariance + covar = StatUtils::GetInvert(fFullCovar); + fDecomp = StatUtils::GetDecomp(fFullCovar); + + // Now scale back + (*fFullCovar) *= 1.0/ScalingFactor; + (*covar) *= ScalingFactor; + (*fDecomp) *= ScalingFactor; + + // Use a TH2D version of the covariance to be able to use the global bin numbering scheme + covar_th2d = new TH2D((fSettings.Title()+"_th2").c_str(), (fSettings.Title()+"_th2").c_str(), covar->GetNrows(), 0, covar->GetNrows(), covar->GetNcols(), 0, covar->GetNcols()); + for (int i = 0; i < covar_th2d->GetXaxis()->GetNbins(); ++i) { + for (int j = 0; j < covar_th2d->GetYaxis()->GetNbins(); ++j) { + covar_th2d->SetBinContent(i+1, j+1, (*covar)(i,j)); + } + } + //std::cout << "covar is " << covar_th2d->GetXaxis()->GetNbins() << " x " << covar_th2d->GetYaxis()->GetNbins() << " = " << covar_th2d->GetXaxis()->GetNbins()*covar_th2d->GetYaxis()->GetNbins() << std::endl; + //std::cout << "data is " << fDataHist->GetXaxis()->GetNbins() << " x " << fDataHist->GetYaxis()->GetNbins() << " = " << fDataHist->GetXaxis()->GetNbins()*fDataHist->GetYaxis()->GetNbins() << std::endl; + + // Let's make our own mapping histogram + // The covariance matrix is dominant in Pt and sub-dominant in Pz and includes all bins with under and overflow + // So we have 13x12 data/MC bins, and including the under/overflow we have 15x14=210 covariance bins + // i.e. need to cut out the first and last bins of covariance matrix + // Mapping histogram will have same dimensions as the data + /* + fMapHist = (TH2I*)(fDataHist->Clone()); + fMapHist->Reset(); + std::string MapTitle = std::string(fDataHist->GetName())+"_MAP"; + fMapHist->SetNameTitle(MapTitle.c_str(), MapTitle.c_str()); + int counter = 1; + for (int i = 0; i <= fDataHist->GetXaxis()->GetNbins()+1; ++i) { + for (int j = 0; j <= fDataHist->GetYaxis()->GetNbins()+1; ++j) { + if (i == 0 || i == fDataHist->GetXaxis()->GetNbins()+1 || j == 0 || j == fDataHist->GetYaxis()->GetNbins()+1) { + fMapHist->SetBinContent(i+1, j+1, 0); + } else { + fMapHist->SetBinContent(i+1, j+1, counter); + counter++; + } + std::cout << fMapHist->GetBinContent(i+1, j+1) << " " << fDataHist->GetBinContent(i+1, j+1) << std::endl; + } + std::cout << std::endl; + } + */ + + // Final setup --------------------------------------------------- + FinaliseMeasurement(); +}; + +//******************************************************************** +void MINERvA_CC0pi_XSec_2D_nu::FillEventVariables(FitEvent *event) { + //******************************************************************** + // Checking to see if there is a Muon + if (event->NumFSParticle(13) == 0) return; + + // Get the muon kinematics + TLorentzVector Pmu = event->GetHMFSParticle(13)->fP; + TLorentzVector Pnu = event->GetNeutrinoIn()->fP; + + Double_t px = Pmu.X()/1000; + Double_t py = Pmu.Y()/1000; + Double_t pt = sqrt(px*px+py*py); + + // y-axis is transverse momentum for both measurements + fYVar = pt; + + // Now we set the x-axis + switch (fDist) { + case kPtPz: + { + // Don't want to assume the event generators all have neutrino coming along z + // pz is muon momentum projected onto the neutrino direction + Double_t pz = Pmu.Vect().Dot(Pnu.Vect()*(1.0/Pnu.Vect().Mag()))/1000.; + // Set Hist Variables + fXVar = pz; + break; + } + case kPtQ2: + { + // 34 MeV binding energy and neutrino mode (true) + double q2qe = FitUtils::Q2QErec(Pmu, Pnu, 34.0, true); + fXVar = q2qe; + break; + } + default: + THROW("DIST NOT FOUND : " << fDist); + break; + } + +}; + +//******************************************************************** +bool MINERvA_CC0pi_XSec_2D_nu::isSignal(FitEvent *event) { + //******************************************************************** + return SignalDef::isCC0pi_MINERvAPTPZ(event, 14, EnuMin, EnuMax); +}; + +//******************************************************************** +// Custom likelihood calculator because binning of covariance matrix +double MINERvA_CC0pi_XSec_2D_nu::GetLikelihood() { + //******************************************************************** + + // The calculated chi2 + double chi2 = 0.0; + + // STRICTLY DEBUGGING WITH DAN + /* + std::vector Names; + //Names.push_back(std::string(std::getenv("NUISANCE"))+"/build/app/DanMC/Genie_MC.root"); + //Names.push_back(std::string(std::getenv("NUISANCE"))+"/build/app/DanMC/Genie2p2hrpa_MC.root"); + //Names.push_back(std::string(std::getenv("NUISANCE"))+"/build/app/DanMC/MnvGENIE_MC.root"); + Names.push_back(std::string(std::getenv("NUISANCE"))+"/data/MINERvA/CC0pi_2D/MINERvA_LE_CCQELike_pzmu.root"); + + // Hack hack hack hack + double scaleF = 0.0; + + for (size_t a = 0; a < Names.size(); ++a) { + std::cout << Names[a] << std::endl; + //TFile *file = new TFile((std::string(std::getenv("NUISANCE"))+"/build/app/DanMC/Genie_MC.root").c_str(), "OPEN"); + //TFile *file = new TFile((std::string(std::getenv("NUISANCE"))+"/build/app/DanMC/Genie2p2hrpa_MC.root").c_str(), "OPEN"); + //TFile *file = new TFile((std::string(std::getenv("NUISANCE"))+"/build/app/DanMC/MnvGENIE_MC.root").c_str(), "OPEN"); + TFile *file = new TFile(Names[a].c_str(), "OPEN"); + //TH2D *mc = (TH2D*)(file->Get("h_pzmu_ptmu_mc_nobck_unfold_effcor_cross_section_CV_WithStatErr")->Clone()); + TH2D *mc = (TH2D*)(file->Get("h_pzmu_ptmu_mc_nobck_unfold_effcor_cross_section_CV_WithErr")->Clone()); + fMCHist = mc; + //fMCHist->Scale(1., "width"); + //fDataHist->Scale(1., "width"); + */ + + // Support shape comparisons + double scaleF = fDataHist->Integral() / fMCHist->Integral(); + if (fIsShape) { + fMCHist->Scale(scaleF); + fMCFine->Scale(scaleF); + //PlotUtils::ScaleNeutModeArray((TH1**)fMCHist_PDG, scaleF); + } + + // Even though this chi2 calculation looks ugly it is _EXACTLY_ what MINERvA used for their measurement + // Can be prettified in due time but for now keep + bool chi2_use_overflow_err = false; + const int lowBin = chi2_use_overflow_err?0:1; // Either they both use underflow, or neither of them does + + int nbinsx=fMCHist->GetNbinsX(); + int nbinsy=fMCHist->GetNbinsY(); + const int highBinX = nbinsx; + const int highBinY = nbinsy; + Int_t Nbins = (highBinX-lowBin+1)* (highBinY-lowBin+1); // from low to high inclusive in each dimension + + TMatrixD covMatrixTmp = (*fFullCovar); + + TMatrixD covMatrix(Nbins, Nbins); + const double scaling = 1e80;//ROOT can't seem to handle small entries with lots of zeros? Suggested scaling the histogram and then rescaling the inverted matrix + for( int i = 0; i != Nbins; ++i ) { + for( int j = 0; j != Nbins; ++j ) { + // I think this is right... it checks out on a small sample + int old_i_bin = (i/nbinsx + 1)* (nbinsx +2) +1 + i%nbinsx; + int old_j_bin = (j/nbinsx + 1)* (nbinsx +2) +1 + j%nbinsx; + covMatrix[i][j] = covMatrixTmp[old_i_bin][old_j_bin]*scaling; + //std::cout << Nbins*Nbins << " (" << Nbins << "*" << Nbins << ")" << std::endl; + //std::cout << i << ", " << j << " = " << old_i_bin << " " << old_j_bin << std::endl; + } + //throw; + } + TDecompSVD error(covMatrix); + TMatrixD errorMatrix(covMatrix); + if( ! error.Invert( errorMatrix ) ) { + std::cout << "Cannot invert total covariance matrix. You could use statistical errors only for Chi2 calculation. But it isn't implemented yet" << std::endl; + } + + covMatrix *= 1/scaling; + errorMatrix *= scaling; + + for( int i = 0; i != Nbins; ++i ) { + int hist_i_bin = chi2_use_overflow_err?i:((i/nbinsx + 1)* (nbinsx +2) +1 + i%nbinsx); // Translate to the histogram bin, if we aren't using the overflow errors, meaning the covariance matrix is smaller than the histogram + const Double_t x_data_i = fDataHist->GetBinContent(hist_i_bin); + const Double_t x_mc_i = fMCHist->GetBinContent(hist_i_bin); + for( int j = 0; j != Nbins; ++j ) { + // Each element of the inverted covariance matrix corresponds to a pair of data and MC + int hist_j_bin = chi2_use_overflow_err?j:((j/nbinsx + 1)* (nbinsx +2) +1 + j%nbinsx); + const Double_t x_data_j = fDataHist->GetBinContent(hist_j_bin); + const Double_t x_mc_j = fMCHist->GetBinContent(hist_j_bin); + const double chi2_ij = (x_data_i - x_mc_i) * errorMatrix[i][j] * (x_data_j - x_mc_j); + std::cout << x_data_i << "\t" << x_mc_i << "\t" << i << "\t" << errorMatrix[i][j] << "\t" << j << "\t" << x_data_j << "\t" << x_mc_j << "\t" << chi2_ij << std::endl; + + chi2 += chi2_ij; + } + } + //std::cout << Names[a] << " chi2 " << chi2 << std::endl; + + /* + * CWRET CALC + // Calculate the test-statistic + for (int i = 0; i < covar_th2d->GetXaxis()->GetNbins()+1; ++i) { + // Get the global bin for x + int xbin1, ybin1, zbin1; + fDataHist->GetBinXYZ(i, xbin1, ybin1, zbin1); + double xlo1 = fDataHist->GetXaxis()->GetBinLowEdge(xbin1); + double xhi1 = fDataHist->GetXaxis()->GetBinLowEdge(xbin1+1); + double ylo1 = fDataHist->GetYaxis()->GetBinLowEdge(ybin1); + double yhi1 = fDataHist->GetYaxis()->GetBinLowEdge(ybin1+1); + if (xlo1 < fDataHist->GetXaxis()->GetBinLowEdge(1) || + ylo1 < fDataHist->GetYaxis()->GetBinLowEdge(1) || + xhi1 > fDataHist->GetXaxis()->GetBinLowEdge(fDataHist->GetXaxis()->GetNbins()+1) || + yhi1 > fDataHist->GetYaxis()->GetBinLowEdge(fDataHist->GetYaxis()->GetNbins()+1)) continue; + + // Get the data + double data1 = fDataHist->GetBinContent(i); + // Get the MC + double mc1 = fMCHist->GetBinContent(i); + + for (int j = 0; j < covar_th2d->GetYaxis()->GetNbins()+1; ++j) { + + int xbin2, ybin2, zbin2; + fDataHist->GetBinXYZ(j, xbin2, ybin2, zbin2); + double xlo2 = fDataHist->GetXaxis()->GetBinLowEdge(xbin2); + double xhi2 = fDataHist->GetXaxis()->GetBinLowEdge(xbin2+1); + double ylo2 = fDataHist->GetYaxis()->GetBinLowEdge(ybin2); + double yhi2 = fDataHist->GetYaxis()->GetBinLowEdge(ybin2+1); + + if (xlo2 < fDataHist->GetXaxis()->GetBinLowEdge(1) || + ylo2 < fDataHist->GetYaxis()->GetBinLowEdge(1) || + xhi2 > fDataHist->GetXaxis()->GetBinLowEdge(fDataHist->GetXaxis()->GetNbins()+1) || + yhi2 > fDataHist->GetYaxis()->GetBinLowEdge(fDataHist->GetYaxis()->GetNbins()+1)) continue; + + + //std::cout << "Correlating: (" << xlo1 << "-" << xhi1 << "," << ylo1 << "-" << yhi1 << ") with (" << xlo2 << "-" << xhi2 << "," << ylo2 << "-" << yhi2 << ")" << std::endl; + + // Get the data + double data2 = fDataHist->GetBinContent(j); + // Get the MC + double mc2 = fMCHist->GetBinContent(j); + //std::cout << data1 << " " << mc1 << std::endl; + //std::cout << data2 << " " << mc2 << std::endl; + //std::cout << std::endl; + + // Get the inverse covariance matrix entry + double coventry = covar_th2d->GetBinContent(i, j); + + //std::cout << fDataHist->GetXaxis()->GetBinLowEdge(i+1) << " - " << fDataHist->GetXaxis()->GetBinLowEdge(i+2) << ", " << fDataHist->GetYaxis()->GetBinLowEdge(j+1) << " - " << fDataHist->GetYaxis()->GetBinLowEdge(j+2) << " = " << coventry << " (global = " << global << ")" << std::endl; + + chi2 += (data1-mc1)*coventry*(data2-mc2); + } + } + */ + //} + + // Normalisation penalty term if included + if (fAddNormPen) { + chi2 += + (1 - (fCurrentNorm)) * (1 - (fCurrentNorm)) / (fNormError * fNormError); + LOG(REC) << "Norm penalty = " + << (1 - (fCurrentNorm)) * (1 - (fCurrentNorm)) / + (fNormError * fNormError) + << std::endl; + } + + // Adjust the shape back to where it was. + if (fIsShape and !FitPar::Config().GetParB("saveshapescaling")) { + fMCHist->Scale(1. / scaleF); + fMCFine->Scale(1. / scaleF); + } + + fLikelihood = chi2; + + return chi2; +}; diff --git a/src/MINERvA/MINERvA_CC0pi_XSec_2D_nu.h b/src/MINERvA/MINERvA_CC0pi_XSec_2D_nu.h new file mode 100755 index 0000000..15cccce --- /dev/null +++ b/src/MINERvA/MINERvA_CC0pi_XSec_2D_nu.h @@ -0,0 +1,80 @@ +// 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 MINERVA_CC0PI_XSEC_2D_NU_H_SEEN +#define MINERVA_CC0PI_XSEC_2D_NU_H_SEEN + +#include "Measurement2D.h" + +//******************************************************************** +class MINERvA_CC0pi_XSec_2D_nu : public Measurement2D { +//******************************************************************** + + public: + + // Constructor + MINERvA_CC0pi_XSec_2D_nu(nuiskey samplekey); + + // Destructor + virtual ~MINERvA_CC0pi_XSec_2D_nu() { + + // Remove all the content histograms * + // for (int i = 0; i < 9; i++) + // delete this->fMCHist_content[i]; + + // delete everything + /* delete difHist; */ + /* delete evtsignalHist; */ + /* delete fluxsignalHist; */ + /* delete fMapHist; */ + /* delete status; */ + /* delete PDGHistogram; */ + + /* // Delete MODE histograms */ + /* for (int i = 0; i < 60; i++) */ + /* delete fMCHist_PDG[i]; */ + + return; + }; + + // Required functions + bool isSignal(FitEvent *nvect); + void FillEventVariables(FitEvent *event); + + protected: + // Converted covariance matrix to provide global binning method in GetLikelihood + TH2D* covar_th2d; + double GetLikelihood(); + + // Set up settings based on distribution + void SetupDataSettings(); + + private: + // The distribution privates + int fDist; + enum Distribution { + // Pt Pz + kPtPz, + // Pt Q2 + kPtQ2 + }; + +}; + +#endif diff --git a/src/MINERvA/MINERvA_CC0pinp_STV_XSec_1D_nu.cxx b/src/MINERvA/MINERvA_CC0pinp_STV_XSec_1D_nu.cxx index 60468c3..b29e4d2 100644 --- a/src/MINERvA/MINERvA_CC0pinp_STV_XSec_1D_nu.cxx +++ b/src/MINERvA/MINERvA_CC0pinp_STV_XSec_1D_nu.cxx @@ -1,291 +1,327 @@ // 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 . *******************************************************************************/ // Implementation of 2018 MINERvA numu CC0pi STV analysis // arxiv 1805.05486.pdf // Clarence Wret // cwret@fnal.gov // Stephen Dolan // Stephen.Dolan@llr.in2p3.fr #include "MINERvA_SignalDef.h" #include "MINERvA_CC0pinp_STV_XSec_1D_nu.h" //******************************************************************** void MINERvA_CC0pinp_STV_XSec_1D_nu::SetupDataSettings() { //******************************************************************** // Set Distribution // See header file for enum and some descriptions std::string name = fSettings.GetS("name"); if (!name.compare("MINERvA_CC0pinp_STV_XSec_1Dpmu_nu")) fDist= kMuonMom; else if (!name.compare("MINERvA_CC0pinp_STV_XSec_1Dthmu_nu")) fDist= kMuonTh; else if (!name.compare("MINERvA_CC0pinp_STV_XSec_1Dpprot_nu")) fDist= kPrMom; else if (!name.compare("MINERvA_CC0pinp_STV_XSec_1Dthprot_nu")) fDist= kPrTh; else if (!name.compare("MINERvA_CC0pinp_STV_XSec_1Dpnreco_nu")) fDist= kNeMom; else if (!name.compare("MINERvA_CC0pinp_STV_XSec_1Ddalphat_nu")) fDist= kDalphaT; else if (!name.compare("MINERvA_CC0pinp_STV_XSec_1Ddpt_nu")) fDist= kDpT; else if (!name.compare("MINERvA_CC0pinp_STV_XSec_1Ddphit_nu")) fDist= kDphiT; // Location of data, correlation matrices and the titles std::string titles = "MINERvA_CC0pinp_STV_XSec_1D"; std::string foldername; std::string distdescript; // Data release is a single file std::string rootfile = "MINERvA_1805.05486.root"; + fMin = -999; + fMax = 999; + switch (fDist) { case (kMuonMom): titles += "pmu"; foldername = "muonmomentum"; distdescript = "Muon momentum in lab frame"; + fMin = 2.0; + fMax = 6.0; break; case (kMuonTh): titles += "thmu"; foldername = "muontheta"; distdescript = "Muon angle relative neutrino in lab frame"; + fMin = 0.0; + fMax = 20.0; break; case (kPrMom): titles += "pprot"; foldername = "protonmomentum"; distdescript = "Proton momentum in lab frame"; + fMin = 0.5; + fMax = 1.2; break; case (kPrTh): titles += "thprot"; foldername = "protontheta"; distdescript = "Proton angle relative neutrino in lab frame"; + fMin = 0.0; + fMax = 70.0; break; case (kNeMom): titles += "pnreco"; foldername = "neutronmomentum"; distdescript = "Neutron momentum in lab frame"; + fMin = 0.0; + fMax = 0.9; break; case (kDalphaT): foldername = "dalphat"; titles += foldername; distdescript = "Delta Alpha_T"; + fMin = 0.0; + fMax = 170; break; case (kDpT): foldername = "dpt"; titles += foldername; distdescript = "Delta p_T"; + fMin = 0.0; + fMax = 170; break; case (kDphiT): foldername = "dphit"; titles += foldername; distdescript = "Delta phi_T"; + fMin = 0.0; + fMax = 60.0; break; default: ERR(FTL) << "Did not find your specified distribution implemented, exiting" << std::endl; ERR(FTL) << "You gave " << fName << std::endl; throw; } titles += "_nu"; // All have the same name std::string dataname = foldername; // Sample overview --------------------------------------------------- std::string descrip = distdescript + \ "Target: CH \n" \ "Flux: MINERvA Forward Horn Current numu ONLY \n" \ "Signal: Any event with 1 muon, and 0pi0 in FS, no mesons, at least one proton with: \n" \ "1.5GeV < p_mu < 10 GeV\n"\ "theta_mu < 20 degrees\n"\ "0.45GeV < p_prot < 1.2 GeV\n"\ "theta_prot < 70 degrees\n" \ "arXiv 1805.05486"; fSettings.SetDescription(descrip); std::string filename = GeneralUtils::GetTopLevelDir()+"/data/MINERvA/CC0pi/CC0pi_STV/MINERvA_1805.05486.root"; // Specify the data fSettings.SetDataInput(filename+";"+dataname); // And the correlations fSettings.SetCovarInput(filename+";"+dataname); // Set titles fSettings.SetTitle(titles); }; //******************************************************************** MINERvA_CC0pinp_STV_XSec_1D_nu::MINERvA_CC0pinp_STV_XSec_1D_nu(nuiskey samplekey) { //******************************************************************** // A few different distributinos // Muon momentum, muon angle, proton momentum, proton angle, neutron momentum, dalphat, dpt, dphit // Setup common settings fSettings = LoadSampleSettings(samplekey); // Load up the data paths and sample descriptions SetupDataSettings(); fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/DIAG"); // No Enu cut - fSettings.SetEnuRange(0.0, 50.0); + fSettings.SetEnuRange(0.0, 100.0); fSettings.DefineAllowedTargets("C,H"); fSettings.DefineAllowedSpecies("numu"); - // Finalise the settings FinaliseSampleSettings(); // Scaling Setup --------------------------------------------------- // ScaleFactor automatically setup for DiffXSec/cm2/Nucleon fScaleFactor = GetEventHistogram()->Integral("width") * double(1E-38) / (double(fNEvents) * TotalIntegratedFlux("width")); // Set the data and covariance matrix - SetCovarianceFromRootFile(fSettings.GetCovarInput() ); SetDataFromRootFile( fSettings.GetDataInput() ); + SetCovarianceFromRootFile(fSettings.GetCovarInput() ); fSettings.SetXTitle(fDataHist->GetXaxis()->GetTitle()); fSettings.SetYTitle(fDataHist->GetYaxis()->GetTitle()); // Final setup --------------------------------------------------- FinaliseMeasurement(); }; //******************************************************************** // Data comes in a TList // Some bins need stripping out because of zero bin content. Why oh why void MINERvA_CC0pinp_STV_XSec_1D_nu::SetDataFromRootFile(std::string filename) { //******************************************************************** std::vector tempfile = GeneralUtils::ParseToStr(filename, ";"); TFile *File = new TFile(tempfile[0].c_str(), "READ"); // First object is the data TH1D* temp = (TH1D*)(((TList*)(File->Get(tempfile[1].c_str())))->At(0)); // Garh, some of the data points are zero in the TH1D (WHY?!) so messes with the covariance entries to data bins check // Skim through the data and check for zero bins std::vector CrossSection; std::vector Error; std::vector BinEdges; - for (int i = 0; i < temp->GetXaxis()->GetNbins()+1; ++i) { - if (temp->GetBinContent(i+1) > 0) { + int lastbin = 0; + startbin = 0; + for (int i = 0; i < temp->GetXaxis()->GetNbins()+2; ++i) { + if (temp->GetBinContent(i+1) > 0 && temp->GetBinLowEdge(i+1) > fMin && temp->GetBinLowEdge(i+1) < fMax) { + if (startbin == 0) startbin = i; + lastbin = i; CrossSection.push_back(temp->GetBinContent(i+1)); BinEdges.push_back(temp->GetXaxis()->GetBinLowEdge(i+1)); Error.push_back(temp->GetBinError(i+1)); } } + BinEdges.push_back(temp->GetXaxis()->GetBinLowEdge(lastbin+2)); fDataHist = new TH1D((fSettings.GetName()+"_data").c_str(), (fSettings.GetFullTitles()).c_str(), BinEdges.size()-1, &BinEdges[0]); fDataHist->SetDirectory(0); - for (unsigned int i = 0; i < BinEdges.size(); ++i) { + for (unsigned int i = 0; i < BinEdges.size()-1; ++i) { fDataHist->SetBinContent(i+1, CrossSection[i]); fDataHist->SetBinError(i+1, Error[i]); } fDataHist->GetXaxis()->SetTitle(temp->GetXaxis()->GetTitle()); fDataHist->GetYaxis()->SetTitle(temp->GetYaxis()->GetTitle()); fDataHist->SetTitle(temp->GetTitle()); File->Close(); } //******************************************************************** // Covariance also needs stripping out // There's padding (two bins...) and overflow (last bin before the two empty bins) void MINERvA_CC0pinp_STV_XSec_1D_nu::SetCovarianceFromRootFile(std::string filename) { //******************************************************************** std::vector tempfile = GeneralUtils::ParseToStr(filename, ";"); TFile *File = new TFile(tempfile[0].c_str(), "READ"); // First object is the data, second is data with statistical error only, third is the covariance matrix TMatrixDSym *tempcov = (TMatrixDSym*)((TList*)File->Get(tempfile[1].c_str()))->At(2); - // Strip out the first two and last two columns - int nbins = tempcov->GetNrows(); // Count the number of zero entries int ngood = 0; int nstart = -1; int nend = -1; // Scan through the middle bin and look for entries int middle = tempcov->GetNrows()/2; + int nbinsdata = fDataHist->GetXaxis()->GetNbins(); + for (int j = 0; j < tempcov->GetNrows(); ++j) { - if ((*tempcov)(middle,j) > 0) { + if ((*tempcov)(middle,j) > 0 && ngood < nbinsdata) { ngood++; if (nstart == -1) nstart = j; if (j > nend) nend = j; } } - ngood--; fFullCovar = new TMatrixDSym(ngood); for (int i = 0; i < fFullCovar->GetNrows(); ++i) { for (int j = 0; j < fFullCovar->GetNrows(); ++j) { - (*fFullCovar)(i,j) = (*tempcov)(i+nstart, j+nstart); + (*fFullCovar)(i,j) = (*tempcov)(i+nstart+startbin-1, j+nstart+startbin-1); } } (*fFullCovar) *= 1E38*1E38; File->Close(); // Fill other covars. covar = StatUtils::GetInvert(fFullCovar); fDecomp = StatUtils::GetDecomp(fFullCovar); } void MINERvA_CC0pinp_STV_XSec_1D_nu::FillEventVariables(FitEvent *event) { fXVar = -999.99; // Need a proton and a muon if (event->NumFSParticle(2212) == 0 || event->NumFSParticle(13) == 0) { return; } TLorentzVector Pnu = event->GetNeutrinoIn()->fP; TLorentzVector Pmu = event->GetHMFSParticle(13)->fP; - TLorentzVector Pprot = event->GetHMFSParticle(2212)->fP; + // Find the highest momentum proton in the event between 450 and 1200 MeV with theta_p < 70 + int HMFSProton = 0; + double HighestMomentum = 0.0; + // Get the stack of protons + std::vector Protons = event->GetAllFSProton(); + for (size_t i = 0; i < Protons.size(); ++i) { + if (Protons[i]->p() > 450 && + Protons[i]->p() < 1200 && + Protons[i]->P3().Angle(Pnu.Vect()) < (M_PI/180.0)*70.0 && + Protons[i]->p() > HighestMomentum) { + HMFSProton = i; + } + } + // Now get the proton + TLorentzVector Pprot = Protons[HMFSProton]->fP; switch (fDist) { case (kMuonMom): fXVar = Pmu.Vect().Mag()/1000.0; break; case (kMuonTh): fXVar = Pmu.Vect().Angle(Pnu.Vect()) * (180.0/M_PI); break; case (kPrMom): fXVar = Pprot.Vect().Mag()/1000.0; break; case (kPrTh): fXVar = Pprot.Vect().Angle(Pnu.Vect()) * (180.0/M_PI); break; // Use Stephen's validated functions case (kNeMom): fXVar = FitUtils::Get_pn_reco_C(event, 14, true); break; case (kDalphaT): fXVar = FitUtils::Get_STV_dalphat(event, 14, true)*(180.0/M_PI); break; case (kDpT): fXVar = FitUtils::Get_STV_dpt(event, 14, true)/1000.0; break; case (kDphiT): fXVar = FitUtils::Get_STV_dphit(event, 14, true)*(180.0/M_PI); break; } return; }; //******************************************************************** bool MINERvA_CC0pinp_STV_XSec_1D_nu::isSignal(FitEvent *event) //******************************************************************** { return SignalDef::isCC0piNp_MINERvA_STV(event, EnuMin, EnuMax); } diff --git a/src/MINERvA/MINERvA_CC0pinp_STV_XSec_1D_nu.h b/src/MINERvA/MINERvA_CC0pinp_STV_XSec_1D_nu.h index 6950d34..5b88e40 100644 --- a/src/MINERvA/MINERvA_CC0pinp_STV_XSec_1D_nu.h +++ b/src/MINERvA/MINERvA_CC0pinp_STV_XSec_1D_nu.h @@ -1,53 +1,58 @@ // 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 MINERvA_CC0PINP_STV_XSEC_1D_NU_H_SEEN #define MINERvA_CC0PINP_STV_XSEC_1D_NU_H_SEEN #include "Measurement1D.h" class MINERvA_CC0pinp_STV_XSec_1D_nu : public Measurement1D { public: MINERvA_CC0pinp_STV_XSec_1D_nu(nuiskey samplekey); virtual ~MINERvA_CC0pinp_STV_XSec_1D_nu() {}; void FillEventVariables(FitEvent *event); bool isSignal(FitEvent *event); private: void SetupDataSettings(); void SetDataFromRootFile(std::string); void SetCovarianceFromRootFile(std::string); enum DataDistribution { kMuonMom, kMuonTh, kPrMom, kPrTh, kNeMom, kDalphaT, kDpT, kDphiT } MINERvA_CC0pi_STV_DataDistributions; DataDistribution fDist; + + // For truncating + double fMin; + double fMax; + int startbin; }; #endif diff --git a/src/MINERvA/MINERvA_SignalDef.cxx b/src/MINERvA/MINERvA_SignalDef.cxx index 1dd6084..2127052 100644 --- a/src/MINERvA/MINERvA_SignalDef.cxx +++ b/src/MINERvA/MINERvA_SignalDef.cxx @@ -1,416 +1,425 @@ // 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 "SignalDef.h" #include "FitUtils.h" #include "MINERvA_SignalDef.h" namespace SignalDef { // ********************************* // MINERvA CC1pi+/- signal definition (2015 release) // Note: There is a 2016 release which is different to this (listed below), but // it is CCNpi+ and has a different W cut // Note2: The W cut is implemented in the class implementation in MINERvA/ // rather than here so we can draw events that don't pass the W cut (W cut is // 1.4 GeV) // Could possibly be changed for slight speed increase since less events // would be used // // MINERvA signal is slightly different to MiniBooNE // // Exactly one negative muon // Exactly one charged pion (both + and -); however, there is a Michel e- // requirement but UNCLEAR IF UNFOLDED OR NOT (so don't know if should be // applied) // Exactly 1 charged pion exits (so + and - charge), however, has Michel // electron requirement, so look for + only here? // No restriction on neutral pions or other mesons // MINERvA has unfolded and not unfolded muon phase space for 2015 // // Possible issues with the data: // 1) pi- is allowed in signal even when Michel cut included; most pi- is efficiency corrected in GENIE // 2) There is a T_pi < 350 MeV cut coming from requiring a stopping pion; this is efficiency corrected in GENIE // 3) There is a 1.5 < Enu < 10.0 cut in signal definition // 4) There is an angular muon cut which is sometimes efficiency corrected (why we have bool isRestricted below) // // Nice things: // Much data given: with and without muon angle cuts and with and without shape // only data + covariance // bool isCC1pip_MINERvA(FitEvent *event, double EnuMin, double EnuMax, bool isRestricted) { // ********************************* // Signal is both pi+ and pi- // WARNING: PI- CONTAMINATION IS FULLY GENIE BECAUSE THE MICHEL TAG // First, make sure it's CCINC if (!isCCINC(event, 14, EnuMin, EnuMax)) return false; // Allow pi+/pi- int piPDG[] = {211, -211}; int nLeptons = event->NumFSLeptons(); int nPion = event->NumFSParticle(piPDG); // Check that the desired pion exists and is the only meson if (nPion != 1) return false; // Check that there is only one final state lepton if (nLeptons != 1) return false; // MINERvA released another CC1pi+ xsec without muon unfolding! // here the muon angle is < 20 degrees (seen in MINOS) TLorentzVector pnu = event->GetHMISParticle(14)->fP; TLorentzVector pmu = event->GetHMFSParticle(13)->fP; if (isRestricted) { double th_nu_mu = FitUtils::th(pmu, pnu) * 180. / M_PI; if (th_nu_mu >= 20) return false; } // Extract Hadronic Mass double hadMass = FitUtils::Wrec(pnu, pmu); // Actual cut is True GENIE Ws! Arg.! Use gNtpcConv definition. #ifdef __GENIE_ENABLED__ if (event->fType == kGENIE){ EventRecord * gevent = static_cast(event->genie_event->event); const Interaction * interaction = gevent->Summary(); const Kinematics & kine = interaction->Kine(); double Ws = kine.W (true); // std::cout << "Ws versus WRec = " << Ws << " vs " << hadMass << " " << kine.W(false) << std::endl; hadMass = Ws * 1000.0; } #endif if (hadMass > 1400.0) return false; return true; }; // Updated MINERvA 2017 Signal using Wexp and no restriction on angle bool isCC1pip_MINERvA_2017(FitEvent *event, double EnuMin, double EnuMax){ // Signal is both pi+ and pi- // WARNING: PI- CONTAMINATION IS FULLY GENIE BECAUSE THE MICHEL TAG // First, make sure it's CCINC if (!isCCINC(event, 14, EnuMin, EnuMax)) return false; // Allow pi+/pi- int piPDG[] = {211, -211}; int nLeptons = event->NumFSLeptons(); int nPion = event->NumFSParticle(piPDG); // Check that the desired pion exists and is the only meson if (nPion != 1) return false; // Check that there is only one final state lepton if (nLeptons != 1) return false; // Get Muon and Lepton Kinematics TLorentzVector pnu = event->GetHMISParticle(14)->fP; TLorentzVector pmu = event->GetHMFSParticle(13)->fP; // Extract Hadronic Mass double hadMass = FitUtils::Wrec(pnu, pmu); if (hadMass > 1400.0) return false; return true; }; // ********************************* // MINERvA CCNpi+/- signal definition from 2016 publication // Different to CC1pi+/- listed above; additional has W < 1.8 GeV // // For notes on strangeness of signal definition, see CC1pip_MINERvA // // One negative muon // At least one charged pion // 1.5 < Enu < 10 // No restrictions on pi0 or other mesons or baryons // W_reconstructed (ignoring initial state motion) cut at 1.8 GeV // // Also writes number of pions (nPions) if studies on this want to be done... bool isCCNpip_MINERvA(FitEvent *event, double EnuMin, double EnuMax, bool isRestricted, bool isWtrue) { // ********************************* // First, make sure it's CCINC if (!isCCINC(event, 14, EnuMin, EnuMax)) return false; int nLeptons = event->NumFSLeptons(); // Write the number of pions to the measurement class... // Maybe better to just do that inside the class? int nPions = event->NumFSParticle(PhysConst::pdg_charged_pions); // Check that there is a pion! if (nPions == 0) return false; // Check that there is only one final state lepton if (nLeptons != 1) return false; // Need the muon and the neutrino to check angles and W TLorentzVector pnu = event->GetNeutrinoIn()->fP; TLorentzVector pmu = event->GetHMFSParticle(13)->fP; // MINERvA released some data with restricted muon angle // Here the muon angle is < 20 degrees (seen in MINOS) if (isRestricted) { double th_nu_mu = FitUtils::th(pmu, pnu) * 180. / M_PI; if (th_nu_mu >= 20.) return false; } // Lastly check the W cut (always at 1.8 GeV) double Wrec = FitUtils::Wrec(pnu, pmu) + 0.; // Actual cut is True GENIE Ws! Arg.! Use gNtpcConv definition. if (isWtrue){ #ifdef __GENIE_ENABLED__ if (event->fType == kGENIE){ GHepRecord* ghep = static_cast(event->genie_event->event); const Interaction * interaction = ghep->Summary(); const Kinematics & kine = interaction->Kine(); double Ws = kine.W (true); Wrec = Ws * 1000.0; // Say Wrec is Ws } #endif } if (Wrec > 1800. || Wrec < 0.0) return false; return true; }; //******************************************************************** bool isCCQEnumu_MINERvA(FitEvent *event, double EnuMin, double EnuMax, bool fullphasespace) { //******************************************************************** if (!isCCQELike(event, 14, EnuMin, EnuMax)) return false; TLorentzVector pnu = event->GetHMISParticle(14)->fP; TLorentzVector pmu = event->GetHMFSParticle(13)->fP; double ThetaMu = pnu.Vect().Angle(pmu.Vect()); double Enu_rec = FitUtils::EnuQErec(pmu, cos(ThetaMu), 34., true); // If Restricted phase space if (!fullphasespace && ThetaMu > 0.34906585) return false; // restrict energy range if (Enu_rec < EnuMin || Enu_rec > EnuMax) return false; return true; }; //******************************************************************** bool isCCQEnumubar_MINERvA(FitEvent *event, double EnuMin, double EnuMax, bool fullphasespace) { //******************************************************************** if (!isCCQELike(event, -14, EnuMin, EnuMax)) return false; TLorentzVector pnu = event->GetHMISParticle(-14)->fP; TLorentzVector pmu = event->GetHMFSParticle(-13)->fP; double ThetaMu = pnu.Vect().Angle(pmu.Vect()); double Enu_rec = FitUtils::EnuQErec(pmu, cos(ThetaMu), 30., true); // If Restricted phase space if (!fullphasespace && ThetaMu > 0.34906585) return false; // restrict energy range if (Enu_rec < EnuMin || Enu_rec > EnuMax) return false; return true; } //******************************************************************** bool isCCincLowRecoil_MINERvA(FitEvent *event, double EnuMin, double EnuMax) { //******************************************************************** if (!isCCINC(event, 14, EnuMin, EnuMax)) return false; // Need at least one muon if (event->NumFSParticle(13) < 1) return false; TLorentzVector pmu = event->GetHMFSParticle(13)->fP; TLorentzVector pnu = event->GetHMISParticle(14)->fP; // Cut on muon angle greated than 20deg if (cos(pnu.Vect().Angle(pmu.Vect())) < 0.93969262078) return false; // Cut on muon energy < 1.5 GeV if (pmu.E()/1000.0 < 1.5) return false; return true; } bool isCC0pi1p_MINERvA(FitEvent *event, double enumin, double enumax) { // Require numu CC0pi event with a proton above threshold bool signal = (isCC0pi(event, 14, enumin, enumax) && HasProtonKEAboveThreshold(event, 110.0)); return signal; } // 2015 analysis just asks for 1pi0 and no pi+/pi- bool isCC1pi0_MINERvA_2015(FitEvent *event, double EnuMin, double EnuMax) { bool CC1pi0_anu = SignalDef::isCC1pi(event, -14, 111, EnuMin, EnuMax); return CC1pi0_anu; } // 2016 analysis just asks for 1pi0 and no other charged tracks. Half-open to interpretation: we go with "charged tracks" meaning pions. You'll be forgiven for thinking proton tracks should be included here too but we checked with MINERvA bool isCC1pi0_MINERvA_2016(FitEvent *event, double EnuMin, double EnuMax) { bool CC1pi0_anu = SignalDef::isCC1pi(event, -14, 111, EnuMin, EnuMax); /* // Additionally look for charged proton track // Comment: This is _NOT_ in the signal definition but was tested bool HasProton = event->HasFSParticle(2212); if (CC1pi0_anu) { if (!HasProton) { return true; } else { return false; } } else { return false; } */ return CC1pi0_anu; } // 2016 analysis just asks for 1pi0 and no other charged tracks bool isCC1pi0_MINERvA_nu(FitEvent *event, double EnuMin, double EnuMax) { bool CC1pi0_nu = SignalDef::isCC1pi(event, 14, 111, EnuMin, EnuMax); return CC1pi0_nu; } //******************************************************************** bool isCC0pi_MINERvAPTPZ(FitEvent* event, int nuPDG, double emin, double emax){ //******************************************************************** // Check it's CCINC if (!SignalDef::isCCINC(event, nuPDG, emin, emax)) return false; // Make Angle Cut > 20.0 TLorentzVector pnu = event->GetHMISParticle(14)->fP; TLorentzVector pmu = event->GetHMFSParticle(13)->fP; double th_nu_mu = FitUtils::th(pmu, pnu) * 180. / M_PI; if (th_nu_mu >= 20.0) return false; int genie_n_muons = 0; int genie_n_mesons = 0; int genie_n_heavy_baryons_plus_pi0s = 0; int genie_n_photons = 0; for(unsigned int i = 0; i < event->NParticles(); ++i) { FitParticle* p = event->GetParticle(i); if (p->Status() != kFinalState) continue; int pdg = p->fPID; double energy = p->fP.E(); if( abs(pdg) == 13 ) genie_n_muons++; else if( pdg == 22 && energy > 10.0 ) genie_n_photons++; else if( abs(pdg) == 211 || abs(pdg) == 321 || abs(pdg) == 323 || pdg == 111 || pdg == 130 || pdg == 310 || pdg == 311 || pdg == 313 ) genie_n_mesons++; else if( pdg == 3112 || pdg == 3122 || pdg == 3212 || pdg == 3222 || pdg == 4112 || pdg == 4122 || pdg == 4212 || pdg == 4222 || pdg == 411 || pdg == 421 || pdg == 111 ) genie_n_heavy_baryons_plus_pi0s++; } if( genie_n_muons == 1 && genie_n_mesons == 0 && genie_n_heavy_baryons_plus_pi0s == 0 && genie_n_photons == 0 ) return true; return false; } bool isCC0pi_anti_MINERvAPTPZ(FitEvent* event, int nuPDG, double emin, double emax){ // Check it's CCINC if (!SignalDef::isCCINC(event, nuPDG, emin, emax)) return false; // Make Angle Cut > 20.0 TLorentzVector pnu = event->GetHMISParticle(-14)->fP; TLorentzVector pmu = event->GetHMFSParticle(-13)->fP; double th_nu_mu = FitUtils::th(pmu, pnu) * 180. / M_PI; if (th_nu_mu >= 20.0) return false; int genie_n_muons = 0; int genie_n_mesons = 0; int genie_n_heavy_baryons_plus_pi0s = 0; int genie_n_photons = 0; for(unsigned int i = 0; i < event->NParticles(); ++i) { FitParticle* p = event->GetParticle(i); if (p->Status() != kFinalState) continue; int pdg = p->fPID; double energy = p->fP.E(); if( abs(pdg) == 13 ) genie_n_muons++; else if( pdg == 22 && energy > 10.0 ) genie_n_photons++; else if( abs(pdg) == 211 || abs(pdg) == 321 || abs(pdg) == 323 || abs(pdg) == 111 || abs(pdg) == 130 || abs(pdg) == 310 || abs(pdg) == 311 || abs(pdg) == 313 ) genie_n_mesons++; else if( abs(pdg) == 3112 || abs(pdg) == 3122 || abs(pdg) == 3212 || abs(pdg) == 3222 || abs(pdg) == 4112 || abs(pdg) == 4122 || abs(pdg) == 4212 || abs(pdg) == 4222 || abs(pdg) == 411 || abs(pdg) == 421 || abs(pdg) == 111 ) genie_n_heavy_baryons_plus_pi0s++; } if( genie_n_muons == 1 && genie_n_mesons == 0 && genie_n_heavy_baryons_plus_pi0s == 0 && genie_n_photons == 0 ) return true; return false; } // MINERvA CC0pi transverse variables signal defintion bool isCC0piNp_MINERvA_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 momentum cuts if (pmu.Vect().Mag() < 1500 || pmu.Vect().Mag() > 10000) return false; // Muon angle cuts if (pmu.Vect().Angle(pnu.Vect()) > (M_PI/180.0)*20.0) return false; + // Since there are emany protons allowed we can't just use the highest proton momentum candidate and place cuts on it + // Get the stack of protons + std::vector Protons = event->GetAllFSProton(); + // Count how many protons pass the threshold + int nProtonsAboveThreshold = 0; + for (size_t i = 0; i < Protons.size(); ++i) { + if (Protons[i]->p() > 450 && Protons[i]->p() < 1200 && Protons[i]->P3().Angle(pnu.Vect()) < (M_PI/180.0)*70.0) nProtonsAboveThreshold++; + } + // Proton momentum cuts - if (pp.Vect().Mag() < 450 || pp.Vect().Mag() > 1200) return false; + //if (pp.Vect().Mag() < 450 || pp.Vect().Mag() > 1200) return false; // Proton angle cuts - if (pp.Vect().Angle(pnu.Vect()) > (M_PI/180.0)*70) return false; + //if (pp.Vect().Angle(pnu.Vect()) > (M_PI/180.0)*70) return false; + if (nProtonsAboveThreshold == 0) return false; return true; }; } diff --git a/src/Utils/FitUtils.cxx b/src/Utils/FitUtils.cxx index 4ffa992..af425c6 100644 --- a/src/Utils/FitUtils.cxx +++ b/src/Utils/FitUtils.cxx @@ -1,978 +1,1048 @@ // 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 "FitUtils.h" /* MISC Functions */ //******************************************************************** double *FitUtils::GetArrayFromMap(std::vector invals, std::map inmap) { //******************************************************************** double *outarr = new double[invals.size()]; int count = 0; for (size_t i = 0; i < invals.size(); i++) { outarr[count++] = inmap[invals.at(i)]; } return outarr; } /* MISC Event */ //******************************************************************** // Returns the kinetic energy of a particle in GeV double FitUtils::T(TLorentzVector part) { //******************************************************************** double E_part = part.E() / 1000.; double p_part = part.Vect().Mag() / 1000.; double m_part = sqrt(E_part * E_part - p_part * p_part); double KE_part = E_part - m_part; return KE_part; }; //******************************************************************** // Returns the momentum of a particle in GeV double FitUtils::p(TLorentzVector part) { //******************************************************************** double p_part = part.Vect().Mag() / 1000.; return p_part; }; double FitUtils::p(FitParticle *part) { return FitUtils::p(part->fP); }; //******************************************************************** // Returns the angle between two particles in radians double FitUtils::th(TLorentzVector part1, TLorentzVector part2) { //******************************************************************** double th = part1.Vect().Angle(part2.Vect()); return th; }; double FitUtils::th(FitParticle *part1, FitParticle *part2) { return FitUtils::th(part1->fP, part2->fP); }; // T2K CC1pi+ helper functions // //******************************************************************** // Returns the angle between q3 and the pion defined in Raquel's CC1pi+ on CH // paper // Uses "MiniBooNE formula" for Enu, here called EnuCC1pip_T2K_MB //******************************************************************** double FitUtils::thq3pi_CC1pip_T2K(TLorentzVector pnu, TLorentzVector pmu, TLorentzVector ppi) { // Want this in GeV TVector3 p_mu = pmu.Vect() * (1. / 1000.); // Get the reconstructed Enu // We are not using Michel e sample, so we have pion kinematic information double Enu = EnuCC1piprec(pnu, pmu, ppi, true); // Get neutrino unit direction, multiply by reconstructed Enu TVector3 p_nu = pnu.Vect() * (1. / (pnu.Vect().Mag())) * Enu; TVector3 p_pi = ppi.Vect() * (1. / 1000.); // This is now in GeV TVector3 q3 = (p_nu - p_mu); // Want this in GeV double th_q3_pi = q3.Angle(p_pi); return th_q3_pi; } //******************************************************************** // Returns the q3 defined in Raquel's CC1pi+ on CH paper // Uses "MiniBooNE formula" for Enu //******************************************************************** double FitUtils::q3_CC1pip_T2K(TLorentzVector pnu, TLorentzVector pmu, TLorentzVector ppi) { // Can't use the true Enu here; need to reconstruct it // We do have Michel e- here so reconstruct Enu by "MiniBooNE formula" without // pion kinematics // The bool false refers to that we don't have pion kinematics // Last bool refers to if we have pion kinematic information or not double Enu = EnuCC1piprec(pnu, pmu, ppi, false); TVector3 p_mu = pmu.Vect() * (1. / 1000.); TVector3 p_nu = pnu.Vect() * (1. / pnu.Vect().Mag()) * Enu; double q3 = (p_nu - p_mu).Mag(); return q3; } //******************************************************************** // Returns the W reconstruction from Raquel CC1pi+ CH thesis // Uses the MiniBooNE formula Enu //******************************************************************** double FitUtils::WrecCC1pip_T2K_MB(TLorentzVector pnu, TLorentzVector pmu, TLorentzVector ppi) { double E_mu = pmu.E() / 1000.; double p_mu = pmu.Vect().Mag() / 1000.; double E_nu = EnuCC1piprec(pnu, pmu, ppi, false); double a1 = (E_nu + PhysConst::mass_neutron) - E_mu; double a2 = E_nu - p_mu; double wrec = sqrt(a1 * a1 - a2 * a2); return wrec; } //******************************************************** double FitUtils::ProtonQ2QErec(double pE, double binding) { //******************************************************** const double V = binding / 1000.; // binding potential const double mn = PhysConst::mass_neutron; // neutron mass const double mp = PhysConst::mass_proton; // proton mass const double mn_eff = mn - V; // effective proton mass const double pki = (pE / 1000.0) - mp; double q2qe = mn_eff * mn_eff - mp * mp + 2 * mn_eff * (pki + mp - mn_eff); return q2qe; }; //******************************************************************** double FitUtils::EnuQErec(TLorentzVector pmu, double costh, double binding, bool neutrino) { //******************************************************************** // Convert all values to GeV const double V = binding / 1000.; // binding potential const double mn = PhysConst::mass_neutron; // neutron mass const double mp = PhysConst::mass_proton; // proton mass double mN_eff = mn - V; double mN_oth = mp; if (!neutrino) { mN_eff = mp - V; mN_oth = mn; } double el = pmu.E() / 1000.; double pl = (pmu.Vect().Mag()) / 1000.; // momentum of lepton double ml = sqrt(el * el - pl * pl); // lepton mass double rEnu = (2 * mN_eff * el - ml * ml + mN_oth * mN_oth - mN_eff * mN_eff) / (2 * (mN_eff - el + pl * costh)); return rEnu; }; double FitUtils::Q2QErec(TLorentzVector pmu, double costh, double binding, bool neutrino) { double el = pmu.E() / 1000.; double pl = (pmu.Vect().Mag()) / 1000.; // momentum of lepton double ml = sqrt(el * el - pl * pl); // lepton mass double rEnu = EnuQErec(pmu, costh, binding, neutrino); double q2 = -ml * ml + 2. * rEnu * (el - pl * costh); return q2; }; double FitUtils::Q2QErec(TLorentzVector Pmu, TLorentzVector Pnu, double binding, bool neutrino) { double q2qe = Q2QErec(Pmu, cos(Pnu.Vect().Angle(Pmu.Vect())), binding, neutrino); return q2qe; } double FitUtils::EnuQErec(double pl, double costh, double binding, bool neutrino) { if (pl < 0) return 0.; // Make sure nobody is silly double mN_eff = PhysConst::mass_neutron - binding / 1000.; double mN_oth = PhysConst::mass_proton; if (!neutrino) { mN_eff = PhysConst::mass_proton - binding / 1000.; mN_oth = PhysConst::mass_neutron; } double ml = PhysConst::mass_muon; double el = sqrt(pl * pl + ml * ml); double rEnu = (2 * mN_eff * el - ml * ml + mN_oth * mN_oth - mN_eff * mN_eff) / (2 * (mN_eff - el + pl * costh)); return rEnu; }; double FitUtils::Q2QErec(double pl, double costh, double binding, bool neutrino) { if (pl < 0) return 0.; // Make sure nobody is silly double ml = PhysConst::mass_muon; double el = sqrt(pl * pl + ml * ml); double rEnu = EnuQErec(pl, costh, binding, neutrino); double q2 = -ml * ml + 2. * rEnu * (el - pl * costh); return q2; }; //******************************************************************** // Reconstructs Enu for CC1pi0 // Very similar for CC1pi+ reconstruction double FitUtils::EnuCC1pi0rec(TLorentzVector pnu, TLorentzVector pmu, TLorentzVector ppi0) { //******************************************************************** double E_mu = pmu.E() / 1000; double p_mu = pmu.Vect().Mag() / 1000; double m_mu = sqrt(E_mu * E_mu - p_mu * p_mu); double th_nu_mu = pnu.Vect().Angle(pmu.Vect()); double E_pi0 = ppi0.E() / 1000; double p_pi0 = ppi0.Vect().Mag() / 1000; double m_pi0 = sqrt(E_pi0 * E_pi0 - p_pi0 * p_pi0); double th_nu_pi0 = pnu.Vect().Angle(ppi0.Vect()); const double m_n = PhysConst::mass_neutron; // neutron mass const double m_p = PhysConst::mass_proton; // proton mass double th_pi0_mu = ppi0.Vect().Angle(pmu.Vect()); double rEnu = (m_mu * m_mu + m_pi0 * m_pi0 + m_n * m_n - m_p * m_p - 2 * m_n * (E_pi0 + E_mu) + 2 * E_pi0 * E_mu - 2 * p_pi0 * p_mu * cos(th_pi0_mu)) / (2 * (E_pi0 + E_mu - p_pi0 * cos(th_nu_pi0) - p_mu * cos(th_nu_mu) - m_n)); return rEnu; }; //******************************************************************** // Reconstruct Q2 for CC1pi0 // Beware: uses true Enu, not reconstructed Enu double FitUtils::Q2CC1pi0rec(TLorentzVector pnu, TLorentzVector pmu, TLorentzVector ppi0) { //******************************************************************** double E_mu = pmu.E() / 1000.; // energy of lepton in GeV double p_mu = pmu.Vect().Mag() / 1000.; // momentum of lepton double m_mu = sqrt(E_mu * E_mu - p_mu * p_mu); // lepton mass double th_nu_mu = pnu.Vect().Angle(pmu.Vect()); // double rEnu = EnuCC1pi0rec(pnu, pmu, ppi0); //reconstructed neutrino energy // Use true neutrino energy double rEnu = pnu.E() / 1000.; double q2 = -m_mu * m_mu + 2. * rEnu * (E_mu - p_mu * cos(th_nu_mu)); return q2; }; //******************************************************************** // Reconstruct Enu for CC1pi+ // pionInfo reflects if we're using pion kinematics or not // In T2K CC1pi+ CH the Michel tag is used for pion in which pion kinematic info // is lost and Enu is reconstructed without pion kinematics double FitUtils::EnuCC1piprec(TLorentzVector pnu, TLorentzVector pmu, TLorentzVector ppi, bool pionInfo) { //******************************************************************** double E_mu = pmu.E() / 1000.; double p_mu = pmu.Vect().Mag() / 1000.; double m_mu = sqrt(E_mu * E_mu - p_mu * p_mu); double E_pi = ppi.E() / 1000.; double p_pi = ppi.Vect().Mag() / 1000.; double m_pi = sqrt(E_pi * E_pi - p_pi * p_pi); const double m_n = PhysConst::mass_neutron; // neutron/proton mass // should really take proton mass for proton interaction, neutron for neutron // interaction. However, difference is pretty much negligable here! // need this angle too double th_nu_pi = pnu.Vect().Angle(ppi.Vect()); double th_nu_mu = pnu.Vect().Angle(pmu.Vect()); double th_pi_mu = ppi.Vect().Angle(pmu.Vect()); double rEnu = -999; // If experiment doesn't have pion kinematic information (T2K CC1pi+ CH // example of this; Michel e sample has no directional information on pion) // We'll need to modify the reconstruction variables if (!pionInfo) { // Assumes that pion angle contribution contributes net zero // Also assumes the momentum of reconstructed pions when using Michel // electrons is 130 MeV // So need to redefine E_pi and p_pi // It's a little unclear what happens to the pmu . ppi term (4-vector); do // we include the angular contribution? // This below doesn't th_nu_pi = M_PI / 2.; p_pi = 0.130; E_pi = sqrt(p_pi * p_pi + m_pi * m_pi); // rEnu = (m_mu*m_mu + m_pi*m_pi - 2*m_n*(E_mu + E_pi) + 2*E_mu*E_pi - // 2*p_mu*p_pi) / (2*(E_mu + E_pi - p_mu*cos(th_nu_mu) - m_n)); } rEnu = (m_mu * m_mu + m_pi * m_pi - 2 * m_n * (E_pi + E_mu) + 2 * E_pi * E_mu - 2 * p_pi * p_mu * cos(th_pi_mu)) / (2 * (E_pi + E_mu - p_pi * cos(th_nu_pi) - p_mu * cos(th_nu_mu) - m_n)); return rEnu; }; //******************************************************************** // Reconstruct neutrino energy from outgoing particles; will differ from the // actual neutrino energy. Here we use assumption of a Delta resonance double FitUtils::EnuCC1piprecDelta(TLorentzVector pnu, TLorentzVector pmu) { //******************************************************************** const double m_Delta = PhysConst::mass_delta; // PDG value for Delta mass in GeV const double m_n = PhysConst::mass_neutron; // neutron/proton mass // should really take proton mass for proton interaction, neutron for neutron // interaction. However, difference is pretty much negligable here! double E_mu = pmu.E() / 1000; double p_mu = pmu.Vect().Mag() / 1000; double m_mu = sqrt(E_mu * E_mu - p_mu * p_mu); double th_nu_mu = pnu.Vect().Angle(pmu.Vect()); double rEnu = (m_Delta * m_Delta - m_n * m_n - m_mu * m_mu + 2 * m_n * E_mu) / (2 * (m_n - E_mu + p_mu * cos(th_nu_mu))); return rEnu; }; // MOVE TO T2K UTILS! //******************************************************************** // Reconstruct Enu using "extended MiniBooNE" as defined in Raquel's T2K TN // // Supposedly includes pion direction and binding energy of target nucleon // I'm not convinced (yet), maybe double FitUtils::EnuCC1piprec_T2K_eMB(TLorentzVector pnu, TLorentzVector pmu, TLorentzVector ppi) { //******************************************************************** // Unit vector for neutrino momentum TVector3 p_nu_vect_unit = pnu.Vect() * (1. / pnu.E()); double E_mu = pmu.E() / 1000.; TVector3 p_mu_vect = pmu.Vect() * (1. / 1000.); double E_pi = ppi.E() / 1000.; TVector3 p_pi_vect = ppi.Vect() * (1. / 1000.); double E_bind = 27. / 1000.; // This should be roughly correct for CH; but not clear! double m_p = PhysConst::mass_proton; // Makes life a little easier, gonna square this one double a1 = m_p - E_bind - E_mu - E_pi; // Just gets the magnitude square of the muon and pion momentum vectors double a2 = (p_mu_vect + p_pi_vect).Mag2(); // Gets the somewhat complicated scalar product between neutrino and // (p_mu+p_pi), scaled to Enu double a3 = p_nu_vect_unit * (p_mu_vect + p_pi_vect); double rEnu = (m_p * m_p - a1 * a1 + a2) / (2 * (m_p - E_bind - E_mu - E_pi + a3)); return rEnu; } //******************************************************************** // Reconstructed Q2 for CC1pi+ // // enuType describes how the neutrino energy is reconstructed // 0 uses true neutrino energy; case for MINERvA and MiniBooNE // 1 uses "extended MiniBooNE" reconstructed (T2K CH) // 2 uses "MiniBooNE" reconstructed (EnuCC1piprec with pionInfo = true) (T2K CH) // "MiniBooNE" reconstructed (EnuCC1piprec with pionInfo = false, the // case for T2K when using Michel tag) (T2K CH) // 3 uses Delta for reconstruction (T2K CH) double FitUtils::Q2CC1piprec(TLorentzVector pnu, TLorentzVector pmu, TLorentzVector ppi, int enuType, bool pionInfo) { //******************************************************************** double E_mu = pmu.E() / 1000.; // energy of lepton in GeV double p_mu = pmu.Vect().Mag() / 1000.; // momentum of lepton double m_mu = sqrt(E_mu * E_mu - p_mu * p_mu); // lepton mass double th_nu_mu = pnu.Vect().Angle(pmu.Vect()); // Different ways of reconstructing the neutrino energy; default is just to // use the truth (case 0) double rEnu = -999; switch (enuType) { case 0: // True neutrino energy, default; this is the case for when Q2 // defined as Q2 true (MiniBooNE, MINERvA) rEnu = pnu.E() / 1000.; break; case 1: // Extended MiniBooNE reconstructed, as defined by Raquel's CC1pi+ // CH T2K analysis // Definitely uses pion info :) rEnu = EnuCC1piprec_T2K_eMB(pnu, pmu, ppi); break; case 2: // MiniBooNE reconstructed, as defined by MiniBooNE and Raquel's // CC1pi+ CH T2K analysis and Linda's CC1pi+ H2O // Can have this with and without pion info, depending on if Michel tag // used (Raquel) rEnu = EnuCC1piprec(pnu, pmu, ppi, pionInfo); break; case 3: rEnu = EnuCC1piprecDelta(pnu, pmu); break; } // No need for default here since default value of enuType = 0, defined in // header double q2 = -m_mu * m_mu + 2. * rEnu * (E_mu - p_mu * cos(th_nu_mu)); return q2; }; //******************************************************************** // Returns the reconstructed W from a nucleon and an outgoing pion // // Could do this a lot more clever (pp + ppi).Mag() would do the job, but this // would be less instructive //******************************************************************** double FitUtils::MpPi(TLorentzVector pp, TLorentzVector ppi) { double E_p = pp.E(); double p_p = pp.Vect().Mag(); double m_p = sqrt(E_p * E_p - p_p * p_p); double E_pi = ppi.E(); double p_pi = ppi.Vect().Mag(); double m_pi = sqrt(E_pi * E_pi - p_pi * p_pi); double th_p_pi = pp.Vect().Angle(ppi.Vect()); // fairly easy thing to derive since bubble chambers measure the proton! double invMass = sqrt(m_p * m_p + m_pi * m_pi + 2 * E_p * E_pi - 2 * p_pi * p_p * cos(th_p_pi)); return invMass; }; //******************************************************** // Reconstruct the hadronic mass using neutrino and muon // Could technically do E_nu = EnuCC1pipRec(pnu,pmu,ppi) too, but this wwill be // reconstructed Enu; so gives reconstructed Wrec which most of the time isn't // used! // // Only MINERvA uses this so far; and the Enu is Enu_true // If we want W_true need to take initial state nucleon motion into account // Return value is in MeV!!! double FitUtils::Wrec(TLorentzVector pnu, TLorentzVector pmu) { //******************************************************** double E_mu = pmu.E(); double p_mu = pmu.Vect().Mag(); double m_mu = sqrt(E_mu * E_mu - p_mu * p_mu); double th_nu_mu = pnu.Vect().Angle(pmu.Vect()); // The factor of 1000 is necessary for downstream functions const double m_p = PhysConst::mass_proton * 1000; // MINERvA cut on W_exp which is tuned to W_true; so use true Enu from // generators double E_nu = pnu.E(); double w_rec = sqrt(m_p * m_p + m_mu * m_mu - 2 * m_p * E_mu + 2 * E_nu * (m_p - E_mu + p_mu * cos(th_nu_mu))); return w_rec; }; //******************************************************** // Reconstruct the true hadronic mass using the initial state and muon // Could technically do E_nu = EnuCC1pipRec(pnu,pmu,ppi) too, but this wwill be // reconstructed Enu; so gives reconstructed Wrec which most of the time isn't // used! // // No one seems to use this because it's fairly MC dependent! // Return value is in MeV!!! double FitUtils::Wtrue(TLorentzVector pnu, TLorentzVector pmu, TLorentzVector pnuc) { //******************************************************** // Could simply do the TLorentzVector operators here but this is more // instructive? // ... and prone to errors ... double E_mu = pmu.E(); double p_mu = pmu.Vect().Mag(); double m_mu = sqrt(E_mu * E_mu - p_mu * p_mu); double th_nu_mu = pnu.Vect().Angle(pmu.Vect()); double E_nuc = pnuc.E(); double p_nuc = pnuc.Vect().Mag(); double m_nuc = sqrt(E_nuc * E_nuc - p_nuc * p_nuc); double th_nuc_mu = pmu.Vect().Angle(pnuc.Vect()); double th_nu_nuc = pnu.Vect().Angle(pnuc.Vect()); double E_nu = pnu.E(); double w_rec = sqrt(m_nuc * m_nuc + m_mu * m_mu - 2 * E_nu * E_mu + 2 * E_nu * p_mu * cos(th_nu_mu) - 2 * E_nuc * E_mu + 2 * p_nuc * p_mu * cos(th_nuc_mu) + 2 * E_nu * E_nuc - 2 * E_nu * p_nuc * cos(th_nu_nuc)); return w_rec; }; double FitUtils::SumKE_PartVect(std::vector const fps) { double sum = 0.0; for (size_t p_it = 0; p_it < fps.size(); ++p_it) { sum += fps[p_it]->KE(); } return sum; } double FitUtils::SumTE_PartVect(std::vector const fps) { double sum = 0.0; for (size_t p_it = 0; p_it < fps.size(); ++p_it) { sum += fps[p_it]->E(); } return sum; } /* E Recoil */ double FitUtils::GetErecoil_TRUE(FitEvent *event) { // Get total energy of hadronic system. double Erecoil = 0.0; for (unsigned int i = 2; i < event->Npart(); i++) { // Only final state if (!event->PartInfo(i)->fIsAlive) continue; if (event->PartInfo(i)->fNEUTStatusCode != 0) continue; // Skip Lepton if (abs(event->PartInfo(i)->fPID) == abs(event->PartInfo(0)->fPID) - 1) continue; // Add Up KE of protons and TE of everything else if (event->PartInfo(i)->fPID == 2212 || event->PartInfo(i)->fPID == 2112) { Erecoil += fabs(event->PartInfo(i)->fP.E()) - fabs(event->PartInfo(i)->fP.Mag()); } else { Erecoil += event->PartInfo(i)->fP.E(); } } return Erecoil; } double FitUtils::GetErecoil_CHARGED(FitEvent *event) { // Get total energy of hadronic system. double Erecoil = 0.0; for (unsigned int i = 2; i < event->Npart(); i++) { // Only final state if (!event->PartInfo(i)->fIsAlive) continue; if (event->PartInfo(i)->fNEUTStatusCode != 0) continue; // Skip Lepton if (abs(event->PartInfo(i)->fPID) == abs(event->PartInfo(0)->fPID) - 1) continue; // Skip Neutral particles if (event->PartInfo(i)->fPID == 2112 || event->PartInfo(i)->fPID == 111 || event->PartInfo(i)->fPID == 22) continue; // Add Up KE of protons and TE of everything else if (event->PartInfo(i)->fPID == 2212) { Erecoil += fabs(event->PartInfo(i)->fP.E()) - fabs(event->PartInfo(i)->fP.Mag()); } else { Erecoil += event->PartInfo(i)->fP.E(); } } return Erecoil; } // MOVE TO MINERVA Utils! double FitUtils::GetErecoil_MINERvA_LowRecoil(FitEvent *event) { // Get total energy of hadronic system. double Erecoil = 0.0; for (unsigned int i = 2; i < event->Npart(); i++) { // Only final state if (!event->PartInfo(i)->fIsAlive) continue; if (event->PartInfo(i)->fNEUTStatusCode != 0) continue; // Skip Lepton if (abs(event->PartInfo(i)->fPID) == 13) continue; // Skip Neutrons particles if (event->PartInfo(i)->fPID == 2112) continue; int PID = event->PartInfo(i)->fPID; // KE of Protons and charged pions if (PID == 2212 or PID == 211 or PID == -211) { // Erecoil += FitUtils::T(event->PartInfo(i)->fP); Erecoil += fabs(event->PartInfo(i)->fP.E()) - fabs(event->PartInfo(i)->fP.Mag()); // Total Energy of non-neutrons // } else if (PID != 2112 and PID < 999 and PID != 22 and abs(PID) != // 14) { } else if (PID == 111 || PID == 11 || PID == -11 || PID == 22) { Erecoil += (event->PartInfo(i)->fP.E()); } } return Erecoil; } TVector3 GetVectorInTPlane(const TVector3 &inp, const TVector3 &planarNormal) { TVector3 pnUnit = planarNormal.Unit(); double inpProjectPN = inp.Dot(pnUnit); TVector3 InPlanarInput = inp - (inpProjectPN * pnUnit); return InPlanarInput; } TVector3 GetUnitVectorInTPlane(const TVector3 &inp, const TVector3 &planarNormal) { return GetVectorInTPlane(inp, planarNormal).Unit(); } Double_t GetDeltaPhiT(TVector3 const &V_lepton, TVector3 const &V_other, TVector3 const &Normal, bool PiMinus = false) { TVector3 V_lepton_T = GetUnitVectorInTPlane(V_lepton, Normal); TVector3 V_other_T = GetUnitVectorInTPlane(V_other, Normal); return PiMinus ? acos(V_lepton_T.Dot(V_other_T)) : (M_PI - acos(V_lepton_T.Dot(V_other_T))); } TVector3 GetDeltaPT(TVector3 const &V_lepton, TVector3 const &V_other, TVector3 const &Normal) { TVector3 V_lepton_T = GetVectorInTPlane(V_lepton, Normal); TVector3 V_other_T = GetVectorInTPlane(V_other, Normal); return V_lepton_T + V_other_T; } Double_t GetDeltaAlphaT(TVector3 const &V_lepton, TVector3 const &V_other, TVector3 const &Normal, bool PiMinus = false) { TVector3 DeltaPT = GetDeltaPT(V_lepton, V_other, Normal); return GetDeltaPhiT(V_lepton, DeltaPT, Normal, PiMinus); } double FitUtils::Get_STV_dpt(FitEvent *event, int ISPDG, bool Is0pi) { // Check that the neutrino exists if (event->NumISParticle(ISPDG) == 0) { return -9999; } // Return 0 if the proton or muon are missing if (event->NumFSParticle(2212) == 0 || event->NumFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1)) == 0) { return -9999; } // Now get the TVector3s for each particle TVector3 const &NuP = event->GetHMISParticle(14)->fP.Vect(); TVector3 const &LeptonP = event->GetHMFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1))->fP.Vect(); - TVector3 HadronP = event->GetHMFSParticle(2212)->fP.Vect(); + // Find the highest momentum proton in the event between 450 and 1200 MeV with theta_p < 70 + TLorentzVector Pnu = event->GetNeutrinoIn()->fP; + int HMFSProton = 0; + double HighestMomentum = 0.0; + // Get the stack of protons + std::vector Protons = event->GetAllFSProton(); + for (size_t i = 0; i < Protons.size(); ++i) { + if (Protons[i]->p() > 450 && + Protons[i]->p() < 1200 && + Protons[i]->P3().Angle(Pnu.Vect()) < (M_PI/180.0)*70.0 && + Protons[i]->p() > HighestMomentum) { + HMFSProton = i; + } + } + // Now get the proton + TLorentzVector Pprot = Protons[HMFSProton]->fP; + // Get highest momentum proton in allowed proton range + TVector3 HadronP = Pprot.Vect(); // If we don't have a CC0pi signal definition we also add in pion momentum if (!Is0pi) { if (event->NumFSParticle(PhysConst::pdg_pions) == 0) { return -9999; } // Count up pion momentum TLorentzVector pp = event->GetHMFSParticle(PhysConst::pdg_pions)->fP; HadronP += pp.Vect(); } return GetDeltaPT(LeptonP, HadronP, NuP).Mag(); } double FitUtils::Get_STV_dphit(FitEvent *event, int ISPDG, bool Is0pi) { // Check that the neutrino exists if (event->NumISParticle(ISPDG) == 0) { return -9999; } // Return 0 if the proton or muon are missing if (event->NumFSParticle(2212) == 0 || event->NumFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1)) == 0) { return -9999; } // Now get the TVector3s for each particle TVector3 const &NuP = event->GetHMISParticle(ISPDG)->fP.Vect(); TVector3 const &LeptonP = event->GetHMFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1))->fP.Vect(); - TVector3 HadronP = event->GetHMFSParticle(2212)->fP.Vect(); + // Find the highest momentum proton in the event between 450 and 1200 MeV with theta_p < 70 + TLorentzVector Pnu = event->GetNeutrinoIn()->fP; + int HMFSProton = 0; + double HighestMomentum = 0.0; + // Get the stack of protons + std::vector Protons = event->GetAllFSProton(); + for (size_t i = 0; i < Protons.size(); ++i) { + if (Protons[i]->p() > 450 && + Protons[i]->p() < 1200 && + Protons[i]->P3().Angle(Pnu.Vect()) < (M_PI/180.0)*70.0 && + Protons[i]->p() > HighestMomentum) { + HMFSProton = i; + } + } + // Now get the proton + TLorentzVector Pprot = Protons[HMFSProton]->fP; + // Get highest momentum proton in allowed proton range + TVector3 HadronP = Pprot.Vect(); if (!Is0pi) { if (event->NumFSParticle(PhysConst::pdg_pions) == 0) { return -9999; } TLorentzVector pp = event->GetHMFSParticle(PhysConst::pdg_pions)->fP; HadronP += pp.Vect(); } return GetDeltaPhiT(LeptonP, HadronP, NuP); } + double FitUtils::Get_STV_dalphat(FitEvent *event, int ISPDG, bool Is0pi) { // Check that the neutrino exists if (event->NumISParticle(ISPDG) == 0) { return -9999; } // Return 0 if the proton or muon are missing if (event->NumFSParticle(2212) == 0 || event->NumFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1)) == 0) { return -9999; } // Now get the TVector3s for each particle TVector3 const &NuP = event->GetHMISParticle(ISPDG)->fP.Vect(); TVector3 const &LeptonP = event->GetHMFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1))->fP.Vect(); - TVector3 HadronP = event->GetHMFSParticle(2212)->fP.Vect(); + // Find the highest momentum proton in the event between 450 and 1200 MeV with theta_p < 70 + TLorentzVector Pnu = event->GetNeutrinoIn()->fP; + int HMFSProton = 0; + double HighestMomentum = 0.0; + // Get the stack of protons + std::vector Protons = event->GetAllFSProton(); + for (size_t i = 0; i < Protons.size(); ++i) { + if (Protons[i]->p() > 450 && + Protons[i]->p() < 1200 && + Protons[i]->P3().Angle(Pnu.Vect()) < (M_PI/180.0)*70.0 && + Protons[i]->p() > HighestMomentum) { + HMFSProton = i; + } + } + // Now get the proton + TLorentzVector Pprot = Protons[HMFSProton]->fP; + // Get highest momentum proton in allowed proton range + TVector3 HadronP = Pprot.Vect(); if (!Is0pi) { if (event->NumFSParticle(PhysConst::pdg_pions) == 0) { return -9999; } TLorentzVector pp = event->GetHMFSParticle(PhysConst::pdg_pions)->fP; HadronP += pp.Vect(); } return GetDeltaAlphaT(LeptonP, HadronP, NuP); } // As defined in PhysRevC.95.065501 // Using prescription from arXiv 1805.05486 // Returns in GeV double FitUtils::Get_pn_reco_C(FitEvent *event, int ISPDG, bool Is0pi) { const double mn = PhysConst::mass_neutron; // neutron mass const double mp = PhysConst::mass_proton; // proton mass // Check that the neutrino exists if (event->NumISParticle(ISPDG) == 0) { return -9999; } // Return 0 if the proton or muon are missing if (event->NumFSParticle(2212) == 0 || event->NumFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1)) == 0) { return -9999; } // Now get the TVector3s for each particle TVector3 const &NuP = event->GetHMISParticle(14)->fP.Vect(); TVector3 const &LeptonP = event->GetHMFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1))->fP.Vect(); - TVector3 HadronP = event->GetHMFSParticle(2212)->fP.Vect(); + + // Find the highest momentum proton in the event between 450 and 1200 MeV with theta_p < 70 + TLorentzVector Pnu = event->GetNeutrinoIn()->fP; + int HMFSProton = 0; + double HighestMomentum = 0.0; + // Get the stack of protons + std::vector Protons = event->GetAllFSProton(); + for (size_t i = 0; i < Protons.size(); ++i) { + if (Protons[i]->p() > 450 && + Protons[i]->p() < 1200 && + Protons[i]->P3().Angle(Pnu.Vect()) < (M_PI/180.0)*70.0 && + Protons[i]->p() > HighestMomentum) { + HMFSProton = i; + } + } + // Now get the proton + TLorentzVector Pprot = Protons[HMFSProton]->fP; + // Get highest momentum proton in allowed proton range + TVector3 HadronP = Pprot.Vect(); double const el = event->GetHMFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1))->E()/1000.; - double const eh = event->GetHMFSParticle(2212)->E()/1000.; + double const eh = Pprot.E(); if (!Is0pi) { if (event->NumFSParticle(PhysConst::pdg_pions) == 0) { return -9999; } - //TLorentzVector pp = event->GetHMFSParticle(PhysConst::pdg_pions)->fP; - //HadronP += pp.Vect(); + TLorentzVector pp = event->GetHMFSParticle(PhysConst::pdg_pions)->fP; + HadronP += pp.Vect(); } TVector3 dpt = GetDeltaPT(LeptonP, HadronP, NuP); double dptMag = dpt.Mag()/1000.; double ma = 6*mn + 6*mp - 0.09216; // target mass (E is from PhysRevC.95.065501) - double map = ma - mn + 0.02713; // reminant mass + double map = ma - mn + 0.02713; // remnant mass double pmul = LeptonP.Dot(NuP.Unit())/1000.; double phl = HadronP.Dot(NuP.Unit())/1000.; //double pmul = GetVectorInTPlane(LeptonP, dpt).Mag()/1000.; //double phl = GetVectorInTPlane(HadronP, dpt).Mag()/1000.; double R = ma + pmul + phl - el - eh; double dpl = 0.5*R - (map*map + dptMag*dptMag)/(2*R); //double dpl = ((R*R)-(dptMag*dptMag)-(map*map))/(2*R); // as in in PhysRevC.95.065501 - gives same result double pn_reco = sqrt((dptMag*dptMag) + (dpl*dpl)); //std::cout << "Diagnostics: " << std::endl; //std::cout << "mn: " << mn << std::endl; //std::cout << "ma: " << ma << std::endl; //std::cout << "map: " << map << std::endl; //std::cout << "pmu: " << LeptonP.Mag()/1000. << std::endl; //std::cout << "ph: " << HadronP.Mag()/1000. << std::endl; //std::cout << "pmul: " << pmul << std::endl; //std::cout << "phl: " << phl << std::endl; //std::cout << "el: " << el << std::endl; //std::cout << "eh: " << eh << std::endl; //std::cout << "R: " << R << std::endl; //std::cout << "dptMag: " << dptMag << std::endl; //std::cout << "dpl: " << dpl << std::endl; //std::cout << "pn_reco: " << pn_reco << std::endl; return pn_reco; } // Get Cos theta with Adler angles double FitUtils::CosThAdler(TLorentzVector Pnu, TLorentzVector Pmu, TLorentzVector Ppi, TLorentzVector Pprot) { // Get the "resonance" lorentz vector (pion proton system) TLorentzVector Pres = Pprot + Ppi; // Boost the particles into the resonance rest frame so we can define the x,y,z axis Pnu.Boost(Pres.BoostVector()); Pmu.Boost(-Pres.BoostVector()); Ppi.Boost(-Pres.BoostVector()); // The z-axis is defined as the axis of three-momentum transfer, \vec{k} // Also unit normalise the axis TVector3 zAxis = (Pnu.Vect()-Pmu.Vect())*(1.0/((Pnu.Vect()-Pmu.Vect()).Mag())); // Then the angle between the pion in the "resonance" rest-frame and the z-axis is the theta Adler angle double costhAdler = cos(Ppi.Vect().Angle(zAxis)); return costhAdler; } // Get phi with Adler angles, a bit more complicated... double FitUtils::PhiAdler(TLorentzVector Pnu, TLorentzVector Pmu, TLorentzVector Ppi, TLorentzVector Pprot) { // Get the "resonance" lorentz vector (pion proton system) TLorentzVector Pres = Pprot + Ppi; // Boost the particles into the resonance rest frame so we can define the x,y,z axis Pnu.Boost(Pres.BoostVector()); Pmu.Boost(-Pres.BoostVector()); Ppi.Boost(-Pres.BoostVector()); // The z-axis is defined as the axis of three-momentum transfer, \vec{k} // Also unit normalise the axis TVector3 zAxis = (Pnu.Vect()-Pmu.Vect())*(1.0/((Pnu.Vect()-Pmu.Vect()).Mag())); // The y-axis is then defined perpendicular to z and muon momentum in the resonance frame TVector3 yAxis = Pnu.Vect().Cross(Pmu.Vect()); yAxis *= 1.0/double(yAxis.Mag()); // And the x-axis is then simply perpendicular to z and x TVector3 xAxis = yAxis.Cross(zAxis); xAxis *= 1.0/double(xAxis.Mag()); double x = Ppi.Vect().Dot(xAxis); double y = Ppi.Vect().Dot(yAxis); //double z = Ppi.Vect().Dot(zAxis); double newphi = atan2(y, x); // Old silly method before atan2 /* // Then finally construct phi as the angle between pion projection and x axis // Get the project of the pion momentum on to the zaxis TVector3 PiVectZ = zAxis*Ppi.Vect().Dot(zAxis); // The subtract the projection off the pion vector to get to get the plane TVector3 PiPlane = Ppi.Vect() - PiVectZ; double phi = -999.99; if (PiPlane.Y() > 0) { phi = (180./M_PI)*PiPlane.Angle(xAxis); } else if (PiPlane.Y() < 0) { phi = (180./M_PI)*(2*M_PI-PiPlane.Angle(xAxis)); } else if (PiPlane.Y() == 0) { TRandom3 rand; double randNo = rand.Rndm(); if (randNo > 0.5) { phi = (180./M_PI)*PiPlane.Angle(xAxis); } else { phi = (180./M_PI)*(2*M_PI-PiPlane.Angle(xAxis)); } } */ return newphi; } //******************************************************************** double FitUtils::ppInfK(TLorentzVector pmu, double costh, double binding, bool neutrino) { //******************************************************************** // Convert all values to GeV //const double V = binding / 1000.; // binding potential //const double mn = PhysConst::mass_neutron; // neutron mass const double mp = PhysConst::mass_proton; // proton mass double el = pmu.E() / 1000.; //double pl = (pmu.Vect().Mag()) / 1000.; // momentum of lepton double enu = EnuQErec(pmu, costh, binding, neutrino); double ep_inf = enu - el + mp; double pp_inf = sqrt(ep_inf * ep_inf - mp * mp); return pp_inf; }; //******************************************************************** TVector3 FitUtils::tppInfK(TLorentzVector pmu, double costh, double binding, bool neutrino) { //******************************************************************** // Convert all values to GeV //const double V = binding / 1000.; // binding potential //const double mn = PhysConst::mass_neutron; // neutron mass //const double mp = PhysConst::mass_proton; // proton mass double pl_x = pmu.X() / 1000.; double pl_y = pmu.Y() / 1000.; double pl_z= pmu.Z() / 1000.; double enu = EnuQErec(pmu, costh, binding, neutrino); TVector3 tpp_inf(-pl_x, -pl_y, -pl_z+enu); return tpp_inf; }; //******************************************************************** double FitUtils::cthpInfK(TLorentzVector pmu, double costh, double binding, bool neutrino) { //******************************************************************** // Convert all values to GeV //const double V = binding / 1000.; // binding potential //const double mn = PhysConst::mass_neutron; // neutron mass const double mp = PhysConst::mass_proton; // proton mass double el = pmu.E() / 1000.; double pl = (pmu.Vect().Mag()) / 1000.; // momentum of lepton double enu = EnuQErec(pmu, costh, binding, neutrino); double ep_inf = enu - el + mp; double pp_inf = sqrt(ep_inf * ep_inf - mp * mp); double cth_inf = (enu*enu + pp_inf*pp_inf - pl*pl)/(2*enu*pp_inf); return cth_inf; };