Page Menu
Home
HEPForge
Search
Configure Global Search
Log In
Files
F9501793
No One
Temporary
Actions
View File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Flag For Later
Size
88 KB
Subscribers
None
View Options
diff --git a/src/Reweight/MINERvAWeightCalcs.cxx b/src/Reweight/MINERvAWeightCalcs.cxx
index 4733a3e..b9228d6 100644
--- a/src/Reweight/MINERvAWeightCalcs.cxx
+++ b/src/Reweight/MINERvAWeightCalcs.cxx
@@ -1,632 +1,1109 @@
#ifdef __MINERVA_RW_ENABLED__
#ifdef __GENIE_ENABLED__
#include "MINERvAWeightCalcs.h"
#include "BaseFitEvt.h"
namespace nuisance {
namespace reweight {
//*******************************************************
MINERvAReWeight_QE::MINERvAReWeight_QE() {
//*******************************************************
fTwk_NormCCQE = 0.0;
fDef_NormCCQE = 1.0;
fCur_NormCCQE = fDef_NormCCQE;
}
//*******************************************************
MINERvAReWeight_QE::~MINERvAReWeight_QE(){};
//*******************************************************
//*******************************************************
double MINERvAReWeight_QE::CalcWeight(BaseFitEvt *evt) {
//*******************************************************
// Check GENIE
- if (evt->fType != kGENIE)
- return 1.0;
+ if (evt->fType != kGENIE) return 1.0;
// Extract the GENIE Record
GHepRecord *ghep = static_cast<GHepRecord *>(evt->genie_event->event);
const Interaction *interaction = ghep->Summary();
// const InitialState& init_state = interaction->InitState();
const ProcessInfo &proc_info = interaction->ProcInfo();
// const Target& tgt = init_state.Tgt();
// If the event is not QE this Calc doesn't handle it
- if (!proc_info.IsQuasiElastic())
- return 1.0;
+ if (!proc_info.IsQuasiElastic()) return 1.0;
// WEIGHT CALCULATIONS -------------
double w = 1.0;
// CCQE Dial
- if (!proc_info.IsWeakCC())
- w *= fCur_NormCCQE;
+ if (!proc_info.IsWeakCC()) w *= fCur_NormCCQE;
// Return Combined Weight
return w;
}
//*******************************************************
void MINERvAReWeight_QE::SetDialValue(std::string name, double val) {
//*******************************************************
SetDialValue(Reweight::ConvDial(name, kCUSTOM), val);
}
//*******************************************************
void MINERvAReWeight_QE::SetDialValue(int rwenum, double val) {
//*******************************************************
// Check Handled
int curenum = rwenum % 1000;
- if (!IsHandled(curenum))
- return;
+ if (!IsHandled(curenum)) return;
// Set Values
if (curenum == Reweight::kMINERvARW_NormCCQE) {
fTwk_NormCCQE = val;
fCur_NormCCQE = fDef_NormCCQE + fTwk_NormCCQE;
}
// Define Tweaked
fTweaked = ((fTwk_NormCCQE != 0.0));
}
//*******************************************************
bool MINERvAReWeight_QE::IsHandled(int rwenum) {
//*******************************************************
int curenum = rwenum % 1000;
switch (curenum) {
case Reweight::kMINERvARW_NormCCQE:
return true;
default:
return false;
}
}
//*******************************************************
MINERvAReWeight_MEC::MINERvAReWeight_MEC() {
//*******************************************************
fTwk_NormCCMEC = 0.0;
fDef_NormCCMEC = 1.0;
fCur_NormCCMEC = fDef_NormCCMEC;
}
//*******************************************************
MINERvAReWeight_MEC::~MINERvAReWeight_MEC(){};
//*******************************************************
//*******************************************************
double MINERvAReWeight_MEC::CalcWeight(BaseFitEvt *evt) {
//*******************************************************
// Check GENIE
if (evt->fType != kGENIE)
return 1.0;
// Extract the GENIE Record
GHepRecord *ghep = static_cast<GHepRecord *>(evt->genie_event->event);
const Interaction *interaction = ghep->Summary();
// const InitialState& init_state = interaction->InitState();
const ProcessInfo &proc_info = interaction->ProcInfo();
// const Target& tgt = init_state.Tgt();
// If the event is not MEC this Calc doesn't handle it
- if (!proc_info.IsMEC())
- return 1.0;
+ if (!proc_info.IsMEC()) return 1.0;
// WEIGHT CALCULATIONS -------------
double w = 1.0;
// CCMEC Dial
- if (!proc_info.IsWeakCC())
- w *= fCur_NormCCMEC;
+ if (!proc_info.IsWeakCC()) w *= fCur_NormCCMEC;
// Return Combined Weight
return w;
}
//*******************************************************
void MINERvAReWeight_MEC::SetDialValue(std::string name, double val) {
//*******************************************************
SetDialValue(Reweight::ConvDial(name, kCUSTOM), val);
}
//*******************************************************
void MINERvAReWeight_MEC::SetDialValue(int rwenum, double val) {
//*******************************************************
// Check Handled
int curenum = rwenum % 1000;
if (!IsHandled(curenum))
return;
// Set Values
if (curenum == Reweight::kMINERvARW_NormCCMEC) {
fTwk_NormCCMEC = val;
fCur_NormCCMEC = fDef_NormCCMEC + fTwk_NormCCMEC;
}
// Define Tweaked
fTweaked = ((fTwk_NormCCMEC != 0.0));
}
//*******************************************************
bool MINERvAReWeight_MEC::IsHandled(int rwenum) {
//*******************************************************
int curenum = rwenum % 1000;
switch (curenum) {
case Reweight::kMINERvARW_NormCCMEC:
return true;
default:
return false;
}
}
//*******************************************************
MINERvAReWeight_RES::MINERvAReWeight_RES() {
//*******************************************************
fTwk_NormCCRES = 0.0;
fDef_NormCCRES = 1.0;
fCur_NormCCRES = fDef_NormCCRES;
}
//*******************************************************
MINERvAReWeight_RES::~MINERvAReWeight_RES(){};
//*******************************************************
//*******************************************************
double MINERvAReWeight_RES::CalcWeight(BaseFitEvt *evt) {
//*******************************************************
// std::cout << "Caculating RES" << std::endl;
// Check GENIE
- if (evt->fType != kGENIE)
- return 1.0;
+ if (evt->fType != kGENIE) return 1.0;
// Extract the GENIE Record
GHepRecord *ghep = static_cast<GHepRecord *>(evt->genie_event->event);
const Interaction *interaction = ghep->Summary();
// const InitialState& init_state = interaction->InitState();
const ProcessInfo &proc_info = interaction->ProcInfo();
// const Target& tgt = init_state.Tgt();
// If the event is not RES this Calc doesn't handle it
- if (!proc_info.IsResonant())
- return 1.0;
+ if (!proc_info.IsResonant()) return 1.0;
// WEIGHT CALCULATIONS -------------
double w = 1.0;
// CCRES Dial
- if (proc_info.IsWeakCC())
- w *= fCur_NormCCRES;
+ if (proc_info.IsWeakCC()) w *= fCur_NormCCRES;
// Return Combined Weight
return w;
}
//*******************************************************
void MINERvAReWeight_RES::SetDialValue(std::string name, double val) {
//*******************************************************
SetDialValue(Reweight::ConvDial(name, kCUSTOM), val);
}
//*******************************************************
void MINERvAReWeight_RES::SetDialValue(int rwenum, double val) {
//*******************************************************
// Check Handled
int curenum = rwenum % 1000;
- if (!IsHandled(curenum))
- return;
+ if (!IsHandled(curenum)) return;
// Set Values
if (curenum == Reweight::kMINERvARW_NormCCRES) {
fTwk_NormCCRES = val;
fCur_NormCCRES = fDef_NormCCRES + fTwk_NormCCRES;
}
// Define Tweaked
fTweaked = ((fTwk_NormCCRES != 0.0));
}
//*******************************************************
bool MINERvAReWeight_RES::IsHandled(int rwenum) {
//*******************************************************
int curenum = rwenum % 1000;
switch (curenum) {
- case Reweight::kMINERvARW_NormCCRES:
- return true;
- default:
- return false;
+ case Reweight::kMINERvARW_NormCCRES:
+ return true;
+ default:
+ return false;
+ }
+}
+
+//*****************************************************************************
+MINOSRPA::MINOSRPA() {
+
+ fApply_MINOSRPA = false;
+
+ fDef_MINOSRPA_A = 1.010;
+ fTwk_MINOSRPA_A = 0.000;
+ fDef_MINOSRPA_B = 0.156;
+ fTwk_MINOSRPA_B = 0.000;
+
+}
+//
+double MINOSRPA::CalcWeight(BaseFitEvt* evt) {
+ if (!fTweaked) return 1.0;
+ double w = 1.0;
+
+ // Extract the GENIE Record
+ GHepRecord* ghep = static_cast<GHepRecord*>(evt->genie_event->event);
+ const Interaction* interaction = ghep->Summary();
+ const InitialState& init_state = interaction->InitState();
+ const ProcessInfo& proc_info = interaction->ProcInfo();
+ const Target& tgt = init_state.Tgt();
+
+ // If not QE return 1.0
+ if (!tgt.IsNucleus()) return 1.0;
+ if (!proc_info.IsQuasiElastic() && !proc_info.IsResonant()) return 1.0;
+
+ // Extract Beam and Target PDG
+ GHepParticle* neutrino = ghep->Probe();
+ // int bpdg = neutrino->Pdg();
+
+ GHepParticle* target = ghep->Particle(1);
+ assert(target);
+ // int tpdg = target->Pdg();
+
+ // Extract Q0-Q3
+ GHepParticle* fsl = ghep->FinalStatePrimaryLepton();
+ const TLorentzVector& k1 = *(neutrino->P4());
+ const TLorentzVector& k2 = *(fsl->P4());
+ //double q0 = fabs((k1 - k2).E());
+ //double q3 = fabs((k1 - k2).Vect().Mag());
+ double Q2 = fabs((k1 - k2).Mag2());
+
+ // Resonant Events
+ if (proc_info.IsResonant() && fApply_MINOSRPA) {
+ w *= GetRPAWeight(
+ Q2,
+ fCur_MINOSRPA_A,
+ fCur_MINOSRPA_B
+ );
+ }
+
+ return w;
+}
+//
+double MINOSRPA::GetRPAWeight(double Q2, double A, double B) {
+ if (Q2 > 0.7) return 1.0;
+ double w = A / (1.0 + TMath::Exp(1.0 - TMath::Sqrt(Q2) / B));
+ return w;
+}
+//
+bool MINOSRPA::IsHandled(int rwenum) {
+ int curenum = rwenum % 1000;
+ switch (curenum) {
+ case Reweight::kMINERvARW_MINOSRPA_Apply:
+ case Reweight::kMINERvARW_MINOSRPA_A:
+ case Reweight::kMINERvARW_MINOSRPA_B:
+ return true;
+ default:
+ return false;
+ }
+}
+//
+void MINOSRPA::SetDialValue(std::string name, double val) {
+ SetDialValue(Reweight::ConvDial(name, kCUSTOM), val);
+}
+//
+void MINOSRPA::SetDialValue(int rwenum, double val) {
+ int curenum = rwenum % 1000;
+
+ // Check Handled
+ if (!IsHandled(curenum)) return;
+ if (curenum == Reweight::kMINERvARW_MINOSRPA_Apply) fApply_MINOSRPA = (val > 0.5);
+ if (curenum == Reweight::kMINERvARW_MINOSRPA_A) fTwk_MINOSRPA_A = val;
+ if (curenum == Reweight::kMINERvARW_MINOSRPA_B) fTwk_MINOSRPA_B = val;
+
+ // Check for changes
+ fTweaked = (fApply_MINOSRPA ||
+ fabs(fTwk_MINOSRPA_A) > 0.0 ||
+ fabs(fTwk_MINOSRPA_B) > 0.0);
+
+ // Update Values
+ fCur_MINOSRPA_A = fDef_MINOSRPA_A * (1.0 + 0.1 * fTwk_MINOSRPA_A);
+ fCur_MINOSRPA_B = fDef_MINOSRPA_B * (1.0 + 0.1 * fTwk_MINOSRPA_B);
+
+}
+//
+//
+//*****************************************************************************
+LagrangeRPA::LagrangeRPA() {
+
+ fApplyRPA = false;
+
+ /*
+ fI1_Def = 4.18962e-01;
+ fI1 = fI1_Def;
+
+ fI2_Def = 7.39927e-01;
+ fI2 = fI2_Def;
+ */
+
+ // Table VIII https://arxiv.org/pdf/1903.01558.pdf
+ fR1_Def = 0.37;
+ fR1 = fR1_Def;
+
+ fR2_Def = 0.60;
+ fR2 = fR2_Def;
+
+}
+//
+double LagrangeRPA::CalcWeight(BaseFitEvt* evt) {
+
+ if (!fTweaked) return 1.0;
+ double w = 1.0;
+
+ // Extract the GENIE Record
+ GHepRecord* ghep = static_cast<GHepRecord*>(evt->genie_event->event);
+ const Interaction* interaction = ghep->Summary();
+ const InitialState& init_state = interaction->InitState();
+ const ProcessInfo& proc_info = interaction->ProcInfo();
+ const Target& tgt = init_state.Tgt();
+
+ // If not QE return 1.0
+ if (!tgt.IsNucleus()) return 1.0;
+ if (!proc_info.IsQuasiElastic() && !proc_info.IsResonant()) return 1.0;
+
+ // Extract Beam and Target PDG
+ GHepParticle* neutrino = ghep->Probe();
+ // int bpdg = neutrino->Pdg();
+
+ GHepParticle* target = ghep->Particle(1);
+ assert(target);
+ // int tpdg = target->Pdg();
+
+ // Extract Q0-Q3
+ GHepParticle* fsl = ghep->FinalStatePrimaryLepton();
+ const TLorentzVector& k1 = *(neutrino->P4());
+ const TLorentzVector& k2 = *(fsl->P4());
+ //double q0 = fabs((k1 - k2).E());
+ //double q3 = fabs((k1 - k2).Vect().Mag());
+ double Q2 = fabs((k1 - k2).Mag2());
+
+ if (proc_info.IsResonant() && fApplyRPA) {
+ w *= GetRPAWeight(Q2);
}
+
+ return w;
}
+//
+//
+double LagrangeRPA::GetRPAWeight(double Q2) {
+ //std::cout << "Getting RPA Weight : " << Q2 << std::endl;
+ if (Q2 > 0.7) return 1.0;
+
+ // Keep original Lagrange RPA for documentation
+ /*
+ double x1 = 0.00;
+ double x2 = 0.30;
+ double x3 = 0.70;
+
+ double y1 = 0.00;
+ double y2 = fI2;
+ double y3 = 1.00;
+
+ double xv = Q2;
+
+ // Automatically 0 because y1 choice
+ double W1 = y1 * (xv-x2)*(xv-x3)/((x1-x2)*(x1-x3));
+ double W2 = y2 * (xv-x1)*(xv-x3)/((x2-x1)*(x2-x3));
+ double W3 = y3 * (xv-x1)*(xv-x2)/((x3-x1)*(x3-x2));
+
+ double P = W1 + W2 + W3;
+ double A1 = (1.0 - sqrt(1.0 - fI1));
+ double R = P * (1.0 - A1) + A1;
+
+ return 1.0 - (1.0-R)*(1.0-R);
+ */
+
+ // Equation 7 https://arxiv.org/pdf/1903.01558.pdf
+ const double x1 = 0.00;
+ const double x2 = 0.35;
+ const double x3 = 0.70;
+
+ // Equation 6 https://arxiv.org/pdf/1903.01558.pdf
+ double RQ2 = fR2 *( (Q2-x1)*(Q2-x3)/((x2-x1)*(x2-x3)) )
+ + (Q2-x1)*(Q2-x2)/((x3-x1)*(x3-x2));
+ double weight = 1-(1-fR1)*(1-RQ2)*(1-RQ2);
+
+ // Check range of R1 and R2
+ // Commented out because this is implementation dependent: user may want strange double peaks
+ /*
+ if (fR1 > 1) return 1;
+ if (fR2 > 1 || fR2 < 0.5) return 1;
+ */
+
+ return weight;
+}
+//
+bool LagrangeRPA::IsHandled(int rwenum) {
+ int curenum = rwenum % 1000;
+ switch (curenum) {
+ case Reweight::kMINERvARW_LagrangeRPA_Apply:
+ case Reweight::kMINERvARW_LagrangeRPA_R1:
+ case Reweight::kMINERvARW_LagrangeRPA_R2:
+ return true;
+ default:
+ return false;
+ }
+}
+//
+void LagrangeRPA::SetDialValue(std::string name, double val) {
+ SetDialValue(Reweight::ConvDial(name, kCUSTOM), val);
+}
+//
+void LagrangeRPA::SetDialValue(int rwenum, double val) {
+ int curenum = rwenum % 1000;
+
+ // Check Handled
+ if (!IsHandled(curenum)) return;
+ if (curenum == Reweight::kMINERvARW_LagrangeRPA_Apply) fApplyRPA = (val > 0.5);
+ if (curenum == Reweight::kMINERvARW_LagrangeRPA_R1) fR1 = val;
+ if (curenum == Reweight::kMINERvARW_LagrangeRPA_R2) fR2 = val;
+
+ // Check for changes
+ fTweaked = (fApplyRPA);
+}
+//
//*******************************************************
RikRPA::RikRPA() {
//*******************************************************
// - Syst : kMINERvA_RikRPA_ApplyRPA
// - Type : Binary
// - Limits : 0.0 (false) -> 1.0 (true)
// - Default : 0.0
fApplyDial_RPACorrection = false;
// - Syst : kMINERvA_RikRPA_LowQ2
// - Type : Absolute
// - Limits : 1.0 -> 1.0
// - Default : 0.0
// - Frac Error : 100%
fDefDial_RPALowQ2 = 0.0;
fCurDial_RPALowQ2 = fDefDial_RPALowQ2;
fErrDial_RPALowQ2 = 0.0;
// - Syst : kMINERvA_RikRPA_HighQ2
// - Type : Absolute
// - Limits : 1.0 -> 1.0
// - Default : 0.0
// - Frac Error : 100%
fDefDial_RPAHighQ2 = 0.0;
fCurDial_RPAHighQ2 = fDefDial_RPAHighQ2;
fErrDial_RPAHighQ2 = 1.0;
// - Syst : kMINERvA_RikRESRPA_ApplyRPA
// - Type : Binary
// - Limits : 0.0 (false) -> 1.0 (true)
// - Default : 0.0
fApplyDial_RESRPACorrection = false;
// - Syst : kMINERvA_RikRESRPA_LowQ2
// - Type : Absolute
// - Limits : 1.0 -> 1.0
// - Default : 0.0
// - Frac Error : 100%
fDefDial_RESRPALowQ2 = 0.0;
fCurDial_RESRPALowQ2 = fDefDial_RESRPALowQ2;
fErrDial_RESRPALowQ2 = 0.0;
// - Syst : kMINERvA_RikRESRPA_HighQ2
// - Type : Absolute
// - Limits : 1.0 -> 1.0
// - Default : 0.0
// - Frac Error : 100%
fDefDial_RESRPAHighQ2 = 0.0;
fCurDial_RESRPAHighQ2 = fDefDial_RESRPAHighQ2;
fErrDial_RESRPAHighQ2 = 1.0;
// Setup Calculators
fEventWeights = new double[5];
for (int i = 0; i < kMaxCalculators; i++) {
fRPACalculators[i] = NULL;
}
fTweaked = false;
}
//*******************************************************
RikRPA::~RikRPA() {
//*******************************************************
// delete fEventWeights;
// for (size_t i = 0; i < kMaxCalculators; i++) {
// if (fRPACalculators[i]) delete fRPACalculators[i];
// fRPACalculators[i] = NULL;
// }
}
//*******************************************************
double RikRPA::CalcWeight(BaseFitEvt *evt) {
//*******************************************************
// LOG(FIT) << "Calculating RikRPA" << std::endl;
// Return 1.0 if not tweaked
if (!fTweaked)
return 1.0;
double w = 1.0;
// Extract the GENIE Record
GHepRecord *ghep = static_cast<GHepRecord *>(evt->genie_event->event);
const Interaction *interaction = ghep->Summary();
const InitialState &init_state = interaction->InitState();
const ProcessInfo &proc_info = interaction->ProcInfo();
// const Kinematics & kine = interaction->Kine();
// const XclsTag & xcls = interaction->ExclTag();
const Target &tgt = init_state.Tgt();
// If not QE return 1.0
// LOG(FIT) << "RikRPA : Event QE = " << proc_info.IsQuasiElastic() <<
// std::endl;
if (!tgt.IsNucleus()) {
return 1.0;
}
if (!proc_info.IsQuasiElastic() && !proc_info.IsResonant())
return 1.0;
// Extract Beam and Target PDG
GHepParticle *neutrino = ghep->Probe();
int bpdg = neutrino->Pdg();
GHepParticle *target = ghep->TargetNucleus();
assert(target);
int tpdg = target->Pdg();
// Find the enum we need
int calcenum = GetRPACalcEnum(bpdg, tpdg);
- if (calcenum == -1)
+ if (calcenum == -1)
return 1.0;
// Check we have the RPA Calc setup for this enum
// if not, set it up at that point
- if (!fRPACalculators[calcenum])
+ if (!fRPACalculators[calcenum])
SetupRPACalculator(calcenum);
weightRPA *rpacalc = fRPACalculators[calcenum];
if (!rpacalc) {
NUIS_ABORT("Failed to grab the RPA Calculator : " << calcenum);
}
// Extract Q0-Q3
GHepParticle *fsl = ghep->FinalStatePrimaryLepton();
const TLorentzVector &k1 = *(neutrino->P4());
const TLorentzVector &k2 = *(fsl->P4());
double q0 = fabs((k1 - k2).E());
double q3 = fabs((k1 - k2).Vect().Mag());
double Q2 = fabs((k1 - k2).Mag2());
// Quasielastic
if (proc_info.IsQuasiElastic()) {
// Now use q0-q3 and RPA Calculator to fill fWeights
rpacalc->getWeight(q0, q3, fEventWeights);
if (fApplyDial_RPACorrection) {
w *= fEventWeights[0]; // CV
}
// Syst Application : kMINERvA_RikRPA_LowQ2
if (fabs(fCurDial_RPALowQ2) > 0.0) {
double interpw = fEventWeights[0];
if (fCurDial_RPALowQ2 > 0.0 && Q2 < 2.0) {
interpw =
fEventWeights[0] - (fEventWeights[0] - fEventWeights[1]) *
- fCurDial_RPALowQ2; // WLow+ } else if
+ fCurDial_RPALowQ2; // WLow+
} else if (fCurDial_RPALowQ2 < 0.0 && Q2 < 2.0) {
interpw = fEventWeights[0] - (fEventWeights[2] - fEventWeights[0]) *
fCurDial_RPALowQ2; // WLow-
}
w *= interpw / fEventWeights[0]; // Div by CV again
}
// Syst Application : kMINERvA_RikRPA_HighQ2
if (fabs(fCurDial_RPAHighQ2) > 0.0) {
double interpw = fEventWeights[0];
if (fCurDial_RPAHighQ2 > 0.0) {
interpw = fEventWeights[0] - (fEventWeights[0] - fEventWeights[3]) *
- fCurDial_RPAHighQ2; // WHigh+
+ fCurDial_RPAHighQ2; // WHigh+
} else if (fCurDial_RPAHighQ2 < 0.0) {
interpw = fEventWeights[0] - (fEventWeights[4] - fEventWeights[0]) *
- fCurDial_RPAHighQ2; // WHigh-
+ fCurDial_RPAHighQ2; // WHigh-
}
w *= interpw / fEventWeights[0]; // Div by CV again
}
}
// Resonant Events
if (proc_info.IsResonant()) {
// Now use Q2 and RESRPA Calculator to fill fWeights
double CV = rpacalc->getWeight(Q2);
if (fApplyDial_RESRPACorrection) {
w *= CV; // fEventWeights[0]; // CVa
}
/*
// Syst Application : kMINERvA_RikRESRPA_LowQ2
if (fabs(fCurDial_RESRPAHighQ2) > 0.0) {
- double interpw = fEventWeights[0];
+ double interpw = fEventWeights[0];
- if (fCurDial_RESRPAHighQ2 > 0.0) {
- interpw = fEventWeights[0] - (fEventWeights[0] - fEventWeights[3]) *
- fCurDial_RESRPAHighQ2; // WHigh+
+ if (fCurDial_RESRPAHighQ2 > 0.0) {
+ interpw = fEventWeights[0] - (fEventWeights[0] - fEventWeights[3]) *
+ fCurDial_RESRPAHighQ2; // WHigh+
- } else if (fCurDial_RESRPAHighQ2 < 0.0) {
- interpw = fEventWeights[0] - (fEventWeights[4] - fEventWeights[0]) *
- fCurDial_RESRPAHighQ2; // WHigh-
- }
- w *= interpw / fEventWeights[0]; // Div by CV again
+ } else if (fCurDial_RESRPAHighQ2 < 0.0) {
+ interpw = fEventWeights[0] - (fEventWeights[4] - fEventWeights[0]) *
+ fCurDial_RESRPAHighQ2; // WHigh-
+ }
+ w *= interpw / fEventWeights[0]; // Div by CV again
}
// Syst Application : kMINERvA_RikRESRPA_HighQ2
if (fabs(fCurDial_RESRPAHighQ2) > 0.0) {
- double interpw = fEventWeights[0];
+ double interpw = fEventWeights[0];
- if (fCurDial_RESRPAHighQ2 > 0.0) {
- interpw = fEventWeights[0] - (fEventWeights[0] - fEventWeights[3]) *
- fCurDial_RESRPAHighQ2; // WHigh+
+ if (fCurDial_RESRPAHighQ2 > 0.0) {
+ interpw = fEventWeights[0] - (fEventWeights[0] - fEventWeights[3]) *
+ fCurDial_RESRPAHighQ2; // WHigh+
- } else if (fCurDial_RESRPAHighQ2 < 0.0) {
- interpw = fEventWeights[0] - (fEventWeights[4] - fEventWeights[0]) *
- fCurDial_RESRPAHighQ2; // WHigh-
- }
- w *= interpw / fEventWeights[0]; // Div by CV again
+ } else if (fCurDial_RESRPAHighQ2 < 0.0) {
+ interpw = fEventWeights[0] - (fEventWeights[4] - fEventWeights[0]) *
+ fCurDial_RESRPAHighQ2; // WHigh-
+ }
+ w *= interpw / fEventWeights[0]; // Div by CV again
}
*/
}
// LOG(FIT) << "RPA Weight = " << w << std::endl;
return w;
} // namespace reweight
//*******************************************************
void RikRPA::SetDialValue(std::string name, double val) {
//*******************************************************
SetDialValue(Reweight::ConvDial(name, kCUSTOM), val);
}
//*******************************************************
void RikRPA::SetDialValue(int rwenum, double val) {
//*******************************************************
int curenum = rwenum % 1000;
// Check Handled
if (!IsHandled(curenum))
return;
if (curenum == Reweight::kMINERvARW_RikRPA_ApplyRPA)
fApplyDial_RPACorrection = (val > 0.5);
if (curenum == Reweight::kMINERvARW_RikRPA_LowQ2)
fCurDial_RPALowQ2 = val;
if (curenum == Reweight::kMINERvARW_RikRPA_HighQ2)
fCurDial_RPAHighQ2 = val;
if (curenum == Reweight::kMINERvARW_RikRESRPA_ApplyRPA)
fApplyDial_RESRPACorrection = (val > 0.5);
if (curenum == Reweight::kMINERvARW_RikRESRPA_LowQ2)
fCurDial_RESRPALowQ2 = val;
if (curenum == Reweight::kMINERvARW_RikRESRPA_HighQ2)
fCurDial_RESRPAHighQ2 = val;
// Assign flag to say stuff has changed
fTweaked = (fApplyDial_RPACorrection ||
- fabs(fCurDial_RPAHighQ2 - fDefDial_RPAHighQ2) > 0.0 ||
- fabs(fCurDial_RPALowQ2 - fDefDial_RPALowQ2) > 0.0 ||
- fApplyDial_RESRPACorrection ||
- fabs(fCurDial_RESRPAHighQ2 - fDefDial_RESRPAHighQ2) > 0.0 ||
- fabs(fCurDial_RESRPALowQ2 - fDefDial_RESRPALowQ2) > 0.0);
+ fabs(fCurDial_RPAHighQ2 - fDefDial_RPAHighQ2) > 0.0 ||
+ fabs(fCurDial_RPALowQ2 - fDefDial_RPALowQ2) > 0.0 ||
+ fApplyDial_RESRPACorrection ||
+ fabs(fCurDial_RESRPAHighQ2 - fDefDial_RESRPAHighQ2) > 0.0 ||
+ fabs(fCurDial_RESRPALowQ2 - fDefDial_RESRPALowQ2) > 0.0);
}
//*******************************************************
bool RikRPA::IsHandled(int rwenum) {
//*******************************************************
int curenum = rwenum % 1000;
switch (curenum) {
- case Reweight::kMINERvARW_RikRESRPA_ApplyRPA:
- return true;
- case Reweight::kMINERvARW_RikRESRPA_LowQ2:
- return true;
- case Reweight::kMINERvARW_RikRESRPA_HighQ2:
- return true;
- case Reweight::kMINERvARW_RikRPA_ApplyRPA:
- return true;
- case Reweight::kMINERvARW_RikRPA_LowQ2:
- return true;
- case Reweight::kMINERvARW_RikRPA_HighQ2:
- return true;
- default:
- return false;
+ case Reweight::kMINERvARW_RikRESRPA_ApplyRPA:
+ return true;
+ case Reweight::kMINERvARW_RikRESRPA_LowQ2:
+ return true;
+ case Reweight::kMINERvARW_RikRESRPA_HighQ2:
+ return true;
+ case Reweight::kMINERvARW_RikRPA_ApplyRPA:
+ return true;
+ case Reweight::kMINERvARW_RikRPA_LowQ2:
+ return true;
+ case Reweight::kMINERvARW_RikRPA_HighQ2:
+ return true;
+ default:
+ return false;
}
}
//*******************************************************
void RikRPA::SetupRPACalculator(int calcenum) {
//*******************************************************
std::string rwdir = FitPar::GetDataBase() + "reweight/MINERvA/RikRPA/";
std::string fidir = "";
switch (calcenum) {
- case kNuMuC12:
- fidir = "outNievesRPAratio-nu12C-20GeV-20170202.root";
- break;
+ case kNuMuC12:
+ fidir = "outNievesRPAratio-nu12C-20GeV-20170202.root";
+ break;
- case kNuMuO16:
- fidir = "outNievesRPAratio-nu16O-20GeV-20170202.root";
- break;
+ case kNuMuO16:
+ fidir = "outNievesRPAratio-nu16O-20GeV-20170202.root";
+ break;
- case kNuMuAr40:
- fidir = "outNievesRPAratio-nu40Ar-20GeV-20170202.root";
- break;
+ case kNuMuAr40:
+ fidir = "outNievesRPAratio-nu40Ar-20GeV-20170202.root";
+ break;
- case kNuMuCa40:
- fidir = "outNievesRPAratio-nu40Ca-20GeV-20170202.root";
- break;
+ case kNuMuCa40:
+ fidir = "outNievesRPAratio-nu40Ca-20GeV-20170202.root";
+ break;
- case kNuMuFe56:
- fidir = "outNievesRPAratio-nu56Fe-20GeV-20170202.root";
- break;
+ case kNuMuFe56:
+ fidir = "outNievesRPAratio-nu56Fe-20GeV-20170202.root";
+ break;
- case kNuMuBarC12:
- fidir = "outNievesRPAratio-anu12C-20GeV-20170202.root";
- break;
+ case kNuMuBarC12:
+ fidir = "outNievesRPAratio-anu12C-20GeV-20170202.root";
+ break;
- case kNuMuBarO16:
- fidir = "outNievesRPAratio-anu16O-20GeV-20170202.root";
- break;
+ case kNuMuBarO16:
+ fidir = "outNievesRPAratio-anu16O-20GeV-20170202.root";
+ break;
- case kNuMuBarAr40:
- fidir = "outNievesRPAratio-anu40Ar-20GeV-20170202.root";
- break;
+ case kNuMuBarAr40:
+ fidir = "outNievesRPAratio-anu40Ar-20GeV-20170202.root";
+ break;
- case kNuMuBarCa40:
- fidir = "outNievesRPAratio-anu40Ca-20GeV-20170202.root";
- break;
+ case kNuMuBarCa40:
+ fidir = "outNievesRPAratio-anu40Ca-20GeV-20170202.root";
+ break;
- case kNuMuBarFe56:
- fidir = "outNievesRPAratio-anu56Fe-20GeV-20170202.root";
- break;
+ case kNuMuBarFe56:
+ fidir = "outNievesRPAratio-anu56Fe-20GeV-20170202.root";
+ break;
}
NUIS_LOG(FIT, "Loading RPA CALC : " << fidir);
TDirectory *olddir = gDirectory;
NUIS_LOG(FIT, "***********************************************");
NUIS_LOG(FIT, "Loading a new weightRPA calculator");
NUIS_LOG(FIT, "Authors: Rik Gran, Heidi Schellman");
NUIS_LOG(FIT, "Citation: arXiv:1705.02932 [hep-ex]");
NUIS_LOG(FIT, "***********************************************");
// Test the file exists
std::ifstream infile((rwdir + fidir).c_str());
if (!infile.good()) {
NUIS_ERR(FTL, "*** ERROR ***");
NUIS_ERR(FTL, "RikRPA file " << rwdir + fidir << " does not exist!");
NUIS_ERR(FTL,
- "These can be found at https://nuisance.hepforge.org/files/RikRPA/");
+ "These can be found at https://nuisance.hepforge.org/files/RikRPA/");
NUIS_ERR(FTL,
- "Please run: wget -r -nH --cut-dirs=2 -np -e robots=off -R "
- "\"index.html*\" https://nuisance.hepforge.org/files/RikRPA/ -P "
- << rwdir);
+ "Please run: wget -r -nH --cut-dirs=2 -np -e robots=off -R "
+ "\"index.html*\" https://nuisance.hepforge.org/files/RikRPA/ -P "
+ << rwdir);
NUIS_ERR(FTL, "And try again");
NUIS_ABORT("*************");
}
fRPACalculators[calcenum] = new weightRPA(rwdir + fidir);
olddir->cd();
return;
}
//*******************************************************
int RikRPA::GetRPACalcEnum(int bpdg, int tpdg) {
//*******************************************************
if (bpdg == 14 && tpdg == 1000060120)
return kNuMuC12;
else if (bpdg == 14 && tpdg == 1000080160)
return kNuMuO16;
else if (bpdg == 14 && tpdg == 1000180400)
return kNuMuAr40;
else if (bpdg == 14 && tpdg == 1000200400)
return kNuMuCa40;
else if (bpdg == 14 && tpdg == 1000280560)
return kNuMuFe56;
else if (bpdg == -14 && tpdg == 1000060120)
return kNuMuBarC12;
else if (bpdg == -14 && tpdg == 1000080160)
return kNuMuBarO16;
else if (bpdg == -14 && tpdg == 1000180400)
return kNuMuBarAr40;
else if (bpdg == -14 && tpdg == 1000200400)
return kNuMuBarCa40;
else if (bpdg == -14 && tpdg == 1000280560)
return kNuMuBarFe56;
else {
// NUIS_ERR(WRN, "Unknown beam and target combination for RPA Calcs! "
//<< bpdg << " " << tpdg);
}
return -1;
}
+
+//*****************************************************************************
+COHBrandon::COHBrandon() {
+
+ fApply_COHNorm = false;
+
+ fDef_COHNorm = 0.5;
+ fCur_COHNorm = fDef_COHNorm;
+ fTwk_COHNorm = 0.0;
+
+ fDef_COHCut = 0.450;
+ fCur_COHCut = fDef_COHCut;
+ fTwk_COHNorm = 0.0;
+
+}
+
+COHBrandon::~COHBrandon() {};
+
+double COHBrandon::CalcWeight(BaseFitEvt* evt) {
+
+ // Check GENIE
+ if (evt->fType != kGENIE) return 1.0;
+
+ // Extract the GENIE Record
+ GHepRecord* ghep = static_cast<GHepRecord*>(evt->genie_event->event);
+ const Interaction* interaction = ghep->Summary();
+ //const InitialState& init_state = interaction->InitState();
+ const ProcessInfo& proc_info = interaction->ProcInfo();
+ //const Target& tgt = init_state.Tgt();
+
+ // If the event is not QE this Calc doesn't handle it
+ if (!proc_info.IsCoherent()) return 1.0;
+
+ // WEIGHT CALCULATIONS -------------
+ double w = 1.0;
+ double pionE = -999.9;
+
+ TObjArrayIter piter(ghep);
+ GHepParticle * p = 0;
+ //int ip = -1;
+ while ( (p = (GHepParticle *) piter.Next()) ) {
+ int pdgc = p->Pdg();
+ int ist = p->Status();
+ // if (ghep->Particle(ip)->FirstMother()==0) continue;
+ if (pdg::IsPseudoParticle(p->Pdg())) continue;
+ if (ist == kIStStableFinalState) {
+ if (pdgc == 211 || pdgc == -211 || pdgc == 111) {
+ pionE = p->Energy();
+ break;
+ }
+ }
+ }
+ // std::cout << "Coherent Pion Energy : " << pionE << std::endl;
+
+ // Apply Weight
+ if (fApply_COHNorm && pionE > 0.0 && pionE < fCur_COHCut) {
+ w *= fCur_COHNorm;
+ }
+
+ // Return Combined Weight
+ return w;
+}
+
+void COHBrandon::SetDialValue(std::string name, double val) {
+ SetDialValue(Reweight::ConvDial(name, kCUSTOM), val);
+}
+
+void COHBrandon::SetDialValue(int rwenum, double val) {
+ // Check Handled
+ int curenum = rwenum % 1000;
+ if (!IsHandled(curenum)) return;
+
+ // Set Values
+ if (curenum == Reweight::kMINERvARW_NormCOH) {
+ fTwk_COHNorm = val;
+ fCur_COHNorm = fDef_COHNorm * (1.0 + 0.1 * fTwk_COHNorm);
+ }
+
+ if (curenum == Reweight::kMINERvARW_CutCOH) {
+ fTwk_COHCut = val;
+ fCur_COHCut = fDef_COHCut * (1.0 + 0.1 * fTwk_COHCut);
+ }
+
+ if (curenum == Reweight::kMINERvARW_ApplyCOH) {
+ fApply_COHNorm = (val > 0.5);
+ }
+
+ // Define Tweaked
+ fTweaked = (fApply_COHNorm);
+}
+
+bool COHBrandon::IsHandled(int rwenum) {
+ int curenum = rwenum % 1000;
+
+ switch (curenum) {
+ case Reweight::kMINERvARW_NormCOH:
+ case Reweight::kMINERvARW_CutCOH:
+ case Reweight::kMINERvARW_ApplyCOH:
+ return true;
+ default:
+ return false;
+ }
+}
+
+//*****************************************************************************
+WEnhancement::WEnhancement() {
+
+ fApply_Enhancement = false;
+
+ fDef_WNorm = 1.0;
+ fCur_WNorm = fDef_WNorm;
+ fTwk_WNorm = 0.0;
+
+ fDef_WMean = 0.95;
+ fCur_WMean = fDef_WMean;
+ fTwk_WMean = 0.0;
+
+ fDef_WSigma = 0.225;
+ fCur_WSigma = fDef_WSigma;
+ fTwk_WSigma = 0.0;
+
+}
+
+WEnhancement::~WEnhancement() {};
+
+double WEnhancement::CalcWeight(BaseFitEvt* evt) {
+
+ // Check GENIE
+ if (evt->fType != kGENIE) return 1.0;
+
+ // Extract the GENIE Record
+ GHepRecord* ghep = static_cast<GHepRecord*>(evt->genie_event->event);
+ const Interaction* interaction = ghep->Summary();
+ //const Kinematics & kine = interaction->Kine();
+ //const InitialState& init_state = interaction->InitState();
+ const ProcessInfo& proc_info = interaction->ProcInfo();
+ //const Target& tgt = init_state.Tgt();
+
+ // If the event is not QE this Calc doesn't handle it
+ if (!proc_info.IsWeakCC()) return 1.0;
+ if (!proc_info.IsResonant()) return 1.0;
+ if (!fApply_Enhancement) return 1.0;
+
+ // WEIGHT CALCULATIONS -------------
+ double w = 1.0;
+ bool isCC1pi0AtVertex = false;
+
+ // Get W
+ GHepParticle* neutrino = ghep->Probe();
+ GHepParticle* fsl = ghep->FinalStatePrimaryLepton();
+ const TLorentzVector& k1 = *(neutrino->P4());
+ const TLorentzVector& k2 = *(fsl->P4());
+
+ double E_mu = k2.E();
+ double p_mu = k2.Vect().Mag();
+ double m_mu = sqrt(E_mu * E_mu - p_mu * p_mu);
+ double th_nu_mu = k1.Vect().Angle(k2.Vect());
+
+ // The factor of 1000 is necessary for downstream functions
+ const double m_p = PhysConst::mass_proton;
+ double E_nu = k1.E();
+ // double hadMass = kine.W (true);
+ double hadMass = 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)));
+
+
+ // Determine if event is CC1pi0 at vertex
+ TObjArrayIter piter(ghep);
+ GHepParticle * p = 0;
+ int pi0 = 0;
+ int piother = 0;
+ while ( (p = (GHepParticle *) piter.Next()) ) {
+ int pdgc = p->Pdg();
+ int ist = p->Status();
+ if (pdg::IsPseudoParticle(p->Pdg())) continue;
+ if (ist == genie::kIStStableFinalState) {
+ if (pdgc == 111) {
+ pi0++;
+ } else if (pdgc == 211 || pdgc == -211) {
+ piother++;
+ }
+ }
+ }
+
+ // std::cout << "Npi0 = " << pi0 << std::endl;
+ isCC1pi0AtVertex = (pi0 == 1 && piother == 0);
+
+ // Apply Weight
+ if (fApply_Enhancement && isCC1pi0AtVertex) {
+
+ double enhancement = 1.0 + fCur_WNorm * (exp( -0.5 * pow((fCur_WMean - hadMass) / (fCur_WSigma), 2) ));
+ w *= enhancement;
+ }
+
+ // Return Combined Weight
+ if (isnan(w) || w < 0.0 || w > 400.0) {
+ w = 1.0;
+ }
+
+ return w;
+}
+
+void WEnhancement::SetDialValue(std::string name, double val) {
+ SetDialValue(Reweight::ConvDial(name, kCUSTOM), val);
+}
+
+void WEnhancement::SetDialValue(int rwenum, double val) {
+ // Check Handled
+ int curenum = rwenum % 1000;
+ if (!IsHandled(curenum)) return;
+
+ // Set Values
+ if (curenum == Reweight::kMINERvARW_ApplyWTune) {
+ fApply_Enhancement = (val > 0.5);
+ }
+
+ if (curenum == Reweight::kMINERvARW_NormWTune) {
+ fTwk_WNorm = val;
+ fCur_WNorm = fDef_WNorm * (1.0 + 0.1 * fTwk_WNorm);
+ }
+
+ if (curenum == Reweight::kMINERvARW_MeanWTune) {
+ fTwk_WMean = val;
+ fCur_WMean = fDef_WMean * (1.0 + 0.1 * fTwk_WMean);
+ }
+
+ if (curenum == Reweight::kMINERvARW_SigmaWTune) {
+ fTwk_WSigma = val;
+ fCur_WSigma = fDef_WSigma * (1.0 + 0.1 * fTwk_WSigma);
+ }
+
+ // Define Tweaked
+ fTweaked = (fApply_Enhancement);
+}
+
+bool WEnhancement::IsHandled(int rwenum) {
+ int curenum = rwenum % 1000;
+
+ switch (curenum) {
+ case Reweight::kMINERvARW_ApplyWTune:
+ case Reweight::kMINERvARW_NormWTune:
+ case Reweight::kMINERvARW_MeanWTune:
+ case Reweight::kMINERvARW_SigmaWTune:
+ return true;
+ default:
+ return false;
+ }
+}
+
} // namespace reweight
} // namespace nuisance
#endif
#endif
diff --git a/src/Reweight/MINERvAWeightCalcs.h b/src/Reweight/MINERvAWeightCalcs.h
index 965010d..0daef1e 100644
--- a/src/Reweight/MINERvAWeightCalcs.h
+++ b/src/Reweight/MINERvAWeightCalcs.h
@@ -1,158 +1,259 @@
#ifndef MINERVA_WEIGHT_CALCS
#define MINERVA_WEIGHT_CALCS
#include <string>
#ifdef __MINERVA_RW_ENABLED__
#ifdef __GENIE_ENABLED__
#ifdef GENIE_PRE_R3
#include "Conventions/Units.h"
#include "EVGCore/EventRecord.h"
#include "GHEP/GHepParticle.h"
#include "GHEP/GHepRecord.h"
#include "GHEP/GHepUtils.h"
#include "Ntuple/NtpMCEventRecord.h"
#include "PDG/PDGUtils.h"
#else
#include "Framework/Conventions/Units.h"
#include "Framework/EventGen/EventRecord.h"
#include "Framework/GHEP/GHepParticle.h"
#include "Framework/GHEP/GHepRecord.h"
#include "Framework/GHEP/GHepUtils.h"
#include "Framework/Ntuple/NtpMCEventRecord.h"
#include "Framework/ParticleData/PDGUtils.h"
#endif
#include "NUISANCEWeightCalcs.h"
#include "GeneralUtils.h"
#include "NUISANCESyst.h"
#include "FitEvent.h"
#include "WeightUtils.h"
#include "weightRPA.h"
using namespace genie;
class BaseFitEvt;
namespace nuisance {
namespace reweight {
- // MEC Dials
+ // QE Dials
class MINERvAReWeight_QE : public NUISANCEWeightCalc {
public:
MINERvAReWeight_QE();
virtual ~MINERvAReWeight_QE();
double CalcWeight(BaseFitEvt* evt);
void SetDialValue(std::string name, double val);
void SetDialValue(int rwenum, double val);
bool IsHandled(int rwenum);
double fTwk_NormCCQE;
double fCur_NormCCQE;
double fDef_NormCCQE;
bool fTweaked;
};
// MEC Dials
class MINERvAReWeight_MEC : public NUISANCEWeightCalc {
public:
MINERvAReWeight_MEC();
virtual ~MINERvAReWeight_MEC();
double CalcWeight(BaseFitEvt* evt);
void SetDialValue(std::string name, double val);
void SetDialValue(int rwenum, double val);
bool IsHandled(int rwenum);
double fTwk_NormCCMEC;
double fCur_NormCCMEC;
double fDef_NormCCMEC;
bool fTweaked;
};
// RES Dials
class MINERvAReWeight_RES : public NUISANCEWeightCalc {
public:
MINERvAReWeight_RES();
virtual ~MINERvAReWeight_RES();
double CalcWeight(BaseFitEvt* evt);
void SetDialValue(std::string name, double val);
void SetDialValue(int rwenum, double val);
bool IsHandled(int rwenum);
double fTwk_NormCCRES;
double fCur_NormCCRES;
double fDef_NormCCRES;
bool fTweaked;
};
+ // MINOS pion tuning, https://arxiv.org/pdf/1903.01558.pdf
+ class MINOSRPA : public NUISANCEWeightCalc {
+ public:
+ MINOSRPA();
+ ~MINOSRPA(){};
+
+ double CalcWeight(BaseFitEvt* evt);
+ void SetDialValue(std::string name, double val);
+ void SetDialValue(int rwenum, double val);
+ bool IsHandled(int rwenum);
+
+ double GetRPAWeight(double Q2, double A, double B);
+
+ bool fTweaked;
+
+ bool fApply_MINOSRPA;
+
+ double fTwk_MINOSRPA_A;
+ double fDef_MINOSRPA_A;
+ double fCur_MINOSRPA_A;
+
+ double fTwk_MINOSRPA_B;
+ double fDef_MINOSRPA_B;
+ double fCur_MINOSRPA_B;
+
+ };
+
+ // MINERvA pion tuning, https://arxiv.org/pdf/1903.01558.pdf
+ class LagrangeRPA : public NUISANCEWeightCalc {
+ public:
+ LagrangeRPA();
+ ~LagrangeRPA(){};
+
+ double CalcWeight(BaseFitEvt* evt);
+ void SetDialValue(std::string name, double val);
+ void SetDialValue(int rwenum, double val);
+ bool IsHandled(int rwenum);
+
+ double GetRPAWeight(double Q2);
+
+ bool fTweaked;
+
+ bool fApplyRPA;
+
+ double fR1;
+ double fR2;
+ double fR1_Def;
+ double fR2_Def;
+ };
+
/// RPA Weight Calculator that applies RPA systematics
/// to GENIE events. GENIE EVENTS ONLY!
class RikRPA : public NUISANCEWeightCalc {
public:
RikRPA();
~RikRPA();
double CalcWeight(BaseFitEvt* evt);
void SetDialValue(std::string name, double val);
void SetDialValue(int rwenum, double val);
bool IsHandled(int rwenum);
void SetupRPACalculator(int calcenum);
int GetRPACalcEnum(int bpdg, int tpdg);
bool fApplyDial_RPACorrection;
double fTwkDial_RPALowQ2;
double fDefDial_RPALowQ2;
double fCurDial_RPALowQ2;
double fErrDial_RPALowQ2;
double fTwkDial_RPAHighQ2;
double fDefDial_RPAHighQ2;
double fCurDial_RPAHighQ2;
double fErrDial_RPAHighQ2;
bool fApplyDial_RESRPACorrection;
double fTwkDial_RESRPALowQ2;
double fDefDial_RESRPALowQ2;
double fCurDial_RESRPALowQ2;
double fErrDial_RESRPALowQ2;
double fTwkDial_RESRPAHighQ2;
double fDefDial_RESRPAHighQ2;
double fCurDial_RESRPAHighQ2;
double fErrDial_RESRPAHighQ2;
double* fEventWeights;
bool fTweaked;
const static int kMaxCalculators = 10;
enum rpacalcenums {
kNuMuC12,
kNuMuO16,
kNuMuAr40,
kNuMuCa40,
kNuMuFe56,
kNuMuBarC12,
kNuMuBarO16,
kNuMuBarAr40,
kNuMuBarCa40,
kNuMuBarFe56
};
weightRPA* fRPACalculators[kMaxCalculators];
};
+ // Custom coherent tune from MINERvA
+ class COHBrandon : public NUISANCEWeightCalc {
+ public:
+ COHBrandon();
+ ~COHBrandon();
+
+ double CalcWeight(BaseFitEvt* evt);
+ void SetDialValue(std::string name, double val);
+ void SetDialValue(int rwenum, double val);
+ bool IsHandled(int rwenum);
+
+ bool fApply_COHNorm;
+
+ double fDef_COHNorm;
+ double fCur_COHNorm;
+ double fTwk_COHNorm;
+
+ double fDef_COHCut;
+ double fCur_COHCut;
+ double fTwk_COHCut;
+
+ bool fTweaked;
+ };
+
+ class WEnhancement : public NUISANCEWeightCalc {
+ public:
+ WEnhancement();
+ ~WEnhancement();
+
+ double CalcWeight(BaseFitEvt* evt);
+ void SetDialValue(std::string name, double val);
+ void SetDialValue(int rwenum, double val);
+ bool IsHandled(int rwenum);
+
+ bool fTweaked;
+
+ bool fApply_Enhancement;
+
+ double fDef_WNorm;
+ double fCur_WNorm;
+ double fTwk_WNorm;
+
+ double fDef_WMean;
+ double fCur_WMean;
+ double fTwk_WMean;
+
+ double fDef_WSigma;
+ double fCur_WSigma;
+ double fTwk_WSigma;
+ };
+
}; // namespace reweight
}; // namespace nuisance
#endif // __GENIE_ENABLED__
#endif //__MINERVA_RW_ENABLED__
#endif
diff --git a/src/Reweight/NUISANCESyst.cxx b/src/Reweight/NUISANCESyst.cxx
index d11114b..def0a68 100644
--- a/src/Reweight/NUISANCESyst.cxx
+++ b/src/Reweight/NUISANCESyst.cxx
@@ -1,80 +1,96 @@
#include "NUISANCESyst.h"
int Reweight::ConvertNUISANCEDial(std::string type) {
- for (int i = kUnkownNUISANCEDial + 1; i < kNUISANCEDial_LAST; i++) {
+ for (int i = kUnknownNUISANCEDial + 1; i < kNUISANCEDial_LAST; i++) {
if (!type.compare(ConvNUISANCEDial(i).c_str())) {
return i;
}
}
- return kUnkownNUISANCEDial;
+ return kUnknownNUISANCEDial;
};
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"; }
+ case kMINERvARW_MINOSRPA_Apply: { return "MINERvARW_MINOSRPA_Apply"; }
+ case kMINERvARW_MINOSRPA_A: { return "MINERvARW_MINOSRPA_A"; }
+ case kMINERvARW_MINOSRPA_B: { return "MINERvARW_MINOSRPA_B"; }
+
+ case kMINERvARW_LagrangeRPA_Apply: { return "MINERvARW_LagrangeRPA_Apply"; }
+ case kMINERvARW_LagrangeRPA_R1: { return "MINERvARW_LagrangeRPA_R1"; }
+ case kMINERvARW_LagrangeRPA_R2: { return "MINERvARW_LagrangeRPA_R2"; }
+
+ case kMINERvARW_NormCOH: { return "MINERvARW_NormCOH"; }
+ case kMINERvARW_CutCOH: { return "MINERvARW_CutCOH"; }
+ case kMINERvARW_ApplyCOH: { return "MINERvARW_ApplyCOH"; }
+
+ case kMINERvARW_ApplyWTune: { return "MINERvARW_ApplyWTune"; }
+ case kMINERvARW_NormWTune: { return "MINERvARW_NormWTune"; }
+ case kMINERvARW_MeanWTune: { return "MINERvARW_MeanWTune"; }
+ case kMINERvARW_SigmaWTune: { return "MINERvARW_SigmaWTune"; }
+
default: return "unknown_nuisance_dial";
}
};
diff --git a/src/Reweight/NUISANCESyst.h b/src/Reweight/NUISANCESyst.h
index 14825ef..94d629e 100644
--- a/src/Reweight/NUISANCESyst.h
+++ b/src/Reweight/NUISANCESyst.h
@@ -1,78 +1,103 @@
#ifndef NUISANCESyst_H
#define NUISANCESyst_H
#include "GeneralUtils.h"
namespace Reweight {
enum NUISANCESyst {
- kUnkownNUISANCEDial = 0,
+ kUnknownNUISANCEDial = 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_CCRES_norm,
+ kGaussianCorr_CCRES_tilt,
+ kGaussianCorr_CCRES_Pq0,
+ kGaussianCorr_CCRES_Wq0,
+ kGaussianCorr_CCRES_Pq3,
+ kGaussianCorr_CCRES_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,
+ kMINERvARW_MINOSRPA_Apply,
+ kMINERvARW_MINOSRPA_A,
+ kMINERvARW_MINOSRPA_B,
+
+ kMINERvARW_LagrangeRPA_Apply,
+ kMINERvARW_LagrangeRPA_R1,
+ kMINERvARW_LagrangeRPA_R2,
+
+ kMINERvARW_NormCOH,
+ kMINERvARW_CutCOH,
+ kMINERvARW_ApplyCOH,
+
+ kMINERvARW_ApplyWTune,
+ kMINERvARW_NormWTune,
+ kMINERvARW_MeanWTune,
+ kMINERvARW_SigmaWTune,
+
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 bdf0a8b..9ae3551 100644
--- a/src/Reweight/NUISANCEWeightCalcs.cxx
+++ b/src/Reweight/NUISANCEWeightCalcs.cxx
@@ -1,538 +1,545 @@
#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))
+ if (!IsHandled(curenum))
return;
- if (curenum == kModeNorm_NormRES)
+ 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<FitEvent *>(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<FitEvent *>(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<FitEvent *>(evt);
double rw_weight = 1.0;
// Get Neutrino
if (!fevt->Npart()) {
NUIS_ABORT("NO particles found in stack!");
}
FitParticle *pnu = fevt->GetHMISAnyLeptons();
if (!pnu) {
NUIS_ABORT("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)
+ 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)
+ if (fevt->PartInfo(j)->fPID == 2212)
npr++;
- else if (fevt->PartInfo(j)->fPID == 2112)
+ 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))) {
+ } 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)
+ 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;
+ 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 (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;
+ 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;
+ 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) {
NUIS_LOG(FIT,
" Using tilt-shift Gaussian parameters for Gaussian enhancement...");
} else {
NUIS_LOG(FIT, " Using Normal Gaussian parameters for Gaussian enhancement...");
}
};
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;
}
+ return w;
+
// 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 z = (q0 - meanq0) * (q0 - meanq0) / (sigmaq0 * sigmaq0) +
+ (q3 - meanq3) * (q3 - meanq3) / (sigmaq3 * sigmaq3) -
+ 2 * corr * (q0 - meanq0) * (q3 - meanq3) / (sigmaq0 * sigmaq3);
+ double ret = 1.0;
+
+ if ( fabs(1 - corr*corr) < 1.E-5 ) {
+ return 1.0;
+ }
+
+ if ( (-0.5 * z / (1 - corr*corr)) > 200 or (-0.5 * z / (1 - corr*corr)) < -200 ) {
+ return 1.0;
+ } else {
+ ret = norm * exp( -0.5 * z / (1 - corr*corr) );
+ }
- double ret = norm * exp(-0.5 * z / (1 - corr * corr));
+ if (ret != ret or ret < 0.0 or isnan(ret)) {
+ return 1.0;
+ }
+
+ if (fAllowSuppression) return ret;
+ return ret + 1.0;
// 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;
+ 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++) {
+ 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++) {
+ 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;
+ 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/NUISANCEWeightEngine.cxx b/src/Reweight/NUISANCEWeightEngine.cxx
index 8bbd716..546cf7d 100644
--- a/src/Reweight/NUISANCEWeightEngine.cxx
+++ b/src/Reweight/NUISANCEWeightEngine.cxx
@@ -1,123 +1,128 @@
#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;
NUIS_LOG(FIT, "Setting up NUISANCE Custom RW : " << fCalcName);
// 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 {
NUIS_ABORT("I do not recognise method "
<< Gaussian_Method
<< " for the Gaussian enhancement, so will die now...");
}
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_QE());
fWeightCalculators.push_back(new nuisance::reweight::MINERvAReWeight_MEC());
fWeightCalculators.push_back(new nuisance::reweight::MINERvAReWeight_RES());
+ fWeightCalculators.push_back(new nuisance::reweight::MINOSRPA());
+ fWeightCalculators.push_back(new nuisance::reweight::LagrangeRPA());
fWeightCalculators.push_back(new nuisance::reweight::RikRPA());
+ fWeightCalculators.push_back(new nuisance::reweight::COHBrandon());
+ fWeightCalculators.push_back(new nuisance::reweight::WEnhancement());
#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<size_t>(0);
fNameIndex[name]; // = std::vector<size_t>(0);
// Split by commas
std::vector<std::string> 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
NUIS_LOG(FIT, "Registering " << singlename << " from " << name);
// 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<size_t> 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<size_t> 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<NUISANCEWeightCalc *>::iterator calciter =
fWeightCalculators.begin();
calciter != fWeightCalculators.end(); calciter++) {
NUISANCEWeightCalc *nuiscalc =
static_cast<NUISANCEWeightCalc *>(*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<NUISANCEWeightCalc *>::iterator iter =
fWeightCalculators.begin();
iter != fWeightCalculators.end(); iter++) {
NUISANCEWeightCalc *nuiscalc = static_cast<NUISANCEWeightCalc *>(*iter);
rw_weight *= nuiscalc->CalcWeight(evt);
}
// Return rw_weight
return rw_weight;
}
diff --git a/src/Reweight/WeightUtils.cxx b/src/Reweight/WeightUtils.cxx
index 8c33b86..59cfa90 100644
--- a/src/Reweight/WeightUtils.cxx
+++ b/src/Reweight/WeightUtils.cxx
@@ -1,614 +1,613 @@
#include "WeightUtils.h"
#include "FitLogger.h"
#ifndef __NO_REWEIGHT__
#ifdef __T2KREW_ENABLED__
#include "T2KGenieReWeight.h"
#include "T2KNIWGReWeight.h"
#include "T2KNIWGUtils.h"
#include "T2KNeutReWeight.h"
#include "T2KNeutUtils.h"
#include "T2KReWeight.h"
using namespace t2krew;
#endif
#ifdef __NIWG_ENABLED__
#include "NIWGReWeight.h"
#include "NIWGSyst.h"
#endif
#ifdef __NEUT_ENABLED__
#include "NReWeight.h"
#include "NSyst.h"
#endif
#ifdef __NUWRO_REWEIGHT_ENABLED__
#include "NuwroReWeight.h"
#include "NuwroSyst.h"
#endif
#ifdef __GENIE_ENABLED__
#ifdef GENIE_PRE_R3
#ifndef __NO_GENIE_REWEIGHT__
#include "ReWeight/GReWeight.h"
#include "ReWeight/GSyst.h"
#endif
#else
using namespace genie;
#ifndef __NO_GENIE_REWEIGHT__
#include "RwFramework/GReWeight.h"
#include "RwFramework/GSyst.h"
using namespace genie::rew;
#endif
#endif
#endif
#ifdef __NOVA_ENABLED__
#include "NOvARwgtEngine.h"
#endif
#ifdef __NUSYST_ENABLED__
#include "nusystematicsWeightEngine.h"
#endif
#endif // end of no reweight
#include "GlobalDialList.h"
#include "ModeNormEngine.h"
#include "NUISANCESyst.h"
#include "OscWeightEngine.h"
//********************************************************************
TF1 FitBase::GetRWConvFunction(std::string const &type,
std::string const &name) {
//********************************************************************
std::string dialfunc = "x";
std::string parType = type;
double low = -10000.0;
double high = 10000.0;
if (parType.find("parameter") == std::string::npos)
parType += "_parameter";
std::string line;
std::ifstream card(
(GeneralUtils::GetTopLevelDir() + "/parameters/dial_conversion.card")
.c_str(),
std::ifstream::in);
while (std::getline(card >> std::ws, line, '\n')) {
std::vector<std::string> inputlist = GeneralUtils::ParseToStr(line, " ");
// Check the line length
if (inputlist.size() < 4)
continue;
// Check whether this is a comment
if (inputlist[0].c_str()[0] == '#')
continue;
// Check whether this is the correct parameter type
if (inputlist[0].compare(parType) != 0)
continue;
// Check the parameter name
if (inputlist[1].compare(name) != 0)
continue;
// inputlist[2] should be the units... ignore for now
dialfunc = inputlist[3];
// High and low are optional, check whether they exist
if (inputlist.size() > 4)
low = GeneralUtils::StrToDbl(inputlist[4]);
if (inputlist.size() > 5)
high = GeneralUtils::StrToDbl(inputlist[5]);
}
TF1 convfunc = TF1((name + "_convfunc").c_str(), dialfunc.c_str(), low, high);
return convfunc;
}
//********************************************************************
std::string FitBase::GetRWUnits(std::string const &type,
std::string const &name) {
//********************************************************************
std::string unit = "sig.";
std::string parType = type;
if (parType.find("parameter") == std::string::npos) {
parType += "_parameter";
}
std::string line;
std::ifstream card(
(GeneralUtils::GetTopLevelDir() + "/parameters/dial_conversion.card")
.c_str(),
std::ifstream::in);
while (std::getline(card >> std::ws, line, '\n')) {
std::vector<std::string> inputlist = GeneralUtils::ParseToStr(line, " ");
// Check the line length
if (inputlist.size() < 3)
continue;
// Check whether this is a comment
if (inputlist[0].c_str()[0] == '#')
continue;
// Check whether this is the correct parameter type
if (inputlist[0].compare(parType) != 0)
continue;
// Check the parameter name
if (inputlist[1].compare(name) != 0)
continue;
unit = inputlist[2];
break;
}
return unit;
}
//********************************************************************
double FitBase::RWAbsToSigma(std::string const &type, std::string const &name,
double val) {
//********************************************************************
TF1 f1 = GetRWConvFunction(type, name);
double conv_val = f1.GetX(val);
if (fabs(conv_val) < 1E-10)
conv_val = 0.0;
NUIS_LOG(FIT, "AbsToSigma(" << name << ") = " << val << " -> " << conv_val);
return conv_val;
}
//********************************************************************
double FitBase::RWSigmaToAbs(std::string const &type, std::string const &name,
double val) {
//********************************************************************
TF1 f1 = GetRWConvFunction(type, name);
double conv_val = f1.Eval(val);
return conv_val;
}
//********************************************************************
double FitBase::RWFracToSigma(std::string const &type, std::string const &name,
double val) {
//********************************************************************
TF1 f1 = GetRWConvFunction(type, name);
double conv_val = f1.GetX((val * f1.Eval(0.0)));
if (fabs(conv_val) < 1E-10)
conv_val = 0.0;
return conv_val;
}
//********************************************************************
double FitBase::RWSigmaToFrac(std::string const &type, std::string const &name,
double val) {
//********************************************************************
TF1 f1 = GetRWConvFunction(type, name);
double conv_val = f1.Eval(val) / f1.Eval(0.0);
return conv_val;
}
int FitBase::ConvDialType(std::string const &type) {
if (!type.compare("neut_parameter"))
return kNEUT;
else if (!type.compare("niwg_parameter"))
return kNIWG;
else if (!type.compare("nuwro_parameter"))
return kNUWRO;
else if (!type.compare("t2k_parameter"))
return kT2K;
else if (!type.compare("genie_parameter"))
return kGENIE;
else if (!type.compare("custom_parameter"))
return kCUSTOM;
else if (!type.compare("norm_parameter"))
return kNORM;
else if (!type.compare("likeweight_parameter"))
return kLIKEWEIGHT;
else if (!type.compare("spline_parameter"))
return kSPLINEPARAMETER;
else if (!type.compare("osc_parameter"))
return kOSCILLATION;
else if (!type.compare("modenorm_parameter"))
return kMODENORM;
else if (!type.compare("nova_parameter"))
return kNOvARWGT;
else if (!type.compare("nusyst_parameter"))
return kNuSystematics;
else
return kUNKNOWN;
}
std::string FitBase::ConvDialType(int type) {
switch (type) {
case kNEUT: {
return "neut_parameter";
}
case kNIWG: {
return "niwg_parameter";
}
case kNUWRO: {
return "nuwro_parameter";
}
case kT2K: {
return "t2k_parameter";
}
case kGENIE: {
return "genie_parameter";
}
case kNORM: {
return "norm_parameter";
}
case kCUSTOM: {
return "custom_parameter";
}
case kLIKEWEIGHT: {
return "likeweight_parameter";
}
case kSPLINEPARAMETER: {
return "spline_parameter";
}
case kOSCILLATION: {
return "osc_parameter";
}
case kMODENORM: {
return "modenorm_parameter";
}
case kNOvARWGT: {
return "nova_parameter";
}
case kNuSystematics: {
return "nusyst_parameter";
}
default:
return "unknown_parameter";
}
}
int FitBase::GetDialEnum(std::string const &type, std::string const &name) {
return FitBase::GetDialEnum(FitBase::ConvDialType(type), name);
}
int FitBase::GetDialEnum(int type, std::string const &name) {
int offset = type * 1000;
int this_enum = Reweight::kNoDialFound; // Not Found
NUIS_LOG(DEB, "Getting dial enum " << type << " " << name);
// Select Types
switch (type) {
// NEUT DIAL TYPE
case kNEUT: {
#if defined(__NEUT_ENABLED__) && !defined(__NO_REWEIGHT__)
int neut_enum = (int)neut::rew::NSyst::FromString(name);
if (neut_enum != 0) {
this_enum = neut_enum + offset;
}
#else
this_enum = Reweight::kNoTypeFound; // Not enabled
#endif
break;
}
// NIWG DIAL TYPE
case kNIWG: {
#if defined(__NIWG_ENABLED__) && !defined(__NO_REWEIGHT__)
int niwg_enum = (int)niwg::rew::NIWGSyst::FromString(name);
if (niwg_enum != 0) {
this_enum = niwg_enum + offset;
}
#else
this_enum = Reweight::kNoTypeFound;
#endif
break;
}
// NUWRO DIAL TYPE
case kNUWRO: {
#if defined(__NUWRO_REWEIGHT_ENABLED__) && !defined(__NO_REWEIGHT__)
int nuwro_enum = (int)nuwro::rew::NuwroSyst::FromString(name);
if (nuwro_enum > 0) {
this_enum = nuwro_enum + offset;
}
#else
this_enum = Reweight::kNoTypeFound;
#endif
}
// GENIE DIAL TYPE
case kGENIE: {
#if defined(__GENIE_ENABLED__) && !defined(__NO_REWEIGHT__)
int genie_enum = (int)genie::rew::GSyst::FromString(name);
if (genie_enum > 0) {
this_enum = genie_enum + offset;
}
#else
this_enum = Reweight::kNoTypeFound;
#endif
break;
}
case kCUSTOM: {
int custom_enum = 0; // PLACEHOLDER
this_enum = custom_enum + offset;
break;
}
// T2K DIAL TYPE
case kT2K: {
#if defined(__T2KREW_ENABLED__) && !defined(__NO_REWEIGHT__)
int t2k_enum = (int)t2krew::T2KSyst::FromString(name);
if (t2k_enum > 0) {
this_enum = t2k_enum + offset;
}
#else
this_enum = Reweight::kNoTypeFound;
#endif
break;
}
case kNORM: {
if (gNormEnums.find(name) == gNormEnums.end()) {
gNormEnums[name] = gNormEnums.size() + 1 + offset;
}
this_enum = gNormEnums[name];
break;
}
case kLIKEWEIGHT: {
if (gLikeWeightEnums.find(name) == gLikeWeightEnums.end()) {
gLikeWeightEnums[name] = gLikeWeightEnums.size() + 1 + offset;
}
this_enum = gLikeWeightEnums[name];
break;
}
case kSPLINEPARAMETER: {
if (gSplineParameterEnums.find(name) == gSplineParameterEnums.end()) {
gSplineParameterEnums[name] = gSplineParameterEnums.size() + 1 + offset;
}
this_enum = gSplineParameterEnums[name];
break;
}
case kOSCILLATION: {
#ifdef __PROB3PP_ENABLED__
int oscEnum = OscWeightEngine::SystEnumFromString(name);
if (oscEnum != 0) {
this_enum = oscEnum + offset;
}
#else
this_enum = Reweight::kNoTypeFound; // Not enabled
#endif
}
case kMODENORM: {
size_t us_pos = name.find_first_of('_');
std::string numstr = name.substr(us_pos + 1);
int mode_num = std::atoi(numstr.c_str());
NUIS_LOG(FTL, "Getting mode num " << mode_num);
if (!mode_num) {
NUIS_ABORT("Attempting to parse dial name: \""
<< name << "\" as a mode norm dial but failed.");
}
this_enum = 60 + mode_num + offset;
break;
}
}
// If Not Enabled
if (this_enum == Reweight::kNoTypeFound) {
NUIS_ERR(FTL,
"RW Engine not supported for " << FitBase::ConvDialType(type));
NUIS_ABORT("Check dial " << name);
}
// If Not Found
if (this_enum == Reweight::kNoDialFound) {
NUIS_ABORT("Dial " << name << " not found.");
}
return this_enum;
}
int Reweight::ConvDialType(std::string const &type) {
return FitBase::ConvDialType(type);
}
std::string Reweight::ConvDialType(int type) {
return FitBase::ConvDialType(type);
}
int Reweight::GetDialType(int type) {
int t = (type / 1000);
return t > kNuSystematics ? Reweight::kNoDialFound : t;
}
int Reweight::RemoveDialType(int type) { return (type % 1000); }
int Reweight::NEUTEnumFromName(std::string const &name) {
#if defined(__NEUT_ENABLED__) && !defined(__NO_REWEIGHT__)
int neutenum = (int)neut::rew::NSyst::FromString(name);
return (neutenum > 0) ? neutenum : Reweight::kNoDialFound;
#else
return Reweight::kGeneratorNotBuilt;
#endif
}
int Reweight::NIWGEnumFromName(std::string const &name) {
#if defined(__NIWG_ENABLED__) && !defined(__NO_REWEIGHT__)
int niwgenum = (int)niwg::rew::NIWGSyst::FromString(name);
return (niwgenum != 0) ? niwgenum : Reweight::kNoDialFound;
#else
return Reweight::kGeneratorNotBuilt;
#endif
}
int Reweight::NUWROEnumFromName(std::string const &name) {
#if defined(__NUWRO_REWEIGHT_ENABLED__) && !defined(__NO_REWEIGHT__)
int nuwroenum = (int)nuwro::rew::NuwroSyst::FromString(name);
return (nuwroenum > 0) ? nuwroenum : Reweight::kNoDialFound;
#else
return Reweight::kGeneratorNotBuilt;
#endif
}
int Reweight::GENIEEnumFromName(std::string const &name) {
#if defined(__GENIE_ENABLED__) && !defined(__NO_REWEIGHT__)
int genieenum = (int)genie::rew::GSyst::FromString(name);
return (genieenum > 0) ? genieenum : Reweight::kNoDialFound;
#else
return Reweight::kGeneratorNotBuilt;
#endif
}
int Reweight::T2KEnumFromName(std::string const &name) {
#if defined(__T2KREW_ENABLED__) && !defined(__NO_REWEIGHT__)
int t2kenum = (int)t2krew::T2KSyst::FromString(name);
return (t2kenum > 0) ? t2kenum : Reweight::kNoDialFound;
#else
return Reweight::kGeneratorNotBuilt;
#endif
}
int Reweight::OscillationEnumFromName(std::string const &name) {
#ifdef __PROB3PP_ENABLED__
int oscEnum = OscWeightEngine::SystEnumFromString(name);
return (oscEnum > 0) ? oscEnum : Reweight::kNoDialFound;
#else
return Reweight::kGeneratorNotBuilt;
#endif
}
int Reweight::NUISANCEEnumFromName(std::string const &name, int type) {
int nuisenum = Reweight::DialList().EnumFromNameAndType(name, type);
return nuisenum;
}
int Reweight::CustomEnumFromName(std::string const &name) {
int custenum = Reweight::ConvertNUISANCEDial(name);
- return custenum;
+ return (custenum != kUnknownNUISANCEDial ? custenum : kNoDialFound);
}
-int Reweight::ConvDial(std::string const &name, std::string const &type,
- bool exceptions) {
+int Reweight::ConvDial(std::string const &name, std::string const &type, bool exceptions) {
return Reweight::ConvDial(name, Reweight::ConvDialType(type), exceptions);
}
int Reweight::ConvDial(std::string const &fullname, int type, bool exceptions) {
std::string name =
GeneralUtils::ParseToStr(fullname, ",")[0]; // Only use first dial given
// Produce offset seperating each type.
int offset = type * 1000;
int genenum = Reweight::kNoDialFound;
switch (type) {
case kNEUT:
genenum = NEUTEnumFromName(name);
break;
case kNIWG:
genenum = NIWGEnumFromName(name);
break;
case kNUWRO:
genenum = NUWROEnumFromName(name);
break;
case kGENIE:
genenum = GENIEEnumFromName(name);
break;
case kT2K:
genenum = T2KEnumFromName(name);
break;
case kCUSTOM:
genenum = CustomEnumFromName(name);
break;
case kNORM:
case kLIKEWEIGHT:
case kSPLINEPARAMETER:
case kNEWSPLINE:
genenum = NUISANCEEnumFromName(name, type);
break;
case kOSCILLATION:
genenum = OscillationEnumFromName(name);
break;
case kMODENORM:
genenum = ModeNormEngine::SystEnumFromString(name);
break;
#ifdef __NOVA_ENABLED__
case kNOvARWGT:
genenum = NOvARwgtEngine::GetWeightGeneratorIndex(name);
break;
#endif
#ifdef __NUSYST_ENABLED__
case kNuSystematics: {
// Super inefficient...
nusystematicsWeightEngine we;
genenum = we.ConvDial(name);
break;
}
#endif
default:
genenum = Reweight::kNoTypeFound;
break;
}
// Throw if required.
if (exceptions) {
// If Not Enabled
if (genenum == Reweight::kGeneratorNotBuilt) {
NUIS_ERR(FTL,
"RW Engine not supported for " << FitBase::ConvDialType(type));
NUIS_ABORT("Check dial " << name);
}
// If no type enabled
if (genenum == Reweight::kNoTypeFound) {
NUIS_ABORT("Type mismatch inside ConvDialEnum");
}
// If Not Found
if (genenum == Reweight::kNoDialFound) {
NUIS_ABORT("Dial " << name << " not found.");
}
}
// Add offset if no issue
int nuisenum = genenum;
if ((genenum != Reweight::kGeneratorNotBuilt) &&
(genenum != Reweight::kNoTypeFound) &&
(genenum != Reweight::kNoDialFound)) {
nuisenum += offset;
}
// Now register dial
Reweight::DialList().RegisterDialEnum(name, type, nuisenum);
return nuisenum;
}
std::string Reweight::ConvDial(int nuisenum) {
for (size_t i = 0; i < Reweight::DialList().fAllDialEnums.size(); i++) {
if (Reweight::DialList().fAllDialEnums[i] == nuisenum) {
return Reweight::DialList().fAllDialNames[i];
}
}
NUIS_LOG(FIT, "Cannot find dial with enum = " << nuisenum);
return "";
}
File Metadata
Details
Attached
Mime Type
text/x-diff
Expires
Sun, Feb 23, 3:04 PM (38 m, 2 s)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4486832
Default Alt Text
(88 KB)
Attached To
rNUISANCEGIT nuisancegit
Event Timeline
Log In to Comment