diff --git a/src/classes/PFitter.cpp b/src/classes/PFitter.cpp index cd1b98a6..db7a6af2 100644 --- a/src/classes/PFitter.cpp +++ b/src/classes/PFitter.cpp @@ -1128,6 +1128,7 @@ Bool_t PFitter::ExecuteSave() else hcorr->Draw("COLZ"); ccorr->Write("ccorr", TObject::kOverwrite, sizeof(ccorr)); + hcorr->Write("hcorr", TObject::kOverwrite, sizeof(hcorr)); ff.Close(); // clean up if (ccorr) { diff --git a/src/external/BMWIntegrator/BMWIntegrator.h b/src/external/BMWIntegrator/BMWIntegrator.h index 5b6b250c..69afe8c2 100644 --- a/src/external/BMWIntegrator/BMWIntegrator.h +++ b/src/external/BMWIntegrator/BMWIntegrator.h @@ -61,10 +61,7 @@ class TIntegrator { }; inline TIntegrator::TIntegrator() : fFunc(0) { - ROOT::Math::GSLIntegrator *integrator = new ROOT::Math::GSLIntegrator(ROOT::Math::Integration::kADAPTIVE,ROOT::Math::Integration::kGAUSS31); - fIntegrator = integrator; - integrator = 0; - delete integrator; + fIntegrator = new ROOT::Math::GSLIntegrator(ROOT::Math::Integration::kADAPTIVE,ROOT::Math::Integration::kGAUSS31); } inline TIntegrator::~TIntegrator(){ @@ -104,10 +101,7 @@ class TMCIntegrator { }; inline TMCIntegrator::TMCIntegrator() : fFunc(0) { - ROOT::Math::GSLMCIntegrator *integrator = new ROOT::Math::GSLMCIntegrator(ROOT::Math::MCIntegration::kMISER, 1.E-6, 1.E-4, 500000); - fMCIntegrator = integrator; - integrator = 0; - delete integrator; + fMCIntegrator = new ROOT::Math::GSLMCIntegrator(ROOT::Math::MCIntegration::kMISER, 1.E-6, 1.E-4, 500000); } inline TMCIntegrator::~TMCIntegrator(){ diff --git a/src/external/TFitPofB-lib/classes/TBulkTriVortexFieldCalc.cpp b/src/external/TFitPofB-lib/classes/TBulkTriVortexFieldCalc.cpp index a8a93032..b00275f4 100644 --- a/src/external/TFitPofB-lib/classes/TBulkTriVortexFieldCalc.cpp +++ b/src/external/TFitPofB-lib/classes/TBulkTriVortexFieldCalc.cpp @@ -286,6 +286,130 @@ void TBulkTriVortexLondonFieldCalc::CalculateGrid() const { } +TBulkSqVortexLondonFieldCalc::TBulkSqVortexLondonFieldCalc(const string& wisdom, const unsigned int steps) { + fWisdom = wisdom; + if (steps % 2) { + fSteps = steps + 1; + } else { + fSteps = steps; + } + fParam.resize(3); + fGridExists = false; + +#ifdef HAVE_LIBFFTW3_THREADS + int init_threads(fftw_init_threads()); + if (init_threads) + fftw_plan_with_nthreads(2); +#endif /* HAVE_LIBFFTW3_THREADS */ + + fFFTin = new fftw_complex[(fSteps/2 + 1) * fSteps]; + fFFTout = new double[fSteps*fSteps]; + +// cout << "Check for the FFT plan..." << endl; + +// Load wisdom from file if it exists and should be used + + fUseWisdom = true; + int wisdomLoaded(0); + + FILE *wordsOfWisdomR; + wordsOfWisdomR = fopen(fWisdom.c_str(), "r"); + if (wordsOfWisdomR == NULL) { + fUseWisdom = false; + } else { + wisdomLoaded = fftw_import_wisdom_from_file(wordsOfWisdomR); + fclose(wordsOfWisdomR); + } + + if (!wisdomLoaded) { + fUseWisdom = false; + } + +// create the FFT plan + + if (fUseWisdom) + fFFTplan = fftw_plan_dft_c2r_2d(fSteps, fSteps, fFFTin, fFFTout, FFTW_EXHAUSTIVE); + else + fFFTplan = fftw_plan_dft_c2r_2d(fSteps, fSteps, fFFTin, fFFTout, FFTW_ESTIMATE); +} + +void TBulkSqVortexLondonFieldCalc::CalculateGrid() const { + // SetParameters - method has to be called from the user before the calculation!! + if (fParam.size() < 3) { + cout << endl << "The SetParameters-method has to be called before B(x,y) can be calculated!" << endl; + return; + } + if (!fParam[0] || !fParam[1] || !fParam[2]) { + cout << endl << "The field, penetration depth and coherence length have to have finite values in order to calculate B(x,y)!" << endl; + return; + } + + double field(fabs(fParam[0])), lambda(fabs(fParam[1])), xi(fabs(fParam[2])); + double Hc2(getHc2(xi)); + + double latConstSq(sqrt(fluxQuantum/field)); + double xisq_2_scaled(2.0*pow(xi*PI/latConstSq,2.0)), lambdasq_scaled(4.0*pow(lambda*PI/latConstSq,2.0)); + + const int NFFT(fSteps); + const int NFFT_2(fSteps/2); + const int NFFTsq(fSteps*fSteps); + + // fill the field Fourier components in the matrix + + // ... but first check that the field is not larger than Hc2 and that we are dealing with a type II SC + if ((field >= Hc2) || (lambda < xi/sqrt(2.0))) { + int m; + #pragma omp parallel for default(shared) private(m) schedule(dynamic) + for (m = 0; m < NFFTsq; m++) { + fFFTout[m] = field; + } + // Set the flag which shows that the calculation has been done + fGridExists = true; + return; + } + + // ... now fill in the Fourier components if everything was okay above + double Gsq, ll; + int k, l, lNFFT_2; + + for (l = 0; l < NFFT_2; ++l) { + lNFFT_2 = l*(NFFT_2 + 1); + ll = static_cast(l*l); + for (k = 0; k <= NFFT_2; ++k) { + Gsq = static_cast(k*k) + ll; + fFFTin[lNFFT_2 + k][0] = exp(-xisq_2_scaled*Gsq)/(1.0+lambdasq_scaled*Gsq); + fFFTin[lNFFT_2 + k][1] = 0.0; + } + } + + for (l = NFFT_2; l < NFFT; ++l) { + lNFFT_2 = l*(NFFT_2 + 1); + ll = static_cast((NFFT-l)*(NFFT-l)); + for (k = 0; k <= NFFT_2; ++k) { + Gsq = static_cast(k*k) + ll; + fFFTin[lNFFT_2 + k][0] = exp(-xisq_2_scaled*Gsq)/(1.0+lambdasq_scaled*Gsq); + fFFTin[lNFFT_2 + k][1] = 0.0; + } + } + + // Do the Fourier transform to get B(x,y) + + fftw_execute(fFFTplan); + + // Multiply by the applied field + #pragma omp parallel for default(shared) private(l) schedule(dynamic) + for (l = 0; l < NFFTsq; l++) { + fFFTout[l] *= field; + } + + // Set the flag which shows that the calculation has been done + + fGridExists = true; + return; + +} + + TBulkTriVortexMLFieldCalc::TBulkTriVortexMLFieldCalc(const string& wisdom, const unsigned int steps) { fWisdom = wisdom; diff --git a/src/external/TFitPofB-lib/classes/TFitPofBStartupHandler.cpp b/src/external/TFitPofB-lib/classes/TFitPofBStartupHandler.cpp index a606ad2a..a007f867 100644 --- a/src/external/TFitPofB-lib/classes/TFitPofBStartupHandler.cpp +++ b/src/external/TFitPofB-lib/classes/TFitPofBStartupHandler.cpp @@ -49,7 +49,7 @@ ClassImpQ(TFitPofBStartupHandler) /** *

*/ -TFitPofBStartupHandler::TFitPofBStartupHandler() : fDeltat(0.), fDeltaB(0.), fNSteps(0) +TFitPofBStartupHandler::TFitPofBStartupHandler() : fDeltat(0.), fDeltaB(0.), fNSteps(0), fGridSteps(0) { } diff --git a/src/external/TFitPofB-lib/classes/TLondon1D.cpp b/src/external/TFitPofB-lib/classes/TLondon1D.cpp index e6c84e9c..a86882e1 100644 --- a/src/external/TFitPofB-lib/classes/TLondon1D.cpp +++ b/src/external/TFitPofB-lib/classes/TLondon1D.cpp @@ -209,18 +209,11 @@ TLondon1DHS::TLondon1DHS() : fCalcNeeded(true), fFirstCall(true) { fParForPofB.push_back(0.005); // Bkg-width fParForPofB.push_back(0.0); // Bkg-weight - TTrimSPData *x = new TTrimSPData(rge_path, energy_vec); - fImpProfile = x; - x = 0; + fImpProfile = new TTrimSPData(rge_path, energy_vec); - TPofBCalc *y = new TPofBCalc(fParForPofB); - fPofB = y; - y = 0; - - TPofTCalc *z = new TPofTCalc(fPofB, fWisdom, fParForPofT); - fPofT = z; - z = 0; + fPofB = new TPofBCalc(fParForPofB); + fPofT = new TPofTCalc(fPofB, fWisdom, fParForPofT); // clean up if (saxParser) { @@ -370,17 +363,11 @@ TLondon1D1L::TLondon1D1L() : fCalcNeeded(true), fFirstCall(true), fCallCounter(0 fParForPofB.push_back(startupHandler->GetDeltaB()); fParForPofB.push_back(0.0); - TTrimSPData *x = new TTrimSPData(rge_path, energy_vec); - fImpProfile = x; - x = 0; + fImpProfile = new TTrimSPData(rge_path, energy_vec); - TPofBCalc *y = new TPofBCalc(fParForPofB); - fPofB = y; - y = 0; + fPofB = new TPofBCalc(fParForPofB); - TPofTCalc *z = new TPofTCalc(fPofB, fWisdom, fParForPofT); - fPofT = z; - z = 0; + fPofT = new TPofTCalc(fPofB, fWisdom, fParForPofT); // clean up if (saxParser) { @@ -524,18 +511,11 @@ TLondon1D2L::TLondon1D2L() : fCalcNeeded(true), fFirstCall(true), fLastTwoChange fParForPofB.push_back(startupHandler->GetDeltaB()); fParForPofB.push_back(0.0); - TTrimSPData *x = new TTrimSPData(rge_path, energy_vec); - fImpProfile = x; - x = 0; + fImpProfile = new TTrimSPData(rge_path, energy_vec); - TPofBCalc *y = new TPofBCalc(fParForPofB); - fPofB = y; - y = 0; - - TPofTCalc *z = new TPofTCalc(fPofB, fWisdom, fParForPofT); - fPofT = z; - z = 0; + fPofB = new TPofBCalc(fParForPofB); + fPofT = new TPofTCalc(fPofB, fWisdom, fParForPofT); // clean up if (saxParser) { @@ -684,17 +664,11 @@ TProximity1D1LHS::TProximity1D1LHS() : fCalcNeeded(true), fFirstCall(true) { fParForPofB.push_back(0.01); // Bkg-width fParForPofB.push_back(0.0); // Bkg-weight - TTrimSPData *x = new TTrimSPData(rge_path, energy_vec); - fImpProfile = x; - x = 0; + fImpProfile = new TTrimSPData(rge_path, energy_vec); - TPofBCalc *y = new TPofBCalc(fParForPofB); - fPofB = y; - y = 0; + fPofB = new TPofBCalc(fParForPofB); - TPofTCalc *z = new TPofTCalc(fPofB, fWisdom, fParForPofT); - fPofT = z; - z = 0; + fPofT = new TPofTCalc(fPofB, fWisdom, fParForPofT); // clean up if (saxParser) { @@ -851,17 +825,11 @@ TProximity1D1LHSGss::TProximity1D1LHSGss() : fCalcNeeded(true), fFirstCall(true) fParForPofB.push_back(0.0); // fParForPofB.push_back(0.0); - TTrimSPData *x = new TTrimSPData(rge_path, energy_vec); - fImpProfile = x; - x = 0; + fImpProfile = new TTrimSPData(rge_path, energy_vec); - TPofBCalc *y = new TPofBCalc(fParForPofB); - fPofB = y; - y = 0; + fPofB = new TPofBCalc(fParForPofB); - TPofTCalc *z = new TPofTCalc(fPofB, fWisdom, fParForPofT); - fPofT = z; - z = 0; + fPofT = new TPofTCalc(fPofB, fWisdom, fParForPofT); // clean up if (saxParser) { @@ -991,18 +959,11 @@ TLondon1D3L::TLondon1D3L() : fCalcNeeded(true), fFirstCall(true), fLastThreeChan fParForPofB.push_back(startupHandler->GetDeltaB()); fParForPofB.push_back(0.0); - TTrimSPData *x = new TTrimSPData(rge_path, energy_vec); - fImpProfile = x; - x = 0; + fImpProfile = new TTrimSPData(rge_path, energy_vec); - TPofBCalc *y = new TPofBCalc(fParForPofB); - fPofB = y; - y = 0; - - TPofTCalc *z = new TPofTCalc(fPofB, fWisdom, fParForPofT); - fPofT = z; - z = 0; + fPofB = new TPofBCalc(fParForPofB); + fPofT = new TPofTCalc(fPofB, fWisdom, fParForPofT); // clean up if (saxParser) { @@ -1162,17 +1123,11 @@ TLondon1D3LS::TLondon1D3LS() : fCalcNeeded(true), fFirstCall(true), fLastThreeCh fParForPofB.push_back(startupHandler->GetDeltaB()); fParForPofB.push_back(0.0); - TTrimSPData *x = new TTrimSPData(rge_path, energy_vec); - fImpProfile = x; - x = 0; + fImpProfile = new TTrimSPData(rge_path, energy_vec); - TPofBCalc *y = new TPofBCalc(fParForPofB); - fPofB = y; - y = 0; + fPofB = new TPofBCalc(fParForPofB); - TPofTCalc *z = new TPofTCalc(fPofB, fWisdom, fParForPofT); - fPofT = z; - z = 0; + fPofT = new TPofTCalc(fPofB, fWisdom, fParForPofT); // clean up if (saxParser) { @@ -1314,15 +1269,9 @@ double TLondon1D3LS::operator()(double t, const vector &par) const { // fParForPofB.push_back(startupHandler->GetDeltaB()); // fParForPofB.push_back(0.0); // -// TTrimSPData *x = new TTrimSPData(rge_path, energy_vec); -// fImpProfile = x; -// x = 0; -// delete x; +// fImpProfile = new TTrimSPData(rge_path, energy_vec); // -// TPofTCalc *y = new TPofTCalc(fWisdom, fParForPofT); -// fPofT = y; -// y = 0; -// delete y; +// fPofT = new TPofTCalc(fWisdom, fParForPofT); // // // clean up // if (saxParser) { @@ -1485,17 +1434,11 @@ TLondon1D3LSub::TLondon1D3LSub() : fCalcNeeded(true), fFirstCall(true), fWeights fParForPofB.push_back(startupHandler->GetDeltaB()); fParForPofB.push_back(0.0); - TTrimSPData *x = new TTrimSPData(rge_path, energy_vec); - fImpProfile = x; - x = 0; + fImpProfile = new TTrimSPData(rge_path, energy_vec); - TPofBCalc *y = new TPofBCalc(fParForPofB); - fPofB = y; - y = 0; + fPofB = new TPofBCalc(fParForPofB); - TPofTCalc *z = new TPofTCalc(fPofB, fWisdom, fParForPofT); - fPofT = z; - z = 0; + fPofT = new TPofTCalc(fPofB, fWisdom, fParForPofT); // clean up if (saxParser) { diff --git a/src/external/TFitPofB-lib/classes/TPofBCalc.cpp b/src/external/TFitPofB-lib/classes/TPofBCalc.cpp index 9162c815..9c3ccaed 100644 --- a/src/external/TFitPofB-lib/classes/TPofBCalc.cpp +++ b/src/external/TFitPofB-lib/classes/TPofBCalc.cpp @@ -388,12 +388,12 @@ void TPofBCalc::Calculate(const TBulkVortexFieldCalc *vortexLattice, const vecto int a3(static_cast(floor(fBmax/fDB))); int a4(static_cast(ceil(fBmax/fDB))); - unsigned int firstZerosEnd ((a1 < a2) ? a1 : ((a1 > 0) ? (a1 - 1) : 0)); + //unsigned int firstZerosEnd ((a1 < a2) ? a1 : ((a1 > 0) ? (a1 - 1) : 0)); unsigned int lastZerosStart ((a3 < a4) ? a4 : (a4 + 1)); unsigned int numberOfSteps(vortexLattice->GetNumberOfSteps()); unsigned int numberOfStepsSq(numberOfSteps*numberOfSteps); unsigned int numberOfSteps_2(numberOfSteps/2); - unsigned int numberOfStepsSq_2(numberOfStepsSq/2); + //unsigned int numberOfStepsSq_2(numberOfStepsSq/2); if (lastZerosStart >= fPBSize) lastZerosStart = fPBSize - 1; @@ -405,30 +405,76 @@ void TPofBCalc::Calculate(const TBulkVortexFieldCalc *vortexLattice, const vecto } double *vortexFields = vortexLattice->DataB(); - unsigned int fill_index, counter(0); + unsigned int fill_index; - for (unsigned int j(0); j < numberOfStepsSq_2; j++) { - fill_index = static_cast(ceil(fabs((vortexFields[j]/fDB)))); - if (fill_index >= fPBSize) - fPB[fPBSize - 1] += 1.0; - else - fPB[fill_index] += 1.0; - counter++; - if (counter == numberOfSteps_2) { // sum only over the first quadrant of B(x,y) - counter = 0; - j += numberOfSteps_2; + if (para.size() == 7 && para[6] == 1.0 && para[5] != 0.0 && vortexLattice->IsTriangular()) { + // weight distribution with Gaussian around vortex-cores + double Rsq1, Rsq2, Rsq3, Rsq4, Rsq5, Rsq6, sigmaSq(-0.5*para[5]*para[5]); + for (unsigned int j(0); j < numberOfSteps_2; ++j) { + for (unsigned int i(0); i < numberOfSteps_2; ++i) { + fill_index = static_cast(ceil(fabs((vortexFields[i + numberOfSteps*j]/fDB)))); + if (fill_index < fPBSize) { + Rsq1 = static_cast(3*i*i + j*j)/static_cast(numberOfStepsSq); + Rsq2 = static_cast(3*(numberOfSteps_2 - i)*(numberOfSteps_2 - i) \ + + (numberOfSteps_2 - j)*(numberOfSteps_2 - j))/static_cast(numberOfStepsSq); + Rsq3 = static_cast(3*(numberOfSteps - i)*(numberOfSteps - i) \ + + j*j)/static_cast(numberOfStepsSq); + Rsq4 = static_cast(3*(numberOfSteps_2 - i)*(numberOfSteps_2 - i) \ + + (numberOfSteps_2 + j)*(numberOfSteps_2 + j))/static_cast(numberOfStepsSq); + Rsq5 = static_cast(3*i*i \ + + (numberOfSteps - j)*(numberOfSteps - j))/static_cast(numberOfStepsSq); + Rsq6 = static_cast(3*(numberOfSteps_2 + i)*(numberOfSteps_2 + i) \ + + (numberOfSteps_2 - j)*(numberOfSteps_2 - j))/static_cast(numberOfStepsSq); + fPB[fill_index] += exp(sigmaSq*Rsq1) + exp(sigmaSq*Rsq2) + exp(sigmaSq*Rsq3) \ + + exp(sigmaSq*Rsq4) + exp(sigmaSq*Rsq5) + exp(sigmaSq*Rsq6); + } + } + } + } else if (para.size() == 7 && para[6] == 2.0 && para[5] != 0.0 && vortexLattice->IsTriangular()) { + // weight distribution with Lorentzian around vortex-cores + double Rsq1, Rsq2, Rsq3, Rsq4, Rsq5, Rsq6, sigmaSq(para[5]*para[5]); + for (unsigned int j(0); j < numberOfSteps_2; ++j) { + for (unsigned int i(0); i < numberOfSteps_2; ++i) { + fill_index = static_cast(ceil(fabs((vortexFields[i + numberOfSteps*j]/fDB)))); + if (fill_index < fPBSize) { + Rsq1 = static_cast(3*i*i + j*j)/static_cast(numberOfStepsSq); + Rsq2 = static_cast(3*(numberOfSteps_2 - i)*(numberOfSteps_2 - i) \ + + (numberOfSteps_2 - j)*(numberOfSteps_2 - j))/static_cast(numberOfStepsSq); + Rsq3 = static_cast(3*(numberOfSteps - i)*(numberOfSteps - i) \ + + j*j)/static_cast(numberOfStepsSq); + Rsq4 = static_cast(3*(numberOfSteps_2 - i)*(numberOfSteps_2 - i) \ + + (numberOfSteps_2 + j)*(numberOfSteps_2 + j))/static_cast(numberOfStepsSq); + Rsq5 = static_cast(3*i*i \ + + (numberOfSteps - j)*(numberOfSteps - j))/static_cast(numberOfStepsSq); + Rsq6 = static_cast(3*(numberOfSteps_2 + i)*(numberOfSteps_2 + i) \ + + (numberOfSteps_2 - j)*(numberOfSteps_2 - j))/static_cast(numberOfStepsSq); + fPB[fill_index] += 1.0/(1.0+sigmaSq*Rsq1) + 1.0/(1.0+sigmaSq*Rsq2) + 1.0/(1.0+sigmaSq*Rsq3) \ + + 1.0/(1.0+sigmaSq*Rsq4) + 1.0/(1.0+sigmaSq*Rsq5) + 1.0/(1.0+sigmaSq*Rsq6); + } + } + } + } else { + for (unsigned int j(0); j < numberOfSteps_2; ++j) { + for (unsigned int i(0); i < numberOfSteps_2; ++i) { + fill_index = static_cast(ceil(fabs((vortexFields[i + numberOfSteps*j]/fDB)))); + if (fill_index < fPBSize) { + fPB[fill_index] += 1.0; + } + } } } vortexFields = 0; - double normalizer(static_cast(numberOfStepsSq_2/2)*fDB); - + // normalize P(B) + double sum(0.0); + for (unsigned int i(0); i < fPBSize; ++i) + sum += fPB[i]; + sum *= fDB; int i; #pragma omp parallel for default(shared) private(i) schedule(dynamic) - for (i = firstZerosEnd; i <= static_cast(lastZerosStart); i++) { - fPB[i] /= normalizer; - } + for (i = 0; i < static_cast(fPBSize); ++i) + fPB[i] /= sum; // end pragma omp parallel if(para.size() == 5) diff --git a/src/external/TFitPofB-lib/classes/TVortex.cpp b/src/external/TFitPofB-lib/classes/TVortex.cpp index 57866511..c24d0b93 100644 --- a/src/external/TFitPofB-lib/classes/TVortex.cpp +++ b/src/external/TFitPofB-lib/classes/TVortex.cpp @@ -39,6 +39,7 @@ using namespace std; #include "TFitPofBStartupHandler.h" ClassImp(TBulkTriVortexLondon) +ClassImp(TBulkSqVortexLondon) ClassImp(TBulkTriVortexML) ClassImp(TBulkTriVortexAGL) ClassImp(TBulkTriVortexNGL) @@ -60,6 +61,23 @@ TBulkTriVortexLondon::~TBulkTriVortexLondon() { fParForPofT.clear(); } +//------------------ +// Destructor of the TBulkSqVortexLondon class -- cleaning up +//------------------ + +TBulkSqVortexLondon::~TBulkSqVortexLondon() { + delete fPofT; + fPofT = 0; + delete fPofB; + fPofT = 0; + delete fVortex; + fVortex = 0; + fPar.clear(); + fParForVortex.clear(); + fParForPofB.clear(); + fParForPofT.clear(); +} + //------------------ // Destructor of the TBulkTriVortexML class -- cleaning up //------------------ @@ -148,18 +166,69 @@ TBulkTriVortexLondon::TBulkTriVortexLondon() : fCalcNeeded(true), fFirstCall(tru fParForPofB.push_back(0.0); // Bkg-Field fParForPofB.push_back(0.005); // Bkg-width fParForPofB.push_back(0.0); // Bkg-weight + fParForPofB.push_back(0.0); // vortex-weighting + fParForPofB.push_back(0.0); // vortex-weighting: 0.0 homogeneous, 1.0 Gaussian, 2.0 Lorentzian - TBulkTriVortexLondonFieldCalc *x = new TBulkTriVortexLondonFieldCalc(fWisdom, fGridSteps); - fVortex = x; - x = 0; + fVortex = new TBulkTriVortexLondonFieldCalc(fWisdom, fGridSteps); - TPofBCalc *y = new TPofBCalc(fParForPofB); - fPofB = y; - y = 0; + fPofB = new TPofBCalc(fParForPofB); - TPofTCalc *z = new TPofTCalc(fPofB, fWisdom, fParForPofT); - fPofT = z; - z = 0; + fPofT = new TPofTCalc(fPofB, fWisdom, fParForPofT); + + // clean up + if (saxParser) { + delete saxParser; + saxParser = 0; + } + if (startupHandler) { + delete startupHandler; + startupHandler = 0; + } +} + +//------------------ +// Constructor of the TBulkSqVortexLondon class +// creates (a pointer to) the TPofTCalc object (with the FFT plan) +//------------------ + +TBulkSqVortexLondon::TBulkSqVortexLondon() : fCalcNeeded(true), fFirstCall(true) { + + // read startup file + string startup_path_name("TFitPofB_startup.xml"); + + TSAXParser *saxParser = new TSAXParser(); + TFitPofBStartupHandler *startupHandler = new TFitPofBStartupHandler(); + saxParser->ConnectToHandler("TFitPofBStartupHandler", startupHandler); + int status (saxParser->ParseFile(startup_path_name.c_str())); + // check for parse errors + if (status) { // error + cerr << endl << "**ERROR** reading/parsing TFitPofB_startup.xml failed." \ + << endl << "**ERROR** Please make sure that the file exists in the local directory and it is set up correctly!" \ + << endl; + assert(false); + } + + fGridSteps = startupHandler->GetGridSteps(); + fWisdom = startupHandler->GetWisdomFile(); + + fParForVortex.resize(3); // field, lambda, xi + + fParForPofT.push_back(0.0); + fParForPofT.push_back(startupHandler->GetDeltat()); + fParForPofT.push_back(startupHandler->GetDeltaB()); + + fParForPofB.push_back(startupHandler->GetDeltat()); + fParForPofB.push_back(startupHandler->GetDeltaB()); + + fParForPofB.push_back(0.0); // Bkg-Field + fParForPofB.push_back(0.005); // Bkg-width + fParForPofB.push_back(0.0); // Bkg-weight + + fVortex = new TBulkSqVortexLondonFieldCalc(fWisdom, fGridSteps); + + fPofB = new TPofBCalc(fParForPofB); + + fPofT = new TPofTCalc(fPofB, fWisdom, fParForPofT); // clean up if (saxParser) { @@ -180,6 +249,91 @@ TBulkTriVortexLondon::TBulkTriVortexLondon() : fCalcNeeded(true), fFirstCall(tru double TBulkTriVortexLondon::operator()(double t, const vector &par) const { + assert(par.size() == 4 || par.size() == 5 || par.size() == 7); // normal, +BkgWeight, +VortexWeighting + + if(t<0.0) + return cos(par[0]*0.017453293); + + // check if the function is called the first time and if yes, read in parameters + + if(fFirstCall){ + fPar = par; + + for (unsigned int i(0); i < 3; i++) { + fParForVortex[i] = fPar[i+1]; + } + fFirstCall = false; + } + + // check if any parameter has changed + + bool par_changed(false); + bool only_phase_changed(false); + + for (unsigned int i(0); iSetParameters(fParForVortex); + fVortex->CalculateGrid(); + fPofB->UnsetPBExists(); + fPofB->Calculate(fVortex, fParForPofB); + fPofT->DoFFT(); + + }/* else { + cout << "Only the phase parameter has changed, (re-)calculating P(t) now..." << endl; + }*/ + + fPofT->CalcPol(fParForPofT); + + fCalcNeeded = false; + } + + return fPofT->Eval(t); + +} + +//------------------ +// TBulkSqVortexLondon-Method that calls the procedures to create B(x,y), p(B) and P(t) +// It finally returns P(t) for a given t. +// Parameters: all the parameters for the function to be fitted through TBulkSqVortexLondon (phase, av.field, lambda, xi, [not implemented: bkg weight]) +//------------------ + +double TBulkSqVortexLondon::operator()(double t, const vector &par) const { + assert(par.size() == 4 || par.size() == 5); if(t<0.0) @@ -252,6 +406,7 @@ double TBulkTriVortexLondon::operator()(double t, const vector &par) con } + //------------------ // Constructor of the TBulkTriVortexML class // creates (a pointer to) the TPofTCalc object (with the FFT plan) @@ -290,17 +445,11 @@ TBulkTriVortexML::TBulkTriVortexML() : fCalcNeeded(true), fFirstCall(true) { fParForPofB.push_back(0.005); // Bkg-width fParForPofB.push_back(0.0); // Bkg-weight - TBulkTriVortexMLFieldCalc *x = new TBulkTriVortexMLFieldCalc(fWisdom, fGridSteps); - fVortex = x; - x = 0; + fVortex = new TBulkTriVortexMLFieldCalc(fWisdom, fGridSteps); - TPofBCalc *y = new TPofBCalc(fParForPofB); - fPofB = y; - y = 0; + fPofB = new TPofBCalc(fParForPofB); - TPofTCalc *z = new TPofTCalc(fPofB, fWisdom, fParForPofT); - fPofT = z; - z = 0; + fPofT = new TPofTCalc(fPofB, fWisdom, fParForPofT); // clean up if (saxParser) { @@ -432,17 +581,11 @@ TBulkTriVortexAGL::TBulkTriVortexAGL() : fCalcNeeded(true), fFirstCall(true) { fParForPofB.push_back(0.005); // Bkg-width fParForPofB.push_back(0.0); // Bkg-weight - TBulkTriVortexAGLFieldCalc *x = new TBulkTriVortexAGLFieldCalc(fWisdom, fGridSteps); - fVortex = x; - x = 0; + fVortex = new TBulkTriVortexAGLFieldCalc(fWisdom, fGridSteps); - TPofBCalc *y = new TPofBCalc(fParForPofB); - fPofB = y; - y = 0; + fPofB = new TPofBCalc(fParForPofB); - TPofTCalc *z = new TPofTCalc(fPofB, fWisdom, fParForPofT); - fPofT = z; - z = 0; + fPofT = new TPofTCalc(fPofB, fWisdom, fParForPofT); // clean up if (saxParser) { @@ -574,17 +717,11 @@ TBulkTriVortexNGL::TBulkTriVortexNGL() : fCalcNeeded(true), fFirstCall(true) { fParForPofB.push_back(0.005); // Bkg-width fParForPofB.push_back(0.0); // Bkg-weight - TBulkTriVortexNGLFieldCalc *x = new TBulkTriVortexNGLFieldCalc(fWisdom, fGridSteps); - fVortex = x; - x = 0; + fVortex = new TBulkTriVortexNGLFieldCalc(fWisdom, fGridSteps); - TPofBCalc *y = new TPofBCalc(fParForPofB); - fPofB = y; - y = 0; + fPofB = new TPofBCalc(fParForPofB); - TPofTCalc *z = new TPofTCalc(fPofB, fWisdom, fParForPofT); - fPofT = z; - z = 0; + fPofT = new TPofTCalc(fPofB, fWisdom, fParForPofT); // clean up if (saxParser) { diff --git a/src/external/TFitPofB-lib/include/TBulkTriVortexFieldCalc.h b/src/external/TFitPofB-lib/include/TBulkTriVortexFieldCalc.h index 08055544..77fda00a 100644 --- a/src/external/TFitPofB-lib/include/TBulkTriVortexFieldCalc.h +++ b/src/external/TFitPofB-lib/include/TBulkTriVortexFieldCalc.h @@ -58,6 +58,7 @@ public: virtual bool GridExists() const {return fGridExists;} virtual void UnsetGridExists() const {fGridExists = false;} virtual unsigned int GetNumberOfSteps() const {return fSteps;} + virtual bool IsTriangular() const = 0; protected: vector fParam; @@ -82,6 +83,23 @@ public: ~TBulkTriVortexLondonFieldCalc() {} void CalculateGrid() const; + bool IsTriangular() const {return true;} + +}; + +//-------------------- +// Class for square vortex lattice, London model with Gaussian Cutoff +//-------------------- + +class TBulkSqVortexLondonFieldCalc : public TBulkVortexFieldCalc { + +public: + + TBulkSqVortexLondonFieldCalc(const string&, const unsigned int steps = 256); + ~TBulkSqVortexLondonFieldCalc() {} + + void CalculateGrid() const; + bool IsTriangular() const {return false;} }; @@ -97,6 +115,7 @@ public: ~TBulkTriVortexMLFieldCalc() {} void CalculateGrid() const; + bool IsTriangular() const {return true;} }; @@ -112,6 +131,7 @@ public: ~TBulkTriVortexAGLFieldCalc() {} void CalculateGrid() const; + bool IsTriangular() const {return true;} }; @@ -134,6 +154,7 @@ public: double* GetBMatrix() const {return fFFTout;} fftw_complex* GetOmegaDiffMatrix() const {return fOmegaDiffMatrix;} fftw_complex* GetQMatrix() const {return fQMatrix;} + bool IsTriangular() const {return true;} private: diff --git a/src/external/TFitPofB-lib/include/TVortex.h b/src/external/TFitPofB-lib/include/TVortex.h index c1cdcac7..2e622aed 100644 --- a/src/external/TFitPofB-lib/include/TVortex.h +++ b/src/external/TFitPofB-lib/include/TVortex.h @@ -59,6 +59,30 @@ private: ClassDef(TBulkTriVortexLondon,1) }; +class TBulkSqVortexLondon : public PUserFcnBase { + +public: + TBulkSqVortexLondon(); + ~TBulkSqVortexLondon(); + + double operator()(double, const vector&) const; + +private: + mutable vector fPar; + TBulkSqVortexLondonFieldCalc *fVortex; + TPofBCalc *fPofB; + TPofTCalc *fPofT; + mutable bool fCalcNeeded; + mutable bool fFirstCall; + mutable vector fParForVortex; + mutable vector fParForPofB; + mutable vector fParForPofT; + string fWisdom; + unsigned int fGridSteps; + + ClassDef(TBulkSqVortexLondon,1) +}; + class TBulkTriVortexML : public PUserFcnBase { public: diff --git a/src/external/TFitPofB-lib/include/TVortexLinkDef.h b/src/external/TFitPofB-lib/include/TVortexLinkDef.h index 8cf7fcd6..b2efbe95 100644 --- a/src/external/TFitPofB-lib/include/TVortexLinkDef.h +++ b/src/external/TFitPofB-lib/include/TVortexLinkDef.h @@ -37,6 +37,7 @@ #pragma link off all functions; #pragma link C++ class TBulkTriVortexLondon+; +#pragma link C++ class TBulkSqVortexLondon+; #pragma link C++ class TBulkTriVortexML+; #pragma link C++ class TBulkTriVortexAGL+; #pragma link C++ class TBulkTriVortexNGL+; diff --git a/src/external/libLFRelaxation/TLFRelaxation.cpp b/src/external/libLFRelaxation/TLFRelaxation.cpp index f46c04ca..b0ffb866 100644 --- a/src/external/libLFRelaxation/TLFRelaxation.cpp +++ b/src/external/libLFRelaxation/TLFRelaxation.cpp @@ -53,10 +53,7 @@ ClassImp(TLFSGInterpolation) // LF Static Gaussian KT TLFStatGssKT::TLFStatGssKT() { - TIntSinGss *intSinGss = new TIntSinGss(); - fIntSinGss = intSinGss; - intSinGss = 0; - delete intSinGss; + fIntSinGss = new TIntSinGss(); } TLFStatGssKT::~TLFStatGssKT() { @@ -89,10 +86,7 @@ double TLFStatGssKT::operator()(double t, const vector &par) const { // LF Static Lorentzian KT TLFStatLorKT::TLFStatLorKT() { - TIntBesselJ0Exp *intBesselJ0Exp = new TIntBesselJ0Exp(); - fIntBesselJ0Exp = intBesselJ0Exp; - intBesselJ0Exp = 0; - delete intBesselJ0Exp; + fIntBesselJ0Exp = new TIntBesselJ0Exp(); } TLFStatLorKT::~TLFStatLorKT() {