diff --git a/src/Reweight/CMakeLists.txt b/src/Reweight/CMakeLists.txt index 55b9fef..01296a8 100644 --- a/src/Reweight/CMakeLists.txt +++ b/src/Reweight/CMakeLists.txt @@ -1,83 +1,84 @@ # 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 ) 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 ) 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/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.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