/*************************************************************************** musrfit.cpp Author: Andreas Suter e-mail: andreas.suter@psi.ch ***************************************************************************/ /*************************************************************************** * Copyright (C) 2007-2023 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 #ifdef HAVE_GOMP #include #endif #include #include //#include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #ifdef HAVE_GIT_REV_H #include "git-revision.h" #endif #include "PMusr.h" #include "PStartupHandler.h" #include "PMsrHandler.h" #include "PRunDataHandler.h" #include "PRunListCollection.h" #include "PFitter.h" //-------------------------------------------------------------------------- static int timeout = 3600; // given in (sec) //-------------------------------------------------------------------------- /** *

Sometimes musrfit is not terminating properly for reasons still not pinned down, hence * this workaround. It should be removed asap. Especially since it is leading to memory leaks. */ void* musrfit_timeout(void *args) { pid_t *pid = (pid_t*)(args); sleep(timeout); std::cerr << std::endl << ">> **FATAL ERROR** musrfit_timeout for task pid=" << *pid << " called! Will kill it!" << std::endl << std::endl; kill(*pid, SIGKILL); return (void*)nullptr; } //-------------------------------------------------------------------------- /** *

Sends the usage description to the standard output. */ void musrfit_syntax() { std::cout << std::endl << "usage: musrfit [ [-k, --keep-mn2-ouput] [-c, --chisq-only] [-t, --title-from-data-file]"; std::cout << std::endl << " [-e, --estimateN0] [-p, --per-run-block-chisq]"; std::cout << std::endl << " [--dump ] [--timeout ] |"; std::cout << std::endl << " -n, --no-of-cores-avail | -u, --use-no-of-threads |"; std::cout << std::endl << " --nexus-support | --show-dynamic-path | --version | --help"; std::cout << std::endl << " : msr input file"; std::cout << std::endl << " 'musrfit ' will execute musrfit"; std::cout << std::endl << " 'musrfit' or 'musrfit --help' will show this help"; std::cout << std::endl << " 'musrfit --version' will print the musrfit version"; std::cout << std::endl << " 'musrfit --nexus-support' will print if NeXus support is available."; std::cout << std::endl << " 'musrfit --show-dynamic-path' will print the internal dynamic library search paths."; std::cout << std::endl << " -k, --keep-mn2-output: will rename the files MINUIT2.OUTPUT and "; std::cout << std::endl << " MINUIT2.root to -mn2.output and -mn2.root,"; std::cout << std::endl << " respectively,"; std::cout << std::endl << " e.g. = 147.msr -> 147-mn2.output, 147-mn2.root"; std::cout << std::endl << " -c, --chisq-only: instead of fitting the data, chisq is just calculated"; std::cout << std::endl << " once and the result is set to the stdout. This feature is useful"; std::cout << std::endl << " to adjust initial parameters."; std::cout << std::endl << " -t, --title-from-data-file: will replace the run title by the"; std::cout << std::endl << " run title of the FIRST run of the run block, if a run title"; std::cout << std::endl << " is present in the data file."; std::cout << std::endl << " -e, --estimateN0: estimate N0 for single histogram fits."; std::cout << std::endl << " -p, --per-run-block-chisq: will write per run block chisq to the msr-file."; std::cout << std::endl << " -n, --no-of-cores-avail: print out how many cores are available (only vaild for OpenMP)"; std::cout << std::endl << " -u, --use-no-of-threads :"; std::cout << std::endl << " : number of threads to be used (OpenMP). Needs to be <= max. number of cores."; std::cout << std::endl << " If OpenMP is enable, the maximal number of cores is used, if it is not limited by this option."; std::cout << std::endl << " --dump is writing a data file with the fit data and the theory"; std::cout << std::endl << " can be 'ascii', 'root'"; std::cout << std::endl << " --timeout : overwrites to predefined timeout of " << timeout << " (sec)."; std::cout << std::endl << " <= 0 means timeout facility is not enabled. = nn"; std::cout << std::endl << " will set the timeout to nn (sec)."; std::cout << std::endl; std::cout << std::endl << " At the end of a fit, musrfit writes the fit results into an and"; std::cout << std::endl << " swaps them, i.e. in the you will find the fit results and in the"; std::cout << std::endl << " your initial guess values."; std::cout << std::endl << std::endl; } //-------------------------------------------------------------------------- /** *

Writes the fitted data- and theory-set in ascii format to disc. * * \param fln output file name * \param data processed-data handler * \param runCounter msr-file run counter needed to form the output file name */ void musrfit_write_ascii(TString fln, PRunData *data, int runCounter) { // generate dump file name TString count("_"); count += runCounter; Ssiz_t index = fln.Last('.'); fln.Insert(index, count); std::ofstream f; // open dump-file f.open(fln.Data(), std::iostream::out); if (!f.is_open()) { std::cout << std::endl << "Couldn't open dump (" << fln.Data() << ") file for writting, sorry ..."; std::cout << std::endl; return; } // dump data f << "% number of data values = " << data->GetValue()->size() << std::endl; f << "% time (us), value, error, theory" << std::endl; double time; for (unsigned int i=0; iGetValue()->size(); i++) { time = data->GetDataTimeStart() + (double)i*data->GetDataTimeStep(); f << time << ", " << data->GetValue()->at(i) << ", " << data->GetError()->at(i) << ", " << data->GetTheory()->at(i) << std::endl; } // close file f.close(); } //-------------------------------------------------------------------------- /** *

Go through all msr-file runs and write each of them in ascii-format to disc * (used for diagnostics). Be aware that this output is different to what * you will get when using musrview! * * \param fileName file name * \param runList run list collection handler */ void musrfit_dump_ascii(char *fileName, PRunListCollection *runList) { TString fln(fileName); fln.ReplaceAll(".msr", ".dat"); // go through the run list, get the data and dump it in a file int runCounter = 0; PRunData *data; // single histos unsigned int size = runList->GetNoOfSingleHisto(); if (size > 0) { for (unsigned int i=0; iGetSingleHisto(i); if (data) { // dump data musrfit_write_ascii(fln, data, runCounter); runCounter++; } } } // single histos size = runList->GetNoOfSingleHistoRRF(); if (size > 0) { for (unsigned int i=0; iGetSingleHistoRRF(i); if (data) { // dump data musrfit_write_ascii(fln, data, runCounter); runCounter++; } } } // asymmetry size = runList->GetNoOfAsymmetry(); if (size > 0) { for (unsigned int i=0; iGetAsymmetry(i); if (data) { // dump data musrfit_write_ascii(fln, data, runCounter); runCounter++; } } } // asymmetry RRF size = runList->GetNoOfAsymmetryRRF(); if (size > 0) { for (unsigned int i=0; iGetAsymmetryRRF(i); if (data) { // dump data musrfit_write_ascii(fln, data, runCounter); runCounter++; } } } // muMinus size = runList->GetNoOfMuMinus(); if (size > 0) { for (unsigned int i=0; iGetMuMinus(i); if (data) { // dump data musrfit_write_ascii(fln, data, runCounter); runCounter++; } } } // nonMusr size = runList->GetNoOfNonMusr(); if (size > 0) { for (unsigned int i=0; iGetNonMusr(i); if (data) { // dump data musrfit_write_ascii(fln, data, runCounter); runCounter++; } } } } //-------------------------------------------------------------------------- /** *

Writes the fitted data- and theory-set in root format to disc. * * \param f root file object * \param fln file name * \param data processed-data handler * \param runCounter msr-file run counter needed to form the output file name */ void musrfit_write_root(TFile &f, TString fln, PRunData *data, int runCounter) { char name[128]; TString title = fln.Copy(); snprintf(name, sizeof(name), "_%d", runCounter); title.ReplaceAll(".root", name); snprintf(name, sizeof(name),"c%d", runCounter); std::unique_ptr c = std::make_unique(name, title.Data(), 10, 10, 800, 600); // create histos Double_t diff = data->GetDataTimeStep(); Double_t start = -diff/2.0; Double_t end = data->GetDataTimeStep()*data->GetValue()->size(); std::unique_ptr hdata = std::make_unique("hdata", "run data", data->GetValue()->size(), start, end); std::unique_ptr htheo = std::make_unique("htheo", "run theory", data->GetValue()->size(), start, end); // fill data for (unsigned int i=0; iGetValue()->size(); i++) { hdata->SetBinContent(i+1, data->GetValue()->at(i)); hdata->SetBinError(i+1, data->GetError()->at(i)); htheo->SetBinContent(i+1, data->GetTheory()->at(i)); } hdata->SetMarkerStyle(20); hdata->Draw("*H HIST E1"); htheo->SetLineColor(2); htheo->SetLineWidth(3); htheo->Draw("C SAME"); f.WriteTObject(c.get()); } //-------------------------------------------------------------------------- /** *

Go through all msr-file runs and write each of them in root format to disc * (used for diagnostics). Be aware that this output is different to what * you will get when using musrview! * * \param fileName file name * \param runList run list connection handler */ void musrfit_dump_root(char *fileName, PRunListCollection *runList) { TString fln(fileName); fln.ReplaceAll(".msr", ".root"); TFile f(fln.Data(), "recreate"); // go through the run list, get the data and dump it in a file int runCounter = 0; PRunData *data; // single histos unsigned int size = runList->GetNoOfSingleHisto(); if (size > 0) { for (unsigned int i=0; iGetSingleHisto(i); if (data) { // dump data musrfit_write_root(f, fln, data, runCounter); runCounter++; } } } // single histo RRF size = runList->GetNoOfSingleHistoRRF(); if (size > 0) { for (unsigned int i=0; iGetSingleHistoRRF(i); if (data) { // dump data musrfit_write_root(f, fln, data, runCounter); runCounter++; } } } // asymmetry size = runList->GetNoOfAsymmetry(); if (size > 0) { for (unsigned int i=0; iGetAsymmetry(i); if (data) { // dump data musrfit_write_root(f, fln, data, runCounter); runCounter++; } } } // asymmetry RRF size = runList->GetNoOfAsymmetryRRF(); if (size > 0) { for (unsigned int i=0; iGetAsymmetryRRF(i); if (data) { // dump data musrfit_write_root(f, fln, data, runCounter); runCounter++; } } } // muMinus size = runList->GetNoOfMuMinus(); if (size > 0) { for (unsigned int i=0; iGetMuMinus(i); if (data) { // dump data musrfit_write_root(f, fln, data, runCounter); runCounter++; } } } // nonMusr size = runList->GetNoOfNonMusr(); if (size > 0) { for (unsigned int i=0; iGetNonMusr(i); if (data) { // dump data musrfit_write_root(f, fln, data, runCounter); runCounter++; } } } f.Close("R"); } //-------------------------------------------------------------------------- /** *

The musrfit program is used to fit muSR data. * For a detailed description/usage of the program, please see * \htmlonly musrfit online help * \endhtmlonly * \latexonly musrfit online help: \texttt{http://lmu.web.psi.ch/musrfit/user/html/user-manual.html#musrfit} * \endlatexonly * * return: * - PMUSR_SUCCESS if everthing went smooth * - PMUSR_WRONG_STARTUP_SYNTAX if syntax error is encountered * - line number if an error in the msr-file was encountered which cannot be handled. * * \param argc number of input arguments * \param argv list of input arguments */ int main(int argc, char *argv[]) { bool show_syntax = false; int status; bool keep_mn2_output = false; bool chisq_only = false; bool title_from_data_file = false; bool timeout_enabled = true; PStartupOptions startup_options; startup_options.writeExpectedChisq = false; startup_options.estimateN0 = false; int number_of_cores=1; #ifdef HAVE_GOMP number_of_cores = omp_get_num_procs(); #endif TString dump(""); char filename[1024]; // check syntax if (argc < 2) { musrfit_syntax(); return PMUSR_WRONG_STARTUP_SYNTAX; } // add default shared library path /usr/local/lib if not already persent const char *dsp = gSystem->GetDynamicPath(); if (strstr(dsp, "/usr/local/lib") == nullptr) gSystem->AddDynamicPath("/usr/local/lib"); if (argc == 2) { if (!strcmp(argv[1], "--version")) { #ifdef HAVE_CONFIG_H #ifdef HAVE_GIT_REV_H std::cout << std::endl << "musrfit version: " << PACKAGE_VERSION << ", git-branch: " << GIT_BRANCH << ", git-rev: " << GIT_CURRENT_SHA1 << " (" << BUILD_TYPE << "), ROOT version: " << ROOT_VERSION_USED << std::endl << std::endl; #else std::cout << std::endl << "musrfit version: " << PACKAGE_VERSION << " (" << BUILD_TYPE << "), ROOT version: " << ROOT_VERSION_USED << std::endl << std::endl; #endif #else #ifdef HAVE_GIT_REV_H std::cout << std::endl << "musrfit git-branch: " << GIT_BRANCH << ", git-rev: " << GIT_CURRENT_SHA1 << std::endl << std::endl; #else std::cout << std::endl << "musrfit version: unknown" << std::endl << std::endl; #endif #endif return PMUSR_SUCCESS; } else if (!strcmp(argv[1], "--nexus-support")) { #ifdef PNEXUS_ENABLED std::cout << std::endl << ">> musrfit: NeXus support enabled." << std::endl << std::endl; #else std::cout << std::endl << "musrfit: NeXus support NOT enabled." << std::endl << std::endl; #endif return PMUSR_SUCCESS; } else if (!strcmp(argv[1], "--show-dynamic-path")) { std::cout << std::endl << "musrfit: internal dynamic search paths for shared libraries/root dictionaries:"; std::cout << std::endl << " '" << gSystem->GetDynamicPath() << "'" << std::endl << std::endl; return PMUSR_SUCCESS; } else if (!strcmp(argv[1], "--help")) { show_syntax = true; } } if (show_syntax) { musrfit_syntax(); return PMUSR_WRONG_STARTUP_SYNTAX; } strcpy(filename, ""); for (int i=1; i> musrfit: **ERROR** found option --dump without " << std::endl; show_syntax = true; break; } } else if (!strcmp(argv[i], "-e") || !strcmp(argv[i], "--estimateN0")) { startup_options.estimateN0 = true; } else if (!strcmp(argv[i], "-p") || !strcmp(argv[i], "--per-run-block-chisq")) { startup_options.writeExpectedChisq = true; } else if (!strcmp(argv[i], "-n") || !strcmp(argv[i], "--no-of-cores-avail")) { #ifdef HAVE_GOMP std::cout << std::endl; std::cout << "musrfit: maxmimal number of cores for OpenMP available: " << omp_get_num_procs() << std::endl; std::cout << std::endl; #else std::cout << std::endl; std::cout << ">> musrfit: this option is only vaild if OpenMP is present. This seems not to be the case here. Sorry!" << std::endl; std::cout << std::endl; #endif return PMUSR_SUCCESS; } else if (!strcmp(argv[i], "-u") || !strcmp(argv[i], "--use-no-of-threads")) { if (i> musrfit: **ERROR** found option --use-no-of-threads with <= 0" << std::endl; std::cerr << " This doesn't make any sense." << std::endl; show_syntax = true; break; } else if (ival > number_of_cores) { #ifdef HAVE_GOMP std::cerr << std::endl << ">> musrfit: **WARNING** found option --use-no-of-threads with =" << ival << " > max available cores=" << number_of_cores << "." << std::endl; std::cerr << " Will set to max available cores." << std::endl; #else std::cerr << std::endl << ">> musrfit: **WARNING** option --use-no-of-threads can only be used if OpenMP is available." << std::endl; std::cerr << " Here it is not the case, and hence this option will be ignored." << std::endl; #endif } else { number_of_cores = ival; } } else { std::cerr << std::endl << ">> musrfit: **ERROR** found option --use-no-of-threads where it not a number: '" << argv[i+1] << "'" << std::endl; show_syntax = true; break; } i++; } else { std::cerr << std::endl << ">> musrfit: **ERROR** found option --use-no-of-threads without " << std::endl; show_syntax = true; break; } } else if (!strcmp(argv[i], "--timeout")) { if (i> musrfit: timeout disabled." << std::endl; } } else { std::cerr << std::endl << ">> musrfit: **ERROR** found option --timeout with unsupported = " << argv[i+1] << std::endl; show_syntax = true; break; } i++; } else { std::cerr << std::endl << ">> musrfit: **ERROR** found option --timeout without " << std::endl; show_syntax = true; break; } } else { show_syntax = true; break; } } // check if a filename is present if (strlen(filename) == 0) { show_syntax = true; std::cout << std::endl << ">> musrfit **ERROR** no msr-file present!" << std::endl; } if (show_syntax) { musrfit_syntax(); return PMUSR_WRONG_STARTUP_SYNTAX; } // check if dump string does make sense if (!dump.IsNull()) { dump.ToLower(); if (!dump.Contains("ascii") && !dump.Contains("root")) { std::cerr << std::endl << ">> musrfit: **ERROR** found option --dump with unsupported = " << dump << std::endl; musrfit_syntax(); return PMUSR_WRONG_STARTUP_SYNTAX; } } // read startup file char startup_path_name[128]; std::unique_ptr saxParser = std::make_unique(); std::unique_ptr startupHandler = std::make_unique(); if (!startupHandler->StartupFileFound()) { std::cerr << std::endl << ">> musrfit **WARNING** couldn't find " << startupHandler->GetStartupFilePath().Data(); std::cerr << std::endl; } else { strcpy(startup_path_name, startupHandler->GetStartupFilePath().Data()); saxParser->ConnectToHandler("PStartupHandler", startupHandler.get()); //status = saxParser->ParseFile(startup_path_name); // parsing the file as above seems to lead to problems in certain environments; // use the parseXmlFile function instead (see PStartupHandler.cpp for the definition) status = parseXmlFile(saxParser.get(), startup_path_name); // check for parse errors if (status) { // error std::cerr << std::endl << ">> musrfit **WARNING** Reading/parsing musrfit_startup.xml failed."; std::cerr << std::endl; } } #ifdef HAVE_GOMP // set omp_set_num_threads omp_set_num_threads(number_of_cores); #endif // read msr-file std::unique_ptr msrHandler; if (startupHandler) msrHandler = std::make_unique(filename, &startup_options); else msrHandler = std::make_unique(filename); status = msrHandler->ReadMsrFile(); if (status != PMUSR_SUCCESS) { switch (status) { case PMUSR_MSR_FILE_NOT_FOUND: std::cout << std::endl << ">> musrfit **ERROR** couldn't find " << filename << std::endl << std::endl; break; case PMUSR_MSR_SYNTAX_ERROR: std::cout << std::endl << ">> musrfit **SYNTAX ERROR** in file " << filename << ", full stop here." << std::endl << std::endl; break; default: std::cout << std::endl << ">> musrfit **UNKOWN ERROR** when trying to read the msr-file" << std::endl << std::endl; break; } return status; } // read all the necessary runs (raw data) std::unique_ptr dataHandler; if (startupHandler) dataHandler = std::make_unique(msrHandler.get(), startupHandler->GetDataPathList()); else dataHandler = std::make_unique(msrHandler.get()); dataHandler->ReadData(); bool success = dataHandler->IsAllDataAvailable(); if (!success) { std::cout << std::endl << ">> musrfit **ERROR** Couldn't read all data files, will quit ..." << std::endl; } // if present, replace the run title of the with the run title of the FIRST run in the run block of the msr-file if (title_from_data_file && success) { PMsrRunList *rl = msrHandler->GetMsrRunList(); PRawRunData *rrd = dataHandler->GetRunData(*(rl->at(0).GetRunName())); if (rrd->GetRunTitle()->Length() > 0) msrHandler->SetMsrTitle(*rrd->GetRunTitle()); } // generate the necessary fit histogramms for the fit std::unique_ptr runListCollection; if (success) { // feed all the necessary histogramms for the fit runListCollection = std::make_unique(msrHandler.get(), dataHandler.get()); for (unsigned int i=0; i < msrHandler->GetMsrRunList()->size(); i++) { success = runListCollection->Add(i, kFit); if (!success) { std::cout << std::endl << ">> musrfit **ERROR** Couldn't handle run no " << i+1 << ": "; std::cout << (*msrHandler->GetMsrRunList())[i].GetRunName()->Data(); break; } } } // start timeout thread std::unique_ptr th; if (timeout_enabled) { pid_t musrfit_pid = getpid(); th = std::make_unique(musrfit_timeout, (void*)&musrfit_pid); if (th) { th->Run(); } } // do fitting std::unique_ptr fitter; if (success) { fitter = std::make_unique(msrHandler.get(), runListCollection.get(), chisq_only); if (fitter->IsValid()) { fitter->DoFit(); if (!fitter->IsScanOnly()) msrHandler->SetMsrStatisticConverged(fitter->HasConverged()); } } // write log file if (success && !chisq_only) { if (!fitter->IsScanOnly()) { status = msrHandler->WriteMsrLogFile(); if (status != PMUSR_SUCCESS) { switch (status) { case PMUSR_MSR_LOG_FILE_WRITE_ERROR: std::cout << std::endl << ">> musrfit **ERROR** couldn't write mlog-file" << std::endl << std::endl; break; case PMUSR_TOKENIZE_ERROR: std::cout << std::endl << ">> musrfit **ERROR** couldn't generate mlog-file name" << std::endl << std::endl; break; default: std::cout << std::endl << ">> musrfit **UNKOWN ERROR** when trying to write the mlog-file" << std::endl << std::endl; break; } } } } // check if dump is wanted if (success && !dump.IsNull()) { std::cout << std::endl << "will write dump file ..." << std::endl; dump.ToLower(); if (dump.Contains("ascii")) musrfit_dump_ascii(filename, runListCollection.get()); else if (dump.Contains("root")) musrfit_dump_root(filename, runListCollection.get()); else std::cout << std::endl << "do not know format " << dump.Data() << ", sorry :-| " << std::endl; } // rename MINUIT2.OUTPUT and MINUIT2.root file if wanted if (success) { if (keep_mn2_output && !chisq_only && !fitter->IsScanOnly()) { // 1st rename MINUIT2.OUTPUT TString fln = TString(filename); char ext[32]; strcpy(ext, "-mn2.output"); fln.ReplaceAll(".msr", 4, ext, strlen(ext)); gSystem->CopyFile("MINUIT2.OUTPUT", fln.Data(), kTRUE); // 2nd rename MINUIT2.ROOT fln = TString(filename); strcpy(ext, "-mn2.root"); fln.ReplaceAll(".msr", 4, ext, strlen(ext)); gSystem->CopyFile("MINUIT2.root", fln.Data(), kTRUE); } } if (success) { if (!chisq_only && !fitter->IsScanOnly()) { // swap msr- and mlog-file std::cout << std::endl << ">> swapping msr-, mlog-file ..." << std::endl; // copy msr-file -> __temp.msr gSystem->CopyFile(filename, "__temp.msr", kTRUE); // copy mlog-file -> msr-file TString fln = TString(filename); char ext[32]; strcpy(ext, ".mlog"); fln.ReplaceAll(".msr", 4, ext, strlen(ext)); gSystem->CopyFile(fln.Data(), filename, kTRUE); // copy __temp.msr -> mlog-file gSystem->CopyFile("__temp.msr", fln.Data(), kTRUE); // delete __temp.msr TSystemFile tmp("__temp.msr", "./"); tmp.Delete(); } } std::cout << std::endl << "done ..." << std::endl; return PMUSR_SUCCESS; } // end ---------------------------------------------------------------------