diff --git a/config/EventGenerator.xml b/config/EventGenerator.xml index d1a418a16..fed181a64 100644 --- a/config/EventGenerator.xml +++ b/config/EventGenerator.xml @@ -830,5 +830,15 @@ XSecModel alg Yes Cross section model used at the thread + + + 3 + genie::InitialStateAppender/Default + genie::VertexGenerator/Default + genie::NormGenerator/Default + genie::NormInteractionListGenerator/Default + genie::NormXSec/Default + + diff --git a/config/EventGeneratorListAssembler.xml b/config/EventGeneratorListAssembler.xml index 30a19530f..b8704ac71 100644 --- a/config/EventGeneratorListAssembler.xml +++ b/config/EventGeneratorListAssembler.xml @@ -562,5 +562,16 @@ Generator-%d alg No genie::EventGenerator/PhotonCOH + + 1 + genie::EventGenerator/NORM + + + + 2 + genie::EventGenerator/QEL-CC + genie::EventGenerator/NORM + + diff --git a/config/NormGenerator.xml b/config/NormGenerator.xml new file mode 100644 index 000000000..04bf65e21 --- /dev/null +++ b/config/NormGenerator.xml @@ -0,0 +1,20 @@ + + + + + + + + + + + + diff --git a/config/NormInteractionListGenerator.xml b/config/NormInteractionListGenerator.xml new file mode 100644 index 000000000..cc322ae6d --- /dev/null +++ b/config/NormInteractionListGenerator.xml @@ -0,0 +1,21 @@ + + + + + + + + + + + + + + diff --git a/config/NormXSec.xml b/config/NormXSec.xml new file mode 100644 index 000000000..f18e0bcbc --- /dev/null +++ b/config/NormXSec.xml @@ -0,0 +1,20 @@ + + + + + + + + + + + diff --git a/config/QELEventGeneratorSM.xml b/config/QELEventGeneratorSM.xml index f60bfb170..0eadb84e1 100644 --- a/config/QELEventGeneratorSM.xml +++ b/config/QELEventGeneratorSM.xml @@ -4,28 +4,28 @@ Configuration for the QELEventGeneratorSM EventRecordVisitorI Configurable Parameters: -............................................................................................................. -Name Type Optional Comment Default -............................................................................................................. -MaxXSec-SafetyFactor double Yes Safety factor for the maximum 1.2 - differential cross section -MaxDiffv-SafetyFactor double Yes Safety factor for the maximum 1.2 - vmax(Q2)-vmin(Q2) -MaxXSec-DiffTolerance double Yes Max allowed 200*(xsec-xsecmax)/(xsec+xsecmax) 999999.0 - if xsec>xsecmax -UniformOverPhaseSpace bool Yes Kinematics uniformly over allowed phase space false - wgt = (phase_space_volume)*(diff_xsec)/(xsec) -IsNucleonInNucleus bool Yes Generate struck nucleon in nucleus true -Threshold-Q2 double Yes Q2-threshold for seeking the second maximum 2.00 -Cache-MinEnergy double Yes min E for which maxxsec is cached - 1.00 - forcing explicit calc +................................................................................................................................ +Name Type Optional Comment Default +................................................................................................................................ +MaxXSec-SafetyFactor double Yes Safety factor: ComputeMaxXSec -> ComputeMaxXSec * MaxXSec-SafetyFactor 1.0 +Safety-Factors vec-double Yes Safety factors: MaxXSec[nkey] -> MaxXSec[nkey] * Safety-Factors[nkey] all 1.0 +Interpolator-Types vec-string Yes Type of interpolator for each key in a branch all "" +MaxXSec-DiffTolerance double Yes Max allowed 200*(xsec-xsecmax)/(xsec+xsecmax) 999999.0 + if xsec>xsecmax +UniformOverPhaseSpace bool Yes Kinematics uniformly over allowed phase space false + wgt = (phase_space_volume)*(diff_xsec)/(xsec) +IsNucleonInNucleus bool Yes Generate struck nucleon in nucleus true +Threshold-Q2 double Yes Q2-threshold for seeking the second maximum 2.00 +Cache-MinEnergy double Yes min E for which maxxsec is cached - 1.00 + forcing explicit calc --> - + 1.06; 1.06; 1.06 + AKIMA; AKIMA; LINEAR diff --git a/config/master_config.xml b/config/master_config.xml index 4ff4ad940..048508c3f 100644 --- a/config/master_config.xml +++ b/config/master_config.xml @@ -264,4 +264,8 @@ EventLibraryInterface.xml EvtLibInteractionListGenerator.xml EvtLibPXSec.xml + + NormGenerator.xml + NormXSec.xml + NormInteractionListGenerator.xml diff --git a/src/Apps/gNtpConv.cxx b/src/Apps/gNtpConv.cxx index b9443ca00..54ed2a16a 100644 --- a/src/Apps/gNtpConv.cxx +++ b/src/Apps/gNtpConv.cxx @@ -324,6 +324,7 @@ void ConvertToGST(void) bool brIsMec = false; // Is MEC? bool brIsDfr = false; // Is Diffractive? bool brIsImd = false; // Is IMD? + bool brIsNrm = false; // Is Norm? bool brIsSingleK = false; // Is single kaon? bool brIsImdAnh = false; // Is IMD annihilation? bool brIsNuEL = false; // Is ve elastic? @@ -438,6 +439,7 @@ void ConvertToGST(void) s_tree->Branch("coh", &brIsCoh, "coh/O" ); s_tree->Branch("dfr", &brIsDfr, "dfr/O" ); s_tree->Branch("imd", &brIsImd, "imd/O" ); + s_tree->Branch("norm", &brIsNrm, "norm/O" ); s_tree->Branch("imdanh", &brIsImdAnh, "imdanh/O" ); s_tree->Branch("singlek", &brIsSingleK, "singlek/O" ); s_tree->Branch("nuel", &brIsNuEL, "nuel/O" ); @@ -657,9 +659,10 @@ void ConvertToGST(void) bool is_mec = proc_info.IsMEC(); bool is_amnugamma = proc_info.IsAMNuGamma(); bool is_hnl = proc_info.IsHNLDecay(); - + bool is_norm = proc_info.IsNorm(); + if (!hitnucl && neutrino) { - assert(is_coh || is_imd || is_imdanh || is_nuel | is_amnugamma || is_coh_el || is_hnl); + assert(is_coh || is_imd || is_imdanh || is_nuel | is_amnugamma || is_coh_el || is_hnl || is_norm); } // Hit quark - set only for DIS events @@ -954,6 +957,7 @@ void ConvertToGST(void) brIsCoh = is_coh; brIsDfr = is_dfr; brIsImd = is_imd; + brIsNrm = is_norm; brIsSingleK = is_singlek; brIsNuEL = is_nuel; brIsEM = is_em; diff --git a/src/Framework/Conventions/KinePhaseSpace.h b/src/Framework/Conventions/KinePhaseSpace.h index 59ad95638..c7806955e 100644 --- a/src/Framework/Conventions/KinePhaseSpace.h +++ b/src/Framework/Conventions/KinePhaseSpace.h @@ -24,7 +24,7 @@ using std::string; namespace genie { - +// Note: please attach new phase space enum element to the end of the list . typedef enum EKinePhaseSpace { kPSNull = 0, kPSfE, @@ -72,7 +72,8 @@ typedef enum EKinePhaseSpace { kPSlog10xlog10Q2fE, kPSEDNufE, // Used for Dark Neutrinos, two body final state kPSn1n2fE, - kPSn1n2n3fE + kPSn1n2n3fE, + kPSQ2vpfE } KinePhaseSpace_t; class KinePhaseSpace @@ -124,6 +125,7 @@ class KinePhaseSpace case(kPSElOlTpifE) : return "<{El,Omega_l,Theta_pi}|E>"; break; case(kPSTkTlctl) : return "<{Tk,Tl,cos(theta_l)}|E>"; break; case(kPSQ2vfE) : return "<{Q2,v}|E>"; break; + case(kPSQ2vpfE) : return "<{Q2,v,p}|E>"; break; // TODO: update this string when the appropriate kinematic variables are known case(kPSQELEvGen) : return ""; break; case(kPSDMELEvGen) : return ""; break; diff --git a/src/Framework/EventGen/GMCJDriver.cxx b/src/Framework/EventGen/GMCJDriver.cxx index f78fb91e8..2793f0625 100644 --- a/src/Framework/EventGen/GMCJDriver.cxx +++ b/src/Framework/EventGen/GMCJDriver.cxx @@ -652,9 +652,12 @@ void GMCJDriver::BootstrapXSecSplineSummation(void) // bit above the maximum beam energy (but below the maximum valid energy) // to avoid the evaluation of the cubic spline around the viscinity of // knots with zero y values (although the GENIE Spline object handles it) - double dE = fEmax/10.; double min = rE.min; - double max = (fEmax+dE < rE.max) ? fEmax+dE : rE.max; + double max = rE.max; + + // Because of edge issue (see GENIE docdb 297) these lines are commented out + // double dE = fEmax/10.; + // double max = (fEmax+dE < rE.max) ? fEmax+dE : rE.max; // in the logaritmic binning is important to have a narrow binning to // describe better the glashow resonance peak. diff --git a/src/Framework/GHEP/GHepRecord.cxx b/src/Framework/GHEP/GHepRecord.cxx index 19f7ba2a0..3306cc170 100644 --- a/src/Framework/GHEP/GHepRecord.cxx +++ b/src/Framework/GHEP/GHepRecord.cxx @@ -1183,8 +1183,8 @@ void GHepRecord::Print(ostream & stream) const case ( kPSQ2fE ) : stream << " dsig(Q2;E)/dQ2 = " << setfill(' ') << setw(13) << fDiffXSec/units::cm2 << " cm^2/GeV^2 |"; break; - case ( kPSQ2vfE ) : - stream << " dsig(Q2,v;E)/dQ2dv = " << setfill(' ') << setw(13) << fDiffXSec/units::cm2 << " cm^2/GeV^3 |"; + case ( kPSQ2vpfE ) : + stream << " dsig(Q2,v,p;E)/dQ2dvdp =" << setfill(' ') << setw(13) << fDiffXSec/units::cm2 << " cm^2/GeV^4 |"; break; case ( kPSWQ2fE ) : stream << " d2sig(W,Q2;E)/dWdQ2 = " << setfill(' ') << setw(13) << fDiffXSec/units::cm2 << " cm^2/GeV^3 |"; diff --git a/src/Framework/Interaction/KPhaseSpace.cxx b/src/Framework/Interaction/KPhaseSpace.cxx index 2b8a75185..f5945f193 100644 --- a/src/Framework/Interaction/KPhaseSpace.cxx +++ b/src/Framework/Interaction/KPhaseSpace.cxx @@ -88,6 +88,8 @@ double KPhaseSpace::Threshold(void) const if( ! pi.IsKnown() ) return 0; + if (pi.IsNorm() ) return 0; + if (pi.IsSingleKaon()) { int kaon_pdgc = xcls.StrangeHadronPdg(); double Mi = tgt.HitNucP4Ptr()->M(); // initial nucleon mass diff --git a/src/Framework/Interaction/ProcessInfo.cxx b/src/Framework/Interaction/ProcessInfo.cxx index 633faf295..297607d5a 100644 --- a/src/Framework/Interaction/ProcessInfo.cxx +++ b/src/Framework/Interaction/ProcessInfo.cxx @@ -128,6 +128,11 @@ bool ProcessInfo::IsIMDAnnihilation(void) const return (fScatteringType == kScIMDAnnihilation); } //____________________________________________________________________________ +bool ProcessInfo::IsNorm(void) const +{ + return (fScatteringType == kScNorm); +} +//____________________________________________________________________________ bool ProcessInfo::IsDarkMatterElectronElastic(void) const { return (fScatteringType == kScDarkMatterElectron); diff --git a/src/Framework/Interaction/ProcessInfo.h b/src/Framework/Interaction/ProcessInfo.h index ccbf5e6a7..48dcbe62d 100644 --- a/src/Framework/Interaction/ProcessInfo.h +++ b/src/Framework/Interaction/ProcessInfo.h @@ -71,6 +71,7 @@ class ProcessInfo : public TObject { bool IsNuElectronElastic (void) const; bool IsInverseMuDecay (void) const; bool IsIMDAnnihilation (void) const; + bool IsNorm (void) const; bool IsDarkMatterElectronElastic (void) const; bool IsInverseBetaDecay (void) const; bool IsGlashowResonance (void) const; diff --git a/src/Framework/Interaction/ScatteringType.h b/src/Framework/Interaction/ScatteringType.h index a51781d2d..ed421f467 100644 --- a/src/Framework/Interaction/ScatteringType.h +++ b/src/Framework/Interaction/ScatteringType.h @@ -54,7 +54,8 @@ typedef enum EScatteringType { kScPhotonResonance, kScDarkMatterElastic = 101, kScDarkMatterDeepInelastic, - kScDarkMatterElectron + kScDarkMatterElectron, + kScNorm } ScatteringType_t; class ScatteringType @@ -80,6 +81,7 @@ class ScatteringType case(kScInverseBetaDecay) : return "IBD"; break; case(kScGlashowResonance) : return "GLR"; break; case(kScIMDAnnihilation) : return "IMDAnh"; break; + case(kScNorm) : return "Norm"; break; case(kScPhotonCoherent) : return "PhotonCOH"; break; case(kScPhotonResonance) : return "PhotonRES"; break; case(kScDarkMatterElastic) : return "DMEL"; break; diff --git a/src/Framework/Numerical/Spline.cxx b/src/Framework/Numerical/Spline.cxx index 778a93c06..09db045a4 100644 --- a/src/Framework/Numerical/Spline.cxx +++ b/src/Framework/Numerical/Spline.cxx @@ -132,7 +132,7 @@ Spline::Spline(const Spline & spline) : TObject(), fInterpolator(0) { LOG("Spline", pDEBUG) << "Spline copy constructor"; - + this->InitSpline(); this->LoadFromTSpline3( *spline.GetAsTSpline(), spline.NKnots() ); } //___________________________________________________________________________ @@ -141,13 +141,15 @@ Spline::Spline(const TSpline3 & spline, int nknots) : { LOG("Spline", pDEBUG) << "Constructing spline from the input TSpline3 object"; - + this->InitSpline(); this->LoadFromTSpline3( spline, nknots ); } //___________________________________________________________________________ Spline::~Spline() { if(fInterpolator) delete fInterpolator; + if(fInterpolator5) delete fInterpolator5; + if(fGSLInterpolator) delete fGSLInterpolator; } //___________________________________________________________________________ bool Spline::LoadFromXmlFile(string filename, string xtag, string ytag) @@ -374,7 +376,12 @@ double Spline::Evaluate(double x) const if(!is0p && !is0n) { // both knots (on the left and right are non-zero) - just interpolate LOG("Spline", pDEBUG) << "Point is between non-zero knots"; - y = fInterpolator->Eval(x); + if (fInterpolatorType == "TSpline3") + y = fInterpolator->Eval(x); + else if (fInterpolatorType == "TSpline5") + y = fInterpolator5->Eval(x); + else + y = fGSLInterpolator->Eval(x); } else { // at least one of the neighboring knots has y=0 if(is0p && is0n) { @@ -721,8 +728,11 @@ void Spline::InitSpline(void) fXMin = 0.0; fXMax = 0.0; fYMax = 0.0; - + fInterpolator = 0; + fInterpolator5 = 0; + fGSLInterpolator = 0; + fInterpolatorType = "TSpline3"; fYCanBeNegative = false; @@ -732,6 +742,8 @@ void Spline::InitSpline(void) void Spline::ResetSpline(void) { if(fInterpolator) delete fInterpolator; + if(fInterpolator5) delete fInterpolator5; + if(fGSLInterpolator) delete fGSLInterpolator; this->InitSpline(); } //___________________________________________________________________________ @@ -750,7 +762,76 @@ void Spline::BuildSpline(int nentries, double x[], double y[]) if(fInterpolator) delete fInterpolator; fInterpolator = new TSpline3("spl3", x, y, nentries, "0"); - + LOG("Spline", pDEBUG) << "...done building spline"; } //___________________________________________________________________________ +void Spline::SetType(string type) +{ + if(!fInterpolator) return; + + fInterpolatorType = genie::utils::str::ToUpper(type); + + ROOT::Math::Interpolation::Type gsltype; + + if ( fInterpolatorType == "TSPLINE3" || fInterpolatorType == "" ) + { + fInterpolatorType == "TSpline3"; + return; + } + else if ( fInterpolatorType == "TSPLINE5" ) + { + fInterpolatorType = "TSpline5"; + if(fInterpolator5) delete fInterpolator5; + double x[fNKnots], y[fNKnots]; + for (int i=0; iGetKnot(i, x[i], y[i]); + } + fInterpolator5 = new TSpline5("spl5", x, y, fNKnots, "0"); + return; + } + else if ( fInterpolatorType == "LINEAR" ) + { + gsltype = ROOT::Math::Interpolation::kLINEAR; + } + else if ( fInterpolatorType == "POLYNOMIAL" ) + { + gsltype = ROOT::Math::Interpolation::kPOLYNOMIAL; + } + else if ( fInterpolatorType == "CSPLINE" ) + { + gsltype = ROOT::Math::Interpolation::kCSPLINE; + } + else if ( fInterpolatorType == "CSPLINE_PERIODIC" ) + { + gsltype = ROOT::Math::Interpolation::kCSPLINE_PERIODIC; + } + else if ( fInterpolatorType == "AKIMA" ) + { + gsltype = ROOT::Math::Interpolation::kAKIMA; + } + else if ( fInterpolatorType == "AKIMA_PERIODIC" ) + { + gsltype = ROOT::Math::Interpolation::kAKIMA_PERIODIC; + } + else + { + fInterpolatorType = "TSpline3"; + LOG("Spline", pWARN) + << "Unknown interpolator type. Setting it to default [TSpline3]."; + return; + } + + if(fGSLInterpolator) delete fGSLInterpolator; + vector x(fNKnots); + vector y(fNKnots); + for (int i=0; iGetKnot(i, x[i], y[i]); + } + fGSLInterpolator = new ROOT::Math::Interpolator(x, y, gsltype); + + +} +//___________________________________________________________________________ diff --git a/src/Framework/Numerical/Spline.h b/src/Framework/Numerical/Spline.h index 7533d3ea0..daf3a3ab3 100644 --- a/src/Framework/Numerical/Spline.h +++ b/src/Framework/Numerical/Spline.h @@ -7,10 +7,20 @@ Uses ROOT's TSpline3 for the actual interpolation and can retrieve function (x,y(x)) pairs from an XML file, a flat ascii file, a - TNtuple, a TTree or an SQL database. + TNtuple, a TTree or an SQL database.\n + + Update May 15, 2022 IK: + Adding as extra interpolators TSpline5 and + ROOT::Math::GSLInterpolator (LINEAR, POLYNOMIAL, CSPLINE, CSPLINE_PERIODIC, + AKIMA, AKIMA_PERIODIC) + + +\ref [1] GENIE docdb 297 \author Costas Andreopoulos - University of Liverpool & STFC Rutherford Appleton Laboratory + University of Liverpool & STFC Rutherford Appleton Laboratory \n + Igor Kakorin + Joint Institute for Nuclear Research \created May 04, 2004 @@ -22,12 +32,14 @@ #ifndef _SPLINE_H_ #define _SPLINE_H_ +#include #include #include #include #include #include +#include class TNtupleD; class TTree; @@ -106,6 +118,8 @@ class Spline : public TObject { void Add (double a); void Multiply (double a); void Divide (double a); + void SetType (string type); + string GetType (void) { return fInterpolatorType; } // Print knots void Print(ostream & stream) const; @@ -128,8 +142,11 @@ class Spline : public TObject { double fYMax; TSpline3 * fInterpolator; bool fYCanBeNegative; + TSpline5 * fInterpolator5; + ROOT::Math::Interpolator * fGSLInterpolator; + string fInterpolatorType; -ClassDef(Spline,1) +ClassDef(Spline,2) }; } diff --git a/src/Framework/Utils/CacheBranchFx.cxx b/src/Framework/Utils/CacheBranchFx.cxx index 612825262..adc67a1c3 100644 --- a/src/Framework/Utils/CacheBranchFx.cxx +++ b/src/Framework/Utils/CacheBranchFx.cxx @@ -65,7 +65,7 @@ void CacheBranchFx::AddValues(double x, double y) fFx.insert(map::value_type(x,y)); } //____________________________________________________________________________ -void CacheBranchFx::CreateSpline(void) +void CacheBranchFx::CreateSpline(string type) { int n = fFx.size(); double * x = new double[n]; @@ -81,6 +81,7 @@ void CacheBranchFx::CreateSpline(void) if(fSpline) delete fSpline; fSpline = new Spline(n,x,y); + fSpline->SetType(type); delete [] x; delete [] y; diff --git a/src/Framework/Utils/CacheBranchFx.h b/src/Framework/Utils/CacheBranchFx.h index 02e6f1807..dee3f2942 100644 --- a/src/Framework/Utils/CacheBranchFx.h +++ b/src/Framework/Utils/CacheBranchFx.h @@ -7,6 +7,18 @@ \author Costas Andreopoulos University of Liverpool & STFC Rutherford Appleton Laboratory + + Update May 15, 2022 IK: + Now type of spline can be: TSpline3, TSpline5 and + ROOT::Math::GSLInterpolator (LINEAR, POLYNOMIAL, CSPLINE, CSPLINE_PERIODIC, + AKIMA, AKIMA_PERIODIC) + +\ref [1] GENIE docdb 297 + +\author Costas Andreopoulos + University of Liverpool & STFC Rutherford Appleton Laboratory \n + Igor Kakorin + Joint Institute for Nuclear Research \created November 26, 2004 @@ -46,7 +58,7 @@ class CacheBranchFx : public CacheBranchI const map & Map (void) const { return fFx; } Spline * Spl (void) const { return fSpline; } - void CreateSpline(void); + void CreateSpline(string type = "TSpline3"); void AddValues(double x, double y); void Reset (void); diff --git a/src/Physics/Common/InitialStateAppender.cxx b/src/Physics/Common/InitialStateAppender.cxx index 631569768..57863c5b5 100644 --- a/src/Physics/Common/InitialStateAppender.cxx +++ b/src/Physics/Common/InitialStateAppender.cxx @@ -86,7 +86,7 @@ void InitialStateAppender::AddNucleus(GHepRecord * evrec) const const InitialState & init_state = interaction->InitState(); const ProcessInfo & proc_info = interaction->ProcInfo(); - bool is_nucleus = init_state.Tgt().IsNucleus(); + bool is_nucleus = init_state.Tgt().IsNucleus() || proc_info.IsNorm(); if(!is_nucleus && !proc_info.IsGlashowResonance() && !proc_info.IsPhotonCoherent()) { LOG("ISApp", pINFO) << "Not an interaction with a nuclear target - no nucleus to add"; diff --git a/src/Physics/Common/KineGeneratorWithCache.cxx b/src/Physics/Common/KineGeneratorWithCache.cxx index 1013b0169..c1e44e31a 100644 --- a/src/Physics/Common/KineGeneratorWithCache.cxx +++ b/src/Physics/Common/KineGeneratorWithCache.cxx @@ -31,20 +31,20 @@ using std::map; using namespace genie; //___________________________________________________________________________ -KineGeneratorWithCache::KineGeneratorWithCache() : -EventRecordVisitorI() +KineGeneratorWithCache::KineGeneratorWithCache() : +EventRecordVisitorI(), fSafetyFactor(1.), fNumOfSafetyFactors(-1), fNumOfInterpolatorTypes(-1) { } //___________________________________________________________________________ -KineGeneratorWithCache::KineGeneratorWithCache(string name) : -EventRecordVisitorI(name) +KineGeneratorWithCache::KineGeneratorWithCache(string name) : +EventRecordVisitorI(name), fSafetyFactor(1.), fNumOfSafetyFactors(-1), fNumOfInterpolatorTypes(-1) { } //___________________________________________________________________________ -KineGeneratorWithCache::KineGeneratorWithCache(string name, string config) : -EventRecordVisitorI(name, config) +KineGeneratorWithCache::KineGeneratorWithCache(string name, string config) : +EventRecordVisitorI(name, config), fSafetyFactor(1.), fNumOfSafetyFactors(-1), fNumOfInterpolatorTypes(-1) { } @@ -54,30 +54,38 @@ KineGeneratorWithCache::~KineGeneratorWithCache() } //___________________________________________________________________________ -double KineGeneratorWithCache::MaxXSec(GHepRecord * event_rec) const +double KineGeneratorWithCache::MaxXSec(GHepRecord * event_rec, const int nkey) const { LOG("Kinematics", pINFO) - << "Getting max. differential xsec for the rejection method"; + << "Getting max. for the rejection method"; double xsec_max = -1; Interaction * interaction = event_rec->Summary(); LOG("Kinematics", pINFO) - << "Attempting to find a cached max{dxsec/dK} value"; - xsec_max = this->FindMaxXSec(interaction); - if(xsec_max>0) return xsec_max; + << "Attempting to find a cached max value"; + xsec_max = this->FindMaxXSec(interaction, nkey); + if(xsec_max>0) return nkey<=fNumOfSafetyFactors-1?vSafetyFactors[nkey]*xsec_max:xsec_max; LOG("Kinematics", pINFO) - << "Attempting to compute the max{dxsec/dK} value"; - xsec_max = this->ComputeMaxXSec(interaction); + << "Attempting to compute the max value"; + if (nkey == 0) + { + xsec_max = this->ComputeMaxXSec(interaction); + } + else + { + xsec_max = this->ComputeMaxXSec(interaction, nkey); + } + if(xsec_max>0) { - LOG("Kinematics", pINFO) << "max{dxsec/dK} = " << xsec_max; - this->CacheMaxXSec(interaction, xsec_max); - return xsec_max; + LOG("Kinematics", pINFO) << "max = " << xsec_max; + this->CacheMaxXSec(interaction, xsec_max, nkey); + return nkey<=fNumOfSafetyFactors-1?vSafetyFactors[nkey]*xsec_max:xsec_max; } LOG("Kinematics", pNOTICE) - << "Can not generate event kinematics {K} (max_xsec({K};E)<=0)"; + << "Can not generate event kinematics max_xsec<=0)"; // xsec for selected kinematics = 0 event_rec->SetDiffXSec(0,kPSNull); // switch on error flag @@ -87,7 +95,7 @@ double KineGeneratorWithCache::MaxXSec(GHepRecord * event_rec) const interaction->ResetBit(kISkipKinematicChk); // throw exception genie::exceptions::EVGThreadException exception; - exception.SetReason("kinematics generation: max_xsec({K};E)<=0"); + exception.SetReason("kinematics generation: max_xsec<=0"); exception.SwitchOnFastForward(); throw exception; @@ -95,7 +103,7 @@ double KineGeneratorWithCache::MaxXSec(GHepRecord * event_rec) const } //___________________________________________________________________________ double KineGeneratorWithCache::FindMaxXSec( - const Interaction * interaction) const + const Interaction * interaction, const int nkey) const { // Find a cached max xsec for the specified xsec algorithm & interaction and // close to the specified energy @@ -111,7 +119,7 @@ double KineGeneratorWithCache::FindMaxXSec( } // access the the cache branch - CacheBranchFx * cb = this->AccessCacheBranch(interaction); + CacheBranchFx * cb = this->AccessCacheBranch(interaction, nkey); // if there are enough points stored in the cache buffer to build a // spline, then intepolate @@ -119,7 +127,7 @@ double KineGeneratorWithCache::FindMaxXSec( if( E >= cb->Spl()->XMin() && E <= cb->Spl()->XMax()) { double spl_max_xsec = cb->Spl()->Evaluate(E); LOG("Kinematics", pINFO) - << "\nInterpolated: max xsec (E=" << E << ") = " << spl_max_xsec; + << "\nInterpolated: max (E=" << E << ") = " << spl_max_xsec; return spl_max_xsec; } LOG("Kinematics", pINFO) @@ -180,22 +188,24 @@ double KineGeneratorWithCache::FindMaxXSec( } //___________________________________________________________________________ void KineGeneratorWithCache::CacheMaxXSec( - const Interaction * interaction, double max_xsec) const + const Interaction * interaction, double max_xsec, const int nkey) const { LOG("Kinematics", pINFO) - << "Adding the computed max{dxsec/dK} value to cache"; - CacheBranchFx * cb = this->AccessCacheBranch(interaction); + << "Adding the computed max value to cache"; + CacheBranchFx * cb = this->AccessCacheBranch(interaction, nkey); double E = this->Energy(interaction); + if (E0) cb->AddValues(E,max_xsec); if(! cb->Spl() ) { - if( cb->Map().size() > 40 ) cb->CreateSpline(); + if( cb->Map().size() > 40 ) + cb->CreateSpline(nkey<=fNumOfInterpolatorTypes-1?vInterpolatorTypes[nkey]:""); } if( cb->Spl() ) { if( E < cb->Spl()->XMin() || E > cb->Spl()->XMax() ) { - cb->CreateSpline(); + cb->CreateSpline(nkey<=fNumOfInterpolatorTypes-1?vInterpolatorTypes[nkey]:""); } } } @@ -212,26 +222,26 @@ double KineGeneratorWithCache::Energy(const Interaction * interaction) const } //___________________________________________________________________________ CacheBranchFx * KineGeneratorWithCache::AccessCacheBranch( - const Interaction * interaction) const + const Interaction * interaction, const int nkey) const { // Returns the cache branch for this algorithm and this interaction. If no // branch is found then one is created. Cache * cache = Cache::Instance(); - // build the cache branch key as: namespace::algorithm/config/interaction + // build the cache branch key as: namespace::algorithm/config/interaction/nkey string algkey = this->Id().Key(); string intkey = interaction->AsString(); - string key = cache->CacheBranchKey(algkey, intkey); + string key = cache->CacheBranchKey(algkey, intkey, std::to_string(nkey)); CacheBranchFx * cache_branch = dynamic_cast (cache->FindCacheBranch(key)); if(!cache_branch) { //-- create the cache branch at the first pass - LOG("Kinematics", pINFO) << "No Max d^nXSec/d{K}^n cache branch found"; + LOG("Kinematics", pINFO) << "No cache branch found"; LOG("Kinematics", pINFO) << "Creating cache branch - key = " << key; - cache_branch = new CacheBranchFx("max[d^nXSec/d^n{K}] over phase space"); + cache_branch = new CacheBranchFx("Max over phase space"); cache->AddCacheBranch(key, cache_branch); } assert(cache_branch); @@ -270,3 +280,15 @@ void KineGeneratorWithCache::AssertXSecLimits( } } //___________________________________________________________________________ +double KineGeneratorWithCache::ComputeMaxXSec (const Interaction * in, const int nkey) const +{ + if (nkey == 0) + { + return this->ComputeMaxXSec(in); + } + else + { + return -1; + } +} +//___________________________________________________________________________ diff --git a/src/Physics/Common/KineGeneratorWithCache.h b/src/Physics/Common/KineGeneratorWithCache.h index a939d1dda..f839236c7 100644 --- a/src/Physics/Common/KineGeneratorWithCache.h +++ b/src/Physics/Common/KineGeneratorWithCache.h @@ -11,10 +11,17 @@ (retrieving, creating, searching, adding to) the cache. The various super-classes should implement the ComputeMaxXSec(...) method for computing the maximum xsec in case it has not already - being pushed into the cache at a previous iteration. + being pushed into the cache at a previous iteration. \n + + Update May 15, 2022 IK: + It makes possible to cache several values having different keys. + The example of using this opportunity see in + the class QELEventGeneratorSM. \author Costas Andreopoulos - University of Liverpool & STFC Rutherford Appleton Laboratory + University of Liverpool & STFC Rutherford Appleton Laboratory \n + Igor Kakorin + Joint Institute for Nuclear Research \created December 15, 2004 @@ -48,21 +55,26 @@ class KineGeneratorWithCache : public EventRecordVisitorI { ~KineGeneratorWithCache(); virtual double ComputeMaxXSec (const Interaction * in) const = 0; - virtual double MaxXSec (GHepRecord * evrec) const; - virtual double FindMaxXSec (const Interaction * in) const; - virtual void CacheMaxXSec (const Interaction * in, double xsec) const; + virtual double ComputeMaxXSec (const Interaction * in, const int nkey) const; + virtual double MaxXSec (GHepRecord * evrec, const int nkey=0) const; + virtual double FindMaxXSec (const Interaction * in, const int nkey=0) const; + virtual void CacheMaxXSec (const Interaction * in, double xsec, const int nkey=0) const; virtual double Energy (const Interaction * in) const; - virtual CacheBranchFx * AccessCacheBranch (const Interaction * in) const; + virtual CacheBranchFx * AccessCacheBranch (const Interaction * in, const int nkey=0) const; virtual void AssertXSecLimits (const Interaction * in, double xsec, double xsec_max) const; mutable const XSecAlgorithmI * fXSecModel; - double fSafetyFactor; ///< maxxsec -> maxxsec * safety_factor - double fMaxXSecDiffTolerance; ///< max{100*(xsec-maxxsec)/.5*(xsec+maxxsec)} if xsec>maxxsec - double fEMin; ///< min E for which maxxsec is cached - forcing explicit calc. - bool fGenerateUniformly; ///< uniform over allowed phase space + event weight? + double fSafetyFactor; ///< ComputeMaxXSec -> ComputeMaxXSec * fSafetyFactor + std::vector vSafetyFactors; ///< MaxXSec -> MaxXSec * fSafetyFactors[nkey] + int fNumOfSafetyFactors; ///< Number of given safety factors + std::vector vInterpolatorTypes; ///< Type of interpolator for each key in a branch + int fNumOfInterpolatorTypes; ///< Number of given interpolators types + double fMaxXSecDiffTolerance; ///< max{100*(xsec-maxxsec)/.5*(xsec+maxxsec)} if xsec>maxxsec + double fEMin; ///< min E for which maxxsec is cached - forcing explicit calc. + bool fGenerateUniformly; ///< uniform over allowed phase space + event weight? }; } // genie namespace diff --git a/src/Physics/Common/LinkDef.h b/src/Physics/Common/LinkDef.h index 8aba58b88..7f12724a5 100644 --- a/src/Physics/Common/LinkDef.h +++ b/src/Physics/Common/LinkDef.h @@ -23,5 +23,9 @@ #pragma link C++ class genie::XSecScaleMap; #pragma link C++ class genie::XSecLinearCombinations; #pragma link C++ class genie::QvalueShifter; +#pragma link C++ class genie::NormXSec; +#pragma link C++ class genie::NormGenerator; +#pragma link C++ class genie::NormInteractionListGenerator; + #endif diff --git a/src/Physics/Common/NormGenerator.cxx b/src/Physics/Common/NormGenerator.cxx new file mode 100644 index 000000000..1107c298f --- /dev/null +++ b/src/Physics/Common/NormGenerator.cxx @@ -0,0 +1,100 @@ +//____________________________________________________________________________ +/* + Copyright (c) 2003-2022, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org + + Igor Kakorin + Joint Institute for Nuclear Research +*/ +//____________________________________________________________________________ + +#include + +#include + + +#include "Framework/Algorithm/AlgConfigPool.h" +#include "Framework/Conventions/GBuild.h" +#include "Framework/EventGen/EventGeneratorI.h" +#include "Framework/EventGen//RunningThreadInfo.h" +#include "Framework/GHEP/GHepRecord.h" +#include "Framework/GHEP/GHepParticle.h" +#include "Physics/Common/NormGenerator.h" + +using namespace genie; + +//___________________________________________________________________________ +NormGenerator::NormGenerator(): +EventRecordVisitorI("genie::NormGenerator") +{ + +} +//___________________________________________________________________________ +NormGenerator::NormGenerator(string config): +EventRecordVisitorI("genie::NormGenerator") +{ + +} +//___________________________________________________________________________ +NormGenerator::~NormGenerator() +{ + +} +//___________________________________________________________________________ +void NormGenerator::ProcessEventRecord(GHepRecord * evrec) const +{ + Interaction * interaction = evrec->Summary(); + const InitialState & init_state = interaction -> InitState(); + + // Access cross section algorithm for running thread + RunningThreadInfo * rtinfo = RunningThreadInfo::Instance(); + const EventGeneratorI * evg = rtinfo->RunningThread(); + fXSecModel = evg->CrossSectionAlg(); + + //bool isHeavyNucleus = tgt->A()>=3; + GHepParticle * probe = evrec->Probe(); + int pdg_probe = probe->Pdg(); + int iprobe = evrec->ProbePosition(); + TLorentzVector * p4v = probe->GetP4(); + TLorentzVector * vtx = probe->X4(); + + GHepParticle * nucltgt = evrec->TargetNucleus(); + int pdg_tgt = nucltgt->Pdg(); + int inucltgt = evrec->TargetNucleusPosition(); + TLorentzVector * p4_tgt = nucltgt->GetP4(); + TLorentzVector * vtx_tgt = nucltgt->X4(); + + + evrec->AddParticle(pdg_probe, kIStStableFinalState, iprobe,-1,-1,-1, *p4v, *vtx); + evrec->AddParticle(pdg_tgt, kIStStableFinalState, inucltgt,-1,-1,-1, *p4_tgt, *vtx_tgt); + + //-- Determine the status code + //const Target & tgt = interaction->InitState().Tgt(); + + + // update the interaction summary + evrec->Summary()->KinePtr()->SetFSLeptonP4(*p4v); + + double xsec = fXSecModel->XSec(interaction, kPSfE); + + evrec->SetDiffXSec(xsec, kPSfE); + + return; +} +//___________________________________________________________________________ +void NormGenerator::Configure(const Registry & config) +{ + Algorithm::Configure(config); + this->LoadConfig(); +} +//____________________________________________________________________________ +void NormGenerator::Configure(string config) +{ + Algorithm::Configure(config); + this->LoadConfig(); +} +//____________________________________________________________________________ +void NormGenerator::LoadConfig(void) +{ +} +//____________________________________________________________________________ diff --git a/src/Physics/Common/NormGenerator.h b/src/Physics/Common/NormGenerator.h new file mode 100644 index 000000000..49cc5388e --- /dev/null +++ b/src/Physics/Common/NormGenerator.h @@ -0,0 +1,58 @@ +//____________________________________________________________________________ +/*! + +\class genie::NormGenerator + +\brief Normalization channel. Its main property is a constant cross section + per nucleon over the whole energy range. For nu/charged probes this produces + NC/EM events with the probe & target "echoed" back as final state particles. + +\ref [1] GENIE docdb 297 + + +\author Igor Kakorin + Joint Institute for Nuclear Research + +\created May 16, 2022 + +\cpright Copyright (c) 2003-2022, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org + +*/ +//____________________________________________________________________________ + +#ifndef _NORM_GENERATOR_H_ +#define _NORM_GENERATOR_H_ + +#include "Framework/EventGen/XSecAlgorithmI.h" +#include "Framework/EventGen/EventRecordVisitorI.h" + + +namespace genie { + + class NormGenerator: public EventRecordVisitorI { + + public : + NormGenerator(); + NormGenerator(string config); + ~NormGenerator(); + + // implement the EventRecordVisitorI interface + void ProcessEventRecord(GHepRecord * event_rec) const; + + // overload the Algorithm::Configure() methods to load private data + // members from configuration options + void Configure(const Registry & config); + void Configure(string config); + + // methods to load sub-algorithms and config data from the Registry + void LoadConfig (void); + + + + private: + mutable const XSecAlgorithmI * fXSecModel; + }; + +} // genie namespace +#endif // _NORM_GENERATOR_H_ diff --git a/src/Physics/Common/NormInteractionListGenerator.cxx b/src/Physics/Common/NormInteractionListGenerator.cxx new file mode 100644 index 000000000..3c4d8b6c5 --- /dev/null +++ b/src/Physics/Common/NormInteractionListGenerator.cxx @@ -0,0 +1,85 @@ +//____________________________________________________________________________ +/* + Copyright (c) 2003-2022, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org + + Igor Kakorin + Joint Institute for Nuclear Research +*/ +//____________________________________________________________________________ + +#include + +#include "Physics/Common/NormInteractionListGenerator.h" +#include "Framework/EventGen/InteractionList.h" +#include "Framework/Interaction/Interaction.h" +#include "Framework/ParticleData/PDGCodes.h" + +using namespace genie; + +//___________________________________________________________________________ +NormInteractionListGenerator::NormInteractionListGenerator() : +InteractionListGeneratorI("genie::NormInteractionListGenerator") +{ + +} +//___________________________________________________________________________ +NormInteractionListGenerator::NormInteractionListGenerator(string config) : +InteractionListGeneratorI("genie::NormInteractionListGenerator", config) +{ + +} +//___________________________________________________________________________ +NormInteractionListGenerator::~NormInteractionListGenerator() +{ + +} +//___________________________________________________________________________ +InteractionList * NormInteractionListGenerator::CreateInteractionList( + const InitialState & init_state) const +{ + InteractionList * intlist = new InteractionList; + int probe = init_state.ProbePdg(); + + bool isNC = probe == kPdgNuE || probe == kPdgAntiNuE || + probe == kPdgNuMu || probe == kPdgAntiNuMu || + probe == kPdgNuTau || probe == kPdgAntiNuTau; + + bool isEM = probe == kPdgElectron || probe == kPdgPositron || + probe == kPdgMuon || probe == kPdgAntiMuon || + probe == kPdgTau || probe == kPdgAntiTau; + + if (!isNC && !isEM) + { + LOG("IntLst", pWARN) + << "Unknown InteractionType! Returning NULL InteractionList " + << "for init-state: " << init_state.AsString(); + return 0; + } + + ProcessInfo proc_info(kScNorm, isNC?kIntWeakNC:kIntEM); + Interaction * interaction = new Interaction(init_state, proc_info); + + + intlist->push_back(interaction); + + return intlist; +} +//___________________________________________________________________________ +void NormInteractionListGenerator::Configure(const Registry & config) +{ + Algorithm::Configure(config); + this->LoadConfigData(); +} +//____________________________________________________________________________ +void NormInteractionListGenerator::Configure(string config) +{ + Algorithm::Configure(config); + this->LoadConfigData(); +} +//____________________________________________________________________________ +void NormInteractionListGenerator::LoadConfigData(void) +{ + +} +//____________________________________________________________________________ diff --git a/src/Physics/Common/NormInteractionListGenerator.h b/src/Physics/Common/NormInteractionListGenerator.h new file mode 100644 index 000000000..b251325ff --- /dev/null +++ b/src/Physics/Common/NormInteractionListGenerator.h @@ -0,0 +1,53 @@ +//____________________________________________________________________________ +/*! + +\class genie::NormInteractionListGenerator + +\brief Concrete implementations of the InteractionListGeneratorI interface. + Generates a list of all the Interaction (= event summary) objects that + can be generated by the NormGenerator. + +\author Igor Kakorin + Joint Institute for Nuclear Research + +\created May 16, 2022 + +\cpright Copyright (c) 2003-2022, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org + +*/ +//____________________________________________________________________________ + +#ifndef _NORM_INTERACTION_LIST_GENERATOR_H_ +#define _NORM_INTERACTION_LIST_GENERATOR_H_ + +#include "Framework/EventGen/InteractionListGeneratorI.h" + +namespace genie { + +class NormInteractionListGenerator : public InteractionListGeneratorI { + +public : + + NormInteractionListGenerator(); + NormInteractionListGenerator(string config); + ~NormInteractionListGenerator(); + + //-- implement the InteractionListGeneratorI interface + InteractionList * CreateInteractionList(const InitialState & init) const; + + //-- overload the Algorithm::Configure() methods to load private data + // members from configuration options + void Configure(const Registry & config); + void Configure(string config); + +private: + + void LoadConfigData(void); + + +}; + +} // genie namespace + +#endif // _NORM_INTERACTION_LIST_GENERATOR_H_ diff --git a/src/Physics/Common/NormXSec.cxx b/src/Physics/Common/NormXSec.cxx new file mode 100644 index 000000000..aee4fa5ad --- /dev/null +++ b/src/Physics/Common/NormXSec.cxx @@ -0,0 +1,74 @@ +//____________________________________________________________________________ +/* + Copyright (c) 2003-2022, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org + + Igor Kakorin + Joint Institute for Nuclear Research +*/ +//____________________________________________________________________________ + +#include "Framework/Algorithm/AlgConfigPool.h" +#include "Framework/Conventions/GBuild.h" +#include "Physics/Common/NormXSec.h" + +using namespace genie; +//____________________________________________________________________________ +NormXSec::NormXSec() : +XSecAlgorithmI("genie::NormXSec") +{ + +} +//____________________________________________________________________________ +NormXSec::NormXSec(string config) : +XSecAlgorithmI("genie::NormXSec", config) +{ + +} +//____________________________________________________________________________ +NormXSec::~NormXSec() +{ + +} +//____________________________________________________________________________ +double NormXSec::XSec( + const Interaction * interaction, KinePhaseSpace_t kps) const +{ + const Target & tgt = interaction->InitState().Tgt(); + double A = tgt.A(); + return A*fNormScale*1e-11; +} +//_________________________________ +double NormXSec::Integral(const Interaction * interaction) const +{ + const Target & tgt = interaction->InitState().Tgt(); + double A = tgt.A(); + return A*fNormScale*1e-11; +} +//____________________________________________________________________________ +bool NormXSec::ValidProcess(const Interaction * interaction) const +{ + return true; +} +//____________________________________________________________________________ +bool NormXSec::ValidKinematics(const Interaction * interaction) const +{ + return true; +} +//____________________________________________________________________________ +void NormXSec::Configure(const Registry & config) +{ + Algorithm::Configure(config); + this->LoadConfig(); +} +//____________________________________________________________________________ +void NormXSec::Configure(string config) +{ + Algorithm::Configure(config); + this->LoadConfig(); +} +//____________________________________________________________________________ +void NormXSec::LoadConfig(void) +{ + GetParamDef( "NormScale", fNormScale, 1.0); +} diff --git a/src/Physics/Common/NormXSec.h b/src/Physics/Common/NormXSec.h new file mode 100644 index 000000000..f1ad12d94 --- /dev/null +++ b/src/Physics/Common/NormXSec.h @@ -0,0 +1,56 @@ +//____________________________________________________________________________ +/*! + +\class genie::NormXSec + +\brief Normalization channel. Its main property is a constant cross section + per nucleon over the whole energy range. + +\ref [1] GENIE docdb 297 + + +\author Igor Kakorin + Joint Institute for Nuclear Research + +\created May 16, 2022 + +\cpright Copyright (c) 2003-2022, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org + + +*/ +//____________________________________________________________________________ + +#ifndef _NORM_CROSS_SECTION_H_ +#define _NORM_CROSS_SECTION_H_ +#include "Framework/EventGen/XSecAlgorithmI.h" + +namespace genie { + + +class NormXSec : public XSecAlgorithmI { + +public: + NormXSec(); + NormXSec(string config); + virtual ~NormXSec(); + + // XSecAlgorithmI interface implementation + double XSec (const Interaction * i, KinePhaseSpace_t k) const; + double Integral (const Interaction * i) const; + bool ValidProcess (const Interaction * i) const; + bool ValidKinematics (const Interaction * i) const; + + // Override the Algorithm::Configure methods to load configuration + // data to private data members + void Configure (const Registry & config); + void Configure (string param_set); + +private: + void LoadConfig (void); + double fNormScale; +}; + +} // genie namespace + +#endif diff --git a/src/Physics/QuasiElastic/EventGen/QELEventGeneratorSM.cxx b/src/Physics/QuasiElastic/EventGen/QELEventGeneratorSM.cxx index d2848dc4c..072f19807 100644 --- a/src/Physics/QuasiElastic/EventGen/QELEventGeneratorSM.cxx +++ b/src/Physics/QuasiElastic/EventGen/QELEventGeneratorSM.cxx @@ -12,9 +12,6 @@ Joint Institute for Nuclear Research, Institute for Theoretical and Experimental Physics - Vladimir Lyubushkin, - Joint Institute for Nuclear Research - Vadim Naumov , Joint Institute for Nuclear Research @@ -25,6 +22,8 @@ //____________________________________________________________________________ #include +#include +#include #include "Framework/Algorithm/AlgFactory.h" #include "Framework/Conventions/GBuild.h" @@ -55,9 +54,12 @@ using namespace genie; using namespace genie::controls; using namespace genie::utils; +using namespace genie::utils::gsl; namespace { // anonymous namespace (file only visibility) const double eps = std::numeric_limits::epsilon(); + const double max_dbl = std::numeric_limits::max(); + const double min_dbl = std::numeric_limits::min(); } //___________________________________________________________________________ QELEventGeneratorSM::QELEventGeneratorSM() : @@ -85,7 +87,6 @@ void QELEventGeneratorSM::ProcessEventRecord(GHepRecord * evrec) const LOG("QELEvent", pNOTICE) << "Generating kinematics uniformly over the allowed phase space"; } - // Get the interaction and set the 'trust' bits Interaction * interaction = evrec->Summary(); const InitialState & init_state = interaction -> InitState(); @@ -115,29 +116,30 @@ void QELEventGeneratorSM::ProcessEventRecord(GHepRecord * evrec) const const EventGeneratorI * evg = rtinfo->RunningThread(); fXSecModel = evg->CrossSectionAlg(); - // heavy nucleus is nucleus that heavier than hydrogen and deuterium - bool isHeavyNucleus = tgt->A()>=3; + // heavy nucleus is nucleus that heavier than tritium or 3He. + bool isHeavyNucleus = tgt->A()>3; sm_utils->SetInteraction(interaction); // phase space for heavy nucleus is different from light one - fkps = isHeavyNucleus?kPSQ2vfE:kPSQ2fE; + fkps = isHeavyNucleus?kPSQ2vpfE:kPSQ2fE; Range1D_t rQ2 = sm_utils->Q2QES_SM_lim(); // Try to calculate the maximum cross-section in kinematical limits // if not pre-computed already - double xsec_max1 = (fGenerateUniformly) ? -1 : this->MaxXSec(evrec); - double xsec_max2 = (fGenerateUniformly) ? -1 : (rQ2.maxMaxXSec2(evrec);// this make correct calculation of probability - double vmax= isHeavyNucleus?this->MaxDiffv(evrec) : 0.; + double xsec_max1 = fGenerateUniformly ? -1 : this->MaxXSec(evrec); + // this make correct calculation of probability + double xsec_max2 = fGenerateUniformly ? -1 : (rQ2.maxMaxXSec(evrec, 1); + double dvmax= isHeavyNucleus ? this->MaxXSec(evrec, 2) : 0.; - // generate Q2, v - double gQ2, v, xsec; + // generate Q2, v, pF + double Q2, v, kF, xsec; unsigned int iter = 0; bool accept = false; TLorentzVector q; while(1) { LOG("QELEvent", pINFO) << "Attempt #: " << iter; - if(iter > kRjMaxIterations) + if(iter > 100*kRjMaxIterations) { LOG("QELEvent", pWARN) << "Couldn't select a valid kinematics after " << iter << " iterations"; @@ -148,47 +150,52 @@ void QELEventGeneratorSM::ProcessEventRecord(GHepRecord * evrec) const throw exception; } - // Pick Q2 and v + // Pick Q2, v and pF double xsec_max = 0.; double pth = rnd->RndKine().Rndm(); //pth < prob1/(prob1+prob2), where prob1,prob2 - probabilities to generate event in area1 (Q2fQ2Min) which are not normalized if (pth <= xsec_max1*(TMath::Min(rQ2.max, fQ2Min)-rQ2.min)/(xsec_max1*(TMath::Min(rQ2.max, fQ2Min)-rQ2.min)+xsec_max2*(rQ2.max-fQ2Min))) { xsec_max = xsec_max1; - gQ2 = (rnd->RndKine().Rndm() * (TMath::Min(rQ2.max, fQ2Min)-rQ2.min)) + rQ2.min; + Q2 = (rnd->RndKine().Rndm() * (TMath::Min(rQ2.max, fQ2Min)-rQ2.min)) + rQ2.min; } else { xsec_max = xsec_max2; - gQ2 = (rnd->RndKine().Rndm() * (rQ2.max-fQ2Min)) + fQ2Min; + Q2 = (rnd->RndKine().Rndm() * (rQ2.max-fQ2Min)) + fQ2Min; } - Range1D_t rv = sm_utils->vQES_SM_lim(gQ2); + Range1D_t rv = sm_utils->vQES_SM_lim(Q2); // for nuclei heavier than deuterium generate energy transfer in defined energy interval - v = 0.; + double dv = 0.; if (isHeavyNucleus) { - v = vmax * rnd->RndKine().Rndm(); - if (v > (rv.max-rv.min)) + dv = dvmax * rnd->RndKine().Rndm(); + if (dv > (rv.max-rv.min)) { continue; } } - v += rv.min; + v = rv.min + dv; + + Range1D_t rkF = sm_utils->kFQES_SM_lim(Q2, v); + // rkF.max = Fermi momentum + kF = rnd->RndKine().Rndm()*sm_utils->GetFermiMomentum(); + if (kF < rkF.min) + { + continue; + } Kinematics * kinematics = interaction->KinePtr(); - kinematics->SetKV(kKVQ2, gQ2); + kinematics->SetKV(kKVQ2, Q2); kinematics->SetKV(kKVv, v); + kinematics->SetKV(kKVPn, kF); xsec = fXSecModel->XSec(interaction, fkps); - - //-- Decide whether to accept the current kinematics + //-- Decide whether to accept the current kinematics if(!fGenerateUniformly) { this->AssertXSecLimits(interaction, xsec, xsec_max); double t = xsec_max * rnd->RndKine().Rndm(); -#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ - LOG("QELEvent", pDEBUG)<< "xsec= " << xsec << ", J= " << J << ", Rnd= " << t; -#endif accept = (t < xsec); } else @@ -207,12 +214,9 @@ void QELEventGeneratorSM::ProcessEventRecord(GHepRecord * evrec) const } // Z-frame is frame where momentum transfer is (v,0,0,qv) - double qv = TMath::Sqrt(v*v+gQ2); + double qv = TMath::Sqrt(v*v+Q2); TLorentzVector transferMom(0, 0, qv, v); - Range1D_t rkF = sm_utils->kFQES_SM_lim(gQ2, v); - double kF = (rnd->RndKine().Rndm() * (rkF.max-rkF.min)) + rkF.min; - // Momentum of initial neutrino in LAB frame TLorentzVector * tempTLV = evrec->Probe()->GetP4(); TLorentzVector neutrinoMom = *tempTLV; @@ -245,7 +249,6 @@ void QELEventGeneratorSM::ProcessEventRecord(GHepRecord * evrec) const // 4-momentum of final nucleon in Z-frame TLorentzVector outNucleonMom = inNucleonMom+transferMom; - // Rotate all vectors from Z frame to LAB frame TVector3 yvec (0.0, 1.0, 0.0); TVector3 zvec (0.0, 0.0, 1.0); @@ -273,20 +276,20 @@ void QELEventGeneratorSM::ProcessEventRecord(GHepRecord * evrec) const int rpdgc = interaction->RecoilNucleonPdg(); assert(rpdgc); - double gW = PDGLibrary::Instance()->Find(rpdgc)->Mass(); - LOG("QELEvent", pNOTICE) << "Selected: W = "<< gW; + double W = PDGLibrary::Instance()->Find(rpdgc)->Mass(); + LOG("QELEvent", pNOTICE) << "Selected: W = "<< W; double M = init_state.Tgt().HitNucP4().M(); double E = init_state.ProbeE(kRfHitNucRest); // (W,Q2) -> (x,y) - double gx=0, gy=0; - kinematics::WQ2toXY(E,M,gW,gQ2,gx,gy); + double x=0, y=0; + kinematics::WQ2toXY(E,M,W,Q2,x,y); // lock selected kinematics & clear running values - interaction->KinePtr()->SetQ2(gQ2, true); - interaction->KinePtr()->SetW (gW, true); - interaction->KinePtr()->Setx (gx, true); - interaction->KinePtr()->Sety (gy, true); + interaction->KinePtr()->SetQ2(Q2, true); + interaction->KinePtr()->SetW (W, true); + interaction->KinePtr()->Setx (x, true); + interaction->KinePtr()->Sety (y, true); interaction->KinePtr()->ClearRunningValues(); // set the cross section for the selected kinematics @@ -315,7 +318,7 @@ void QELEventGeneratorSM::ProcessEventRecord(GHepRecord * evrec) const if (!fGenerateNucleonInNucleus) ist = kIStStableFinalState; else - ist = (tgt->IsNucleus()) ? kIStHadronInTheNucleus : kIStStableFinalState; + ist = (tgt->IsNucleus() && isHeavyNucleus) ? kIStHadronInTheNucleus : kIStStableFinalState; GHepParticle outNucleon(interaction->RecoilNucleonPdg(), ist, evrec->HitNucleonPosition(),-1,-1,-1, outNucleonMom , x4l); evrec->AddParticle(outNucleon); @@ -330,7 +333,6 @@ void QELEventGeneratorSM::ProcessEventRecord(GHepRecord * evrec) const // add a recoiled nucleus remnant this->AddTargetNucleusRemnant(evrec); - LOG("QELEvent", pINFO) << "Done generating QE event kinematics!"; } //___________________________________________________________________________ @@ -423,9 +425,12 @@ void QELEventGeneratorSM::Configure(string config) void QELEventGeneratorSM::LoadConfig(void) { - // Safety factor for the maximum differential cross section - GetParamDef( "MaxXSec-SafetyFactor", fSafetyFactor, 1.2) ; - GetParamDef( "MaxDiffv-SafetyFactor",fSafetyFacor_nu, 1.2); + // Safety factors for the maximum differential cross section + fNumOfSafetyFactors = GetParamVect("Safety-Factors", vSafetyFactors, false); + + // Interpolator types for the maximum differential cross section + fNumOfInterpolatorTypes = GetParamVect("Interpolator-Types", vInterpolatorTypes, false); + // Minimum energy for which max xsec would be cached, forcing explicit // calculation for lower eneries @@ -451,374 +456,289 @@ void QELEventGeneratorSM::LoadConfig(void) double QELEventGeneratorSM::ComputeMaxXSec(const Interaction * interaction) const { double xsec_max = -1; - const int N_Q2 = 8; - const int N_v = 8; + double tmp_xsec_max = -1; + double Q20, v0; + const int N_Q2 = 32; + const InitialState & init_state = interaction -> InitState(); + Target * tgt = init_state.TgtPtr(); + bool isHeavyNucleus = tgt->A()>3; + int N_v = isHeavyNucleus?32:0; Range1D_t rQ2 = sm_utils->Q2QES_SM_lim(); const double logQ2min = TMath::Log(TMath::Max(rQ2.min, eps)); const double logQ2max = TMath::Log(TMath::Min(rQ2.max, fQ2Min)); - - double tmp_xsec_max = -1; - // Now scan through kinematical variables Q2,v - for (int Q2_n=0; Q2_n < N_Q2; Q2_n++) + Kinematics * kinematics = interaction->KinePtr(); + const double pFmax = sm_utils->GetFermiMomentum(); + // Now scan through kinematical variables Q2,v,kF + for (int Q2_n=0; Q2_n <= N_Q2; Q2_n++) { // Scan around Q2 - double Q2 = TMath::Exp(Q2_n * (logQ2max-logQ2min)/N_Q2 + logQ2min); + double Q2 = TMath::Exp(Q2_n*(logQ2max-logQ2min)/N_Q2 + logQ2min); + kinematics->SetKV(kKVQ2, Q2); Range1D_t rv = sm_utils->vQES_SM_lim(Q2); const double logvmin = TMath::Log(TMath::Max(rv.min, eps)); const double logvmax = TMath::Log(TMath::Max(rv.max, TMath::Max(rv.min, eps))); - for (int v_n=0; v_n < N_v; v_n++) + for (int v_n=0; v_n <= N_v; v_n++) { // Scan around v - double v = TMath::Exp(v_n * (logvmax-logvmin)/N_v + logvmin); - Kinematics * kinematics = interaction->KinePtr(); - kinematics->SetKV(kKVQ2, Q2); + double v = TMath::Exp(v_n*(logvmax-logvmin)/N_v + logvmin); kinematics->SetKV(kKVv, v); + kinematics->SetKV(kKVPn, pFmax); // Compute the QE cross section for the current kinematics double xs = fXSecModel->XSec(interaction, fkps); if (xs > tmp_xsec_max) - tmp_xsec_max = xs; + { + tmp_xsec_max = xs; + Q20 = Q2; + v0 = v; + } } // Done with v scan }// Done with Q2 scan - - xsec_max = tmp_xsec_max; - // Apply safety factor, since value retrieved from the cache might - // correspond to a slightly different value - xsec_max *= fSafetyFactor; - return xsec_max; - -} -//___________________________________________________________________________ -double QELEventGeneratorSM::ComputeMaxXSec2(const Interaction * interaction) const -{ - double xsec_max = -1; - const int N_Q2 = 8; - const int N_v = 8; - - Range1D_t rQ2 = sm_utils->Q2QES_SM_lim(); - if (rQ2.maxvQES_SM_min(Q2min, Q2max); + const double vmax = sm_utils->vQES_SM_max(Q2min, Q2max); + + ROOT::Math::Minimizer * min = ROOT::Math::Factory::CreateMinimizer("Minuit", "Minimize"); + ROOT::Math::IBaseFunctionMultiDim * f = isHeavyNucleus?static_cast(new d3XSecSM_dQ2dvdkF_E(fXSecModel, interaction, pFmax)): + static_cast(new d1XSecSM_dQ2_E(fXSecModel, interaction)); + min->SetFunction( *f ); + min->SetMaxFunctionCalls(10000); // for Minuit/Minuit2 + min->SetMaxIterations(10000); // for GSL + min->SetTolerance(0.001); + min->SetPrintLevel(0); + double step = 1e-7; + min->SetVariable(0, "Q2", Q20, step); + min->SetVariableLimits(0, Q2min, Q2max); + if (isHeavyNucleus) { - // Scan around Q2 - double Q2 = TMath::Exp(Q2_n * (logQ2max-logQ2min)/N_Q2 + logQ2min); - Range1D_t rv = sm_utils->vQES_SM_lim(Q2); - const double logvmin = TMath::Log(TMath::Max(rv.min, eps)); - const double logvmax = TMath::Log(TMath::Max(rv.max, TMath::Max(rv.min, eps))); - for (int v_n=0; v_n < N_v; v_n++) - { - // Scan around v - double v = TMath::Exp(v_n * (logvmax-logvmin)/N_v + logvmin); - Kinematics * kinematics = interaction->KinePtr(); - kinematics->SetKV(kKVQ2, Q2); - kinematics->SetKV(kKVv, v); - // Compute the QE cross section for the current kinematics - double xs = fXSecModel->XSec(interaction, fkps); - if (xs > tmp_xsec_max) - tmp_xsec_max = xs; - } // Done with v scan - }// Done with Q2 scan + min->SetVariable(1, "v", v0, step); + min->SetVariableLimits(1, vmin, vmax); + } + min->Minimize(); + xsec_max = -min->MinValue(); + if (tmp_xsec_max > xsec_max) + { + xsec_max = tmp_xsec_max; + } - xsec_max = tmp_xsec_max; - // Apply safety factor, since value retrieved from the cache might - // correspond to a slightly different value - xsec_max *= fSafetyFactor; return xsec_max; } //___________________________________________________________________________ -double QELEventGeneratorSM::MaxXSec2(GHepRecord * event_rec) const +double QELEventGeneratorSM::ComputeMaxXSec(const Interaction * interaction, const int nkey) const { - LOG("Kinematics", pINFO) - << "Getting max. differential xsec for the rejection method"; - - double xsec_max = -1; - Interaction * interaction = event_rec->Summary(); - - LOG("Kinematics", pINFO) - << "Attempting to find a cached max{dxsec/dK} value"; - xsec_max = this->FindMaxXSec2(interaction); - if(xsec_max>0) return xsec_max; - - LOG("Kinematics", pINFO) - << "Attempting to compute the max{dxsec/dK} value"; - xsec_max = this->ComputeMaxXSec2(interaction); - if(xsec_max>0) { - LOG("Kinematics", pINFO) << "max{dxsec/dK} = " << xsec_max; - this->CacheMaxXSec2(interaction, xsec_max); + switch (nkey){ + case 0: + return this->ComputeMaxXSec(interaction); + + case 1: + { + Range1D_t rQ2 = sm_utils->Q2QES_SM_lim(); + if (rQ2.max InitState(); + Target * tgt = init_state.TgtPtr(); + bool isHeavyNucleus = tgt->A()>3; + int N_v = isHeavyNucleus?32:0; + + const double logQ2min = TMath::Log(fQ2Min); + const double logQ2max = TMath::Log(rQ2.max); + Kinematics * kinematics = interaction->KinePtr(); + const double pFmax = sm_utils->GetFermiMomentum(); + // Now scan through kinematical variables Q2,v,kF + for (int Q2_n=0; Q2_n <= N_Q2; Q2_n++) + { + // Scan around Q2 + double Q2 = TMath::Exp(Q2_n*(logQ2max-logQ2min)/N_Q2 + logQ2min); + kinematics->SetKV(kKVQ2, Q2); + Range1D_t rv = sm_utils->vQES_SM_lim(Q2); + const double logvmin = TMath::Log(TMath::Max(rv.min, eps)); + const double logvmax = TMath::Log(TMath::Max(rv.max, TMath::Max(rv.min, eps))); + for (int v_n=0; v_n <= N_v; v_n++) + { + // Scan around v + double v = TMath::Exp(v_n*(logvmax-logvmin)/N_v + logvmin); + kinematics->SetKV(kKVv, v); + kinematics->SetKV(kKVPn, pFmax); + // Compute the QE cross section for the current kinematics + double xs = fXSecModel->XSec(interaction, fkps); + if (xs > tmp_xsec_max) + { + tmp_xsec_max = xs; + Q20 = Q2; + v0 = v; + } + } // Done with v scan + }// Done with Q2 scan + + const double Q2min = fQ2Min; + const double Q2max = rQ2.max; + const double vmin = sm_utils->vQES_SM_min(Q2min, Q2max); + const double vmax = sm_utils->vQES_SM_max(Q2min, Q2max); + ROOT::Math::Minimizer * min = ROOT::Math::Factory::CreateMinimizer("Minuit", "Minimize"); + ROOT::Math::IBaseFunctionMultiDim * f = isHeavyNucleus?static_cast(new d3XSecSM_dQ2dvdkF_E(fXSecModel, interaction, pFmax)): + static_cast(new d1XSecSM_dQ2_E(fXSecModel, interaction)); + min->SetFunction( *f ); + min->SetMaxFunctionCalls(10000); // for Minuit/Minuit2 + min->SetMaxIterations(10000); // for GSL + min->SetTolerance(0.001); + min->SetPrintLevel(0); + double step = 1e-7; + min->SetVariable(0, "Q2", Q20, step); + min->SetVariableLimits(0, Q2min, Q2max); + if (isHeavyNucleus) + { + min->SetVariable(1, "v", v0, step); + min->SetVariableLimits(1, vmin, vmax); + } + min->Minimize(); + xsec_max = -min->MinValue(); + if (tmp_xsec_max > xsec_max) + { + xsec_max = tmp_xsec_max; + } return xsec_max; } - - LOG("Kinematics", pNOTICE) - << "Can not generate event kinematics {K} (max_xsec({K};E)<=0)"; - // xsec for selected kinematics = 0 - event_rec->SetDiffXSec(0,kPSNull); - // switch on error flag - event_rec->EventFlags()->SetBitNumber(kKineGenErr, true); - // reset 'trust' bits - interaction->ResetBit(kISkipProcessChk); - interaction->ResetBit(kISkipKinematicChk); - // throw exception - genie::exceptions::EVGThreadException exception; - exception.SetReason("kinematics generation: max_xsec({K};E)<=0"); - exception.SwitchOnFastForward(); - throw exception; - - return 0; -} -//___________________________________________________________________________ -double QELEventGeneratorSM::FindMaxXSec2( - const Interaction * interaction) const -{ -// Find a cached max xsec for the specified xsec algorithm & interaction and -// close to the specified energy - - // get neutrino energy - double E = this->Energy(interaction); - LOG("Kinematics", pINFO) << "E = " << E; - - if(E < fEMin) { - LOG("Kinematics", pINFO) - << "Below minimum energy - Forcing explicit calculation"; - return -1.; - } - - // access the the cache branch - CacheBranchFx * cb = this->AccessCacheBranch2(interaction); - - // if there are enough points stored in the cache buffer to build a - // spline, then intepolate - if( cb->Spl() ) { - if( E >= cb->Spl()->XMin() && E <= cb->Spl()->XMax()) { - double spl_max_xsec = cb->Spl()->Evaluate(E); - LOG("Kinematics", pINFO) - << "\nInterpolated: max xsec (E=" << E << ") = " << spl_max_xsec; - return spl_max_xsec; + case 2: + { + double diffv_max = -1; + double tmp_diffv_max = -1; + const int N_Q2 = 100; + double Q20; + Range1D_t rQ2 = sm_utils->Q2QES_SM_lim(); + for (int Q2_n = 0; Q2_n<=N_Q2; Q2_n++) // Scan around Q2 + { + double Q2 = rQ2.min + 1.*Q2_n * (rQ2.max-rQ2.min)/N_Q2; + Range1D_t rv = sm_utils->vQES_SM_lim(Q2); + if (rv.max-rv.min > tmp_diffv_max) + { + tmp_diffv_max = rv.max-rv.min; + Q20 = Q2; + } + } // Done with Q2 scan + + ROOT::Math::Minimizer * min = ROOT::Math::Factory::CreateMinimizer("Minuit", "Minimize"); + ROOT::Math::IBaseFunctionMultiDim * f = new dv_dQ2_E(interaction); + min->SetFunction( *f ); + min->SetMaxFunctionCalls(10000); // for Minuit/Minuit2 + min->SetMaxIterations(10000); // for GSL + min->SetTolerance(0.001); + min->SetPrintLevel(0); + double step = 1e-7; + min->SetVariable(0, "Q2", Q20, step); + min->SetVariableLimits(0, rQ2.min, rQ2.max); + min->Minimize(); + diffv_max = -min->MinValue(); + + if (tmp_diffv_max > diffv_max) + { + diffv_max = tmp_diffv_max; } - LOG("Kinematics", pINFO) - << "Outside spline boundaries - Forcing explicit calculation"; - return -1.; + return diffv_max; } - - // if there are not enough points at the cache buffer to have a spline, - // look whether there is another point that is sufficiently close - double dE = TMath::Min(0.25, 0.05*E); - const map & fmap = cb->Map(); - map::const_iterator iter = fmap.lower_bound(E); - if(iter != fmap.end()) { - if(TMath::Abs(E - iter->first) < dE) return iter->second; + default: + return -1.; } - - return -1; - } //___________________________________________________________________________ -void QELEventGeneratorSM::CacheMaxXSec2( - const Interaction * interaction, double max_xsec) const +// GSL wrappers +//____________________________________________________________________________ +d3XSecSM_dQ2dvdkF_E::d3XSecSM_dQ2dvdkF_E( + const XSecAlgorithmI * m, + const Interaction * i, + double pF) : ROOT::Math::IBaseFunctionMultiDim(), + fModel(m), + fInteraction(i), + fpF(pF) { - LOG("Kinematics", pINFO) - << "Adding the computed max{dxsec/dK} value to cache"; - CacheBranchFx * cb = this->AccessCacheBranch2(interaction); - - double E = this->Energy(interaction); - if(max_xsec>0) cb->AddValues(E,max_xsec); - - if(! cb->Spl() ) { - if( cb->Map().size() > 40 ) cb->CreateSpline(); - } - - if( cb->Spl() ) { - if( E < cb->Spl()->XMin() || E > cb->Spl()->XMax() ) { - cb->CreateSpline(); - } - } } -//___________________________________________________________________________ -CacheBranchFx * QELEventGeneratorSM::AccessCacheBranch2( - const Interaction * interaction) const +d3XSecSM_dQ2dvdkF_E::~d3XSecSM_dQ2dvdkF_E() { -// Returns the cache branch for this algorithm and this interaction. If no -// branch is found then one is created. - - Cache * cache = Cache::Instance(); - - // build the cache branch key as: namespace::algorithm/config/interaction - string algkey = this->Id().Key(); - string intkey = interaction->AsString(); - string key = cache->CacheBranchKey(algkey, intkey, "2nd"); - - CacheBranchFx * cache_branch = - dynamic_cast (cache->FindCacheBranch(key)); - if(!cache_branch) { - //-- create the cache branch at the first pass - LOG("Kinematics", pINFO) << "No Max d^nXSec/d{K}^n cache branch found"; - LOG("Kinematics", pINFO) << "Creating cache branch - key = " << key; - - cache_branch = new CacheBranchFx("max[d^nXSec/d^n{K}] over phase space"); - cache->AddCacheBranch(key, cache_branch); - } - assert(cache_branch); - - return cache_branch; } -//___________________________________________________________________________ -double QELEventGeneratorSM::ComputeMaxDiffv(const Interaction *) const +unsigned int d3XSecSM_dQ2dvdkF_E::NDim(void) const { - double max_diffv = -1; - const int N_Q2 = 10; - - Range1D_t rQ2 = sm_utils->Q2QES_SM_lim(); - for (int Q2_n = 0; Q2_nvQES_SM_lim(Q2); - if (rv.max-rv.min > max_diffv) - max_diffv = rv.max-rv.min; - } // Done with Q2 scan - max_diffv *= fSafetyFactor; - return max_diffv; - + return 2; } -//___________________________________________________________________________ -double QELEventGeneratorSM::MaxDiffv(GHepRecord * event_rec) const +double d3XSecSM_dQ2dvdkF_E::DoEval(const double * xin) const { - LOG("Kinematics", pINFO) - << "Getting max. vmax(Q2)-vmin(Q2) for the rejection method"; - - double max_diffv = -1; - Interaction * interaction = event_rec->Summary(); - - LOG("Kinematics", pINFO) - << "Attempting to find a cached max{vmax(Q2)-vmin(Q2)} value"; - max_diffv = this->FindMaxDiffv(interaction); - if(max_diffv>0) return max_diffv; - - LOG("Kinematics", pINFO) - << "Attempting to compute the max{vmax(Q2)-vmin(Q2)} value"; - max_diffv = this->ComputeMaxDiffv(interaction); - if(max_diffv>0) { - LOG("Kinematics", pINFO) << "max{vmax(Q2)-vmin(Q2)} = " << max_diffv; - this->CacheMaxDiffv(interaction, max_diffv); - return max_diffv; - } - - LOG("Kinematics", pNOTICE) - << "Can not generate event kinematics (max{vmax(Q2)-vmin(Q2);E}<=0)"; - // xsec for selected kinematics = 0 - event_rec->SetDiffXSec(0,kPSNull); - // switch on error flag - event_rec->EventFlags()->SetBitNumber(kKineGenErr, true); - // reset 'trust' bits - interaction->ResetBit(kISkipProcessChk); - interaction->ResetBit(kISkipKinematicChk); - // throw exception - genie::exceptions::EVGThreadException exception; - exception.SetReason("kinematics generation: max{vmax(Q2)-vmin(Q2);E}<=0"); - exception.SwitchOnFastForward(); - throw exception; - - return 0; +// outputs: +// differential cross section +// + fInteraction->KinePtr()->SetKV(kKVQ2, xin[0]); + fInteraction->KinePtr()->SetKV(kKVv, xin[1]); + fInteraction->KinePtr()->SetKV(kKVPn, fpF); + double xsec = -fModel->XSec(fInteraction, kPSQ2vpfE); + return xsec; } -//___________________________________________________________________________ -double QELEventGeneratorSM::FindMaxDiffv(const Interaction * interaction) const +ROOT::Math::IBaseFunctionMultiDim * + d3XSecSM_dQ2dvdkF_E::Clone() const { -// Find a cached maximum of vmax(Q2)-vmin(Q2) for xsec algorithm & interaction and -// close to the specified energy - - // get neutrino energy - double E = this->Energy(interaction); - LOG("Kinematics", pINFO) << "E = " << E; - - if(E < fEMin) { - LOG("Kinematics", pINFO) - << "Below minimum energy - Forcing explicit calculation"; - return -1.; - } - - // access the the cache branch - CacheBranchFx * cb = this->AccessCacheBranchDiffv(interaction); - - // if there are enough points stored in the cache buffer to build a - // spline, then intepolate - if( cb->Spl() ) { - if( E >= cb->Spl()->XMin() && E <= cb->Spl()->XMax()) - { - double spl_maxdiffv = cb->Spl()->Evaluate(E); - LOG("Kinematics", pINFO) - << "\nInterpolated: max vmax(Q2)-vmin(Q2) (E=" << E << ") = " << spl_maxdiffv; - return spl_maxdiffv; - } - LOG("Kinematics", pINFO) - << "Outside spline boundaries - Forcing explicit calculation"; - return -1.; - } - - // if there are not enough points at the cache buffer to have a spline, - // look whether there is another point that is sufficiently close - double dE = TMath::Min(0.25, 0.05*E); - const map & fmap = cb->Map(); - map::const_iterator iter = fmap.lower_bound(E); - if(iter != fmap.end()) - { - if(TMath::Abs(E - iter->first) < dE) return iter->second; - } - - return -1; - + return new d3XSecSM_dQ2dvdkF_E(fModel, fInteraction, fpF); } -//___________________________________________________________________________ -void QELEventGeneratorSM::CacheMaxDiffv(const Interaction * interaction, double max_diffv) const +//____________________________________________________________________________ +d1XSecSM_dQ2_E::d1XSecSM_dQ2_E( + const XSecAlgorithmI * m, + const Interaction * i) : ROOT::Math::IBaseFunctionMultiDim(), + fModel(m), + fInteraction(i) { - LOG("Kinematics", pINFO) - << "Adding the computed max{vmax(Q2)-vmin(Q2)} value to cache"; - CacheBranchFx * cb = this->AccessCacheBranchDiffv(interaction); - - double E = this->Energy(interaction); - if(max_diffv>0) cb->AddValues(E,max_diffv); - - if(! cb->Spl() ) - { - if( cb->Map().size() > 40 ) cb->CreateSpline(); - } - - if( cb->Spl() ) - { - if( E < cb->Spl()->XMin() || E > cb->Spl()->XMax() ) - { - cb->CreateSpline(); - } - } } -//___________________________________________________________________________ -CacheBranchFx * QELEventGeneratorSM::AccessCacheBranchDiffv(const Interaction * interaction) const +d1XSecSM_dQ2_E::~d1XSecSM_dQ2_E() { -// Returns the cache branch for this algorithm and this interaction. If no -// branch is found then one is created. - - Cache * cache = Cache::Instance(); - - // build the cache branch key as: namespace::algorithm/config/interaction - string algkey = this->Id().Key(); - string intkey = interaction->AsString(); - string key = cache->CacheBranchKey(algkey, intkey, "diffv"); - - CacheBranchFx * cache_branch = - dynamic_cast (cache->FindCacheBranch(key)); - if(!cache_branch) - { - //-- create the cache branch at the first pass - LOG("Kinematics", pINFO) << "No Max vmax(Q2)-vmin(Q2) cache branch found"; - LOG("Kinematics", pINFO) << "Creating cache branch - key = " << key; - - cache_branch = new CacheBranchFx("max[vmax(Q2)-vmin(Q2)] over phase space"); - cache->AddCacheBranch(key, cache_branch); - } - assert(cache_branch); - - return cache_branch; } -//___________________________________________________________________________ +unsigned int d1XSecSM_dQ2_E::NDim(void) const +{ + return 1; +} +double d1XSecSM_dQ2_E::DoEval(const double * xin) const +{ +// outputs: +// differential cross section +// + fInteraction->KinePtr()->SetKV(kKVQ2, xin[0]); + double xsec = -fModel->XSec(fInteraction, kPSQ2fE); + return xsec; +} +ROOT::Math::IBaseFunctionMultiDim * + d1XSecSM_dQ2_E::Clone() const +{ + return new d1XSecSM_dQ2_E(fModel, fInteraction); +} +//____________________________________________________________________________ +dv_dQ2_E::dv_dQ2_E(const Interaction * i) : ROOT::Math::IBaseFunctionMultiDim(), + fInteraction(i) +{ + AlgFactory * algf = AlgFactory::Instance(); + sm_utils = const_cast(dynamic_cast(algf->GetAlgorithm("genie::SmithMonizUtils","Default"))); + sm_utils->SetInteraction(fInteraction); +} +dv_dQ2_E::~dv_dQ2_E() +{ +} +unsigned int dv_dQ2_E::NDim(void) const +{ + return 1; +} +double dv_dQ2_E::DoEval(const double * xin) const +{ +// outputs: +// differential cross section +// + double Q2 = xin[0]; + Range1D_t rv = sm_utils->vQES_SM_lim(Q2); + return rv.min-rv.max; +} +ROOT::Math::IBaseFunctionMultiDim * + dv_dQ2_E::Clone() const +{ + return new dv_dQ2_E(fInteraction); +} +//____________________________________________________________________________ diff --git a/src/Physics/QuasiElastic/EventGen/QELEventGeneratorSM.h b/src/Physics/QuasiElastic/EventGen/QELEventGeneratorSM.h index 7d4792204..97b23e7ae 100644 --- a/src/Physics/QuasiElastic/EventGen/QELEventGeneratorSM.h +++ b/src/Physics/QuasiElastic/EventGen/QELEventGeneratorSM.h @@ -19,9 +19,6 @@ Joint Institute for Nuclear Research, Institute for Theoretical and Experimental Physics \n - Vladimir Lyubushkin, \n - Joint Institute for Nuclear Research \n - Vadim Naumov , \n Joint Institute for Nuclear Research \n @@ -39,6 +36,8 @@ #ifndef _QEL_EVENT_GENERATORSM_H_ #define _QEL_EVENT_GENERATORSM_H_ +#include + #include "Physics/Common/KineGeneratorWithCache.h" #include "Framework/Conventions/Controls.h" #include "Physics/QuasiElastic/XSection/SmithMonizUtils.h" @@ -66,19 +65,10 @@ public : void LoadConfig (void); double ComputeMaxXSec(const Interaction * in) const; + double ComputeMaxXSec (const Interaction * in, const int nkey) const; void AddTargetNucleusRemnant (GHepRecord * evrec) const; ///< add a recoiled nucleus remnant - double ComputeMaxXSec2 (const Interaction * in) const; - double MaxXSec2 (GHepRecord * evrec) const; - double FindMaxXSec2 (const Interaction * in) const; - void CacheMaxXSec2 (const Interaction * in, double xsec) const; - CacheBranchFx * AccessCacheBranch2 (const Interaction * in) const; - - double ComputeMaxDiffv (const Interaction * in) const; - double MaxDiffv (GHepRecord * evrec) const; - double FindMaxDiffv (const Interaction * in) const; - void CacheMaxDiffv (const Interaction * in, double xsec) const; - CacheBranchFx * AccessCacheBranchDiffv (const Interaction * in) const; + mutable KinePhaseSpace_t fkps; @@ -86,11 +76,77 @@ public : bool fGenerateNucleonInNucleus; ///< generate struck nucleon in nucleus double fQ2Min; ///< Q2-threshold for seeking the second maximum - double fSafetyFacor_nu; - }; // class definition +class XSecAlgorithmI; +class Interaction; + +namespace utils { +namespace gsl { +//..................................................................................... +// +// genie::utils::gsl::d3XSecSM_dQ2dvdkF_E +// A 3-D cross section function: d3XSecSM_dQ2dvdkF_E = f(Q2, v, kF=fixed)|(fixed E) +// +class d3XSecSM_dQ2dvdkF_E: public ROOT::Math::IBaseFunctionMultiDim +{ +public: + d3XSecSM_dQ2dvdkF_E(const XSecAlgorithmI *, const Interaction *, double pF); + ~d3XSecSM_dQ2dvdkF_E(); + + // ROOT::Math::IBaseFunctionMultiDim interface + unsigned int NDim (void) const; + double DoEval (const double *) const; + ROOT::Math::IBaseFunctionMultiDim * Clone (void) const; + +private: + const XSecAlgorithmI * fModel; + const Interaction * fInteraction; + const double fpF; +}; +// +// genie::utils::gsl::d1XSecSM_dQ2_E +// A 1-D cross section function: d1XSecSM_dQ2_E = f(Q2)|(fixed E) +// +class d1XSecSM_dQ2_E: public ROOT::Math::IBaseFunctionMultiDim +{ +public: + d1XSecSM_dQ2_E(const XSecAlgorithmI *, const Interaction *); + ~d1XSecSM_dQ2_E(); + + // ROOT::Math::IBaseFunctionMultiDim interface + unsigned int NDim (void) const; + double DoEval (const double *) const; + ROOT::Math::IBaseFunctionMultiDim * Clone (void) const; + +private: + const XSecAlgorithmI * fModel; + const Interaction * fInteraction; +}; +// +// genie::utils::gsl::dv_dQ2_E=f(Q2)|(fixed E) +// A 1-D dependence of allowable \nu-range from Q2 +// +class dv_dQ2_E: public ROOT::Math::IBaseFunctionMultiDim +{ +public: + dv_dQ2_E(const Interaction *); + ~dv_dQ2_E(); + + // ROOT::Math::IBaseFunctionMultiDim interface + unsigned int NDim (void) const; + double DoEval (const double *) const; + ROOT::Math::IBaseFunctionMultiDim * Clone (void) const; + +private: + const Interaction * fInteraction; + mutable SmithMonizUtils * sm_utils; +}; +} // gsl namespace +} // utils namespace + + } // genie namespace #endif // _QEL_EVENT_GENERATORSM_H_ diff --git a/src/Physics/QuasiElastic/XSection/SmithMonizQELCCPXSec.cxx b/src/Physics/QuasiElastic/XSection/SmithMonizQELCCPXSec.cxx index 55fcec099..bf3baee6c 100644 --- a/src/Physics/QuasiElastic/XSection/SmithMonizQELCCPXSec.cxx +++ b/src/Physics/QuasiElastic/XSection/SmithMonizQELCCPXSec.cxx @@ -75,7 +75,7 @@ double SmithMonizQELCCPXSec::XSec( double xsec; // dimension of kine phase space std::string s = KinePhaseSpace::AsString(kps); - int kpsdim = 1 + std::count(s.begin(), s.end(), ','); + int kpsdim = s!="<|E>"?1 + std::count(s.begin(), s.begin()+s.find('}'), ','):0; if(!this -> ValidProcess (interaction) ) { @@ -98,9 +98,13 @@ double SmithMonizQELCCPXSec::XSec( xsec = this->d2sQES_dQ2dv_SM(interaction); } + if(kpsdim == 3) + { + xsec = this->d3sQES_dQ2dvdkF_SM(interaction); + } - // The algorithm computes d^1xsec/dQ2 or d^2xsec/dQ2dv + // The algorithm computes d^1xsec/dQ2, d^2xsec/dQ2dv or d^3xsec/dQ2dvdp // Check whether variable tranformation is needed if ( kps != kPSQ2fE && kps != kPSQ2vfE ) { @@ -109,6 +113,8 @@ double SmithMonizQELCCPXSec::XSec( J = utils::kinematics::Jacobian(interaction, kPSQ2fE, kps); else if (kpsdim == 2) J = utils::kinematics::Jacobian(interaction, kPSQ2vfE, kps); + else if (kpsdim == 3) + J = utils::kinematics::Jacobian(interaction, kPSQ2vpfE, kps); xsec *= J; } @@ -192,25 +198,73 @@ void SmithMonizQELCCPXSec::LoadConfig(void) //____________________________________________________________________________ double SmithMonizQELCCPXSec::d3sQES_dQ2dvdkF_SM(const Interaction * interaction) const { + // Assuming that variables E_nu, Q2, \nu and kF are within allowable kinematic region + // which are specified in methods: genie::utils::gsl::d2Xsec_dQ2dv::DoEval and QELEventGeneratorSM::ProcessEventRecord // Get kinematics & init-state parameters - const Kinematics & kinematics = interaction -> Kine(); - + const Kinematics & kinematics = interaction -> Kine(); + sm_utils->SetInteraction(interaction); + const InitialState & init_state = interaction -> InitState(); + const Target & target = init_state.Tgt(); + double E_nu = init_state.ProbeE(kRfLab); + double Q2 = kinematics.GetKV(kKVQ2); + double v = kinematics.GetKV(kKVv); double kF = kinematics.GetKV(kKVPn); double kkF = kF*kF; - double P_Fermi, E_nuBIN; + int nucl_pdg_ini = target.HitNucPdg(); + int nucl_pdg_fin = genie::pdg::SwitchProtonNeutron(nucl_pdg_ini); + + PDGLibrary * pdglib = PDGLibrary::Instance(); + TParticlePDG * nucl_fin = pdglib->Find( nucl_pdg_fin ); - E_nuBIN = sm_utils->GetBindingEnergy(); + double E_BIN = sm_utils->GetBindingEnergy(); + double m_ini = target.HitNucMass(); + double mm_ini = m_ini*m_ini; + double m_fin = nucl_fin -> Mass(); // Mass of final hadron or hadron system (GeV) + double mm_fin = m_fin*m_fin; + double m_tar = target.Mass(); // Mass of target nucleus (GeV) + double mm_tar = m_tar*m_tar; + + // One of the xsec terms changes sign for antineutrinos + bool is_neutrino = pdg::IsNeutrino(init_state.ProbePdg()); + int n_NT = (is_neutrino) ? +1 : -1; + + double E_p = TMath::Sqrt(mm_ini+kkF)-E_BIN; + //|\vec{q}| + double qqv = v*v+Q2; + double qv = TMath::Sqrt(qqv); + double cosT_p = ((v-E_BIN)*(2*E_p+v+E_BIN)-qqv+mm_ini-mm_fin)/(2*kF*qv); //\cos\theta_p + if (cosT_p < -1.0 || cosT_p > 1.0 ) + { + return 0.0; + } + + double pF = TMath::Sqrt(kkF+(2*kF*qv)*cosT_p+qqv); + + double E_lep = E_nu-v; + double m_lep = interaction->FSPrimLepton()->Mass(); + double mm_lep = m_lep*m_lep; + if (E_lep < m_lep) + { + return 0.0; + } + double P_lep = TMath::Sqrt(E_lep*E_lep-mm_lep); + double k6 = (Q2+mm_lep)/(2*E_nu); + double cosT_lep= (E_lep-k6)/P_lep; + if (cosT_lep < -1.0 || cosT_lep > 1.0 ) return 0.0; + + double cosT_k = (v+k6)/qv; + if (cosT_k < -1.0 || cosT_k > 1.0 ) return 0.0; - double E_p = TMath::Sqrt(fmm_ini+kkF)-E_nuBIN; - double cosT_p = ((fv-E_nuBIN)*(2*E_p+fv+E_nuBIN)-fqqv+fmm_ini-fmm_fin)/(2*kF*fqv); //\cos\theta_p - if (cosT_p < -1.0 || cosT_p > 1.0 ) return 0.0; - double pF = TMath::Sqrt(kkF+(2*kF*fqv)*cosT_p+fqqv); - double b2_flux = (E_p-kF*fcosT_k*cosT_p)*(E_p-kF*fcosT_k*cosT_p); - double c2_flux = kkF*(1-cosT_p*cosT_p)*(1-fcosT_k*fcosT_k); + double b2_flux = (E_p-kF*cosT_k*cosT_p)*(E_p-kF*cosT_k*cosT_p); + double c2_flux = kkF*(1-cosT_p*cosT_p)*(1-cosT_k*cosT_k); + + double k1 = fVud2*kNucleonMass2*kPi; + double k2 = mm_lep/(2*mm_tar); + double k7 = P_lep*cosT_lep; - P_Fermi = sm_utils->GetFermiMomentum(); + double P_Fermi = sm_utils->GetFermiMomentum(); double FV_SM = 4.0*TMath::Pi()/3*TMath::Power(P_Fermi, 3); - double factor = fk1*(fm_tar*kF/(FV_SM*fqv*TMath::Sqrt(b2_flux-c2_flux)))*SmithMonizUtils::rho(P_Fermi, 0.0, kF)*(1-SmithMonizUtils::rho(P_Fermi, 0.01, pF)); + double factor = k1*(m_tar*kF/(FV_SM*qv*TMath::Sqrt(b2_flux-c2_flux)))*SmithMonizUtils::rho(P_Fermi, 0.0, kF)*(1-SmithMonizUtils::rho(P_Fermi, 0.01, pF)); double a2 = kkF/kNucleonMass2; double a3 = a2*cosT_p*cosT_p; @@ -219,19 +273,36 @@ double SmithMonizQELCCPXSec::d3sQES_dQ2dvdkF_SM(const Interaction * interaction) double a4 = a7*a7; double a5 = 2*a7*a6; - double k3 = fv/fqv; - double k4 = (3*a3-a2)/fqqv; - double k5 = (a7-a6*k3)*fm_tar/kNucleonMass; - - double T_1 = 1.0*fW_1+(a2-a3)*0.5*fW_2; //Ref.[1], W_1 - double T_2 = ((a2-a3)*fQ2/(2*fqqv)+a4-k3*(a5-k3*a3))*fW_2; //Ref.[1], W_2 - double T_3 = k5*fW_3; //Ref.[1], W_8 - double T_4 = fmm_tar*(0.5*fW_2*k4+1.0*fW_4/kNucleonMass2+a6*fW_5/(kNucleonMass*fqv)); //Ref.[1], W_\alpha - double T_5 = k5*fW_5+fm_tar*(a5/fqv-fv*k4)*fW_2; - - double xsec = kGF2*factor*((fE_lep-fk7)*(T_1+fk2*T_4)/fm_tar+(fE_lep+fk7)*T_2/(2*fm_tar) - +fn_NT*T_3*((fE_nu+fE_lep)*(fE_lep-fk7)/(2*fmm_tar)-fk2)-fk2*T_5) - *(kMw2/(kMw2+fQ2))*(kMw2/(kMw2+fQ2))/fE_nu/kPi; + double k3 = v/qv; + double k4 = (3*a3-a2)/qqv; + double k5 = (a7-a6*k3)*m_tar/kNucleonMass; + + // Calculate the QEL form factors + fFormFactors.Calculate(interaction); + double F_V = fFormFactors.F1V(); + double F_M = fFormFactors.xiF2V(); + double F_A = fFormFactors.FA(); + double F_P = fFormFactors.Fp(); + double FF_V = F_V*F_V; + double FF_M = F_M*F_M; + double FF_A = F_A*F_A; + + double t = Q2/(4*kNucleonMass2); + double W_1 = FF_A*(1+t)+t*(F_V+F_M)*(F_V+F_M); //Ref.[1], \tilde{T}_1 + double W_2 = FF_A+FF_V+t*FF_M; //Ref.[1], \tilde{T}_2 + double W_3 =-2*F_A*(F_V+F_M); //Ref.[1], \tilde{T}_8 + double W_4 =-0.5*F_V*F_M-F_A*F_P+t*F_P*F_P-0.25*(1-t)*FF_M; //Ref.[1], \tilde{T}_\alpha + double W_5 = FF_V+t*FF_M+FF_A; + + double T_1 = 1.0*W_1+(a2-a3)*0.5*W_2; //Ref.[1], W_1 + double T_2 = ((a2-a3)*Q2/(2*qqv)+a4-k3*(a5-k3*a3))*W_2; //Ref.[1], W_2 + double T_3 = k5*W_3; //Ref.[1], W_8 + double T_4 = mm_tar*(0.5*W_2*k4+1.0*W_4/kNucleonMass2+a6*W_5/(kNucleonMass*qv)); //Ref.[1], W_\alpha + double T_5 = k5*W_5+m_tar*(a5/qv-v*k4)*W_2; + + double xsec = kGF2*factor*((E_lep-k7)*(T_1+k2*T_4)/m_tar+(E_lep+k7)*T_2/(2*m_tar) + +n_NT*T_3*((E_nu+E_lep)*(E_lep-k7)/(2*mm_tar)-k2)-k2*T_5) + *(kMw2/(kMw2+Q2))*(kMw2/(kMw2+Q2))/E_nu/kPi; return xsec; @@ -242,64 +313,20 @@ double SmithMonizQELCCPXSec::d2sQES_dQ2dv_SM(const Interaction * interaction) co Kinematics * kinematics = interaction -> KinePtr(); sm_utils->SetInteraction(interaction); const InitialState & init_state = interaction -> InitState(); - fE_nu = init_state.ProbeE(kRfLab); - if (fE_nu < sm_utils->E_nu_thr_SM()) return 0; - fQ2 = kinematics->GetKV(kKVQ2); - fv = kinematics->GetKV(kKVv); - Range1D_t rkF = sm_utils->kFQES_SM_lim(fQ2,fv); + // Assuming that the energy is greater of threshold. + // See condition in method SmithMonizQELCCXSec::Integrate + // interaction->InitState().ProbeE(kRfLab)E_nu_thr_SM() + // of SmithMonizQELCCXSec.cxx + // if (E_nu < sm_utils->E_nu_thr_SM()) return 0; + // Assuming that variables Q2 and \nu are within allowable kinematic region + // which are specified in method: genie::utils::gsl::d2Xsec_dQ2dv::DoEval + double Q2 = kinematics->GetKV(kKVQ2); + double v = kinematics->GetKV(kKVv); + Range1D_t rkF = sm_utils->kFQES_SM_lim(Q2,v); const Target & target = init_state.Tgt(); - PDGLibrary * pdglib = PDGLibrary::Instance(); - - // One of the xsec terms changes sign for antineutrinos - bool is_neutrino = pdg::IsNeutrino(init_state.ProbePdg()); - fn_NT = (is_neutrino) ? +1 : -1; - - int nucl_pdg_ini = target.HitNucPdg(); - double m_ini = target.HitNucMass(); - fmm_ini = TMath::Power(m_ini, 2); - int nucl_pdg_fin = genie::pdg::SwitchProtonNeutron(nucl_pdg_ini); - TParticlePDG * nucl_fin = pdglib->Find( nucl_pdg_fin ); - double m_fin = nucl_fin -> Mass(); // Mass of final hadron or hadron system (GeV) - fmm_fin = TMath::Power(m_fin, 2); - fm_tar = target.Mass(); // Mass of target nucleus (GeV) - fmm_tar = TMath::Power(fm_tar, 2); - - fE_lep = fE_nu-fv; - double m_lep = interaction->FSPrimLepton()->Mass(); - double mm_lep = m_lep*m_lep; - if (fE_lep < m_lep) return 0.0; - double P_lep = TMath::Sqrt(fE_lep*fE_lep-mm_lep); - double k6 = (fQ2+mm_lep)/(2*fE_nu); - double cosT_lep= (fE_lep-k6)/P_lep; - if (cosT_lep < -1.0 || cosT_lep > 1.0 ) return 0.0; - //|\vec{q}| - fqqv = fv*fv+fQ2; - fqv = TMath::Sqrt(fqqv); - fcosT_k = (fv+k6)/fqv; - if (fcosT_k < -1.0 || fcosT_k > 1.0 ) return 0.0; - - fk1 = fVud2*kNucleonMass2*kPi; - fk2 = mm_lep/(2*fmm_tar); - fk7 = P_lep*cosT_lep; - + - // Calculate the QEL form factors - fFormFactors.Calculate(interaction); - fF_V = fFormFactors.F1V(); - fF_M = fFormFactors.xiF2V(); - fF_A = fFormFactors.FA(); - fF_P = fFormFactors.Fp(); - fFF_V = fF_V*fF_V; - fFF_M = fF_M*fF_M; - fFF_A = fF_A*fF_A; - - double t = fQ2/(4*kNucleonMass2); - fW_1 = fFF_A*(1+t)+t*(fF_V+fF_M)*(fF_V+fF_M); //Ref.[1], \tilde{T}_1 - fW_2 = fFF_A+fFF_V+t*fFF_M; //Ref.[1], \tilde{T}_2 - fW_3 =-2*fF_A*(fF_V+fF_M); //Ref.[1], \tilde{T}_8 - fW_4 =-0.5*fF_V*fF_M-fF_A*fF_P+t*fF_P*fF_P-0.25*(1-t)*fFF_M; //Ref.[1], \tilde{T}_\alpha - fW_5 = fFF_V+t*fFF_M+fFF_A; // Gaussian quadratures integrate over Fermi momentum double R[48]= { 0.16276744849602969579e-1,0.48812985136049731112e-1, @@ -428,8 +455,9 @@ double SmithMonizQELCCPXSec::dsQES_dQ2_SM(const Interaction * interaction) const // Apply given scaling factor xsec *= fXSecScale; - // Deuterium and tritium is a special case - if (target.A()>1 && target.A()<4) + // Pauli-correction factor for deuterium, we formally apply this factor for He-3 and tritium, + // because RFG model is not applicable for them. + if (1 Joint Institute for Nuclear Research - Vladimir Lyubushkin - Joint Institute for Nuclear Research - Vadim Naumov Joint Institute for Nuclear Research @@ -70,7 +67,7 @@ double SmithMonizQELCCXSec::Integrate( const InitialState & init_state = in -> InitState(); const Target & target = init_state.Tgt(); - if (target.A()<3) + if (target.A()<4) { const KPhaseSpace & kps = in->PhaseSpace(); if(!kps.IsAboveThreshold()) { @@ -95,7 +92,10 @@ double SmithMonizQELCCXSec::Integrate( { Interaction * interaction = new Interaction(*in); sm_utils->SetInteraction(in); - if (interaction->InitState().ProbeE(kRfLab)E_nu_thr_SM()) return 0; + if (interaction->InitState().ProbeE(kRfLab)E_nu_thr_SM()) + { + return 0; + } interaction->SetBit(kISkipProcessChk); interaction->SetBit(kISkipKinematicChk); double xsec = 0; @@ -194,7 +194,6 @@ double genie::utils::gsl::d2Xsec_dQ2dv::DoEval(const double * xin) const // differential cross section [10^-38 cm^2] // - Range1D_t rQ2 = sm_utils->Q2QES_SM_lim(); double Q2 = (rQ2.max-rQ2.min)*xin[0]+rQ2.min; Range1D_t rv = sm_utils->vQES_SM_lim(Q2); double v = (rv.max-rv.min)*xin[1]+rv.min; diff --git a/src/Physics/QuasiElastic/XSection/SmithMonizQELCCXSec.h b/src/Physics/QuasiElastic/XSection/SmithMonizQELCCXSec.h index 77f6e39fa..107725cd8 100644 --- a/src/Physics/QuasiElastic/XSection/SmithMonizQELCCXSec.h +++ b/src/Physics/QuasiElastic/XSection/SmithMonizQELCCXSec.h @@ -14,9 +14,6 @@ Konstantin Kuzmin Joint Institute for Nuclear Research \n - Vladimir Lyubushkin - Joint Institute for Nuclear Research \n - Vadim Naumov Joint Institute for Nuclear Research \n diff --git a/src/Physics/QuasiElastic/XSection/SmithMonizUtils.cxx b/src/Physics/QuasiElastic/XSection/SmithMonizUtils.cxx index ca0f2c3df..c49d46419 100644 --- a/src/Physics/QuasiElastic/XSection/SmithMonizUtils.cxx +++ b/src/Physics/QuasiElastic/XSection/SmithMonizUtils.cxx @@ -22,9 +22,6 @@ University of Liverpool & STFC Rutherford Appleton Laboratory */ //____________________________________________________________________________ - -#include - #include #include #include @@ -83,9 +80,9 @@ void SmithMonizUtils::LoadConfig(void) GetParam( "RFG-UseParametrization", fUseParametrization); - // Load removal energy for specific nuclei from either the algorithm's + // load removal energy for specific nuclei from either the algorithm's // configuration file or the UserPhysicsOptions file. - // If none is used use Wapstra's semi-empirical formula. + // if none is used use Wapstra's semi-empirical formula. const std::string keyStart = "RFG-NucRemovalE@Pdg="; RgIMap entries = GetConfig().GetItemMap(); @@ -120,30 +117,34 @@ void SmithMonizUtils::SetInteraction(const Interaction * interaction) { fInteraction = interaction; - // Get kinematics & init-state parameters - // unused // const Kinematics & kinematics = interaction -> Kine(); + // get kinematics & init-state parameters const InitialState & init_state = interaction -> InitState(); const Target & target = init_state.Tgt(); PDGLibrary * pdglib = PDGLibrary::Instance(); - - E_nu = interaction->InitState().ProbeE(kRfLab); // Neutrino energy (GeV) + + // neutrino energy (GeV) + E_nu = interaction->InitState().ProbeE(kRfLab); assert(target.HitNucIsSet()); // get lepton&nuclear masses (init & final state nucleus) - m_lep = interaction->FSPrimLepton()->Mass(); // Mass of final charged lepton (GeV) + + // mass of final charged lepton (GeV) + m_lep = interaction->FSPrimLepton()->Mass(); mm_lep = TMath::Power(m_lep, 2); int nucl_pdg_ini = target.HitNucPdg(); m_ini = target.HitNucMass(); mm_ini = TMath::Power(m_ini, 2); int nucl_pdg_fin = genie::pdg::SwitchProtonNeutron(nucl_pdg_ini); TParticlePDG * nucl_fin = pdglib->Find( nucl_pdg_fin ); - m_fin = nucl_fin -> Mass(); // Mass of final hadron or hadron system (GeV) + // mass of final hadron or hadron system (GeV) + m_fin = nucl_fin -> Mass(); mm_fin = TMath::Power(m_fin, 2); - m_tar = target.Mass(); // Mass of target nucleus (GeV) + // mass of target nucleus (GeV) + m_tar = target.Mass(); mm_tar = TMath::Power(m_tar, 2); - // For hydrogen and deuterium RFG is not applied - if (target.A()<3) + // RFG is not applied for A<4 + if (target.A()<4) { E_BIN = P_Fermi = m_rnu = mm_rnu = 0; return; @@ -162,8 +163,8 @@ void SmithMonizUtils::SetInteraction(const Interaction * interaction) << pdg::IonPdgCode(Af, Zf) << " in PDGLibrary!"; exit(1); } - - m_rnu = nucl_f -> Mass(); // Mass of residual nucleus (GeV) + // mass of residual nucleus (GeV) + m_rnu = nucl_f -> Mass(); mm_rnu = TMath::Power(m_rnu, 2); int Z = target.Z(); @@ -171,26 +172,26 @@ void SmithMonizUtils::SetInteraction(const Interaction * interaction) int N = A-Z; - // Maximum value of Fermi momentum of target nucleon (GeV) + // maximum value of Fermi momentum of target nucleon (GeV) if (A < 6 || !fUseParametrization) { - //-- look up the Fermi momentum for this Target + // look up the Fermi momentum for this Target FermiMomentumTablePool * kftp = FermiMomentumTablePool::Instance(); const FermiMomentumTable * kft = kftp->GetTable(fKFTable); P_Fermi = kft->FindClosestKF(pdg::IonPdgCode(A, Z), nucl_pdg_ini); } else { - //-- define the Fermi momentum for this Target + // define the Fermi momentum for this Target // P_Fermi = utils::nuclear::FermiMomentumForIsoscalarNucleonParametrization(target); - //-- correct the Fermi momentum for the struck nucleon + // correct the Fermi momentum for the struck nucleon if(is_p) P_Fermi *= TMath::Power( 2.*Z/A, 1./3); else P_Fermi *= TMath::Power( 2.*N/A, 1./3); } - // Neutrino binding energy (GeV) + // neutrino binding energy (GeV) if (target.A() < 6 || !fUseParametrization) { map::const_iterator it = fNucRmvE.find(Z); @@ -211,30 +212,33 @@ double SmithMonizUtils::E_nu_thr_SM(void) const Func1D QEL_EnuMin_SM_(*this, &SmithMonizUtils::QEL_EnuMin_SM); - - const int MFC = 10000; // Maximum of function call - const double EPSABS = 0; + // maximum of function call + const int MFC = 10000; + const double EPSABS = 0.; const double EPSREL = 1.0e-08; - double E_min = ((m_lep + m_rnu + m_fin)*(m_lep + m_rnu + m_fin) - mm_tar)/(2*m_tar); - double Enu_2 = 5.0e+00; - double Enu_rf; + + // Energy threshold of scattering on nucleus (Eq. (1) of Ref. 3) + double E_min = ((m_lep + m_rnu + m_fin)*(m_lep + m_rnu + m_fin) - mm_tar)/2/m_tar; + + // Energy threshold of scattering on bound nucleon (Eq. (2) of Ref. 3) + double E_min2 = ((m_lep + m_fin)*(m_lep + m_fin)-mm_ini-E_BIN*E_BIN+2*E_BIN*TMath::Sqrt(mm_ini+P_Fermi*P_Fermi))/2/(TMath::Sqrt(mm_ini+P_Fermi*P_Fermi)-E_BIN+P_Fermi); + + E_min = TMath::Max(E_min, E_min2); + + // if nu_1>nu_max at E_min then we try to find energy in range (E_min, 50.) where nu_1=nu_max and put E_min equal to it. + // nu_1 -- minimal energy transfer for bound nucleon (see Eqs. (11) of Ref. 3), + // nu_max -- maximal energy transfer on nucleus (see Eq. (9) of Ref.3) if (QEL_EnuMin_SM(E_min) > 0) { // C++ analog of fortran function Enu_rf= DZEROX(E_min,Enu_2,EPS,MFC,QEL_EnuMin_SM,1) ROOT::Math::RootFinder rfgb(ROOT::Math::RootFinder::kGSL_BRENT); - rfgb.Solve(QEL_EnuMin_SM_, E_min, Enu_2, MFC, EPSABS, EPSREL); //convergence is reached using tolerance = 2 *( epsrel * abs(x) + epsabs) - Enu_rf = rfgb.Root(); - } - else - { - Enu_rf = -1.0e+01; - } - E_min = TMath::Max(E_min,Enu_rf); - if(E_min < 0) - { - E_min = 0; - LOG("SmithMoniz", pDEBUG) << "E_min = " << E_min; + //convergence is reached using tolerance = 2 *( epsrel * abs(x) + epsabs) + if ( rfgb.Solve(QEL_EnuMin_SM_, E_min, 50., MFC, EPSABS, EPSREL) ) + { + E_min = rfgb.Root(); + } } + E_min = TMath::Max(E_min, 0.); return E_min; } @@ -242,33 +246,39 @@ double SmithMonizUtils::E_nu_thr_SM(void) const // The auxiliary function for determining energy threshold double SmithMonizUtils::QEL_EnuMin_SM(double Enu) const { + // return the minimum of nu_1-nu_max as function of Q2 in range ( Q2_lim1(Enu), Q2_lim2(Enu) ) const double EPS = 1.0e-06; const double Delta= 1.0e-14; const double Precision = std::numeric_limits::epsilon(); - Enu_in = Enu; double s = 2*Enu*m_tar+mm_tar; double W2 = (m_rnu+m_fin)*(m_rnu+m_fin); - double c = 0.5*(W2+mm_lep-mm_tar*(W2-mm_lep)/s); - double sqrtD = TMath::Sqrt(TMath::Max(Precision,LambdaFUNCTION(1.0, mm_lep/s, W2/s))); - double tmp = 0.5*(s-mm_tar); - double Q2_lim1= tmp*(1.0-sqrtD)-c; - double Q2_lim2= tmp*(1.0+sqrtD)-c; + // neutrino energy in CMS (see Eq. (4) of Ref.3) + double E_nu_CM = (s-mm_tar)/2/TMath::Sqrt(s); + // final lepton energy and momentum at W2_min (see Eqs. (5) and (6) of Ref.3) + double E_l_CM = (s+mm_lep-W2)/2/TMath::Sqrt(s); + double P_l_CM = E_l_CM>m_lep?TMath::Sqrt(E_l_CM*E_l_CM-mm_lep):Precision; + // minimal and maximal allowed Q2 for scattering on nucleus (see Eqs. (3) of Ref.3) + double Q2_lim1 = 2*E_nu_CM*(E_l_CM-P_l_CM)-mm_lep; + double Q2_lim2 = 2*E_nu_CM*(E_l_CM+P_l_CM)-mm_lep; // C++ analog of fortran function DMINFC(Q2lim1_SM,Q2_lim1,Q2_lim2,EPS,Delta,Q2_0,F_MIN,LLM) - Func1D Q2lim1_SM_(*this, &SmithMonizUtils::Q2lim1_SM); + Func1Dpar Q2lim1_SM_(*this, &SmithMonizUtils::Q2lim1_SM, Enu); double Q2_0,F_MIN; bool LLM; + // find minimum of nu_1-nu_max as function of Q2 in range (Q2_lim1,Q2_lim2) DMINFC(Q2lim1_SM_,Q2_lim1,Q2_lim2,EPS,Delta,Q2_0,F_MIN,LLM); return F_MIN; } //____________________________________________________________________________ // The auxiliary function for determining Q2-range -double SmithMonizUtils::Q2lim1_SM(double Q2) const +double SmithMonizUtils::Q2lim1_SM(double Q2, double Enu) const { - double s = 2*Enu_in*m_tar+mm_tar; - double nu_max = (s-mm_tar-mm_lep*(s-mm_tar)/(Q2+mm_lep)-mm_tar*(Q2+mm_lep)/(s-mm_tar))/(2*m_tar); + // maximal energy transfer (see Eq. (9) of Ref.3) + double nu_max = Enu*Q2/(Q2+mm_lep)-(Q2+mm_lep)/4/Enu; + double E = sqrt(P_Fermi*P_Fermi+mm_ini); double b = (E-E_BIN)*(E-E_BIN)-P_Fermi*P_Fermi; double a = 0.5*(Q2+mm_fin-b); + // minimal energy transfer for bound nucleon (see Eqs. (11) of Ref. 3), double nu_1 = (a*(E-E_BIN)-P_Fermi*TMath::Sqrt(a*a+Q2*b))/b; return nu_1-nu_max; @@ -277,10 +287,13 @@ double SmithMonizUtils::Q2lim1_SM(double Q2) const // The auxiliary function for determining Q2-range double SmithMonizUtils::Q2lim2_SM(double Q2) const { + // minimal energy transfer for scattering on nucleus (see Eq. (7) of Ref.3) double nu_min = ((m_rnu+m_fin)*(m_rnu+m_fin)+Q2-mm_tar)/(2*m_tar); + double E = sqrt(P_Fermi*P_Fermi+mm_ini); double b = (E-E_BIN)*(E-E_BIN)-P_Fermi*P_Fermi; double a = (Q2+mm_fin-b)*0.5; + // maximal energy transfer for bound nucleon (see Eqs. (11) of Ref. 3) double nu_2 = (a*(E-E_BIN)+P_Fermi*TMath::Sqrt(a*a+Q2*b))/b; return nu_min-nu_2; @@ -290,76 +303,117 @@ double SmithMonizUtils::Q2lim2_SM(double Q2) const Range1D_t SmithMonizUtils::Q2QES_SM_lim(void) const { - - const int MFC = 1000; // Maximum of function call + // here we find Q2_min and Q2_max such that + // 0. Q2_min>=0 + // 1. nu_1(Q2_min)<=nu_max(Q2_min) + // 2. nu_1(Q2_max)<=nu_max(Q2_max) + // 3. nu_min(Q2_min)<=nu_2(Q2_min) + // 4. Q2_min>=the value of minimal Q2 defined for scattering on nucleus + // 5. Q2_max<=the value of maximal Q2 defined for scattering on nucleus (see Eqs. (3) of Ref.3) + // nu_1 -- minimal energy transfer for bound nucleon (see Eqs. (11) of Ref. 3), + // nu_max -- maximal energy transfer on nucleus (see Eq. (9) of Ref.3), + // nu_2 -- maximal energy transfer for bound nucleon (see Eqs. (11) of Ref. 3), + // nu_min -- minimal energy transfer on nucleus (see Eq. (7) of Ref.3) + + // maximum of function call + const int MFC = 1000; const double EPS = 1.0e-08; const double Delta= 1.0e-14; const double EPSABS = 0; const double EPSREL = 1.0e-08; - - Enu_in = E_nu; - double s = 2*E_nu*m_tar+mm_tar; - double W2 = (m_rnu+m_fin)*(m_rnu+m_fin); - double c = 0.5*(W2+mm_lep-mm_tar*(W2-mm_lep)/s); - double sqrtD = TMath::Sqrt(LambdaFUNCTION(1.0, mm_lep/s, W2/s)); - double tmp = 0.5*(s-mm_tar); - double Q2_min = tmp*(1.0-sqrtD)-c; - double Q2_max = tmp*(1.0+sqrtD)-c; - // if the nucleus is hydrogen or deuterium then there is no need in further calculation + const double Precision = std::numeric_limits::epsilon(); + // if the nucleus mass is less than 4 then this is a special case if (E_BIN == 0 && P_Fermi == 0) { + double s = 2*E_nu*m_ini+mm_ini; + // minimal W2 for scattering on nucleus (see Eq. (6) of Ref.3) + double W2 = mm_fin; + // neutrino energy in CMS (see Eq. (4) of Ref.3) + double E_nu_CM = (s-mm_ini)/2/TMath::Sqrt(s); + // final lepton energy and momentum at W2_min (see Eqs. (5) of Ref.3) + double E_l_CM = (s+mm_lep-W2)/2/TMath::Sqrt(s); + double P_l_CM = E_l_CM>m_lep?TMath::Sqrt(E_l_CM*E_l_CM-mm_lep):Precision; + // minimal and maximal allowed Q2 for scattering on nucleus (see Eqs. (3) of Ref.3) + double Q2_min = 2*E_nu_CM*(E_l_CM-P_l_CM)-mm_lep; + double Q2_max = 2*E_nu_CM*(E_l_CM+P_l_CM)-mm_lep; Q2_min= TMath::Max(Q2_min,0.0); Range1D_t R(Q2_min,Q2_max); return R; } + double s = 2*E_nu*m_tar+mm_tar; + // minimal W2 for scattering on nucleus (see Eq. (6) of Ref.3) + double W2 = (m_rnu+m_fin)*(m_rnu+m_fin); + // neutrino energy in CMS (see Eq. (4) of Ref.3) + double E_nu_CM = (s-mm_tar)/2/TMath::Sqrt(s); + // final lepton energy and momentum at W2_min (see Eqs. (5) of Ref.3) + double E_l_CM = (s+mm_lep-W2)/2/TMath::Sqrt(s); + double P_l_CM = E_l_CM>m_lep?TMath::Sqrt(E_l_CM*E_l_CM-mm_lep):Precision; + // minimal and maximal allowed Q2 for scattering on nucleus (see Eqs. (3) of Ref.3) + double Q2_min = 2*E_nu_CM*(E_l_CM-P_l_CM)-mm_lep; + double Q2_max = 2*E_nu_CM*(E_l_CM+P_l_CM)-mm_lep; double F_MIN, Q2_0; bool LLM; // C++ analog of fortran function DMINFC(Q2lim1_SM,Q2_min,Q2_max,EPS,Delta,Q2_0,F_MIN,LLM) - Func1D Q2lim1_SM_(*this, &SmithMonizUtils::Q2lim1_SM); + Func1Dpar Q2lim1_SM_(*this, &SmithMonizUtils::Q2lim1_SM, E_nu); + // if minimum of nu_1-nu_max>0 then exit with error, because it's impossible DMINFC(Q2lim1_SM_,Q2_min,Q2_max,EPS,Delta,Q2_0,F_MIN,LLM); if (F_MIN>0) { - LOG("SmithMoniz", pFATAL) - << "No overlapped area for energy " << E_nu << "\n" << - "Q2_min=" << Q2_min << " Q2_max=" << Q2_max << "\n" << - "Q2_0=" << Q2_0 << " F_MIN=" << F_MIN; - exit(1); + LOG("SmithMoniz", pFATAL) + << "No overlapped area for energy " << E_nu << "\n" << + "Q2_min=" << Q2_min << " Q2_max=" << Q2_max << "\n" << + "Q2_0=" << Q2_0 << " F_MIN=" << F_MIN; + exit(1); } - if (Q2lim1_SM(Q2_min)>0) + // at Q2_0 here we have: nu_1(Q2_0)-nu_max(Q2_0)<0 + // if nu_1(Q2_min)-nu_max(Q2_min)>0 we find corrected Q2_min_cor>Q2_min where nu_1(Q2_min_cor)-nu_max(Q2_min_cor)=0 + // (it is always possible because of above conditions) + if (Q2lim1_SM(Q2_min, E_nu)>0) { //C++ analog of fortran function Q2_RF=DZEROX(Q2_min,Q2_0,EPS,MFC,Q2lim1_SM,1) ROOT::Math::RootFinder rfgb(ROOT::Math::RootFinder::kGSL_BRENT); - rfgb.Solve(Q2lim1_SM_, Q2_min, Q2_0, MFC, EPSABS, EPSREL); //convergence is reached using tolerance = 2 *( epsrel * abs(x) + epsabs) - double Q2_RF = rfgb.Root(); - Q2_min= TMath::Max(Q2_min,Q2_RF); + // convergence is reached using tolerance = 2 *( epsrel * abs(x) + epsabs) + if (rfgb.Solve(Q2lim1_SM_, Q2_min, Q2_0, MFC, EPSABS, EPSREL)) + { + Q2_min= rfgb.Root(); + } } - if(Q2lim1_SM(Q2_max)>0) + // if nu_1(Q2_max)-nu_max(Q2_max)>0 we find Q2_max_cor0) { // C++ analog of fortran function Q2_RF=DZEROX(Q2_0,Q2_max,Eps,MFC,Q2lim1_SM,1) ROOT::Math::RootFinder rfgb(ROOT::Math::RootFinder::kGSL_BRENT); - rfgb.Solve(Q2lim1_SM_, Q2_0, Q2_max, MFC, EPSABS, EPSREL); //convergence is reached using tolerance = 2 *( epsrel * abs(x) + epsabs) - double Q2_RF = rfgb.Root(); - Q2_max= TMath::Min(Q2_max,Q2_RF); + //convergence is reached using tolerance = 2 *( epsrel * abs(x) + epsabs) + if (rfgb.Solve(Q2lim1_SM_, Q2_0, Q2_max, MFC, EPSABS, EPSREL)) + { + Q2_max= rfgb.Root(); + } } Func1D Q2lim2_SM_(*this, &SmithMonizUtils::Q2lim2_SM); + // if nu_min(Q2_min)-nu_2(Q2_min)>0 and nu_min(Q2_max)-nu_2(Q2_max)>0 then set Q2_min=Q2_max (it makes xsec equal to zero). if (Q2lim2_SM(Q2_min)>0) { if(Q2lim2_SM(Q2_max)>0) { LOG("SmithMoniz", pWARN) << "The RFG model is not applicable! The cross section is set zero!"; - Q2_min = Q2_max; + Q2_min = Q2_max; } + // here we have nu_min(Q2_min)-nu_2(Q2_min)>0 and nu_min(Q2_max)-nu_2(Q2_max)<0 or vice versa + // so we always can find Q2_min_cor where nu_min(Q2_min_cor)-nu_2(Q2_min_cor)=0 else { // C++ analog of fortran function Q2_RF = DZEROX(Q2_min,Q2_max,Eps,MFC,Q2lim2_SM,1) ROOT::Math::RootFinder rfgb(ROOT::Math::RootFinder::kGSL_BRENT); - rfgb.Solve(Q2lim2_SM_, Q2_min,Q2_max, MFC, EPSABS, EPSREL); //convergence is reached using tolerance = 2 *( epsrel * abs(x) + epsabs) - double Q2_RF = rfgb.Root(); - Q2_min= TMath::Max(Q2_min,Q2_RF); + // convergence is reached using tolerance = 2 *( epsrel * abs(x) + epsabs) + if (rfgb.Solve(Q2lim2_SM_, Q2_min,Q2_max, MFC, EPSABS, EPSREL)) + { + Q2_min= rfgb.Root(); + } } } + Q2_min = TMath::Max(Q2_min,0.0); - Q2_min= TMath::Max(Q2_min,0.0); Range1D_t R(Q2_min,Q2_max); return R; @@ -368,37 +422,118 @@ Range1D_t SmithMonizUtils::Q2QES_SM_lim(void) const // Return allowed v-range for given Q2 Range1D_t SmithMonizUtils::vQES_SM_lim(double Q2) const { - double s = 2*E_nu*m_tar+mm_tar; - double nu_min= ((m_rnu+m_fin)*(m_rnu+m_fin)+Q2-mm_tar)/(2*m_tar); - double nu_max= (s-mm_tar-mm_lep*(s-mm_tar)/(Q2+mm_lep)-mm_tar*(Q2+mm_lep)/(s-mm_tar))/(2*m_tar); + + // minimal energy transfer for scattering on nucleus (see Eq. (7) of Ref.3) + double nu_min= ((m_rnu+m_fin)*(m_rnu+m_fin)+Q2-mm_tar)/2/m_tar; + + // if the target is nucleon then nu_min=nu_max=(m_fin^2+Q^2-m_ini^2)/(2*m_ini) + if (E_BIN == 0 && P_Fermi == 0) + return Range1D_t(nu_min, nu_min); + + // maximal energy transfer (see Eq. (9) of Ref.3) + double nu_max = E_nu*Q2/(Q2+mm_lep)-(Q2+mm_lep)/4/E_nu; + + // now we find limits for bound nucleon double E = TMath::Sqrt(P_Fermi*P_Fermi+mm_ini); double b = (E-E_BIN)*(E-E_BIN)-P_Fermi*P_Fermi; double a = (Q2+mm_fin-b)*0.5; double tmp1 = a*(E-E_BIN); double tmp2 = P_Fermi*TMath::Sqrt(a*a+Q2*b); + // minimal and maximal energy transfer for bound nucleon (see Eqs. (11) of Ref. 3) double nu_1 = (tmp1-tmp2)/b; double nu_2 = (tmp1+tmp2)/b; + // for minimal energy transfer we take maximum of corresponding values on nucleus and bound nucleon nu_min= TMath::Max(nu_min,nu_1); + // for maximal energy transfer we take minimum of corresponding values on nucleus and bound nucleon nu_max= TMath::Min(nu_max,nu_2); - Range1D_t R; - if (E_BIN == 0 && P_Fermi == 0) - nu_max=nu_min; + if (nu_min<=nu_max) - R = Range1D_t(nu_min,nu_max); - else - R = Range1D_t(0.5*(nu_min+nu_max),0.5*(nu_min+nu_max)); - return R; + return Range1D_t(nu_min,nu_max); + else + // to avoid machine precision errors + return Range1D_t(0.5*(nu_min+nu_max),0.5*(nu_min+nu_max)); } //____________________________________________________________________________ +// Return minimum of low edge of v-range for given neutrino energy +double SmithMonizUtils::vQES_SM_min(double Q2_min, double Q2_max) const +{ + // maximum of function call + const double EPS = 1.0e-08; + const double Delta= 1.0e-14; + + if (E_BIN == 0 && P_Fermi == 0) + return vQES_SM_lim(Q2_min).min; + + double F_MIN, Q2_0; + bool LLM; + // C++ analog of fortran function DMINFC(Q2lim1_SM,Q2_min,Q2_max,EPS,Delta,Q2_0,F_MIN,LLM) + Func1D vlim1_SM_(*this, &SmithMonizUtils::vlim1_SM); + DMINFC(vlim1_SM_,Q2_min,Q2_max,EPS,Delta,Q2_0,F_MIN,LLM); + + return F_MIN; + +} +//____________________________________________________________________________ +// Return maximum of low edge of v-range for given neutrino energy +double SmithMonizUtils::vQES_SM_max(double Q2_min, double Q2_max) const +{ + // maximum of function call + const double EPS = 1.0e-08; + const double Delta= 1.0e-14; + + if (E_BIN == 0 && P_Fermi == 0) + return vQES_SM_lim(Q2_max).min; + + double F_MIN, Q2_0; + bool LLM; + // C++ analog of fortran function DMINFC(Q2lim1_SM,Q2_min,Q2_max,EPS,Delta,Q2_0,F_MIN,LLM) + Func1D vlim2_SM_(*this, &SmithMonizUtils::vlim2_SM); + DMINFC(vlim2_SM_,Q2_min,Q2_max,EPS,Delta,Q2_0,F_MIN,LLM); + + return -F_MIN; + +} +//____________________________________________________________________________ +// The auxiliary function for determining minimum of low edge of v-range +double SmithMonizUtils::vlim1_SM(double Q2) const +{ + // minimal energy transfer for scattering on nucleus (see Eq. (7) of Ref.3) + double nu_min = ((m_rnu+m_fin)*(m_rnu+m_fin)+Q2-mm_tar)/(2*m_tar); + + double E = sqrt(P_Fermi*P_Fermi+mm_ini); + double b = (E-E_BIN)*(E-E_BIN)-P_Fermi*P_Fermi; + double a = (Q2+mm_fin-b)*0.5; + // minimal energy transfer for bound nucleon (see Eqs. (11) of Ref. 3) + double nu_1 = (a*(E-E_BIN)-P_Fermi*TMath::Sqrt(a*a+Q2*b))/b; + nu_min= TMath::Max(nu_min,nu_1); + return nu_min; +} +//____________________________________________________________________________ +// The auxiliary function for determining maximum of up edge of v-range +double SmithMonizUtils::vlim2_SM(double Q2) const +{ + // maximal energy transfer (see Eq. (9) of Ref.3) + double nu_max = E_nu*Q2/(Q2+mm_lep)-(Q2+mm_lep)/4/E_nu; + + double E = sqrt(P_Fermi*P_Fermi+mm_ini); + double b = (E-E_BIN)*(E-E_BIN)-P_Fermi*P_Fermi; + double a = 0.5*(Q2+mm_fin-b); + // maximal energy transfer for bound nucleon (see Eqs. (11) of Ref. 3) + double nu_2 = (a*(E-E_BIN)+P_Fermi*TMath::Sqrt(a*a+Q2*b))/b; + nu_max= TMath::Min(nu_max,nu_2); + return -nu_max; +} +//____________________________________________________________________________ // Return allowed Fermi momentum range for given Q2 and v Range1D_t SmithMonizUtils::kFQES_SM_lim(double Q2, double nu) const { double qv = TMath::Sqrt(nu*nu+Q2); double c_f = (nu-E_BIN)/qv; double d_f = (E_BIN*E_BIN-2*nu*E_BIN-Q2+mm_ini-mm_fin)/(2*qv*m_ini); - double Ef_min= m_ini*(c_f*d_f+TMath::Sqrt(1.0-c_f*c_f+d_f*d_f))/(1.0-c_f*c_f); - double kF_min= TMath::Sqrt(TMath::Max(Ef_min*Ef_min-mm_ini,0.0)); + // minimal energy of initial nucleon (see Eq. (13) of Ref.3) + double Ef_min= TMath::Max(m_ini*(c_f*d_f+TMath::Sqrt(1.0-c_f*c_f+d_f*d_f))/(1.0-c_f*c_f), TMath::Sqrt(P_Fermi*P_Fermi+mm_ini)-nu); + double kF_min= P_Fermi!=0?TMath::Sqrt(TMath::Max(Ef_min*Ef_min-mm_ini,0.0)):0.; double kF_max= P_Fermi; Range1D_t R; if (kF_min<=kF_max) @@ -409,13 +544,8 @@ Range1D_t SmithMonizUtils::kFQES_SM_lim(double Q2, double nu) const } //____________________________________________________________________________ -double SmithMonizUtils::LambdaFUNCTION(double a, double b, double c) const -{ - return a*a+b*b+c*c-2*(a*b+b*c+a*c); -} -//____________________________________________________________________________ // C++ implementation of DMINC function from CERN library -void SmithMonizUtils::DMINFC(Func1D &F, double A,double B, double EPS, double DELTA, double &X, double &Y, bool &LLM) const +void SmithMonizUtils::DMINFC(Functor1D &F, double A,double B, double EPS, double DELTA, double &X, double &Y, bool &LLM) const { const double W5 = TMath::Sqrt(5); const double HV = (3-W5)/2; @@ -479,13 +609,19 @@ void SmithMonizUtils::DMINFC(Func1D &F, double A,double B, doub double SmithMonizUtils::rho(double P_Fermi, double T_Fermi, double p) { - if (T_Fermi==0) //Pure Fermi gaz with T_Fermi=0 + if (T_Fermi==0) + { + //Pure Fermi gaz with T_Fermi=0 if(p<=P_Fermi) return 1.0; else return 0.0; - else //Fermi-Dirac distribution + } + else + { + //Fermi-Dirac distribution return 1.0/(1.0 + TMath::Exp(-(P_Fermi-p)/T_Fermi)); + } } @@ -514,16 +650,6 @@ double SmithMonizUtils::GetTheta_p(double pv, double v, double qv, double &E_p) return 0; } //____________________________________________________________________________ -double SmithMonizUtils::vQES_SM_low_bound(double Q2) const -{ - return vQES_SM_lim(Q2).min; -} -//____________________________________________________________________________ -double SmithMonizUtils::vQES_SM_upper_bound(double Q2) const -{ - return -1.0*vQES_SM_lim(Q2).max; -} -//____________________________________________________________________________ double SmithMonizUtils::PhaseSpaceVolume(KinePhaseSpace_t ps) const { double vol = 0.0; @@ -533,17 +659,24 @@ double SmithMonizUtils::PhaseSpaceVolume(KinePhaseSpace_t ps) const vol = rQ2.max - rQ2.min; return vol; } - else if (ps == kPSQ2vfE) + else if (ps == kPSQ2vpfE) { - Range1D_t rQ2 = Q2QES_SM_lim(); const int kNQ2 = 101; + const int kNv = 101; + Range1D_t rQ2 = Q2QES_SM_lim(); double dQ2 = (rQ2.max-rQ2.min)/(kNQ2-1); - for(int iq2=0; iq2 Joint Institute for Nuclear Research \n @@ -65,15 +67,26 @@ class SmithMonizUtils : public Algorithm { Range1D_t kFQES_SM_lim(double nu, double Q2) const; static double rho(double P_Fermi, double T_Fermi, double p); double PhaseSpaceVolume(KinePhaseSpace_t ps) const; - + double vQES_SM_min(double Q2min, double Q2max) const; + double vQES_SM_max(double Q2min, double Q2max) const; //! methods overloading the Algorithm() interface implementation //! to build the fragmentation function from configuration data void Configure(const Registry & config); void Configure(string config); + +double QEL_EnuMin_SM(double E_nu) const; private: + class Functor1D + { + public: + virtual ~Functor1D() {} + virtual double operator()(double d) = 0; + }; + + template - class Func1D + class Func1D : public Functor1D { public: Func1D(const C &obj, double (C::*f)(double) const):obj_(obj), f_(f){} @@ -83,15 +96,27 @@ class SmithMonizUtils : public Algorithm { const C &obj_; double (C::*f_)(double) const; }; + + template + class Func1Dpar : public Functor1D + { + public: + Func1Dpar(const C &obj, double (C::*f)(double,double) const, double par):obj_(obj), f_(f), parameter(par){} + ~Func1Dpar(){} + double operator()(double d) {return (obj_.*f_)( d, parameter);} + private: + const C &obj_; + double (C::*f_)(double, double) const; + double parameter; + }; void LoadConfig (void); - double QEL_EnuMin_SM(double E_nu) const; - double Q2lim1_SM(double Q2) const; + // double QEL_EnuMin_SM(double E_nu) const; + double Q2lim1_SM(double Q2, double Enu) const; double Q2lim2_SM(double Q2) const; - double LambdaFUNCTION(double a, double b, double c) const; - void DMINFC(Func1D &F, double A,double B, double EPS, double DELTA, double &X, double &Y, bool &LLM) const; - double vQES_SM_low_bound (double Q2) const; - double vQES_SM_upper_bound(double Q2) const; + void DMINFC(Functor1D &F, double A,double B, double EPS, double DELTA, double &X, double &Y, bool &LLM) const; + double vlim1_SM(double Q2) const; + double vlim2_SM(double Q2) const; map fNucRmvE; string fKFTable; @@ -115,7 +140,6 @@ class SmithMonizUtils : public Algorithm { double mm_rnu; ///< Squared mass of residual nucleus (GeV) double P_Fermi; ///< Maximum value of Fermi momentum of target nucleon (GeV) double E_BIN; ///< Binding energy (GeV) -mutable double Enu_in; ///< Running neutrino energy (GeV) };