ROOTPWA
fitResult.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 #include <algorithm>
39 
40 #include "Math/SpecFuncMathCore.h"
41 #include "Math/ProbFuncMathCore.h"
42 #include "TDecompChol.h"
43 #include "TMath.h"
44 #include "TMatrixDSym.h"
45 #include "TRandom3.h"
46 
47 // PWA2000 classes
48 #include "integral.h"
49 
50 #include "fitResult.h"
51 
52 
53 using namespace std;
54 using namespace rpwa;
55 
56 
57 ClassImp(fitResult);
58 
59 
60 fitResult::fitResult()
61  : _nmbEvents (0),
62  _normNmbEvents (0),
63  _massBinCenter (0),
64  _logLikelihood (0),
65  _rank (0),
66  _covMatrixValid(false),
67  _converged (false),
68  _hasHessian (false)
69 { }
70 
71 
72 fitResult::fitResult(const TFitBin& fitBin)
73  : _nmbEvents (fitBin.nmbEvents()),
74  _normNmbEvents (fitBin.normNmbEvents()),
75  _massBinCenter (fitBin.massBinCenter()),
76  _logLikelihood (fitBin.logLikelihood()),
77  _rank (fitBin.rank()),
78  _covMatrixValid (fitBin.fitParCovMatrixValid()),
79  _fitParCovMatrix (fitBin.fitParCovMatrix()),
80  _fitParCovMatrixIndices(fitBin.fitParCovIndices()),
81  _normIntegral (fitBin.normIntegral()),
82  _normIntIndexMap (fitBin.prodAmpIndexMap()),
83  _converged (true),
84  _hasHessian (false)
85 {
86  _prodAmps = fitBin.prodAmps();
87  {
88  const unsigned int nmbAmpNames = fitBin.prodAmpNames().size();
89  _prodAmpNames.resize(nmbAmpNames, "");
90  for (unsigned int i = 0; i < nmbAmpNames; ++i)
91  _prodAmpNames[i] = fitBin.prodAmpNames()[i].Data();
92  }
93  {
94  const unsigned int nmbWaveNames = fitBin.waveNames().size();
95  _waveNames.resize(nmbWaveNames, "");
96  for (unsigned int i = 0; i < nmbWaveNames; ++i)
97  _waveNames[i] = fitBin.waveNames()[i].Data();
98  }
99 }
100 
101 
102 fitResult::fitResult(const fitResult& result)
103  : _nmbEvents (result.nmbEvents()),
104  _normNmbEvents (result.normNmbEvents()),
105  _massBinCenter (result.massBinCenter()),
106  _logLikelihood (result.logLikelihood()),
107  _rank (result.rank()),
108  _prodAmps (result.prodAmps()),
109  _prodAmpNames (result.prodAmpNames()),
110  _waveNames (result.waveNames()),
111  _covMatrixValid (result.covMatrixValid()),
112  _fitParCovMatrix (result.fitParCovMatrix()),
113  _fitParCovMatrixIndices(result.fitParCovIndices()),
114  _normIntegral (result.normIntegralMatrix()),
115  _normIntIndexMap (result.normIntIndexMap()),
116  _phaseSpaceIntegral (result._phaseSpaceIntegral),
117  _converged (result._converged),
118  _hasHessian (result._hasHessian)
119 { }
120 
121 
122 // enable copying from TFitResult for older ROOT versions
123 #ifdef USE_TFITRESULT
124 fitResult::fitResult(const TFitResult& result)
125  : _nmbEvents (result.nmbEvents()),
126  _normNmbEvents (result.normNmbEvents()),
127  _massBinCenter (result.massBinCenter()),
128  _logLikelihood (result.logLikelihood()),
129  _rank (result.rank()),
130  _prodAmps (result.prodAmps()),
131  _prodAmpNames (result.prodAmpNames()),
132  _waveNames (result.waveNames()),
133  _covMatrixValid (result.covMatrixValid()),
134  _fitParCovMatrix (result.fitParCovMatrix()),
135  _fitParCovMatrixIndices(result.fitParCovIndices()),
136  _normIntegral (result.normIntegralMatrix()),
137  _normIntIndexMap (result.normIntIndexMap())
138 { }
139 #endif
140 
141 
142 fitResult::~fitResult()
143 { }
144 
145 
146 fitResult*
147 fitResult::cloneVar(){
148  // copy complete info
149  fitResult* result = new fitResult(*this);
150 
151  // Cholesky decomposition (hoepfully this will not slow us down too much
152  // if it does we have to cache
153  cerr << "Starting Cholesky Decomposition" << endl;
154  TDecompChol decomp(_fitParCovMatrix);
155  decomp.Decompose();
156  TMatrixT<Double_t> C(decomp.GetU());
157  cerr << "...Done" << endl;
158 
159  unsigned int npar=C.GetNrows();
160  TMatrixD x(npar,1);
161  // generate npar independent random numbers
162  for(unsigned int ipar=0;ipar<npar;++ipar){
163  x(ipar,0)=gRandom->Gaus();
164  }
165 
166  // this tells us how to permute the parameters taking into account all
167  // correlations in the covariance matrix
168  TMatrixD y(C,TMatrixD::kMult,x);
169 
170  // now we need mapping from parameters to production amp re and im
171  // loop over production amps
172  unsigned int nt=_prodAmps.size();
173  for(unsigned int it=0;it<nt;++it){
174  Int_t jre=_fitParCovMatrixIndices[it].first;
175  Int_t jim=_fitParCovMatrixIndices[it].second;
176  double re=result->_prodAmps[it].Re(); double im=result->_prodAmps[it].Im();
177  if(jre>-1)re+=y(jre,0);
178  if(jim>-1)im+=y(jim,0);
179  result->_prodAmps[it]=TComplex(re,im);
180  cerr << "Modifying amp " << it << " by (" << (jre>-1 ? y(jre,0) : 0)
181  << "," << (jim>-1 ? y(jim,0) : 0) <<")" << endl;
182  }
183  return result;
184 }
185 
186 
192 double
193 fitResult::evidence() const
194 {
195 
196  // REMOVE CONSTRAINT TO NUMBER OF EVENTS!
197  double l = -logLikelihood();// - intensity(".*");
198 
199  double det = _fitParCovMatrix.Determinant();
200  // simple determinant neglecting all off-diagonal entries
201  // unsigned int n= _fitParCovMatrix.GetNcols();
202  // double det2=1;
203  // for(unsigned int i=0;i<n;++i){
204  // det2*=_fitParCovMatrix[i][i];
205  // }
206 
207  double d = (double)_fitParCovMatrix.GetNcols();
208  double sum = 0;
209  unsigned int ni = _normIntegral.ncols();
210  for (unsigned int i = 0; i < ni ; ++i)
211  sum += 1. / _normIntegral(i, i).Re();
212  // parameter-volume after observing data
213  //double vad = TMath::Power(2*TMath::Pi(), d * 0.5) * TMath::Sqrt(det);
214  //double lvad = TMath::Log(vad);
215 
216  double lvad=0.5*(d*1.837877066+TMath::Log(det));
217 
218  // parameter volume prior to observing the data
219  // n-Sphere:
220  double lva= TMath::Log(d) + 0.5*(d*1.144729886+(d-1)*TMath::Log(_nmbEvents))-ROOT::Math::lgamma(0.5*d+1);
221 
222  // n-ball:
223  // double lva = 0.5*(d*1.144729886+d*TMath::Log(_nmbEvents))-ROOT::Math::lgamma(0.5*d+1);
224 
225  // finally we calculate the probability of single waves being negligible and
226  // take these reults into account
227 
228  unsigned int nwaves=nmbWaves();
229  double logprob=0;
230  for(unsigned int iwaves=0;iwaves<nwaves;++iwaves){
231  double val=intensity(iwaves);
232  double err=intensityErr(iwaves);
233  // P(val>0); (assuming gaussian...) dirty!
234  // require 3simga significance!
235  double prob=ROOT::Math::normal_cdf_c(5.*err,err,val);
236  //cerr << "val="<<val<<" err="<<err<<" prob(val>0)="<<prob<<endl;
237  logprob+=TMath::Log(prob);
238  }
239 
240 
241 
242  // cerr << "fitResult::evidence()" << endl
243  // << " det : " << det << endl
244  // //<< " detsimple : " << det2 << endl
245  // << " LogLikeli : " << l << endl
246  // << " logVA : " << lva << endl
247  // // << " logVASphere : " << lvaS << endl
248  // << " logVA|D : " << lvad << endl
249  // // << " logVA|D2 : " << lvad2 << endl
250  // << " Occamfactor : " << -lva+lvad << endl
251  // << " LogProb : " << logprob << endl
252  // << " evidence : " << l + lvad - lva + logprob<< endl;
253 
254  return l + lvad - lva + logprob;
255 }
256 
257 
261 complex<double>
262 fitResult::spinDensityMatrixElem(const unsigned int waveIndexA,
263  const unsigned int waveIndexB) const
264 {
265  // get pairs of amplitude indices with the same rank for waves A and B
266  const vector<pair<unsigned int, unsigned int> > prodAmpIndexPairs
267  = prodAmpIndexPairsForWaves(waveIndexA, waveIndexB);
268  if (prodAmpIndexPairs.size() == 0)
269  return 0;
270  // sum up amplitude products
271  complex<double> spinDens = 0;
272  for (unsigned int i = 0; i < prodAmpIndexPairs.size(); ++i)
273  spinDens += prodAmp(prodAmpIndexPairs[i].first) * conj(prodAmp(prodAmpIndexPairs[i].second));
274  return spinDens;
275 }
276 
277 
279 double
280 fitResult::fitParameter(const string& parName) const
281 {
282  // check if parameter corresponds to real or imaginary part of production amplitude
283  TString name(parName);
284  const bool realPart = (name.Contains("RE") || name.Contains("flat"));
285  // find corresponding production amplitude
286  if (realPart)
287  name.ReplaceAll("_RE", "");
288  else
289  name.ReplaceAll("_IM", "");
290  const int index = prodAmpIndex(name.Data());
291  if (index >= 0) {
292  if (realPart)
293  return prodAmp(index).real();
294  else
295  return prodAmp(index).imag();
296  }
297  return 0; // not found
298 }
299 
300 
302 // where n is the number of amplitudes
303 // layout:
304 // cov(A0.re, A0.re) cov(A0.re, A0.im) ... cov(A0.re, A(n - 1).re) cov(A0.re, A(n - 1).im)
305 // cov(A0.im, A0.re) cov(A0.im, A0.im) ... cov(A0.im, A(n - 1).re) cov(A0.im, A(n - 1).im)
306 // . . . .
307 // . . . .
308 // . . . .
309 // 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)
310 // 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)
311 // !!! possible optimization: exploit symmetry of cov matrix
312 TMatrixT<double>
313 fitResult::prodAmpCov(const vector<unsigned int>& prodAmpIndices) const
314 {
315  const unsigned int dim = 2 * prodAmpIndices.size();
316  TMatrixT<double> prodAmpCov(dim, dim);
317  if (!_covMatrixValid) {
318  printWarn << "fitResult does not have a valid error matrix. Returning zero covariance matrix." << endl;
319  return prodAmpCov;
320  }
321  // get corresponding indices for parameter covariances
322  vector<int> parCovIndices;
323  for (unsigned int i = 0; i < prodAmpIndices.size(); ++i) {
324  parCovIndices.push_back(_fitParCovMatrixIndices[prodAmpIndices[i]].first); // real part
325  parCovIndices.push_back(_fitParCovMatrixIndices[prodAmpIndices[i]].second); // imaginary part
326  }
327  // build covariance matrix
328  for (unsigned int row = 0; row < dim; ++row)
329  for (unsigned int col = 0; col < dim; ++col) {
330  const int i = parCovIndices[row];
331  const int j = parCovIndices[col];
332  if ((i >= 0) && (j >= 0))
333  prodAmpCov[row][col] = fitParameterCov(i, j);
334  }
335  return prodAmpCov;
336 }
337 
338 
342 // !!! possible optimization: make special case for waveIndexA == waveIndexB
343 TMatrixT<double>
344 fitResult::spinDensityMatrixElemCov(const unsigned int waveIndexA,
345  const unsigned int waveIndexB) const
346 {
347  // get pairs of amplitude indices with the same rank for waves A and B
348  const vector<pair<unsigned int, unsigned int> > prodAmpIndexPairs
349  = prodAmpIndexPairsForWaves(waveIndexA, waveIndexB);
350  if (!_covMatrixValid || (prodAmpIndexPairs.size() == 0)) {
351  TMatrixT<double> spinDensCov(2, 2);
352  return spinDensCov;
353  }
354  // build covariance matrix for amplitudes
355  const TMatrixT<double> prodAmpCov = this->prodAmpCov(prodAmpIndexPairs);
356  // build Jacobian for rho_AB, which is a 2 x 4m matrix composed of 2m sub-Jacobians:
357  // J = (JA0, ..., JA(m - 1), JB0, ..., JB(m - 1))
358  // m is the number of production amplitudes for waves A and B that have the same rank
359  const unsigned int dim = prodAmpCov.GetNcols();
360  TMatrixT<double> jacobian(2, dim);
361  // build m sub-Jacobians for d rho_AB / d V_Ar = M(V_Br^*)
362  for (unsigned int i = 0; i < prodAmpIndexPairs.size(); ++i) {
363  const unsigned int ampIndexB = prodAmpIndexPairs[i].second;
364  const TMatrixT<double> subJacobianA = matrixRepr(conj(prodAmp(ampIndexB)));
365  jacobian.SetSub(0, 2 * i, subJacobianA);
366  }
367  // build m sub-Jacobian for d rho_AB / d V_Br = M(V_Ar) {{1, 0},
368  // {0, -1}}
369  TMatrixT<double> M(2, 2); // complex conjugation of V_Br is non-analytic operation
370  M[0][0] = 1;
371  M[1][1] = -1;
372  const unsigned int colOffset = 2 * prodAmpIndexPairs.size();
373  for (unsigned int i = 0; i < prodAmpIndexPairs.size(); ++i) {
374  const unsigned int ampIndexA = prodAmpIndexPairs[i].first;
375  TMatrixT<double> subJacobianB = matrixRepr(prodAmp(ampIndexA));
376  subJacobianB *= M;
377  jacobian.SetSub(0, colOffset + 2 * i, subJacobianB);
378  }
379  // calculate spin density covariance matrix
380  // cov(rho_AB) = J cov(V_A0, ..., V_A(m - 1), V_B0, ..., V_B(m - 1)) J^T
381  const TMatrixT<double> jacobianT(TMatrixT<double>::kTransposed, jacobian);
382  // binary operations are unavoidable, since matrices are not squared
383  // !!! possible optimaztion: use special TMatrixT constructors to perform the multiplication
384  const TMatrixT<double> prodAmpCovJT = prodAmpCov * jacobianT;
385  const TMatrixT<double> spinDensCov = jacobian * prodAmpCovJT;
386  return spinDensCov;
387 }
388 
389 
393 double
394 fitResult::intensity(const char* waveNamePattern) const
395 {
396  vector<unsigned int> waveIndices = waveIndicesMatchingPattern(waveNamePattern);
397  double intensity = 0;
398  for (unsigned int i = 0; i < waveIndices.size(); ++i) {
399  intensity += this->intensity(waveIndices[i]);
400  // cout << " contribution from " << _waveNames[waveIndices[i]]
401  // << " = " << this->intensity(waveIndices[i]) << endl;
402  for (unsigned int j = 0; j < i; ++j) {
403  intensity += overlap(waveIndices[i], waveIndices[j]);
404  // cout << " overlap with " << _waveNames[waveIndices[j]]
405  // << " = " << overlap(waveIndices[i], waveIndices[j]) << endl;
406  }
407  }
408  return intensity;
409 }
410 
411 
413 complex<double>
414 fitResult::normIntegralForProdAmp(const unsigned int prodAmpIndexA,
415  const unsigned int prodAmpIndexB) const
416 {
417  // treat special case of flat wave which has no normalization integral
418  const bool flatWaveA = prodAmpName(prodAmpIndexA).Contains("flat");
419  const bool flatWaveB = prodAmpName(prodAmpIndexB).Contains("flat");
420  if (flatWaveA && flatWaveB)
421  return 1;
422  else if (flatWaveA || flatWaveB)
423  return 0;
424  else {
425  map<int, int>::const_iterator indexA = _normIntIndexMap.find(prodAmpIndexA);
426  map<int, int>::const_iterator indexB = _normIntIndexMap.find(prodAmpIndexB);
427  if ((indexA == _normIntIndexMap.end()) || (indexB == _normIntIndexMap.end())) {
428  printWarn << "Amplitude index " << prodAmpIndexA << " or " << prodAmpIndexB
429  << " is out of bound." << endl;
430  return 0;
431  }
432  return normIntegral(indexA->second, indexB->second);
433  }
434 }
435 
436 
440 double
441 fitResult::intensityErr(const char* waveNamePattern) const
442 {
443  // get amplitudes that correspond to wave name pattern
444  const vector<unsigned int> prodAmpIndices = prodAmpIndicesMatchingPattern(waveNamePattern);
445  const unsigned int nmbAmps = prodAmpIndices.size();
446  if (!_covMatrixValid || (nmbAmps == 0))
447  return 0;
448  // build Jacobian for intensity, which is a 1 x 2n matrix composed of n sub-Jacobians:
449  // J = (JA_0, ..., JA_{n - 1}), where n is the number of production amplitudes
450  TMatrixT<double> jacobian(1, 2 * nmbAmps);
451  for (unsigned int i = 0; i < nmbAmps; ++i) {
452  // build sub-Jacobian for each amplitude; intensity is real valued function, so J has only one row
453  // JA_ir = 2 * sum_j (A_jr Norm_ji)
454  complex<double> ampNorm = 0; // sum_j (A_jr Norm_ji)
455  const int currentRank = rankOfProdAmp(prodAmpIndices[i]);
456  for (unsigned int j = 0; j < nmbAmps; ++j) {
457  if (rankOfProdAmp(prodAmpIndices[j]) != currentRank)
458  continue;
459  ampNorm += prodAmp(prodAmpIndices[j]) * normIntegralForProdAmp(j, i); // order of indices is essential
460  }
461  jacobian[0][2 * i ] = ampNorm.real();
462  jacobian[0][2 * i + 1] = ampNorm.imag();
463  }
464  jacobian *= 2;
465  const TMatrixT<double> prodAmpCov = this->prodAmpCov(prodAmpIndices); // 2n x 2n matrix
466  const TMatrixT<double> jacobianT(TMatrixT<double>::kTransposed, jacobian); // 2n x 1 matrix
467  const TMatrixT<double> prodAmpCovJT = prodAmpCov * jacobianT; // 2n x 1 matrix
468  const TMatrixT<double> intensityCov = jacobian * prodAmpCovJT; // 1 x 1 matrix
469  return sqrt(intensityCov[0][0]);
470 }
471 
472 
474 double
475 fitResult::phase(const unsigned int waveIndexA,
476  const unsigned int waveIndexB) const
477 {
478  if (waveIndexA == waveIndexB)
479  return 0;
480  return arg(spinDensityMatrixElem(waveIndexA, waveIndexB)) * TMath::RadToDeg();
481 }
482 
483 
485 double
486 fitResult::phaseErr(const unsigned int waveIndexA,
487  const unsigned int waveIndexB) const
488 {
489  if (!_covMatrixValid || (waveIndexA == waveIndexB))
490  return 0;
491  // construct Jacobian for phi_AB = +- arctan(Im[rho_AB] / Re[rho_AB])
492  const complex<double> spinDens = spinDensityMatrixElem(waveIndexA, waveIndexB);
493  TMatrixT<double> jacobian(1, 2); // phase is real valued function, so J has only one row
494  {
495  const double x = spinDens.real();
496  const double y = spinDens.imag();
497  if ((x != 0) || (y != 0)) {
498  jacobian[0][0] = 1 / (x + y * y / x);
499  jacobian[0][1] = -y / (x * x + y * y);
500  }
501  }
502  // calculate variance
503  const double phaseVariance = realValVariance(waveIndexA, waveIndexB, jacobian);
504  return sqrt(phaseVariance) * TMath::RadToDeg();
505 }
506 
507 
509 double
510 fitResult::coherence(const unsigned int waveIndexA,
511  const unsigned int waveIndexB) const
512 {
513  const double rhoAA = spinDensityMatrixElem(waveIndexA, waveIndexA).real(); // rho_AA is real by definition
514  const double rhoBB = spinDensityMatrixElem(waveIndexB, waveIndexB).real(); // rho_BB is real by definition
515  const complex<double> rhoAB = spinDensityMatrixElem(waveIndexA, waveIndexB);
516  return sqrt(std::norm(rhoAB) / (rhoAA * rhoBB));
517 }
518 
519 
521 double
522 fitResult::coherenceErr(const unsigned int waveIndexA,
523  const unsigned int waveIndexB) const
524 {
525  // get amplitude indices for waves A and B
526  const vector<unsigned int> prodAmpIndices[2] = {prodAmpIndicesForWave(waveIndexA),
527  prodAmpIndicesForWave(waveIndexB)};
528  if (!_covMatrixValid || (prodAmpIndices[0].size() == 0) || (prodAmpIndices[1].size() == 0))
529  return 0;
530  // build Jacobian for coherence, which is a 1 x 2(n + m) matrix composed of (n + m) sub-Jacobians:
531  // J = (JA_0, ..., JA_{n - 1}, JB_0, ..., JB_{m - 1})
532  const unsigned int nmbAmps = prodAmpIndices[0].size() + prodAmpIndices[1].size();
533  TMatrixT<double> jacobian(1, 2 * nmbAmps);
534  // precalculate some variables
535  const double rhoAA = spinDensityMatrixElem(waveIndexA, waveIndexA).real(); // rho_AA is real by definition
536  const double rhoBB = spinDensityMatrixElem(waveIndexB, waveIndexB).real(); // rho_BB is real by definition
537  const complex<double> rhoAB = spinDensityMatrixElem(waveIndexA, waveIndexB);
538  const double rhoABRe = rhoAB.real();
539  const double rhoABIm = rhoAB.imag();
540  const double rhoABNorm = std::norm(rhoAB);
541  const double coh = sqrt(rhoABNorm / (rhoAA * rhoBB));
542  if (coh == 0)
543  return 0;
544  // build m sub-Jacobians for JA_r = coh_AB / d V_Ar
545  for (unsigned int i = 0; i < prodAmpIndices[0].size(); ++i) {
546  const unsigned int prodAmpIndexA = prodAmpIndices[0][i];
547  const complex<double> prodAmpA = prodAmp(prodAmpIndexA);
548  const int prodAmpRankA = rankOfProdAmp(prodAmpIndexA);
549  // find production amplitude of wave B with same rank
550  complex<double> prodAmpB = 0;
551  for (unsigned int j = 0; j < prodAmpIndices[1].size(); ++j) {
552  const unsigned int prodAmpIndexB = prodAmpIndices[1][j];
553  if (rankOfProdAmp(prodAmpIndexB) == prodAmpRankA) {
554  prodAmpB = prodAmp(prodAmpIndexB);
555  break;
556  }
557  }
558  jacobian[0][2 * i ] = rhoABRe * prodAmpB.real() - rhoABIm * prodAmpB.imag() - (rhoABNorm / rhoAA) * prodAmpA.real();
559  jacobian[0][2 * i + 1] = rhoABRe * prodAmpB.imag() + rhoABIm * prodAmpB.real() - (rhoABNorm / rhoAA) * prodAmpA.imag();
560  }
561  // !!! possible optimization: join the loops for JA_r and JB_r
562  // build m sub-Jacobian for JB_r = d coh_AB / d V_Br
563  const unsigned int colOffset = 2 * prodAmpIndices[0].size();
564  for (unsigned int i = 0; i < prodAmpIndices[1].size(); ++i) {
565  const unsigned int prodAmpIndexB = prodAmpIndices[1][i];
566  const complex<double> prodAmpB = prodAmp(prodAmpIndexB);
567  const int prodAmpRankB = rankOfProdAmp(prodAmpIndexB);
568  // find production amplitude of wave A with same rank
569  complex<double> prodAmpA = 0;
570  for (unsigned int j = 0; j < prodAmpIndices[0].size(); ++j) {
571  const unsigned int prodAmpIndexA = prodAmpIndices[0][j];
572  if (rankOfProdAmp(prodAmpIndexA) == prodAmpRankB) {
573  prodAmpA = prodAmp(prodAmpIndexA);
574  break;
575  }
576  }
577  jacobian[0][colOffset + 2 * i ] = rhoABRe * prodAmpA.real() + rhoABIm * prodAmpA.imag() - (rhoABNorm / rhoBB) * prodAmpB.real();
578  jacobian[0][colOffset + 2 * i + 1] = rhoABRe * prodAmpA.imag() - rhoABIm * prodAmpA.real() - (rhoABNorm / rhoBB) * prodAmpB.imag();
579  }
580  jacobian *= 1 / (coh * rhoAA * rhoBB);
581  // build covariance matrix for amplitudes and calculate coherence covariance matrix
582  const TMatrixT<double> prodAmpCov = this->prodAmpCov(prodAmpIndices[0], prodAmpIndices[1]); // 2(n + m) x 2(n + m) matrix
583  const TMatrixT<double> jacobianT(TMatrixT<double>::kTransposed, jacobian); // 2(n + m) x 1 matrix
584  const TMatrixT<double> prodAmpCovJT = prodAmpCov * jacobianT; // 2(n + m) x 1 matrix
585  const TMatrixT<double> cohCov = jacobian * prodAmpCovJT; // 1 x 1 matrix
586  return sqrt(cohCov[0][0]);
587 }
588 
589 
591 double
592 fitResult::overlap(const unsigned int waveIndexA,
593  const unsigned int waveIndexB) const
594 {
595  const complex<double> spinDens = spinDensityMatrixElem(waveIndexA, waveIndexB);
596  const complex<double> normInt = normIntegral (waveIndexA, waveIndexB);
597  return 2 * (spinDens * normInt).real();
598 }
599 
600 
602 double
603 fitResult::overlapErr(const unsigned int waveIndexA,
604  const unsigned int waveIndexB) const
605 {
606  if (!_covMatrixValid)
607  return 0;
608  const complex<double> normInt = normIntegral(waveIndexA, waveIndexB);
609  TMatrixT<double> jacobian(1, 2); // overlap is real valued function, so J has only one row
610  jacobian[0][0] = 2 * normInt.real();
611  jacobian[0][1] = -2 * normInt.imag();
612  const double overlapVariance = realValVariance(waveIndexA, waveIndexB, jacobian);
613  return sqrt(overlapVariance);
614 }
615 
616 
617 void
618 fitResult::reset()
619 {
620  _nmbEvents = 0;
621  _normNmbEvents = 0;
622  _massBinCenter = 0;
623  _logLikelihood = 0;
624  _rank = 0;
625  _prodAmps.clear();
626  _prodAmpNames.clear();
627  _waveNames.clear();
628  _covMatrixValid = false;
629  _fitParCovMatrix.ResizeTo(0, 0);
630  _fitParCovMatrixIndices.clear();
631  _normIntegral.ResizeTo(0, 0);
632  _normIntIndexMap.clear();
633  _phaseSpaceIntegral.clear();
634  _converged = false;
635  _hasHessian = false;
636 }
637 
638 
639 void
640 fitResult::fill
641 (const unsigned int nmbEvents, // number of events in bin
642  const unsigned int normNmbEvents, // number of events to normalize to
643  const double massBinCenter, // center value of mass bin
644  const double logLikelihood, // log(likelihood) at maximum
645  const int rank, // rank of fit
646  const vector<complex<double> >& prodAmps, // production amplitudes
647  const vector<string>& prodAmpNames, // names of production amplitudes used in fit
648  const TMatrixT<double>& fitParCovMatrix, // covariance matrix of fit parameters
649  const vector<pair<int, int> >& fitParCovMatrixIndices, // indices of fit parameters for real and imaginary part in covariance matrix matrix
650  const TCMatrix& normIntegral, // normalization integral matrix
651  const vector<double>& phaseSpaceIntegral, // normalization integral over full phase space without acceptance
652  const bool converged,
653  const bool hasHessian)
654 {
655  _converged = converged;
656  _hasHessian = hasHessian;
657  _nmbEvents = nmbEvents;
658  _normNmbEvents = normNmbEvents;
659  _massBinCenter = massBinCenter;
660  _logLikelihood = logLikelihood;
661  _rank = rank;
662  _prodAmps.resize(prodAmps.size());
663  for (unsigned int i = 0; i < prodAmps.size(); ++i)
664  _prodAmps[i] = TComplex(prodAmps[i].real(), prodAmps[i].imag());
665  _prodAmpNames = prodAmpNames;
666  _fitParCovMatrix.ResizeTo(fitParCovMatrix.GetNrows(), fitParCovMatrix.GetNcols());
667  _fitParCovMatrix = fitParCovMatrix;
668  _fitParCovMatrixIndices = fitParCovMatrixIndices;
669  // check whether there really is an error matrix
670  if (!(fitParCovMatrix.GetNrows() == 0) && !(fitParCovMatrix.GetNcols() == 0))
671  _covMatrixValid = true;
672  else
673  _covMatrixValid = false;
674  _normIntegral.ResizeTo(normIntegral.nrows(), normIntegral.ncols());
675  _normIntegral = normIntegral;
676  _phaseSpaceIntegral = phaseSpaceIntegral;
677 
678  buildWaveMap();
679 
680  // check consistency
681  if (_prodAmps.size() != _prodAmpNames.size())
682  cout << "fitResult::fill(): warning: number of production amplitudes "
683  << "(" << _prodAmps.size() << ") does not match number of "
684  << "production amplitude names (" << _prodAmpNames.size() << ")." << endl;
685  if (_prodAmps.size() != _fitParCovMatrixIndices.size())
686  cout << "fitResult::fill(): warning: number of production amplitudes "
687  << "(" << _prodAmps.size() << ") does not match number of "
688  << "covariance matrix indices (" << _fitParCovMatrixIndices.size() << ")." << endl;
689  if ( ((int)_waveNames.size() != _normIntegral.nrows())
690  || ((int)_waveNames.size() != _normIntegral.ncols()))
691  cout << "fitResult::fill(): warning: number of waves (" << _waveNames.size()
692  << ") does not match size of normalization integral "
693  << "(" << _normIntegral.nrows() << ", " << _normIntegral.ncols() << ")." << endl;
694 
695  if (0) {
696  // print debug information that allows to check ordering of arrays
697  map<TString, complex<double> > V;
698  for (unsigned int i = 0; i < nmbProdAmps(); ++i)
699  V[prodAmpName(i)] = prodAmp(i);
700  printInfo << "production amplitudes:" << endl;
701  for (map<TString, complex<double> >::const_iterator i = V.begin(); i != V.end(); ++i)
702  cout << " " << i->first << " = " << i->second << endl;
703  printInfo << "normalization integral matrix:" << endl;
704  map<TString, unsigned int> waveIndexMap;
705  for (unsigned int i = 0; i < nmbWaves(); ++i)
706  waveIndexMap[waveName(i)] = i;
707  for (map<TString, unsigned int>::const_iterator i = waveIndexMap.begin();
708  i != waveIndexMap.end(); ++i)
709  for (map<TString, unsigned int>::const_iterator j = waveIndexMap.begin();
710  j != waveIndexMap.end(); ++j) {
711  cout << " [" << i->first << "][" << j->first << "] = "
712  << this->normIntegral(i->second, j->second) << endl;
713  }
714  printInfo << "covariance matrix:" << endl;
715  map<TString, unsigned int> prodAmpIndexMap;
716  for (unsigned int i = 0; i < nmbProdAmps(); ++i)
717  prodAmpIndexMap[prodAmpName(i)] = i;
718  for (map<TString, unsigned int>::const_iterator i = prodAmpIndexMap.begin();
719  i != prodAmpIndexMap.end(); ++i) {
720  const int iIdx[2] = { _fitParCovMatrixIndices[i->second].first, // real part index
721  _fitParCovMatrixIndices[i->second].second}; // imaginary part index
722  for (map<TString, unsigned int>::const_iterator j = prodAmpIndexMap.begin();
723  j != prodAmpIndexMap.end(); ++j) {
724  const int jIdx[2] = { _fitParCovMatrixIndices[j->second].first, // real part index
725  _fitParCovMatrixIndices[j->second].second}; // imaginary part index
726  cout << " [" << i->first << "][" << j->first << "]: " << flush;
727  if (iIdx[0] >= 0) {
728  if (jIdx[0] >= 0)
729  cout << "ReRe = " << fitParameterCov(iIdx[0], jIdx[0]) << ", " << flush;
730  if (jIdx[1] >= 0)
731  cout << "ReIm = " << fitParameterCov(iIdx[0], jIdx[1]) << ", " << flush;
732  }
733  if (iIdx[1] >= 0) {
734  if (jIdx[0] >= 0)
735  cout << "ImRe = " << fitParameterCov(iIdx[1], jIdx[0]) << ", " << flush;
736  if (jIdx[1] >= 0)
737  cout << "ImIm = " << fitParameterCov(iIdx[1], jIdx[1]) << flush;
738  }
739  cout << endl;
740  }
741  }
742  }
743 }
744 
745 
746 int
747 fitResult::waveIndex(const string& waveName) const
748 {
749  int index = -1;
750  for (unsigned int i = 0; i < nmbWaves(); ++i)
751  if (waveName == _waveNames[i]) {
752  index = i;
753  break; // assumes that waves have unique names
754  }
755  if (index == -1)
756  printWarn << "could not find any wave named '" << waveName << "'." << endl;
757  return index;
758 }
759 
760 
761 int
762 fitResult::prodAmpIndex(const string& prodAmpName) const
763 {
764  int index = -1;
765  for (unsigned int i = 0; i < nmbProdAmps(); ++i)
766  if (prodAmpName == _prodAmpNames[i]) {
767  index = i;
768  break; // assumes that production amplitudes have unique names
769  }
770  if (index == -1)
771  printWarn << "could not find any production amplitude named '" << prodAmpName << "'." << endl;
772  return index;
773 }
774 
775 
776 void
777 fitResult::buildWaveMap()
778 {
779  const int n = _prodAmpNames.size();
780  for (int i = 0; i < n; ++i) {
781  // strip rank
782  const TString title = wavetitle(i);
783  if (find(_waveNames.begin(), _waveNames.end(), title.Data()) == _waveNames.end())
784  _waveNames.push_back(title.Data());
785 
786  // look for index of first occurence
787  int j;
788  for (j = 0; j < n; ++j)
789  if(prodAmpName(j).Contains(title))
790  break;
791  _normIntIndexMap[i] = j;
792  }
793 }
Generated by   doxygen 1.8.1.2