From 2c514a881cdd1d0b78466e1dfb14a8c2f2e286be Mon Sep 17 00:00:00 2001 From: "Bastian M. Wojek" Date: Sun, 10 Apr 2011 16:27:36 +0000 Subject: [PATCH] Tried to fix the ASCII export from musrview in the case of a Fourier-power-difference (some more tests needed); additionally minor changes to the BMWlibs --- ChangeLog | 1 + src/classes/PFunction.cpp | 1 + src/classes/PMusrCanvas.cpp | 39 ++-- src/external/BMWtools/TTrimSPDataHandler.cpp | 45 ++-- src/external/BMWtools/TTrimSPDataHandler.h | 2 + .../classes/TBulkTriVortexFieldCalc.cpp | 197 +++++++++++++++++- src/external/libFitPofB/classes/TPofBCalc.cpp | 26 ++- src/external/libFitPofB/classes/TPofTCalc.cpp | 3 + src/external/libFitPofB/classes/TVortex.cpp | 166 +++++++++++++++ .../include/TBulkTriVortexFieldCalc.h | 16 ++ src/external/libFitPofB/include/TPofBCalc.h | 1 + src/external/libFitPofB/include/TVortex.h | 33 +++ .../libFitPofB/include/TVortexLinkDef.h | 1 + 13 files changed, 492 insertions(+), 39 deletions(-) diff --git a/ChangeLog b/ChangeLog index 219955a2..69c1040b 100644 --- a/ChangeLog +++ b/ChangeLog @@ -6,6 +6,7 @@ changes since 0.9.0 =================================== +FIXED ASCII export from musrview in case of a Fourier-Power- or Fourier-Phase-difference plot FIXED bug in asymmetry fit with fixed background musrfit 0.9.0 - changes since 0.8.0 diff --git a/src/classes/PFunction.cpp b/src/classes/PFunction.cpp index d4c42ff2..629dfc67 100644 --- a/src/classes/PFunction.cpp +++ b/src/classes/PFunction.cpp @@ -490,6 +490,7 @@ Double_t PFunction::EvalNode(PFuncTreeNode &node) Double_t denominator = EvalNode(node.children[1]); if (denominator == 0.0) { cerr << endl << "**PANIC ERROR**: PFunction::EvalNode: division by 0.0"; + cerr << endl << "**PANIC ERROR**: PFunction::EvalNode: requested operation: " << EvalNode(node.children[0]) << "/" << EvalNode(node.children[1]); cerr << endl; assert(0); } diff --git a/src/classes/PMusrCanvas.cpp b/src/classes/PMusrCanvas.cpp index 762b650d..a26c9dca 100644 --- a/src/classes/PMusrCanvas.cpp +++ b/src/classes/PMusrCanvas.cpp @@ -1433,6 +1433,7 @@ void PMusrCanvas::InitMusrCanvas(const Char_t* title, Int_t wtopx, Int_t wtopy, fValid = true; fMainCanvas->cd(); + fMainCanvas->Show(); fMainCanvas->Connect("ProcessedEvent(Int_t,Int_t,Int_t,TObject*)", "PMusrCanvas", @@ -4437,10 +4438,10 @@ void PMusrCanvas::SaveDataAscii() break; case PV_FOURIER_PWR: // get current x-range - xminBin = fData[0].dataFourierPwr->GetXaxis()->GetFirst(); // first bin of the zoomed range - xmaxBin = fData[0].dataFourierPwr->GetXaxis()->GetLast(); // last bin of the zoomed range - xmin = fData[0].dataFourierPwr->GetXaxis()->GetBinCenter(xminBin); - xmax = fData[0].dataFourierPwr->GetXaxis()->GetBinCenter(xmaxBin); + xminBin = fData[0].diffFourierPwr->GetXaxis()->GetFirst(); // first bin of the zoomed range + xmaxBin = fData[0].diffFourierPwr->GetXaxis()->GetLast(); // last bin of the zoomed range + xmin = fData[0].diffFourierPwr->GetXaxis()->GetBinCenter(xminBin); + xmax = fData[0].diffFourierPwr->GetXaxis()->GetBinCenter(xmaxBin); // fill ascii dump data for (UInt_t i=0; iGetNbinsX(); j++) { + for (Int_t j=1; jGetNbinsX(); j++) { // get frequency - freq = fData[i].dataFourierPwr->GetBinCenter(j); + freq = fData[i].diffFourierPwr->GetBinCenter(j); // check if time is in the current range if ((freq >= xmin) && (freq <= xmax)) { dump.dataX.push_back(freq); - dump.data.push_back(fData[i].dataFourierPwr->GetBinContent(j)); + dump.data.push_back(fData[i].diffFourierPwr->GetBinContent(j)); } } @@ -4469,10 +4470,10 @@ void PMusrCanvas::SaveDataAscii() break; case PV_FOURIER_PHASE: // get current x-range - xminBin = fData[0].dataFourierPhase->GetXaxis()->GetFirst(); // first bin of the zoomed range - xmaxBin = fData[0].dataFourierPhase->GetXaxis()->GetLast(); // last bin of the zoomed range - xmin = fData[0].dataFourierPhase->GetXaxis()->GetBinCenter(xminBin); - xmax = fData[0].dataFourierPhase->GetXaxis()->GetBinCenter(xmaxBin); + xminBin = fData[0].diffFourierPhase->GetXaxis()->GetFirst(); // first bin of the zoomed range + xmaxBin = fData[0].diffFourierPhase->GetXaxis()->GetLast(); // last bin of the zoomed range + xmin = fData[0].diffFourierPhase->GetXaxis()->GetBinCenter(xminBin); + xmax = fData[0].diffFourierPhase->GetXaxis()->GetBinCenter(xmaxBin); // fill ascii dump data for (UInt_t i=0; iGetNbinsX(); j++) { + for (Int_t j=1; jGetNbinsX(); j++) { // get frequency - freq = fData[i].dataFourierPhase->GetBinCenter(j); + freq = fData[i].diffFourierPhase->GetBinCenter(j); // check if time is in the current range if ((freq >= xmin) && (freq <= xmax)) { dump.dataX.push_back(freq); - dump.data.push_back(fData[i].dataFourierPhase->GetBinContent(j)); + dump.data.push_back(fData[i].diffFourierPhase->GetBinContent(j)); } } @@ -5005,8 +5006,18 @@ void PMusrCanvas::SaveDataAscii() fout << "freq" << dumpVector.size()/2-1 << ", F_diffIm" << dumpVector.size()/2-1 << endl; break; case PV_FOURIER_PWR: + fout << "% "; + for (UInt_t i=0; i &energies, bool fEnergyIter = fEnergy.end(); } +// Method checking if an implantation profile is available for a given energy +// The behavior is the similar to the find-algorithm but more robust (tiny deviations in the energies are allowed). +// If the given energy is found the methods sets the internal energy iterator to the element of the energy vector. +// If it is not found the energy iterator will point to the end() of the energy vector. + +void TTrimSPData::FindEnergy(double e) const { + for(fEnergyIter = fEnergy.begin(); fEnergyIter != fEnergy.end(); ++fEnergyIter) { + if(fabs(*fEnergyIter - e) < 0.05) + return; + } + return; +} + void TTrimSPData::UseHighResolution(double e) { - fEnergyIter = find(fEnergy.begin(), fEnergy.end(), e); + FindEnergy(e); if(fEnergyIter != fEnergy.end()) { unsigned int i(fEnergyIter - fEnergy.begin()); @@ -149,7 +162,7 @@ void TTrimSPData::UseHighResolution(double e) { vector TTrimSPData::DataZ(double e) const { - fEnergyIter = find(fEnergy.begin(), fEnergy.end(), e); + FindEnergy(e); if(fEnergyIter != fEnergy.end()) { unsigned int i(fEnergyIter - fEnergy.begin()); @@ -167,7 +180,7 @@ vector TTrimSPData::DataZ(double e) const { vector TTrimSPData::DataNZ(double e) const { - fEnergyIter = find(fEnergy.begin(), fEnergy.end(), e); + FindEnergy(e); if(fEnergyIter != fEnergy.end()) { unsigned int i(fEnergyIter - fEnergy.begin()); @@ -185,7 +198,7 @@ vector TTrimSPData::DataNZ(double e) const { vector TTrimSPData::OrigDataNZ(double e) const { - fEnergyIter = find(fEnergy.begin(), fEnergy.end(), e); + FindEnergy(e); if(fEnergyIter != fEnergy.end()) { unsigned int i(fEnergyIter - fEnergy.begin()); @@ -199,7 +212,7 @@ vector TTrimSPData::OrigDataNZ(double e) const { double TTrimSPData::DataDZ(double e) const { - fEnergyIter = find(fEnergy.begin(), fEnergy.end(), e); + FindEnergy(e); if(fEnergyIter != fEnergy.end()) { unsigned int i(fEnergyIter - fEnergy.begin()); @@ -223,9 +236,9 @@ double TTrimSPData::LayerFraction(double e, unsigned int layno, const vector& interface, const } } - fEnergyIter = find(fEnergy.begin(), fEnergy.end(), e); + FindEnergy(e); // If all weights are equal to one, use the original n(z) vector for(unsigned int i(0); i z, nz; - fEnergyIter = find(fEnergy.begin(), fEnergy.end(), e); + FindEnergy(e); if(fEnergyIter != fEnergy.end()) { unsigned int i(fEnergyIter - fEnergy.begin()); z = fDataZ[i]; nz = fDataNZ[i]; } else { - cout << "TTrimSPData::GetNofZ: No implantation profile available for the specified energy... Quitting!" << endl; + cout << "TTrimSPData::GetNofZ: No implantation profile available for the specified energy " << e << " keV... Quitting!" << endl; exit(-1); } @@ -379,7 +392,7 @@ double TTrimSPData::GetNofZ(double zz, double e) const { void TTrimSPData::Normalize(double e) const { - fEnergyIter = find(fEnergy.begin(), fEnergy.end(), e); + FindEnergy(e); if(fEnergyIter != fEnergy.end()) { unsigned int i(fEnergyIter - fEnergy.begin()); @@ -404,7 +417,7 @@ void TTrimSPData::Normalize(double e) const { //--------------------- bool TTrimSPData::IsNormalized(double e) const { - fEnergyIter = find(fEnergy.begin(), fEnergy.end(), e); + FindEnergy(e); if(fEnergyIter != fEnergy.end()) { unsigned int i(fEnergyIter - fEnergy.begin()); @@ -420,7 +433,7 @@ bool TTrimSPData::IsNormalized(double e) const { //--------------------- double TTrimSPData::MeanRange(double e) const { - fEnergyIter = find(fEnergy.begin(), fEnergy.end(), e); + FindEnergy(e); if(fEnergyIter != fEnergy.end()) { unsigned int i(fEnergyIter - fEnergy.begin()); @@ -444,7 +457,7 @@ double TTrimSPData::MeanRange(double e) const { double TTrimSPData::PeakRange(double e) const { - fEnergyIter = find(fEnergy.begin(), fEnergy.end(), e); + FindEnergy(e); if(fEnergyIter != fEnergy.end()) { unsigned int i(fEnergyIter - fEnergy.begin()); @@ -476,7 +489,7 @@ void TTrimSPData::ConvolveGss(double w, double e) const { vector z, nz, gss; double nn; - fEnergyIter = find(fEnergy.begin(), fEnergy.end(), e); + FindEnergy(e); if(fEnergyIter != fEnergy.end()) { unsigned int i(fEnergyIter - fEnergy.begin()); diff --git a/src/external/BMWtools/TTrimSPDataHandler.h b/src/external/BMWtools/TTrimSPDataHandler.h index 091f5da3..5adecba3 100644 --- a/src/external/BMWtools/TTrimSPDataHandler.h +++ b/src/external/BMWtools/TTrimSPDataHandler.h @@ -70,6 +70,8 @@ public: double PeakRange(double) const; private: + void FindEnergy(double) const; + vector fEnergy; ///< vector holding all available muon energies vector fDZ; ///< vector holding the spatial resolution of the TRIM.SP output for all energies vector< vector > fDataZ; ///< discrete points in real space for which n(z) has been calculated for all energies diff --git a/src/external/libFitPofB/classes/TBulkTriVortexFieldCalc.cpp b/src/external/libFitPofB/classes/TBulkTriVortexFieldCalc.cpp index de0ee522..04f8d99e 100644 --- a/src/external/libFitPofB/classes/TBulkTriVortexFieldCalc.cpp +++ b/src/external/libFitPofB/classes/TBulkTriVortexFieldCalc.cpp @@ -804,6 +804,195 @@ void TBulkTriVortexAGLFieldCalc::CalculateGrid() const { +TBulkTriVortexAGLIIFieldCalc::TBulkTriVortexAGLIIFieldCalc(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) { +#ifdef HAVE_GOMP + fftw_plan_with_nthreads(omp_get_num_procs()); +#else + fftw_plan_with_nthreads(2); +#endif /* HAVE_GOMP */ + } +#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 TBulkTriVortexAGLIIFieldCalc::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 vortex-core radius have to have finite values in order to calculate B(x,y)!" << endl; + return; + } + + double field(fabs(fParam[0])), lambda(fabs(fParam[1])), xiV(fabs(fParam[2])); + double Hc2(getHc2(xiV)); // use the vortex-core radius for Hc2-calculation (which is wrong since xi_GL should be used instead: one would therefore need one more parameter) + + 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 < xiV/sqrt(2.0))) { + int m; + #ifdef HAVE_GOMP + #pragma omp parallel for default(shared) private(m) schedule(dynamic) + #endif + for (m = 0; m < NFFTsq; ++m) { + fFFTout[m] = field; + } + // Set the flag which shows that the calculation has been done + fGridExists = true; + return; + } + + double latConstTr(sqrt(2.0*fluxQuantum/(field*sqrt3))); + double fInf(1.0); + + double lambdasq_scaled(4.0/3.0*pow(lambda*PI/latConstTr,2.0)); + + // ... now fill in the Fourier components if everything was okay above + double Gsq, sqrtfInfSqPlusLsqGsq, ll; + int k, l, lNFFT_2; + + for (l = 0; l < NFFT_2; l += 2) { + lNFFT_2 = l*(NFFT_2 + 1); + ll = 3.0*static_cast(l*l); + for (k = 0; k < NFFT_2; k += 2) { + Gsq = static_cast(k*k) + ll; + sqrtfInfSqPlusLsqGsq = sqrt(fInf*fInf + lambdasq_scaled*Gsq); + fFFTin[lNFFT_2 + k][0] = ((!k && !l) ? 1.0 : \ + fInf*TMath::BesselK1(xiV/lambda*sqrtfInfSqPlusLsqGsq)/(sqrtfInfSqPlusLsqGsq*TMath::BesselK1(xiV/lambda*fInf))); + fFFTin[lNFFT_2 + k][1] = 0.0; + fFFTin[lNFFT_2 + k + 1][0] = 0.0; + fFFTin[lNFFT_2 + k + 1][1] = 0.0; + } + k = NFFT_2; + Gsq = static_cast(k*k) + ll; + sqrtfInfSqPlusLsqGsq = sqrt(fInf*fInf + lambdasq_scaled*Gsq); + fFFTin[lNFFT_2 + k][0] = fInf*TMath::BesselK1(xiV/lambda*sqrtfInfSqPlusLsqGsq)/(sqrtfInfSqPlusLsqGsq*TMath::BesselK1(xiV/lambda*fInf)); + fFFTin[lNFFT_2 + k][1] = 0.0; + } + + + for (l = NFFT_2; l < NFFT; l += 2) { + lNFFT_2 = l*(NFFT_2 + 1); + ll = 3.0*static_cast((NFFT-l)*(NFFT-l)); + for (k = 0; k < NFFT_2; k += 2) { + Gsq = static_cast(k*k) + ll; + sqrtfInfSqPlusLsqGsq = sqrt(fInf*fInf + lambdasq_scaled*Gsq); + fFFTin[lNFFT_2 + k][0] = fInf*TMath::BesselK1(xiV/lambda*sqrtfInfSqPlusLsqGsq)/(sqrtfInfSqPlusLsqGsq*TMath::BesselK1(xiV/lambda*fInf)); + fFFTin[lNFFT_2 + k][1] = 0.0; + fFFTin[lNFFT_2 + k + 1][0] = 0.0; + fFFTin[lNFFT_2 + k + 1][1] = 0.0; + } + k = NFFT_2; + Gsq = static_cast(k*k) + ll; + sqrtfInfSqPlusLsqGsq = sqrt(fInf*fInf + lambdasq_scaled*Gsq); + fFFTin[lNFFT_2 + k][0] = fInf*TMath::BesselK1(xiV/lambda*sqrtfInfSqPlusLsqGsq)/(sqrtfInfSqPlusLsqGsq*TMath::BesselK1(xiV/lambda*fInf)); + fFFTin[lNFFT_2 + k][1] = 0.0; + } + + // intermediate rows + + for (l = 1; l < NFFT_2; l += 2) { + lNFFT_2 = l*(NFFT_2 + 1); + ll = 3.0*static_cast(l*l); + for (k = 0; k < NFFT_2; k += 2) { + Gsq = static_cast((k + 1)*(k + 1)) + ll; + sqrtfInfSqPlusLsqGsq = sqrt(fInf*fInf + lambdasq_scaled*Gsq); + fFFTin[lNFFT_2 + k][0] = 0.0; + fFFTin[lNFFT_2 + k][1] = 0.0; + fFFTin[lNFFT_2 + k + 1][0] = fInf*TMath::BesselK1(xiV/lambda*sqrtfInfSqPlusLsqGsq)/(sqrtfInfSqPlusLsqGsq*TMath::BesselK1(xiV/lambda*fInf)); + fFFTin[lNFFT_2 + k + 1][1] = 0.0; + } + k = NFFT_2; + fFFTin[lNFFT_2 + k][0] = 0.0; + fFFTin[lNFFT_2 + k][1] = 0.0; + } + + for (l = NFFT_2 + 1; l < NFFT; l += 2) { + lNFFT_2 = l*(NFFT_2 + 1); + ll = 3.0*static_cast((NFFT-l)*(NFFT-l)); + for (k = 0; k < NFFT_2; k += 2) { + Gsq = static_cast((k+1)*(k+1)) + ll; + sqrtfInfSqPlusLsqGsq = sqrt(fInf*fInf + lambdasq_scaled*Gsq); + fFFTin[lNFFT_2 + k][0] = 0.0; + fFFTin[lNFFT_2 + k][1] = 0.0; + fFFTin[lNFFT_2 + k + 1][0] = fInf*TMath::BesselK1(xiV/lambda*sqrtfInfSqPlusLsqGsq)/(sqrtfInfSqPlusLsqGsq*TMath::BesselK1(xiV/lambda*fInf)); + fFFTin[lNFFT_2 + k + 1][1] = 0.0; + } + k = NFFT_2; + fFFTin[lNFFT_2 + k][0] = 0.0; + fFFTin[lNFFT_2 + k][1] = 0.0; + } + + // Do the Fourier transform to get B(x,y) + + fftw_execute(fFFTplan); + + // Multiply by the applied field + #ifdef HAVE_GOMP + #pragma omp parallel for default(shared) private(l) schedule(dynamic) + #endif + for (l = 0; l < NFFTsq; ++l) { + fFFTout[l] *= field; + } + + // Set the flag which shows that the calculation has been done + + fGridExists = true; + return; + +} + + + + TBulkTriVortexNGLFieldCalc::TBulkTriVortexNGLFieldCalc(const string& wisdom, const unsigned int steps) : fLatticeConstant(0.0), fKappa(0.0), fSumAk(0.0), fSumOmegaSq(0.0), fSumSum(0.0) { @@ -1728,11 +1917,11 @@ void TBulkTriVortexNGLFieldCalc::CalculateGrid() const { for (l = 0; l < NFFT; l++) { if (fFFTin[l][0]){ - if (((fabs(fFFTin[l][0]) > 1.0E-6) && (fabs(fCheckAkConvergence[l] - fFFTin[l][0])/fFFTin[l][0] > 1.0E-6)) || \ + if (((fabs(fFFTin[l][0]) > 1.0E-6) && (fabs(fCheckAkConvergence[l] - fFFTin[l][0])/fFFTin[l][0] > 1.0E-3)) || \ (fCheckAkConvergence[l]/fFFTin[l][0] < 0.0)) { - //cout << "old: " << fCheckAkConvergence[l] << ", new: " << fFFTin[l][0] << endl; + cout << "old: " << fCheckAkConvergence[l] << ", new: " << fFFTin[l][0] << endl; akConverged = false; - //cout << "index = " << l << endl; + cout << "index = " << l << endl; break; } } @@ -1808,7 +1997,7 @@ void TBulkTriVortexNGLFieldCalc::CalculateGrid() const { for (l = 0; l < NFFT; l++) { if (fBkMatrix[l][0]) { - if (((fabs(fBkMatrix[l][0]) > 1.0E-6) && (fabs(fCheckBkConvergence[l] - fBkMatrix[l][0])/fabs(fBkMatrix[l][0]) > 1.0E-6)) || \ + if (((fabs(fBkMatrix[l][0]) > 1.0E-6) && (fabs(fCheckBkConvergence[l] - fBkMatrix[l][0])/fabs(fBkMatrix[l][0]) > 1.0E-3)) || \ (fCheckBkConvergence[l]/fBkMatrix[l][0] < 0.0)) { // cout << "old: " << fCheckBkConvergence[l] << ", new: " << fBkMatrix[l][0] << endl; bkConverged = false; diff --git a/src/external/libFitPofB/classes/TPofBCalc.cpp b/src/external/libFitPofB/classes/TPofBCalc.cpp index 3f0fa70d..42b55019 100644 --- a/src/external/libFitPofB/classes/TPofBCalc.cpp +++ b/src/external/libFitPofB/classes/TPofBCalc.cpp @@ -153,6 +153,23 @@ void TPofBCalc::Normalize(unsigned int minFilledIndex = 0, unsigned int maxFille fPB[i] /= pBsum; } +// Do not actually calculate P(B) but take it from a B and a P(B) vector of the same size + +void TPofBCalc::SetPB(const vector &pb) const { + assert(fPBSize == pb.size()); + + int i; + #ifdef HAVE_GOMP + #pragma omp parallel for default(shared) private(i) schedule(dynamic) + #endif + for (i = 0; i < static_cast(fPBSize); ++i) { + fPB[i] = pb[i]; + } + + fPBExists = true; + return; +} + void TPofBCalc::Calculate(const string &type, const vector ¶) { if (type == "skg"){ // skewed Gaussian @@ -659,7 +676,7 @@ void TPofBCalc::ConvolveGss(double w) { TBin = 1.0/(gBar*static_cast(NFFT-1)*fDB); - FFTout = new fftw_complex[NFFT/2 + 1]; //(fftw_complex *)fftw_malloc(sizeof(fftw_complex) * (NFFT/2+1)); + FFTout = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * (NFFT/2+1));//new fftw_complex[NFFT/2 + 1]; // do the FFT to time domain @@ -674,10 +691,10 @@ void TPofBCalc::ConvolveGss(double w) { int i; #ifdef HAVE_GOMP - #pragma omp parallel for default(shared) private(GssInTimeDomain, i) schedule(dynamic) + #pragma omp parallel for default(shared) private(i, GssInTimeDomain) schedule(dynamic) #endif for (i = 0; i < static_cast(NFFT/2+1); ++i) { - GssInTimeDomain = exp(expo*static_cast(i*i)); + GssInTimeDomain = exp(expo*static_cast(i)*static_cast(i)); FFTout[i][0] *= GssInTimeDomain; FFTout[i][1] *= GssInTimeDomain; } @@ -692,8 +709,7 @@ void TPofBCalc::ConvolveGss(double w) { fftw_destroy_plan(FFTplanToTimeDomain); fftw_destroy_plan(FFTplanToFieldDomain); - delete[] FFTout; // fftw_free(FFTout); - FFTout = 0; + fftw_free(FFTout);//delete[] FFTout; FFTout = 0; // fftw_cleanup(); // normalize p(B) diff --git a/src/external/libFitPofB/classes/TPofTCalc.cpp b/src/external/libFitPofB/classes/TPofTCalc.cpp index 52b6c4c3..30a61b62 100644 --- a/src/external/libFitPofB/classes/TPofTCalc.cpp +++ b/src/external/libFitPofB/classes/TPofTCalc.cpp @@ -415,9 +415,12 @@ void TPofTCalc::FakeData(const string &rootOutputFileName, const vector histoDataPPC.clear(); fakeHisto = 0; + histoFolder->Clear(); delete histoFolder; histoFolder = 0; + decayAnaModule->Clear(); delete decayAnaModule; decayAnaModule = 0; + runInfoFolder->Clear(); delete runInfoFolder; runInfoFolder = 0; delete runHeader; runHeader = 0; diff --git a/src/external/libFitPofB/classes/TVortex.cpp b/src/external/libFitPofB/classes/TVortex.cpp index 1cf0e533..3a713ec9 100644 --- a/src/external/libFitPofB/classes/TVortex.cpp +++ b/src/external/libFitPofB/classes/TVortex.cpp @@ -41,6 +41,7 @@ ClassImp(TBulkTriVortexLondon) ClassImp(TBulkSqVortexLondon) ClassImp(TBulkTriVortexML) ClassImp(TBulkTriVortexAGL) +ClassImp(TBulkTriVortexAGLII) ClassImp(TBulkTriVortexNGL) ClassImp(TBulkAnisotropicTriVortexLondonGlobal) ClassImp(TBulkAnisotropicTriVortexLondon) @@ -117,6 +118,23 @@ TBulkTriVortexAGL::~TBulkTriVortexAGL() { fParForPofT.clear(); } +//------------------ +// Destructor of the TBulkTriVortexAGLII class -- cleaning up +//------------------ + +TBulkTriVortexAGLII::~TBulkTriVortexAGLII() { + delete fPofT; + fPofT = 0; + delete fPofB; + fPofT = 0; + delete fVortex; + fVortex = 0; + fPar.clear(); + fParForVortex.clear(); + fParForPofB.clear(); + fParForPofT.clear(); +} + //------------------ // Destructor of the TBulkTriVortexNGL class -- cleaning up //------------------ @@ -689,6 +707,153 @@ double TBulkTriVortexAGL::operator()(double t, const vector &par) const } +//------------------ +// Constructor of the TBulkTriVortexAGLII class +// creates (a pointer to) the TPofTCalc object (with the FFT plan) +//------------------ + +TBulkTriVortexAGLII::TBulkTriVortexAGLII() : fCalcNeeded(true), fFirstCall(true) { + + // read startup file + string startup_path_name("BMW_startup.xml"); + + TSAXParser *saxParser = new TSAXParser(); + BMWStartupHandler *startupHandler = new BMWStartupHandler(); + saxParser->ConnectToHandler("BMWStartupHandler", startupHandler); + int status (saxParser->ParseFile(startup_path_name.c_str())); + // check for parse errors + if (status) { // error + cerr << endl << "**ERROR** reading/parsing " << startup_path_name << " 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 + fParForPofB.push_back(0.0); // vortex-weighting || antiferromagnetic field + fParForPofB.push_back(0.0); // vortex-weighting: 0.0 homogeneous, 1.0 Gaussian, 2.0 Lorentzian || 3.0 antiferromagnetic vortex-cores + + fVortex = new TBulkTriVortexAGLIIFieldCalc(fWisdom, fGridSteps); + + fPofB = new TPofBCalc(fParForPofB); + + fPofT = new TPofTCalc(fPofB, fWisdom, fParForPofT); + + // clean up + if (saxParser) { + delete saxParser; + saxParser = 0; + } + if (startupHandler) { + delete startupHandler; + startupHandler = 0; + } +} + +//------------------ +// TBulkTriVortexAGLII-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 TBulkTriVortexAGLII (phase, av.field, lambda, core-radius, [not implemented: bkg weight]) +//------------------ + +double TBulkTriVortexAGLII::operator()(double t, const vector &par) const { + + assert(par.size() == 4 || par.size() == 5 || par.size() == 7 || par.size() == 8); + + 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); + +} //------------------ // Constructor of the TBulkTriVortexNGL class @@ -1454,6 +1619,7 @@ void TBulkAnisotropicTriVortexAGLGlobal::Calc(const vector &par) const { return; } + //------------------------------------------------------------------------------------ /** *

Constructor. diff --git a/src/external/libFitPofB/include/TBulkTriVortexFieldCalc.h b/src/external/libFitPofB/include/TBulkTriVortexFieldCalc.h index 8af4af28..63d488e5 100644 --- a/src/external/libFitPofB/include/TBulkTriVortexFieldCalc.h +++ b/src/external/libFitPofB/include/TBulkTriVortexFieldCalc.h @@ -166,6 +166,22 @@ public: }; +/** + *

Class for the calculation of the spatial field distribution B(x,y) within a 2D triangular vortex lattice + * using an analytical Ginzburg-Landau approximation for tiny applied fields + */ +class TBulkTriVortexAGLIIFieldCalc : public TBulkVortexFieldCalc { + +public: + + TBulkTriVortexAGLIIFieldCalc(const string&, const unsigned int steps = 256); + ~TBulkTriVortexAGLIIFieldCalc() {} + + void CalculateGrid() const; + bool IsTriangular() const {return true;} + +}; + /** *

Class for the calculation of the spatial field distribution B(x,y) within a 2D triangular vortex lattice * using the analytical Ginzburg-Landau approximation for an anisotropic superconductor with the field applied diff --git a/src/external/libFitPofB/include/TPofBCalc.h b/src/external/libFitPofB/include/TPofBCalc.h index 08476bd1..0c52cc18 100644 --- a/src/external/libFitPofB/include/TPofBCalc.h +++ b/src/external/libFitPofB/include/TPofBCalc.h @@ -64,6 +64,7 @@ public: double GetBmin() const {return fBmin;} double GetBmax() const {return fBmax;} unsigned int GetPBSize() const {return fPBSize;} + void SetPB(const vector&) const; void ConvolveGss(double); void AddBackground(double, double, double); double GetFirstMoment() const; diff --git a/src/external/libFitPofB/include/TVortex.h b/src/external/libFitPofB/include/TVortex.h index 482cf209..ba3d2768 100644 --- a/src/external/libFitPofB/include/TVortex.h +++ b/src/external/libFitPofB/include/TVortex.h @@ -166,6 +166,39 @@ private: ClassDef(TBulkTriVortexAGL,1) }; +/** + *

Class implementing the musrfit user function interface for calculating muon-spin depolarization functions + * originating from the spatial field distribution B(x,y) within a 2D triangular vortex lattice + * calculated using an analytical Ginzburg-Landau approximation for very small applied fields + */ +class TBulkTriVortexAGLII : public PUserFcnBase { + +public: + TBulkTriVortexAGLII(); + ~TBulkTriVortexAGLII(); + + virtual Bool_t NeedGlobalPart() const { return false; } + virtual void SetGlobalPart(vector &globalPart, UInt_t idx) { } + virtual Bool_t GlobalPartIsValid() const { return true; } + + double operator()(double, const vector&) const; + +private: + mutable vector fPar; ///< parameters of the model + TBulkTriVortexAGLIIFieldCalc *fVortex; ///< spatial field distribution B(x,y) within the vortex lattice + TPofBCalc *fPofB; ///< static field distribution P(B) + TPofTCalc *fPofT; ///< muon spin polarization p(t) + mutable bool fCalcNeeded; ///< tag needed to avoid unnecessary calculations if the core parameters were unchanged + mutable bool fFirstCall; ///< tag for checking if the function operator is called the first time + mutable vector fParForVortex; ///< parameters for the calculation of B(x,y) + mutable vector fParForPofB; ///< parameters for the calculation of P(B) + mutable vector fParForPofT; ///< parameters for the calculation of p(t) + string fWisdom; ///< file name of the FFTW wisdom file + unsigned int fGridSteps; ///< number of points in x- and y-direction for which B(x,y) is calculated + + ClassDef(TBulkTriVortexAGLII,1) +}; + /** *

Class implementing the musrfit user function interface for calculating muon spin depolarization functions * originating from the spatial field distribution B(x,y) within a 2D triangular vortex lattice diff --git a/src/external/libFitPofB/include/TVortexLinkDef.h b/src/external/libFitPofB/include/TVortexLinkDef.h index 59412579..6a27e6fe 100644 --- a/src/external/libFitPofB/include/TVortexLinkDef.h +++ b/src/external/libFitPofB/include/TVortexLinkDef.h @@ -39,6 +39,7 @@ #pragma link C++ class TBulkSqVortexLondon+; #pragma link C++ class TBulkTriVortexML+; #pragma link C++ class TBulkTriVortexAGL+; +#pragma link C++ class TBulkTriVortexAGLII+; #pragma link C++ class TBulkTriVortexNGL+; #pragma link C++ class TBulkAnisotropicTriVortexLondon+; #pragma link C++ class TBulkAnisotropicTriVortexLondonGlobal+;