musrfit 1.10.0
musrfit.cpp
Go to the documentation of this file.
1/***************************************************************************
2
3 musrfit.cpp
4
5 Author: Andreas Suter
6 e-mail: andreas.suter@psi.ch
7
8***************************************************************************/
9
10/***************************************************************************
11 * Copyright (C) 2007-2026 by Andreas Suter *
12 * andreas.suter@psi.ch *
13 * *
14 * This program is free software; you can redistribute it and/or modify *
15 * it under the terms of the GNU General Public License as published by *
16 * the Free Software Foundation; either version 2 of the License, or *
17 * (at your option) any later version. *
18 * *
19 * This program is distributed in the hope that it will be useful, *
20 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
21 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
22 * GNU General Public License for more details. *
23 * *
24 * You should have received a copy of the GNU General Public License *
25 * along with this program; if not, write to the *
26 * Free Software Foundation, Inc., *
27 * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
28 ***************************************************************************/
29
30#ifdef HAVE_CONFIG_H
31#include "config.h"
32#endif
33
34#ifdef HAVE_GOMP
35#include <omp.h>
36#endif
37
38#include <cstdio>
39#include <cstdlib>
40#include <cstring>
41#include <sys/types.h>
42#include <unistd.h>
43#include <signal.h>
44
45#include <iostream>
46#include <fstream>
47#include <memory>
48#include <string>
49
50#include <TSAXParser.h>
51#include <TString.h>
52#include <TFile.h>
53#include <TCanvas.h>
54#include <TH1.h>
55#include <TSystem.h>
56#include <TSystemFile.h>
57#include <TThread.h>
58
59#ifdef HAVE_GIT_REV_H
60#include "git-revision.h"
61#endif
62
63#include "PMusr.h"
64#include "PStartupHandler.h"
65#include "PMsrHandler.h"
66#include "PRunDataHandler.h"
67#include "PRunListCollection.h"
68#include "PFitter.h"
69
70//--------------------------------------------------------------------------
71
72static int timeout = 3600; // given in (sec)
73
74//--------------------------------------------------------------------------
79void* musrfit_timeout(void *args)
80{
81 pid_t *pid = (pid_t*)(args);
82
83 sleep(timeout);
84
85 std::cerr << std::endl << ">> **FATAL ERROR** musrfit_timeout for task pid=" << *pid << " called! Will kill it!" << std::endl << std::endl;
86
87 kill(*pid, SIGKILL);
88
89 return (void*)nullptr;
90}
91
92//--------------------------------------------------------------------------
97{
98 std::cout << std::endl << "usage: musrfit [<msr-file> [-k, --keep-mn2-ouput] [-c, --chisq-only] [-t, --title-from-data-file]";
99 std::cout << std::endl << " [-e, --estimateN0] [-p, --per-run-block-chisq]";
100 std::cout << std::endl << " [--dump <type>] [--timeout <timeout_tag>] |";
101 std::cout << std::endl << " -n, --no-of-cores-avail | -u, --use-no-of-threads <number> |";
102 std::cout << std::endl << " --nexus-support | --show-dynamic-path | --version | --help";
103 std::cout << std::endl << " <msr-file>: msr input file";
104 std::cout << std::endl << " 'musrfit <msr-file>' will execute musrfit";
105 std::cout << std::endl << " 'musrfit' or 'musrfit --help' will show this help";
106 std::cout << std::endl << " 'musrfit --version' will print the musrfit version";
107 std::cout << std::endl << " 'musrfit --nexus-support' will print if NeXus support is available.";
108 std::cout << std::endl << " 'musrfit --show-dynamic-path' will print the internal dynamic library search paths.";
109 std::cout << std::endl << " -k, --keep-mn2-output: will rename the files MINUIT2.OUTPUT and ";
110 std::cout << std::endl << " MINUIT2.root to <msr-file>-mn2.output and <msr-file>-mn2.root,";
111 std::cout << std::endl << " respectively,";
112 std::cout << std::endl << " e.g. <msr-file> = 147.msr -> 147-mn2.output, 147-mn2.root";
113 std::cout << std::endl << " -c, --chisq-only: instead of fitting the data, chisq is just calculated";
114 std::cout << std::endl << " once and the result is set to the stdout. This feature is useful";
115 std::cout << std::endl << " to adjust initial parameters.";
116 std::cout << std::endl << " -t, --title-from-data-file: will replace the <msr-file> run title by the";
117 std::cout << std::endl << " run title of the FIRST run of the <msr-file> run block, if a run title";
118 std::cout << std::endl << " is present in the data file.";
119 std::cout << std::endl << " -e, --estimateN0: estimate N0 for single histogram fits.";
120 std::cout << std::endl << " -p, --per-run-block-chisq: will write per run block chisq to the msr-file.";
121 std::cout << std::endl << " -n, --no-of-cores-avail: print out how many cores are available (only vaild for OpenMP)";
122 std::cout << std::endl << " -u, --use-no-of-threads <number>:";
123 std::cout << std::endl << " <number>: number of threads to be used (OpenMP). Needs to be <= max. number of cores.";
124 std::cout << std::endl << " If OpenMP is enable, the maximal number of cores is used, if it is not limited by this option.";
125 std::cout << std::endl << " -r, --reset: reset startup musrfit_startup.xml, i.e. rewrite a default, and quit.";
126 std::cout << std::endl << " The order of which musrfit_startup.xml is reset is:";
127 std::cout << std::endl << " (i) if present in the current dir.";
128 std::cout << std::endl << " (ii) if present under $HOME/.musrfit/";
129 std::cout << std::endl << " (iii) if present under $MUSRFITPATH/";
130 std::cout << std::endl << " (iv) if present under $ROOTSYS/";
131 std::cout << std::endl << " -y, --yaml: write fit results (MINUIT2.OUTPUT) into a yaml-file. Output <msr-file>.yaml";
132 std::cout << std::endl << " --dump <type> is writing a data file with the fit data and the theory";
133 std::cout << std::endl << " <type> can be 'ascii', 'root'";
134 std::cout << std::endl << " --timeout <timeout_tag>: overwrites to predefined timeout of " << timeout << " (sec).";
135 std::cout << std::endl << " <timeout_tag> <= 0 means timeout facility is not enabled. <timeout_tag> = nn";
136 std::cout << std::endl << " will set the timeout to nn (sec).";
137 std::cout << std::endl;
138 std::cout << std::endl << " At the end of a fit, musrfit writes the fit results into an <mlog-file> and";
139 std::cout << std::endl << " swaps them, i.e. in the <msr-file> you will find the fit results and in the";
140 std::cout << std::endl << " <mlog-file> your initial guess values.";
141 std::cout << std::endl << std::endl;
142}
143
144//--------------------------------------------------------------------------
152void musrfit_write_ascii(TString fln, PRunData *data, int runCounter)
153{
154 // generate dump file name
155 TString count("_");
156 count += runCounter;
157 Ssiz_t index = fln.Last('.');
158 fln.Insert(index, count);
159
160 std::ofstream f;
161
162 // open dump-file
163 f.open(fln.Data(), std::iostream::out);
164 if (!f.is_open()) {
165 std::cout << std::endl << "Couldn't open dump (" << fln.Data() << ") file for writting, sorry ...";
166 std::cout << std::endl;
167 return;
168 }
169
170 // dump data
171 f << "% number of data values = " << data->GetValue()->size() << std::endl;
172 f << "% time (us), value, error, theory" << std::endl;
173 double time;
174 for (unsigned int i=0; i<data->GetValue()->size(); i++) {
175 time = data->GetDataTimeStart() + (double)i*data->GetDataTimeStep();
176 f << time << ", " << data->GetValue()->at(i) << ", " << data->GetError()->at(i) << ", " << data->GetTheory()->at(i) << std::endl;
177 }
178
179 // close file
180 f.close();
181}
182
183//--------------------------------------------------------------------------
192void musrfit_dump_ascii(const char *fileName, PRunListCollection *runList)
193{
194 TString fln(fileName);
195 fln.ReplaceAll(".msr", ".dat");
196
197 // go through the run list, get the data and dump it in a file
198
199 int runCounter = 0;
200 PRunData *data;
201
202 // single histos
203 unsigned int size = runList->GetNoOfSingleHisto();
204 if (size > 0) {
205 for (unsigned int i=0; i<size; i++) {
206 data = runList->GetSingleHisto(i);
207 if (data) {
208 // dump data
209 musrfit_write_ascii(fln, data, runCounter);
210 runCounter++;
211 }
212 }
213 }
214
215 // single histos
216 size = runList->GetNoOfSingleHistoRRF();
217 if (size > 0) {
218 for (unsigned int i=0; i<size; i++) {
219 data = runList->GetSingleHistoRRF(i);
220 if (data) {
221 // dump data
222 musrfit_write_ascii(fln, data, runCounter);
223 runCounter++;
224 }
225 }
226 }
227
228 // asymmetry
229 size = runList->GetNoOfAsymmetry();
230 if (size > 0) {
231 for (unsigned int i=0; i<size; i++) {
232 data = runList->GetAsymmetry(i);
233 if (data) {
234 // dump data
235 musrfit_write_ascii(fln, data, runCounter);
236 runCounter++;
237 }
238 }
239 }
240
241 // asymmetry RRF
242 size = runList->GetNoOfAsymmetryRRF();
243 if (size > 0) {
244 for (unsigned int i=0; i<size; i++) {
245 data = runList->GetAsymmetryRRF(i);
246 if (data) {
247 // dump data
248 musrfit_write_ascii(fln, data, runCounter);
249 runCounter++;
250 }
251 }
252 }
253
254 // muMinus
255 size = runList->GetNoOfMuMinus();
256 if (size > 0) {
257 for (unsigned int i=0; i<size; i++) {
258 data = runList->GetMuMinus(i);
259 if (data) {
260 // dump data
261 musrfit_write_ascii(fln, data, runCounter);
262 runCounter++;
263 }
264 }
265 }
266
267 // nonMusr
268 size = runList->GetNoOfNonMusr();
269 if (size > 0) {
270 for (unsigned int i=0; i<size; i++) {
271 data = runList->GetNonMusr(i);
272 if (data) {
273 // dump data
274 musrfit_write_ascii(fln, data, runCounter);
275 runCounter++;
276 }
277 }
278 }
279}
280
281//--------------------------------------------------------------------------
290void musrfit_write_root(TFile &f, TString fln, PRunData *data, int runCounter)
291{
292 char name[128];
293
294 TString title = fln.Copy();
295 snprintf(name, sizeof(name), "_%d", runCounter);
296 title.ReplaceAll(".root", name);
297
298 snprintf(name, sizeof(name),"c%d", runCounter);
299
300 std::unique_ptr<TCanvas> c = std::make_unique<TCanvas>(name, title.Data(), 10, 10, 800, 600);
301
302 // create histos
303 Double_t diff = data->GetDataTimeStep();
304 Double_t start = -diff/2.0;
305 Double_t end = data->GetDataTimeStep()*data->GetValue()->size();
306 std::unique_ptr<TH1F> hdata = std::make_unique<TH1F>("hdata", "run data", data->GetValue()->size(), start, end);
307 std::unique_ptr<TH1F> htheo = std::make_unique<TH1F>("htheo", "run theory", data->GetValue()->size(), start, end);
308
309 // fill data
310 for (unsigned int i=0; i<data->GetValue()->size(); i++) {
311 hdata->SetBinContent(i+1, data->GetValue()->at(i));
312 hdata->SetBinError(i+1, data->GetError()->at(i));
313 htheo->SetBinContent(i+1, data->GetTheory()->at(i));
314 }
315
316 hdata->SetMarkerStyle(20);
317 hdata->Draw("*H HIST E1");
318
319 htheo->SetLineColor(2);
320 htheo->SetLineWidth(3);
321 htheo->Draw("C SAME");
322
323 f.WriteTObject(c.get());
324}
325
326//--------------------------------------------------------------------------
335void musrfit_dump_root(const char *fileName, PRunListCollection *runList)
336{
337 TString fln(fileName);
338 fln.ReplaceAll(".msr", ".root");
339
340 TFile f(fln.Data(), "recreate");
341
342 // go through the run list, get the data and dump it in a file
343
344 int runCounter = 0;
345 PRunData *data;
346
347 // single histos
348 unsigned int size = runList->GetNoOfSingleHisto();
349 if (size > 0) {
350 for (unsigned int i=0; i<size; i++) {
351 data = runList->GetSingleHisto(i);
352 if (data) {
353 // dump data
354 musrfit_write_root(f, fln, data, runCounter);
355 runCounter++;
356 }
357 }
358 }
359
360 // single histo RRF
361 size = runList->GetNoOfSingleHistoRRF();
362 if (size > 0) {
363 for (unsigned int i=0; i<size; i++) {
364 data = runList->GetSingleHistoRRF(i);
365 if (data) {
366 // dump data
367 musrfit_write_root(f, fln, data, runCounter);
368 runCounter++;
369 }
370 }
371 }
372
373 // asymmetry
374 size = runList->GetNoOfAsymmetry();
375 if (size > 0) {
376 for (unsigned int i=0; i<size; i++) {
377 data = runList->GetAsymmetry(i);
378 if (data) {
379 // dump data
380 musrfit_write_root(f, fln, data, runCounter);
381 runCounter++;
382 }
383 }
384 }
385
386 // asymmetry RRF
387 size = runList->GetNoOfAsymmetryRRF();
388 if (size > 0) {
389 for (unsigned int i=0; i<size; i++) {
390 data = runList->GetAsymmetryRRF(i);
391 if (data) {
392 // dump data
393 musrfit_write_root(f, fln, data, runCounter);
394 runCounter++;
395 }
396 }
397 }
398
399 // muMinus
400 size = runList->GetNoOfMuMinus();
401 if (size > 0) {
402 for (unsigned int i=0; i<size; i++) {
403 data = runList->GetMuMinus(i);
404 if (data) {
405 // dump data
406 musrfit_write_root(f, fln, data, runCounter);
407 runCounter++;
408 }
409 }
410 }
411
412 // nonMusr
413 size = runList->GetNoOfNonMusr();
414 if (size > 0) {
415 for (unsigned int i=0; i<size; i++) {
416 data = runList->GetNonMusr(i);
417 if (data) {
418 // dump data
419 musrfit_write_root(f, fln, data, runCounter);
420 runCounter++;
421 }
422 }
423 }
424
425 f.Close("R");
426}
427
428//--------------------------------------------------------------------------
445int main(int argc, char *argv[])
446{
447
448 bool show_syntax = false;
449
450 // check syntax (root unrelated)
451 if (argc < 2) {
454 }
455
456 if (argc == 2) {
457 if (!strcmp(argv[1], "--version")) {
458#ifdef HAVE_CONFIG_H
459#ifdef HAVE_GIT_REV_H
460 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;
461#else
462 std::cout << std::endl << "musrfit version: " << PACKAGE_VERSION << " (" << BUILD_TYPE << "), ROOT version: " << ROOT_VERSION_USED << std::endl << std::endl;
463#endif
464#else
465#ifdef HAVE_GIT_REV_H
466 std::cout << std::endl << "musrfit git-branch: " << GIT_BRANCH << ", git-rev: " << GIT_CURRENT_SHA1 << std::endl << std::endl;
467#else
468 std::cout << std::endl << "musrfit version: unknown" << std::endl << std::endl;
469#endif
470#endif
471 return PMUSR_SUCCESS;
472 } else if (!strcmp(argv[1], "--nexus-support")) {
473#ifdef PNEXUS_ENABLED
474 std::cout << std::endl << ">> musrfit: NeXus support enabled. ";
475 std::cout << " HDF5 ready. ";
476#ifdef HAVE_HDF4
477 std::cout << " HDF4 ready." << std::endl << std::endl;
478#else
479 std::cout << " HDF4 not ready." << std::endl << std::endl;
480#endif
481#else
482 std::cout << std::endl << "musrfit: NeXus support NOT enabled." << std::endl << std::endl;
483#endif
484 return PMUSR_SUCCESS;
485 } else if (!strcmp(argv[1], "--help")) {
486 show_syntax = true;
487 }
488 }
489
490 if (show_syntax) {
493 }
494
495 int status;
496 bool keep_mn2_output = false;
497 bool chisq_only = false;
498 bool yaml_out = false;
499 bool title_from_data_file = false;
500 bool timeout_enabled = true;
501 bool reset_startup_file = false;
502 PStartupOptions startup_options;
503 startup_options.writeExpectedChisq = false;
504 startup_options.estimateN0 = false;
505 int number_of_cores=1;
506
507#ifdef HAVE_GOMP
508 number_of_cores = omp_get_num_procs();
509#endif
510
511 TString dump("");
512 std::string filename{};
513
514 // add default shared library path /usr/local/lib if not already persent
515 const char *dsp = gSystem->GetDynamicPath();
516 if (strstr(dsp, "/usr/local/lib") == nullptr)
517 gSystem->AddDynamicPath("/usr/local/lib");
518
519 // check syntax (root related)
520 if (argc == 2) {
521 if (!strcmp(argv[1], "--show-dynamic-path")) {
522 std::cout << std::endl << "musrfit: internal dynamic search paths for shared libraries/root dictionaries:";
523 std::cout << std::endl << " '" << gSystem->GetDynamicPath() << "'" << std::endl << std::endl;
524 return PMUSR_SUCCESS;
525 }
526 }
527
528 for (int i=1; i<argc; i++) {
529 if (strstr(argv[i], ".msr")) {
530 filename = argv[i];
531 } else if (!strcmp(argv[i], "-k") || !strcmp(argv[i], "--keep-mn2-output")) {
532 keep_mn2_output = true;
533 } else if (!strcmp(argv[i], "-c") || !strcmp(argv[i], "--chisq-only")) {
534 chisq_only = true;
535 } else if (!strcmp(argv[i], "-t") || !strcmp(argv[i], "--title-from-data-file")) {
536 title_from_data_file = true;
537 } else if (!strcmp(argv[i], "--dump")) {
538 if (i<argc-1) {
539 dump = TString(argv[i+1]);
540 i++;
541 } else {
542 std::cerr << std::endl << ">> musrfit: **ERROR** found option --dump without <type>" << std::endl;
543 show_syntax = true;
544 break;
545 }
546 } else if (!strcmp(argv[i], "-e") || !strcmp(argv[i], "--estimateN0")) {
547 startup_options.estimateN0 = true;
548 } else if (!strcmp(argv[i], "-p") || !strcmp(argv[i], "--per-run-block-chisq")) {
549 startup_options.writeExpectedChisq = true;
550 } else if (!strcmp(argv[i], "-y") || !strcmp(argv[i], "--yaml")) {
551 yaml_out = true;
552 } else if (!strcmp(argv[i], "-n") || !strcmp(argv[i], "--no-of-cores-avail")) {
553#ifdef HAVE_GOMP
554 std::cout << std::endl;
555 std::cout << "musrfit: maxmimal number of cores for OpenMP available: " << omp_get_num_procs() << std::endl;
556 std::cout << std::endl;
557#else
558 std::cout << std::endl;
559 std::cout << ">> musrfit: this option is only vaild if OpenMP is present. This seems not to be the case here. Sorry!" << std::endl;
560 std::cout << std::endl;
561#endif
562 return PMUSR_SUCCESS;
563 } else if (!strcmp(argv[i], "-u") || !strcmp(argv[i], "--use-no-of-threads")) {
564 if (i<argc-1) {
565 TString str(argv[i+1]);
566 if (str.IsFloat()) {
567 int ival = str.Atoi();
568 if (ival <= 0) {
569 std::cerr << std::endl << ">> musrfit: **ERROR** found option --use-no-of-threads with <number> <= 0" << std::endl;
570 std::cerr << " This doesn't make any sense." << std::endl;
571 show_syntax = true;
572 break;
573 } else if (ival > number_of_cores) {
574#ifdef HAVE_GOMP
575 std::cerr << std::endl << ">> musrfit: **WARNING** found option --use-no-of-threads with <number>=" << ival << " > max available cores=" << number_of_cores << "." << std::endl;
576 std::cerr << " Will set <number> to max available cores." << std::endl;
577#else
578 std::cerr << std::endl << ">> musrfit: **WARNING** option --use-no-of-threads can only be used if OpenMP is available." << std::endl;
579 std::cerr << " Here it is not the case, and hence this option will be ignored." << std::endl;
580#endif
581 } else {
582 number_of_cores = ival;
583 }
584 } else {
585 std::cerr << std::endl << ">> musrfit: **ERROR** found option --use-no-of-threads where <number> it not a number: '" << argv[i+1] << "'" << std::endl;
586 show_syntax = true;
587 break;
588 }
589 i++;
590 } else {
591 std::cerr << std::endl << ">> musrfit: **ERROR** found option --use-no-of-threads without <number>" << std::endl;
592 show_syntax = true;
593 break;
594 }
595 } else if (!strcmp(argv[i], "-r") || !strcmp(argv[i], "--reset")) {
596 reset_startup_file = true;
597 } else if (!strcmp(argv[i], "--timeout")) {
598 if (i<argc-1) {
599 TString str(argv[i+1]);
600 if (str.IsFloat()) {
601 timeout = str.Atoi();
602 if (timeout <= 0) {
603 timeout_enabled = false;
604 std::cout << std::endl << ">> musrfit: timeout disabled." << std::endl;
605 }
606 } else {
607 std::cerr << std::endl << ">> musrfit: **ERROR** found option --timeout with unsupported <timeout_tag> = " << argv[i+1] << std::endl;
608 show_syntax = true;
609 break;
610 }
611 i++;
612 } else {
613 std::cerr << std::endl << ">> musrfit: **ERROR** found option --timeout without <timeout_tag>" << std::endl;
614 show_syntax = true;
615 break;
616 }
617 } else {
618 show_syntax = true;
619 break;
620 }
621 }
622
623 // check if a filename is present
624 if (filename.empty() && !reset_startup_file) {
625 show_syntax = true;
626 std::cout << std::endl << ">> musrfit **ERROR** no msr-file present!" << std::endl;
627 }
628
629 if (show_syntax) {
632 }
633
634 // check if dump string does make sense
635 if (!dump.IsNull()) {
636 dump.ToLower();
637 if (!dump.Contains("ascii") && !dump.Contains("root")) {
638 std::cerr << std::endl << ">> musrfit: **ERROR** found option --dump with unsupported <type> = " << dump << std::endl;
641 }
642 }
643
644 // read startup file
645 std::string startup_path_name{};
646 std::unique_ptr<TSAXParser> saxParser = std::make_unique<TSAXParser>();
647 std::unique_ptr<PStartupHandler> startupHandler = std::make_unique<PStartupHandler>(reset_startup_file);
648 if (reset_startup_file) // only rewrite musrfit_startup.xml has been requested
649 return PMUSR_SUCCESS;
650 if (!startupHandler->StartupFileFound()) {
651 std::cerr << std::endl << ">> musrfit **WARNING** couldn't find " << startupHandler->GetStartupFilePath().Data();
652 std::cerr << std::endl;
653 } else {
654 startup_path_name = startupHandler->GetStartupFilePath().Data();
655 saxParser->ConnectToHandler("PStartupHandler", startupHandler.get());
656 //status = saxParser->ParseFile(startup_path_name);
657 // parsing the file as above seems to lead to problems in certain environments;
658 // use the parseXmlFile function instead (see PStartupHandler.cpp for the definition)
659 status = parseXmlFile(saxParser.get(), startup_path_name.c_str());
660 // check for parse errors
661 if (status) { // error
662 std::cerr << std::endl << ">> musrfit **WARNING** Reading/parsing musrfit_startup.xml failed.";
663 std::cerr << std::endl;
664 }
665 }
666
667#ifdef HAVE_GOMP
668 // set omp_set_num_threads
669 omp_set_num_threads(number_of_cores);
670#endif
671
672 // read msr-file
673 std::unique_ptr<PMsrHandler> msrHandler;
674 if (startupHandler)
675 msrHandler = std::make_unique<PMsrHandler>(filename.c_str(), &startup_options);
676 else
677 msrHandler = std::make_unique<PMsrHandler>(filename.c_str());
678 status = msrHandler->ReadMsrFile();
679 if (status != PMUSR_SUCCESS) {
680 switch (status) {
682 std::cout << std::endl << ">> musrfit **ERROR** couldn't find " << filename << std::endl << std::endl;
683 break;
685 std::cout << std::endl << ">> musrfit **SYNTAX ERROR** in file " << filename << ", full stop here." << std::endl << std::endl;
686 break;
687 default:
688 std::cout << std::endl << ">> musrfit **UNKOWN ERROR** when trying to read the msr-file" << std::endl << std::endl;
689 break;
690 }
691 return status;
692 }
693
694 // read all the necessary runs (raw data)
695 std::unique_ptr<PRunDataHandler> dataHandler;
696 if (startupHandler)
697 dataHandler = std::make_unique<PRunDataHandler>(msrHandler.get(), startupHandler->GetDataPathList());
698 else
699 dataHandler = std::make_unique<PRunDataHandler>(msrHandler.get());
700
701 dataHandler->ReadData();
702
703 bool success = dataHandler->IsAllDataAvailable();
704 if (!success) {
705 std::cerr << std::endl << ">> musrfit **ERROR** Couldn't read all data files, will quit ..." << std::endl;
706 }
707
708 // if present, replace the run title of the <msr-file> with the run title of the FIRST run in the run block of the msr-file
709 if (title_from_data_file && success) {
710 PMsrRunList *rl = msrHandler->GetMsrRunList();
711 if (rl->empty()) {
712 success = false;
713 std::cerr << std::endl << ">> musrfit **ERROR** no run list present." << std::endl;
714 }
715 if (success) {
716 TString *name = rl->at(0).GetRunName();
717 if (name == nullptr) {
718 std::cerr << std::endl << ">> musrfit **ERROR** to obtain run list name." << std::endl;
719 success = false;
720 }
721 PRawRunData *rrd = nullptr;
722 if (success) {
723 rrd = dataHandler->GetRunData(*(rl->at(0).GetRunName()));
724 if (rrd == nullptr) {
725 std::cerr << std::endl << ">> musrfit **ERROR** no raw run data avaliable." << std::endl;
726 success = false;
727 }
728 }
729 if (success) {
730 if (rrd->GetRunTitle()->Length() > 0)
731 msrHandler->SetMsrTitle(*rrd->GetRunTitle());
732 }
733 }
734 }
735
736 // generate the necessary fit histogramms for the fit
737 std::unique_ptr<PRunListCollection> runListCollection;
738 if (success) {
739 // feed all the necessary histogramms for the fit
740 runListCollection = std::make_unique<PRunListCollection>(msrHandler.get(), dataHandler.get());
741 for (unsigned int i=0; i < msrHandler->GetMsrRunList()->size(); i++) {
742 success = runListCollection->Add(i, kFit);
743 if (!success) {
744 std::cout << std::endl << ">> musrfit **ERROR** Couldn't handle run no " << i+1 << ": ";
745 std::cout << (*msrHandler->GetMsrRunList())[i].GetRunName()->Data();
746 break;
747 }
748 }
749 }
750
751 // start timeout thread
752 std::unique_ptr<TThread> th;
753 if (timeout_enabled) {
754 static pid_t musrfit_pid = getpid();
755 th = std::make_unique<TThread>(musrfit_timeout, (void*)&musrfit_pid);
756 if (th) {
757 th->Run();
758 }
759 }
760
761 // do fitting
762 std::unique_ptr<PFitter> fitter;
763 if (success) {
764 fitter = std::make_unique<PFitter>(msrHandler.get(), runListCollection.get(), chisq_only, yaml_out);
765 if (fitter->IsValid()) {
766 fitter->DoFit();
767 if (!fitter->IsScanOnly())
768 msrHandler->SetMsrStatisticConverged(fitter->HasConverged());
769 }
770 }
771
772 // write log file
773 if (success && !chisq_only) {
774 if (!fitter->IsScanOnly()) {
775 status = msrHandler->WriteMsrLogFile();
776 if (status != PMUSR_SUCCESS) {
777 switch (status) {
779 std::cout << std::endl << ">> musrfit **ERROR** couldn't write mlog-file" << std::endl << std::endl;
780 break;
782 std::cout << std::endl << ">> musrfit **ERROR** couldn't generate mlog-file name" << std::endl << std::endl;
783 break;
784 default:
785 std::cout << std::endl << ">> musrfit **UNKOWN ERROR** when trying to write the mlog-file" << std::endl << std::endl;
786 break;
787 }
788 }
789 }
790 }
791
792 // check if dump is wanted
793 if (success && !dump.IsNull()) {
794 std::cout << std::endl << "will write dump file ..." << std::endl;
795 dump.ToLower();
796 if (dump.Contains("ascii"))
797 musrfit_dump_ascii(filename.c_str(), runListCollection.get());
798 else if (dump.Contains("root"))
799 musrfit_dump_root(filename.c_str(), runListCollection.get());
800 else
801 std::cout << std::endl << "do not know format " << dump.Data() << ", sorry :-| " << std::endl;
802 }
803
804 // rename MINUIT2.OUTPUT and MINUIT2.root file if wanted
805 if (success) {
806 if (keep_mn2_output && !chisq_only && !fitter->IsScanOnly()) {
807 // 1st rename MINUIT2.OUTPUT
808 TString fln = TString(filename.c_str());
809 fln.ReplaceAll(".msr", "-mn2.output");
810 gSystem->CopyFile("MINUIT2.OUTPUT", fln.Data(), kTRUE);
811
812 // 2nd rename MINUIT2.ROOT
813 fln = TString(filename.c_str());
814 fln.ReplaceAll(".msr", "-mn2.root");
815 gSystem->CopyFile("MINUIT2.root", fln.Data(), kTRUE);
816 }
817 }
818
819 if (success) {
820 if (!chisq_only && !fitter->IsScanOnly()) {
821 // swap msr- and mlog-file
822 std::cout << std::endl << ">> swapping msr-, mlog-file ..." << std::endl;
823 // copy msr-file -> __temp.msr
824 gSystem->CopyFile(filename.c_str(), "__temp.msr", kTRUE);
825 // copy mlog-file -> msr-file
826 TString fln = TString(filename.c_str());
827 fln.ReplaceAll(".msr", ".mlog");
828 gSystem->CopyFile(fln.Data(), filename.c_str(), kTRUE);
829 // copy __temp.msr -> mlog-file
830 gSystem->CopyFile("__temp.msr", fln.Data(), kTRUE);
831 // delete __temp.msr
832 TSystemFile tmp("__temp.msr", "./");
833 tmp.Delete();
834 }
835 }
836
837 if (th && timeout_enabled) {
838 th->Kill();
839 }
840
841 std::cout << std::endl << "done ..." << std::endl;
842
843 return PMUSR_SUCCESS;
844}
845
846// end ---------------------------------------------------------------------
847
#define PMUSR_MSR_LOG_FILE_WRITE_ERROR
Failed to write to MSR log file.
Definition PMusr.h:72
#define PMUSR_TOKENIZE_ERROR
Error during tokenization/parsing of input.
Definition PMusr.h:70
@ kFit
Fitting mode - perform least-squares fit to data.
Definition PMusr.h:429
#define PMUSR_SUCCESS
Successful operation completion.
Definition PMusr.h:58
std::vector< PMsrRunBlock > PMsrRunList
Definition PMusr.h:1263
#define PMUSR_MSR_FILE_NOT_FOUND
MSR file could not be found at specified path.
Definition PMusr.h:64
#define PMUSR_MSR_SYNTAX_ERROR
Syntax error detected in MSR file content.
Definition PMusr.h:68
#define PMUSR_WRONG_STARTUP_SYNTAX
Incorrect startup command syntax provided.
Definition PMusr.h:62
const char * startup_path_name
return status
int parseXmlFile(TSAXParser *, const char *)
Replacement function for TSAXParser::ParseFile().
virtual const TString * GetRunTitle()
Definition PMusr.h:859
virtual const PDoubleVector * GetValue()
Returns pointer to data value vector (asymmetry, counts, or y-data)
Definition PMusr.h:482
virtual const PDoubleVector * GetTheory()
Returns pointer to theory value vector.
Definition PMusr.h:488
virtual Double_t GetDataTimeStep()
Returns the time bin width for data in microseconds (μs)
Definition PMusr.h:473
virtual const PDoubleVector * GetError()
Returns pointer to data error vector (statistical uncertainties)
Definition PMusr.h:484
virtual Double_t GetDataTimeStart()
Returns the start time of the data set in microseconds (μs)
Definition PMusr.h:471
Manager class for all processed μSR run data during fitting.
virtual PRunData * GetSingleHisto(UInt_t index, EDataSwitch tag=kIndex)
Retrieves processed data for a single histogram run.
virtual UInt_t GetNoOfMuMinus() const
Returns the number of μ⁻ runs in the collection.
virtual UInt_t GetNoOfSingleHisto() const
Returns the number of single histogram runs in the collection.
virtual UInt_t GetNoOfAsymmetryRRF() const
Returns the number of asymmetry RRF runs in the collection.
virtual PRunData * GetSingleHistoRRF(UInt_t index, EDataSwitch tag=kIndex)
Retrieves processed data for a single histogram RRF run.
virtual PRunData * GetNonMusr(UInt_t index, EDataSwitch tag=kIndex)
Retrieves processed data for a non-μSR run.
virtual PRunData * GetAsymmetry(UInt_t index, EDataSwitch tag=kIndex)
Retrieves processed data for an asymmetry run.
virtual UInt_t GetNoOfNonMusr() const
Returns the number of non-μSR runs in the collection.
virtual UInt_t GetNoOfAsymmetry() const
Returns the number of asymmetry runs in the collection.
virtual UInt_t GetNoOfSingleHistoRRF() const
Returns the number of single histogram RRF runs in the collection.
virtual PRunData * GetMuMinus(UInt_t index, EDataSwitch tag=kIndex)
Retrieves processed data for a μ⁻ run.
virtual PRunData * GetAsymmetryRRF(UInt_t index, EDataSwitch tag=kIndex)
Retrieves processed data for an asymmetry RRF run.
int main(int argc, char *argv[])
Definition musrfit.cpp:445
void musrfit_write_ascii(TString fln, PRunData *data, int runCounter)
Definition musrfit.cpp:152
static int timeout
Definition musrfit.cpp:72
void musrfit_dump_ascii(const char *fileName, PRunListCollection *runList)
Definition musrfit.cpp:192
void * musrfit_timeout(void *args)
Definition musrfit.cpp:79
void musrfit_write_root(TFile &f, TString fln, PRunData *data, int runCounter)
Definition musrfit.cpp:290
void musrfit_dump_root(const char *fileName, PRunListCollection *runList)
Definition musrfit.cpp:335
void musrfit_syntax()
Definition musrfit.cpp:96
Bool_t writeExpectedChisq
if set to true, expected chisq and chisq per block will be written
Definition PMusr.h:1397
Bool_t estimateN0
if set to true, for single histogram fits N0 will be estimated
Definition PMusr.h:1398