diff --git a/src/external/TFitPofB-lib/classes/Makefile.TFitPofB b/src/external/TFitPofB-lib/classes/Makefile.TFitPofB deleted file mode 100644 index be2e475e..00000000 --- a/src/external/TFitPofB-lib/classes/Makefile.TFitPofB +++ /dev/null @@ -1,82 +0,0 @@ -#--------------------------------------------------- -# get compilation flags from root-config - -ROOTCFLAGS = $(shell $(ROOTSYS)/bin/root-config --cflags) - -#--------------------------------------------------- - -OS = LINUX -CXX = g++-4.4.0 -CXXFLAGS = -O3 -Wall -fopenmp -Wno-trigraphs -fPIC -MUSRFITINCLUDE = ../../../include -#MUSRFITINCLUDE = /home/l_wojek/rep/analysis/musrfit/src/include -LOCALINCLUDE = ../include -ROOTINCLUDE = $(ROOTSYS)/include -INCLUDES = -I$(LOCALINCLUDE) -I$(MUSRFITINCLUDE) -I$(ROOTINCLUDE) -LD = g++-4.4.0 -LDFLAGS = -O3 -SOFLAGS = -O3 -shared -fopenmp -lfftw3_threads -lfftw3 -lm -lpthread - -# the output from the root-config script: -CXXFLAGS += $(ROOTCFLAGS) -LDFLAGS += - -# some definitions: headers (used to generate *Dict* stuff), sources, objects,... -OBJS = -OBJS += TTrimSPDataHandler.o -OBJS += TBulkVortexFieldCalc.o -OBJS += TBofZCalc.o -OBJS += TPofBCalc.o -OBJS += TPofTCalc.o -OBJS += TFitPofBStartupHandler.o TFitPofBStartupHandlerDict.o -OBJS += TVortex.o TVortexDict.o -OBJS += TLondon1D.o TLondon1DDict.o -OBJS += TSkewedGss.o TSkewedGssDict.o - -SHLIB = libTFitPofB.so - -# make the shared lib: -# -all: $(SHLIB) - -$(SHLIB): $(OBJS) - @echo "---> Building shared library $(SHLIB) ..." - /bin/rm -f $(SHLIB) - $(LD) $(OBJS) $(SOFLAGS) -o $(SHLIB) - @echo "done" - -# clean up: remove all object file (and core files) -# semicolon needed to tell make there is no source -# for this target! -# -clean:; @rm -f $(OBJS) *Dict* core* - @echo "---> removing $(OBJS)" - -# -$(OBJS): %.o: %.cpp - $(CXX) $(INCLUDES) $(CXXFLAGS) -c $< - -# Generate the ROOT CINT dictionary - -TVortexDict.cpp: ../include/TVortex.h ../include/TVortexLinkDef.h - @echo "Generating dictionary $@..." - rootcint -f $@ -c -p -I$(MUSRFITINCLUDE) $^ - -TLondon1DDict.cpp: ../include/TLondon1D.h ../include/TLondon1DLinkDef.h - @echo "Generating dictionary $@..." - rootcint -f $@ -c -p -I$(MUSRFITINCLUDE) $^ - -TSkewedGssDict.cpp: ../include/TSkewedGss.h ../include/TSkewedGssLinkDef.h - @echo "Generating dictionary $@..." - rootcint -f $@ -c -p -I$(MUSRFITINCLUDE) $^ - -TFitPofBStartupHandlerDict.cpp: ../include/TFitPofBStartupHandler.h ../include/TFitPofBStartupHandlerLinkDef.h - @echo "Generating dictionary $@..." - rootcint -f $@ -c -p $^ - -install: all - @echo "Installing shared lib: libTFitPofB.so" -ifeq ($(OS),LINUX) - cp -pv $(SHLIB) $(ROOTSYS)/lib - cp -pv $(LOCALINCLUDE)/*.h $(ROOTSYS)/include -endif diff --git a/src/external/TFitPofB-lib/classes/Makefile.libTFitPofB b/src/external/TFitPofB-lib/classes/Makefile.libTFitPofB new file mode 100644 index 00000000..f359df5c --- /dev/null +++ b/src/external/TFitPofB-lib/classes/Makefile.libTFitPofB @@ -0,0 +1,187 @@ +#--------------------------------------------------- +# Makefile.libTFitPofB +# +# Author: Bastian M. Wojek +# e-mail: bastian.wojek@psi.ch +# +# $Id$ +# +#--------------------------------------------------- + +#--------------------------------------------------- +# get compilation and library flags from root-config + +ROOTCFLAGS = $(shell $(ROOTSYS)/bin/root-config --cflags) +ROOTLIBS = $(shell $(ROOTSYS)/bin/root-config --libs) +ROOTGLIBS = $(shell $(ROOTSYS)/bin/root-config --glibs) + +#--------------------------------------------------- +# depending on the architecture, choose the compiler, +# linker, and the flags to use +# + +ARCH = $(shell $(ROOTSYS)/bin/root-config --arch) + +ifeq ($(ARCH),linux) +OS = LINUX +endif +ifeq ($(ARCH),linuxx8664gcc) +OS = LINUX +endif +ifeq ($(ARCH),win32gcc) +OS = WIN32GCC +endif +ifeq ($(ARCH),macosx) +OS = DARWIN +endif + +# -- Linux +ifeq ($(OS),LINUX) +CXX = g++-4.4.0 +CXXFLAGS = -O3 -fopenmp -Wall -Wno-trigraphs -fPIC +PMUSRPATH = ../../../include +MNPATH = $(ROOTSYS)/include +LOCALPATH = ../include +FFTW3PATH = /usr/include +INCLUDES = -I$(PMUSRPATH) -I$(MNPATH) -I$(LOCALPATH) -I$(FFTW3PATH) +LD = g++-4.4.0 +LDFLAGS = -O +SOFLAGS = -shared -fopenmp +SHLIB = libTFitPofB.so +endif + +# -- Windows/Cygwin +ifeq ($(OS),WIN32GCC) +CXX = g++ +CXXFLAGS = -O3 -Wall -Wno-trigraphs -D_DLL +PMUSRPATH = ../../../include +MNPATH = $(ROOTSYS)/include +LOCALPATH = ../include +FFTW3PATH = /usr/include +INCLUDES = -I$(PMUSRPATH) -I$(MNPATH) -I$(LOCALPATH) -I$(FFTW3PATH) +LD = g++ +LDFLAGS = -O -Wl,--enable-auto-import -Wl,--enable-runtime-pseudo-reloc +SOFLAGS = -shared -Wl,--export-all-symbols +SHLIB = libTFitPofB.dll +endif + +# -- MacOSX/Darwin +ifeq ($(OS),DARWIN) +CXX = g++ +CXXFLAGS = -O3 -Wall -Wno-trigraphs -fPIC +PMUSRPATH = ../../../include +MNPATH = $(ROOTSYS)/include +FINKPATH = /sw/include +LOCALPATH = ../include +FFTW3PATH = /sw/include +INCLUDES = -I$(PMUSRPATH) -I$(MNPATH) -I$(LOCALPATH) -I$(FINKPATH) -I$(FFTW3PATH) +LD = g++ +LDFLAGS = -O -Xlinker -bind_at_load +SOFLAGS = -dynamiclib -flat_namespace -undefined suppress -Wl,-x +SHLIB = libTFitPofB.dylib +endif + +# the output from the root-config script: +CXXFLAGS += $(ROOTCFLAGS) +LDFLAGS += + +# the ROOT libraries (G = graphic) +LIBS = $(ROOTLIBS) +GLIBS = $(ROOTGLIBS) + +# PMusr lib +PMUSRLIB = -L$(ROOTSYS)/lib -lPMusr +# FFTW lib +FFTW3LIB = -lfftw3_threads -lfftw3 -lm -lpthread + +ifeq ($(OS),WIN32GCC) +# PMusr lib +PMUSRLIB = -L$(ROOTSYS)/lib -lPMusr -lMathMore +endif + +ifeq ($(OS),DARWIN) +# FFTW lib +FFTW3LIB = -L/sw/lib -lfftw3_threads -lfftw3 -lm -lpthread +# PMusr lib +PMUSRLIB = -L$(ROOTSYS)/lib -lPMusr -lMathMore +endif + +# some definitions: headers (used to generate *Dict* stuff), sources, objects,... +OBJS = +OBJS += TTrimSPDataHandler.o +OBJS += TBulkVortexFieldCalc.o +OBJS += TBofZCalc.o +OBJS += TPofBCalc.o +OBJS += TPofTCalc.o +OBJS += TFitPofBStartupHandler.o TFitPofBStartupHandlerDict.o +OBJS += TVortex.o TVortexDict.o +OBJS += TLondon1D.o TLondon1DDict.o +OBJS += TSkewedGss.o TSkewedGssDict.o + +INST_HEADER = +INST_HEADER += ../include/TBofZCalc.h +INST_HEADER += ../include/TBulkVortexFieldCalc.h +INST_HEADER += ../include/TFitPofBStartupHandler.h +INST_HEADER += ../include/TLondon1D.h +INST_HEADER += ../include/TPofBCalc.h +INST_HEADER += ../include/TPofTCalc.h +INST_HEADER += ../include/TSkewedGss.h +INST_HEADER += ../include/TTrimSPDataHandler.h +INST_HEADER += ../include/TVortex.h + +# make the shared libs: + +all: $(SHLIB) + +$(SHLIB): $(EXTOBJS) $(OBJS) + @echo "---> Building shared library $(SHLIB) ..." + /bin/rm -f $(SHLIB) + $(LD) $(SOFLAGS) $(LDFLAGS) -o $(SHLIB) $(OBJS) $(LIBS) $(PMUSRLIB) $(FFTW3LIB) + @echo "done" + +# clean up: remove all object file (and core files) +# semicolon needed to tell make there is no source +# for this target! +# +clean:; @rm -f $(OBJS) $(EXTOBJS) $(SHLIB) *Dict* core* *~ + @echo "---> removing $(OBJS) $(SHLIB)" + + +$(OBJS): %.o: %.cpp + $(CXX) $(INCLUDES) $(CXXFLAGS) -c $< + +# Generate the ROOT CINT dictionary + +TVortexDict.cpp: ../include/TVortex.h ../include/TVortexLinkDef.h + @echo "Generating dictionary $@..." + rootcint -f $@ -c -p -I$(PMUSRPATH) $^ + +TLondon1DDict.cpp: ../include/TLondon1D.h ../include/TLondon1DLinkDef.h + @echo "Generating dictionary $@..." + rootcint -f $@ -c -p -I$(PMUSRPATH) $^ + +TSkewedGssDict.cpp: ../include/TSkewedGss.h ../include/TSkewedGssLinkDef.h + @echo "Generating dictionary $@..." + rootcint -f $@ -c -p -I$(PMUSRPATH) $^ + +TFitPofBStartupHandlerDict.cpp: ../include/TFitPofBStartupHandler.h ../include/TFitPofBStartupHandlerLinkDef.h + @echo "Generating dictionary $@..." + rootcint -f $@ -c -p $^ + +install: all + @echo "Installing shared lib: $(SHLIB)" +ifeq ($(OS),LINUX) + cp -pv $(SHLIB) $(ROOTSYS)/lib + cp -pv $(INST_HEADER) $(ROOTSYS)/include +endif +ifeq ($(OS),WIN32GCC) + cp -pv $(SHLIB) $(ROOTSYS)/bin + ln -sf $(ROOTSYS)/bin/$(SHLIB) $(ROOTSYS)/lib/$(SHLIB) + cp -pv $(INST_HEADER) $(ROOTSYS)/include +endif +ifeq ($(OS),DARWIN) + cp -pv $(SHLIB) $(ROOTSYS)/lib + cp -pv $(INST_HEADER) $(ROOTSYS)/include +endif + +cleaninstall: clean install diff --git a/src/external/TFitPofB-lib/classes/TBofZCalc.cpp b/src/external/TFitPofB-lib/classes/TBofZCalc.cpp index 7c23855a..1218d1f2 100644 --- a/src/external/TFitPofB-lib/classes/TBofZCalc.cpp +++ b/src/external/TFitPofB-lib/classes/TBofZCalc.cpp @@ -9,6 +9,26 @@ ***************************************************************************/ +/*************************************************************************** + * Copyright (C) 2009 by Bastian M. Wojek * + * * + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program; if not, write to the * + * Free Software Foundation, Inc., * + * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. * + ***************************************************************************/ + #include "TBofZCalc.h" #include #include diff --git a/src/external/TFitPofB-lib/classes/TBulkVortexFieldCalc.cpp b/src/external/TFitPofB-lib/classes/TBulkVortexFieldCalc.cpp index 08613310..b25cc743 100644 --- a/src/external/TFitPofB-lib/classes/TBulkVortexFieldCalc.cpp +++ b/src/external/TFitPofB-lib/classes/TBulkVortexFieldCalc.cpp @@ -9,16 +9,40 @@ ***************************************************************************/ +/*************************************************************************** + * Copyright (C) 2009 by Bastian M. Wojek * + * * + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program; if not, write to the * + * Free Software Foundation, Inc., * + * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. * + ***************************************************************************/ + #include "TBulkVortexFieldCalc.h" +#include #include #include +#include +using namespace std; #define PI 3.14159265358979323846 #define TWOPI 6.28318530717958647692 const double fluxQuantum(2.067833667e7); // 10e14 times CGS units %% in CGS units should be 10^-7 // in this case this is Gauss per square nm +const double sqrt3(sqrt(3.0)); double getXi(const double hc2) { // get xi given Hc2 in Gauss if (hc2) @@ -34,68 +58,295 @@ double getHc2(const double xi) { // get Hc2 given xi in nm return 0.0; } -TBulkTriVortexLondonFieldCalc::TBulkTriVortexLondonFieldCalc(const vector& param, const unsigned int steps, const unsigned int Ncomp) - : fNFourierComp(Ncomp) { - fSteps = steps; - fParam = param; - fB.clear(); -} - -double TBulkTriVortexLondonFieldCalc::GetBatPoint(double relX, double relY) const { - - double latConstTr(sqrt(fluxQuantum/fParam[2])*sqrt(sqrt(4.0/3.0))); - double Hc2(getHc2(fParam[1])); - - double xCoord(latConstTr*relX); // in nanometers - double yCoord(latConstTr*relY); - - if ((fParam[2] < Hc2) && (fParam[0] > fParam[1]/sqrt(2.0))) { - double KLatTr(4.0*PI/(latConstTr*sqrt(3.0))); - double fourierSum(0.0); - for(int mm = -fNFourierComp; mm <= static_cast(fNFourierComp); mm++) { - for(int nn = -fNFourierComp; nn <= static_cast(fNFourierComp); nn++) { - fourierSum += cos(KLatTr*(xCoord*mm*(0.5*sqrt(3.0)) + yCoord*mm*0.5 + yCoord*nn))*exp(-(0.5*fParam[1]*fParam[1]*KLatTr*KLatTr)* \ - (0.75*mm*mm + (nn + 0.5*mm)*(nn + 0.5*mm)))/(1.0+(fParam[0]*KLatTr*fParam[0]*KLatTr)*(0.75*mm*mm + (nn + 0.5*mm)*(nn + 0.5*mm))); - } +TBulkVortexFieldCalc::~TBulkVortexFieldCalc() { + // if a wisdom file is used export the wisdom so it has not to be checked for the FFT-plan next time + if (fUseWisdom) { + FILE *wordsOfWisdomW; + wordsOfWisdomW = fopen(fWisdom.c_str(), "w"); + if (wordsOfWisdomW == NULL) { + cout << "TBulkVortexFieldCalc::~TBulkVortexFieldCalc(): Could not open file ... No wisdom is exported..." << endl; + } else { + fftw_export_wisdom_to_file(wordsOfWisdomW); + fclose(wordsOfWisdomW); } - return fParam[2]*fourierSum; } - else - return fParam[2]; + // clean up + + fftw_destroy_plan(fFFTplan); + delete[] fFFTin; fFFTin = 0; + delete[] fFFTout; fFFTout = 0; + //fftw_cleanup(); + //fftw_cleanup_threads(); } +TBulkTriVortexLondonFieldCalc::TBulkTriVortexLondonFieldCalc(const string& wisdom, const unsigned int steps) { + fWisdom = wisdom; + fSteps = steps; + fParam.resize(3); + fGridExists = false; + if (fSteps%2) + fSteps++; + + int init_threads(fftw_init_threads()); + if (init_threads) + fftw_plan_with_nthreads(2); + + 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); +} + +// double TBulkTriVortexLondonFieldCalc::GetBatPoint(double relX, double relY) const { +// +// double field(fParam[0]), lambda(fParam[1]), xi(fParam[2]); +// double xisq_2(0.5*xi*xi), lambdasq(lambda*lambda); +// +// double latConstTr(sqrt(fluxQuantum/field)*sqrt(sqrt(4.0/3.0))); +// double Hc2(getHc2(xi)); +// +// double xCoord(latConstTr*relX); // in nanometers +// double yCoord(latConstTr*relY); +// +// if ((field < Hc2) && (lambda > xi/sqrt(2.0))) { +// double KLatTr(4.0*PI/(latConstTr*sqrt3)); +// double fourierSum(1.0), fourierAdd(0.0), expo; +// +// // k = 0, l = 0 gives 1.0, already in fourierSum +// +// // l = 0, only integer k +// for (double k (1.0); k < static_cast(fNFourierComp); k+=1.0){ +// expo = 3.0*pow(KLatTr*k, 2.0); +// fourierAdd += exp(-xisq_2*expo)/(lambdasq*expo + 1.0)*cos(sqrt3*KLatTr*k*xCoord); +// } +// +// fourierSum += 2.0*fourierAdd; +// fourierAdd = 0.0; +// +// // k = 0, only integer l +// for (double l (1.0); l < static_cast(fNFourierComp); l+=1.0){ +// expo = pow(KLatTr*l, 2.0); +// fourierAdd += exp(-xisq_2*expo)/(lambdasq*expo + 1.0)*cos(KLatTr*l*yCoord); +// } +// +// fourierSum += 2.0*fourierAdd; +// fourierAdd = 0.0; +// +// // k != 0, l != 0 both integers +// for (double k (1.0); k < static_cast(fNFourierComp); k+=1.0){ +// for (double l (1.0); l < static_cast(fNFourierComp); l+=1.0){ +// expo = KLatTr*KLatTr*(3.0*k*k + l*l); +// fourierAdd += exp(-xisq_2*expo)/(lambdasq*expo + 1.0)*cos(sqrt3*KLatTr*k*xCoord)*cos(KLatTr*l*yCoord); +// } +// } +// +// fourierSum += 4.0*fourierAdd; +// fourierAdd = 0.0; +// +// // k != 0, l != 0 both half-integral numbers +// for (double k (0.5); k < static_cast(fNFourierComp); k+=1.0){ +// for (double l (0.5); l < static_cast(fNFourierComp); l+=1.0){ +// expo = KLatTr*KLatTr*(3.0*k*k + l*l); +// fourierAdd += exp(-xisq_2*expo)/(lambdasq*expo + 1.0)*cos(sqrt3*KLatTr*k*xCoord)*cos(KLatTr*l*yCoord); +// } +// } +// +// fourierSum += 4.0*fourierAdd; +// +// // for(int mm = -fNFourierComp; mm <= static_cast(fNFourierComp); mm++) { +// // for(int nn = -fNFourierComp; nn <= static_cast(fNFourierComp); nn++) { +// // fourierSum += cos(KLatTr*(xCoord*mm*(0.5*sqrt(3.0)) + yCoord*mm*0.5 + yCoord*nn))*exp(-(0.5*fParam[1]*fParam[1]*KLatTr*KLatTr)* +// // (0.75*mm*mm + (nn + 0.5*mm)*(nn + 0.5*mm)))/(1.0+(fParam[0]*KLatTr*fParam[0]*KLatTr)*(0.75*mm*mm + (nn + 0.5*mm)*(nn + 0.5*mm))); +// // } +// // } +// // cout << " " << fourierSum << ", "; +// return field*fourierSum; +// } +// else +// return field; +// +// } + void TBulkTriVortexLondonFieldCalc::CalculateGrid() const { + // SetParameters - method has to be called from the user before the calculation!! + double field(abs(fParam[0])), lambda(abs(fParam[1])), xi(abs(fParam[2])); + double Hc2(getHc2(xi)); - double bb(cos(PI/6.0)*0.5); - double yStep(bb/static_cast(fSteps)); - double xStep(0.5/static_cast(fSteps)); + double latConstTr(sqrt(fluxQuantum/field*sqrt(4.0/3.0))); + double xisq_2_scaled(2.0/3.0*pow(xi*PI/latConstTr,2.0)), lambdasq_scaled(4.0/3.0*pow(lambda*PI/latConstTr,2.0)); - fB.resize(fSteps); - for (unsigned int i(0); i < fSteps; i++) - fB[i].resize(fSteps); + int NFFT(fSteps); + int NFFT_2(fSteps/2); + int NFFTsq(fSteps*fSteps); - for(unsigned int i(0); i < fSteps; i++) { - int j; -#pragma omp parallel for default(shared) private(j) schedule(dynamic) - for(j=0; j < static_cast(fSteps); j++) { - fB[i][j] = GetBatPoint(static_cast(i)*xStep, static_cast(j)*yStep); + // fill the field Fourier components in the matrix + + int m; + + // ... 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))) { +#pragma omp parallel for default(shared) private(m) schedule(dynamic) + for (m = 0; m < NFFTsq; m++) { + fFFTout[m] = field; } -// end pragma omp parallel + // 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; + int k, l, lNFFT_2; + +// omp causes problems with the fftw_complex*... comment it out for the moment +//#pragma omp parallel default(shared) private(l,lNFFT_2,k,Gsq) + { +// #pragma omp sections nowait + { +// #pragma omp section +// #pragma omp parallel for default(shared) private(l,lNFFT_2,k) schedule(dynamic) + for (l = 0; l < NFFT_2; l += 2) { + lNFFT_2 = l*(NFFT_2 + 1); + for (k = 0; k < NFFT_2; k += 2) { + Gsq = 3.0*static_cast(k*k) + static_cast(l*l); + fFFTin[lNFFT_2 + k][0] = exp(-xisq_2_scaled*Gsq)/(1.0+lambdasq_scaled*Gsq); + 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 = 3.0*static_cast(k*k) + static_cast(l*l); + fFFTin[lNFFT_2 + k][0] = exp(-xisq_2_scaled*Gsq)/(1.0+lambdasq_scaled*Gsq); + fFFTin[lNFFT_2 + k][1] = 0.0; + } + +// #pragma omp section +// #pragma omp parallel for default(shared) private(l,lNFFT_2,k) schedule(dynamic) + for (l = NFFT_2; l < NFFT; l += 2) { + lNFFT_2 = l*(NFFT_2 + 1); + for (k = 0; k < NFFT_2; k += 2) { + Gsq = 3.0*static_cast(k*k) + static_cast((NFFT-l)*(NFFT-l)); + fFFTin[lNFFT_2 + k][0] = exp(-xisq_2_scaled*Gsq)/(1.0+lambdasq_scaled*Gsq); + 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 = 3.0*static_cast(k*k) + static_cast((NFFT-l)*(NFFT-l)); + fFFTin[lNFFT_2 + k][0] = exp(-xisq_2_scaled*Gsq)/(1.0+lambdasq_scaled*Gsq); + fFFTin[lNFFT_2 + k][1] = 0.0; + } + + // intermediate rows +// #pragma omp section +// #pragma omp parallel for default(shared) private(l,lNFFT_2,k) schedule(dynamic) + for (l = 1; l < NFFT_2; l += 2) { + lNFFT_2 = l*(NFFT_2 + 1); + for (k = 0; k < NFFT_2; k += 2) { + Gsq = 3.0*static_cast((k + 1)*(k + 1)) + static_cast(l*l); + fFFTin[lNFFT_2 + k][0] = 0.0; + fFFTin[lNFFT_2 + k][1] = 0.0; + fFFTin[lNFFT_2 + k + 1][0] = exp(-xisq_2_scaled*Gsq)/(1.0+lambdasq_scaled*Gsq); + 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; + } + +// #pragma omp section +// #pragma omp parallel for default(shared) private(l,lNFFT_2,k) schedule(dynamic) + for (l = NFFT_2 + 1; l < NFFT; l += 2) { + lNFFT_2 = l*(NFFT_2 + 1); + for (k = 0; k < NFFT_2; k += 2) { + Gsq = 3.0*static_cast((k+1)*(k+1)) + static_cast((NFFT-l)*(NFFT-l)); + fFFTin[lNFFT_2 + k][0] = 0.0; + fFFTin[lNFFT_2 + k][1] = 0.0; + fFFTin[lNFFT_2 + k + 1][0] = exp(-xisq_2_scaled*Gsq)/(1.0+lambdasq_scaled*Gsq); + 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; + } + } /* end of sections */ + + } /* end of parallel section */ + + // 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; } double TBulkTriVortexLondonFieldCalc::GetBmin() const { - return GetBatPoint(0.5, 0.5*tan(PI/6.0)); + if (fGridExists) { + double min(fFFTout[0]); + unsigned int minindex(0), counter(0); + for (unsigned int j(0); j < fSteps * fSteps / 2; j++) { + if (fFFTout[j] <= 0.0) { + return 0.0; + } + if (fFFTout[j] < min) { + min = fFFTout[j]; + minindex = j; + } + counter++; + if (counter == fSteps/2) { // check only the first quadrant of B(x,y) + counter = 0; + j += fSteps/2; + } + } + return fFFTout[minindex]; + } else { + CalculateGrid(); + return GetBmin(); + } } double TBulkTriVortexLondonFieldCalc::GetBmax() const { - if (!fB.empty()) - return fB[0][0]; - else - return GetBatPoint(0.0, 0.0); + if (fGridExists) + return fFFTout[0]; + else { + CalculateGrid(); + return GetBmax(); + } } diff --git a/src/external/TFitPofB-lib/classes/TFitPofBStartupHandler.cpp b/src/external/TFitPofB-lib/classes/TFitPofBStartupHandler.cpp index 454d66ae..8f9fef62 100644 --- a/src/external/TFitPofB-lib/classes/TFitPofBStartupHandler.cpp +++ b/src/external/TFitPofB-lib/classes/TFitPofBStartupHandler.cpp @@ -1,6 +1,6 @@ /*************************************************************************** - TFitPofBStartupHandler.h + TFitPofBStartupHandler.cpp Author: Bastian M. Wojek e-mail: bastian.wojek@psi.ch @@ -14,7 +14,7 @@ ***************************************************************************/ /*************************************************************************** - * Copyright (C) 2007 by Andreas Suter, Bastian M. Wojek * + * Copyright (C) 2007-2008 by Andreas Suter, Bastian M. Wojek * * * * * * This program is free software; you can redistribute it and/or modify * diff --git a/src/external/TFitPofB-lib/classes/TLondon1D.cpp b/src/external/TFitPofB-lib/classes/TLondon1D.cpp index 15cf6fc8..d97133c0 100644 --- a/src/external/TFitPofB-lib/classes/TLondon1D.cpp +++ b/src/external/TFitPofB-lib/classes/TLondon1D.cpp @@ -9,6 +9,26 @@ ***************************************************************************/ +/*************************************************************************** + * Copyright (C) 2009 by Bastian M. Wojek * + * * + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program; if not, write to the * + * Free Software Foundation, Inc., * + * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. * + ***************************************************************************/ + #include "TLondon1D.h" #include #include @@ -43,6 +63,8 @@ TLondon1DHS::~TLondon1DHS() { fParForPofT.clear(); delete fImpProfile; fImpProfile = 0; + delete fPofB; + fPofB = 0; delete fPofT; fPofT = 0; } @@ -54,6 +76,8 @@ TLondon1D1L::~TLondon1D1L() { fParForPofT.clear(); delete fImpProfile; fImpProfile = 0; + delete fPofB; + fPofB = 0; delete fPofT; fPofT = 0; } @@ -65,6 +89,8 @@ TLondon1D2L::~TLondon1D2L() { fParForPofT.clear(); delete fImpProfile; fImpProfile = 0; + delete fPofB; + fPofB = 0; delete fPofT; fPofT = 0; } @@ -76,6 +102,8 @@ TProximity1D1LHS::~TProximity1D1LHS() { fParForPofT.clear(); delete fImpProfile; fImpProfile = 0; + delete fPofB; + fPofB = 0; delete fPofT; fPofT = 0; } @@ -87,6 +115,8 @@ TProximity1D1LHSGss::~TProximity1D1LHSGss() { fParForPofT.clear(); delete fImpProfile; fImpProfile = 0; + delete fPofB; + fPofB = 0; delete fPofT; fPofT = 0; } @@ -98,6 +128,8 @@ TLondon1D3L::~TLondon1D3L() { fParForPofT.clear(); delete fImpProfile; fImpProfile = 0; + delete fPofB; + fPofB = 0; delete fPofT; fPofT = 0; } @@ -109,6 +141,8 @@ TLondon1D3LS::~TLondon1D3LS() { fParForPofT.clear(); delete fImpProfile; fImpProfile = 0; + delete fPofB; + fPofB = 0; delete fPofT; fPofT = 0; } @@ -131,6 +165,8 @@ TLondon1D3LSub::~TLondon1D3LSub() { fParForPofT.clear(); delete fImpProfile; fImpProfile = 0; + delete fPofB; + fPofB = 0; delete fPofT; fPofT = 0; } @@ -159,7 +195,7 @@ TLondon1DHS::TLondon1DHS() : fCalcNeeded(true), fFirstCall(true) { string rge_path(startupHandler->GetDataPath()); vector energy_vec(startupHandler->GetEnergyList()); - fParForPofT.push_back(0.0); + fParForPofT.push_back(0.0); // phase fParForPofT.push_back(startupHandler->GetDeltat()); fParForPofT.push_back(startupHandler->GetDeltaB()); @@ -173,12 +209,15 @@ TLondon1DHS::TLondon1DHS() : fCalcNeeded(true), fFirstCall(true) { TTrimSPData *x = new TTrimSPData(rge_path, energy_vec); fImpProfile = x; x = 0; - delete x; - TPofTCalc *y = new TPofTCalc(fWisdom, fParForPofT); - fPofT = y; + TPofBCalc *y = new TPofBCalc(fParForPofB); + fPofB = y; y = 0; - delete y; + + TPofTCalc *z = new TPofTCalc(fPofB, fWisdom, fParForPofT); + fPofT = z; + z = 0; + // clean up if (saxParser) { @@ -212,13 +251,8 @@ double TLondon1DHS::operator()(double t, const vector &par) const { if(fFirstCall){ fPar = par; -// for (unsigned int i(0); i &par) const { fParForPofB[2] = par[1]; // energy TLondon1D_1L BofZ1(fParForBofZ); - TPofBCalc PofB1(BofZ1, *fImpProfile, fParForPofB); - fPofT->DoFFT(PofB1); + fPofB->UnsetPBExists(); + fPofB->Calculate(&BofZ1, fImpProfile, fParForPofB); + fPofT->DoFFT(); }/* else { cout << "Only the phase parameter has changed, (re-)calculating P(t) now..." << endl; @@ -480,12 +518,15 @@ TLondon1D2L::TLondon1D2L() : fCalcNeeded(true), fFirstCall(true), fLastTwoChange TTrimSPData *x = new TTrimSPData(rge_path, energy_vec); fImpProfile = x; x = 0; - delete x; - TPofTCalc *y = new TPofTCalc(fWisdom, fParForPofT); - fPofT = y; + TPofBCalc *y = new TPofBCalc(fParForPofB); + fPofB = y; y = 0; - delete y; + + TPofTCalc *z = new TPofTCalc(fPofB, fWisdom, fParForPofT); + fPofT = z; + z = 0; + // clean up if (saxParser) { @@ -577,8 +618,10 @@ double TLondon1D2L::operator()(double t, const vector &par) const { } TLondon1D_2L BofZ2(fParForBofZ); - TPofBCalc PofB2(BofZ2, *fImpProfile, fParForPofB); - fPofT->DoFFT(PofB2); + fPofB->UnsetPBExists(); + fPofB->Calculate(&BofZ2, fImpProfile, fParForPofB); + fPofT->DoFFT(); + }/* else { cout << "Only the phase parameter has changed, (re-)calculating P(t) now..." << endl; @@ -633,10 +676,14 @@ TProximity1D1LHS::TProximity1D1LHS() : fCalcNeeded(true), fFirstCall(true) { fImpProfile = x; x = 0; - TPofTCalc *y = new TPofTCalc(fWisdom, fParForPofT); - fPofT = y; + TPofBCalc *y = new TPofBCalc(fParForPofB); + fPofB = y; y = 0; + TPofTCalc *z = new TPofTCalc(fPofB, fWisdom, fParForPofT); + fPofT = z; + z = 0; + // clean up if (saxParser) { delete saxParser; @@ -738,8 +785,10 @@ double TProximity1D1LHS::operator()(double t, const vector &par) const { } TProximity1D_1LHS BofZ(fParForBofZ); - TPofBCalc PofB(BofZ, *fImpProfile, fParForPofB); - fPofT->DoFFT(PofB); + fPofB->UnsetPBExists(); + fPofB->Calculate(&BofZ, fImpProfile, fParForPofB); + fPofT->DoFFT(); + }/* else { cout << "Only the phase parameter has changed, (re-)calculating P(t) now..." << endl; @@ -791,10 +840,14 @@ TProximity1D1LHSGss::TProximity1D1LHSGss() : fCalcNeeded(true), fFirstCall(true) fImpProfile = x; x = 0; - TPofTCalc *y = new TPofTCalc(fWisdom, fParForPofT); - fPofT = y; + TPofBCalc *y = new TPofBCalc(fParForPofB); + fPofB = y; y = 0; + TPofTCalc *z = new TPofTCalc(fPofB, fWisdom, fParForPofT); + fPofT = z; + z = 0; + // clean up if (saxParser) { delete saxParser; @@ -872,8 +925,9 @@ double TProximity1D1LHSGss::operator()(double t, const vector &par) cons fParForPofB[2] = par[1]; // energy TProximity1D_1LHSGss BofZ(fParForBofZ); - TPofBCalc PofB(BofZ, *fImpProfile, fParForPofB); - fPofT->DoFFT(PofB); + fPofB->UnsetPBExists(); + fPofB->Calculate(&BofZ, fImpProfile, fParForPofB); + fPofT->DoFFT(); }/* else { cout << "Only the phase parameter has changed, (re-)calculating P(t) now..." << endl; @@ -922,12 +976,15 @@ TLondon1D3L::TLondon1D3L() : fCalcNeeded(true), fFirstCall(true), fLastThreeChan TTrimSPData *x = new TTrimSPData(rge_path, energy_vec); fImpProfile = x; x = 0; - delete x; - TPofTCalc *y = new TPofTCalc(fWisdom, fParForPofT); - fPofT = y; + TPofBCalc *y = new TPofBCalc(fParForPofB); + fPofB = y; y = 0; - delete y; + + TPofTCalc *z = new TPofTCalc(fPofB, fWisdom, fParForPofT); + fPofT = z; + z = 0; + // clean up if (saxParser) { @@ -1034,8 +1091,10 @@ double TLondon1D3L::operator()(double t, const vector &par) const { } TLondon1D_3L BofZ3(fParForBofZ); - TPofBCalc PofB3(BofZ3, *fImpProfile, fParForPofB); - fPofT->DoFFT(PofB3); + fPofB->UnsetPBExists(); + fPofB->Calculate(&BofZ3, fImpProfile, fParForPofB); + fPofT->DoFFT(); + }/* else { cout << "Only the phase parameter has changed, (re-)calculating P(t) now..." << endl; @@ -1085,12 +1144,14 @@ TLondon1D3LS::TLondon1D3LS() : fCalcNeeded(true), fFirstCall(true), fLastThreeCh TTrimSPData *x = new TTrimSPData(rge_path, energy_vec); fImpProfile = x; x = 0; - delete x; - TPofTCalc *y = new TPofTCalc(fWisdom, fParForPofT); - fPofT = y; + TPofBCalc *y = new TPofBCalc(fParForPofB); + fPofB = y; y = 0; - delete y; + + TPofTCalc *z = new TPofTCalc(fPofB, fWisdom, fParForPofT); + fPofT = z; + z = 0; // clean up if (saxParser) { @@ -1183,8 +1244,9 @@ double TLondon1D3LS::operator()(double t, const vector &par) const { } TLondon1D_3LS BofZ3S(fParForBofZ); - TPofBCalc PofB3S(BofZ3S, *fImpProfile, fParForPofB); - fPofT->DoFFT(PofB3S); + fPofB->UnsetPBExists(); + fPofB->Calculate(&BofZ3S, fImpProfile, fParForPofB); + fPofT->DoFFT(); }/* else { cout << "Only the phase parameter has changed, (re-)calculating P(t) now..." << endl; @@ -1402,12 +1464,14 @@ TLondon1D3LSub::TLondon1D3LSub() : fCalcNeeded(true), fFirstCall(true), fWeights TTrimSPData *x = new TTrimSPData(rge_path, energy_vec); fImpProfile = x; x = 0; - delete x; - TPofTCalc *y = new TPofTCalc(fWisdom, fParForPofT); - fPofT = y; + TPofBCalc *y = new TPofBCalc(fParForPofB); + fPofB = y; y = 0; - delete y; + + TPofTCalc *z = new TPofTCalc(fPofB, fWisdom, fParForPofT); + fPofT = z; + z = 0; // clean up if (saxParser) { @@ -1519,13 +1583,14 @@ double TLondon1D3LSub::operator()(double t, const vector &par) const { } TLondon1D_3L BofZ3(fParForBofZ); - TPofBCalc PofB3(BofZ3, *fImpProfile, fParForPofB); + fPofB->UnsetPBExists(); + fPofB->Calculate(&BofZ3, fImpProfile, fParForPofB); // Add background contribution from the substrate - PofB3.AddBackground(par[2], par[14], fImpProfile->LayerFraction(par[1], 4, interfaces)); + fPofB->AddBackground(par[2], par[14], fImpProfile->LayerFraction(par[1], 4, interfaces)); // FourierTransform of P(B) - fPofT->DoFFT(PofB3); + fPofT->DoFFT(); }/* else { cout << "Only the phase parameter has changed, (re-)calculating P(t) now..." << endl; diff --git a/src/external/TFitPofB-lib/classes/TPofBCalc.cpp b/src/external/TFitPofB-lib/classes/TPofBCalc.cpp index 6b37f88a..18208ffc 100644 --- a/src/external/TFitPofB-lib/classes/TPofBCalc.cpp +++ b/src/external/TFitPofB-lib/classes/TPofBCalc.cpp @@ -9,6 +9,26 @@ ***************************************************************************/ +/*************************************************************************** + * Copyright (C) 2009 by Bastian M. Wojek * + * * + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program; if not, write to the * + * Free Software Foundation, Inc., * + * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. * + ***************************************************************************/ + #include "TPofTCalc.h" #include #include @@ -23,13 +43,39 @@ #include /-------------------------------------------------------*/ +TPofBCalc::TPofBCalc(const vector ¶) : fBmin(0.0), fBmax(0.0), fDT(para[0]), fDB(para[1]), fPBExists(false) { + fPBSize = static_cast(1.0/(gBar*fDT*fDB)); + if (fPBSize % 2) { + fPBSize += 1; + } else { + fPBSize += 2; + } + + fB = new double[fPBSize]; + fPB = new double[fPBSize]; + + int i; + + for (i = 0; i < static_cast(fPBSize); i++) { + fB[i] = static_cast(i)*fDB; + fPB[i] = 0.0; + } +} + + // Do not actually calculate P(B) but take it from a B and a P(B) vector of the same size -TPofBCalc::TPofBCalc( const vector& b, const vector& pb, double dt) { +TPofBCalc::TPofBCalc(const vector& b, const vector& pb, double dt) { assert(b.size() == pb.size() && b.size() >= 2); + fPBSize = pb.size(); - fB = b; - fPB = pb; + fB = new double[fPBSize]; + fPB = new double[fPBSize]; + + for (unsigned int i(0); i < fPBSize; i++) { + fB[i] = b[i]; + fPB[i] = pb[i]; + } vector::const_iterator iter, iterB; iterB = b.begin(); @@ -37,7 +83,7 @@ TPofBCalc::TPofBCalc( const vector& b, const vector& pb, double for(iter = pb.begin(); iter != pb.end(); iter++){ if(*iter != 0.0) { fBmin = *iterB; - cout << fBmin << endl; +// cout << fBmin << endl; break; } iterB++; @@ -46,7 +92,7 @@ TPofBCalc::TPofBCalc( const vector& b, const vector& pb, double for( ; iter != b.end(); iter++){ if(*iter == 0.0) { fBmax = *(iterB-1); - cout << fBmax << endl; +// cout << fBmax << endl; break; } iterB++; @@ -54,130 +100,93 @@ TPofBCalc::TPofBCalc( const vector& b, const vector& pb, double fDT = dt; // needed if a convolution should be done fDB = b[1]-b[0]; - - cout << fDB << endl; + +// cout << fDB << endl; + + fPBExists = true; } - -TPofBCalc::TPofBCalc( const string &type, const vector ¶ ) : fDT(para[0]), fDB(para[1]) { - -if (type == "skg"){ // skewed Gaussian - - fBmin = 0.0; - fBmax = para[2]/gBar+10.0*fabs(para[4])/(2.0*pi*gBar); - - double BB; - double expominus(para[3]*para[3]/(2.0*pi*pi*gBar*gBar)); - double expoplus(para[4]*para[4]/(2.0*pi*pi*gBar*gBar)); - double B0(para[2]/gBar); - double BmaxFFT(1.0/gBar/fDT); - - for ( BB = 0.0 ; BB < B0 ; BB += fDB ) { - fB.push_back(BB); - fPB.push_back(exp(-(BB-B0)*(BB-B0)/expominus)); +void TPofBCalc::UnsetPBExists() { + int i; +#pragma omp parallel for default(shared) private(i) schedule(dynamic) + for (i = 0; i < static_cast(fPBSize); i++) { + fPB[i] = 0.0; } - - for ( ; BB < fBmax ; BB += fDB ) { - fB.push_back(BB); - fPB.push_back(exp(-(BB-B0)*(BB-B0)/expoplus)); - } - - unsigned int lastZerosStart(fB.size()); - - for ( ; BB <= BmaxFFT ; BB += fDB ) { - fB.push_back(BB); - fPB.push_back(0.0); - } - - // make sure that we have an even number of elements in p(B) for FFT, so we do not have to care later - - if (fB.size() % 2) { - fB.push_back(BB); - fPB.push_back(0.0); - } - - // normalize p(B) - - double pBsum = 0.0; - for (unsigned int i(0); i ¶) { + + if (type == "skg"){ // skewed Gaussian + + fBmin = 0.0; + fBmax = para[2]/gBar+10.0*fabs(para[4])/(2.0*pi*gBar); + + int a3(static_cast(floor(fBmax/fDB))); + int a4(static_cast(ceil(fBmax/fDB))); + + unsigned int BmaxIndex((a3 < a4) ? a4 : (a4 + 1)); + unsigned int B0Index(static_cast(ceil(para[2]/(gBar*fDB)))); + + double expominus(para[3]*para[3]/(2.0*pi*pi*gBar*gBar)); + double expoplus(para[4]*para[4]/(2.0*pi*pi*gBar*gBar)); + double B0(para[2]/(gBar)); + + for (unsigned int i(0); i < B0Index; i++) { + fPB[i] = exp(-(fB[i]-B0)*(fB[i]-B0)/expominus); + } + + for (unsigned int i(B0Index); i <= BmaxIndex; i++) { + fPB[i] = exp(-(fB[i]-B0)*(fB[i]-B0)/expoplus); + } + + // normalize p(B) + + double pBsum = 0.0; + for (unsigned int i(0); i <= BmaxIndex; i++) + pBsum += fPB[i]; + pBsum *= fDB; + for (unsigned int i(0); i <= BmaxIndex; i++) + fPB[i] /= pBsum; + + } + + fPBExists = true; + return; } //----------- -// Constructor that does the P(B) calculation for given analytical inverse of B(z) and its derivative and n(z) +// Calculate-method that does the P(B) calculation for given analytical inverse of B(z) and its derivative and n(z) // Parameters: dt[us], dB[G], Energy[keV], Bbg[G], width[us^{-1}], weight[1] //----------- -TPofBCalc::TPofBCalc( const TBofZCalcInverse &BofZ, const TTrimSPData &dataTrimSP, const vector ¶ ) : fDT(para[0]), fDB(para[1]) -{ +void TPofBCalc::Calculate(const TBofZCalcInverse *BofZ, const TTrimSPData *dataTrimSP, const vector ¶) { - fBmin = BofZ.GetBmin(); - fBmax = BofZ.GetBmax(); + if(fPBExists) + return; - double BB; + fBmin = BofZ->GetBmin(); + fBmax = BofZ->GetBmax(); - // fill not used Bs before Bmin with 0.0 + int a1(static_cast(floor(fBmin/fDB))); + int a2(static_cast(ceil(fBmin/fDB))); + int a3(static_cast(floor(fBmax/fDB))); + int a4(static_cast(ceil(fBmax/fDB))); - for (BB = 0.0; BB < fBmin; BB += fDB) - fB.push_back(BB); - fPB.resize(fB.size(),0.0); + unsigned int firstZerosEnd ((a1 < a2) ? a1 : ((a1 > 0) ? (a1 - 1) : 0)); + unsigned int lastZerosStart ((a3 < a4) ? a4 : (a4 + 1)); - unsigned int firstZerosEnd(fB.size()); + unsigned int i; // calculate p(B) from the inverse of B(z) - for ( ; BB <= fBmax ; BB += fDB) { - fB.push_back(BB); - fPB.push_back(0.0); + for (i = firstZerosEnd; i <= lastZerosStart; i++) { vector< pair > inv; - inv = BofZ.GetInverseAndDerivative(BB); + inv = BofZ->GetInverseAndDerivative(fB[i]); - for (unsigned int i(0); iGetNofZ(inv[j].first, para[2])*fabs(inv[j].second); } // normalize p(B) @@ -194,35 +203,35 @@ TPofBCalc::TPofBCalc( const TBofZCalcInverse &BofZ, const TTrimSPData &dataTrimS } - - //----------- -// Constructor that does the P(B) calculation for given B(z) and n(z) +// Calculate-method that does the P(B) calculation for given B(z) and n(z) // Parameters: dt[us], dB[G], Energy[keV] //----------- -TPofBCalc::TPofBCalc( const TBofZCalc &BofZ, const TTrimSPData &dataTrimSP, const vector ¶, unsigned int zonk ) : fDT(para[0]), fDB(para[1]) { +void TPofBCalc::Calculate(const TBofZCalc *BofZ, const TTrimSPData *dataTrimSP, const vector ¶, unsigned int zonk) { - fBmin = BofZ.GetBmin(); - fBmax = BofZ.GetBmax(); + if(fPBExists) + return; + + fBmin = BofZ->GetBmin(); + fBmax = BofZ->GetBmax(); + + int a1(static_cast(floor(fBmin/fDB))); + int a2(static_cast(ceil(fBmin/fDB))); + 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 lastZerosStart ((a3 < a4) ? a4 : (a4 + 1)); double BB, BBnext; double zm, zp, zNext, dz; - // fill not used Bs before Bmin with 0.0 - - for (BB = 0.0; BB < fBmin; BB += fDB) { - fB.push_back(BB); - fPB.push_back(0.0); - } - - unsigned int firstZerosEnd(fB.size()); - // calculate p(B) from B(z) - vector bofzZ(BofZ.DataZ()); - vector bofzBZ(BofZ.DataBZ()); - double ddZ(BofZ.GetDZ()); + vector *bofzZ = BofZ->DataZ(); + vector *bofzBZ = BofZ->DataBZ(); + double ddZ(BofZ->GetDZ()); /* USED FOR DEBUGGING----------------------------------- cout << "Bmin = " << fBmin << ", Bmax = " << fBmax << endl; @@ -281,19 +290,20 @@ TPofBCalc::TPofBCalc( const TBofZCalc &BofZ, const TTrimSPData &dataTrimSP, cons double nn; bool zNextFound(false); - for ( ; BB <= fBmax ; BB += fDB) { - BBnext = BB + fDB; - fB.push_back(BB); - fPB.push_back(0.0); + unsigned int i; - for ( unsigned int j(0); j < bofzZ.size() - 1; j++ ) { - if ( bofzBZ[j] >= BB && bofzBZ[j+1] <= BB ) { - zm = (BB-bofzBZ[j])*ddZ/(bofzBZ[j+1]-bofzBZ[j]) + bofzZ[j]; + for (i = 0; i <= lastZerosStart; i++) { + BB = fB[i]; + BBnext = fB[i+1]; + + for ( unsigned int j(0); j < bofzZ->size() - 1; j++ ) { + if ( (*bofzBZ)[j] >= BB && (*bofzBZ)[j+1] <= BB ) { + zm = (BB-(*bofzBZ)[j])*ddZ/((*bofzBZ)[j+1]-(*bofzBZ)[j]) + (*bofzZ)[j]; for (unsigned int k(0); k < j; k++) { - if ( ( bofzBZ[j-k] <= BBnext && bofzBZ[j-k-1] >= BBnext ) ) { + if ( ( (*bofzBZ)[j-k] <= BBnext && (*bofzBZ)[j-k-1] >= BBnext ) ) { // cout << "1 " << j << " " << k << endl; - zNext = (BBnext-bofzBZ[j-k-1])*ddZ/(bofzBZ[j-k]-bofzBZ[j-k-1]) + bofzZ[j-k-1]; + zNext = (BBnext-(*bofzBZ)[j-k-1])*ddZ/((*bofzBZ)[j-k]-(*bofzBZ)[j-k-1]) + (*bofzZ)[j-k-1]; zNextFound = true; break; } @@ -303,20 +313,20 @@ TPofBCalc::TPofBCalc( const TBofZCalc &BofZ, const TTrimSPData &dataTrimSP, cons zNextFound = false; dz = zNext-zm; - nn = dataTrimSP.GetNofZ(zm, para[2]); + nn = dataTrimSP->GetNofZ(zm, para[2]); if (nn != -1.0) { // cout << "zNext = " << zNextm << ", zm = " << zm << ", dz = " << dz << endl; - *(fPB.end()-1) += nn*fabs(dz/fDB); + fPB[i] += nn*fabs(dz/fDB); } } - } else if (bofzBZ[j] <= BB && bofzBZ[j+1] >= BB ) { - zp = (BB-bofzBZ[j])*ddZ/(bofzBZ[j+1]-bofzBZ[j]) + bofzZ[j]; + } else if ((*bofzBZ)[j] <= BB && (*bofzBZ)[j+1] >= BB ) { + zp = (BB-(*bofzBZ)[j])*ddZ/((*bofzBZ)[j+1]-(*bofzBZ)[j]) + (*bofzZ)[j]; - for (unsigned int k(0); k < bofzZ.size() - j - 1; k++) { - if ( ( bofzBZ[j+k] <= BBnext && bofzBZ[j+k+1] >= BBnext ) ) { + for (unsigned int k(0); k < bofzZ->size() - j - 1; k++) { + if ( ( (*bofzBZ)[j+k] <= BBnext && (*bofzBZ)[j+k+1] >= BBnext ) ) { // cout << "2 " << j << " " << k << endl; - zNext = (BBnext-bofzBZ[j+k])*ddZ/(bofzBZ[j+k+1]-bofzBZ[j+k]) + bofzZ[j+k]; + zNext = (BBnext-(*bofzBZ)[j+k])*ddZ/((*bofzBZ)[j+k+1]-(*bofzBZ)[j+k]) + (*bofzZ)[j+k]; zNextFound = true; break; } @@ -326,10 +336,10 @@ TPofBCalc::TPofBCalc( const TBofZCalc &BofZ, const TTrimSPData &dataTrimSP, cons zNextFound = false; dz = zNext-zp; - nn = dataTrimSP.GetNofZ(zp, para[2]); + nn = dataTrimSP->GetNofZ(zp, para[2]); if (nn != -1.0) { // cout << "zNext = " << zNextp << ", zp = " << zp << ", dz = " << dz << endl; - *(fPB.end()-1) += nn*fabs(dz/fDB); + fPB[i] += nn*fabs(dz/fDB); } } } @@ -337,114 +347,90 @@ TPofBCalc::TPofBCalc( const TBofZCalc &BofZ, const TTrimSPData &dataTrimSP, cons } - unsigned int lastZerosStart(fB.size()); - - // fill not used Bs after Bext with 0.0 - - double BmaxFFT(1.0/gBar/fDT); - -// cout << "N = " << int(BmaxFFT/para[1]+1.0) << endl; - - for ( ; BB <= BmaxFFT ; BB += fDB ) { - fB.push_back(BB); - fPB.push_back(0.0); - } - - // make sure that we have an even number of elements in p(B) for FFT, so we do not have to care later - - if (fB.size() % 2) { - fB.push_back(BB); - fPB.push_back(0.0); - } - -// cout << "size of B = " << fB.size() << ", size of p(B) = " << fPB.size() << endl; + bofzZ = 0; + bofzBZ = 0; // normalize p(B) double pBsum = 0.0; - for (unsigned int i(firstZerosEnd); i ¶ ) : fDT(para[0]), fDB(para[1]) -{ +void TPofBCalc::Calculate(const TBulkVortexFieldCalc *vortexLattice, const vector ¶) { - fBmin = vortexLattice.GetBmin(); - fBmax = vortexLattice.GetBmax(); + if(fPBExists) + return; - unsigned int a1(static_cast(floor(fBmin/fDB))); - unsigned int a2(static_cast(ceil(fBmin/fDB))); - unsigned int a3(static_cast(floor(fBmax/fDB))); - unsigned int a4(static_cast(ceil(fBmax/fDB))); + fBmin = vortexLattice->GetBmin(); + fBmax = vortexLattice->GetBmax(); - unsigned int firstZerosEnd ((a1 != a2) ? a1 : (a1 - 1)); - unsigned int lastZerosStart ((a3 != a4) ? a4 : (a4 + 1)); + int a1(static_cast(floor(fBmin/fDB))); + int a2(static_cast(ceil(fBmin/fDB))); + int a3(static_cast(floor(fBmax/fDB))); + int a4(static_cast(ceil(fBmax/fDB))); - // fill the B vector - unsigned int indexBmaxFFT(static_cast(ceil(1.0/(gBar*fDT*fDB)))); - fB.resize(indexBmaxFFT, 0.0); - fPB.resize(indexBmaxFFT, 0.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); - int i; -#pragma omp parallel for default(shared) private(i) schedule(dynamic) - for (i = 0; i < static_cast(indexBmaxFFT); i++) - fB[i] = fDB*static_cast(i); -// end pragma omp parallel + if (lastZerosStart >= fPBSize) + lastZerosStart = fPBSize - 1; - // make sure that we have an even number of elements in p(B) for FFT, so we do not have to care later - if (fB.size() % 2) { - fB.push_back(*(fB.end()-1)+fDB); - fPB.push_back(0.0); +// cout << endl << fBmin << " " << fBmax << " " << firstZerosEnd << " " << lastZerosStart << " " << numberOfPoints << endl; + + if (!vortexLattice->GridExists()) { + vortexLattice->CalculateGrid(); } - vector< vector > vortexFields(vortexLattice.DataB()); + double *vortexFields = vortexLattice->DataB(); + unsigned int fill_index, counter(0); - if (vortexFields.empty()) { - vortexLattice.CalculateGrid(); - vortexFields = vortexLattice.DataB(); - } - - unsigned int numberOfPoints(vortexFields.size()); - - for (unsigned int j(0); j < numberOfPoints; j++){ - for (unsigned int k(0); k < numberOfPoints; k++){ - fPB[static_cast(ceil(vortexFields[j][k]/fDB))] += 1.0; + for (unsigned int j(0); j < numberOfStepsSq_2; j++) { + fill_index = static_cast(ceil(abs((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; } } - double normalizer(static_cast(numberOfPoints*numberOfPoints)*fDB); + vortexFields = 0; + double normalizer(static_cast(numberOfStepsSq_2/2)*fDB); + + int i; #pragma omp parallel for default(shared) private(i) schedule(dynamic) for (i = firstZerosEnd; i <= static_cast(lastZerosStart); i++) { fPB[i] /= normalizer; } // end pragma omp parallel -// // normalize p(B) -// -// double pBsum = 0.0; -// for (unsigned int i(firstZerosEnd); i1.0 || B<0.0) return; - unsigned int sizePB(fPB.size()); double BsSq(s*s/(gBar*gBar*4.0*pi*pi)); // calculate Gaussian background - vector bg; - for(unsigned int i(0); i < sizePB; i++) { - bg.push_back(exp(-(fB[i]-B)*(fB[i]-B)/(2.0*BsSq))); + double bg[fPBSize]; + for(unsigned int i(0); i < fPBSize; i++) { + bg[i] = exp(-(fB[i]-B)*(fB[i]-B)/(2.0*BsSq)); } // normalize background double bgsum(0.0); - for (unsigned int i(0); i < sizePB; i++) + for (unsigned int i(0); i < fPBSize; i++) bgsum += bg[i]; bgsum *= fDB; - for (unsigned int i(0); i < sizePB; i++) + for (unsigned int i(0); i < fPBSize; i++) bg[i] /= bgsum; // add background to P(B) - for (unsigned int i(0); i < sizePB; i++) + for (unsigned int i(0); i < fPBSize; i++) fPB[i] = (1.0 - w)*fPB[i] + w*bg[i]; // // check if normalization is still valid @@ -487,30 +472,19 @@ void TPofBCalc::ConvolveGss(double w) { if(!w) return; - unsigned int NFFT; + unsigned int NFFT(fPBSize); double TBin; fftw_plan FFTplanToTimeDomain; fftw_plan FFTplanToFieldDomain; fftw_complex *FFTout; - double *FFTin; - NFFT = ( int(1.0/gBar/fDT/fDB+1.0) % 2 ) ? int(1.0/gBar/fDT/fDB+2.0) : int(1.0/gBar/fDT/fDB+1.0); - TBin = 1.0/gBar/double(NFFT-1)/fDB; + TBin = 1.0/(gBar*double(NFFT-1)*fDB); - FFTin = (double *)malloc(sizeof(double) * NFFT); - FFTout = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * (NFFT/2+1)); + FFTout = new fftw_complex[NFFT/2 + 1]; //(fftw_complex *)fftw_malloc(sizeof(fftw_complex) * (NFFT/2+1)); // do the FFT to time domain - FFTplanToTimeDomain = fftw_plan_dft_r2c_1d(NFFT, FFTin, FFTout, FFTW_ESTIMATE); - - for (unsigned int i(0); i(i*i)); FFTout[i][0] *= GssInTimeDomain; FFTout[i][1] *= GssInTimeDomain; } // FFT back to the field domain - FFTplanToFieldDomain = fftw_plan_dft_c2r_1d(NFFT, FFTout, FFTin, FFTW_ESTIMATE); + FFTplanToFieldDomain = fftw_plan_dft_c2r_1d(NFFT, FFTout, fPB, FFTW_ESTIMATE); fftw_execute(FFTplanToFieldDomain); - for (unsigned int i(0); i &par) : fWisdom(wisdom) { +TPofTCalc::TPofTCalc (const TPofBCalc *PofB, const string &wisdom, const vector &par) : fWisdom(wisdom) { int init_threads(fftw_init_threads()); if (!init_threads) @@ -71,88 +71,95 @@ TPofTCalc::TPofTCalc (const string &wisdom, const vector &par) : fWisdom else fftw_plan_with_nthreads(2); - fNFFT = ( int(1.0/gBar/par[1]/par[2]+1.0) % 2 ) ? int(1.0/gBar/par[1]/par[2]+2.0) : int(1.0/gBar/par[1]/par[2]+1.0); - fTBin = 1.0/gBar/double(fNFFT-1)/par[2]; + fNFFT = static_cast(1.0/(gBar*par[1]*par[2])); + if (fNFFT % 2) { + fNFFT += 1; + } else { + fNFFT += 2; + } - fT.resize(fNFFT/2+1); - fPT.resize(fNFFT/2+1); + fTBin = 1.0/(gBar*double(fNFFT-1)*par[2]); + + int NFFT_2p1(fNFFT/2 + 1); + + // allocating memory for the time- and polarisation vectors + + fT = new double[NFFT_2p1]; //static_cast(malloc(sizeof(double) * NFFT_2p1)); + fPT = new double[NFFT_2p1]; //static_cast(malloc(sizeof(double) * NFFT_2p1)); int i; #pragma omp parallel for default(shared) private(i) schedule(dynamic) - for (i=0; i(i)*fTBin; } - fFFTin = (double *)malloc(sizeof(double) * fNFFT); - fFFTout = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * (fNFFT/2+1)); + fFFTin = PofB->DataPB(); + fFFTout = new fftw_complex[NFFT_2p1]; //static_cast(fftw_malloc(sizeof(fftw_complex) * NFFT_2p1)); - cout << "TPofTCalc::TPofTCalc: Check for the FFT plan..." << endl; - - // Load wisdom from file + // Load wisdom from file if it exists and should be used + fUseWisdom = true; int wisdomLoaded(0); FILE *wordsOfWisdomR; - wordsOfWisdomR = fopen(fWisdom.c_str(), "r"); + wordsOfWisdomR = fopen(wisdom.c_str(), "r"); if (wordsOfWisdomR == NULL) { - cout << "TPofTCalc::TPofTCalc: Couldn't open wisdom file ..." << endl; + fUseWisdom = false; } else { wisdomLoaded = fftw_import_wisdom_from_file(wordsOfWisdomR); fclose(wordsOfWisdomR); } if (!wisdomLoaded) { - cout << "TPofTCalc::TPofTCalc: No wisdom is imported..." << endl; + fUseWisdom = false; } - fFFTplan = fftw_plan_dft_r2c_1d(fNFFT, fFFTin, fFFTout, FFTW_EXHAUSTIVE); +// create the FFT plan -// cout << &fFFTplan << endl; + if (fUseWisdom) + fFFTplan = fftw_plan_dft_r2c_1d(fNFFT, fFFTin, fFFTout, FFTW_EXHAUSTIVE); + else + fFFTplan = fftw_plan_dft_r2c_1d(fNFFT, fFFTin, fFFTout, FFTW_ESTIMATE); + +} + +//--------------------- +// Destructor of the TPofTCalc class - it saves the FFT plan and cleans up +//--------------------- + +TPofTCalc::~TPofTCalc() { + // if a wisdom file is used export the wisdom so it has not to be checked for the FFT-plan next time + if (fUseWisdom) { + FILE *wordsOfWisdomW; + wordsOfWisdomW = fopen(fWisdom.c_str(), "w"); + if (wordsOfWisdomW == NULL) { + cout << "TPofTCalc::~TPofTCalc(): Could not open file ... No wisdom is exported..." << endl; + } else { + fftw_export_wisdom_to_file(wordsOfWisdomW); + fclose(wordsOfWisdomW); + } + } + + // clean up + + fftw_destroy_plan(fFFTplan); + delete[] fFFTout; //fftw_free(fFFTout); + fFFTout = 0; +// fftw_cleanup(); +// fftw_cleanup_threads(); + + delete[] fT; + fT = 0; + delete[] fPT; + fPT = 0; } //-------------- // Method that does the FFT of a given p(B) //-------------- -void TPofTCalc::DoFFT(const TPofBCalc &PofB) { - - vector pB(PofB.DataPB()); - -/* USED FOR DEBUGGING ----------------------- - - time_t seconds; - seconds = time(NULL); - - vector B(PofB.DataB()); - double Bmin(PofB.GetBmin()); - - char debugfile[50]; - int n = sprintf (debugfile, "test_PB_%ld_%f.dat", seconds, Bmin); - - if (n > 0) { - ofstream of(debugfile); - - for (unsigned int i(0); i &par) { #pragma omp parallel for default(shared) private(i) schedule(dynamic) for (i=0; i(t/fTBin)); + if (i < fNFFT/2){ + return fPT[i]+(fPT[i+1]-fPT[i])/(fT[i+1]-fT[i])*(t-fT[i]); + } + cout << "TPofTCalc::Eval: No data for the time " << t << " us available! Returning -999.0 ..." << endl; + return -999.0; +} + + + //--------------------- // Method for generating fake LEM decay histograms from p(B) // Parameters: output filename, par(dt, dB, timeres, channels, asyms, phases, t0s, N0s, bgs) @@ -381,49 +403,5 @@ void TPofTCalc::FakeData(const string &rootOutputFileName, const vector return; } -//--------------------- -// Method for evaluating P(t) at a given t -//--------------------- -double TPofTCalc::Eval(double t) const { - - unsigned int i(int(t/fTBin)); - if (i #include @@ -27,6 +47,8 @@ TSkewedGss::~TSkewedGss() { fPar.clear(); fParForPofB.clear(); fParForPofT.clear(); + delete fPofB; + fPofB = 0; delete fPofT; fPofT = 0; } @@ -61,10 +83,13 @@ TSkewedGss::TSkewedGss() : fCalcNeeded(true), fFirstCall(true) { fParForPofB.push_back(0.0); // s- fParForPofB.push_back(0.0); // s+ - TPofTCalc *y = new TPofTCalc(fWisdom, fParForPofT); + TPofBCalc *x = new TPofBCalc(fParForPofB); + fPofB = x; + x = 0; + + TPofTCalc *y = new TPofTCalc(fPofB, fWisdom, fParForPofT); fPofT = y; y = 0; - delete y; // clean up if (saxParser) { @@ -131,8 +156,8 @@ double TSkewedGss::operator()(double t, const vector &par) const { fParForPofB[3] = par[2]; // sigma- fParForPofB[4] = par[3]; // sigma+ - TPofBCalc PofB("skg", fParForPofB); - fPofT->DoFFT(PofB); + fPofB->Calculate("skg", fParForPofB); + fPofT->DoFFT(); }/* else { cout << "Only the phase parameter has changed, (re-)calculating P(t) now..." << endl; diff --git a/src/external/TFitPofB-lib/classes/TTrimSPDataHandler.cpp b/src/external/TFitPofB-lib/classes/TTrimSPDataHandler.cpp index a707ecfb..360bee81 100644 --- a/src/external/TFitPofB-lib/classes/TTrimSPDataHandler.cpp +++ b/src/external/TFitPofB-lib/classes/TTrimSPDataHandler.cpp @@ -9,6 +9,26 @@ ***************************************************************************/ +/*************************************************************************** + * Copyright (C) 2009 by Bastian M. Wojek * + * * + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program; if not, write to the * + * Free Software Foundation, Inc., * + * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. * + ***************************************************************************/ + #include "TTrimSPDataHandler.h" #include #include diff --git a/src/external/TFitPofB-lib/classes/TVortex.cpp b/src/external/TFitPofB-lib/classes/TVortex.cpp index 22c7ac8d..4aaa7e0a 100644 --- a/src/external/TFitPofB-lib/classes/TVortex.cpp +++ b/src/external/TFitPofB-lib/classes/TVortex.cpp @@ -9,6 +9,26 @@ ***************************************************************************/ +/*************************************************************************** + * Copyright (C) 2009 by Bastian M. Wojek * + * * + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program; if not, write to the * + * Free Software Foundation, Inc., * + * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. * + ***************************************************************************/ + #include "TVortex.h" #include #include @@ -25,12 +45,16 @@ ClassImp(TBulkTriVortexLondon) //------------------ TBulkTriVortexLondon::~TBulkTriVortexLondon() { + delete fPofT; + fPofT = 0; + delete fPofB; + fPofT = 0; + delete fVortex; + fVortex = 0; fPar.clear(); fParForVortex.clear(); fParForPofB.clear(); fParForPofT.clear(); - delete fPofT; - fPofT = 0; } //------------------ @@ -54,9 +78,8 @@ TBulkTriVortexLondon::TBulkTriVortexLondon() : fCalcNeeded(true), fFirstCall(tru fGridSteps = startupHandler->GetGridSteps(); fWisdom = startupHandler->GetWisdomFile(); - fVortexFourierComp = startupHandler->GetVortexFourierComp(); - fParForVortex.clear(); + fParForVortex.resize(3); // field, lambda, xi fParForPofT.push_back(0.0); fParForPofT.push_back(startupHandler->GetDeltat()); @@ -69,10 +92,18 @@ TBulkTriVortexLondon::TBulkTriVortexLondon() : fCalcNeeded(true), fFirstCall(tru fParForPofB.push_back(0.005); // Bkg-width fParForPofB.push_back(0.0); // Bkg-weight - TPofTCalc *y = new TPofTCalc(fWisdom, fParForPofT); - fPofT = y; + TBulkTriVortexLondonFieldCalc *x = new TBulkTriVortexLondonFieldCalc(fWisdom, fGridSteps); + fVortex = x; + x = 0; + + TPofBCalc *y = new TPofBCalc(fParForPofB); + fPofB = y; y = 0; + TPofTCalc *z = new TPofTCalc(fPofB, fWisdom, fParForPofT); + fPofT = z; + z = 0; + // clean up if (saxParser) { delete saxParser; @@ -87,7 +118,7 @@ TBulkTriVortexLondon::TBulkTriVortexLondon() : fCalcNeeded(true), fFirstCall(tru //------------------ // TBulkTriVortexLondon-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 TBulkTriVortexLondon (phase, appl.field, lambda, xi) +// Parameters: all the parameters for the function to be fitted through TBulkTriVortexLondon (phase, appl.field, lambda, xi, [not implemented: bkg weight]) //------------------ double TBulkTriVortexLondon::operator()(double t, const vector &par) const { @@ -102,8 +133,8 @@ double TBulkTriVortexLondon::operator()(double t, const vector &par) con if(fFirstCall){ fPar = par; - for (unsigned int i(1); i < 4; i++){ - fParForVortex.push_back(fPar[i]); + for (unsigned int i(0); i < 3; i++){ + fParForVortex[i] = fPar[i+1]; } fFirstCall = false; } @@ -138,16 +169,18 @@ double TBulkTriVortexLondon::operator()(double t, const vector &par) con // cout << " Parameters have changed, (re-)calculating p(B) and P(t) now..." << endl; - fParForVortex[0] = par[2]; // lambda - fParForVortex[1] = par[3]; // xi - fParForVortex[2] = par[1]; // field + for (unsigned int i(0); i < 3; i++) { + fParForVortex[i] = par[i+1]; + } fParForPofB[2] = par[1]; // Bkg-Field //fParForPofB[3] = 0.005; // Bkg-width (in principle zero) - TBulkTriVortexLondonFieldCalc vortexLattice(fParForVortex, fGridSteps, fVortexFourierComp); - TPofBCalc PofB(vortexLattice, fParForPofB); - fPofT->DoFFT(PofB); + fVortex->SetParameters(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; diff --git a/src/external/TFitPofB-lib/include/TBofZCalc.h b/src/external/TFitPofB-lib/include/TBofZCalc.h index edb2fc3d..81d5b382 100644 --- a/src/external/TFitPofB-lib/include/TBofZCalc.h +++ b/src/external/TFitPofB-lib/include/TBofZCalc.h @@ -9,6 +9,26 @@ ***************************************************************************/ +/*************************************************************************** + * Copyright (C) 2009 by Bastian M. Wojek * + * * + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program; if not, write to the * + * Free Software Foundation, Inc., * + * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. * + ***************************************************************************/ + #ifndef _TBofZCalc_H_ #define _TBofZCalc_H_ @@ -31,8 +51,8 @@ public: fParam.clear(); } - virtual vector DataZ() const {return fZ;} - virtual vector DataBZ() const {return fBZ;} + virtual vector* DataZ() const {return &fZ;} + virtual vector* DataBZ() const {return &fBZ;} virtual void Calculate(); virtual double GetBofZ(double) const = 0; virtual double GetBmin() const = 0; @@ -43,8 +63,8 @@ protected: int fSteps; double fDZ; vector fParam; - vector fZ; - vector fBZ; + mutable vector fZ; + mutable vector fBZ; }; //-------------------- diff --git a/src/external/TFitPofB-lib/include/TBulkVortexFieldCalc.h b/src/external/TFitPofB-lib/include/TBulkVortexFieldCalc.h index 552e982b..c5db18ec 100644 --- a/src/external/TFitPofB-lib/include/TBulkVortexFieldCalc.h +++ b/src/external/TFitPofB-lib/include/TBulkVortexFieldCalc.h @@ -9,12 +9,35 @@ ***************************************************************************/ +/*************************************************************************** + * Copyright (C) 2009 by Bastian M. Wojek * + * * + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program; if not, write to the * + * Free Software Foundation, Inc., * + * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. * + ***************************************************************************/ + #ifndef _TBulkVortexFieldCalc_H_ #define _TBulkVortexFieldCalc_H_ #include +#include using namespace std; +#include + //-------------------- // Base class for any kind of vortex symmetry //-------------------- @@ -25,43 +48,43 @@ public: TBulkVortexFieldCalc() {} - virtual ~TBulkVortexFieldCalc() { - for(unsigned int i(0); i < fB.size(); i++) - fB[i].clear(); - fB.clear(); - fParam.clear(); - } + virtual ~TBulkVortexFieldCalc(); - virtual vector< vector > DataB() const {return fB;} + virtual double* DataB() const {return fFFTout;} + virtual void SetParameters(const vector& par) {fParam = par; fGridExists = false;} virtual void CalculateGrid() const = 0; - virtual double GetBatPoint(double, double) const = 0; virtual double GetBmin() const = 0; virtual double GetBmax() const = 0; + virtual bool GridExists() const {return fGridExists;} + virtual void UnsetGridExists() const {fGridExists = false;} + virtual unsigned int GetNumberOfSteps() const {return fSteps;} protected: + vector fParam; unsigned int fSteps; // number of steps, the "unit cell" of the vortex lattice is devided in (in x- and y- direction) - vector fParam; // lambda, xi, appl.field - mutable vector< vector > fB; // fields in a "unit cell" of the vortex lattice + fftw_complex *fFFTin; // Fourier components of the field + double *fFFTout; // fields in a "unit cell" of the vortex lattice + fftw_plan fFFTplan; + bool fUseWisdom; + string fWisdom; + mutable bool fGridExists; }; //-------------------- -// Class for triangular vortex lattice, London model +// Class for triangular vortex lattice, London model with Gaussian Cutoff //-------------------- class TBulkTriVortexLondonFieldCalc : public TBulkVortexFieldCalc { public: - TBulkTriVortexLondonFieldCalc(const vector&, const unsigned int steps = 200, const unsigned int Ncomp = 17); + TBulkTriVortexLondonFieldCalc(const string&, const unsigned int steps = 256); ~TBulkTriVortexLondonFieldCalc() {} void CalculateGrid() const; - double GetBatPoint(double, double) const; double GetBmin() const; double GetBmax() const; -private: - unsigned int fNFourierComp; }; #endif // _TBulkVortexFieldCalc_H_ diff --git a/src/external/TFitPofB-lib/include/TFitPofBStartupHandlerLinkDef.h b/src/external/TFitPofB-lib/include/TFitPofBStartupHandlerLinkDef.h index 86d09b90..4218241e 100644 --- a/src/external/TFitPofB-lib/include/TFitPofBStartupHandlerLinkDef.h +++ b/src/external/TFitPofB-lib/include/TFitPofBStartupHandlerLinkDef.h @@ -1,6 +1,6 @@ /*************************************************************************** - TLondon1DLinkDef.h + TFitPofBStartupHandlerLinkDef.h Author: Bastian M. Wojek e-mail: bastian.wojek@psi.ch @@ -9,6 +9,26 @@ ***************************************************************************/ +/*************************************************************************** + * Copyright (C) 2009 by Bastian M. Wojek * + * * + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program; if not, write to the * + * Free Software Foundation, Inc., * + * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. * + ***************************************************************************/ + // root dictionary stuff -------------------------------------------------- #ifdef __CINT__ diff --git a/src/external/TFitPofB-lib/include/TLondon1D.h b/src/external/TFitPofB-lib/include/TLondon1D.h index 014281d3..87032ce7 100644 --- a/src/external/TFitPofB-lib/include/TLondon1D.h +++ b/src/external/TFitPofB-lib/include/TLondon1D.h @@ -9,6 +9,26 @@ ***************************************************************************/ +/*************************************************************************** + * Copyright (C) 2009 by Bastian M. Wojek * + * * + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program; if not, write to the * + * Free Software Foundation, Inc., * + * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. * + ***************************************************************************/ + #ifndef _TLondon1D_H_ #define _TLondon1D_H_ @@ -27,6 +47,7 @@ public: private: mutable vector fPar; TTrimSPData *fImpProfile; + TPofBCalc *fPofB; TPofTCalc *fPofT; mutable bool fCalcNeeded; mutable bool fFirstCall; @@ -51,6 +72,7 @@ public: private: mutable vector fPar; TTrimSPData *fImpProfile; + TPofBCalc *fPofB; TPofTCalc *fPofT; mutable bool fCalcNeeded; mutable bool fFirstCall; @@ -76,6 +98,7 @@ public: private: mutable vector fPar; TTrimSPData *fImpProfile; + TPofBCalc *fPofB; TPofTCalc *fPofT; mutable bool fCalcNeeded; mutable bool fFirstCall; @@ -101,6 +124,7 @@ public: private: mutable vector fPar; TTrimSPData *fImpProfile; + TPofBCalc *fPofB; TPofTCalc *fPofT; mutable bool fCalcNeeded; mutable bool fFirstCall; @@ -125,6 +149,7 @@ public: private: mutable vector fPar; TTrimSPData *fImpProfile; + TPofBCalc *fPofB; TPofTCalc *fPofT; mutable bool fCalcNeeded; mutable bool fFirstCall; @@ -150,6 +175,7 @@ public: private: mutable vector fPar; TTrimSPData *fImpProfile; + TPofBCalc *fPofB; TPofTCalc *fPofT; mutable bool fCalcNeeded; mutable bool fFirstCall; @@ -176,6 +202,7 @@ public: private: mutable vector fPar; TTrimSPData *fImpProfile; + TPofBCalc *fPofB; TPofTCalc *fPofT; mutable bool fCalcNeeded; mutable bool fFirstCall; @@ -226,6 +253,7 @@ public: private: mutable vector fPar; TTrimSPData *fImpProfile; + TPofBCalc *fPofB; TPofTCalc *fPofT; mutable bool fCalcNeeded; mutable bool fFirstCall; diff --git a/src/external/TFitPofB-lib/include/TLondon1DLinkDef.h b/src/external/TFitPofB-lib/include/TLondon1DLinkDef.h index f8ff47f6..bd816267 100644 --- a/src/external/TFitPofB-lib/include/TLondon1DLinkDef.h +++ b/src/external/TFitPofB-lib/include/TLondon1DLinkDef.h @@ -9,6 +9,26 @@ ***************************************************************************/ +/*************************************************************************** + * Copyright (C) 2009 by Bastian M. Wojek * + * * + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program; if not, write to the * + * Free Software Foundation, Inc., * + * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. * + ***************************************************************************/ + // root dictionary stuff -------------------------------------------------- #ifdef __CINT__ diff --git a/src/external/TFitPofB-lib/include/TPofBCalc.h b/src/external/TFitPofB-lib/include/TPofBCalc.h index ecd0ee0a..ab2ada53 100644 --- a/src/external/TFitPofB-lib/include/TPofBCalc.h +++ b/src/external/TFitPofB-lib/include/TPofBCalc.h @@ -9,6 +9,26 @@ ***************************************************************************/ +/*************************************************************************** + * Copyright (C) 2009 by Bastian M. Wojek * + * * + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program; if not, write to the * + * Free Software Foundation, Inc., * + * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. * + ***************************************************************************/ + #ifndef _TPofBCalc_H_ #define _TPofBCalc_H_ @@ -23,30 +43,38 @@ class TPofBCalc { public: - TPofBCalc( const string&, const vector& ); - TPofBCalc( const TBofZCalc&, const TTrimSPData&, const vector&, unsigned int ); - TPofBCalc( const TBofZCalcInverse&, const TTrimSPData&, const vector& ); - TPofBCalc( const TBulkVortexFieldCalc&, const vector& ); - TPofBCalc( const vector&, const vector& , double dt = 0.01 ); + // standard constructor: allocates memory for B and PB, the arrays are filled later by one of the Calculate-methods + TPofBCalc(const vector&); + // alternative constructor: PB is not actually calculated but copied from a vector + TPofBCalc(const vector&, const vector& , double dt = 0.01); ~TPofBCalc() { - fB.clear(); - fPB.clear(); + delete[] fB; fB = 0; + delete[] fPB; fPB = 0; } - vector DataB() const {return fB;} - vector DataPB() const {return fPB;} + void Calculate(const string&, const vector&); + void Calculate(const TBofZCalc*, const TTrimSPData*, const vector&, unsigned int); + void Calculate(const TBofZCalcInverse*, const TTrimSPData*, const vector&); + void Calculate(const TBulkVortexFieldCalc*, const vector&); + + const double* DataB() const {return fB;} + double* DataPB() const {return fPB;} double GetBmin() const {return fBmin;} double GetBmax() const {return fBmax;} + unsigned int GetPBSize() const {return fPBSize;} void ConvolveGss(double); void AddBackground(double, double, double); + void UnsetPBExists(); private: - vector fB; - vector fPB; + double *fB; + mutable double *fPB; double fBmin; double fBmax; double fDT; double fDB; + mutable bool fPBExists; + unsigned int fPBSize; }; #endif // _TPofBCalc_H_ diff --git a/src/external/TFitPofB-lib/include/TPofTCalc.h b/src/external/TFitPofB-lib/include/TPofTCalc.h index a81d86b8..a2104ede 100644 --- a/src/external/TFitPofB-lib/include/TPofTCalc.h +++ b/src/external/TFitPofB-lib/include/TPofTCalc.h @@ -9,11 +9,31 @@ ***************************************************************************/ +/*************************************************************************** + * Copyright (C) 2009 by Bastian M. Wojek * + * * + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program; if not, write to the * + * Free Software Foundation, Inc., * + * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. * + ***************************************************************************/ + #ifndef _TPofTCalc_H_ #define _TPofTCalc_H_ #include "TPofBCalc.h" -#include "fftw3.h" +#include #include #define PI 3.14159265358979323846 @@ -23,25 +43,26 @@ class TPofTCalc { public: - TPofTCalc(const string&, const vector&); + TPofTCalc(const TPofBCalc*, const string&, const vector&); ~TPofTCalc(); - vector DataT() const {return fT;} - vector DataPT() const {return fPT;} - void DoFFT(const TPofBCalc&); + const double* DataT() const {return fT;} + const double* DataPT() const {return fPT;} + void DoFFT(); void CalcPol(const vector&); void FakeData(const string&, const vector&); double Eval(double) const; private: fftw_plan fFFTplan; - fftw_complex *fFFTout; double *fFFTin; - vector fT; - vector fPT; + fftw_complex *fFFTout; + double *fT; + double *fPT; double fTBin; int fNFFT; const string fWisdom; + bool fUseWisdom; }; diff --git a/src/external/TFitPofB-lib/include/TSkewedGss.h b/src/external/TFitPofB-lib/include/TSkewedGss.h index b7074ffd..891a2612 100644 --- a/src/external/TFitPofB-lib/include/TSkewedGss.h +++ b/src/external/TFitPofB-lib/include/TSkewedGss.h @@ -9,6 +9,26 @@ ***************************************************************************/ +/*************************************************************************** + * Copyright (C) 2009 by Bastian M. Wojek * + * * + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program; if not, write to the * + * Free Software Foundation, Inc., * + * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. * + ***************************************************************************/ + #ifndef _TSkewedGss_H_ #define _TSkewedGss_H_ @@ -26,6 +46,7 @@ public: private: mutable vector fPar; + TPofBCalc *fPofB; TPofTCalc *fPofT; mutable bool fCalcNeeded; mutable bool fFirstCall; diff --git a/src/external/TFitPofB-lib/include/TSkewedGssLinkDef.h b/src/external/TFitPofB-lib/include/TSkewedGssLinkDef.h index bd87c505..23fe4c47 100644 --- a/src/external/TFitPofB-lib/include/TSkewedGssLinkDef.h +++ b/src/external/TFitPofB-lib/include/TSkewedGssLinkDef.h @@ -9,6 +9,26 @@ ***************************************************************************/ +/*************************************************************************** + * Copyright (C) 2009 by Bastian M. Wojek * + * * + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program; if not, write to the * + * Free Software Foundation, Inc., * + * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. * + ***************************************************************************/ + // root dictionary stuff -------------------------------------------------- #ifdef __CINT__ diff --git a/src/external/TFitPofB-lib/include/TTrimSPDataHandler.h b/src/external/TFitPofB-lib/include/TTrimSPDataHandler.h index 49f3b0b1..2331acd1 100644 --- a/src/external/TFitPofB-lib/include/TTrimSPDataHandler.h +++ b/src/external/TFitPofB-lib/include/TTrimSPDataHandler.h @@ -9,6 +9,26 @@ ***************************************************************************/ +/*************************************************************************** + * Copyright (C) 2009 by Bastian M. Wojek * + * * + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program; if not, write to the * + * Free Software Foundation, Inc., * + * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. * + ***************************************************************************/ + #ifndef _TTrimSPDataHandler_H_ #define _TTrimSPDataHandler_H_ diff --git a/src/external/TFitPofB-lib/include/TVortex.h b/src/external/TFitPofB-lib/include/TVortex.h index 46d7312e..77960fa1 100644 --- a/src/external/TFitPofB-lib/include/TVortex.h +++ b/src/external/TFitPofB-lib/include/TVortex.h @@ -9,6 +9,26 @@ ***************************************************************************/ +/*************************************************************************** + * Copyright (C) 2009 by Bastian M. Wojek * + * * + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program; if not, write to the * + * Free Software Foundation, Inc., * + * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. * + ***************************************************************************/ + #ifndef _TVortex_H_ #define _TVortex_H_ @@ -25,15 +45,16 @@ public: private: mutable vector fPar; + TBulkTriVortexLondonFieldCalc *fVortex; + TPofBCalc *fPofB; TPofTCalc *fPofT; mutable bool fCalcNeeded; mutable bool fFirstCall; - mutable vector fParForPofT; mutable vector fParForVortex; mutable vector fParForPofB; + mutable vector fParForPofT; string fWisdom; unsigned int fGridSteps; - unsigned int fVortexFourierComp; ClassDef(TBulkTriVortexLondon,1) }; diff --git a/src/external/TFitPofB-lib/include/TVortexLinkDef.h b/src/external/TFitPofB-lib/include/TVortexLinkDef.h index 8e9fdfc4..d43b7308 100644 --- a/src/external/TFitPofB-lib/include/TVortexLinkDef.h +++ b/src/external/TFitPofB-lib/include/TVortexLinkDef.h @@ -9,6 +9,26 @@ ***************************************************************************/ +/*************************************************************************** + * Copyright (C) 2009 by Bastian M. Wojek * + * * + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program; if not, write to the * + * Free Software Foundation, Inc., * + * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. * + ***************************************************************************/ + // root dictionary stuff -------------------------------------------------- #ifdef __CINT__ diff --git a/src/external/TFitPofB-lib/test/Makefile.testVortex-v2 b/src/external/TFitPofB-lib/test/Makefile.testVortex-v2 new file mode 100644 index 00000000..68a0042e --- /dev/null +++ b/src/external/TFitPofB-lib/test/Makefile.testVortex-v2 @@ -0,0 +1,49 @@ +#--------------------------------------------------- +# get compilation and library flags from root-config + +ROOTCFLAGS = $(shell $(ROOTSYS)/bin/root-config --cflags) +ROOTLIBS = $(shell $(ROOTSYS)/bin/root-config --libs) + +#--------------------------------------------------- + +CXX = g++-4.4.0 +CXXFLAGS = -O3 -Wall +LOCALINCLUDE = ../include +ROOTINCLUDE = $(ROOTSYS)/include +INCLUDES = -I$(LOCALINCLUDE) -I$(ROOTINCLUDE) +LD = g++-4.4.0 +LDFLAGS = -O3 -L../classes -lTFitPofB -lfftw3_threads -lfftw3 -lm -lpthread -fopenmp -lPMusr + +# the output from the root-config script: +CXXFLAGS += $(ROOTCFLAGS) +LDFLAGS += + +# the ROOT libraries +LIBS = $(ROOTLIBS) -lXMLParser -lMathMore + +EXEC = testVortex-v2 + +# some definitions: headers, sources, objects,... +OBJS = +OBJS += $(EXEC).o + +# make the executable: +# +all: $(EXEC) + +$(EXEC): $(OBJS) + @echo "---> Building $(EXEC) ..." + $(LD) $(LDFLAGS) $(OBJS) -o $(EXEC) $(LIBS) + @echo "done" + +# clean up: remove all object file (and core files) +# semicolon needed to tell make there is no source +# for this target! +# +clean:; @rm -f $(OBJS) + @echo "---> removing $(OBJS)" + +# +$(OBJS): %.o: %.cpp + $(CXX) $(INCLUDES) $(CXXFLAGS) -c $< + diff --git a/src/external/TFitPofB-lib/test/testVortex-v2.cpp b/src/external/TFitPofB-lib/test/testVortex-v2.cpp new file mode 100644 index 00000000..799df076 --- /dev/null +++ b/src/external/TFitPofB-lib/test/testVortex-v2.cpp @@ -0,0 +1,635 @@ +#include "TPofTCalc.h" +#include +#include + +using namespace std; + +#include +#include + +int main(){ + + unsigned int NFFT(256); + + vector parForVortex; + parForVortex.resize(3); + +// parForVortex[0] = 100.0; //app.field +// parForVortex[1] = 200.0; //lambda +// parForVortex[2] = 4.0; //xi + + vector parForPofB; + parForPofB.push_back(0.01); //dt + parForPofB.push_back(0.5); //dB + + vector parForPofT; + parForPofT.push_back(0.0); //phase + parForPofT.push_back(0.01); //dt + parForPofT.push_back(0.5); //dB + + TBulkTriVortexLondonFieldCalc *vortexLattice = new TBulkTriVortexLondonFieldCalc("/home/l_wojek/analysis/WordsOfWisdom.dat", NFFT); + + parForVortex[0] = 500.0; //app.field + parForVortex[1] = 300.0; //lambda + parForVortex[2] = 3.0; //xi + + vortexLattice->SetParameters(parForVortex); + vortexLattice->CalculateGrid(); + +// ofstream ofy("testVortex-B.dat"); +// for (unsigned int j(0); j < NFFT * NFFT; j++) { +// ofy << vortexLattice->DataB()[j] << " "; +// if (!((j+1)%(NFFT))) +// ofy << endl; +// } +// ofy.close(); + + TPofBCalc *PofB = new TPofBCalc(parForPofB); + PofB->Calculate(vortexLattice, parForPofB); + +// const double *b(PofB->DataB()); +// const double *pb(PofB->DataPB()); +// unsigned int s(PofB->GetPBSize()); +// + double test(0.0); +// +// ofstream ofx("testVortex.dat"); +// for (unsigned int i(0); i < s; i++) { +// ofx << b[i] << " " << pb[i] << endl; +// test+=pb[i]; +// } +// ofx.close(); + +// cout << test << endl; + + TPofTCalc poft(PofB, "/home/l_wojek/analysis/WordsOfWisdom.dat", parForPofT); + + poft.DoFFT(); + poft.CalcPol(parForPofT); + + + +// ofstream of8("testVortex-Pt.dat"); + for (double i(0.); i<12.0; i+=0.003) { + test = poft.Eval(i); +// of8 << i << " " << poft.Eval(i) << endl; + } +// of8.close(); + +// parForVortex[0] = 500.0; //app.field +// parForVortex[1] = 100.0; //lambda +// parForVortex[2] = 3.0; //xi +// +// vortexLattice->SetParameters(parForVortex); +// vortexLattice->CalculateGrid(); +// PofB->UnsetPBExists(); +// PofB->Calculate(vortexLattice, parForPofB); +// +// poft.DoFFT(); +// poft.CalcPol(parForPofT); +// +// ofstream of9("testVortex-Pt1.dat"); +// for (double i(0.); i<12.0; i+=0.003) { +// // test = poft.Eval(i); +// of9 << i << " " << poft.Eval(i) << endl; +// } +// of9.close(); + + delete vortexLattice; + vortexLattice = 0; + delete PofB; + PofB = 0; + + parForPofB.clear(); + parForPofT.clear(); + parForVortex.clear(); + +/* + double par_arr[] = {24.4974, E, 94.6, 10.0, 130.0}; + vector par_vec(par_arr, par_arr+(sizeof(par_arr)/sizeof(par_arr[0]))); + + + vector parForBofZ; + for (unsigned int i(2); i parForPofB; + parForPofB.push_back(0.01); //dt + parForPofB.push_back(0.1); //dB + parForPofB.push_back(par_vec[1]); + parForPofB.push_back(par_vec[2]); // Bkg-Field + parForPofB.push_back(0.005); // Bkg-width (in principle zero) + vector interfaces; + interfaces.push_back(par_vec[3]); // dead layer + parForPofB.push_back(calcData.LayerFraction(par_vec[1], 1, interfaces)); // Fraction of muons in the deadlayer + + + vector zonk; + vector donk; + vector hurgaB; + vector hurgaPB; + + for (unsigned int u(0); u<10 ; u++) { + + parForBofZ[par_vec.size()-3] = par_vec[4] + double(u)*5.; + + TLondon1D_HS *BofZ = new TLondon1D_HS(parForBofZ); + TPofBCalc *PofB = new TPofBCalc(*BofZ, calcData, parForPofB); + + if(u == 0){ + hurgaB = PofB->DataB(); + hurgaPB = PofB->DataPB(); + PofB->ConvolveGss(7.0); + donk = PofB->DataPB(); + } + + delete BofZ; + delete PofB; + } + +// TPofBCalc sumPB(hurgaB, zonk); + +// sumPB.AddBackground(par_vec[2], 0.15, 0.1); +// sumPB.ConvolveGss(7.0); +// donk = sumPB.DataPB(); + + char debugfile[50]; + int n = sprintf (debugfile, "testInverseExp-BpB_B-%.4f_l-%.3f_E-%.1f_normal.dat", par_vec[2], par_vec[4], par_vec[1]); + + if (n > 0) { + ofstream of7(debugfile); + for (unsigned int i(0); i parForPofT; +// parForPofT.push_back(par_vec[0]); //phase +// parForPofT.push_back(0.01); //dt +// parForPofT.push_back(0.02); //dB + +// TLondon1D_HS BofZ(parForBofZ); + +// vector zz; +// vector Bz; +// for(double i(0.5); i<2001; i+=1.0) { +// zz.push_back(i/10.0); +// if(!(i/10.0 < parForBofZ[1])) { +// Bz.push_back(BofZ.GetBofZ(i/10.0)); +// } +// else { +// Bz.push_back(parForBofZ[0]); +// } +// } +/* + char debugfile1[50]; + int nn = sprintf (debugfile1, "testInverseExp-Bz_z-%.4f_l-%.3f_E-%.1f_normal.dat", par_vec[2], par_vec[9], par_vec[1]); + if (nn > 0) { + ofstream of01(debugfile1); + for (unsigned int i(0); i<2000; i++) { + of01 << zz[i] << " " << Bz[i] << endl; + } + of01.close(); + } + */ +// TPofBCalc PofB(BofZ, calcData, parForPofB); + + +// double t1(0.0); +// get start time +// struct timeval tv_start, tv_stop; +// gettimeofday(&tv_start, 0); +// +// TPofBCalc PofB(BofZ, calcData, parForPofB); +// +// gettimeofday(&tv_stop, 0); +// t1 = (tv_stop.tv_sec - tv_start.tv_sec)*1000.0 + (tv_stop.tv_usec - tv_start.tv_usec)/1000.0; + +// cout << "p(B) calculation time with derivatives of the inverse function: " << t1 << " ms" << endl; + +// gettimeofday(&tv_start, 0); +// +// BofZ.Calculate(); +// TPofBCalc PofB1(BofZ, calcData, parForPofB, 1); +// +// gettimeofday(&tv_stop, 0); +// t1 = (tv_stop.tv_sec - tv_start.tv_sec)*1000.0 + (tv_stop.tv_usec - tv_start.tv_usec)/1000.0; +// +// cout << "p(B) calculation time without derivatives of the inverse function: " << t1 << " ms" << endl; + +// PofB.ConvolveGss(1.17); + +// PofB.AddBackground(par_vec[2], 0.2, calcData.LayerFraction(E, 4, interfaces)); +// PofB.ConvolveGss(1.17); + +// vector hurgaB(PofB.DataB()); +// vector hurgaPB(PofB.DataPB()); +// vector hurgaPB1(PofB1.DataPB()); +// +// char debugfile[50]; +// int n = sprintf (debugfile, "testInverseExp-BpB_B-%.4f_l-%.3f_E-%.1f_normal.dat", par_vec[2], par_vec[4], par_vec[1]); +// +// if (n > 0) { +// ofstream of7(debugfile); +// for (unsigned int i(0); i 0) { +// ofstream of8(debugfilex); +// for (double i(0.); i<12.0; i+=0.003) { +// of8 << i << " " << poft.Eval(i) << endl; +// } +// of8.close(); +// } +/* + PofB.ConvolveGss(8.8); + + vector hurgaB1(PofB.DataB()); + vector hurgaPB1(PofB.DataPB()); + + n = sprintf (debugfile, "BpB_B-%.4f_l-%.3f_E-%.1f_broadend8.8G.dat", par_vec[2], par_vec[4], par_vec[1]); + + if (n > 0) { + ofstream of1(debugfile); + for (unsigned int i(0); i z2(calcData.DataZ(par_vec[1])); + vector nz2(calcData.DataNZ(par_vec[1])); + + ofstream of2("Implantation-profile-broad.dat"); + for (unsigned int i(0); i hurgaB2(PofB2.DataB()); + vector hurgaPB2(PofB2.DataPB()); + + n = sprintf (debugfile, "BpB_B-%.4f_l-%.3f_E-%.1f_broadend10nm.dat", par_vec[2], par_vec[4], par_vec[1]); + + if (n>0) { + ofstream of8(debugfile); + for (unsigned int i(0); i hurgaB3(PofB2.DataB()); + vector hurgaPB3(PofB2.DataPB()); + + n = sprintf (debugfile, "BpB_B-%.4f_l-%.3f_E-%.1f_broadend10nm+7.5G.dat", par_vec[2], par_vec[4], par_vec[1]); + + if (n > 0) { + ofstream of9(debugfile); + for (unsigned int i(0); i parameter(param,param+8); + vector param_for_BofZ; + vector param_for_PofB; + vector param_for_PofT; + + for (unsigned int i(0); i<4; i++) + param_for_BofZ.push_back(parameter[i]); + + for (unsigned int i(5); i<8; i++) + param_for_PofB.push_back(parameter[i]); + + for (unsigned int i(4); i<7; i++) + param_for_PofT.push_back(parameter[i]); + + TLondon1D_1L bofz(param_for_BofZ); + + cout << "Bmin = " << bofz.GetBmin() << endl; + cout << "Bmax = " << bofz.GetBmax() << endl; + cout << "dZ = " << bofz.GetdZ() << endl; + + ofstream of5("test_Bz.dat"); + for (double i(0); i<50.; i+=0.1) { + of5 << i << " " << bofz.GetBofZ(i) << endl; + } + of5.close(); + + ofstream of6("test_zBz.dat"); + for (unsigned int i(0); i hurgaB(pofb.DataB()); + vector hurgaPB(pofb.DataPB()); + + ofstream of7("test_BpB.dat"); + for (unsigned int i(0); i parNo_vec(parNo_arr, parNo_arr+(sizeof(parNo_arr)/sizeof(parNo_arr[0]))); + vector par_vec(par_arr, par_arr+(sizeof(par_arr)/sizeof(par_arr[0]))); + +// vector par_vec_sub; + +// for(unsigned int i(0); i parNo_vec(parNo_arr, parNo_arr+(sizeof(parNo_arr)/sizeof(parNo_arr[0]))); + vector par_vec(par_arr, par_arr+(sizeof(par_arr)/sizeof(par_arr[0]))); + + vector par_vec_sub; + + for(unsigned int i(0); i parNo_vec(parNo_arr, parNo_arr+(sizeof(parNo_arr)/sizeof(parNo_arr[0]))); + vector par_vec(par_arr, par_arr+(sizeof(par_arr)/sizeof(par_arr[0]))); + + vector par_vec_sub; + + for(unsigned int i(0); i parNo_vec(parNo_arr, parNo_arr+(sizeof(parNo_arr)/sizeof(parNo_arr[0]))); + vector par_vec(par_arr, par_arr+(sizeof(par_arr)/sizeof(par_arr[0]))); + + vector par_vec_sub; + + for(unsigned int i(0); i parNo_vec(parNo_arr, parNo_arr+(sizeof(parNo_arr)/sizeof(parNo_arr[0]))); + vector par_vec(par_arr, par_arr+(sizeof(par_arr)/sizeof(par_arr[0]))); + + vector par_vec_sub; + + for(unsigned int i(0); i data01, data02, data03, data04, data05, data06, data07, data08, data09, data10; + + for (double i(0.); i<12.0; i+=0.003) + data01.push_back(fitter.Eval(i, par_vec_sub)); + + + par_vec_sub[0] += 10.0; + par_vec_sub[10] -= 10.0; + + for (double i(0.); i<12.0; i+=0.003) + data02.push_back(fitter.Eval(i, par_vec_sub)); + + + par_vec_sub[0] += 10.0; + par_vec_sub[10] -= 10.0; + + for (double i(0.); i<12.0; i+=0.003) + data03.push_back(fitter.Eval(i, par_vec_sub)); + + par_vec_sub[0] += 10.0; + + for (double i(0.); i<12.0; i+=0.003) + data04.push_back(fitter.Eval(i, par_vec_sub)); + + + + par_vec_sub[0] += 10.0; + par_vec_sub[10] -= 10.0; + + for (double i(0.); i<12.0; i+=0.003) + data05.push_back(fitter.Eval(i, par_vec_sub)); + + par_vec_sub[0] += 10.0; + par_vec_sub[10] -= 10.0; + + for (double i(0.); i<12.0; i+=0.003) + data06.push_back(fitter.Eval(i, par_vec_sub)); + + par_vec_sub[0] += 10.0; + + for (double i(0.); i<12.0; i+=0.003) + data07.push_back(fitter.Eval(i, par_vec_sub)); + + par_vec_sub[0] += 10.0; + par_vec_sub[10] = 190.0; + + for (double i(0.); i<12.0; i+=0.003) + data08.push_back(fitter.Eval(i, par_vec_sub)); + + par_vec_sub[0] += 10.0; + par_vec_sub[10] -= 10.0; + + for (double i(0.); i<12.0; i+=0.003) + data09.push_back(fitter.Eval(i, par_vec_sub)); + + par_vec_sub[0] += 10.0; + + for (double i(0.); i<12.0; i+=0.003) + data10.push_back(fitter.Eval(i, par_vec_sub)); + + */ + + +// par_vec.clear(); + + + + + + return 0; +}