musrfit  1.9.2
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-2023 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 
71 static int timeout = 3600; // given in (sec)
72 
73 //--------------------------------------------------------------------------
78 void* 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 << " --dump <type> is writing a data file with the fit data and the theory";
125  std::cout << std::endl << " <type> can be 'ascii', 'root'";
126  std::cout << std::endl << " --timeout <timeout_tag>: overwrites to predefined timeout of " << timeout << " (sec).";
127  std::cout << std::endl << " <timeout_tag> <= 0 means timeout facility is not enabled. <timeout_tag> = nn";
128  std::cout << std::endl << " will set the timeout to nn (sec).";
129  std::cout << std::endl;
130  std::cout << std::endl << " At the end of a fit, musrfit writes the fit results into an <mlog-file> and";
131  std::cout << std::endl << " swaps them, i.e. in the <msr-file> you will find the fit results and in the";
132  std::cout << std::endl << " <mlog-file> your initial guess values.";
133  std::cout << std::endl << std::endl;
134 }
135 
136 //--------------------------------------------------------------------------
144 void musrfit_write_ascii(TString fln, PRunData *data, int runCounter)
145 {
146  // generate dump file name
147  TString count("_");
148  count += runCounter;
149  Ssiz_t index = fln.Last('.');
150  fln.Insert(index, count);
151 
152  std::ofstream f;
153 
154  // open dump-file
155  f.open(fln.Data(), std::iostream::out);
156  if (!f.is_open()) {
157  std::cout << std::endl << "Couldn't open dump (" << fln.Data() << ") file for writting, sorry ...";
158  std::cout << std::endl;
159  return;
160  }
161 
162  // dump data
163  f << "% number of data values = " << data->GetValue()->size() << std::endl;
164  f << "% time (us), value, error, theory" << std::endl;
165  double time;
166  for (unsigned int i=0; i<data->GetValue()->size(); i++) {
167  time = data->GetDataTimeStart() + (double)i*data->GetDataTimeStep();
168  f << time << ", " << data->GetValue()->at(i) << ", " << data->GetError()->at(i) << ", " << data->GetTheory()->at(i) << std::endl;
169  }
170 
171  // close file
172  f.close();
173 }
174 
175 //--------------------------------------------------------------------------
184 void musrfit_dump_ascii(char *fileName, PRunListCollection *runList)
185 {
186  TString fln(fileName);
187  fln.ReplaceAll(".msr", ".dat");
188 
189  // go through the run list, get the data and dump it in a file
190 
191  int runCounter = 0;
192  PRunData *data;
193 
194  // single histos
195  unsigned int size = runList->GetNoOfSingleHisto();
196  if (size > 0) {
197  for (unsigned int i=0; i<size; i++) {
198  data = runList->GetSingleHisto(i);
199  if (data) {
200  // dump data
201  musrfit_write_ascii(fln, data, runCounter);
202  runCounter++;
203  }
204  }
205  }
206 
207  // single histos
208  size = runList->GetNoOfSingleHistoRRF();
209  if (size > 0) {
210  for (unsigned int i=0; i<size; i++) {
211  data = runList->GetSingleHistoRRF(i);
212  if (data) {
213  // dump data
214  musrfit_write_ascii(fln, data, runCounter);
215  runCounter++;
216  }
217  }
218  }
219 
220  // asymmetry
221  size = runList->GetNoOfAsymmetry();
222  if (size > 0) {
223  for (unsigned int i=0; i<size; i++) {
224  data = runList->GetAsymmetry(i);
225  if (data) {
226  // dump data
227  musrfit_write_ascii(fln, data, runCounter);
228  runCounter++;
229  }
230  }
231  }
232 
233  // asymmetry RRF
234  size = runList->GetNoOfAsymmetryRRF();
235  if (size > 0) {
236  for (unsigned int i=0; i<size; i++) {
237  data = runList->GetAsymmetryRRF(i);
238  if (data) {
239  // dump data
240  musrfit_write_ascii(fln, data, runCounter);
241  runCounter++;
242  }
243  }
244  }
245 
246  // muMinus
247  size = runList->GetNoOfMuMinus();
248  if (size > 0) {
249  for (unsigned int i=0; i<size; i++) {
250  data = runList->GetMuMinus(i);
251  if (data) {
252  // dump data
253  musrfit_write_ascii(fln, data, runCounter);
254  runCounter++;
255  }
256  }
257  }
258 
259  // nonMusr
260  size = runList->GetNoOfNonMusr();
261  if (size > 0) {
262  for (unsigned int i=0; i<size; i++) {
263  data = runList->GetNonMusr(i);
264  if (data) {
265  // dump data
266  musrfit_write_ascii(fln, data, runCounter);
267  runCounter++;
268  }
269  }
270  }
271 }
272 
273 //--------------------------------------------------------------------------
282 void musrfit_write_root(TFile &f, TString fln, PRunData *data, int runCounter)
283 {
284  char name[128];
285 
286  TString title = fln.Copy();
287  snprintf(name, sizeof(name), "_%d", runCounter);
288  title.ReplaceAll(".root", name);
289 
290  snprintf(name, sizeof(name),"c%d", runCounter);
291 
292  std::unique_ptr<TCanvas> c = std::make_unique<TCanvas>(name, title.Data(), 10, 10, 800, 600);
293 
294  // create histos
295  Double_t diff = data->GetDataTimeStep();
296  Double_t start = -diff/2.0;
297  Double_t end = data->GetDataTimeStep()*data->GetValue()->size();
298  std::unique_ptr<TH1F> hdata = std::make_unique<TH1F>("hdata", "run data", data->GetValue()->size(), start, end);
299  std::unique_ptr<TH1F> htheo = std::make_unique<TH1F>("htheo", "run theory", data->GetValue()->size(), start, end);
300 
301  // fill data
302  for (unsigned int i=0; i<data->GetValue()->size(); i++) {
303  hdata->SetBinContent(i+1, data->GetValue()->at(i));
304  hdata->SetBinError(i+1, data->GetError()->at(i));
305  htheo->SetBinContent(i+1, data->GetTheory()->at(i));
306  }
307 
308  hdata->SetMarkerStyle(20);
309  hdata->Draw("*H HIST E1");
310 
311  htheo->SetLineColor(2);
312  htheo->SetLineWidth(3);
313  htheo->Draw("C SAME");
314 
315  f.WriteTObject(c.get());
316 }
317 
318 //--------------------------------------------------------------------------
327 void musrfit_dump_root(char *fileName, PRunListCollection *runList)
328 {
329  TString fln(fileName);
330  fln.ReplaceAll(".msr", ".root");
331 
332  TFile f(fln.Data(), "recreate");
333 
334  // go through the run list, get the data and dump it in a file
335 
336  int runCounter = 0;
337  PRunData *data;
338 
339  // single histos
340  unsigned int size = runList->GetNoOfSingleHisto();
341  if (size > 0) {
342  for (unsigned int i=0; i<size; i++) {
343  data = runList->GetSingleHisto(i);
344  if (data) {
345  // dump data
346  musrfit_write_root(f, fln, data, runCounter);
347  runCounter++;
348  }
349  }
350  }
351 
352  // single histo RRF
353  size = runList->GetNoOfSingleHistoRRF();
354  if (size > 0) {
355  for (unsigned int i=0; i<size; i++) {
356  data = runList->GetSingleHistoRRF(i);
357  if (data) {
358  // dump data
359  musrfit_write_root(f, fln, data, runCounter);
360  runCounter++;
361  }
362  }
363  }
364 
365  // asymmetry
366  size = runList->GetNoOfAsymmetry();
367  if (size > 0) {
368  for (unsigned int i=0; i<size; i++) {
369  data = runList->GetAsymmetry(i);
370  if (data) {
371  // dump data
372  musrfit_write_root(f, fln, data, runCounter);
373  runCounter++;
374  }
375  }
376  }
377 
378  // asymmetry RRF
379  size = runList->GetNoOfAsymmetryRRF();
380  if (size > 0) {
381  for (unsigned int i=0; i<size; i++) {
382  data = runList->GetAsymmetryRRF(i);
383  if (data) {
384  // dump data
385  musrfit_write_root(f, fln, data, runCounter);
386  runCounter++;
387  }
388  }
389  }
390 
391  // muMinus
392  size = runList->GetNoOfMuMinus();
393  if (size > 0) {
394  for (unsigned int i=0; i<size; i++) {
395  data = runList->GetMuMinus(i);
396  if (data) {
397  // dump data
398  musrfit_write_root(f, fln, data, runCounter);
399  runCounter++;
400  }
401  }
402  }
403 
404  // nonMusr
405  size = runList->GetNoOfNonMusr();
406  if (size > 0) {
407  for (unsigned int i=0; i<size; i++) {
408  data = runList->GetNonMusr(i);
409  if (data) {
410  // dump data
411  musrfit_write_root(f, fln, data, runCounter);
412  runCounter++;
413  }
414  }
415  }
416 
417  f.Close("R");
418 }
419 
420 //--------------------------------------------------------------------------
437 int main(int argc, char *argv[])
438 {
439  bool show_syntax = false;
440  int status;
441  bool keep_mn2_output = false;
442  bool chisq_only = false;
443  bool title_from_data_file = false;
444  bool timeout_enabled = true;
445  PStartupOptions startup_options;
446  startup_options.writeExpectedChisq = false;
447  startup_options.estimateN0 = false;
448  int number_of_cores=1;
449 
450 #ifdef HAVE_GOMP
451  number_of_cores = omp_get_num_procs();
452 #endif
453 
454  TString dump("");
455  char filename[1024];
456 
457  // check syntax
458  if (argc < 2) {
459  musrfit_syntax();
461  }
462 
463  // add default shared library path /usr/local/lib if not already persent
464  const char *dsp = gSystem->GetDynamicPath();
465  if (strstr(dsp, "/usr/local/lib") == nullptr)
466  gSystem->AddDynamicPath("/usr/local/lib");
467 
468  if (argc == 2) {
469  if (!strcmp(argv[1], "--version")) {
470 #ifdef HAVE_CONFIG_H
471 #ifdef HAVE_GIT_REV_H
472  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;
473 #else
474  std::cout << std::endl << "musrfit version: " << PACKAGE_VERSION << " (" << BUILD_TYPE << "), ROOT version: " << ROOT_VERSION_USED << std::endl << std::endl;
475 #endif
476 #else
477 #ifdef HAVE_GIT_REV_H
478  std::cout << std::endl << "musrfit git-branch: " << GIT_BRANCH << ", git-rev: " << GIT_CURRENT_SHA1 << std::endl << std::endl;
479 #else
480  std::cout << std::endl << "musrfit version: unknown" << std::endl << std::endl;
481 #endif
482 #endif
483  return PMUSR_SUCCESS;
484  } else if (!strcmp(argv[1], "--nexus-support")) {
485 #ifdef PNEXUS_ENABLED
486  std::cout << std::endl << ">> musrfit: NeXus support enabled." << std::endl << std::endl;
487 #else
488  std::cout << std::endl << "musrfit: NeXus support NOT enabled." << std::endl << std::endl;
489 #endif
490  return PMUSR_SUCCESS;
491  } else if (!strcmp(argv[1], "--show-dynamic-path")) {
492  std::cout << std::endl << "musrfit: internal dynamic search paths for shared libraries/root dictionaries:";
493  std::cout << std::endl << " '" << gSystem->GetDynamicPath() << "'" << std::endl << std::endl;
494  return PMUSR_SUCCESS;
495  } else if (!strcmp(argv[1], "--help")) {
496  show_syntax = true;
497  }
498  }
499 
500  if (show_syntax) {
501  musrfit_syntax();
503  }
504 
505  strcpy(filename, "");
506  for (int i=1; i<argc; i++) {
507  if (strstr(argv[i], ".msr")) {
508  strncpy(filename, argv[i], sizeof(filename));
509  } else if (!strcmp(argv[i], "-k") || !strcmp(argv[i], "--keep-mn2-output")) {
510  keep_mn2_output = true;
511  } else if (!strcmp(argv[i], "-c") || !strcmp(argv[i], "--chisq-only")) {
512  chisq_only = true;
513  } else if (!strcmp(argv[i], "-t") || !strcmp(argv[i], "--title-from-data-file")) {
514  title_from_data_file = true;
515  } else if (!strcmp(argv[i], "--dump")) {
516  if (i<argc-1) {
517  dump = TString(argv[i+1]);
518  i++;
519  } else {
520  std::cerr << std::endl << ">> musrfit: **ERROR** found option --dump without <type>" << std::endl;
521  show_syntax = true;
522  break;
523  }
524  } else if (!strcmp(argv[i], "-e") || !strcmp(argv[i], "--estimateN0")) {
525  startup_options.estimateN0 = true;
526  } else if (!strcmp(argv[i], "-p") || !strcmp(argv[i], "--per-run-block-chisq")) {
527  startup_options.writeExpectedChisq = true;
528  } else if (!strcmp(argv[i], "-n") || !strcmp(argv[i], "--no-of-cores-avail")) {
529 #ifdef HAVE_GOMP
530  std::cout << std::endl;
531  std::cout << "musrfit: maxmimal number of cores for OpenMP available: " << omp_get_num_procs() << std::endl;
532  std::cout << std::endl;
533 #else
534  std::cout << std::endl;
535  std::cout << ">> musrfit: this option is only vaild if OpenMP is present. This seems not to be the case here. Sorry!" << std::endl;
536  std::cout << std::endl;
537 #endif
538  return PMUSR_SUCCESS;
539  } else if (!strcmp(argv[i], "-u") || !strcmp(argv[i], "--use-no-of-threads")) {
540  if (i<argc-1) {
541  TString str(argv[i+1]);
542  if (str.IsFloat()) {
543  int ival = str.Atoi();
544  if (ival <= 0) {
545  std::cerr << std::endl << ">> musrfit: **ERROR** found option --use-no-of-threads with <number> <= 0" << std::endl;
546  std::cerr << " This doesn't make any sense." << std::endl;
547  show_syntax = true;
548  break;
549  } else if (ival > number_of_cores) {
550 #ifdef HAVE_GOMP
551  std::cerr << std::endl << ">> musrfit: **WARNING** found option --use-no-of-threads with <number>=" << ival << " > max available cores=" << number_of_cores << "." << std::endl;
552  std::cerr << " Will set <number> to max available cores." << std::endl;
553 #else
554  std::cerr << std::endl << ">> musrfit: **WARNING** option --use-no-of-threads can only be used if OpenMP is available." << std::endl;
555  std::cerr << " Here it is not the case, and hence this option will be ignored." << std::endl;
556 #endif
557  } else {
558  number_of_cores = ival;
559  }
560  } else {
561  std::cerr << std::endl << ">> musrfit: **ERROR** found option --use-no-of-threads where <number> it not a number: '" << argv[i+1] << "'" << std::endl;
562  show_syntax = true;
563  break;
564  }
565  i++;
566  } else {
567  std::cerr << std::endl << ">> musrfit: **ERROR** found option --use-no-of-threads without <number>" << std::endl;
568  show_syntax = true;
569  break;
570  }
571  } else if (!strcmp(argv[i], "--timeout")) {
572  if (i<argc-1) {
573  TString str(argv[i+1]);
574  if (str.IsFloat()) {
575  timeout = str.Atoi();
576  if (timeout <= 0) {
577  timeout_enabled = false;
578  std::cout << std::endl << ">> musrfit: timeout disabled." << std::endl;
579  }
580  } else {
581  std::cerr << std::endl << ">> musrfit: **ERROR** found option --timeout with unsupported <timeout_tag> = " << argv[i+1] << std::endl;
582  show_syntax = true;
583  break;
584  }
585  i++;
586  } else {
587  std::cerr << std::endl << ">> musrfit: **ERROR** found option --timeout without <timeout_tag>" << std::endl;
588  show_syntax = true;
589  break;
590  }
591  } else {
592  show_syntax = true;
593  break;
594  }
595  }
596 
597  // check if a filename is present
598  if (strlen(filename) == 0) {
599  show_syntax = true;
600  std::cout << std::endl << ">> musrfit **ERROR** no msr-file present!" << std::endl;
601  }
602 
603  if (show_syntax) {
604  musrfit_syntax();
606  }
607 
608  // check if dump string does make sense
609  if (!dump.IsNull()) {
610  dump.ToLower();
611  if (!dump.Contains("ascii") && !dump.Contains("root")) {
612  std::cerr << std::endl << ">> musrfit: **ERROR** found option --dump with unsupported <type> = " << dump << std::endl;
613  musrfit_syntax();
615  }
616  }
617 
618  // read startup file
619  char startup_path_name[128];
620  std::unique_ptr<TSAXParser> saxParser = std::make_unique<TSAXParser>();
621  std::unique_ptr<PStartupHandler> startupHandler = std::make_unique<PStartupHandler>();
622  if (!startupHandler->StartupFileFound()) {
623  std::cerr << std::endl << ">> musrfit **WARNING** couldn't find " << startupHandler->GetStartupFilePath().Data();
624  std::cerr << std::endl;
625  } else {
626  strcpy(startup_path_name, startupHandler->GetStartupFilePath().Data());
627  saxParser->ConnectToHandler("PStartupHandler", startupHandler.get());
628  //status = saxParser->ParseFile(startup_path_name);
629  // parsing the file as above seems to lead to problems in certain environments;
630  // use the parseXmlFile function instead (see PStartupHandler.cpp for the definition)
631  status = parseXmlFile(saxParser.get(), startup_path_name);
632  // check for parse errors
633  if (status) { // error
634  std::cerr << std::endl << ">> musrfit **WARNING** Reading/parsing musrfit_startup.xml failed.";
635  std::cerr << std::endl;
636  }
637  }
638 
639 #ifdef HAVE_GOMP
640  // set omp_set_num_threads
641  omp_set_num_threads(number_of_cores);
642 #endif
643 
644  // read msr-file
645  std::unique_ptr<PMsrHandler> msrHandler;
646  if (startupHandler)
647  msrHandler = std::make_unique<PMsrHandler>(filename, &startup_options);
648  else
649  msrHandler = std::make_unique<PMsrHandler>(filename);
650  status = msrHandler->ReadMsrFile();
651  if (status != PMUSR_SUCCESS) {
652  switch (status) {
654  std::cout << std::endl << ">> musrfit **ERROR** couldn't find " << filename << std::endl << std::endl;
655  break;
657  std::cout << std::endl << ">> musrfit **SYNTAX ERROR** in file " << filename << ", full stop here." << std::endl << std::endl;
658  break;
659  default:
660  std::cout << std::endl << ">> musrfit **UNKOWN ERROR** when trying to read the msr-file" << std::endl << std::endl;
661  break;
662  }
663  return status;
664  }
665 
666  // read all the necessary runs (raw data)
667  std::unique_ptr<PRunDataHandler> dataHandler;
668  if (startupHandler)
669  dataHandler = std::make_unique<PRunDataHandler>(msrHandler.get(), startupHandler->GetDataPathList());
670  else
671  dataHandler = std::make_unique<PRunDataHandler>(msrHandler.get());
672 
673  dataHandler->ReadData();
674 
675  bool success = dataHandler->IsAllDataAvailable();
676  if (!success) {
677  std::cout << std::endl << ">> musrfit **ERROR** Couldn't read all data files, will quit ..." << std::endl;
678  }
679 
680  // 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
681  if (title_from_data_file && success) {
682  PMsrRunList *rl = msrHandler->GetMsrRunList();
683  PRawRunData *rrd = dataHandler->GetRunData(*(rl->at(0).GetRunName()));
684  if (rrd->GetRunTitle()->Length() > 0)
685  msrHandler->SetMsrTitle(*rrd->GetRunTitle());
686  }
687 
688  // generate the necessary fit histogramms for the fit
689  std::unique_ptr<PRunListCollection> runListCollection;
690  if (success) {
691  // feed all the necessary histogramms for the fit
692  runListCollection = std::make_unique<PRunListCollection>(msrHandler.get(), dataHandler.get());
693  for (unsigned int i=0; i < msrHandler->GetMsrRunList()->size(); i++) {
694  success = runListCollection->Add(i, kFit);
695  if (!success) {
696  std::cout << std::endl << ">> musrfit **ERROR** Couldn't handle run no " << i+1 << ": ";
697  std::cout << (*msrHandler->GetMsrRunList())[i].GetRunName()->Data();
698  break;
699  }
700  }
701  }
702 
703  // start timeout thread
704  std::unique_ptr<TThread> th;
705  if (timeout_enabled) {
706  pid_t musrfit_pid = getpid();
707  th = std::make_unique<TThread>(musrfit_timeout, (void*)&musrfit_pid);
708  if (th) {
709  th->Run();
710  }
711  }
712 
713  // do fitting
714  std::unique_ptr<PFitter> fitter;
715  if (success) {
716  fitter = std::make_unique<PFitter>(msrHandler.get(), runListCollection.get(), chisq_only);
717  if (fitter->IsValid()) {
718  fitter->DoFit();
719  if (!fitter->IsScanOnly())
720  msrHandler->SetMsrStatisticConverged(fitter->HasConverged());
721  }
722  }
723 
724  // write log file
725  if (success && !chisq_only) {
726  if (!fitter->IsScanOnly()) {
727  status = msrHandler->WriteMsrLogFile();
728  if (status != PMUSR_SUCCESS) {
729  switch (status) {
731  std::cout << std::endl << ">> musrfit **ERROR** couldn't write mlog-file" << std::endl << std::endl;
732  break;
734  std::cout << std::endl << ">> musrfit **ERROR** couldn't generate mlog-file name" << std::endl << std::endl;
735  break;
736  default:
737  std::cout << std::endl << ">> musrfit **UNKOWN ERROR** when trying to write the mlog-file" << std::endl << std::endl;
738  break;
739  }
740  }
741  }
742  }
743 
744  // check if dump is wanted
745  if (success && !dump.IsNull()) {
746  std::cout << std::endl << "will write dump file ..." << std::endl;
747  dump.ToLower();
748  if (dump.Contains("ascii"))
749  musrfit_dump_ascii(filename, runListCollection.get());
750  else if (dump.Contains("root"))
751  musrfit_dump_root(filename, runListCollection.get());
752  else
753  std::cout << std::endl << "do not know format " << dump.Data() << ", sorry :-| " << std::endl;
754  }
755 
756  // rename MINUIT2.OUTPUT and MINUIT2.root file if wanted
757  if (success) {
758  if (keep_mn2_output && !chisq_only && !fitter->IsScanOnly()) {
759  // 1st rename MINUIT2.OUTPUT
760  TString fln = TString(filename);
761  char ext[32];
762  strcpy(ext, "-mn2.output");
763  fln.ReplaceAll(".msr", 4, ext, strlen(ext));
764  gSystem->CopyFile("MINUIT2.OUTPUT", fln.Data(), kTRUE);
765 
766  // 2nd rename MINUIT2.ROOT
767  fln = TString(filename);
768  strcpy(ext, "-mn2.root");
769  fln.ReplaceAll(".msr", 4, ext, strlen(ext));
770  gSystem->CopyFile("MINUIT2.root", fln.Data(), kTRUE);
771  }
772  }
773 
774  if (success) {
775  if (!chisq_only && !fitter->IsScanOnly()) {
776  // swap msr- and mlog-file
777  std::cout << std::endl << ">> swapping msr-, mlog-file ..." << std::endl;
778  // copy msr-file -> __temp.msr
779  gSystem->CopyFile(filename, "__temp.msr", kTRUE);
780  // copy mlog-file -> msr-file
781  TString fln = TString(filename);
782  char ext[32];
783  strcpy(ext, ".mlog");
784  fln.ReplaceAll(".msr", 4, ext, strlen(ext));
785  gSystem->CopyFile(fln.Data(), filename, kTRUE);
786  // copy __temp.msr -> mlog-file
787  gSystem->CopyFile("__temp.msr", fln.Data(), kTRUE);
788  // delete __temp.msr
789  TSystemFile tmp("__temp.msr", "./");
790  tmp.Delete();
791  }
792  }
793 
794  std::cout << std::endl << "done ..." << std::endl;
795 
796  return PMUSR_SUCCESS;
797 }
798 
799 // end ---------------------------------------------------------------------
800 
void musrfit_write_root(TFile &f, TString fln, PRunData *data, int runCounter)
Definition: musrfit.cpp:282
virtual PRunData * GetAsymmetry(UInt_t index, EDataSwitch tag=kIndex)
virtual UInt_t GetNoOfAsymmetryRRF() const
returns the number of asymmetry RRF data sets present in the msr-file
virtual const PDoubleVector * GetError()
Definition: PMusr.h:249
void musrfit_dump_root(char *fileName, PRunListCollection *runList)
Definition: musrfit.cpp:327
Bool_t writeExpectedChisq
if set to true, expected chisq and chisq per block will be written
Definition: PMusr.h:849
std::vector< PMsrRunBlock > PMsrRunList
Definition: PMusr.h:754
Bool_t estimateN0
if set to true, for single histogram fits N0 will be estimated
Definition: PMusr.h:850
virtual UInt_t GetNoOfSingleHistoRRF() const
returns the number of single histogram RRF data sets present in the msr-file
virtual const PDoubleVector * GetTheory()
Definition: PMusr.h:251
virtual UInt_t GetNoOfMuMinus() const
returns the number of mu minus data sets present in the msr-file
void musrfit_syntax()
Definition: musrfit.cpp:95
int parseXmlFile(TSAXParser *, const char *)
#define PMUSR_MSR_FILE_NOT_FOUND
Definition: PMusr.h:47
virtual PRunData * GetMuMinus(UInt_t index, EDataSwitch tag=kIndex)
virtual PRunData * GetAsymmetryRRF(UInt_t index, EDataSwitch tag=kIndex)
void musrfit_dump_ascii(char *fileName, PRunListCollection *runList)
Definition: musrfit.cpp:184
virtual PRunData * GetSingleHisto(UInt_t index, EDataSwitch tag=kIndex)
void musrfit_write_ascii(TString fln, PRunData *data, int runCounter)
Definition: musrfit.cpp:144
int main(int argc, char *argv[])
Definition: musrfit.cpp:437
const char * startup_path_name
static int timeout
Definition: musrfit.cpp:71
virtual const PDoubleVector * GetValue()
Definition: PMusr.h:248
virtual PRunData * GetNonMusr(UInt_t index, EDataSwitch tag=kIndex)
#define PMUSR_MSR_SYNTAX_ERROR
Definition: PMusr.h:49
#define PMUSR_TOKENIZE_ERROR
Definition: PMusr.h:50
virtual PRunData * GetSingleHistoRRF(UInt_t index, EDataSwitch tag=kIndex)
#define PMUSR_WRONG_STARTUP_SYNTAX
Definition: PMusr.h:46
#define PMUSR_MSR_LOG_FILE_WRITE_ERROR
Definition: PMusr.h:51
virtual UInt_t GetNoOfNonMusr() const
returns the number of non-muSR data sets present in the msr-file
#define PMUSR_SUCCESS
Definition: PMusr.h:44
virtual Double_t GetDataTimeStart()
Definition: PMusr.h:242
virtual const TString * GetRunTitle()
Definition: PMusr.h:408
return status
void * musrfit_timeout(void *args)
Definition: musrfit.cpp:78
virtual Double_t GetDataTimeStep()
Definition: PMusr.h:243
virtual UInt_t GetNoOfSingleHisto() const
returns the number of single histogram data sets present in the msr-file
Definition: PMusr.h:220
virtual UInt_t GetNoOfAsymmetry() const
returns the number of asymmetry data sets present in the msr-file