LMU/musrfit
LMU/musrfit
5
0
Fork 0
Code Issues 6 Pull Requests Actions Packages Projects Releases 3 Activity

added a magnetic proximity user function to the ASlibs

Browse Source
This commit is contained in:
nemu 2011-03-10 11:45:02 +00:00
parent 34e462e05f
commit d634a9286c
15 changed files with 1283 additions and 3 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
ChangeLog
View File

@ -6,6 +6,7 @@
changes since 0.8.0
===================================
NEW added a magnetic proximity user function to the ASlibs
NEW added the command SCALE_N0_BKG TRUE | FALSE to the command-block. This can be used to force a single histogram fit
to use either 1/ns scaling for N0 and background or 1/bins one.
NEW any2many: some more work, including the PSI-BIN write routines which are officially not released yet.

3
configure.ac
View File

@ -622,6 +622,7 @@ AC_CONFIG_FILES([Makefile \
src/external/libLFRelaxation/Makefile \
src/external/libGapIntegrals/Makefile \
src/external/libCalcMeanFieldsLEM/Makefile \
src/external/Nonlocal/Makefile])
src/external/Nonlocal/Makefile \
src/external/MagProximity/Makefile])
AC_OUTPUT

64
src/external/MagProximity/Makefile.am vendored Normal file
View File

@ -0,0 +1,64 @@
## Process this file with automake to create Makefile.in
h_sources = \
PMagProximityFitter.h \
PMPRgeHandler.h \
PMPStartupHandler.h \
PMagProximity.h
h_linkdef = \
PMagProximityFitterLinkDef.h \
PMPStartupHandlerLinkDef.h
dict_h_sources = \
PMagProximityFitterDict.h \
PMPStartupHandlerDict.h
cpp_sources = \
PMagProximityFitter.cpp \
PMPRgeHandler.cpp \
PMPStartupHandler.cpp
dict_cpp_sources = \
PMagProximityFitterDict.cpp \
PMPStartupHandlerDict.cpp
include_HEADERS = $(h_sources)
noinst_HEADERS = $(h_linkdef) $(dict_h_sources)
INCLUDES = -I$(top_srcdir)/src/include -I../include $(PMUSR_CFLAGS) $(FFTW3_CFLAGS) $(ROOT_CFLAGS)
AM_CXXFLAGS = $(LOCAL_LIB_CXXFLAGS)
BUILT_SOURCES = $(dict_cpp_sources) $(dict_h_sources)
AM_LDFLAGS = $(LOCAL_LIB_LDFLAGS) -L@ROOTLIBDIR@
CLEANFILES = *Dict.cpp *Dict.h *~ core
%Dict.cpp %Dict.h: %.h %LinkDef.h
@ROOTCINT@ -v -f $*Dict.cpp -c -p $(INCLUDES) $^
lib_LTLIBRARIES = libPMagProximityFitter.la
libPMagProximityFitter_la_SOURCES = $(h_sources) $(cpp_sources) $(dict_h_sources) $(dict_cpp_sources)
libPMagProximityFitter_la_LIBADD = $(PMUSR_LIBS) $(FFTW3_LIBS) $(ROOT_LIBS)
libPMagProximityFitter_la_LDFLAGS = -version-info $(PLUGIN_LIBRARY_VERSION) -release $(PLUGIN_RELEASE) $(AM_LDFLAGS)
## For the moment do not build pkgconfig files for musrfit plug-ins...
## pkgconfigdir = $(libdir)/pkgconfig
## pkgconfig_DATA = PMagProximityFitter.pc
## However, create some symbolic links to the shared library
## in order to unify the function call on different operating systems
if IS_DARWIN
install-exec-hook:
$(LN_S) $(libdir)/libPMagProximityFitter.dylib $(libdir)/libPMagProximityFitter.so
uninstall-hook:
rm -f $(libdir)/libPMagProximityFitter.so
endif
if IS_CYGWIN
install-exec-hook:
$(LN_S) $(bindir)/cygPMagProximityFitter-$(PLUGIN_MAJOR_VERSION)-$(PLUGIN_MINOR_VERSION)-$(PLUGIN_MAJOR_VERSION).dll $(libdir)/libPMagProximityFitter.so
uninstall-hook:
rm -f $(libdir)/libPMagProximityFitter.so
endif

215
src/external/MagProximity/PMPRgeHandler.cpp vendored Normal file
View File

@ -0,0 +1,215 @@
/***************************************************************************
PMPRgeHandler.cpp
Author: Andreas Suter
e-mail: andreas.suter@psi.ch
$Id$
***************************************************************************/
/***************************************************************************
* Copyright (C) 2011 by Andreas Suter *
* andreas.suter@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 <iostream>
#include <fstream>
using namespace std;
#include "PMPRgeHandler.h"
//--------------------------------------------------------------------------
// Constructor
//--------------------------------------------------------------------------
/**
*
*/
PMPRgeHandler::PMPRgeHandler(const PStringVector &rgeDataPathList, const PDoubleVector &rgeDataEnergyList)
{
fIsValid = false;
fIsValid = LoadRgeData(rgeDataPathList, rgeDataEnergyList);
}
//--------------------------------------------------------------------------
// Destructor
//--------------------------------------------------------------------------
/**
*
*/
PMPRgeHandler::~PMPRgeHandler()
{
fRgeDataList.clear();
}
//--------------------------------------------------------------------------
// GetRgeEnergyIndex
//--------------------------------------------------------------------------
/**
*
* \param energy in (keV)
*/
Int_t PMPRgeHandler::GetRgeEnergyIndex(const Double_t energy)
{
Int_t idx = -1;
for (UInt_t i=0; i<fRgeDataList.size(); i++) {
if (energy == fRgeDataList[i].energy) {
idx = i;
break;
}
}
return idx;
}
//--------------------------------------------------------------------------
// GetRgeValue
//--------------------------------------------------------------------------
/**
*
* \param energy in (keV)
* \param dist in (nm)
*/
Double_t PMPRgeHandler::GetRgeValue(const Int_t index, const Double_t dist)
{
Double_t rgeVal = 0.0;
UInt_t distIdx = (UInt_t)(dist/(fRgeDataList[index].stoppingDistance[1]-fRgeDataList[index].stoppingDistance[0]));
if ((distIdx >= fRgeDataList[index].stoppingDistance.size()) || (distIdx < 0)) {
rgeVal = 0.0;
} else {
rgeVal = fRgeDataList[index].stoppingAmplitude[distIdx] +
(fRgeDataList[index].stoppingAmplitude[distIdx+1] - fRgeDataList[index].stoppingAmplitude[distIdx]) *
(dist-fRgeDataList[index].stoppingDistance[distIdx])/(fRgeDataList[index].stoppingDistance[distIdx+1]-fRgeDataList[index].stoppingDistance[distIdx]);
}
return rgeVal;
}
//--------------------------------------------------------------------------
// GetRgeValue
//--------------------------------------------------------------------------
/**
*
* \param energy in (keV)
* \param dist in (nm)
*/
Double_t PMPRgeHandler::GetRgeValue(const Double_t energy, const Double_t dist)
{
// check if energy is present in rge data list
Int_t idx = -1;
for (UInt_t i=0; i<fRgeDataList.size(); i++) {
if (energy == fRgeDataList[i].energy) {
idx = i;
break;
}
}
// energy not found
if (idx == -1)
return -1.0;
return GetRgeValue(idx, dist);
}
//--------------------------------------------------------------------------
// LoadRgeData
//--------------------------------------------------------------------------
/**
*
*/
Bool_t PMPRgeHandler::LoadRgeData(const PStringVector &rgeDataPathList, const PDoubleVector &rgeDataEnergyList)
{
ifstream fin;
PMPRgeData data;
Int_t idx=0;
TString dataName, tstr;
char line[512];
int result;
double dist, val;
for (UInt_t i=0; i<rgeDataPathList.size(); i++) {
// open rge-file for reading
fin.open(rgeDataPathList[i].Data(), iostream::in);
if (!fin.is_open()) {
cout << endl << "PMPRgeHandler::LoadRgeData **ERROR**";
cout << endl << " Could not open file " << rgeDataPathList[i].Data();
cout << endl;
return false;
}
// keep energy (in keV)
data.energy = rgeDataEnergyList[i]/1000.0;
// read msr-file
idx = 0;
while (!fin.eof()) {
// read a line
fin.getline(line, sizeof(line));
idx++;
// ignore first line
if (idx <= 1)
continue;
// ignore empty lines
tstr = line;
if (tstr.IsWhitespace())
continue;
// get values
result = sscanf(line, "%lf %lf", &dist, &val);
// check if data are valid
if (result != 2) {
fin.close();
return false;
}
// feed fRgeDataList
data.stoppingDistance.push_back(dist/10.0); // keep distancies in (nm)
data.stoppingAmplitude.push_back(val);
}
// normalize stopping distribution
Double_t norm = 0.0;
for (UInt_t j=0; j<data.stoppingAmplitude.size(); j++)
norm += data.stoppingAmplitude[j];
norm *= (data.stoppingDistance[1] - data.stoppingDistance[0]);
for (UInt_t j=0; j<data.stoppingAmplitude.size(); j++)
data.stoppingAmplitude[j] /= norm;
// keep data
fRgeDataList.push_back(data);
// clean up
data.stoppingAmplitude.clear();
data.stoppingDistance.clear();
fin.close();
}
return true;
}

55
src/external/MagProximity/PMPRgeHandler.h vendored Normal file
View File

@ -0,0 +1,55 @@
/***************************************************************************
PMPRgeHandler.h
Author: Andreas Suter
e-mail: andreas.suter@psi.ch
$Id$
***************************************************************************/
/***************************************************************************
* Copyright (C) 2011 by Andreas Suter *
* andreas.suter@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 _PMPRGEHANDLER_H_
#define _PMPRGEHANDLER_H_
#include "PMagProximity.h"
class PMPRgeHandler
{
public:
PMPRgeHandler(const PStringVector &rgeDataPathList, const PDoubleVector &rgeDataEnergyList);
virtual ~PMPRgeHandler();
virtual Bool_t IsValid() { return fIsValid; }
virtual Int_t GetRgeEnergyIndex(const Double_t energy);
virtual Double_t GetRgeValue(const Double_t energy, const Double_t dist);
virtual Double_t GetRgeValue(const Int_t index, const Double_t dist);
private:
Bool_t fIsValid;
PMPRgeDataList fRgeDataList;
virtual Bool_t LoadRgeData(const PStringVector &rgeDataPathList, const PDoubleVector &rgeDataEnergyList);
};
#endif // _PMPRGEHANDLER_H_

268
src/external/MagProximity/PMPStartupHandler.cpp vendored Normal file
View File

@ -0,0 +1,268 @@
/***************************************************************************
PMPStartupHandler.cpp
Author: Andreas Suter
e-mail: andreas.suter@psi.ch
$Id$
***************************************************************************/
/***************************************************************************
* Copyright (C) 2011 by Andreas Suter *
* andreas.suter@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 <cmath>
#include <iostream>
#include <fstream>
using namespace std;
#include "PMPStartupHandler.h"
ClassImpQ(PMPStartupHandler)
//--------------------------------------------------------------------------
// Constructor
//--------------------------------------------------------------------------
/**
*
*/
PMPStartupHandler::PMPStartupHandler()
{
fIsValid = true;
fStartupFileFound = false;
fStartupFilePath = "";
fTrimSpDataPath = TString("");
// get default path (for the moment only linux like)
char startup_path_name[128];
// check if the startup file is found in the current directory
strcpy(startup_path_name, "./mag_proximity_startup.xml");
if (StartupFileExists(startup_path_name)) {
fStartupFileFound = true;
fStartupFilePath = TString(startup_path_name);
} else { // startup file is not found in the current directory
cout << endl << "PMPStartupHandler(): **WARNING** Couldn't find mag_proximity_startup.xml in the current directory, will try default one." << endl;
strncpy(startup_path_name, "/home/nemu/analysis/musrfit/src/external/MagProximity/mag_proximity_startup.xml", sizeof(startup_path_name));
if (StartupFileExists(startup_path_name)) {
fStartupFileFound = true;
fStartupFilePath = TString(startup_path_name);
}
}
}
//--------------------------------------------------------------------------
// Destructor
//--------------------------------------------------------------------------
/**
*
*/
PMPStartupHandler::~PMPStartupHandler()
{
fTrimSpDataPathList.clear();
fTrimSpDataEnergyList.clear();
}
//--------------------------------------------------------------------------
// OnStartDocument
//--------------------------------------------------------------------------
/**
* <p>
*/
void PMPStartupHandler::OnStartDocument()
{
fKey = eEmpty;
}
//--------------------------------------------------------------------------
// OnEndDocument
//--------------------------------------------------------------------------
/**
* <p>
*/
void PMPStartupHandler::OnEndDocument()
{
// nothing to be done for now
}
//--------------------------------------------------------------------------
// OnStartElement
//--------------------------------------------------------------------------
/**
* <p>
*
* \param str
* \param attributes
*/
void PMPStartupHandler::OnStartElement(const char *str, const TList *attributes)
{
if (!strcmp(str, "data_path")) {
fKey = eDataPath;
} else if (!strcmp(str, "energy")) {
fKey = eEnergy;
}
}
//--------------------------------------------------------------------------
// OnEndElement
//--------------------------------------------------------------------------
/**
* <p>
*
* \param str
*/
void PMPStartupHandler::OnEndElement(const char *str)
{
fKey = eEmpty;
}
//--------------------------------------------------------------------------
// OnCharacters
//--------------------------------------------------------------------------
/**
* <p>
*
* \param str
*/
void PMPStartupHandler::OnCharacters(const char *str)
{
TString tstr;
switch (fKey) {
case eDataPath:
fTrimSpDataPath = str;
break;
case eEnergy:
tstr = str;
if (tstr.IsFloat()) {
fTrimSpDataEnergyList.push_back(tstr.Atof());
tstr = fTrimSpDataPath;
tstr += str;
tstr += ".rge";
fTrimSpDataPathList.push_back(tstr);
} else {
cout << endl << "PMPStartupHandler::OnCharacters: **ERROR** when finding energy:";
cout << endl << "\"" << str << "\" is not a floating point number, will ignore it and use the default value.";
cout << endl;
}
break;
default:
break;
}
}
//--------------------------------------------------------------------------
// OnComment
//--------------------------------------------------------------------------
/**
* <p>
*
* \param str
*/
void PMPStartupHandler::OnComment(const char *str)
{
// nothing to be done for now
}
//--------------------------------------------------------------------------
// OnWarning
//--------------------------------------------------------------------------
/**
* <p>
*
* \param str
*/
void PMPStartupHandler::OnWarning(const char *str)
{
cout << endl << "PMPStartupHandler **WARNING** " << str;
cout << endl;
}
//--------------------------------------------------------------------------
// OnError
//--------------------------------------------------------------------------
/**
* <p>
*
* \param str
*/
void PMPStartupHandler::OnError(const char *str)
{
cout << endl << "PMPStartupHandler **ERROR** " << str;
cout << endl;
}
//--------------------------------------------------------------------------
// OnFatalError
//--------------------------------------------------------------------------
/**
* <p>
*
* \param str
*/
void PMPStartupHandler::OnFatalError(const char *str)
{
cout << endl << "PMPStartupHandler **FATAL ERROR** " << str;
cout << endl;
}
//--------------------------------------------------------------------------
// OnCdataBlock
//--------------------------------------------------------------------------
/**
* <p>
*
* \param str
*/
void PMPStartupHandler::OnCdataBlock(const char *str, Int_t len)
{
// nothing to be done for now
}
//--------------------------------------------------------------------------
// StartupFileExists
//--------------------------------------------------------------------------
/**
* <p>
*
*/
bool PMPStartupHandler::StartupFileExists(char *fln)
{
bool result = false;
ifstream ifile(fln);
if (ifile.fail()) {
result = false;
} else {
result = true;
ifile.close();
}
return result;
}
// -------------------------------------------------------------------------
// end
// -------------------------------------------------------------------------

83
src/external/MagProximity/PMPStartupHandler.h vendored Normal file
View File

@ -0,0 +1,83 @@
/***************************************************************************
PMPStartupHandler.h
Author: Andreas Suter
e-mail: andreas.suter@psi.ch
$Id$
***************************************************************************/
/***************************************************************************
* Copyright (C) 2011 by Andreas Suter *
* andreas.suter@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 _PMPSTARTUPHANDLER_H_
#define _PMPSTARTUPHANDLER_H_
#include <TObject.h>
#include <TQObject.h>
#include <TString.h>
#include "PMagProximity.h"
class PMPStartupHandler : public TObject
{
public:
PMPStartupHandler();
virtual ~PMPStartupHandler();
virtual void OnStartDocument(); // SLOT
virtual void OnEndDocument(); // SLOT
virtual void OnStartElement(const char*, const TList*); // SLOT
virtual void OnEndElement(const char*); // SLOT
virtual void OnCharacters(const char*); // SLOT
virtual void OnComment(const char*); // SLOT
virtual void OnWarning(const char*); // SLOT
virtual void OnError(const char*); // SLOT
virtual void OnFatalError(const char*); // SLOT
virtual void OnCdataBlock(const char*, Int_t); // SLOT
virtual bool IsValid() { return fIsValid; }
virtual TString GetStartupFilePath() { return fStartupFilePath; }
virtual const PStringVector GetTrimSpDataPathList() const { return fTrimSpDataPathList; }
virtual const PDoubleVector GetTrimSpDataVectorList() const { return fTrimSpDataEnergyList; }
virtual bool StartupFileFound() { return fStartupFileFound; }
private:
enum EKeyWords {eEmpty, eComment, eDataPath, eEnergy};
EKeyWords fKey;
bool fIsValid;
bool fStartupFileFound;
TString fStartupFilePath;
TString fTrimSpDataPath;
PStringVector fTrimSpDataPathList;
PDoubleVector fTrimSpDataEnergyList;
bool StartupFileExists(char *fln);
ClassDef(PMPStartupHandler, 1)
};
#endif // _PMPSTARTUPHANDLER_H_

40
src/external/MagProximity/PMPStartupHandlerLinkDef.h vendored Normal file
View File

@ -0,0 +1,40 @@
/***************************************************************************
PMPStartupHandlerLinkDef.h
Author: Andreas Suter
e-mail: andreas.suter@psi.ch
$Id$
***************************************************************************/
/***************************************************************************
* Copyright (C) 2011 by Andreas Suter *
* andreas.suter@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. *
***************************************************************************/
#ifdef __CINT__
#pragma link off all globals;
#pragma link off all classes;
#pragma link off all functions;
#pragma link C++ class PMPStartupHandler+;
#endif

67
src/external/MagProximity/PMagProximity.h vendored Normal file
View File

@ -0,0 +1,67 @@
/***************************************************************************
PMagProximity.h
Author: Andreas Suter
e-mail: andreas.suter@psi.ch
$Id$
***************************************************************************/
/***************************************************************************
* Copyright (C) 2011 by Andreas Suter *
* andreas.suter@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 _PMAGPROXIMITY_H_
#define _PMAGPROXIMITY_H_
#include <vector>
#include <TString.h>
//-------------------------------------------------------------
/**
* <p> typedef to make to code more readable.
*/
typedef std::vector<TString> PStringVector;
//-------------------------------------------------------------
/**
* <p> typedef to make to code more readable.
*/
typedef std::vector<Double_t> PDoubleVector;
//-------------------------------------------------------------
/**
* <p>
*/
typedef struct {
Double_t energy;
PDoubleVector stoppingDistance;
PDoubleVector stoppingAmplitude;
} PMPRgeData;
//-------------------------------------------------------------
/**
* <p>
*/
typedef std::vector<PMPRgeData> PMPRgeDataList;
#endif // _PMAGPROXIMITY_H_

340
src/external/MagProximity/PMagProximityFitter.cpp vendored Normal file
View File

@ -0,0 +1,340 @@
/***************************************************************************
PMagProximityFitter.cpp
Author: Andreas Suter
e-mail: andreas.suter@psi.ch
$Id$
***************************************************************************/
/***************************************************************************
* Copyright (C) 2011 by Andreas Suter *
* andreas.suter@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 <cmath>
#include <iostream>
using namespace std;
#include <TSAXParser.h>
#include <TMath.h>
#include "PMagProximityFitter.h"
#define GAMMA_MU 0.0851615503527
#define DEGREE2RAD 0.0174532925199
ClassImp(PMagProximityFitterGlobal)
//--------------------------------------------------------------------------
// Constructor
//--------------------------------------------------------------------------
/**
*
*/
PMagProximityFitterGlobal::PMagProximityFitterGlobal()
{
fValid = true;
fStartupHandler = 0;
fRgeHandler = 0;
fDz = 0.01; // 0.01nm
// read XML startup file
char startup_path_name[128];
TSAXParser *saxParser = new TSAXParser();
PMPStartupHandler *fStartupHandler = new PMPStartupHandler();
strcpy(startup_path_name, fStartupHandler->GetStartupFilePath().Data());
saxParser->ConnectToHandler("PMPStartupHandler", fStartupHandler);
Int_t status = saxParser->ParseFile(startup_path_name);
// check for parse errors
if (status) { // error
cout << endl << ">> PMagProximityFitterGlobal::PMagProximityFitterGlobal: **WARNING** reading/parsing mag_proximity_startup.xml.";
cout << endl;
fValid = false;
}
// clean up
if (saxParser) {
delete saxParser;
saxParser = 0;
}
// check if everything went fine with the startup handler
if (!fStartupHandler->IsValid()) {
cout << endl << ">> PMagProximityFitterGlobal::PMagProximityFitterGlobal **PANIC ERROR**";
cout << endl << ">> startup handler too unhappy. Will terminate unfriendly, sorry.";
cout << endl;
fValid = false;
}
// load all the TRIM.SP rge-files
fRgeHandler = new PMPRgeHandler(fStartupHandler->GetTrimSpDataPathList(), fStartupHandler->GetTrimSpDataVectorList());
if (!fRgeHandler->IsValid()) {
cout << endl << ">> PMagProximityFitterGlobal::PMagProximityFitterGlobal **PANIC ERROR**";
cout << endl << ">> rge data handler too unhappy. Will terminate unfriendly, sorry.";
cout << endl;
fValid = false;
}
}
//--------------------------------------------------------------------------
// Destructor
//--------------------------------------------------------------------------
/**
*
*/
PMagProximityFitterGlobal::~PMagProximityFitterGlobal()
{
fPreviousParam.clear();
fField.clear();
if (fRgeHandler) {
delete fRgeHandler;
fRgeHandler = 0;
}
if (fStartupHandler) {
delete fStartupHandler;
fStartupHandler = 0;
}
}
//--------------------------------------------------------------------------
// CalculateField (public)
//--------------------------------------------------------------------------
/**
*
*/
void PMagProximityFitterGlobal::CalculateField(const std::vector<Double_t> &param) const
{
// param: [0] energy (keV), [1] zStart (nm), [2] zEnd (nm), [3] B0 (G), [4] B1 (G), [5] zeta (nm), [6] phase (°)
assert(param.size() == 7);
// check that param are new and hence a calculation is needed (except the phase, which is not needed here)
Bool_t newParams = false;
if (fPreviousParam.size() == 0) {
for (UInt_t i=0; i<param.size(); i++)
fPreviousParam.push_back(param[i]);
newParams = true;
} else {
assert(param.size() == fPreviousParam.size());
for (UInt_t i=0; i<param.size(); i++) {
if (i == 6) // phase, nothing to be done
continue;
if (param[i] != fPreviousParam[i]) {
newParams = true;
break;
}
}
}
if (!newParams)
return;
// keep parameters
for (UInt_t i=0; i<param.size(); i++)
fPreviousParam[i] = param[i];
// delete previous field vector
fField.clear();
// calculate the magnetic field with a 0.1nm step resolution
Double_t z = param[1];
Double_t dval = 0.0;
do {
if (param[5] == 0.0) { // zeta == 0
dval = 0.0;
} else {
dval = param[4] * exp(-(param[2]-z)/param[5]) + param[3];
}
fField.push_back(dval);
z += fDz;
} while (z<=param[2]);
}
//--------------------------------------------------------------------------
// GetMagneticField
//--------------------------------------------------------------------------
/**
*
*/
Double_t PMagProximityFitterGlobal::GetMagneticField(const Double_t z) const
{
Double_t result = -1.0;
// check if z is out of range [zStart, zEnd]
UInt_t idx=0;
if ((z<fPreviousParam[1]) || (z>fPreviousParam[2])) {
result = 0.0;
} else {
idx = (UInt_t)((z-fPreviousParam[1])/fDz);
if (idx == fField.size()-1) {
result = fField[idx];
} else {
result = fField[idx] + (fField[idx+1]-fField[idx])*((z-fPreviousParam[1])-idx*fDz)/fDz;
}
}
/*
static UInt_t count=0;
if (count < 10) {
count++;
if (idx == fField.size()-1)
cout << endl << "debug> z=" << z << ", idx=" << idx << ", fDz=" << fDz << ", fField[idx]=" << fField[idx] << ", result=" << result;
else
cout << endl << "debug> z=" << z << ", idx=" << idx << ", fDz=" << fDz << ", fField[idx]=" << fField[idx] << ", fField[idx+1]=" << fField[idx+1] << ", result=" << result;
cout << endl;
}
*/
// cout << endl << z << ", " << result;
return result;
}
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
ClassImp(PMagProximityFitter)
//--------------------------------------------------------------------------
// Constructor
//--------------------------------------------------------------------------
/**
*
*/
PMagProximityFitter::PMagProximityFitter()
{
fValid = false;
fInvokedGlobal = false;
fMagProximityFitterGlobal = 0;
}
//--------------------------------------------------------------------------
// Destructor
//--------------------------------------------------------------------------
/**
*
*/
PMagProximityFitter::~PMagProximityFitter()
{
if ((fMagProximityFitterGlobal != 0) && fInvokedGlobal) {
delete fMagProximityFitterGlobal;
fMagProximityFitterGlobal = 0;
}
}
//--------------------------------------------------------------------------
// SetGlobalPart (public)
//--------------------------------------------------------------------------
/**
* <p>
*
* <b>return:</b>
*
* \param globalPart
* \param idx
*/
void PMagProximityFitter::SetGlobalPart(vector<void*> &globalPart, UInt_t idx)
{
fIdxGlobal = static_cast<Int_t>(idx);
if ((Int_t)globalPart.size() <= fIdxGlobal) {
fMagProximityFitterGlobal = new PMagProximityFitterGlobal();
if (fMagProximityFitterGlobal == 0) {
fValid = false;
cerr << endl << ">> PMagProximityFitter::SetGlobalPart(): **ERROR** Couldn't invoke global user function object, sorry ..." << endl;
} else if (!fMagProximityFitterGlobal->IsValid()) {
fValid = false;
cerr << endl << ">> PMagProximityFitter::SetGlobalPart(): **ERROR** initialization of global user function object failed, sorry ..." << endl;
} else {
fValid = true;
fInvokedGlobal = true;
globalPart.resize(fIdxGlobal+1);
globalPart[fIdxGlobal] = dynamic_cast<PMagProximityFitterGlobal*>(fMagProximityFitterGlobal);
}
} else {
fValid = true;
fMagProximityFitterGlobal = (PMagProximityFitterGlobal*)globalPart[fIdxGlobal];
}
}
//--------------------------------------------------------------------------
// GlobalPartIsValid (public)
//--------------------------------------------------------------------------
/**
* <p>
*
* <b>return:</b>
*/
Bool_t PMagProximityFitter::GlobalPartIsValid() const
{
return (fValid && fMagProximityFitterGlobal->IsValid());
}
//--------------------------------------------------------------------------
// operator()
//--------------------------------------------------------------------------
/**
*
*/
Double_t PMagProximityFitter::operator()(Double_t t, const std::vector<Double_t> &param) const
{
// param: [0] energy (keV), [1] zStart (nm), [2] zEnd (nm), [3] B0 (G), [4] B1 (G), [5] zeta (nm), [6] phase (°)
assert(param.size() == 7);
// for negative time return polarization == 1
if (t <= 0.0)
return 1.0;
// if zeta=param[5] has a negative value, return 1.0. This can be used to switch off the proximity effect but still being able to fit asymmetries etc.
if (param[5] <= 0.0)
return 1.0;
// calculate field if parameter have changed
fMagProximityFitterGlobal->CalculateField(param);
Int_t energyIndex = fMagProximityFitterGlobal->GetEnergyIndex(param[0]);
// calculate polarization
Bool_t done = false;
Double_t pol = 0.0, dPol = 0.0;
Double_t z=param[1]; // set to zStart
Int_t terminate = 0;
Double_t dz = 0.1; // setp width == 0.1nm
do {
dPol = fMagProximityFitterGlobal->GetMuonStoppingDensity(energyIndex, z) * cos(GAMMA_MU * fMagProximityFitterGlobal->GetMagneticField(z) * t + param[6] * DEGREE2RAD);
z += dz;
pol += dPol;
// change in polarization is very small hence start termination counting
if (fabs(dPol) < 1.0e-5) {
terminate++;
} else {
terminate = 0;
}
if (terminate > 10) // polarization died out hence one can stop
done = true;
} while (!done);
return pol*dz;
}

89
src/external/MagProximity/PMagProximityFitter.h vendored Normal file
View File

@ -0,0 +1,89 @@
/***************************************************************************
PMagProximity.h
Author: Andreas Suter
e-mail: andreas.suter@psi.ch
$Id$
***************************************************************************/
/***************************************************************************
* Copyright (C) 2009 by Andreas Suter *
* andreas.suter@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 _PMAGPROXIMITYFITTER_H_
#define _PMAGPROXIMITYFITTER_H_
#include "PUserFcnBase.h"
#include "PMPStartupHandler.h"
#include "PMPRgeHandler.h"
class PMagProximityFitterGlobal
{
public:
PMagProximityFitterGlobal();
virtual ~PMagProximityFitterGlobal();
Bool_t IsValid() { return fValid; }
virtual void CalculateField(const std::vector<Double_t> &param) const;
virtual Int_t GetEnergyIndex(const Double_t energy) { return fRgeHandler->GetRgeEnergyIndex(energy); }
virtual Double_t GetMuonStoppingDensity(const Int_t energyIndex, const Double_t z) const { return fRgeHandler->GetRgeValue(energyIndex, z); }
virtual Double_t GetMagneticField(const Double_t z) const;
private:
Bool_t fValid;
PMPStartupHandler *fStartupHandler;
PMPRgeHandler *fRgeHandler;
mutable std::vector<Double_t> fPreviousParam;
mutable Int_t fEnergyIndex; // keeps the proper index to select n(z)
mutable Double_t fDz;
mutable std::vector<Double_t> fField;
ClassDef(PMagProximityFitterGlobal, 1)
};
class PMagProximityFitter : public PUserFcnBase
{
public:
PMagProximityFitter();
virtual ~PMagProximityFitter();
virtual Bool_t NeedGlobalPart() const { return true; }
virtual void SetGlobalPart(vector<void*> &globalPart, UInt_t idx);
virtual Bool_t GlobalPartIsValid() const;
virtual Double_t operator()(Double_t t, const std::vector<Double_t> &param) const;
private:
Bool_t fValid;
Bool_t fInvokedGlobal;
Int_t fIdxGlobal;
PMagProximityFitterGlobal *fMagProximityFitterGlobal;
ClassDef(PMagProximityFitter, 1)
};
#endif // _PMAGPROXIMITYFITTER_H_

41
src/external/MagProximity/PMagProximityFitterLinkDef.h vendored Normal file
View File

@ -0,0 +1,41 @@
/***************************************************************************
PMagProximityFitterLinkDef.h
Author: Andreas Suter
e-mail: andreas.suter@psi.ch
$Id$
***************************************************************************/
/***************************************************************************
* Copyright (C) 2011 by Andreas Suter *
* andreas.suter@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. *
***************************************************************************/
#ifdef __CINT__
#pragma link off all globals;
#pragma link off all classes;
#pragma link off all functions;
#pragma link C++ class PMagProximityFitterGlobal+;
#pragma link C++ class PMagProximityFitter+;
#endif

15
src/external/MagProximity/mag_proximity_startup.xml vendored Normal file
View File

@ -0,0 +1,15 @@
<?xml version="1.0" encoding="UTF-8"?>
<mag_proximity xmlns="http://nemu.web.psi.ch/musrfit/mag_proximity">
<comment>
$Id$
</comment>
<trim_sp_part>
<data_path>/mnt/home/nemu/analysis/2010/EuS-Co/trimsp2/EuS-Co_E</data_path>
<energy_list>
<energy>5030</energy>
<energy>6330</energy>
<energy>7530</energy>
<energy>8730</energy>
</energy_list>
</trim_sp_part>
</mag_proximity>

3
src/external/Makefile.am vendored
View File

@ -1,5 +1,6 @@
if BUILD_ASLIBS
ASDIRS = Nonlocal
ASDIRS = Nonlocal \
MagProximity
endif
if BUILD_CUBALIB

2
src/external/Nonlocal/PNL_RgeHandler.cpp vendored
View File

@ -98,7 +98,7 @@ Double_t PNL_RgeHandler::GetRgeValue(const Int_t index, const Double_t dist)
} else {
rgeVal = fRgeDataList[index].stoppingAmplitude[distIdx] +
(fRgeDataList[index].stoppingAmplitude[distIdx+1] - fRgeDataList[index].stoppingAmplitude[distIdx]) *
(fRgeDataList[index].stoppingDistance[distIdx+1]-dist)/(fRgeDataList[index].stoppingDistance[distIdx+1]-fRgeDataList[index].stoppingDistance[distIdx]);
(dist-fRgeDataList[index].stoppingDistance[distIdx])/(fRgeDataList[index].stoppingDistance[distIdx+1]-fRgeDataList[index].stoppingDistance[distIdx]);
}
/*