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

Merged muonspin/musrfit:root6 into master

Browse Source
This commit is contained in:
Zaher Salman 2021-04-06 15:28:19 +02:00
commit 6b32dcaa38
31 changed files with 846 additions and 1435 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

9
CMakeLists.txt
View File

@ -5,7 +5,7 @@ if (CMAKE_VERSION GREATER_EQUAL 3.12)
cmake_policy(SET CMP0075 NEW)
endif (CMAKE_VERSION GREATER_EQUAL 3.12)
project(musrfit VERSION 1.6.5 LANGUAGES C CXX)
project(musrfit VERSION 1.7.0 LANGUAGES C CXX)
#--- musrfit specific options -------------------------------------------------
option(nexus "build optional NeXus support. Needed for ISIS" OFF)
@ -83,7 +83,12 @@ if (ROOT_mathmore_FOUND)
endif (ROOT_mathmore_FOUND)
#--- check for boost ----------------------------------------------------------
find_package(Boost REQUIRED)
find_package(Boost REQUIRED
COMPONENTS
system
filesystem
)
message(STATUS "Boost libs: ${Boost_LIBRARIES}")
#--- check for gsl ------------------------------------------------------------
find_package(GSL REQUIRED)

5
ChangeLog
View File

@ -12,6 +12,11 @@ or
https://bitbucket.org/muonspin/musrfit/commits/all
Release of V1.7.0, 2021/03/30
=============================
Centralize the rge-handling for user functions.
Release of V1.6.5, 2021/01/22
=============================

9
doc/examples/ASlibs/nonlocal_startup.xml
View File

@ -1,13 +1,14 @@
<?xml version="1.0" encoding="UTF-8"?>
<nonlocal xmlns="http://nemu.web.psi.ch/musrfit/nonlocal">
<comment>
$Id: nonlocal_startup.xml 4047 2009-07-02 13:36:18Z nemu $
nonlocal_startup.xml
</comment>
<nonlocal_par>
<fourier_points>262144</fourier_points>
</nonlocal_par>
<trim_sp_part>
<data_path>profiles/Sn_E</data_path>
<trim_sp>
<data_path>./profiles/</data_path>
<rge_fln_pre>Sn_E</rge_fln_pre>
<energy_list>
<energy>1000</energy>
<energy>2000</energy>
@ -22,5 +23,5 @@
<energy>25000</energy>
<energy>27300</energy>
</energy_list>
</trim_sp_part>
</trim_sp>
</nonlocal>

32
doc/examples/ASlibs/test-nonlocal.msr
View File

@ -4,21 +4,21 @@ FITPARAMETER
# Nr. Name Value Step Pos_Error Boundaries
1 one 1 0 none
2 zero 0 0 none
3 Asy 0.2235 -0.0066 0.0072 0 0.33
4 energy 22 0 none
3 Asy 0.2235 -0.0065 0.0072 0 0.33
4 energy 22000 0 none
5 redTemp 0.8683 0 none
6 thickness 5000 0 none
7 ell 12000 0 none
8 xi 94 0 none
9 lambdaL 52.5 -1.4 1.2
9 lambdaL 52.5 -1.3 1.2
10 Bext 47.11 0 none
11 deadLayer 0.00052 -0.00052 6.89563 0 none
11 deadLayer 0.0004 0.0037 none 0 none
12 RateSmear 0.418 -0.015 0.015
13 N0_L 301.61 0.62 none
14 Bkg_L 24.442 0.076 none
13 N0_L 301.6 -1.1 1.2
14 Bkg_L 24.441 -0.092 0.092
15 Phase_L 42.0 -2.2 2.2
16 Alpha_LR 1.0614 -0.0066 0.0064
17 Bkg_R 27.646 -0.090 0.090
16 Alpha_LR 1.0613 -0.0066 0.0064
17 Bkg_R 27.646 -0.089 0.090
18 RelPhase_R 154.4 -1.1 1.1
###############################################################
@ -37,7 +37,6 @@ RUN data/lem10_his_0825 MUE4 PSI ROOT-NPP (name beamline institute data-file-f
fittype 0 (single histogram fit)
norm 13
backgr.fit 14
lifetimecorrection
map 1 2 0 0 0 0 0 0 0 0
forward 1
data 3289 65000
@ -49,7 +48,6 @@ RUN data/lem10_his_0825 MUE4 PSI ROOT-NPP (name beamline institute data-file-f
fittype 0 (single histogram fit)
norm fun1
backgr.fit 17
lifetimecorrection
map 16 18 0 0 0 0 0 0 0 0
forward 3
data 3289 65000
@ -69,23 +67,27 @@ SAVE
###############################################################
PLOT 0 (single histo plot)
lifetimecorrection
runs 1 2
range 0 9 -0.3 0.3
view_packing 500
###############################################################
FOURIER
units Gauss # units either 'Gauss', 'MHz', or 'Mc/s'
units Gauss # units either 'Gauss', 'Tesla', 'MHz', or 'Mc/s'
fourier_power 10
apodization STRONG # NONE, WEAK, MEDIUM, STRONG
plot POWER # REAL, IMAG, REAL_AND_IMAG, POWER, PHASE
plot POWER # REAL, IMAG, REAL_AND_IMAG, POWER, PHASE, PHASE_OPT_REAL
phase 8.5
#range_for_phase_correction 50.0 70.0
range 0.0 200.0
range 0 200
###############################################################
STATISTIC --- 2014-11-05 12:34:54
maxLH = 415.6, NDF = 358, maxLH/NDF = 1.160933
STATISTIC --- 2021-03-30 15:04:30
maxLH = 415.6, NDF = 358, maxLH/NDF = 1.160942
expected maxLH = 418.1, NDF = 358, expected maxLH/NDF = 1.168011
run block 1: (NDF/red.maxLH/red.maxLH_e) = (177/1.072511/1.087460)
run block 2: (NDF/red.maxLH/red.maxLH_e) = (175/1.290186/1.289528)

38
src/classes/CMakeLists.txt
View File

@ -1,4 +1,4 @@
#--- build the PMusr and PUserFcnBase libs ------------------------------------
#--- build the PMusr, PUserFcnBase, and PRgeHandler libs ----------------------
#--- generate necessary dictionaries ------------------------------------------
set(MUSRFIT_INC ${CMAKE_SOURCE_DIR}/src/include)
@ -53,6 +53,18 @@ root_generate_dictionary(
-I${FFTW3_INCLUDE_DIR} -I${MUSRFIT_INC}
MODULE PUserFcnBase
)
root_generate_dictionary(
PRgeHandlerDict
PRgeHandler.h
OPTIONS
-I${NONLOCAL_INC}
-I${MUSRFIT_INC}
-I${CMAKE_CURRENT_SOURCE_DIR}
-inlineInputHeader
LINKDEF ${MUSRFIT_INC}/PRgeHandlerLinkDef.h
MODULE PRgeHandler
)
#--- create pkg-config info ---------------------------------------------------
set(prefix "${CMAKE_INSTALL_PREFIX}")
@ -64,6 +76,8 @@ set(MUSR_LIBRARY_NAME "PMusr")
configure_file("PMusr.pc.in" "PMusr.pc" @ONLY)
set(USERFCN_LIBRARY_NAME "PUserFcnBase")
configure_file("PUserFcnBase.pc.in" "PUserFcnBase.pc" @ONLY)
set(PRGEHANDLER_LIBRARY_NAME "PRgeHandler")
configure_file("PRgeHandler.pc.in" "PRgeHandler.pc" @ONLY)
#--- lib creation -------------------------------------------------------------
add_library(PMusr SHARED
@ -117,6 +131,11 @@ add_library(PUserFcnBase SHARED
PUserFcnBaseDict.cxx
)
add_library(PRgeHandler SHARED
PRgeHandler.cpp
PRgeHandlerDict.cxx
)
#--- set target properties, e.g. version --------------------------------------
set_target_properties(PMusr
PROPERTIES
@ -126,11 +145,20 @@ set_target_properties(PUserFcnBase
PROPERTIES
VERSION ${MUSR_VERSION}
)
set_target_properties(PRgeHandler
PROPERTIES
VERSION ${MUSR_VERSION}
)
#--- make sure that the include directory is found ----------------------------
target_include_directories(
PUserFcnBase BEFORE PRIVATE $<BUILD_INTERFACE:${MUSRFIT_INC}>
)
target_include_directories(
PRgeHandler BEFORE PRIVATE $<BUILD_INTERFACE:${MUSRFIT_INC}>
)
#--- add OpenMP compile options if needed -------------------------------------
if (OpenMP_FOUND)
target_compile_options(PMusr PUBLIC ${OpenMP_CXX_FLAGS})
@ -157,11 +185,15 @@ if (OpenMP_FOUND)
endif (OpenMP_FOUND)
target_link_libraries(PUserFcnBase ${ROOT_LIBRARIES})
target_link_libraries(PRgeHandler ${Boost_LIBRARIES} ${ROOT_LIBRARIES})
target_link_libraries(PMusr ${DependOnLibs})
#--- install PUserFcnBase solib -----------------------------------------------
install(TARGETS PUserFcnBase DESTINATION lib)
#--- install PRgeHandler solib -----------------------------------------------
install(TARGETS PRgeHandler DESTINATION lib)
#--- install PMusr solib ------------------------------------------------------
install(TARGETS PMusr DESTINATION lib)
@ -173,6 +205,8 @@ install(
${CMAKE_CURRENT_BINARY_DIR}/libPMusrCanvas.rootmap
${CMAKE_CURRENT_BINARY_DIR}/libPMusrT0_rdict.pcm
${CMAKE_CURRENT_BINARY_DIR}/libPMusrT0.rootmap
${CMAKE_CURRENT_BINARY_DIR}/libPRgeHandler_rdict.pcm
${CMAKE_CURRENT_BINARY_DIR}/libPRgeHandler.rootmap
${CMAKE_CURRENT_BINARY_DIR}/libPStartupHandler_rdict.pcm
${CMAKE_CURRENT_BINARY_DIR}/libPStartupHandler.rootmap
${CMAKE_CURRENT_BINARY_DIR}/libPUserFcnBase_rdict.pcm
@ -195,6 +229,7 @@ install(
${MUSRFIT_INC}/PMusr.h
${MUSRFIT_INC}/PMusrT0.h
${MUSRFIT_INC}/PPrepFourier.h
${MUSRFIT_INC}/PRgeHandler.h
${MUSRFIT_INC}/PRunAsymmetry.h
${MUSRFIT_INC}/PRunAsymmetryBNMR.h
${MUSRFIT_INC}/PRunAsymmetryRRF.h
@ -214,6 +249,7 @@ install(
#--- install pkg-config info --------------------------------------------------
install(
FILES ${CMAKE_CURRENT_BINARY_DIR}/PUserFcnBase.pc
${CMAKE_CURRENT_BINARY_DIR}/PRgeHandler.pc
${CMAKE_CURRENT_BINARY_DIR}/PMusr.pc
DESTINATION lib/pkgconfig
)

603
src/classes/PRgeHandler.cpp Normal file
View File

@ -0,0 +1,603 @@
/***************************************************************************
PRgeHandler.cpp
Author: Andreas Suter
e-mail: andreas.suter@psi.ch
***************************************************************************/
/***************************************************************************
* Copyright (C) 2007-2021 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 <fstream>
#include <iostream>
#include <vector>
#include <string>
#include <boost/filesystem.hpp>
#include <boost/algorithm/string.hpp>
#include "PRgeHandler.h"
ClassImpQ(PXmlRgeHandler)
//--------------------------------------------------------------------------
// OnStartDocument
//--------------------------------------------------------------------------
/**
* <p>Called on start of the XML file reading. Initializes all necessary variables.
*/
void PXmlRgeHandler::OnStartDocument()
{
fKey = eEmpty;
}
//--------------------------------------------------------------------------
// OnEndDocument
//--------------------------------------------------------------------------
/**
* <p>Called on end of XML file reading.
*/
void PXmlRgeHandler::OnEndDocument()
{
if (!fIsValid)
return;
// check if anything was set
if (fTrimSpDataPath.empty()) {
std::string err = "<data_path> content is missing!";
fIsValid = false;
OnError(err.c_str());
}
if (fTrimSpFlnPre.empty()) {
std::string err = "<rge_fln_pre> content is missing!";
fIsValid = false;
OnError(err.c_str());
}
if (fTrimSpDataEnergyList.size()==0) {
std::string err = "no implantation energies present!";
fIsValid = false;
OnError(err.c_str());
}
}
//--------------------------------------------------------------------------
// OnStartElement
//--------------------------------------------------------------------------
/**
* <p>Called when a XML start element is found. Filters out the needed elements
* and sets a proper key.
*
* \param str XML element name
* \param attributes not used
*/
void PXmlRgeHandler::OnStartElement(const Char_t *str, const TList *attributes)
{
if (!strcmp(str, "trim_sp")) {
isTrimSp=true;
} else if (!strcmp(str, "data_path") && isTrimSp) {
fKey = eDataPath;
} else if (!strcmp(str, "rge_fln_pre") && isTrimSp) {
fKey = eFlnPre;
} else if (!strcmp(str, "energy") && isTrimSp) {
fKey = eEnergy;
}
}
//--------------------------------------------------------------------------
// OnEndElement
//--------------------------------------------------------------------------
/**
* <p>Called when a XML end element is found. Resets the handler key.
*
* \param str not used
*/
void PXmlRgeHandler::OnEndElement(const Char_t *str)
{
if (!strcmp(str, "trim_sp")) {
isTrimSp=false;
}
fKey = eEmpty;
}
//--------------------------------------------------------------------------
// OnCharacters
//--------------------------------------------------------------------------
/**
* <p>Content of a given XML element. Filters out the data and feeds them to
* the internal variables.
*
* \param str XML element string
*/
void PXmlRgeHandler::OnCharacters(const Char_t *str)
{
int ival;
std::string msg(""), sstr(str);
size_t pos;
switch(fKey) {
case eDataPath:
fTrimSpDataPath=str;
break;
case eFlnPre:
fTrimSpFlnPre=str;
break;
case eEnergy:
try {
ival = std::stoi(str, &pos);
} catch(std::invalid_argument& e) {
fIsValid = false;
msg = "The found energy '" + sstr + "' which is not a number";
OnError(msg.c_str());
} catch(std::out_of_range& e) {
fIsValid = false;
msg = "The found energy '" + sstr + "' which is out-of-range.";
OnError(msg.c_str());
} catch(...) {
fIsValid = false;
msg = "The found energy '" + sstr + "' which generates an error.";
OnError(msg.c_str());
}
if (pos != sstr.length()) {
fIsValid = false;
msg = "The found energy '" + sstr + "' which is not an integer";
OnError(msg.c_str());
}
fTrimSpDataEnergyList.push_back(ival);
break;
default:
break;
}
}
//--------------------------------------------------------------------------
// OnComment
//--------------------------------------------------------------------------
/**
* <p>Called when a XML comment is found. Not used.
*
* \param str not used.
*/
void PXmlRgeHandler::OnComment(const Char_t *str)
{
// nothing to be done for now
}
//--------------------------------------------------------------------------
// OnWarning
//--------------------------------------------------------------------------
/**
* <p>Called when the XML parser emits a warning.
*
* \param str warning string
*/
void PXmlRgeHandler::OnWarning(const Char_t *str)
{
std::cerr << std::endl << "PXmlRgeHandler **WARNING** " << str;
std::cerr << std::endl;
}
//--------------------------------------------------------------------------
// OnError
//--------------------------------------------------------------------------
/**
* <p>Called when the XML parser emits an error.
*
* \param str error string
*/
void PXmlRgeHandler::OnError(const Char_t *str)
{
std::cerr << std::endl << "PXmlRgeHandler **ERROR** " << str;
std::cerr << std::endl;
}
//--------------------------------------------------------------------------
// OnFatalError
//--------------------------------------------------------------------------
/**
* <p>Called when the XML parser emits a fatal error.
*
* \param str fatal error string
*/
void PXmlRgeHandler::OnFatalError(const Char_t *str)
{
std::cerr << std::endl << "PXmlRgeHandler **FATAL ERROR** " << str;
std::cerr << std::endl;
}
//--------------------------------------------------------------------------
// OnCdataBlock
//--------------------------------------------------------------------------
/**
* <p>Not used.
*
* \param str not used
* \param len not used
*/
void PXmlRgeHandler::OnCdataBlock(const Char_t *str, Int_t len)
{
// nothing to be done for now
}
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
ClassImp(PRgeHandler)
//--------------------------------------------------------------------------
// Ctor
//--------------------------------------------------------------------------
/**
* @brief PRgeHandler::PRgeHandler
* @param fln
*/
PRgeHandler::PRgeHandler(const std::string fln)
{
// make sure there is an xml-startup file name
if (fln.empty()) { // fln not given
std::cerr << std::endl;
std::cerr << "**ERROR** NO xml file name provided." << std::endl;
std::cerr << std::endl;
fValid=false;
return;
}
// read the startup xml-file to extract the necessary rge infos.
if (!boost::filesystem::exists(fln)) {
std::cerr << std::endl;
std::cerr << "**ERROR** xml file named: " << fln << " does not exist." << std::endl;
std::cerr << std::endl;
fValid=false;
return;
}
// create the rge xml handler
PXmlRgeHandler *xmlRge = new PXmlRgeHandler();
if (xmlRge == nullptr) {
std::cerr << std::endl;
std::cerr << "**ERROR** couldn't invoke PXmlRgeHandler." << std::endl;
std::cerr << std::endl;
fValid=false;
return;
}
// create the SAX parser
TSAXParser *saxParser = new TSAXParser();
if (saxParser == nullptr) {
std::cerr << std::endl;
std::cerr << "**ERROR** couldn't invoke TSAXParser." << std::endl;
std::cerr << std::endl;
fValid=false;
return;
}
saxParser->SetStopOnError();
// connect SAX parser and rge handler
saxParser->ConnectToHandler("PXmlRgeHandler", xmlRge);
int status = saxParser->ParseFile(fln.c_str());
if (status != 0) {
std::cerr << std::endl;
std::cerr << "**ERROR** parsing the xml-file: " << fln << "." << std::endl;
std::cerr << std::endl;
fValid=false;
return;
}
if (xmlRge->IsValid())
fValid = true;
else
fValid = false;
// read the rge-file(s) content if everything went fine so far
std::string rgeFln;
PRgeData dataSet;
if (fValid) {
const PIntVector energy = xmlRge->GetTrimSpDataVectorList();
for (int i=0; i<energy.size(); i++) {
// construct the file name
rgeFln = xmlRge->GetTrimSpDataPath();
if (rgeFln[rgeFln.size()-1] != '/')
rgeFln += "/";
rgeFln += xmlRge->GetTrimSpFlnPre();
rgeFln += std::to_string(energy[i]);
rgeFln += ".rge";
if (!boost::filesystem::exists(rgeFln)) {
fValid = false;
std::cerr << std::endl;
std::cerr << "**ERROR** rge-file: " << rgeFln << " not found." << std::endl;
std::cerr << std::endl;
break;
}
// init data set
dataSet.energy = energy[i];
dataSet.depth.clear();
dataSet.amplitude.clear();
if (!ReadRgeFile(rgeFln, dataSet)) {
fValid = false;
std::cerr << std::endl;
std::cerr << "**ERROR** read error in rge-file: " << rgeFln << " not found." << std::endl;
std::cerr << std::endl;
break;
}
// get the total number of particles
double tot=0.0;
for (int j=0; j<dataSet.amplitude.size(); j++)
tot += dataSet.amplitude[j];
dataSet.noOfParticles = tot;
// sum_j nn dzz == 1, (dzz_j = depth[j]-depth[j-1])
dataSet.nn.resize(dataSet.amplitude.size());
tot = 0.0;
double zz=0.0; // "previous" zz (since depth[j]-depth[j-1] != const) it needs to be done this way.
for (int j=0; j<dataSet.nn.size(); j++) {
tot += dataSet.amplitude[j] * (dataSet.depth[j] - zz);
zz = dataSet.depth[j];
}
for (int j=0; j<dataSet.nn.size(); j++) {
dataSet.nn[j] = dataSet.amplitude[j] / tot;
}
fData.push_back(dataSet);
}
}
delete saxParser;
delete xmlRge;
}
//--------------------------------------------------------------------------
// ReadRgeFile
//--------------------------------------------------------------------------
/**
* <p>Read the content of a rge-file.
*
* @param fln file name of the rge-file
* @return true on success.
*/
bool PRgeHandler::ReadRgeFile(const std::string fln, PRgeData &data)
{
std::ifstream fin(fln);
if (!fin.is_open())
return false;
std::string line, msg;
std::vector<std::string> tok;
Double_t zz, nn;
size_t pos;
int lineNo=0;
while (fin.good() && fValid) {
std::getline(fin, line);
lineNo++;
boost::algorithm::trim(line);
if (line.empty())
continue;
if (!std::isdigit(line[0]))
continue;
tok.clear();
boost::algorithm::split(tok, line, boost::algorithm::is_any_of(" \t"), boost::algorithm::token_compress_on);
if (tok.size() != 2) {
std::cerr << std::endl;
std::cerr << "**ERROR** in rge-file: " << fln << ", unexpected number of tokens (" << tok.size() << ")." << std::endl;
std::cerr << std::endl;
fin.close();
return false;
}
// check distance
try {
zz = std::stod(tok[0], &pos);
} catch(std::invalid_argument& e) {
fValid = false;
msg = "The found depth '" + tok[0] + "' which is not a number";
std::cerr << std::endl;
std::cerr << "**ERROR** in rge-file: " << fln << ": lineNo: " << lineNo << std::endl;
std::cerr << " " << msg << std::endl;
std::cerr << std::endl;
} catch(std::out_of_range& e) {
fValid = false;
msg = "The found depth '" + tok[0] + "' which is out-of-range.";
std::cerr << std::endl;
std::cerr << "**ERROR** in rge-file: " << fln << ": lineNo: " << lineNo << std::endl;
std::cerr << " " << msg << std::endl;
std::cerr << std::endl;
} catch(...) {
fValid = false;
msg = "The found depth '" + tok[0] + "' which generates an error.";
std::cerr << std::endl;
std::cerr << "**ERROR** in rge-file: " << fln << ": lineNo: " << lineNo << std::endl;
std::cerr << " " << msg << std::endl;
std::cerr << std::endl;
}
if (pos != tok[0].length()) {
fValid = false;
msg = "The found depth '" + tok[0] + "' which is not an number.";
std::cerr << std::endl;
std::cerr << "**ERROR** in rge-file: " << fln << ": lineNo: " << lineNo << std::endl;
std::cerr << " " << msg << std::endl;
std::cerr << std::endl;
}
// check number of implanted particles
try {
nn = std::stod(tok[1], &pos);
} catch(std::invalid_argument& e) {
fValid = false;
msg = "The found #particles '" + tok[1] + "' which is not a number";
std::cerr << std::endl;
std::cerr << "**ERROR** in rge-file: " << fln << ": lineNo: " << lineNo << std::endl;
std::cerr << " " << msg << std::endl;
std::cerr << std::endl;
} catch(std::out_of_range& e) {
fValid = false;
msg = "The found #particles '" + tok[1] + "' which is out-of-range.";
std::cerr << std::endl;
std::cerr << "**ERROR** in rge-file: " << fln << ": lineNo: " << lineNo << std::endl;
std::cerr << " " << msg << std::endl;
std::cerr << std::endl;
} catch(...) {
fValid = false;
msg = "The found #particles '" + tok[1] + "' which generates an error.";
std::cerr << std::endl;
std::cerr << "**ERROR** in rge-file: " << fln << ": lineNo: " << lineNo << std::endl;
std::cerr << " " << msg << std::endl;
std::cerr << std::endl;
}
if (pos != tok[1].length()) {
fValid = false;
msg = "The found #particles '" + tok[1] + "' which is not an integer.";
std::cerr << std::endl;
std::cerr << "**ERROR** in rge-file: " << fln << ": lineNo: " << lineNo << std::endl;
std::cerr << " " << msg << std::endl;
std::cerr << std::endl;
}
data.depth.push_back(zz/10.0); // distance in nm
data.amplitude.push_back(nn);
}
fin.close();
return true;
}
//--------------------------------------------------------------------------
// GetZmax via energy
//--------------------------------------------------------------------------
/**
* <p>Get maximal depth for a given energy.
*
* @param energy energy in (eV)
* @return zMax if energy is found, -1 otherwise.
*/
Double_t PRgeHandler::GetZmax(const Double_t energy)
{
int idx=-1;
for (int i=0; i<fData.size(); i++) {
if (fabs(fData[i].energy-energy) < 1.0) {
idx = i;
break;
}
}
if (idx != -1)
return GetZmax(idx);
return -1.0;
}
//--------------------------------------------------------------------------
// GetZmax via index
//--------------------------------------------------------------------------
/**
* <p>Get maximal depth for a given index.
*
* @param idx index for which zMax is requested.
* @return zMax if idx is in range, -1 otherwise.
*/
Double_t PRgeHandler::GetZmax(const Int_t idx)
{
if ((idx < 0) || (idx >= fData.size()))
return -1.0;
return fData[idx].depth[fData[idx].depth.size()-1];
}
//--------------------------------------------------------------------------
// Get_n via energy
//--------------------------------------------------------------------------
/**
* <p>Get the normalized n(E,z) value.
*
* @param energy (eV)
* @param z (nm)
* @return n(E,z) if energy and z are in proper range, -1.0 otherwise.
*/
Double_t PRgeHandler::Get_n(const Double_t energy, const Double_t z)
{
int idx=-1;
for (int i=0; i<fData.size(); i++) {
if (fabs(fData[i].energy-energy) < 1.0) {
idx = i;
break;
}
}
if (idx == -1)
return 0.0;
return Get_n(idx, z);
}
//--------------------------------------------------------------------------
// Get_n via index
//--------------------------------------------------------------------------
/**
* <p>Get the normalized n(idx,z) value.
*
* @param idx index of the rge-dataset
* @param z (nm)
* @return n(idx,z) if idx and z are in proper range, -1.0 otherwise.
*/
Double_t PRgeHandler::Get_n(const Int_t idx, const Double_t z)
{
if ((idx < 0) || (idx >= fData.size()))
return 0.0;
if ((z < 0.0) || (z > GetZmax(idx)))
return 0.0;
int pos=0;
for (int i=0; i<fData[idx].depth.size(); i++) {
if (z <= fData[idx].depth[i]) {
pos = i-1;
break;
}
}
Double_t nn=0.0;
if (pos < 0) {
nn = fData[idx].nn[0] * z/(fData[idx].depth[1]-fData[idx].depth[0]);
} else { // linear interpolation
nn = fData[idx].nn[pos] +
(fData[idx].nn[pos+1] - fData[idx].nn[pos]) *
(z-fData[idx].depth[pos])/(fData[idx].depth[pos+1]-fData[idx].depth[pos]);
}
return nn;
}
//--------------------------------------------------------------------------
// GetEnergyIndex
//--------------------------------------------------------------------------
/**
* <p>Get the energy index by providing an energy in (eV).
*
* @param energy in (eV).
* @return energy index if energy was found, -1 otherwise.
*/
Int_t PRgeHandler::GetEnergyIndex(const Double_t energy)
{
int idx=-1;
for (int i=0; i<fData.size(); i++) {
if (fabs(fData[i].energy-energy) < 1.0) {
idx = i;
break;
}
}
return idx;
}

10
src/classes/PRgeHandler.pc.in Normal file
View File

@ -0,0 +1,10 @@
prefix=@prefix@
exec_prefix=@exec_prefix@
libdir=@libdir@
includedir=@includedir@
Name: PRgeHandler
Description: C++ shared library providing the support for handling of TrimSP rge-files
Version: @MUSR_VERSION@
Libs: -L${libdir} -l@PRGEHANDLER_LIBRARY_NAME@
Cflags: -I${includedir}

2
src/classes/PRunAsymmetry.cpp
View File

@ -1297,7 +1297,7 @@ Bool_t PRunAsymmetry::PrepareViewData(PRawRunData* runData, UInt_t histoNo[2])
// calculate theory
Double_t time;
UInt_t size = runData->GetDataBin(histoNo[0])->size();
UInt_t size = runData->GetDataBin(histoNo[0])->size()/packing;
Int_t factor = 8; // 8 times more points for the theory (if fTheoAsData == false)
fData.SetTheoryTimeStart(fData.GetDataTimeStart());

2
src/classes/PRunSingleHisto.cpp
View File

@ -1315,7 +1315,7 @@ Bool_t PRunSingleHisto::PrepareViewData(PRawRunData* runData, const UInt_t histo
// calculate theory
Double_t theoryValue;
UInt_t size = fForward.size();
UInt_t size = fForward.size()/packing;
Int_t factor = 8; // 8 times more points for the theory (if fTheoAsData == false)
UInt_t rebinRRF = 0;

4
src/external/MagProximity/CMakeLists.txt vendored
View File

@ -38,7 +38,6 @@ configure_file("PMagProximityFitter.pc.in" "PMagProximityFitter.pc" @ONLY)
add_library(PMagProximityFitter SHARED
PMagProximityFitter.cpp
PMagProximityFitterDict.cxx
PMPRgeHandler.cpp
PMPStartupHandler.cpp
PMPStartupHandlerDict.cxx
)
@ -57,7 +56,7 @@ set_target_properties(PMagProximityFitter
)
#--- add library dependencies -------------------------------------------------
target_link_libraries(PMagProximityFitter ${FFTW3_LIBRARY} ${ROOT_LIBRARIES} PUserFcnBase)
target_link_libraries(PMagProximityFitter ${FFTW3_LIBRARY} ${ROOT_LIBRARIES} PRgeHandler PUserFcnBase)
#--- install PMagProximityFitter solib ----------------------------------------
install(TARGETS PMagProximityFitter DESTINATION lib)
@ -76,7 +75,6 @@ install(
FILES
PMagProximity.h
PMagProximityFitter.h
PMPRgeHandler.h
PMPStartupHandler.h
DESTINATION
include

212
src/external/MagProximity/PMPRgeHandler.cpp vendored
View File

@ -1,212 +0,0 @@
/***************************************************************************
PMPRgeHandler.cpp
Author: Andreas Suter
e-mail: andreas.suter@psi.ch
***************************************************************************/
/***************************************************************************
* Copyright (C) 2011-2021 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>
#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 = static_cast<UInt_t>(dist/(fRgeDataList[index].stoppingDistance[1]-fRgeDataList[index].stoppingDistance[0]));
if (distIdx >= fRgeDataList[index].stoppingDistance.size()) {
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)
{
std::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(), std::iostream::in);
if (!fin.is_open()) {
std::cout << std::endl << "PMPRgeHandler::LoadRgeData **ERROR**";
std::cout << std::endl << " Could not open file " << rgeDataPathList[i].Data();
std::cout << std::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;
}

53
src/external/MagProximity/PMPRgeHandler.h vendored
View File

@ -1,53 +0,0 @@
/***************************************************************************
PMPRgeHandler.h
Author: Andreas Suter
e-mail: andreas.suter@psi.ch
***************************************************************************/
/***************************************************************************
* Copyright (C) 2011-2021 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_

2
src/external/MagProximity/PMagProximityFitter.cpp vendored
View File

@ -89,7 +89,7 @@ PMagProximityFitterGlobal::PMagProximityFitterGlobal()
}
// load all the TRIM.SP rge-files
fRgeHandler = new PMPRgeHandler(fStartupHandler->GetTrimSpDataPathList(), fStartupHandler->GetTrimSpDataVectorList());
fRgeHandler = new PRgeHandler();
if (!fRgeHandler->IsValid()) {
std::cout << std::endl << ">> PMagProximityFitterGlobal::PMagProximityFitterGlobal **PANIC ERROR**";
std::cout << std::endl << ">> rge data handler too unhappy. Will terminate unfriendly, sorry.";

8
src/external/MagProximity/PMagProximityFitter.h vendored
View File

@ -32,7 +32,7 @@
#include "PUserFcnBase.h"
#include "PMPStartupHandler.h"
#include "PMPRgeHandler.h"
#include "PRgeHandler.h"
class PMagProximityFitterGlobal
{
@ -42,15 +42,15 @@ class 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 Int_t GetEnergyIndex(const Double_t energy) { return fRgeHandler->GetEnergyIndex(energy); }
virtual Double_t GetMuonStoppingDensity(const Int_t energyIndex, const Double_t z) const { return fRgeHandler->Get_n(energyIndex, z); }
virtual Double_t GetMagneticField(const Double_t z) const;
private:
Bool_t fValid;
PMPStartupHandler *fStartupHandler;
PMPRgeHandler *fRgeHandler;
PRgeHandler *fRgeHandler;
mutable std::vector<Double_t> fPreviousParam;

7
src/external/MagProximity/mag_proximity_startup.xml vendored
View File

@ -3,13 +3,14 @@
<comment>
contains meta file infromation of the needed trim.sp files.
</comment>
<trim_sp_part>
<data_path>/mnt/home/nemu/analysis/2010/EuS-Co/trimsp2/EuS-Co_E</data_path>
<trim_sp>
<data_path>/afs/psi.ch/project/nemu/analysis/2010/EuS-Co/trimsp2</data_path>
<rge_fln_pre>EuS-Co_E</rge_fln_pre>
<energy_list>
<energy>5030</energy>
<energy>6330</energy>
<energy>7530</energy>
<energy>8730</energy>
</energy_list>
</trim_sp_part>
</trim_sp>
</mag_proximity>

7
src/external/Nonlocal/CMakeLists.txt vendored
View File

@ -36,7 +36,7 @@ set(prefix "${CMAKE_INSTALL_PREFIX}")
set(exec_prefix "\$\{prefix\}")
set(libdir "\$\{exec_prefix\}/lib")
set(includedir "\$\{prefix\}/include")
set(PNL_PIPPARDFITTER_VERSION "1.0.0")
set(PNL_PIPPARDFITTER_VERSION "1.1.0")
set(PNL_PIPPARDFITTER_LIBRARY_NAME "PNL_PippardFitter")
configure_file("PNL_PippardFitter.pc.in" "PNL_PippardFitter.pc" @ONLY)
@ -44,7 +44,6 @@ configure_file("PNL_PippardFitter.pc.in" "PNL_PippardFitter.pc" @ONLY)
add_library(PNL_PippardFitter SHARED
PNL_PippardFitter.cpp
PNL_PippardFitterDict.cxx
PNL_RgeHandler.cpp
PNL_StartupHandler.cpp
PNL_StartupHandlerDict.cxx
)
@ -63,7 +62,7 @@ set_target_properties(PNL_PippardFitter
)
#--- add library dependencies -------------------------------------------------
target_link_libraries(PNL_PippardFitter ${FFTW3_LIBRARY} ${ROOT_LIBRARIES} PUserFcnBase)
target_link_libraries(PNL_PippardFitter ${Boost_LIBRARIES} ${FFTW3_LIBRARY} ${ROOT_LIBRARIES} PRgeHandler PUserFcnBase)
#--- install PNL_PippardFitter solib ------------------------------------------
install(TARGETS PNL_PippardFitter DESTINATION lib)
@ -80,9 +79,7 @@ install(
#--- install PNL_PippardFitter header -----------------------------------------
install(
FILES
PNonlocal.h
PNL_PippardFitter.h
PNL_RgeHandler.h
PNL_StartupHandler.h
DESTINATION
include

54
src/external/Nonlocal/PNL_PippardFitter.cpp vendored
View File

@ -31,7 +31,6 @@
#include <cmath>
#include <iostream>
using namespace std;
#include <TSAXParser.h>
#include <TMath.h>
@ -51,10 +50,6 @@ ClassImp(PNL_PippardFitterGlobal)
*/
PNL_PippardFitterGlobal::PNL_PippardFitterGlobal()
{
fValid = true;
fStartupHandler = 0;
fRgeHandler = 0;
// read XML startup file
char startup_path_name[128];
TSAXParser *saxParser = new TSAXParser();
@ -67,8 +62,8 @@ PNL_PippardFitterGlobal::PNL_PippardFitterGlobal()
Int_t status = parseXmlFile(saxParser, startup_path_name);
// check for parse errors
if (status) { // error
cout << endl << ">> PNL_PippardFitterGlobal::PNL_PippardFitterGlobal: **WARNING** Reading/parsing nonlocal_startup.xml failed.";
cout << endl;
std::cout << std::endl << ">> PNL_PippardFitterGlobal::PNL_PippardFitterGlobal: **WARNING** Reading/parsing nonlocal_startup.xml failed.";
std::cout << std::endl;
fValid = false;
}
@ -80,20 +75,20 @@ PNL_PippardFitterGlobal::PNL_PippardFitterGlobal()
// check if everything went fine with the startup handler
if (!fStartupHandler->IsValid()) {
cout << endl << ">> PNL_PippardFitterGlobal::PNL_PippardFitterGlobal **PANIC ERROR**";
cout << endl << ">> startup handler too unhappy. Will terminate unfriendly, sorry.";
cout << endl;
std::cout << std::endl << ">> PNL_PippardFitterGlobal::PNL_PippardFitterGlobal **PANIC ERROR**";
std::cout << std::endl << ">> startup handler too unhappy. Will terminate unfriendly, sorry.";
std::cout << std::endl;
fValid = false;
}
fFourierPoints = fStartupHandler->GetFourierPoints();
// load all the TRIM.SP rge-files
fRgeHandler = new PNL_RgeHandler(fStartupHandler->GetTrimSpDataPathList(), fStartupHandler->GetTrimSpDataVectorList());
fRgeHandler = new PRgeHandler("./nonlocal_startup.xml");
if (!fRgeHandler->IsValid()) {
cout << endl << ">> PNL_PippardFitterGlobal::PNL_PippardFitterGlobal **PANIC ERROR**";
cout << endl << ">> rge data handler too unhappy. Will terminate unfriendly, sorry.";
cout << endl;
std::cout << std::endl << ">> PNL_PippardFitterGlobal::PNL_PippardFitterGlobal **PANIC ERROR**";
std::cout << std::endl << ">> rge data handler too unhappy. Will terminate unfriendly, sorry.";
std::cout << std::endl;
fValid = false;
}
@ -182,16 +177,16 @@ void PNL_PippardFitterGlobal::CalculateField(const std::vector<Double_t> &param)
if (fFieldq == 0) {
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;
std::cout << std::endl << "PPippard::CalculateField(): **ERROR** couldn't allocate memory for fFieldq";
std::cout << std::endl;
return;
}
}
if (fFieldB == 0) {
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;
std::cout << std::endl << "PPippard::CalculateField(): **ERROR** couldn't allocate memory for fFieldB";
std::cout << std::endl;
return;
}
}
@ -349,19 +344,6 @@ Double_t PNL_PippardFitterGlobal::XiP_T(const Double_t xi0, const Double_t meanF
ClassImp(PNL_PippardFitter)
//--------------------------------------------------------------------------
// Constructor
//--------------------------------------------------------------------------
/**
*
*/
PNL_PippardFitter::PNL_PippardFitter()
{
fValid = false;
fInvokedGlobal = false;
fPippardFitterGlobal = 0;
}
//--------------------------------------------------------------------------
// Destructor
//--------------------------------------------------------------------------
@ -387,18 +369,18 @@ PNL_PippardFitter::~PNL_PippardFitter()
* \param globalPart
* \param idx
*/
void PNL_PippardFitter::SetGlobalPart(vector<void*> &globalPart, UInt_t idx)
void PNL_PippardFitter::SetGlobalPart(std::vector<void*> &globalPart, UInt_t idx)
{
fIdxGlobal = static_cast<Int_t>(idx);
if ((Int_t)globalPart.size() <= fIdxGlobal) {
fPippardFitterGlobal = new PNL_PippardFitterGlobal();
if (fPippardFitterGlobal == 0) {
if (fPippardFitterGlobal == nullptr) {
fValid = false;
cerr << endl << ">> PNL_PippardFitter::SetGlobalPart(): **ERROR** Couldn't invoke global user function object, sorry ..." << endl;
std::cerr << std::endl << ">> PNL_PippardFitter::SetGlobalPart(): **ERROR** Couldn't invoke global user function object, sorry ..." << std::endl;
} else if (!fPippardFitterGlobal->IsValid()) {
fValid = false;
cerr << endl << ">> PNL_PippardFitter::SetGlobalPart(): **ERROR** initialization of global user function object failed, sorry ..." << endl;
std::cerr << std::endl << ">> PNL_PippardFitter::SetGlobalPart(): **ERROR** initialization of global user function object failed, sorry ..." << std::endl;
} else {
fValid = true;
fInvokedGlobal = true;
@ -447,7 +429,7 @@ Double_t PNL_PippardFitter::operator()(Double_t t, const std::vector<Double_t> &
fPippardFitterGlobal->CalculateField(param);
Int_t energyIndex = fPippardFitterGlobal->GetEnergyIndex(param[0]);
if (energyIndex == -1) { // energy not found
cerr << endl << ">> PNL_PippardFitter::operator() energy " << param[0] << " not found. Will terminate." << endl;
std::cerr << std::endl << ">> PNL_PippardFitter::operator() energy " << param[0] << " not found. Will terminate." << std::endl;
assert(0);
}

23
src/external/Nonlocal/PNL_PippardFitter.h vendored
View File

@ -36,9 +36,10 @@
//#endif
#include <fftw3.h>
#include "PMusr.h"
#include "PUserFcnBase.h"
#include "PNL_StartupHandler.h"
#include "PNL_RgeHandler.h"
#include "PRgeHandler.h"
class PNL_PippardFitterGlobal
{
@ -49,15 +50,15 @@ class PNL_PippardFitterGlobal
Bool_t IsValid() { return fValid; }
virtual void SetTempExponent(const Double_t nn) { f_nn = nn; }
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 Int_t GetEnergyIndex(const Double_t energy) { return fRgeHandler->GetEnergyIndex(energy); }
virtual Double_t GetMuonStoppingDensity(const Int_t energyIndex, const Double_t z) const { return fRgeHandler->Get_n(energyIndex, z); }
virtual Double_t GetMagneticField(const Double_t z) const;
private:
Bool_t fValid;
Bool_t fValid{true};
PNL_StartupHandler *fStartupHandler;
PNL_RgeHandler *fRgeHandler;
PNL_StartupHandler *fStartupHandler{nullptr};
PRgeHandler *fRgeHandler{nullptr};
mutable std::vector<Double_t> fPreviousParam;
@ -85,21 +86,21 @@ class PNL_PippardFitterGlobal
class PNL_PippardFitter : public PUserFcnBase
{
public:
PNL_PippardFitter();
PNL_PippardFitter() {}
virtual ~PNL_PippardFitter();
virtual Bool_t NeedGlobalPart() const { return true; }
virtual void SetGlobalPart(vector<void*> &globalPart, UInt_t idx);
virtual void SetGlobalPart(std::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;
Bool_t fValid{true};
Bool_t fInvokedGlobal{false};
Int_t fIdxGlobal;
PNL_PippardFitterGlobal *fPippardFitterGlobal;
PNL_PippardFitterGlobal *fPippardFitterGlobal{nullptr};
ClassDef(PNL_PippardFitter, 1)
};

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

@ -1,208 +0,0 @@
/***************************************************************************
PNL_RgeHandler.cpp
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. *
***************************************************************************/
#include <cassert>
#include <iostream>
#include <fstream>
#include "PNL_RgeHandler.h"
//--------------------------------------------------------------------------
// Constructor
//--------------------------------------------------------------------------
/**
*
*/
PNL_RgeHandler::PNL_RgeHandler(const PStringVector &rgeDataPathList, const PDoubleVector &rgeDataEnergyList)
{
fIsValid = false;
fIsValid = LoadRgeData(rgeDataPathList, rgeDataEnergyList);
}
//--------------------------------------------------------------------------
// Destructor
//--------------------------------------------------------------------------
/**
*
*/
PNL_RgeHandler::~PNL_RgeHandler()
{
fRgeDataList.clear();
}
//--------------------------------------------------------------------------
// GetRgeEnergyIndex
//--------------------------------------------------------------------------
/**
*
*/
Int_t PNL_RgeHandler::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
//--------------------------------------------------------------------------
/**
*
*/
Double_t PNL_RgeHandler::GetRgeValue(const Int_t index, const Double_t dist)
{
Double_t rgeVal = 0.0;
// find the proper position index
Int_t distIdx = -2;
for (UInt_t i=0; i<fRgeDataList[index].stoppingDistance.size(); i++) {
if (fRgeDataList[index].stoppingDistance[i] > dist) {
distIdx = i-1;
break;
}
}
// get the proper n(z) value
if (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
//--------------------------------------------------------------------------
/**
*
*/
Double_t PNL_RgeHandler::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 PNL_RgeHandler::LoadRgeData(const PStringVector &rgeDataPathList, const PDoubleVector &rgeDataEnergyList)
{
std::ifstream fin;
PNL_RgeData data;
TString 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(), std::iostream::in);
if (!fin.is_open()) {
std::cout << std::endl << "PNL_RgeHandler::LoadRgeData **ERROR**";
std::cout << std::endl << " Could not open file " << rgeDataPathList[i].Data();
std::cout << std::endl;
return false;
}
// keep energy (in keV)
data.energy = rgeDataEnergyList[i]/1000.0;
// read msr-file
while (!fin.eof()) {
// read a line
fin.getline(line, sizeof(line));
// ignore empty lines
tstr = line;
if (tstr.IsWhitespace())
continue;
// get values
result = sscanf(line, "%lf %lf", &dist, &val);
// check if data are valid, otherwise ignore it
if (result != 2) {
continue;
}
// feed fRgeDataList
data.stoppingDistance.push_back(dist/10.0);
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/Nonlocal/PNL_RgeHandler.h vendored
View File

@ -1,55 +0,0 @@
/***************************************************************************
PNL_RgeHandler.h
Author: Andreas Suter
e-mail: andreas.suter@psi.ch
$Id$
***************************************************************************/
/***************************************************************************
* Copyright (C) 2009-2021 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 _PNL_RGEHANDLER_H_
#define _PNL_RGEHANDLER_H_
#include "PNonlocal.h"
class PNL_RgeHandler
{
public:
PNL_RgeHandler(const PStringVector &rgeDataPathList, const PDoubleVector &rgeDataEnergyList);
virtual ~PNL_RgeHandler();
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;
PNL_RgeDataList fRgeDataList;
virtual Bool_t LoadRgeData(const PStringVector &rgeDataPathList, const PDoubleVector &rgeDataEnergyList);
};
#endif // _PNL_RGEHANDLER_H_

70
src/external/Nonlocal/PNL_StartupHandler.cpp vendored
View File

@ -32,6 +32,8 @@
#include <iostream>
#include <fstream>
#include <boost/filesystem.hpp>
#include "PNL_StartupHandler.h"
ClassImpQ(PNL_StartupHandler)
@ -44,46 +46,26 @@ ClassImpQ(PNL_StartupHandler)
*/
PNL_StartupHandler::PNL_StartupHandler()
{
fIsValid = true;
fStartupFileFound = false;
fStartupFilePath = "";
fFourierPoints = 0;
fTrimSpDataPath = TString("");
// get default path (for the moment only linux like)
char startup_path_name[512];
char *home_str=0;
// check if the startup file is found in the current directory
strcpy(startup_path_name, "./nonlocal_startup.xml");
if (StartupFileExists(startup_path_name)) {
if (boost::filesystem::exists(startup_path_name)) {
fStartupFileFound = true;
fStartupFilePath = TString(startup_path_name);
} else { // startup file is not found in the current directory
std::cout << std::endl << "PNL_StartupHandler(): **WARNING** Couldn't find nonlocal_startup.xml in the current directory, will try default one." << std::endl;
home_str = getenv("HOME");
snprintf(startup_path_name, sizeof(startup_path_name), "%s/.musrfit/nonlocal_startup.xml", home_str);
if (StartupFileExists(startup_path_name)) {
if (boost::filesystem::exists(startup_path_name)) {
fStartupFileFound = true;
fStartupFilePath = TString(startup_path_name);
}
}
}
//--------------------------------------------------------------------------
// Destructor
//--------------------------------------------------------------------------
/**
*
*/
PNL_StartupHandler::~PNL_StartupHandler()
{
fTrimSpDataPathList.clear();
fTrimSpDataEnergyList.clear();
}
//--------------------------------------------------------------------------
// OnStartDocument
//--------------------------------------------------------------------------
@ -124,10 +106,6 @@ void PNL_StartupHandler::OnStartElement(const char *str, const TList *attributes
{
if (!strcmp(str, "fourier_points")) {
fKey = eFourierPoints;
} else if (!strcmp(str, "data_path")) {
fKey = eDataPath;
} else if (!strcmp(str, "energy")) {
fKey = eEnergy;
}
}
@ -167,23 +145,6 @@ void PNL_StartupHandler::OnCharacters(const char *str)
std::cout << std::endl;
}
break;
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 {
std::cout << std::endl << "PNL_StartupHandler::OnCharacters: **ERROR** when finding energy:";
std::cout << std::endl << "\"" << str << "\" is not a floating point number, will ignore it and use the default value.";
std::cout << std::endl;
}
break;
default:
break;
}
@ -257,29 +218,6 @@ void PNL_StartupHandler::OnCdataBlock(const char *str, Int_t len)
// nothing to be done for now
}
//--------------------------------------------------------------------------
// StartupFileExists
//--------------------------------------------------------------------------
/**
* <p>
*
*/
bool PNL_StartupHandler::StartupFileExists(char *fln)
{
bool result = false;
std::ifstream ifile(fln);
if (ifile.fail()) {
result = false;
} else {
result = true;
ifile.close();
}
return result;
}
// -------------------------------------------------------------------------
// end
// -------------------------------------------------------------------------

24
src/external/Nonlocal/PNL_StartupHandler.h vendored
View File

@ -5,8 +5,6 @@
Author: Andreas Suter
e-mail: andreas.suter@psi.ch
$Id$
***************************************************************************/
/***************************************************************************
@ -36,13 +34,11 @@
#include <TQObject.h>
#include <TString.h>
#include "PNonlocal.h"
class PNL_StartupHandler : public TObject
{
public:
PNL_StartupHandler();
virtual ~PNL_StartupHandler();
virtual ~PNL_StartupHandler() {}
virtual void OnStartDocument(); // SLOT
virtual void OnEndDocument(); // SLOT
@ -58,26 +54,18 @@ class PNL_StartupHandler : public TObject
virtual bool IsValid() { return fIsValid; }
virtual TString GetStartupFilePath() { return fStartupFilePath; }
virtual const Int_t GetFourierPoints() const { return fFourierPoints; }
virtual const PStringVector GetTrimSpDataPathList() const { return fTrimSpDataPathList; }
virtual const PDoubleVector GetTrimSpDataVectorList() const { return fTrimSpDataEnergyList; }
virtual bool StartupFileFound() { return fStartupFileFound; }
private:
enum EKeyWords {eEmpty, eComment, eFourierPoints, eDataPath, eEnergy};
enum EKeyWords {eEmpty, eFourierPoints};
EKeyWords fKey;
bool fIsValid;
bool fIsValid{true};
bool fStartupFileFound;
TString fStartupFilePath;
bool fStartupFileFound{false};
TString fStartupFilePath{""};
Int_t fFourierPoints;
TString fTrimSpDataPath;
PStringVector fTrimSpDataPathList;
PDoubleVector fTrimSpDataEnergyList;
bool StartupFileExists(char *fln);
Int_t fFourierPoints{0};
ClassDef(PNL_StartupHandler, 1)
};

65
src/external/Nonlocal/PNonlocal.h vendored
View File

@ -1,65 +0,0 @@
/***************************************************************************
PNonlocal.h
Author: Andreas Suter
e-mail: andreas.suter@psi.ch
***************************************************************************/
/***************************************************************************
* Copyright (C) 2009-2021 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 _PNONLOCAL_H_
#define _PNONLOCAL_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;
} PNL_RgeData;
//-------------------------------------------------------------
/**
* <p>
*/
typedef std::vector<PNL_RgeData> PNL_RgeDataList;
#endif // _PNONLOCAL_H_

32
src/external/Nonlocal/nonlocal_startup.xml vendored
View File

@ -1,25 +1,27 @@
<?xml version="1.0" encoding="UTF-8"?>
<nonlocal xmlns="http://nemu.web.psi.ch/musrfit/nonlocal">
<comment>
trim.sp meta information
Fourier and TrimSp information
</comment>
<nonlocal_par>
<fourier_points>262144</fourier_points>
</nonlocal_par>
<trim_sp_part>
<data_path>/afs/psi.ch/project/nemu/analysis/2009/Nonlocal/trimsp/InSne</data_path>
<trim_sp>
<data_path>./profiles/</data_path>
<rge_fln_pre>Sn_E</rge_fln_pre>
<energy_list>
<energy>2.5</energy>
<energy>4.0</energy>
<energy>6.0</energy>
<energy>8.0</energy>
<energy>10.0</energy>
<energy>12.5</energy>
<energy>14.1</energy>
<energy>17.5</energy>
<energy>22.0</energy>
<energy>25.0</energy>
<energy>28.2</energy>
<energy>1000</energy>
<energy>2000</energy>
<energy>4000</energy>
<energy>6000</energy>
<energy>8000</energy>
<energy>10000</energy>
<energy>12000</energy>
<energy>14100</energy>
<energy>18000</energy>
<energy>22000</energy>
<energy>25000</energy>
<energy>27300</energy>
</energy_list>
</trim_sp_part>
</trim_sp>
</nonlocal>

17
src/external/libPhotoMeissner/classes/CMakeLists.txt vendored
View File

@ -21,16 +21,6 @@ root_generate_dictionary(
LINKDEF ${PHOTO_MEISSNER_INC}/PPhotoMeissnerLinkDef.h
MODULE PPhotoMeissner
)
root_generate_dictionary(
PStartupHandler_PMDict
PStartupHandler_PM.h
OPTIONS
-I${MUSRFIT_INC}
-I${PHOTO_MEISSNER_INC}
-inlineInputHeader
LINKDEF ${PHOTO_MEISSNER_INC}/PStartupHandler_PMLinkDef.h
MODULE PStartupHandler_PM
)
#--- create pkg-config info ---------------------------------------------------
set(prefix "${CMAKE_INSTALL_PREFIX}")
@ -45,8 +35,6 @@ configure_file("PPhotoMeissner.pc.in" "PPhotoMeissner.pc" @ONLY)
add_library(PPhotoMeissner SHARED
PPhotoMeissner.cpp
PPhotoMeissnerDict.cxx
PStartupHandler_PM.cpp
PStartupHandler_PMDict.cxx
)
#--- set target properties, e.g. version --------------------------------------
@ -67,7 +55,7 @@ target_include_directories(
#--- add library dependencies -------------------------------------------------
target_link_libraries(PPhotoMeissner
${FFTW3_LIBRARY} ${GSL_LIBRARY} ${ROOT_LIBRARIES} PUserFcnBase
${FFTW3_LIBRARY} ${GSL_LIBRARY} ${ROOT_LIBRARIES} PRgeHandler PUserFcnBase
)
#--- install PPhotoMeissner solib ---------------------------------------------
@ -77,8 +65,6 @@ install(TARGETS PPhotoMeissner DESTINATION lib)
install(
FILES ${CMAKE_CURRENT_BINARY_DIR}/libPPhotoMeissner_rdict.pcm
${CMAKE_CURRENT_BINARY_DIR}/libPPhotoMeissner.rootmap
${CMAKE_CURRENT_BINARY_DIR}/libPStartupHandler_PM_rdict.pcm
${CMAKE_CURRENT_BINARY_DIR}/libPStartupHandler_PM.rootmap
DESTINATION lib
)
@ -86,7 +72,6 @@ install(
install(
FILES
${CMAKE_CURRENT_SOURCE_DIR}/../include/PPhotoMeissner.h
${CMAKE_CURRENT_SOURCE_DIR}/../include/PStartupHandler_PM.h
DESTINATION
include
)

68
src/external/libPhotoMeissner/classes/PPhotoMeissner.cpp vendored
View File

@ -8,7 +8,7 @@
***************************************************************************/
/***************************************************************************
* Copyright (C) 2013 by Andreas Suter *
* Copyright (C) 2013-2021 by Andreas Suter *
* andreas.suter@psi.ch *
* *
* This program is free software; you can redistribute it and/or modify *
@ -30,15 +30,11 @@
#include <cmath>
#include <cassert>
#include <iostream>
using namespace std;
#include <gsl/gsl_sf_bessel.h>
#include <TSAXParser.h>
#include "PMusr.h"
#include "../include/PPhotoMeissner.h"
#include "PPhotoMeissner.h"
ClassImp(PPhotoMeissner) // for ROOT dictionary
@ -51,38 +47,14 @@ ClassImp(PPhotoMeissner) // for ROOT dictionary
*/
PPhotoMeissner::PPhotoMeissner()
{
// init
fStartupHandler = 0;
fValid = true;
// read XML startup file
char startup_path_name[128];
TSAXParser *saxParser = new TSAXParser();
fStartupHandler = new PStartupHandler_PM();
strcpy(startup_path_name, fStartupHandler->GetStartupFilePath().Data());
saxParser->ConnectToHandler("PStartupHandler_PM", fStartupHandler);
//Int_t status = saxParser->ParseFile(startup_path_name);
// parsing the file as above seems to lead to problems in certain environments;
// use the parseXmlFile function instead (see PUserFcnBase.cpp for the definition)
Int_t status = parseXmlFile(saxParser, startup_path_name);
// check for parse errors
if (status) { // error
cout << endl << ">> PPhotoMeissner::PPhotoMeissner: **WARNING** Reading/parsing photoMeissner_startup.xml failed.";
cout << endl;
fValid = false;
}
// clean up
if (saxParser) {
delete saxParser;
saxParser = 0;
}
fRgeHandler = new PRgeHandler("./photoMeissner_startup.xml");
// check if everything went fine with the startup handler
if (!fStartupHandler->IsValid()) {
cout << endl << ">> PPhotoMeissner::PPhotoMeissner **PANIC ERROR**";
cout << endl << ">> startup handler too unhappy. Will terminate unfriendly, sorry.";
cout << endl;
if (!fRgeHandler->IsValid()) {
std::cout << std::endl << ">> PPhotoMeissner::PPhotoMeissner **PANIC ERROR**";
std::cout << std::endl << ">> startup handler too unhappy. Will terminate unfriendly, sorry.";
std::cout << std::endl;
fValid = false;
}
}
@ -95,8 +67,8 @@ PPhotoMeissner::PPhotoMeissner()
*/
PPhotoMeissner::~PPhotoMeissner()
{
if (fStartupHandler)
delete fStartupHandler;
if (fRgeHandler)
delete fRgeHandler;
}
//--------------------------------------------------------------------------
@ -108,9 +80,9 @@ PPhotoMeissner::~PPhotoMeissner()
* \param t
* \param par
*/
Double_t PPhotoMeissner::operator()(Double_t t, const vector<Double_t> &par) const
Double_t PPhotoMeissner::operator()(Double_t t, const std::vector<Double_t> &par) const
{
// expected parameters: (0) implantation energy (kV), (1) dead layer (nm), (2) Bext (G), (3) lambdaLondon (nm),
// expected parameters: (0) implantation energy (eV), (1) dead layer (nm), (2) Bext (G), (3) lambdaLondon (nm),
// (4) lambdaPhoto (nm), (5) z0 (nm), (6) detector phase (°), [(7) layer thickness]
assert((par.size() == 7) || (par.size() == 8));
@ -122,39 +94,33 @@ Double_t PPhotoMeissner::operator()(Double_t t, const vector<Double_t> &par) con
Double_t dz = 1.0; // go in 1A steps
Double_t zz = 0.0;
Double_t nn = 0.0;
Double_t sum = 0.0;
Double_t BB = 0.0;
Double_t gamma = fTwoPi*GAMMA_BAR_MUON;
int idx = fRgeHandler->GetEnergyIndex(par[0]);
int done = 0;
if (par.size() == 7) { // semi-infinite sample
do {
nn = fStartupHandler->GetRgeValue(par[0], zz);
nn = fRgeHandler->Get_n(idx, zz);
BB = FieldHalfSpace(zz, par);
result += nn*cos(gamma*BB*t+fDegToRad*par[6]);
zz += dz;
sum += nn;
if (nn == 0.0)
done++;
} while (done < 5);
} else { // film
do {
nn = fStartupHandler->GetRgeValue(par[0], zz);
nn = fRgeHandler->Get_n(idx, zz);
BB = FieldFilm(zz, par);
result += nn*cos(gamma*BB*t+fDegToRad*par[6]);
zz += dz;
sum += nn;
if (nn == 0.0)
done++;
} while (done < 5);
}
if (sum == 0.0) {
cerr << endl << ">> PPhotoMeissner::operator(): **PANIC ERROR** sum of RGE values == 0!" << endl;
assert(0);
}
return result/sum;
return result;
}
//--------------------------------------------------------------------------
@ -180,7 +146,7 @@ Double_t PPhotoMeissner::InuMinus(const Double_t nu, const Double_t x) const
*
* \param par
*/
Double_t PPhotoMeissner::FieldFilm(const Double_t z, const vector<Double_t> &par) const
Double_t PPhotoMeissner::FieldFilm(const Double_t z, const std::vector<Double_t> &par) const
{
// prevent dividing by zero
double lamL = par[3];
@ -217,7 +183,7 @@ Double_t PPhotoMeissner::FieldFilm(const Double_t z, const vector<Double_t> &par
*
* \param par
*/
Double_t PPhotoMeissner::FieldHalfSpace(const Double_t z, const vector<Double_t> &par) const
Double_t PPhotoMeissner::FieldHalfSpace(const Double_t z, const std::vector<Double_t> &par) const
{
// prevent dividing by zero
double lamL = par[3];

507
src/external/libPhotoMeissner/classes/PStartupHandler_PM.cpp vendored
View File

@ -1,507 +0,0 @@
/***************************************************************************
PStartupHandler_PM.cpp
Author: Andreas Suter
e-mail: andreas.suter@psi.ch
***************************************************************************/
/***************************************************************************
* Copyright (C) 2013-2021 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>
#include "PStartupHandler_PM.h"
//--------------------------------------------------------------------------
// Constructor
//--------------------------------------------------------------------------
/**
*
*/
PRgeHandler_PM::PRgeHandler_PM(const PStringVector &rgeDataPathList, const PDoubleVector &rgeDataEnergyList)
{
fIsValid = false;
fIsValid = LoadRgeData(rgeDataPathList, rgeDataEnergyList);
}
//--------------------------------------------------------------------------
// Destructor
//--------------------------------------------------------------------------
/**
*
*/
PRgeHandler_PM::~PRgeHandler_PM()
{
fRgeDataList.clear();
}
//--------------------------------------------------------------------------
// GetRgeEnergyIndex
//--------------------------------------------------------------------------
/**
*
* \param energy in (keV)
*/
Int_t PRgeHandler_PM::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 PRgeHandler_PM::GetRgeValue(const Int_t index, const Double_t dist)
{
Double_t rgeVal = 0.0;
UInt_t distIdx = static_cast<UInt_t>(dist/(fRgeDataList[index].stoppingDistance[1]-fRgeDataList[index].stoppingDistance[0]));
if (distIdx >= fRgeDataList[index].stoppingDistance.size()) {
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 PRgeHandler_PM::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 PRgeHandler_PM::LoadRgeData(const PStringVector &rgeDataPathList, const PDoubleVector &rgeDataEnergyList)
{
std::ifstream fin;
PRgeData_PM 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(), std::iostream::in);
if (!fin.is_open()) {
std::cerr << std::endl << "PRgeHandler_PM::LoadRgeData **ERROR**";
std::cerr << std::endl << " Could not open file " << rgeDataPathList[i].Data();
std::cerr << std::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;
}
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
ClassImpQ(PStartupHandler_PM)
//--------------------------------------------------------------------------
// Constructor
//--------------------------------------------------------------------------
/**
*
*/
PStartupHandler_PM::PStartupHandler_PM()
{
fIsValid = true;
fStartupFileFound = false;
fStartupFilePath = "";
// get default path (for the moment only linux like)
char startup_path_name[512];
char *home_path=0;
// check if the startup file is found in the current directory
strcpy(startup_path_name, "./photoMeissner_startup.xml");
if (StartupFileExists(startup_path_name)) {
fStartupFileFound = true;
fStartupFilePath = TString(startup_path_name);
} else { // startup file is not found in the current directory
std::cout << std::endl << ">> PStartupHandler_PM(): **WARNING** Couldn't find photoMeissner_startup.xml in the current directory, will try default one." << std::endl;
home_path = getenv("HOME");
snprintf(startup_path_name, sizeof(startup_path_name), "%s/.musrfit/external/photoMeissner_startup.xml", home_path);
if (StartupFileExists(startup_path_name)) {
fStartupFileFound = true;
fStartupFilePath = TString(startup_path_name);
}
}
// init RGE handler
fRgeHandler = nullptr;
}
//--------------------------------------------------------------------------
// OnStartDocument
//--------------------------------------------------------------------------
/**
* <p>
*/
void PStartupHandler_PM::OnStartDocument()
{
fKey = eEmpty;
}
//--------------------------------------------------------------------------
// OnEndDocument
//--------------------------------------------------------------------------
/**
* <p>
*/
void PStartupHandler_PM::OnEndDocument()
{
if (!fIsValid)
return;
// generate the file path list
TString str;
for (unsigned int i=0; i<fRgeDataEnergyList.size(); i++) {
str = fRgePath;
str += fRgeDataEnergyList[i];
str += ".rge";
fRgeFilePathList.push_back(str);
}
fRgeHandler = new PRgeHandler_PM(fRgeFilePathList, fRgeDataEnergyList);
if (fRgeHandler == nullptr) { // severe problem
fIsValid = false;
std::cerr << std::endl << ">> PStartupHandler_PM::OnEndDocument(): **ERROR** couldn't invoke RGE handler." << std::endl << std::endl;
return;
}
if (!fRgeHandler->IsValid()) { // severe problem
fIsValid = false;
std::cerr << std::endl << ">> PStartupHandler_PM::OnEndDocument(): **ERROR** RGE handler not valid." << std::endl << std::endl;
return;
}
}
//--------------------------------------------------------------------------
// OnStartElement
//--------------------------------------------------------------------------
/**
* <p>
*
* \param str
* \param attributes
*/
void PStartupHandler_PM::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 PStartupHandler_PM::OnEndElement(const char *str)
{
fKey = eEmpty;
}
//--------------------------------------------------------------------------
// OnCharacters
//--------------------------------------------------------------------------
/**
* <p>
*
* \param str
*/
void PStartupHandler_PM::OnCharacters(const char *str)
{
TString tstr;
Double_t dval;
switch (fKey) {
case eDataPath:
fRgePath = str;
break;
case eEnergy:
tstr = str;
if (tstr.IsFloat()) {
dval = tstr.Atof();
fRgeDataEnergyList.push_back(dval);
} else {
std::cerr << std::endl << "PStartupHandler_PM::OnCharacters: **ERROR** when finding energy:";
std::cerr << std::endl << "\"" << str << "\" is not a floating point number, will ignore it.";
std::cerr << std::endl;
}
break;
default:
break;
}
}
//--------------------------------------------------------------------------
// OnComment
//--------------------------------------------------------------------------
/**
* <p>
*
* \param str
*/
void PStartupHandler_PM::OnComment(const char *str)
{
// nothing to be done for now
}
//--------------------------------------------------------------------------
// OnWarning
//--------------------------------------------------------------------------
/**
* <p>
*
* \param str
*/
void PStartupHandler_PM::OnWarning(const char *str)
{
std::cerr << std::endl << "PStartupHandler_PM **WARNING** " << str;
std::cerr << std::endl;
}
//--------------------------------------------------------------------------
// OnError
//--------------------------------------------------------------------------
/**
* <p>
*
* \param str
*/
void PStartupHandler_PM::OnError(const char *str)
{
std::cerr << std::endl << "PStartupHandler_PM **ERROR** " << str;
std::cerr << std::endl;
}
//--------------------------------------------------------------------------
// OnFatalError
//--------------------------------------------------------------------------
/**
* <p>
*
* \param str
*/
void PStartupHandler_PM::OnFatalError(const char *str)
{
std::cerr << std::endl << "PStartupHandler_PM **FATAL ERROR** " << str;
std::cerr << std::endl;
}
//--------------------------------------------------------------------------
// OnCdataBlock
//--------------------------------------------------------------------------
/**
* <p>
*
* \param str
*/
void PStartupHandler_PM::OnCdataBlock(const char *str, Int_t len)
{
// nothing to be done for now
}
//--------------------------------------------------------------------------
// GetRgeEnergyIndex
//--------------------------------------------------------------------------
/**
* <p>
*
* \param energy
*/
Int_t PStartupHandler_PM::GetRgeEnergyIndex(const Double_t energy)
{
Int_t result = -1;
if (fIsValid)
result = fRgeHandler->GetRgeEnergyIndex(energy);
return result;
}
//--------------------------------------------------------------------------
// GetRgeValue
//--------------------------------------------------------------------------
/**
* <p>
*
* \param energy
* \param dist
*/
Double_t PStartupHandler_PM::GetRgeValue(const Double_t energy, const Double_t dist)
{
Double_t result = 0.0;
if (fIsValid)
result = fRgeHandler->GetRgeValue(energy, dist);
return result;
}
//--------------------------------------------------------------------------
// GetRgeValue
//--------------------------------------------------------------------------
/**
* <p>
*
* \param index
* \param dist
*/
Double_t PStartupHandler_PM::GetRgeValue(const Int_t index, const Double_t dist)
{
Double_t result = 0.0;
if (fIsValid)
result = fRgeHandler->GetRgeValue(index, dist);
return result;
}
//--------------------------------------------------------------------------
// StartupFileExists
//--------------------------------------------------------------------------
/**
* <p>
*
*/
bool PStartupHandler_PM::StartupFileExists(char *fln)
{
bool result = false;
std::ifstream ifile(fln);
if (ifile.fail()) {
result = false;
} else {
result = true;
ifile.close();
}
return result;
}
// -------------------------------------------------------------------------
// end
// -------------------------------------------------------------------------

8
src/external/libPhotoMeissner/include/PPhotoMeissner.h vendored
View File

@ -33,7 +33,7 @@
#include <vector>
#include "PUserFcnBase.h"
#include "PStartupHandler_PM.h"
#include "PRgeHandler.h"
class PPhotoMeissner : public PUserFcnBase
{
@ -44,12 +44,12 @@ class PPhotoMeissner : public PUserFcnBase
virtual Bool_t IsValid() { return fValid; }
// function operator
Double_t operator()(Double_t, const std::vector<Double_t>&) const;
Double_t operator()(Double_t t, const std::vector<Double_t> &param) const;
private:
PStartupHandler_PM *fStartupHandler;
PRgeHandler *fRgeHandler{nullptr};
Bool_t fValid; ///< flag indicating if initialization went through smoothly
Bool_t fValid{true}; ///< flag indicating if initialization went through smoothly
constexpr static const Double_t fDegToRad = 0.0174532925199432955;
constexpr static const Double_t fTwoPi = 6.28318530717958623;

7
src/external/libPhotoMeissner/test/photoMeissner_startup.xml vendored
View File

@ -3,8 +3,9 @@
<comment>
contains the needed trim.sp meta-file information
</comment>
<trim_sp_part>
<data_path>./trimsp/YBCO_E</data_path>
<trim_sp>
<data_path>./trimsp</data_path>
<rge_fln_pre>YBCO_E</rge_fln_pre>
<energy_list>
<energy>1000</energy>
<energy>3000</energy>
@ -16,5 +17,5 @@
<energy>20300</energy>
<energy>24300</energy>
</energy_list>
</trim_sp_part>
</trim_sp>
</photo_meissner>

133
src/external/libPhotoMeissner/include/PStartupHandler_PM.h → src/include/PRgeHandler.h
View File

@ -1,6 +1,6 @@
/***************************************************************************
PStartupHandler_PM.h
PRgeHandler.h
Author: Andreas Suter
e-mail: andreas.suter@psi.ch
@ -8,7 +8,7 @@
***************************************************************************/
/***************************************************************************
* Copyright (C) 2013-2021 by Andreas Suter *
* Copyright (C) 2007-2021 by Andreas Suter *
* andreas.suter@psi.ch *
* *
* This program is free software; you can redistribute it and/or modify *
@ -27,74 +27,46 @@
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
***************************************************************************/
#ifndef _PSTARTUPHANDLER_PM_H_
#define _PSTARTUPHANDLER_PM_H_
#ifndef _PRGEHANDLER_H_
#define _PRGEHANDLER_H_
#include <string>
#include <vector>
#include <TObject.h>
#include <TQObject.h>
#include <TString.h>
#include <TSAXParser.h>
//-------------------------------------------------------------
/**
* <p> typedef to make to code more readable.
*/
typedef std::vector<TString> PStringVector;
#include "PMusr.h"
//-------------------------------------------------------------
//-----------------------------------------------------------------------------
/**
* <p> typedef to make to code more readable.
*/
typedef std::vector<Double_t> PDoubleVector;
//-------------------------------------------------------------
/**
* <p>
* <p>Keep a single rge table from TrimSP for a given energy.
*/
typedef struct {
Double_t energy;
PDoubleVector stoppingDistance;
PDoubleVector stoppingAmplitude;
} PRgeData_PM;
PDoubleVector depth;
PDoubleVector amplitude;
PDoubleVector nn; // normalized int n(z) dz = 1 amplitudes
Double_t noOfParticles;
} PRgeData;
//-------------------------------------------------------------
//-----------------------------------------------------------------------------
/**
* <p>
* <p>Keep all rge tables from TrimSP.
*/
typedef std::vector<PRgeData_PM> PRgeData_PM_List;
typedef std::vector<PRgeData> PRgeDataList;
//-------------------------------------------------------------
//-----------------------------------------------------------------------------
/**
* <p>
* <p>Reads the xml-startup file in order to extract all necessary information
* about the RGE files (TrimSP) needed.
*/
class PRgeHandler_PM
class PXmlRgeHandler : public TObject, public TQObject
{
public:
PRgeHandler_PM(const PStringVector &rgeDataPathList, const PDoubleVector &rgeDataEnergyList);
virtual ~PRgeHandler_PM();
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;
PRgeData_PM_List fRgeDataList;
virtual Bool_t LoadRgeData(const PStringVector &rgeDataPathList, const PDoubleVector &rgeDataEnergyList);
};
//-------------------------------------------------------------
/**
* <p>
*/
class PStartupHandler_PM : public TObject
{
public:
PStartupHandler_PM();
virtual ~PStartupHandler_PM() {}
PXmlRgeHandler() {}
virtual ~PXmlRgeHandler() {}
virtual void OnStartDocument(); // SLOT
virtual void OnEndDocument(); // SLOT
@ -108,32 +80,49 @@ class PStartupHandler_PM : public TObject
virtual void OnCdataBlock(const char*, Int_t); // SLOT
virtual bool IsValid() { return fIsValid; }
virtual TString GetStartupFilePath() { return fStartupFilePath; }
virtual bool StartupFileFound() { return fStartupFileFound; }
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);
virtual std::string GetTrimSpDataPath() { return fTrimSpDataPath; }
virtual std::string GetTrimSpFlnPre() { return fTrimSpFlnPre; }
virtual const PIntVector GetTrimSpDataVectorList() const { return fTrimSpDataEnergyList; }
private:
enum EKeyWords {eEmpty, eDataPath, eEnergy};
enum EKeyWords {eEmpty, eDataPath, eFlnPre, eEnergy};
EKeyWords fKey;
bool fIsValid;
bool isTrimSp{false};
bool fIsValid{true};
bool fStartupFileFound;
TString fStartupFilePath;
std::string fTrimSpDataPath{""}; //< where to find the rge-files
std::string fTrimSpFlnPre{""}; //< keeps the preface of the rge file name, e.g. LCCO_E
PIntVector fTrimSpDataEnergyList; //< TrimSP implantation energy list (in eV)
TString fRgePath;
PStringVector fRgeFilePathList;
PDoubleVector fRgeDataEnergyList;
PRgeHandler_PM *fRgeHandler;
bool StartupFileExists(char *fln);
ClassDef(PStartupHandler_PM, 1)
ClassDef(PXmlRgeHandler, 1)
};
#endif // _PSTARTUPHANDLER_PM_H_
//-----------------------------------------------------------------------------
/**
* @brief The PRegHandler class
*/
class PRgeHandler : public TObject
{
public:
PRgeHandler(std::string fln="");
virtual ~PRgeHandler() {}
virtual bool IsValid() { return fValid; }
virtual PRgeDataList GetRgeData() { return fData; }
virtual Double_t GetZmax(const Double_t energy);
virtual Double_t GetZmax(const Int_t idx);
virtual Double_t Get_n(const Double_t energy, const Double_t z);
virtual Double_t Get_n(const Int_t idx, const Double_t z);
virtual Int_t GetEnergyIndex(const Double_t energy);
private:
bool fValid{false};
PRgeDataList fData;
virtual bool ReadRgeFile(const std::string fln, PRgeData &data);
ClassDef(PRgeHandler, 1)
};
#endif // _PRGEHANDLER_H_

7
src/external/libPhotoMeissner/include/PStartupHandler_PMLinkDef.h → src/include/PRgeHandlerLinkDef.h
View File

@ -1,6 +1,6 @@
/***************************************************************************
PStartupHandler_PMLinkDef.h
PRgeHandlerLinkDef.h
Author: Andreas Suter
e-mail: andreas.suter@psi.ch
@ -8,7 +8,7 @@
***************************************************************************/
/***************************************************************************
* Copyright (C) 2013-2021 by Andreas Suter *
* Copyright (C) 2007-2021 by Andreas Suter *
* andreas.suter@psi.ch *
* *
* This program is free software; you can redistribute it and/or modify *
@ -33,6 +33,7 @@
#pragma link off all classes;
#pragma link off all functions;
#pragma link C++ class PStartupHandler_PM+;
#pragma link C++ class PXmlRgeHandler+;
#pragma link C++ class PRgeHandler+;
#endif