From cd041649438b04ef0dd89f46736137653515b12b Mon Sep 17 00:00:00 2001 From: Thomas Prokscha Date: Thu, 3 Mar 2016 17:28:50 +0100 Subject: [PATCH 1/3] Fixed problems with run header --- src/tests/MuTransition/7700.msr | 78 ++++++++++++++++ src/tests/MuTransition/runMuSimulation.C | 113 +++++++++++++---------- 2 files changed, 144 insertions(+), 47 deletions(-) create mode 100644 src/tests/MuTransition/7700.msr diff --git a/src/tests/MuTransition/7700.msr b/src/tests/MuTransition/7700.msr new file mode 100644 index 00000000..ee61c73f --- /dev/null +++ b/src/tests/MuTransition/7700.msr @@ -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 diff --git a/src/tests/MuTransition/runMuSimulation.C b/src/tests/MuTransition/runMuSimulation.C index 96b3f5d4..07d66475 100644 --- a/src/tests/MuTransition/runMuSimulation.C +++ b/src/tests/MuTransition/runMuSimulation.C @@ -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; icd(); + 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; From 47580b66d877e6c8104b007b911f7e3be92a72a4 Mon Sep 17 00:00:00 2001 From: Suter Andreas Date: Fri, 4 Mar 2016 09:43:05 +0100 Subject: [PATCH 2/3] added a GbG LF calculator (standalone for now) --- src/tests/GbG_LF/CMakeLists.txt | 11 +++ src/tests/GbG_LF/GbG_LF.cpp | 116 ++++++++++++++++++++++++++++++++ 2 files changed, 127 insertions(+) create mode 100644 src/tests/GbG_LF/CMakeLists.txt create mode 100644 src/tests/GbG_LF/GbG_LF.cpp diff --git a/src/tests/GbG_LF/CMakeLists.txt b/src/tests/GbG_LF/CMakeLists.txt new file mode 100644 index 00000000..af870c2e --- /dev/null +++ b/src/tests/GbG_LF/CMakeLists.txt @@ -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 gslcblas gsl) + diff --git a/src/tests/GbG_LF/GbG_LF.cpp b/src/tests/GbG_LF/GbG_LF.cpp new file mode 100644 index 00000000..7b25806f --- /dev/null +++ b/src/tests/GbG_LF/GbG_LF.cpp @@ -0,0 +1,116 @@ +#include +#include +#include +using namespace std; + +#include + +#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 time; + vector 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 Date: Mon, 14 Mar 2016 11:30:06 +0100 Subject: [PATCH 3/3] added user function for GbG LF --- configure.ac | 1 + src/classes/PTheory.cpp | 1 + src/external/Makefile.am | 3 +- src/external/libGbGLF/Makefile.am | 56 +++++++++ src/external/libGbGLF/PGbGLF.cpp | 173 ++++++++++++++++++++++++++ src/external/libGbGLF/PGbGLF.h | 76 +++++++++++ src/external/libGbGLF/PGbGLFLinkDef.h | 40 ++++++ src/tests/GbG_LF/CMakeLists.txt | 2 +- 8 files changed, 350 insertions(+), 2 deletions(-) create mode 100644 src/external/libGbGLF/Makefile.am create mode 100644 src/external/libGbGLF/PGbGLF.cpp create mode 100644 src/external/libGbGLF/PGbGLF.h create mode 100644 src/external/libGbGLF/PGbGLFLinkDef.h diff --git a/configure.ac b/configure.ac index 4931abfc..86f073d6 100644 --- a/configure.ac +++ b/configure.ac @@ -1239,6 +1239,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 \ diff --git a/src/classes/PTheory.cpp b/src/classes/PTheory.cpp index efa1c7d8..30c0003c 100644 --- a/src/classes/PTheory.cpp +++ b/src/classes/PTheory.cpp @@ -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()); } diff --git a/src/external/Makefile.am b/src/external/Makefile.am index 755d9210..c2d05957 100644 --- a/src/external/Makefile.am +++ b/src/external/Makefile.am @@ -3,7 +3,8 @@ if BUILD_ASLIBS ASDIRS = Nonlocal \ MagProximity \ libSpinValve \ - libPhotoMeissner + libPhotoMeissner \ + libGbGLF endif if BUILD_CUBALIB diff --git a/src/external/libGbGLF/Makefile.am b/src/external/libGbGLF/Makefile.am new file mode 100644 index 00000000..2c77bbd1 --- /dev/null +++ b/src/external/libGbGLF/Makefile.am @@ -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 diff --git a/src/external/libGbGLF/PGbGLF.cpp b/src/external/libGbGLF/PGbGLF.cpp new file mode 100644 index 00000000..edfafba5 --- /dev/null +++ b/src/external/libGbGLF/PGbGLF.cpp @@ -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 +#include +#include +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() +//-------------------------------------------------------------------------- +/** + *

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 &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) +//-------------------------------------------------------------------------- +/** + *

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); +} diff --git a/src/external/libGbGLF/PGbGLF.h b/src/external/libGbGLF/PGbGLF.h new file mode 100644 index 00000000..5a07f748 --- /dev/null +++ b/src/external/libGbGLF/PGbGLF.h @@ -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 +using namespace std; + +#include + +#include "PUserFcnBase.h" + +typedef struct { + double wExt; + double s0; + double s1; +} gslFunParam; + + +double pz_GbG_2(double x, void *param); + +//-------------------------------------------------------------------------------------------- +/** + *

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 &globalPart, UInt_t idx) {} + virtual Bool_t GlobalPartIsValid() const { return true; } + + virtual Double_t operator()(Double_t t, const std::vector ¶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_ diff --git a/src/external/libGbGLF/PGbGLFLinkDef.h b/src/external/libGbGLF/PGbGLFLinkDef.h new file mode 100644 index 00000000..30ebc943 --- /dev/null +++ b/src/external/libGbGLF/PGbGLFLinkDef.h @@ -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 -------------------------------------------------- + diff --git a/src/tests/GbG_LF/CMakeLists.txt b/src/tests/GbG_LF/CMakeLists.txt index af870c2e..56dd3d7e 100644 --- a/src/tests/GbG_LF/CMakeLists.txt +++ b/src/tests/GbG_LF/CMakeLists.txt @@ -7,5 +7,5 @@ set(DependOnLibs gslcblas) add_executable(GbG_LF GbG_LF.cpp) -target_link_libraries(GbG_LF gslcblas gsl) +target_link_libraries(GbG_LF gsl gslcblas)