diff --git a/cmake/cacheVariables.cmake b/cmake/cacheVariables.cmake
index 70e0763..c106104 100644
--- a/cmake/cacheVariables.cmake
+++ b/cmake/cacheVariables.cmake
@@ -1,225 +1,241 @@
# 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 .
################################################################################
function(CheckAndSetDefaultEnv VARNAME DEFAULT CACHETYPE DOCSTRING ENVNAME)
#cmessage(DEBUG "Trying to assign variable ${VARNAME} into the cache.")
if(NOT DEFINED ${VARNAME})
if(DEFINED ENV{${ENVNAME}} AND NOT $ENV{${ENVNAME}} STREQUAL "")
set(${VARNAME} $ENV{${ENVNAME}} CACHE ${CACHETYPE} ${DOCSTRING})
cmessage(DEBUG " Read ${VARNAME} from ENVVAR ${ENVNAME} as $ENV{${ENVNAME}}.")
else()
set(${VARNAME} ${DEFAULT} CACHE ${CACHETYPE} ${DOCSTRING})
endif()
else()
set(${VARNAME} ${${VARNAME}} CACHE ${CACHETYPE} ${DOCSTRING})
unset(${VARNAME})
endif()
cmessage(CACHE "--Set cache variable: \"${VARNAME}\" to \"${${VARNAME}}\", in cache ${CACHETYPE}.")
endfunction()
function(CheckAndSetDefaultCache VARNAME DEFAULT CACHETYPE DOCSTRING)
# cmessage(DEBUG "Trying to assign variable ${VARNAME} into the cache.")
if(NOT DEFINED ${VARNAME})
set(${VARNAME} ${DEFAULT} CACHE ${CACHETYPE} ${DOCSTRING})
else()
set(${VARNAME} ${${VARNAME}} CACHE ${CACHETYPE} ${DOCSTRING})
unset(${VARNAME})
endif()
cmessage(CACHE "--Set cache variable: \"${VARNAME}\" to \"${${VARNAME}}\", in cache ${CACHETYPE}.")
endfunction()
function(CheckAndSetDefault VARNAME DEFAULT)
# cmessage(DEBUG "Trying to assign variable ${VARNAME}.")
if(NOT DEFINED ${VARNAME})
set(${VARNAME} ${DEFAULT} PARENT_SCOPE)
set(${VARNAME} ${DEFAULT})
endif()
cmessage(CACHE "--Set variable: \"${VARNAME}\" to \"${${VARNAME}}\".")
endfunction()
CheckAndSetDefaultCache(VERBOSE TRUE BOOL "Whether to configure loudly.")
set (CMAKE_SKIP_BUILD_RPATH TRUE)
#Changes default install path to be a subdirectory of the build dir.
#Can set build dir at configure time with -DCMAKE_INSTALL_PREFIX=/install/path
if(CMAKE_INSTALL_PREFIX STREQUAL "" OR CMAKE_INSTALL_PREFIX STREQUAL
"/usr/local")
set(CMAKE_INSTALL_PREFIX "${CMAKE_BINARY_DIR}/${CMAKE_SYSTEM_NAME}")
elseif(NOT DEFINED CMAKE_INSTALL_PREFIX)
set(CMAKE_INSTALL_PREFIX "${CMAKE_BINARY_DIR}/${CMAKE_SYSTEM_NAME}")
endif()
if(CMAKE_BUILD_TYPE STREQUAL "")
set(CMAKE_BUILD_TYPE DEBUG)
elseif(NOT DEFINED CMAKE_BUILD_TYPE)
set(CMAKE_BUILD_TYPE DEBUG)
endif()
CheckAndSetDefaultCache(EXTRA_SETUP_SCRIPT "" PATH "The path to an extra script to inject into the NUISANCE setup script. <>")
CheckAndSetDefaultCache(USE_MINIMIZER TRUE INTERNAL "Whether we are using the ROOT minimization libraries. ")
CheckAndSetDefaultCache(USE_ROOT6 FALSE INTERNAL "Whether we are using the ROOT 6. ")
CheckAndSetDefaultCache(USE_HEPMC FALSE BOOL "Whether to enable HepMC input support. ")
CheckAndSetDefaultEnv(HEPMC "" PATH "Path to HepMC source tree root directory. Overrides environment variable \$HEPMC <>" HEPMC)
CheckAndSetDefaultCache(HEPMC_MOMUNIT "GEV" STRING "HepMC momentum units [MEV|GEV]. ")
CheckAndSetDefaultCache(HEPMC_LENUNIT "CM" STRING "HepMC momentum units [MM|CM]. ")
CheckAndSetDefaultCache(HEPMC_USED_EP FALSE INTERNAL "Whether we built HepMC or not. ")
CheckAndSetDefaultCache(USE_NEUT FALSE BOOL "Whether to enable NEUT (reweight) support. Requires external libraries. ")
CheckAndSetDefaultEnv(NEUT_VERSION FALSE STRING "NEUT version string, e.g. 5.4.0. <5.4.0>" NEUT_VERSION)
CheckAndSetDefaultEnv(NEUT_ROOT "" PATH "Path to NEUT source tree root directory. Overrides environment variable \$NEUT_ROOT <>" NEUT_ROOT)
CheckAndSetDefaultEnv(CERN "" PATH "Path to CERNLIB source tree root directory that NEUT was built against. Overrides environment variable \$CERN <>" CERN)
CheckAndSetDefaultEnv(CERN_LEVEL "" STRING "CERNLIB Library version. Overrides environment variable \$CERN_LEVEL <>" CERN_LEVEL)
CheckAndSetDefaultCache(USE_NuWro FALSE BOOL "Whether to enable NuWro support. ")
CheckAndSetDefaultEnv(NUWRO "" PATH "Path to NuWro source tree root directory. Overrides environment variable \$NUWRO <>" NUWRO)
CheckAndSetDefaultEnv(NUWRO_INC "" PATH "Path to NuWro installed includes directory, needs to contain \"params_all.h\". Overrides environment variable \$NUWRO_INC <>" NUWRO_INC)
CheckAndSetDefaultCache(NUWRO_INPUT_FILE "" FILEPATH "Path to an input NuWro event vector, which can be used to build NuWro i/o libraries. <>")
CheckAndSetDefaultCache(NUWRO_BUILT_FROM_FILE FALSE INTERNAL "Whether the NuWro libraries were built by NUISANCE. ")
CheckAndSetDefaultCache(USE_NuWro_RW FALSE BOOL "Whether to try and build support for NuWro reweighting. ")
CheckAndSetDefaultCache(USE_NuWro_SRW_Event FALSE BOOL "Whether to use cut down NuWro reweight event format. Requires NuWro reweight. ")
CheckAndSetDefaultCache(USE_GENIE FALSE BOOL "Whether to enable GENIE (reweight) support. Requires external libraries. ")
CheckAndSetDefaultCache(GENIE_VERSION "AUTO" STRING "GENIE Version ")
CheckAndSetDefaultEnv(GENIE "" PATH "Path to GENIE source tree root directory. Overrides environment variable \$GENIE <>" GENIE)
CheckAndSetDefaultEnv(GENIE_REWEIGHT "" PATH "Path to GENIE ReWeight directory. Only relevant for GENIE v3+. Overrides environment variable \$GENIE_REWEIGHT <>" GENIE_REWEIGHT)
CheckAndSetDefaultCache(GENIE_EMPMEC_REWEIGHT FALSE BOOL "Whether to use GENIE EMP MEC reweight (requires custom GENIE) ")
CheckAndSetDefaultEnv(LHAPDF_LIB "" PATH "Path to pre-built LHAPDF libraries. Overrides environment variable \$LHAPDF_LIB. <>" LHAPDF_LIB)
CheckAndSetDefaultEnv(LHAPDF_INC "" PATH "Path to installed LHAPDF headers. Overrides environment variable \$LHAPDF_INC. <>" LHAPDF_INC)
CheckAndSetDefaultEnv(LHAPATH "" PATH "Path to LHA PDF inputs. Overrides environment variable \$LHAPATH. <>" LHAPATH)
CheckAndSetDefaultEnv(LIBXML2_LIB "" PATH "Path to pre-built LIBXML2 libraries. Overrides environment variable \$LIBXML2_LIB. <>" LIBXML2_LIB)
CheckAndSetDefaultEnv(LIBXML2_INC "" PATH "Path to installed LIBXML2 headers. Overrides environment variable \$LIBXML2_INC. <>" LIBXML2_INC)
CheckAndSetDefaultEnv(LOG4CPP_LIB "" PATH "Path to pre-built LOG4CPP libraries. Overrides environment variable \$LOG4CPP_LIB. <>" LOG4CPP_LIB)
CheckAndSetDefaultEnv(LOG4CPP_INC "" PATH "Path to installed LOG4CPP headers. Overrides environment variable \$LOG4CPP_INC. <>" LOG4CPP_INC)
CheckAndSetDefaultEnv(GSL_LIB "" PATH "Path to pre-built gsl libraries. Overrides environment variable \$GSL_LIB. <>" GSL_LIB)
CheckAndSetDefaultEnv(GSL_INC "" PATH "Path to installed gsl headers. Overrides environment variable \$GSL_INC. <>" GSL_INC)
CheckAndSetDefaultCache(BUILD_GEVGEN FALSE BOOL "Whether to build nuisance_gevgen app.")
CheckAndSetDefaultCache(USE_T2K FALSE BOOL "Whether to enable T2KReWeight support. Requires external libraries. ")
CheckAndSetDefaultEnv(T2KREWEIGHT "" PATH "Path to installed T2KREWEIGHTReWeight. Overrides environment variable \$T2KREWEIGHT. <>" T2KREWEIGHT)
CheckAndSetDefaultCache(USE_NIWG FALSE BOOL "Whether to enable (T2K) NIWG ReWeight support. Requires external libraries. ")
CheckAndSetDefaultEnv(NIWG_ROOT "" PATH "Path to installed NIWGReWeight. Overrides environment variable \$NIWG. <>" NIWG)
CheckAndSetDefaultCache(USE_MINERvA_RW FALSE BOOL "Whether to enable MINERvA ReWeight support. ")
-CheckAndSetDefaultCache(USE_NOvARwgt FALSE BOOL "Whether to enable MINERvA ReWeight support. ")
+CheckAndSetDefaultCache(USE_NOvARwgt FALSE BOOL "Whether to enable NOvA ReWeight support. ")
CheckAndSetDefaultEnv(NOVARWGT "" PATH "Path to directory containing libPythia6.so. Overrides environment variable \$NOVARWGT <>" NOVARWGT)
CheckAndSetDefaultEnv(PYTHIA6 "" PATH "Path to directory containing libPythia6.so. Overrides environment variable \$PYTHIA6 <>" PYTHIA6)
CheckAndSetDefaultEnv(PYTHIA8 "" PATH "Path to directory containing libPythia8.so. Overrides environment variable \$PYTHIA8 <>" PYTHIA8)
CheckAndSetDefaultCache(USE_PYTHIA8 FALSE BOOL "Whether to enable PYTHIA8 event support. ")
CheckAndSetDefaultCache(USE_GiBUU TRUE BOOL "Whether to enable GiBUU event support. ")
CheckAndSetDefaultCache(BUILD_GiBUU FALSE BOOL "Whether to build supporting GiBUU event tools along with a patched version of GiBUU. ")
CheckAndSetDefaultCache(USE_NUANCE TRUE BOOL "Whether to enable NUANCE event support. ")
CheckAndSetDefaultCache(USE_PROB3PP FALSE BOOL "Whether to download and compile in Prob3++ support. ")
CheckAndSetDefaultCache(NO_EXTERNAL_UPDATE FALSE BOOL "Whether to perform the update target for external dependencies. Note this may produce errors for CMake < 3.8 where a bug was fixed for the feature that this option invokes. ")
CheckAndSetDefaultCache(USE_GPERFTOOLS FALSE BOOL "Whether to compile in google performance tools. ")
CheckAndSetDefault(NEED_PYTHIA6 FALSE)
CheckAndSetDefault(NEED_PYTHIA8 FALSE)
CheckAndSetDefault(NEED_ROOTEVEGEN FALSE)
CheckAndSetDefault(NEED_ROOTPYTHIA6 FALSE)
CheckAndSetDefaultCache(USE_OMP FALSE BOOL "Whether to enable multicore features (there currently are none...). ")
CheckAndSetDefaultCache(USE_DYNSAMPLES TRUE BOOL "Whether to enable the dynamic sample loader. ")
CheckAndSetDefault(NO_EXPERIMENTS FALSE)
cmessage(STATUS "NO_EXPERIMENTS: ${NO_EXPERIMENTS}")
CheckAndSetDefaultCache(NO_ANL ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build ANL samples. <-DNO_EXPERIMENTS=FALSE>")
CheckAndSetDefaultCache(NO_ArgoNeuT ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build ArgoNeuT samples. <-DNO_EXPERIMENTS=FALSE>")
CheckAndSetDefaultCache(NO_BEBC ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build BEBC samples. <-DNO_EXPERIMENTS=FALSE>")
CheckAndSetDefaultCache(NO_BNL ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build BNL samples. <-DNO_EXPERIMENTS=FALSE>")
CheckAndSetDefaultCache(NO_FNAL ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build FNAL samples. <-DNO_EXPERIMENTS=FALSE>")
CheckAndSetDefaultCache(NO_GGM ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build GGM samples. <-DNO_EXPERIMENTS=FALSE>")
CheckAndSetDefaultCache(NO_K2K ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build K2K samples. <-DNO_EXPERIMENTS=FALSE>")
CheckAndSetDefaultCache(NO_MINERvA ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build MINERvA samples. <-DNO_EXPERIMENTS=FALSE>")
CheckAndSetDefaultCache(NO_MiniBooNE ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build MiniBooNE samples. <-DNO_EXPERIMENTS=FALSE>")
CheckAndSetDefaultCache(NO_T2K ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build T2K samples. <-DNO_EXPERIMENTS=FALSE>")
CheckAndSetDefaultCache(NO_SciBooNE ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build SciBooNE samples. <-DNO_EXPERIMENTS=FALSE>")
function(SAYVARS)
LIST(APPEND VARS
USE_HEPMC
HEPMC
HEPMC_MOMUNIT
HEPMC_LENUNIT
HEPMC_USED_EP
USE_NEUT
+ NEUT_VERSION
NEUT_ROOT
CERN
CERN_LEVEL
USE_NuWro
NUWRO
NUWRO_INC
NUWRO_INPUT_FILE
NUWRO_BUILT_FROM_FILE
USE_GENIE
+ GENIE_VERSION
GENIE
+ GENIE_REWEIGHT
LHAPDF_LIB
LHAPDF_INC
+ LHAPATH
LIBXML2_LIB
LIBXML2_INC
LOG4CPP_LIB
- GENIE_LOG4CPP_INC
+ LOG4CPP_INC
+ GSL_LIB
+ GSL_INC
+ PYTHIA6
+ PYTHIA8
+ USE_PYTHIA8
BUILD_GEVGEN
USE_T2K
+ T2KREWEIGHT
USE_NIWG
+ NIWG_ROOT
+ USE_MINERvA_RW
+ USE_NOvARwgt
+ NOVARWGT
USE_GiBUU
BUILD_GiBUU
USE_NUANCE
+ USE_PROB3PP
NO_EXTERNAL_UPDATE
USE_GPERFTOOLS
+ NO_EXPERIMENTS
NO_ANL
NO_ArgoNeuT
NO_BEBC
NO_BNL
NO_FNAL
NO_GGM
NO_K2K
NO_MINERvA
NO_MiniBooNE
NO_T2K
NO_SciBooNE)
foreach(v ${VARS})
if(DEFINED ${v})
cmessage(DEBUG "VARIABLE: \"${v}\" = \"${${v}}\"")
endif()
endforeach(v)
endfunction()
diff --git a/parameters/dial_conversion.card b/parameters/dial_conversion.card
index 0227a07..7a50e83 100644
--- a/parameters/dial_conversion.card
+++ b/parameters/dial_conversion.card
@@ -1,38 +1,46 @@
# par Name Units Nominal FracErr ConvFunc
-neut_parameter MaCCQE GeV 1.21*(1.0+x*0.16)
+neut_parameter MaCCQE GeV 1.21*(1.0+x*0.16)
neut_parameter MaNFFRES GeV 0.95*(1.0+x*0.157894737)
neut_parameter CA5RES x 1.01*(1.0+x*0.247524752)
neut_parameter BgSclRES x 1.30*(1.0+x*0.153846154)
neut_parameter FrAbs_pi x 1.1*(1.0+x*0.5)
-neut_parameter FrInelLow_pi x 1*(1.0+x*0.5)
+neut_parameter FrInelLow_pi x 1.0+x*0.5
neut_parameter FrInelHigh_pi x 1.8*(1.0+x*0.3)
-neut_parameter FrPiProd_pi x 1*(1.0+x*0.5)
-neut_parameter FrCExLow_pi x 1*(1.0+x*0.5)
+neut_parameter FrPiProd_pi x 1.0+x*0.5
+neut_parameter FrCExLow_pi x 1.0+x*0.5
neut_parameter FrCExHigh_pi x 1.8*(1.0+x*0.3)
niwg_parameter NIWGMEC_Norm_C12 % 100.0*(1.0+x)
niwg_parameter VecFFCCQE MDLQE x
niwg_parameter CCQEFermiSurfMom MeV 217*(1.0+x*0.15)
t2k_parameter NIWG2014a_pF_C12 MeV 217*(1.0+x*0.073733)
t2k_parameter NIWG2014a_pF_O16 MeV 225*(1.0+x*0.071111)
t2k_parameter NIWGMEC_PDDWeight_C12 x 0.5*(1.0+x*1.0)
t2k_parameter NIWGMEC_PDDWeight_O16 x 0.5*(1.0+x*1.0)
-t2k_parameter NXSec_MaCCQE GeV 1.21*(1.0+x*0.165289256)
+t2k_parameter NIWG2012a_dismpishp x 1.0*(1.0+x*0.4)
+t2k_parameter NXSec_MaCCQE GeV 1.21*(1.0+x*0.165289256)
t2k_parameter NXSec_MaNFFRES GeV 0.95*(1.0+x*0.157894737)
t2k_parameter NXSec_CA5RES x 1.01*(1.0+x*0.247524752)
-t2k_parameter NXSec_BgSclCCRES x 1.30*(1.0+x*0.153846154)
+t2k_parameter NXSec_BgSclRES x 1.30*(1.0+x*0.153846154)
+
+t2k_parameter NCasc_FrAbs_pi x 1.1*(1.0+x*0.5)
+t2k_parameter NCasc_FrInelLow_pi x 1.0+x*0.5
+t2k_parameter NCasc_FrInelHigh_pi x 1.8*(1.0+x*0.3)
+t2k_parameter NCasc_FrPiProd_pi x 1.0+x*0.5
+t2k_parameter NCasc_FrCExLow_pi x 1.0+x*0.5
+t2k_parameter NCasc_FrCExHigh_pi x 1.8*(1.0+x*0.3)
t2k_parameter NIWG_Effective_rpaCCQE_A x 1.0*(1.0 + x*0.1)
t2k_parameter NIWG_Effective_rpaCCQE_B x 1.0*(1.0 + x*0.1)
t2k_parameter NIWG_Effective_rpaCCQE_C x 1.0*(1.0 + x*0.1)
t2k_parameter NIWG_Effective_rpaCCQE_D x 1.0*(1.0 + x*0.1)
t2k_parameter NIWG_Effective_rpaCCQE_U x 1.0*(1.0 + x*0.1)
nuwro_parameter kNuwro_Ma_CCQE MeV 1200*(1.0+x*0.160)
nuwro_parameter kNuwro_MaRES GeV 0.940*(1.0+x*0.1)
nuwro_parameter kNuwro_CA5 x 1.19*(1.0+x*0.1)
nuwro_parameter kNuwro_SPPBkgScale x 1.0*(1.0+x*0.26)
diff --git a/src/MCStudies/GenericFlux_Vectors.cxx b/src/MCStudies/GenericFlux_Vectors.cxx
index 045f0c0..0543dc5 100644
--- a/src/MCStudies/GenericFlux_Vectors.cxx
+++ b/src/MCStudies/GenericFlux_Vectors.cxx
@@ -1,339 +1,363 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#include "GenericFlux_Vectors.h"
GenericFlux_Vectors::GenericFlux_Vectors(std::string name,
std::string inputfile, FitWeight *rw,
std::string type,
std::string fakeDataFile) {
// Measurement Details
fName = name;
eventVariables = NULL;
// Define our energy range for flux calcs
EnuMin = 0.;
EnuMax = 1E10; // Arbritrarily high energy limit
if (Config::HasPar("EnuMin")) {
EnuMin = Config::GetParD("EnuMin");
}
if (Config::HasPar("EnuMax")) {
EnuMax = Config::GetParD("EnuMax");
}
// Set default fitter flags
fIsDiag = true;
fIsShape = false;
fIsRawEvents = false;
// This function will sort out the input files automatically and parse all the
// inputs,flags,etc.
// There may be complex cases where you have to do this by hand, but usually
// this will do.
Measurement1D::SetupMeasurement(inputfile, type, rw, fakeDataFile);
eventVariables = NULL;
// Setup fDataHist as a placeholder
this->fDataHist = new TH1D(("empty_data"), ("empty-data"), 1, 0, 1);
this->SetupDefaultHist();
fFullCovar = StatUtils::MakeDiagonalCovarMatrix(fDataHist);
covar = StatUtils::GetInvert(fFullCovar);
// 1. The generator is organised in SetupMeasurement so it gives the
// cross-section in "per nucleon" units.
// So some extra scaling for a specific measurement may be required. For
// Example to get a "per neutron" measurement on carbon
// which we do here, we have to multiple by the number of nucleons 12 and
// divide by the number of neutrons 6.
// N.B. MeasurementBase::PredictedEventRate includes the 1E-38 factor that is
// often included here in other classes that directly integrate the event
// histogram. This method is used here as it now respects EnuMin and EnuMax
// correctly.
this->fScaleFactor =
(this->PredictedEventRate("width", 0, EnuMax) / double(fNEvents)) /
- this->TotalIntegratedFlux();
+ this->TotalIntegratedFlux("width");
LOG(SAM) << " Generic Flux Scaling Factor = " << fScaleFactor
<< " [= " << (GetEventHistogram()->Integral("width") * 1E-38) << "/("
- << (fNEvents + 0.) << "*" << this->TotalIntegratedFlux() << ")]"
+ << (fNEvents + 0.) << "*" << TotalIntegratedFlux("width") << ")]"
<< std::endl;
-
if (fScaleFactor <= 0.0) {
ERR(WRN) << "SCALE FACTOR TOO LOW " << std::endl;
throw;
}
// Setup our TTrees
this->AddEventVariablesToTree();
this->AddSignalFlagsToTree();
}
void GenericFlux_Vectors::AddEventVariablesToTree() {
// Setup the TTree to save everything
if (!eventVariables) {
Config::Get().out->cd();
eventVariables = new TTree((this->fName + "_VARS").c_str(),
(this->fName + "_VARS").c_str());
}
LOG(SAM) << "Adding Event Variables" << std::endl;
eventVariables->Branch("Mode", &Mode, "Mode/I");
eventVariables->Branch("cc", &cc, "cc/B");
eventVariables->Branch("PDGnu", &PDGnu, "PDGnu/I");
eventVariables->Branch("Enu_true", &Enu_true, "Enu_true/F");
eventVariables->Branch("tgt", &tgt, "tgt/I");
eventVariables->Branch("PDGLep", &PDGLep, "PDGLep/I");
eventVariables->Branch("ELep", &ELep, "ELep/F");
eventVariables->Branch("CosLep", &CosLep, "CosLep/F");
// Basic interaction kinematics
eventVariables->Branch("Q2", &Q2, "Q2/F");
eventVariables->Branch("q0", &q0, "q0/F");
eventVariables->Branch("q3", &q3, "q3/F");
eventVariables->Branch("Enu_QE", &Enu_QE, "Enu_QE/F");
eventVariables->Branch("Q2_QE", &Q2_QE, "Q2_QE/F");
eventVariables->Branch("W_nuc_rest", &W_nuc_rest, "W_nuc_rest/F");
eventVariables->Branch("W", &W, "W/F");
+ eventVariables->Branch("W_genie", &W_genie, "W_genie/F");
eventVariables->Branch("x", &x, "x/F");
eventVariables->Branch("y", &y, "y/F");
eventVariables->Branch("Eav", &Eav, "Eav/F");
eventVariables->Branch("EavAlt", &EavAlt, "EavAlt/F");
+ eventVariables->Branch("dalphat", &dalphat, "dalphat/F");
+ eventVariables->Branch("dpt", &dpt, "dpt/F");
+ eventVariables->Branch("dphit", &dphit, "dphit/F");
+ eventVariables->Branch("pnreco_C", &pnreco_C, "pnreco_C/F");
+
// Save outgoing particle vectors
eventVariables->Branch("nfsp", &nfsp, "nfsp/I");
eventVariables->Branch("px", px, "px[nfsp]/F");
eventVariables->Branch("py", py, "py[nfsp]/F");
eventVariables->Branch("pz", pz, "pz[nfsp]/F");
eventVariables->Branch("E", E, "E[nfsp]/F");
eventVariables->Branch("pdg", pdg, "pdg[nfsp]/I");
// Event Scaling Information
eventVariables->Branch("Weight", &Weight, "Weight/F");
eventVariables->Branch("InputWeight", &InputWeight, "InputWeight/F");
eventVariables->Branch("RWWeight", &RWWeight, "RWWeight/F");
// Should be a double because may be 1E-39 and less
eventVariables->Branch("fScaleFactor", &fScaleFactor, "fScaleFactor/D");
// The customs
eventVariables->Branch("CustomWeight", &CustomWeight, "CustomWeight/F");
eventVariables->Branch("CustomWeightArray", CustomWeightArray, "CustomWeightArray[6]/F");
return;
}
void GenericFlux_Vectors::FillEventVariables(FitEvent *event) {
ResetVariables();
// Fill Signal Variables
FillSignalFlags(event);
LOG(DEB) << "Filling signal" << std::endl;
// Now fill the information
Mode = event->Mode;
cc = (abs(event->Mode) < 30);
// Get the incoming neutrino and outgoing lepton
FitParticle *nu = event->GetNeutrinoIn();
FitParticle *lep = event->GetHMFSAnyLepton();
PDGnu = nu->fPID;
Enu_true = nu->fP.E() / 1E3;
tgt = event->fTargetPDG;
if (lep != NULL) {
PDGLep = lep->fPID;
ELep = lep->fP.E() / 1E3;
CosLep = cos(nu->fP.Vect().Angle(lep->fP.Vect()));
// Basic interaction kinematics
Q2 = -1 * (nu->fP - lep->fP).Mag2() / 1E6;
q0 = (nu->fP - lep->fP).E() / 1E3;
q3 = (nu->fP - lep->fP).Vect().Mag() / 1E3;
// These assume C12 binding from MINERvA... not ideal
Enu_QE = FitUtils::EnuQErec(lep->fP, CosLep, 34., true);
Q2_QE = FitUtils::Q2QErec(lep->fP, CosLep, 34., true);
Eav = FitUtils::GetErecoil_MINERvA_LowRecoil(event)/1.E3;
EavAlt = FitUtils::Eavailable(event)/1.E3;
// Get W_true with assumption of initial state nucleon at rest
float m_n = (float)PhysConst::mass_proton;
// Q2 assuming nucleon at rest
W_nuc_rest = sqrt(-Q2 + 2 * m_n * q0 + m_n * m_n);
// True Q2
W = sqrt(-Q2 + 2 * m_n * q0 + m_n * m_n);
x = Q2 / (2 * m_n * q0);
y = 1 - ELep / Enu_true;
+
+ dalphat = FitUtils::Get_STV_dalphat(event, PDGnu, true);
+ dpt = FitUtils::Get_STV_dpt(event, PDGnu, true);
+ dphit = FitUtils::Get_STV_dphit(event, PDGnu, true);
+ pnreco_C = FitUtils::Get_pn_reco_C(event, PDGnu, true);
}
// Loop over the particles and store all the final state particles in a vector
for (UInt_t i = 0; i < event->Npart(); ++i) {
bool part_alive = event->PartInfo(i)->fIsAlive &&
- event->PartInfo(i)->Status() == kFinalState;
+ event->PartInfo(i)->Status() == kFinalState;
if (!part_alive) continue;
partList.push_back(event->PartInfo(i));
}
// Save outgoing particle vectors
nfsp = (int)partList.size();
for (int i = 0; i < nfsp; ++i) {
px[i] = partList[i]->fP.X() / 1E3;
py[i] = partList[i]->fP.Y() / 1E3;
pz[i] = partList[i]->fP.Z() / 1E3;
E[i] = partList[i]->fP.E() / 1E3;
pdg[i] = partList[i]->fPID;
}
+#ifdef __GENIE_ENABLED__
+ if (event->fType == kGENIE) {
+ EventRecord * gevent = static_cast(event->genie_event->event);
+ const Interaction * interaction = gevent->Summary();
+ const Kinematics & kine = interaction->Kine();
+ double W_genie = kine.W();
+ }
+#endif
+
// Fill event weights
Weight = event->RWWeight * event->InputWeight;
RWWeight = event->RWWeight;
InputWeight = event->InputWeight;
// And the Customs
CustomWeight = event->CustomWeight;
for (int i = 0; i < 6; ++i) {
CustomWeightArray[i] = event->CustomWeightArray[i];
}
// Fill the eventVariables Tree
eventVariables->Fill();
return;
};
//********************************************************************
void GenericFlux_Vectors::ResetVariables() {
-//********************************************************************
+ //********************************************************************
cc = false;
// Reset all Function used to extract any variables of interest to the event
Mode = PDGnu = tgt = PDGLep = 0;
Enu_true = ELep = CosLep = Q2 = q0 = q3 = Enu_QE = Q2_QE = W_nuc_rest = W = x = y = Eav = EavAlt = -999.9;
+ W_genie = -999;
+ // Other fun variables
+ // MINERvA-like ones
+ dalphat = dpt = dphit = pnreco_C = -999.99;
+
nfsp = 0;
for (int i = 0; i < kMAX; ++i){
px[i] = py[i] = pz[i] = E[i] = -999;
pdg[i] = 0;
}
Weight = InputWeight = RWWeight = 0.0;
CustomWeight = 0.0;
for (int i = 0; i < 6; ++i) CustomWeightArray[i] = 0.0;
partList.clear();
flagCCINC = flagNCINC = flagCCQE = flagCC0pi = flagCCQELike = flagNCEL = flagNC0pi = flagCCcoh = flagNCcoh = flagCC1pip = flagNC1pip = flagCC1pim = flagNC1pim = flagCC1pi0 = flagNC1pi0 = false;
}
//********************************************************************
void GenericFlux_Vectors::FillSignalFlags(FitEvent *event) {
//********************************************************************
// Some example flags are given from SignalDef.
// See src/Utils/SignalDef.cxx for more.
int nuPDG = event->PartInfo(0)->fPID;
// Generic signal flags
flagCCINC = SignalDef::isCCINC(event, nuPDG);
flagNCINC = SignalDef::isNCINC(event, nuPDG);
flagCCQE = SignalDef::isCCQE(event, nuPDG);
flagCCQELike = SignalDef::isCCQELike(event, nuPDG);
flagCC0pi = SignalDef::isCC0pi(event, nuPDG);
flagNCEL = SignalDef::isNCEL(event, nuPDG);
flagNC0pi = SignalDef::isNC0pi(event, nuPDG);
flagCCcoh = SignalDef::isCCCOH(event, nuPDG, 211);
flagNCcoh = SignalDef::isNCCOH(event, nuPDG, 111);
flagCC1pip = SignalDef::isCC1pi(event, nuPDG, 211);
flagNC1pip = SignalDef::isNC1pi(event, nuPDG, 211);
flagCC1pim = SignalDef::isCC1pi(event, nuPDG, -211);
flagNC1pim = SignalDef::isNC1pi(event, nuPDG, -211);
flagCC1pi0 = SignalDef::isCC1pi(event, nuPDG, 111);
flagNC1pi0 = SignalDef::isNC1pi(event, nuPDG, 111);
}
void GenericFlux_Vectors::AddSignalFlagsToTree() {
if (!eventVariables) {
Config::Get().out->cd();
eventVariables = new TTree((this->fName + "_VARS").c_str(),
- (this->fName + "_VARS").c_str());
+ (this->fName + "_VARS").c_str());
}
LOG(SAM) << "Adding signal flags" << std::endl;
// Signal Definitions from SignalDef.cxx
eventVariables->Branch("flagCCINC", &flagCCINC, "flagCCINC/O");
eventVariables->Branch("flagNCINC", &flagNCINC, "flagNCINC/O");
eventVariables->Branch("flagCCQE", &flagCCQE, "flagCCQE/O");
eventVariables->Branch("flagCC0pi", &flagCC0pi, "flagCC0pi/O");
eventVariables->Branch("flagCCQELike", &flagCCQELike, "flagCCQELike/O");
eventVariables->Branch("flagNCEL", &flagNCEL, "flagNCEL/O");
eventVariables->Branch("flagNC0pi", &flagNC0pi, "flagNC0pi/O");
eventVariables->Branch("flagCCcoh", &flagCCcoh, "flagCCcoh/O");
eventVariables->Branch("flagNCcoh", &flagNCcoh, "flagNCcoh/O");
eventVariables->Branch("flagCC1pip", &flagCC1pip, "flagCC1pip/O");
eventVariables->Branch("flagNC1pip", &flagNC1pip, "flagNC1pip/O");
eventVariables->Branch("flagCC1pim", &flagCC1pim, "flagCC1pim/O");
eventVariables->Branch("flagNC1pim", &flagNC1pim, "flagNC1pim/O");
eventVariables->Branch("flagCC1pi0", &flagCC1pi0, "flagCC1pi0/O");
eventVariables->Branch("flagNC1pi0", &flagNC1pi0, "flagNC1pi0/O");
};
void GenericFlux_Vectors::Write(std::string drawOpt) {
// First save the TTree
eventVariables->Write();
// Save Flux and Event Histograms too
GetInput()->GetFluxHistogram()->Write();
GetInput()->GetEventHistogram()->Write();
return;
}
// Override functions which aren't really necessary
bool GenericFlux_Vectors::isSignal(FitEvent *event) {
(void)event;
return true;
};
void GenericFlux_Vectors::ScaleEvents() { return; }
void GenericFlux_Vectors::ApplyNormScale(float norm) {
this->fCurrentNorm = norm;
return;
}
void GenericFlux_Vectors::FillHistograms() { return; }
void GenericFlux_Vectors::ResetAll() {
eventVariables->Reset();
return;
}
float GenericFlux_Vectors::GetChi2() { return 0.0; }
diff --git a/src/MCStudies/GenericFlux_Vectors.h b/src/MCStudies/GenericFlux_Vectors.h
index 921b247..101335b 100644
--- a/src/MCStudies/GenericFlux_Vectors.h
+++ b/src/MCStudies/GenericFlux_Vectors.h
@@ -1,127 +1,133 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#ifndef GenericFlux_Vectors_H_SEEN
#define GenericFlux_Vectors_H_SEEN
#include "Measurement1D.h"
+#include "FitEvent.h"
class GenericFlux_Vectors : public Measurement1D {
public:
GenericFlux_Vectors(std::string name, std::string inputfile, FitWeight *rw, std::string type, std::string fakeDataFile);
virtual ~GenericFlux_Vectors() {};
//! Grab info from event
void FillEventVariables(FitEvent *event);
//! Fill signal flags
void FillSignalFlags(FitEvent *event);
void ResetVariables();
//! Fill Custom Histograms
void FillHistograms();
//! ResetAll
void ResetAll();
//! Scale
void ScaleEvents();
//! Norm
void ApplyNormScale(float norm);
//! Define this samples signal
bool isSignal(FitEvent *nvect);
//! Write Files
void Write(std::string drawOpt);
//! Get Chi2
float GetChi2();
void AddEventVariablesToTree();
void AddSignalFlagsToTree();
private:
TTree* eventVariables;
std::vector partList;
int Mode;
bool cc;
int PDGnu;
int tgt;
int PDGLep;
float ELep;
float CosLep;
// Basic interaction kinematics
float Q2;
float q0;
float q3;
float Enu_QE;
float Enu_true;
float Q2_QE;
float W_nuc_rest;
float W;
float x;
float y;
float Eav;
float EavAlt;
+ float dalphat;
+ float W_genie;
+ float dpt;
+ float dphit;
+ float pnreco_C;
// Save outgoing particle vectors
int nfsp;
static const int kMAX = 200;
float px[kMAX];
float py[kMAX];
float pz[kMAX];
float E[kMAX];
int pdg[kMAX];
// Basic event info
float Weight;
float InputWeight;
float RWWeight;
double fScaleFactor;
// Custom weights
float CustomWeight;
float CustomWeightArray[6];
// Generic signal flags
bool flagCCINC;
bool flagNCINC;
bool flagCCQE;
bool flagCC0pi;
bool flagCCQELike;
bool flagNCEL;
bool flagNC0pi;
bool flagCCcoh;
bool flagNCcoh;
bool flagCC1pip;
bool flagNC1pip;
bool flagCC1pim;
bool flagNC1pim;
bool flagCC1pi0;
bool flagNC1pi0;
};
#endif
diff --git a/src/Reweight/BeRPA.h b/src/Reweight/BeRPA.h
new file mode 100644
index 0000000..d6e85c8
--- /dev/null
+++ b/src/Reweight/BeRPA.h
@@ -0,0 +1,93 @@
+#ifndef _ERPA_H_SEEN_
+#define _ERPA_H_SEEN_
+
+#include
+#include
+
+// =====================
+//
+// Roughly a copy/paste job from NIWGReWeight/NIWGReWeightEffectiveRPA.cc/h
+// Written by Asmita Redij, University of Bern
+// Implemented in MaCh3 by Clarence Wret: c.wret14@imperial.ac.uk
+//
+// =====================
+//
+// Purpose:
+//
+// eRPA attempts to mimic Nieves RPA calculation in functional form (cubic, switching to exponential)
+// Currently, there are no RPA shape parameters to include as RPA systematics
+// eRPA remedies this by ffitting a functional form to the Nieves calculation
+//
+// Defaults have been fit to Nieves RPA calculation
+// +/- 1 sigma in eRPA are also set after a fit to the +/- 1 sigma from Nieves RPA
+//
+// The eRPA calculation is a function of Q2 and return a RPA/no-RPA correction factor for a given value of Q2
+// It should currently only be applied to CCQE (or 2p2h?) interactions on nuclear targets (i.e. ignore H interactions!)
+//
+//
+// The eRPA calculation is a function which depends on five parameters. The naming of the parameters has changed as work
+// has developed on eRPA. Originally, eRPA was discussed in terms of five parameters: A, B, C, D, and U. The interpretation
+// of these parameters is:
+// Polynomial terms:
+// A = RPA correction at Q2 = 0
+// B = Slope of RPA correction at Q2 = 0
+// C = RPA correction when polynominal switches to exponential (when Q2 = U)
+// Exponential terms:
+// D = Controls the damping constant of the exponential when Q2 > U
+// Control terms:
+// U = 1.20, the value of Q2 where we switch from polynomial to exponential.
+// Can be varied but is in Jan 2016 not recommended use.
+//
+// We have now changed to describing eRPA in terms of the five parameters A, B, D, E, and U. The formula is identical
+// to the one which depends on A, B, C, D, and U. The relation between the new and old parameters is:
+// A (old) = A (new)
+// B (old) = B (new)
+// C (old) = D (new)
+// D (old) = E (new)
+// U (old) = U (new)
+//
+// To make thing extra confusing, the new parameterisation also includes a parameter C given by
+// C = D+U*E*(D-1)/3
+// This is *not* related to C in the old parameterisation, and is *not* directly used as a parameter by us (it's indirectly
+// determined by the values of D, U, and E)
+//
+// More info about the eRPA parameterisation (and anything else you could possibly want to know!) is in T2K-TN-309.
+// Some more info about eRPA is also available in the slides at http://www.t2k.org/asg/oagroup/meeting/2016/banffoa27seppm/ffsanderpa/
+// (beware: they use the "old" A,B,C,D,U parameter names).
+//
+// =====================
+
+// Calculates the absolute eRPA for given Q2, BERNSTEIN STYLE
+inline double const calcRPA(const double Q2, const double A, const double B, const double D, const double E, const double U = 1.20) {
+
+ // Callum's eye-balled nominals
+ // A = 0.6
+ // B = 1.0
+ // D = 1.2
+ // E = 1.0
+ // U = 1.2
+ //
+ // Callum's fitted nominals
+ // A = 0.59 +/- 20%
+ // B = 1.05 +/- 20%
+ // D = 1.13 +/- 15%
+ // E = 0.88 +/- 40%
+ // U = 1.2
+
+ // Kept for convenience
+ double eRPA = 1.;
+
+ // Q2 transition; less than U -> polynominal
+ if (Q2 < U) {
+ // xprime as prescribed by Callum
+ const double xprime = Q2/U;
+ const double C = D + U*E*(D-1)/3.;
+ eRPA = A*(1-xprime)*(1-xprime)*(1-xprime) + 3*B*(1-xprime)*(1-xprime)*xprime + 3*C*(1-xprime)*xprime*xprime + D*xprime*xprime*xprime;
+ } else {
+ eRPA = 1 + (D-1)*exp(-E*(Q2-U));
+ }
+
+ return eRPA;
+};
+
+#endif
diff --git a/src/Reweight/CMakeLists.txt b/src/Reweight/CMakeLists.txt
index 5ca22b8..00f4bda 100644
--- a/src/Reweight/CMakeLists.txt
+++ b/src/Reweight/CMakeLists.txt
@@ -1,91 +1,92 @@
# Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
################################################################################
# This file is part of NUISANCE.
#
# NUISANCE is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# NUISANCE is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with NUISANCE. If not, see .
################################################################################
set(IMPLFILES
GlobalDialList.cxx
FitWeight.cxx
WeightEngineBase.cxx
NEUTWeightEngine.cxx
NuWroWeightEngine.cxx
GENIEWeightEngine.cxx
WeightUtils.cxx
SampleNormEngine.cxx
LikelihoodWeightEngine.cxx
SplineWeightEngine.cxx
NUISANCESyst.cxx
T2KWeightEngine.cxx
NUISANCEWeightEngine.cxx
NUISANCEWeightCalcs.cxx
NIWGWeightEngine.cxx
OscWeightEngine.cxx
MINERvAWeightCalcs.cxx
weightRPA.h
)
if(USE_NOvARwgt)
LIST(APPEND IMPLFILES NOvARwgtEngine.cxx)
endif()
set(HEADERFILES
GlobalDialList.h
FitWeight.h
WeightEngineBase.h
NEUTWeightEngine.h
NuWroWeightEngine.h
GENIEWeightEngine.h
WeightUtils.h
SampleNormEngine.h
LikelihoodWeightEngine.h
SplineWeightEngine.h
NUISANCESyst.h
T2KWeightEngine.h
NUISANCEWeightEngine.h
NUISANCEWeightCalcs.h
NIWGWeightEngine.h
OscWeightEngine.h
MINERvAWeightCalcs.h
weightRPA.h
+BeRPA.h
)
if(USE_NOvARwgt)
LIST(APPEND HEADERFILES NOvARwgtEngine.h)
endif()
set(LIBNAME Reweight)
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
"${ExtFit_VERSION_MAJOR}.${ExtFit_VERSION_MINOR}.${ExtFit_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/Reweight/FitWeight.cxx b/src/Reweight/FitWeight.cxx
index 863d143..2794e7a 100644
--- a/src/Reweight/FitWeight.cxx
+++ b/src/Reweight/FitWeight.cxx
@@ -1,284 +1,288 @@
#include "FitWeight.h"
#include "GENIEWeightEngine.h"
#include "LikelihoodWeightEngine.h"
#include "ModeNormEngine.h"
#include "NEUTWeightEngine.h"
#include "NIWGWeightEngine.h"
#include "NUISANCEWeightEngine.h"
#include "NuWroWeightEngine.h"
#include "OscWeightEngine.h"
#include "SampleNormEngine.h"
#include "SplineWeightEngine.h"
#include "T2KWeightEngine.h"
#ifdef __NOVA_ENABLED__
#include "NOvARwgtEngine.h"
#endif
void FitWeight::AddRWEngine(int type) {
LOG(FIT) << "Adding reweight engine " << type << std::endl;
switch (type) {
case kNEUT:
fAllRW[type] = new NEUTWeightEngine("neutrw");
break;
case kNUWRO:
fAllRW[type] = new NuWroWeightEngine("nuwrorw");
break;
case kGENIE:
fAllRW[type] = new GENIEWeightEngine("genierw");
break;
case kNORM:
fAllRW[type] = new SampleNormEngine("normrw");
break;
case kLIKEWEIGHT:
fAllRW[type] = new LikelihoodWeightEngine("likerw");
break;
case kT2K:
fAllRW[type] = new T2KWeightEngine("t2krw");
break;
case kCUSTOM:
fAllRW[type] = new NUISANCEWeightEngine("nuisrw");
break;
case kSPLINEPARAMETER:
fAllRW[type] = new SplineWeightEngine("splinerw");
break;
case kNIWG:
fAllRW[type] = new NIWGWeightEngine("niwgrw");
break;
case kOSCILLATION:
fAllRW[type] = new OscWeightEngine();
break;
case kMODENORM:
fAllRW[type] = new ModeNormEngine();
break;
+#ifdef __NOVA_ENABLED__
case kNOvARWGT:
fAllRW[type] = new NOvARwgtEngine();
break;
+#endif
default:
THROW("CANNOT ADD RW Engine for unknown dial type: " << type);
break;
}
}
WeightEngineBase *FitWeight::GetRWEngine(int type) {
if (HasRWEngine(type)) {
return fAllRW[type];
}
THROW("CANNOT get RW Engine for dial type: " << type);
}
bool FitWeight::HasRWEngine(int type) {
switch (type) {
case kNEUT:
case kNUWRO:
case kGENIE:
case kNORM:
case kLIKEWEIGHT:
case kT2K:
case kCUSTOM:
case kSPLINEPARAMETER:
case kNIWG:
case kOSCILLATION:
+#ifdef __NOVA_ENABLED__
case kNOvARWGT: {
return fAllRW.count(type);
}
+#endif
default: { THROW("CANNOT get RW Engine for dial type: " << type); }
}
}
void FitWeight::IncludeDial(std::string name, std::string type, double val) {
// Should register the dial here.
int typeenum = Reweight::ConvDialType(type);
IncludeDial(name, typeenum, val);
}
void FitWeight::IncludeDial(std::string name, int dialtype, double val) {
// Get the dial enum
int nuisenum = Reweight::ConvDial(name, dialtype);
if (nuisenum == -1) {
THROW("NUISENUM == " << nuisenum << " for " << name);
}
// Setup RW Engine Pointer
if (fAllRW.find(dialtype) == fAllRW.end()) {
AddRWEngine(dialtype);
}
WeightEngineBase *rw = fAllRW[dialtype];
// Include the dial
rw->IncludeDial(name, val);
// Set Dial Value
if (val != -9999.9) {
rw->SetDialValue(name, val);
}
// Sort Maps
fAllEnums[name] = nuisenum;
fAllValues[nuisenum] = val;
// Sort Lists
fNameList.push_back(name);
fEnumList.push_back(nuisenum);
fValueList.push_back(val);
}
void FitWeight::Reconfigure(bool silent) {
// Reconfigure all added RW engines
for (std::map::iterator iter = fAllRW.begin();
iter != fAllRW.end(); iter++) {
(*iter).second->Reconfigure(silent);
}
}
void FitWeight::SetDialValue(std::string name, double val) {
// Add extra check, if name not found look for one with name in it.
int nuisenum = fAllEnums[name];
SetDialValue(nuisenum, val);
}
// Allow for name aswell using GlobalList to determine sample name.
void FitWeight::SetDialValue(int nuisenum, double val) {
// Conv dial type
int dialtype = Reweight::GetDialType(nuisenum);
if (fAllRW.find(dialtype) == fAllRW.end()) {
THROW("Cannot find RW Engine for dialtype = "
<< dialtype << ", " << Reweight::RemoveDialType(nuisenum));
}
// Get RW Engine for this dial
fAllRW[dialtype]->SetDialValue(nuisenum, val);
fAllValues[nuisenum] = val;
// Update ValueList
for (size_t i = 0; i < fEnumList.size(); i++) {
if (fEnumList[i] == nuisenum) {
fValueList[i] = val;
}
}
}
void FitWeight::SetAllDials(const double *x, int n) {
for (size_t i = 0; i < (UInt_t)n; i++) {
int rwenum = fEnumList[i];
SetDialValue(rwenum, x[i]);
}
Reconfigure();
}
double FitWeight::GetDialValue(std::string name) {
// Add extra check, if name not found look for one with name in it.
int nuisenum = fAllEnums[name];
return GetDialValue(nuisenum);
}
double FitWeight::GetDialValue(int nuisenum) { return fAllValues[nuisenum]; }
int FitWeight::GetDialPos(std::string name) {
int rwenum = fAllEnums[name];
return GetDialPos(rwenum);
}
int FitWeight::GetDialPos(int nuisenum) {
for (size_t i = 0; i < fEnumList.size(); i++) {
if (fEnumList[i] == nuisenum) {
return i;
}
}
ERR(FTL) << "No Dial Found! " << std::endl;
throw;
return -1;
}
bool FitWeight::DialIncluded(std::string name) {
return (fAllEnums.find(name) != fAllEnums.end());
}
bool FitWeight::DialIncluded(int rwenum) {
return (fAllValues.find(rwenum) != fAllValues.end());
}
double FitWeight::CalcWeight(BaseFitEvt *evt) {
double rwweight = 1.0;
for (std::map::iterator iter = fAllRW.begin();
iter != fAllRW.end(); iter++) {
double w = (*iter).second->CalcWeight(evt);
rwweight *= w;
}
return rwweight;
}
void FitWeight::UpdateWeightEngine(const double *x) {
size_t count = 0;
for (std::vector::iterator iter = fEnumList.begin();
iter != fEnumList.end(); iter++) {
SetDialValue((*iter), x[count]);
count++;
}
}
void FitWeight::GetAllDials(double *x, int n) {
for (int i = 0; i < n; i++) {
x[i] = GetDialValue(fEnumList[i]);
}
}
// bool FitWeight::NeedsEventReWeight(const double* x) {
// bool haschange = false;
// size_t count = 0;
// // Compare old to new and decide if RW needed.
// for (std::vector::iterator iter = fEnumList.begin();
// iter != fEnumList.end(); iter++) {
// int nuisenum = (*iter);
// int type = (nuisenum / 1000) - (nuisenum % 1000);
// // Compare old to new
// double oldval = GetDialValue(nuisenum);
// double newval = x[count];
// if (oldval != newval) {
// if (fAllRW[type]->NeedsEventReWeight()) {
// haschange = true;
// }
// }
// count++;
// }
// return haschange;
// }
double FitWeight::GetSampleNorm(std::string name) {
if (name.empty())
return 1.0;
// Find norm dial
if (fAllEnums.find(name + "_norm") != fAllEnums.end()) {
if (fAllValues.find(fAllEnums[name + "_norm"]) != fAllValues.end()) {
return fAllValues[fAllEnums[name + "_norm"]];
} else {
return 1.0;
}
} else {
return 1.0;
}
}
void FitWeight::Print() {
LOG(REC) << "Fit Weight State: " << std::endl;
for (size_t i = 0; i < fNameList.size(); i++) {
LOG(REC) << " -> Par " << i << ". " << fNameList[i] << " " << fValueList[i]
<< std::endl;
}
}
diff --git a/src/Reweight/NUISANCESyst.cxx b/src/Reweight/NUISANCESyst.cxx
index 04c425c..d11114b 100644
--- a/src/Reweight/NUISANCESyst.cxx
+++ b/src/Reweight/NUISANCESyst.cxx
@@ -1,74 +1,80 @@
#include "NUISANCESyst.h"
int Reweight::ConvertNUISANCEDial(std::string type) {
for (int i = kUnkownNUISANCEDial + 1; i < kNUISANCEDial_LAST; i++) {
if (!type.compare(ConvNUISANCEDial(i).c_str())) {
return i;
}
}
return kUnkownNUISANCEDial;
};
std::string Reweight::ConvNUISANCEDial(int type) {
switch (type) {
case kGaussianCorr_CCQE_norm: { return "GaussianCorr_CCQE_norm"; }
case kGaussianCorr_CCQE_tilt: { return "GaussianCorr_CCQE_tilt"; }
case kGaussianCorr_CCQE_Pq0: { return "GaussianCorr_CCQE_Pq0"; }
case kGaussianCorr_CCQE_Wq0: { return "GaussianCorr_CCQE_Wq0"; }
case kGaussianCorr_CCQE_Pq3: { return "GaussianCorr_CCQE_Pq3"; }
case kGaussianCorr_CCQE_Wq3: { return "GaussianCorr_CCQE_Wq3"; }
case kGaussianCorr_2p2h_norm: { return "GaussianCorr_2p2h_norm"; }
case kGaussianCorr_2p2h_tilt: { return "GaussianCorr_2p2h_tilt"; }
case kGaussianCorr_2p2h_Pq0: { return "GaussianCorr_2p2h_Pq0"; }
case kGaussianCorr_2p2h_Wq0: { return "GaussianCorr_2p2h_Wq0"; }
case kGaussianCorr_2p2h_Pq3: { return "GaussianCorr_2p2h_Pq3"; }
case kGaussianCorr_2p2h_Wq3: { return "GaussianCorr_2p2h_Wq3"; }
case kGaussianCorr_2p2h_PPandNN_norm: { return "GaussianCorr_2p2h_PPandNN_norm"; }
case kGaussianCorr_2p2h_PPandNN_tilt: { return "GaussianCorr_2p2h_PPandNN_tilt"; }
case kGaussianCorr_2p2h_PPandNN_Pq0: { return "GaussianCorr_2p2h_PPandNN_Pq0"; }
case kGaussianCorr_2p2h_PPandNN_Wq0: { return "GaussianCorr_2p2h_PPandNN_Wq0"; }
case kGaussianCorr_2p2h_PPandNN_Pq3: { return "GaussianCorr_2p2h_PPandNN_Pq3"; }
case kGaussianCorr_2p2h_PPandNN_Wq3: { return "GaussianCorr_2p2h_PPandNN_Wq3"; }
case kGaussianCorr_2p2h_NP_norm: { return "GaussianCorr_2p2h_NP_norm"; }
case kGaussianCorr_2p2h_NP_tilt: { return "GaussianCorr_2p2h_NP_tilt"; }
case kGaussianCorr_2p2h_NP_Pq0: { return "GaussianCorr_2p2h_NP_Pq0"; }
case kGaussianCorr_2p2h_NP_Wq0: { return "GaussianCorr_2p2h_NP_Wq0"; }
case kGaussianCorr_2p2h_NP_Pq3: { return "GaussianCorr_2p2h_NP_Pq3"; }
case kGaussianCorr_2p2h_NP_Wq3: { return "GaussianCorr_2p2h_NP_Wq3"; }
case kGaussianCorr_CC1pi_norm: { return "GaussianCorr_CC1pi_norm"; }
case kGaussianCorr_CC1pi_tilt: { return "GaussianCorr_CC1pi_tilt"; }
case kGaussianCorr_CC1pi_Pq0: { return "GaussianCorr_CC1pi_Pq0"; }
case kGaussianCorr_CC1pi_Wq0: { return "GaussianCorr_CC1pi_Wq0"; }
case kGaussianCorr_CC1pi_Pq3: { return "GaussianCorr_CC1pi_Pq3"; }
case kGaussianCorr_CC1pi_Wq3: { return "GaussianCorr_CC1pi_Wq3"; }
case kGaussianCorr_AllowSuppression: { return "GaussianCorr_AllowSuppression"; }
+ case kBeRPA_A: { return "BeRPA_A"; }
+ case kBeRPA_B: { return "BeRPA_B"; }
+ case kBeRPA_D: { return "BeRPA_D"; }
+ case kBeRPA_E: { return "BeRPA_E"; }
+ case kBeRPA_U: { return "BeRPA_U"; }
+
case kModeNorm_NormRES: { return "NormRES"; }
case kMINERvARW_NormCCQE: { return "MINERvARW_NormCCQE"; }
case kMINERvARW_NormCCMEC: { return "MINERvARW_NormCCMEC"; }
case kMINERvARW_NormCCRES: { return "MINERvARW_NormCCRES"; }
case kMINERvARW_RikRPA_ApplyRPA: { return "MINERvARW_RikRPA_ApplyRPA"; }
case kMINERvARW_RikRPA_LowQ2: { return "MINERvARW_RikRPA_LowQ2"; }
case kMINERvARW_RikRPA_HighQ2: { return "MINERvARW_RikRPA_HighQ2"; }
case kMINERvARW_RikRESRPA_ApplyRPA: { return "MINERvARW_RikRESRPA_ApplyRPA"; }
case kMINERvARW_RikRESRPA_LowQ2: { return "MINERvARW_RikRESRPA_LowQ2"; }
case kMINERvARW_RikRESRPA_HighQ2: { return "MINERvARW_RikRESRPA_HighQ2"; }
case kSBLOsc_Distance: { return "SBLOsc_Distance"; }
case kSBLOsc_MassSplitting: { return "SBLOsc_MassSplitting"; }
case kSBLOsc_Sin2Theta: { return "SBLOsc_Sin2Theta"; }
default: return "unknown_nuisance_dial";
}
};
diff --git a/src/Reweight/NUISANCESyst.h b/src/Reweight/NUISANCESyst.h
index 652e2ff..14825ef 100644
--- a/src/Reweight/NUISANCESyst.h
+++ b/src/Reweight/NUISANCESyst.h
@@ -1,73 +1,78 @@
#ifndef NUISANCESyst_H
#define NUISANCESyst_H
#include "GeneralUtils.h"
namespace Reweight {
enum NUISANCESyst {
kUnkownNUISANCEDial = 0,
kGaussianCorr_CCQE_norm,
kGaussianCorr_CCQE_tilt,
kGaussianCorr_CCQE_Pq0,
kGaussianCorr_CCQE_Wq0,
kGaussianCorr_CCQE_Pq3,
kGaussianCorr_CCQE_Wq3,
kGaussianCorr_2p2h_norm,
kGaussianCorr_2p2h_tilt,
kGaussianCorr_2p2h_Pq0,
kGaussianCorr_2p2h_Wq0,
kGaussianCorr_2p2h_Pq3,
kGaussianCorr_2p2h_Wq3,
kGaussianCorr_2p2h_PPandNN_norm,
kGaussianCorr_2p2h_PPandNN_tilt,
kGaussianCorr_2p2h_PPandNN_Pq0,
kGaussianCorr_2p2h_PPandNN_Wq0,
kGaussianCorr_2p2h_PPandNN_Pq3,
kGaussianCorr_2p2h_PPandNN_Wq3,
kGaussianCorr_2p2h_NP_norm,
kGaussianCorr_2p2h_NP_tilt,
kGaussianCorr_2p2h_NP_Pq0,
kGaussianCorr_2p2h_NP_Wq0,
kGaussianCorr_2p2h_NP_Pq3,
kGaussianCorr_2p2h_NP_Wq3,
kGaussianCorr_CC1pi_norm,
kGaussianCorr_CC1pi_tilt,
kGaussianCorr_CC1pi_Pq0,
kGaussianCorr_CC1pi_Wq0,
kGaussianCorr_CC1pi_Pq3,
kGaussianCorr_CC1pi_Wq3,
kGaussianCorr_AllowSuppression,
+ kBeRPA_A,
+ kBeRPA_B,
+ kBeRPA_D,
+ kBeRPA_E,
+ kBeRPA_U,
+
kModeNorm_NormRES,
kMINERvARW_NormCCQE,
kMINERvARW_NormCCMEC,
kMINERvARW_NormCCRES,
//kMINERvARW_NormCCNonRES1pi,
kMINERvARW_RikRPA_ApplyRPA,
kMINERvARW_RikRPA_LowQ2,
kMINERvARW_RikRPA_HighQ2,
kMINERvARW_RikRESRPA_ApplyRPA,
kMINERvARW_RikRESRPA_LowQ2,
kMINERvARW_RikRESRPA_HighQ2,
kSBLOsc_Distance,
kSBLOsc_MassSplitting,
kSBLOsc_Sin2Theta,
kNUISANCEDial_LAST
};
int ConvertNUISANCEDial(std::string type);
std::string ConvNUISANCEDial(int type);
-
};
#endif
diff --git a/src/Reweight/NUISANCEWeightCalcs.cxx b/src/Reweight/NUISANCEWeightCalcs.cxx
index be3c9ae..8c61eef 100644
--- a/src/Reweight/NUISANCEWeightCalcs.cxx
+++ b/src/Reweight/NUISANCEWeightCalcs.cxx
@@ -1,444 +1,503 @@
#include "NUISANCEWeightCalcs.h"
#include "FitEvent.h"
#include "GeneralUtils.h"
#include "NUISANCESyst.h"
#include "WeightUtils.h"
using namespace Reweight;
ModeNormCalc::ModeNormCalc() { fNormRES = 1.0; }
double ModeNormCalc::CalcWeight(BaseFitEvt* evt) {
int mode = abs(evt->Mode);
double w = 1.0;
if (mode == 11 or mode == 12 or mode == 13) {
w *= fNormRES;
}
return w;
}
void ModeNormCalc::SetDialValue(std::string name, double val) {
SetDialValue(Reweight::ConvDial(name, kCUSTOM), val);
}
void ModeNormCalc::SetDialValue(int rwenum, double val) {
int curenum = rwenum % 1000;
// Check Handled
if (!IsHandled(curenum)) return;
if (curenum == kModeNorm_NormRES) fNormRES = val;
}
bool ModeNormCalc::IsHandled(int rwenum) {
int curenum = rwenum % 1000;
switch (curenum) {
case kModeNorm_NormRES:
return true;
default:
return false;
}
}
+BeRPACalc::BeRPACalc() : fBeRPA_A(0.59), fBeRPA_B(1.05), fBeRPA_D(1.13), fBeRPA_E(0.88), fBeRPA_U(1.2), nParams(0) {
+ // A = 0.59 +/- 20%
+ // B = 1.05 +/- 20%
+ // D = 1.13 +/- 15%
+ // E = 0.88 +/- 40%
+ // U = 1.2
+}
+
+double BeRPACalc::CalcWeight(BaseFitEvt* evt) {
+ FitEvent* fevt = static_cast(evt);
+ int mode = abs(evt->Mode);
+ double w = 1.0;
+ if (nParams == 0) {
+ return w;
+ }
+
+ // Get Q2
+ // Get final state lepton
+ if (mode == 1) {
+ double Q2 = -1.0*(fevt->GetHMFSAnyLeptons()->P4()-fevt->GetNeutrinoIn()->P4())*(fevt->GetHMFSAnyLeptons()->P4()-fevt->GetNeutrinoIn()->P4())/1.E6;
+ // Only CCQE events
+ w *= calcRPA(Q2, fBeRPA_A, fBeRPA_B, fBeRPA_D, fBeRPA_E, fBeRPA_U);
+ }
+
+ return w;
+}
+
+void BeRPACalc::SetDialValue(std::string name, double val) {
+ SetDialValue(Reweight::ConvDial(name, kCUSTOM), val);
+}
+
+void BeRPACalc::SetDialValue(int rwenum, double val) {
+ int curenum = rwenum % 1000;
+
+ // Check Handled
+ if (!IsHandled(curenum)) return;
+ // Need 4 or 5 reconfigures
+ if (curenum == kBeRPA_A) fBeRPA_A = val;
+ else if (curenum == kBeRPA_B) fBeRPA_B = val;
+ else if (curenum == kBeRPA_D) fBeRPA_D = val;
+ else if (curenum == kBeRPA_E) fBeRPA_E = val;
+ else if (curenum == kBeRPA_U) fBeRPA_U = val;
+ nParams++;
+}
+
+bool BeRPACalc::IsHandled(int rwenum) {
+ int curenum = rwenum % 1000;
+ switch (curenum) {
+ case kBeRPA_A:
+ case kBeRPA_B:
+ case kBeRPA_D:
+ case kBeRPA_E:
+ case kBeRPA_U:
+ return true;
+ default:
+ return false;
+ }
+}
+
SBLOscWeightCalc::SBLOscWeightCalc() {
fDistance = 0.0;
fMassSplitting = 0.0;
fSin2Theta = 0.0;
}
double SBLOscWeightCalc::CalcWeight(BaseFitEvt* evt) {
FitEvent* fevt = static_cast(evt);
FitParticle* pnu = fevt->PartInfo(0);
double E = pnu->fP.E() / 1.E3;
// Extract energy
return GetSBLOscWeight(E);
}
void SBLOscWeightCalc::SetDialValue(std::string name, double val) {
SetDialValue(Reweight::ConvDial(name, kCUSTOM), val);
}
void SBLOscWeightCalc::SetDialValue(int rwenum, double val) {
int curenum = rwenum % 1000;
if (!IsHandled(curenum)) return;
if (curenum == kSBLOsc_Distance) fDistance = val;
if (curenum == kSBLOsc_MassSplitting) fMassSplitting = val;
if (curenum == kSBLOsc_Sin2Theta) fSin2Theta = val;
}
bool SBLOscWeightCalc::IsHandled(int rwenum) {
int curenum = rwenum % 1000;
switch (curenum) {
case kSBLOsc_Distance:
return true;
case kSBLOsc_MassSplitting:
return true;
case kSBLOsc_Sin2Theta:
return true;
default:
return false;
}
}
double SBLOscWeightCalc::GetSBLOscWeight(double E) {
if (E <= 0.0) return 1.0 - 0.5 * fSin2Theta;
return 1.0 - fSin2Theta * pow(sin(1.267 * fMassSplitting * fDistance / E), 2);
}
GaussianModeCorr::GaussianModeCorr() {
// Apply the tilt-shift Gauss by Patrick
// Alternatively set in config
fMethod = true;
// Init
fApply_CCQE = false;
fGausVal_CCQE[kPosNorm] = 0.0;
fGausVal_CCQE[kPosTilt] = 0.0;
fGausVal_CCQE[kPosPq0] = 1.0;
fGausVal_CCQE[kPosWq0] = 1.0;
fGausVal_CCQE[kPosPq3] = 1.0;
fGausVal_CCQE[kPosWq3] = 1.0;
fApply_2p2h = false;
fGausVal_2p2h[kPosNorm] = 0.0;
fGausVal_2p2h[kPosTilt] = 0.0;
fGausVal_2p2h[kPosPq0] = 1.0;
fGausVal_2p2h[kPosWq0] = 1.0;
fGausVal_2p2h[kPosPq3] = 1.0;
fGausVal_2p2h[kPosWq3] = 1.0;
fApply_2p2h_PPandNN = false;
fGausVal_2p2h_PPandNN[kPosNorm] = 0.0;
fGausVal_2p2h_PPandNN[kPosTilt] = 0.0;
fGausVal_2p2h_PPandNN[kPosPq0] = 1.0;
fGausVal_2p2h_PPandNN[kPosWq0] = 1.0;
fGausVal_2p2h_PPandNN[kPosPq3] = 1.0;
fGausVal_2p2h_PPandNN[kPosWq3] = 1.0;
fApply_2p2h_NP = false;
fGausVal_2p2h_NP[kPosNorm] = 0.0;
fGausVal_2p2h_NP[kPosTilt] = 0.0;
fGausVal_2p2h_NP[kPosPq0] = 1.0;
fGausVal_2p2h_NP[kPosWq0] = 1.0;
fGausVal_2p2h_NP[kPosPq3] = 1.0;
fGausVal_2p2h_NP[kPosWq3] = 1.0;
fApply_CC1pi = false;
fGausVal_CC1pi[kPosNorm] = 0.0;
fGausVal_CC1pi[kPosTilt] = 0.0;
fGausVal_CC1pi[kPosPq0] = 1.0;
fGausVal_CC1pi[kPosWq0] = 1.0;
fGausVal_CC1pi[kPosPq3] = 1.0;
fGausVal_CC1pi[kPosWq3] = 1.0;
fAllowSuppression = false;
fDebugStatements = FitPar::Config().GetParB("GaussianModeCorr_DEBUG");
//fDebugStatements = true;
}
double GaussianModeCorr::CalcWeight(BaseFitEvt* evt) {
FitEvent* fevt = static_cast(evt);
double rw_weight = 1.0;
// Get Neutrino
if (!fevt->Npart()) {
THROW("NO particles found in stack!");
}
FitParticle* pnu = fevt->GetHMISAnyLeptons();
if (!pnu) {
THROW("NO Starting particle found in stack!");
}
int pdgnu = pnu->fPID;
int expect_fsleppdg = 0;
if (pdgnu & 1) {
expect_fsleppdg = pdgnu;
} else {
expect_fsleppdg = abs(pdgnu) - 1;
}
FitParticle* plep = fevt->GetHMFSParticle(expect_fsleppdg);
if (!plep) return 1.0;
TLorentzVector q = pnu->fP - plep->fP;
// Extra q0,q3
double q0 = fabs(q.E()) / 1.E3;
double q3 = fabs(q.Vect().Mag()) / 1.E3;
int initialstate = -1; // Undef
if (abs(fevt->Mode) == 2) {
int npr = 0;
int nne = 0;
for (UInt_t j = 0; j < fevt->Npart(); j++) {
if ((fevt->PartInfo(j))->fIsAlive) continue;
if (fevt->PartInfo(j)->fPID == 2212)
npr++;
else if (fevt->PartInfo(j)->fPID == 2112)
nne++;
}
if (fevt->Mode == 2 && npr == 1 && nne == 1) {
initialstate = 2;
} else if (fevt->Mode == 2 && ((npr == 0 && nne == 2) || (npr == 2 && nne == 0))) {
initialstate = 1;
}
}
// Apply weighting
if (fApply_CCQE && abs(fevt->Mode) == 1) {
if (fDebugStatements) std::cout << "Getting CCQE Weight" << std::endl;
double g = GetGausWeight(q0, q3, fGausVal_CCQE);
if (g < 1.0) g = 1.0;
rw_weight *= g;
}
if (fApply_2p2h && abs(fevt->Mode) == 2) {
if (fDebugStatements) std::cout << "Getting 2p2h Weight" << std::endl;
if (fDebugStatements) std::cout << "Got q0 q3 = " << q0 << " " << q3 << " mode = " << fevt->Mode << std::endl;
rw_weight *= GetGausWeight(q0, q3, fGausVal_2p2h);
if (fDebugStatements) std::cout << "Returning Weight " << rw_weight << std::endl;
}
if (fApply_2p2h_PPandNN && abs(fevt->Mode) == 2 && initialstate == 1) {
if (fDebugStatements) std::cout << "Getting 2p2h PPandNN Weight" << std::endl;
rw_weight *= GetGausWeight(q0, q3, fGausVal_2p2h_PPandNN);
}
if (fApply_2p2h_NP && abs(fevt->Mode) == 2 && initialstate == 2) {
if (fDebugStatements) std::cout << "Getting 2p2h NP Weight" << std::endl;
rw_weight *= GetGausWeight(q0, q3, fGausVal_2p2h_NP);
}
if (fApply_CC1pi && abs(fevt->Mode) >= 11 && abs(fevt->Mode) <= 13) {
if (fDebugStatements) std::cout << "Getting CC1pi Weight" << std::endl;
rw_weight *= GetGausWeight(q0, q3, fGausVal_CC1pi);
}
return rw_weight;
}
void GaussianModeCorr::SetMethod(bool method) {
fMethod = method;
if (fMethod == true) {
LOG(FIT) << " Using tilt-shift Gaussian parameters for Gaussian enhancement..." << std::endl;
} else {
LOG(FIT) << " Using Normal Gaussian parameters for Gaussian enhancement..." << std::endl;
}
};
double GaussianModeCorr::GetGausWeight(double q0, double q3, double vals[]) {
// The weight
double w = 1.0;
// Use tilt-shift method by Patrick
if (fMethod) {
if (fDebugStatements) {
std::cout << "Using Patrick gaussian" << std::endl;
}
// // CCQE Without Suppression
// double Norm = 4.82788679036;
// double Tilt = 2.3501416116;
// double Pq0 = 0.363964889702;
// double Wq0 = 0.133976806938;
// double Pq3 = 0.431769740224;
// double Wq3 = 0.207666663434;
// // Also add for CCQE at the end
// return (w > 1.0) ? w : 1.0;
// // 2p2h with suppression
// double Norm = 15.967;
// double Tilt = -0.455655;
// double Pq0 = 0.214598;
// double Wq0 = 0.0291061;
// double Pq3 = 0.480194;
// double Wq3 = 0.134588;
double Norm = vals[kPosNorm];
double Tilt = vals[kPosTilt];
double Pq0 = vals[kPosPq0];
double Wq0 = vals[kPosWq0];
double Pq3 = vals[kPosPq3];
double Wq3 = vals[kPosWq3];
double a = cos(Tilt) * cos(Tilt) / (2 * Wq0 * Wq0);
a += sin(Tilt) * sin(Tilt) / (2 * Wq3 * Wq3);
double b = -sin(2 * Tilt) / (4 * Wq0 * Wq0);
b += sin(2 * Tilt) / (4 * Wq3 * Wq3);
double c = sin(Tilt) * sin(Tilt) / (2 * Wq0 * Wq0);
c += cos(Tilt) * cos(Tilt) / (2 * Wq3 * Wq3);
w = Norm;
w *= exp(-a * (q0 - Pq0) * (q0 - Pq0));
w *= exp(+2.0 * b * (q0 - Pq0) * (q3 - Pq3));
w *= exp(-c * (q3 - Pq3) * (q3 - Pq3));
if (fDebugStatements) {
std::cout << "Applied Tilt " << Tilt << " " << cos(Tilt) << " " << sin(Tilt) << std::endl;
std::cout << "abc = " << a << " " << b << " " << c << std::endl;
std::cout << "Returning " << Norm << " " << Pq0 << " " << Wq0 << " " << Pq3 << " " << Wq3 << " " << w << std::endl;
}
if (w != w || std::isnan(w) || w < 0.0) {
w = 0.0;
}
if (w < 1.0 and !fAllowSuppression) {
w = 1.0;
}
// Use the MINERvA Gaussian method
} else {
/*
* From MINERvA and Daniel Ruterbories:
* Old notes here: * http://cdcvs.fnal.gov/cgi-bin/public-cvs/cvsweb-public.cgi/AnalysisFramework/Ana/CCQENu/ana_common/data/?cvsroot=mnvsoft
* These parameters are slightly altered
*
* FRESH:
* 10.5798
* 0.254032
* 0.50834
* 0.0571035
* 0.129051
* 0.875287
*/
if (fDebugStatements) {
std::cout << "Using MINERvA Gaussian" << std::endl;
}
double norm = vals[kPosNorm];
double meanq0 = vals[kPosTilt];
double meanq3 = vals[kPosPq0];
double sigmaq0 = vals[kPosWq0];
double sigmaq3 = vals[kPosPq3];
double corr = vals[kPosWq3];
double z = (q0 - meanq0)*(q0 - meanq0) /sigmaq0/sigmaq0
+ (q3 - meanq3)*(q3 - meanq3) / sigmaq3/sigmaq3
- 2*corr*(q0-meanq0)*(q3-meanq3)/ (sigmaq0 * sigmaq3);
double ret = norm*exp( -0.5 * z / (1 - corr*corr) );
//Need to add 1 to the results
w = 1.0 + ret;
}
return w;
}
void GaussianModeCorr::SetDialValue(std::string name, double val) {
SetDialValue(Reweight::ConvDial(name, kCUSTOM), val);
}
void GaussianModeCorr::SetDialValue(int rwenum, double val) {
int curenum = rwenum % 1000;
// Check Handled
if (!IsHandled(curenum)) return;
// CCQE Setting
for (int i = kGaussianCorr_CCQE_norm; i <= kGaussianCorr_CCQE_Wq3; i++) {
if (i == curenum) {
int index = i - kGaussianCorr_CCQE_norm;
fGausVal_CCQE[index] = val;
fApply_CCQE = true;
}
}
// 2p2h Setting
for (int i = kGaussianCorr_2p2h_norm; i <= kGaussianCorr_2p2h_Wq3; i++) {
if (i == curenum) {
int index = i - kGaussianCorr_2p2h_norm;
fGausVal_2p2h[index] = val;
fApply_2p2h = true;
}
}
// 2p2h_PPandNN Setting
for (int i = kGaussianCorr_2p2h_PPandNN_norm; i <= kGaussianCorr_2p2h_PPandNN_Wq3; i++) {
if (i == curenum) {
int index = i - kGaussianCorr_2p2h_PPandNN_norm;
fGausVal_2p2h_PPandNN[index] = val;
fApply_2p2h_PPandNN = true;
}
}
// 2p2h_NP Setting
for (int i = kGaussianCorr_2p2h_NP_norm; i <= kGaussianCorr_2p2h_NP_Wq3; i++) {
if (i == curenum) {
int index = i - kGaussianCorr_2p2h_NP_norm;
fGausVal_2p2h_NP[index] = val;
fApply_2p2h_NP = true;
}
}
// CC1pi Setting
for (int i = kGaussianCorr_CC1pi_norm; i <= kGaussianCorr_CC1pi_Wq3; i++) {
if (i == curenum) {
int index = i - kGaussianCorr_CC1pi_norm;
fGausVal_CC1pi[index] = val;
fApply_CC1pi = true;
}
}
if (curenum == kGaussianCorr_AllowSuppression) {
fAllowSuppression = (val > 0.5);
}
}
bool GaussianModeCorr::IsHandled(int rwenum) {
int curenum = rwenum % 1000;
switch (curenum) {
case kGaussianCorr_CCQE_norm:
case kGaussianCorr_CCQE_tilt:
case kGaussianCorr_CCQE_Pq0:
case kGaussianCorr_CCQE_Wq0:
case kGaussianCorr_CCQE_Pq3:
case kGaussianCorr_CCQE_Wq3:
case kGaussianCorr_2p2h_norm:
case kGaussianCorr_2p2h_tilt:
case kGaussianCorr_2p2h_Pq0:
case kGaussianCorr_2p2h_Wq0:
case kGaussianCorr_2p2h_Pq3:
case kGaussianCorr_2p2h_Wq3:
case kGaussianCorr_2p2h_PPandNN_norm:
case kGaussianCorr_2p2h_PPandNN_tilt:
case kGaussianCorr_2p2h_PPandNN_Pq0:
case kGaussianCorr_2p2h_PPandNN_Wq0:
case kGaussianCorr_2p2h_PPandNN_Pq3:
case kGaussianCorr_2p2h_PPandNN_Wq3:
case kGaussianCorr_2p2h_NP_norm:
case kGaussianCorr_2p2h_NP_tilt:
case kGaussianCorr_2p2h_NP_Pq0:
case kGaussianCorr_2p2h_NP_Wq0:
case kGaussianCorr_2p2h_NP_Pq3:
case kGaussianCorr_2p2h_NP_Wq3:
case kGaussianCorr_CC1pi_norm:
case kGaussianCorr_CC1pi_tilt:
case kGaussianCorr_CC1pi_Pq0:
case kGaussianCorr_CC1pi_Wq0:
case kGaussianCorr_CC1pi_Pq3:
case kGaussianCorr_CC1pi_Wq3:
case kGaussianCorr_AllowSuppression:
return true;
default:
return false;
}
}
diff --git a/src/Reweight/NUISANCEWeightCalcs.h b/src/Reweight/NUISANCEWeightCalcs.h
index c09b424..a1d2b9e 100644
--- a/src/Reweight/NUISANCEWeightCalcs.h
+++ b/src/Reweight/NUISANCEWeightCalcs.h
@@ -1,106 +1,129 @@
#ifndef NUISANCE_WEIGHT_CALCS
#define NUISANCE_WEIGHT_CALCS
#include "BaseFitEvt.h"
+#include "BeRPA.h"
class NUISANCEWeightCalc {
public:
NUISANCEWeightCalc() {};
virtual ~NUISANCEWeightCalc() {};
virtual double CalcWeight(BaseFitEvt* evt){return 1.0;};
virtual void SetDialValue(std::string name, double val){};
virtual void SetDialValue(int rwenum, double val){};
virtual bool IsHandled(int rwenum){return false;};
virtual void Print(){};
std::map fDialNameIndex;
std::map fDialEnumIndex;
std::vector fDialValues;
std::string fName;
};
class ModeNormCalc : public NUISANCEWeightCalc {
public:
ModeNormCalc();
~ModeNormCalc(){};
double CalcWeight(BaseFitEvt* evt);
void SetDialValue(std::string name, double val);
void SetDialValue(int rwenum, double val);
bool IsHandled(int rwenum);
double fNormRES;
};
+class BeRPACalc : public NUISANCEWeightCalc {
+ public:
+ BeRPACalc();
+ ~BeRPACalc(){};
+
+ double CalcWeight(BaseFitEvt* evt);
+ void SetDialValue(std::string name, double val);
+ void SetDialValue(int rwenum, double val);
+ bool IsHandled(int rwenum);
+
+ private:
+ // Parameter values
+ double fBeRPA_A;
+ double fBeRPA_B;
+ double fBeRPA_D;
+ double fBeRPA_E;
+ double fBeRPA_U;
+ // Counts of enabled parameters
+ int nParams;
+
+};
+
class SBLOscWeightCalc : public NUISANCEWeightCalc {
public:
SBLOscWeightCalc();
~SBLOscWeightCalc(){};
double CalcWeight(BaseFitEvt* evt);
void SetDialValue(std::string name, double val);
void SetDialValue(int rwenum, double val);
bool IsHandled(int rwenum);
double GetSBLOscWeight(double E);
double fDistance;
double fMassSplitting;
double fSin2Theta;
};
class GaussianModeCorr : public NUISANCEWeightCalc {
public:
GaussianModeCorr();
~GaussianModeCorr(){};
double CalcWeight(BaseFitEvt* evt);
void SetDialValue(std::string name, double val);
void SetDialValue(int rwenum, double val);
bool IsHandled(int rwenum);
double GetGausWeight(double q0, double q3, double vals[]);
// Set the Gaussian method (tilt-shift or normal Gaussian parameters)
void SetMethod(bool method);
// 5 pars describe the Gaussain
// 0 norm.
// 1 q0 pos
// 2 q0 width
// 3 q3 pos
// 4 q3 width
static const int kPosNorm = 0;
static const int kPosTilt = 1;
static const int kPosPq0 = 2;
static const int kPosWq0 = 3;
static const int kPosPq3 = 4;
static const int kPosWq3 = 5;
bool fApply_CCQE;
double fGausVal_CCQE[6];
bool fApply_2p2h;
double fGausVal_2p2h[6];
bool fApply_2p2h_PPandNN;
double fGausVal_2p2h_PPandNN[6];
bool fApply_2p2h_NP;
double fGausVal_2p2h_NP[6];
bool fApply_CC1pi;
double fGausVal_CC1pi[6];
bool fAllowSuppression;
bool fDebugStatements;
bool fMethod;
};
#endif
diff --git a/src/Reweight/NUISANCEWeightEngine.cxx b/src/Reweight/NUISANCEWeightEngine.cxx
index ebb62fc..39be11c 100644
--- a/src/Reweight/NUISANCEWeightEngine.cxx
+++ b/src/Reweight/NUISANCEWeightEngine.cxx
@@ -1,118 +1,119 @@
#include "NUISANCEWeightEngine.h"
#include "NUISANCEWeightCalcs.h"
#ifdef __MINERVA_RW_ENABLED__
#ifdef __GENIE_ENABLED__
#include "MINERvAWeightCalcs.h"
#endif
#endif
NUISANCEWeightEngine::NUISANCEWeightEngine(std::string name) {
// Setup the NUISANCE Reweight engine
fCalcName = name;
LOG(FIT) << "Setting up NUISANCE Custom RW : " << fCalcName << std::endl;
// Load in all Weight Calculations
GaussianModeCorr* GaussianMode = new GaussianModeCorr();
std::string Gaussian_Method = FitPar::Config().GetParS("Gaussian_Enhancement");
if (Gaussian_Method == "Tilt-Shift") {
GaussianMode->SetMethod(true);
} else if (Gaussian_Method == "Normal") {
GaussianMode->SetMethod(false);
} else {
ERR(FTL) << "I do not recognise method " << Gaussian_Method << " for the Gaussian enhancement, so will die now..." << std::endl;
throw;
}
fWeightCalculators.push_back(GaussianMode);
fWeightCalculators.push_back(new ModeNormCalc());
fWeightCalculators.push_back(new SBLOscWeightCalc());
+ fWeightCalculators.push_back(new BeRPACalc());
#ifdef __MINERVA_RW_ENABLED__
#ifdef __GENIE_ENABLED__
fWeightCalculators.push_back(new nuisance::reweight::MINERvAReWeight_MEC());
fWeightCalculators.push_back(new nuisance::reweight::MINERvAReWeight_RES());
fWeightCalculators.push_back(new nuisance::reweight::RikRPA());
#endif
#endif
// Set Abs Twk Config
fIsAbsTwk = true;
};
void NUISANCEWeightEngine::IncludeDial(std::string name, double startval) {
// Get First enum
int nuisenum = Reweight::ConvDial(name, kCUSTOM);
// Setup Maps
fEnumIndex[nuisenum]; // = std::vector(0);
fNameIndex[name]; // = std::vector(0);
// Split by commas
std::vector allnames = GeneralUtils::ParseToStr(name, ",");
for (uint i = 0; i < allnames.size(); i++) {
std::string singlename = allnames[i];
// Get RW
int singleenum = Reweight::ConvDial(singlename, kCUSTOM);
// Fill Maps
int index = fValues.size();
fValues.push_back(0.0);
fNUISANCEEnums.push_back(singleenum);
// Initialize dial
LOG(FIT) << "Registering " << singlename << " from " << name << std::endl;
// Setup index
fEnumIndex[nuisenum].push_back(index);
fNameIndex[name].push_back(index);
}
// Set Value if given
if (startval != -999.9) {
SetDialValue(nuisenum, startval);
}
};
void NUISANCEWeightEngine::SetDialValue(int nuisenum, double val) {
std::vector indices = fEnumIndex[nuisenum];
for (uint i = 0; i < indices.size(); i++) {
fValues[indices[i]] = val;
}
}
void NUISANCEWeightEngine::SetDialValue(std::string name, double val) {
std::vector indices = fNameIndex[name];
for (uint i = 0; i < indices.size(); i++) {
fValues[indices[i]] = val;
}
}
void NUISANCEWeightEngine::Reconfigure(bool silent) {
for (size_t i = 0; i < fNUISANCEEnums.size(); i++) {
for (std::vector::iterator calciter = fWeightCalculators.begin(); calciter != fWeightCalculators.end(); calciter++) {
NUISANCEWeightCalc* nuiscalc = static_cast(*calciter);
if (nuiscalc->IsHandled(fNUISANCEEnums[i])) {
nuiscalc->SetDialValue(fNUISANCEEnums[i], fValues[i]);
}
}
}
}
double NUISANCEWeightEngine::CalcWeight(BaseFitEvt* evt) {
double rw_weight = 1.0;
// Cast as usable class
for (std::vector::iterator iter =
fWeightCalculators.begin();
iter != fWeightCalculators.end(); iter++) {
NUISANCEWeightCalc* nuiscalc = static_cast(*iter);
rw_weight *= nuiscalc->CalcWeight(evt);
}
// Return rw_weight
return rw_weight;
}
diff --git a/src/Reweight/T2KWeightEngine.cxx b/src/Reweight/T2KWeightEngine.cxx
index ac0123b..6fadd76 100644
--- a/src/Reweight/T2KWeightEngine.cxx
+++ b/src/Reweight/T2KWeightEngine.cxx
@@ -1,149 +1,140 @@
#include "T2KWeightEngine.h"
T2KWeightEngine::T2KWeightEngine(std::string name) {
#ifdef __T2KREW_ENABLED__
- // Setup the NEUT Reweight engien
- fCalcName = name;
- LOG(FIT) << "Setting up T2K RW : " << fCalcName << std::endl;
+ // Setup the NEUT Reweight engien
+ fCalcName = name;
+ LOG(FIT) << "Setting up T2K RW : " << fCalcName << std::endl;
- // Create RW Engine suppressing cout
- StopTalking();
+ // Create RW Engine suppressing cout
+ StopTalking();
- // Create Main RW Engine
- fT2KRW = new t2krew::T2KReWeight();
+ // Create Main RW Engine
+ fT2KRW = new t2krew::T2KReWeight();
- // Setup Sub RW Engines (Only activated for neut and niwg)
- fT2KNeutRW = new t2krew::T2KNeutReWeight();
- fT2KNIWGRW = new t2krew::T2KNIWGReWeight();
+ // Setup Sub RW Engines (Only activated for neut and niwg)
+ fT2KNeutRW = new t2krew::T2KNeutReWeight();
+ fT2KNIWGRW = new t2krew::T2KNIWGReWeight();
- fT2KRW->AdoptWghtEngine("fNeutRW", fT2KNeutRW);
- fT2KRW->AdoptWghtEngine("fNIWGRW", fT2KNIWGRW);
+ fT2KRW->AdoptWghtEngine("fNeutRW", fT2KNeutRW);
+ fT2KRW->AdoptWghtEngine("fNIWGRW", fT2KNIWGRW);
- fT2KRW->Reconfigure();
+ fT2KRW->Reconfigure();
- // allow cout again
- StartTalking();
+ // allow cout again
+ StartTalking();
- // Set Abs Twk Config
- fIsAbsTwk = (FitPar::Config().GetParB("setabstwk"));
+ // Set Abs Twk Config
+ fIsAbsTwk = (FitPar::Config().GetParB("setabstwk"));
#else
- ERR(FTL) << "T2K RW NOT ENABLED" << std::endl;
- throw;
+ ERR(FTL) << "T2K RW NOT ENABLED" << std::endl;
+ throw;
#endif
};
void T2KWeightEngine::IncludeDial(std::string name, double startval) {
#ifdef __T2KREW_ENABLED__
- // Get First enum
- int nuisenum = Reweight::ConvDial(name, kT2K);
+ // Get First enum
+ int nuisenum = Reweight::ConvDial(name, kT2K);
- // Setup Maps
- fEnumIndex[nuisenum];// = std::vector(0);
- fNameIndex[name]; // = std::vector(0);
+ // Setup Maps
+ fEnumIndex[nuisenum];// = std::vector(0);
+ fNameIndex[name]; // = std::vector(0);
- // Split by commas
- std::vector allnames = GeneralUtils::ParseToStr(name, ",");
- for (uint i = 0; i < allnames.size(); i++) {
- std::string singlename = allnames[i];
+ // Split by commas
+ std::vector allnames = GeneralUtils::ParseToStr(name, ",");
+ for (uint i = 0; i < allnames.size(); i++) {
+ std::string singlename = allnames[i];
- // Get RW
- t2krew::T2KSyst_t gensyst = t2krew::T2KSyst::FromString(name);
+ // Get RW
+ t2krew::T2KSyst_t gensyst = t2krew::T2KSyst::FromString(name);
- // Fill Maps
- int index = fValues.size();
- fValues.push_back(0.0);
- fT2KSysts.push_back(gensyst);
+ // Fill Maps
+ int index = fValues.size();
+ fValues.push_back(0.0);
+ fT2KSysts.push_back(gensyst);
- // Initialize dial
- std::cout << "Registering " << singlename << " from " << name << std::endl;
- fT2KRW->Systematics().Include(gensyst);
+ // Initialize dial
+ std::cout << "Registering " << singlename << " from " << name << std::endl;
+ fT2KRW->Systematics().Include(gensyst);
- // If Absolute
- if (fIsAbsTwk) {
- fT2KRW->Systematics().SetAbsTwk(gensyst);
- }
+ // If Absolute
+ if (fIsAbsTwk) {
+ fT2KRW->Systematics().SetAbsTwk(gensyst);
+ }
- // Setup index
- fEnumIndex[nuisenum].push_back(index);
- fNameIndex[name].push_back(index);
+ // Setup index
+ fEnumIndex[nuisenum].push_back(index);
+ fNameIndex[name].push_back(index);
- }
+ }
- // Set Value if given
- if (startval != -999.9) {
- SetDialValue(nuisenum, startval);
- }
+ // Set Value if given
+ if (startval != -999.9) {
+ SetDialValue(nuisenum, startval);
+ }
#endif
}
void T2KWeightEngine::SetDialValue(int nuisenum, double val) {
#ifdef __T2KREW_ENABLED__
- std::vector indices = fEnumIndex[nuisenum];
- for (uint i = 0; i < indices.size(); i++) {
- fValues[indices[i]] = val;
- fT2KRW->Systematics().SetTwkDial(fT2KSysts[indices[i]], val);
- }
+ std::vector indices = fEnumIndex[nuisenum];
+ for (uint i = 0; i < indices.size(); i++) {
+ fValues[indices[i]] = val;
+ fT2KRW->Systematics().SetTwkDial(fT2KSysts[indices[i]], val);
+ }
#endif
}
void T2KWeightEngine::SetDialValue(std::string name, double val) {
#ifdef __T2KREW_ENABLED__
- std::vector indices = fNameIndex[name];
- for (uint i = 0; i < indices.size(); i++) {
- fValues[indices[i]] = val;
- fT2KRW->Systematics().SetTwkDial(fT2KSysts[indices[i]], val);
- }
+ std::vector indices = fNameIndex[name];
+ for (uint i = 0; i < indices.size(); i++) {
+ fValues[indices[i]] = val;
+ fT2KRW->Systematics().SetTwkDial(fT2KSysts[indices[i]], val);
+ }
#endif
}
void T2KWeightEngine::Reconfigure(bool silent) {
#ifdef __T2KREW_ENABLED__
- // Hush now...
- if (silent) StopTalking();
+ // Hush now...
+ if (silent) StopTalking();
- // Reconf
- StopTalking();
- fT2KRW->Reconfigure();
- StartTalking();
+ // Reconf
+ StopTalking();
+ fT2KRW->Reconfigure();
+ StartTalking();
- // Shout again
- if (silent) StartTalking();
+ // Shout again
+ if (silent) StartTalking();
#endif
}
double T2KWeightEngine::CalcWeight(BaseFitEvt* evt) {
- double rw_weight = 1.0;
+ double rw_weight = 1.0;
#ifdef __T2KREW_ENABLED__
- // Skip Non GENIE
- if (evt->fType != kNEUT) return 1.0;
+ // Skip Non GENIE
+ if (evt->fType != kNEUT) return 1.0;
- // Hush now
- StopTalking();
+ // Hush now
+ StopTalking();
- // Get Weight For NEUT
- rw_weight = fT2KRW->CalcWeight(evt->fNeutVect);
+ // Get Weight For NEUT
+ rw_weight = fT2KRW->CalcWeight(evt->fNeutVect);
- // Speak Now
- StartTalking();
+ // Speak Now
+ StartTalking();
#endif
- // Return rw_weight
- return rw_weight;
+ // Return rw_weight
+ return rw_weight;
}
-
-
-
-
-
-
-
-
-
diff --git a/src/Utils/FitUtils.cxx b/src/Utils/FitUtils.cxx
index 225c98c..87b1efe 100644
--- a/src/Utils/FitUtils.cxx
+++ b/src/Utils/FitUtils.cxx
@@ -1,1095 +1,1186 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#include "FitUtils.h"
/*
MISC Functions
*/
//********************************************************************
double *FitUtils::GetArrayFromMap(std::vector invals,
std::map inmap) {
//********************************************************************
double *outarr = new double[invals.size()];
int count = 0;
for (size_t i = 0; i < invals.size(); i++) {
outarr[count++] = inmap[invals.at(i)];
}
return outarr;
}
/*
MISC Event
*/
//********************************************************************
// Returns the kinetic energy of a particle in GeV
double FitUtils::T(TLorentzVector part) {
//********************************************************************
double E_part = part.E() / 1000.;
double p_part = part.Vect().Mag() / 1000.;
double m_part = sqrt(E_part * E_part - p_part * p_part);
double KE_part = E_part - m_part;
return KE_part;
};
//********************************************************************
// Returns the momentum of a particle in GeV
double FitUtils::p(TLorentzVector part) {
//********************************************************************
double p_part = part.Vect().Mag() / 1000.;
return p_part;
};
double FitUtils::p(FitParticle *part) { return FitUtils::p(part->fP); };
//********************************************************************
// Returns the angle between two particles in radians
double FitUtils::th(TLorentzVector part1, TLorentzVector part2) {
//********************************************************************
double th = part1.Vect().Angle(part2.Vect());
return th;
};
double FitUtils::th(FitParticle *part1, FitParticle *part2) {
return FitUtils::th(part1->fP, part2->fP);
};
// T2K CC1pi+ helper functions
//
//********************************************************************
// Returns the angle between q3 and the pion defined in Raquel's CC1pi+ on CH
// paper
// Uses "MiniBooNE formula" for Enu, here called EnuCC1pip_T2K_MB
//********************************************************************
double FitUtils::thq3pi_CC1pip_T2K(TLorentzVector pnu, TLorentzVector pmu,
TLorentzVector ppi) {
// Want this in GeV
TVector3 p_mu = pmu.Vect() * (1. / 1000.);
// Get the reconstructed Enu
// We are not using Michel e sample, so we have pion kinematic information
double Enu = EnuCC1piprec(pnu, pmu, ppi, true);
// Get neutrino unit direction, multiply by reconstructed Enu
TVector3 p_nu = pnu.Vect() * (1. / (pnu.Vect().Mag())) * Enu;
TVector3 p_pi = ppi.Vect() * (1. / 1000.);
// This is now in GeV
TVector3 q3 = (p_nu - p_mu);
// Want this in GeV
double th_q3_pi = q3.Angle(p_pi);
return th_q3_pi;
}
//********************************************************************
// Returns the q3 defined in Raquel's CC1pi+ on CH paper
// Uses "MiniBooNE formula" for Enu
//********************************************************************
double FitUtils::q3_CC1pip_T2K(TLorentzVector pnu, TLorentzVector pmu,
TLorentzVector ppi) {
// Can't use the true Enu here; need to reconstruct it
// We do have Michel e- here so reconstruct Enu by "MiniBooNE formula" without
// pion kinematics
// The bool false refers to that we don't have pion kinematics
// Last bool refers to if we have pion kinematic information or not
double Enu = EnuCC1piprec(pnu, pmu, ppi, false);
TVector3 p_mu = pmu.Vect() * (1. / 1000.);
TVector3 p_nu = pnu.Vect() * (1. / pnu.Vect().Mag()) * Enu;
double q3 = (p_nu - p_mu).Mag();
return q3;
}
//********************************************************************
// Returns the W reconstruction from Raquel CC1pi+ CH thesis
// Uses the MiniBooNE formula Enu
//********************************************************************
double FitUtils::WrecCC1pip_T2K_MB(TLorentzVector pnu, TLorentzVector pmu,
TLorentzVector ppi) {
double E_mu = pmu.E() / 1000.;
double p_mu = pmu.Vect().Mag() / 1000.;
double E_nu = EnuCC1piprec(pnu, pmu, ppi, false);
double a1 = (E_nu + PhysConst::mass_neutron) - E_mu;
double a2 = E_nu - p_mu;
double wrec = sqrt(a1 * a1 - a2 * a2);
return wrec;
}
//********************************************************
double FitUtils::ProtonQ2QErec(double pE, double binding) {
//********************************************************
const double V = binding / 1000.; // binding potential
const double mn = PhysConst::mass_neutron; // neutron mass
const double mp = PhysConst::mass_proton; // proton mass
const double mn_eff = mn - V; // effective proton mass
const double pki = (pE / 1000.0) - mp;
double q2qe = mn_eff * mn_eff - mp * mp + 2 * mn_eff * (pki + mp - mn_eff);
return q2qe;
};
//********************************************************************
double FitUtils::EnuQErec(TLorentzVector pmu, double costh, double binding,
bool neutrino) {
//********************************************************************
// Convert all values to GeV
const double V = binding / 1000.; // binding potential
const double mn = PhysConst::mass_neutron; // neutron mass
const double mp = PhysConst::mass_proton; // proton mass
double mN_eff = mn - V;
double mN_oth = mp;
if (!neutrino) {
mN_eff = mp - V;
mN_oth = mn;
}
double el = pmu.E() / 1000.;
double pl = (pmu.Vect().Mag()) / 1000.; // momentum of lepton
double ml = sqrt(el * el - pl * pl); // lepton mass
double rEnu =
(2 * mN_eff * el - ml * ml + mN_oth * mN_oth - mN_eff * mN_eff) /
(2 * (mN_eff - el + pl * costh));
return rEnu;
};
double FitUtils::Q2QErec(TLorentzVector pmu, double costh, double binding, bool neutrino) {
double el = pmu.E() / 1000.;
double pl = (pmu.Vect().Mag()) / 1000.; // momentum of lepton
double ml = sqrt(el * el - pl * pl); // lepton mass
double rEnu = EnuQErec(pmu, costh, binding, neutrino);
double q2 = -ml * ml + 2. * rEnu * (el - pl * costh);
return q2;
};
double FitUtils::Q2QErec(TLorentzVector Pmu, TLorentzVector Pnu, double binding, bool neutrino) {
double q2qe = Q2QErec(Pmu, cos(Pnu.Vect().Angle(Pmu.Vect())), binding, neutrino);
return q2qe;
}
double FitUtils::EnuQErec(double pl, double costh, double binding,
bool neutrino) {
if (pl < 0) return 0.; // Make sure nobody is silly
double mN_eff = PhysConst::mass_neutron - binding / 1000.;
double mN_oth = PhysConst::mass_proton;
if (!neutrino) {
mN_eff = PhysConst::mass_proton - binding / 1000.;
mN_oth = PhysConst::mass_neutron;
}
double ml = PhysConst::mass_muon;
double el = sqrt(pl * pl + ml * ml);
double rEnu =
(2 * mN_eff * el - ml * ml + mN_oth * mN_oth - mN_eff * mN_eff) /
(2 * (mN_eff - el + pl * costh));
return rEnu;
};
double FitUtils::Q2QErec(double pl, double costh, double binding,
bool neutrino) {
if (pl < 0) return 0.; // Make sure nobody is silly
double ml = PhysConst::mass_muon;
double el = sqrt(pl * pl + ml * ml);
double rEnu = EnuQErec(pl, costh, binding, neutrino);
double q2 = -ml * ml + 2. * rEnu * (el - pl * costh);
return q2;
};
//********************************************************************
// Reconstructs Enu for CC1pi0
// Very similar for CC1pi+ reconstruction
double FitUtils::EnuCC1pi0rec(TLorentzVector pnu, TLorentzVector pmu,
TLorentzVector ppi0) {
//********************************************************************
double E_mu = pmu.E() / 1000;
double p_mu = pmu.Vect().Mag() / 1000;
double m_mu = sqrt(E_mu * E_mu - p_mu * p_mu);
double th_nu_mu = pnu.Vect().Angle(pmu.Vect());
double E_pi0 = ppi0.E() / 1000;
double p_pi0 = ppi0.Vect().Mag() / 1000;
double m_pi0 = sqrt(E_pi0 * E_pi0 - p_pi0 * p_pi0);
double th_nu_pi0 = pnu.Vect().Angle(ppi0.Vect());
const double m_n = PhysConst::mass_neutron; // neutron mass
const double m_p = PhysConst::mass_proton; // proton mass
double th_pi0_mu = ppi0.Vect().Angle(pmu.Vect());
double rEnu = (m_mu * m_mu + m_pi0 * m_pi0 + m_n * m_n - m_p * m_p -
2 * m_n * (E_pi0 + E_mu) + 2 * E_pi0 * E_mu -
2 * p_pi0 * p_mu * cos(th_pi0_mu)) /
(2 * (E_pi0 + E_mu - p_pi0 * cos(th_nu_pi0) -
p_mu * cos(th_nu_mu) - m_n));
return rEnu;
};
//********************************************************************
// Reconstruct Q2 for CC1pi0
// Beware: uses true Enu, not reconstructed Enu
double FitUtils::Q2CC1pi0rec(TLorentzVector pnu, TLorentzVector pmu,
TLorentzVector ppi0) {
//********************************************************************
double E_mu = pmu.E() / 1000.; // energy of lepton in GeV
double p_mu = pmu.Vect().Mag() / 1000.; // momentum of lepton
double m_mu = sqrt(E_mu * E_mu - p_mu * p_mu); // lepton mass
double th_nu_mu = pnu.Vect().Angle(pmu.Vect());
// double rEnu = EnuCC1pi0rec(pnu, pmu, ppi0); //reconstructed neutrino energy
// Use true neutrino energy
double rEnu = pnu.E() / 1000.;
double q2 = -m_mu * m_mu + 2. * rEnu * (E_mu - p_mu * cos(th_nu_mu));
return q2;
};
//********************************************************************
// Reconstruct Enu for CC1pi+
// pionInfo reflects if we're using pion kinematics or not
// In T2K CC1pi+ CH the Michel tag is used for pion in which pion kinematic info
// is lost and Enu is reconstructed without pion kinematics
double FitUtils::EnuCC1piprec(TLorentzVector pnu, TLorentzVector pmu,
TLorentzVector ppi, bool pionInfo) {
//********************************************************************
double E_mu = pmu.E() / 1000.;
double p_mu = pmu.Vect().Mag() / 1000.;
double m_mu = sqrt(E_mu * E_mu - p_mu * p_mu);
double E_pi = ppi.E() / 1000.;
double p_pi = ppi.Vect().Mag() / 1000.;
double m_pi = sqrt(E_pi * E_pi - p_pi * p_pi);
const double m_n = PhysConst::mass_neutron; // neutron/proton mass
// should really take proton mass for proton interaction, neutron for neutron
// interaction. However, difference is pretty much negligable here!
// need this angle too
double th_nu_pi = pnu.Vect().Angle(ppi.Vect());
double th_nu_mu = pnu.Vect().Angle(pmu.Vect());
double th_pi_mu = ppi.Vect().Angle(pmu.Vect());
double rEnu = -999;
// If experiment doesn't have pion kinematic information (T2K CC1pi+ CH
// example of this; Michel e sample has no directional information on pion)
// We'll need to modify the reconstruction variables
if (!pionInfo) {
// Assumes that pion angle contribution contributes net zero
// Also assumes the momentum of reconstructed pions when using Michel
// electrons is 130 MeV
// So need to redefine E_pi and p_pi
// It's a little unclear what happens to the pmu . ppi term (4-vector); do
// we include the angular contribution?
// This below doesn't
th_nu_pi = M_PI / 2.;
p_pi = 0.130;
E_pi = sqrt(p_pi * p_pi + m_pi * m_pi);
// rEnu = (m_mu*m_mu + m_pi*m_pi - 2*m_n*(E_mu + E_pi) + 2*E_mu*E_pi -
// 2*p_mu*p_pi) / (2*(E_mu + E_pi - p_mu*cos(th_nu_mu) - m_n));
}
rEnu =
(m_mu * m_mu + m_pi * m_pi - 2 * m_n * (E_pi + E_mu) + 2 * E_pi * E_mu -
2 * p_pi * p_mu * cos(th_pi_mu)) /
(2 * (E_pi + E_mu - p_pi * cos(th_nu_pi) - p_mu * cos(th_nu_mu) - m_n));
return rEnu;
};
//********************************************************************
// Reconstruct neutrino energy from outgoing particles; will differ from the
// actual neutrino energy. Here we use assumption of a Delta resonance
double FitUtils::EnuCC1piprecDelta(TLorentzVector pnu, TLorentzVector pmu) {
//********************************************************************
const double m_Delta =
PhysConst::mass_delta; // PDG value for Delta mass in GeV
const double m_n = PhysConst::mass_neutron; // neutron/proton mass
// should really take proton mass for proton interaction, neutron for neutron
// interaction. However, difference is pretty much negligable here!
double E_mu = pmu.E() / 1000;
double p_mu = pmu.Vect().Mag() / 1000;
double m_mu = sqrt(E_mu * E_mu - p_mu * p_mu);
double th_nu_mu = pnu.Vect().Angle(pmu.Vect());
double rEnu = (m_Delta * m_Delta - m_n * m_n - m_mu * m_mu + 2 * m_n * E_mu) /
(2 * (m_n - E_mu + p_mu * cos(th_nu_mu)));
return rEnu;
};
// MOVE TO T2K UTILS!
//********************************************************************
// Reconstruct Enu using "extended MiniBooNE" as defined in Raquel's T2K TN
//
// Supposedly includes pion direction and binding energy of target nucleon
// I'm not convinced (yet), maybe
double FitUtils::EnuCC1piprec_T2K_eMB(TLorentzVector pnu, TLorentzVector pmu,
TLorentzVector ppi) {
//********************************************************************
// Unit vector for neutrino momentum
TVector3 p_nu_vect_unit = pnu.Vect() * (1. / pnu.E());
double E_mu = pmu.E() / 1000.;
TVector3 p_mu_vect = pmu.Vect() * (1. / 1000.);
double E_pi = ppi.E() / 1000.;
TVector3 p_pi_vect = ppi.Vect() * (1. / 1000.);
double E_bind =
27. / 1000.; // This should be roughly correct for CH; but not clear!
double m_p = PhysConst::mass_proton;
// Makes life a little easier, gonna square this one
double a1 = m_p - E_bind - E_mu - E_pi;
// Just gets the magnitude square of the muon and pion momentum vectors
double a2 = (p_mu_vect + p_pi_vect).Mag2();
// Gets the somewhat complicated scalar product between neutrino and
// (p_mu+p_pi), scaled to Enu
double a3 = p_nu_vect_unit * (p_mu_vect + p_pi_vect);
double rEnu =
(m_p * m_p - a1 * a1 + a2) / (2 * (m_p - E_bind - E_mu - E_pi + a3));
return rEnu;
}
//********************************************************************
// Reconstructed Q2 for CC1pi+
//
// enuType describes how the neutrino energy is reconstructed
// 0 uses true neutrino energy; case for MINERvA and MiniBooNE
// 1 uses "extended MiniBooNE" reconstructed (T2K CH)
// 2 uses "MiniBooNE" reconstructed (EnuCC1piprec with pionInfo = true) (T2K CH)
// "MiniBooNE" reconstructed (EnuCC1piprec with pionInfo = false, the
// case for T2K when using Michel tag) (T2K CH)
// 3 uses Delta for reconstruction (T2K CH)
double FitUtils::Q2CC1piprec(TLorentzVector pnu, TLorentzVector pmu,
TLorentzVector ppi, int enuType, bool pionInfo) {
//********************************************************************
double E_mu = pmu.E() / 1000.; // energy of lepton in GeV
double p_mu = pmu.Vect().Mag() / 1000.; // momentum of lepton
double m_mu = sqrt(E_mu * E_mu - p_mu * p_mu); // lepton mass
double th_nu_mu = pnu.Vect().Angle(pmu.Vect());
// Different ways of reconstructing the neutrino energy; default is just to
// use the truth (case 0)
double rEnu = -999;
switch (enuType) {
case 0: // True neutrino energy, default; this is the case for when Q2
// defined as Q2 true (MiniBooNE, MINERvA)
rEnu = pnu.E() / 1000.;
break;
case 1: // Extended MiniBooNE reconstructed, as defined by Raquel's CC1pi+
// CH T2K analysis
// Definitely uses pion info :)
rEnu = EnuCC1piprec_T2K_eMB(pnu, pmu, ppi);
break;
case 2: // MiniBooNE reconstructed, as defined by MiniBooNE and Raquel's
// CC1pi+ CH T2K analysis and Linda's CC1pi+ H2O
// Can have this with and without pion info, depending on if Michel tag
// used (Raquel)
rEnu = EnuCC1piprec(pnu, pmu, ppi, pionInfo);
break;
case 3:
rEnu = EnuCC1piprecDelta(pnu, pmu);
break;
} // No need for default here since default value of enuType = 0, defined in
// header
double q2 = -m_mu * m_mu + 2. * rEnu * (E_mu - p_mu * cos(th_nu_mu));
return q2;
};
//********************************************************************
// Returns the reconstructed W from a nucleon and an outgoing pion
//
// Could do this a lot more clever (pp + ppi).Mag() would do the job, but this
// would be less instructive
//********************************************************************
double FitUtils::MpPi(TLorentzVector pp, TLorentzVector ppi) {
double E_p = pp.E();
double p_p = pp.Vect().Mag();
double m_p = sqrt(E_p * E_p - p_p * p_p);
double E_pi = ppi.E();
double p_pi = ppi.Vect().Mag();
double m_pi = sqrt(E_pi * E_pi - p_pi * p_pi);
double th_p_pi = pp.Vect().Angle(ppi.Vect());
// fairly easy thing to derive since bubble chambers measure the proton!
double invMass = sqrt(m_p * m_p + m_pi * m_pi + 2 * E_p * E_pi -
2 * p_pi * p_p * cos(th_p_pi));
return invMass;
};
//********************************************************
// Reconstruct the hadronic mass using neutrino and muon
// Could technically do E_nu = EnuCC1pipRec(pnu,pmu,ppi) too, but this wwill be
// reconstructed Enu; so gives reconstructed Wrec which most of the time isn't
// used!
//
// Only MINERvA uses this so far; and the Enu is Enu_true
// If we want W_true need to take initial state nucleon motion into account
// Return value is in MeV!!!
double FitUtils::Wrec(TLorentzVector pnu, TLorentzVector pmu) {
//********************************************************
double E_mu = pmu.E();
double p_mu = pmu.Vect().Mag();
double m_mu = sqrt(E_mu * E_mu - p_mu * p_mu);
double th_nu_mu = pnu.Vect().Angle(pmu.Vect());
// The factor of 1000 is necessary for downstream functions
const double m_p = PhysConst::mass_proton * 1000;
// MINERvA cut on W_exp which is tuned to W_true; so use true Enu from
// generators
double E_nu = pnu.E();
double w_rec = sqrt(m_p * m_p + m_mu * m_mu - 2 * m_p * E_mu +
2 * E_nu * (m_p - E_mu + p_mu * cos(th_nu_mu)));
return w_rec;
};
//********************************************************
// Reconstruct the true hadronic mass using the initial state and muon
// Could technically do E_nu = EnuCC1pipRec(pnu,pmu,ppi) too, but this wwill be
// reconstructed Enu; so gives reconstructed Wrec which most of the time isn't
// used!
//
// No one seems to use this because it's fairly MC dependent!
// Return value is in MeV!!!
double FitUtils::Wtrue(TLorentzVector pnu, TLorentzVector pmu,
TLorentzVector pnuc) {
//********************************************************
// Could simply do the TLorentzVector operators here but this is more
// instructive?
// ... and prone to errors ...
double E_mu = pmu.E();
double p_mu = pmu.Vect().Mag();
double m_mu = sqrt(E_mu * E_mu - p_mu * p_mu);
double th_nu_mu = pnu.Vect().Angle(pmu.Vect());
double E_nuc = pnuc.E();
double p_nuc = pnuc.Vect().Mag();
double m_nuc = sqrt(E_nuc * E_nuc - p_nuc * p_nuc);
double th_nuc_mu = pmu.Vect().Angle(pnuc.Vect());
double th_nu_nuc = pnu.Vect().Angle(pnuc.Vect());
double E_nu = pnu.E();
double w_rec = sqrt(m_nuc * m_nuc + m_mu * m_mu - 2 * E_nu * E_mu +
2 * E_nu * p_mu * cos(th_nu_mu) - 2 * E_nuc * E_mu +
2 * p_nuc * p_mu * cos(th_nuc_mu) + 2 * E_nu * E_nuc -
2 * E_nu * p_nuc * cos(th_nu_nuc));
return w_rec;
};
double FitUtils::SumKE_PartVect(std::vector const fps) {
double sum = 0.0;
for (size_t p_it = 0; p_it < fps.size(); ++p_it) {
sum += fps[p_it]->KE();
}
return sum;
}
double FitUtils::SumTE_PartVect(std::vector const fps) {
double sum = 0.0;
for (size_t p_it = 0; p_it < fps.size(); ++p_it) {
sum += fps[p_it]->E();
}
return sum;
}
/*
E Recoil
*/
double FitUtils::GetErecoil_TRUE(FitEvent *event) {
// Get total energy of hadronic system.
double Erecoil = 0.0;
for (unsigned int i = 2; i < event->Npart(); i++) {
// Only final state
if (!event->PartInfo(i)->fIsAlive) continue;
if (event->PartInfo(i)->fNEUTStatusCode != 0) continue;
// Skip Lepton
if (abs(event->PartInfo(i)->fPID) == abs(event->PartInfo(0)->fPID) - 1)
continue;
// Add Up KE of protons and TE of everything else
if (event->PartInfo(i)->fPID == 2212 || event->PartInfo(i)->fPID == 2112) {
Erecoil +=
fabs(event->PartInfo(i)->fP.E()) - fabs(event->PartInfo(i)->fP.Mag());
} else {
Erecoil += event->PartInfo(i)->fP.E();
}
}
return Erecoil;
}
double FitUtils::GetErecoil_CHARGED(FitEvent *event) {
// Get total energy of hadronic system.
double Erecoil = 0.0;
for (unsigned int i = 2; i < event->Npart(); i++) {
// Only final state
if (!event->PartInfo(i)->fIsAlive) continue;
if (event->PartInfo(i)->fNEUTStatusCode != 0) continue;
// Skip Lepton
if (abs(event->PartInfo(i)->fPID) == abs(event->PartInfo(0)->fPID) - 1)
continue;
// Skip Neutral particles
if (event->PartInfo(i)->fPID == 2112 || event->PartInfo(i)->fPID == 111 ||
event->PartInfo(i)->fPID == 22)
continue;
// Add Up KE of protons and TE of everything else
if (event->PartInfo(i)->fPID == 2212) {
Erecoil +=
fabs(event->PartInfo(i)->fP.E()) - fabs(event->PartInfo(i)->fP.Mag());
} else {
Erecoil += event->PartInfo(i)->fP.E();
}
}
return Erecoil;
}
// MOVE TO MINERVA Utils!
double FitUtils::GetErecoil_MINERvA_LowRecoil(FitEvent *event) {
// Get total energy of hadronic system.
double Erecoil = 0.0;
for (unsigned int i = 2; i < event->Npart(); i++) {
// Only final state
if (!event->PartInfo(i)->fIsAlive) continue;
if (event->PartInfo(i)->fNEUTStatusCode != 0) continue;
// Skip Lepton
if (abs(event->PartInfo(i)->fPID) == 13) continue;
// Skip Neutrons particles
if (event->PartInfo(i)->fPID == 2112) continue;
int PID = event->PartInfo(i)->fPID;
// KE of Protons and charged pions
if (PID == 2212 or PID == 211 or PID == -211) {
// Erecoil += FitUtils::T(event->PartInfo(i)->fP);
Erecoil +=
fabs(event->PartInfo(i)->fP.E()) - fabs(event->PartInfo(i)->fP.Mag());
// Total Energy of non-neutrons
// } else if (PID != 2112 and PID < 999 and PID != 22 and abs(PID) !=
// 14) {
} else if (PID == 111 || PID == 11 || PID == -11 || PID == 22) {
Erecoil += (event->PartInfo(i)->fP.E());
}
}
return Erecoil;
}
// MOVE TO MINERVA Utils!
// The alternative Eavailble definition takes true q0 and subtracts the kinetic energy of neutrons and pion masses
// returns in MeV
double FitUtils::Eavailable(FitEvent *event) {
double Eav = 0.0;
// Now take q0 and subtract Eav
double q0 = event->GetNeutrinoIn()->fP.E();
if (event->GetHMFSParticle(13)) {
q0 -= event->GetHMFSParticle(13)->fP.E();
} else if (event->GetHMFSParticle(-13)) {
q0 -= event->GetHMFSParticle(-13)->fP.E();
} else if (event->GetHMFSParticle(14)) {
q0 -= event->GetHMFSParticle(14)->fP.E();
} else if (event->GetHMFSParticle(-14)) {
q0 -= event->GetHMFSParticle(-14)->fP.E();
}else {
std::cerr << "Found no Muon or Muon Neutrino" << std::endl;
}
for (unsigned int i = 2; i < event->Npart(); i++) {
// Only final state
if (!event->PartInfo(i)->fIsAlive) continue;
if (event->PartInfo(i)->fNEUTStatusCode != 0) continue;
int PID = event->PartInfo(i)->fPID;
// Neutrons
if (PID == 2112) {
// Adding kinetic energy of neutron
Eav += FitUtils::T(event->PartInfo(i)->fP)*1000.;
// All pion masses
} else if (abs(PID) == 211 || PID == 111) {
Eav += event->PartInfo(i)->fP.M();
}
}
return q0-Eav;
}
TVector3 GetVectorInTPlane(const TVector3 &inp, const TVector3 &planarNormal) {
TVector3 pnUnit = planarNormal.Unit();
double inpProjectPN = inp.Dot(pnUnit);
TVector3 InPlanarInput = inp - (inpProjectPN * pnUnit);
return InPlanarInput;
}
TVector3 GetUnitVectorInTPlane(const TVector3 &inp,
const TVector3 &planarNormal) {
return GetVectorInTPlane(inp, planarNormal).Unit();
}
Double_t GetDeltaPhiT(TVector3 const &V_lepton, TVector3 const &V_other,
TVector3 const &Normal, bool PiMinus = false) {
TVector3 V_lepton_T = GetUnitVectorInTPlane(V_lepton, Normal);
TVector3 V_other_T = GetUnitVectorInTPlane(V_other, Normal);
return PiMinus ? acos(V_lepton_T.Dot(V_other_T))
: (M_PI - acos(V_lepton_T.Dot(V_other_T)));
}
TVector3 GetDeltaPT(TVector3 const &V_lepton, TVector3 const &V_other,
TVector3 const &Normal) {
TVector3 V_lepton_T = GetVectorInTPlane(V_lepton, Normal);
TVector3 V_other_T = GetVectorInTPlane(V_other, Normal);
return V_lepton_T + V_other_T;
}
Double_t GetDeltaAlphaT(TVector3 const &V_lepton, TVector3 const &V_other,
TVector3 const &Normal, bool PiMinus = false) {
TVector3 DeltaPT = GetDeltaPT(V_lepton, V_other, Normal);
return GetDeltaPhiT(V_lepton, DeltaPT, Normal, PiMinus);
}
double FitUtils::Get_STV_dpt(FitEvent *event, int ISPDG, bool Is0pi) {
// Check that the neutrino exists
if (event->NumISParticle(ISPDG) == 0) {
return -9999;
}
// Return 0 if the proton or muon are missing
if (event->NumFSParticle(2212) == 0 ||
event->NumFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1)) == 0) {
return -9999;
}
// Now get the TVector3s for each particle
TVector3 const &NuP = event->GetHMISParticle(14)->fP.Vect();
TVector3 const &LeptonP =
event->GetHMFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1))->fP.Vect();
// Find the highest momentum proton in the event between 450 and 1200 MeV with theta_p < 70
TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
int HMFSProton = 0;
double HighestMomentum = 0.0;
// Get the stack of protons
std::vector Protons = event->GetAllFSProton();
for (size_t i = 0; i < Protons.size(); ++i) {
if (Protons[i]->p() > 450 &&
Protons[i]->p() < 1200 &&
Protons[i]->P3().Angle(Pnu.Vect()) < (M_PI/180.0)*70.0 &&
Protons[i]->p() > HighestMomentum) {
HighestMomentum = Protons[i]->p();
HMFSProton = i;
}
}
// Now get the proton
TLorentzVector Pprot = Protons[HMFSProton]->fP;
// Get highest momentum proton in allowed proton range
TVector3 HadronP = Pprot.Vect();
// If we don't have a CC0pi signal definition we also add in pion momentum
if (!Is0pi) {
if (event->NumFSParticle(PhysConst::pdg_pions) == 0) {
return -9999;
}
// Count up pion momentum
TLorentzVector pp = event->GetHMFSParticle(PhysConst::pdg_pions)->fP;
HadronP += pp.Vect();
}
return GetDeltaPT(LeptonP, HadronP, NuP).Mag();
}
double FitUtils::Get_STV_dphit(FitEvent *event, int ISPDG, bool Is0pi) {
// Check that the neutrino exists
if (event->NumISParticle(ISPDG) == 0) {
return -9999;
}
// Return 0 if the proton or muon are missing
if (event->NumFSParticle(2212) == 0 ||
event->NumFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1)) == 0) {
return -9999;
}
// Now get the TVector3s for each particle
TVector3 const &NuP = event->GetHMISParticle(ISPDG)->fP.Vect();
TVector3 const &LeptonP =
event->GetHMFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1))->fP.Vect();
// Find the highest momentum proton in the event between 450 and 1200 MeV with theta_p < 70
TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
int HMFSProton = 0;
double HighestMomentum = 0.0;
// Get the stack of protons
std::vector Protons = event->GetAllFSProton();
for (size_t i = 0; i < Protons.size(); ++i) {
if (Protons[i]->p() > 450 &&
Protons[i]->p() < 1200 &&
Protons[i]->P3().Angle(Pnu.Vect()) < (M_PI/180.0)*70.0 &&
Protons[i]->p() > HighestMomentum) {
HighestMomentum = Protons[i]->p();
HMFSProton = i;
}
}
// Now get the proton
TLorentzVector Pprot = Protons[HMFSProton]->fP;
// Get highest momentum proton in allowed proton range
TVector3 HadronP = Pprot.Vect();
if (!Is0pi) {
if (event->NumFSParticle(PhysConst::pdg_pions) == 0) {
return -9999;
}
TLorentzVector pp = event->GetHMFSParticle(PhysConst::pdg_pions)->fP;
HadronP += pp.Vect();
}
return GetDeltaPhiT(LeptonP, HadronP, NuP);
}
double FitUtils::Get_STV_dalphat(FitEvent *event, int ISPDG, bool Is0pi) {
// Check that the neutrino exists
if (event->NumISParticle(ISPDG) == 0) {
return -9999;
}
// Return 0 if the proton or muon are missing
if (event->NumFSParticle(2212) == 0 ||
event->NumFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1)) == 0) {
return -9999;
}
// Now get the TVector3s for each particle
TVector3 const &NuP = event->GetHMISParticle(ISPDG)->fP.Vect();
TVector3 const &LeptonP =
event->GetHMFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1))->fP.Vect();
// Find the highest momentum proton in the event between 450 and 1200 MeV with theta_p < 70
TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
int HMFSProton = 0;
double HighestMomentum = 0.0;
// Get the stack of protons
std::vector Protons = event->GetAllFSProton();
for (size_t i = 0; i < Protons.size(); ++i) {
if (Protons[i]->p() > 450 &&
Protons[i]->p() < 1200 &&
Protons[i]->P3().Angle(Pnu.Vect()) < (M_PI/180.0)*70.0 &&
Protons[i]->p() > HighestMomentum) {
HighestMomentum = Protons[i]->p();
HMFSProton = i;
}
}
// Now get the proton
TLorentzVector Pprot = Protons[HMFSProton]->fP;
// Get highest momentum proton in allowed proton range
TVector3 HadronP = Pprot.Vect();
if (!Is0pi) {
if (event->NumFSParticle(PhysConst::pdg_pions) == 0) {
return -9999;
}
TLorentzVector pp = event->GetHMFSParticle(PhysConst::pdg_pions)->fP;
HadronP += pp.Vect();
}
return GetDeltaAlphaT(LeptonP, HadronP, NuP);
}
// As defined in PhysRevC.95.065501
// Using prescription from arXiv 1805.05486
// Returns in GeV
double FitUtils::Get_pn_reco_C(FitEvent *event, int ISPDG, bool Is0pi) {
const double mn = PhysConst::mass_neutron; // neutron mass
const double mp = PhysConst::mass_proton; // proton mass
// Check that the neutrino exists
if (event->NumISParticle(ISPDG) == 0) {
return -9999;
}
// Return 0 if the proton or muon are missing
if (event->NumFSParticle(2212) == 0 ||
event->NumFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1)) == 0) {
return -9999;
}
// Now get the TVector3s for each particle
TVector3 const &NuP = event->GetHMISParticle(14)->fP.Vect();
TVector3 const &LeptonP =
event->GetHMFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1))->fP.Vect();
// Find the highest momentum proton in the event between 450 and 1200 MeV with theta_p < 70
TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
int HMFSProton = 0;
double HighestMomentum = 0.0;
// Get the stack of protons
std::vector Protons = event->GetAllFSProton();
for (size_t i = 0; i < Protons.size(); ++i) {
// Update the highest momentum particle
if (Protons[i]->p() > 450 &&
Protons[i]->p() < 1200 &&
Protons[i]->P3().Angle(Pnu.Vect()) < (M_PI/180.0)*70.0 &&
Protons[i]->p() > HighestMomentum) {
HighestMomentum = Protons[i]->p();
HMFSProton = i;
}
}
// Now get the proton
TLorentzVector Pprot = Protons[HMFSProton]->fP;
// Get highest momentum proton in allowed proton range
TVector3 HadronP = Pprot.Vect();
//TVector3 HadronP = event->GetHMFSParticle(2212)->fP.Vect();
double const el = event->GetHMFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1))->E()/1000.;
double const eh = Pprot.E()/1000.;
if (!Is0pi) {
if (event->NumFSParticle(PhysConst::pdg_pions) == 0) {
return -9999;
}
TLorentzVector pp = event->GetHMFSParticle(PhysConst::pdg_pions)->fP;
HadronP += pp.Vect();
}
TVector3 dpt = GetDeltaPT(LeptonP, HadronP, NuP);
double dptMag = dpt.Mag()/1000.;
double ma = 6*mn + 6*mp - 0.09216; // target mass (E is from PhysRevC.95.065501)
double map = ma - mn + 0.02713; // remnant mass
double pmul = LeptonP.Dot(NuP.Unit())/1000.;
double phl = HadronP.Dot(NuP.Unit())/1000.;
//double pmul = GetVectorInTPlane(LeptonP, dpt).Mag()/1000.;
//double phl = GetVectorInTPlane(HadronP, dpt).Mag()/1000.;
double R = ma + pmul + phl - el - eh;
double dpl = 0.5*R - (map*map + dptMag*dptMag)/(2*R);
//double dpl = ((R*R)-(dptMag*dptMag)-(map*map))/(2*R); // as in in PhysRevC.95.065501 - gives same result
double pn_reco = sqrt((dptMag*dptMag) + (dpl*dpl));
//std::cout << "Diagnostics: " << std::endl;
//std::cout << "mn: " << mn << std::endl;
//std::cout << "ma: " << ma << std::endl;
//std::cout << "map: " << map << std::endl;
//std::cout << "pmu: " << LeptonP.Mag()/1000. << std::endl;
//std::cout << "ph: " << HadronP.Mag()/1000. << std::endl;
//std::cout << "pmul: " << pmul << std::endl;
//std::cout << "phl: " << phl << std::endl;
//std::cout << "el: " << el << std::endl;
//std::cout << "eh: " << eh << std::endl;
//std::cout << "R: " << R << std::endl;
//std::cout << "dptMag: " << dptMag << std::endl;
//std::cout << "dpl: " << dpl << std::endl;
//std::cout << "pn_reco: " << pn_reco << std::endl;
return pn_reco;
}
+double FitUtils::Get_pn_reco_Ar(FitEvent *event, int ISPDG, bool Is0pi) {
+
+ const double mn = PhysConst::mass_neutron; // neutron mass
+ const double mp = PhysConst::mass_proton; // proton mass
+
+ // Check that the neutrino exists
+ if (event->NumISParticle(ISPDG) == 0) {
+ return -9999;
+ }
+ // Return 0 if the proton or muon are missing
+ if (event->NumFSParticle(2212) == 0 ||
+ event->NumFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1)) == 0) {
+ return -9999;
+ }
+
+ // Now get the TVector3s for each particle
+ TVector3 const &NuP = event->GetHMISParticle(14)->fP.Vect();
+ TVector3 const &LeptonP =
+ event->GetHMFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1))->fP.Vect();
+
+ // Find the highest momentum proton in the event between 450 and 1200 MeV with theta_p < 70
+ TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
+ int HMFSProton = 0;
+ double HighestMomentum = 0.0;
+ // Get the stack of protons
+ std::vector Protons = event->GetAllFSProton();
+ for (size_t i = 0; i < Protons.size(); ++i) {
+ // Update the highest momentum particle
+ if (Protons[i]->p() > 450 &&
+ Protons[i]->p() < 1200 &&
+ Protons[i]->P3().Angle(Pnu.Vect()) < (M_PI/180.0)*70.0 &&
+ Protons[i]->p() > HighestMomentum) {
+ HighestMomentum = Protons[i]->p();
+ HMFSProton = i;
+ }
+ }
+ // Now get the proton
+ TLorentzVector Pprot = Protons[HMFSProton]->fP;
+ // Get highest momentum proton in allowed proton range
+ TVector3 HadronP = Pprot.Vect();
+ //TVector3 HadronP = event->GetHMFSParticle(2212)->fP.Vect();
+
+ double const el = event->GetHMFSParticle(ISPDG + ((ISPDG < 0) ? 1 : -1))->E()/1000.;
+ double const eh = Pprot.E()/1000.;
+
+ if (!Is0pi) {
+ if (event->NumFSParticle(PhysConst::pdg_pions) == 0) {
+ return -9999;
+ }
+ TLorentzVector pp = event->GetHMFSParticle(PhysConst::pdg_pions)->fP;
+ HadronP += pp.Vect();
+ }
+ TVector3 dpt = GetDeltaPT(LeptonP, HadronP, NuP);
+ double dptMag = dpt.Mag()/1000.;
+
+ //double ma = 6*mn + 6*mp - 0.09216; // target mass (E is from PhysRevC.95.065501)
+ //double map = ma - mn + 0.02713; // remnant mass
+ double ma = 6*mn + 6*mp - 0.34381; // target mass (E is from PhysRevC.95.065501)
+ double map = ma - mn + 0.02713; // remnant mass
+
+ double pmul = LeptonP.Dot(NuP.Unit())/1000.;
+ double phl = HadronP.Dot(NuP.Unit())/1000.;
+
+ //double pmul = GetVectorInTPlane(LeptonP, dpt).Mag()/1000.;
+ //double phl = GetVectorInTPlane(HadronP, dpt).Mag()/1000.;
+
+ double R = ma + pmul + phl - el - eh;
+
+ double dpl = 0.5*R - (map*map + dptMag*dptMag)/(2*R);
+ //double dpl = ((R*R)-(dptMag*dptMag)-(map*map))/(2*R); // as in in PhysRevC.95.065501 - gives same result
+
+ double pn_reco = sqrt((dptMag*dptMag) + (dpl*dpl));
+
+ //std::cout << "Diagnostics: " << std::endl;
+ //std::cout << "mn: " << mn << std::endl;
+ //std::cout << "ma: " << ma << std::endl;
+ //std::cout << "map: " << map << std::endl;
+ //std::cout << "pmu: " << LeptonP.Mag()/1000. << std::endl;
+ //std::cout << "ph: " << HadronP.Mag()/1000. << std::endl;
+ //std::cout << "pmul: " << pmul << std::endl;
+ //std::cout << "phl: " << phl << std::endl;
+ //std::cout << "el: " << el << std::endl;
+ //std::cout << "eh: " << eh << std::endl;
+ //std::cout << "R: " << R << std::endl;
+ //std::cout << "dptMag: " << dptMag << std::endl;
+ //std::cout << "dpl: " << dpl << std::endl;
+ //std::cout << "pn_reco: " << pn_reco << std::endl;
+
+ return pn_reco;
+}
+
// Get Cos theta with Adler angles
double FitUtils::CosThAdler(TLorentzVector Pnu, TLorentzVector Pmu, TLorentzVector Ppi, TLorentzVector Pprot) {
// Get the "resonance" lorentz vector (pion proton system)
TLorentzVector Pres = Pprot + Ppi;
// Boost the particles into the resonance rest frame so we can define the x,y,z axis
Pnu.Boost(Pres.BoostVector());
Pmu.Boost(-Pres.BoostVector());
Ppi.Boost(-Pres.BoostVector());
// The z-axis is defined as the axis of three-momentum transfer, \vec{k}
// Also unit normalise the axis
TVector3 zAxis = (Pnu.Vect()-Pmu.Vect())*(1.0/((Pnu.Vect()-Pmu.Vect()).Mag()));
// Then the angle between the pion in the "resonance" rest-frame and the z-axis is the theta Adler angle
double costhAdler = cos(Ppi.Vect().Angle(zAxis));
return costhAdler;
}
// Get phi with Adler angles, a bit more complicated...
double FitUtils::PhiAdler(TLorentzVector Pnu, TLorentzVector Pmu, TLorentzVector Ppi, TLorentzVector Pprot) {
// Get the "resonance" lorentz vector (pion proton system)
TLorentzVector Pres = Pprot + Ppi;
// Boost the particles into the resonance rest frame so we can define the x,y,z axis
Pnu.Boost(Pres.BoostVector());
Pmu.Boost(-Pres.BoostVector());
Ppi.Boost(-Pres.BoostVector());
// The z-axis is defined as the axis of three-momentum transfer, \vec{k}
// Also unit normalise the axis
TVector3 zAxis = (Pnu.Vect()-Pmu.Vect())*(1.0/((Pnu.Vect()-Pmu.Vect()).Mag()));
// The y-axis is then defined perpendicular to z and muon momentum in the resonance frame
TVector3 yAxis = Pnu.Vect().Cross(Pmu.Vect());
yAxis *= 1.0/double(yAxis.Mag());
// And the x-axis is then simply perpendicular to z and x
TVector3 xAxis = yAxis.Cross(zAxis);
xAxis *= 1.0/double(xAxis.Mag());
double x = Ppi.Vect().Dot(xAxis);
double y = Ppi.Vect().Dot(yAxis);
//double z = Ppi.Vect().Dot(zAxis);
double newphi = atan2(y, x)*(180./M_PI);
// Convert negative angles to positive
if (newphi < 0.0) newphi += 360.0;
// Old silly method before atan2
/*
// Then finally construct phi as the angle between pion projection and x axis
// Get the project of the pion momentum on to the zaxis
TVector3 PiVectZ = zAxis*Ppi.Vect().Dot(zAxis);
// The subtract the projection off the pion vector to get to get the plane
TVector3 PiPlane = Ppi.Vect() - PiVectZ;
double phi = -999.99;
if (PiPlane.Y() > 0) {
phi = (180./M_PI)*PiPlane.Angle(xAxis);
} else if (PiPlane.Y() < 0) {
phi = (180./M_PI)*(2*M_PI-PiPlane.Angle(xAxis));
} else if (PiPlane.Y() == 0) {
TRandom3 rand;
double randNo = rand.Rndm();
if (randNo > 0.5) {
phi = (180./M_PI)*PiPlane.Angle(xAxis);
} else {
phi = (180./M_PI)*(2*M_PI-PiPlane.Angle(xAxis));
}
}
*/
return newphi;
}
//********************************************************************
double FitUtils::ppInfK(TLorentzVector pmu, double costh, double binding,
bool neutrino) {
//********************************************************************
// Convert all values to GeV
//const double V = binding / 1000.; // binding potential
//const double mn = PhysConst::mass_neutron; // neutron mass
const double mp = PhysConst::mass_proton; // proton mass
double el = pmu.E() / 1000.;
//double pl = (pmu.Vect().Mag()) / 1000.; // momentum of lepton
double enu = EnuQErec(pmu, costh, binding, neutrino);
double ep_inf = enu - el + mp;
double pp_inf = sqrt(ep_inf * ep_inf - mp * mp);
return pp_inf;
};
//********************************************************************
TVector3 FitUtils::tppInfK(TLorentzVector pmu, double costh, double binding,
bool neutrino) {
//********************************************************************
// Convert all values to GeV
//const double V = binding / 1000.; // binding potential
//const double mn = PhysConst::mass_neutron; // neutron mass
//const double mp = PhysConst::mass_proton; // proton mass
double pl_x = pmu.X() / 1000.;
double pl_y = pmu.Y() / 1000.;
double pl_z= pmu.Z() / 1000.;
double enu = EnuQErec(pmu, costh, binding, neutrino);
TVector3 tpp_inf(-pl_x, -pl_y, -pl_z+enu);
return tpp_inf;
};
//********************************************************************
double FitUtils::cthpInfK(TLorentzVector pmu, double costh, double binding,
bool neutrino) {
//********************************************************************
// Convert all values to GeV
//const double V = binding / 1000.; // binding potential
//const double mn = PhysConst::mass_neutron; // neutron mass
const double mp = PhysConst::mass_proton; // proton mass
double el = pmu.E() / 1000.;
double pl = (pmu.Vect().Mag()) / 1000.; // momentum of lepton
double enu = EnuQErec(pmu, costh, binding, neutrino);
double ep_inf = enu - el + mp;
double pp_inf = sqrt(ep_inf * ep_inf - mp * mp);
double cth_inf = (enu*enu + pp_inf*pp_inf - pl*pl)/(2*enu*pp_inf);
return cth_inf;
};
diff --git a/src/Utils/FitUtils.h b/src/Utils/FitUtils.h
index 9510e20..ccfe8bd 100644
--- a/src/Utils/FitUtils.h
+++ b/src/Utils/FitUtils.h
@@ -1,186 +1,187 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#ifndef FITUTILS_H_SEEN
#define FITUTILS_H_SEEN
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include "FitEvent.h"
#include "TGraph.h"
#include "TH2Poly.h"
#include "FitEvent.h"
#include "FitLogger.h"
/*!
* \addtogroup Utils
* @{
*/
/// Functions needed by individual samples for calculating kinematic quantities.
namespace FitUtils {
/// Return a vector of all values saved in map
double *GetArrayFromMap(std::vector invals,
std::map inmap);
/// Returns kinetic energy of particle
double T(TLorentzVector part);
/// Returns momentum of particle
double p(TLorentzVector part);
double p(FitParticle* part);
/// Returns angle between particles (_NOT_ cosine!)
double th(TLorentzVector part, TLorentzVector part2);
double th(FitParticle* part1, FitParticle* part2);
/// Hadronic mass reconstruction
double Wrec(TLorentzVector pnu, TLorentzVector pmu);
/// Hadronic mass true from initial state particles and muon; useful if the full
/// FSI vectors aren't not saved and we for some reasons need W_true
double Wtrue(TLorentzVector pnu, TLorentzVector pmu, TLorentzVector pnuc);
double SumKE_PartVect(std::vector const fps);
double SumTE_PartVect(std::vector const fps);
/// Return E Hadronic for all FS Particles in Hadronic System
double GetErecoil_TRUE(FitEvent *event);
/// Return E Hadronic for all Charged FS Particles in Hadronic System
double GetErecoil_CHARGED(FitEvent *event);
double Eavailable(FitEvent *event);
/*
CCQE MiniBooNE/MINERvA
*/
/// Function to calculate the reconstructed Q^{2}_{QE}
double Q2QErec(TLorentzVector pmu, double costh, double binding, bool neutrino = true);
/// Function returns the reconstructed E_{nu} values
double EnuQErec(TLorentzVector pmu, double costh, double binding, bool neutrino = true);
//! Function to calculate the reconstructed Q^{2}_{QE}
double Q2QErec(double pl, double costh, double binding, bool neutrino = true);
//! Function to calculate the reconstructed Q^{2}_{QE}
double Q2QErec(TLorentzVector Pmu, TLorentzVector Pnu, double binding, bool neutrino = true);
//! Function returns the reconstructed E_{nu} values
double EnuQErec(double pl, double costh, double binding,
bool neutrino = true);
/*
CCQE1p MINERvA
*/
/// Reconstruct Q2QE given just the maximum energy proton.
double ProtonQ2QErec(double pE, double binding);
/*
E Recoil MINERvA
*/
double GetErecoil_MINERvA_LowRecoil(FitEvent *event);
/*
CC1pi0 MiniBooNE
*/
/// Reconstruct Enu from CCpi0 vectors and binding energy
double EnuCC1pi0rec(TLorentzVector pnu, TLorentzVector pmu,
TLorentzVector ppi0 = TLorentzVector(0, 0, 0, 0));
/// Reconstruct Q2 from CCpi0 vectors and binding energy
double Q2CC1pi0rec(TLorentzVector pnu, TLorentzVector pmu,
TLorentzVector ppi0 = TLorentzVector(0, 0, 0, 0));
/*
CC1pi+ MiniBooNE
*/
/// returns reconstructed Enu a la MiniBooNE CCpi+
/// returns reconstructed Enu a la MiniBooNE CCpi+
// Also for when not having pion info (so when we have a Michel tag in T2K)
double EnuCC1piprec(TLorentzVector pnu, TLorentzVector pmu, TLorentzVector ppip,
bool pionInfo = true);
/// returns reconstructed Enu assumming resonance interaction where intermediate
/// resonance was a Delta
double EnuCC1piprecDelta(TLorentzVector pnu, TLorentzVector pmu);
/// returns reconstructed in a variety of flavours
double Q2CC1piprec(TLorentzVector pnu, TLorentzVector pmu, TLorentzVector ppip,
int enuType = 0, bool pionInfo = true);
/*
T2K CC1pi+ on CH
*/
double thq3pi_CC1pip_T2K(TLorentzVector pnu, TLorentzVector pmu,
TLorentzVector ppi);
double q3_CC1pip_T2K(TLorentzVector pnu, TLorentzVector pmu,
TLorentzVector ppi);
double WrecCC1pip_T2K_MB(TLorentzVector pnu, TLorentzVector pmu,
TLorentzVector ppip);
double EnuCC1piprec_T2K_eMB(TLorentzVector pnu, TLorentzVector pmu,
TLorentzVector ppi);
/*
nucleon single pion
*/
double MpPi(TLorentzVector pp, TLorentzVector ppi);
/// Gets delta p T as defined in Phys.Rev. C94 (2016) no.1, 015503
double Get_STV_dpt(FitEvent *event, int ISPDG, bool Is0pi);
/// Gets delta phi T as defined in Phys.Rev. C94 (2016) no.1, 015503
double Get_STV_dphit(FitEvent *event, int ISPDG, bool Is0pi);
/// Gets delta alpha T as defined in Phys.Rev. C94 (2016) no.1, 015503
double Get_STV_dalphat(FitEvent *event, int ISPDG, bool Is0pi);
// As defined in PhysRevC.95.065501
// Using prescription from arXiv 1805.05486
double Get_pn_reco_C(FitEvent *event, int ISPDG, bool Is0pi);
+double Get_pn_reco_Ar(FitEvent *event, int ISPDG, bool Is0pi);
//For T2K inferred kinematics analyis - variables defined as on page 7 of T2K TN287v11 (and now arXiv 1802.05078)
double ppInfK(TLorentzVector pmu, double costh, double binding, bool neutrino);
TVector3 tppInfK(TLorentzVector pmu, double costh, double binding, bool neutrino);
double cthpInfK(TLorentzVector pmu, double costh, double binding, bool neutrino);
double CosThAdler(TLorentzVector Pnu, TLorentzVector Pmu, TLorentzVector Ppi, TLorentzVector Pprot);
double PhiAdler(TLorentzVector Pnu, TLorentzVector Pmu, TLorentzVector Ppi, TLorentzVector Pprot);
}
/*! @} */
#endif