From c990453194dfe4c01c70d907ec4431a816ee2e47 Mon Sep 17 00:00:00 2001 From: "Bastian M. Wojek" Date: Wed, 12 Nov 2008 18:55:50 +0000 Subject: [PATCH] Implemented some more London-multilayer-scenarios for testing purposes --- .../TFitPofB-lib/classes/TBofZCalc.cpp | 46 ++ .../TFitPofB-lib/classes/TLondon1D.cpp | 395 +++++++++++++++++- .../TFitPofB-lib/classes/TPofBCalc.cpp | 36 ++ .../classes/TTrimSPDataHandler.cpp | 46 ++ 4 files changed, 512 insertions(+), 11 deletions(-) diff --git a/src/external/TFitPofB-lib/classes/TBofZCalc.cpp b/src/external/TFitPofB-lib/classes/TBofZCalc.cpp index ce117a36..99332251 100644 --- a/src/external/TFitPofB-lib/classes/TBofZCalc.cpp +++ b/src/external/TFitPofB-lib/classes/TBofZCalc.cpp @@ -195,3 +195,49 @@ TLondon1D_3LS::TLondon1D_3LS(unsigned int steps, const vector ¶m) { fBZ.push_back(BBz); } } + +//------------------ +// Constructor of the TLondon1D_4L class +// 1D-London screening in a thin superconducting film, four layers, four lambdas +// Parameters: Bext[G], deadlayer[nm], thickness1[nm], thickness2[nm], thickness3[nm], thickness4[nm], +// lambda1[nm], lambda2[nm], lambda3[nm], lambda4[nm] +//------------------ + +TLondon1D_4L::TLondon1D_4L(unsigned int steps, const vector ¶m) { + + double N1((param[6]+exp(2.0*param[2]/param[6])*(param[6]-param[7])+param[7])*(-1.0*param[9]*cosh(param[5]/param[9])*(param[8]*cosh(param[4]/param[8])-param[7]*sinh(param[4]/param[8]))+param[8]*(param[7]*cosh(param[4]/param[8])-param[8]*sinh(param[4]/param[8]))*sinh(param[5]/param[9]))+exp(2.0*param[3]/param[7])*(param[6]-param[7]+exp(2.0*param[2]/param[6])*(param[6]+param[7]))*(param[9]*cosh(param[5]/param[9])*(param[8]*cosh(param[4]/param[8])+param[7]*sinh(param[4]/param[8]))+param[8]*(param[7]*cosh(param[4]/param[8])+param[8]*sinh(param[4]/param[8]))*sinh(param[5]/param[9]))); + + double N11(param[9]*cosh(param[5]/param[9])*(param[8]*cosh(param[4]/param[8])*(param[7]*cosh(param[3]/param[7])-param[6]*sinh(param[3]/param[7]))+param[7]*(-1.0*param[6]*cosh(param[3]/param[7])+param[7]*sinh(param[3]/param[7]))*sinh(param[4]/param[8]))+param[8]*(param[7]*cosh(param[4]/param[8])*(-1.0*param[6]*cosh(param[3]/param[7])+param[7]*sinh(param[3]/param[7]))+param[8]*(param[7]*cosh(param[3]/param[7])-param[6]*sinh(param[3]/param[7]))*sinh(param[4]/param[8]))*sinh(param[5]/param[9])); + + double N12(exp(2.0*(param[2]+param[1])/param[6])*(param[9]*cosh(param[5]/param[9])*(param[8]*cosh(param[4]/param[8])*(param[7]*cosh(param[3]/param[7])+param[6]*sinh(param[3]/param[7]))+param[7]*(param[6]*cosh(param[3]/param[7])+param[7]*sinh(param[3]/param[7]))*sinh(param[4]/param[8]))+param[8]*(param[7]*cosh(param[4]/param[8])*(param[6]*cosh(param[3]/param[7])+param[7]*sinh(param[3]/param[7]))+param[8]*(param[7]*cosh(param[3]/param[7])+param[6]*sinh(param[3]/param[7]))*sinh(param[4]/param[8]))*sinh(param[5]/param[9]))); + + double N2((param[6]+param[7]+(param[6]-param[7])*exp(2.0*param[2]/param[6]))*(-1.0*param[9]*cosh(param[5]/param[9])*(param[8]*cosh(param[4]/param[8])-param[7]*sinh(param[4]/param[8]))+param[8]*(param[7]*cosh(param[4]/param[8])-param[8]*sinh(param[4]/param[8]))*sinh(param[5]/param[9]))+exp(2.0*param[3]/param[7])*(param[6]-param[7]+exp(2.0*param[2]/param[6])*(param[6]+param[7]))*(param[9]*cosh(param[5]/param[9])*(param[8]*cosh(param[4]/param[8])+param[7]*sinh(param[4]/param[8]))+param[8]*(param[7]*cosh(param[4]/param[8])+param[8]*sinh(param[4]/param[8]))*sinh(param[5]/param[9]))); + + double N21(exp(param[3]/param[7])*param[8]*param[9]*(param[6]*cosh(param[2]/param[6])+param[7]*sinh(param[2]/param[6]))+param[6]*param[9]*cosh(param[5]/param[9])*(-1.0*param[8]*cosh(param[4]/param[8])+param[7]*sinh(param[4]/param[8]))+param[6]*param[8]*(param[7]*cosh(param[4]/param[8])-param[8]*sinh(param[4]/param[8]))*sinh(param[5]/param[9])); + + double N22(exp((2.0*param[2]+param[3]+2.0*param[1])/param[7])*(-1.0*param[6]*param[8]*param[9]*cosh(param[2]/param[6])+param[7]*param[8]*param[9]*sinh(param[2]/param[6])+exp(param[3]/param[7])*param[6]*param[9]*cosh(param[5]/param[9])*(param[8]*cosh(param[4]/param[8])+param[7]*sinh(param[4]/param[8]))+exp(param[3]/param[7])*param[6]*param[8]*(param[7]*cosh(param[4]/param[8])+param[8]*sinh(param[4]/param[8]))*sinh(param[5]/param[9]))); + + double N3(2.0*((param[6]+exp(2.0*param[2]/param[6])*(param[6]-param[7])+param[7])*(-1.0*param[9]*cosh(param[5]/param[9])*(param[8]*cosh(param[4]/param[8])-param[7]*sinh(param[4]/param[8]))+param[8]*(param[7]*cosh(param[4]/param[8])-param[8]*sinh(param[4]/param[8]))*sinh(param[5]/param[9]))+exp(2.0*param[3]/param[7])*(param[6]-param[7]+exp(2.0*param[2]/param[6])*(param[6]+param[7]))*(param[9]*cosh(param[5]/param[9])*(param[8]*cosh(param[4]/param[8])+param[7]*sinh(param[4]/param[8]))+param[8]*(param[7]*cosh(param[4]/param[8])+param[8]*sinh(param[4]/param[8]))*sinh(param[5]/param[9])))); + + double N4(4.0*((param[6]+exp(2.0*param[2]/param[6])*(param[6]-param[7])+param[7])*(-1.0*param[9]*cosh(param[5]/param[9])*(param[8]*cosh(param[4]/param[8])-param[7]*sinh(param[4]/param[8]))+param[8]*(param[7]*cosh(param[4]/param[8])-param[8]*sinh(param[4]/param[8]))*sinh(param[5]/param[9]))+exp(2.0*param[3]/param[7])*(param[6]-param[7]+exp(2.0*param[2]/param[6])*(param[6]+param[7]))*(param[9]*cosh(param[5]/param[9])*(param[8]*cosh(param[4]/param[8])+param[7]*sinh(param[4]/param[8]))+param[8]*(param[7]*cosh(param[4]/param[8])+param[8]*sinh(param[4]/param[8]))*sinh(param[5]/param[9])))); + + double N42(-1.0*param[6]*param[7]*param[8]+exp(param[5]/param[9])*(param[6]*cosh(param[2]/param[6])*(param[8]*sinh(param[3]/param[7])*(param[9]*cosh(param[4]/param[8])+param[8]*sinh(param[4]/param[8]))+param[7]*cosh(param[3]/param[7])*(param[8]*cosh(param[4]/param[8])+param[9]*sinh(param[4]/param[8])))+param[7]*sinh(param[2]/param[6])*(param[8]*cosh(param[3]/param[7])*(param[9]*cosh(param[4]/param[8])+param[8]*sinh(param[4]/param[8]))+param[7]*sinh(param[3]/param[7])*(param[8]*cosh(param[4]/param[8])+param[9]*sinh(param[4]/param[8]))))); + + fDZ = (param[2]+param[3]+param[4]+param[5])/double(steps); + double ZZ, BBz; + + for (unsigned int j(0); j +#include using namespace std; #include @@ -21,6 +22,8 @@ ClassImp(TLondon1D1L) ClassImp(TLondon1D2L) ClassImp(TLondon1D3L) ClassImp(TLondon1D3LS) +ClassImp(TLondon1D4L) +ClassImp(TLondon1D3LSub) //------------------ @@ -82,6 +85,28 @@ TLondon1D3LS::~TLondon1D3LS() { fPofT = 0; } +TLondon1D4L::~TLondon1D4L() { + fPar.clear(); + fParForBofZ.clear(); + fParForPofB.clear(); + fParForPofT.clear(); + delete fImpProfile; + fImpProfile = 0; + delete fPofT; + fPofT = 0; +} + +TLondon1D3LSub::~TLondon1D3LSub() { + fPar.clear(); + fParForBofZ.clear(); + fParForPofB.clear(); + fParForPofT.clear(); + delete fImpProfile; + fImpProfile = 0; + delete fPofT; + fPofT = 0; +} + //------------------ // Constructor of the TLondon1DHS class -- reading available implantation profiles and // creates (a pointer to) the TPofTCalc object (with the FFT plan) @@ -143,6 +168,9 @@ TLondon1DHS::TLondon1DHS() : fCalcNeeded(true), fFirstCall(true) { //------------------ double TLondon1DHS::operator()(double t, const vector &par) const { + + assert(par.size() == 5); + if(t<0.0) return 0.0; @@ -278,9 +306,11 @@ TLondon1D1L::TLondon1D1L() : fCalcNeeded(true), fFirstCall(true), fCallCounter(0 double TLondon1D1L::operator()(double t, const vector &par) const { + assert(par.size() == 6); + // Debugging // Count the number of function calls - fCallCounter++; +// fCallCounter++; if(t<0.0) return 0.0; @@ -350,15 +380,15 @@ double TLondon1D1L::operator()(double t, const vector &par) const { fCalcNeeded = false; } - // Debugging start - if (!(fCallCounter%10000)){ - cout << fCallCounter-1 << "\t"; - for (unsigned int i(0); iEval(t); @@ -423,6 +453,9 @@ TLondon1D2L::TLondon1D2L() : fCalcNeeded(true), fFirstCall(true), fLastTwoChange //------------------ double TLondon1D2L::operator()(double t, const vector &par) const { + + assert(par.size() == 10); + if(t<0.0) return 0.0; @@ -568,6 +601,9 @@ TLondon1D3L::TLondon1D3L() : fCalcNeeded(true), fFirstCall(true), fLastThreeChan //------------------ double TLondon1D3L::operator()(double t, const vector &par) const { + + assert(par.size() == 13); + if(t<0.0) return 0.0; @@ -728,6 +764,9 @@ TLondon1D3LS::TLondon1D3LS() : fCalcNeeded(true), fFirstCall(true), fLastThreeCh //------------------ double TLondon1D3LS::operator()(double t, const vector &par) const { + + assert(par.size() == 12); + if(t<0.0) return 0.0; @@ -814,3 +853,337 @@ double TLondon1D3LS::operator()(double t, const vector &par) const { return fPofT->Eval(t); } + +//------------------ +// Constructor of the TLondon1D4L class -- reading available implantation profiles and +// creates (a pointer to) the TPofTCalc object (with the FFT plan) +//------------------ + +TLondon1D4L::TLondon1D4L() : fCalcNeeded(true), fFirstCall(true), fLastFourChanged(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; + } + + fNSteps = startupHandler->GetNSteps(); + fWisdom = startupHandler->GetWisdomFile(); + string rge_path(startupHandler->GetDataPath()); + vector energy_vec(startupHandler->GetEnergyList()); + + 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); + + TTrimSPData *x = new TTrimSPData(rge_path, energy_vec); + fImpProfile = x; + x = 0; + delete x; + + TPofTCalc *y = new TPofTCalc(fWisdom, fParForPofT); + fPofT = y; + y = 0; + delete y; + + // clean up + if (saxParser) { + delete saxParser; + saxParser = 0; + } + if (startupHandler) { + delete startupHandler; + startupHandler = 0; + } +} + +//------------------ +// TLondon1D4L-Method that calls the procedures to create B(z), 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 TLondon1D4L +//------------------ + +double TLondon1D4L::operator()(double t, const vector &par) const { + + assert(par.size() == 16); + + if(t<0.0) + return 0.0; + + // check if the function is called the first time and if yes, read in parameters + + if(fFirstCall){ + fPar = par; + +/* for (unsigned int i(0); iWeightLayers(par[1], interfaces, weights); + } + + TLondon1D_4L BofZ4(fNSteps, fParForBofZ); + TPofBCalc PofB4(BofZ4, *fImpProfile, fParForPofB); + fPofT->DoFFT(PofB4); + + }/* else { + cout << "Only the phase parameter has changed, (re-)calculating P(t) now..." << endl; + }*/ + + fPofT->CalcPol(fParForPofT); + + fCalcNeeded = false; + fLastFourChanged = false; + } + + return fPofT->Eval(t); + +} + +//------------------ +// Constructor of the TLondon1D3LSub class -- reading available implantation profiles and +// creates (a pointer to) the TPofTCalc object (with the FFT plan) +//------------------ + +TLondon1D3LSub::TLondon1D3LSub() : fCalcNeeded(true), fFirstCall(true), fWeightsChanged(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; + } + + fNSteps = startupHandler->GetNSteps(); + fWisdom = startupHandler->GetWisdomFile(); + string rge_path(startupHandler->GetDataPath()); + vector energy_vec(startupHandler->GetEnergyList()); + + 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); + + TTrimSPData *x = new TTrimSPData(rge_path, energy_vec); + fImpProfile = x; + x = 0; + delete x; + + TPofTCalc *y = new TPofTCalc(fWisdom, fParForPofT); + fPofT = y; + y = 0; + delete y; + + // clean up + if (saxParser) { + delete saxParser; + saxParser = 0; + } + if (startupHandler) { + delete startupHandler; + startupHandler = 0; + } +} + +//------------------ +// TLondon1D3L-Method that calls the procedures to create B(z), 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 TLondon1D3L +//------------------ + +double TLondon1D3LSub::operator()(double t, const vector &par) const { + + assert(par.size() == 15); + + if(t<0.0) + return 0.0; + + // check if the function is called the first time and if yes, read in parameters + + if(fFirstCall){ + fPar = par; + +/* for (unsigned int i(0); iWeightLayers(par[1], interfaces, weights); + } + + TLondon1D_3L BofZ3(fNSteps, fParForBofZ); + TPofBCalc PofB3(BofZ3, *fImpProfile, fParForPofB); + + // Add background contribution from the substrate + PofB3.AddBackground(par[2], par[14], fImpProfile->LayerFraction(par[1], 4, interfaces)); + + // FourierTransform of P(B) + fPofT->DoFFT(PofB3); + + }/* else { + cout << "Only the phase parameter has changed, (re-)calculating P(t) now..." << endl; + }*/ + + fPofT->CalcPol(fParForPofT); + + fCalcNeeded = false; + fWeightsChanged = false; + } + + return fPofT->Eval(t); + +} + diff --git a/src/external/TFitPofB-lib/classes/TPofBCalc.cpp b/src/external/TFitPofB-lib/classes/TPofBCalc.cpp index 0f63063c..1739fdd4 100644 --- a/src/external/TFitPofB-lib/classes/TPofBCalc.cpp +++ b/src/external/TFitPofB-lib/classes/TPofBCalc.cpp @@ -193,6 +193,42 @@ TPofBCalc::TPofBCalc( const TBofZCalc &BofZ, const TTrimSPData &dataTrimSP, cons } +// TPofBCalc::AddBackground, Parameters: field B[G], width s[us], weight w + +void TPofBCalc::AddBackground(double B, double s, double w) { + + if(!s || w<0.0 || w>1.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)); + } + + // normalize background + double bgsum(0.0); + for (unsigned int i(0); i < sizePB; i++) + bgsum += bg[i]; + bgsum *= fDB; + for (unsigned int i(0); i < sizePB; i++) + bg[i] /= bgsum; + + // add background to P(B) + for (unsigned int i(0); i < sizePB; i++) + fPB[i] = (1.0 - w)*fPB[i] + w*bg[i]; + +// // check if normalization is still valid +// double pBsum(0.0); +// for (unsigned int i(0); i < sizePB; i++) +// pBsum += fPB[i]; +// +// cout << "pBsum = " << pBsum << endl; + +} void TPofBCalc::ConvolveGss(double w) { diff --git a/src/external/TFitPofB-lib/classes/TTrimSPDataHandler.cpp b/src/external/TFitPofB-lib/classes/TTrimSPDataHandler.cpp index 42122ca8..86e4db89 100644 --- a/src/external/TFitPofB-lib/classes/TTrimSPDataHandler.cpp +++ b/src/external/TFitPofB-lib/classes/TTrimSPDataHandler.cpp @@ -127,6 +127,52 @@ vector TTrimSPData::OrigDataNZ(double e) const { } +//--------------------- +// Method returning fraction of muons implanted in the specified layer for a given energy[keV] +// Parameters: Energy[keV], LayerNumber[1], Interfaces[nm] +//--------------------- + +double TTrimSPData::LayerFraction(double e, unsigned int layno, const vector& interface) const { + + if(layno < 1 && layno > (interface.size()+1)) { + cout << "No such layer available according to your specified interfaces... Returning 0.0!" << endl; + return 0.0; + } + + for(unsigned int i(0); i= *(interface.end()-1)*10.0) + layerNumber += fDataNZ[i][j]; + } else { + for(unsigned int j(0); j= interface[layno-2]*10.0 && fDataZ[i][j] < interface[layno-1]*10.0) + layerNumber += fDataNZ[i][j]; + } + // fraction of muons in layer layno + cout << "Fraction of muons in layer " << layno << ": " << layerNumber/totalNumber << endl; + return layerNumber/totalNumber; + } + } + + // default + cout << "No implantation profile available for the specified energy... Returning 0.0" << endl; + return 0.0; + +} + //--------------------- // Method putting different weight to different layers of your thin film // Parameters: Implantation Energy[keV], Interfaces[nm], Weights [0.0 <= w[i] <= 1.0]