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 <dirent.h>
// linux
#include <dlfcn.h>
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<std::string> 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<std::string> 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<DSF_NSamples_ptr>(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<DSF_GetSampleName_ptr>(
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<DSF_GetSample_ptr>(
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<DSF_DestroySample_ptr>(
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<std::string, std::vector<std::string> > ManifestSamples;
for (std::map<std::string, std::pair<std::string, int> >::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<std::string>();
}
ManifestSamples[smp_it->second.first].push_back(smp_it->first);
}
QLOG(FIT, "Dynamic sample manifest: ");
for (std::map<std::string, std::vector<std::string> >::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<std::string, int> 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 <http://www.gnu.org/licenses/>.
################################################################################
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 <http://www.gnu.org/licenses/>.
+*******************************************************************************/
+
+/*
+ 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<std::string> 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 <http://www.gnu.org/licenses/>.
+*******************************************************************************/
+
+#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 <http://www.gnu.org/licenses/>.
*******************************************************************************/
// 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<std::string> 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<double> CrossSection;
std::vector<double> Error;
std::vector<double> 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<std::string> 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<FitParticle*> 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 <http://www.gnu.org/licenses/>.
*******************************************************************************/
#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 <http://www.gnu.org/licenses/>.
*******************************************************************************/
#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<EventRecord*>(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<GHepRecord*>(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<FitParticle*> 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 <http://www.gnu.org/licenses/>.
*******************************************************************************/
#include "FitUtils.h"
/*
MISC Functions
*/
//********************************************************************
double *FitUtils::GetArrayFromMap(std::vector<std::string> invals,
std::map<std::string, double> 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<FitParticle *> 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<FitParticle *> 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<FitParticle*> 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<FitParticle*> 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<FitParticle*> 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<FitParticle*> 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;
};