Merge branch 'master' into root6
This commit is contained in:
commit
b1ba14ec76
@ -1242,6 +1242,7 @@ AC_CONFIG_FILES([Makefile \
|
||||
src/external/libSpinValve/classes/Makefile \
|
||||
src/external/libPhotoMeissner/Makefile \
|
||||
src/external/libPhotoMeissner/classes/Makefile \
|
||||
src/external/libGbGLF/Makefile \
|
||||
src/external/libBNMR/Makefile \
|
||||
src/musredit_qt5/Makefile \
|
||||
src/musredit/Makefile \
|
||||
|
||||
@ -301,6 +301,7 @@ PTheory::PTheory(PMsrHandler *msrInfo, UInt_t runNo, const Bool_t hasParent) : f
|
||||
cerr << endl << ">> PTheory::PTheory: **ERROR** user function object could not be invoked. See line no " << line->fLineNo;
|
||||
cerr << endl;
|
||||
fValid = false;
|
||||
return;
|
||||
} else { // user function valid, hence expand the fUserParam vector to the proper size
|
||||
fUserParam.resize(fParamNo.size());
|
||||
}
|
||||
|
||||
3
src/external/Makefile.am
vendored
3
src/external/Makefile.am
vendored
@ -3,7 +3,8 @@ if BUILD_ASLIBS
|
||||
ASDIRS = Nonlocal \
|
||||
MagProximity \
|
||||
libSpinValve \
|
||||
libPhotoMeissner
|
||||
libPhotoMeissner \
|
||||
libGbGLF
|
||||
endif
|
||||
|
||||
if BUILD_CUBALIB
|
||||
|
||||
56
src/external/libGbGLF/Makefile.am
vendored
Normal file
56
src/external/libGbGLF/Makefile.am
vendored
Normal file
@ -0,0 +1,56 @@
|
||||
## Process this file with automake to create Makefile.in
|
||||
|
||||
h_sources = \
|
||||
PGbGLF.h
|
||||
|
||||
h_linkdef = \
|
||||
PGbGLFLinkDef.h
|
||||
|
||||
dict_h_sources = \
|
||||
PGbGLFDict.h
|
||||
|
||||
cpp_sources = \
|
||||
PGbGLF.cpp
|
||||
|
||||
dict_cpp_sources = \
|
||||
PGbGLFDict.cpp
|
||||
|
||||
include_HEADERS = $(h_sources)
|
||||
noinst_HEADERS = $(h_linkdef) $(dict_h_sources)
|
||||
|
||||
AM_CPPFLAGS = -I$(top_srcdir)/src/include $(PMUSR_CFLAGS) $(FFTW3_CFLAGS) $(GSL_CFLAGS) -I$(ROOTINCDIR)
|
||||
AM_CXXFLAGS = $(LOCAL_LIB_CXXFLAGS)
|
||||
|
||||
BUILT_SOURCES = $(dict_cpp_sources) $(dict_h_sources)
|
||||
AM_LDFLAGS = $(LOCAL_LIB_LDFLAGS) -L@ROOTLIBDIR@
|
||||
CLEANFILES = *Dict.cpp *Dict.h *~ core
|
||||
|
||||
%Dict.cpp %Dict.h: %.h %LinkDef.h
|
||||
@ROOTCINT@ -v -f $*Dict.cpp -c -p $(AM_CPPFLAGS) $^
|
||||
|
||||
lib_LTLIBRARIES = libGbGLF.la
|
||||
|
||||
libGbGLF_la_SOURCES = $(h_sources) $(cpp_sources) $(dict_h_sources) $(dict_cpp_sources)
|
||||
libGbGLF_la_LIBADD = $(USERFCN_LIBS) $(GSL_LIBS)
|
||||
libGbGLF_la_LDFLAGS = -version-info $(PLUGIN_LIBRARY_VERSION) -release $(PLUGIN_RELEASE) $(AM_LDFLAGS)
|
||||
|
||||
## For the moment do not build pkgconfig files for musrfit plug-ins...
|
||||
## pkgconfigdir = $(libdir)/pkgconfig
|
||||
## pkgconfig_DATA = LFRelaxation.pc
|
||||
|
||||
## However, create some symbolic links to the shared library
|
||||
## in order to unify the function call on different operating systems
|
||||
|
||||
if IS_DARWIN
|
||||
install-exec-hook:
|
||||
$(LN_S) -f $(libdir)/libGbGLF.dylib $(libdir)/libGbGLF.so
|
||||
uninstall-hook:
|
||||
rm -f $(libdir)/libGbGLF.so
|
||||
endif
|
||||
|
||||
if IS_CYGWIN
|
||||
install-exec-hook:
|
||||
$(LN_S) -f $(bindir)/cygGbGLF-$(PLUGIN_MAJOR_VERSION)-$(PLUGIN_MINOR_VERSION)-$(PLUGIN_MAJOR_VERSION).dll $(libdir)/libGbGLF.so
|
||||
uninstall-hook:
|
||||
rm -f $(libdir)/libGbGLF.so
|
||||
endif
|
||||
173
src/external/libGbGLF/PGbGLF.cpp
vendored
Normal file
173
src/external/libGbGLF/PGbGLF.cpp
vendored
Normal file
@ -0,0 +1,173 @@
|
||||
/***************************************************************************
|
||||
|
||||
PGbGLF.cpp
|
||||
|
||||
Author: Andreas Suter
|
||||
e-mail: andreas.suter@psi.ch
|
||||
|
||||
***************************************************************************/
|
||||
|
||||
/***************************************************************************
|
||||
* Copyright (C) 2007-2016 by Andreas Suter *
|
||||
* andreas.suter@psi.ch *
|
||||
* *
|
||||
* This program is free software; you can redistribute it and/or modify *
|
||||
* it under the terms of the GNU General Public License as published by *
|
||||
* the Free Software Foundation; either version 2 of the License, or *
|
||||
* (at your option) any later version. *
|
||||
* *
|
||||
* This program is distributed in the hope that it will be useful, *
|
||||
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
|
||||
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
|
||||
* GNU General Public License for more details. *
|
||||
* *
|
||||
* You should have received a copy of the GNU General Public License *
|
||||
* along with this program; if not, write to the *
|
||||
* Free Software Foundation, Inc., *
|
||||
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
|
||||
***************************************************************************/
|
||||
|
||||
#ifdef HAVE_CONFIG_H
|
||||
#include "config.h"
|
||||
#endif
|
||||
|
||||
#include <cassert>
|
||||
#include <iostream>
|
||||
#include <cmath>
|
||||
using namespace std;
|
||||
|
||||
#include "PMusr.h"
|
||||
#include "PGbGLF.h"
|
||||
|
||||
#define TWO_PI 6.28318530717958647692528676656
|
||||
|
||||
ClassImp(PGbGLF)
|
||||
|
||||
//--------------------------------------------------------------------------
|
||||
// Constructor
|
||||
//--------------------------------------------------------------------------
|
||||
/**
|
||||
*
|
||||
*/
|
||||
PGbGLF::PGbGLF()
|
||||
{
|
||||
fWs = 0;
|
||||
fWs = gsl_integration_workspace_alloc(65536);
|
||||
if (fWs == 0) {
|
||||
cerr << "debug> PGbGLF::PGbGLF(): **ERROR** couldn't allocate fWs ..." << endl;
|
||||
}
|
||||
|
||||
fIntTab = 0;
|
||||
|
||||
fGslFunParam.wExt = 0.0;
|
||||
fGslFunParam.s0 = 0.0;
|
||||
fGslFunParam.s1 = 0.0;
|
||||
|
||||
fPrevParam[0] = 0.0;
|
||||
fPrevParam[1] = 0.0;
|
||||
fPrevParam[2] = 0.0;
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------
|
||||
// Destructor
|
||||
//--------------------------------------------------------------------------
|
||||
/**
|
||||
*
|
||||
*/
|
||||
PGbGLF::~PGbGLF()
|
||||
{
|
||||
if (fWs)
|
||||
gsl_integration_workspace_free(fWs);
|
||||
if (fIntTab)
|
||||
gsl_integration_qawo_table_free(fIntTab);
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------
|
||||
// operator()
|
||||
//--------------------------------------------------------------------------
|
||||
/**
|
||||
* <p>Method returning the function value at a given time for a given set of parameters.
|
||||
*
|
||||
* \param t time
|
||||
* \param par
|
||||
*/
|
||||
Double_t PGbGLF::operator()(Double_t t, const std::vector<Double_t> &par) const
|
||||
{
|
||||
// param: [0] Bext (G), [1] Delta0 (1/us), [2] DeltaGbG (1/us)
|
||||
assert(par.size()==3);
|
||||
|
||||
if (t<0.0)
|
||||
return 1.0;
|
||||
if (t>fTmax)
|
||||
t = fTmax;
|
||||
|
||||
Double_t dval = 0.0;
|
||||
Double_t wExt = TWO_PI * GAMMA_BAR_MUON * par[0];
|
||||
Double_t s0 = par[1];
|
||||
Double_t s1 = par[2];
|
||||
|
||||
Double_t wExt2 = wExt*wExt;
|
||||
Double_t s02 = s0*s0;
|
||||
Double_t s12 = s1*s1;
|
||||
Double_t t2 = t*t;
|
||||
|
||||
Double_t result = 0.0;
|
||||
Double_t err = 0.0;
|
||||
|
||||
// if integration table is not already setup, do it
|
||||
if (fIntTab == 0) {
|
||||
fIntTab = gsl_integration_qawo_table_alloc(wExt, fTmax, GSL_INTEG_SINE, 16);
|
||||
if (fIntTab == 0) {
|
||||
cerr << ">> PGbGLF::operator(): **ERROR** couldn't invoke qawo table." << endl;
|
||||
return 0.0;
|
||||
}
|
||||
fFun.function = &pz_GbG_2;
|
||||
fGslFunParam.wExt = wExt;
|
||||
fGslFunParam.s0 = s0;
|
||||
fGslFunParam.s1 = s1;
|
||||
fFun.params = &fGslFunParam;
|
||||
}
|
||||
// if parameter set has changed, adopted fFun parameters
|
||||
if (fPrevParam[1] != par[1]) {
|
||||
fPrevParam[0] = par[0];
|
||||
fPrevParam[1] = par[1];
|
||||
fPrevParam[2] = par[2];
|
||||
fGslFunParam.wExt = wExt;
|
||||
fGslFunParam.s0 = s0;
|
||||
fGslFunParam.s1 = s1;
|
||||
fFun.params = &fGslFunParam;
|
||||
}
|
||||
|
||||
// P_z^LF (GbG, 1st part)
|
||||
Double_t aa = 1.0+t2*s12;
|
||||
dval = 1.0 - 2.0*(s02+s12)/wExt2 + 2.0*(s02+s12*aa)/(wExt2*pow(aa,2.5))*cos(wExt*t)*exp(-0.5*s02*t2/aa);
|
||||
|
||||
// P_z^LF (GbG, 2nd part)
|
||||
gsl_integration_qawo_table_set_length(fIntTab, t);
|
||||
gsl_integration_qawo(&fFun, 0.0, 1.0e-9, 1.0e-9, 1024, fWs, fIntTab, &result, &err);
|
||||
dval += result;
|
||||
|
||||
return dval;
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------
|
||||
// pz_GbG_2 (private)
|
||||
//--------------------------------------------------------------------------
|
||||
/**
|
||||
* <p>Integrand of the non-analytic part of the GbG-LF polarization
|
||||
*
|
||||
* \param x
|
||||
* \param param
|
||||
*/
|
||||
|
||||
double pz_GbG_2(double x, void *param)
|
||||
{
|
||||
gslFunParam *p = (gslFunParam*) param;
|
||||
|
||||
double s02 = p->s0 * p->s0;
|
||||
double s12 = p->s1 * p->s1;
|
||||
double x2 = x*x;
|
||||
double aa = 1.0+x2*s12;
|
||||
|
||||
return 2.0*(s02*s02+3.0*s12*s12*aa*aa+6.0*s02*s12*aa)/(pow(p->wExt, 3.0)*pow(aa, 4.5))*exp(-0.5*s02*x2/aa);
|
||||
}
|
||||
76
src/external/libGbGLF/PGbGLF.h
vendored
Normal file
76
src/external/libGbGLF/PGbGLF.h
vendored
Normal file
@ -0,0 +1,76 @@
|
||||
/***************************************************************************
|
||||
|
||||
PGbGLF.h
|
||||
|
||||
Author: Andreas Suter
|
||||
e-mail: andreas.suter@psi.ch
|
||||
|
||||
***************************************************************************/
|
||||
|
||||
/***************************************************************************
|
||||
* Copyright (C) 2007-2016 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 _PGBGLF_H_
|
||||
#define _PGBGLF_H_
|
||||
|
||||
#include <vector>
|
||||
using namespace std;
|
||||
|
||||
#include <gsl/gsl_integration.h>
|
||||
|
||||
#include "PUserFcnBase.h"
|
||||
|
||||
typedef struct {
|
||||
double wExt;
|
||||
double s0;
|
||||
double s1;
|
||||
} gslFunParam;
|
||||
|
||||
|
||||
double pz_GbG_2(double x, void *param);
|
||||
|
||||
//--------------------------------------------------------------------------------------------
|
||||
/**
|
||||
* <p>Interface class for the user function.
|
||||
*/
|
||||
class PGbGLF : public PUserFcnBase
|
||||
{
|
||||
public:
|
||||
PGbGLF();
|
||||
virtual ~PGbGLF();
|
||||
|
||||
virtual Bool_t NeedGlobalPart() const { return false; }
|
||||
virtual void SetGlobalPart(vector<void *> &globalPart, UInt_t idx) {}
|
||||
virtual Bool_t GlobalPartIsValid() const { return true; }
|
||||
|
||||
virtual Double_t operator()(Double_t t, const std::vector<Double_t> ¶m) const;
|
||||
|
||||
private:
|
||||
static const Double_t fTmax = 12.0;
|
||||
mutable Double_t fPrevParam[3];
|
||||
mutable gslFunParam fGslFunParam;
|
||||
gsl_integration_workspace *fWs;
|
||||
mutable gsl_function fFun;
|
||||
mutable gsl_integration_qawo_table *fIntTab;
|
||||
|
||||
ClassDef(PGbGLF, 1)
|
||||
};
|
||||
|
||||
#endif // _PGBGLF_H_
|
||||
40
src/external/libGbGLF/PGbGLFLinkDef.h
vendored
Normal file
40
src/external/libGbGLF/PGbGLFLinkDef.h
vendored
Normal file
@ -0,0 +1,40 @@
|
||||
/***************************************************************************
|
||||
|
||||
PGbGLFLinkDef.h
|
||||
|
||||
Author: Andreas Suter
|
||||
|
||||
***************************************************************************/
|
||||
|
||||
/***************************************************************************
|
||||
* Copyright (C) 2016 by Andreas Suter *
|
||||
* *
|
||||
* *
|
||||
* 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. *
|
||||
***************************************************************************/
|
||||
|
||||
// root dictionary stuff --------------------------------------------------
|
||||
#ifdef __CINT__
|
||||
|
||||
#pragma link off all globals;
|
||||
#pragma link off all classes;
|
||||
#pragma link off all functions;
|
||||
|
||||
#pragma link C++ class PGbGLF+;
|
||||
|
||||
#endif //__CINT__
|
||||
// root dictionary stuff --------------------------------------------------
|
||||
|
||||
11
src/tests/GbG_LF/CMakeLists.txt
Normal file
11
src/tests/GbG_LF/CMakeLists.txt
Normal file
@ -0,0 +1,11 @@
|
||||
cmake_minimum_required(VERSION 2.8)
|
||||
|
||||
project(GbG_LF)
|
||||
|
||||
set(DependOnLibs gsl)
|
||||
set(DependOnLibs gslcblas)
|
||||
|
||||
add_executable(GbG_LF GbG_LF.cpp)
|
||||
|
||||
target_link_libraries(GbG_LF gsl gslcblas)
|
||||
|
||||
116
src/tests/GbG_LF/GbG_LF.cpp
Normal file
116
src/tests/GbG_LF/GbG_LF.cpp
Normal file
@ -0,0 +1,116 @@
|
||||
#include <iostream>
|
||||
#include <fstream>
|
||||
#include <vector>
|
||||
using namespace std;
|
||||
|
||||
#include <gsl/gsl_integration.h>
|
||||
|
||||
#define GAMMA_BAR_MUON 1.35538817e-2
|
||||
#define TWO_PI 6.28318530717958647692528676656
|
||||
|
||||
typedef struct {
|
||||
double wExt;
|
||||
double s0;
|
||||
double s1;
|
||||
} fun_params;
|
||||
|
||||
void syntax()
|
||||
{
|
||||
cout << endl;
|
||||
cout << "usage: GbG_LF Bext sigma0 sigmaGbG tmax fln" << endl;
|
||||
cout << " Bext in (G)" << endl;
|
||||
cout << " sigma0 in (1/us)" << endl;
|
||||
cout << " sigmaGbG in (1/us)" << endl;
|
||||
cout << " tmax in (us)" << endl;
|
||||
cout << " fln : output file name" << endl;
|
||||
cout << endl << endl;
|
||||
}
|
||||
|
||||
double pz_GbG_2(double x, void *param)
|
||||
{
|
||||
fun_params *par = (fun_params*) param;
|
||||
double s02 = par->s0 * par->s0;
|
||||
double s12 = par->s1 * par->s1;
|
||||
double x2 = x*x;
|
||||
double aa = 1.0+x2*s12;
|
||||
|
||||
return 2.0*(s02*s02+3.0*s12*s12*aa*aa+6.0*s02*s12*aa)/(pow(par->wExt, 3.0)*pow(aa, 4.5))*exp(-0.5*s02*x2/aa);
|
||||
}
|
||||
|
||||
int main(int argc, char *argv[])
|
||||
{
|
||||
if (argc != 6) {
|
||||
syntax();
|
||||
return -1;
|
||||
}
|
||||
|
||||
double Bext = atof(argv[1]);
|
||||
double sigma0 = atof(argv[2]);
|
||||
double sigma1 = atof(argv[3]);
|
||||
double tmax = atof(argv[4]);
|
||||
|
||||
// calculate needed time step size
|
||||
double dt = 1.0 / (GAMMA_BAR_MUON*Bext) / 10.0;
|
||||
if (10.0/dt < 500.0)
|
||||
dt = 0.02;
|
||||
cout << "dt = " << dt << " (us)" << endl;
|
||||
|
||||
double wExt = TWO_PI * GAMMA_BAR_MUON * Bext;
|
||||
double s0 = sigma0;
|
||||
double s1 = sigma1;
|
||||
|
||||
double wExt2 = wExt*wExt;
|
||||
double s02 = s0*s0;
|
||||
double s12 = s1*s1;
|
||||
|
||||
gsl_integration_workspace *w = gsl_integration_workspace_alloc(65536);
|
||||
gsl_function fun;
|
||||
fun.function = &pz_GbG_2;
|
||||
fun_params param = { wExt, s0, s1 };
|
||||
fun.params = ¶m;
|
||||
|
||||
gsl_integration_qawo_table *tab = gsl_integration_qawo_table_alloc(wExt, tmax, GSL_INTEG_SINE, 16);
|
||||
|
||||
double result = 0.0, err = 0.0;
|
||||
double t = 0.0, t2 = 0.0, dval=0.0, aa=0.0;
|
||||
vector<double> time;
|
||||
vector<double> pol, pol1, pol2;
|
||||
do {
|
||||
time.push_back(t);
|
||||
t2 = t*t;
|
||||
|
||||
dval = 0.0;
|
||||
|
||||
// P_z^LF (GbG, 1st part)
|
||||
aa = 1.0+t2*s12;
|
||||
dval = 1.0 - 2.0*(s02+s12)/wExt2 + 2.0*(s02+s12*aa)/(wExt2*pow(aa,2.5))*cos(wExt*t)*exp(-0.5*s02*t2/aa);
|
||||
pol1.push_back(dval);
|
||||
|
||||
// P_z^LF (GbG, 2nd part)
|
||||
gsl_integration_qawo_table_set_length(tab, t);
|
||||
gsl_integration_qawo(&fun, 0.0, 1.0e-9, 1.0e-9, 1024, w, tab, &result, &err);
|
||||
pol2.push_back(result);
|
||||
dval += result;
|
||||
|
||||
pol.push_back(dval);
|
||||
|
||||
t += dt;
|
||||
} while (t <= tmax);
|
||||
|
||||
gsl_integration_workspace_free(w);
|
||||
gsl_integration_qawo_table_free(tab);
|
||||
|
||||
ofstream fout(argv[5], ofstream::out);
|
||||
fout << "% Bext (G) : " << Bext << endl;
|
||||
fout << "% sigma0 (1/us) : " << sigma0 << endl;
|
||||
fout << "% sigma1 (1/us) : " << sigma1 << endl;
|
||||
fout << "%----" << endl;
|
||||
fout << "% time (us), pol1, pol2, pol" << endl;
|
||||
for (unsigned int i=0; i<time.size(); i++) {
|
||||
fout << time[i] << ", " << pol1[i] << ", " << pol2[i] << ", " << pol[i] << endl;
|
||||
}
|
||||
|
||||
fout.close();
|
||||
|
||||
return 0;
|
||||
}
|
||||
78
src/tests/MuTransition/7700.msr
Normal file
78
src/tests/MuTransition/7700.msr
Normal file
@ -0,0 +1,78 @@
|
||||
07700- complexMuPol, A0 100MHz, Mu-frac 1.00, Mu12 40.43 MHz(0.49), Mu23 256.25 MHz(0.01), ionRate 0.000 MHz, capRate 0.000 MHz, SF rate 0.010 MHz, 106.5 G
|
||||
###############################################################
|
||||
FITPARAMETER
|
||||
# Nr. Name Value Step Pos_Error Boundaries
|
||||
1 alpha 0.99961 -0.00091 0.00091 0 none
|
||||
2 asy 0 0 none 0 0.33
|
||||
3 phase 0 0 none
|
||||
4 field 106.5 0 none 0 none
|
||||
5 rate 0 0 none 0 150
|
||||
6 asyMu12 0.135 0 none
|
||||
7 phaseMu12 0.68 -0.42 0.42
|
||||
8 freqMu12 40.43248 -0.00045 0.00045
|
||||
9 rateMu12 0.0104 -0.0019 0.0019 0 100
|
||||
10 asyMu34 0.135 0 none
|
||||
11 phaseMu34 -0.01 -0.43 0.43
|
||||
12 freqMu34 59.56737 -0.00046 0.00046
|
||||
13 rateMu34 0.0164 -0.0019 0.0019 0 100
|
||||
|
||||
###############################################################
|
||||
THEORY
|
||||
asymmetry 2
|
||||
TFieldCos 3 fun1 (phase frequency)
|
||||
simplExpo 5 (rate)
|
||||
+
|
||||
asymmetry 6
|
||||
TFieldCos 7 8 (phase frequency)
|
||||
simplExpo 9 (rate)
|
||||
+
|
||||
asymmetry 10
|
||||
TFieldCos 11 12 (phase frequency)
|
||||
simplExpo 13 (rate)
|
||||
|
||||
###############################################################
|
||||
FUNCTIONS
|
||||
fun1 = par4 * gamma_mu
|
||||
fun2 = par8 * gamma_mu
|
||||
|
||||
###############################################################
|
||||
GLOBAL
|
||||
fittype 2 (asymmetry fit)
|
||||
data 1 18000 1 18000
|
||||
fit 0.005 8
|
||||
packing 1
|
||||
|
||||
###############################################################
|
||||
RUN 07700 MUE4 PSI MUSR-ROOT (name beamline institute data-file-format)
|
||||
#ADDRUN 07701 MUE4 PSI MUSR-ROOT (name beamline institute data-file-format)
|
||||
alpha 1
|
||||
map 0 0 0 0 0 0 0 0 0 0
|
||||
forward 1
|
||||
backward 2
|
||||
backgr.fix 0.000000 0.000000
|
||||
|
||||
###############################################################
|
||||
COMMANDS
|
||||
MINIMIZE
|
||||
MINOS
|
||||
SAVE
|
||||
|
||||
###############################################################
|
||||
FOURIER
|
||||
units MHz # units either 'Gauss', 'Tesla', 'MHz', or 'Mc/s'
|
||||
fourier_power 12
|
||||
apodization STRONG # NONE, WEAK, MEDIUM, STRONG
|
||||
plot POWER # REAL, IMAG, REAL_AND_IMAG, POWER, PHASE
|
||||
phase 8.000000
|
||||
#range_for_phase_correction 50.0 70.0
|
||||
#range 0 200
|
||||
|
||||
###############################################################
|
||||
PLOT 2 (asymmetry plot)
|
||||
runs 1
|
||||
range 0 8
|
||||
view_packing 2
|
||||
|
||||
###############################################################
|
||||
STATISTIC --- 2016-03-03 17:25:40
|
||||
chisq = 8091.8, NDF = 7988, chisq/NDF = 1.012996
|
||||
@ -43,38 +43,38 @@ void runMuSimulation()
|
||||
|
||||
char titleStr[256];
|
||||
TFolder *histosFolder;
|
||||
TFolder *decayAnaModule;
|
||||
TFolder *decayAnaModule, *scAnaModule;
|
||||
TFolder *gRunHeader;
|
||||
TString runTitle;
|
||||
TString histogramFileName;
|
||||
TObjArray Slist(0);
|
||||
TMusrRunPhysicalQuantity prop;
|
||||
|
||||
//prepare to run simulation; here: isotropic Mu in Germanium
|
||||
UInt_t runNo = 9900;
|
||||
//prepare to run simulation; here: isotropic Mu with A0 = 100.0 MHz
|
||||
UInt_t runNo = 7701;
|
||||
Double_t T = 300.; //temperature
|
||||
Double_t spinFlipRate = 0.001; //if spinFlipRate > 0.001 only spin-flip processes will be simulated
|
||||
Double_t capRate = 1.0;//*sqrt(T/200.); //assume that capture rate varies as sqrt(T), capRate = sigma*v*p , v ~ sqrt(T)
|
||||
Double_t EA = 100; //activation energy (meV)
|
||||
Double_t spinFlipRate = 0.01; //if spinFlipRate > 0.001 only spin-flip processes will be simulated
|
||||
Double_t capRate = 0.0001;//*sqrt(T/200.); //assume that capture rate varies as sqrt(T), capRate = sigma*v*p , v ~ sqrt(T)
|
||||
Double_t ionRate; //assume Arrhenius behaviour ionRate = preFac*exp(-EA/kT)
|
||||
Double_t EA = 100.; //activation energy (meV)
|
||||
ionRate = 2.9e7 * exp(-EA/(0.08625*T)); // Ge: 2.9*10^7MHz "attempt" frequency; 1K = 0.08625 meV
|
||||
Double_t B = 100.; //field in G
|
||||
ionRate = 0.0001; //2.9e7 * exp(-EA/(0.08625*T)); // Ge: 2.9*10^7MHz "attempt" frequency; 1K = 0.08625 meV
|
||||
Double_t B = 106.5; //field in G
|
||||
Double_t Bvar = 0.; //field variance
|
||||
Double_t Freq12 = 134.858; //Mu freq of the 12 transition
|
||||
Double_t Freq34 = 4328.142; //Mu freq of the 34 transition
|
||||
Double_t Freq23 = 143.713; //Mu freq of the 23 transition
|
||||
Double_t Freq14 = 4606.713; //Mu freq of the 14 transition
|
||||
Double_t Freq12 = 40.433; //Mu freq of the 12 transition
|
||||
Double_t Freq34 = 59.567; //Mu freq of the 34 transition
|
||||
Double_t Freq23 = 256.245; //Mu freq of the 23 transition
|
||||
Double_t Freq14 = 356.245; //Mu freq of the 14 transition
|
||||
Double_t MuFrac = 1.0; //total Mu fraction
|
||||
Double_t MuFrac12 = 2*0.266; //Mu in states 12 and 34
|
||||
Double_t MuFrac23 = 2*0.234; //Mu in states 23 and 14
|
||||
Int_t Nmuons = 1e6; //number of muons
|
||||
Double_t MuFrac12 = 2*0.487; //Mu in states 12 and 34
|
||||
Double_t MuFrac23 = 2*0.013; //Mu in states 23 and 14
|
||||
Int_t Nmuons = 5e6; //number of muons
|
||||
Double_t Asym = 0.27; //muon decay asymmetry
|
||||
|
||||
histogramFileName = TString("0");
|
||||
histogramFileName += runNo;
|
||||
histogramFileName += TString(".root");
|
||||
|
||||
sprintf(titleStr,"- Mu-frac %3.2f, Mu12 %6.2f MHz(%3.2f), Mu23 %6.2f MHz(%3.2f), ionRate %8.2f MHz, capRate %6.2f MHz, SF rate %6.2f MHz, %5.0f G", MuFrac, Freq12, MuFrac12/2, Freq23, MuFrac23/2, ionRate, capRate, spinFlipRate, B);
|
||||
sprintf(titleStr,"- complexMuPol, A0 100MHz, Mu-frac %3.2f, Mu12 %6.2f MHz(%3.2f), Mu23 %6.2f MHz(%3.2f), ionRate %8.3f MHz, capRate %6.3f MHz, SF rate %6.3f MHz, %5.1f G", MuFrac, Freq12, MuFrac12/2, Freq23, MuFrac23/2, ionRate, capRate, spinFlipRate, B);
|
||||
runTitle = TString("0");
|
||||
runTitle += runNo;
|
||||
runTitle += TString(titleStr);
|
||||
@ -103,76 +103,92 @@ void runMuSimulation()
|
||||
gRunHeader = gROOT->GetRootFolder()->AddFolder("RunHeader", "MuTransition Simulation Header Info");
|
||||
gROOT->GetListOfBrowsables()->Add(gRunHeader, "RunHeader");
|
||||
// header = new TLemRunHeader();
|
||||
header = new TMusrRunHeader(true);
|
||||
header->FillFolder(gRunHeader);
|
||||
gRunHeader->Add(&Slist);
|
||||
Slist.SetName("RunSummary");
|
||||
TMusrRunHeader *header = new TMusrRunHeader(true);
|
||||
|
||||
header->Set("RunInfo/Generic Validator URL", "http://lmu.web.psi.ch/facilities/software/MusrRoot/validation/MusrRoot.xsd");
|
||||
header->Set("RunInfo/Specific Validator URL", "http://lmu.web.psi.ch/facilities/software/MusrRoot/validation/MusrRootLEM.xsd");
|
||||
header->Set("RunInfo/Generator", "runMuSimulation");
|
||||
|
||||
|
||||
header->Set("RunInfo/File Name", histogramFileName.Data());
|
||||
header->Set("RunInfo/Run Title", runTitle.Data());
|
||||
header->Set("RunInfo/Run Number", (Int_t) runNo);
|
||||
header->Set("RunInfo/Run Start Time", "0");
|
||||
header->Set("RunInfo/Run Stop Time", "1");
|
||||
prop.Set("Run Duration", 1.0, "sec");
|
||||
header->Set("RunInfo/Run Start Time", "2016-03-01 06:20:00");
|
||||
header->Set("RunInfo/Run Stop Time", "2016-03-01 06:20:11");
|
||||
|
||||
prop.Set("Run Duration", 11.0, "sec");
|
||||
header->Set("RunInfo/Run Duration", prop);
|
||||
|
||||
header->Set("RunInfo/Laboratory", "PSI");
|
||||
|
||||
header->Set("RunInfo/Instrument", "MC-Simulation");
|
||||
prop.Set("Muon Beam Momentum", 0.0, "MeV/c");
|
||||
|
||||
header->Set("RunInfo/Muon Beam Momentum", prop);
|
||||
header->Set("RunInfo/Muon Species", "positive muon and muonium");
|
||||
header->Set("RunInfo/Muon Source", "Simulation");
|
||||
header->Set("RunInfo/Setup", "Monte-Carlo setup");
|
||||
header->Set("RunInfo/Comment", "Testing effect of charge-exchange or Mu0 spin flip processes on uSR signal");
|
||||
header->Set("RunInfo/Sample Name", "Monte-Carlo");
|
||||
header->Set("RunInfo/Sample Temperature", 300);
|
||||
|
||||
prop.Set("Sample Temperature", MRH_UNDEFINED, T, 0.01, "K");
|
||||
header->Set("RunInfo/Sample Temperature", prop);
|
||||
|
||||
prop.Set("Sample Magnetic Field", MRH_UNDEFINED, B, Bvar, "G");
|
||||
header->Set("RunInfo/Sample Magnetic Field", prop);
|
||||
|
||||
header->Set("RunInfo/No of Histos", 2);
|
||||
|
||||
prop.Set("Time Resolution", 1.0, "ns", "Simulation");
|
||||
header->Set("RunInfo/Time Resolution", prop);
|
||||
|
||||
prop.Set("Implantation Energy", 0, "keV");
|
||||
header->Set("RunInfo/Implantation Energy", prop);
|
||||
|
||||
prop.Set("Muon Spin Angle", 0, "degree along x");
|
||||
header->Set("RunInfo/Muon Spin Angle", prop);
|
||||
|
||||
header->Set("DetectorInfo/Detector001/Name", "hDecay001");
|
||||
header->Set("DetectorInfo/Detector001/Name", "e+ forward");
|
||||
header->Set("DetectorInfo/Detector001/Histo Number", 1);
|
||||
header->Set("DetectorInfo/Detector001/Histo Length", 12001);
|
||||
header->Set("DetectorInfo/Detector001/Time Zero Bin", 0);
|
||||
header->Set("DetectorInfo/Detector001/Histo Length", 18001);
|
||||
header->Set("DetectorInfo/Detector001/Time Zero Bin", 0.001); //doesn't like 0.0 as time zero
|
||||
header->Set("DetectorInfo/Detector001/First Good Bin", 1);
|
||||
header->Set("DetectorInfo/Detector001/Last Good Bin", 12001);
|
||||
|
||||
header->Set("DetectorInfo/Detector002/Name", "hDecay002");
|
||||
header->Set("DetectorInfo/Detector002/Histo Number", 1);
|
||||
header->Set("DetectorInfo/Detector002/Histo Length", 12001);
|
||||
header->Set("DetectorInfo/Detector002/Time Zero Bin", 0);
|
||||
header->Set("DetectorInfo/Detector001/Last Good Bin", 18000);
|
||||
|
||||
header->Set("DetectorInfo/Detector002/Name", "e+ backward");
|
||||
header->Set("DetectorInfo/Detector002/Histo Number", 2);
|
||||
header->Set("DetectorInfo/Detector002/Histo Length", 18001);
|
||||
header->Set("DetectorInfo/Detector002/Time Zero Bin", 0.001);
|
||||
header->Set("DetectorInfo/Detector002/First Good Bin", 1);
|
||||
header->Set("DetectorInfo/Detector002/Last Good Bin", 12001);
|
||||
|
||||
header->Set("DetectorInfo/Detector002/Last Good Bin", 18000);
|
||||
|
||||
// simulation parameters
|
||||
header->Set("Simulation/Mu Precession frequency 12", Freq12);
|
||||
header->Set("Simulation/Mu Precession frequency 34", Freq34);
|
||||
header->Set("Simulation/Mu Precession frequency 23", Freq23);
|
||||
header->Set("Simulation/Mu Precession frequency 14", Freq14);
|
||||
header->Set("Simulation/Mu Fraction", MuFrac);
|
||||
header->Set("Simulation/Mu Fraction 12", MuFrac12);
|
||||
header->Set("Simulation/Mu Fraction 23", MuFrac23);
|
||||
header->Set("Simulation/Mu+ Capture Rate", capRate);
|
||||
header->Set("Simulation/Mu0 Precession frequency 12", Freq12);
|
||||
header->Set("Simulation/Mu0 Precession frequency 34", Freq34);
|
||||
header->Set("Simulation/Mu0 Precession frequency 23", Freq23);
|
||||
header->Set("Simulation/Mu0 Precession frequency 14", Freq14);
|
||||
header->Set("Simulation/Mu0 Fraction", MuFrac);
|
||||
header->Set("Simulation/Mu0 Fraction 12", MuFrac12);
|
||||
header->Set("Simulation/Mu0 Fraction 23", MuFrac23);
|
||||
header->Set("Simulation/muon Capture Rate", capRate);
|
||||
header->Set("Simulation/Mu0 Ionization Rate", ionRate);
|
||||
header->Set("Simulation/Mu0 Spin Flip Rate", spinFlipRate);
|
||||
header->Set("Simulation/Number of Muons", Nmuons);
|
||||
header->Set("Simulation/Decay Asymmetry", Asym);
|
||||
|
||||
header->Set("SampleEnvironmentInfo/Cryo", "no cryostat");
|
||||
header->Set("MagneticFieldEnvironmentInfo/Magnet Name", "Field along z");
|
||||
header->Set("BeamlineInfo/Name", "Monte-Carlo setup");
|
||||
|
||||
histosFolder = gROOT->GetRootFolder()->AddFolder("histos", "Histograms");
|
||||
gROOT->GetListOfBrowsables()->Add(histosFolder, "histos");
|
||||
decayAnaModule = histosFolder->AddFolder("DecayAnaModule", "muSR decay histograms");
|
||||
decayAnaModule = histosFolder->AddFolder("DecayAnaModule", "muSR decay histograms");
|
||||
scAnaModule = histosFolder->AddFolder("SCAnaModule", "SlowControl histograms");
|
||||
|
||||
TH1F *histo[NDECAYHISTS];
|
||||
char str[128];
|
||||
for (UInt_t i=0; i<NDECAYHISTS; i++) {
|
||||
sprintf(str, "hDecay00%d", (Int_t)i+1);
|
||||
histo[i] = new TH1F(str, str, 12001, -0.5, 12000.5);
|
||||
histo[i] = new TH1F(str, str, 18001, -0.5, 18000.5);
|
||||
}
|
||||
|
||||
for (i=0; i<NDECAYHISTS; i++)
|
||||
@ -191,9 +207,12 @@ void runMuSimulation()
|
||||
}
|
||||
|
||||
fout->cd();
|
||||
|
||||
header->FillFolder(gRunHeader);
|
||||
gRunHeader->Write();
|
||||
gRunHeader->Add(&Slist);
|
||||
Slist.SetName("RunSummary");
|
||||
histosFolder->Write();
|
||||
gRunHeader->Write();
|
||||
fout->Close();
|
||||
cout << "Histograms written to " << histogramFileName.Data() << endl;
|
||||
delete fout;
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user