diff --git a/src/external/BMWtools/BMWIntegrator.h b/src/external/BMWtools/BMWIntegrator.h index a99c7092..370f8ca3 100644 --- a/src/external/BMWtools/BMWIntegrator.h +++ b/src/external/BMWtools/BMWIntegrator.h @@ -481,11 +481,13 @@ inline double TAnSWaveGapIntegral::FuncAtX(double *x) const // x = {E, phi}, fPa *

Class for the 1D integration of j0(a*x)*exp(-b*x) * The integration uses the GSL integration routines. */ -class TIntBesselJ0Exp : public TIntegrator { +class TIntBesselJ0Exp : public T2Integrator { public: TIntBesselJ0Exp() {} ~TIntBesselJ0Exp() {} - double FuncAtX(double) const; + + private: + double FuncAtX(double, const vector&) const; }; /** @@ -496,27 +498,29 @@ class TIntBesselJ0Exp : public TIntegrator { * * \param x point where the function should be evaluated */ -inline double TIntBesselJ0Exp::FuncAtX(double x) const +inline double TIntBesselJ0Exp::FuncAtX(double x, const vector &par) const { - double w0t(TMath::TwoPi()*fPar[0]*x); + double w0t(TMath::TwoPi()*par[0]*x); double j0; if (fabs(w0t) < 0.001) { // check zero time limits of the spherical bessel functions j0(x) and j1(x) j0 = 1.0; } else { j0 = TMath::Sin(w0t)/w0t; } - return j0 * TMath::Exp(-fPar[1]*x); + return j0 * TMath::Exp(-par[1]*x); } /** *

Class for the 1D integration of sin(a*x)*exp(-b*x*x) * The integration uses the GSL integration routines. */ -class TIntSinGss : public TIntegrator { +class TIntSinGss : public T2Integrator { public: TIntSinGss() {} ~TIntSinGss() {} - double FuncAtX(double) const; + + private: + double FuncAtX(double, const vector&) const; }; /** @@ -527,9 +531,9 @@ class TIntSinGss : public TIntegrator { * * \param x point where the function should be evaluated */ -inline double TIntSinGss::FuncAtX(double x) const +inline double TIntSinGss::FuncAtX(double x, const vector &par) const { - return TMath::Sin(TMath::TwoPi()*fPar[0]*x) * TMath::Exp(-0.5*fPar[1]*fPar[1]*x*x); + return TMath::Sin(TMath::TwoPi()*par[0]*x) * TMath::Exp(-0.5*par[1]*par[1]*x*x); } /** @@ -538,11 +542,13 @@ inline double TIntSinGss::FuncAtX(double x) const * doi:10.1016/S0921-4526(00)00368-9 * The integration uses the GSL integration routines. */ -class TIntSGInterpolation : public TIntegrator { +class TIntSGInterpolation : public T2Integrator { public: TIntSGInterpolation() {} ~TIntSGInterpolation() {} - double FuncAtX(double) const; + + private: + double FuncAtX(double, const vector&) const; }; /** @@ -553,12 +559,12 @@ class TIntSGInterpolation : public TIntegrator { * * \param x point where the function should be evaluated */ -inline double TIntSGInterpolation::FuncAtX(double x) const +inline double TIntSGInterpolation::FuncAtX(double x, const vector &par) const { // Parameters: nu_L [MHz], a [1/us], lambda [1/us], beta [1], t [us] - double wt(TMath::TwoPi()*fPar[0]*x); - double expo(0.5*fPar[1]*fPar[1]*x*x/fPar[3]+fPar[2]*fPar[4]); - return (wt*TMath::Cos(wt)-TMath::Sin(wt))/(wt*wt)*TMath::Exp(-TMath::Power(expo,fPar[3]))/TMath::Power(expo,(1.0-fPar[3])); + double wt(TMath::TwoPi()*par[0]*x); + double expo(0.5*par[1]*par[1]*x*x/par[3]+par[2]*par[4]); + return (wt*TMath::Cos(wt)-TMath::Sin(wt))/(wt*wt)*TMath::Exp(-TMath::Power(expo,par[3]))/TMath::Power(expo,(1.0-par[3])); } /** @@ -593,8 +599,10 @@ inline double TGapIntegral::FuncAtX(double e) const class TFirstUniaxialGssKTIntegral : public T2Integrator { public: TFirstUniaxialGssKTIntegral() {} - ~TFirstUniaxialGssKTIntegral() {} - double FuncAtX(double, const vector&) const; // variable: x + virtual ~TFirstUniaxialGssKTIntegral() {} + + private: + virtual double FuncAtX(double, const vector&) const; // variable: x }; /** @@ -621,8 +629,10 @@ inline double TFirstUniaxialGssKTIntegral::FuncAtX(double x, const vector&) const; // variable: x + virtual ~TSecondUniaxialGssKTIntegral() {} + + private: + virtual double FuncAtX(double, const vector&) const; // variable: x }; /** diff --git a/src/external/libLFRelaxation/TLFRelaxation.cpp b/src/external/libLFRelaxation/TLFRelaxation.cpp index 2da09e40..00c6bc65 100644 --- a/src/external/libLFRelaxation/TLFRelaxation.cpp +++ b/src/external/libLFRelaxation/TLFRelaxation.cpp @@ -66,9 +66,8 @@ ClassImp(TLFSGInterpolation) /** *

Constructor: Initialize variables and allocate memory for the integrator. */ -TLFStatGssKT::TLFStatGssKT() : fCalcNeeded(true), fLastTime(0.) { +TLFStatGssKT::TLFStatGssKT() { fIntSinGss = new TIntSinGss(); - fPar.resize(2, 0.); } //-------------------------------------------------------------------------- @@ -80,8 +79,6 @@ TLFStatGssKT::TLFStatGssKT() : fCalcNeeded(true), fLastTime(0.) { TLFStatGssKT::~TLFStatGssKT() { delete fIntSinGss; fIntSinGss = 0; - fIntValues.clear(); - fPar.clear(); } //-------------------------------------------------------------------------- @@ -103,39 +100,17 @@ double TLFStatGssKT::operator()(double t, const vector &par) const { double sigsq(par[1]*par[1]); double sigsqtt(sigsq*t*t); - if(TMath::Abs(par[0])<0.00135538817){ + if(TMath::Abs(par[0])<0.00135538817) { return (0.33333333333333333333+0.66666666666666666667*(1.0-sigsqtt)*TMath::Exp(-0.5*sigsqtt)); } - for (unsigned int i(0); i<2; ++i) { - if( fPar[i]-par[i] ) { - fPar[i] = par[i]; - fCalcNeeded = true; - } - } - - if (fCalcNeeded) { - fIntValues.clear(); - fIntValues[0.] = 0.; - fLastTime = 0.; - } else { - if (t > fIntValues.rbegin()->first) { - fCalcNeeded = true; - fLastTime = fIntValues.rbegin()->first; - } - } - double omegaL(TMath::TwoPi()*par[0]); double coeff1(2.0*sigsq/(omegaL*omegaL)); double coeff2(coeff1*sigsq/omegaL); - if (fCalcNeeded) { - fIntSinGss->SetParameters(par); - fIntValues[t] = fIntValues[fLastTime] + fIntSinGss->IntegrateFunc(fLastTime,t); - fCalcNeeded = false; - } + double intValue = fIntSinGss->IntegrateFunc(0.0, t, par); - return (1.0-(coeff1*(1.0-TMath::Exp(-0.5*sigsqtt)*TMath::Cos(omegaL*t)))+(coeff2*fIntValues[t])); + return (1.0-(coeff1*(1.0-TMath::Exp(-0.5*sigsqtt)*TMath::Cos(omegaL*t)))+(coeff2*intValue)); } //-------------------------------------------------------------------------- @@ -144,9 +119,8 @@ double TLFStatGssKT::operator()(double t, const vector &par) const { /** *

Constructor: Initialize variables and allocate memory for the integrator. */ -TLFStatExpKT::TLFStatExpKT() : fCalcNeeded(true), fLastTime(0.) { +TLFStatExpKT::TLFStatExpKT() { fIntBesselJ0Exp = new TIntBesselJ0Exp(); - fPar.resize(2, 0.); } //-------------------------------------------------------------------------- @@ -158,8 +132,6 @@ TLFStatExpKT::TLFStatExpKT() : fCalcNeeded(true), fLastTime(0.) { TLFStatExpKT::~TLFStatExpKT() { delete fIntBesselJ0Exp; fIntBesselJ0Exp = 0; - fIntValues.clear(); - fPar.clear(); } //-------------------------------------------------------------------------- @@ -178,28 +150,10 @@ double TLFStatExpKT::operator()(double t, const vector &par) const { if(t<0.0) return 1.0; - if(TMath::Abs(par[0])<0.00135538817){ + if(TMath::Abs(par[0])<0.00135538817) { return (0.33333333333333333333+0.66666666666666666667*(1.0-par[1]*t)*TMath::Exp(-par[1]*t)); } - for (unsigned int i(0); i<2; ++i) { - if( fPar[i]-par[i] ) { - fPar[i] = par[i]; - fCalcNeeded = true; - } - } - - if (fCalcNeeded) { - fIntValues.clear(); - fIntValues[0.] = 0.; - fLastTime = 0.; - } else { - if (t > fIntValues.rbegin()->first) { - fCalcNeeded = true; - fLastTime = fIntValues.rbegin()->first; - } - } - double omegaL(TMath::TwoPi()*par[0]); double coeff1(par[1]/omegaL); double coeff2(coeff1*coeff1); @@ -215,13 +169,9 @@ double TLFStatExpKT::operator()(double t, const vector &par) const { j1 = (TMath::Sin(w0t)-w0t*TMath::Cos(w0t))/(w0t*w0t); } - if (fCalcNeeded) { - fIntBesselJ0Exp->SetParameters(par); - fIntValues[t] = fIntValues[fLastTime] + fIntBesselJ0Exp->IntegrateFunc(fLastTime,t); - fCalcNeeded = false; - } + double intValue = fIntBesselJ0Exp->IntegrateFunc(0.0, t, par); - return (1.0-(coeff1*TMath::Exp(-par[1]*t)*j1)-(coeff2*(j0*TMath::Exp(-par[1]*t)-1.0))-coeff3*fIntValues[t]); + return (1.0-(coeff1*TMath::Exp(-par[1]*t)*j1)-(coeff2*(j0*TMath::Exp(-par[1]*t)-1.0))-coeff3*intValue); } //-------------------------------------------------------------------------- @@ -245,7 +195,7 @@ TLFDynGssKT::TLFDynGssKT() : fCalcNeeded(true), fFirstCall(true), fCounter(0) { #ifdef HAVE_GOMP fftwf_plan_with_nthreads(omp_get_num_procs()); #else - fftwf_plan_with_nthreads(2); + fftwf_plan_with_nthreads(1); #endif /* HAVE_GOMP */ } #endif /* HAVE_LIBFFTW3F_THREADS */ @@ -395,11 +345,16 @@ double TLFDynGssKT::operator()(double t, const vector &par) const { */ double tt(0.), sigsqtsq(0.); int i; + #ifdef HAVE_GOMP + int chunk = fNSteps/omp_get_num_procs(); + if (chunk < 10) + chunk = 10; + #endif if(fabs(par[0])<0.00135538817) { double mcplusnu(-(fC+par[2])); #ifdef HAVE_GOMP - #pragma omp parallel for default(shared) private(i,tt) schedule(dynamic) + #pragma omp parallel for default(shared) private(i,tt) schedule(dynamic, chunk) #endif for(i = 0; i < static_cast(fNSteps); ++i) { tt=static_cast(i)*fDt; @@ -414,7 +369,7 @@ double TLFDynGssKT::operator()(double t, const vector &par) const { fFFTtime[0] = fDt; #ifdef HAVE_GOMP - #pragma omp parallel for default(shared) private(i,tt) schedule(dynamic) + #pragma omp parallel for default(shared) private(i,tt) schedule(dynamic, chunk) #endif for(i = 1; i < static_cast(fNSteps); ++i) { tt=(static_cast(i)-0.5)*fDt; @@ -442,7 +397,10 @@ double TLFDynGssKT::operator()(double t, const vector &par) const { double denom(1.0), imagsq(0.0), oneMINrealnu(1.0); #ifdef HAVE_GOMP - #pragma omp parallel for default(shared) private(i, imagsq, oneMINrealnu, denom) schedule(dynamic) + chunk = (fNSteps/2+1)/omp_get_num_procs(); + if (chunk < 10) + chunk = 10; + #pragma omp parallel for default(shared) private(i, imagsq, oneMINrealnu, denom) schedule(dynamic, chunk) #endif for (i = 0; i < static_cast(fNSteps)/2+1; ++i) { imagsq=fFFTfreq[i][1]*fFFTfreq[i][1]; @@ -564,7 +522,7 @@ TLFDynExpKT::TLFDynExpKT() : fCalcNeeded(true), fFirstCall(true), fCounter(0), f #ifdef HAVE_GOMP fftwf_plan_with_nthreads(omp_get_num_procs()); #else - fftwf_plan_with_nthreads(2); + fftwf_plan_with_nthreads(1); #endif /* HAVE_GOMP */ } #endif /* HAVE_LIBFFTW3F_THREADS */ @@ -739,10 +697,15 @@ double TLFDynExpKT::operator()(double t, const vector &par) const { double tt(0.); int i; + #ifdef HAVE_GOMP + int chunk = fNSteps/omp_get_num_procs(); + if (chunk < 10) + chunk = 10; + #endif if(TMath::Abs(par[0])<0.00135538817){ #ifdef HAVE_GOMP - #pragma omp parallel for default(shared) private(i, tt) schedule(dynamic) + #pragma omp parallel for default(shared) private(i, tt) schedule(dynamic, chunk) #endif for(i = 0; i < static_cast(fNSteps); ++i) { tt=static_cast(i)*fDt; @@ -757,7 +720,7 @@ double TLFDynExpKT::operator()(double t, const vector &par) const { fFFTtime[0] = 1.0*fDt; #ifdef HAVE_GOMP - #pragma omp parallel for default(shared) private(i, tt) schedule(dynamic) + #pragma omp parallel for default(shared) private(i, tt) schedule(dynamic, chunk) #endif for(i = 1; i < static_cast(fNSteps); ++i) { tt=(double(i)-0.5)*fDt; @@ -785,7 +748,10 @@ double TLFDynExpKT::operator()(double t, const vector &par) const { double denom(1.0), imagsq(0.0), oneMINrealnu(1.0); #ifdef HAVE_GOMP - #pragma omp parallel for default(shared) private(i, imagsq, oneMINrealnu, denom) schedule(dynamic) + chunk = (fNSteps/2+1)/omp_get_num_procs(); + if (chunk < 10) + chunk = 10; + #pragma omp parallel for default(shared) private(i, imagsq, oneMINrealnu, denom) schedule(dynamic, chunk) #endif for (i = 0; i < static_cast(fNSteps)/2+1; ++i) { imagsq=fFFTfreq[i][1]*fFFTfreq[i][1]; @@ -815,9 +781,8 @@ double TLFDynExpKT::operator()(double t, const vector &par) const { /** *

Constructor: Initialize variables and allocate memory for the integrator. */ -TLFSGInterpolation::TLFSGInterpolation() : fCalcNeeded(true), fLastTime(0.) { +TLFSGInterpolation::TLFSGInterpolation() { fIntegral = new TIntSGInterpolation(); - fPar.resize(4, 0.); } //-------------------------------------------------------------------------- @@ -829,8 +794,6 @@ TLFSGInterpolation::TLFSGInterpolation() : fCalcNeeded(true), fLastTime(0.) { TLFSGInterpolation::~TLFSGInterpolation() { delete fIntegral; fIntegral = 0; - fPar.clear(); - fIntValues.clear(); } //-------------------------------------------------------------------------- @@ -856,35 +819,14 @@ double TLFSGInterpolation::operator()(double t, const vector &par) const 0.66666666666666666667*(1.0-par[1]*par[1]*t*t/TMath::Power(expo,(1.0-par[3])))*TMath::Exp(-TMath::Power(expo,par[3]))); } - for (unsigned int i(0); i<4; ++i) { - if( fPar[i]-par[i] ) { - fPar[i] = par[i]; - fCalcNeeded = true; - } - } - - if (fCalcNeeded) { - fIntValues.clear(); - fIntValues[0.] = 0.; - fLastTime = 0.; - } else { - if (t > fIntValues.rbegin()->first) { - fCalcNeeded = true; - fLastTime = fIntValues.rbegin()->first; - } - } - double omegaL(TMath::TwoPi()*par[0]); - if (fCalcNeeded) { - vector intpar(par); - intpar.push_back(t); - fIntegral->SetParameters(intpar); - fIntValues[t] = fIntValues[fLastTime] + fIntegral->IntegrateFunc(fLastTime,t); - intpar.clear(); - } + vector intpar(par); + intpar.push_back(t); + double intValue = fIntegral->IntegrateFunc(0.0, t, intpar); + intpar.clear(); return (TMath::Exp(-TMath::Power(par[2]*t,par[3])) + 2.0*par[1]*par[1]/(omegaL*omegaL) * \ ((TMath::Cos(omegaL*t)-TMath::Sin(omegaL*t)/(omegaL*t))*TMath::Exp(-TMath::Power(expo,par[3]))/TMath::Power(expo,(1.0-par[3])) + \ - fIntValues[t])); + intValue)); } diff --git a/src/external/libLFRelaxation/TLFRelaxation.h b/src/external/libLFRelaxation/TLFRelaxation.h index 10ceb9dd..bf2370c0 100644 --- a/src/external/libLFRelaxation/TLFRelaxation.h +++ b/src/external/libLFRelaxation/TLFRelaxation.h @@ -65,12 +65,8 @@ public: private: TIntSinGss *fIntSinGss; ///< integrator - mutable map fIntValues; ///< previously calculated integral values - mutable vector fPar; ///< parameters of the function [\htmlonly νL=Bγμ/2π (MHz), σ (μs-1)\endhtmlonly \latexonly $\nu_{\mathrm{L}}=B\gamma_{\mu}/2\pi~(\mathrm{MHz})$, $\sigma~(\mu\mathrm{s}^{-1})$ \endlatexonly] - mutable bool fCalcNeeded; ///< flag specifying if an integration has to be done - mutable double fLastTime; ///< time of the last calculated function value - ClassDef(TLFStatGssKT,2) + ClassDef(TLFStatGssKT,3) }; //----------------------------------------------------------------------------------------------------------------- @@ -96,12 +92,8 @@ public: private: TIntBesselJ0Exp *fIntBesselJ0Exp; ///< integrator - mutable map fIntValues; ///< previously calculated integral values - mutable vector fPar; ///< parameters of the function [\htmlonly νL=Bγμ/2π (MHz), a (μs-1)\endhtmlonly \latexonly $\nu_{\mathrm{L}}=B\gamma_{\mu}/2\pi~(\mathrm{MHz})$, $a~(\mu\mathrm{s}^{-1})$ \endlatexonly] - mutable bool fCalcNeeded; ///< flag specifying if an integration has to be done - mutable double fLastTime; ///< time of the last calculated function value - ClassDef(TLFStatExpKT,2) + ClassDef(TLFStatExpKT,3) }; //----------------------------------------------------------------------------------------------------------------- @@ -270,12 +262,8 @@ public: private: TIntSGInterpolation *fIntegral; ///< integrator - mutable map fIntValues; ///< previously calculated integral values - mutable vector fPar; ///< parameters of the function [\htmlonly νL=Bγμ/2π (MHz), a (μs-1), λ (μs-1), β (1) \endhtmlonly \latexonly $\nu_{\mathrm{L}}=B\gamma_{\mu}/2\pi~(\mathrm{MHz})$, $a~(\mu\mathrm{s}^{-1})$, $\lambda~(\mu\mathrm{s}^{-1})$, $\beta~(1)$ \endlatexonly] - mutable bool fCalcNeeded; ///< flag specifying if an integration has to be done - mutable double fLastTime; ///< time of the last calculated function value - ClassDef(TLFSGInterpolation,2) + ClassDef(TLFSGInterpolation,3) };