ROOTPWA
pwafit.cc
Go to the documentation of this file.
1 
2 //
3 // Copyright 2009 Sebastian Neubert
4 //
5 // This file is part of rootpwa
6 //
7 // rootpwa is free software: you can redistribute it and/or modify
8 // it under the terms of the GNU General Public License as published by
9 // the Free Software Foundation, either version 3 of the License, or
10 // (at your option) any later version.
11 //
12 // rootpwa is distributed in the hope that it will be useful,
13 // but WITHOUT ANY WARRANTY; without even the implied warranty of
14 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 // GNU General Public License for more details.
16 //
17 // You should have received a copy of the GNU General Public License
18 // along with rootpwa. If not, see <http://www.gnu.org/licenses/>.
19 //
21 //-----------------------------------------------------------
22 // File and Version Information:
23 // $Id$
24 //
25 // Description:
26 // fitting program for rootpwa
27 // minimizes TPWALikelihood function
28 //
29 //
30 // Author List:
31 // Sebastian Neubert TUM (original author)
32 //
33 //
34 //-----------------------------------------------------------
35 
36 
37 #include <iostream>
38 #include <iomanip>
39 #include <vector>
40 #include <string>
41 #include <complex>
42 #include <cassert>
43 #include <time.h>
44 
45 #include "TROOT.h"
46 #include "TTree.h"
47 #include "TFile.h"
48 #include "TString.h"
49 #include "TComplex.h"
50 #include "TRandom3.h"
51 #include "Math/Minimizer.h"
52 #include "Math/Factory.h"
53 #include "TStopwatch.h"
54 
55 #include "conversionUtils.hpp"
56 #include "TPWALikelihood.h"
57 #include "TFitBin.h"
58 #include "fitResult.h"
59 #ifdef USE_CUDA
60 #include "complex.cuh"
61 #include "likelihoodInterface.cuh"
62 #endif
63 #include "amplitudeTreeLeaf.h"
64 
65 
66 using namespace std;
67 using namespace ROOT::Math;
68 using namespace rpwa;
69 
70 
71 void
72 usage(const string& progName,
73  const int errCode = 0)
74 {
75 
76  cerr << "performs PWA fit for given mass bin and list of waves" << endl
77  << endl
78  << "usage:" << endl
79  << progName
80  << " -l # -u # -w wavelist [-d amplitude directory -R -o outfile -S start value file -N -n normfile"
81  << " [-a normfile] -r rank -M minimizer [-m algorithm -g strategy -t #] -q -h]" << endl
82  << " where:" << endl
83  << " -l # lower edge of mass bin [MeV/c^2]" << endl
84  << " -u # upper edge of mass bin [MeV/c^2]" << endl
85  << " -w file path to wavelist file" << endl
86  << " -d dir path to directory with decay amplitude files (default: '.')" << endl
87 #ifdef USE_STD_COMPLEX_TREE_LEAFS
88  << " -R use .root amplitude files (default: false)" << endl
89 #else
90  << " -R use .root amplitude files [not supported; ROOT version too low]" << endl
91 #endif
92  << " -o file path to output file (default: 'fitresult.root')" << endl
93  << " -S file path to file with start values (default: none; highest priority)" << endl
94  << " -s # seed for random start values (default: 1234567)" << endl
95  << " -x # use fixed instead of random start values (default: 0.01)" << endl
96 
97  << " -N use normalization of decay amplitudes (default: false)" << endl
98  << " -n file path to normalization integral file (default: 'norm.int')" << endl
99  << " -a file path to acceptance integral file (default: 'norm.int')" << endl
100  << " -A # number of input events to normalize acceptance to" << endl
101  << " -r # rank of spin density matrix (default: 1)" << endl
102  << " -M name minimizer (default: Minuit2)" << endl
103  << " -m name minimization algorithm (optional, default: Migrad)" << endl
104  << " available minimizers: Minuit: Migrad, Simplex, Minimize, Migrad_imp" << endl
105  << " Minuit2: Migrad, Simplex, Combined, Scan, Fumili" << endl
106  << " GSLMultiMin: ConjugateFR, ConjugatePR, BFGS, BFGS2, SteepestDescent" << endl
107  << " GSLMultiFit: -" << endl
108  << " GSLSimAn: -" << endl
109  << " Linear: Robust" << endl
110  << " Fumili: -" << endl
111  << " -g # minimizer strategy: 0 = low, 1 = medium, 2 = high effort (default: 1)" << endl
112  << " -t # minimizer tolerance (default: 1e-10)" << endl
113 #ifdef USE_CUDA
114  << " -c enable CUDA acceleration (default: off)" << endl
115 #else
116  << " -c enable CUDA acceleration [not supported by your platform]" << endl
117 #endif
118  << " -q run quietly (default: false)" << endl
119  << " -h print help" << endl
120  << endl;
121  exit(errCode);
122 }
123 
124 
125 ostream&
126 printMinimizerStatus(ostream& out,
127  Minimizer& minimizer)
128 {
129  out << "minimization stopped after " << minimizer.NCalls() << " function calls. "
130  << "minimizer status summary:" << endl
131  << " total number of parameters .......................... " << minimizer.NDim() << endl
132  << " number of free parameters ........................... " << minimizer.NFree() << endl
133  << " maximum allowed number of iterations ................ " << minimizer.MaxIterations() << endl
134  << " maximum allowed number of function calls ............ " << minimizer.MaxFunctionCalls() << endl
135  << " minimizer status .................................... " << minimizer.Status() << endl
136  << " minimizer provides error and error matrix ........... " << yesNo(minimizer.ProvidesError()) << endl
137  << " minimizer has performed detailed error validation ... " << yesNo(minimizer.IsValidError()) << endl
138  << " estimated distance to minimum ....................... " << minimizer.Edm() << endl
139  << " statistical scale used for error calculation ........ " << minimizer.ErrorDef() << endl
140  << " strategy ............................................ " << minimizer.Strategy() << endl
141  << " absolute tolerance .................................. " << minimizer.Tolerance() << endl
142  << " precision ........................................... " << minimizer.Precision() << endl;
143  return out;
144 }
145 
146 
147 ostream&
148 operator <<(ostream& out,
149  Minimizer& minimizer)
150 {
151  return printMinimizerStatus(out, minimizer);
152 }
153 
154 
155 int
156 main(int argc,
157  char** argv)
158 {
159  printCompilerInfo();
160  printLibraryInfo ();
161  printSvnVersion ();
162  cout << endl;
163 
164  // force loading predefined std::complex dictionary
165  // see http://root.cern.ch/phpBB3/viewtopic.php?f=5&t=9618&p=50164
166  gROOT->ProcessLine("#include <complex>");
167 
168  // ---------------------------------------------------------------------------
169  // internal parameters
170  const string valTreeName = "pwa";
171  const string valBranchName = "fitResult_v2";
172  double defaultStartValue = 0.01;
173  bool useFixedStartValues = false;
174  double startValStep = 0.0005;
175  const unsigned int maxNmbOfIterations = 20000;
176  const unsigned int maxNmbOfFunctionCalls = 40000;
177  const bool runHesse = true;
178  const bool runMinos = false;
179  int startValSeed = 1234567;
180 
181  // ---------------------------------------------------------------------------
182  // parse command line options
183  const string progName = argv[0];
184  double massBinMin = 0; // [MeV/c^2]
185  double massBinMax = 0; // [MeV/c^2]
186  string waveListFileName = ""; // wavelist filename
187  string ampDirName = "."; // decay amplitude directory name
188  bool useRootAmps = false; // if true .root amplitude files are read
189  //bool useRootAmps = true; // if true .root amplitude files are read
190  string outFileName = "fitresult.root"; // output filename
191  string startValFileName = ""; // file with start values
192  bool useNormalizedAmps = false; // if true normalized amplitudes are used
193  string normIntFileName = ""; // file with normalization integrals
194  string accIntFileName = ""; // file with acceptance integrals
195  unsigned int numbAccEvents = 0; // number of events used for acceptance integrals
196  unsigned int rank = 1; // rank of fit
197  string minimizerType[2] = {"Minuit2", "Migrad"}; // minimizer, minimization algorithm
198  int minimizerStrategy = 1; // minimizer strategy
199  double minimizerTolerance = 1e-10; // minimizer tolerance
200  bool cudaEnabled = false; // if true CUDA kernels are activated
201  bool quiet = false;
202  extern char* optarg;
203  // extern int optind;
204  int c;
205  while ((c = getopt(argc, argv, "l:u:w:d:Ro:S:s:x::Nn:a:A:r:M:m:g:t:cqh")) != -1)
206  switch (c) {
207  case 'l':
208  massBinMin = atof(optarg);
209  break;
210  case 'u':
211  massBinMax = atof(optarg);
212  break;
213  case 'w':
214  waveListFileName = optarg;
215  break;
216  case 'd':
217  ampDirName = optarg;
218  break;
219  case 'R':
220 #ifdef USE_STD_COMPLEX_TREE_LEAFS
221  useRootAmps = true;
222 #endif
223  break;
224  case 'o':
225  outFileName = optarg;
226  break;
227  case 'S':
228  startValFileName = optarg;
229  break;
230  case 's':
231  startValSeed = atoi(optarg);
232  break;
233  case 'x':
234  if (optarg)
235  defaultStartValue = atof(optarg);
236  useFixedStartValues = true;
237  break;
238  case 'N':
239  useNormalizedAmps = true;
240  break;
241  case 'n':
242  normIntFileName = optarg;
243  break;
244  case 'a':
245  accIntFileName = optarg;
246  break;
247  case 'A':
248  numbAccEvents = atoi(optarg);
249  break;
250  case 'r':
251  rank = atoi(optarg);
252  break;
253  case 'M':
254  minimizerType[0] = optarg;
255  break;
256  case 'm':
257  minimizerType[1] = optarg;
258  break;
259  case 'g':
260  minimizerStrategy = atoi(optarg);
261  break;
262  case 't':
263  minimizerTolerance = atof(optarg);
264  break;
265  case 'c':
266 #ifdef USE_CUDA
267  cudaEnabled = true;
268 #endif
269  break;
270  case 'q':
271  quiet = true;
272  break;
273  case 'h':
274  usage(progName);
275  break;
276  }
277  if (normIntFileName.length() <= 1) {
278  normIntFileName = "norm.int";
279  printWarn << "using default normalization integral file '" << normIntFileName << "'" << endl;
280  }
281  if (accIntFileName.length() <= 1) {
282  accIntFileName = "norm.int";
283  printWarn << "using default acceptance normalization integral file "
284  << "'" << accIntFileName << "'" << endl;
285  }
286  if (waveListFileName.length() <= 1) {
287  printErr << "no wavelist file specified. aborting." << endl;
288  usage(progName, 1);
289  }
290  // report parameters
291  printInfo << "running " << progName << " with the following parameters:" << endl;
292  cout << " mass bin [" << massBinMin << ", " << massBinMax << "] MeV/c^2" << endl
293  << " path to wave list file ......................... '" << waveListFileName << "'" << endl
294  << " path to amplitude directory .................... '" << ampDirName << "'" << endl
295  << " use .root amplitude files ...................... " << yesNo(useRootAmps) << endl
296  << " path to output file ............................ '" << outFileName << "'" << endl
297  << " path to file with start values ................. '" << startValFileName << "'" << endl
298  << " seed for random start values ................... " << startValSeed << endl;
299  if (useFixedStartValues)
300  cout << " using fixed instead of random start values ..... " << defaultStartValue << endl;
301  cout << " use normalization .............................. " << yesNo(useNormalizedAmps) << endl
302  << " path to file with normalization integral ... '" << normIntFileName << "'" << endl
303  << " path to file with acceptance integral ...... '" << accIntFileName << "'" << endl
304  << " number of acceptance norm. events .......... " << numbAccEvents << endl
305  << " rank of spin density matrix .................... " << rank << endl
306  << " minimizer ...................................... " << minimizerType[0] << ", " << minimizerType[1] << endl
307  << " minimizer strategy ............................. " << minimizerStrategy << endl
308  << " minimizer tolerance ............................ " << minimizerTolerance << endl
309  << " CUDA acceleration .............................. " << enDisabled(cudaEnabled) << endl
310  << " quiet .......................................... " << yesNo(quiet) << endl;
311 
312  // ---------------------------------------------------------------------------
313  // setup likelihood function
314  printInfo << "creating and setting up likelihood function" << endl;
315  TPWALikelihood<complex<double> > L;
316  if (quiet)
317  L.setQuiet();
318  L.useNormalizedAmps(useNormalizedAmps);
319 #ifdef USE_CUDA
320  L.enableCuda(cudaEnabled);
321 #endif
322  L.init(rank, waveListFileName, normIntFileName, accIntFileName,
323  ampDirName, numbAccEvents, useRootAmps);
324  if (not quiet)
325  cout << L << endl;
326  const unsigned int nmbPar = L.NDim();
327  const unsigned int nmbEvts = L.nmbEvents();
328 
329  // ---------------------------------------------------------------------------
330  // setup minimizer
331  printInfo << "creating and setting up minimizer '" << minimizerType[0] << "' "
332  << "using algorithm '" << minimizerType[1] << "'" << endl;
333  Minimizer* minimizer = Factory::CreateMinimizer(minimizerType[0], minimizerType[1]);
334  if (not minimizer) {
335  printErr << "could not create minimizer. exiting." << endl;
336  throw;
337  }
338  minimizer->SetFunction (L);
339  minimizer->SetStrategy (minimizerStrategy);
340  minimizer->SetTolerance (minimizerTolerance);
341 
342  // setting the ErrorDef to 1 since the ROOT interface does not
343  // Propagate the value. Will do the error rescaling by hand below.
344  minimizer->SetErrorDef(1);
345  minimizer->SetPrintLevel ((quiet) ? 0 : 3);
346  minimizer->SetMaxIterations (maxNmbOfIterations);
347  minimizer->SetMaxFunctionCalls(maxNmbOfFunctionCalls);
348 
349  // ---------------------------------------------------------------------------
350  // read in fitResult with start values
351  printInfo << "reading start values from '" << startValFileName << "'" << endl;
352  const double massBinCenter = (massBinMin + massBinMax) / 2;
353  fitResult* startFitResult = NULL;
354  bool startValValid = false;
355  TFile* startValFile = NULL;
356  if (startValFileName.length() <= 2)
357  printWarn << "start value file name '" << startValFileName << "' is invalid. "
358  << "using default start values." << endl;
359  else {
360  // open root file
361  startValFile = TFile::Open(startValFileName.c_str(), "READ");
362  if (not startValFile or startValFile->IsZombie())
363  printWarn << "cannot open start value file '" << startValFileName << "'. "
364  << "using default start values." << endl;
365  else {
366  // get tree with start values
367  TTree* tree;
368  startValFile->GetObject(valTreeName.c_str(), tree);
369  if (not tree)
370  printWarn << "cannot find start value tree '"<< valTreeName << "' in file "
371  << "'" << startValFileName << "'" << endl;
372  else {
373  //startBin = new TFitBin();
374  //tree->SetBranchAddress("fitbin", &startBin);
375  startFitResult = new fitResult();
376  tree->SetBranchAddress(valBranchName.c_str(), &startFitResult);
377  // find tree entry which is closest to mass bin center
378  unsigned int bestIndex = 0;
379  double bestMass = 0;
380  for (unsigned int i = 0; i < tree->GetEntriesFast(); ++i) {
381  tree->GetEntry(i);
382  if (fabs(massBinCenter - startFitResult->massBinCenter()) <= fabs(massBinCenter - bestMass)) {
383  bestIndex = i;
384  bestMass = startFitResult->massBinCenter();
385  }
386  }
387  tree->GetEntry(bestIndex);
388  startValValid = true;
389  }
390  }
391  }
392 
393  // ---------------------------------------------------------------------------
394  // set start parameter values
395  printInfo << "setting start values for " << nmbPar << " parameters" << endl
396  << " parameter naming scheme is: V[rank index]_[IGJPCMR][isobar spec]" << endl;
397  unsigned int maxParNameLength = 0; // maximum length of parameter names
398  vector<bool> parIsFixed(nmbPar, false); // memorizes state of variables; ROOT::Math::Minimizer has no corresponding accessor
399  {
400  // use local instance of random number generator so that other
401  // code has no chance of tampering with gRandom and thus cannot
402  // affect the reproducability of the start values
403  TRandom3 random(startValSeed);
404  bool success = true;
405  for (unsigned int i = 0; i < nmbPar; ++i) {
406  const string parName = L.parName(i);
407  if (parName.length() > maxParNameLength)
408  maxParNameLength = parName.length();
409  // workaround, because Minuit2Minimizer::SetVariable() expects
410  // that variables are set consecutively. how stupid is that?
411  // so we prepare variables here and set values below
412  if ((L.parThreshold(i) == 0) or (L.parThreshold(i) < massBinCenter)) {
413  if (not minimizer->SetVariable(i, parName, 0, startValStep))
414  success = false;
415  } else {
416  if (not minimizer->SetFixedVariable(i, parName, 0.)) // fix this parameter to 0
417  success = false;
418  parIsFixed[i] = true;
419  }
420  }
421  const double sqrtNmbEvts = sqrt((double)nmbEvts);
422  vector<unsigned int> parIndices = L.orderedParIndices();
423  for (unsigned int i = 0; i < parIndices.size(); ++i) {
424  const unsigned int parIndex = parIndices[i];
425  double startVal;
426  const string parName = minimizer->VariableName(parIndex);
427  if (parName != L.parName(parIndex)) {
428  printWarn << "parameter name in minimizer and likelihood is inconsistent "
429  << "(" << parName << " vs. " << L.parName(parIndex) << ")" << endl;
430  success = false;
431  }
432  if (startValValid) {
433  // get parameter value from fitResult
434  assert(startFitResult);
435  startVal = startFitResult->fitParameter(parName.c_str());
436  } else
437  startVal = (useFixedStartValues) ? defaultStartValue
438  : random.Uniform(defaultStartValue, sqrtNmbEvts);
439  // check if parameter needs to be fixed because of threshold
440  if ((L.parThreshold(parIndex) == 0) or (L.parThreshold(parIndex) < massBinCenter)) {
441  if (startVal == 0) {
442  cout << " read start value 0 for parameter " << parName << ". "
443  << "using default start value." << endl;
444  startVal = (useFixedStartValues) ? defaultStartValue
445  : random.Uniform(defaultStartValue, sqrtNmbEvts);
446  }
447  cout << " setting parameter [" << setw(3) << i << "] "
448  << setw(maxParNameLength) << parName << " = " << maxPrecisionAlign(startVal) << endl;
449  if (not minimizer->SetVariableValue(parIndex, startVal))
450  success = false;
451  } else {
452  cout << " fixing parameter [" << setw(3) << i << "] "
453  << setw(maxParNameLength) << parName << " = 0" << endl;
454  if (not minimizer->SetVariableValue(parIndex, 0.)) // fix this parameter to 0
455  success = false;
456  parIsFixed[parIndex] = true;
457  }
458  if (not success) {
459  printErr << "something went wrong when setting log likelihood parameters. aborting." << endl;
460  throw;
461  }
462  }
463  // cleanup
464  if(startValFile) {
465  startValFile->Close();
466  delete startValFile;
467  startValFile = NULL;
468  }
469  }
470 
471  // ---------------------------------------------------------------------------
472  // find minimum of likelihood function
473  bool converged = false;
474  bool hasHesse = false;
475  printInfo << "performing minimization" << endl;
476  {
477  TStopwatch timer;
478  timer.Start();
479  bool success = minimizer->Minimize();
480  timer.Stop();
481  if (success)
482  printInfo << "minimization finished successfully. " << flush;
483  else
484  printWarn << "minimization failed. " << flush;
485  cout << "used " << flush;
486  timer.Print();
487  converged = success;
488  printInfo << *minimizer;
489  // printInfo << "covariance matrix:" <<endl;
490  // for(unsigned int i = 0; i < nmbPar; ++i)
491  // for(unsigned int j = 0; j < nmbPar; ++j)
492  // cout << " [" << i << "][" << j << "] = " << minimizer->CovMatrix(i, j) << endl;
493  if (runHesse) {
494  printInfo << "calculating Hessian matrix" << endl;
495  timer.Start();
496  success = minimizer->Hesse();
497  timer.Stop();
498  if (success)
499  printInfo << "successfully calculated Hessian matrix. " << flush;
500  else
501  printWarn << "calculation of Hessian matrix failed. " << flush;
502  cout << "used " << flush;
503  timer.Print();
504  hasHesse = success;
505  printInfo << *minimizer;
506  // printInfo << "covariance matrix:" <<endl;
507  // for(unsigned int i = 0; i < nmbPar; ++i)
508  // for(unsigned int j = 0; j < nmbPar; ++j)
509  // cout << " [" << i << "][" << j << "] = " << minimizer->CovMatrix(i, j) << endl;
510  }
511  }
512 
513  // ---------------------------------------------------------------------------
514  // print results
515  printInfo << "minimization result:" << endl;
516  vector<unsigned int> parIndices = L.orderedParIndices();
517  for (unsigned int i = 0; i< parIndices.size(); ++i) {
518  const unsigned int parIndex = parIndices[i];
519  cout << " parameter [" << setw(3) << i << "] "
520  << setw(maxParNameLength) << L.parName(parIndex) << " = ";
521  if (parIsFixed[parIndex])
522  cout << minimizer->X()[parIndex] << " (fixed)" << endl;
523  else {
524  cout << setw(12) << maxPrecisionAlign(minimizer->X() [parIndex]) << " +- "
525  << setw(12) << maxPrecisionAlign(minimizer->Errors()[parIndex]);
526  if (runMinos) {
527  double minosErrLow = 0;
528  double minosErrUp = 0;
529  const bool success = minimizer->GetMinosError(parIndex, minosErrLow, minosErrUp);
530  if (success)
531  cout << " Minos: " << "[" << minosErrLow << ", +" << minosErrUp << "]" << endl;
532  } else
533  cout << endl;
534  }
535  }
536  printInfo << "function call summary:" << endl;
537  L.printFuncInfo(cout);
538 #ifdef USE_CUDA
539  printInfo << "total CUDA kernel time: "
540  << cuda::likelihoodInterface<cuda::complex<double> >::kernelTime() << " sec" << endl;
541 #endif
542 
543  // ---------------------------------------------------------------------------
544  // write out result
545  printInfo << "writing result to '" << outFileName << "'" << endl;
546  {
547  // open output file and create tree for writing
548  TFile* outFile = new TFile(outFileName.c_str(), "UPDATE");
549  if ((not outFile) or outFile->IsZombie())
550  printWarn << "cannot open output file '" << outFileName << "'. "
551  << "no results will be written." << endl;
552  else {
553  // check whether output tree already exists
554  TTree* tree;
555  TFitBin* fitBinResult = new TFitBin();
556  fitResult* result = new fitResult();
557  outFile->GetObject(valTreeName.c_str(), tree);
558  if (not tree) {
559  printInfo << "file '" << outFileName << "' is empty. "
560  << "creating new tree '" << valTreeName << "' for PWA result." << endl;
561  tree = new TTree(valTreeName.c_str(), valTreeName.c_str());
562  tree->Branch("fitbin", &fitBinResult); // depricated legacy branch
563  tree->Branch(valBranchName.c_str(), &result);
564  } else {
565  tree->SetBranchAddress("fitbin", &fitBinResult);
566  tree->SetBranchAddress(valBranchName.c_str(), &result);
567  }
568 
569  {
570  // get data structures to construct fitResult
571  vector<std::complex<double> > prodAmps; // production amplitudes
572  vector<string> prodAmpNames; // names of production amplitudes used in fit
573  vector<pair<int,int> > fitParCovMatrixIndices; // indices of fit parameters for real and imaginary part in covariance matrix matrix
574  L.buildProdAmpArrays(minimizer->X(), prodAmps, fitParCovMatrixIndices, prodAmpNames, true);
575  TMatrixT<double> fitParCovMatrix(nmbPar, nmbPar); // covariance matrix of fit parameters
576  for(unsigned int i = 0; i < nmbPar; ++i)
577  for(unsigned int j = 0; j < nmbPar; ++j)
578  // The factor 0.5 is needed because
579  // MINUIT by default assumes a Chi2
580  // function and not a loglikeli
581  // (see Minuit manual!)
582  // Note: SetErrorDef in ROOT does not work
583  fitParCovMatrix[i][j] = 0.5* minimizer->CovMatrix(i, j);
584  const unsigned int nmbWaves = L.nmbWaves() + 1; // flat wave is not included in L.nmbWaves()
585  TCMatrix normIntegral(nmbWaves, nmbWaves); // normalization integral over full phase space without acceptance
586  TCMatrix accIntegral (nmbWaves, nmbWaves); // normalization integral over full phase space with acceptance
587  vector<double> phaseSpaceIntegral;
588  L.getIntegralMatrices(normIntegral, accIntegral, phaseSpaceIntegral);
589  const int normNmbEvents = (useNormalizedAmps) ? 1 : L.nmbEvents(); // number of events to normalize to
590 
591  cout << "filling fitResult:" << endl
592  << " number of fit parameters ............... " << nmbPar << endl
593  << " number of production amplitudes ........ " << prodAmps.size() << endl
594  << " number of production amplitude names ... " << prodAmpNames.size() << endl
595  << " number of wave names ................... " << nmbWaves << endl
596  << " number of cov. matrix indices .......... " << fitParCovMatrixIndices.size() << endl
597  << " dimension of covariance matrix ......... " << fitParCovMatrix.GetNrows() << " x " << fitParCovMatrix.GetNcols() << endl
598  << " dimension of normalization matrix ...... " << normIntegral.nrows() << " x " << normIntegral.ncols() << endl
599  << " dimension of acceptance matrix ......... " << accIntegral.nrows() << " x " << accIntegral.ncols() << endl;
600  result->fill(L.nmbEvents(),
601  normNmbEvents,
602  massBinCenter,
603  minimizer->MinValue(),
604  rank,
605  prodAmps,
606  prodAmpNames,
607  fitParCovMatrix,
608  fitParCovMatrixIndices,
609  normIntegral, // contains the sqrt of the integral matrix diagonal elements!!!
610  phaseSpaceIntegral,
611  converged,
612  hasHesse);
613  }
614 
615  if (1) {
616  // get data structures to construct TFitBin
617  vector<TComplex> prodAmplitudes; // production amplitudes
618  vector<pair<int,int> > indices; // indices for error matrix access
619  vector<TString> waveNames; // contains rank information
620  {
621  vector<std::complex<double> > V;
622  vector<string> names;
623  L.buildProdAmpArrays(minimizer->X(), V, indices, names, true);
624  // convert to TComplex;
625  for (unsigned int i = 0; i < V.size(); ++i)
626  prodAmplitudes.push_back(TComplex(V[i].real(), V[i].imag()));
627  assert(indices.size() == prodAmplitudes.size());
628  for(unsigned int i = 0; i < names.size(); ++i)
629  waveNames.push_back(TString(names[i].c_str()));
630  }
631  vector<TString> waveTitles; // without rank
632  {
633  const vector<string>& titles = L.waveNames();
634  for(unsigned int i = 0; i < titles.size(); ++i)
635  waveTitles.push_back(TString(titles[i].c_str()));
636  waveTitles.push_back("flat");
637  }
638  // error matrix
639  TMatrixD errMatrix(nmbPar, nmbPar);
640  for(unsigned int i = 0; i < nmbPar; ++i)
641  for(unsigned int j = 0; j < nmbPar; ++j)
642  errMatrix[i][j] = minimizer->CovMatrix(i, j);
643  // normalization integral and acceptance matrix
644  cout << " setting up integrals" << endl;
645  const unsigned int n = waveTitles.size();
646  TCMatrix integralMatrix(n, n);
647  TCMatrix accMatrix (n, n);
648  vector<double> phaseSpaceIntegral;
649  L.getIntegralMatrices(integralMatrix, accMatrix, phaseSpaceIntegral);
650  //integralMatrix.Print();
651  // representation of number of events depends on whether normalization was done
652  const int nmbEvt = (useNormalizedAmps) ? 1 : L.nmbEvents();
653 
654  cout << "filling TFitBin:" << endl
655  << " number of fit parameters ........... " << nmbPar << endl
656  << " number of production amplitudes .... " << prodAmplitudes.size() << endl
657  << " number of indices .................. " << indices.size() << endl
658  << " number of wave names (with rank) ... " << waveNames.size() << endl
659  << " number of wave titles (w/o rank) ... " << waveTitles.size() << endl
660  << " dimension of error matrix .......... " << errMatrix.GetNrows() << endl
661  << " dimension of integral matrix ....... " << integralMatrix.nrows() << endl
662  << " dimension of acceptance matrix ..... " << accMatrix.nrows() << endl;
663  fitBinResult->fill(prodAmplitudes,
664  indices,
665  waveNames,
666  nmbEvt,
667  L.nmbEvents(), // raw number of data events
668  massBinCenter,
669  integralMatrix,
670  errMatrix,
671  minimizer->MinValue(),
672  rank);
673  //fitBinResult->PrintParameters();
674  }
675 
676  // write result to file
677  tree->Fill();
678  tree->Write("", TObject::kOverwrite);
679  outFile->Close();
680  }
681  }
682 
683  if (minimizer)
684  delete minimizer;
685  return 0;
686 }
Generated by   doxygen 1.8.1.2