Page MenuHomeHEPForge
Files F7879237

No OneTemporary
Actions

View Options

diff --git a/src/FCN/SampleList.cxx b/src/FCN/SampleList.cxx
index 3db2359..5fb0a68 100644
--- a/src/FCN/SampleList.cxx
+++ b/src/FCN/SampleList.cxx
@@ -1,1016 +1,1023 @@
#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"
// 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"
// 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 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_nu.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_1Dth_antinu.h"
#include "MINERvA_CCCOHPI_XSec_1Dth_nu.h"
#include "MINERvA_CCCOHPI_XSec_1DQ2_nu.h"
#include "MINERvA_CC0pi_XSec_1DQ2_TgtRatio_nu.h"
#include "MINERvA_CC0pi_XSec_1DQ2_Tgt_nu.h"
#include "MINERvA_CC1pip_XSec_1D_2017Update.h"
-
+#include "MINERvA_CC0pi_XSec_2Dptpx_nu.h"
+#include "MINERvA_CC0pi_XSec_2Dptpx_antinu.h"
#endif
#ifndef __NO_T2K__
// T2K CC0pi
#include "T2K_CC0pi_XSec_2DPcos_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"
// T2K STV CC0pi
#include "T2K_CC0pinp_STV_XSec_1Ddpt_nu.h"
#include "T2K_CC0pi_XSec_2DPcos_nu_nonuniform.h"
#endif
#ifndef __NO_SciBooNE__
// SciBooNE COH studies
#include "SciBooNE_CCCOH_1TRK_1DQ2_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_MuPiVA_1DQ2_nu.h"
#include "SciBooNE_CCCOH_MuPr_1DQ2_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 "FitWeight.h"
#include "NuisConfig.h"
#include "NuisKey.h"
//! 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.AddS("name", name);
samplekey.AddS("input", file);
samplekey.AddS("type", type);
return CreateSample(samplekey);
}
MeasurementBase* CreateSample(nuiskey samplekey) {
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_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_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));
/*
CC1pi0
*/
// 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));
/*
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));
/*
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(file, rw, type, fkdt));
} else if (!name.compare("T2K_CC1pip_CH_XSec_1Dppi_nu")) {
return (new T2K_CC1pip_CH_XSec_1Dppi_nu(file, rw, type, fkdt));
} 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));
// 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_MuPr_1DQ2_nu")) {
return (new SciBooNE_CCCOH_MuPr_1DQ2_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_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_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 {
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/FitBase/ParamPull.h b/src/FitBase/ParamPull.h
index 7da9eac..aa99977 100644
--- a/src/FitBase/ParamPull.h
+++ b/src/FitBase/ParamPull.h
@@ -1,193 +1,194 @@
// 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 PARAM_PULL_H_SEEN
#define PARAM_PULL_H_SEEN
/*!
* \addtogroup FitBase
* @{
*/
#include <stdlib.h>
#include <numeric>
#include <math.h>
#include <iostream>
#include <sstream>
#include <fstream>
#include <list>
#include <vector>
// ROOT includes
#include <TROOT.h>
#include <TH1D.h>
#include <TTree.h>
#include <TFile.h>
#include <TMatrixDSym.h>
#include <TDecompSVD.h>
#include <TVectorD.h>
+#include "TRandom3.h"
// Fit Includes
#include "PlotUtils.h"
#include "StatUtils.h"
#include "FitWeight.h"
#include "FitLogger.h"
#include "EventManager.h"
#include "TVector.h"
using namespace std;
//! Enums to allow Non Guassian Pulls in Future
enum FitPullTypes {
kUnknownPull = -1,
kNoPull = 0,
kGausPull = 1
};
enum FitThrowTypes {
kUnknownThrow = -1,
kNoThrow = 0,
kGausThrow = 1,
kFlatThrow = 2
};
//! Used to produce gaussian penalty terms in the fit.
class ParamPull {
public:
//! Default Constructor
ParamPull(std::string name, std::string inputfile, std::string type, std::string dials="");
//! Default destructor
virtual ~ParamPull(void) {};
// Set dial types (DEFAULT,ABS,FRAC)
void SetType(std::string type);
// Setup Histogram inputs (from previous fit file, or ROOT file)
void SetupHistograms(std::string input);
// Read a previous NUISANCE file
void ReadFitFile(std::string input);
// Read a ROOT file with any histograms in
void ReadRootFile(std::string input);
// Read Text file
void ReadVectFile(std::string input);
// Read a single dial central value and error
void ReadDialInput(std::string input);
//! Reconfigure function reads in current RW engine dials and sets their value to MC
void Reconfigure(void);
//! Get likelihood given the current values
double GetLikelihood(void);
//! Get NDOF if used in likelihoods
int GetNDOF(void);
// Get Covariance Matrices as plots
TH2D GetCovar(void);
TH2D GetFullCovar(void);
TH2D GetDecompCovar(void);
// Get Covariance Matrices
inline TMatrixDSym GetCovarMatrix (void) const { return *fInvCovar; };
inline TMatrixDSym GetFullCovarMatrix (void) const { return *fCovar; };
inline TMatrixDSym GetDecompCovarMatrix (void) const { return *fDecomp; };
//! Save the histograms
void Write(std::string writeopt="");
//! Throw the dial values using the current covariance. Useful for parameter throws.
void ThrowCovariance(void);
//! Compare dials to RW
bool CheckDialsValid(void);
//! Reset toy data back to original data
void ResetToy(void);
//! Read fake data from MC
void SetFakeData(std::string fakeinput);
//! Reset fake data back to original data (before fake data or throws)
void RemoveFakeData();
// Get Functions
inline std::string GetName (void) const { return fName; };
inline std::string GetInput (void) const { return fInput; };
inline std::string GetType (void) const { return fType; };
inline std::string GetFileType (void) const { return fFileType; };
inline std::string GetDialOptions (void) const { return fDialOptions; };
std::map<std::string, int> GetAllDials();
inline TH1D GetDataHist (void) const { return *fDataHist; };
inline TH1D GetDataTrue (void) const { return *fDataTrue; };
inline TH1D GetDataOrig (void) const { return *fDataOrig; };
inline TH1D GetMCHist (void) const { return *fMCHist; };
inline TH1D GetMaxHist (void) const { return *fMaxHist; };
inline TH1D GetMinHist (void) const { return *fMinHist; };
inline TH1I GetDialTypes (void) const { return *fTypeHist; };
inline TH1D GetLimitHist (void) const { return *fLimitHist; };
private:
TH1D RemoveBinsNotInString(TH1D hist, std::string mystr);
TH1I RemoveBinsNotInString(TH1I hist, std::string mystr);
std::string fName; //!< Pull Name
std::string fInput; //!< Pull input string
std::string fType; //!< Pull options type
std::string fFileType; //!< Pull input file types
std::string fPlotTitles; //! Axis format
std::string fDialOptions; //!< Dial handling options
std::string fDialSelection; //!< Dial Selection
TH1D* fMCHist; //!< Current MC Histogram
TH1D* fDataHist; //!< Current data Histogram
TH1D* fDataTrue; //!< True Data (before histogram throws)
TH1D* fDataOrig; //!< Orig Data (without toys or fake data)
TH1D* fMaxHist; //!< Maximum limit on MC/Data
TH1D* fMinHist; //!< Maximum limit on MC/Data
TH1I* fTypeHist; //!< Dial Types
int fCalcType; //!< Method to calculate likelihood
int fThrowType; //!< Method to calculate throws
TMatrixDSym* fCovar; //!< Covariance
TMatrixDSym* fInvCovar; //!< Inverted Covariance
TMatrixDSym* fDecomp; //!< Decomposition
TH1D* fLimitHist;
};
// Class TypeDefs
typedef std::list<ParamPull*>::const_iterator PullListConstIter;
typedef std::list<ParamPull*>::iterator PullListIter;
typedef std::vector<ParamPull*>::const_iterator PullVectConstIter;
typedef std::vector<ParamPull*>::iterator PullVectIter;
/*! @} */
#endif
diff --git a/src/MINERvA/CMakeLists.txt b/src/MINERvA/CMakeLists.txt
index d41a492..18ee3ba 100644
--- a/src/MINERvA/CMakeLists.txt
+++ b/src/MINERvA/CMakeLists.txt
@@ -1,156 +1,164 @@
# 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_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_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_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_CC1pip_XSec_1D_2017Update.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
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_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_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_1DEnu_antinu.h
MINERvA_CCCOHPI_XSec_1DEpi_antinu.h
MINERvA_CCCOHPI_XSec_1Dth_antinu.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_2Dptpx_antinu.cxx b/src/MINERvA/MINERvA_CC0pi_XSec_2Dptpx_antinu.cxx
new file mode 100755
index 0000000..36d7e17
--- /dev/null
+++ b/src/MINERvA/MINERvA_CC0pi_XSec_2Dptpx_antinu.cxx
@@ -0,0 +1,112 @@
+//Adrian Orea
+//I used the file MINERvA_CCinc_XSec_2DEavq3_nu.cxx as a template
+//Also, I am fully aware of the naming typo (should be ptpz), but Everything is already named the same way so...
+
+//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/>.
+*******************************************************************************/
+/*
+ Author : Adrian Orea
+*/
+
+#include "MINERvA_SignalDef.h"
+#include "MINERvA_CC0pi_XSec_2Dptpx_antinu.h"
+
+//********************************************************************
+MINERvA_CC0pi_XSec_2Dptpx_antinu::MINERvA_CC0pi_XSec_2Dptpx_antinu(nuiskey samplekey) {
+//********************************************************************
+
+ // Sample overview ---------------------------------------------------
+ std::string descrip = "MINERvA_CC0pi_XSec_2Dptpx_antinu sample. \n" \
+ "Target: CH \n" \
+ "Flux: MINERvA Medium Energy FHC numu \n" \
+ "Signal: CC-0pi \n";
+
+ // Setup common settings
+ fSettings = LoadSampleSettings(samplekey);
+ fSettings.SetDescription(descrip);
+ fSettings.SetXTitle("p_{t} (GeV)");
+ fSettings.SetYTitle("p_{z} (GeV)");
+ fSettings.SetZTitle("d^{2}#sigma/dP_{t}dP_{z} (cm^{2}/GeV^{2}/nucleon)");
+ fSettings.SetAllowedTypes("FIX,FREE,SHAPE/FULL,DIAG/MASK", "FIX/FULL");
+ fSettings.SetEnuRange(0.0, 100.0);
+ fSettings.DefineAllowedTargets("C,H");
+
+ // CCQELike plot information
+ fSettings.SetTitle("MINERvA_CC0pi_XSec_2Dptpx_antinu");
+
+ fSettings.SetDataInput( FitPar::GetDataBase() + "MINERvA/CC0pi/data_2D.txt" );
+ fSettings.SetErrorInput( FitPar::GetDataBase() + "MINERvA/CC0pi/error_2D.txt" );
+ fSettings.SetCovarInput( FitPar::GetDataBase() + "MINERvA/CC0pi/covar_2D.txt" );
+ fSettings.SetMapInput( FitPar::GetDataBase() + "MINERvA/CC0pi/map_2D.txt" );
+ fSettings.DefineAllowedSpecies("antinumu");
+ FinaliseSampleSettings();
+
+ // Scaling Setup ---------------------------------------------------
+ // ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
+ fScaleFactor = (GetEventHistogram()->Integral("width") * 1E-38 / (fNEvents + 0.)) / this->TotalIntegratedFlux();
+
+ // Plot Setup -------------------------------------------------------
+ Double_t P_t[7] = {0,0.15,0.25,0.4,0.7,1.0,1.5};
+ Double_t P_z[12] = {1.5,2.0,2.5,3.0,3.5,4.0,4.5,5.0,6.0,8.0,10.0,15.0};
+
+ CreateDataHistogram(7, P_t, 12, P_z);
+ SetDataValuesFromTextFile( fSettings.GetDataInput() );
+ ScaleData(1E-41);
+
+ SetMapValuesFromText( fSettings.GetMapInput() );
+
+ SetDataErrorsFromTextFile( fSettings.GetErrorInput() );
+ ScaleDataErrors(1E-41);
+
+ //SetCholDecompFromTextFile( fSettings.GetCovarInput() );
+ //ScaleCovar(1E-16);
+
+ //StatUtils::SetDataErrorFromCov(fDataHist, fFullCovar, fMapHist, 1E-38);
+
+ // Final setup ---------------------------------------------------
+ FinaliseMeasurement();
+
+};
+
+//********************************************************************
+void MINERvA_CC0pi_XSec_2Dptpx_antinu::FillEventVariables(FitEvent *event) {
+//********************************************************************
+
+// Checking to see if there is a Muon
+ if (event->NumFSParticle(-13) == 0)
+ return;
+
+ TLorentzVector Pmu = event->GetHMFSParticle(-13)->fP; //I Added this part
+ Double_t px = Pmu.X()/1000;
+ Double_t py = Pmu.Y()/1000;
+ Double_t pz = Pmu.Z()/1000;
+ Double_t pt = sqrt(px*px+py*py);
+
+// Set Hist Variables
+ fYVar = pz;
+ fXVar = pt;
+
+ return;
+};
+
+//********************************************************************
+bool MINERvA_CC0pi_XSec_2Dptpx_antinu::isSignal(FitEvent *event) {
+//********************************************************************
+ return (SignalDef::isCC0pi_anti_MINERvAPTPZ(event, -14, EnuMin, EnuMax)/* && SignalDef::HasProtonMomBelowThreshold(event, 120)*/);
+};
diff --git a/src/MINERvA/MINERvA_CC0pi_XSec_2Dptpx_antinu.h b/src/MINERvA/MINERvA_CC0pi_XSec_2Dptpx_antinu.h
new file mode 100755
index 0000000..1cbaac6
--- /dev/null
+++ b/src/MINERvA/MINERvA_CC0pi_XSec_2Dptpx_antinu.h
@@ -0,0 +1,69 @@
+//Adrian Orea
+//NEED TO MODIFY
+//I used the file MINERvA_CCinc_XSec_2DEavq3_nu.h as a template
+
+// 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_2DPTPX_ANTINU_H_SEEN
+#define MINERVA_CC0PI_XSEC_2DPTPX_ANTINU_H_SEEN
+
+#include "Measurement2D.h"
+
+//********************************************************************
+class MINERvA_CC0pi_XSec_2Dptpx_antinu : public Measurement2D {
+//********************************************************************
+
+ public:
+
+ // Constructor
+ MINERvA_CC0pi_XSec_2Dptpx_antinu(nuiskey samplekey);
+
+ // Destructor
+ virtual ~MINERvA_CC0pi_XSec_2Dptpx_antinu() {
+
+ // 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:
+
+ // Cuts
+};
+
+#endif
diff --git a/src/MINERvA/MINERvA_CC0pi_XSec_2Dptpx_nu.cxx b/src/MINERvA/MINERvA_CC0pi_XSec_2Dptpx_nu.cxx
new file mode 100755
index 0000000..0d88e2c
--- /dev/null
+++ b/src/MINERvA/MINERvA_CC0pi_XSec_2Dptpx_nu.cxx
@@ -0,0 +1,110 @@
+//Adrian Orea
+//I used the file MINERvA_CCinc_XSec_2DEavq3_nu.cxx as a template
+//Also, I am fully aware of the naming typo (should be ptpz), but Everything is already named the same way so...
+
+//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/>.
+*******************************************************************************/
+
+/*
+ Author : Adrian Orea
+ */
+
+
+#include "MINERvA_SignalDef.h"
+#include "MINERvA_CC0pi_XSec_2Dptpx_nu.h"
+
+//********************************************************************
+MINERvA_CC0pi_XSec_2Dptpx_nu::MINERvA_CC0pi_XSec_2Dptpx_nu(nuiskey samplekey) {
+//********************************************************************
+
+ // Sample overview ---------------------------------------------------
+ std::string descrip = "MINERvA_CC0pi_XSec_2Dptpx_nu sample. \n" \
+ "Target: CH \n" \
+ "Flux: MINERvA Medium Energy FHC numu \n" \
+ "Signal: CC-0pi \n";
+
+ // Setup common settings
+ fSettings = LoadSampleSettings(samplekey);
+ fSettings.SetDescription(descrip);
+ fSettings.SetXTitle("p_{z} (GeV)");
+ fSettings.SetYTitle("p_{t} (GeV)");
+ fSettings.SetZTitle("d^{2}#sigma/dP_{t}dP_{z} (cm^{2}/GeV^{2}/nucleon)");
+ fSettings.SetAllowedTypes("FIX,FREE,SHAPE/FULL,DIAG/MASK", "FIX/FULL");
+ fSettings.SetEnuRange(0.0, 100.0);
+ fSettings.DefineAllowedTargets("C,H");
+
+ // CCQELike plot information
+ fSettings.SetTitle("MINERvA_CC0pi_XSec_2Dptpx_nu");
+
+ fSettings.SetDataInput( FitPar::GetDataBase() + "MINERvA/CC0pi/CCQE_Neutrino_CV_Result.txt" );
+ fSettings.SetCovarInput( FitPar::GetDataBase() + "MINERvA/CC0pi/covar_2D.txt" );
+ fSettings.SetMapInput( FitPar::GetDataBase() + "MINERvA/CC0pi/map_2D.txt" );
+ fSettings.DefineAllowedSpecies("numu");
+ FinaliseSampleSettings();
+
+ // Scaling Setup ---------------------------------------------------
+ // ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
+ fScaleFactor = (GetEventHistogram()->Integral("width") * 1E-37 / (fNEvents + 0.)) / this->TotalIntegratedFlux();
+
+ // Plot Setup -------------------------------------------------------
+ Double_t P_t[14] = {0,0.075,0.15,0.25,0.325,0.4,0.475,0.55,0.7,0.85,1.0,1.25,1.5,2.5};
+ Double_t P_z[13] = {1.5,2.0,2.5,3.0,3.5,4.0,4.5,5.0,6.0,8.0,10.0,15.0,20.0};
+ CreateDataHistogram(13, P_z, 14, P_t);
+ SetDataValuesFromTextFile( fSettings.GetDataInput() );
+ ScaleData(1E-41);
+
+ SetMapValuesFromText( fSettings.GetMapInput() );
+
+ SetCholDecompFromTextFile( fSettings.GetCovarInput() );
+ ScaleCovar(1E-16);
+
+ StatUtils::SetDataErrorFromCov(fDataHist, fFullCovar, fMapHist, 1E-38);
+
+ // Final setup ---------------------------------------------------
+ FinaliseMeasurement();
+
+};
+
+//********************************************************************
+void MINERvA_CC0pi_XSec_2Dptpx_nu::FillEventVariables(FitEvent *event) {
+//********************************************************************
+
+// Checking to see if there is a Muon
+ if (event->NumFSParticle(13) == 0)
+ return;
+
+ TLorentzVector Pmu = event->GetHMFSParticle(13)->fP; //I Added this part
+ Double_t px = Pmu.X()/1000;
+ Double_t py = Pmu.Y()/1000;
+ Double_t pz = Pmu.Z()/1000;
+ Double_t pt = sqrt(px*px+py*py);
+
+// Set Hist Variables
+ fYVar = pt;
+ fXVar = pz;
+
+ return;
+};
+
+//********************************************************************
+bool MINERvA_CC0pi_XSec_2Dptpx_nu::isSignal(FitEvent *event) {
+//********************************************************************
+ return SignalDef::isCC0pi_MINERvAPTPZ(event, 14, EnuMin, EnuMax);
+
+};
diff --git a/src/MINERvA/MINERvA_CC0pi_XSec_2Dptpx_nu.h b/src/MINERvA/MINERvA_CC0pi_XSec_2Dptpx_nu.h
new file mode 100755
index 0000000..40bb3d0
--- /dev/null
+++ b/src/MINERvA/MINERvA_CC0pi_XSec_2Dptpx_nu.h
@@ -0,0 +1,69 @@
+//Adrian Orea
+//NEED TO MODIFY
+//I used the file MINERvA_CCinc_XSec_2DEavq3_nu.h as a template
+
+// 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_2DPTPX_NU_H_SEEN
+#define MINERVA_CC0PI_XSEC_2DPTPX_NU_H_SEEN
+
+#include "Measurement2D.h"
+
+//********************************************************************
+class MINERvA_CC0pi_XSec_2Dptpx_nu : public Measurement2D {
+//********************************************************************
+
+ public:
+
+ // Constructor
+ MINERvA_CC0pi_XSec_2Dptpx_nu(nuiskey samplekey);
+
+ // Destructor
+ virtual ~MINERvA_CC0pi_XSec_2Dptpx_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:
+
+ // Cuts
+};
+
+#endif
diff --git a/src/MINERvA/MINERvA_CCCOHPI_XSec_joint.cxx b/src/MINERvA/MINERvA_CCCOHPI_XSec_joint.cxx
new file mode 100644
index 0000000..5eca830
--- /dev/null
+++ b/src/MINERvA/MINERvA_CCCOHPI_XSec_joint.cxx
@@ -0,0 +1,203 @@
+// 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 "MINERvA_SignalDef.h"
+#include "MINERvA_CCCOHPI_XSec_joint.h"
+#include "MINERvA_CCCOHPI_XSec_1DEpi_nu.h"
+#include "MINERvA_CCCOHPI_XSec_1DEpi_antinu.h"
+#include "MINERvA_CCCOHPI_XSec_1Dth_nu.h"
+#include "MINERvA_CCCOHPI_XSec_1Dth_antinu.h"
+#include "MINERvA_CCCOHPI_XSec_1DEnu_nu.h"
+#include "MINERvA_CCCOHPI_XSec_1DEnu_antinu.h"
+#include "MINERvA_CCCOHPI_XSec_1DQ2_nu.h"
+#include "MINERvA_CCCOHPI_XSec_1DQ2_antinu.h"
+
+//********************************************************************
+void MINERvA_CCCOHPI_XSec_joint::SetupDataSettings(){
+//********************************************************************
+
+ // Set Distribution
+ std::string name = fSettings.GetS("name");
+ if (!name.compare("MINERvA_CCCOHPI_XSec_1DEpi_joint")) fDist = kEpi;
+ else if (!name.compare("MINERvA_CCCOHPI_XSec_1Dth_joint")) fDist = kth;
+ else if (!name.compare("MINERvA_CCCOHPI_XSec_1DQ2_joint")) fDist = kQ2;
+ else if (!name.compare("MINERvA_CCCOHPI_XSec_1DEnu_joint")) fDist = kEnu;
+
+ // Define what files to use from the dist
+ std::string datafile = "";
+ std::string covarfile = "";
+ std::string titles = "";
+ std::string distdescript = "";
+ switch(fDist){
+
+ case (kEpi):
+ datafile = "Epi_joint_data.csv";
+ covarfile = "Epi_joint_cov.csv";
+ titles = "CC-Coherent #pi;E_{#pi} Bins;d#sigma/dE_{#pi} (cm^{2}/GeV/C12)";
+ break;
+
+ case (kth):
+ datafile = "Thpi_joint_data.csv";
+ covarfile = "Thpi_joint_cov.csv";
+ titles = "CC-Coherent #pi;#theta_{#pi} Bins;d#sigma/d#theta_{#pi} (cm^{2}/C12)";
+ break;
+
+ case (kQ2):
+ datafile = "Q2_joint_data.csv";
+ covarfile = "Q2_joint_cov.csv";
+ titles ="CC-Coherent #pi;Q^{2} Bins;d#sigma/dQ^{2} (cm^{2}/GeV^{2}/C12)";
+ break;
+
+ case (kEnu):
+ datafile = "Enu_joint_data.csv";
+ covarfile = "Enu_joint_cov.csv";
+ titles ="CC-Coherent #pi;E_{#nu} Bins;#sigma(E_{#nu}) (cm^{2}/C12";
+ break;
+
+ default:
+ THROW("Unknown Analysis Distribution : " << fDist);
+ }
+
+ // Now setup each data distribution and description.
+ std::string descrip = distdescript + \
+ "Target: CH \n" \
+ "Flux: 2 Files, MINERvA FHC numu;MINERvA FHC numubar \n" \
+ "Signal: Any event with 1 electron, any nucleons, and no other FS particles \n";
+
+ fSettings.SetDescription(descrip);
+ fSettings.SetDataInput( GeneralUtils::GetTopLevelDir()+"/data/MINERvA/CCcoh/" + datafile );
+ fSettings.SetCovarInput( GeneralUtils::GetTopLevelDir()+"/data/MINERvA/CCcoh/" + covarfile );
+ fSettings.SetTitle( GeneralUtils::ParseToStr(titles,";")[0] );
+ fSettings.SetXTitle( GeneralUtils::ParseToStr(titles,";")[1] );
+ fSettings.SetYTitle( GeneralUtils::ParseToStr(titles,";")[2] );
+
+ return;
+}
+
+//********************************************************************
+MINERvA_CCCOHPI_XSec_joint::MINERvA_CCCOHPI_XSec_joint(nuiskey samplekey) {
+//********************************************************************
+
+ // Define Sample Settings common to all data distributions
+ fSettings = LoadSampleSettings(samplekey);
+ fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/FULL");
+ fSettings.SetEnuRange(0.0, 20.0);
+ fSettings.DefineAllowedTargets("C,H");
+ fSettings.DefineAllowedSpecies("numu");
+ SetupDataSettings();
+ FinaliseSampleSettings();
+
+ // Scaling Setup ---------------------------------------------------
+ fScaleFactor = 0.0;
+
+ // Plot Setup -------------------------------------------------------
+ SetDataFromTextFile( fSettings.GetDataInput() );
+ SetCovarFromMultipleTextFiles(fSettings.GetCovarInput());
+
+ // Setup Sub Measurements -------------------------------------------------------
+ SetupSubMeasurements();
+
+ // Final setup ---------------------------------------------------
+ FinaliseMeasurement();
+
+};
+
+//********************************************************************
+void MINERvA_CCCOHPI_XSec_joint::SetupSubMeasurements(){
+//********************************************************************
+
+ // Get parsed input files
+ if (fSubInFiles.size() != 2) {
+ ERR(FTL) << "MINERvA Joint requires input files in format: nu;antinu" << std::endl;
+ }
+ std::string inFileNeutrino = fSubInFiles.at(0);
+ std::string inFileAntineutrino = fSubInFiles.at(1);
+
+ // Create some config keys
+ nuiskey nukey = Config::CreateKey("sample");
+ nukey.SetS("input", inFileNeutrino);
+ nukey.SetS("type", fSettings.GetS("type"));
+
+ nuiskey antinukey = Config::CreateKey("sample");
+ antinukey.SetS("input", inFileAntineutrino);
+ antinukey.SetS("type", fSettings.GetS("type"));
+
+ // Create samples for given DIST
+ switch(fDist){
+
+ case kEpi:
+ nukey.SetS("name","MINERvA_CCCOHPI_XSec_1DEpi_nu");
+ MIN_nu = new MINERvA_CCCOHPI_XSec_1DEpi_nu(nukey);
+ antinukey.SetS("name","MINERvA_CCCOHPI_XSec_1DEpi_antinu");
+ MIN_anu = new MINERvA_CCCOHPI_XSec_1DEpi_antinu(antinukey);
+ break;
+
+ case kth:
+ nukey.SetS("name","MINERvA_CCCOHPI_XSec_1Dth_nu");
+ MIN_nu = new MINERvA_CCCOHPI_XSec_1Dth_nu(nukey);
+ antinukey.SetS("name","MINERvA_CCCOHPI_XSec_1Dth_antinu");
+ MIN_anu = new MINERvA_CCCOHPI_XSec_1Dth_antinu(antinukey);
+ break;
+
+ case kEnu:
+ nukey.SetS("name","MINERvA_CCCOHPI_XSec_1DEnu_nu");
+ MIN_nu = new MINERvA_CCCOHPI_XSec_1DEnu_nu(nukey);
+ antinukey.SetS("name","MINERvA_CCCOHPI_XSec_1DEnu_antinu");
+ MIN_anu = new MINERvA_CCCOHPI_XSec_1DEnu_antinu(antinukey);
+ break;
+
+ case kQ2:
+ nukey.SetS("name","MINERvA_CCCOHPI_XSec_1DQ2_nu");
+ MIN_nu = new MINERvA_CCCOHPI_XSec_1DQ2_nu(nukey);
+ antinukey.SetS("name","MINERvA_CCCOHPI_XSec_1DQ2_antinu");
+ MIN_anu = new MINERvA_CCCOHPI_XSec_1DQ2_antinu(antinukey);
+ break;
+ }
+
+ // Add to chain for processing
+ fSubChain.clear();
+ fSubChain.push_back(MIN_anu);
+ fSubChain.push_back(MIN_nu);
+}
+
+//********************************************************************
+void MINERvA_CCCOHPI_XSec_joint::MakePlots() {
+//********************************************************************
+
+ TH1D* MIN_nu_mc = (TH1D*) MIN_nu->GetMCHistogram();
+ TH1D* MIN_anu_mc = (TH1D*) MIN_anu->GetMCHistogram();
+ int count = 0;
+ for (int i = 0; i < MIN_nu_mc->GetNbinsX(); i++){
+ fMCHist->SetBinContent( count+1, MIN_nu_mc->GetBinContent(i+1) );
+ fMCHist->SetBinError( count+1, MIN_nu_mc->GetBinError(i+1) );
+ fMCHist->GetXaxis()->SetBinLabel(count+1, GeneralUtils::IntToStr(count+1).c_str());
+ count++;
+ }
+ for (int i = 0; i < MIN_anu_mc->GetNbinsX(); i++){
+ fMCHist->SetBinContent(count+1, MIN_anu_mc->GetBinContent(i+1) );
+ fMCHist->SetBinError(count+1, MIN_anu_mc->GetBinError(i+1) );
+ fMCHist->GetXaxis()->SetBinLabel(count+1, GeneralUtils::IntToStr(count+1).c_str());
+ count++;
+ }
+
+ return;
+}
+
+
+
diff --git a/src/MINERvA/MINERvA_CCCOHPI_XSec_joint.h b/src/MINERvA/MINERvA_CCCOHPI_XSec_joint.h
new file mode 100644
index 0000000..e64205b
--- /dev/null
+++ b/src/MINERvA/MINERvA_CCCOHPI_XSec_joint.h
@@ -0,0 +1,61 @@
+// 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_CCCOHPI_XSec_joint_H_SEEN
+#define MINERVA_CCCOHPI_XSec_joint_H_SEEN
+
+// Fit Includes
+#include "MeasurementBase.h"
+#include "JointMeas1D.h"
+#include "Measurement1D.h"
+#include "MINERvA_CCCOHPI_XSec_1DEpi_nu.h"
+#include "MINERvA_CCCOHPI_XSec_1DEpi_antinu.h"
+#include "MINERvA_CCCOHPI_XSec_1Dth_nu.h"
+#include "MINERvA_CCCOHPI_XSec_1Dth_antinu.h"
+#include "MINERvA_CCCOHPI_XSec_1DEnu_nu.h"
+#include "MINERvA_CCCOHPI_XSec_1DEnu_antinu.h"
+#include "MINERvA_CCCOHPI_XSec_1DQ2_nu.h"
+#include "MINERvA_CCCOHPI_XSec_1DQ2_antinu.h"
+
+class MINERvA_CCCOHPI_XSec_joint : public JointMeas1D {
+public:
+
+ MINERvA_CCCOHPI_XSec_joint(nuiskey samplekey);
+ virtual ~MINERvA_CCCOHPI_XSec_joint() {};
+ void MakePlots();
+ void SetupDataSettings();
+ void SetupSubMeasurements();
+
+ enum DataDistribution {
+ kEpi,
+ kth,
+ kQ2,
+ kEnu
+ } CCCOH_DataDistributions;
+
+ private:
+
+ // This is a dummy, the signal is defined separately for each sample!
+ bool isSignal(){return false;};
+ Measurement1D* MIN_nu;
+ Measurement1D* MIN_anu;
+ int fDist;
+};
+
+#endif
diff --git a/src/MINERvA/MINERvA_SignalDef.cxx b/src/MINERvA/MINERvA_SignalDef.cxx
index 8c2be8c..7d62827 100644
--- a/src/MINERvA/MINERvA_SignalDef.cxx
+++ b/src/MINERvA/MINERvA_SignalDef.cxx
@@ -1,305 +1,384 @@
// 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
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
bool HasProton = event->HasFSParticle(2212);
if (CC1pi0_anu) {
if (!HasProton) {
return true;
} else {
return false;
}
} else {
return false;
}
*/
return CC1pi0_anu;
}
+ //********************************************************************
+ 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(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(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;
+ }
}
diff --git a/src/MINERvA/MINERvA_SignalDef.h b/src/MINERvA/MINERvA_SignalDef.h
index 2b1f03e..66f2799 100644
--- a/src/MINERvA/MINERvA_SignalDef.h
+++ b/src/MINERvA/MINERvA_SignalDef.h
@@ -1,94 +1,98 @@
// 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_SIGNALDEF_H_SEEN
#define MINERVA_SIGNALDEF_H_SEEN
#include "FitEvent.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 problems:
/// 1) Should there be a pi+ only cut implemented due to Michel requirement, or
/// is pi- events filled from MC?
/// 2) There is a T_pi < 350 MeV cut coming from requiring a stopping pion so
/// the
/// Michel e is seen, this is also unclear if it's unfolded so any pion is OK
///
/// 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 = false);
bool isCC1pip_MINERvA_2017(FitEvent *event, double EnuMin, double EnuMax);
// *********************************
/// MINERvA CCNpi+/- signal definition from 2016 publication
/// Different to CC1pi+/- listed above; additional has W < 1.8 GeV
///
/// Still asks for a Michel e and still unclear if this is unfolded or not
/// Says stuff like "requirement that a Michel e isolates a subsample that is
/// more nearly a pi+ prodution", yet the signal definition is both pi+ and pi-?
///
/// One negative muon
/// At least one charged pion
/// 1.5 < Enu < 10
/// No restrictions on pi0 or other mesons or baryons
///
/// 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 = false, bool isWtrue=false);
bool isCCQEnumu_MINERvA(FitEvent *event, double EnuMin, double EnuMax,
bool fullphasespace = true);
bool isCCQEnumubar_MINERvA(FitEvent *event, double EnuMin, double EnuMax,
bool fullphasespace = true);
bool isCCincLowRecoil_MINERvA(FitEvent *event, double EnuMin, double EnuMax);
bool isCC0pi1p_MINERvA(FitEvent *event, double enumin, double enumax);
bool isCC1pi0_MINERvA_2015(FitEvent *event, double EnuMin, double EnuMax);
bool isCC1pi0_MINERvA_2016(FitEvent *event, double EnuMin, double EnuMax);
+
+ bool isCC0pi_MINERvAPTPZ(FitEvent *event, int nuPDG, double EnuMin = 0, double EnuMax = 0);
+ bool isCC0pi_anti_MINERvAPTPZ(FitEvent* event, int nuPDG, double EnuMin = 0, double EnuMax = 0);
+
}
#endif
diff --git a/src/Routines/SystematicRoutines.cxx b/src/Routines/SystematicRoutines.cxx
index d33b35e..1afc9bd 100755
--- a/src/Routines/SystematicRoutines.cxx
+++ b/src/Routines/SystematicRoutines.cxx
@@ -1,1517 +1,1517 @@
// 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 "SystematicRoutines.h"
void SystematicRoutines::Init(){
fInputFile = "";
fInputRootFile = NULL;
fOutputFile = "";
fOutputRootFile = NULL;
fCovar = fCovarFree = NULL;
fCorrel = fCorrelFree = NULL;
fDecomp = fDecompFree = NULL;
fStrategy = "ErrorBands";
fRoutines.clear();
fRoutines.push_back("ErrorBands");
fCardFile = "";
fFakeDataInput = "";
fSampleFCN = NULL;
fAllowedRoutines = ("ErrorBands,PlotLimits");
};
SystematicRoutines::~SystematicRoutines(){
};
SystematicRoutines::SystematicRoutines(int argc, char* argv[]){
// Initialise Defaults
Init();
nuisconfig configuration = Config::Get();
// Default containers
std::string cardfile = "";
std::string maxevents = "-1";
int errorcount = 0;
int verbocount = 0;
std::vector<std::string> xmlcmds;
std::vector<std::string> configargs;
fNThrows = 250;
fStartThrows = 0;
fThrowString = "";
// Make easier to handle arguments.
std::vector<std::string> args = GeneralUtils::LoadCharToVectStr(argc, argv);
ParserUtils::ParseArgument(args, "-c", fCardFile, true);
ParserUtils::ParseArgument(args, "-o", fOutputFile, false, false);
ParserUtils::ParseArgument(args, "-n", maxevents, false, false);
ParserUtils::ParseArgument(args, "-f", fStrategy, false, false);
ParserUtils::ParseArgument(args, "-d", fFakeDataInput, false, false);
ParserUtils::ParseArgument(args, "-s", fStartThrows, false, false);
ParserUtils::ParseArgument(args, "-t", fNThrows, false, false);
ParserUtils::ParseArgument(args, "-p", fThrowString, false, false);
ParserUtils::ParseArgument(args, "-i", xmlcmds);
ParserUtils::ParseArgument(args, "-q", configargs);
ParserUtils::ParseCounter(args, "e", errorcount);
ParserUtils::ParseCounter(args, "v", verbocount);
ParserUtils::CheckBadArguments(args);
// Add extra defaults if none given
if (fCardFile.empty() and xmlcmds.empty()) {
ERR(FTL) << "No input supplied!" << std::endl;
throw;
}
if (fOutputFile.empty() and !fCardFile.empty()) {
fOutputFile = fCardFile + ".root";
ERR(WRN) << "No output supplied so saving it to: " << fOutputFile << std::endl;
} else if (fOutputFile.empty()) {
ERR(FTL) << "No output file or cardfile supplied!" << std::endl;
throw;
}
// Configuration Setup =============================
// Check no comp key is available
if (Config::Get().GetNodes("nuiscomp").empty()) {
fCompKey = Config::Get().CreateNode("nuiscomp");
} else {
fCompKey = Config::Get().GetNodes("nuiscomp")[0];
}
if (!fCardFile.empty()) fCompKey.AddS("cardfile", fCardFile);
if (!fOutputFile.empty()) fCompKey.AddS("outputfile", fOutputFile);
if (!fStrategy.empty()) fCompKey.AddS("strategy", fStrategy);
// Load XML Cardfile
configuration.LoadConfig( fCompKey.GetS("cardfile"), "");
// Add CMD XML Structs
for (size_t i = 0; i < xmlcmds.size(); i++) {
configuration.AddXMLLine(xmlcmds[i]);
}
// Add Config Args
for (size_t i = 0; i < configargs.size(); i++) {
configuration.OverrideConfig(configargs[i]);
}
if (maxevents.compare("-1")){
configuration.OverrideConfig("MAXEVENTS=" + maxevents);
}
// Finish configuration XML
configuration.FinaliseConfig(fCompKey.GetS("outputfile") + ".xml");
// Add Error Verbo Lines
verbocount += Config::Get().GetParI("VERBOSITY");
errorcount += Config::Get().GetParI("ERROR");
std::cout << "[ NUISANCE ]: Setting VERBOSITY=" << verbocount << std::endl;
std::cout << "[ NUISANCE ]: Setting ERROR=" << errorcount << std::endl;
SETVERBOSITY(verbocount);
// Proper Setup
if (fStrategy.find("ErrorBands") != std::string::npos ||
fStrategy.find("MergeErrors") != std::string::npos){
fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "RECREATE");
}
// fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "RECREATE");
SetupSystematicsFromXML();
SetupCovariance();
SetupRWEngine();
SetupFCN();
GetCovarFromFCN();
// Run();
return;
};
void SystematicRoutines::SetupSystematicsFromXML(){
LOG(FIT) << "Setting up nuismin" << std::endl;
// Setup Parameters ------------------------------------------
std::vector<nuiskey> parkeys = Config::QueryKeys("parameter");
if (!parkeys.empty()) {
LOG(FIT) << "Number of parameters : " << parkeys.size() << std::endl;
}
for (size_t i = 0; i < parkeys.size(); i++) {
nuiskey key = parkeys.at(i);
// Check for type,name,nom
if (!key.Has("type")) {
ERR(FTL) << "No type given for parameter " << i << std::endl;
throw;
} else if (!key.Has("name")) {
ERR(FTL) << "No name given for parameter " << i << std::endl;
throw;
} else if (!key.Has("nominal")) {
ERR(FTL) << "No nominal given for parameter " << i << std::endl;
throw;
}
// Get Inputs
std::string partype = key.GetS("type");
std::string parname = key.GetS("name");
double parnom = key.GetD("nominal");
double parlow = parnom - 1;
double parhigh = parnom + 1;
double parstep = 1;
// Override if state not given
if (!key.Has("state")){
key.SetS("state","FIX");
}
std::string parstate = key.GetS("state");
// Extra limits
if (key.Has("low")) {
parlow = key.GetD("low");
parhigh = key.GetD("high");
parstep = key.GetD("step");
LOG(FIT) << "Read " << partype << " : "
<< parname << " = "
<< parnom << " : "
<< parlow << " < p < " << parhigh
<< " : " << parstate << std::endl;
} else {
LOG(FIT) << "Read " << partype << " : "
<< parname << " = "
<< parnom << " : "
<< parstate << std::endl;
}
// Run Parameter Conversion if needed
if (parstate.find("ABS") != std::string::npos) {
parnom = FitBase::RWAbsToSigma( partype, parname, parnom );
parlow = FitBase::RWAbsToSigma( partype, parname, parlow );
parhigh = FitBase::RWAbsToSigma( partype, parname, parhigh );
parstep = FitBase::RWAbsToSigma( partype, parname, parstep );
} else if (parstate.find("FRAC") != std::string::npos) {
parnom = FitBase::RWFracToSigma( partype, parname, parnom );
parlow = FitBase::RWFracToSigma( partype, parname, parlow );
parhigh = FitBase::RWFracToSigma( partype, parname, parhigh );
parstep = FitBase::RWFracToSigma( partype, parname, parstep );
}
// Push into vectors
fParams.push_back(parname);
fTypeVals[parname] = FitBase::ConvDialType(partype);;
fStartVals[parname] = parnom;
fCurVals[parname] = parnom;
fErrorVals[parname] = 0.0;
fStateVals[parname] = parstate;
bool fixstate = parstate.find("FIX") != std::string::npos;
fFixVals[parname] = fixstate;
fStartFixVals[parname] = fFixVals[parname];
fMinVals[parname] = parlow;
fMaxVals[parname] = parhigh;
fStepVals[parname] = parstep;
}
// Setup Samples ----------------------------------------------
std::vector<nuiskey> samplekeys = Config::QueryKeys("sample");
if (!samplekeys.empty()) {
LOG(FIT) << "Number of samples : " << samplekeys.size() << std::endl;
}
for (size_t i = 0; i < samplekeys.size(); i++) {
nuiskey key = samplekeys.at(i);
// Get Sample Options
std::string samplename = key.GetS("name");
std::string samplefile = key.GetS("input");
std::string sampletype =
key.Has("type") ? key.GetS("type") : "DEFAULT";
double samplenorm =
key.Has("norm") ? key.GetD("norm") : 1.0;
// Print out
LOG(FIT) << "Read sample info " << i << " : "
<< samplename << std::endl
<< "\t\t input -> " << samplefile << std::endl
<< "\t\t state -> " << sampletype << std::endl
<< "\t\t norm -> " << samplenorm << std::endl;
// If FREE add to parameters otherwise continue
if (sampletype.find("FREE") == std::string::npos) {
continue;
}
// Form norm dial from samplename + sampletype + "_norm";
std::string normname = samplename + "_norm";
// Check normname not already present
if (fTypeVals.find(normname) != fTypeVals.end()) {
continue;
}
// Add new norm dial to list if its passed above checks
fParams.push_back(normname);
fTypeVals[normname] = kNORM;
fStateVals[normname] = sampletype;
fCurVals[normname] = samplenorm;
fErrorVals[normname] = 0.0;
fMinVals[normname] = 0.1;
fMaxVals[normname] = 10.0;
fStepVals[normname] = 0.5;
bool state = sampletype.find("FREE") == std::string::npos;
fFixVals[normname] = state;
fStartFixVals[normname] = state;
}
// Setup Fake Parameters -----------------------------
std::vector<nuiskey> fakekeys = Config::QueryKeys("fakeparameter");
if (!fakekeys.empty()) {
LOG(FIT) << "Number of fake parameters : " << fakekeys.size() << std::endl;
}
for (size_t i = 0; i < fakekeys.size(); i++) {
nuiskey key = fakekeys.at(i);
// Check for type,name,nom
if (!key.Has("name")) {
ERR(FTL) << "No name given for fakeparameter " << i << std::endl;
throw;
} else if (!key.Has("nom")) {
ERR(FTL) << "No nominal given for fakeparameter " << i << std::endl;
throw;
}
// Get Inputs
std::string parname = key.GetS("name");
double parnom = key.GetD("nom");
// Push into vectors
fFakeVals[parname] = parnom;
}
}
void SystematicRoutines::ReadCard(std::string cardfile){
// Read cardlines into vector
std::vector<std::string> cardlines = GeneralUtils::ParseFileToStr(cardfile,"\n");
FitPar::Config().cardLines = cardlines;
// Read Samples first (norm params can be overridden)
int linecount = 0;
for (std::vector<std::string>::iterator iter = cardlines.begin();
iter != cardlines.end(); iter++){
std::string line = (*iter);
linecount++;
// Skip Empties
if (line.empty()) continue;
if (line.c_str()[0] == '#') continue;
// Read Valid Samples
int samstatus = ReadSamples(line);
// Show line if bad to help user
if (samstatus == kErrorStatus) {
ERR(FTL) << "Bad Input in cardfile " << fCardFile
<< " at line " << linecount << "!" << std::endl;
LOG(FIT) << line << std::endl;
throw;
}
}
// Read Parameters second
linecount = 0;
for (std::vector<std::string>::iterator iter = cardlines.begin();
iter != cardlines.end(); iter++){
std::string line = (*iter);
linecount++;
// Skip Empties
if (line.empty()) continue;
if (line.c_str()[0] == '#') continue;
// Try Parameter Reads
int parstatus = ReadParameters(line);
int fakstatus = ReadFakeDataPars(line);
// Show line if bad to help user
if (parstatus == kErrorStatus ||
fakstatus == kErrorStatus ){
ERR(FTL) << "Bad Parameter Input in cardfile " << fCardFile
<< " at line " << linecount << "!" << std::endl;
LOG(FIT) << line << std::endl;
throw;
}
}
return;
};
int SystematicRoutines::ReadParameters(std::string parstring){
std::string inputspec = "RW Dial Inputs Syntax \n"
"free input w/ limits: TYPE NAME START MIN MAX STEP [STATE] \n"
"fix input: TYPE NAME VALUE [STATE] \n"
"free input w/o limits: TYPE NAME START FREE,[STATE] \n"
"Allowed Types: \n"
"neut_parameter,niwg_parameter,t2k_parameter,"
"nuwro_parameter,gibuu_parameter";
// Check sample input
if (parstring.find("parameter") == std::string::npos) return kGoodStatus;
// Parse inputs
std::vector<std::string> strvct = GeneralUtils::ParseToStr(parstring, " ");
// Skip if comment or parameter somewhere later in line
if (strvct[0].c_str()[0] == '#' ||
strvct[0].find("parameter") == std::string::npos){
return kGoodStatus;
}
// Check length
if (strvct.size() < 3){
ERR(FTL) << "Input rw dials need to provide at least 3 inputs." << std::endl;
std::cout << inputspec << std::endl;
return kErrorStatus;
}
// Setup default inputs
std::string partype = strvct[0];
std::string parname = strvct[1];
double parval = GeneralUtils::StrToDbl(strvct[2]);
double minval = parval - 1.0;
double maxval = parval + 1.0;
double stepval = 1.0;
std::string state = "FIX"; //[DEFAULT]
// Check Type
if (FitBase::ConvDialType(partype) == kUNKNOWN){
ERR(FTL) << "Unknown parameter type! " << partype << std::endl;
std::cout << inputspec << std::endl;
return kErrorStatus;
}
// Check Parameter Name
if (FitBase::GetDialEnum(partype, parname) == -1){
ERR(FTL) << "Bad RW parameter name! " << partype << " " << parname << std::endl;
std::cout << inputspec << std::endl;
return kErrorStatus;
}
// Option Extra (No Limits)
if (strvct.size() == 4){
state = strvct[3];
}
// Check for weirder inputs
if (strvct.size() > 4 && strvct.size() < 6){
ERR(FTL) << "Provided incomplete limits for " << parname << std::endl;
std::cout << inputspec << std::endl;
return kErrorStatus;
}
// Option Extra (With limits and steps)
if (strvct.size() >= 6){
minval = GeneralUtils::StrToDbl(strvct[3]);
maxval = GeneralUtils::StrToDbl(strvct[4]);
stepval = GeneralUtils::StrToDbl(strvct[5]);
}
// Option Extra (dial state after limits)
if (strvct.size() == 7){
state = strvct[6];
}
// Run Parameter Conversion if needed
if (state.find("ABS") != std::string::npos){
parval = FitBase::RWAbsToSigma( partype, parname, parval );
minval = FitBase::RWAbsToSigma( partype, parname, minval );
maxval = FitBase::RWAbsToSigma( partype, parname, maxval );
stepval = FitBase::RWAbsToSigma( partype, parname, stepval );
} else if (state.find("FRAC") != std::string::npos){
parval = FitBase::RWFracToSigma( partype, parname, parval );
minval = FitBase::RWFracToSigma( partype, parname, minval );
maxval = FitBase::RWFracToSigma( partype, parname, maxval );
stepval = FitBase::RWFracToSigma( partype, parname, stepval );
}
// Check no repeat params
if (std::find(fParams.begin(), fParams.end(), parname) != fParams.end()){
ERR(FTL) << "Duplicate parameter names given for " << parname << std::endl;
throw;
}
// Setup Containers
fParams.push_back(parname);
fTypeVals[parname] = FitBase::ConvDialType(partype);
fStartVals[parname] = parval;
fCurVals[parname] = fStartVals[parname];
fErrorVals[parname] = 0.0;
fStateVals[parname] = state;
bool fixstate = state.find("FIX") != std::string::npos;
fFixVals[parname] = fixstate;
fStartFixVals[parname] = fFixVals[parname];
fMinVals[parname] = minval;
fMaxVals[parname] = maxval;
fStepVals[parname] = stepval;
// Print the parameter
LOG(MIN) << "Read Parameter " << parname << " " << parval << " "
<< minval << " " << maxval << " "
<< stepval << " " << state << std::endl;
// Tell reader its all good
return kGoodStatus;
}
//*******************************************
int SystematicRoutines::ReadFakeDataPars(std::string parstring){
//******************************************
std::string inputspec = "Fake Data Dial Inputs Syntax \n"
"fake value: fake_parameter NAME VALUE \n"
"Name should match dialnames given in actual dial specification.";
// Check sample input
if (parstring.find("fake_parameter") == std::string::npos)
return kGoodStatus;
// Parse inputs
std::vector<std::string> strvct = GeneralUtils::ParseToStr(parstring, " ");
// Skip if comment or parameter somewhere later in line
if (strvct[0].c_str()[0] == '#' ||
strvct[0] == "fake_parameter"){
return kGoodStatus;
}
// Check length
if (strvct.size() < 3){
ERR(FTL) << "Fake dials need to provide at least 3 inputs." << std::endl;
std::cout << inputspec << std::endl;
return kErrorStatus;
}
// Read Inputs
std::string parname = strvct[1];
double parval = GeneralUtils::StrToDbl(strvct[2]);
// Setup Container
fFakeVals[parname] = parval;
// Print the fake parameter
LOG(MIN) << "Read Fake Parameter " << parname << " " << parval << std::endl;
// Tell reader its all good
return kGoodStatus;
}
//******************************************
int SystematicRoutines::ReadSamples(std::string samstring){
//******************************************
const static std::string inputspec =
"\tsample <sample_name> <input_type>:inputfile.root [OPTS] "
"[norm]\nsample_name: Name "
"of sample to include. e.g. MiniBooNE_CCQE_XSec_1DQ2_nu\ninput_type: The "
"input event format. e.g. NEUT, GENIE, EVSPLN, ...\nOPTS: Additional, "
"optional sample options.\nnorm: Additional, optional sample "
"normalisation factor.";
// Check sample input
if (samstring.find("sample") == std::string::npos)
return kGoodStatus;
// Parse inputs
std::vector<std::string> strvct = GeneralUtils::ParseToStr(samstring, " ");
// Skip if comment or parameter somewhere later in line
if (strvct[0].c_str()[0] == '#' ||
strvct[0] != "sample"){
return kGoodStatus;
}
// Check length
if (strvct.size() < 3){
ERR(FTL) << "Sample need to provide at least 3 inputs." << std::endl;
return kErrorStatus;
}
// Setup default inputs
std::string samname = strvct[1];
std::string samfile = strvct[2];
if (samfile == "FIX") {
ERR(FTL) << "Input filename was \"FIX\", this line is probably malformed "
"in the input card file. Line:\'"
<< samstring << "\'" << std::endl;
ERR(FTL) << "Expect sample lines to look like:\n\t" << inputspec
<< std::endl;
throw;
}
std::string samtype = "DEFAULT";
double samnorm = 1.0;
// Optional Type
if (strvct.size() > 3) {
samtype = strvct[3];
// samname += "_"+samtype;
// Also get rid of the / and replace it with underscore because it might not be supported character
// while (samname.find("/") != std::string::npos) {
// samname.replace(samname.find("/"), 1, std::string("_"));
// }
}
// Optional Norm
if (strvct.size() > 4) samnorm = GeneralUtils::StrToDbl(strvct[4]);
// Add Sample Names as Norm Dials
std::string normname = samname + "_norm";
// Check no repeat params
if (std::find(fParams.begin(), fParams.end(), normname) != fParams.end()){
ERR(FTL) << "Duplicate samples given for " << samname << std::endl;
throw;
}
fParams.push_back(normname);
fTypeVals[normname] = kNORM;
fStartVals[normname] = samnorm;
fCurVals[normname] = fStartVals[normname];
fErrorVals[normname] = 0.0;
fMinVals[normname] = 0.1;
fMaxVals[normname] = 10.0;
fStepVals[normname] = 0.5;
bool state = samtype.find("FREE") == std::string::npos;
fFixVals[normname] = state;
fStartFixVals[normname] = state;
// Print read in
LOG(MIN) << "Read sample " << samname << " "
<< samfile << " " << samtype << " "
<< samnorm << std::endl;
// Tell reader its all good
return kGoodStatus;
}
/*
Setup Functions
*/
//*************************************
void SystematicRoutines::SetupRWEngine(){
//*************************************
for (UInt_t i = 0; i < fParams.size(); i++){
std::string name = fParams[i];
FitBase::GetRW() -> IncludeDial(name, fTypeVals.at(name) );
}
UpdateRWEngine(fStartVals);
return;
}
//*************************************
void SystematicRoutines::SetupFCN(){
//*************************************
LOG(FIT)<<"Making the jointFCN"<<std::endl;
if (fSampleFCN) delete fSampleFCN;
fSampleFCN = new JointFCN(fOutputRootFile);
SetFakeData();
return;
}
//*************************************
void SystematicRoutines::SetFakeData(){
//*************************************
if (fFakeDataInput.empty()) return;
if (fFakeDataInput.compare("MC") == 0){
LOG(FIT)<<"Setting fake data from MC starting prediction." <<std::endl;
UpdateRWEngine(fFakeVals);
FitBase::GetRW()->Reconfigure();
fSampleFCN->ReconfigureAllEvents();
fSampleFCN->SetFakeData("MC");
UpdateRWEngine(fCurVals);
LOG(FIT)<<"Set all data to fake MC predictions."<<std::endl;
} else {
fSampleFCN->SetFakeData(fFakeDataInput);
}
return;
}
//*****************************************
void SystematicRoutines::GetCovarFromFCN(){
//*****************************************
LOG(FIT) << "Loading ParamPull objects from FCN to build covar" << std::endl;
// Make helperstring
std::ostringstream helperstr;
// Keep track of what is being thrown
std::map<std::string, std::string> dialthrowhandle;
// Get Covariance Objects from FCN
std::list<ParamPull*> inputpulls = fSampleFCN->GetPullList();
for (PullListConstIter iter = inputpulls.begin();
iter != inputpulls.end(); iter++){
ParamPull* pull = (*iter);
if (pull->GetType().find("THROW")){
fInputThrows.push_back(pull);
fInputCovar.push_back(pull->GetFullCovarMatrix());
fInputDials.push_back(pull->GetDataHist());
LOG(FIT) << "Read ParamPull: " << pull->GetName() << " " << pull->GetType() << std::endl;
}
TH1D dialhist = pull->GetDataHist();
TH1D minhist = pull->GetMinHist();
TH1D maxhist = pull->GetMaxHist();
TH1I typehist = pull->GetDialTypes();
for (int i = 0; i < dialhist.GetNbinsX(); i++){
std::string name = std::string(dialhist.GetXaxis()->GetBinLabel(i+1));
dialthrowhandle[name] = pull->GetName();
if (fCurVals.find(name) == fCurVals.end()){
// Add to Containers
fParams.push_back(name);
fCurVals[name] = dialhist.GetBinContent(i+1);
fStartVals[name] = dialhist.GetBinContent(i+1);
fMinVals[name] = minhist.GetBinContent(i+1);
fMaxVals[name] = maxhist.GetBinContent(i+1);
fStepVals[name] = 1.0;
fFixVals[name] = false;
fStartFixVals[name] = false;
fTypeVals[name] = typehist.GetBinContent(i+1);
fStateVals[name] = "FREE" + pull->GetType();
// Maker Helper
helperstr << std::string(16, ' ' ) << FitBase::ConvDialType(fTypeVals[name]) << " "
<< name << " " << fMinVals[name] << " "
<< fMaxVals[name] << " " << fStepVals[name] << " " << fStateVals[name]
<< std::endl;
}
}
}
// Check if no throws given
if (fInputThrows.empty()){
ERR(WRN) << "No covariances given to nuissyst" << std::endl;
ERR(WRN) << "Pushing back an uncorrelated gaussian throw error for each free parameter using step size" << std::endl;
for (UInt_t i = 0; i < fParams.size(); i++){
std::string syst = fParams[i];
if (fFixVals[syst]) continue;
// Make Terms
std::string name = syst + "_pull";
std::ostringstream pullterm;
pullterm << "DIAL:" << syst << ";"
<< fStartVals[syst] << ";"
<< fStepVals[syst];
std::string type = "GAUSTHROW/NEUT";
// Push Back Pulls
ParamPull* pull = new ParamPull( name, pullterm.str(), type );
fInputThrows.push_back(pull);
fInputCovar.push_back(pull->GetFullCovarMatrix());
fInputDials.push_back(pull->GetDataHist());
// Print Whats added
ERR(WRN) << "Added ParamPull : " << name << " " << pullterm.str() << " " << type << std::endl;
// Add helper string for future fits
helperstr << std::string(16, ' ' ) << "covar " << name << " " << pullterm.str() << " " << type << std::endl;
// Keep Track of Throws
dialthrowhandle[syst] = pull->GetName();
}
}
// Print Helper String
if (!helperstr.str().empty()){
LOG(FIT) << "To remove these statements in future studies, add the lines below to your card:" << std::endl;
// Can't use the logger properly because this can be multi-line. Use cout and added spaces to look better!
std::cout << helperstr.str();
sleep(2);
}
// Print Throw State
for (UInt_t i = 0; i < fParams.size(); i++){
std::string syst = fParams[i];
if (dialthrowhandle.find(syst) != dialthrowhandle.end()){
LOG(FIT) << "Dial " << i << ". " << setw(40) << syst << " = THROWING with " << dialthrowhandle[syst] << std::endl;
} else {
LOG(FIT) << "Dial " << i << ". " << setw(40) << syst << " = FIXED" << std::endl;
}
}
// Pause anyway
sleep(1);
return;
}
/*
Fitting Functions
*/
//*************************************
void SystematicRoutines::UpdateRWEngine(std::map<std::string,double>& updateVals){
//*************************************
for (UInt_t i = 0; i < fParams.size(); i++){
std::string name = fParams[i];
if (updateVals.find(name) == updateVals.end()) continue;
FitBase::GetRW()->SetDialValue(name,updateVals.at(name));
}
FitBase::GetRW()->Reconfigure();
return;
}
//*************************************
void SystematicRoutines::PrintState(){
//*************************************
LOG(FIT)<<"------------"<<std::endl;
// Count max size
int maxcount = 0;
for (UInt_t i = 0; i < fParams.size(); i++){
maxcount = max(int(fParams[i].size()), maxcount);
}
// Header
LOG(FIT) << " # " << left << setw(maxcount) << "Parameter "
<< " = "
<< setw(10) << "Value" << " +- "
<< setw(10) << "Error" << " "
<< setw(8) << "(Units)" << " "
<< setw(10) << "Conv. Val" << " +- "
<< setw(10) << "Conv. Err" << " "
<< setw(8) << "(Units)" << std::endl;
// Parameters
for (UInt_t i = 0; i < fParams.size(); i++){
std::string syst = fParams.at(i);
std::string typestr = FitBase::ConvDialType(fTypeVals[syst]);
std::string curunits = "(sig.)";
double curval = fCurVals[syst];
double curerr = fErrorVals[syst];
if (fStateVals[syst].find("ABS") != std::string::npos){
curval = FitBase::RWSigmaToAbs(typestr, syst, curval);
curerr = (FitBase::RWSigmaToAbs(typestr, syst, curerr) -
FitBase::RWSigmaToAbs(typestr, syst, 0.0));
curunits = "(Abs.)";
} else if (fStateVals[syst].find("FRAC") != std::string::npos){
curval = FitBase::RWSigmaToFrac(typestr, syst, curval);
curerr = (FitBase::RWSigmaToFrac(typestr, syst, curerr) -
FitBase::RWSigmaToFrac(typestr, syst, 0.0));
curunits = "(Frac)";
}
std::string convunits = "(" + FitBase::GetRWUnits(typestr, syst) + ")";
double convval = FitBase::RWSigmaToAbs(typestr, syst, curval);
double converr = (FitBase::RWSigmaToAbs(typestr, syst, curerr) -
FitBase::RWSigmaToAbs(typestr, syst, 0.0));
std::ostringstream curparstring;
curparstring << " " << setw(3) << left
<< i << ". "
<< setw(maxcount) << syst << " = "
<< setw(10) << curval << " +- "
<< setw(10) << curerr << " "
<< setw(8) << curunits << " "
<< setw(10) << convval << " +- "
<< setw(10) << converr << " "
<< setw(8) << convunits;
LOG(FIT) << curparstring.str() << std::endl;
}
LOG(FIT)<<"------------"<<std::endl;
double like = fSampleFCN->GetLikelihood();
LOG(FIT) << std::left << std::setw(46) << "Likelihood for JointFCN: " << like << std::endl;
LOG(FIT)<<"------------"<<std::endl;
}
/*
Write Functions
*/
//*************************************
void SystematicRoutines::SaveResults(){
//*************************************
if (!fOutputRootFile)
fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "RECREATE");
fOutputRootFile->cd();
SaveCurrentState();
}
//*************************************
void SystematicRoutines::SaveCurrentState(std::string subdir){
//*************************************
LOG(FIT)<<"Saving current full FCN predictions" <<std::endl;
// Setup DIRS
TDirectory* curdir = gDirectory;
if (!subdir.empty()){
TDirectory* newdir =(TDirectory*) gDirectory->mkdir(subdir.c_str());
newdir->cd();
}
FitBase::GetRW()->Reconfigure();
fSampleFCN->ReconfigureAllEvents();
fSampleFCN->Write();
// Change back to current DIR
curdir->cd();
return;
}
//*************************************
void SystematicRoutines::SaveNominal(){
//*************************************
if (!fOutputRootFile)
fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "RECREATE");
fOutputRootFile->cd();
LOG(FIT)<<"Saving Nominal Predictions (be cautious with this)" <<std::endl;
FitBase::GetRW()->Reconfigure();
SaveCurrentState("nominal");
};
//*************************************
void SystematicRoutines::SavePrefit(){
//*************************************
if (!fOutputRootFile)
fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "RECREATE");
fOutputRootFile->cd();
LOG(FIT)<<"Saving Prefit Predictions"<<std::endl;
UpdateRWEngine(fStartVals);
SaveCurrentState("prefit");
UpdateRWEngine(fCurVals);
};
/*
MISC Functions
*/
//*************************************
int SystematicRoutines::GetStatus(){
//*************************************
return 0;
}
//*************************************
void SystematicRoutines::SetupCovariance(){
//*************************************
// Remove covares if they exist
if (fCovar) delete fCovar;
if (fCovarFree) delete fCovarFree;
if (fCorrel) delete fCorrel;
if (fCorrelFree) delete fCorrelFree;
if (fDecomp) delete fDecomp;
if (fDecompFree) delete fDecompFree;
int NFREE = 0;
int NDIM = 0;
// Get NFREE from min or from vals (for cases when doing throws)
NDIM = fParams.size();
for (UInt_t i = 0; i < fParams.size(); i++){
if (!fFixVals[fParams[i]]) NFREE++;
}
if (NDIM == 0) return;
fCovar = new TH2D("covariance","covariance",NDIM,0,NDIM,NDIM,0,NDIM);
if (NFREE > 0){
fCovarFree = new TH2D("covariance_free",
"covariance_free",
NFREE,0,NFREE,
NFREE,0,NFREE);
}
// Set Bin Labels
int countall = 0;
int countfree = 0;
for (UInt_t i = 0; i < fParams.size(); i++){
fCovar->GetXaxis()->SetBinLabel(countall+1,fParams[i].c_str());
fCovar->GetYaxis()->SetBinLabel(countall+1,fParams[i].c_str());
countall++;
if (!fFixVals[fParams[i]] and NFREE > 0){
fCovarFree->GetXaxis()->SetBinLabel(countfree+1,fParams[i].c_str());
fCovarFree->GetYaxis()->SetBinLabel(countfree+1,fParams[i].c_str());
countfree++;
}
}
fCorrel = PlotUtils::GetCorrelationPlot(fCovar,"correlation");
fDecomp = PlotUtils::GetDecompPlot(fCovar,"decomposition");
if (NFREE > 0)fCorrelFree = PlotUtils::GetCorrelationPlot(fCovarFree, "correlation_free");
if (NFREE > 0)fDecompFree = PlotUtils::GetDecompPlot(fCovarFree,"decomposition_free");
return;
};
//*************************************
void SystematicRoutines::ThrowCovariance(bool uniformly){
//*************************************
// Set fThrownVals to all values in currentVals
for (UInt_t i = 0; i < fParams.size(); i++){
std::string name = fParams.at(i);
fThrownVals[name] = fCurVals[name];
}
for (PullListConstIter iter = fInputThrows.begin();
iter != fInputThrows.end(); iter++){
ParamPull* pull = *iter;
pull->ThrowCovariance();
TH1D dialhist = pull->GetDataHist();
for (int i = 0; i < dialhist.GetNbinsX(); i++){
std::string name = std::string(dialhist.GetXaxis()->GetBinLabel(i+1));
if (fCurVals.find(name) != fCurVals.end()){
fThrownVals[name] = dialhist.GetBinContent(i+1);
}
}
// Reset throw incase pulls are calculated.
pull->ResetToy();
}
return;
};
//*************************************
void SystematicRoutines::PlotLimits(){
//*************************************
std::cout << "Plotting Limits" << std::endl;
if (!fOutputRootFile)
fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "RECREATE");
TDirectory* limfolder = (TDirectory*) fOutputRootFile->mkdir("Limits");
limfolder->cd();
// Set all parameters at their starting values
for (UInt_t i = 0; i < fParams.size(); i++){
fCurVals[fParams[i]] = fStartVals[fParams[i]];
}
TDirectory* nomfolder = (TDirectory*) limfolder->mkdir("nominal");
nomfolder->cd();
UpdateRWEngine(fCurVals);
fSampleFCN->ReconfigureAllEvents();
fSampleFCN->Write();
limfolder->cd();
std::vector<std::string> allfolders;
// Loop through each parameter
for (UInt_t i = 0; i < fParams.size(); i++){
std::string syst = fParams[i];
std::cout << "Starting Param " << syst << std::endl;
if (fFixVals[syst]) continue;
// Loop Downwards
while (fCurVals[syst] > fMinVals[syst]){
fCurVals[syst] = fCurVals[syst] - fStepVals[syst];
// Check Limit
if (fCurVals[syst] < fMinVals[syst])
fCurVals[syst] = fMinVals[syst];
// Check folder exists
std::string curvalstring = std::string( Form( (syst + "_%f").c_str(), fCurVals[syst] ) );
if (std::find(allfolders.begin(), allfolders.end(), curvalstring) != allfolders.end())
break;
// Make new folder for variation
TDirectory* minfolder = (TDirectory*) limfolder->mkdir(Form( (syst + "_%f").c_str(), fCurVals[syst] ) );
minfolder->cd();
allfolders.push_back(curvalstring);
// Update Iterations
double *vals = FitUtils::GetArrayFromMap( fParams, fCurVals );
fSampleFCN->DoEval( vals );
delete vals;
// Save to folder
fSampleFCN->Write();
}
// Reset before next loop
fCurVals[syst] = fStartVals[syst];
// Loop Upwards now
while (fCurVals[syst] < fMaxVals[syst]){
fCurVals[syst] = fCurVals[syst] + fStepVals[syst];
// Check Limit
if (fCurVals[syst] > fMaxVals[syst])
fCurVals[syst] = fMaxVals[syst];
// Check folder exists
std::string curvalstring = std::string( Form( (syst + "_%f").c_str(), fCurVals[syst] ) );
if (std::find(allfolders.begin(), allfolders.end(), curvalstring) != allfolders.end())
break;
// Make new folder
TDirectory* maxfolder = (TDirectory*) limfolder->mkdir(Form( (syst + "_%f").c_str(), fCurVals[syst] ) );
maxfolder->cd();
allfolders.push_back(curvalstring);
// Update Iterations
double *vals = FitUtils::GetArrayFromMap( fParams, fCurVals );
fSampleFCN->DoEval( vals );
delete vals;
// Save to file
fSampleFCN->Write();
}
// Reset before leaving
fCurVals[syst] = fStartVals[syst];
UpdateRWEngine(fCurVals);
}
return;
}
//*************************************
void SystematicRoutines::Run(){
//*************************************
std::cout << "Running routines "<< std::endl;
fRoutines = GeneralUtils::ParseToStr(fStrategy,",");
for (UInt_t i = 0; i < fRoutines.size(); i++){
std::string routine = fRoutines.at(i);
int fitstate = kFitUnfinished;
LOG(FIT)<<"Running Routine: "<<routine<<std::endl;
if (routine.compare("PlotLimits") == 0) PlotLimits();
else if (routine.compare("ErrorBands") == 0) GenerateErrorBands();
else if (routine.compare("ThrowErrors") == 0) GenerateThrows();
else if (routine.compare("MergeErrors") == 0) MergeThrows();
else {
std::cout << "Unknown ROUTINE : " << routine << std::endl;
}
// If ending early break here
if (fitstate == kFitFinished || fitstate == kNoChange){
LOG(FIT) << "Ending fit routines loop." << std::endl;
break;
}
}
return;
}
void SystematicRoutines::GenerateErrorBands(){
GenerateThrows();
MergeThrows();
}
//*************************************
void SystematicRoutines::GenerateThrows(){
//*************************************
TFile* tempfile = new TFile((fOutputFile + ".throws.root").c_str(),"RECREATE");
tempfile->cd();
int nthrows = fNThrows;
int startthrows = fStartThrows;
int endthrows = startthrows + nthrows;
if (nthrows < 0) nthrows = endthrows;
if (startthrows < 0) startthrows = 0;
if (endthrows < 0) endthrows = startthrows + nthrows;
- int seed = (gRandom->Uniform(0.0,1.0)*100000 + 100000000*(startthrows + endthrows) + time(NULL))/35;
+ int seed = (gRandom->Uniform(0.0,1.0)*100000 + 100000000*(startthrows + endthrows) + time(NULL) + int(getpid()) );
gRandom->SetSeed(seed);
LOG(FIT) << "Using Seed : " << seed << std::endl;
LOG(FIT) << "nthrows = " << nthrows << std::endl;
LOG(FIT) << "startthrows = " << startthrows << std::endl;
LOG(FIT) << "endthrows = " << endthrows << std::endl;
UpdateRWEngine(fCurVals);
fSampleFCN->ReconfigureAllEvents();
if (startthrows == 0){
LOG(FIT) << "Making nominal " << std::endl;
TDirectory* nominal = (TDirectory*) tempfile->mkdir("nominal");
nominal->cd();
fSampleFCN->Write();
}
LOG(SAM) << "nthrows = " << nthrows << std::endl;
LOG(SAM) << "startthrows = " << startthrows << std::endl;
LOG(SAM) << "endthrows = " << endthrows << std::endl;
TTree* parameterTree = new TTree("throws","throws");
double chi2;
for (UInt_t i = 0; i < fParams.size(); i++)
parameterTree->Branch(fParams[i].c_str(), &fThrownVals[fParams[i]], (fParams[i] + "/D").c_str());
parameterTree->Branch("chi2",&chi2,"chi2/D");
fSampleFCN->CreateIterationTree("error_iterations", FitBase::GetRW());
// Would anybody actually want to do uniform throws of any parameter??
bool uniformly = FitPar::Config().GetParB("error_uniform");
// Run Throws and save
for (Int_t i = 0; i < endthrows+1; i++){
LOG(FIT) << "Loop " << i << std::endl;
ThrowCovariance(uniformly);
if (i < startthrows) continue;
if (i == 0) continue;
LOG(FIT) << "Throw " << i << " ================================" << std::endl;
// Generate Random Parameter Throw
// ThrowCovariance(uniformly);
TDirectory* throwfolder = (TDirectory*)tempfile->mkdir(Form("throw_%i",i));
throwfolder->cd();
// Run Eval
double *vals = FitUtils::GetArrayFromMap( fParams, fThrownVals );
chi2 = fSampleFCN->DoEval( vals );
delete vals;
// Save the FCN
fSampleFCN->Write();
parameterTree->Fill();
}
tempfile->cd();
fSampleFCN->WriteIterationTree();
tempfile->Close();
}
void SystematicRoutines::MergeThrows(){
fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "RECREATE");
fOutputRootFile->cd();
// Make a container folder
TDirectory* errorDIR = (TDirectory*) fOutputRootFile->mkdir("error_bands");
errorDIR->cd();
TDirectory* outnominal = (TDirectory*) fOutputRootFile->mkdir("nominal_throw");
outnominal->cd();
// Split Input Files
if (!fThrowString.empty()) fThrowList = GeneralUtils::ParseToStr(fThrowString,",");
// Add default if no throwlist given
if (fThrowList.size() < 1) fThrowList.push_back( fOutputFile + ".throws.root" );
/// Save location of file containing nominal
std::string nominalfile;
bool nominalfound;
// Loop over files and check they exist.
for (uint i = 0; i < fThrowList.size(); i++){
std::string file = fThrowList[i];
bool found = false;
// normal
std::string newfile = file;
TFile* throwfile = new TFile(file.c_str(),"READ");
if (throwfile and !throwfile->IsZombie()){
found = true;
}
// normal.throws.root
if (!found){
newfile = file + ".throws.root";
throwfile = new TFile((file + ".throws.root").c_str(),"READ");
if (throwfile and !throwfile->IsZombie()) {
found = true;
}
}
// If its found save to throwlist, else save empty.
// Also search for nominal
if (found){
fThrowList[i] = newfile;
LOG(FIT) << "Throws File :" << newfile << std::endl;
// Find input which contains nominal
if (throwfile->Get("nominal")){
nominalfound = true;
nominalfile = newfile;
}
throwfile->Close();
} else {
fThrowList[i] = "";
}
delete throwfile;
}
// Make sure we have a nominal file
if (!nominalfound or nominalfile.empty()){
ERR(FTL) << "No nominal found when mergining! Exiting!" << std::endl;
throw;
}
// Get the nominal throws file
TFile* tempfile = new TFile((nominalfile).c_str(),"READ");
tempfile->cd();
TDirectory* nominal = (TDirectory*)tempfile->Get("nominal");
// int nthrows = FitPar::Config().GetParI("error_throws");
bool uniformly = FitPar::Config().GetParB("error_uniform");
// Check percentage of bad files is okay.
int badfilecount = 0;
for (uint i = 0; i < fThrowList.size(); i++){
if (!fThrowList[i].empty()){
LOG(FIT) << "Loading Throws From File " << i << " : "
<< fThrowList[i] << std::endl;
} else {
badfilecount++;
}
}
// Check we have at least one good file
if ((uint)badfilecount == fThrowList.size()){
ERR(FTL) << "Found no good throw files for MergeThrows" << std::endl;
throw;
} else if (badfilecount > fThrowList.size()*0.25){
ERR(WRN) << "Over 25% of your throw files are dodgy. Please check this is okay!" << std::endl;
ERR(WRN) << "Will continue for the time being..." << std::endl;
sleep(5);
}
// Now go through the keys in the temporary file and look for TH1D, and TH2D plots
TIter next(nominal->GetListOfKeys());
TKey *key;
while ((key = (TKey*)next())) {
TClass *cl = gROOT->GetClass(key->GetClassName());
if (!cl->InheritsFrom("TH1D") and !cl->InheritsFrom("TH2D")) continue;
TH1* baseplot = (TH1D*)key->ReadObj();
std::string plotname = std::string(baseplot->GetName());
LOG(FIT) << "Creating error bands for " << plotname;
if (LOG_LEVEL(FIT)){
if (!uniformly) std::cout << " : Using COVARIANCE Throws! " << std::endl;
else std::cout << " : Using UNIFORM THROWS!!! " << std::endl;
}
int nbins = 0;
if (cl->InheritsFrom("TH1D")) nbins = ((TH1D*)baseplot)->GetNbinsX();
else nbins = ((TH1D*)baseplot)->GetNbinsX()* ((TH1D*)baseplot)->GetNbinsY();
// Setup TProfile with RMS option
TProfile* tprof = new TProfile((plotname + "_prof").c_str(),(plotname + "_prof").c_str(),nbins, 0, nbins, "S");
// Setup The TTREE
double* bincontents;
bincontents = new double[nbins];
double* binlowest;
binlowest = new double[nbins];
double* binhighest;
binhighest = new double[nbins];
errorDIR->cd();
TTree* bintree = new TTree((plotname + "_tree").c_str(), (plotname + "_tree").c_str());
for (Int_t i = 0; i < nbins; i++){
bincontents[i] = 0.0;
binhighest[i] = 0.0;
binlowest[i] = 0.0;
bintree->Branch(Form("content_%i",i),&bincontents[i],Form("content_%i/D",i));
}
// Make new throw plot
TH1* newplot;
// Run Throw Merging.
for (UInt_t i = 0; i < fThrowList.size(); i++){
TFile* throwfile = new TFile(fThrowList[i].c_str(), "READ");
// Loop over all throws in a folder
TIter nextthrow(throwfile->GetListOfKeys());
TKey *throwkey;
while ((throwkey = (TKey*)nextthrow())) {
// Skip non throw folders
if (std::string(throwkey->GetName()).find("throw_") == std::string::npos) continue;
// Get Throw DIR
TDirectory* throwdir = (TDirectory*)throwkey->ReadObj();
// Get Plot From Throw
newplot = (TH1*)throwdir->Get(plotname.c_str());
if (!newplot) continue;
// Loop Over Plot
for (Int_t j = 0; j < nbins; j++){
tprof->Fill(j+0.5, newplot->GetBinContent(j+1));
bincontents[j] = newplot->GetBinContent(j+1);
if (bincontents[j] < binlowest[j] or i == 0) binlowest[j] = bincontents[j];
if (bincontents[j] > binhighest[j] or i == 0) binhighest[j] = bincontents[j];
}
errorDIR->cd();
bintree->Fill();
}
}
errorDIR->cd();
if (uniformly){
LOG(FIT) << "Uniformly Calculating Plot Errors!" << std::endl;
}
TH1* statplot = (TH1*) baseplot->Clone();
for (Int_t j = 0; j < nbins; j++){
if (!uniformly){
// if ((baseplot->GetBinError(j+1)/baseplot->GetBinContent(j+1)) < 1.0) {
// baseplot->SetBinError(j+1,sqrt(pow(tprof->GetBinError(j+1),2) + pow(baseplot->GetBinError(j+1),2)));
// } else {
baseplot->SetBinContent(j+1,tprof->GetBinContent(j+1));
baseplot->SetBinError(j+1,tprof->GetBinError(j+1));
// }
} else {
baseplot->SetBinContent(j+1, 0.0);//(binlowest[j] + binhighest[j]) / 2.0);
baseplot->SetBinError(j+1, 0.0); //(binhighest[j] - binlowest[j])/2.0);
}
}
errorDIR->cd();
baseplot->Write();
tprof->Write();
bintree->Write();
outnominal->cd();
for (int i = 0; i < nbins; i++){
baseplot->SetBinError(i+1, sqrt(pow(statplot->GetBinError(i+1),2) + pow(baseplot->GetBinError(i+1),2)));
}
baseplot->Write();
delete statplot;
delete baseplot;
delete tprof;
delete bintree;
delete [] bincontents;
}
return;
};

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 7:46 PM (1 d, 7 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3804326
Default Alt Text
(139 KB)

Event Timeline