From 121496f5e570d0ebfd8a9204778b7196934c1043 Mon Sep 17 00:00:00 2001 From: nemu Date: Thu, 2 Jul 2009 13:36:18 +0000 Subject: [PATCH] some more work towards nonlocal fitting --- src/external/Nonlocal/PNL_PippardFitter.cpp | 116 +++++++++++++++----- src/external/Nonlocal/PNL_PippardFitter.h | 7 +- src/external/Nonlocal/PNL_RgeHandler.cpp | 19 +++- src/external/Nonlocal/nonlocal_startup.xml | 2 +- 4 files changed, 111 insertions(+), 33 deletions(-) diff --git a/src/external/Nonlocal/PNL_PippardFitter.cpp b/src/external/Nonlocal/PNL_PippardFitter.cpp index bc859bbb..aa13c44a 100644 --- a/src/external/Nonlocal/PNL_PippardFitter.cpp +++ b/src/external/Nonlocal/PNL_PippardFitter.cpp @@ -40,6 +40,9 @@ using namespace std; #include "PNL_PippardFitter.h" +#define GAMMA_MU 0.0851615503527 +#define DEGREE2RAD 0.0174532925199 + ClassImp(PNL_PippardFitter) //-------------------------------------------------------------------------- @@ -87,6 +90,8 @@ PNL_PippardFitter::PNL_PippardFitter() assert(false); } + fFourierPoints = fStartupHandler->GetFourierPoints(); + // load all the TRIM.SP rge-files fRgeHandler = new PNL_RgeHandler(fStartupHandler->GetTrimSpDataPathList()); if (!fRgeHandler->IsValid()) { @@ -112,15 +117,6 @@ PNL_PippardFitter::PNL_PippardFitter() */ PNL_PippardFitter::~PNL_PippardFitter() { - if (fStartupHandler) { - delete fStartupHandler; - fStartupHandler = 0; - } - if (fRgeHandler) { - delete fRgeHandler; - fRgeHandler = 0; - } - fPreviousParam.clear(); if (fPlanPresent) { @@ -135,6 +131,17 @@ PNL_PippardFitter::~PNL_PippardFitter() fftw_free(fFieldq); fFieldB = 0; } + + if (fRgeHandler) { + delete fRgeHandler; + fRgeHandler = 0; + } +/* + if (fStartupHandler) { + delete fStartupHandler; + fStartupHandler = 0; + } +*/ } //-------------------------------------------------------------------------- @@ -145,18 +152,52 @@ PNL_PippardFitter::~PNL_PippardFitter() */ Double_t PNL_PippardFitter::operator()(Double_t t, const std::vector ¶m) const { - // expected parameters: energy, temp, thickness, meanFreePath, xi0, lambdaL, phase - assert(param.size() == 7); + // param: [0] energy, [1] temp, [2] thickness, [3] meanFreePath, [4] xi0, [5] lambdaL, [6] Bext, [7] phase, [8] dead-layer + assert(param.size() == 9); + // for negative time return polarization == 1 + if (t <= 0.0) + return 1.0; + + // calculate field if parameter have changed if (NewParameters(param)) { // new parameters, hence B(z), P(t), ..., needs to be calculated // keep parameters for (UInt_t i=0; iGetRgeEnergyIndex(param[0]); CalculateField(param); - CalculatePolarization(param); } - return 0.0; + // calcualte polarization + Bool_t done = false; + Double_t pol = 0.0, dPol = 0.0; + Double_t z=0.0; + Int_t terminate = 0; + Double_t dz = 1.0; + do { + + if (z < param[8]) { // z < dead-layer + dPol = fRgeHandler->GetRgeValue(fEnergyIndex, z); + } else { + dPol = fRgeHandler->GetRgeValue(fEnergyIndex, z) * cos(GAMMA_MU * param[6] * GetMagneticField(z-param[8]) * t + param[7] * DEGREE2RAD); + } + z += dz; + pol += dPol; + + // change in polarization is very small hence start termination counting + if (fabs(dPol) < 1.0e-7) { + terminate++; + } else { + terminate = 0; + } + + if (terminate > 10) // polarization died out hence one can stop + done = true; + } while (!done); + +// cout << endl << "t = " << t << ", pol = " << pol*dz; + + return pol*dz; } //-------------------------------------------------------------------------- @@ -170,6 +211,7 @@ Bool_t PNL_PippardFitter::NewParameters(const std::vector ¶m) cons if (fPreviousParam.size() == 0) { for (UInt_t i=0; i ¶m) cons */ void PNL_PippardFitter::CalculateField(const std::vector ¶m) const { - // param: [0] energy, [1] temp, [2] thickness, [3] meanFreePath, [4] xi0, [5] lambdaL, [6] phase + // param: [0] energy, [1] temp, [2] thickness, [3] meanFreePath, [4] xi0, [5] lambdaL, [6] Bext, [7] phase, [8] dead-layer - Int_t pippardFourierPoints = fStartupHandler->GetFourierPoints(); +//cout << endl << "in CalculateField ..." << endl; +//cout << endl << "fFourierPoints = " << fFourierPoints; - f_dz = XiP_T(param[4], param[3], param[1])*TMath::TwoPi()/pippardFourierPoints/f_dx; // see lab-book p.137, used for specular reflection boundary conditions (default) + + f_dz = XiP_T(param[4], param[3], param[1])*TMath::TwoPi()/fFourierPoints/f_dx; // see lab-book p.137, used for specular reflection boundary conditions (default) +//cout << endl << "f_dz = " << f_dz; // check if it is necessary to allocate memory if (fFieldq == 0) { - fFieldq = (fftw_complex *) fftw_malloc(sizeof(fftw_complex) * pippardFourierPoints); + fFieldq = (fftw_complex *) fftw_malloc(sizeof(fftw_complex) * fFourierPoints); if (fFieldq == 0) { cout << endl << "PPippard::CalculateField(): **ERROR** couldn't allocate memory for fFieldq"; cout << endl; @@ -210,7 +255,7 @@ void PNL_PippardFitter::CalculateField(const std::vector ¶m) const } } if (fFieldB == 0) { - fFieldB = (fftw_complex *) fftw_malloc(sizeof(fftw_complex) * pippardFourierPoints); + fFieldB = (fftw_complex *) fftw_malloc(sizeof(fftw_complex) * fFourierPoints); if (fFieldB == 0) { cout << endl << "PPippard::CalculateField(): **ERROR** couldn't allocate memory for fFieldB"; cout << endl; @@ -228,7 +273,7 @@ void PNL_PippardFitter::CalculateField(const std::vector ¶m) const Double_t x; fFieldq[0][0] = 0.0; fFieldq[0][1] = 0.0; - for (Int_t i=1; i ¶m) const // Fourier transform if (!fPlanPresent) { - fPlan = fftw_plan_dft_1d(pippardFourierPoints, fFieldq, fFieldB, FFTW_FORWARD, FFTW_EXHAUSTIVE); +// fPlan = fftw_plan_dft_1d(fFourierPoints, fFieldq, fFieldB, FFTW_FORWARD, FFTW_EXHAUSTIVE); + fPlan = fftw_plan_dft_1d(fFourierPoints, fFieldq, fFieldB, FFTW_FORWARD, FFTW_ESTIMATE); fPlanPresent = true; } @@ -245,41 +291,57 @@ void PNL_PippardFitter::CalculateField(const std::vector ¶m) const // normalize fFieldB Double_t norm = 0.0; fShift=0; - for (Int_t i=0; i fabs(norm)) { norm = fFieldB[i][1]; fShift = i; } } - for (Int_t i=0; i ¶m) const +Double_t PNL_PippardFitter::GetMagneticField(const Double_t z) const { + Double_t result = -1.0; + + if (fFieldB == 0) + return -1.0; + + if (z <= 0.0) + return 1.0; + + if (z > f_dz*fFourierPoints/2.0) + return 0.0; + + Int_t bin = (Int_t)(z/f_dz); + + result = fFieldB[bin+fShift][1]; + + return result; } //-------------------------------------------------------------------------- diff --git a/src/external/Nonlocal/PNL_PippardFitter.h b/src/external/Nonlocal/PNL_PippardFitter.h index 73b165c2..947515b7 100644 --- a/src/external/Nonlocal/PNL_PippardFitter.h +++ b/src/external/Nonlocal/PNL_PippardFitter.h @@ -53,17 +53,20 @@ class PNL_PippardFitter : public PUserFcnBase mutable std::vector fPreviousParam; Double_t f_dx; // dx = xiPT dq - mutable Double_t f_dz; // spatial step size + mutable Double_t f_dz; // spatial step size mutable Bool_t fPlanPresent; + mutable Int_t fFourierPoints; mutable fftw_plan fPlan; mutable fftw_complex *fFieldq; // (xiPT x)/(x^2 + xiPT^2 K(x,T)), x = q xiPT mutable fftw_complex *fFieldB; // field calculated for specular boundary conditions mutable Int_t fShift; // shift needed to pick up fFieldB at the maximum for B->0 + mutable Int_t fEnergyIndex; // keeps the proper index to select n(z) + virtual Bool_t NewParameters(const std::vector ¶m) const; virtual void CalculateField(const std::vector ¶m) const; - virtual void CalculatePolarization(const std::vector ¶m) const; + virtual Double_t GetMagneticField(const Double_t z) const; virtual Double_t DeltaBCS(const Double_t t) const; virtual Double_t LambdaL_T(const Double_t lambdaL, const Double_t t) const; diff --git a/src/external/Nonlocal/PNL_RgeHandler.cpp b/src/external/Nonlocal/PNL_RgeHandler.cpp index 60344bcd..25e0bd6c 100644 --- a/src/external/Nonlocal/PNL_RgeHandler.cpp +++ b/src/external/Nonlocal/PNL_RgeHandler.cpp @@ -29,6 +29,8 @@ * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. * ***************************************************************************/ +#include + #include #include using namespace std; @@ -95,10 +97,21 @@ Double_t PNL_RgeHandler::GetRgeValue(const Int_t index, const Double_t dist) rgeVal = 0.0; } else { rgeVal = fRgeDataList[index].stoppingAmplitude[distIdx] + - (fRgeDataList[index].stoppingAmplitude[distIdx] - fRgeDataList[index].stoppingAmplitude[distIdx+1]) * + (fRgeDataList[index].stoppingAmplitude[distIdx+1] - fRgeDataList[index].stoppingAmplitude[distIdx]) * (fRgeDataList[index].stoppingDistance[distIdx+1]-dist)/(fRgeDataList[index].stoppingDistance[distIdx+1]-fRgeDataList[index].stoppingDistance[distIdx]); } +/* +cout << endl << "distIdx = " << distIdx << ", fRgeDataList[index].stoppingDistance.size() = " << fRgeDataList[index].stoppingDistance.size(); +cout << endl << "fRgeDataList[index].stoppingAmplitude[distIdx] = " << fRgeDataList[index].stoppingAmplitude[distIdx]; +cout << endl << "fRgeDataList[index].stoppingAmplitude[distIdx+1] = " << fRgeDataList[index].stoppingAmplitude[distIdx+1]; +cout << endl << "dist = " << dist; +cout << endl << "fRgeDataList[index].stoppingDistance[distIdx] = " << fRgeDataList[index].stoppingDistance[distIdx]; +cout << endl << "fRgeDataList[index].stoppingDistance[distIdx+1] = " << fRgeDataList[index].stoppingDistance[distIdx+1]; +cout << endl << "rgeVal = " << rgeVal; +cout << endl << "----" << endl; +*/ + return rgeVal; } @@ -198,10 +211,10 @@ Bool_t PNL_RgeHandler::LoadRgeData(const PStringVector &rgeDataPathList) // normalize stopping distribution Double_t norm = 0.0; for (UInt_t j=0; j - 200000 + 262144 /afs/psi.ch/project/nemu/analysis/2009/Nonlocal/trimsp/InSne