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