Fixed a Bessel-function-bug in the class calculating LF relaxation functions through Laplace transforms
This commit is contained in:
Bastian M. Wojek
parent
c3aebf04d6
commit
ccaba04d51
@ -24,10 +24,12 @@ 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
|
||||
|
||||
@ -56,6 +58,10 @@ $(OBJS): %.o: %.cpp
|
||||
|
||||
# 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) $^
|
||||
|
||||
101
src/external/TFitPofB-lib/classes/TBulkVortexFieldCalc.cpp
vendored
Normal file
101
src/external/TFitPofB-lib/classes/TBulkVortexFieldCalc.cpp
vendored
Normal file
@ -0,0 +1,101 @@
|
||||
/***************************************************************************
|
||||
|
||||
TBulkVortexFieldCalc.cpp
|
||||
|
||||
Author: Bastian M. Wojek, Alexander Maisuradze
|
||||
e-mail: bastian.wojek@psi.ch
|
||||
|
||||
2009/10/17
|
||||
|
||||
***************************************************************************/
|
||||
|
||||
#include "TBulkVortexFieldCalc.h"
|
||||
|
||||
#include <cmath>
|
||||
#include <omp.h>
|
||||
|
||||
#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
|
||||
|
||||
double getXi(const double hc2) { // get xi given Hc2 in Gauss
|
||||
if (hc2)
|
||||
return sqrt(fluxQuantum/(TWOPI*hc2));
|
||||
else
|
||||
return 0.0;
|
||||
}
|
||||
|
||||
double getHc2(const double xi) { // get Hc2 given xi in nm
|
||||
if (xi)
|
||||
return fluxQuantum/(TWOPI*xi*xi);
|
||||
else
|
||||
return 0.0;
|
||||
}
|
||||
|
||||
TBulkTriVortexLondonFieldCalc::TBulkTriVortexLondonFieldCalc(const vector<double>& 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<int>(fNFourierComp); mm++) {
|
||||
for(int nn = -fNFourierComp; nn <= static_cast<int>(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)));
|
||||
}
|
||||
}
|
||||
return fParam[2]*fourierSum;
|
||||
}
|
||||
else
|
||||
return fParam[2];
|
||||
|
||||
}
|
||||
|
||||
void TBulkTriVortexLondonFieldCalc::CalculateGrid() const {
|
||||
|
||||
double bb(cos(PI/6.0)*0.5);
|
||||
double yStep(bb/static_cast<double>(fSteps));
|
||||
double xStep(0.5/static_cast<double>(fSteps));
|
||||
|
||||
fB.resize(fSteps);
|
||||
for (unsigned int i(0); i < fSteps; i++)
|
||||
fB[i].resize(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<int>(fSteps); j++) {
|
||||
fB[i][j] = GetBatPoint(static_cast<double>(i)*xStep, static_cast<double>(j)*yStep);
|
||||
}
|
||||
// end pragma omp parallel
|
||||
}
|
||||
|
||||
return;
|
||||
|
||||
}
|
||||
|
||||
double TBulkTriVortexLondonFieldCalc::GetBmin() const {
|
||||
return GetBatPoint(0.5, 0.5*tan(PI/6.0));
|
||||
}
|
||||
|
||||
double TBulkTriVortexLondonFieldCalc::GetBmax() const {
|
||||
if (!fB.empty())
|
||||
return fB[0][0];
|
||||
else
|
||||
return GetBatPoint(0.0, 0.0);
|
||||
}
|
||||
|
||||
@ -107,8 +107,12 @@ void TFitPofBStartupHandler::OnStartElement(const char *str, const TList *attrib
|
||||
fKey = eDeltaB;
|
||||
} else if (!strcmp(str, "wisdom")) {
|
||||
fKey = eWisdomFile;
|
||||
} if (!strcmp(str, "N_theory")) {
|
||||
} else if (!strcmp(str, "N_theory")) {
|
||||
fKey = eNSteps;
|
||||
} else if (!strcmp(str, "N_VortexGrid")) {
|
||||
fKey = eGridSteps;
|
||||
} else if (!strcmp(str, "N_VortexFourier")) {
|
||||
fKey = eVortexFourierComp;
|
||||
}
|
||||
}
|
||||
|
||||
@ -157,9 +161,17 @@ void TFitPofBStartupHandler::OnCharacters(const char *str)
|
||||
fWisdomFile = str;
|
||||
break;
|
||||
case eNSteps:
|
||||
// convert str to int and assign it to the deltat-member
|
||||
// convert str to int and assign it to the NSteps-member
|
||||
fNSteps = atoi(str);
|
||||
break;
|
||||
case eGridSteps:
|
||||
// convert str to int and assign it to the GridSteps-member
|
||||
fGridSteps = atoi(str);
|
||||
break;
|
||||
case eVortexFourierComp:
|
||||
// convert str to int and assign it to the VortexFourierComp-member
|
||||
fVortexFourierComp = atoi(str);
|
||||
break;
|
||||
default:
|
||||
break;
|
||||
}
|
||||
@ -292,6 +304,20 @@ void TFitPofBStartupHandler::CheckLists()
|
||||
fNSteps = 3000;
|
||||
}
|
||||
|
||||
// check if any number of steps for the theory function is specified
|
||||
cout << endl << "TFitPofBStartupHandler::CheckLists: check number of steps for Vortex grid ..." << endl;
|
||||
if (!fGridSteps) {
|
||||
cout << endl << "TFitPofBStartupHandler::CheckLists: You did not specify the number of steps for the grid. Setting the default." << endl;
|
||||
fGridSteps = 200;
|
||||
}
|
||||
|
||||
// check if any number of steps for the theory function is specified
|
||||
cout << endl << "TFitPofBStartupHandler::CheckLists: check number of steps for Vortex lattice Fourier components ..." << endl;
|
||||
if (!fVortexFourierComp) {
|
||||
cout << endl << "TFitPofBStartupHandler::CheckLists: You did not specify the number of Vortex lattice Fourier components. Setting the default." << endl;
|
||||
fVortexFourierComp = 17;
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
// end ---------------------------------------------------------------------
|
||||
|
||||
77
src/external/TFitPofB-lib/classes/TPofBCalc.cpp
vendored
77
src/external/TFitPofB-lib/classes/TPofBCalc.cpp
vendored
@ -16,6 +16,8 @@
|
||||
#include <fstream>
|
||||
#include <cassert>
|
||||
|
||||
#include <omp.h>
|
||||
|
||||
/* USED FOR DEBUGGING-----------------------------------
|
||||
#include <cstdio>
|
||||
#include <ctime>
|
||||
@ -368,6 +370,81 @@ TPofBCalc::TPofBCalc( const TBofZCalc &BofZ, const TTrimSPData &dataTrimSP, cons
|
||||
|
||||
}
|
||||
|
||||
//-----------
|
||||
// Constructor that does the P(B) calculation for a bulk vortex lattice
|
||||
// Parameters: dt[us], dB[G] [, Bbg[G], width[us^{-1}], weight[1] ]
|
||||
//-----------
|
||||
|
||||
TPofBCalc::TPofBCalc( const TBulkVortexFieldCalc &vortexLattice, const vector<double> ¶ ) : fDT(para[0]), fDB(para[1])
|
||||
{
|
||||
|
||||
fBmin = vortexLattice.GetBmin();
|
||||
fBmax = vortexLattice.GetBmax();
|
||||
|
||||
unsigned int a1(static_cast<int>(floor(fBmin/fDB)));
|
||||
unsigned int a2(static_cast<int>(ceil(fBmin/fDB)));
|
||||
unsigned int a3(static_cast<int>(floor(fBmax/fDB)));
|
||||
unsigned int a4(static_cast<int>(ceil(fBmax/fDB)));
|
||||
|
||||
unsigned int firstZerosEnd ((a1 != a2) ? a1 : (a1 - 1));
|
||||
unsigned int lastZerosStart ((a3 != a4) ? a4 : (a4 + 1));
|
||||
|
||||
// fill the B vector
|
||||
unsigned int indexBmaxFFT(static_cast<int>(ceil(1.0/(gBar*fDT*fDB))));
|
||||
fB.resize(indexBmaxFFT, 0.0);
|
||||
fPB.resize(indexBmaxFFT, 0.0);
|
||||
|
||||
int i;
|
||||
#pragma omp parallel for default(shared) private(i) schedule(dynamic)
|
||||
for (i = 0; i < static_cast<int>(indexBmaxFFT); i++)
|
||||
fB[i] = fDB*static_cast<double>(i);
|
||||
// end pragma omp parallel
|
||||
|
||||
// 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);
|
||||
}
|
||||
|
||||
vector< vector<double> > vortexFields(vortexLattice.DataB());
|
||||
|
||||
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<int>(ceil(vortexFields[j][k]/fDB))] += 1.0;
|
||||
}
|
||||
}
|
||||
|
||||
double normalizer(static_cast<double>(numberOfPoints*numberOfPoints)*fDB);
|
||||
|
||||
#pragma omp parallel for default(shared) private(i) schedule(dynamic)
|
||||
for (i = firstZerosEnd; i <= static_cast<int>(lastZerosStart); i++) {
|
||||
fPB[i] /= normalizer;
|
||||
}
|
||||
// end pragma omp parallel
|
||||
|
||||
// // normalize p(B)
|
||||
//
|
||||
// double pBsum = 0.0;
|
||||
// for (unsigned int i(firstZerosEnd); i<lastZerosStart; i++)
|
||||
// pBsum += fPB[i];
|
||||
// pBsum *= fDB;
|
||||
// for (unsigned int i(firstZerosEnd); i<lastZerosStart; i++)
|
||||
// fPB[i] /= pBsum;
|
||||
|
||||
if(para.size() == 5)
|
||||
AddBackground(para[2], para[3], para[4]);
|
||||
|
||||
}
|
||||
|
||||
|
||||
|
||||
// TPofBCalc::AddBackground, Parameters: field B[G], width s[us], weight w
|
||||
|
||||
void TPofBCalc::AddBackground(double B, double s, double w) {
|
||||
|
||||
163
src/external/TFitPofB-lib/classes/TVortex.cpp
vendored
Normal file
163
src/external/TFitPofB-lib/classes/TVortex.cpp
vendored
Normal file
@ -0,0 +1,163 @@
|
||||
/***************************************************************************
|
||||
|
||||
TVortex.cpp
|
||||
|
||||
Author: Bastian M. Wojek
|
||||
e-mail: bastian.wojek@psi.ch
|
||||
|
||||
2009/10/17
|
||||
|
||||
***************************************************************************/
|
||||
|
||||
#include "TVortex.h"
|
||||
#include <iostream>
|
||||
#include <cassert>
|
||||
#include <cmath>
|
||||
using namespace std;
|
||||
|
||||
#include <TSAXParser.h>
|
||||
#include "TFitPofBStartupHandler.h"
|
||||
|
||||
ClassImp(TBulkTriVortexLondon)
|
||||
|
||||
//------------------
|
||||
// Destructor of the TBulkTriVortexLondon class -- cleaning up
|
||||
//------------------
|
||||
|
||||
TBulkTriVortexLondon::~TBulkTriVortexLondon() {
|
||||
fPar.clear();
|
||||
fParForVortex.clear();
|
||||
fParForPofB.clear();
|
||||
fParForPofT.clear();
|
||||
delete fPofT;
|
||||
fPofT = 0;
|
||||
}
|
||||
|
||||
//------------------
|
||||
// Constructor of the TBulkTriVortexLondon class
|
||||
// creates (a pointer to) the TPofTCalc object (with the FFT plan)
|
||||
//------------------
|
||||
|
||||
TBulkTriVortexLondon::TBulkTriVortexLondon() : fCalcNeeded(true), fFirstCall(true) {
|
||||
|
||||
// read startup file
|
||||
string startup_path_name("TFitPofB_startup.xml");
|
||||
|
||||
TSAXParser *saxParser = new TSAXParser();
|
||||
TFitPofBStartupHandler *startupHandler = new TFitPofBStartupHandler();
|
||||
saxParser->ConnectToHandler("TFitPofBStartupHandler", startupHandler);
|
||||
int status (saxParser->ParseFile(startup_path_name.c_str()));
|
||||
// check for parse errors
|
||||
if (status) { // error
|
||||
cout << endl << "**WARNING** reading/parsing TFitPofB_startup.xml failed." << endl;
|
||||
}
|
||||
|
||||
fGridSteps = startupHandler->GetGridSteps();
|
||||
fWisdom = startupHandler->GetWisdomFile();
|
||||
fVortexFourierComp = startupHandler->GetVortexFourierComp();
|
||||
|
||||
fParForVortex.clear();
|
||||
|
||||
fParForPofT.push_back(0.0);
|
||||
fParForPofT.push_back(startupHandler->GetDeltat());
|
||||
fParForPofT.push_back(startupHandler->GetDeltaB());
|
||||
|
||||
fParForPofB.push_back(startupHandler->GetDeltat());
|
||||
fParForPofB.push_back(startupHandler->GetDeltaB());
|
||||
|
||||
fParForPofB.push_back(0.0); // Bkg-Field
|
||||
fParForPofB.push_back(0.005); // Bkg-width
|
||||
fParForPofB.push_back(0.0); // Bkg-weight
|
||||
|
||||
TPofTCalc *y = new TPofTCalc(fWisdom, fParForPofT);
|
||||
fPofT = y;
|
||||
y = 0;
|
||||
|
||||
// clean up
|
||||
if (saxParser) {
|
||||
delete saxParser;
|
||||
saxParser = 0;
|
||||
}
|
||||
if (startupHandler) {
|
||||
delete startupHandler;
|
||||
startupHandler = 0;
|
||||
}
|
||||
}
|
||||
|
||||
//------------------
|
||||
// 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)
|
||||
//------------------
|
||||
|
||||
double TBulkTriVortexLondon::operator()(double t, const vector<double> &par) const {
|
||||
|
||||
assert(par.size() == 4 || par.size() == 5);
|
||||
|
||||
if(t<0.0)
|
||||
return cos(par[0]*0.017453293);
|
||||
|
||||
// check if the function is called the first time and if yes, read in parameters
|
||||
|
||||
if(fFirstCall){
|
||||
fPar = par;
|
||||
|
||||
for (unsigned int i(1); i < 4; i++){
|
||||
fParForVortex.push_back(fPar[i]);
|
||||
}
|
||||
fFirstCall = false;
|
||||
}
|
||||
|
||||
// check if any parameter has changed
|
||||
|
||||
bool par_changed(false);
|
||||
bool only_phase_changed(false);
|
||||
|
||||
for (unsigned int i(0); i<fPar.size(); i++) {
|
||||
if( fPar[i]-par[i] ) {
|
||||
fPar[i] = par[i];
|
||||
par_changed = true;
|
||||
if (i == 0) {
|
||||
only_phase_changed = true;
|
||||
} else {
|
||||
only_phase_changed = false;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if (par_changed)
|
||||
fCalcNeeded = true;
|
||||
|
||||
// if model parameters have changed, recalculate B(x,y), P(B) and P(t)
|
||||
|
||||
if (fCalcNeeded) {
|
||||
|
||||
fParForPofT[0] = par[0]; // phase
|
||||
|
||||
if(!only_phase_changed) {
|
||||
|
||||
// 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
|
||||
|
||||
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);
|
||||
|
||||
}/* else {
|
||||
cout << "Only the phase parameter has changed, (re-)calculating P(t) now..." << endl;
|
||||
}*/
|
||||
|
||||
fPofT->CalcPol(fParForPofT);
|
||||
|
||||
fCalcNeeded = false;
|
||||
}
|
||||
|
||||
return fPofT->Eval(t);
|
||||
|
||||
}
|
||||
67
src/external/TFitPofB-lib/include/TBulkVortexFieldCalc.h
vendored
Normal file
67
src/external/TFitPofB-lib/include/TBulkVortexFieldCalc.h
vendored
Normal file
@ -0,0 +1,67 @@
|
||||
/***************************************************************************
|
||||
|
||||
TBulkVortexFieldCalc.h
|
||||
|
||||
Author: Bastian M. Wojek, Alexander Maisuradze
|
||||
e-mail: bastian.wojek@psi.ch
|
||||
|
||||
2009/10/17
|
||||
|
||||
***************************************************************************/
|
||||
|
||||
#ifndef _TBulkVortexFieldCalc_H_
|
||||
#define _TBulkVortexFieldCalc_H_
|
||||
|
||||
#include <vector>
|
||||
using namespace std;
|
||||
|
||||
//--------------------
|
||||
// Base class for any kind of vortex symmetry
|
||||
//--------------------
|
||||
|
||||
class TBulkVortexFieldCalc {
|
||||
|
||||
public:
|
||||
|
||||
TBulkVortexFieldCalc() {}
|
||||
|
||||
virtual ~TBulkVortexFieldCalc() {
|
||||
for(unsigned int i(0); i < fB.size(); i++)
|
||||
fB[i].clear();
|
||||
fB.clear();
|
||||
fParam.clear();
|
||||
}
|
||||
|
||||
virtual vector< vector<double> > DataB() const {return fB;}
|
||||
virtual void CalculateGrid() const = 0;
|
||||
virtual double GetBatPoint(double, double) const = 0;
|
||||
virtual double GetBmin() const = 0;
|
||||
virtual double GetBmax() const = 0;
|
||||
|
||||
protected:
|
||||
unsigned int fSteps; // number of steps, the "unit cell" of the vortex lattice is devided in (in x- and y- direction)
|
||||
vector<double> fParam; // lambda, xi, appl.field
|
||||
mutable vector< vector<double> > fB; // fields in a "unit cell" of the vortex lattice
|
||||
};
|
||||
|
||||
//--------------------
|
||||
// Class for triangular vortex lattice, London model
|
||||
//--------------------
|
||||
|
||||
class TBulkTriVortexLondonFieldCalc : public TBulkVortexFieldCalc {
|
||||
|
||||
public:
|
||||
|
||||
TBulkTriVortexLondonFieldCalc(const vector<double>&, const unsigned int steps = 200, const unsigned int Ncomp = 17);
|
||||
~TBulkTriVortexLondonFieldCalc() {}
|
||||
|
||||
void CalculateGrid() const;
|
||||
double GetBatPoint(double, double) const;
|
||||
double GetBmin() const;
|
||||
double GetBmax() const;
|
||||
|
||||
private:
|
||||
unsigned int fNFourierComp;
|
||||
};
|
||||
|
||||
#endif // _TBulkVortexFieldCalc_H_
|
||||
@ -65,9 +65,11 @@ class TFitPofBStartupHandler : public TQObject {
|
||||
virtual const double GetDeltaB() const { return fDeltaB; }
|
||||
virtual const string GetWisdomFile() const { return fWisdomFile; }
|
||||
virtual const unsigned int GetNSteps() const { return fNSteps; }
|
||||
virtual const unsigned int GetGridSteps() const { return fGridSteps; }
|
||||
virtual const unsigned int GetVortexFourierComp() const { return fVortexFourierComp; }
|
||||
|
||||
private:
|
||||
enum EKeyWords {eEmpty, eComment, eDataPath, eEnergy, eEnergyList, eDeltat, eDeltaB, eWisdomFile, eNSteps};
|
||||
enum EKeyWords {eEmpty, eComment, eDataPath, eEnergy, eEnergyList, eDeltat, eDeltaB, eWisdomFile, eNSteps, eGridSteps, eVortexFourierComp};
|
||||
|
||||
EKeyWords fKey;
|
||||
|
||||
@ -77,6 +79,8 @@ class TFitPofBStartupHandler : public TQObject {
|
||||
double fDeltaB;
|
||||
string fWisdomFile;
|
||||
unsigned int fNSteps;
|
||||
unsigned int fGridSteps;
|
||||
unsigned int fVortexFourierComp;
|
||||
|
||||
ClassDef(TFitPofBStartupHandler, 1)
|
||||
};
|
||||
|
||||
@ -14,6 +14,7 @@
|
||||
|
||||
#include "TBofZCalc.h"
|
||||
#include "TTrimSPDataHandler.h"
|
||||
#include "TBulkVortexFieldCalc.h"
|
||||
|
||||
#define gBar 0.0135538817
|
||||
#define pi 3.14159265358979323846
|
||||
@ -25,6 +26,7 @@ public:
|
||||
TPofBCalc( const string&, const vector<double>& );
|
||||
TPofBCalc( const TBofZCalc&, const TTrimSPData&, const vector<double>&, unsigned int );
|
||||
TPofBCalc( const TBofZCalcInverse&, const TTrimSPData&, const vector<double>& );
|
||||
TPofBCalc( const TBulkVortexFieldCalc&, const vector<double>& );
|
||||
TPofBCalc( const vector<double>&, const vector<double>& , double dt = 0.01 );
|
||||
~TPofBCalc() {
|
||||
fB.clear();
|
||||
|
||||
41
src/external/TFitPofB-lib/include/TVortex.h
vendored
Normal file
41
src/external/TFitPofB-lib/include/TVortex.h
vendored
Normal file
@ -0,0 +1,41 @@
|
||||
/***************************************************************************
|
||||
|
||||
TVortex.h
|
||||
|
||||
Author: Bastian M. Wojek
|
||||
e-mail: bastian.wojek@psi.ch
|
||||
|
||||
2009/10/17
|
||||
|
||||
***************************************************************************/
|
||||
|
||||
#ifndef _TVortex_H_
|
||||
#define _TVortex_H_
|
||||
|
||||
#include "PUserFcnBase.h"
|
||||
#include "TPofTCalc.h"
|
||||
|
||||
class TBulkTriVortexLondon : public PUserFcnBase {
|
||||
|
||||
public:
|
||||
TBulkTriVortexLondon();
|
||||
~TBulkTriVortexLondon();
|
||||
|
||||
double operator()(double, const vector<double>&) const;
|
||||
|
||||
private:
|
||||
mutable vector<double> fPar;
|
||||
TPofTCalc *fPofT;
|
||||
mutable bool fCalcNeeded;
|
||||
mutable bool fFirstCall;
|
||||
mutable vector<double> fParForPofT;
|
||||
mutable vector<double> fParForVortex;
|
||||
mutable vector<double> fParForPofB;
|
||||
string fWisdom;
|
||||
unsigned int fGridSteps;
|
||||
unsigned int fVortexFourierComp;
|
||||
|
||||
ClassDef(TBulkTriVortexLondon,1)
|
||||
};
|
||||
|
||||
#endif //_TLondon1D_H_
|
||||
23
src/external/TFitPofB-lib/include/TVortexLinkDef.h
vendored
Normal file
23
src/external/TFitPofB-lib/include/TVortexLinkDef.h
vendored
Normal file
@ -0,0 +1,23 @@
|
||||
/***************************************************************************
|
||||
|
||||
TVortexLinkDef.h
|
||||
|
||||
Author: Bastian M. Wojek
|
||||
e-mail: bastian.wojek@psi.ch
|
||||
|
||||
2009/10/17
|
||||
|
||||
***************************************************************************/
|
||||
|
||||
// root dictionary stuff --------------------------------------------------
|
||||
#ifdef __CINT__
|
||||
|
||||
#pragma link off all globals;
|
||||
#pragma link off all classes;
|
||||
#pragma link off all functions;
|
||||
|
||||
#pragma link C++ class TBulkTriVortexLondon+;
|
||||
|
||||
#endif //__CINT__
|
||||
// root dictionary stuff --------------------------------------------------
|
||||
|
||||
49
src/external/TFitPofB-lib/test/Makefile.testVortex
vendored
Normal file
49
src/external/TFitPofB-lib/test/Makefile.testVortex
vendored
Normal file
@ -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 -lm -fopenmp -lPMusr
|
||||
|
||||
# the output from the root-config script:
|
||||
CXXFLAGS += $(ROOTCFLAGS)
|
||||
LDFLAGS +=
|
||||
|
||||
# the ROOT libraries
|
||||
LIBS = $(ROOTLIBS) -lXMLParser -lMathMore
|
||||
|
||||
EXEC = testVortex
|
||||
|
||||
# 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 $<
|
||||
|
||||
@ -5,12 +5,16 @@
|
||||
Defines path/prefix and energies (keV, format: %02u_%1u) of TrimSP-rge-files,
|
||||
path/name to the FFTW-wisdom-file and time/field binning (us/G)
|
||||
N_theory determines the number of points in "real space" where the theory function will be calculated
|
||||
N_VortexGrid determines the number of points in which each axis of the "real space" will be divided when calculating the vortex B(x,y)
|
||||
N_VortexFourier defines the number of Fourier components which should be summed up in order to get B(x,y)
|
||||
</comment>
|
||||
<data_path>/home/l_wojek/TrimSP/YBCOxtal/YBCOxtal-500000-</data_path>
|
||||
<wisdom>/home/l_wojek/analysis/WordsOfWisdom.dat</wisdom>
|
||||
<delta_t>0.01</delta_t>
|
||||
<delta_B>0.1</delta_B>
|
||||
<N_theory>5000</N_theory>
|
||||
<N_VortexGrid>200</N_VortexGrid>
|
||||
<N_VortexFourier>17</N_VortexFourier>
|
||||
<energy_list>
|
||||
<energy>03_0</energy>
|
||||
<energy>03_6</energy>
|
||||
|
||||
565
src/external/TFitPofB-lib/test/testVortex.cpp
vendored
Normal file
565
src/external/TFitPofB-lib/test/testVortex.cpp
vendored
Normal file
@ -0,0 +1,565 @@
|
||||
#include "TPofTCalc.h"
|
||||
#include <iostream>
|
||||
#include <fstream>
|
||||
|
||||
using namespace std;
|
||||
|
||||
#include <cstdio>
|
||||
#include <sys/time.h>
|
||||
|
||||
int main(){
|
||||
|
||||
vector<double> parForVortex;
|
||||
parForVortex.push_back(150.0); //lambda
|
||||
parForVortex.push_back(3.0); //xi
|
||||
parForVortex.push_back(500.0); //app.field
|
||||
|
||||
vector<double> parForPofB;
|
||||
parForPofB.push_back(0.01); //dt
|
||||
parForPofB.push_back(1.0); //dB
|
||||
|
||||
TBulkTriVortexLondonFieldCalc *vortexLattice = new TBulkTriVortexLondonFieldCalc(parForVortex, 250, 17);
|
||||
TPofBCalc *PofB = new TPofBCalc(*vortexLattice, parForPofB);
|
||||
|
||||
vector<double> b(PofB->DataB());
|
||||
vector<double> pb(PofB->DataPB());
|
||||
|
||||
double test(0.0);
|
||||
|
||||
ofstream ofx("testVector.dat");
|
||||
for (unsigned int i(0); i<b.size(); i++) {
|
||||
ofx << b[i] << " " << pb[i] << endl;
|
||||
test+=pb[i];
|
||||
}
|
||||
ofx.close();
|
||||
|
||||
cout << test << endl;
|
||||
/*
|
||||
double par_arr[] = {24.4974, E, 94.6, 10.0, 130.0};
|
||||
vector<double> par_vec(par_arr, par_arr+(sizeof(par_arr)/sizeof(par_arr[0])));
|
||||
|
||||
|
||||
vector<double> parForBofZ;
|
||||
for (unsigned int i(2); i<par_vec.size(); i++)
|
||||
parForBofZ.push_back(par_vec[i]);
|
||||
|
||||
vector<double> 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<double> 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<double> zonk;
|
||||
vector<double> donk;
|
||||
vector<double> hurgaB;
|
||||
vector<double> 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<hurgaB.size(); i++) {
|
||||
of7 << hurgaB[i] << " " << donk[i] << " " << hurgaPB[i] << endl;
|
||||
}
|
||||
of7.close();
|
||||
}
|
||||
*/
|
||||
// vector<double> 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<double> zz;
|
||||
// vector<double> 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<double> hurgaB(PofB.DataB());
|
||||
// vector<double> hurgaPB(PofB.DataPB());
|
||||
// vector<double> 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<hurgaB.size(); i++) {
|
||||
// of7 << hurgaB[i] << " " << hurgaPB[i] << " " << hurgaPB1[i] << endl;
|
||||
// }
|
||||
// of7.close();
|
||||
// }
|
||||
|
||||
// TPofTCalc poft("/home/l_wojek/analysis/WordsOfWisdom.dat", parForPofT);
|
||||
//
|
||||
// poft.DoFFT(PofB);
|
||||
// poft.CalcPol(parForPofT);
|
||||
//
|
||||
// char debugfilex[50];
|
||||
// int nx = sprintf (debugfilex, "4Ltest-P_t-%.4f_l-%.3f_E-%.1f_normal.dat", par_vec[2], par_vec[9], par_vec[1]);
|
||||
//
|
||||
// if (nx > 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<double> hurgaB1(PofB.DataB());
|
||||
vector<double> 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<hurgaB1.size(); i++) {
|
||||
of1 << hurgaB1[i] << " " << hurgaPB1[i] << endl;
|
||||
}
|
||||
of1.close();
|
||||
}
|
||||
|
||||
calcData.ConvolveGss(10.0, par_vec[1]);
|
||||
calcData.Normalize(par_vec[1]);
|
||||
|
||||
TPofBCalc PofB2(BofZ, calcData, parForPofB);
|
||||
|
||||
vector<double> z2(calcData.DataZ(par_vec[1]));
|
||||
vector<double> nz2(calcData.DataNZ(par_vec[1]));
|
||||
|
||||
ofstream of2("Implantation-profile-broad.dat");
|
||||
for (unsigned int i(0); i<z2.size(); i++) {
|
||||
of2 << z2[i] << " " << nz2[i] << endl;
|
||||
}
|
||||
of2.close();
|
||||
|
||||
vector<double> hurgaB2(PofB2.DataB());
|
||||
vector<double> 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<hurgaB2.size(); i++) {
|
||||
of8 << hurgaB2[i] << " " << hurgaPB2[i] << endl;
|
||||
}
|
||||
of8.close();
|
||||
}
|
||||
|
||||
PofB2.ConvolveGss(7.5);
|
||||
|
||||
vector<double> hurgaB3(PofB2.DataB());
|
||||
vector<double> 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<hurgaB3.size(); i++) {
|
||||
of9 << hurgaB3[i] << " " << hurgaPB3[i] << endl;
|
||||
}
|
||||
of9.close();
|
||||
}
|
||||
|
||||
z.clear();
|
||||
nz.clear();
|
||||
z2.clear();
|
||||
nz2.clear();
|
||||
*/
|
||||
/*
|
||||
|
||||
double param[8] = {100.0, 5.0, 50.0, 100.0, 40.0, 0.01, 0.1, 15.0};
|
||||
vector<double> parameter(param,param+8);
|
||||
vector<double> param_for_BofZ;
|
||||
vector<double> param_for_PofB;
|
||||
vector<double> 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<bofz.DataZ().size(); i++) {
|
||||
of6 << bofz.DataZ()[i] << " " << bofz.DataBZ()[i] << endl;
|
||||
}
|
||||
of6.close();
|
||||
|
||||
TPofBCalc pofb(bofz, calcData, param_for_PofB);
|
||||
|
||||
cout << "Output to file now..." << endl;
|
||||
|
||||
vector<double> hurgaB(pofb.DataB());
|
||||
vector<double> hurgaPB(pofb.DataPB());
|
||||
|
||||
ofstream of7("test_BpB.dat");
|
||||
for (unsigned int i(0); i<hurgaB.size(); i++) {
|
||||
of7 << hurgaB[i] << " " << hurgaPB[i] << endl;
|
||||
}
|
||||
of7.close();
|
||||
|
||||
TPofTCalc poft(param_for_PofT);
|
||||
|
||||
poft.DoFFT(pofb, param_for_PofT);
|
||||
|
||||
ofstream of8("test_tpt4.dat");
|
||||
for (double i(0.); i<12.0; i+=0.003) {
|
||||
of8 << i << " " << poft(i) << endl;
|
||||
}
|
||||
of8.close();
|
||||
|
||||
*/
|
||||
/**************** Test TLondon1DHS *********************************
|
||||
|
||||
// unsigned int parNo_arr[] = {1, 2, 5, 7, 9, 10, 11, 12};
|
||||
double par_arr[] = {24.4974, 22.0, 95.8253, 7.62096, 143.215};
|
||||
|
||||
// vector<unsigned int> parNo_vec(parNo_arr, parNo_arr+(sizeof(parNo_arr)/sizeof(parNo_arr[0])));
|
||||
vector<double> par_vec(par_arr, par_arr+(sizeof(par_arr)/sizeof(par_arr[0])));
|
||||
|
||||
// vector<double> par_vec_sub;
|
||||
|
||||
// for(unsigned int i(0); i<parNo_vec.size(); i++) {
|
||||
// par_vec_sub.push_back(par_vec[parNo_vec[i]-1]);
|
||||
// }
|
||||
|
||||
TLondon1DHS fitter;
|
||||
|
||||
************************************************************************/
|
||||
|
||||
|
||||
/**************** Test TLondon1D1L *********************************
|
||||
|
||||
unsigned int parNo_arr[] = {1, 3, 5, 7, 9, 10, 11, 12};
|
||||
double par_arr[] = {0.0, 999.0, 0.01, 999.0, 0.01, 999.0, 4.6, 999.0, 100.0, 5.0, 190.0, 180.0};
|
||||
|
||||
vector<unsigned int> parNo_vec(parNo_arr, parNo_arr+(sizeof(parNo_arr)/sizeof(parNo_arr[0])));
|
||||
vector<double> par_vec(par_arr, par_arr+(sizeof(par_arr)/sizeof(par_arr[0])));
|
||||
|
||||
vector<double> par_vec_sub;
|
||||
|
||||
for(unsigned int i(0); i<parNo_vec.size(); i++) {
|
||||
par_vec_sub.push_back(par_vec[parNo_vec[i]-1]);
|
||||
}
|
||||
|
||||
TLondon1D1L fitter(parNo_vec, par_vec);
|
||||
|
||||
************************************************************************/
|
||||
|
||||
/**************** Test TLondon1D2L **********************************
|
||||
|
||||
unsigned int parNo_arr[] = {1, 3, 5, 7, 9, 10, 11, 12, 13, 14, 15, 16};
|
||||
double par_arr[] = {0.0, 999.0, 0.01, 999.0, 0.01, 999.0, 21.6, 999.0, 100.0, 5.0, 70.0, 70.0, 180.0, 500.0, 1.0, 0.3};
|
||||
|
||||
vector<unsigned int> parNo_vec(parNo_arr, parNo_arr+(sizeof(parNo_arr)/sizeof(parNo_arr[0])));
|
||||
vector<double> par_vec(par_arr, par_arr+(sizeof(par_arr)/sizeof(par_arr[0])));
|
||||
|
||||
vector<double> par_vec_sub;
|
||||
|
||||
for(unsigned int i(0); i<parNo_vec.size(); i++) {
|
||||
par_vec_sub.push_back(par_vec[parNo_vec[i]-1]);
|
||||
}
|
||||
|
||||
TLondon1D2L fitter(parNo_vec, par_vec);
|
||||
|
||||
************************************************************************/
|
||||
|
||||
/**************** Test TLondon1D3L ********************************
|
||||
|
||||
unsigned int parNo_arr[] = {1, 7, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19};
|
||||
double par_arr[] = {0.0, 999.0, 0.01, 999.0, 0.01, 999.0, 4.6, 999.0, 100.0, 5.0, 70.0, 50.0, 70.0, 180.0, 180.0, 180.0, 1.0, 1.0, 1.0};
|
||||
|
||||
vector<unsigned int> parNo_vec(parNo_arr, parNo_arr+(sizeof(parNo_arr)/sizeof(parNo_arr[0])));
|
||||
vector<double> par_vec(par_arr, par_arr+(sizeof(par_arr)/sizeof(par_arr[0])));
|
||||
|
||||
vector<double> par_vec_sub;
|
||||
|
||||
for(unsigned int i(0); i<parNo_vec.size(); i++) {
|
||||
par_vec_sub.push_back(par_vec[parNo_vec[i]-1]);
|
||||
}
|
||||
|
||||
TLondon1D3L fitter(parNo_vec, par_vec);
|
||||
|
||||
************************************************************************/
|
||||
|
||||
|
||||
/**************** Test TLondon1D3LS *********************************
|
||||
|
||||
unsigned int parNo_arr[] = {1, 3, 5, 7, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18};
|
||||
double par_arr[] = {0.0, 999.0, 0.01, 999.0, 0.01, 999.0, 21.6, 999.0, 100.0, 5.0, 70.0, 50.0, 70.0, 180.0, 500.0, 1.0, 0.3, 1.0};
|
||||
|
||||
vector<unsigned int> parNo_vec(parNo_arr, parNo_arr+(sizeof(parNo_arr)/sizeof(parNo_arr[0])));
|
||||
vector<double> par_vec(par_arr, par_arr+(sizeof(par_arr)/sizeof(par_arr[0])));
|
||||
|
||||
vector<double> par_vec_sub;
|
||||
|
||||
for(unsigned int i(0); i<parNo_vec.size(); i++) {
|
||||
par_vec_sub.push_back(par_vec[parNo_vec[i]-1]);
|
||||
}
|
||||
|
||||
TLondon1D3LS fitter(parNo_vec, par_vec);
|
||||
|
||||
************************************************************************/
|
||||
|
||||
// ofstream of01("test_fitter01.dat");
|
||||
// ofstream of02("test_fitter02.dat");
|
||||
// ofstream of03("test_fitter03.dat");
|
||||
// ofstream of04("test_fitter04.dat");
|
||||
// ofstream of05("test_fitter05.dat");
|
||||
// ofstream of06("test_fitter06.dat");
|
||||
// ofstream of07("test_fitter07.dat");
|
||||
// ofstream of08("test_fitter08.dat");
|
||||
// ofstream of09("test_fitter09.dat");
|
||||
// ofstream of10("test_fitter10.dat");
|
||||
|
||||
|
||||
// for (double i(0.); i<12.0; i+=0.003) {
|
||||
// of01 << i << " " << fitter(i, par_vec) << endl;
|
||||
// }
|
||||
// of01.close();
|
||||
|
||||
// par_vec[1] = 7.7;
|
||||
//
|
||||
// for (double i(0.); i<12.0; i+=0.003) {
|
||||
// of02 << i << " " << fitter(i, par_vec) << endl;
|
||||
// }
|
||||
// of02.close();
|
||||
//
|
||||
// par_vec[0] = 0.0;
|
||||
//
|
||||
// for (double i(0.); i<12.0; i+=0.003) {
|
||||
// of03 << i << " " << fitter(i, par_vec) << endl;
|
||||
// }
|
||||
// of03.close();
|
||||
//
|
||||
// par_vec[2] = 200.0;
|
||||
//
|
||||
// for (double i(0.); i<12.0; i+=0.003) {
|
||||
// of04 << i << " " << fitter(i, par_vec) << endl;
|
||||
// }
|
||||
// of04.close();
|
||||
//
|
||||
// par_vec[4] = 100.0;
|
||||
//
|
||||
// for (double i(0.); i<12.0; i+=0.003) {
|
||||
// of05 << i << " " << fitter(i, par_vec) << endl;
|
||||
// }
|
||||
// of05.close();
|
||||
//
|
||||
// par_vec[0] = 20.0;
|
||||
// par_vec[1] = 24.6;
|
||||
// par_vec[2] = 96.5;
|
||||
// par_vec[3] = 15.0;
|
||||
// par_vec[4] = 130.0;
|
||||
//
|
||||
// for (double i(0.); i<12.0; i+=0.003) {
|
||||
// of06 << i << " " << fitter(i, par_vec) << endl;
|
||||
// }
|
||||
// of06.close();
|
||||
//
|
||||
// par_vec[0] = 20.0;
|
||||
// par_vec[1] = 24.6;
|
||||
// par_vec[2] = 96.5;
|
||||
// par_vec[3] = 15.0;
|
||||
// par_vec[4] = 140.0;
|
||||
//
|
||||
// for (double i(0.); i<12.0; i+=0.003) {
|
||||
// of07 << i << " " << fitter(i, par_vec) << endl;
|
||||
// }
|
||||
// of07.close();
|
||||
//
|
||||
// par_vec[0] = 20.0;
|
||||
// par_vec[1] = 24.6;
|
||||
// par_vec[2] = 96.5;
|
||||
// par_vec[3] = 20.0;
|
||||
// par_vec[4] = 130.0;
|
||||
//
|
||||
// for (double i(0.); i<12.0; i+=0.003) {
|
||||
// of08 << i << " " << fitter(i, par_vec) << endl;
|
||||
// }
|
||||
// of08.close();
|
||||
/*
|
||||
par_vec_sub[0] = 0.0;
|
||||
par_vec_sub[7] = 1000.0;
|
||||
|
||||
for (double i(0.); i<12.0; i+=0.003) {
|
||||
of09 << i << " " << fitter.Eval(i, par_vec_sub) << endl;
|
||||
}
|
||||
of09.close();
|
||||
|
||||
par_vec_sub[0] = 0.0;
|
||||
par_vec_sub[7] = 10.0;
|
||||
|
||||
for (double i(0.); i<12.0; i+=0.003) {
|
||||
of10 << i << " " << fitter.Eval(i, par_vec_sub) << endl;
|
||||
}
|
||||
of10.close();
|
||||
*/
|
||||
/*
|
||||
vector<double> 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;
|
||||
}
|
||||
@ -1,63 +1,169 @@
|
||||
#---------------------------------------------------
|
||||
# get compilation flags from root-config
|
||||
# Makefile.libGapIntegrals
|
||||
#
|
||||
# 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)
|
||||
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
|
||||
#
|
||||
|
||||
OS = LINUX
|
||||
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 -Wall -Wno-trigraphs -fPIC
|
||||
MUSRFITINCLUDE = ../../include
|
||||
#MUSRFITINCLUDE = /home/l_wojek/rep/analysis/musrfit/src/include
|
||||
LOCALINCLUDE = .
|
||||
ROOTINCLUDE = $(ROOTSYS)/include
|
||||
INCLUDES = -I$(LOCALINCLUDE) -I$(MUSRFITINCLUDE) -I$(ROOTINCLUDE)
|
||||
PMUSRPATH = ../../include
|
||||
MNPATH = $(ROOTSYS)/include
|
||||
GSLPATH = /usr/include/gsl
|
||||
CUBAPATH = /usr/local/include
|
||||
INCLUDES = -I$(PMUSRPATH) -I$(MNPATH) -I$(GSLPATH) -I$(CUBAPATH)
|
||||
LD = g++-4.4.0
|
||||
LDFLAGS = -O3
|
||||
LDFLAGS = -O
|
||||
SOFLAGS = -shared
|
||||
SHLIB = libLFRelaxation.so
|
||||
endif
|
||||
|
||||
# -- Windows/Cygwin
|
||||
ifeq ($(OS),WIN32GCC)
|
||||
CXX = g++
|
||||
CXXFLAGS = -O3 -Wall -Wno-trigraphs -D_DLL
|
||||
PMUSRPATH = ../../include
|
||||
MNPATH = $(ROOTSYS)/include
|
||||
GSLPATH = /usr/include/gsl
|
||||
CUBAPATH = /usr/local/include
|
||||
INCLUDES = -I$(PMUSRPATH) -I$(MNPATH) -I$(GSLPATH) -I$(CUBAPATH)
|
||||
LD = g++
|
||||
LDFLAGS = -O -Wl,--enable-auto-import -Wl,--enable-runtime-pseudo-reloc
|
||||
SOFLAGS = -shared -Wl,--export-all-symbols
|
||||
SHLIB = libLFRelaxation.dll
|
||||
endif
|
||||
|
||||
# -- MacOSX/Darwin
|
||||
ifeq ($(OS),DARWIN)
|
||||
CXX = g++
|
||||
CXXFLAGS = -O3 -Wall -Wno-trigraphs -fPIC
|
||||
PMUSRPATH = ../../include
|
||||
MNPATH = $(ROOTSYS)/include
|
||||
FINKPATH = /sw/include
|
||||
GSLPATH = $(FINKPATH)/gsl
|
||||
CUBAPATH = $(FINKPATH)
|
||||
INCLUDES = -I$(PMUSRPATH) -I$(MNPATH) -I$(GSLPATH) -I$(CUBAPATH)
|
||||
LD = g++
|
||||
LDFLAGS = -O -Xlinker -bind_at_load
|
||||
SOFLAGS = -dynamiclib -flat_namespace -undefined suppress -Wl,-x
|
||||
SHLIB = libLFRelaxation.dylib
|
||||
endif
|
||||
|
||||
# the output from the root-config script:
|
||||
CXXFLAGS += $(ROOTCFLAGS)
|
||||
LDFLAGS +=
|
||||
|
||||
# the ROOT libraries (G = graphic)
|
||||
LIBS = $(ROOTLIBS)
|
||||
GLIBS = $(ROOTGLIBS)
|
||||
|
||||
# GSL lib
|
||||
GSLLIB = -lgslcblas -lgsl
|
||||
# Cuba lib
|
||||
CUBALIB = -L/usr/local/lib -lcuba -lm
|
||||
# PMusr lib
|
||||
PMUSRLIB = -L$(ROOTSYS)/lib -lPMusr
|
||||
|
||||
ifeq ($(OS),WIN32GCC)
|
||||
# GSL lib
|
||||
GSLLIB = -L/usr/lib -lgslcblas -lgsl
|
||||
# Cuba lib
|
||||
CUBALIB = -L/usr/local/lib -lcuba -lm
|
||||
# PMusr lib
|
||||
PMUSRLIB = -L$(ROOTSYS)/lib -lPMusr -lMathMore
|
||||
endif
|
||||
|
||||
ifeq ($(OS),DARWIN)
|
||||
# GSL lib
|
||||
GSLLIB = -L/sw/lib -lgslcblas -lgsl
|
||||
# Cuba lib
|
||||
CUBALIB = -L/sw/lib -lcuba -lm
|
||||
# PMusr lib
|
||||
PMUSRLIB = -L$(ROOTSYS)/lib -lPMusr -lMathMore
|
||||
endif
|
||||
|
||||
# some definitions: headers (used to generate *Dict* stuff), sources, objects,...
|
||||
OBJS =
|
||||
OBJS += TLFRelaxation.o TLFRelaxationDict.o
|
||||
|
||||
SHLIB = libLFRelaxation.so
|
||||
EXTOBJS = BMWIntegrator.o
|
||||
|
||||
# make the shared libs:
|
||||
|
||||
# make the shared lib:
|
||||
#
|
||||
all: $(SHLIB)
|
||||
|
||||
$(SHLIB): $(OBJS)
|
||||
$(SHLIB): $(EXTOBJS) $(OBJS)
|
||||
@echo "---> Building shared library $(SHLIB) ..."
|
||||
/bin/rm -f $(SHLIB)
|
||||
$(LD) $(LDFLAGS) $(OBJS) $(SOFLAGS) -o $(SHLIB)
|
||||
$(LD) $(SOFLAGS) $(LDFLAGS) -o $(SHLIB) $(EXTOBJS) $(OBJS) $(LIBS) $(PMUSRLIB) $(GSLLIB) $(CUBALIB)
|
||||
@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) $(SHLIB) *Dict* core* *~
|
||||
@echo "---> removing $(OBJS) $(SHLIB)"
|
||||
clean:; @rm -f $(OBJS) $(EXTOBJS) $(SHLIB) *Dict* core* *~
|
||||
@echo "---> removing $(OBJS) $(EXTOBJS) $(SHLIB)"
|
||||
|
||||
#
|
||||
BMWIntegrator.o: ../BMWIntegrator/BMWIntegrator.cpp
|
||||
$(CXX) $(INCLUDES) $(CXXFLAGS) -c $<
|
||||
|
||||
$(OBJS): %.o: %.cpp
|
||||
$(CXX) $(INCLUDES) $(CXXFLAGS) -c $<
|
||||
|
||||
# Generate the ROOT CINT dictionary
|
||||
|
||||
TLFRelaxationDict.cpp: ./TLFRelaxation.h ./TLFRelaxationLinkDef.h
|
||||
TLFRelaxationDict.cpp: ../BMWIntegrator/BMWIntegrator.h ./TLFRelaxation.h ./TLFRelaxationLinkDef.h
|
||||
@echo "Generating dictionary $@..."
|
||||
rootcint -f $@ -c -p -I$(MUSRFITINCLUDE) $^
|
||||
rootcint -f $@ -c -p $^
|
||||
|
||||
install: all
|
||||
@echo "Installing shared lib: libLFRelaxation.so"
|
||||
@echo "Installing shared lib: $(SHLIB)"
|
||||
ifeq ($(OS),LINUX)
|
||||
cp -pv $(SHLIB) $(ROOTSYS)/lib
|
||||
cp -pv $(LOCALINCLUDE)/*.h $(ROOTSYS)/include
|
||||
cp -pv ../BMWIntegrator/BMWIntegrator.h TLFRelaxation.h $(ROOTSYS)/include
|
||||
endif
|
||||
ifeq ($(OS),WIN32GCC)
|
||||
cp -pv $(SHLIB) $(ROOTSYS)/bin
|
||||
ln -sf $(ROOTSYS)/bin/$(SHLIB) $(ROOTSYS)/lib/$(SHLIB)
|
||||
cp -pv ../BMWIntegrator/BMWIntegrator.h TLFRelaxation.h $(ROOTSYS)/include
|
||||
endif
|
||||
ifeq ($(OS),DARWIN)
|
||||
cp -pv $(SHLIB) $(ROOTSYS)/lib
|
||||
cp -pv ../BMWIntegrator/BMWIntegrator.h TLFRelaxation.h $(ROOTSYS)/include
|
||||
endif
|
||||
|
||||
cleaninstall: clean install
|
||||
|
||||
154
src/external/libLFRelaxation/TLFRelaxation.cpp
vendored
154
src/external/libLFRelaxation/TLFRelaxation.cpp
vendored
@ -9,11 +9,31 @@
|
||||
|
||||
***************************************************************************/
|
||||
|
||||
/***************************************************************************
|
||||
* Copyright (C) 2009 by Bastian M. Wojek *
|
||||
* bastian.wojek@psi.ch *
|
||||
* *
|
||||
* 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 <cassert>
|
||||
#include <sys/time.h>
|
||||
#include <iostream>
|
||||
|
||||
#include "TIntegrator.h"
|
||||
#include "../BMWIntegrator/BMWIntegrator.h"
|
||||
#include "TLFRelaxation.h"
|
||||
|
||||
using namespace std;
|
||||
@ -26,6 +46,7 @@ ClassImp(TLFStatGssKT)
|
||||
ClassImp(TLFStatLorKT)
|
||||
ClassImp(TLFDynGssKT)
|
||||
ClassImp(TLFDynLorKT)
|
||||
ClassImp(TLFSGInterpolation)
|
||||
|
||||
|
||||
// LF Static Gaussian KT
|
||||
@ -96,7 +117,17 @@ double TLFStatLorKT::operator()(double t, const vector<double> &par) const {
|
||||
double coeff2(coeff1*coeff1);
|
||||
double coeff3((1.0+coeff2)*par[1]);
|
||||
|
||||
return (1.0-(coeff1*TMath::Exp(-par[1]*t)*TMath::BesselJ1(omegaL*t))-(coeff2*(TMath::BesselJ0(omegaL*t)*TMath::Exp(-par[1]*t)-1.0))-coeff3*fIntBesselJ0Exp->IntegrateFunc(0.,t));
|
||||
double w0t(omegaL*t);
|
||||
double j1, j0;
|
||||
if (fabs(w0t) < 0.001) { // check zero time limits of the spherical bessel functions j0(x) and j1(x)
|
||||
j0 = 1.0;
|
||||
j1 = 0.0;
|
||||
} else {
|
||||
j0 = TMath::Sin(w0t)/w0t;
|
||||
j1 = (TMath::Sin(w0t)-w0t*TMath::Cos(w0t))/(w0t*w0t);
|
||||
}
|
||||
|
||||
return (1.0-(coeff1*TMath::Exp(-par[1]*t)*j1)-(coeff2*(j0*TMath::Exp(-par[1]*t)-1.0))-coeff3*fIntBesselJ0Exp->IntegrateFunc(0.,t));
|
||||
}
|
||||
|
||||
// LF Dynamic Gaussian KT
|
||||
@ -188,7 +219,7 @@ double TLFDynGssKT::operator()(double t, const vector<double> &par) const {
|
||||
double sigsqtsq(sigsq*t*t);
|
||||
return 0.33333333333333333333+0.66666666666666666667*(1.0-sigsqtsq)*exp(-0.5*sigsqtsq); // ZF static Gauss KT
|
||||
}
|
||||
return exp(-2.0*sigsq/nusq*(enut-1.0+nut)); // ZF Abragam Delta^2->2*Delta^2
|
||||
// return exp(-2.0*sigsq/nusq*(enut-1.0+nut)); // ZF Abragam Delta^2->2*Delta^2 (only here for TESTING - DO NOT return this value in an production environment!!!!!!)
|
||||
}
|
||||
|
||||
if((par[2] >= 5.0*par[1]) || (omegaL >= 10.0*par[1])) {
|
||||
@ -272,7 +303,7 @@ double TLFDynGssKT::operator()(double t, const vector<double> &par) const {
|
||||
|
||||
// LF Dynamic Lorentz KT
|
||||
|
||||
TLFDynLorKT::TLFDynLorKT() : fCalcNeeded(true), fFirstCall(true), fWisdom("/home/l_wojek/analysis/WordsOfWisdom.dat"), fNSteps(524288), fDt(0.000040), fCounter(0), fA(1.0), fL1(0.0), fL2(0.0) {
|
||||
TLFDynLorKT::TLFDynLorKT() : fCalcNeeded(true), fFirstCall(true), fWisdom("/home/l_wojek/analysis/WordsOfWisdom.dat"), fNSteps(524288), fDt(0.000040), fCounter(0), fL1(0.0), fL2(0.0) {
|
||||
// Calculate d_omega and C for given NFFT and dt
|
||||
fDw = TMath::Pi()/fNSteps/fDt;
|
||||
fC = 2.0*TMath::Log(double(fNSteps))/(double(fNSteps-1)*fDt);
|
||||
@ -283,14 +314,14 @@ TLFDynLorKT::TLFDynLorKT() : fCalcNeeded(true), fFirstCall(true), fWisdom("/home
|
||||
FILE *wordsOfWisdomR;
|
||||
wordsOfWisdomR = fopen(fWisdom.c_str(), "r");
|
||||
if (wordsOfWisdomR == NULL) {
|
||||
cout << "TLFDynLorKT::TLFDynGssKT: Couldn't open wisdom file ..." << endl;
|
||||
cout << "TLFDynLorKT::TLFDynLorKT: Couldn't open wisdom file ..." << endl;
|
||||
} else {
|
||||
wisdomLoaded = fftw_import_wisdom_from_file(wordsOfWisdomR);
|
||||
fclose(wordsOfWisdomR);
|
||||
}
|
||||
|
||||
if (!wisdomLoaded) {
|
||||
cout << "TLFDynLorKT::TLFDynGssKT: No wisdom is imported..." << endl;
|
||||
cout << "TLFDynLorKT::TLFDynLorKT: No wisdom is imported..." << endl;
|
||||
}
|
||||
// END of WisdomLoading
|
||||
|
||||
@ -307,7 +338,7 @@ TLFDynLorKT::~TLFDynLorKT() {
|
||||
FILE *wordsOfWisdomW;
|
||||
wordsOfWisdomW = fopen(fWisdom.c_str(), "w");
|
||||
if (wordsOfWisdomW == NULL) {
|
||||
cout << "TLFDynGssKT::~TLFDynGssKT: Could not open file ... No wisdom is exported..." << endl;
|
||||
cout << "TLFDynLorKT::~TLFDynLorKT: Could not open file ... No wisdom is exported..." << endl;
|
||||
} else {
|
||||
fftw_export_wisdom_to_file(wordsOfWisdomW);
|
||||
fclose(wordsOfWisdomW);
|
||||
@ -322,10 +353,6 @@ TLFDynLorKT::~TLFDynLorKT() {
|
||||
cout << "TLFDynLorKT::~TLFDynLorKT(): " << fCounter << " full FFT cyles needed..." << endl;
|
||||
}
|
||||
|
||||
const double TLFDynLorKT::fX[16]={1.61849, 1.27059, 0.890315, 6.72608, 0.327258, 0.981975, 1.5964, 2.31023, 2.37689, 5.63945, 0.235734, 1.10617, 1.1664, 0.218402, 0.266562, 0.471331};
|
||||
|
||||
const double TLFDynLorKT::fY[20]={1.54699, 0.990321, 0.324923, 0.928399, 2.11155, 0.615106, 1.49907, 0.679423, 6.17621, 1.38907, 0.979608, 0.0593631, 0.806427, 3.33144, 1.05059, 1.29922, 0.254661, 0.284348, 0.991419, 0.375213};
|
||||
|
||||
double TLFDynLorKT::operator()(double t, const vector<double> &par) const {
|
||||
|
||||
assert(par.size() == 3); // three parameters nuL=gbar*B,sigma,fluct.rate nu
|
||||
@ -349,6 +376,8 @@ double TLFDynLorKT::operator()(double t, const vector<double> &par) const {
|
||||
}
|
||||
|
||||
double omegaL(TWOPI*par[0]); // Larmor frequency
|
||||
double a(par[1]); // static width
|
||||
double nu(par[2]); // hopping rate
|
||||
|
||||
if(fabs(par[0])<0.00135538817){ // Zero Field
|
||||
if(!par[2]){ // nu=0
|
||||
@ -357,37 +386,43 @@ double TLFDynLorKT::operator()(double t, const vector<double> &par) const {
|
||||
}
|
||||
}
|
||||
|
||||
// semi-analytical approximations for {nu >= 4a} and {omegaL >= 10a}
|
||||
// semi-analytical approximation for {nu >= 5a} and {omegaL >= 30a}
|
||||
|
||||
if(par[2]){ // nu!=0
|
||||
if(par[2] >= 4.0*par[1]){ // nu >= 4a
|
||||
if(fCalcNeeded){
|
||||
double wLnu(omegaL/par[2]); // omegaL/nu
|
||||
double wLnusq(wLnu*wLnu); // (omegaL/nu)^2
|
||||
double anu(par[1]/par[2]); // a/nu
|
||||
// check if hopping > 5 * damping, of Larmor angular frequency is > 30 * damping (BMW limit)
|
||||
if((par[2] > 5.0*par[1]) || (omegaL > 30.0*par[1])){
|
||||
if(fCalcNeeded){
|
||||
|
||||
fA=1.0+wLnu*anu*(-fX[0]/(wLnusq+fX[1])+fX[2]/(wLnusq*wLnu+fX[3])+fX[4]/(wLnusq*wLnusq+fX[5]));
|
||||
fL1=par[1]*(fX[6]/(wLnu+fX[7])+fX[8]/(wLnusq+fX[9])+fX[10]/(wLnusq*wLnusq+fX[11]));
|
||||
fL2=par[2]*(fX[12]*tanh(fX[13]*wLnu)+fX[14]*wLnu+fX[15]);
|
||||
|
||||
fCalcNeeded=false;
|
||||
// 'c' and 'd' are parameters BMW obtained by fitting large parameter space LF-curves to the model below
|
||||
const double c[7] = {1.15331, 1.64826, -0.71763, 3.0, 0.386683, -5.01876, 2.41854};
|
||||
const double d[4] = {2.44056, 2.92063, 1.69581, 0.667277};
|
||||
double w0N[4];
|
||||
double nuN[4];
|
||||
w0N[0] = omegaL;
|
||||
nuN[0] = nu;
|
||||
for (unsigned int i=1; i<4; i++) {
|
||||
w0N[i] = omegaL * w0N[i-1];
|
||||
nuN[i] = nu * nuN[i-1];
|
||||
}
|
||||
return fA*exp(-fL1*t)+(1.0-fA)*exp(-fL2*t)*cos(omegaL*t);
|
||||
}
|
||||
if(omegaL >= 10.0*par[1]){ // omegaL >= 10a
|
||||
if(fCalcNeeded){
|
||||
double wLnu(omegaL/par[2]); // omegaL/nu
|
||||
double wLnusq(wLnu*wLnu); // (omegaL/nu)^2
|
||||
double anu(par[1]/par[2]); // a/nu
|
||||
double denom(w0N[3]+d[0]*w0N[2]*nuN[0]+d[1]*w0N[1]*nuN[1]+d[2]*w0N[0]*nuN[2]+d[3]*nuN[3]);
|
||||
fL1 = (c[0]*w0N[2]+c[1]*w0N[1]*nuN[0]+c[2]*w0N[0]*nuN[1])/denom;
|
||||
fL2 = (c[3]*w0N[2]+c[4]*w0N[1]*nuN[0]+c[5]*w0N[0]*nuN[1]+c[6]*nuN[2])/denom;
|
||||
|
||||
fA=1.0-fY[0]*pow(anu,fY[1])/(wLnu+fY[2])+fY[3]*pow(anu,fY[4])/(wLnusq-fY[5])+fY[6]*pow(anu,fY[7])/(wLnusq*wLnu+fY[8]);
|
||||
fL1=par[1]*(fY[9]/(pow(wLnu,fY[10])+fY[11])-fY[12]/(pow(wLnu,fY[13])+fY[14]));
|
||||
fL2=par[2]*(fY[15]*tanh(fY[16]*wLnu)+fY[17]*pow(wLnu,fY[18])+fY[19]);
|
||||
|
||||
fCalcNeeded=false;
|
||||
}
|
||||
return fA*exp(-fL1*t)+(1.0-fA)*exp(-fL2*t)*cos(omegaL*t);
|
||||
fCalcNeeded=false;
|
||||
}
|
||||
|
||||
double w0t(omegaL*t);
|
||||
double j1, j0;
|
||||
if (fabs(w0t) < 0.001) { // check zero time limits of the spherical bessel functions j0(x) and j1(x)
|
||||
j0 = 1.0;
|
||||
j1 = 0.0;
|
||||
} else {
|
||||
j0 = TMath::Sin(w0t)/w0t;
|
||||
j1 = (TMath::Sin(w0t)-w0t*TMath::Cos(w0t))/(w0t*w0t);
|
||||
}
|
||||
|
||||
double Gamma_t = -4.0/3.0*a*(fL1*(1.0-j0*TMath::Exp(-nu*t))+fL2*j1*TMath::Exp(-nu*t)+(1.0-fL2*omegaL/3.0-fL1*nu)*t);
|
||||
|
||||
return TMath::Exp(Gamma_t);
|
||||
}
|
||||
|
||||
// if no approximation can be used -> Laplace transform
|
||||
@ -411,15 +446,17 @@ double TLFDynLorKT::operator()(double t, const vector<double> &par) const {
|
||||
|
||||
for(unsigned int i(1); i<fNSteps; i++) {
|
||||
tt=(double(i)-0.5)*fDt;
|
||||
fFFTtime[i]=TMath::BesselJ0(omegaL*tt)*TMath::Exp(-par[1]*tt)*fDt;
|
||||
fFFTtime[i]=TMath::Sin(omegaL*tt)/(omegaL*tt)*TMath::Exp(-par[1]*tt)*fDt;
|
||||
}
|
||||
|
||||
double totoIntegrale(0.);
|
||||
double w0t(1.0);
|
||||
|
||||
for(unsigned int i(1); i<fNSteps; i++) {
|
||||
tt=double(i)*fDt;
|
||||
totoIntegrale+=fFFTtime[i];
|
||||
fFFTtime[i]=(1.0-(coeff1*TMath::Exp(-par[1]*tt)*TMath::BesselJ1(omegaL*tt))-(coeff2*(TMath::BesselJ0(omegaL*tt)*TMath::Exp(-par[1]*tt)-1.0))-coeff3*totoIntegrale)*TMath::Exp(-(fC+par[2])*tt)*fDt;
|
||||
w0t = omegaL*tt;
|
||||
fFFTtime[i]=(1.0-(coeff1*TMath::Exp(-par[1]*tt)*(TMath::Sin(w0t)-w0t*TMath::Cos(w0t))/(w0t*w0t))-(coeff2*(TMath::Sin(w0t)/(w0t)*TMath::Exp(-par[1]*tt)-1.0))-coeff3*totoIntegrale)*TMath::Exp(-(fC+par[2])*tt)*fDt;
|
||||
}
|
||||
}
|
||||
|
||||
@ -453,3 +490,42 @@ double TLFDynLorKT::operator()(double t, const vector<double> &par) const {
|
||||
// return fFFTtime[int(t/fDt)];
|
||||
return fDw*TMath::Exp(fC*t)/TMath::Pi()*fFFTtime[int(t/fDt)];
|
||||
}
|
||||
|
||||
|
||||
// De Renzi Spin Glass Interpolation
|
||||
|
||||
TLFSGInterpolation::TLFSGInterpolation() {
|
||||
TIntSGInterpolation *integral = new TIntSGInterpolation();
|
||||
fIntegral = integral;
|
||||
integral = 0;
|
||||
}
|
||||
|
||||
TLFSGInterpolation::~TLFSGInterpolation() {
|
||||
delete fIntegral;
|
||||
fIntegral = 0;
|
||||
}
|
||||
|
||||
double TLFSGInterpolation::operator()(double t, const vector<double> &par) const {
|
||||
|
||||
assert(par.size() == 4); // two parameters nu=gbar*B [MHz], a [1/us], lambda [1/us], beta [1]
|
||||
|
||||
if((t <= 0.0) || (par[3] <= 0.0))
|
||||
return 1.0;
|
||||
|
||||
double expo(0.5*par[1]*par[1]*t*t/par[3]+par[2]*t);
|
||||
|
||||
if(TMath::Abs(par[0])<0.00135538817){
|
||||
return (0.33333333333333333333*TMath::Exp(-TMath::Power(par[2]*t,par[3])) + \
|
||||
0.66666666666666666667*(1.0-par[1]*par[1]*t*t/TMath::Power(expo,(1.0-par[3])))*TMath::Exp(-TMath::Power(expo,par[3])));
|
||||
}
|
||||
|
||||
vector<double> intpar(par);
|
||||
intpar.push_back(t);
|
||||
fIntegral->SetParameters(intpar);
|
||||
|
||||
double omegaL(TMath::TwoPi()*par[0]);
|
||||
|
||||
return (TMath::Exp(-TMath::Power(par[2]*t,par[3])) + 2.0*par[1]*par[1]/(omegaL*omegaL) * \
|
||||
((TMath::Cos(omegaL*t)-TMath::Sin(omegaL*t)/(omegaL*t))*TMath::Exp(-TMath::Power(expo,par[3]))/TMath::Power(expo,(1.0-par[3])) + \
|
||||
fIntegral->IntegrateFunc(0.,t)));
|
||||
}
|
||||
|
||||
39
src/external/libLFRelaxation/TLFRelaxation.h
vendored
39
src/external/libLFRelaxation/TLFRelaxation.h
vendored
@ -9,6 +9,26 @@
|
||||
|
||||
***************************************************************************/
|
||||
|
||||
/***************************************************************************
|
||||
* Copyright (C) 2009 by Bastian M. Wojek *
|
||||
* bastian.wojek@psi.ch *
|
||||
* *
|
||||
* 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 _TLFRelaxation_H_
|
||||
#define _TLFRelaxation_H_
|
||||
|
||||
@ -34,7 +54,7 @@ using namespace std;
|
||||
//#include "TMath.h"
|
||||
#include "PUserFcnBase.h"
|
||||
#include "fftw3.h"
|
||||
#include "TIntegrator.h"
|
||||
#include "../BMWIntegrator/BMWIntegrator.h"
|
||||
|
||||
class TLFStatGssKT : public PUserFcnBase {
|
||||
|
||||
@ -112,14 +132,25 @@ private:
|
||||
double *fFFTtime;
|
||||
fftw_complex *fFFTfreq;
|
||||
mutable unsigned int fCounter;
|
||||
static const double fX[16];
|
||||
static const double fY[20];
|
||||
mutable double fA;
|
||||
mutable double fL1;
|
||||
mutable double fL2;
|
||||
|
||||
ClassDef(TLFDynLorKT,1)
|
||||
};
|
||||
|
||||
class TLFSGInterpolation : public PUserFcnBase {
|
||||
|
||||
public:
|
||||
TLFSGInterpolation();
|
||||
~TLFSGInterpolation();
|
||||
|
||||
double operator()(double, const vector<double>&) const;
|
||||
|
||||
private:
|
||||
TIntSGInterpolation *fIntegral;
|
||||
|
||||
ClassDef(TLFSGInterpolation,1)
|
||||
};
|
||||
|
||||
|
||||
#endif //_LFRelaxation_H_
|
||||
|
||||
@ -9,6 +9,26 @@
|
||||
|
||||
***************************************************************************/
|
||||
|
||||
/***************************************************************************
|
||||
* Copyright (C) 2009 by Bastian M. Wojek *
|
||||
* bastian.wojek@psi.ch *
|
||||
* *
|
||||
* 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__
|
||||
|
||||
@ -20,6 +40,7 @@
|
||||
#pragma link C++ class TLFStatLorKT+;
|
||||
#pragma link C++ class TLFDynGssKT+;
|
||||
#pragma link C++ class TLFDynLorKT+;
|
||||
#pragma link C++ class TLFSGInterpolation+;
|
||||
|
||||
#endif //__CINT__
|
||||
// root dictionary stuff --------------------------------------------------
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user