diff --git a/config/GTEST23_10k/CommonParam.xml b/config/GTEST23_10k/CommonParam.xml new file mode 100644 index 000000000..ce485ff9a --- /dev/null +++ b/config/GTEST23_10k/CommonParam.xml @@ -0,0 +1,366 @@ + + + + + + + + + + + + + 2.30 + 3.00 + + + + + 1.065047 + 1.019031 + 0.877255 + + + + + 0.501716712132 + + + + + 0.804 + 0.463 + + + + + + + 0.97417 + 0.2248 + 0.220 + 0.995 + + 0.227780466682 + + + + + + + 1.4 + 3.0 + + + + + 1.430 + + + + + + true + 1.927862 + + + 0.007713 + 0.787999 + 0.100 + 1.000 + 0.100 + 1.000 + 0.127858 + 2.115230 + 0.300 + 1.000 + 0.300 + 1.000 + 0.127858 + 2.115230 + 0.300 + 1.000 + 0.300 + 1.000 + 0.007713 + 0.787999 + 0.100 + 1.000 + 0.100 + 1.000 + + + + + + + 0.5000 + 0.0170 + 0.0250 + 0.0270 + 0.0320 + 0.0295 + 0.0280 + 0.0300 + 0.0360 + 0.0360 + 0.0440 + false + + + Default + + + + + + + + 0.00 + 0.05 + 1.0 + 1.0 + 1.0 + 0.041 + 1.0 + 0.250 + true + true + true + + + + 0.5 + 0.7 + + + + 0.0 + 0.0 + true + + + + + 1.000 + + + 0.000 + 1.000 + + + 0.250 + + + + + + + + 7.0 + + + 0.35 + + + + + + -1.2670 + + + 1 + + + 2.793 + + + 0 + + + -1.913042 + + + 0.961242 + 0.840 + + + + + + 2.7930 + -1.913042 + + + + + + 1.2 + + + + + genie::ZExpELFormFactorModel/Default + + + false + genie::TransverseEnhancementFFModel/Default + + + + + + + P33(1232),S11(1535),D13(1520),S11(1650),D13(1700),D15(1675), + S31(1620),D33(1700),P11(1440),P33(1600),P13(1720),F15(1680), + P31(1910),P33(1920),F35(1905),F37(1950),P11(1710) + + + + + + + + + + false + + + + + + + genie_pdg_table.txt + + + + + + 0.010 + 1000.000 + + + + diff --git a/config/GTEST23_10k/MECInteractionListGenerator.xml b/config/GTEST23_10k/MECInteractionListGenerator.xml new file mode 100644 index 000000000..a6ea945a6 --- /dev/null +++ b/config/GTEST23_10k/MECInteractionListGenerator.xml @@ -0,0 +1,47 @@ + + + + + + + + + + false + true + + + + + + + true + false + false + + + + false + true + false + + + + false + false + true + + + + diff --git a/config/GTEST23_10k/ModelConfiguration.xml b/config/GTEST23_10k/ModelConfiguration.xml new file mode 100644 index 000000000..5e91814f9 --- /dev/null +++ b/config/GTEST23_10k/ModelConfiguration.xml @@ -0,0 +1,146 @@ + + + + + + + + + + + + + + genie::LocalFGM/Default + + + + + + + + true + genie::HAIntranuke2018/Default + + + + + + + genie::NievesQELCCPXSec/ZExp + + + genie::AhrensNCELPXSec/Default + genie::RosenbluthPXSec/Default + + + genie::BergerSehgalRESPXSec2014/NoPauliBlock + genie::BergerSehgalRESPXSec2014/NoPauliBlock + genie::BergerSehgalRESPXSec2014/EM-NoPauliBlock + + + + + + + + + genie::KNOTunedQPMDISPXSec/Default + genie::KNOTunedQPMDISPXSec/Default + genie::KNOTunedQPMDISPXSec/Default + + + + + + + + + genie::BergerSehgalCOHPiPXSec2015/Default + genie::BergerSehgalCOHPiPXSec2015/Default + + + + + genie::BardinIMDRadCorPXSec/Default + genie::IMDAnnihilationPXSec/Default + genie::NuElectronPXSec/Default + genie::KovalenkoQELCharmPXSec/Default + genie::AivazisCharmPXSecLO/CC-Default + genie::AlamSimoAtharVacasSKPXSec2014/Default + genie::PaisQELLambdaPXSec/Default + genie::H3AMNuGammaPXSec/Default + genie::ReinDFRPXSec/Default + genie::ReinDFRPXSec/Default + + + genie::NievesSimoVacasMECPXSec2016/Default + + genie::EmpiricalMECPXSec2015/Default + genie::EmpiricalMECPXSec2015/Default + genie::DummyPXSec/Default + genie::NNBarOscDummyPXSec/Default + + + + + + diff --git a/config/GTEST23_10k/TuneGeneratorList.xml b/config/GTEST23_10k/TuneGeneratorList.xml new file mode 100644 index 000000000..938e22406 --- /dev/null +++ b/config/GTEST23_10k/TuneGeneratorList.xml @@ -0,0 +1,72 @@ + + + + + + + + 18 + genie::EventGenerator/QEL-CC + genie::EventGenerator/QEL-NC + genie::EventGenerator/RES-CC + genie::EventGenerator/RES-NC + genie::EventGenerator/DIS-CC + genie::EventGenerator/DIS-NC + genie::EventGenerator/COH-CC-PION + genie::EventGenerator/COH-NC-PION + genie::EventGenerator/DIS-CC-CHARM + genie::EventGenerator/QEL-CC-CHARM + genie::EventGenerator/NUE-EL + genie::EventGenerator/IMD + genie::EventGenerator/IMD-ANH + genie::EventGenerator/DFR-CC + genie::EventGenerator/DFR-NC + genie::EventGenerator/QEL-CC-LAMBDA + genie::EventGenerator/MEC-CC + genie::EventGenerator/MEC-NC + + + + diff --git a/config/Messenger.xml b/config/Messenger.xml index 20607f9dd..cf8bc114f 100644 --- a/config/Messenger.xml +++ b/config/Messenger.xml @@ -206,6 +206,7 @@ NOTICE NOTICE NOTICE + NOTICE WARN WARN WARN @@ -223,7 +224,6 @@ WARN ERROR FATAL - NOTICE NOTICE NOTICE diff --git a/config/ZExpELFormFactorModel.xml b/config/ZExpELFormFactorModel.xml new file mode 100644 index 000000000..644be8559 --- /dev/null +++ b/config/ZExpELFormFactorModel.xml @@ -0,0 +1,64 @@ + + + + + + + + + QuasiElastic + + true + 4 + + -0.21 + 0.0779198 + + -1.486 + -0.096 + 1.82 + 1.29 + -4.08895 + 0.175959 + 4.85982 + -0.97755 + 0.084 + -0.279 + -0.15 + 0.35 + 2.70696 + -0.420869 + -2.65913 + 0.0 + + + + + + + + + diff --git a/config/master_config.xml b/config/master_config.xml index 048508c3f..36930ec6a 100644 --- a/config/master_config.xml +++ b/config/master_config.xml @@ -134,6 +134,7 @@ BBA07ELFormFactorsModel.xml DipoleAxialFormFactorModel.xml ZExpAxialFormFactorModel.xml + ZExpELFormFactorModel.xml KuzminNaumov2016AxialFormFactorModel.xml LwlynSmithFFCC.xml LwlynSmithFFDeltaS.xml diff --git a/src/Physics/QuasiElastic/XSection/LinkDef.h b/src/Physics/QuasiElastic/XSection/LinkDef.h index 05ffb8045..00f591a2a 100644 --- a/src/Physics/QuasiElastic/XSection/LinkDef.h +++ b/src/Physics/QuasiElastic/XSection/LinkDef.h @@ -27,6 +27,7 @@ #pragma link C++ class genie::AxialFormFactor; #pragma link C++ class genie::DipoleAxialFormFactorModel; #pragma link C++ class genie::ZExpAxialFormFactorModel; +#pragma link C++ class genie::ZExpELFormFactorModel; #pragma link C++ class genie::KuzminNaumov2016AxialFormFactorModel; #pragma link C++ class genie::NievesQELCCPXSec; #pragma link C++ class genie::SuSAv2QELPXSec; diff --git a/src/Physics/QuasiElastic/XSection/ZExpELFormFactorModel.cxx b/src/Physics/QuasiElastic/XSection/ZExpELFormFactorModel.cxx new file mode 100644 index 000000000..ee8c45511 --- /dev/null +++ b/src/Physics/QuasiElastic/XSection/ZExpELFormFactorModel.cxx @@ -0,0 +1,736 @@ +//____________________________________________________________________________ +/* + Copyright (c) 2003-2019, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org + or see $GENIE/LICENSE + +\author Kaushik Borah + University of Kentucky, Lexington, KY 40506, USA + based off ZExpAxialFormFactorModel by + Aaron Meyer + University of Chicago, Chicago, Illinois, 60637, USA + + For the class documentation see the corresponding header file. + + + +*/ +//____________________________________________________________________________ + +#include +#include + +#include "Framework/Algorithm/AlgConfigPool.h" +#include "Framework/Conventions/Constants.h" +#include "Framework/Interaction/Interaction.h" +#include "Physics/QuasiElastic/XSection/ZExpELFormFactorModel.h" +#include "Framework/Messenger/Messenger.h" + +using std::ostringstream; + +using namespace genie; + +//____________________________________________________________________________ +ZExpELFormFactorModel::ZExpELFormFactorModel() : +ELFormFactorsModelI("genie::ZExpELFormFactorModel") +{ + +} +//____________________________________________________________________________ +ZExpELFormFactorModel::ZExpELFormFactorModel(string config) : +ELFormFactorsModelI("genie::ZExpELFormFactorModel", config) +{ + +} +//____________________________________________________________________________ +ZExpELFormFactorModel::~ZExpELFormFactorModel() +{ + delete[] fZ_APn; + delete[] fZ_BPn; + delete[] fZ_ANn; + delete[] fZ_BNn; +} +//____________________________________________________________________________ +double ZExpELFormFactorModel::Gep(const Interaction * interaction) const +{ + // calculate and return Gep + double q2 = interaction->KinePtr()->q2(); + double zparam = this->CalculateZ(q2); + if (zparam != zparam) // checks for nan + { + LOG("ZExpELFormFactorModel",pWARN) << "Undefined expansion parameter"; + return 0.; + } + double gep = 0.; + for (int ki=0;ki<=fKmax+(fQ4limit ? 4 : 0);ki++) + { + gep = gep + TMath::Power(zparam,ki) * fZ_APn[ki]; + } + + return gep; +} +//____________________________________________________________________________ +double ZExpELFormFactorModel::Gmp(const Interaction * interaction) const +{ + // calculate and return Gmp + double q2 = interaction->KinePtr()->q2(); + double zparam = this->CalculateZ(q2); + if (zparam != zparam) // checks for nan + { + LOG("ZExpELFormFactorModel",pWARN) << "Undefined expansion parameter"; + return 0.; + } + double gmp = 0.; + for (int ki=0;ki<=fKmax+(fQ4limit ? 4 : 0);ki++) + { + gmp = gmp + TMath::Power(zparam,ki) * fZ_BPn[ki]; + } + + return gmp; +} + +//____________________________________________________________________________ +double ZExpELFormFactorModel::Gen(const Interaction * interaction) const +{ + // calculate and return Gen + double q2 = interaction->KinePtr()->q2(); + double zparam = this->CalculateZ(q2); + if (zparam != zparam) // checks for nan + { + LOG("ZExpELFormFactorModel",pWARN) << "Undefined expansion parameter"; + return 0.; + } + double gen = 0.; + for (int ki=0;ki<=fKmax+(fQ4limit ? 4 : 0);ki++) + { + gen = gen + TMath::Power(zparam,ki) * fZ_ANn[ki]; + } + + return gen; +} +//____________________________________________________________________________ +double ZExpELFormFactorModel::Gmn(const Interaction * interaction) const +{ + // calculate and return Gmn + double q2 = interaction->KinePtr()->q2(); + double zparam = this->CalculateZ(q2); + if (zparam != zparam) // checks for nan + { + LOG("ZExpELFormFactorModel",pWARN) << "Undefined expansion parameter"; + return 0.; + } + double gmn = 0.; + for (int ki=0;ki<=fKmax+(fQ4limit ? 4 : 0);ki++) + { + gmn = gmn + TMath::Power(zparam,ki) * fZ_BNn[ki]; + } + + return gmn; +} +//____________________________________________________________________________ +double ZExpELFormFactorModel::CalculateZ(double q2) const +{ + + // calculate z expansion parameter + double znum = TMath::Sqrt(fTcut - q2) - TMath::Sqrt(fTcut - fT0); + double zden = TMath::Sqrt(fTcut - q2) + TMath::Sqrt(fTcut - fT0); + + return znum/zden; +} +//____________________________________________________________________________ +void ZExpELFormFactorModel::FixCoeffs(void) +{ + //if (fKmax < 1 ) { fKmax = 1; } + //else if (fKmax > 10) { fKmax = 10; } + + if (fQ4limit) this->FixQ4Limit(); + else this->FixEL0(); +} +//____________________________________________________________________________ +void ZExpELFormFactorModel::FixEL0(void) +{ + // Function to fix form factor such that Gep(q2=0) = 1 + // For T0 = 0, this will set Gep0 = 1 + double zparam = this->CalculateZ(0.); + double gep = 0.; + for (int ki=1;ki<=fKmax;ki++) + { + gep = gep + TMath::Power(zparam,ki) * fZ_APn[ki]; + } + fZ_APn[0] = fGep0 - gep; + + // Function to fix form factor such that Gmp(q2=0) = 2.7928 + // For T0 = 0, this will set Gmp0 = 2.7928 + + double gmp = 0.; + for (int ki=1;ki<=fKmax;ki++) + { + gmp = gmp + TMath::Power(zparam,ki) * fZ_BPn[ki]; + } + fZ_BPn[0] = fGmp0 - gmp; + + // Function to fix form factor such that Gen(q2=0) = 0 + // For T0 = 0, this will set Gen0 = 0 + + double gen = 0.; + for (int ki=1;ki<=fKmax;ki++) + { + gen = gen + TMath::Power(zparam,ki) * fZ_ANn[ki]; + } + fZ_ANn[0] = fGen0 - gen; + + // Function to fix form factor such that Gmn(q2=0) = -1.9130 + // For T0 = 0, this will set Gmn0 = -1.9130 + + double gmn = 0.; + for (int ki=1;ki<=fKmax;ki++) + { + gmn = gmn + TMath::Power(zparam,ki) * fZ_BNn[ki]; + } + fZ_BNn[0] = fGmn0 - gmn; + +} +//____________________________________________________________________________ +void ZExpELFormFactorModel::FixQ4Limit(void) +{ + // fixes parameters such that the model limits to 1/Q^4 at large t + // -- requires at least 5 parameters to do so -- + // 4 parameters for Q^4 behavior + + // will use AP_0 and AP_Kmax through AP_Kmax-3 to do the fitting + // calculate some useful numbers (redundancy for readability) + double kp4 = (double)fKmax+4; + double kp3 = (double)fKmax+3; + double kp2 = (double)fKmax+2; + double kp1 = (double)fKmax+1; + double kp0 = (double)fKmax ; + //double km5 = (double)fKmax-5; + double z0 = this->CalculateZ(0.); + double zkp4 = TMath::Power(z0,(int)kp4); + double zkp3 = TMath::Power(z0,(int)kp3); + double zkp2 = TMath::Power(z0,(int)kp2); + double zkp1 = TMath::Power(z0,(int)kp1); + + // denominator + double denom = \ + 6. - kp4*kp3*kp2*zkp1 + 3.*kp4*kp3*kp1*zkp2 \ + - 3.*kp4*kp2*kp1*zkp3 + kp3*kp2*kp1*zkp4; + + // extra parameters (effectively constants) + // number refers to the number of derivatives + double b0 = 0.; + for (int ki = 1;ki <= fKmax;ki++) + { + b0 = b0 + fZ_APn[ki]; + } + + double b0z = -fGep0; + for (int ki = 1;ki <= fKmax;ki++) + { + b0z = b0z + fZ_APn[ki]*TMath::Power(z0,ki); + } + + double b1 = 0.; + for (int ki = 1;ki <= fKmax;ki++) + { + b1 = b1 + ki*fZ_APn[ki]; + } + + double b2 = 0.; + for (int ki = 1;ki <= fKmax;ki++) + { + b2 = b2 + ki*(ki-1)*fZ_APn[ki]; + } + + double b3 = 0.; + for (int ki = 1;ki <= fKmax;ki++) + { + b3 = b3 + ki*(ki-1)*(ki-2)*fZ_APn[ki]; + } + + // Assign new parameters + fZ_APn[(int)kp4] = (1./denom) * \ + ( (b0-b0z)*kp3*kp2*kp1 \ + + b3*( -1. + .5*kp3*kp2*zkp1 - kp3*kp1*zkp2 \ + +.5*kp2*kp1*zkp3 ) \ + + b2*( 3.*kp1 - kp3*kp2*kp1*zkp1 \ + +kp3*kp1*(2*fKmax+1)*zkp2 - kp2*kp1*kp0*zkp3 ) \ + + b1*( -3.*kp2*kp1 + .5*kp3*kp2*kp2*kp1*zkp1 \ + -kp3*kp2*kp1*kp0*zkp2 + .5*kp2*kp1*kp1*kp0*zkp3) ); + + fZ_APn[(int)kp3] = (1./denom) * \ + ( -3.*(b0-b0z)*kp4*kp2*kp1 \ + + b3*( 3. - kp4*kp2*zkp1 + (3./2.)*kp4*kp1*zkp2 \ + -.5*kp2*kp1*zkp4 ) \ + + b2*( -3.*(3*fKmax+4) + kp4*kp2*(2*fKmax+3)*zkp1 \ + -3.*kp4*kp1*kp1*zkp2 + kp2*kp1*kp0*zkp4 ) \ + + b1*( 3.*kp1*(3*fKmax+8) - kp4*kp3*kp2*kp1*zkp1 \ + +(3./2.)*kp4*kp3*kp1*kp0*zkp2 - .5*kp2*kp1*kp1*kp0*zkp4) ); + + fZ_APn[(int)kp2] = (1./denom) * \ + ( 3.*(b0-b0z)*kp4*kp3*kp1 \ + + b3*( -3. + .5*kp4*kp3*zkp1 - (3./2.)*kp4*kp1*zkp3 \ + +kp3*kp1*zkp4 ) \ + + b2*( 3.*(3*fKmax+5) - kp4*kp3*kp2*zkp1 \ + +3.*kp4*kp1*kp1*zkp3 - kp3*kp1*(2*fKmax+1)*zkp4) \ + + b1*( -3.*kp3*(3*fKmax+4) + .5*kp4*kp3*kp3*kp2*zkp1 \ + -(3./2.)*kp4*kp3*kp1*kp0*zkp3 + kp3*kp2*kp1*kp0*zkp4) ); + + fZ_APn[(int)kp1] = (1./denom) * \ + ( -(b0-b0z)*kp4*kp3*kp2 \ + + b3*( 1. - .5*kp4*kp3*zkp2 + kp4*kp2*zkp3 \ + -.5*kp3*kp2*zkp4 ) \ + + b2*( -3.*kp2 + kp4*kp3*kp2*zkp2 \ + -kp4*kp2*(2*fKmax+3)*zkp3 + kp3*kp2*kp1*zkp4) \ + + b1*( 3.*kp3*kp2 - .5*kp4*kp3*kp3*kp2*zkp2 \ + +kp4*kp3*kp2*kp1*zkp3 - .5*kp3*kp2*kp2*kp1*zkp4) ); + + fZ_APn[0] = (1./denom) * \ + ( -6.*b0z \ + + b0*( kp4*kp3*kp2*zkp1 - 3.*kp4*kp3*kp1*zkp2 \ + +3.*kp4*kp2*kp1*zkp3 - kp3*kp2*kp1*zkp4 ) \ + + b3*( -zkp1 + 3.*zkp2 - 3.*zkp3 + zkp4 ) \ + + b2*( 3.*kp2*zkp1 - 3.*(3*fKmax+5)*zkp2 \ + +3.*(3*fKmax+4)*zkp3 - 3.*kp1*zkp4 ) \ + + b1*( -3.*kp3*kp2*zkp1 + 3.*kp3*(3*fKmax+4)*zkp2 \ + -3.*kp1*(3*fKmax+8)*zkp3 + 3.*kp2*kp1*zkp4 ) ); + + + + + + + // fixes parameters such that the model limits to 1/Q^4 at large t + // -- requires at least 5 parameters to do so -- + // 4 parameters for Q^4 behavior + + // will use BP_0 and BP_Kmax through BP_Kmax-3 to do the fitting + // calculate some useful numbers (redundancy for readability) + double kkp4 = (double)fKmax+4; + double kkp3 = (double)fKmax+3; + double kkp2 = (double)fKmax+2; + double kkp1 = (double)fKmax+1; + double kkp0 = (double)fKmax ; + //double km5 = (double)fKmax-5; + double kz0 = this->CalculateZ(0.); + double zkkp4 = TMath::Power(kz0,(int)kkp4); + double zkkp3 = TMath::Power(kz0,(int)kkp3); + double zkkp2 = TMath::Power(kz0,(int)kkp2); + double zkkp1 = TMath::Power(kz0,(int)kkp1); + + // denominator + double kdenom = \ + 6. - kkp4*kkp3*kkp2*zkkp1 + 3.*kkp4*kkp3*kkp1*zkkp2 \ + - 3.*kkp4*kkp2*kkp1*zkkp3 + kkp3*kkp2*kkp1*zkkp4; + + // extra parameters (effectively constants) + // number refers to the number of derivatives + double kb0 = 0.; + for (int ki = 1;ki <= fKmax;ki++) + { + kb0 = kb0 + fZ_BPn[ki]; + } + + double kb0z = -fGmp0; + for (int ki = 1;ki <= fKmax;ki++) + { + kb0z = kb0z + fZ_BPn[ki]*TMath::Power(kz0,ki); + } + + double kb1 = 0.; + for (int ki = 1;ki <= fKmax;ki++) + { + kb1 = kb1 + ki*fZ_BPn[ki]; + } + + double kb2 = 0.; + for (int ki = 1;ki <= fKmax;ki++) + { + kb2 = kb2 + ki*(ki-1)*fZ_BPn[ki]; + } + + double kb3 = 0.; + for (int ki = 1;ki <= fKmax;ki++) + { + kb3 = kb3 + ki*(ki-1)*(ki-2)*fZ_BPn[ki]; + } + + // Assign new parameters + fZ_BPn[(int)kkp4] = (1./kdenom) * \ + ( (kb0-kb0z)*kkp3*kkp2*kkp1 \ + + kb3*( -1. + .5*kkp3*kkp2*zkkp1 - kkp3*kkp1*zkkp2 \ + +.5*kkp2*kkp1*zkkp3 ) \ + + kb2*( 3.*kkp1 - kkp3*kkp2*kkp1*zkkp1 \ + +kkp3*kkp1*(2*fKmax+1)*zkkp2 - kkp2*kkp1*kkp0*zkkp3 ) \ + + kb1*( -3.*kkp2*kkp1 + .5*kkp3*kkp2*kkp2*kkp1*zkkp1 \ + -kkp3*kkp2*kkp1*kkp0*zkkp2 + .5*kkp2*kkp1*kkp1*kkp0*zkkp3) ); + + fZ_BPn[(int)kkp3] = (1./kdenom) * \ + ( -3.*(kb0-kb0z)*kkp4*kkp2*kkp1 \ + + kb3*( 3. - kkp4*kkp2*zkkp1 + (3./2.)*kkp4*kkp1*zkkp2 \ + -.5*kkp2*kkp1*zkkp4 ) \ + + kb2*( -3.*(3*fKmax+4) + kkp4*kkp2*(2*fKmax+3)*zkkp1 \ + -3.*kkp4*kkp1*kkp1*zkkp2 + kkp2*kkp1*kkp0*zkkp4 ) \ + + kb1*( 3.*kkp1*(3*fKmax+8) - kkp4*kkp3*kkp2*kkp1*zkkp1 \ + +(3./2.)*kkp4*kkp3*kkp1*kkp0*zkkp2 - .5*kkp2*kkp1*kkp1*kkp0*zkkp4) ); + + fZ_BPn[(int)kkp2] = (1./kdenom) * \ + ( 3.*(kb0-kb0z)*kkp4*kkp3*kkp1 \ + + kb3*( -3. + .5*kkp4*kkp3*zkkp1 - (3./2.)*kkp4*kkp1*zkkp3 \ + +kkp3*kkp1*zkkp4 ) \ + + kb2*( 3.*(3*fKmax+5) - kkp4*kkp3*kkp2*zkkp1 \ + +3.*kkp4*kkp1*kkp1*zkkp3 - kkp3*kkp1*(2*fKmax+1)*zkkp4) \ + + kb1*( -3.*kkp3*(3*fKmax+4) + .5*kkp4*kkp3*kkp3*kkp2*zkkp1 \ + -(3./2.)*kkp4*kkp3*kkp1*kkp0*zkkp3 + kkp3*kkp2*kkp1*kkp0*zkkp4) ); + + fZ_BPn[(int)kkp1] = (1./kdenom) * \ + ( -(kb0-kb0z)*kkp4*kkp3*kkp2 \ + + kb3*( 1. - .5*kkp4*kkp3*zkkp2 + kkp4*kkp2*zkkp3 \ + -.5*kkp3*kkp2*zkkp4 ) \ + + kb2*( -3.*kkp2 + kkp4*kkp3*kkp2*zkkp2 \ + -kkp4*kkp2*(2*fKmax+3)*zkkp3 + kkp3*kkp2*kkp1*zkkp4) \ + + kb1*( 3.*kkp3*kkp2 - .5*kkp4*kkp3*kkp3*kkp2*zkkp2 \ + +kkp4*kkp3*kkp2*kkp1*zkkp3 - .5*kkp3*kkp2*kkp2*kkp1*zkkp4) ); + + fZ_BPn[0] = (1./kdenom) * \ + ( -6.*kb0z \ + + kb0*( kkp4*kkp3*kkp2*zkkp1 - 3.*kkp4*kkp3*kkp1*zkkp2 \ + +3.*kkp4*kkp2*kkp1*zkkp3 - kkp3*kkp2*kkp1*zkkp4 ) \ + + kb3*( -zkkp1 + 3.*zkkp2 - 3.*zkkp3 + zkkp4 ) \ + + kb2*( 3.*kkp2*zkkp1 - 3.*(3*fKmax+5)*zkkp2 \ + +3.*(3*fKmax+4)*zkkp3 - 3.*kkp1*zkkp4 ) \ + + kb1*( -3.*kkp3*kkp2*zkkp1 + 3.*kkp3*(3*fKmax+4)*zkkp2 \ + -3.*kkp1*(3*fKmax+8)*zkkp3 + 3.*kkp2*kkp1*zkkp4 ) ); + + + + + + + // fixes parameters such that the model limits to 1/Q^4 at large t + // -- requires at least 5 parameters to do so -- + // 4 parameters for Q^4 behavior + + // will use AN_0 and AN_Kmax through AN_Kmax-3 to do the fitting + // calculate some useful numbers (redundancy for readability) + double lp4 = (double)fKmax+4; + double lp3 = (double)fKmax+3; + double lp2 = (double)fKmax+2; + double lp1 = (double)fKmax+1; + double lp0 = (double)fKmax ; + //double km5 = (double)fKmax-5; + double lz0 = this->CalculateZ(0.); + double zlp4 = TMath::Power(lz0,(int)lp4); + double zlp3 = TMath::Power(lz0,(int)lp3); + double zlp2 = TMath::Power(lz0,(int)lp2); + double zlp1 = TMath::Power(lz0,(int)lp1); + + // ldenominator + double ldenom = \ + 6. - lp4*lp3*lp2*zlp1 + 3.*lp4*lp3*lp1*zlp2 \ + - 3.*lp4*lp2*lp1*zlp3 + lp3*lp2*lp1*zlp4; + + // extra parameters (effectively constants) + // number refers to the number of derivatives + double lb0 = 0.; + for (int ki = 1;ki <= fKmax;ki++) + { + lb0 = lb0 + fZ_ANn[ki]; + } + + double lb0z = -fGen0; + for (int ki = 1;ki <= fKmax;ki++) + { + lb0z = lb0z + fZ_ANn[ki]*TMath::Power(lz0,ki); + } + + double lb1 = 0.; + for (int ki = 1;ki <= fKmax;ki++) + { + lb1 = lb1 + ki*fZ_ANn[ki]; + } + + double lb2 = 0.; + for (int ki = 1;ki <= fKmax;ki++) + { + lb2 = lb2 + ki*(ki-1)*fZ_ANn[ki]; + } + + double lb3 = 0.; + for (int ki = 1;ki <= fKmax;ki++) + { + lb3 = lb3 + ki*(ki-1)*(ki-2)*fZ_ANn[ki]; + } + + // Assign new parameters + fZ_ANn[(int)lp4] = (1./ldenom) * \ + ( (lb0-lb0z)*lp3*lp2*lp1 \ + + lb3*( -1. + .5*lp3*lp2*zlp1 - lp3*lp1*zlp2 \ + +.5*lp2*lp1*zlp3 ) \ + + lb2*( 3.*lp1 - lp3*lp2*lp1*zlp1 \ + +lp3*lp1*(2*fKmax+1)*zlp2 - lp2*lp1*lp0*zlp3 ) \ + + lb1*( -3.*lp2*lp1 + .5*lp3*lp2*lp2*lp1*zlp1 \ + -lp3*lp2*lp1*lp0*zlp2 + .5*lp2*lp1*lp1*lp0*zlp3) ); + + fZ_ANn[(int)lp3] = (1./ldenom) * \ + ( -3.*(lb0-lb0z)*lp4*lp2*lp1 \ + + lb3*( 3. - lp4*lp2*zlp1 + (3./2.)*lp4*lp1*zlp2 \ + -.5*lp2*lp1*zlp4 ) \ + + lb2*( -3.*(3*fKmax+4) + lp4*lp2*(2*fKmax+3)*zlp1 \ + -3.*lp4*lp1*lp1*zlp2 + lp2*lp1*lp0*zlp4 ) \ + + lb1*( 3.*lp1*(3*fKmax+8) - lp4*lp3*lp2*lp1*zlp1 \ + +(3./2.)*lp4*lp3*lp1*lp0*zlp2 - .5*lp2*lp1*lp1*lp0*zlp4) ); + + fZ_ANn[(int)lp2] = (1./ldenom) * \ + ( 3.*(lb0-lb0z)*lp4*lp3*lp1 \ + + lb3*( -3. + .5*lp4*lp3*zlp1 - (3./2.)*lp4*lp1*zlp3 \ + +lp3*lp1*zlp4 ) \ + + lb2*( 3.*(3*fKmax+5) - lp4*lp3*lp2*zlp1 \ + +3.*lp4*lp1*lp1*zlp3 - lp3*lp1*(2*fKmax+1)*zlp4) \ + + lb1*( -3.*lp3*(3*fKmax+4) + .5*lp4*lp3*lp3*lp2*zlp1 \ + -(3./2.)*lp4*lp3*lp1*lp0*zlp3 + lp3*lp2*lp1*lp0*zlp4) ); + + fZ_ANn[(int)lp1] = (1./ldenom) * \ + ( -(lb0-lb0z)*lp4*lp3*lp2 \ + + lb3*( 1. - .5*lp4*lp3*zlp2 + lp4*lp2*zlp3 \ + -.5*lp3*lp2*zlp4 ) \ + + lb2*( -3.*lp2 + lp4*lp3*lp2*zlp2 \ + -lp4*lp2*(2*fKmax+3)*zlp3 + lp3*lp2*lp1*zlp4) \ + + lb1*( 3.*lp3*lp2 - .5*lp4*lp3*lp3*lp2*zlp2 \ + +lp4*lp3*lp2*lp1*zlp3 - .5*lp3*lp2*lp2*lp1*zlp4) ); + + fZ_ANn[0] = (1./ldenom) * \ + ( -6.*lb0z \ + + lb0*( lp4*lp3*lp2*zlp1 - 3.*lp4*lp3*lp1*zlp2 \ + +3.*lp4*lp2*lp1*zlp3 - lp3*lp2*lp1*zlp4 ) \ + + lb3*( -zlp1 + 3.*zlp2 - 3.*zlp3 + zlp4 ) \ + + lb2*( 3.*lp2*zlp1 - 3.*(3*fKmax+5)*zlp2 \ + +3.*(3*fKmax+4)*zlp3 - 3.*lp1*zlp4 ) \ + + lb1*( -3.*lp3*lp2*zlp1 + 3.*lp3*(3*fKmax+4)*zlp2 \ + -3.*lp1*(3*fKmax+8)*zlp3 + 3.*lp2*lp1*zlp4 ) ); + + + + + + + + // fixes parameters such that the model limits to 1/Q^4 at large t + // -- requires at least 5 parameters to do so -- + // 4 parameters for Q^4 behavior + + // will use BN_0 and BN_Kmax through BN_Kmax-3 to do the fitting + // calculate some useful numbers (redundancy for readability) + double llp4 = (double)fKmax+4; + double llp3 = (double)fKmax+3; + double llp2 = (double)fKmax+2; + double llp1 = (double)fKmax+1; + double llp0 = (double)fKmax ; + //double km5 = (double)fKmax-5; + double llz0 = this->CalculateZ(0.); + double zllp4 = TMath::Power(llz0,(int)llp4); + double zllp3 = TMath::Power(llz0,(int)llp3); + double zllp2 = TMath::Power(llz0,(int)llp2); + double zllp1 = TMath::Power(llz0,(int)llp1); + + // denominator + double lldenom = \ + 6. - llp4*llp3*llp2*zllp1 + 3.*llp4*llp3*llp1*zllp2 \ + - 3.*llp4*llp2*llp1*zllp3 + llp3*llp2*llp1*zllp4; + + // extra parameters (effectively constants) + // number refers to the number of derivatives + double llb0 = 0.; + for (int ki = 1;ki <= fKmax;ki++) + { + llb0 = llb0 + fZ_BNn[ki]; + } + + double llb0z = -fGmn0; + for (int ki = 1;ki <= fKmax;ki++) + { + llb0z = llb0z + fZ_BNn[ki]*TMath::Power(llz0,ki); + } + + double llb1 = 0.; + for (int ki = 1;ki <= fKmax;ki++) + { + llb1 = llb1 + ki*fZ_BNn[ki]; + } + + double llb2 = 0.; + for (int ki = 1;ki <= fKmax;ki++) + { + llb2 = llb2 + ki*(ki-1)*fZ_BNn[ki]; + } + + double llb3 = 0.; + for (int ki = 1;ki <= fKmax;ki++) + { + llb3 = llb3 + ki*(ki-1)*(ki-2)*fZ_BNn[ki]; + } + + // Assign new parameters + fZ_BNn[(int)llp4] = (1./lldenom) * \ + ( (llb0-llb0z)*llp3*llp2*llp1 \ + + llb3*( -1. + .5*llp3*llp2*zllp1 - llp3*llp1*zllp2 \ + +.5*llp2*llp1*zllp3 ) \ + + llb2*( 3.*llp1 - llp3*llp2*llp1*zllp1 \ + +llp3*llp1*(2*fKmax+1)*zllp2 - llp2*llp1*llp0*zllp3 ) \ + + llb1*( -3.*llp2*llp1 + .5*llp3*llp2*llp2*llp1*zllp1 \ + -llp3*llp2*llp1*llp0*zllp2 + .5*llp2*llp1*llp1*llp0*zllp3) ); + + fZ_BNn[(int)llp3] = (1./lldenom) * \ + ( -3.*(llb0-llb0z)*llp4*llp2*llp1 \ + + llb3*( 3. - llp4*llp2*zllp1 + (3./2.)*llp4*llp1*zllp2 \ + -.5*llp2*llp1*zllp4 ) \ + + llb2*( -3.*(3*fKmax+4) + llp4*llp2*(2*fKmax+3)*zllp1 \ + -3.*llp4*llp1*llp1*zllp2 + llp2*llp1*llp0*zllp4 ) \ + + llb1*( 3.*llp1*(3*fKmax+8) - llp4*llp3*llp2*llp1*zllp1 \ + +(3./2.)*llp4*llp3*llp1*llp0*zllp2 - .5*llp2*llp1*llp1*llp0*zllp4) ); + + fZ_BNn[(int)llp2] = (1./lldenom) * \ + ( 3.*(llb0-llb0z)*llp4*llp3*llp1 \ + + llb3*( -3. + .5*llp4*llp3*zllp1 - (3./2.)*llp4*llp1*zllp3 \ + +llp3*llp1*zllp4 ) \ + + llb2*( 3.*(3*fKmax+5) - llp4*llp3*llp2*zllp1 \ + +3.*llp4*llp1*llp1*zllp3 - llp3*llp1*(2*fKmax+1)*zllp4) \ + + llb1*( -3.*llp3*(3*fKmax+4) + .5*llp4*llp3*llp3*llp2*zllp1 \ + -(3./2.)*llp4*llp3*llp1*llp0*zllp3 + llp3*llp2*llp1*llp0*zllp4) ); + + fZ_BNn[(int)llp1] = (1./lldenom) * \ + ( -(llb0-llb0z)*llp4*llp3*llp2 \ + + llb3*( 1. - .5*llp4*llp3*zllp2 + llp4*llp2*zllp3 \ + -.5*llp3*llp2*zllp4 ) \ + + llb2*( -3.*llp2 + llp4*llp3*llp2*zllp2 \ + -llp4*llp2*(2*fKmax+3)*zllp3 + llp3*llp2*llp1*zllp4) \ + + llb1*( 3.*llp3*llp2 - .5*llp4*llp3*llp3*llp2*zllp2 \ + +llp4*llp3*llp2*llp1*zllp3 - .5*llp3*llp2*llp2*llp1*zllp4) ); + + fZ_BNn[0] = (1./lldenom) * \ + ( -6.*llb0z \ + + llb0*( llp4*llp3*llp2*zllp1 - 3.*llp4*llp3*llp1*zllp2 \ + +3.*llp4*llp2*llp1*zllp3 - llp3*llp2*llp1*zllp4 ) \ + + llb3*( -zllp1 + 3.*zllp2 - 3.*zllp3 + zllp4 ) \ + + llb2*( 3.*llp2*zllp1 - 3.*(3*fKmax+5)*zllp2 \ + +3.*(3*fKmax+4)*zllp3 - 3.*llp1*zllp4 ) \ + + llb1*( -3.*llp3*llp2*zllp1 + 3.*llp3*(3*fKmax+4)*zllp2 \ + -3.*llp1*(3*fKmax+8)*zllp3 + 3.*llp2*llp1*zllp4 ) ); + + + + +} +//____________________________________________________________________________ +void ZExpELFormFactorModel::Configure(const Registry & config) +{ + Algorithm::Configure(config); + this->LoadConfig(); +} +//____________________________________________________________________________ +void ZExpELFormFactorModel::Configure(string param_set) +{ + Algorithm::Configure(param_set); + this->LoadConfig(); +} +//____________________________________________________________________________ +void ZExpELFormFactorModel::LoadConfig(void ) +{ +// get config options from the configuration registry or set defaults +// from the global parameter list + + GetParam( "QEL-Q4limit", fQ4limit ) ; + GetParam( "QEL-Kmax", fKmax ) ; + + GetParam( "QEL-T0", fT0 ) ; + GetParam( "QEL-Tcut", fTcut ) ; + + GetParam( "QEL-Gep0", fGep0 ) ; + GetParam( "QEL-Gmp0", fGmp0 ) ; + GetParam( "QEL-Gen0", fGen0 ) ; + GetParam( "QEL-Gmn0", fGmn0 ) ; + assert(fKmax > 0); + + // z expansion coefficients + if (fQ4limit) fZ_APn = new double [fKmax+5]; + else fZ_APn = new double [fKmax+1]; + + // load the user-defined coefficient values + // -- AP0 and APn for nFixCoeffs(); +Interaction * interaction = new Interaction(); +for (int i=0;i<10;i++) { + double Q2 = i*0.1; + interaction->KinePtr()->SetQ2( Q2); + LOG("some_identifying_string",pWARN) + << "Q2=" <