diff --git a/src/Reweight/NUISANCESyst.cxx b/src/Reweight/NUISANCESyst.cxx index 6063aa3..e6cde59 100644 --- a/src/Reweight/NUISANCESyst.cxx +++ b/src/Reweight/NUISANCESyst.cxx @@ -1,61 +1,65 @@ #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 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 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 dfc84e2..7d7e1c1 100644 --- a/src/Reweight/NUISANCESyst.h +++ b/src/Reweight/NUISANCESyst.h @@ -1,62 +1,66 @@ #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, kModeNorm_NormRES, kMINERvARW_NormCCQE, kMINERvARW_NormCCMEC, kMINERvARW_NormCCRES, kMINERvARW_RikRPA_ApplyRPA, kMINERvARW_RikRPA_LowQ2, kMINERvARW_RikRPA_HighQ2, kMINERvARW_RikRESRPA_ApplyRPA, - kMINERvARW_RikRESRPA_LowQ2, + 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 4a70d6b..d21637b 100644 --- a/src/Reweight/NUISANCEWeightCalcs.cxx +++ b/src/Reweight/NUISANCEWeightCalcs.cxx @@ -1,346 +1,394 @@ #include "NUISANCEWeightCalcs.h" #include "GeneralUtils.h" #include "FitEvent.h" #include "WeightUtils.h" #include "NUISANCESyst.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; } } + + +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() { // 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"); } 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->PartInfo(0); if (!pnu){ THROW("NO Starting particle found in stack!"); } int pdgnu = pnu->fPID; FitParticle* plep = fevt->GetHMFSParticle(abs(pdgnu) - 1); 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++; } // std::cout << "PN State = " << npr << " " << nne << std::endl; if (fevt->Mode == 2 and npr == 1 and nne == 1) { initialstate = 2; } else if (fevt->Mode == 2 and ((npr == 0 and nne == 2) or (npr == 2 and nne == 0))) { initialstate = 1; } } // std::cout << "Got q0 q3 = " << q0 << " " << q3 << std::endl; // Apply weighting if (fApply_CCQE and 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 and abs(fevt->Mode) == 2) { if (fDebugStatements) std::cout << "Getting 2p2h Weight" << std::endl; rw_weight *= GetGausWeight(q0, q3, fGausVal_2p2h); } if (fApply_2p2h_PPandNN and abs(fevt->Mode) == 2 and initialstate == 1) { if (fDebugStatements) std::cout << "Getting 2p2h PPandNN Weight" << std::endl; rw_weight *= GetGausWeight(q0, q3, fGausVal_2p2h_PPandNN); } if (fApply_2p2h_NP and abs(fevt->Mode) == 2 and initialstate == 2) { if (fDebugStatements) std::cout << "Getting 2p2h NP Weight" << std::endl; rw_weight *= GetGausWeight(q0, q3, fGausVal_2p2h_NP); } if (fApply_CC1pi and abs(fevt->Mode) >= 11 and abs(fevt->Mode) <= 13) { if (fDebugStatements) std::cout << "Getting CC1pi Weight" << std::endl; rw_weight *= GetGausWeight(q0, q3, fGausVal_CC1pi); } // if (fDebugStatements) std::cout << "Returning Weight " << rw_weight << std::endl; return rw_weight; } double GaussianModeCorr::GetGausWeight(double q0, double q3, double vals[]) { // // 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); double 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 || isnan(w) || w < 0.0){ w = 0.0; } if (w < 1.0 and !fAllowSuppression){ w = 1.0; } 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 15e2c6b..640842e 100644 --- a/src/Reweight/NUISANCEWeightCalcs.h +++ b/src/Reweight/NUISANCEWeightCalcs.h @@ -1,90 +1,108 @@ #ifndef NUISANCE_WEIGHT_CALCS #define NUISANCE_WEIGHT_CALCS #include "BaseFitEvt.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 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[]); // 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; }; #endif diff --git a/src/Reweight/NUISANCEWeightEngine.cxx b/src/Reweight/NUISANCEWeightEngine.cxx index 116385c..39eb1f9 100644 --- a/src/Reweight/NUISANCEWeightEngine.cxx +++ b/src/Reweight/NUISANCEWeightEngine.cxx @@ -1,110 +1,111 @@ #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 fWeightCalculators.push_back(new GaussianModeCorr()); fWeightCalculators.push_back(new ModeNormCalc()); + fWeightCalculators.push_back(new SBLOscWeightCalc()); #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 std::cout << "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; }