ROOTPWA
fitResult.h
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 // Data storage class for PWA fit result of one kinematic bin
27 //
28 // Environment:
29 // Software developed for the COMPASS experiment at CERN
30 //
31 // Author List:
32 // Sebastian Neubert TUM (original author)
33 //
34 //
35 //-----------------------------------------------------------
36 
37 
38 #ifndef FITRESULT_H
39 #define FITRESULT_H
40 
41 
42 #include <iostream>
43 #include <iomanip>
44 #include <vector>
45 #include <map>
46 #include <string>
47 #include <complex>
48 
49 #include "TObject.h"
50 #include "TComplex.h"
51 #include "TMatrixT.h"
52 #include "TString.h"
53 #include "TRegexp.h"
54 
55 #include "reportingUtils.hpp"
56 #include "reportingUtilsRoot.hpp"
57 #include "TCMatrix.h"
58 #include "TFitBin.h"
59 
60 // enable copying from TFitResult for older ROOT versions
61 #include "TFitResult.h"
62 
63 
64 namespace rpwa {
65 
66  inline
67  TString
68  escapeRegExpSpecialChar(const std::string& s)
69  {
70  TString escapedS(s);
71  const char specialChars[] = {'^', '$', '.', '[', ']', '^', '*', '+', '?'};
72  for (unsigned int i = 0; i < sizeof(specialChars) / sizeof(specialChars[0]); ++i) {
73  TString escapedChar = TString("\\").Append(specialChars[i]);
74  escapedS.ReplaceAll(specialChars[i], escapedChar);
75  }
76  return escapedS;
77  }
78 
79 
80  inline
81  TString
82  unescapeRegExpSpecialChar(const std::string& s)
83  {
84  TString escapedS(s);
85  const char specialChars[] = {'^', '$', '.', '[', ']', '^', '*', '+', '?'};
86  for (unsigned int i = 0; i < sizeof(specialChars) / sizeof(specialChars[0]); ++i) {
87  TString escapedChar = TString("\\").Append(specialChars[i]);
88  escapedS.ReplaceAll(escapedChar, specialChars[i]);
89  }
90  return escapedS;
91  }
92 
93 
95  class fitResult : public TObject {
96 
97  public:
98 
99  fitResult();
100  fitResult(const fitResult& result);
101  fitResult(const TFitBin& fitBin);
102  // enable copying from TFitResult for older ROOT versions
103 #ifdef USE_TFITRESULT
104  fitResult(const TFitResult& result);
105 #endif
106  virtual ~fitResult();
107 
108  // create a copy of the current object with slightly modified
109  // parameters according to the distribution described by the
110  // covariance matrix
111  fitResult* cloneVar();
112 
113  void reset();
114  void fill(const unsigned int nmbEvents, // number of events in bin
115  const unsigned int normNmbEvents, // number of events to normalize to
116  const double massBinCenter, // center value of mass bin
117  const double logLikelihood, // log(likelihood) at maximum
118  const int rank, // rank of fit
119  const std::vector<std::complex<double> >& prodAmps, // production amplitudes
120  const std::vector<std::string>& prodAmpNames, // names of production amplitudes used in fit
121  const TMatrixT<double>& fitParCovMatrix, // covariance matrix of fit parameters
122  const std::vector<std::pair<int, int> >& fitParCovMatrixIndices, // indices of fit parameters for real and imaginary part in covariance matrix matrix
123  const TCMatrix& normIntegral, // normalization integral matrix
124  const std::vector<double>& phaseSpaceIntegral, // normalization integral over full phase space without acceptance
125  const bool converged, // indicates whether fit has converged (according to minimizer)
126  const bool hasHessian); // indicates whether Hessian matrix has been calculated successfully
127 
128  double massBinCenter() const { return _massBinCenter; }
129  double logLikelihood() const { return _logLikelihood; }
130  double evidence () const;
131  bool converged () const { return _converged; }
132  bool hasHessian () const { return _hasHessian; }
133  unsigned int rank () const { return _rank; }
134  unsigned int nmbEvents () const { return _nmbEvents; }
135  unsigned int normNmbEvents() const { return _normNmbEvents; }
136  unsigned int nmbWaves () const { return _waveNames.size(); }
137  unsigned int nmbProdAmps () const { return _prodAmps.size(); }
138 
139  TString waveName (const unsigned int waveIdx) const { return _waveNames[waveIdx]; }
140  TString waveNameEsc (const unsigned int waveIdx) const { return escapeRegExpSpecialChar(_waveNames[waveIdx]); }
141  TString prodAmpName (const unsigned int prodAmpIdx) const { return _prodAmpNames[prodAmpIdx]; }
142  TString prodAmpNameEsc(const unsigned int prodAmpIdx) const { return escapeRegExpSpecialChar(_prodAmpNames[prodAmpIdx]); }
143 
144  int waveIndex (const std::string& waveName ) const;
145  int prodAmpIndex(const std::string& prodAmpName) const;
146 
147  double fitParameter (const std::string& parName) const;
148 
149  double fitParameterCov(const unsigned int parIndexA,
150  const unsigned int parIndexB) const { return _fitParCovMatrix[parIndexA][parIndexB]; }
151  inline void fitParameters (double* parArray) const;
152 
154  std::complex<double> prodAmp (const unsigned int prodAmpIdx) const
155  { return std::complex<double>(_prodAmps[prodAmpIdx].Re(), _prodAmps[prodAmpIdx].Im()); }
156  inline TMatrixT<double> prodAmpCov(const unsigned int prodAmpIdx) const;
157 
158  TMatrixT<double> prodAmpCov(const std::vector<unsigned int>& prodAmpIndices) const;
160  inline TMatrixT<double> prodAmpCov(const std::vector<std::pair<unsigned int, unsigned int> >& prodAmpIndexPairs) const;
162  inline TMatrixT<double> prodAmpCov(const std::vector<unsigned int>& prodAmpIndicesA,
163  const std::vector<unsigned int>& prodAmpIndicesB) const;
164  bool covMatrixValid() const { return _covMatrixValid; }
165 
167  inline std::complex<double> normIntegral(const unsigned int waveIndexA,
168  const unsigned int waveIndexB) const;
170  double phaseSpaceIntegral(const unsigned int waveIdx) const { return _phaseSpaceIntegral[waveIdx]; }
171  double phaseSpaceIntegral(const std::string& waveName) const { return _phaseSpaceIntegral[waveIndex(waveName)]; }
172 
174  std::complex<double> spinDensityMatrixElem (const unsigned int waveIndexA,
175  const unsigned int waveIndexB) const;
177  TMatrixT<double> spinDensityMatrixElemCov(const unsigned int waveIndexA,
178  const unsigned int waveIndexB) const;
179 
181  double intensity (const unsigned int waveIdx) const { return spinDensityMatrixElem(waveIdx, waveIdx).real(); }
183  double intensityErr(const unsigned int waveIdx) const { return sqrt(spinDensityMatrixElemCov(waveIdx, waveIdx)[0][0]); }
184  double intensity (const char* waveNamePattern) const;
185  double intensityErr(const char* waveNamePattern) const;
186  double intensity () const { return intensity (".*"); }
187  double intensityErr() const { return intensityErr(".*"); }
188 
189  double phase (const unsigned int waveIndexA,
190  const unsigned int waveIndexB) const;
191  double phaseErr (const unsigned int waveIndexA,
192  const unsigned int waveIndexB) const;
193  double phase (const std::string waveNameA,
194  const std::string waveNameB) const
195  { return phase(waveIndex(waveNameA), waveIndex(waveNameB)); }
196  double phaseErr (const std::string waveNameA,
197  const std::string waveNameB) const
198  { return phaseErr(waveIndex(waveNameA), waveIndex(waveNameB)); }
199 
200  double coherence (const unsigned int waveIndexA,
201  const unsigned int waveIndexB) const;
202  double coherenceErr(const unsigned int waveIndexA,
203  const unsigned int waveIndexB) const;
204  double overlap (const unsigned int waveIndexA,
205  const unsigned int waveIndexB) const;
206  double overlapErr (const unsigned int waveIndexA,
207  const unsigned int waveIndexB) const;
208 
209  // low level interface to make copying easier
210  const std::vector<TComplex>& prodAmps () const { return _prodAmps; }
211  const std::vector<std::string>& prodAmpNames () const { return _prodAmpNames; }
212  const std::vector<std::string>& waveNames () const { return _waveNames; }
213  const TMatrixT<Double_t>& fitParCovMatrix () const { return _fitParCovMatrix; }
214  const std::vector<std::pair<Int_t, Int_t> >& fitParCovIndices () const { return _fitParCovMatrixIndices; }
215  const TCMatrix& normIntegralMatrix() const { return _normIntegral; }
216  const std::map<Int_t, Int_t>& normIntIndexMap () const { return _normIntIndexMap; }
217 
218  inline std::ostream& printProdAmpNames(std::ostream& out = std::cout) const;
219  inline std::ostream& printWaveNames (std::ostream& out = std::cout) const;
220  inline std::ostream& printProdAmps (std::ostream& out = std::cout) const;
221  inline std::ostream& printWaves (std::ostream& out = std::cout) const;
222  inline std::ostream& printAmpsGenPW (std::ostream& out = std::cout) const;
223 
224  virtual inline std::ostream& print(std::ostream& out = std::cout) const;
225  friend std::ostream& operator << (std::ostream& out,
226  const fitResult& result) { return result.print(out); }
227 
228  private:
229 
230  UInt_t _nmbEvents;
231  UInt_t _normNmbEvents;
232  Double_t _massBinCenter;
233  Double_t _logLikelihood;
234  Int_t _rank;
235  std::vector<TComplex> _prodAmps;
236  std::vector<std::string> _prodAmpNames;
237  std::vector<std::string> _waveNames;
238  Bool_t _covMatrixValid;
239  TMatrixT<Double_t> _fitParCovMatrix;
240  std::vector<std::pair<Int_t, Int_t> > _fitParCovMatrixIndices;
241  TCMatrix _normIntegral;
242  std::map<Int_t, Int_t> _normIntIndexMap;
243  std::vector<double> _phaseSpaceIntegral;
244  bool _converged;
245  bool _hasHessian;
246  // add more info about fit: quality of fit information, ndf, list of fixed parameters, ...
247 
248  // helper functions
249  inline static TMatrixT<double> matrixRepr(const std::complex<double>& c);
250 
251  inline int rankOfProdAmp(const unsigned int prodAmpIndex) const;
252 
253  std::complex<double> normIntegralForProdAmp(const unsigned int prodAmpIndexA,
254  const unsigned int prodAmpIndexB) const;
255 
256  inline std::vector<unsigned int> waveIndicesMatchingPattern (const std::string& waveNamePattern ) const;
257  inline std::vector<unsigned int> prodAmpIndicesMatchingPattern(const std::string& prodAmpNamePattern) const;
258  std::vector<unsigned int> prodAmpIndicesForWave (const unsigned int waveIdx ) const
259  { return prodAmpIndicesMatchingPattern(waveNameEsc(waveIdx).Data()); }
260  inline std::vector<std::pair<unsigned int, unsigned int> > prodAmpIndexPairsForWaves
261  (const unsigned int waveIndexA,
262  const unsigned int waveIndexB) const;
263 
264  inline double realValVariance(const unsigned int waveIndexA,
265  const unsigned int waveIndexB,
266  const TMatrixT<double>& jacobian) const;
267 
268  TString wavetitle(int i) const
269  {
270  const TString ampName = _prodAmpNames[i];
271  if (ampName.Contains("V_"))
272  return ampName(2, 1000);
273  else
274  return ampName(3, 1000);
275  }
276  void buildWaveMap();
277 
278  public:
279 
280  ClassDef(fitResult,4)
281 
282  }; // class fitResult
283 
284 
286  inline
287  void
288  fitResult::fitParameters(double* parArray) const
289  {
290  unsigned int countPar = 0;
291  for (unsigned int i = 0; i < nmbProdAmps(); ++i) {
292  parArray[countPar++] = prodAmp(i).real();
293  // check if amplitude has imaginary part
294  if (_fitParCovMatrixIndices[i].second >= 0)
295  parArray[countPar++] = prodAmp(i).imag();
296  }
297  }
298 
299 
301  inline
302  TMatrixT<double>
303  fitResult::prodAmpCov(const unsigned int prodAmpIdx) const {
304  // get parameter indices
305  const int i = _fitParCovMatrixIndices[prodAmpIdx].first;
306  const int j = _fitParCovMatrixIndices[prodAmpIdx].second;
307  TMatrixT<double> cov(2, 2);
308  cov[0][0] = _fitParCovMatrix[i][i];
309  if (j >= 0) {
310  cov[0][1] = _fitParCovMatrix[i][j];
311  cov[1][0] = _fitParCovMatrix[j][i];
312  cov[1][1] = _fitParCovMatrix[j][j];
313  }
314  return cov;
315  }
316 
317 
323  inline
324  TMatrixT<double>
325  fitResult::prodAmpCov(const std::vector<std::pair<unsigned int, unsigned int> >& prodAmpIndexPairs) const
326  {
327  std::vector<unsigned int> prodAmpIndices;
328  // copy first production amplitude indices
329  for (unsigned int i = 0; i < prodAmpIndexPairs.size(); ++i)
330  prodAmpIndices.push_back(prodAmpIndexPairs[i].first);
331  // copy second production amplitude indices
332  for (unsigned int i = 0; i < prodAmpIndexPairs.size(); ++i)
333  prodAmpIndices.push_back(prodAmpIndexPairs[i].second);
334  return prodAmpCov(prodAmpIndices);
335  }
336 
337 
343  inline
344  TMatrixT<double>
345  fitResult::prodAmpCov(const std::vector<unsigned int>& prodAmpIndicesA,
346  const std::vector<unsigned int>& prodAmpIndicesB) const
347  {
348  std::vector<unsigned int> prodAmpIndices;
349  // copy wave A production amplitude indices
350  prodAmpIndices.assign(prodAmpIndicesA.begin(), prodAmpIndicesA.end());
351  // copy wave B production amplitude indices
352  prodAmpIndices.insert(prodAmpIndices.end(), prodAmpIndicesB.begin(), prodAmpIndicesB.end());
353  return prodAmpCov(prodAmpIndices);
354  }
355 
356 
357  // returns normalization integral for pair of waves at index A and B
358  inline
359  std::complex<double>
360  fitResult::normIntegral(const unsigned int waveIndexA,
361  const unsigned int waveIndexB) const
362  {
363  const TComplex norm = _normIntegral(waveIndexA, waveIndexB);
364  return std::complex<double>(norm.Re(), norm.Im());
365  }
366 
367 
368  // prints all production amplitude names
369  inline
370  std::ostream&
371  fitResult::printProdAmpNames(std::ostream& out) const
372  {
373  out << " Production amplitude names:" << std::endl;
374  for (unsigned int i = 0; i < nmbProdAmps(); ++i)
375  out << " " << std::setw(3) << i << " " << _prodAmpNames[i] << std::endl;
376  return out;
377  }
378 
379 
380  // prints all wave names
381  inline
382  std::ostream&
383  fitResult::printWaveNames(std::ostream& out) const
384  {
385  out << " Wave names:" << std::endl;
386  for (unsigned int i = 0; i < nmbWaves(); ++i)
387  out << " " << std::setw(3) << i << " " << _waveNames[i] << std::endl;
388  return out;
389  }
390 
391 
392  // prints all production amplitudes and their covariance matrix
393  inline
394  std::ostream&
395  fitResult::printProdAmps(std::ostream& out) const
396  {
397  out << "Production amplitudes:" << std::endl;
398  for (unsigned int i = 0; i < nmbProdAmps(); ++i) {
399  out << " " << std::setw(3) << i << " " << _prodAmpNames[i] << " = " << prodAmp(i)
400  << ", cov = " << prodAmpCov(i) << std::endl;
401  }
402  return out;
403  }
404 
405 
406  // prints all wave intensities and their errors
407  inline
408  std::ostream&
409  fitResult::printWaves(std::ostream& out) const
410  {
411  out << "Waves:" << std::endl;
412  for (unsigned int i = 0; i < nmbWaves(); ++i) {
413  out << " " << std::setw(3) << i << " " << _waveNames[i] << " = " << intensity(i)
414  << " +- " << intensityErr(i) << std::endl;
415  }
416  return out;
417  }
418 
419 
420  // prints all amplitudes in format used for genpw
421  // only supports rank 1 at the moment!!!
422  inline
423  std::ostream&
424  fitResult::printAmpsGenPW(std::ostream& s) const {
425  for(unsigned int i=0;i<_waveNames.size();++i){
426  s << _waveNames[i] << " "
427  << prodAmp(i).real() << " "
428  << prodAmp(i).imag() << std::endl;
429  }
430  return s;
431  }
432 
433 
435  inline
436  std::ostream&
437  fitResult::print(std::ostream& out) const
438  {
439  out << "fitResult dump:" << std::endl
440  << " number of events .................... " << _nmbEvents << std::endl
441  << " number of events to normalize to .... " << _normNmbEvents << std::endl
442  << " center value of mass bin ............ " << _massBinCenter << std::endl
443  << " log(likelihood) at maximum .......... " << _logLikelihood << std::endl
444  << " rank of fit ......................... " << _rank << std::endl
445  << " bin has a valid covariance matrix ... " << _covMatrixValid << std::endl;
446  printProdAmps(out);
447  printWaveNames(out);
448  out << " covariance matrix:" << std::endl << _fitParCovMatrix << std::endl;
449  out << " covariance matrix indices:" << std::endl;
450  for (unsigned int i = 0; i < _fitParCovMatrixIndices.size(); ++i)
451  out << " index " << std::setw(3) << i << " = (" << std::setw(3) << _fitParCovMatrixIndices[i].first
452  << ", " << std::setw(3) << _fitParCovMatrixIndices[i].second << ")" << std::endl;
453  out << " normalization integral (w/o acceptance):" << std::endl << _normIntegral << std::endl;
454  out << " map of production amplitude indices to indices in normalization integral:" << std::endl;
455  for (std::map<Int_t, Int_t>::const_iterator i = _normIntIndexMap.begin(); i != _normIntIndexMap.end(); ++i)
456  out << " prod. amp [" << std::setw(3) << i->first << "] "
457  << "-> norm. int. [" << std::setw(3) << i->second << "]" << std::endl;
458  return out;
459  }
460 
461 
466  // !!! possible optimization: return TMatrixTLazy
467  inline
468  TMatrixT<double>
469  fitResult::matrixRepr(const std::complex<double>& c)
470  {
471  TMatrixT<double> m(2, 2);
472  m[0][0] = m[1][1] = c.real();
473  m[0][1] = -c.imag();
474  m[1][0] = c.imag();
475  return m;
476  }
477 
478 
480  inline
481  int
482  fitResult::rankOfProdAmp(const unsigned int prodAmpIdx) const
483  {
484  const TString ampName = _prodAmpNames[prodAmpIdx];
485  if (ampName.Contains("flat"))
486  return -1;
487  else
488  // the rank is encoded in second character of parameter name
489  return TString(ampName(1, 1)).Atoi();
490  }
491 
492 
494  inline
495  std::vector<unsigned int>
496  fitResult::waveIndicesMatchingPattern(const std::string& waveNamePattern) const
497  {
498  // escape special characters:
499  TString Pattern(waveNamePattern);
500  Pattern.ReplaceAll("+","\\+");
501  Pattern.ReplaceAll("\\\\+","\\+");
502  std::vector<unsigned int> waveIndices;
503  for (unsigned int waveIdx = 0; waveIdx < nmbWaves(); ++waveIdx){
504  //std::cout<<waveName(waveIdx)<<std::endl;
505  if (waveName(waveIdx).Contains(TRegexp(Pattern)))
506  waveIndices.push_back(waveIdx);
507  }
508  return waveIndices;
509  }
510 
511 
513  inline
514  std::vector<unsigned int>
515  fitResult::prodAmpIndicesMatchingPattern(const std::string& ampNamePattern) const
516  {
517  // escape special characters:
518  TString Pattern(ampNamePattern);
519  Pattern.ReplaceAll("+","\\+");
520  Pattern.ReplaceAll("\\\\+","\\+");
521  std::vector<unsigned int> prodAmpIndices;
522  for (unsigned int prodAmpIdx = 0; prodAmpIdx < nmbProdAmps(); ++prodAmpIdx)
523  if (prodAmpName(prodAmpIdx).Contains(TRegexp(Pattern)))
524  prodAmpIndices.push_back(prodAmpIdx);
525  return prodAmpIndices;
526  }
527 
528 
530  inline
531  std::vector<std::pair<unsigned int, unsigned int> >
532  fitResult::prodAmpIndexPairsForWaves(const unsigned int waveIndexA,
533  const unsigned int waveIndexB) const
534  {
535  std::vector<std::pair<unsigned int, unsigned int> > prodAmpIndexPairs;
536  const std::vector<unsigned int> prodAmpIndicesA = prodAmpIndicesForWave(waveIndexA);
537  const std::vector<unsigned int> prodAmpIndicesB = prodAmpIndicesForWave(waveIndexB);
538  for (unsigned int countAmpA = 0; countAmpA < prodAmpIndicesA.size(); ++countAmpA) {
539  const unsigned int ampIndexA = prodAmpIndicesA[countAmpA];
540  const int ampRankA = rankOfProdAmp(ampIndexA);
541  // find production amplitude of wave B with same rank
542  int ampIndexB = -1;
543  for (unsigned int countAmpB = 0; countAmpB < prodAmpIndicesB.size(); ++countAmpB)
544  if (rankOfProdAmp(prodAmpIndicesB[countAmpB]) == ampRankA) {
545  ampIndexB = prodAmpIndicesB[countAmpB];
546  break;
547  }
548  if (ampIndexB >=0)
549  prodAmpIndexPairs.push_back(std::make_pair(ampIndexA, (unsigned int)ampIndexB));
550  }
551  return prodAmpIndexPairs;
552  }
553 
554 
556  inline
557  double
558  fitResult::realValVariance(const unsigned int waveIndexA,
559  const unsigned int waveIndexB,
560  const TMatrixT<double>& jacobian) const // Jacobian of real valued function (d f/ d Re[rho] d f / d Im[rho])
561  {
562  if (!_covMatrixValid)
563  return 0;
564  const TMatrixT<double> spinDensCov = spinDensityMatrixElemCov(waveIndexA, waveIndexB); // 2 x 2 matrix
565  const TMatrixT<double> jacobianT(TMatrixT<double>::kTransposed, jacobian); // 2 x 1 matrix
566  const TMatrixT<double> spinDensCovJT = spinDensCov * jacobianT; // 2 x 1 matrix
567  const TMatrixT<double> cov = jacobian * spinDensCovJT; // 1 x 1 matrix
568  return cov[0][0];
569  }
570 
571 
572 } // namespace rpwa
573 
574 
575 #endif // FITRESULT_H
Generated by   doxygen 1.8.1.2