diff --git a/config/GSystUncertaintyTable.xml b/config/GSystUncertaintyTable.xml index 56c2033a..903ded2b 100644 --- a/config/GSystUncertaintyTable.xml +++ b/config/GSystUncertaintyTable.xml @@ -294,8 +294,8 @@ values for each Reweight tweak dial. All of these have default values of zero. 41.4534 41.4534 - - - data/covariances/ELZExpFF.dat + 1.0 + 1.0 + diff --git a/src/Apps/gRwght1Param.cxx b/src/Apps/gRwght1Param.cxx index d975d7d9..412c22d6 100644 --- a/src/Apps/gRwght1Param.cxx +++ b/src/Apps/gRwght1Param.cxx @@ -110,6 +110,7 @@ #include "RwCalculators/GReWeightNuXSecNCEL.h" #include "RwCalculators/GReWeightNuXSecCCQE.h" #include "RwCalculators/GReWeightNuXSecCCQEELFF.h" +#include "RwCalculators/GReWeightNuXSecCCQEZAFF.h" #include "RwCalculators/GReWeightNuXSecCCRES.h" #include "RwCalculators/GReWeightNuXSecCOH.h" #include "RwCalculators/GReWeightNonResonanceBkg.h" @@ -235,6 +236,7 @@ int main(int argc, char ** argv) rw.AdoptWghtCalc( "xsec_ncel", new GReWeightNuXSecNCEL ); rw.AdoptWghtCalc( "xsec_ccqe", new GReWeightNuXSecCCQE ); rw.AdoptWghtCalc( "xsec_ccqe_elff", new GReWeightNuXSecCCQEELFF ); + rw.AdoptWghtCalc( "xsec_ccqe_zaff", new GReWeightNuXSecCCQEZAFF ); rw.AdoptWghtCalc( "xsec_ccqe_axial", new GReWeightNuXSecCCQEaxial ); //rwh - xsec_ccqe_vec is problematic for various tunes rw.AdoptWghtCalc( "xsec_ccqe_vec", new GReWeightNuXSecCCQEvec ); @@ -261,11 +263,14 @@ int main(int argc, char ** argv) rw.AdoptWghtCalc( "delta_rad", new GReWeightDeltaradAngle); // Get GSystSet and include the (single) input systematic parameter + GSystSet & syst = rw.Systematics(); syst.Init(gOptSyst); // Fine-tune weight calculators + + if ( gOptSyst == kXSecTwkDial_MaCCQE ) { // By default GReWeightNuXSecCCQE is in `NormAndMaShape' mode diff --git a/src/Apps/gRwghtNCorrelatedParams.cxx b/src/Apps/gRwghtNCorrelatedParams.cxx index 9c2acfae..0ddc2417 100644 --- a/src/Apps/gRwghtNCorrelatedParams.cxx +++ b/src/Apps/gRwghtNCorrelatedParams.cxx @@ -895,6 +895,7 @@ void AdoptWeightCalcs (vector lsyst, GReWeight & rw) rw.AdoptWghtCalc( "xsec_mec", new GReWeightXSecMEC ); } break; + /* case kXSecTwkDial_ZExpELFF_AN1: case kXSecTwkDial_ZExpELFF_AN2: case kXSecTwkDial_ZExpELFF_AN3: @@ -916,6 +917,7 @@ void AdoptWeightCalcs (vector lsyst, GReWeight & rw) rw.AdoptWghtCalc( "xsec_ccqe_elff", new GReWeightNuXSecCCQEELFF ); } break; + */ default: // no fine-tuning needed break; } diff --git a/src/RwCalculators/GReWeightNuXSecCCQE.cxx b/src/RwCalculators/GReWeightNuXSecCCQE.cxx index 427fb6cb..a1e827ee 100644 --- a/src/RwCalculators/GReWeightNuXSecCCQE.cxx +++ b/src/RwCalculators/GReWeightNuXSecCCQE.cxx @@ -321,7 +321,7 @@ void GReWeightNuXSecCCQE::Reconfigure(void) for (int i=0;iGetDouble(alg_key.str()); } else diff --git a/src/RwCalculators/GReWeightNuXSecCCQEELFF.cxx b/src/RwCalculators/GReWeightNuXSecCCQEELFF.cxx index 52b67b92..1b9fab8e 100644 --- a/src/RwCalculators/GReWeightNuXSecCCQEELFF.cxx +++ b/src/RwCalculators/GReWeightNuXSecCCQEELFF.cxx @@ -66,7 +66,6 @@ GReWeightNuXSecCCQEELFF::GReWeightNuXSecCCQEELFF(std::string model, std::string GReWeightNuXSecCCQEELFF::~GReWeightNuXSecCCQEELFF() { if ( fXSecModelConfig ) delete fXSecModelConfig; - if ( fXSecModel ) delete fXSecModel; if ( fXSecModelDef ) delete fXSecModelDef; } @@ -78,28 +77,6 @@ bool GReWeightNuXSecCCQEELFF::IsHandled(GSyst_t syst) const switch(syst) { // add ZExp vector parameters - case ( kXSecTwkDial_ZExpELFF_AN1 ) : - case ( kXSecTwkDial_ZExpELFF_AN2 ) : - case ( kXSecTwkDial_ZExpELFF_AN3 ) : - case ( kXSecTwkDial_ZExpELFF_AN4 ) : - case ( kXSecTwkDial_ZExpELFF_AP1 ) : - case ( kXSecTwkDial_ZExpELFF_AP2 ) : - case ( kXSecTwkDial_ZExpELFF_AP3 ) : - case ( kXSecTwkDial_ZExpELFF_AP4 ) : - case ( kXSecTwkDial_ZExpELFF_BN1 ) : - case ( kXSecTwkDial_ZExpELFF_BN2 ) : - case ( kXSecTwkDial_ZExpELFF_BN3 ) : - case ( kXSecTwkDial_ZExpELFF_BN4 ) : - case ( kXSecTwkDial_ZExpELFF_BP1 ) : - case ( kXSecTwkDial_ZExpELFF_BP2 ) : - case ( kXSecTwkDial_ZExpELFF_BP3 ) : - case ( kXSecTwkDial_ZExpELFF_BP4 ) : - if(fMode==kModeZExp && fModelIsZExp){ - handle = true; - }else { - handle = false; - } - break; case ( kXSecTwkDial_ZExpELFF ) : if(fMode==kModeZExp && fModelIsZExp){ handle = true; @@ -136,71 +113,6 @@ void GReWeightNuXSecCCQEELFF::SetSystematic(GSyst_t syst, double twk_dial) switch(syst) { case (kXSecTwkDial_ZExpELFF): fZExpTwkDial = twk_dial; - fIsAllPara = true; - break; - case ( kXSecTwkDial_ZExpELFF_AP1 ) : - fZExpParaTwkDial.fZ_APn[0] = twk_dial; - fIsSinglePara = true; - break; - case ( kXSecTwkDial_ZExpELFF_AP2 ) : - fZExpParaTwkDial.fZ_APn[1] = twk_dial; - fIsSinglePara = true; - break; - case ( kXSecTwkDial_ZExpELFF_AP3 ) : - fZExpParaTwkDial.fZ_APn[2] = twk_dial; - fIsSinglePara = true; - break; - case ( kXSecTwkDial_ZExpELFF_AP4 ) : - fZExpParaTwkDial.fZ_APn[3] = twk_dial; - fIsSinglePara = true; - break; - case ( kXSecTwkDial_ZExpELFF_AN1 ) : - fZExpParaTwkDial.fZ_ANn[0] = twk_dial; - fIsSinglePara = true; - break; - case ( kXSecTwkDial_ZExpELFF_AN2 ) : - fZExpParaTwkDial.fZ_ANn[1] = twk_dial; - fIsSinglePara = true; - break; - case ( kXSecTwkDial_ZExpELFF_AN3 ) : - fZExpParaTwkDial.fZ_ANn[2] = twk_dial; - fIsSinglePara = true; - break; - case ( kXSecTwkDial_ZExpELFF_AN4 ) : - fZExpParaTwkDial.fZ_ANn[3] = twk_dial; - fIsSinglePara = true; - break; - case ( kXSecTwkDial_ZExpELFF_BP1 ) : - fZExpParaTwkDial.fZ_BPn[0] = twk_dial; - fIsSinglePara = true; - break; - case ( kXSecTwkDial_ZExpELFF_BP2 ) : - fZExpParaTwkDial.fZ_BPn[1] = twk_dial; - fIsSinglePara = true; - break; - case ( kXSecTwkDial_ZExpELFF_BP3 ) : - fZExpParaTwkDial.fZ_BPn[2] = twk_dial; - fIsSinglePara = true; - break; - case ( kXSecTwkDial_ZExpELFF_BP4 ) : - fZExpParaTwkDial.fZ_BPn[3] = twk_dial; - fIsSinglePara = true; - break; - case ( kXSecTwkDial_ZExpELFF_BN1 ) : - fZExpParaTwkDial.fZ_BNn[0] = twk_dial; - fIsSinglePara = true; - break; - case ( kXSecTwkDial_ZExpELFF_BN2 ) : - fZExpParaTwkDial.fZ_BNn[1] = twk_dial; - fIsSinglePara = true; - break; - case ( kXSecTwkDial_ZExpELFF_BN3 ) : - fZExpParaTwkDial.fZ_BNn[2] = twk_dial; - fIsSinglePara = true; - break; - case ( kXSecTwkDial_ZExpELFF_BN4 ) : - fZExpParaTwkDial.fZ_BNn[3] = twk_dial; - fIsSinglePara = true; break; default: break; @@ -238,92 +150,13 @@ void GReWeightNuXSecCCQEELFF::Reset(void) //_______________________________________________________________________________________ void GReWeightNuXSecCCQEELFF::Reconfigure(void) { - if(fIsSinglePara && fIsAllPara){ - LOG("ReW",pERROR) << "can't not set kXSecTwkDial_ZExpELFF and kXSecTwkDial_ZExpELFFXXX at the same time"; - abort(); - } GSystUncertainty * fracerr = GSystUncertainty::Instance(); if(fMode==kModeZExp && fModelIsZExp) { - if(fIsSinglePara){ - int sign_twk = 0; - double fracerr_zexp = 0.; - GSyst_t syst; - // loop over all indices and update each - for (int i=0;iOneSigmaErr(syst, sign_twk); - fZExpPara.fZ_APn[i] = fZExpParaDef.fZ_APn[i] * ( 1.0 + fZExpParaTwkDial.fZ_APn[i] * fracerr_zexp); - switch(i){ - case 0: syst = kXSecTwkDial_ZExpELFF_BP1; break; - case 1: syst = kXSecTwkDial_ZExpELFF_BP2; break; - case 2: syst = kXSecTwkDial_ZExpELFF_BP3; break; - case 3: syst = kXSecTwkDial_ZExpELFF_BP4; break; - default: return; break; - } - sign_twk = utils::rew::Sign(fZExpParaTwkDial.fZ_BPn[i]); - fracerr_zexp = fracerr->OneSigmaErr(syst, sign_twk); - fZExpPara.fZ_BPn[i] = fZExpParaDef.fZ_BPn[i] * (1.0 + fZExpParaTwkDial.fZ_BPn[i] * fracerr_zexp); - switch(i){ - case 0: syst = kXSecTwkDial_ZExpELFF_AN1; break; - case 1: syst = kXSecTwkDial_ZExpELFF_AN2; break; - case 2: syst = kXSecTwkDial_ZExpELFF_AN3; break; - case 3: syst = kXSecTwkDial_ZExpELFF_AN4; break; - default: return; break; - } - sign_twk = utils::rew::Sign(fZExpParaTwkDial.fZ_ANn[i]); - fracerr_zexp = fracerr->OneSigmaErr(syst, sign_twk); - fZExpPara.fZ_ANn[i] = fZExpParaDef.fZ_ANn[i] * ( 1.0 + fZExpParaTwkDial.fZ_ANn[i] * fracerr_zexp); - switch(i){ - case 0: syst = kXSecTwkDial_ZExpELFF_BN1; break; - case 1: syst = kXSecTwkDial_ZExpELFF_BN2; break; - case 2: syst = kXSecTwkDial_ZExpELFF_BN3; break; - case 3: syst = kXSecTwkDial_ZExpELFF_BN4; break; - default: return; break; - } - sign_twk = utils::rew::Sign(fZExpParaTwkDial.fZ_BNn[i]); - fracerr_zexp = fracerr->OneSigmaErr(syst, sign_twk); - fZExpPara.fZ_BNn[i] = fZExpParaDef.fZ_BNn[i] * ( 1.0 + fZExpParaTwkDial.fZ_BNn[i] * fracerr_zexp); - } - - - Registry r("GReWeightNuXSecCCQEELFF",false); - //~ Registry r(fXSecModel->GetConfig()); - if (fMode==kModeZExp) - { - ostringstream alg_key; - for (int i=0;iConfigure(r); - } - else if(fIsAllPara){ - int sign_twk = 0; - double fracerr_zexp = 0.; - sign_twk = utils::rew::Sign(fZExpTwkDial); - fracerr_zexp = fracerr->OneSigmaErr(kXSecTwkDial_ZExpELFF, sign_twk); - fZExp_Scale = fZExpTwkDial * fracerr_zexp; - } + int sign_twk = 0; + double fracerr_zexp = 0.; + sign_twk = utils::rew::Sign(fZExpTwkDial); + fracerr_zexp = fracerr->OneSigmaErr(kXSecTwkDial_ZExpELFF, sign_twk); + fZExp_Scale = fZExpTwkDial * fracerr_zexp; } else { return; @@ -355,6 +188,7 @@ double GReWeightNuXSecCCQEELFF::CalcWeight(const genie::EventRecord & event) if ( nupdg==kPdgAntiNuE && !fRewNuebar ) return 1.; double wght = 1.0; + if ( fMode==kModeZExp && fModelIsZExp ) { wght *= this->CalcWeightZExp( event ); return wght; @@ -366,11 +200,14 @@ void GReWeightNuXSecCCQEELFF::Init(void) { AlgConfigPool * conf_pool = AlgConfigPool::Instance(); Registry * gpl = conf_pool->GlobalParameterList(); + + // --- get the model configuration from current tune conf_pool->Print(std::cout); gpl->Print(std::cout); - Registry * rgcovmat = conf_pool->FindRegistry("genie::rew::GSystUncertaintyTable/CovarianceMatrix"); - int n_row = rgcovmat->GetInt(Algorithm::BuildParamMatRowSizeKey("ZExpELFF@CovarianceMatrix")); - int n_col = rgcovmat->GetInt(Algorithm::BuildParamMatColSizeKey("ZExpELFF@CovarianceMatrix")); + RgAlg elff_id = conf_pool->FindRegistry("CommonParamList/ElasticFF")->GetAlg("ElasticFormFactorsModel"); + Registry * elff_model = conf_pool->FindRegistry(elff_id); + int n_row = elff_model->GetInt(Algorithm::BuildParamMatRowSizeKey("ZExpELFF@CovarianceMatrix")); + int n_col = elff_model->GetInt(Algorithm::BuildParamMatColSizeKey("ZExpELFF@CovarianceMatrix")); if(n_row != n_col){ LOG( "GReWeightNuXSecCCQEELFF", pFATAL ) << "Non-square covariance matrix" << "encountered in GReWeightNuXSecCCQEELFF::Init()"; @@ -381,13 +218,11 @@ void GReWeightNuXSecCCQEELFF::Init(void) A_f.resize(n_row); for(int i = 0; i < n_row; i++){ for(int j = 0; j < n_row; j++){ - error_mat[i][j] = rgcovmat->GetDouble(Algorithm::BuildParamMatKey("ZExpELFF@CovarianceMatrix", i, j)); + error_mat[i][j] = elff_model->GetDouble(Algorithm::BuildParamMatKey("ZExpELFF@CovarianceMatrix", i, j)); } } - rgcovmat->Print(std::cout); + error_mat.Print(); RgAlg xsec_alg = gpl->GetAlg("XSecModel@genie::EventGenerator/QEL-CC"); - - AlgId id(xsec_alg); AlgId twk_id(id); @@ -414,8 +249,6 @@ void GReWeightNuXSecCCQEELFF::Init(void) fModelIsZExp = (strcmp(fFFModel.c_str(),kModelZExp ) == 0); - fIsSinglePara = false; - fIsAllPara = false; this->RewNue (true); this->RewNuebar (true); @@ -481,60 +314,12 @@ void GReWeightNuXSecCCQEELFF::Init(void) fZExpPara.fZ_BNn[i] = 0.; fZExpPara.fZ_BPn[i] = 0.; } + fZExpTwkDial = 0.; } } //_______________________________________________________________________________________ double GReWeightNuXSecCCQEELFF::CalcWeightZExp(const genie::EventRecord & event) { - if(fIsSinglePara){ - bool tweaked = false; - for (int i=0;i controls::kASmallNum); - tweaked = tweaked || (TMath::Abs(fZExpParaTwkDial.fZ_ANn[i]) > controls::kASmallNum); - tweaked = tweaked || (TMath::Abs(fZExpParaTwkDial.fZ_BPn[i]) > controls::kASmallNum); - tweaked = tweaked || (TMath::Abs(fZExpParaTwkDial.fZ_BNn[i]) > controls::kASmallNum); - } - if(!tweaked) { return 1.0; } - - Interaction * interaction = event.Summary(); - - interaction->KinePtr()->UseSelectedKinematics(); - - // Retrieve the kinematic phase space used to generate the event - const KinePhaseSpace_t phase_space = event.DiffXSecVars(); - double old_xsec = event.DiffXSec(); - - if (phase_space == kPSQ2fE) { - interaction->SetBit(kIAssumeFreeNucleon); - } - - if (!fUseOldWeightFromFile || fNWeightChecksDone < fNWeightChecksToDo) { - double calc_old_xsec = fXSecModelDef->XSec(interaction, phase_space); - if (fNWeightChecksDone < fNWeightChecksToDo) { - if (std::abs(calc_old_xsec - old_xsec)/old_xsec > controls::kASmallNum) { - LOG("ReW",pWARN) << "Warning - default dxsec does not match dxsec saved in tree. Does the config match?"; - } - fNWeightChecksDone++; - } - if(!fUseOldWeightFromFile) { - old_xsec = calc_old_xsec; - } - } - double old_weight = event.Weight(); - double new_xsec = fXSecModel->XSec(interaction, phase_space); - double new_weight = old_weight * (new_xsec/old_xsec); - - - interaction->KinePtr()->ClearRunningValues(); - - if (phase_space == kPSQ2fE) { - interaction->ResetBit(kIAssumeFreeNucleon); - } - - return new_weight; - } - else{ // uncertainty propagation bool tweaked = false; tweaked = tweaked || (TMath::Abs(fZExpTwkDial) > controls::kASmallNum); if(!tweaked) { return 1.0; } @@ -543,7 +328,6 @@ double GReWeightNuXSecCCQEELFF::CalcWeightZExp(const genie::EventRecord & event) double new_weight = (fZExp_Scale * oneSigma + old_xsec) / old_xsec; // LOG("ReW", pNOTICE) << fZExp_Scale << " " << oneSigma << " " << old_xsec; return new_weight; - } } //_______________________________________________________________________________________ diff --git a/src/RwCalculators/GReWeightNuXSecCCQEZAFF.cxx b/src/RwCalculators/GReWeightNuXSecCCQEZAFF.cxx new file mode 100644 index 00000000..aa56fb5a --- /dev/null +++ b/src/RwCalculators/GReWeightNuXSecCCQEZAFF.cxx @@ -0,0 +1,419 @@ +//____________________________________________________________________________ +/* + Copyright (c) 2003-2025, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org + + Liang Liu + Fermi National Accelerator Laboratory + +*/ +//____________________________________________________________________________ + +#include +#include +#include +#include +#include + +// GENIE/Generator includes +#include "Framework/Algorithm/AlgFactory.h" +#include "Framework/Algorithm/AlgConfigPool.h" +#include "Framework/EventGen/XSecAlgorithmI.h" +#include "Framework/Conventions/Units.h" +#include "Framework/Conventions/Controls.h" +#include "Framework/EventGen/EventRecord.h" +#include "Framework/GHEP/GHepParticle.h" +#include "Framework/Interaction/Interaction.h" +#include "Framework/Messenger/Messenger.h" +#include "Framework/ParticleData/PDGCodes.h" +#include "Framework/Registry/Registry.h" +#include "Physics/XSectionIntegration/XSecIntegratorI.h" + +// GENIE/Reweight includes +#include "GReWeightNuXSecCCQEZAFF.h" +#include "RwCalculators/GReWeightUtils.h" +#include "RwFramework/GSystSet.h" +#include "RwFramework/GSystUncertainty.h" +#include "TRandom3.h" +#include "TVectorD.h" +#include "TDecompChol.h" + +using namespace genie; +using namespace genie::rew; +using std::ostringstream; + +static const char* kModelZExp = "genie::ZExpAxialFormFactorModel"; + +const int GReWeightNuXSecCCQEZAFF::kModeZExp; + +//_______________________________________________________________________________________ +GReWeightNuXSecCCQEZAFF::GReWeightNuXSecCCQEZAFF() : + GReWeightModel("CCQE"), + fManualModelName(), + fManualModelType() +{ + this->Init(); +} +//_______________________________________________________________________________________ +GReWeightNuXSecCCQEZAFF::GReWeightNuXSecCCQEZAFF(std::string model, std::string type) : + GReWeightModel("CCQE"), + fManualModelName(model), + fManualModelType(type) +{ + this->Init(); +} +//_______________________________________________________________________________________ +GReWeightNuXSecCCQEZAFF::~GReWeightNuXSecCCQEZAFF() +{ + if ( fXSecModelConfig ) delete fXSecModelConfig; + if ( fXSecModel ) delete fXSecModel; + if ( fXSecModelDef ) delete fXSecModelDef; +} +//_______________________________________________________________________________________ +bool GReWeightNuXSecCCQEZAFF::IsHandled(GSyst_t syst) const +{ + // read form factor model and compare to mode + bool handle; + + switch(syst) { + case ( kXSecTwkDial_ZExpZAFF ) : + if(fMode==kModeZExp && fModelIsZExp){ + handle = true; + }else { + handle = false; + } + break; + default: + handle = false; + break; + } + return handle; +} +//_______________________________________________________________________________________ +bool GReWeightNuXSecCCQEZAFF::AppliesTo(const EventRecord &event) const +{ + auto type = event.Summary()->ProcInfo().ScatteringTypeId(); + bool is_cc = event.Summary()->ProcInfo().IsWeakCC(); + if (type==kScQuasiElastic && is_cc) { + return true; + } + return false; +} +//_______________________________________________________________________________________ +void GReWeightNuXSecCCQEZAFF::SetSystematic(GSyst_t syst, double twk_dial) +{ + if(!this->IsHandled(syst)) + { + LOG("ReW",pWARN) << "Systematic " << GSyst::AsString(syst) << " is not handled for algorithm " + << fFFModel << " and mode " << fMode; + return; + } + switch(syst) { + case (kXSecTwkDial_ZExpZAFF): + fZExpTwkDial = twk_dial; + break; + default: + break; + } +} +//_______________________________________________________________________________________ +void GReWeightNuXSecCCQEZAFF::Reset(void) +{ + fZExpPara.fQ4limit = fZExpParaDef.fQ4limit; + fZExpPara.fKmax = fZExpParaDef.fKmax; + fZExpPara.fT0 = fZExpParaDef.fT0; + fZExpParaTwkDial.fT0 = 0.; + fZExpPara.fTcut = fZExpParaDef.fTcut; + fZExpParaTwkDial.fTcut = 0.; + fZExpPara.fGep0 = fZExpParaDef.fGep0; + fZExpParaTwkDial.fGep0 = 0.; + fZExpPara.fGmp0 = fZExpParaDef.fGmp0; + fZExpParaTwkDial.fGmp0 = 0.; + fZExpPara.fGen0 = fZExpParaDef.fGen0; + fZExpParaTwkDial.fGen0 = 0.; + fZExpPara.fGmn0 = fZExpParaDef.fGmn0; + fZExpParaTwkDial.fGmn0 = 0.; + for(int i = 0; i < fZExpPara.fKmax; i++){ + fZExpPara.fZ_An[i] = fZExpParaDef.fZ_An[i]; + fZExpParaTwkDial.fZ_An[i] = 0.; + } + this->Reconfigure(); +} +//_______________________________________________________________________________________ +void GReWeightNuXSecCCQEZAFF::Reconfigure(void) +{ + GSystUncertainty * fracerr = GSystUncertainty::Instance(); + if(fMode==kModeZExp && fModelIsZExp) { + int sign_twk = 0; + double fracerr_zexp = 0.; + sign_twk = utils::rew::Sign(fZExpTwkDial); + fracerr_zexp = fracerr->OneSigmaErr(kXSecTwkDial_ZExpZAFF, sign_twk); + fZExp_Scale = fZExpTwkDial * fracerr_zexp; + } + else { + return; + } +} +//_______________________________________________________________________________________ +double GReWeightNuXSecCCQEZAFF::CalcWeight(const genie::EventRecord & event) +{ + bool is_qe = event.Summary()->ProcInfo().IsQuasiElastic(); + bool is_cc = event.Summary()->ProcInfo().IsWeakCC(); + if ( !is_qe || !is_cc ) return 1.; + + bool charm = event.Summary()->ExclTag().IsCharmEvent(); // skip CCQE charm + if ( charm ) return 1.; + + // Skip CCQE strange + bool strange = event.Summary()->ExclTag().IsStrangeEvent(); + if ( strange ) return 1.; + + // Skip any other CCQE channels that do not produce a final-state nucleon + int final_nucleon_pdgc = event.Summary()->RecoilNucleonPdg(); + if ( !pdg::IsNucleon(final_nucleon_pdgc) ) return 1.; + + int nupdg = event.Probe()->Pdg(); + + if ( nupdg==kPdgNuMu && !fRewNumu ) return 1.; + if ( nupdg==kPdgAntiNuMu && !fRewNumubar) return 1.; + if ( nupdg==kPdgNuE && !fRewNue ) return 1.; + if ( nupdg==kPdgAntiNuE && !fRewNuebar ) return 1.; + + double wght = 1.0; + if ( fMode==kModeZExp && fModelIsZExp ) { + wght *= this->CalcWeightZExp( event ); + return wght; + } + return 1.; +} +//_______________________________________________________________________________________ +void GReWeightNuXSecCCQEZAFF::Init(void) +{ + + // Get the model and parameters of axial form factor from current tune + AlgConfigPool * conf_pool = AlgConfigPool::Instance(); + Registry * gpl = conf_pool->GlobalParameterList(); + // get axial form factor tune from current CCQE model + RgAlg cc_qel_id = gpl->GetAlg( "XSecModel@genie::EventGenerator/QEL-CC" ); + RgAlg ff_id = conf_pool->FindRegistry(cc_qel_id)->GetAlg("FormFactorsAlg"); + RgAlg aff_model_id = conf_pool->FindRegistry(ff_id)->GetAlg("AxialFormFactorModel"); + // get the covariance matrix + Registry * zexp_axial_model = conf_pool->FindRegistry(aff_model_id); + int n_row = zexp_axial_model->GetInt(Algorithm::BuildParamMatRowSizeKey("ZExpZAFF@CovarianceMatrix")); + int n_col = zexp_axial_model->GetInt(Algorithm::BuildParamMatColSizeKey("ZExpZAFF@CovarianceMatrix")); + + if(n_row != n_col){ + LOG( "GReWeightNuXSecCCQEZAFF", pFATAL ) << "Non-square covariance matrix" + << "encountered in GReWeightNuXSecCCQEZAFF::Init()"; + std::exit(1); + } + error_mat.ResizeTo(n_row, n_row); + errors.resize(n_row); + A_f.resize(n_row); + for(int i = 0; i < n_row; i++){ + for(int j = 0; j < n_row; j++){ + error_mat[i][j] = zexp_axial_model->GetDouble(Algorithm::BuildParamMatKey("ZExpZAFF@CovarianceMatrix", i, j)); + } + } + LOG( "GReWeightNuXSecCCQEZAFF", pINFO ) << "CovarianceMatrix of z expansion:"; + error_mat.Print(); + + AlgId id(cc_qel_id); + + AlgId twk_id(id); + if (fManualModelName.size()) { + twk_id = AlgId(fManualModelName,fManualModelType); + } + + AlgFactory * algf = AlgFactory::Instance(); + + Algorithm * alg_def = algf->AdoptAlgorithm(id); + fXSecModelDef = dynamic_cast(alg_def); + fXSecModelDef->AdoptSubstructure(); + + Algorithm * alg_twk = algf->AdoptAlgorithm(twk_id); + fXSecModel = dynamic_cast(alg_twk); + fXSecModel->AdoptSubstructure(); + + + // Check what kind of form factors we're using in the tweaked cross section + // model + fXSecModelConfig = new Registry(fXSecModel->GetConfig()); + fFFModel = fXSecModelConfig->GetAlg("FormFactorsAlg/AxialFormFactorModel").name; + fXSecModelConfig->Print(std::cout); + + fModelIsZExp = (strcmp(fFFModel.c_str(), kModelZExp ) == 0); + + + this->RewNue (true); + this->RewNuebar (true); + this->RewNumu (true); + this->RewNumubar(true); + + this->SetZExpPath("FormFactorsAlg/AxialFormFactorModel/"); + + if (fModelIsZExp) + { + this->SetMode(kModeZExp); + fZExpParaDef.fQ4limit = fXSecModelConfig->GetBool(fZExpPath + "QEL-Q4limit"); + fZExpParaDef.fKmax = fXSecModelConfig->GetInt(fZExpPath + "QEL-Kmax"); + fZExpParaDef.fT0 = fXSecModelConfig->GetDouble(fZExpPath + "QEL-T0"); + fZExpParaDef.fTcut = fXSecModelConfig->GetDouble(fZExpPath + "QEL-Tcut"); + fZExpParaDef.fGep0 = fXSecModelConfig->GetDouble(fZExpPath + "QEL-Gep0"); + fZExpParaDef.fGmp0 = fXSecModelConfig->GetDouble(fZExpPath + "QEL-Gmp0"); + fZExpParaDef.fGen0 = fXSecModelConfig->GetDouble(fZExpPath + "QEL-Gen0"); + fZExpParaDef.fGmn0 = fXSecModelConfig->GetDouble(fZExpPath + "QEL-Gmn0"); + + fZExpParaDef.fZ_An.resize(fZExpParaDef.fKmax); + fZExpPara.fZ_An.resize(fZExpParaDef.fKmax); + fZExpParaTwkDial.fZ_An.resize(fZExpParaDef.fKmax); + + fZExpParaTwkDial.fQ4limit = 0.; + fZExpParaTwkDial.fKmax = 0.; + fZExpParaTwkDial.fT0 = 0.; + fZExpParaTwkDial.fTcut = 0.; + fZExpParaTwkDial.fGep0 = 0.; + fZExpParaTwkDial.fGmp0 = 0.; + fZExpParaTwkDial.fGen0 = 0.; + fZExpParaTwkDial.fGmn0 = 0.; + + ostringstream alg_key; + for(int i = 0; i < fZExpParaDef.fKmax; i++){ + alg_key.str(""); + alg_key << fZExpPath << "QEL-Z_A-" << i; + fZExpParaDef.fZ_An[i] = fXSecModelConfig->GetDouble(alg_key.str()); + fZExpParaTwkDial.fZ_An[i] = 0.; + fZExpPara.fZ_An[i] = 0.; + } + fZExpTwkDial = 0.; + } +} +//_______________________________________________________________________________________ +double GReWeightNuXSecCCQEZAFF::CalcWeightZExp(const genie::EventRecord & event) +{ + bool tweaked = false; + tweaked = tweaked || (TMath::Abs(fZExpTwkDial) > controls::kASmallNum); + if(!tweaked) { return 1.0; } + double oneSigma = GetOneSigma(event); + double old_xsec = event.DiffXSec(); + double new_weight = (fZExp_Scale * oneSigma + old_xsec) / old_xsec; + return new_weight; +} + +//_______________________________________________________________________________________ +// Calculate the partial derivative of the cross section +// A_f is the partial derivative +// A_f = \patial XSec() / \patial p_i = (XSec(p_i + \delta * error_i ) - XSec(p_i - \delta * error_i )) / ( 2.0 * \delta * error_i ) +// + +void GReWeightNuXSecCCQEZAFF::XSecPartialDerivative(const EventRecord & event){ + // Get the uncertainties from the error matrix + // ap1, ap2, ap3, ap4, + // bp1, bp2, bp3, bp4, + // an1, an2, an3, an4, + // bn1, bn2, bn3, bn4 + + for(int i = 0; i < fZExpParaDef.fKmax; i++){ + errors[i] = TMath::Sqrt(error_mat[i][i]); + A_f[i] = 0.0; + } + + double delta = 0.1; + + for(int index = 0; index < fZExpParaDef.fKmax; index++){ + double xsec_tmp_0 = 0.0; + double xsec_tmp_1 = 0.0; + + for(int sign = -1; sign < 2; sign++){ + if(sign == 0) continue; + + for(int ipara = 0; ipara < fZExpParaDef.fKmax; ipara++){ + fZExpPara.fZ_An[ipara] = fZExpParaDef.fZ_An[ipara]; + } + + fZExpPara.fZ_An[index] = fZExpParaDef.fZ_An[index] + errors[index] * delta * sign; + + Registry r("GReWeightNuXSecCCQEZAFF",false); + //~ Registry r(fXSecModel->GetConfig()); + if (fMode==kModeZExp) + { + ostringstream alg_key; + for (int i=0;iConfigure(r); + if(sign == -1){ + xsec_tmp_0 = UpdateXSec(event) ; + } + else if(sign == 1){ + xsec_tmp_1 = UpdateXSec(event) ; + } + } + + double delta_xsec = xsec_tmp_1 - xsec_tmp_0; + A_f[index] = delta_xsec/(errors[index] * delta * 2.0); + // LOG("ReW", pNOTICE) << A_f[index] ; + } +} + +//_______________________________________________________________________________________ +// Uncertainty propagation +// +// \sigma_{XSec}^2 = A_f[i] *A_f[j] *M_ij +double GReWeightNuXSecCCQEZAFF::GetOneSigma(const EventRecord & event){ + XSecPartialDerivative(event); + double OneSigma2 = 0; + for(int i = 0; i < fZExpParaDef.fKmax; i++){ + for(int j = 0; j < fZExpParaDef.fKmax; j++){ + OneSigma2+= A_f[i]*A_f[j]*error_mat[i][j]; + } + } + return TMath::Sqrt(OneSigma2); +} + +//_______________________________________________________________________________________ + + +double GReWeightNuXSecCCQEZAFF::UpdateXSec(const EventRecord & event){ + + Interaction * interaction = event.Summary(); + + interaction->KinePtr()->UseSelectedKinematics(); + + // Retrieve the kinematic phase space used to generate the event + const KinePhaseSpace_t phase_space = event.DiffXSecVars(); + double old_xsec = event.DiffXSec(); + + if (phase_space == kPSQ2fE) { + interaction->SetBit(kIAssumeFreeNucleon); + } + + if (!fUseOldWeightFromFile || fNWeightChecksDone < fNWeightChecksToDo) { + double calc_old_xsec = fXSecModelDef->XSec(interaction, phase_space); + if (fNWeightChecksDone < fNWeightChecksToDo) { + if (std::abs(calc_old_xsec - old_xsec)/old_xsec > controls::kASmallNum) { + LOG("ReW",pWARN) << "Warning - default dxsec does not match dxsec saved in tree. Does the config match?"; + } + fNWeightChecksDone++; + } + if(!fUseOldWeightFromFile) { + old_xsec = calc_old_xsec; + } + } + double new_xsec = fXSecModel->XSec(interaction, phase_space); + + interaction->KinePtr()->ClearRunningValues(); + + if (phase_space == kPSQ2fE) { + interaction->ResetBit(kIAssumeFreeNucleon); + } + + return new_xsec; +} + + diff --git a/src/RwCalculators/GReWeightNuXSecCCQEZAFF.h b/src/RwCalculators/GReWeightNuXSecCCQEZAFF.h new file mode 100644 index 00000000..52eacf0e --- /dev/null +++ b/src/RwCalculators/GReWeightNuXSecCCQEZAFF.h @@ -0,0 +1,130 @@ +//____________________________________________________________________________ +/*! + + \class genie::rew::GReWeightNuXSecCCQEZAFF + + \brief Reweighting CCQE GENIE neutrino cross sections + Z expansion axial form factor model + + \author Liang Liu + Fermi National Accelerator Laboratory + + \created March 25, 2024 + + \cpright Copyright (c) 2003-2025, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org + */ +//____________________________________________________________________________ + +#ifndef _G_REWEIGHT_NU_XSEC_CCQE_ZAFF_H_ +#define _G_REWEIGHT_NU_XSEC_CCQE_ZAFF_H_ + +#include + +// GENIE/Reweight includes +#include "RwCalculators/GReWeightModel.h" +#include "TRandom3.h" +#include "TMatrixDSym.h" + +class TFile; +class TNtupleD; + +namespace genie { + + class XSecAlgorithmI; + class XSecIntegratorI; + class Registry; + + namespace rew { + + class GReWeightNuXSecCCQEZAFF : public GReWeightModel + { + public: + static const int kModeZExp = 0; + + GReWeightNuXSecCCQEZAFF(); + GReWeightNuXSecCCQEZAFF(std::string model, std::string type); + ~GReWeightNuXSecCCQEZAFF(); + + // implement the GReWeightI interface + bool AppliesTo (const EventRecord &event) const; + bool IsHandled (GSyst_t syst) const; + void SetSystematic (GSyst_t syst, double val); + void Reset (void); + void Reconfigure (void); + double CalcWeight (const EventRecord & event); + + // various config options + void SetMode (int mode) { fMode = mode; } + void RewNue (bool tf ) { fRewNue = tf; } + void RewNuebar (bool tf ) { fRewNuebar = tf; } + void RewNumu (bool tf ) { fRewNumu = tf; } + void RewNumubar (bool tf ) { fRewNumubar = tf; } + // z-expansion specific options + void SetZExpPath (string p){ fZExpPath = p; } + + private: + void Init (void); + double CalcWeightZExp (const EventRecord & event); + + double GetOneSigma(const EventRecord & event); + double UpdateXSec(const EventRecord & event); + void XSecPartialDerivative(const EventRecord & event); + + XSecAlgorithmI * fXSecModelDef; ///< default model + XSecAlgorithmI * fXSecModel; ///< tweaked model + Registry * fXSecModelConfig; ///< config in tweaked model + string fFFModel; ///< String name of form factor model + bool fModelIsZExp; ///< Using Zexp form factors? + + + std::string fManualModelName; ///< If using a tweaked model that isn't the same as default, name + std::string fManualModelType; ///< If using a tweaked model that isn't the same as default, type + + int fMode; ///< 0: ZExp; TODO: currently there is only one model implemented. + bool fRewNue; ///< reweight nu_e CC? + bool fRewNuebar; ///< reweight nu_e_bar CC? + bool fRewNumu; ///< reweight nu_mu CC? + bool fRewNumubar; ///< reweight nu_mu_bar CC? + + // unused // int fZExpCurrIdx; ///< current coefficient index + int fZExpMaxCoef; ///< max number of coefficients to use + string fZExpPath; ///< algorithm path to get coefficients + + // parameters in z expansion vector model + // define the default, current and tweak dial + struct ZExpPara { + bool fQ4limit; + int fKmax; + double fT0; + double fTcut; + double fGep0; + double fGmp0; + double fGen0; + double fGmn0; + std::vector fZ_An; + } fZExpParaDef, fZExpPara, fZExpParaTwkDial; + // tweek dial and scale factor in propagation method + double fZExpTwkDial; + double fZExp_Scale; + + // Two methods are provided to calculate the uncertainties of XSec + // 1. propagation of errors: it is based on grwght1p + // 2. Cholesky decomposition: it is based on grwghtnp + bool fIsSinglePara; // it will be used for Cholesky decomposition + bool fIsAllPara; // flag of propagation method + std::vector A_f; + + // List of the uncertainties of parameters from Kaushik + // ap1, ap2, ap3, ap4, + // bp1, bp2, bp3, bp4, + // an1, an2, an3, an4, + // bn1, bn2, bn3, bn4 + std::vector errors; + TMatrixDSym error_mat; + }; + + } // rew namespace +} // genie namespace + +#endif diff --git a/src/RwCalculators/LinkDef.h b/src/RwCalculators/LinkDef.h index 124683af..b32d26ea 100644 --- a/src/RwCalculators/LinkDef.h +++ b/src/RwCalculators/LinkDef.h @@ -21,6 +21,8 @@ #pragma link C++ class genie::rew::GReWeightDISNuclMod; #pragma link C++ class genie::rew::GReWeightNuXSecNCEL; #pragma link C++ class genie::rew::GReWeightNuXSecCCQE; +#pragma link C++ class genie::rew::GReWeightNuXSecCCQEELFF; +#pragma link C++ class genie::rew::GReWeightNuXSecCCQEZAFF; #pragma link C++ class genie::rew::GReWeightNuXSecCCQEvec; #pragma link C++ class genie::rew::GReWeightNuXSecCCQEaxial; #pragma link C++ class genie::rew::GReWeightNuXSecCCRES; diff --git a/src/RwFramework/GSyst.cxx b/src/RwFramework/GSyst.cxx index d2cfd0cf..567524a0 100644 --- a/src/RwFramework/GSyst.cxx +++ b/src/RwFramework/GSyst.cxx @@ -122,22 +122,7 @@ std::map GSyst::BuildGSystToStringMap() { temp_map[ kXSecTwkDial_NormCCCOHpi ] = "NormCCCOH"; temp_map[ kXSecTwkDial_NormNCCOHpi ] = "NormNCCOH"; temp_map[ kXSecTwkDial_ZExpELFF ] = "ZExpELFFCCQE"; - temp_map[ kXSecTwkDial_ZExpELFF_AP1 ] = "ZExpELFFAP1CCQE"; - temp_map[ kXSecTwkDial_ZExpELFF_AP2 ] = "ZExpELFFAP2CCQE"; - temp_map[ kXSecTwkDial_ZExpELFF_AP3 ] = "ZExpELFFAP3CCQE"; - temp_map[ kXSecTwkDial_ZExpELFF_AP4 ] = "ZExpELFFAP4CCQE"; - temp_map[ kXSecTwkDial_ZExpELFF_AN1 ] = "ZExpELFFAN1CCQE"; - temp_map[ kXSecTwkDial_ZExpELFF_AN2 ] = "ZExpELFFAN2CCQE"; - temp_map[ kXSecTwkDial_ZExpELFF_AN3 ] = "ZExpELFFAN3CCQE"; - temp_map[ kXSecTwkDial_ZExpELFF_AN4 ] = "ZExpELFFAN4CCQE"; - temp_map[ kXSecTwkDial_ZExpELFF_BP1 ] = "ZExpELFFBP1CCQE"; - temp_map[ kXSecTwkDial_ZExpELFF_BP2 ] = "ZExpELFFBP2CCQE"; - temp_map[ kXSecTwkDial_ZExpELFF_BP3 ] = "ZExpELFFBP3CCQE"; - temp_map[ kXSecTwkDial_ZExpELFF_BP4 ] = "ZExpELFFBP4CCQE"; - temp_map[ kXSecTwkDial_ZExpELFF_BN1 ] = "ZExpELFFBN1CCQE"; - temp_map[ kXSecTwkDial_ZExpELFF_BN2 ] = "ZExpELFFBN2CCQE"; - temp_map[ kXSecTwkDial_ZExpELFF_BN3 ] = "ZExpELFFBN3CCQE"; - temp_map[ kXSecTwkDial_ZExpELFF_BN4 ] = "ZExpELFFBN4CCQE"; + temp_map[ kXSecTwkDial_ZExpZAFF ] = "ZExpZAFFCCQE"; return temp_map; } diff --git a/src/RwFramework/GSyst.h b/src/RwFramework/GSyst.h index bec7ddda..de9cd342 100644 --- a/src/RwFramework/GSyst.h +++ b/src/RwFramework/GSyst.h @@ -223,22 +223,10 @@ typedef enum EGSyst { // Alternative approach to CCQE form factors (z-expansion) vector form factor // kXSecTwkDial_ZExpELFF, - kXSecTwkDial_ZExpELFF_AP1, - kXSecTwkDial_ZExpELFF_AP2, - kXSecTwkDial_ZExpELFF_AP3, - kXSecTwkDial_ZExpELFF_AP4, - kXSecTwkDial_ZExpELFF_AN1, - kXSecTwkDial_ZExpELFF_AN2, - kXSecTwkDial_ZExpELFF_AN3, - kXSecTwkDial_ZExpELFF_AN4, - kXSecTwkDial_ZExpELFF_BP1, - kXSecTwkDial_ZExpELFF_BP2, - kXSecTwkDial_ZExpELFF_BP3, - kXSecTwkDial_ZExpELFF_BP4, - kXSecTwkDial_ZExpELFF_BN1, - kXSecTwkDial_ZExpELFF_BN2, - kXSecTwkDial_ZExpELFF_BN3, - kXSecTwkDial_ZExpELFF_BN4, + // + // Alternative approach to CCQE form factors (z-expansion) axial form factor + // + kXSecTwkDial_ZExpZAFF, // // Misc