ROOTPWA
TFitResult.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 //
39 // !!! deprecated class --- use fitResult instead !!!
40 //
41 
42 
43 #ifndef TFITRESULT_HH
44 #define TFITRESULT_HH
45 
46 
47 #ifdef USE_TFITRESULT
48 
49 
50 #include <iostream>
51 #include <iomanip>
52 #include <vector>
53 #include <map>
54 #include <string>
55 #include <complex>
56 
57 #include "TObject.h"
58 #include "TComplex.h"
59 #include "TMatrixT.h"
60 #include "TString.h"
61 
62 #include "reportingUtils.hpp"
63 #include "reportingUtilsRoot.hpp"
64 #include "TCMatrix.h"
65 #include "TFitBin.h"
66 
67 
69 class TFitResult : public TObject {
70 
71 public:
72 
73  TFitResult();
74  TFitResult(const TFitResult& result);
75  TFitResult(const TFitBin& fitBin);
76  virtual ~TFitResult();
77 
78  void reset();
79  void fill(const unsigned int nmbEvents, // number of events in bin
80  const unsigned int normNmbEvents, // number of events to normalize to
81  const double massBinCenter, // center value of mass bin
82  const double logLikelihood, // log(likelihood) at maximum
83  const int rank, // rank of fit
84  const std::vector<std::complex<double> >& prodAmps, // production amplitudes
85  const std::vector<std::string>& prodAmpNames, // names of production amplitudes used in fit
86  const TMatrixT<double>& fitParCovMatrix, // covariance matrix of fit parameters
87  const std::vector<std::pair<int, int> >& fitParCovMatrixIndices, // indices of fit parameters for real and imaginary part in covariance matrix matrix
88  const TCMatrix& normIntegral); // normalization integral over full phase space without acceptance
89 
90  double massBinCenter() const { return _massBinCenter; }
91  double logLikelihood() const { return _logLikelihood; }
92  double evidence () const ;
93  unsigned int rank () const { return _rank; }
94  unsigned int nmbEvents () const { return _nmbEvents; }
95  unsigned int normNmbEvents() const { return _normNmbEvents; }
96  unsigned int nmbWaves () const { return _waveNames.size(); }
97  unsigned int nmbProdAmps () const { return _prodAmps.size(); }
98 
99  TString waveName (const unsigned int waveIndex) const { return _waveNames[waveIndex]; }
100  TString prodAmpName (const unsigned int prodAmpIndex) const { return _prodAmpNames[prodAmpIndex]; }
101 
102  double fitParameter (const std::string& parName) const;
103 
104  double fitParameterCov(const unsigned int parIndexA,
105  const unsigned int parIndexB) const { return _fitParCovMatrix[parIndexA][parIndexB]; }
106  inline void fitParameters (double* parArray) const;
107 
109  std::complex<double> prodAmp (const unsigned int prodAmpIndex) const { return std::complex<double>(_prodAmps[prodAmpIndex].Re(), _prodAmps[prodAmpIndex].Im()); }
110  inline TMatrixT<double> prodAmpCov(const unsigned int prodAmpIndex) const;
111 
112  TMatrixT<double> prodAmpCov(const std::vector<unsigned int>& prodAmpIndices) const;
114  inline TMatrixT<double> prodAmpCov(const std::vector<std::pair<unsigned int, unsigned int> >& prodAmpIndexPairs) const;
116  inline TMatrixT<double> prodAmpCov(const std::vector<unsigned int>& prodAmpIndicesA,
117  const std::vector<unsigned int>& prodAmpIndicesB) const;
118  bool covMatrixValid() const { return _hasErrors; }
119 
121  inline std::complex<double> normIntegral(const unsigned int waveIndexA,
122  const unsigned int waveIndexB) const;
123 
125  std::complex<double> spinDensityMatrixElem (const unsigned int waveIndexA,
126  const unsigned int waveIndexB) const;
128  TMatrixT<double> spinDensityMatrixElemCov(const unsigned int waveIndexA,
129  const unsigned int waveIndexB) const;
130 
132  double intensity (const unsigned int waveIndex) const { return spinDensityMatrixElem(waveIndex, waveIndex).real(); }
134  double intensityErr(const unsigned int waveIndex) const { return sqrt(spinDensityMatrixElemCov(waveIndex, waveIndex)[0][0]); }
135  double intensity (const char* waveNamePattern) const;
136  double intensityErr(const char* waveNamePattern) const;
137  double intensity () const { return intensity(""); }
138  double intensityErr() const { return intensityErr(""); }
139 
140  double phase (const unsigned int waveIndexA,
141  const unsigned int waveIndexB) const;
142  double phaseErr (const unsigned int waveIndexA,
143  const unsigned int waveIndexB) const;
144  double coherence (const unsigned int waveIndexA,
145  const unsigned int waveIndexB) const;
146  double coherenceErr(const unsigned int waveIndexA,
147  const unsigned int waveIndexB) const;
148  double overlap (const unsigned int waveIndexA,
149  const unsigned int waveIndexB) const;
150  double overlapErr (const unsigned int waveIndexA,
151  const unsigned int waveIndexB) const;
152 
153  // low level interface to make copying easier
154  const std::vector<TComplex>& prodAmps () const { return _prodAmps; }
155  const std::vector<std::string>& prodAmpNames () const { return _prodAmpNames; }
156  const std::vector<std::string>& waveNames () const { return _waveNames; }
157  const TMatrixT<Double_t>& fitParCovMatrix () const { return _fitParCovMatrix; }
158  const std::vector<std::pair<Int_t, Int_t> >& fitParCovIndices () const { return _fitParCovMatrixIndices; }
159  const TCMatrix& normIntegralMatrix() const { return _normIntegral; }
160  const std::map<Int_t, Int_t>& normIntIndexMap () const { return _normIntIndexMap; }
161 
162  inline std::ostream& printProdAmpNames(std::ostream& out = std::cout) const;
163  inline std::ostream& printWaveNames (std::ostream& out = std::cout) const;
164  inline std::ostream& printProdAmps (std::ostream& out = std::cout) const;
165  inline std::ostream& printWaves (std::ostream& out = std::cout) const;
166 
167  virtual inline std::ostream& print(std::ostream& out = std::cout) const;
168  friend std::ostream& operator << (std::ostream& out,
169  const TFitResult& fitResult) { return fitResult.print(out); }
170 
171 private:
172 
173  UInt_t _nmbEvents;
174  UInt_t _normNmbEvents;
175  Double_t _massBinCenter;
176  Double_t _logLikelihood;
177  Int_t _rank;
178  std::vector<TComplex> _prodAmps;
179  std::vector<std::string> _prodAmpNames;
180  std::vector<std::string> _waveNames;
181  Bool_t _hasErrors;
182  TMatrixT<Double_t> _fitParCovMatrix;
183  std::vector<std::pair<Int_t, Int_t> > _fitParCovMatrixIndices;
184  TCMatrix _normIntegral;
185  std::map<Int_t, Int_t> _normIntIndexMap;
186 
187  // add more info about fit: quality of fit information, ndf, list of fixed parameters, ...
188 
189  // helper functions
190  inline static TMatrixT<double> matrixRepr(const std::complex<double>& c);
191 
192  inline int rankOfProdAmp(const unsigned int prodAmpIndex) const;
193 
194  std::complex<double> normIntegralForProdAmp(const unsigned int prodAmpIndexA,
195  const unsigned int prodAmpIndexB) const;
196 
197  inline std::vector<unsigned int> waveIndicesMatchingPattern (const std::string& waveNamePattern) const;
198  inline std::vector<unsigned int> prodAmpIndicesMatchingPattern(const std::string& ampNamePattern ) const;
199  std::vector<unsigned int> prodAmpIndicesForWave (const unsigned int waveIndex) const
200  { return prodAmpIndicesMatchingPattern(_waveNames[waveIndex]); }
201  inline std::vector<std::pair<unsigned int, unsigned int> > prodAmpIndexPairsForWaves(const unsigned int waveIndexA,
202  const unsigned int waveIndexB) const;
203 
204  inline double realValVariance(const unsigned int waveIndexA,
205  const unsigned int waveIndexB,
206  const TMatrixT<double>& jacobian) const;
207 
208 
209  TString wavetitle(int i) const
210  {
211  const TString ampName = _prodAmpNames[i];
212  if (ampName.Contains("V_"))
213  return ampName(2, 1000);
214  else
215  return ampName(3, 1000);
216  }
217  void buildWaveMap();
218 
219 public:
220 
221  ClassDef(TFitResult,1)
222 
223 };
224 
225 
227 inline
228 void
229 TFitResult::fitParameters(double* parArray) const
230 {
231  unsigned int countPar = 0;
232  for (unsigned int i = 0; i < nmbProdAmps(); ++i) {
233  parArray[countPar++] = prodAmp(i).real();
234  // check if amplitude has imaginary part
235  if (_fitParCovMatrixIndices[i].second >= 0)
236  parArray[countPar++] = prodAmp(i).imag();
237  }
238 }
239 
240 
242 inline
243 TMatrixT<double>
244 TFitResult::prodAmpCov(const unsigned int prodAmpIndex) const {
245  // get parameter indices
246  const int i = _fitParCovMatrixIndices[prodAmpIndex].first;
247  const int j = _fitParCovMatrixIndices[prodAmpIndex].second;
248  TMatrixT<double> cov(2, 2);
249  cov[0][0] = _fitParCovMatrix[i][i];
250  if (j >= 0) {
251  cov[0][1] = _fitParCovMatrix[i][j];
252  cov[1][0] = _fitParCovMatrix[j][i];
253  cov[1][1] = _fitParCovMatrix[j][j];
254  }
255  return cov;
256 }
257 
258 
264 inline
265 TMatrixT<double>
266 TFitResult::prodAmpCov(const std::vector<std::pair<unsigned int, unsigned int> >& prodAmpIndexPairs) const
267 {
268  std::vector<unsigned int> prodAmpIndices;
269  // copy first production amplitude indices
270  for (unsigned int i = 0; i < prodAmpIndexPairs.size(); ++i)
271  prodAmpIndices.push_back(prodAmpIndexPairs[i].first);
272  // copy second production amplitude indices
273  for (unsigned int i = 0; i < prodAmpIndexPairs.size(); ++i)
274  prodAmpIndices.push_back(prodAmpIndexPairs[i].second);
275  return prodAmpCov(prodAmpIndices);
276 }
277 
278 
284 inline
285 TMatrixT<double>
286 TFitResult::prodAmpCov(const std::vector<unsigned int>& prodAmpIndicesA,
287  const std::vector<unsigned int>& prodAmpIndicesB) const
288 {
289  std::vector<unsigned int> prodAmpIndices;
290  // copy wave A production amplitude indices
291  prodAmpIndices.assign(prodAmpIndicesA.begin(), prodAmpIndicesA.end());
292  // copy wave B production amplitude indices
293  prodAmpIndices.insert(prodAmpIndices.end(), prodAmpIndicesB.begin(), prodAmpIndicesB.end());
294  return prodAmpCov(prodAmpIndices);
295 }
296 
297 
298 // returns normalization integral for pair of waves at index A and B
299 inline
300 std::complex<double>
301 TFitResult::normIntegral(const unsigned int waveIndexA,
302  const unsigned int waveIndexB) const
303 {
304  const TComplex norm = _normIntegral(waveIndexA, waveIndexB);
305  return std::complex<double>(norm.Re(), norm.Im());
306 }
307 
308 
309 // prints all production amplitude names
310 inline
311 std::ostream&
312 TFitResult::printProdAmpNames(std::ostream& out) const
313 {
314  out << " Production amplitude names:" << std::endl;
315  for (unsigned int i = 0; i < nmbProdAmps(); ++i)
316  out << " " << std::setw(3) << i << " " << _prodAmpNames[i] << std::endl;
317  return out;
318 }
319 
320 
321 // prints all wave names
322 inline
323 std::ostream&
324 TFitResult::printWaveNames(std::ostream& out) const
325 {
326  out << " Wave names:" << std::endl;
327  for (unsigned int i = 0; i < nmbWaves(); ++i)
328  out << " " << std::setw(3) << i << " " << _waveNames[i] << std::endl;
329  return out;
330 }
331 
332 
333 // prints all production amplitudes and their covariance matrix
334 inline
335 std::ostream&
336 TFitResult::printProdAmps(std::ostream& out) const
337 {
338  out << "Production amplitudes:" << std::endl;
339  for (unsigned int i = 0; i < nmbProdAmps(); ++i) {
340  out << " " << std::setw(3) << i << " " << _prodAmpNames[i] << " = " << prodAmp(i)
341  << ", cov = " << prodAmpCov(i) << std::endl;
342  }
343  return out;
344 }
345 
346 
347 // prints all wave intensities and their errors
348 inline
349 std::ostream&
350 TFitResult::printWaves(std::ostream& out) const
351 {
352  out << "Waves:" << std::endl;
353  for (unsigned int i = 0; i < nmbWaves(); ++i) {
354  out << " " << std::setw(3) << i << " " << _waveNames[i] << " = " << intensity(i)
355  << " +- " << intensityErr(i) << std::endl;
356  }
357  return out;
358 }
359 
360 
362 inline
363 std::ostream&
364 TFitResult::print(std::ostream& out) const
365 {
366  out << "TFitResult dump:" << std::endl
367  << " number of events .................... " << _nmbEvents << std::endl
368  << " number of events to normalize to .... " << _normNmbEvents << std::endl
369  << " center value of mass bin ............ " << _massBinCenter << std::endl
370  << " log(likelihood) at maximum .......... " << _logLikelihood << std::endl
371  << " rank of fit ......................... " << _rank << std::endl
372  << " bin has a valid covariance matrix ... " << _hasErrors << std::endl;
373  printProdAmps(out);
374  printWaveNames(out);
375  out << " covariance matrix:" << std::endl << _fitParCovMatrix << std::endl;
376  out << " covariance matrix indices:" << std::endl;
377  for (unsigned int i = 0; i < _fitParCovMatrixIndices.size(); ++i)
378  out << " index " << std::setw(3) << i << " = (" << std::setw(3) << _fitParCovMatrixIndices[i].first
379  << ", " << std::setw(3) << _fitParCovMatrixIndices[i].second << ")" << std::endl;
380  out << " normalization integral (w/o acceptance):" << std::endl << _normIntegral << std::endl;
381  out << " map of production amplitude indices to indices in normalization integral:" << std::endl;
382  for (std::map<Int_t, Int_t>::const_iterator i = _normIntIndexMap.begin(); i != _normIntIndexMap.end(); ++i)
383  out << " prod. amp [" << std::setw(3) << i->first << "] "
384  << "-> norm. int. [" << std::setw(3) << i->second << "]" << std::endl;
385  return out;
386 }
387 
388 
393 // !!! possible optimization: return TMatrixTLazy
394 inline
395 TMatrixT<double>
396 TFitResult::matrixRepr(const std::complex<double>& c)
397 {
398  TMatrixT<double> m(2, 2);
399  m[0][0] = m[1][1] = c.real();
400  m[0][1] = -c.imag();
401  m[1][0] = c.imag();
402  return m;
403 }
404 
405 
407 inline
408 int
409 TFitResult::rankOfProdAmp(const unsigned int prodAmpIndex) const
410 {
411  const TString ampName = _prodAmpNames[prodAmpIndex];
412  if (ampName.Contains("flat"))
413  return -1;
414  else
415  // the rank is encoded in second character of parameter name
416  return TString(ampName(1, 1)).Atoi();
417 }
418 
419 
421 inline
422 std::vector<unsigned int>
423 TFitResult::waveIndicesMatchingPattern(const std::string& waveNamePattern) const
424 {
425  std::vector<unsigned int> waveIndices;
426  for (unsigned int waveIndex = 0; waveIndex < nmbWaves(); ++waveIndex)
427  if (waveName(waveIndex).Contains(waveNamePattern))
428  waveIndices.push_back(waveIndex);
429  return waveIndices;
430 }
431 
432 
434 inline
435 std::vector<unsigned int>
436 TFitResult::prodAmpIndicesMatchingPattern(const std::string& ampNamePattern) const
437 {
438  std::vector<unsigned int> prodAmpIndices;
439  for (unsigned int prodAmpIndex = 0; prodAmpIndex < nmbProdAmps(); ++prodAmpIndex)
440  if (prodAmpName(prodAmpIndex).Contains(ampNamePattern))
441  prodAmpIndices.push_back(prodAmpIndex);
442  return prodAmpIndices;
443 }
444 
445 
447 inline
448 std::vector<std::pair<unsigned int, unsigned int> >
449 TFitResult::prodAmpIndexPairsForWaves(const unsigned int waveIndexA,
450  const unsigned int waveIndexB) const
451 {
452  std::vector<std::pair<unsigned int, unsigned int> > prodAmpIndexPairs;
453  const std::vector<unsigned int> prodAmpIndicesA = prodAmpIndicesForWave(waveIndexA);
454  const std::vector<unsigned int> prodAmpIndicesB = prodAmpIndicesForWave(waveIndexB);
455  for (unsigned int countAmpA = 0; countAmpA < prodAmpIndicesA.size(); ++countAmpA) {
456  const unsigned int ampIndexA = prodAmpIndicesA[countAmpA];
457  const int ampRankA = rankOfProdAmp(ampIndexA);
458  // find production amplitude of wave B with same rank
459  int ampIndexB = -1;
460  for (unsigned int countAmpB = 0; countAmpB < prodAmpIndicesB.size(); ++countAmpB)
461  if (rankOfProdAmp(prodAmpIndicesB[countAmpB]) == ampRankA) {
462  ampIndexB = prodAmpIndicesB[countAmpB];
463  break;
464  }
465  if (ampIndexB >=0)
466  prodAmpIndexPairs.push_back(std::make_pair(ampIndexA, (unsigned int)ampIndexB));
467  }
468  return prodAmpIndexPairs;
469 }
470 
471 
473 inline
474 double
475 TFitResult::realValVariance(const unsigned int waveIndexA,
476  const unsigned int waveIndexB,
477  const TMatrixT<double>& jacobian) const // Jacobian of real valued function (d f/ d Re[rho] d f / d Im[rho])
478 {
479  if (!_hasErrors)
480  return 0;
481  const TMatrixT<double> spinDensCov = spinDensityMatrixElemCov(waveIndexA, waveIndexB); // 2 x 2 matrix
482  const TMatrixT<double> jacobianT(TMatrixT<double>::kTransposed, jacobian); // 2 x 1 matrix
483  const TMatrixT<double> spinDensCovJT = spinDensCov * jacobianT; // 2 x 1 matrix
484  const TMatrixT<double> cov = jacobian * spinDensCovJT; // 1 x 1 matrix
485  return cov[0][0];
486 }
487 
488 
489 #endif // USE_TFITRESULT
490 
491 
492 #endif // TFITRESULT_HH
493 
494 
495 //--------------------------------------------------------------
496 // $Log$
497 //--------------------------------------------------------------
Generated by   doxygen 1.8.1.2