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

added the option to split a user function into a global and run-block related part (see MUSR-134)

Browse Source
This commit is contained in:
nemu
2010-11-10 13:13:00 +00:00
parent c4dfc3cbce
commit f4d6e349fe
41 changed files with 2401 additions and 161 deletions
Download Patch File Download Diff File Expand all files Collapse all files

85
src/tests/userFcnGlobal/1951_tf_h13.msr Normal file
View File

@@ -0,0 +1,85 @@
YBCO(110)/PrBCO/STO on Ag, T=4.96 K, E=19.12 keV, B=~98(G)/3.04(A), Tr/Sa=15.02/-5.00 kV, RA-LR 0.05 kV, RA-TB 0.14 kV
###############################################################
FITPARAMETER
# No Name Value Step Pos_Error Boundaries
1 phase 25.1328 -0.979855 0.987925 -180 180
2 field 84.2201 -0.243991 0.264511 0 none
3 asym 0.11616 -0.00735602 0.00786111 0 0.3
4 rate 0.297241 -0.028157 0.0289587 0 2
5 fieldBKG 94.9012 -0.142071 0.143479 0 none
6 asymBKG 0.139904 -0.00933264 0.00905566 0 0.3
7 rateBKG 0.0623834 -0.0189223 0.0188684 0 100
8 Norm_L 272.022 -0.565339 0.566487 0 1000
9 BG_L 25.7891 -0.120907 0.120769 0 100
10 Norm_R 320.245 -0.608216 0.60492 0 1000
11 BG_R 30.7947 -0.130858 0.131547 0 100
12 relasy_R 0.982832 -0.0159124 0.0162089 0 1.2
13 relphase_R 159.584 -0.898672 0.898169 150 210
14 one 1 0 none
15 zero 0 0 none
###############################################################
THEORY
asymmetry fun1
userFcn libPRelax.so PGauss 4 (rate)
TFieldCos fun3 fun4 (phase frequency)
+
asymmetry fun2
userFcn libPRelax.so PExp 7 (rate)
TFieldCos fun3 fun5 (phase frequency)
###############################################################
FUNCTIONS
fun1 = par3 * map1
fun2 = par6 * map1
fun3 = par1 + map2
fun4 = par2 * gamma_mu
fun5 = par5 * gamma_mu
###############################################################
RUN data/lem10_his_1951 MUE4 PSI ROOT-NPP (name beamline institute data-file-format)
fittype 0 (single histogram fit)
norm 8
backgr.fit 9
lifetimecorrection
map 14 15 0 0 0 0 0 0 0 0
forward 1
background 65000 66500
data 3289 63000
fit 0.15 8
packing 200
RUN data/lem10_his_1951 MUE4 PSI ROOT-NPP (name beamline institute data-file-format)
fittype 0 (single histogram fit)
norm 10
backgr.fit 11
lifetimecorrection
map 12 13 0 0 0 0 0 0 0 0
forward 3
background 65000 66500
data 3289 63000
fit 0.15 8
packing 200
###############################################################
COMMANDS
MINIMIZE
MINOS
SAVE
###############################################################
FOURIER
units Gauss # units either 'Gauss', 'MHz', or 'Mc/s'
fourier_power 10
apodization STRONG # NONE, WEAK, MEDIUM, STRONG
plot POWER # REAL, IMAG, REAL_AND_IMAG, POWER, PHASE
range 0 250
###############################################################
PLOT 0 (single histo plot)
runs 1 2
range 0.00 8.00 -0.30 0.30
###############################################################
STATISTIC --- 2010-11-09 20:13:52
chisq = 401.71526823978365, NDF = 389, chisq/NDF = 1.0326870648837627

99
src/tests/userFcnGlobal/Makefile Normal file
View File

@@ -0,0 +1,99 @@
#---------------------------------------------------
# Makefile
#
# Author: Andreas Suter
# e-mail: andreas.suter@psi.ch
#
# $Id: Makefile.libPNL_PippardFitter 4314 2009-12-22 17:03:48Z l_wojek $
#
#---------------------------------------------------
#---------------------------------------------------
# get compilation and library flags from root-config
ROOTCFLAGS = $(shell $(ROOTSYS)/bin/root-config --cflags)
ROOTLIBS = $(shell $(ROOTSYS)/bin/root-config --libs)
ROOTGLIBS = $(shell $(ROOTSYS)/bin/root-config --glibs)
#---------------------------------------------------
# depending on the architecture, choose the compiler,
# linker, and the flags to use
#
OSTYPE = $(shell uname)
ifeq ($(OSTYPE),Linux)
OS = LINUX
endif
ifeq ($(OSTYPE),Linux-gnu)
OS = LINUX
endif
ifeq ($(OSTYPE),Darwin)
OS = DARWIN
endif
# -- Linux
ifeq ($(OS),LINUX)
CXX = g++
CXXFLAGS = -Wall -Wno-trigraphs -fPIC
INCLUDES = -I../include
LD = g++
LDFLAGS = -g
SOFLAGS = -O -shared
endif
# -- Darwin
ifeq ($(OS),DARWIN)
CXX = g++
CXXFLAGS = -Wall -Wno-trigraphs -fPIC
INCLUDES = -I../include
LD = g++
LDFLAGS = -g
SOFLAGS = -dynamic
endif
# the output from the root-config script:
CXXFLAGS += $(ROOTCFLAGS)
LDFLAGS +=
# the ROOT libraries (G = graphic)
LIBS = $(ROOTLIBS) -lXMLParser
GLIBS = $(ROOTGLIBS) -lXMLParser
# some definitions: headers (used to generate *Dict* stuff), sources, objects,...
OBJS =
OBJS += PRelax.o PRelaxDict.o
SHLIB = libPRelax.so
# make the shared lib:
#
all: $(SHLIB)
$(SHLIB): $(OBJS)
@echo "---> Building shared library $(SHLIB) ..."
/bin/rm -f $(SHLIB)
$(LD) $(OBJS) $(SOFLAGS) -o $(SHLIB) $(LIBS)
@echo "done"
# clean up: remove all object file (and core files)
# semicolon needed to tell make there is no source
# for this target!
#
clean:; @rm -f $(OBJS) *Dict* core*
@echo "---> removing $(OBJS)"
#
$(OBJS): %.o: %.cpp
$(CXX) $(INCLUDES) $(CXXFLAGS) -c $<
PRelaxDict.cpp: PRelax.h PRelaxLinkDef.h
@echo "Generating dictionary $@..."
rootcint -f $@ -c -p $^
install: all
@echo "Installing shared lib: libPRelax.so ( you must be root ;-) )"
ifeq ($(OS),LINUX)
cp -pv $(SHLIB) $(ROOTSYS)/lib
cp -pv PRelax.h $(ROOTSYS)/include
endif

338
src/tests/userFcnGlobal/PRelax.cpp Normal file
View File

@@ -0,0 +1,338 @@
/***************************************************************************
PRelax.cpp
Author: Andreas Suter
e-mail: andreas.suter@psi.ch
$Id$
***************************************************************************/
/***************************************************************************
* Copyright (C) 2007 by Andreas Suter *
* andreas.suter@psi.c *
* *
* 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 <iostream>
using namespace std;
#include <cassert>
#include <cstdlib>
#include <cmath>
#include "PRelax.h"
ClassImp(PExp)
//------------------------------------------------------------------------------------
/**
* <p> Constructor.
*/
PExp::PExp()
{
fValid = false;
fInvokedGlobal = false;
fIdxGlobal = -1;
fGlobalUserFcn = 0;
}
//------------------------------------------------------------------------------------
/**
* <p> Destructor.
*/
PExp::~PExp()
{
if ((fGlobalUserFcn != 0) && fInvokedGlobal) {
delete fGlobalUserFcn;
fGlobalUserFcn = 0;
}
}
//------------------------------------------------------------------------------------
/**
* <p> Used to invoke/retrieve the proper global user function
*
* \param globalPart reference to the global user function object vector
* \param idx global user function index within the theory tree
*/
void PExp::SetGlobalPart(vector<void *> &globalPart, UInt_t idx)
{
fIdxGlobal = static_cast<Int_t>(idx);
if ((Int_t)globalPart.size() <= fIdxGlobal) { // global user function not present, invoke it
fGlobalUserFcn = new PExpGlobal();
if (fGlobalUserFcn == 0) { // global user function object couldn't be invoked -> error
fValid = false;
cerr << endl << ">> PExp::SetGlobalPart(): **ERROR** Couldn't invoke global user function object, sorry ..." << endl;
} else { // global user function object could be invoked -> resize to global user function vector and keep the pointer to the corresponding object
globalPart.resize(fIdxGlobal+1);
globalPart[fIdxGlobal] = dynamic_cast<PExpGlobal*>(fGlobalUserFcn);
fValid = true;
fInvokedGlobal = true;
}
} else { // global user function already present hence just retrieve a pointer to it
fValid = true;
fGlobalUserFcn = (PExpGlobal*)globalPart[fIdxGlobal];
}
}
//------------------------------------------------------------------------------------
/**
* <p> Used to check if the user function is OK.
*
* <b>return:</b> true if both, the user function and the global user function object are valid
*/
Bool_t PExp::GlobalPartIsValid() const
{
return (fValid && fGlobalUserFcn->IsValid());
}
//------------------------------------------------------------------------------------
/**
* <p> user function example: global exponential
*
* \f[ = exp(-t \lambda) \f]
*
* <b>meaning of paramValues:</b> \f$\lambda\f$ depol. rate
*
* <b>return:</b> function value
*
* \param t time in \f$(\mu\mathrm{s})\f$, or x-axis value for non-muSR fit
* \param param parameter vector
*/
Double_t PExp::operator()(Double_t t, const std::vector<Double_t> &param) const
{
// expected parameters: lambda
assert(param.size() == 1);
assert(fGlobalUserFcn);
// call the global user function object
fGlobalUserFcn->Calc(param);
return exp(-t*param[0]);
}
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
ClassImp(PExpGlobal)
//------------------------------------------------------------------------------------
/**
* <p> Constructor
*/
PExpGlobal::PExpGlobal()
{
fValid = true;
fRandom = 0.0;
}
//------------------------------------------------------------------------------------
/**
* <p> Destructor
*/
PExpGlobal::~PExpGlobal()
{
fPrevParam.clear();
}
//------------------------------------------------------------------------------------
/**
* <p> Check if the parameter vector has changed to the previous one, and if so,
* do some global user function calculation (here calculate only a stupid random
* number).
*
* \param param fit parameter vector
*/
void PExpGlobal::Calc(const std::vector<Double_t> &param)
{
// check if previous parameters where ever fed
if (fPrevParam.size() == 0) {
fPrevParam = param;
fRandom = rand() / (RAND_MAX+1.0);
cout << endl << "debug> PExpGlobal::Calc(): fRandom = " << fRandom << " (first time)";
}
// check if param has changed, i.e. something needs to be done
Bool_t newParam = false;
for (UInt_t i=0; i<param.size(); i++) {
if (param[i] != fPrevParam[i]) {
newParam = true;
break;
}
}
if (newParam) {
// do the global calculations (here only a stupid random number)
fRandom = rand() / (RAND_MAX+1.0);
// feed the previous parameter vector
for (UInt_t i=0; i<param.size(); i++)
fPrevParam[i] = param[i];
}
}
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
ClassImp(PGauss)
//------------------------------------------------------------------------------------
/**
* <p> Constructor.
*/
PGauss::PGauss()
{
fValid = false;
fInvokedGlobal = false;
fIdxGlobal = -1;
fGlobalUserFcn = 0;
}
//------------------------------------------------------------------------------------
/**
* <p> Destructor.
*/
PGauss::~PGauss()
{
if ((fGlobalUserFcn != 0) && fInvokedGlobal) {
delete fGlobalUserFcn;
fGlobalUserFcn = 0;
}
}
//------------------------------------------------------------------------------------
/**
* <p> Used to invoke/retrieve the proper global user function
*
* \param globalPart reference to the global user function object vector
* \param idx global user function index within the theory tree
*/
void PGauss::SetGlobalPart(vector<void *> &globalPart, UInt_t idx)
{
fIdxGlobal = static_cast<Int_t>(idx);
if ((Int_t)globalPart.size() <= fIdxGlobal) {
fGlobalUserFcn = new PGaussGlobal();
if (fGlobalUserFcn == 0) {
fValid = false;
cerr << endl << ">> PGauss::SetGlobalPart(): **ERROR** Couldn't invoke global user function object, sorry ..." << endl;
} else {
globalPart.resize(fIdxGlobal+1);
globalPart[fIdxGlobal] = dynamic_cast<PGaussGlobal*>(fGlobalUserFcn);
fValid = true;
fInvokedGlobal = true;
}
} else {
fValid = true;
fGlobalUserFcn = (PGaussGlobal*)globalPart[fIdxGlobal];
}
}
//------------------------------------------------------------------------------------
/**
* <p> Used to check if the user function is OK.
*
* <b>return:</b> true if both, the user function and the global user function object are valid
*/
Bool_t PGauss::GlobalPartIsValid() const
{
return (fValid && fGlobalUserFcn->IsValid());
}
//------------------------------------------------------------------------------------
/**
* <p> user function example: global Gaussian
*
* \f[ = exp(-1/2 (t sigma)^2) \f]
*
* <b>meaning of paramValues:</b> \f$\sigma\f$ depol. rate
*
* <b>return:</b> function value
*
* \param t time in \f$(\mu\mathrm{s})\f$, or x-axis value for non-muSR fit
* \param param parameter vector
*/
Double_t PGauss::operator()(Double_t t, const std::vector<Double_t> &param) const
{
// expected parameters: sigma
assert(param.size() == 1);
assert(fGlobalUserFcn);
fGlobalUserFcn->Calc(param);
return exp(-0.5*pow(t*param[0],2.0));
}
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
ClassImp(PGaussGlobal)
//------------------------------------------------------------------------------------
/**
* <p> Constructor
*/
PGaussGlobal::PGaussGlobal()
{
fValid = true;
fRandom = 0.0;
}
//------------------------------------------------------------------------------------
/**
* <p> Destructor
*/
PGaussGlobal::~PGaussGlobal()
{
fPrevParam.clear();
}
//------------------------------------------------------------------------------------
/**
* <p> Check if the parameter vector has changed to the previous one, and if so,
* do some global user function calculation (here calculate only a stupid random
* number).
*
* \param param fit parameter vector
*/
void PGaussGlobal::Calc(const std::vector<Double_t> &param)
{
// check if previous parameters where ever fed
if (fPrevParam.size() == 0) {
fPrevParam = param;
fRandom = rand() / (RAND_MAX+1.0);
cout << endl << "debug> PGaussGlobal::Calc(): fRandom = " << fRandom << " (first time)";
}
// check if param has changed, i.e. something needs to be done
Bool_t newParam = false;
for (UInt_t i=0; i<param.size(); i++) {
if (param[i] != fPrevParam[i]) {
newParam = true;
break;
}
}
if (newParam) {
// do the global calculations (here only a stupid random number)
fRandom = rand() / (RAND_MAX+1.0);
// feed the previous parameter vector
for (UInt_t i=0; i<param.size(); i++)
fPrevParam[i] = param[i];
}
}

135
src/tests/userFcnGlobal/PRelax.h Normal file
View File

@@ -0,0 +1,135 @@
/***************************************************************************
PRelax.h
Author: Andreas Suter
e-mail: andreas.suter@psi.ch
$Id$
***************************************************************************/
/***************************************************************************
* Copyright (C) 2007 by Andreas Suter *
* andreas.suter@psi.c *
* *
* 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 _PRELAX_H_
#define _PRELAX_H_
#include <vector>
#include <iostream>
#include "PUserFcnBase.h"
class PExpGlobal;
class PGaussGlobal;
//-----------------------------------------------------------------------------------
/**
* <p>User function example class. Global exponential function.
*/
class PExp : public PUserFcnBase
{
public:
PExp();
virtual ~PExp();
virtual Bool_t NeedGlobalPart() const { return true; }
virtual void SetGlobalPart(vector<void*> &globalPart, UInt_t idx);
virtual Bool_t GlobalPartIsValid() const;
Double_t operator()(Double_t t, const std::vector<Double_t> &param) const;
private:
Bool_t fValid;
Bool_t fInvokedGlobal;
Int_t fIdxGlobal;
PExpGlobal *fGlobalUserFcn;
ClassDef(PExp, 1)
};
//-----------------------------------------------------------------------------------
/**
* <p>User function example class (global part). Global exponential function.
*/
class PExpGlobal
{
public:
PExpGlobal();
virtual ~PExpGlobal();
Bool_t IsValid() { return fValid; }
void Calc(const std::vector<Double_t> &param);
private:
Bool_t fValid;
std::vector<Double_t> fPrevParam;
Double_t fRandom;
ClassDef(PExpGlobal, 1)
};
//-----------------------------------------------------------------------------------
/**
* <p>User function example class. Global Gaussian function.
*/
class PGauss : public PUserFcnBase
{
public:
PGauss();
virtual ~PGauss();
virtual Bool_t NeedGlobalPart() const { return true; }
virtual void SetGlobalPart(vector<void*> &globalPart, UInt_t idx);
virtual Bool_t GlobalPartIsValid() const;
Double_t operator()(Double_t t, const std::vector<Double_t> &param) const;
private:
Bool_t fValid;
Bool_t fInvokedGlobal;
Int_t fIdxGlobal;
PGaussGlobal *fGlobalUserFcn;
ClassDef(PGauss, 1)
};
//-----------------------------------------------------------------------------------
/**
* <p>User function example class (global part). Global Gaussian function.
*/
class PGaussGlobal
{
public:
PGaussGlobal();
virtual ~PGaussGlobal();
Bool_t IsValid() { return fValid; }
void Calc(const std::vector<Double_t> &param);
private:
Bool_t fValid;
std::vector<Double_t> fPrevParam;
Double_t fRandom;
ClassDef(PGaussGlobal, 1)
};
#endif // _PRELAX_H_

12
src/tests/userFcnGlobal/PRelaxLinkDef.h Normal file
View File

@@ -0,0 +1,12 @@
#ifdef __CINT__
#pragma link off all globals;
#pragma link off all classes;
#pragma link off all functions;
#pragma link C++ class PExp+;
#pragma link C++ class PExpGlobal+;
#pragma link C++ class PGauss+;
#pragma link C++ class PGaussGlobal+;
#endif

BIN
src/tests/userFcnGlobal/data/lem10_his_1951.root Normal file
View File

Binary file not shown.