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

Merged master into beta-NMR

Browse Source
This commit is contained in:
Zaher Salman
2018-12-18 10:01:28 +01:00
parent 8207aea183 91e2eca915
commit d1cf374354
49 changed files with 1603 additions and 266 deletions
Download Patch File Download Diff File Expand all files Collapse all files

5
src/CMakeLists.txt
View File

@@ -49,6 +49,7 @@ add_executable(any2many git-revision.h any2many.cpp)
target_compile_options(any2many BEFORE PRIVATE "-DHAVE_CONFIG_H")
target_include_directories(any2many
BEFORE PRIVATE
$<BUILD_INTERFACE:${Boost_INCLUDE_DIR}>
$<BUILD_INTERFACE:${FFTW3_INCLUDE_DIR}>
$<BUILD_INTERFACE:${CMAKE_BINARY_DIR}>
$<BUILD_INTERFACE:${CMAKE_BINARY_DIR}/src>
@@ -93,6 +94,7 @@ add_executable(musrfit git-revision.h musrfit.cpp)
target_compile_options(musrfit BEFORE PRIVATE "-DHAVE_CONFIG_H")
target_include_directories(musrfit
BEFORE PRIVATE
$<BUILD_INTERFACE:${Boost_INCLUDE_DIR}>
$<BUILD_INTERFACE:${FFTW3_INCLUDE_DIR}>
$<BUILD_INTERFACE:${CMAKE_BINARY_DIR}>
$<BUILD_INTERFACE:${CMAKE_BINARY_DIR}/src>
@@ -104,6 +106,7 @@ add_executable(musrFT git-revision.h musrFT.cpp)
target_compile_options(musrFT BEFORE PRIVATE "-DHAVE_CONFIG_H")
target_include_directories(musrFT
BEFORE PRIVATE
$<BUILD_INTERFACE:${Boost_INCLUDE_DIR}>
$<BUILD_INTERFACE:${FFTW3_INCLUDE_DIR}>
$<BUILD_INTERFACE:${CMAKE_BINARY_DIR}>
$<BUILD_INTERFACE:${CMAKE_BINARY_DIR}/src>
@@ -127,6 +130,7 @@ add_executable(musrt0 git-revision.h musrt0.cpp)
target_compile_options(musrt0 BEFORE PRIVATE "-DHAVE_CONFIG_H")
target_include_directories(musrt0
BEFORE PRIVATE
$<BUILD_INTERFACE:${Boost_INCLUDE_DIR}>
$<BUILD_INTERFACE:${FFTW3_INCLUDE_DIR}>
$<BUILD_INTERFACE:${CMAKE_BINARY_DIR}>
$<BUILD_INTERFACE:${CMAKE_BINARY_DIR}/src>
@@ -138,6 +142,7 @@ add_executable(musrview git-revision.h musrview.cpp)
target_compile_options(musrview BEFORE PRIVATE "-DHAVE_CONFIG_H")
target_include_directories(musrview
BEFORE PRIVATE
$<BUILD_INTERFACE:${Boost_INCLUDE_DIR}>
$<BUILD_INTERFACE:${FFTW3_INCLUDE_DIR}>
$<BUILD_INTERFACE:${CMAKE_BINARY_DIR}>
$<BUILD_INTERFACE:${CMAKE_BINARY_DIR}/src>

5
src/classes/CMakeLists.txt
View File

@@ -11,13 +11,13 @@ root_generate_dictionary(
)
root_generate_dictionary(
PMusrCanvasDict
-I${FFTW3_INCLUDE_DIR} ${MUSRFIT_INC}/PMusrCanvas.h
-I${Boost_INCLUDE_DIR} -I${FFTW3_INCLUDE_DIR} ${MUSRFIT_INC}/PMusrCanvas.h
LINKDEF ${MUSRFIT_INC}/PMusrCanvasLinkDef.h
MODULE PMusrCanvas
)
root_generate_dictionary(
PMusrT0Dict
-I${FFTW3_INCLUDE_DIR} ${MUSRFIT_INC}/PMusrT0.h
-I${Boost_INCLUDE_DIR} -I${FFTW3_INCLUDE_DIR} ${MUSRFIT_INC}/PMusrT0.h
LINKDEF ${MUSRFIT_INC}/PMusrT0LinkDef.h
MODULE PMusrT0
)
@@ -81,6 +81,7 @@ add_library(PMusr SHARED
#--- make sure that the include directory is found ----------------------------
target_include_directories(
PMusr BEFORE PRIVATE
$<BUILD_INTERFACE:${Boost_INCLUDE_DIR}>
$<BUILD_INTERFACE:${FFTW3_INCLUDE_DIR}>
$<BUILD_INTERFACE:${NEXUS_INCLUDE_DIR}>
$<BUILD_INTERFACE:${MUSRFIT_INC}>

10
src/classes/PMsr2Data.cpp
View File

@@ -1025,12 +1025,11 @@ bool PMsr2Data::PrepareGlobalInputFile(unsigned int tempRun, const string &msrOu
}
// FOURIER block - in case a parameter is used for the phase
tempPar = fMsrHandler->GetMsrFourierList()->fPhaseParamNo;
if (tempPar > 0) {
if (fMsrHandler->GetMsrFourierList()->fPhaseParamNo.size() > 0) {
// go through the whole parameter list ...
for (unsigned int k(0); k < msrParamList->size(); ++k) {
if (tempPar == msrParamList->at(k).fNo) {
fMsrHandler->GetMsrFourierList()->fPhaseParamNo = k + 1;
fMsrHandler->GetMsrFourierList()->fPhaseParamNo.push_back(k + 1);
break;
}
}
@@ -2199,7 +2198,8 @@ int PMsr2Data::WriteOutput(const string &outfile, const vector<unsigned int>& pa
WriteValue(outFile, (*msrParamList)[i].fValue, (*msrParamList)[i].fPosError, outFile.width(), db);
outFile << ", ";
} else {
outFile << (*msrParamList)[i].fValue << ", ";
WriteValue(outFile, (*msrParamList)[i].fValue, fabs((*msrParamList)[i].fStep), outFile.width(), db);
outFile << ", ";
}
if ((*msrParamList)[i].fPosErrorPresent) {
WriteValue(outFile, (*msrParamList)[i].fPosError, (*msrParamList)[i].fPosError, outFile.width(), db);
@@ -2417,7 +2417,7 @@ int PMsr2Data::WriteOutput(const string &outfile, const vector<unsigned int>& pa
if ((*msrParamList)[i].fPosErrorPresent)
WriteValue(outFile, (*msrParamList)[i].fValue, (*msrParamList)[i].fPosError, maxlength, db);
else
WriteValue(outFile, (*msrParamList)[i].fValue, maxlength);
WriteValue(outFile, (*msrParamList)[i].fValue, fabs((*msrParamList)[i].fStep), maxlength, db);
if ((*msrParamList)[i].fPosErrorPresent)
WriteValue(outFile, (*msrParamList)[i].fPosError, (*msrParamList)[i].fPosError, maxlength, db);

493
src/classes/PMsrHandler.cpp
View File

@@ -8,7 +8,7 @@
***************************************************************************/
/***************************************************************************
* Copyright (C) 2007-2016 by Andreas Suter *
* Copyright (C) 2007-2018 by Andreas Suter *
* andreas.suter@psi.ch *
* *
* This program is free software; you can redistribute it and/or modify *
@@ -324,6 +324,15 @@ Int_t PMsrHandler::ReadMsrFile()
CheckLegacyLifetimecorrection(); // check if lifetimecorrection is found in RUN blocks, if yes transfer it to PLOT blocks
}
// check if the given phases in the Fourier block are in agreement with the Plot block settings
if ((fFourier.fPhase.size() > 1) && (fPlots.size() > 0)) {
if (fFourier.fPhase.size() != fPlots[0].fRuns.size()) {
cerr << endl << ">> PMsrHandler::ReadMsrFile: **ERROR** if more than one phase is given in the Fourier block,";
cerr << endl << ">> it needs to correspond to the number of runs in the Plot block!" << endl;
result = PMUSR_MSR_SYNTAX_ERROR;
}
}
// clean up
fit_parameter.clear();
theory.clear();
@@ -1122,12 +1131,15 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t messages)
fout << " # REAL, IMAG, REAL_AND_IMAG, POWER, PHASE, PHASE_OPT_REAL";
fout << endl;
} else if (sstr.BeginsWith("phase")) {
if (fFourier.fPhaseParamNo > 0) {
fout << "phase par" << fFourier.fPhaseParamNo << endl;
} else {
if (fFourier.fPhase != -999.0) {
fout << "phase " << fFourier.fPhase << endl;
if (fFourier.fPhaseParamNo.size() > 0) {
TString phaseParamStr = BeautifyFourierPhaseParameterString();
fout << "phase " << phaseParamStr << endl;
} else if (fFourier.fPhase.size() > 0) {
fout << "phase ";
for (UInt_t i=0; i<fFourier.fPhase.size()-1; i++) {
fout << fFourier.fPhase[i] << ", ";
}
fout << fFourier.fPhase[fFourier.fPhase.size()-1] << endl;
}
} else if (sstr.BeginsWith("range_for_phase_correction")) {
fout << "range_for_phase_correction " << fFourier.fRangeForPhaseCorrection[0] << " " << fFourier.fRangeForPhaseCorrection[1] << endl;
@@ -2213,11 +2225,16 @@ Int_t PMsrHandler::WriteMsrFile(const Char_t *filename, map<UInt_t, TString> *co
fout << endl;
}
// phase
if (fFourier.fPhaseParamNo > 0) {
fout << "phase par" << fFourier.fPhaseParamNo << endl;
} else if (fFourier.fPhase != -999.0) {
fout << "phase " << fFourier.fPhase << endl;
// phase
if (fFourier.fPhaseParamNo.size() > 0) {
TString phaseParamStr = BeautifyFourierPhaseParameterString();
fout << "phase " << phaseParamStr << endl;
} else if (fFourier.fPhase.size() > 0) {
fout << "phase ";
for (UInt_t i=0; i<fFourier.fPhase.size()-1; i++) {
fout << fFourier.fPhase[i] << ", ";
}
fout << fFourier.fPhase[fFourier.fPhase.size()-1] << endl;
}
// range_for_phase_correction
@@ -3879,14 +3896,306 @@ void PMsrHandler::InitFourierParameterStructure(PMsrFourierStructure &fourier)
fourier.fDCCorrected = false; // dc-corrected FFT, default: false
fourier.fApodization = FOURIER_APOD_NOT_GIVEN; // apodization, default: NOT GIVEN
fourier.fPlotTag = FOURIER_PLOT_NOT_GIVEN; // initial plot tag, default: NOT GIVEN
fourier.fPhaseParamNo = 0; // initial parameter no = 0 means not a parameter
fourier.fPhase = -999.0; // fourier phase: -999 = NOT GIVEN
fourier.fPhaseRef = -1; // initial phase reference -1 means: use absolute phases
fourier.fPhaseParamNo.clear(); // initial phase parameter no vector is empty
fourier.fPhase.clear(); // initial phase vector is empty
for (UInt_t i=0; i<2; i++) {
fourier.fRangeForPhaseCorrection[i] = -1.0; // frequency range for phase correction, default: {-1, -1} = NOT GIVEN
fourier.fPlotRange[i] = -1.0; // fourier plot range, default: {-1, -1} = NOT GIVEN
fourier.fRangeForPhaseCorrection[i] = -1.0; // frequency range for phase correction, default: {-1, -1} = NOT GIVEN
fourier.fPlotRange[i] = -1.0; // fourier plot range, default: {-1, -1} = NOT GIVEN
}
}
//--------------------------------------------------------------------------
// RemoveComment (private)
//--------------------------------------------------------------------------
/**
* <p>Removes a potentially present comment from str and returns the truncated
* string in truncStr. A comment starts with '#'
*
* @param str original string which might contain a comment
* @param truncStr string from which the comment has been removed
*/
void PMsrHandler::RemoveComment(const TString &str, TString &truncStr)
{
truncStr = str;
Ssiz_t idx = str.First('#'); // find the index of the comment character
// truncate string if comment is found
if (idx > 0) {
truncStr.Resize(idx-1);
}
}
//--------------------------------------------------------------------------
// ParseFourierPhaseValueVector (private)
//--------------------------------------------------------------------------
/**
* <p>examines if str has the form 'phase val0 [sep val1 ... sep valN]'.
* If this form is found, fill in val0 ... valN to fFourier.fPhase
* vector.
*
* @param fourier msr-file Fourier structure
* @param str string to be analyzed
* @param error flag needed to propagate a fatal error
*
* @return true if a phase value form is found, otherwise return false
*/
Bool_t PMsrHandler::ParseFourierPhaseValueVector(PMsrFourierStructure &fourier, const TString &str, Bool_t &error)
{
Bool_t result = true;
TObjArray *tok = str.Tokenize(" ,;\t");
if (tok == 0) {
cerr << endl << ">> PMsrHandler::ParseFourierPhaseValueVector: **ERROR** couldn't tokenize Fourier phase line." << endl << endl;
return false;
}
// make sure there are enough tokens
if (tok->GetEntries() < 2) {
error = true;
return false;
}
// convert all acceptable tokens
TObjString *ostr=0;
TString sstr("");
for (Int_t i=1; i<tok->GetEntries(); i++) {
ostr = dynamic_cast<TObjString*>(tok->At(i));
sstr = ostr->GetString();
if (sstr.IsFloat()) {
fourier.fPhase.push_back(sstr.Atof());
} else {
result = false;
if (i>1) { // make sure that no 'phase val, parX' mixture is present
cerr << endl << ">> PMsrHandler::ParseFourierPhaseValueVector: **ERROR** in Fourier phase line.";
cerr << endl << ">> Attempt to mix val, parX? This is currently not supported." << endl << endl;
error = true;
}
break;
}
}
// clean up
if (tok) {
delete tok;
}
return result;
}
//--------------------------------------------------------------------------
// ParseFourierPhaseParVector (private)
//--------------------------------------------------------------------------
/**
* <p> examines if str has the form 'phase parX0 [sep parX1 ... sep parXN]'.
* Also allowed is that instead of parXn only one of the parameters could have the
* form parRn which markes a reference phase for relative phase fittings.
* If this form is found, fill in parX0 ... parXN to fourier.fPhaseParamNo, and
* in case a parR is present, set the fourier.fPhaseRef accordingly.
*
* @param fourier msr-file Fourier structure
* @param str string to be analyzed
* @param error flag needed to propagate a fatal error
*
* @return true if a phase parameter form is found, otherwise return false
*/
Bool_t PMsrHandler::ParseFourierPhaseParVector(PMsrFourierStructure &fourier, const TString &str, Bool_t &error)
{
Bool_t result = true;
Int_t refCount = 0;
TObjArray *tok = str.Tokenize(" ,;\t");
if (tok == 0) {
cerr << endl << ">> PMsrHandler::ParseFourierPhaseParVector: **ERROR** couldn't tokenize Fourier phase line." << endl << endl;
return false;
}
// make sure there are enough tokens
if (tok->GetEntries() < 2) {
error = true;
return false;
}
// check that all tokens start with par
TString sstr;
for (Int_t i=1; i<tok->GetEntries(); i++) {
TObjString *ostr = dynamic_cast<TObjString*>(tok->At(i));
sstr = ostr->GetString();
if (!sstr.BeginsWith("par")) {
cerr << ">> PMsrHandler::ParseFourierPhaseParVector: **ERROR** found unhandable token '" << sstr << "'" << endl;
error = true;
result = false;
break;
}
if (sstr.BeginsWith("parR")) {
refCount++;
}
// rule out par(X, offset, #Param) syntax
if (sstr.BeginsWith("par(")) {
result = false;
break;
}
}
if (refCount > 1) {
cerr << ">> PMsrHandler::ParseFourierPhaseParVector: **ERROR** found multiple parR's! Only one reference phase is accepted." << endl;
result = false;
}
// check that token has the form parX, where X is an int
Int_t rmNoOf = 3;
if (result != false) {
for (Int_t i=1; i<tok->GetEntries(); i++) {
TObjString *ostr = dynamic_cast<TObjString*>(tok->At(i));
sstr = ostr->GetString();
rmNoOf = 3;
if (sstr.BeginsWith("parR")) {
rmNoOf++;
}
sstr.Remove(0, rmNoOf); // remove 'par' of 'parR' part. Rest should be an integer
if (sstr.IsDigit()) {
if (rmNoOf == 4) // parR
fourier.fPhaseRef = sstr.Atoi();
fourier.fPhaseParamNo.push_back(sstr.Atoi());
} else {
cerr << ">> PMsrHandler::ParseFourierPhaseParVector: **ERROR** found token '" << ostr->GetString() << "' which is not parX with X an integer." << endl;
fourier.fPhaseParamNo.clear();
error = true;
break;
}
}
}
if (fourier.fPhaseParamNo.size() == tok->GetEntries()-1) { // everything as expected
result = true;
} else {
result = false;
}
// clean up
if (tok) {
delete tok;
}
return result;
}
//--------------------------------------------------------------------------
// ParseFourierPhaseParIterVector (private)
//--------------------------------------------------------------------------
/**
* <p> examines if str has the form 'phase par(X0, offset, #params)' or 'phase parR(X0, offset, #params)'.
* If this form is found, fill in parX0 ... parXN to fourier.fPhaseParamNo, and
* in case of 'parR' also set the fourier.fPhaseRef accordingly.
*
* @param fourier msr-file Fourier structure
* @param str string to be analyzed
* @param error flag needed to propagate a fatal error
*
* @return true if a phase parameter iterator form is found, otherwise return false
*/
Bool_t PMsrHandler::ParseFourierPhaseParIterVector(PMsrFourierStructure &fourier, const TString &str, Bool_t &error)
{
TString wstr = str;
// remove 'phase' from string
wstr.Remove(0, 5);
wstr = wstr.Strip(TString::kLeading, ' ');
// remove 'par(' from string if present, otherwise and error is issued
if (!wstr.BeginsWith("par(") && !wstr.BeginsWith("parR(")) {
cout << ">> PMsrHandler::ParseFourierPhaseParIterVector: **ERROR** token should start with 'par(' or 'parR(', found: '" << wstr << "' -> ERROR" << endl;
error = true;
return false;
}
Int_t noOf = 4; // number of characters to be removed
Bool_t relativePhase = false; // relative phase handling wished
if (wstr.BeginsWith("parR(")) {
noOf += 1;
relativePhase = true;
}
wstr.Remove(0, noOf);
// remove trailing white spaces
wstr = wstr.Strip(TString::kTrailing, ' ');
// remove last ')'
Ssiz_t idx=wstr.Last(')');
wstr.Remove(idx, wstr.Length()-idx);
// tokenize rest which should have the form 'X0, offset, #Param'
TObjArray *tok = wstr.Tokenize(",;");
if (tok == 0) {
cerr << ">> PMsrHandler::ParseFourierPhaseParIterVector: **ERROR** tokenize failed." << endl;
error = true;
return false;
}
// check for proper number of expected elements
if (tok->GetEntries() != 3) {
cerr << ">> PMsrHandler::ParseFourierPhaseParIterVector: **ERROR** wrong syntax for the expected par(X0, offset, #param)." << endl;
error = true;
delete tok;
return false;
}
Int_t x0, offset, noParam;
// get X0
TObjString *ostr = dynamic_cast<TObjString*>(tok->At(0));
wstr = ostr->GetString();
if (wstr.IsDigit()) {
x0 = wstr.Atoi();
} else {
cerr << ">> PMsrHandler::ParseFourierPhaseParIterVector: **ERROR** X0='" << wstr << "' is not an integer." << endl;
error = true;
delete tok;
return false;
}
// get offset
ostr = dynamic_cast<TObjString*>(tok->At(1));
wstr = ostr->GetString();
if (wstr.IsDigit()) {
offset = wstr.Atoi();
} else {
cerr << ">> PMsrHandler::ParseFourierPhaseParIterVector: **ERROR** offset='" << wstr << "' is not an integer." << endl;
error = true;
delete tok;
return false;
}
// get noParam
ostr = dynamic_cast<TObjString*>(tok->At(2));
wstr = ostr->GetString();
if (wstr.IsDigit()) {
noParam = wstr.Atoi();
} else {
cerr << ">> PMsrHandler::ParseFourierPhaseParIterVector: **ERROR** #Param='" << wstr << "' is not an integer." << endl;
error = true;
delete tok;
return false;
}
// set the reference phase parameter number for 'parR'
if (relativePhase)
fourier.fPhaseRef = x0;
else
fourier.fPhaseRef = -1;
for (Int_t i=0; i<noParam; i++)
fourier.fPhaseParamNo.push_back(x0 + i*offset);
// clean up
if (tok) {
delete tok;
}
return true;
}
//--------------------------------------------------------------------------
// HandleFourierEntry (private)
//--------------------------------------------------------------------------
@@ -4031,39 +4340,66 @@ Bool_t PMsrHandler::HandleFourierEntry(PMsrLines &lines)
}
}
} else if (iter->fLine.BeginsWith("phase", TString::kIgnoreCase)) { // phase
if (tokens->GetEntries() < 2) { // phase value is missing
if (tokens->GetEntries() < 2) { // phase value(s)/par(s) is(are) missing
error = true;
continue;
} else {
ostr = dynamic_cast<TObjString*>(tokens->At(1));
str = ostr->GetString();
if (str.BeginsWith("par", TString::kIgnoreCase)) { // parameter value
if (fFourierOnly) {
cerr << endl << ">> PMsrHandler::HandleFourierEntry: **WARNING** Found phase parameter for Fourier only.";
cerr << endl << ">> This is currently not supported. Will set the phase to 0." << endl;
fourier.fPhase = 0.0;
} else {
Int_t no = 0;
if (FilterNumber(str, "par", 0, no)) {
// check that the parameter is in range
if ((Int_t)fParam.size() < no) {
error = true;
continue;
}
// keep the parameter number
fourier.fPhaseParamNo = no;
// get parameter value
fourier.fPhase = fParam[no-1].fValue;
} else {
// allowed phase parameter patterns:
// (i) phase val [sep val sep val ...] [# comment], val=double, sep=' ,;\t'
// (ii) phase parX0 [sep parX1 sep parX2 ...] [# comment], val=double, sep=' ,;\t'
// (iii) phase par(X0 sep1 offset sep1 #param) [# comment], sep1= ',;'
// remove potential comment
TString wstr("");
RemoveComment(iter->fLine, wstr);
// check for 'phase val ...'
Bool_t result = ParseFourierPhaseValueVector(fourier, wstr, error);
if (error)
continue;
// check for 'phase parX0 ...' if not already val are found
if (!result) {
result = ParseFourierPhaseParVector(fourier, wstr, error);
if (error)
continue;
}
// check for 'phase par(X0, offset, #param)' if not already covered by the previous ones
if (!result) {
result = ParseFourierPhaseParIterVector(fourier, wstr, error);
}
if (!result || error) {
continue;
}
// if parameter vector is given: check that all parameters are within range
if (fourier.fPhaseParamNo.size() > 0) {
for (UInt_t i=0; i<fourier.fPhaseParamNo.size(); i++) {
if (fourier.fPhaseParamNo[i] > fParam.size()) {
cerr << ">> PMsrHandler::HandleFourierEntry: found Fourier parameter entry par" << fourier.fPhaseParamNo[i] << " > #Param = " << fParam.size() << endl;
error = true;
--iter;
continue;
}
}
} else if (str.IsFloat()) { // phase value
fourier.fPhase = str.Atof();
} else {
error = true;
continue;
}
// if parameter vector is given -> fill corresponding phase values
Double_t phaseRef = 0.0;
if (fourier.fPhaseParamNo.size() > 0) {
// check if a relative parameter phase number is set
if (fourier.fPhaseRef != -1) {
phaseRef = fParam[fourier.fPhaseRef-1].fValue;
}
fourier.fPhase.clear();
for (UInt_t i=0; i<fourier.fPhaseParamNo.size(); i++) {
if (fourier.fPhaseRef == fourier.fPhaseParamNo[i]) // reference phase
fourier.fPhase.push_back(fParam[fourier.fPhaseParamNo[i]-1].fValue);
else
fourier.fPhase.push_back(fParam[fourier.fPhaseParamNo[i]-1].fValue+phaseRef);
}
}
}
} else if (iter->fLine.BeginsWith("range_for_phase_correction", TString::kIgnoreCase)) {
@@ -4162,7 +4498,11 @@ Bool_t PMsrHandler::HandleFourierEntry(PMsrLines &lines)
cerr << endl << ">> [dc-corrected true | false]";
cerr << endl << ">> [apodization none | weak | medium | strong]";
cerr << endl << ">> [plot real | imag | real_and_imag | power | phase | phase_opt_real]";
cerr << endl << ">> [phase value]";
cerr << endl << ">> [phase valList | parList | parIterList [# comment]]";
cerr << endl << ">> valList : val [sep val ... sep val]. sep=' ,;\\t'";
cerr << endl << ">> parList : parX0 [sep parX1 ... sep parX1]";
cerr << endl << ">> parIterList : par(X0,offset,#param), with X0=first parameter number";
cerr << endl << ">> offset=parameter offset, #param=number of phase parameters.";
cerr << endl << ">> [range_for_phase_correction min max | all]";
cerr << endl << ">> [range min max]";
cerr << endl;
@@ -6294,6 +6634,75 @@ void PMsrHandler::MakeDetectorGroupingString(TString str, PIntVector &group, TSt
} while (i<group.size());
}
//--------------------------------------------------------------------------
// BeautifyFourierPhaseParameterString (private)
//--------------------------------------------------------------------------
/**
* <p>Returns the Fourier phase string if the phase is either of type
* phase parX0 sep ... sep parXn where sep = ',' or
* phase par(X0, offset, #param)
*
* @return Fourier phase parameter string if phase parameter(s) is(are) given, "??" otherwise
*/
TString PMsrHandler::BeautifyFourierPhaseParameterString()
{
TString str("??");
TString formatStr("par%d, par%d");
if (fFourier.fPhaseParamNo.size() == 0)
return str;
Int_t phaseRef = fFourier.fPhaseRef;
if (fFourier.fPhaseParamNo.size() == 1) {
str = TString::Format("par%d", fFourier.fPhaseParamNo[0]);
} else if (fFourier.fPhaseParamNo.size() == 2) {
if (phaseRef == fFourier.fPhaseParamNo[0])
formatStr = "parR%d, par%d";
if (phaseRef == fFourier.fPhaseParamNo[1])
formatStr = "par%d, parR%d";
str = TString::Format(formatStr, fFourier.fPhaseParamNo[0], fFourier.fPhaseParamNo[1]);
} else {
Bool_t phaseIter = true;
// first check if fPhaseParamNo vector can be compacted into par(X0, offset, #param) form
Int_t offset = fFourier.fPhaseParamNo[1] - fFourier.fPhaseParamNo[0];
for (Int_t i=2; i<fFourier.fPhaseParamNo.size(); i++) {
if (fFourier.fPhaseParamNo[i]-fFourier.fPhaseParamNo[i-1] != offset) {
phaseIter = false;
break;
}
}
if (phaseIter) {
if (phaseRef != -1) {
str = TString::Format("parR(%d, %d, %d)", fFourier.fPhaseParamNo[0], offset, fFourier.fPhaseParamNo.size());
} else {
str = TString::Format("par(%d, %d, %d)", fFourier.fPhaseParamNo[0], offset, fFourier.fPhaseParamNo.size());
}
} else {
str = TString("");
for (Int_t i=0; i<fFourier.fPhaseParamNo.size()-1; i++) {
if (phaseRef == fFourier.fPhaseParamNo[i]) {
str += "parR";
} else {
str += "par";
}
str += fFourier.fPhaseParamNo[i];
str += ", ";
}
if (phaseRef == fFourier.fPhaseParamNo[fFourier.fPhaseParamNo.size()-1]) {
str += "parR";
} else {
str += "par";
}
str += fFourier.fPhaseParamNo[fFourier.fPhaseParamNo.size()-1];
}
}
return str;
}
//--------------------------------------------------------------------------
// CheckLegacyLifetimecorrection (private)
//--------------------------------------------------------------------------

65
src/classes/PMusrCanvas.cpp
View File

@@ -28,6 +28,7 @@
***************************************************************************/
#include <iostream>
#include <iomanip>
#include <fstream>
using namespace std;
@@ -152,7 +153,6 @@ PMusrCanvas::PMusrCanvas()
InitFourier();
InitAverage();
fCurrentFourierPhase = fFourier.fPhaseIncrement;
fCurrentFourierPhaseText = 0;
fRRFText = 0;
@@ -205,7 +205,6 @@ PMusrCanvas::PMusrCanvas(const Int_t number, const Char_t* title,
CreateStyle();
InitMusrCanvas(title, wtopx, wtopy, ww, wh);
fCurrentFourierPhase = 0.0;
fCurrentFourierPhaseText = 0;
fRRFText = 0;
@@ -262,7 +261,6 @@ PMusrCanvas::PMusrCanvas(const Int_t number, const Char_t* title,
CreateStyle();
InitMusrCanvas(title, wtopx, wtopy, ww, wh);
fCurrentFourierPhase = 0.0;
fCurrentFourierPhaseText = 0;
fRRFText = 0;
@@ -322,10 +320,6 @@ PMusrCanvas::~PMusrCanvas()
CleanupDataSet(fNonMusrData[i]);
fNonMusrData.clear();
}
if (fCurrentFourierPhaseText) {
delete fCurrentFourierPhaseText;
fCurrentFourierPhaseText = 0;
}
if (fMultiGraphLegend) {
fMultiGraphLegend->Clear();
delete fMultiGraphLegend;
@@ -394,9 +388,8 @@ void PMusrCanvas::SetMsrHandler(PMsrHandler *msrHandler)
if (fMsrHandler->GetMsrFourierList()->fPlotTag != FOURIER_PLOT_NOT_GIVEN) {
fFourier.fPlotTag = fMsrHandler->GetMsrFourierList()->fPlotTag;
}
if (fMsrHandler->GetMsrFourierList()->fPhase != -999.0) {
fFourier.fPhase = fMsrHandler->GetMsrFourierList()->fPhase;
}
fFourier.fPhase = fMsrHandler->GetMsrFourierList()->fPhase;
if ((fMsrHandler->GetMsrFourierList()->fRangeForPhaseCorrection[0] != -1.0) &&
(fMsrHandler->GetMsrFourierList()->fRangeForPhaseCorrection[1] != -1.0)) {
fFourier.fRangeForPhaseCorrection[0] = fMsrHandler->GetMsrFourierList()->fRangeForPhaseCorrection[0];
@@ -2221,10 +2214,10 @@ void PMusrCanvas::ExportData(const Char_t *fileName)
// write data/theory
for (UInt_t j=0; j<dumpVector.size()-1; j++) {
if (i<dumpVector[j].dataX.size()) {
fout << dumpVector[j].dataX[i] << ", ";
fout << dumpVector[j].data[i] << ", ";
fout << setprecision(9) << dumpVector[j].dataX[i] << ", ";
fout << setprecision(9) << dumpVector[j].data[i] << ", ";
if (dumpVector[j].dataErr.size() > 0)
fout << dumpVector[j].dataErr[i] << ", ";
fout << setprecision(9) << dumpVector[j].dataErr[i] << ", ";
} else {
if (dumpVector[j].dataErr.size() > 0)
fout << " , , , ";
@@ -2234,8 +2227,8 @@ void PMusrCanvas::ExportData(const Char_t *fileName)
}
// write last data/theory entry
if (i<dumpVector[dumpVector.size()-1].dataX.size()) {
fout << dumpVector[dumpVector.size()-1].dataX[i] << ", ";
fout << dumpVector[dumpVector.size()-1].data[i];
fout << setprecision(9) << dumpVector[dumpVector.size()-1].dataX[i] << ", ";
fout << setprecision(9) << dumpVector[dumpVector.size()-1].data[i];
} else {
fout << " , ";
}
@@ -2326,7 +2319,8 @@ void PMusrCanvas::InitFourier()
fFourier.fFourierPower = 0; // no zero padding
fFourier.fApodization = FOURIER_APOD_NONE; // no apodization
fFourier.fPlotTag = FOURIER_PLOT_REAL_AND_IMAG; // initial plot tag, plot real and imaginary part
fFourier.fPhase = 0.0; // fourier phase 0°
fFourier.fPhaseParamNo.clear();
fFourier.fPhase.clear();
for (UInt_t i=0; i<2; i++) {
fFourier.fRangeForPhaseCorrection[i] = -1.0; // frequency range for phase correction, default: {-1, -1} = NOT GIVEN
fFourier.fPlotRange[i] = -1.0; // fourier plot range, default: {-1, -1} = NOT GIVEN
@@ -3480,13 +3474,20 @@ void PMusrCanvas::HandleFourier()
}
// apply global phase if present
if (fFourier.fPhase != 0.0) {
const double cp = TMath::Cos(fFourier.fPhase/180.0*TMath::Pi());
const double sp = TMath::Sin(fFourier.fPhase/180.0*TMath::Pi());
if (fFourier.fPhase.size() != 0.0) {
double cp;
double sp;
fCurrentFourierPhase = fFourier.fPhase;
for (UInt_t i=0; i<fData.size(); i++) { // loop over all data sets
if (fFourier.fPhase.size() == 1) {
cp = TMath::Cos(fFourier.fPhase[0]/180.0*TMath::Pi());
sp = TMath::Sin(fFourier.fPhase[0]/180.0*TMath::Pi());
} else {
cp = TMath::Cos(fFourier.fPhase[i]/180.0*TMath::Pi());
sp = TMath::Sin(fFourier.fPhase[i]/180.0*TMath::Pi());
}
if ((fData[i].dataFourierRe != 0) && (fData[i].dataFourierIm != 0)) {
for (Int_t j=0; j<fData[i].dataFourierRe->GetNbinsX(); j++) { // loop over a fourier data set
// calculate new fourier data set value
@@ -3583,16 +3584,23 @@ void PMusrCanvas::HandleDifferenceFourier()
fData[i].diffFourierTag = 1; // d-f
}
// apply global phase
if (fFourier.fPhase != 0.0) {
// apply phase
if (fFourier.fPhase.size() != 0.0) {
double re, im;
const double cp = TMath::Cos(fFourier.fPhase/180.0*TMath::Pi());
const double sp = TMath::Sin(fFourier.fPhase/180.0*TMath::Pi());
double cp;
double sp;
fCurrentFourierPhase = fFourier.fPhase;
for (UInt_t i=0; i<fData.size(); i++) { // loop over all data sets
if ((fData[i].diffFourierRe != 0) && (fData[i].diffFourierIm != 0)) {
if (fFourier.fPhase.size() == 1) {
cp = TMath::Cos(fFourier.fPhase[0]/180.0*TMath::Pi());
sp = TMath::Sin(fFourier.fPhase[0]/180.0*TMath::Pi());
} else {
cp = TMath::Cos(fFourier.fPhase[i]/180.0*TMath::Pi());
sp = TMath::Sin(fFourier.fPhase[i]/180.0*TMath::Pi());
}
for (Int_t j=0; j<fData[i].diffFourierRe->GetNbinsX(); j++) { // loop over a fourier data set
// calculate new fourier data set value
re = fData[i].diffFourierRe->GetBinContent(j) * cp + fData[i].diffFourierIm->GetBinContent(j) * sp;
@@ -6074,7 +6082,10 @@ void PMusrCanvas::PlotFourierPhaseValue(Bool_t unzoom)
// plot Fourier phase
str = TString("phase = ");
str += fCurrentFourierPhase;
str += fCurrentFourierPhase[0];
if (fFourier.fPhase.size() > 1) { // if more than one phase is present, do NOT plot phase info
str = TString("");
}
x = 0.7;
y = 0.85;
fCurrentFourierPhaseText = new TLatex();
@@ -6336,7 +6347,8 @@ void PMusrCanvas::IncrementFourierPhase()
const double cp = TMath::Cos(fFourier.fPhaseIncrement/180.0*TMath::Pi());
const double sp = TMath::Sin(fFourier.fPhaseIncrement/180.0*TMath::Pi());
fCurrentFourierPhase += fFourier.fPhaseIncrement;
for (UInt_t i=0; i<fCurrentFourierPhase.size(); i++)
fCurrentFourierPhase[i] += fFourier.fPhaseIncrement;
PlotFourierPhaseValue();
for (UInt_t i=0; i<fData.size(); i++) { // loop over all data sets
@@ -6388,7 +6400,8 @@ void PMusrCanvas::DecrementFourierPhase()
const double cp = TMath::Cos(fFourier.fPhaseIncrement/180.0*TMath::Pi());
const double sp = TMath::Sin(fFourier.fPhaseIncrement/180.0*TMath::Pi());
fCurrentFourierPhase -= fFourier.fPhaseIncrement;
for (UInt_t i=0; i<fCurrentFourierPhase.size(); i++)
fCurrentFourierPhase[i] -= fFourier.fPhaseIncrement;
PlotFourierPhaseValue();
for (UInt_t i=0; i<fData.size(); i++) { // loop over all data sets

4
src/classes/PRunAsymmetry.cpp
View File

@@ -1640,8 +1640,8 @@ Bool_t PRunAsymmetry::GetProperT0(PRawRunData* runData, PMsrGlobalBlock *globalB
}
// fill in the T0's from the msr-file (if present)
for (Int_t j=0; j<fRunInfo->GetAddT0BinSize(i); j++) {
fAddT0s[i-1][j] = fRunInfo->GetAddT0Bin(i, j);
for (Int_t j=0; j<fRunInfo->GetAddT0BinSize(i-1); j++) {
fAddT0s[i-1][j] = fRunInfo->GetAddT0Bin(i-1, j);
}
// fill in the T0's from the data file, if not already present in the msr-file

4
src/classes/PRunAsymmetryRRF.cpp
View File

@@ -1278,8 +1278,8 @@ Bool_t PRunAsymmetryRRF::GetProperT0(PRawRunData* runData, PMsrGlobalBlock *glob
}
// fill in the T0's from the msr-file (if present)
for (Int_t j=0; j<fRunInfo->GetAddT0BinSize(i); j++) {
fAddT0s[i-1][j] = fRunInfo->GetAddT0Bin(i, j);
for (Int_t j=0; j<fRunInfo->GetAddT0BinSize(i-1); j++) {
fAddT0s[i-1][j] = fRunInfo->GetAddT0Bin(i-1, j);
}
// fill in the T0's from the data file, if not already present in the msr-file

3
src/classes/PStartupHandler.cpp
View File

@@ -192,7 +192,6 @@ void PStartupHandler::OnStartDocument()
fFourierDefaults.fFourierPower = 0;
fFourierDefaults.fApodization = FOURIER_APOD_NONE;
fFourierDefaults.fPlotTag = FOURIER_PLOT_REAL_AND_IMAG;
fFourierDefaults.fPhase = 0.0;
fFourierDefaults.fRangeForPhaseCorrection[0] = -1.0;
fFourierDefaults.fRangeForPhaseCorrection[1] = -1.0;
fFourierDefaults.fPlotRange[0] = -1.0;
@@ -412,7 +411,7 @@ void PStartupHandler::OnCharacters(const Char_t *str)
case ePhase:
tstr = TString(str);
if (tstr.IsFloat()) {
fFourierDefaults.fPhase = tstr.Atof();
fFourierDefaults.fPhase.push_back(tstr.Atof());
} else {
cerr << endl << "PStartupHandler **WARNING** '" << str << "' is not a valid phase, will ignore it.";
cerr << endl;

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

@@ -2,9 +2,10 @@
#--- generate necessary dictionaries ------------------------------------------
set(MUSRFIT_INC ${CMAKE_SOURCE_DIR}/src/include)
root_generate_dictionary(
PMagProximityFitterDict
-I${FFTW3_INCLUDE}
-I${FFTW3_INCLUDE_DIR}
-I${MUSRFIT_INC}
PMagProximityFitter.h
LINKDEF PMagProximityFitterLinkDef.h

14
src/external/libBNMR/libLineProfile/CMakeLists.txt vendored
View File

@@ -4,18 +4,18 @@
set(MUSRFIT_INC ${CMAKE_SOURCE_DIR}/src/include)
root_generate_dictionary(
libLineProfileDict
LineProfileDict
-I${FFTW3_INCLUDE_DIR}
-I${MUSRFIT_INC}
libLineProfile.h
LINKDEF libLineProfileLinkDef.h
MODULE libLineProfile
LineProfile.h
LINKDEF LineProfileLinkDef.h
MODULE LineProfile
)
#--- lib creation -------------------------------------------------------------
add_library(LineProfile SHARED
libLineProfile.cpp
libLineProfileDict.cxx
LineProfile.cpp
LineProfileDict.cxx
)
#--- make sure that the include directory is found ----------------------------
@@ -48,7 +48,7 @@ install(
#--- install libLineProfile header ---------------------------------------------------
install(
FILES
libLineProfile.h
LineProfile.h
DESTINATION
include
)

2
src/external/libBNMR/libLineProfile/libLineProfile.cpp → src/external/libBNMR/libLineProfile/LineProfile.cpp vendored
View File

@@ -20,7 +20,7 @@
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
***************************************************************************/
#include "libLineProfile.h"
#include "LineProfile.h"
#include <iostream> //for testing purposes
//Implement helperfunctions

0
src/external/libBNMR/libLineProfile/libLineProfile.h → src/external/libBNMR/libLineProfile/LineProfile.h vendored
View File

0
src/external/libBNMR/libLineProfile/libLineProfileLinkDef.h → src/external/libBNMR/libLineProfile/LineProfileLinkDef.h vendored
View File

10
src/external/libBNMR/libLineProfile/Makefile.am vendored
View File

@@ -1,19 +1,19 @@
## Process this file with automake to create Makefile.in
h_sources = \
libLineProfile.h
LineProfile.h
h_linkdef = \
libLineProfile.h
LineProfile.h
dict_h_sources = \
libLineProfileDict.h
LineProfileDict.h
cpp_sources = \
libLineProfile.cpp
LineProfile.cpp
dict_cpp_sources = \
libLineProfileDict.cpp
LineProfileDict.cpp
include_HEADERS = $(h_sources)
noinst_HEADERS = $(h_linkdef) $(dict_h_sources)

5
src/include/PMsrHandler.h
View File

@@ -146,6 +146,10 @@ class PMsrHandler
virtual void FillParameterInUse(PMsrLines &theory, PMsrLines &funcs, PMsrLines &run);
virtual void InitFourierParameterStructure(PMsrFourierStructure &fourier);
virtual void RemoveComment(const TString &str, TString &truncStr);
virtual Bool_t ParseFourierPhaseValueVector(PMsrFourierStructure &fourier, const TString &str, Bool_t &error);
virtual Bool_t ParseFourierPhaseParVector(PMsrFourierStructure &fourier, const TString &str, Bool_t &error);
virtual Bool_t ParseFourierPhaseParIterVector(PMsrFourierStructure &fourier, const TString &str, Bool_t &error);
virtual Bool_t FilterNumber(TString str, const Char_t *filter, Int_t offset, Int_t &no);
@@ -153,6 +157,7 @@ class PMsrHandler
virtual UInt_t LastSignificant(Double_t dval, UInt_t precLimit=6);
virtual void MakeDetectorGroupingString(TString str, PIntVector &group, TString &result, Bool_t includeDetector = true);
virtual TString BeautifyFourierPhaseParameterString();
virtual void CheckLegacyLifetimecorrection();
};

7
src/include/PMusr.h
View File

@@ -8,7 +8,7 @@
***************************************************************************/
/***************************************************************************
* Copyright (C) 2007-2016 by Andreas Suter *
* Copyright (C) 2007-2018 by Andreas Suter *
* andreas.suter@psi.ch *
* *
* This program is free software; you can redistribute it and/or modify *
@@ -745,8 +745,9 @@ typedef struct {
Int_t fFourierPower; ///< i.e. zero padding up to 2^fFourierPower, default = 0 which means NO zero padding
Int_t fApodization; ///< tag indicating the kind of apodization wished, 0=no appodization (default), 1=weak, 2=medium, 3=strong (for details see the docu)
Int_t fPlotTag; ///< tag used for initial plot. 0=real, 1=imaginary, 2=real & imaginary (default), 3=power, 4=phase
Int_t fPhaseParamNo; ///< parameter number if used instead of a phase value
Double_t fPhase; ///< phase
Int_t fPhaseRef; ///< phase reference for relative phase(s)
PIntVector fPhaseParamNo; ///< parameter number(s) if used instead of a phase value
PDoubleVector fPhase; ///< phase(s)
Double_t fRangeForPhaseCorrection[2]; ///< field/frequency range for automatic phase correction
Double_t fPlotRange[2]; ///< field/frequency plot range
Double_t fPhaseIncrement; ///< phase increment for manual phase optimization

2
src/include/PMusrCanvas.h
View File

@@ -252,7 +252,7 @@ class PMusrCanvas : public TObject, public TQObject
Bool_t fXRangePresent, fYRangePresent; ///< flag indicating if x-/y-range is present
Double_t fXmin, fXmax, fYmin, fYmax; ///< data/theory frame range
Double_t fCurrentFourierPhase; ///< holds the current Fourier phase
PDoubleVector fCurrentFourierPhase; ///< holds the current Fourier phase(s)
TLatex *fCurrentFourierPhaseText; ///< used in Re/Im Fourier to show the current phase in the pad
TString *fRRFText; ///< RRF information
TLatex *fRRFLatexText; ///< used to display RRF info

3
src/musredit_qt5/musrStep/CMakeLists.txt
View File

@@ -21,6 +21,9 @@ set(CMAKE_AUTOUIC OFF)
# add qt/rcc
qt5_add_resources(musrStep_rcc musrStep.qrc)
# remove generated files from automoc and autouic
set_property(SOURCE qrc_musrStep.cpp PROPERTY SKIP_AUTOMOC ON)
set(macosx_icon icons/musrStep.icns)
if (APPLE)
add_executable(musrStep MACOSX_BUNDLE

3
src/musredit_qt5/musrWiz/CMakeLists.txt
View File

@@ -25,6 +25,9 @@ set(CMAKE_AUTOUIC OFF)
# add qt/rcc
qt5_add_resources(musrWiz_rcc musrWiz.qrc)
# remove generated files from automoc and autouic
set_property(SOURCE qrc_musrWiz.cpp PROPERTY SKIP_AUTOMOC ON)
set(macosx_icon icons/musrWiz.icns)
if (APPLE)
add_executable(musrWiz MACOSX_BUNDLE

40
src/musredit_qt5/musredit/CMakeLists.txt
View File

@@ -5,16 +5,20 @@ find_package(Qt5WebEngine QUIET)
find_package(Qt5WebKit QUIET)
set(qt_libs Qt5::Core Qt5::Widgets Qt5::Network Qt5::Xml Qt5::Svg Qt5::PrintSupport)
set(Qt5NoWeb 0)
if (Qt5WebEngine_FOUND)
message("-- Qt5WebEngine is present.")
find_package(Qt5WebEngineWidgets QUIET CONFIG REQUIRED)
set(qt_libs ${qt_libs} Qt5::WebEngine Qt5::WebEngineWidgets)
# unset a potentially found Qt5Webkit
unset(Qt5WebKit_FOUND)
elseif (Qt5WebKit_FOUND)
message("-- Qt5WebKit is present.")
find_package(Qt5WebKitWidgets QUIET CONFIG REQUIRED)
set(qt_libs ${qt_libs} Qt5::WebKit Qt5::WebKitWidgets)
else (Qt5WebEngine_FOUND)
message("-- Neither Qt5WebEngine nor Qt5WebKit found.")
set(Qt5NoWeb 1)
endif (Qt5WebEngine_FOUND)
set(musredit_src
@@ -80,6 +84,26 @@ qt5_wrap_ui(out_ui ${musredit_ui})
# add qt/rcc
qt5_add_resources(musredit_rcc musredit.qrc)
# remove generated files from automoc and autouic
set_property(SOURCE ui_PFindDialog.h PROPERTY SKIP_AUTOMOC ON)
set_property(SOURCE ui_PReplaceDialog.h PROPERTY SKIP_AUTOMOC ON)
set_property(SOURCE ui_PReplaceConfirmationDialog.h PROPERTY SKIP_AUTOMOC ON)
set_property(SOURCE ui_PMusrEditAbout.h PROPERTY SKIP_AUTOMOC ON)
set_property(SOURCE ui_PPrefsDialog.h PROPERTY SKIP_AUTOMOC ON)
set_property(SOURCE ui_PGetMusrFTOptionsDialog.h PROPERTY SKIP_AUTOMOC ON)
set_property(SOURCE ui_PGetTitleBlockDialog.h PROPERTY SKIP_AUTOMOC ON)
set_property(SOURCE ui_PGetParameterBlockDialog.h PROPERTY SKIP_AUTOMOC ON)
set_property(SOURCE ui_PGetTheoryBlockDialog.h PROPERTY SKIP_AUTOMOC ON)
set_property(SOURCE ui_PGetFunctionsBlockDialog.h PROPERTY SKIP_AUTOMOC ON)
set_property(SOURCE ui_PGetAsymmetryRunBlockDialog.h PROPERTY SKIP_AUTOMOC ON)
set_property(SOURCE ui_PGetSingleHistoRunBlockDialog.h PROPERTY SKIP_AUTOMOC ON)
set_property(SOURCE ui_PGetNonMusrRunBlockDialog.h PROPERTY SKIP_AUTOMOC ON)
set_property(SOURCE ui_PGetFourierBlockDialog.h PROPERTY SKIP_AUTOMOC ON)
set_property(SOURCE ui_PGetPlotBlockDialog.h PROPERTY SKIP_AUTOMOC ON)
set_property(SOURCE ui_PMsr2DataDialog.h PROPERTY SKIP_AUTOMOC ON)
set_property(SOURCE ui_PChangeDefaultPathsDialog.h PROPERTY SKIP_AUTOMOC ON)
set_property(SOURCE qrc_musredit.cpp PROPERTY SKIP_AUTOMOC ON)
set(macosx_icon icons/musredit.icns)
if (APPLE)
add_executable(musredit MACOSX_BUNDLE
@@ -96,6 +120,15 @@ else (APPLE)
)
endif (APPLE)
# set necessary tags depending if QtWebEngine, QtWebKit,
# or none of both are given
if (Qt5WebKit_FOUND)
target_compile_options(musredit
BEFORE PRIVATE
-DHAVE_QT_WEB_KIT
)
endif (Qt5WebKit_FOUND)
if (Qt5WebEngine_FOUND)
target_compile_options(musredit
BEFORE PRIVATE
@@ -103,6 +136,13 @@ if (Qt5WebEngine_FOUND)
)
endif (Qt5WebEngine_FOUND)
if (Qt5NoWeb)
target_compile_options(musredit
BEFORE PRIVATE
-DHAVE_QT_NO_WEB
)
endif (Qt5NoWeb)
target_include_directories(musredit
BEFORE PRIVATE
$<BUILD_INTERFACE:${CMAKE_CURRENT_BINARY_DIR}>

17
src/musredit_qt5/musredit/PChangeDefaultPathsDialog.cpp
View File

@@ -302,6 +302,17 @@ void PChangeDefaultPathsDialog::saveDefaultPathList()
}
fileIn.close();
// check if there is any data_path is present in the musrfit_startup.xml
bool dataPathPresent = false;
QString str;
for (int i=0; i<fileContent.count(); i++) {
str = fileContent[i];
if (str.trimmed().startsWith("<data_path>")) {
dataPathPresent = true;
break;
}
}
// write the new musrfit_startup.xml
QFile fileOut(fDefaultPath->getPrefPathName());
if (!fileOut.open(QIODevice::WriteOnly | QIODevice::Text)) {
@@ -309,11 +320,15 @@ void PChangeDefaultPathsDialog::saveDefaultPathList()
return;
}
QTextStream out(&fileOut);
QString str;
bool first = true;
for (int i=0; i<fileContent.count(); i++) {
str = fileContent[i];
if (!str.trimmed().startsWith("<data_path>")) {
// if not data_path was present, add the new data_paths just before the end of the musrfit_start.xml close tag
if ((dataPathPresent == false) && (str.trimmed().startsWith("</musrfit>"))) {
for (int j=0; j<fSearchPath_listWidget->count(); j++)
out << " <data_path>" << fSearchPath_listWidget->item(j)->text() << "</data_path>" << endl;
}
out << fileContent[i] << endl;
} else {
if (first) {

31
src/musredit_qt5/musredit/PHelp.cpp
View File

@@ -30,9 +30,13 @@
#include <QtWidgets>
#ifdef HAVE_QT_WEB_ENGINE
#include <QWebEngineView>
#else
#endif
#ifdef HAVE_QT_WEB_KIT
#include <QtWebKitWidgets>
#endif
#ifdef HAVE_QT_NO_WEB
#include <QPlainTextEdit>
#endif
#include <QNetworkProxyFactory>
#include <QtDebug>
@@ -61,14 +65,24 @@ PHelp::PHelp(const QString &url, const bool isDarkTheme) :
#ifdef HAVE_QT_WEB_ENGINE
fView = new QWebEngineView(this);
#else
#endif
#ifdef HAVE_QT_WEB_KIT
fView = new QWebView(this);
#endif
#ifdef HAVE_QT_NO_WEB
fView = new QPlainTextEdit(this);
#endif
#ifndef HAVE_QT_NO_WEB
fView->load(QUrl(url));
connect(fView, SIGNAL(loadFinished(bool)), SLOT(adjustLocation()));
connect(fView, SIGNAL(titleChanged(QString)), SLOT(adjustTitle()));
connect(fView, SIGNAL(loadProgress(int)), SLOT(setProgress(int)));
connect(fView, SIGNAL(loadFinished(bool)), SLOT(finishLoading(bool)));
#else
fView->setPlainText("Within the current setup there is NO Help available.");
fView->appendPlainText("The necessary Qt web libs where not found when setting up musredit.");
fView->setReadOnly(true);
#endif
fLocationEdit = new QLineEdit(this);
fLocationEdit->setSizePolicy(QSizePolicy::Expanding, fLocationEdit->sizePolicy().verticalPolicy());
@@ -80,12 +94,13 @@ PHelp::PHelp(const QString &url, const bool isDarkTheme) :
toolBar->addAction(fView->pageAction(QWebEnginePage::Forward));
toolBar->addAction(fView->pageAction(QWebEnginePage::Reload));
toolBar->addAction(fView->pageAction(QWebEnginePage::Stop));
#else
#endif
#ifdef HAVE_QT_WEB_KIT
toolBar->addAction(fView->pageAction(QWebPage::Back));
toolBar->addAction(fView->pageAction(QWebPage::Forward));
toolBar->addAction(fView->pageAction(QWebPage::Reload));
toolBar->addAction(fView->pageAction(QWebPage::Stop));
#endif
#endif
toolBar->addWidget(fLocationEdit);
QMenu *exitMenu = menuBar()->addMenu(tr("&File"));
@@ -127,7 +142,9 @@ void PHelp::done()
*/
void PHelp::adjustLocation()
{
#ifndef HAVE_QT_NO_WEB
fLocationEdit->setText(fView->url().toString());
#endif
}
//---------------------------------------------------------------------------
@@ -138,8 +155,10 @@ void PHelp::adjustLocation()
void PHelp::changeLocation()
{
QUrl url = QUrl(fLocationEdit->text());
#ifndef HAVE_QT_NO_WEB
fView->load(url);
fView->setFocus();
#endif
}
//---------------------------------------------------------------------------
@@ -149,10 +168,12 @@ void PHelp::changeLocation()
*/
void PHelp::adjustTitle()
{
#ifndef HAVE_QT_NO_WEB
if (fProgress <= 0 || fProgress >= 100)
setWindowTitle(fView->title());
setWindowTitle(fView->title());
else
setWindowTitle(QString("%1 (%2%)").arg(fView->title()).arg(fProgress));
#endif
}
//---------------------------------------------------------------------------

9
src/musredit_qt5/musredit/PHelp.h
View File

@@ -34,7 +34,8 @@
#ifdef HAVE_QT_WEB_ENGINE
class QWebEngineView;
#else
#endif
#ifdef HAVE_QT_WEB_KIT
class QWebView;
#endif
QT_BEGIN_NAMESPACE
@@ -66,8 +67,12 @@ class PHelp : public QMainWindow
bool fDarkTheme;
#ifdef HAVE_QT_WEB_ENGINE
QWebEngineView *fView; ///< web viewer
#else
#endif
#ifdef HAVE_QT_WEB_KIT
QWebView *fView; ///< web viewer
#endif
#ifdef HAVE_QT_NO_WEB
QPlainTextEdit *fView; ///< dialog stating that there is NO web viewer
#endif
QLineEdit *fLocationEdit; ///< url address line edit
int fProgress; ///< progress value (0-100) while loading an url

25
src/tests/MuTransition/PSimulateMuTransition.cpp
View File

@@ -100,10 +100,13 @@ PSimulateMuTransition::PSimulateMuTransition(UInt_t seed)
}
fNmuons = 100; // number of muons to simulate
fNshowProgress = 100; // print progress on screen every fNshowProgress events
fMuPrecFreq34 = 4463.; // vacuum Mu hyperfine coupling constant
fMuPrecFreq12 = 0.; // Mu precession frequency of a 12 transition
fMuPrecFreq23 = 0.; // Mu precession frequency of a 23 transition
fMuPrecFreq14 = 0.; // Mu precession frequency of a 14 transition
fMuPrecFreq13 = 0.; // Mu precession frequency of a 13 transition
fMuPrecFreq24 = 0.; // Mu precession frequency of a 24 transition
fMuonPrecFreq = 0.; // muon precession frequency
fBfield = 0.01; // magnetic field (T)
fCaptureRate = 0.01; // Mu+ capture rate (MHz)
@@ -118,6 +121,9 @@ PSimulateMuTransition::PSimulateMuTransition(UInt_t seed)
fMuFractionState34 = 0.25;
fMuFractionState23 = 0.25;
fMuFractionState14 = 0.25;
fMuFractionState13 = 0.;
fMuFractionState24 = 0.;
fDebugFlag = kFALSE;
}
@@ -147,6 +153,8 @@ void PSimulateMuTransition::PrintSettings() const
cout << endl << "Mu0 precession frequency 34 (MHz) = " << fMuPrecFreq34;
cout << endl << "Mu0 precession frequency 23 (MHz) = " << fMuPrecFreq23;
cout << endl << "Mu0 precession frequency 14 (MHz) = " << fMuPrecFreq14;
cout << endl << "Mu0 precession frequency 13 (MHz) = " << fMuPrecFreq13;
cout << endl << "Mu0 precession frequency 24 (MHz) = " << fMuPrecFreq24;
cout << endl << "Mu+ precession frequency (MHz) = " << fMuonGyroRatio * fBfield;
cout << endl << "B field (T) = " << fBfield;
cout << endl << "Mu+ electron capture rate (MHz) = " << fCaptureRate;
@@ -164,7 +172,10 @@ void PSimulateMuTransition::PrintSettings() const
cout << endl << "Muonium fraction state34 = " << fMuFractionState34;
cout << endl << "Muonium fraction state23 = " << fMuFractionState23;
cout << endl << "Muonium fraction state14 = " << fMuFractionState14;
cout << endl << "Muonium fraction state13 = " << fMuFractionState13;
cout << endl << "Muonium fraction state24 = " << fMuFractionState24;
cout << endl << "Number of particles to simulate = " << fNmuons;
cout << endl << "Print progress on screen frequency = " << fNshowProgress;
cout << endl << "Initial muon spin phase (degree) = " << fInitialPhase;
cout << endl << "Debug flag = " << fDebugFlag;
cout << endl << endl;
@@ -224,7 +235,7 @@ void PSimulateMuTransition::Run(TH1F *histoForward, TH1F *histoBackward)
else
histoBackward->Fill(fMuonDecayTime*1000., 1. - fAsymmetry*TMath::Cos(fMuonPhase));
if ( (i%100000) == 0) cout << "number of events processed: " << i << endl;
if ( (i%fNshowProgress) == 0) cout << "number of events processed: " << i << endl;
}
cout << "number of events processed: " << i << endl;
return;
@@ -291,11 +302,14 @@ TComplex PSimulateMuTransition::GTFunction(const Double_t &time, const TString c
complexPol = TComplex::Exp(-TComplex::I()*twoPi*fMuonPrecFreq*time);
else{
complexPol =
(fMuFractionState12 * TComplex::Exp(TComplex::I()*twoPi*fMuPrecFreq12*time) +
(fMuFractionState12 * TComplex::Exp(+TComplex::I()*twoPi*fMuPrecFreq12*time) +
fMuFractionState34 * TComplex::Exp(-TComplex::I()*twoPi*fMuPrecFreq34*time))
+
(fMuFractionState23 * TComplex::Exp(TComplex::I()*twoPi*fMuPrecFreq23*time) +
fMuFractionState14 * TComplex::Exp(TComplex::I()*twoPi*fMuPrecFreq14*time));
(fMuFractionState23 * TComplex::Exp(+TComplex::I()*twoPi*fMuPrecFreq23*time) +
fMuFractionState14 * TComplex::Exp(+TComplex::I()*twoPi*fMuPrecFreq14*time))
+
(fMuFractionState13 * TComplex::Exp(+TComplex::I()*twoPi*fMuPrecFreq13*time) +
fMuFractionState24 * TComplex::Exp(+TComplex::I()*twoPi*fMuPrecFreq24*time));
}
return complexPol;
@@ -354,9 +368,6 @@ Double_t PSimulateMuTransition::GTSpinFlip(const Double_t &time)
* after ionization process is given by Px(t_i+1)*Px(t_i).
* 4) get the next electron capture time, continue until t_d is reached.
*
* <p> For isotropic muonium, TF:
* nu_12 and nu_34 with equal probabilities, probability for both states fMuFractionState12
* ni_23 and nu_14 with equal probabilities, probability for both states fMuFractionState23
*
* <p>Calculates Mu0 polarization in x direction during cyclic charge exchange.
* See M. Senba, J.Phys. B23, 1545 (1990), equations (9), (11)

20
src/tests/MuTransition/PSimulateMuTransition.h
View File

@@ -47,13 +47,16 @@ class PSimulateMuTransition : public TObject
virtual ~PSimulateMuTransition();
virtual void PrintSettings() const;
virtual void SetNmuons(Int_t value) { fNmuons = value; } //!< number of muons
virtual void SetDebugFlag(Bool_t value) { fDebugFlag = value; } //!< debug flag
virtual void SetBfield(Double_t value) { fBfield = value; } //!< sets magnetic field (T)
virtual void SetNmuons(Int_t value) { fNmuons = value; } //!< number of muons
virtual void SetNshowProgress(Int_t value) { fNshowProgress = value; } //!< frequency of output on screen how many muons have been processed
virtual void SetDebugFlag(Bool_t value) { fDebugFlag = value; } //!< debug flag
virtual void SetBfield(Double_t value) { fBfield = value; } //!< sets magnetic field (T)
virtual void SetMuPrecFreq12(Double_t value) { fMuPrecFreq12 = value; } //!< sets Mu transition frequency (MHz)
virtual void SetMuPrecFreq34(Double_t value) { fMuPrecFreq34 = value; } //!< sets Mu transition frequency (MHz)
virtual void SetMuPrecFreq23(Double_t value) { fMuPrecFreq23 = value; } //!< sets Mu transition frequency (MHz)
virtual void SetMuPrecFreq14(Double_t value) { fMuPrecFreq14 = value; } //!< sets Mu transition frequency (MHz)
virtual void SetMuPrecFreq13(Double_t value) { fMuPrecFreq13 = value; } //!< sets Mu transition frequency (MHz)
virtual void SetMuPrecFreq24(Double_t value) { fMuPrecFreq24 = value; } //!< sets Mu transition frequency (MHz)
virtual void SetCaptureRate(Double_t value){ fCaptureRate = value; } //!< sets Mu+ electron capture rate (MHz)
virtual void SetIonizationRate(Double_t value){ fIonizationRate = value; } //!< sets Mu0 ionization rate (MHz)
virtual void SetSpinFlipRate(Double_t value){ fSpinFlipRate = value; } //!< sets Mu0 spin flip rate (MHz)
@@ -63,6 +66,8 @@ class PSimulateMuTransition : public TObject
virtual void SetMuFractionState34(Double_t value){ fMuFractionState34 = value; }
virtual void SetMuFractionState23(Double_t value){ fMuFractionState23 = value; }
virtual void SetMuFractionState14(Double_t value){ fMuFractionState14 = value; }
virtual void SetMuFractionState13(Double_t value){ fMuFractionState13 = value; }
virtual void SetMuFractionState24(Double_t value){ fMuFractionState24 = value; }
virtual Bool_t IsValid() { return fValid; }
virtual void SetSeed(UInt_t seed);
@@ -81,6 +86,8 @@ class PSimulateMuTransition : public TObject
Double_t fMuPrecFreq34; //!< Mu transition frequency 34 (MHz)
Double_t fMuPrecFreq23; //!< Mu transition frequency 23 (MHz)
Double_t fMuPrecFreq14; //!< Mu transition frequency 14 (MHz)
Double_t fMuPrecFreq13; //!< Mu transition frequency 13 (MHz)
Double_t fMuPrecFreq24; //!< Mu transition frequency 24 (MHz)
Double_t fMuonPrecFreq; //!< muon precession frequency (MHz)
Double_t fCaptureRate; //!< Mu+ electron capture rate (MHz)
Double_t fIonizationRate; //!< Mu0 ionization rate (MHz)
@@ -94,8 +101,11 @@ class PSimulateMuTransition : public TObject
Double_t fMuFractionState34; //!< fraction of Mu in state 34
Double_t fMuFractionState23; //!< fraction of Mu in state 23
Double_t fMuFractionState14; //!< fraction of Mu in state 14
Int_t fNmuons; //!< number of muons to simulate
Bool_t fDebugFlag; //!< debug flag
Double_t fMuFractionState13; //!< fraction of Mu in state 13
Double_t fMuFractionState24; //!< fraction of Mu in state 24
Int_t fNmuons; //!< number of muons to simulate
Int_t fNshowProgress; //!< output on screen how many muons have been processed
Bool_t fDebugFlag; //!< debug flag
virtual Double_t NextEventTime(const Double_t &EventRate);
// virtual Double_t PrecessionPhase(const Double_t &time, const TString chargeState);

71
src/tests/MuTransition/runMuSimulation.C
View File

@@ -63,32 +63,47 @@ void runMuSimulation()
Double_t T = 300.; //temperature
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 capRate = 0.001;//*sqrt(T/200.); //assume that capture rate varies as sqrt(T), capRate = sigma*v*p , v ~ sqrt(T)
Double_t preFac = 6.7e7;
Double_t ionRate; //assume Arrhenius behaviour ionRate = preFac*exp(-EA/kT)
ionRate = 0.1; //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 = 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 = 0.487; //weight of transition 12
Double_t MuFrac34 = 0.487; //weight of transition 34
Double_t MuFrac23 = 0.013; //weight of transition 23
Double_t MuFrac14 = 0.013; //weight of transition 14
Int_t Nmuons = 5e6; //number of muons
Double_t Asym = 0.27; //muon decay asymmetry
Int_t debugFlag = 0; //print debug information on screen
ionRate = 0.001; //preFac * exp(-EA/(0.08625*T)); // Ge: 2.9*10^7MHz "attempt" frequency; 1K = 0.08625 meV
Double_t B = 100.0; //field in G
Double_t Bvar = 0.; //field variance
Double_t Freq12 = 40.023; //Mu freq of the 12 transition
Double_t Freq34 = 59.977; //Mu freq of the 34 transition
Double_t Freq23 = 238.549; //Mu freq of the 23 transition
Double_t Freq14 = 338.549; //Mu freq of the 14 transition
Double_t Freq13 = 278.571; //Mu freq of the 23 transition
Double_t Freq24 = 325.165; //Mu freq of the 14 transition
Double_t MuFrac = 1.0; //total Mu fraction
Double_t MuFrac12 = 0.486; //weight of transition 12
Double_t MuFrac34 = 0.486; //weight of transition 34
Double_t MuFrac23 = 0.014; //weight of transition 23
Double_t MuFrac14 = 0.014; //weight of transition 14
Double_t MuFrac13 = 0.0; //weight of transition 13
Double_t MuFrac24 = 0.0; //weight of transition 24
Int_t Nmuons = 5e6; //number of muons
Int_t NshowProgress = 1e4; //frequency to show progress on screen
Double_t Asym = 0.27; //muon decay asymmetry
Int_t debugFlag = 0; //print debug information on screen
TTimeStamp *timeStampStart = new TTimeStamp();
cout << endl << "Simulation started on:" << endl;
timeStampStart->Print("l);
cout << endl;
histogramFileName = TString("0");
histogramFileName += runNo;
histogramFileName += TString(".root");
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);
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, Freq23, MuFrac23, ionRate, capRate, spinFlipRate, B);
runTitle = TString("0");
runTitle += runNo;
runTitle += TString(titleStr);
cout << runTitle << endl << endl;
PSimulateMuTransition *simulateMuTransition = new PSimulateMuTransition();
if (!simulateMuTransition->IsValid()){
@@ -99,15 +114,20 @@ void runMuSimulation()
simulateMuTransition->SetMuPrecFreq34(Freq34); // MHz
simulateMuTransition->SetMuPrecFreq23(Freq23); // MHz
simulateMuTransition->SetMuPrecFreq14(Freq14); // MHz
simulateMuTransition->SetMuPrecFreq13(Freq13); // MHz
simulateMuTransition->SetMuPrecFreq24(Freq24); // MHz
simulateMuTransition->SetMuFraction(MuFrac); // initial Mu fraction
simulateMuTransition->SetMuFractionState12(MuFrac12);
simulateMuTransition->SetMuFractionState34(MuFrac34);
simulateMuTransition->SetMuFractionState23(MuFrac23);
simulateMuTransition->SetMuFractionState14(MuFrac14);
simulateMuTransition->SetMuFractionState14(MuFrac14);
simulateMuTransition->SetMuFractionState13(MuFrac13);
simulateMuTransition->SetMuFractionState24(MuFrac24);
simulateMuTransition->SetBfield(B/10000.); // Tesla
simulateMuTransition->SetCaptureRate(capRate); // MHz
simulateMuTransition->SetIonizationRate(ionRate); // MHz
simulateMuTransition->SetSpinFlipRate(spinFlipRate); // MHz
simulateMuTransition->SetNshowProgress(NshowProgress);
simulateMuTransition->SetNmuons(Nmuons);
simulateMuTransition->SetDecayAsymmetry(Asym);
simulateMuTransition->SetDebugFlag(debugFlag); // to print time and phase during charge-changing cycle
@@ -179,12 +199,18 @@ void runMuSimulation()
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 Precession frequency 13", Freq13);
header->Set("Simulation/Mu0 Precession frequency 24", Freq24);
header->Set("Simulation/Mu0 Fraction", MuFrac);
header->Set("Simulation/Mu0 Fraction 12", MuFrac12);
header->Set("Simulation/Mu0 Fraction 34", MuFrac34);
header->Set("Simulation/Mu0 Fraction 23", MuFrac23);
header->Set("Simulation/Mu0 Fraction 14", MuFrac14);
header->Set("Simulation/muon Capture Rate", capRate);
header->Set("Simulation/Mu0 Fraction 13", MuFrac13);
header->Set("Simulation/Mu0 Fraction 24", MuFrac24);
header->Set("Simulation/Mu0 Activation Energy", EA);
header->Set("Simulation/Mu0 Activation PreFactor", preFac);
header->Set("Simulation/Mux 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);
@@ -206,7 +232,7 @@ void runMuSimulation()
histo[i] = new TH1F(str, str, 18001, -0.5, 18000.5);
}
for (i=0; i<NDECAYHISTS; i++)
for (UInt_t i=0; i<NDECAYHISTS; i++)
decayAnaModule->Add(histo[i]);
// run simulation
@@ -230,6 +256,11 @@ void runMuSimulation()
gRunHeader->Write();
fout->Close();
cout << "Histograms written to " << histogramFileName.Data() << endl;
cout << endl << "Simulation stopped on:" << endl;
TTimeStamp *timeStampEnd = new TTimeStamp();
timeStampEnd->Print("l");
cout << endl;
// delete fout;
// delete header;