ROOTPWA
TFitResult.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 // 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 #include <algorithm>
44 
45 #include "integral.h" // PWA2000 classes
46 
47 #include "TFitResult.h"
48 
49 
50 #ifdef USE_TFITRESULT
51 
52 
53 using namespace std;
54 
55 
56 ClassImp(TFitResult);
57 
58 
59 TFitResult::TFitResult()
60  : _nmbEvents (0),
61  _normNmbEvents(0),
62  _massBinCenter(0),
63  _logLikelihood(0),
64  _rank (0),
65  _hasErrors (false)
66 { }
67 
68 
69 TFitResult::TFitResult(const TFitBin& fitBin)
70  : _nmbEvents (fitBin.nmbEvents()),
71  _normNmbEvents (fitBin.normNmbEvents()),
72  _massBinCenter (fitBin.massBinCenter()),
73  _logLikelihood (fitBin.logLikelihood()),
74  _rank (fitBin.rank()),
75  _hasErrors (fitBin.fitParCovMatrixValid()),
76  _fitParCovMatrix (fitBin.fitParCovMatrix()),
77  _fitParCovMatrixIndices(fitBin.fitParCovIndices()),
78  _normIntegral (fitBin.normIntegral()),
79  _normIntIndexMap (fitBin.prodAmpIndexMap())
80 {
81  _prodAmps = fitBin.prodAmps();
82  {
83  const unsigned int nmbAmpNames = fitBin.prodAmpNames().size();
84  _prodAmpNames.resize(nmbAmpNames, "");
85  for (unsigned int i = 0; i < nmbAmpNames; ++i)
86  _prodAmpNames[i] = fitBin.prodAmpNames()[i].Data();
87  }
88  {
89  const unsigned int nmbWaveNames = fitBin.waveNames().size();
90  _waveNames.resize(nmbWaveNames, "");
91  for (unsigned int i = 0; i < nmbWaveNames; ++i)
92  _waveNames[i] = fitBin.waveNames()[i].Data();
93  }
94 }
95 
96 
97 TFitResult::TFitResult(const TFitResult& result)
98  : _nmbEvents (result.nmbEvents()),
99  _normNmbEvents (result.normNmbEvents()),
100  _massBinCenter (result.massBinCenter()),
101  _logLikelihood (result.logLikelihood()),
102  _rank (result.rank()),
103  _prodAmps (result.prodAmps()),
104  _prodAmpNames (result.prodAmpNames()),
105  _waveNames (result.waveNames()),
106  _hasErrors (result.covMatrixValid()),
107  _fitParCovMatrix (result.fitParCovMatrix()),
108  _fitParCovMatrixIndices(result.fitParCovIndices()),
109  _normIntegral (result.normIntegralMatrix()),
110  _normIntIndexMap (result.normIntIndexMap())
111 { }
112 
113 
114 TFitResult::~TFitResult()
115 { }
116 
117 
123 double
124 TFitResult::evidence() const {
125  double l = -logLikelihood();
126  double det = _fitParCovMatrix.Determinant();
127  double d = (double)_fitParCovMatrix.GetNcols();
128  double sum = 0;
129  unsigned int ni = _normIntegral.ncols();
130  for(unsigned int i=0;i<ni;++i){
131  sum += 1. / _normIntegral(i, i).Re();
132  }
133  double occ = TMath::Power(TMath::Pi(), d * 0.5 - 1.) * TMath::Sqrt(2 * det) / sum;
134  double locc = TMath::Log(occ);
135  cout << "TFitResult::evidence" << endl;
136  cout << " LogLikeli: " << l;
137  cout << " Occamfactor: " << locc;
138  return l + locc;
139 }
140 
141 
145 complex<double>
146 TFitResult::spinDensityMatrixElem(const unsigned int waveIndexA,
147  const unsigned int waveIndexB) const
148 {
149  // get pairs of amplitude indices with the same rank for waves A and B
150  const vector<pair<unsigned int, unsigned int> > prodAmpIndexPairs = prodAmpIndexPairsForWaves(waveIndexA, waveIndexB);
151  if (prodAmpIndexPairs.size() == 0)
152  return 0;
153  // sum up amplitude products
154  complex<double> spinDens = 0;
155  for (unsigned int i = 0; i < prodAmpIndexPairs.size(); ++i)
156  spinDens += prodAmp(prodAmpIndexPairs[i].first) * conj(prodAmp(prodAmpIndexPairs[i].second));
157  return spinDens;
158 }
159 
160 
162 double
163 TFitResult::fitParameter(const string& parName) const
164 {
165  // check if parameter corresponds to real or imaginary part of production amplitude
166  TString name(parName);
167  const bool realPart = (name.Contains("RE") || name.Contains("flat"));
168  // find corresponding production amplitude
169  if (realPart)
170  name.ReplaceAll("_RE", "");
171  else
172  name.ReplaceAll("_IM", "");
173  for (unsigned int i = 0; i < nmbProdAmps(); ++i)
174  if (name == _prodAmpNames[i]) {
175  if (realPart)
176  return prodAmp(i).real();
177  else
178  return prodAmp(i).imag();
179  }
180  // not found
181  printWarn << "Could not find any parameter named '" << parName << "'." << endl;
182  return 0;
183 }
184 
185 
187 // where n is the number of amplitudes
188 // layout:
189 // cov(A0.re, A0.re) cov(A0.re, A0.im) ... cov(A0.re, A(n - 1).re) cov(A0.re, A(n - 1).im)
190 // cov(A0.im, A0.re) cov(A0.im, A0.im) ... cov(A0.im, A(n - 1).re) cov(A0.im, A(n - 1).im)
191 // . . . .
192 // . . . .
193 // . . . .
194 // cov(A(n - 1).re, A0.re) cov(A(n - 1).re, A0.im) ... cov(A(n - 1).re, A(n - 1).re) cov(A(n - 1).re, A(n - 1).im)
195 // cov(A(n - 1).im, A0.re) cov(A(n - 1).im, A0.im) ... cov(A(n - 1).im, A(n - 1).re) cov(A(n - 1).im, A(n - 1).im)
196 // !!! possible optimization: exploit symmetry of cov matrix
197 TMatrixT<double>
198 TFitResult::prodAmpCov(const vector<unsigned int>& prodAmpIndices) const
199 {
200  const unsigned int dim = 2 * prodAmpIndices.size();
201  TMatrixT<double> prodAmpCov(dim, dim);
202  if (!_hasErrors) {
203  printWarn << "TFitResult does not have a valid error matrix. Returning zero covariance matrix." << endl;
204  return prodAmpCov;
205  }
206  // get corresponding indices for parameter covariances
207  vector<int> parCovIndices;
208  for (unsigned int i = 0; i < prodAmpIndices.size(); ++i) {
209  parCovIndices.push_back(_fitParCovMatrixIndices[prodAmpIndices[i]].first); // real part
210  parCovIndices.push_back(_fitParCovMatrixIndices[prodAmpIndices[i]].second); // imaginary part
211  }
212  // build covariance matrix
213  for (unsigned int row = 0; row < dim; ++row)
214  for (unsigned int col = 0; col < dim; ++col) {
215  const int i = parCovIndices[row];
216  const int j = parCovIndices[col];
217  if ((i >= 0) && (j >= 0))
218  prodAmpCov[row][col] = fitParameterCov(i, j);
219  }
220  return prodAmpCov;
221 }
222 
223 
227 // !!! possible optimization: make special case for waveIndexA == waveIndexB
228 TMatrixT<double>
229 TFitResult::spinDensityMatrixElemCov(const unsigned int waveIndexA,
230  const unsigned int waveIndexB) const
231 {
232  // get pairs of amplitude indices with the same rank for waves A and B
233  const vector<pair<unsigned int, unsigned int> > prodAmpIndexPairs = prodAmpIndexPairsForWaves(waveIndexA, waveIndexB);
234  if (!_hasErrors || (prodAmpIndexPairs.size() == 0)) {
235  TMatrixT<double> spinDensCov(2, 2);
236  return spinDensCov;
237  }
238  // build covariance matrix for amplitudes
239  const TMatrixT<double> prodAmpCov = this->prodAmpCov(prodAmpIndexPairs);
240  // build Jacobian for rho_AB, which is a 2 x 4m matrix composed of 2m sub-Jacobians:
241  // J = (JA0, ..., JA(m - 1), JB0, ..., JB(m - 1))
242  // m is the number of production amplitudes for waves A and B that have the same rank
243  const unsigned int dim = prodAmpCov.GetNcols();
244  TMatrixT<double> jacobian(2, dim);
245  // build m sub-Jacobians for d rho_AB / d V_Ar = M(V_Br^*)
246  for (unsigned int i = 0; i < prodAmpIndexPairs.size(); ++i) {
247  const unsigned int ampIndexB = prodAmpIndexPairs[i].second;
248  const TMatrixT<double> subJacobianA = matrixRepr(conj(prodAmp(ampIndexB)));
249  jacobian.SetSub(0, 2 * i, subJacobianA);
250  }
251  // build m sub-Jacobian for d rho_AB / d V_Br = M(V_Ar) {{1, 0},
252  // {0, -1}}
253  TMatrixT<double> M(2, 2); // complex conjugation of V_Br is non-analytic operation
254  M[0][0] = 1;
255  M[1][1] = -1;
256  const unsigned int colOffset = 2 * prodAmpIndexPairs.size();
257  for (unsigned int i = 0; i < prodAmpIndexPairs.size(); ++i) {
258  const unsigned int ampIndexA = prodAmpIndexPairs[i].first;
259  TMatrixT<double> subJacobianB = matrixRepr(prodAmp(ampIndexA));
260  subJacobianB *= M;
261  jacobian.SetSub(0, colOffset + 2 * i, subJacobianB);
262  }
263  // calculate spin density covariance matrix cov(rho_AB) = J cov(V_A0, ..., V_A(m - 1), V_B0, ..., V_B(m - 1)) J^T
264  const TMatrixT<double> jacobianT(TMatrixT<double>::kTransposed, jacobian);
265  // binary operations are unavoidable, since matrices are not squared
266  // !!! possible optimaztion: use special TMatrixT constructors to perform the multiplication
267  const TMatrixT<double> prodAmpCovJT = prodAmpCov * jacobianT;
268  const TMatrixT<double> spinDensCov = jacobian * prodAmpCovJT;
269  return spinDensCov;
270 }
271 
272 
276 double
277 TFitResult::intensity(const char* waveNamePattern) const
278 {
279  vector<unsigned int> waveIndices = waveIndicesMatchingPattern(waveNamePattern);
280  double intensity = 0;
281  for (unsigned int i = 0; i < waveIndices.size(); ++i) {
282  intensity += this->intensity(waveIndices[i]);
283  for (unsigned int j = 0; j < i; ++j)
284  intensity += overlap(waveIndices[i], waveIndices[j]);
285  }
286  return intensity;
287 }
288 
289 
291 complex<double>
292 TFitResult::normIntegralForProdAmp(const unsigned int prodAmpIndexA,
293  const unsigned int prodAmpIndexB) const
294 {
295  // treat special case of flat wave which has no normalization integral
296  const bool flatWaveA = prodAmpName(prodAmpIndexA).Contains("flat");
297  const bool flatWaveB = prodAmpName(prodAmpIndexB).Contains("flat");
298  if (flatWaveA && flatWaveB)
299  return 1;
300  else if (flatWaveA || flatWaveB)
301  return 0;
302  else {
303  map<int, int>::const_iterator indexA = _normIntIndexMap.find(prodAmpIndexA);
304  map<int, int>::const_iterator indexB = _normIntIndexMap.find(prodAmpIndexB);
305  if ((indexA == _normIntIndexMap.end()) || (indexB == _normIntIndexMap.end())) {
306  printWarn << "Amplitude index " << prodAmpIndexA << " or " << prodAmpIndexB << " is out of bound." << endl;
307  return 0;
308  }
309  return normIntegral(indexA->second, indexB->second);
310  }
311 }
312 
313 
317 double
318 TFitResult::intensityErr(const char* waveNamePattern) const
319 {
320  // get amplitudes that correspond to wave name pattern
321  const vector<unsigned int> prodAmpIndices = prodAmpIndicesMatchingPattern(waveNamePattern);
322  const unsigned int nmbAmps = prodAmpIndices.size();
323  if (!_hasErrors || (nmbAmps == 0))
324  return 0;
325  // build Jacobian for intensity, which is a 1 x 2n matrix composed of n sub-Jacobians:
326  // J = (JA_0, ..., JA_{n - 1}), where n is the number of production amplitudes
327  TMatrixT<double> jacobian(1, 2 * nmbAmps);
328  for (unsigned int i = 0; i < nmbAmps; ++i) {
329  // build sub-Jacobian for each amplitude; intensity is real valued function, so J has only one row
330  // JA_ir = 2 * sum_j (A_jr Norm_ji)
331  complex<double> ampNorm = 0; // sum_j (A_jr Norm_ji)
332  const int currentRank = rankOfProdAmp(prodAmpIndices[i]);
333  for (unsigned int j = 0; j < nmbAmps; ++j) {
334  if (rankOfProdAmp(prodAmpIndices[j]) != currentRank)
335  continue;
336  ampNorm += prodAmp(prodAmpIndices[j]) * normIntegralForProdAmp(j, i); // order of indices is essential
337  }
338  jacobian[0][2 * i ] = ampNorm.real();
339  jacobian[0][2 * i + 1] = ampNorm.imag();
340  }
341  jacobian *= 2;
342  const TMatrixT<double> prodAmpCov = this->prodAmpCov(prodAmpIndices); // 2n x 2n matrix
343  const TMatrixT<double> jacobianT(TMatrixT<double>::kTransposed, jacobian); // 2n x 1 matrix
344  const TMatrixT<double> prodAmpCovJT = prodAmpCov * jacobianT; // 2n x 1 matrix
345  const TMatrixT<double> intensityCov = jacobian * prodAmpCovJT; // 1 x 1 matrix
346  return sqrt(intensityCov[0][0]);
347 }
348 
349 
351 double
352 TFitResult::phase(const unsigned int waveIndexA,
353  const unsigned int waveIndexB) const
354 {
355  if (waveIndexA == waveIndexB)
356  return 0;
357  return arg(spinDensityMatrixElem(waveIndexA, waveIndexB)) * TMath::RadToDeg();
358 }
359 
360 
362 double
363 TFitResult::phaseErr(const unsigned int waveIndexA,
364  const unsigned int waveIndexB) const
365 {
366  if (!_hasErrors || (waveIndexA == waveIndexB))
367  return 0;
368  // construct Jacobian for phi_AB = +- arctan(Im[rho_AB] / Re[rho_AB])
369  const complex<double> spinDens = spinDensityMatrixElem(waveIndexA, waveIndexB);
370  TMatrixT<double> jacobian(1, 2); // phase is real valued function, so J has only one row
371  {
372  const double x = spinDens.real();
373  const double y = spinDens.imag();
374  if ((x != 0) || (y != 0)) {
375  jacobian[0][0] = 1 / (x + y * y / x);
376  jacobian[0][1] = -y / (x * x + y * y);
377  }
378  }
379  // calculate variance
380  const double phaseVariance = realValVariance(waveIndexA, waveIndexB, jacobian);
381  return sqrt(phaseVariance) * TMath::RadToDeg();
382 }
383 
384 
386 double
387 TFitResult::coherence(const unsigned int waveIndexA,
388  const unsigned int waveIndexB) const
389 {
390  const double rhoAA = spinDensityMatrixElem(waveIndexA, waveIndexA).real(); // rho_AA is real by definition
391  const double rhoBB = spinDensityMatrixElem(waveIndexB, waveIndexB).real(); // rho_BB is real by definition
392  const complex<double> rhoAB = spinDensityMatrixElem(waveIndexA, waveIndexB);
393 
394  cout << "coherence("<<waveIndexA<<","<<waveIndexB<<"): "
395  <<rhoAA<<" "
396  <<rhoBB<<" "
397  <<rhoAB<<" "
398  <<sqrt(std::norm(rhoAB))<<" "
399  <<sqrt((rhoAA * rhoBB))<<endl;
400 
401 
402 
403  return sqrt(std::norm(rhoAB) / (rhoAA * rhoBB));
404 }
405 
406 
408 double
409 TFitResult::coherenceErr(const unsigned int waveIndexA,
410  const unsigned int waveIndexB) const
411 {
412  // get amplitude indices for waves A and B
413  const vector<unsigned int> prodAmpIndices[2] = {prodAmpIndicesForWave(waveIndexA),
414  prodAmpIndicesForWave(waveIndexB)};
415  if (!_hasErrors || (prodAmpIndices[0].size() == 0) || (prodAmpIndices[1].size() == 0))
416  return 0;
417  // build Jacobian for coherence, which is a 1 x 2(n + m) matrix composed of (n + m) sub-Jacobians:
418  // J = (JA_0, ..., JA_{n - 1}, JB_0, ..., JB_{m - 1})
419  const unsigned int nmbAmps = prodAmpIndices[0].size() + prodAmpIndices[1].size();
420  TMatrixT<double> jacobian(1, 2 * nmbAmps);
421  // precalculate some variables
422  const double rhoAA = spinDensityMatrixElem(waveIndexA, waveIndexA).real(); // rho_AA is real by definition
423  const double rhoBB = spinDensityMatrixElem(waveIndexB, waveIndexB).real(); // rho_BB is real by definition
424  const complex<double> rhoAB = spinDensityMatrixElem(waveIndexA, waveIndexB);
425  const double rhoABRe = rhoAB.real();
426  const double rhoABIm = rhoAB.imag();
427  const double rhoABNorm = std::norm(rhoAB);
428  const double coh = sqrt(rhoABNorm / (rhoAA * rhoBB));
429  if (coh == 0)
430  return 0;
431  // build m sub-Jacobians for JA_r = coh_AB / d V_Ar
432  for (unsigned int i = 0; i < prodAmpIndices[0].size(); ++i) {
433  const unsigned int prodAmpIndexA = prodAmpIndices[0][i];
434  const complex<double> prodAmpA = prodAmp(prodAmpIndexA);
435  const int prodAmpRankA = rankOfProdAmp(prodAmpIndexA);
436  // find production amplitude of wave B with same rank
437  complex<double> prodAmpB = 0;
438  for (unsigned int j = 0; j < prodAmpIndices[1].size(); ++j) {
439  const unsigned int prodAmpIndexB = prodAmpIndices[1][j];
440  if (rankOfProdAmp(prodAmpIndexB) == prodAmpRankA) {
441  prodAmpB = prodAmp(prodAmpIndexB);
442  break;
443  }
444  }
445  jacobian[0][2 * i ] = rhoABRe * prodAmpB.real() - rhoABIm * prodAmpB.imag() - (rhoABNorm / rhoAA) * prodAmpA.real();
446  jacobian[0][2 * i + 1] = rhoABRe * prodAmpB.imag() + rhoABIm * prodAmpB.real() - (rhoABNorm / rhoAA) * prodAmpA.imag();
447  }
448  // !!! possible optimization: join the loops for JA_r and JB_r
449  // build m sub-Jacobian for JB_r = d coh_AB / d V_Br
450  const unsigned int colOffset = 2 * prodAmpIndices[0].size();
451  for (unsigned int i = 0; i < prodAmpIndices[1].size(); ++i) {
452  const unsigned int prodAmpIndexB = prodAmpIndices[1][i];
453  const complex<double> prodAmpB = prodAmp(prodAmpIndexB);
454  const int prodAmpRankB = rankOfProdAmp(prodAmpIndexB);
455  // find production amplitude of wave A with same rank
456  complex<double> prodAmpA = 0;
457  for (unsigned int j = 0; j < prodAmpIndices[0].size(); ++j) {
458  const unsigned int prodAmpIndexA = prodAmpIndices[0][j];
459  if (rankOfProdAmp(prodAmpIndexA) == prodAmpRankB) {
460  prodAmpA = prodAmp(prodAmpIndexA);
461  break;
462  }
463  }
464  jacobian[0][colOffset + 2 * i ] = rhoABRe * prodAmpA.real() + rhoABIm * prodAmpA.imag() - (rhoABNorm / rhoBB) * prodAmpB.real();
465  jacobian[0][colOffset + 2 * i + 1] = rhoABRe * prodAmpA.imag() - rhoABIm * prodAmpA.real() - (rhoABNorm / rhoBB) * prodAmpB.imag();
466  }
467  jacobian *= 1 / (coh * rhoAA * rhoBB);
468  // build covariance matrix for amplitudes and calculate coherence covariance matrix
469  const TMatrixT<double> prodAmpCov = this->prodAmpCov(prodAmpIndices[0], prodAmpIndices[1]); // 2(n + m) x 2(n + m) matrix
470  const TMatrixT<double> jacobianT(TMatrixT<double>::kTransposed, jacobian); // 2(n + m) x 1 matrix
471  const TMatrixT<double> prodAmpCovJT = prodAmpCov * jacobianT; // 2(n + m) x 1 matrix
472  const TMatrixT<double> cohCov = jacobian * prodAmpCovJT; // 1 x 1 matrix
473  return sqrt(cohCov[0][0]);
474 }
475 
476 
478 double
479 TFitResult::overlap(const unsigned int waveIndexA,
480  const unsigned int waveIndexB) const
481 {
482  const complex<double> spinDens = spinDensityMatrixElem(waveIndexA, waveIndexB);
483  const complex<double> normInt = normIntegral (waveIndexA, waveIndexB);
484  return 2 * (spinDens * normInt).real();
485 }
486 
487 
489 double
490 TFitResult::overlapErr(const unsigned int waveIndexA,
491  const unsigned int waveIndexB) const
492 {
493  if (!_hasErrors)
494  return 0;
495  const complex<double> normInt = normIntegral(waveIndexA, waveIndexB);
496  TMatrixT<double> jacobian(1, 2); // overlap is real valued function, so J has only one row
497  jacobian[0][0] = 2 * normInt.real();
498  jacobian[0][1] = -2 * normInt.imag();
499  const double overlapVariance = realValVariance(waveIndexA, waveIndexB, jacobian);
500  return sqrt(overlapVariance);
501 }
502 
503 
504 void
505 TFitResult::reset()
506 {
507  _nmbEvents = 0;
508  _normNmbEvents = 0;
509  _massBinCenter = 0;
510  _logLikelihood = 0;
511  _rank = 0;
512  _prodAmps.clear();
513  _prodAmpNames.clear();
514  _waveNames.clear();
515  _hasErrors = false;
516  _fitParCovMatrix.ResizeTo(0, 0);
517  _fitParCovMatrixIndices.clear();
518  _normIntegral.ResizeTo(0, 0);
519  _normIntIndexMap.clear();
520 }
521 
522 
523 void
524 TFitResult::fill(const unsigned int nmbEvents, // number of events in bin
525  const unsigned int normNmbEvents, // number of events to normalize to
526  const double massBinCenter, // center value of mass bin
527  const double logLikelihood, // log(likelihood) at maximum
528  const int rank, // rank of fit
529  const vector<complex<double> >& prodAmps, // production amplitudes
530  const vector<string>& prodAmpNames, // names of production amplitudes used in fit
531  const TMatrixT<double>& fitParCovMatrix, // covariance matrix of fit parameters
532  const vector<pair<int, int> >& fitParCovMatrixIndices, // indices of fit parameters for real and imaginary part in covariance matrix matrix
533  const TCMatrix& normIntegral) // normalization integral over full phase space without acceptance
534 {
535  _nmbEvents = nmbEvents;
536  _normNmbEvents = normNmbEvents;
537  _massBinCenter = massBinCenter;
538  _logLikelihood = logLikelihood;
539  _rank = rank;
540  _prodAmps.resize(prodAmps.size());
541  for (unsigned int i = 0; i < prodAmps.size(); ++i)
542  _prodAmps[i] = TComplex(prodAmps[i].real(), prodAmps[i].imag());
543  _prodAmpNames = prodAmpNames;
544  _fitParCovMatrix.ResizeTo(fitParCovMatrix.GetNrows(), fitParCovMatrix.GetNcols());
545  _fitParCovMatrix = fitParCovMatrix;
546  _fitParCovMatrixIndices = fitParCovMatrixIndices;
547  if (!(fitParCovMatrix.GetNrows() == 0) && !(fitParCovMatrix.GetNcols() == 0)) // check whether there really is an error matrix
548  _hasErrors = true;
549  else
550  _hasErrors = false;
551  _normIntegral.ResizeTo(normIntegral.nrows(), normIntegral.ncols());
552  _normIntegral = normIntegral;
553 
554  buildWaveMap();
555 
556  // check consistency
557  if (_prodAmps.size() != _prodAmpNames.size())
558  cout << "TFitResult::fill(): warning: number of production amplitudes (" << _prodAmps.size()
559  << ") does not match number of production amplitude names (" << _prodAmpNames.size() << ")." << endl;
560  if (_prodAmps.size() != _fitParCovMatrixIndices.size())
561  cout << "TFitResult::fill(): warning: number of production amplitudes (" << _prodAmps.size()
562  << ") does not match number of covariance matrix indices (" << _fitParCovMatrixIndices.size() << ")." << endl;
563  if (((int)_waveNames.size() != _normIntegral.nrows()) || ((int)_waveNames.size() != _normIntegral.ncols()))
564  cout << "TFitResult::fill(): warning: number of waves (" << _waveNames.size()
565  << ") does not match size of normalization integral (" << _normIntegral.nrows() << ", " << _normIntegral.ncols() << ")." << endl;
566 }
567 
568 
569 void
570 TFitResult::buildWaveMap() {
571  int n=_prodAmpNames.size();
572  for(int i=0;i<n;++i){
573  // strip rank
574  TString title=wavetitle(i);
575  if(find(_waveNames.begin(),_waveNames.end(),title.Data())==_waveNames.end())
576  _waveNames.push_back(title.Data());
577 
578  // look for index of first occurence
579  int j;
580  for(j=0;j<n;++j)
581  if(prodAmpName(j).Contains(title))
582  break;
583  _normIntIndexMap[i]=j;
584  }
585 }
586 
587 
588 #endif // USE_TFITRESULT
Generated by   doxygen 1.8.1.2