ROOTPWA
pwaConstraintFit.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 //fitting program for rootpwa
22 
23 #include <iostream>
24 #include <iomanip>
25 #include <vector>
26 #include <string>
27 #include <complex>
28 #include <cassert>
29 
30 #include "TTree.h"
31 #include "TFile.h"
32 #include "TString.h"
33 #include "TComplex.h"
34 #include "Math/Minimizer.h"
35 #include "Math/Factory.h"
36 
37 #include "reportingUtils.hpp"
38 #include "TFitBin.h"
39 #include "TPWALikelihoodC.h"
40 
41 
42 using namespace std;
43 using namespace ROOT::Math;
44 using namespace rpwa;
45 
46 
47 //int lineno = 1; // global variables needed for lex (not understood)
48 //char *progname;
49 
50 
51 void usage(const char* progName,
52  const int errCode = 0)
53 {
54  cerr << "usage:" << endl
55  << progName << " -l # -u # [-i -h -q -N] [-n normfile -s start values -r rank] -w wavelist [-o outfile] " << endl
56  << " where:" << endl
57  << " -l # lower edge of mass bin [MeV/c^2]" << endl
58  << " -u # upper edge of mass bin [MeV/c^2]" << endl
59  << " -w file path to wavelist file" << endl
60  << " -c file path to constraint config file" << endl
61  << " -o file path to output file (default: 'fitresult.root')" << endl
62  << " -a file path to acceptance file" << endl
63  << " -n file path to normalization file" << endl
64  << " -N use normalization of decay amplitudes (default: false)" << endl
65  << " -S file path to file with start values" <<endl
66  << " -r # rank of spin density matrix (default: 1)" <<endl
67  << " -q run quietly" << endl
68  << " -h print help" << endl
69  << endl;
70  exit(errCode);
71 }
72 
73 
74 int main(int argc, char** argv)
75 {
76  // ---------------------------------------------------------------------------
77  // internal parameters
78  const TString valTreeName = "pwa";
79  const TString valBranchName = "fitbin";
80  const double defaultStartValue = 0.01;
81  double defaultMinStep = 0.0005;
82  const unsigned int maxNmbOfIterations = 20000;
83  const bool runHesse = false;
84  const bool runMinos = false;
85  const string minimizerType[2] = {"Minuit2", "Migrad"};
86 // const string minimizerType[2] = {"Minuit2", "Simplex"};
87 // const string minimizerType[2] = {"Minuit2", "Combined"};
88 // const string minimizerType[2] = {"Minuit2", "Scan"};
89 // const string minimizerType[2] = {"Minuit2", "Fumili"};
90 // const string minimizerType[2] = {"GSLMultiMin", "ConjugateFR"};
91 // const string minimizerType[2] = {"GSLMultiMin", "ConjugatePR"};
92 // const string minimizerType[2] = {"GSLMultiMin", "BFGS"};
93 // const string minimizerType[2] = {"GSLMultiMin", "BFGS2"};
94 // const string minimizerType[2] = {"GSLMultiMin", "SteepestDescent"};
95 // const string minimizerType[2] = {"GSLMultiFit", ""};
96 // const string minimizerType[2] = {"GSLSimAn", ""};
97 
98  // ---------------------------------------------------------------------------
99  // parse command line options
100  char* progName = argv[0];
101  double binLowMass = 0;
102  double binHighMass = 0;
103  TString waveListFileName = ""; // wavelist filename
104  TString waveConstraintConfigFileName = ""; // constraint config file name
105  TString outFileName = "fitresult.root"; // output filename
106  TString startValFileName = ""; // file with start values
107  TString normFileName; // file with normalization integrals
108  TString accFileName; // file with acceptance integrals
109  bool quiet = false;
110  //bool useStartVal = false; // are there start values?
111  bool useNorm = false;
112  unsigned int rank = 1; // rank
113  extern char* optarg;
114  // extern int optind;
115  int c;
116  while ((c = getopt(argc, argv, "l:u:iw:o:c:a:S:n:Nr:qh")) != -1)
117  switch (c) {
118  case 'N':
119  useNorm = true;
120  break;
121  case 'l':
122  binLowMass = atof(optarg);
123  break;
124  case 'u':
125  binHighMass = atof(optarg);
126  break;
127  case 'o':
128  outFileName = optarg;
129  break;
130  case 'c':
131  waveConstraintConfigFileName = optarg;
132  break;
133  case 'S':
134  startValFileName = optarg;
135  //useStartVal = 1;
136  break;
137  case 'n':
138  normFileName = optarg;
139  break;
140  case 'a':
141  accFileName = optarg;
142  break;
143  case 'w':
144  waveListFileName = optarg;
145  break;
146  case 'r':
147  rank = atoi(optarg);
148  break;
149  case 'q':
150  quiet = true;
151  break;
152  case 'h':
153  usage(progName);
154  break;
155  }
156  if (normFileName.Length() <= 1) {
157  normFileName="norm.int";
158  printWarn << "Using default normalization integral file '" << normFileName << "'." << endl;
159  }
160  if (accFileName.Length() <= 1) {
161  accFileName="norm.int";
162  printWarn << "Using default acceptance normalization integral file '" << accFileName << "'." << endl;
163  }
164  if (waveListFileName.Length() <= 1) {
165  printErr << "No wavelist file specified! Aborting!" << endl;
166  usage(progName, 1);
167  }
168 
169  // ---------------------------------------------------------------------------
170  // setup likelihood function
171  printInfo << "Creating and setting up likelihood function" << endl;
172  TPWALikelihoodC L;
173  if(quiet)L.SetQuiet();
174  L.UseNormalizedAmps(useNorm);
175 
176  L.Init(waveListFileName,
177  rank,
178  normFileName,
179  accFileName,
180  20000,
181  waveConstraintConfigFileName);
182  L.LoadAmplitudes();
183 
184  const unsigned int nmbPar = L.NDim();
185 
186  // 12 parameters + flat
187  //double x[13]={0.52707,0.21068,-0.604365,0.17596,-0.216668,-0.0990815,-0.348459,0.208961,0,0,0,0,0};
188  //double x[13]; for(int i=0;i<13;++i)x[i]=0.001;
189  //string a[13]={"a","b","c","d","e","f","g","h","i","j","k","l","flat"};
190  //cout << L.DoEval(x) << endl;
191  //return 0;
192 
193  // ---------------------------------------------------------------------------
194  // setup minimizer
195  printInfo << "Creating and setting up minimizer " << minimizerType[0] << " " << minimizerType[1] << endl;
196  Minimizer* minimizer = Factory::CreateMinimizer(minimizerType[0], minimizerType[1]);
197  if (!minimizer) {
198  printErr << "Error creating minimizer. Exiting." << endl;
199  throw;
200  }
201  minimizer->SetFunction(L);
202  minimizer->SetPrintLevel((quiet) ? 0 : 3);
203 
204  // ---------------------------------------------------------------------------
205  // read in object with start values
206  printInfo << "Reading start values from '" << startValFileName << "'." << endl;
207  double binCenter = 0.5 * (binLowMass + binHighMass);
208  bool hasStartVal = false;
209  TFile* startValFile = NULL;
210  TFitBin* startBin = NULL;
211  if (startValFileName.Length() > 2) {
212  // open root file
213  startValFile = TFile::Open(startValFileName, "READ");
214  if (!startValFile || startValFile->IsZombie())
215  printWarn << "Cannot open start value file '" << startValFileName << "'. Using default start values." << endl;
216  else {
217  // get tree with start values
218  TTree* tree;
219  startValFile->GetObject(valTreeName, tree);
220  if (!tree)
221  printWarn << "Cannot find start value tree '"<< valTreeName << "' in file '" << startValFileName << "'." << endl;
222  else {
223  startBin = new TFitBin();
224  tree->SetBranchAddress(valBranchName, &startBin);
225  // find entry which is closest to mass bin center
226  unsigned int iBest = 0;
227  double mBest = 0;
228  for (unsigned int i = 0; i < tree->GetEntriesFast(); ++i) {
229  tree->GetEntry(i);
230  if (fabs(binCenter - startBin->mass()) <= fabs(binCenter - mBest)) {
231  iBest = i;
232  mBest = startBin->mass();
233  }
234  } // end loop over TFitBins
235  tree->GetEntry(iBest);
236  hasStartVal = true;
237  }
238  }
239  } // endif startValFileName non-empty
240 
241  // ---------------------------------------------------------------------------
242  // set start parameter values
243  printInfo << "Setting start values for " << nmbPar << " parameters." << endl
244  << " Parameter naming scheme is: V[rank index]_[IGJPCME][isobar spec]" << endl;
245  unsigned int maxParNameLength = 0;
246  for (unsigned int i = 0; i < nmbPar; ++i){
247  if (L.parname(i).length() > maxParNameLength)
248  maxParNameLength = L.parname(i).length();
249  cout<<L.parname(i)<<endl;
250  }
251  bool success = true;
252  vector<bool> parIsFixed(nmbPar, false); // memorizes state of variables; ROOT::Math::Minimizer has no corresponding accessor
253  for (unsigned int i = 0; i < nmbPar; ++i) {
254  cout << i << ": " ;
255  cout.flush();
256  double startVal = defaultStartValue;
257  string parName = L.parname(i);
258  cout<<parName<< endl;
259  if (hasStartVal){
260  // look into start FitBin if we find a suitable parameter
261  assert(startBin);
262  startVal = startBin->getParameter(parName.c_str());
263  }
264  // check if parameter needs to be fixed because of threshold
265  if ((L.parthreshold(i) == 0) || (L.parthreshold(i) < binHighMass)) {
266  if (!minimizer->SetVariable(i, parName, startVal, defaultMinStep))
267  success = false;
268  if (startVal == 0) {
269  cout << " Read start value 0 for parameter " << parName << ". Setting default start value." << endl;
270  startVal = defaultStartValue;
271  }
272  cout << " Setting parameter [" << setw(3) << i << "] "
273  << setw(maxParNameLength) << parName << " = " << maxPrecisionAlign(startVal) << endl;
274  } else {
275  if (!minimizer->SetFixedVariable(i, parName, 0.)) // fix this parameter to 0
276  success = false;
277  cout << " Fixing parameter [" << setw(3) << i << "] "
278  << setw(maxParNameLength) << parName << " = 0" << endl;
279  parIsFixed[i] = true;
280  }
281  if (!success) {
282  printErr << "Something went wrong when setting minimizer parameters. Exiting." << endl;
283  throw;
284  }
285  }
286  // cleanup
287  if(startValFile) {
288  startValFile->Close();
289  delete startValFile;
290  startValFile = NULL;
291  }
292 
293  // ---------------------------------------------------------------------------
294  // find minimum of likelihood function
295  printInfo << "Performing minimization." << endl;
296  //L.SetMaxSampDL(1000);
297  // Pre-Minimizaton:
298  //minimizer->SetMaxIterations(10);
299  //minimizer->Minimize();
300  //L.SetMaxSampDL(20000);
301  // Minimizaton:
302  minimizer->SetMaxIterations(maxNmbOfIterations);
303  success = minimizer->Minimize();
304  if (!success)
305  printWarn << "Minimization failed for mass bin " << binCenter << endl;
306  if (runHesse) {
307  printInfo << "Calculating Hessian matrix." << endl;
308  //success = minimizer->Hesse(); // comes only with ROOT 5.24+
309  if (!success)
310  printWarn << "Calculation of Hessian matrix failed." << endl;
311  }
312  cout << "Minimization stopped after " << minimizer->NCalls() << " function calls." << endl
313  << " Total number of mimimzer parameters ................. " << minimizer->NDim() << endl
314  << " Number of free minimzer parameters .................. " << minimizer->NFree() << endl
315  << " Maximum allowed number of iterations ................ " << minimizer->MaxIterations() << endl
316  << " Maximum allowed number of function calls ............ " << minimizer->MaxFunctionCalls() << endl
317  << " Minimizer status .................................... " << minimizer->Status() << endl
318  << " Minimizer provides error and error matrix ........... " << minimizer->ProvidesError() << endl
319  << " Minimizer has performed detailed error validation ... " << minimizer->IsValidError() << endl
320  << " Estimated Distance to Minimum ....................... " << minimizer->Edm() << endl
321  << " Statistical scale used for error calculation ........ " << minimizer->ErrorDef() << endl
322  << " Minimizer strategy .................................. " << minimizer->Strategy() << endl
323  << " Absolute tolerance .................................. " << minimizer->Tolerance() << endl;
324 
325  // ---------------------------------------------------------------------------
326  // print results
327  printInfo << "Minimization result:" << endl;
328  //double x[nmbPar];
329  for (unsigned int i = 0; i< nmbPar; ++i) {
330  //x[i] = minimizer->X()[i]; // step 1% from minimum for derivative
331  cout << " Parameter [" << setw(3) << i << "] "
332  << setw(maxParNameLength) << L.parname(i) << " = ";
333  if (parIsFixed[i])
334  cout << minimizer->X()[i] << " (fixed)" << endl;
335  else {
336  cout << setw(12) << maxPrecisionAlign(minimizer->X()[i]) << " +- "
337  << setw(12) << maxPrecisionAlign(minimizer->Errors()[i]);
338  if (runMinos && (i == 156)) { // does not work for all parameters
339  double minosErrLow = 0;
340  double minosErrUp = 0;
341  success = minimizer->GetMinosError(i, minosErrLow, minosErrUp);
342  if (success)
343  cout << " Minos: " << "[" << minosErrLow << ", +" << minosErrUp << "]" << endl;
344  } else
345  cout << endl;
346  }
347  }
348  cout << "Number of calls to likelihood function FdF() ... " << L.ncalls() << endl
349  << "Total time spent for likelihood calculation .... " << L.Ltime() << endl
350  << "Total time spent for normalization ............. " << L.Ntime() << endl;
351 
352  // check derivatives:
353  // first numerical:
354  /*
355  double h=1E-8;
356  double dxNum[nmbPar];
357  double dxAna[nmbPar];
358  for(unsigned int i=0; i<nmbPar;++i){
359  x[i]+=h;
360  double L2=L.DoEval(x);
361  x[i]-=2*h;
362  double L1=L.DoEval(x);
363  dxNum[i]=0.5*(L2-L1)/h;
364  x[i]+=h;
365  }
366  double F;
367  L.FdF(x,F,dxAna);
368  for(unsigned int i=0; i<nmbPar;++i){
369  cout<< "dL/d"<<i<<"(num)="<<dxNum[i]<<endl;
370  cout<< "dL/d"<<i<<"(ana)="<<dxAna[i]<<endl;
371  }
372  */
373 
374  // ---------------------------------------------------------------------------
375  // write out result
376  printInfo << "Writing result to '" << outFileName << "'." << endl;
377  {
378  // open output file and create tree for writing
379  TFile* outFile = new TFile(outFileName, "UPDATE");
380  if (!outFile || outFile->IsZombie())
381  printWarn << "Cannot open output file '" << outFileName << "'. No results will be written." << endl;
382  else {
383  // check whether output tree already exists
384  TTree* tree;
385  TFitBin* result = new TFitBin();
386  outFile->GetObject(valTreeName, tree);
387  if (!tree) {
388  printInfo << "File '" << outFileName << "' is empty. Creating new tree for PWA result." << endl;
389  tree = new TTree(valTreeName, valTreeName);
390  tree->Branch(valBranchName, &result);
391  } else
392  tree->SetBranchAddress(valBranchName, &result);
393 
394  // get data structures to construct result TFitBin
395  vector<TComplex> prodAmplitudes; // production amplitudes
396  vector<pair<int,int> > indices; // indices for error matrix access
397  vector<TString> waveNames; // contains rank information
398  vector<complex<double> > V;
399  vector<string> names;
400 
401  // error matrix of fit parameters
402  TMatrixD errMatrix(nmbPar, nmbPar);
403  for(unsigned int i = 0; i < nmbPar; ++i){
404  for(unsigned int j = 0; j < nmbPar; ++j) {
405  errMatrix[i][j] = minimizer->CovMatrix(i,j);
406  //if(i==j)cout << "Sig"<< i << "=" << sqrt(errMatrix[i][i]) << endl;
407  }
408  }
409  TMatrixD covMatrix(2,2); // covariants of amplitudes
410  // These two matrices can differ in the case when there are
411  // constraints present such that the fit parameters do not
412  // coincide with amplitude Re/Im anymore.
413 
414  L.buildCAmps(minimizer->X(), V, indices, names, errMatrix,covMatrix,true);
415 
416  // convert to TComplex;
417  for (unsigned int i = 0; i < V.size(); ++i)
418  prodAmplitudes.push_back(TComplex(V[i].real(), V[i].imag()));
419  assert(indices.size() == prodAmplitudes.size());
420  for(unsigned int i = 0; i < names.size(); ++i)
421  waveNames.push_back(TString(names[i].c_str()));
422 
423  vector<TString> waveTitles; // without rank
424  {
425  const vector<string>& titles = L.wavetitles();
426  for(unsigned int i = 0; i < titles.size(); ++i)
427  waveTitles.push_back(TString(titles[i].c_str()));
428  waveTitles.push_back("flat");
429  }
430 
431  // // normalization integral
432 // integral norm = L.normInt();
433  // integral and acceptance Matrix
434  cout << " Setup integrals" << endl;
435  unsigned int n = waveTitles.size();
436  TCMatrix integralMatrix(n, n);
437  TCMatrix accMatrix(n, n);
438 
439 
440  L.getIntCMatrix(integralMatrix, accMatrix);
441  //integralMatrix.Print();
442  // representation of number of events depends on whether normalization was done
443  int nmbEvt = useNorm ? 1 : L.nevents();
444 
445  cout << "Filling TFitBin:" << endl;
446  cout << " Number of fit parameters ........... " << nmbPar << endl;
447  cout << " Number of production amplitudes .... " << prodAmplitudes.size() << endl;
448  cout << " Number of indices .................. " << indices.size() << endl;
449  cout << " Number of wave names (with rank) ... " << waveNames.size() << endl;
450  cout << " Number of wave titles (w/o rank) ... " << waveTitles.size() << endl;
451  cout << " Dimension of error matrix .......... " << errMatrix.GetNrows() << endl;
452  cout << " Dimension of integral matrix ....... " << integralMatrix.nrows() << endl;
453  cout << " Dimension of acceptance matrix ..... " << accMatrix.nrows() << endl;
454  result->fill(prodAmplitudes,
455  indices,
456  waveNames,
457  nmbEvt,
458  L.nevents(), // raw number of data events
459  binCenter,
460  integralMatrix,
461  covMatrix,
462  minimizer->MinValue(),
463  rank);
464  //result->PrintParameters();
465 
466  // write result to file
467  TString binName = valBranchName;
468  binName += binLowMass;
469  binName += "_";
470  binName += binHighMass;
471  binName.ReplaceAll(" ","");
472  tree->Fill();
473  tree->Write("", TObject::kOverwrite);
474  outFile->Close();
475  }
476  }
477 
478  if (minimizer)
479  delete minimizer;
480  return 0;
481 }
482 
483 
484 // dummy function needed since we link to but do not use minuit.o
485 int mnparm(int, string, double, double, double, double)
486 {
487  printErr << "this is impossible" << endl;
488  throw "aFit";
489  return 0;
490 }
Generated by   doxygen 1.8.1.2