Page Menu
Home
HEPForge
Search
Configure Global Search
Log In
Files
F10275468
No One
Temporary
Actions
View File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Flag For Later
Size
36 KB
Subscribers
None
View Options
diff --git a/src/Reweight/NUISANCEWeightCalcs.cxx b/src/Reweight/NUISANCEWeightCalcs.cxx
index c3f919a..acdf04c 100644
--- a/src/Reweight/NUISANCEWeightCalcs.cxx
+++ b/src/Reweight/NUISANCEWeightCalcs.cxx
@@ -1,1006 +1,1037 @@
#include "NUISANCEWeightCalcs.h"
#include "FitEvent.h"
#include "GeneralUtils.h"
#include "NUISANCESyst.h"
#include "WeightUtils.h"
#include "WeightEngineBase.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 % NUIS_DIAL_OFFSET;
// Check Handled
if (!IsHandled(curenum))
return;
if (curenum == kModeNorm_NormRES)
fNormRES = val;
}
bool ModeNormCalc::IsHandled(int rwenum) {
int curenum = rwenum % NUIS_DIAL_OFFSET;
switch (curenum) {
case kModeNorm_NormRES:
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 || !fApply_MINOSRPA) return 1.0;
double w = 1.0;
// If GENIE is enabled, use old code
#ifdef __GENIE_ENABLED__
// 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 on nucleus, not resonant, or NC
if (!tgt.IsNucleus() || !proc_info.IsResonant() || proc_info.IsWeakNC()) 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());
w *= GetRPAWeight(Q2);
#else
// Get the Q2 from NUISANCE if not GENIE
FitEvent *fevt = static_cast<FitEvent*>(evt);
// Check the event is resonant
if (!fevt->IsResonant() || fevt->IsNC()) return 1.0;
int targeta = fevt->GetTargetA();
int targetz = fevt->GetTargetZ();
// Apply only to nuclear targets, ignore free protons
if (targeta == 1 || targetz == 1) return 1.0;
// Q2 in GeV2
double Q2 = fevt->GetQ2();
w *= GetRPAWeight(Q2);
#endif
return w;
}
// Do the actual weight calculation
double MINOSRPA::GetRPAWeight(double Q2) {
if (Q2 > 0.7) return 1.0;
double w = fCur_MINOSRPA_A / (1.0 + TMath::Exp(1.0 - TMath::Sqrt(Q2) / fCur_MINOSRPA_B));
return w;
}
bool MINOSRPA::IsHandled(int rwenum) {
int curenum = rwenum % NUIS_DIAL_OFFSET;
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 % NUIS_DIAL_OFFSET;
// 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 || !fApplyRPA) return 1.0;
double w = 1.0;
// If GENIE is enabled, use old code
#ifdef __GENIE_ENABLED__
// 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 on nucleus, not resonant, or NC
if (!tgt.IsNucleus() || !proc_info.IsResonant() || proc_info.IsWeakNC()) 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());
w *= GetRPAWeight(Q2);
#else
// Get the Q2 from NUISANCE if not GENIE
FitEvent *fevt = static_cast<FitEvent*>(evt);
// Check the event is resonant
if (!fevt->IsResonant() || fevt->IsNC()) return 1.0;
int targeta = fevt->GetTargetA();
int targetz = fevt->GetTargetZ();
// Apply only to nuclear targets, ignore free protons
if (targeta == 1 || targetz == 1) return 1.0;
// Q2 in GeV2
double Q2 = fevt->GetQ2();
w *= GetRPAWeight(Q2);
#endif
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 % NUIS_DIAL_OFFSET;
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 % NUIS_DIAL_OFFSET;
// 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);
}
//
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) {
int mode = abs(evt->Mode);
double w = 1.0;
if (nParams == 0) {
return w;
}
// Get Q2
// Get final state lepton
if (mode == 1) {
FitEvent *fevt = static_cast<FitEvent*>(evt);
double Q2 = fevt->GetQ2();
// 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 % NUIS_DIAL_OFFSET;
// 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 % NUIS_DIAL_OFFSET;
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 % NUIS_DIAL_OFFSET;
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 % NUIS_DIAL_OFFSET;
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->GetNeutrinoIn();
if (!pnu) {
NUIS_ABORT("NO Starting particle found in stack!");
}
FitParticle *plep = fevt->GetLeptonOut();
if (!plep) return 1.0;
TLorentzVector q = pnu->fP - plep->fP;
// Extra q0,q3
double q0 = fabs(q.E()) / 1.E3;
double q3 = fabs(q.Vect().Mag()) / 1.E3;
int initialstate = -1; // Undef
if (abs(fevt->Mode) == 2) {
int npr = 0;
int nne = 0;
for (UInt_t j = 0; j < fevt->Npart(); j++) {
if ((fevt->PartInfo(j))->fIsAlive)
continue;
if (fevt->PartInfo(j)->fPID == 2212)
npr++;
else if (fevt->PartInfo(j)->fPID == 2112)
nne++;
}
if (fevt->Mode == 2 && npr == 1 && nne == 1) {
initialstate = 2;
} else if (fevt->Mode == 2 &&
((npr == 0 && nne == 2) || (npr == 2 && nne == 0))) {
initialstate = 1;
}
}
// Apply weighting
if (fApply_CCQE && abs(fevt->Mode) == 1) {
if (fDebugStatements)
std::cout << "Getting CCQE Weight" << std::endl;
double g = GetGausWeight(q0, q3, fGausVal_CCQE);
if (g < 1.0)
g = 1.0;
rw_weight *= g;
}
if (fApply_2p2h && abs(fevt->Mode) == 2) {
if (fDebugStatements) std::cout << "Getting 2p2h Weight" << std::endl;
if (fDebugStatements) std::cout << "Got q0 q3 = " << q0 << " " << q3 << " mode = " << fevt->Mode << std::endl;
rw_weight *= GetGausWeight(q0, q3, fGausVal_2p2h);
if (fDebugStatements) std::cout << "Returning Weight " << rw_weight << std::endl;
}
if (fApply_2p2h_PPandNN && abs(fevt->Mode) == 2 && initialstate == 1) {
if (fDebugStatements) std::cout << "Getting 2p2h PPandNN Weight" << std::endl;
rw_weight *= GetGausWeight(q0, q3, fGausVal_2p2h_PPandNN);
}
if (fApply_2p2h_NP && abs(fevt->Mode) == 2 && initialstate == 2) {
if (fDebugStatements) std::cout << "Getting 2p2h NP Weight" << std::endl;
rw_weight *= GetGausWeight(q0, q3, fGausVal_2p2h_NP);
}
if (fApply_CC1pi && abs(fevt->Mode) >= 11 && abs(fevt->Mode) <= 13) {
if (fDebugStatements)
std::cout << "Getting CC1pi Weight" << std::endl;
rw_weight *= GetGausWeight(q0, q3, fGausVal_CC1pi);
}
return rw_weight;
}
void GaussianModeCorr::SetMethod(bool method) {
fMethod = method;
if (fMethod == true) {
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 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) );
}
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
}
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 % NUIS_DIAL_OFFSET;
// 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 % NUIS_DIAL_OFFSET;
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;
}
}
// Radcorr weight calculator
RadCorrQ2::RadCorrQ2() {
// Input files for the splines
TFile *fInputs[kNumuBar];
std::string fileloc = std::string(std::getenv("NUISANCE"));
fileloc += "/data/radcorr";
// Two files for numu and numubar
fInputs[kNumu] = new TFile((fileloc+"/Elspectrum_muon_neutrino_merge.root").c_str());
fInputs[kNumuBar] = new TFile((fileloc+"/Elspectrum_muon_antineutrino_merge.root").c_str());
//EnuRange[9] = {0.2, 0.3, 0.5, 1., 2., 3., 5., 10., 20.};
EnuRange[0] = 0.2;
EnuRange[1] = 0.3;
EnuRange[2] = 0.5;
- EnuRange[3] = 1.0;
- EnuRange[4] = 2.0;
- EnuRange[5] = 3.0;
- EnuRange[6] = 5.0;
- EnuRange[7] = 10.;
- EnuRange[8] = 20.;
-
- nEnu = 9;
+ EnuRange[3] = 0.6;
+ EnuRange[4] = 0.75;
+ EnuRange[5] = 0.9;
+ EnuRange[6] = 1.0;
+ EnuRange[7] = 2.0;
+ EnuRange[8] = 3.0;
+ EnuRange[9] = 5.0;
+ EnuRange[10] = 10.;
+ EnuRange[11] = 20.;
+
+ nEnu = 12;
if (nEnu != sizeof(EnuRange)/sizeof(double)) {
std::cerr << "Something wrong with Enu range" << std::endl;
throw;
}
// The spline base names
std::string basename_numu = "Elspectrum_muon_neutrino_neutron";
std::string basename_numub = "Elspectrum_muon_antineutrino_proton";
Graphs = new TGraph**[kNumuBar+1];
for (int i = 0; i < kNumuBar+1; ++i) {
Graphs[i] = new TGraph*[nEnu];
std::string basename;
if (i == kNumu) basename = basename_numu;
else basename = basename_numub;
Graphs[i][0] = (TGraph*)fInputs[i]->Get((basename+"_02GeV").c_str())->Clone();
Graphs[i][1] = (TGraph*)fInputs[i]->Get((basename+"_03GeV").c_str())->Clone();
Graphs[i][2] = (TGraph*)fInputs[i]->Get((basename+"_05GeV").c_str())->Clone();
- Graphs[i][3] = (TGraph*)fInputs[i]->Get((basename+"_1GeV").c_str())->Clone();
- Graphs[i][4] = (TGraph*)fInputs[i]->Get((basename+"_2GeV").c_str())->Clone();
- Graphs[i][5] = (TGraph*)fInputs[i]->Get((basename+"_3GeV").c_str())->Clone();
- Graphs[i][6] = (TGraph*)fInputs[i]->Get((basename+"_5GeV").c_str())->Clone();
- Graphs[i][7] = (TGraph*)fInputs[i]->Get((basename+"_10GeV").c_str())->Clone();
- Graphs[i][8] = (TGraph*)fInputs[i]->Get((basename+"_20GeV").c_str())->Clone();
+ Graphs[i][3] = (TGraph*)fInputs[i]->Get((basename+"_06GeV").c_str())->Clone();
+ Graphs[i][4] = (TGraph*)fInputs[i]->Get((basename+"_075GeV").c_str())->Clone();
+ Graphs[i][5] = (TGraph*)fInputs[i]->Get((basename+"_09GeV").c_str())->Clone();
+ Graphs[i][6] = (TGraph*)fInputs[i]->Get((basename+"_1GeV").c_str())->Clone();
+ Graphs[i][7] = (TGraph*)fInputs[i]->Get((basename+"_2GeV").c_str())->Clone();
+ Graphs[i][8] = (TGraph*)fInputs[i]->Get((basename+"_3GeV").c_str())->Clone();
+ Graphs[i][9] = (TGraph*)fInputs[i]->Get((basename+"_5GeV").c_str())->Clone();
+ Graphs[i][10] = (TGraph*)fInputs[i]->Get((basename+"_10GeV").c_str())->Clone();
+ Graphs[i][11] = (TGraph*)fInputs[i]->Get((basename+"_20GeV").c_str())->Clone();
}
fInputs[0]->Close();
fInputs[1]->Close();
+
+ // Default to not use
+ type = 0;
}
// Calculate the weight from the radiative correction in Q2
// Function of Enu and Q2, using linear interpolation
double RadCorrQ2::CalcWeight(BaseFitEvt *evt) {
- if (!fUse) return 1.0;
+
+ // If we don't want this calculation
+ if (type == 0) return 1.0;
// Check interaction mode is CCQE
- if (abs(evt->Mode) != 1) return 1.0;
+ // If CCQE only requested, return if
+ if (type == 1 && abs(evt->Mode) != 1) return 1.0;
// Only apply to muon (anti)neutrinos
if (abs(evt->probe_pdg) != 14) return 1.0;
// Get the neutrino type
NuType nutype;
if (evt->Mode > 0) nutype = kNumu;
else nutype = kNumuBar;
// Get the event Q2 and Enu
FitEvent *fevt = static_cast<FitEvent*>(evt);
double Q2 = fevt->GetQ2(); // Get in GeV2
double Enu = fevt->GetNeutrinoIn()->E()/1.E3; // Convert to GeV
+ // Apply only to events with outgoing muon
+ if (abs(fevt->GetLeptonOutPDG()) != 13) return 1.0;
+ if (fevt->GetLeptonOut() == NULL) return 1.0;
+
// Since all happens on the nucleon, need to boost into nucleon frame and use Enu there
TLorentzVector initnu = fevt->GetBeamNeutrinoP4();
// Then get struck nucleon
const int pdg[] = {2112, 2212};
FitParticle* initialstate = fevt->GetHMISParticle(pdg);
+ if (initialstate == NULL) return 1.0; // Don't apply to coherent?
TLorentzVector initnuc = initialstate->P4();
// Boost the neutrino
initnu.Boost(-1*initnuc.BoostVector());
//std::cout << Enu << " " << initnu.E()/1E3 << std::endl;
Enu = initnu.E()/1E3; // update to be enu in nucleon frame
// Pick some random Q2 and Enu
//double Q2 = 1;
//double Enu = 2;
// Find the nearest point in Enu
// Get the true neutrino energy and interpolate between nearest points
// This is always the point below our actual point
int nearest = 0;
for (int j = 0; j < nEnu; ++j) {
if (Enu > EnuRange[j]) nearest = j;
}
// Then get the maximum Q2 for each of the Enu points
// to make sure we're not extrapolating unphysically
double mumass = fevt->GetLeptonOut()->M()/1.E3; //convert to GeV
double Q2max_near = GetQ2max(EnuRange[nearest], mumass);
// The allowed Q2 is always going to be lower at the lower energy bin
// If this is the case, evaluate the spline just before the drop off instead
double low = 0;
// Sometimes the Q2 will be right on the edge of allowed
if (Q2-Q2max_near > -0.02) {
// Try to find the spot where we're no longer dropping abruptly
TGraph *g = Graphs[nutype][nearest];
double *y = g->GetY();
int npoints = g->GetN();
double prevy = y[0];
// Which point does the job occur at
int jump = 0;
// Find at what point the jump happens
for (int j = 0; j < npoints; ++j) {
// Look for jumps in the points
if (prevy-y[j] > 0.9) {
jump = j;
break;
}
prevy = y[j];
}
low = y[jump-1];
} else {
low = Graphs[nutype][nearest]->Eval(Q2, 0, "");
}
// This is the point above
int nextbin = nearest+1;
// The Enu might be beyond our last range
if (Enu > EnuRange[nEnu-1]) nextbin = nearest;
// Now also need to check the high Q2
double Q2max_far = GetQ2max(EnuRange[nextbin], mumass);
double high = 0;
// Sometimes the Q2 will be right on the edge of allowed
if (Q2max_far-Q2 < 0.02) {
// Try to find the spot where we're no longer dropping abruptly
TGraph *g = Graphs[nutype][nextbin];
double *y = g->GetY();
int npoints = g->GetN();
double prevy = y[0];
// Which point does the job occur at
int jump = 0;
// Find at what point the jump happens
for (int j = 0; j < npoints; ++j) {
// Look for jumps in the points
if (prevy-y[j] > 0.9) {
jump = j;
break;
}
prevy = y[j];
}
high = y[jump-1];
} else {
high = Graphs[nutype][nextbin]->Eval(Q2, 0, "");
}
// linear intepolation
double weight = (high-low)*(Enu-EnuRange[nearest])/(EnuRange[nearest+1]-EnuRange[nearest])+low;
//std::cout << Q2 << " " << Enu << " " << weight << std::endl;
// Put in a weight cap
if (weight > 10) weight = 10;
if (weight < 0) weight = 0;
return weight;
}
// Get the max Q2 for a given Enu to check interpolaton
double RadCorrQ2::GetQ2max(double Enu, double lepmass) {
// Nucleon mass
const double Mn = 0.93956542052;
const double Mp = 0.93827208816;
const double M = (Mn+Mp)/2.;
double val = -(M+Enu)*lepmass*lepmass+2*M*Enu*Enu+Enu*sqrt(pow(2*M*Enu-lepmass*lepmass, 2)-4*lepmass*lepmass*M*M);
val /= (M+2*Enu);
return val;
}
bool RadCorrQ2::IsHandled(int rwenum) {
int curenum = rwenum % NUIS_DIAL_OFFSET;
switch (curenum) {
case Reweight::kRadCorrQ2:
return true;
default:
return false;
}
}
void RadCorrQ2::SetDialValue(std::string name, double val) {
SetDialValue(Reweight::ConvDial(name, kCUSTOM), val);
}
void RadCorrQ2::SetDialValue(int rwenum, double val) {
int curenum = rwenum % NUIS_DIAL_OFFSET;
// Check Handled
if (!IsHandled(curenum)) return;
// If greater than 0.5, use
if (curenum == kRadCorrQ2) {
- if (val > 0.5) fUse = true;
- else fUse = false;
+ if (val >= 0 && val < 0.5) type = 0;
+ else if (val > 0 && val < 1.5) type = 1;
+ else if (val > 0 && val < 2.5) type = 2;
+ else {
+ std::cerr << "wrong value given to radiative correction. Please give 0, 1, 2. I will round to nearest integer" << std::endl;
+ throw;
+ }
+ }
+}
+
+RadCorrQ2::~RadCorrQ2() {
+ for (int i = 0; i < kNumuBar+1; ++i) {
+ for (int j = 0; j < 12; ++j) {
+ delete Graphs[i][j];
+ }
+ delete[] Graphs[i];
}
}
diff --git a/src/Reweight/NUISANCEWeightCalcs.h b/src/Reweight/NUISANCEWeightCalcs.h
index 0b73f23..279ad7a 100644
--- a/src/Reweight/NUISANCEWeightCalcs.h
+++ b/src/Reweight/NUISANCEWeightCalcs.h
@@ -1,226 +1,226 @@
#ifndef NUISANCE_WEIGHT_CALCS
#define NUISANCE_WEIGHT_CALCS
#include "BaseFitEvt.h"
#include "BeRPA.h"
#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"
using namespace genie;
#endif // End GENIE pre v3
#endif // End GENIE enabled
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<std::string, int> fDialNameIndex;
std::map<int, int> fDialEnumIndex;
std::vector<double> 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;
};
// 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);
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;
};
class BeRPACalc : public NUISANCEWeightCalc {
public:
BeRPACalc();
~BeRPACalc(){};
double CalcWeight(BaseFitEvt* evt);
void SetDialValue(std::string name, double val);
void SetDialValue(int rwenum, double val);
bool IsHandled(int rwenum);
private:
// Parameter values
double fBeRPA_A;
double fBeRPA_B;
double fBeRPA_D;
double fBeRPA_E;
double fBeRPA_U;
// Counts of enabled parameters
int nParams;
};
class SBLOscWeightCalc : public NUISANCEWeightCalc {
public:
SBLOscWeightCalc();
~SBLOscWeightCalc(){};
double CalcWeight(BaseFitEvt* evt);
void SetDialValue(std::string name, double val);
void SetDialValue(int rwenum, double val);
bool IsHandled(int rwenum);
double GetSBLOscWeight(double E);
double fDistance;
double fMassSplitting;
double fSin2Theta;
};
class GaussianModeCorr : public NUISANCEWeightCalc {
public:
GaussianModeCorr();
~GaussianModeCorr(){};
double CalcWeight(BaseFitEvt* evt);
void SetDialValue(std::string name, double val);
void SetDialValue(int rwenum, double val);
bool IsHandled(int rwenum);
double GetGausWeight(double q0, double q3, double vals[]);
// Set the Gaussian method (tilt-shift or normal Gaussian parameters)
void SetMethod(bool method);
// 5 pars describe the Gaussain
// 0 norm.
// 1 q0 pos
// 2 q0 width
// 3 q3 pos
// 4 q3 width
static const int kPosNorm = 0;
static const int kPosTilt = 1;
static const int kPosPq0 = 2;
static const int kPosWq0 = 3;
static const int kPosPq3 = 4;
static const int kPosWq3 = 5;
bool fApply_CCQE;
double fGausVal_CCQE[6];
bool fApply_2p2h;
double fGausVal_2p2h[6];
bool fApply_2p2h_PPandNN;
double fGausVal_2p2h_PPandNN[6];
bool fApply_2p2h_NP;
double fGausVal_2p2h_NP[6];
bool fApply_CC1pi;
double fGausVal_CC1pi[6];
bool fAllowSuppression;
bool fDebugStatements;
bool fMethod;
};
class RadCorrQ2 : public NUISANCEWeightCalc {
public:
RadCorrQ2(); // Constructor sets up files
- ~RadCorrQ2() {};
+ ~RadCorrQ2();
double CalcWeight(BaseFitEvt *evt);
void SetDialValue(std::string name, double val);
void SetDialValue(int rwenum, double val);
bool IsHandled(int rwenum);
- bool fUse;
+ int type;
private:
// Just a handy enum
enum NuType {kNumu = 0, kNumuBar = 1};
// Input TGraphs, first index is numu or numubar, second index is the fixed enu
TGraph ***Graphs;
// Enu Range that the inputs come in
- double EnuRange[9];
+ double EnuRange[12];
int nEnu;
double GetQ2max(double, double);
};
#endif
File Metadata
Details
Attached
Mime Type
text/x-diff
Expires
Fri, Apr 4, 9:36 PM (16 h, 2 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4729312
Default Alt Text
(36 KB)
Attached To
rNUISANCEGIT nuisancegit
Event Timeline
Log In to Comment