ROOTPWA
TPWALikelihood.cc
Go to the documentation of this file.
1 
2 //
3 // Copyright 2009 Sebastian Neubert, Boris Grube
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 // Implementation of class TPWALikelihood
27 // see TPWALikelihood.h for details
28 //
29 //
30 // Author List:
31 // Sebastian Neubert TUM (original author)
32 // Boris Grube Universe Cluster Munich
33 //
34 //
35 //-----------------------------------------------------------
36 
37 
38 #include <iomanip>
39 #include <fstream>
40 #include <complex>
41 #include <cassert>
42 #include <limits>
43 
44 #include "TString.h"
45 #include "TSystem.h"
46 #include "TCMatrix.h"
47 #include "TStopwatch.h"
48 #include "TTree.h"
49 
50 #include "reportingUtils.hpp"
51 #include "conversionUtils.hpp"
52 #include "fileUtils.hpp"
53 #ifdef USE_CUDA
54 #include "complex.cuh"
55 #include "likelihoodInterface.cuh"
56 #endif
57 #include "amplitudeTreeLeaf.h"
58 #include "TPWALikelihood.h"
59 
60 
61 // #define USE_FDF
62 
63 
64 using namespace std;
65 using namespace rpwa;
66 using namespace boost;
67 using namespace boost::accumulators;
68 
69 
70 template<typename complexT> bool TPWALikelihood<complexT>::_debug = true;
71 
72 
73 template<typename complexT>
74 TPWALikelihood<complexT>::TPWALikelihood()
75  : _nmbEvents (0),
76  _rank (1),
77  _nmbWaves (0),
78  _nmbWavesReflMax (0),
79  _nmbPars (0),
80 #ifdef USE_CUDA
81  _cudaEnabled (false),
82 #endif
83  _useNormalizedAmps(true),
84  _numbAccEvents (0)
85 {
86  _nmbWavesRefl[0] = 0;
87  _nmbWavesRefl[1] = 0;
88  resetFuncCallInfo();
89 #ifdef USE_FDF
90  printInfo << "using FdF() to calculate likelihood" << endl;
91 #else
92  printInfo << "using DoEval() to calculate likelihood" << endl;
93 #endif
94 }
95 
96 
97 template<typename complexT>
98 TPWALikelihood<complexT>::~TPWALikelihood()
99 {
100  clear();
101 }
102 
103 
104 template<typename complexT>
105 void
106 TPWALikelihood<complexT>::FdF
107 (const double* par, // parameter array; reduced by rank conditions
108  double& funcVal, // function value
109  double* gradient) const // array of derivatives
110 {
111  ++(_funcCallInfo[FDF].nmbCalls);
112 
113  // timer for total time
114  TStopwatch timerTot;
115  timerTot.Start();
116 
117  // copy arguments into parameter cache
118  for (unsigned int i = 0; i < _nmbPars; ++i)
119  _parCache[i] = par[i];
120 
121  // build complex production amplitudes from function parameters taking into account rank restrictions
122  value_type prodAmpFlat;
123  ampsArrayType prodAmps;
124  copyFromParArray(par, prodAmps, prodAmpFlat);
125  const value_type prodAmpFlat2 = prodAmpFlat * prodAmpFlat;
126 
127  // create array of likelihood derivative w.r.t. real and imaginary
128  // parts of the production amplitudes
129  // !NOTE! although stored as and constructed from complex values,
130  // the dL themselves are _not_ well defined complex numbers!
131  value_type derivativeFlat = 0;
132  array<typename ampsArrayType::index, 3> derivShape = {{ _rank, 2, _nmbWavesReflMax }};
133  ampsArrayType derivatives(derivShape);
134 
135  // loop over events and calculate real-data term of log likelihood
136  // as well as derivatives with respect to parameters
137  TStopwatch timer;
138  timer.Start();
139  accumulator_set<value_type, stats<tag::sum(compensated)> > logLikelihoodAcc;
140  accumulator_set<value_type, stats<tag::sum(compensated)> > derivativeFlatAcc;
141  multi_array<accumulator_set<complexT, stats<tag::sum(compensated)> >, 3>
142  derivativesAcc(derivShape);
143  for (unsigned int iEvt = 0; iEvt < _nmbEvents; ++iEvt) {
144  accumulator_set<value_type, stats<tag::sum(compensated)> > likelihoodAcc;
145  ampsArrayType derivative(derivShape); // likelihood derivative for this event
146  for (unsigned int iRank = 0; iRank < _rank; ++iRank) { // incoherent sum over ranks
147  for (unsigned int iRefl = 0; iRefl < 2; ++iRefl) { // incoherent sum over reflectivities
148  accumulator_set<complexT, stats<tag::sum(compensated)> > ampProdAcc;
149  for (unsigned int iWave = 0; iWave < _nmbWavesRefl[iRefl]; ++iWave) // coherent sum over waves
150  ampProdAcc(prodAmps[iRank][iRefl][iWave] * _decayAmps[iEvt][iRefl][iWave]);
151  const complexT ampProdSum = sum(ampProdAcc);
152  likelihoodAcc(norm(ampProdSum));
153  // set derivative term that is independent on derivative wave index
154  for (unsigned int jWave = 0; jWave < _nmbWavesRefl[iRefl]; ++jWave)
155  // amplitude sums for current rank and for waves with same reflectivity
156  derivative[iRank][iRefl][jWave] = ampProdSum;
157  }
158  // loop again over waves for current rank and multiply with complex conjugate
159  // of decay amplitude of the wave with the derivative wave index
160  for (unsigned int iRefl = 0; iRefl < 2; ++iRefl)
161  for (unsigned int iWave = 0; iWave < _nmbWavesRefl[iRefl]; ++iWave)
162  derivative[iRank][iRefl][iWave] *= conj(_decayAmps[iEvt][iRefl][iWave]);
163  } // end loop over rank
164  likelihoodAcc (prodAmpFlat2 );
165  logLikelihoodAcc(-log(sum(likelihoodAcc)));
166  // incorporate factor 2 / sigma
167  const value_type factor = 2. / sum(likelihoodAcc);
168  for (unsigned int iRank = 0; iRank < _rank; ++iRank)
169  for (unsigned int iRefl = 0; iRefl < 2; ++iRefl)
170  for (unsigned int iWave = 0; iWave < _nmbWavesRefl[iRefl]; ++iWave)
171  derivativesAcc[iRank][iRefl][iWave](-factor * derivative[iRank][iRefl][iWave]);
172  derivativeFlatAcc(-factor * prodAmpFlat);
173  } // end loop over events
174  for (unsigned int iRank = 0; iRank < _rank; ++iRank)
175  for (unsigned int iRefl = 0; iRefl < 2; ++iRefl)
176  for (unsigned int iWave = 0; iWave < _nmbWavesRefl[iRefl]; ++iWave)
177  derivatives[iRank][iRefl][iWave] = sum(derivativesAcc[iRank][iRefl][iWave]);
178  derivativeFlat = sum(derivativeFlatAcc);
179  // log time needed for likelihood calculation
180  timer.Stop();
181  _funcCallInfo[FDF].funcTime(timer.RealTime());
182 
183  // compute normalization term of log likelihood and normalize derivatives w.r.t. parameters
184  timer.Start();
185  accumulator_set<value_type, stats<tag::sum(compensated)> > normFactorAcc;
186  const value_type nmbEvt = (_useNormalizedAmps) ? 1 : _nmbEvents;
187  const value_type twiceNmbEvt = 2 * nmbEvt;
188  for (unsigned int iRank = 0; iRank < _rank; ++iRank)
189  for (unsigned int iRefl = 0; iRefl < 2; ++iRefl)
190  for (unsigned int iWave = 0; iWave < _nmbWavesRefl[iRefl]; ++iWave) {
191  accumulator_set<complexT, stats<tag::sum(compensated)> > normFactorDerivAcc;
192  for (unsigned int jWave = 0; jWave < _nmbWavesRefl[iRefl]; ++jWave) { // inner loop over waves with same reflectivity
193  const complexT I = _accMatrix[iRefl][iWave][iRefl][jWave];
194  normFactorAcc(real((prodAmps[iRank][iRefl][iWave] * conj(prodAmps[iRank][iRefl][jWave]))
195  * I));
196  normFactorDerivAcc(prodAmps[iRank][iRefl][jWave] * conj(I));
197  }
198  derivatives[iRank][iRefl][iWave] += sum(normFactorDerivAcc) * twiceNmbEvt; // account for 2 * nmbEvents
199  }
200  // take care of flat wave
201  normFactorAcc(prodAmpFlat2 * _totAcc);
202  derivativeFlat += prodAmpFlat * twiceNmbEvt * _totAcc;
203  // log time needed for normalization
204  timer.Stop();
205  _funcCallInfo[FDF].normTime(timer.RealTime());
206 
207  // sort derivative results into output array and cache
208  copyToParArray(derivatives, derivativeFlat, gradient);
209  copyToParArray(derivatives, derivativeFlat, toArray(_derivCache));
210 
211  // set function return value
212  funcVal = sum(logLikelihoodAcc) + nmbEvt * sum(normFactorAcc);
213 
214  // log total consumed time
215  timerTot.Stop();
216  _funcCallInfo[FDF].totalTime(timerTot.RealTime());
217 
218  if (_debug)
219  printDebug << "log likelihood = " << maxPrecisionAlign(sum(logLikelihoodAcc)) << ", "
220  << "normalization = " << maxPrecisionAlign(sum(normFactorAcc) ) << ", "
221  << "normalized likelihood = " << maxPrecisionAlign(funcVal ) << endl;
222 }
223 
224 
225 template<typename complexT>
226 double
227 TPWALikelihood<complexT>::DoEval(const double* par) const
228 {
229  ++(_funcCallInfo[DOEVAL].nmbCalls);
230 
231 #ifdef USE_FDF
232 
233  // call FdF
234  double logLikelihood;
235  double gradient[_nmbPars];
236  FdF(par, logLikelihood, gradient);
237  return logLikelihood;
238 
239 #else // USE_FDF
240 
241  // timer for total time
242  TStopwatch timerTot;
243  timerTot.Start();
244 
245  // build complex production amplitudes from function parameters taking into account rank restrictions
246  value_type prodAmpFlat;
247  ampsArrayType prodAmps;
248  copyFromParArray(par, prodAmps, prodAmpFlat);
249  const value_type prodAmpFlat2 = prodAmpFlat * prodAmpFlat;
250 
251  // loop over events and calculate real-data term of log likelihood
252  TStopwatch timer;
253  timer.Start();
254  value_type logLikelihood = 0;
255 #ifdef USE_CUDA
256  if (_cudaEnabled) {
257  logLikelihood = cuda::likelihoodInterface<cuda::complex<value_type> >::logLikelihood
258  (reinterpret_cast<cuda::complex<value_type>*>(prodAmps.data()),
259  prodAmps.num_elements(), prodAmpFlat, _rank);
260  } else
261 #endif
262  {
263  accumulator_set<value_type, stats<tag::sum(compensated)> > logLikelihoodAcc;
264  for (unsigned int iEvt = 0; iEvt < _nmbEvents; ++iEvt) {
265  accumulator_set<value_type, stats<tag::sum(compensated)> > likelihoodAcc;
266  for (unsigned int iRank = 0; iRank < _rank; ++iRank) { // incoherent sum over ranks
267  for (unsigned int iRefl = 0; iRefl < 2; ++iRefl) { // incoherent sum over reflectivities
268  accumulator_set<complexT, stats<tag::sum(compensated)> > ampProdAcc;
269  for (unsigned int iWave = 0; iWave < _nmbWavesRefl[iRefl]; ++iWave) { // coherent sum over waves
270  ampProdAcc(prodAmps[iRank][iRefl][iWave] * _decayAmps[iEvt][iRefl][iWave]);
271  // cout << "prodAmps[" << iRank << "][" << iRefl << "][" << iWave<< "] = "
272  // << maxPrecisionDouble(prodAmps[iRank][iRefl][iWave]) << "; "
273  // << "decayAmps[" << iEvt << "][" << iRefl << "][" << iWave << "] = "
274  // << maxPrecisionDouble(_decayAmps[iEvt][iRefl][iWave]) << endl;
275  }
276  likelihoodAcc(norm(sum(ampProdAcc)));
277  }
278  }
279  likelihoodAcc(prodAmpFlat2);
280  logLikelihoodAcc(-log(sum(likelihoodAcc)));
281  // cout << endl;
282  }
283  logLikelihood = sum(logLikelihoodAcc);
284  assert(logLikelihood==logLikelihood);
285  if(numeric_limits<double>::has_infinity)assert(logLikelihood != numeric_limits<double>::infinity());
286  }
287  // log time needed for likelihood calculation
288  timer.Stop();
289  _funcCallInfo[DOEVAL].funcTime(timer.RealTime());
290 
291  // compute normalization term of log likelihood
292  timer.Start();
293  accumulator_set<value_type, stats<tag::sum(compensated)> > normFactorAcc;
294  const value_type nmbEvt = (_useNormalizedAmps) ? 1 : _nmbEvents;
295  for (unsigned int iRank = 0; iRank < _rank; ++iRank)
296  for (unsigned int iRefl = 0; iRefl < 2; ++iRefl)
297  for (unsigned int iWave = 0; iWave < _nmbWavesRefl[iRefl]; ++iWave)
298  for (unsigned int jWave = 0; jWave < _nmbWavesRefl[iRefl]; ++jWave) // loop over waves with same reflectivity
299  normFactorAcc(real((prodAmps[iRank][iRefl][iWave] * conj(prodAmps[iRank][iRefl][jWave]))
300  * _accMatrix[iRefl][iWave][iRefl][jWave]));
301  // take care of flat wave
302  normFactorAcc(prodAmpFlat2 * _totAcc);
303  // log time needed for normalization
304  timer.Stop();
305  _funcCallInfo[DOEVAL].normTime(timer.RealTime());
306 
307  // calculate and return log likelihood value
308  const double funcVal = logLikelihood + nmbEvt * sum(normFactorAcc);
309  assert(funcVal==funcVal);
310  if(numeric_limits<double>::has_infinity)assert(funcVal != numeric_limits<double>::infinity());
311 
312 
313  // log total consumed time
314  timerTot.Stop();
315  _funcCallInfo[DOEVAL].totalTime(timerTot.RealTime());
316 
317  if (_debug)
318  printDebug << "raw log likelihood = " << maxPrecisionAlign(logLikelihood ) << ", "
319  << "normalization = " << maxPrecisionAlign(sum(normFactorAcc)) << ", "
320  << "normalized log likelihood = " << maxPrecisionAlign(funcVal ) << endl;
321 
322  return funcVal;
323 
324 #endif // USE_FDF
325 }
326 
327 
328 template<typename complexT>
329 double
330 TPWALikelihood<complexT>::DoDerivative(const double* par,
331  unsigned int derivativeIndex) const
332 {
333  ++(_funcCallInfo[DODERIVATIVE].nmbCalls);
334 
335  // timer for total time
336  TStopwatch timerTot;
337  timerTot.Start();
338 
339  // check whether parameter is in cache
340  bool samePar = true;
341  for (unsigned int i = 0; i < _nmbPars; ++i)
342  if (_parCache[i] != par[i]) {
343  samePar = false;
344  break;
345  }
346  if (samePar) {
347  //cout << "using cached derivative! " << endl;
348  timerTot.Stop();
349  _funcCallInfo[DODERIVATIVE].totalTime(timerTot.RealTime());
350  _funcCallInfo[DODERIVATIVE].funcTime (sum(_funcCallInfo[DODERIVATIVE].totalTime));
351  return _derivCache[derivativeIndex];
352  }
353  // call FdF
354  double logLikelihood;
355  double gradient[_nmbPars];
356  FdF(par, logLikelihood, gradient);
357  return gradient[derivativeIndex];
358 }
359 
360 
361 // calculate derivatives with respect to parameters
362 template<typename complexT>
363 void
364 TPWALikelihood<complexT>::Gradient
365 (const double* par, // parameter array; reduced by rank conditions
366  double* gradient) const // array of derivatives
367 {
368  ++(_funcCallInfo[GRADIENT].nmbCalls);
369 
370  // timer for total time
371  TStopwatch timerTot;
372  timerTot.Start();
373 
374 #ifdef USE_FDF
375 
376  // check whether parameter is in cache
377  bool samePar = true;
378  for (unsigned int i = 0; i < _nmbPars; ++i)
379  if (_parCache[i] != par[i]) {
380  samePar = false;
381  break;
382  }
383  if (samePar) {
384  for (unsigned int i = 0; i < _nmbPars ; ++i)
385  gradient[i] = _derivCache[i];
386  timerTot.Stop();
387  _funcCallInfo[GRADIENT].totalTime(timerTot.RealTime());
388  _funcCallInfo[GRADIENT].funcTime (sum(_funcCallInfo[GRADIENT].totalTime));
389  return;
390  }
391  // call FdF
392  double logLikelihood;
393  FdF(par, logLikelihood, gradient);
394 
395 #else // USE_FDF
396 
397  // copy arguments into parameter cache
398  for (unsigned int i = 0; i < _nmbPars; ++i)
399  _parCache[i] = par[i];
400 
401  // build complex production amplitudes from function parameters taking into account rank restrictions
402  value_type prodAmpFlat;
403  ampsArrayType prodAmps;
404  copyFromParArray(par, prodAmps, prodAmpFlat);
405 
406  // create array of likelihood derivative w.r.t. real and imaginary
407  // parts of the production amplitudes
408  // !NOTE! although stored as and constructed from complex values,
409  // the dL themselves are _not_ well defined complex numbers!
410  value_type derivativeFlat = 0;
411  array<typename ampsArrayType::index, 3> derivShape = {{ _rank, 2, _nmbWavesReflMax }};
412  ampsArrayType derivatives(derivShape);
413 
414  // loop over events and calculate derivatives with respect to parameters
415  TStopwatch timer;
416  timer.Start();
417 #ifdef USE_CUDA
418  if (_cudaEnabled) {
419  cuda::likelihoodInterface<cuda::complex<value_type> >::logLikelihoodDeriv
420  (reinterpret_cast<cuda::complex<value_type>*>(prodAmps.data()),
421  prodAmps.num_elements(), prodAmpFlat, _rank,
422  reinterpret_cast<cuda::complex<value_type>*>(derivatives.data()),
423  derivativeFlat);
424  } else
425 #endif
426  {
427  accumulator_set<value_type, stats<tag::sum(compensated)> > derivativeFlatAcc;
428  multi_array<accumulator_set<complexT, stats<tag::sum(compensated)> >, 3>
429  derivativesAcc(derivShape);
430  const value_type prodAmpFlat2 = prodAmpFlat * prodAmpFlat;
431  for (unsigned int iEvt = 0; iEvt < _nmbEvents; ++iEvt) {
432  accumulator_set<value_type, stats<tag::sum(compensated)> > likelihoodAcc;
433  ampsArrayType derivative(derivShape); // likelihood derivatives for this event
434  for (unsigned int iRank = 0; iRank < _rank; ++iRank) { // incoherent sum over ranks
435  for (unsigned int iRefl = 0; iRefl < 2; ++iRefl) { // incoherent sum over reflectivities
436  accumulator_set<complexT, stats<tag::sum(compensated)> > ampProdAcc;
437  for (unsigned int iWave = 0; iWave < _nmbWavesRefl[iRefl]; ++iWave) // coherent sum over waves
438  ampProdAcc(prodAmps[iRank][iRefl][iWave] * _decayAmps[iEvt][iRefl][iWave]);
439  const complexT ampProdSum = sum(ampProdAcc);
440  likelihoodAcc(norm(ampProdSum));
441  // set derivative term that is independent on derivative wave index
442  for (unsigned int jWave = 0; jWave < _nmbWavesRefl[iRefl]; ++jWave)
443  // amplitude sums for current rank and for waves with same reflectivity
444  derivative[iRank][iRefl][jWave] = ampProdSum;
445  }
446  // loop again over waves for current rank and multiply with complex conjugate
447  // of decay amplitude of the wave with the derivative wave index
448  for (unsigned int iRefl = 0; iRefl < 2; ++iRefl)
449  for (unsigned int jWave = 0; jWave < _nmbWavesRefl[iRefl]; ++jWave)
450  derivative[iRank][iRefl][jWave] *= conj(_decayAmps[iEvt][iRefl][jWave]);
451  } // end loop over rank
452  likelihoodAcc(prodAmpFlat2);
453  // incorporate factor 2 / sigma
454  const value_type factor = 2. / sum(likelihoodAcc);
455  for (unsigned int iRank = 0; iRank < _rank; ++iRank)
456  for (unsigned int iRefl = 0; iRefl < 2; ++iRefl)
457  for (unsigned int jWave = 0; jWave < _nmbWavesRefl[iRefl]; ++jWave)
458  derivativesAcc[iRank][iRefl][jWave](-factor * derivative[iRank][iRefl][jWave]);
459  derivativeFlatAcc(-factor * prodAmpFlat);
460  } // end loop over events
461  for (unsigned int iRank = 0; iRank < _rank; ++iRank)
462  for (unsigned int iRefl = 0; iRefl < 2; ++iRefl)
463  for (unsigned int iWave = 0; iWave < _nmbWavesRefl[iRefl]; ++iWave)
464  derivatives[iRank][iRefl][iWave] = sum(derivativesAcc[iRank][iRefl][iWave]);
465  derivativeFlat = sum(derivativeFlatAcc);
466  }
467  // log time needed for gradient calculation
468  timer.Stop();
469  _funcCallInfo[GRADIENT].funcTime(timer.RealTime());
470 
471  // normalize derivatives w.r.t. parameters
472  timer.Start();
473  const value_type nmbEvt = (_useNormalizedAmps) ? 1 : _nmbEvents;
474  const value_type twiceNmbEvt = 2 * nmbEvt;
475  for (unsigned int iRank = 0; iRank < _rank; ++iRank)
476  for (unsigned int iRefl = 0; iRefl < 2; ++iRefl)
477  for (unsigned int iWave = 0; iWave < _nmbWavesRefl[iRefl]; ++iWave) {
478  accumulator_set<complexT, stats<tag::sum(compensated)> > normFactorDerivAcc;
479  for (unsigned int jWave = 0; jWave < _nmbWavesRefl[iRefl]; ++jWave) { // inner loop over waves with same reflectivity
480  const complexT I = _accMatrix[iRefl][iWave][iRefl][jWave];
481  normFactorDerivAcc(prodAmps[iRank][iRefl][jWave] * conj(I));
482  }
483  derivatives[iRank][iRefl][iWave] += sum(normFactorDerivAcc) * twiceNmbEvt; // account for 2 * nmbEvents
484  }
485  // take care of flat wave
486  derivativeFlat += prodAmpFlat * twiceNmbEvt * _totAcc;
487  // log time needed for normalization
488  timer.Stop();
489  _funcCallInfo[GRADIENT].normTime(timer.RealTime());
490 
491  // set return gradient values
492  copyToParArray(derivatives, derivativeFlat, gradient);
493  copyToParArray(derivatives, derivativeFlat, toArray(_derivCache));
494 
495  // log total consumed time
496  timerTot.Stop();
497  _funcCallInfo[GRADIENT].totalTime(timerTot.RealTime());
498 
499 #endif // USE_FDF
500 }
501 
502 
503 template<typename complexT>
504 unsigned int
505 TPWALikelihood<complexT>::nmbWaves(const int reflectivity) const
506 {
507  if (reflectivity == 0)
508  return _nmbWaves;
509  else if (reflectivity > 0)
510  return _nmbWavesRefl[1]; // positive reflectivity
511  else
512  return _nmbWavesRefl[0]; // negative reflectivity
513 }
514 
515 
516 template<typename complexT>
517 void
518 #ifdef USE_CUDA
519 TPWALikelihood<complexT>::enableCuda(const bool enableCuda)
520 {
521  _cudaEnabled = enableCuda;
522 }
523 #else
524 TPWALikelihood<complexT>::enableCuda(const bool) { }
525 #endif
526 
527 
528 template<typename complexT>
529 bool
530 TPWALikelihood<complexT>::cudaEnabled() const
531 {
532 #ifdef USE_CUDA
533  return _cudaEnabled;
534 #else
535  return false;
536 #endif
537 }
538 
539 
540 template<typename complexT>
541 void
542 TPWALikelihood<complexT>::init(const unsigned int rank,
543  const std::string& waveListFileName,
544  const std::string& normIntFileName,
545  const std::string& accIntFileName,
546  const std::string& ampDirName,
547  const unsigned int numbAccEvents,
548  const bool useRootAmps)
549 {
550  _numbAccEvents = numbAccEvents;
551  readWaveList(waveListFileName);
552  buildParDataStruct(rank);
553  readIntegrals(normIntFileName, accIntFileName);
554  readDecayAmplitudes(ampDirName, useRootAmps);
555 #ifdef USE_CUDA
556  if (_cudaEnabled)
557  cuda::likelihoodInterface<cuda::complex<value_type> >::init
558  (reinterpret_cast<cuda::complex<value_type>*>(_decayAmps.data()),
559  _decayAmps.num_elements(), _nmbEvents, _nmbWavesRefl, true);
560 #endif
561 }
562 
563 
564 template<typename complexT>
565 void
566 TPWALikelihood<complexT>::readWaveList(const string& waveListFileName)
567 {
568  printInfo << "reading amplitude names and thresholds from wave list file "
569  << "'" << waveListFileName << "'." << endl;
570  ifstream waveListFile(waveListFileName.c_str());
571  if (not waveListFile) {
572  printErr << "cannot open file '" << waveListFileName << "'. aborting." << endl;
573  throw;
574  }
575  vector<string> waveNames [2];
576  vector<double> waveThresholds[2];
577  vector<unsigned int> waveIndices [2];
578  unsigned int countWave = 0;
579  unsigned int lineNmb = 0;
580  string line;
581  while (getline(waveListFile, line)) {
582  if (line[0] == '#') // comments start with #
583  continue;
584  stringstream lineStream;
585  lineStream.str(line);
586  string waveName;
587  if (lineStream >> waveName) {
588  double threshold;
589  // !!! it would be safer to make the threshold value in the wave list file mandatory
590  if (not (lineStream >> threshold))
591  threshold = 0;
592  if (_debug)
593  printDebug << "reading line " << setw(3) << lineNmb + 1 << ": " << waveName<< ", "
594  << "threshold = " << setw(4) << threshold << " MeV/c^2" << endl;
595  if (getReflectivity(waveName) > 0) {
596  ++_nmbWavesRefl[1]; // positive reflectivity
597  waveNames [1].push_back(waveName);
598  waveThresholds[1].push_back(threshold);
599  waveIndices [1].push_back(countWave);
600  } else {
601  ++_nmbWavesRefl[0]; // negative reflectivity
602  waveNames [0].push_back(waveName);
603  waveThresholds[0].push_back(threshold);
604  waveIndices [0].push_back(countWave);
605  }
606  ++countWave;
607  } else
608  printWarn << "cannot parse line '" << line << "' in wave list file "
609  << "'" << waveListFileName << "'" << endl;
610  ++lineNmb;
611  }
612  waveListFile.close();
613  printInfo << "read " << lineNmb << " lines from wave list file "
614  << "'" << waveListFileName << "'" << endl;
615  _nmbWaves = _nmbWavesRefl[0] + _nmbWavesRefl[1];
616  _nmbWavesReflMax = max(_nmbWavesRefl[0], _nmbWavesRefl[1]);
617  _waveNames.resize (extents[2][_nmbWavesReflMax]);
618  _waveThresholds.resize (extents[2][_nmbWavesReflMax]);
619  _waveToWaveIndex.resize(extents[2][_nmbWavesReflMax]);
620  for (unsigned int iRefl = 0; iRefl < 2; ++iRefl)
621  for (unsigned int iWave = 0; iWave < _nmbWavesRefl[iRefl]; ++iWave) {
622  _waveNames [iRefl][iWave] = waveNames [iRefl][iWave];
623  _waveThresholds [iRefl][iWave] = waveThresholds[iRefl][iWave];
624  _waveToWaveIndex[iRefl][iWave] = waveIndices [iRefl][iWave];
625  }
626 }
627 
628 
629 template<typename complexT>
630 void
631 TPWALikelihood<complexT>::buildParDataStruct(const unsigned int rank)
632 {
633  if ((_nmbWavesRefl[0] + _nmbWavesRefl[1] == 0) or (_waveThresholds.size() == 0)) {
634  printErr << "no wave info. was readWaveList() executed successfully? aborting.";
635  throw;
636  }
637  _rank = rank;
638  // calculate dimension of function taking into account rank restrictions and flat wave
639  _nmbPars = 0;
640  for (unsigned int iRank = 0; iRank < _rank; ++iRank) {
641  const int nmbProdAmpsPos = _nmbWavesRefl[1] - iRank; // number non-zero production amplitudes in this rank with positive reflectivity
642  int nmbProdAmpsPosC = nmbProdAmpsPos - 1; // number of complex-valued production amplitudes in this rank with positive reflectivity
643  if (nmbProdAmpsPosC < 0)
644  nmbProdAmpsPosC = 0;
645  const int nmbProdAmpsNeg = _nmbWavesRefl[0] - iRank; // number non-zero production amplitudes in this rank with negative reflectivity
646  int nmbProdAmpsNegC = nmbProdAmpsNeg - 1; // number of complex-valued production amplitudes in this rank with negative reflectivity
647  if (nmbProdAmpsNegC < 0)
648  nmbProdAmpsNegC = 0;
649  _nmbPars += 2 * (nmbProdAmpsPosC + nmbProdAmpsNegC); // 2 parameters for each complex production amplitude
650  // 1 real production amplitude per rank and reflectivity
651  if (nmbProdAmpsPos > 0)
652  ++_nmbPars;
653  if (nmbProdAmpsNeg > 0)
654  ++_nmbPars;
655  }
656  _nmbPars += 1; // additonal flat wave
657  printInfo << "dimension of likelihood function is " << _nmbPars << "." << endl;
658  _parNames.resize (_nmbPars, "");
659  _parThresholds.resize(_nmbPars, 0);
660  _parCache.resize (_nmbPars, 0);
661  _derivCache.resize (_nmbPars, 0);
662  _prodAmpToFuncParMap.resize(extents[_rank][2][_nmbWavesReflMax]);
663  // build parameter names
664  unsigned int parIndex = 0;
665  for (unsigned int iRank = 0; iRank < _rank; ++iRank)
666  for (unsigned int iRefl = 0; iRefl < 2; ++iRefl)
667  for (unsigned int iWave = 0; iWave < _nmbWavesRefl[iRefl]; ++iWave) {
668  ostringstream parName;
669  if (iWave < iRank) // production amplitude is zero
670  _prodAmpToFuncParMap[iRank][iRefl][iWave] = make_tuple(-1, -1);
671  else if (iWave == iRank) { // production amplitude is real
672  parName << "V" << iRank << "_" << _waveNames[iRefl][iWave] << "_RE";
673  _parNames [parIndex] = parName.str();
674  _parThresholds[parIndex] = _waveThresholds[iRefl][iWave];
675  _prodAmpToFuncParMap[iRank][iRefl][iWave] = make_tuple(parIndex, -1);
676  ++parIndex;
677  } else { // production amplitude is complex
678  parName << "V" << iRank << "_" << _waveNames[iRefl][iWave];
679  _parNames [parIndex] = parName.str() + "_RE";
680  _parThresholds[parIndex] = _waveThresholds[iRefl][iWave];
681  _parNames [parIndex + 1] = parName.str() + "_IM";
682  _parThresholds[parIndex + 1] = _waveThresholds[iRefl][iWave];
683  _prodAmpToFuncParMap[iRank][iRefl][iWave] = make_tuple(parIndex, parIndex + 1);
684  parIndex += 2;
685  }
686  }
687  // flat wave
688  _parNames [parIndex] = "V_flat";
689  _parThresholds[parIndex] = 0;
690 }
691 
692 
693 // returns integral matrix reordered according to _waveNames array
694 template<typename complexT>
695 void
696 TPWALikelihood<complexT>::reorderIntegralMatrix(integral& integral,
697  normMatrixArrayType& reorderedMatrix) const
698 {
699  // get original matrix and list of wave names
700  const matrix<complex<double> > intMatrix = integral.mat();
701  list<string> intWaveNames = integral.files();
702  // "int" saves filenames with path, this must be treated here
703  for (list<string>::iterator it = intWaveNames.begin(); it != intWaveNames.end(); ++it) {
704  // find the slash, if not available -> take the first position of the string
705  const size_t slashPos = it->rfind('/');
706  if (slashPos != string::npos)
707  it->erase(0, slashPos + 1);
708  }
709  // build index lookup-table [reflectivity][wave index] to index in normalization integral
710  waveToIntMapType indexLookUp(extents[2][_nmbWavesReflMax]);
711  for (unsigned int iRefl = 0; iRefl < 2; ++iRefl)
712  for (unsigned int iWave = 0; iWave < _nmbWavesRefl[iRefl]; ++iWave) {
713  if (find(intWaveNames.begin(), intWaveNames.end(), _waveNames[iRefl][iWave])
714  == intWaveNames.end()) {
715  printErr << "wave " << _waveNames[iRefl][iWave] << " is not in integral. aborting." << endl;
716  throw;
717  }
718  indexLookUp[iRefl][iWave] = integral.index(_waveNames[iRefl][iWave]);
719  if (_debug)
720  printDebug << " mapping wave [" << sign((int)iRefl * 2 - 1) << ", "
721  << setw(3) << iWave << "] '" << _waveNames[iRefl][iWave] << "' "
722  << "to index " << setw(3) << indexLookUp[iRefl][iWave] << " in integral." << endl;
723  }
724  // create reordered matrix
725  reorderedMatrix.resize(extents[2][_nmbWavesReflMax][2][_nmbWavesReflMax]);
726  for (unsigned int iRefl = 0; iRefl < 2; ++iRefl)
727  for (unsigned int iWave = 0; iWave < _nmbWavesRefl[iRefl]; ++iWave)
728  for (unsigned int jRefl = 0; jRefl < 2; ++jRefl)
729  for (unsigned int jWave = 0; jWave < _nmbWavesRefl[jRefl]; ++jWave) {
730  const complex<double> val = intMatrix.element(indexLookUp[iRefl][iWave],
731  indexLookUp[jRefl][jWave]);
732  reorderedMatrix[iRefl][iWave][jRefl][jWave] = complexT(val.real(), val.imag());
733  }
734 }
735 
736 
737 // returns integral matrix reordered according to _waveNames array
738 template<typename complexT>
739 void
740 TPWALikelihood<complexT>::reorderIntegralMatrix(const ampIntegralMatrix& integral,
741  normMatrixArrayType& reorderedMatrix) const
742 {
743  // create reordered matrix
744  reorderedMatrix.resize(extents[2][_nmbWavesReflMax][2][_nmbWavesReflMax]);
745  for (unsigned int iRefl = 0; iRefl < 2; ++iRefl)
746  for (unsigned int iWave = 0; iWave < _nmbWavesRefl[iRefl]; ++iWave)
747  for (unsigned int jRefl = 0; jRefl < 2; ++jRefl)
748  for (unsigned int jWave = 0; jWave < _nmbWavesRefl[jRefl]; ++jWave) {
749  const complex<double> val = integral.element(_waveNames[iRefl][iWave],
750  _waveNames[jRefl][jWave]);
751  reorderedMatrix[iRefl][iWave][jRefl][jWave] = complexT(val.real(), val.imag());
752  }
753 }
754 
755 
756 template<typename complexT>
757 void
758 TPWALikelihood<complexT>::readIntegrals
759 (const string& normIntFileName, // name of file with normalization integrals
760  const string& accIntFileName, // name of file with acceptance integrals
761  const string& integralTKeyName) // name of TKey which stores integral in .root file
762 {
763  printInfo << "loading normalization integral from '" << normIntFileName << "'" << endl;
764  const string normIntFileExt = extensionFromPath(normIntFileName);
765 #ifdef USE_STD_COMPLEX_TREE_LEAFS
766  if (normIntFileExt == "root") {
767  TFile* intFile = TFile::Open(normIntFileName.c_str(), "READ");
768  if (not intFile or intFile->IsZombie()) {
769  printErr << "could open normalization integral file '" << normIntFileName << "'. "
770  << "aborting." << endl;
771  throw;
772  }
773  ampIntegralMatrix* integral = 0;
774  intFile->GetObject(integralTKeyName.c_str(), integral);
775  if (not integral) {
776  printErr << "cannot find integral object in TKey '" << integralTKeyName << "' in file "
777  << "'" << normIntFileName << "'. aborting." << endl;
778  throw;
779  }
780  reorderIntegralMatrix(*integral, _normMatrix);
781  intFile->Close();
782  } else
783 #endif // USE_STD_COMPLEX_TREE_LEAFS
784  if (normIntFileExt == "int") {
785  ifstream intFile(normIntFileName.c_str());
786  if (not intFile) {
787  printErr << "cannot open file '" << normIntFileName << "'. aborting." << endl;
788  throw;
789  }
790  integral integral;
791  // !!! integral.scan() performs no error checks!
792  integral.scan(intFile);
793  intFile.close();
794  reorderIntegralMatrix(integral, _normMatrix);
795  } else {
796  printErr << "unknown file type '" << normIntFileName << "'. "
797  << "only .int and .root files are supported. aborting." << endl;
798  throw;
799  }
800 
801  printInfo << "loading acceptance integral from '" << accIntFileName << "'" << endl;
802  const string accIntFileExt = extensionFromPath(accIntFileName);
803 #ifdef USE_STD_COMPLEX_TREE_LEAFS
804  if (accIntFileExt == "root") {
805  TFile* intFile = TFile::Open(accIntFileName.c_str(), "READ");
806  if (not intFile or intFile->IsZombie()) {
807  printErr << "could open normalization integral file '" << accIntFileName << "'. "
808  << "aborting." << endl;
809  throw;
810  }
811  ampIntegralMatrix* integral = 0;
812  intFile->GetObject(integralTKeyName.c_str(), integral);
813  if (not integral) {
814  printErr << "cannot find integral object in TKey '" << integralTKeyName << "' in file "
815  << "'" << accIntFileName << "'. aborting." << endl;
816  throw;
817  }
818  if (_numbAccEvents != 0) {
819  _totAcc = ((double)integral->nmbEvents()) / (double)_numbAccEvents;
820  printInfo << "total acceptance in this bin: " << _totAcc << endl;
821  integral->setNmbEvents(_numbAccEvents);
822  } else
823  _totAcc = 1;
824  reorderIntegralMatrix(*integral, _accMatrix);
825  intFile->Close();
826  } else
827 #endif // USE_STD_COMPLEX_TREE_LEAFS
828  if (accIntFileExt == "int") {
829  ifstream intFile(accIntFileName.c_str());
830  if (not intFile) {
831  printErr << "cannot open file '" << accIntFileName << "'. aborting." << endl;
832  throw;
833  }
834  integral integral;
835  // !!! integral.scan() performs no error checks!
836  integral.scan(intFile);
837  intFile.close();
838  if (_numbAccEvents != 0) {
839  _totAcc = ((double)integral.nevents()) / (double)_numbAccEvents;
840  printInfo << "total acceptance in this bin: " << _totAcc << endl;
841  integral.events(_numbAccEvents);
842  } else
843  _totAcc = 1;
844  reorderIntegralMatrix(integral, _accMatrix);
845  } else {
846  printErr << "unknown file type '" << accIntFileName << "'. "
847  << "only .int and .root files are supported. exiting." << endl;
848  throw;
849  }
850 }
851 
852 
853 template<typename complexT>
854 void
855 TPWALikelihood<complexT>::readDecayAmplitudes(const string& ampDirName,
856  const bool useRootAmps,
857  const string& ampLeafName)
858 {
859  // check that normalization integrals are loaded
860  if (_normMatrix.num_elements() == 0) {
861  printErr << "normalization integrals have to be loaded before loading the amplitudes. "
862  << "aborting." << endl;
863  throw;
864  }
865  clear();
866 
867  printInfo << "loading amplitude data from " << ((useRootAmps) ? ".root" : ".amp")
868  << " files" << endl;
869  // loop over amplitudes and read in data
870  unsigned int nmbEvents = 0;
871  bool firstWave = true;
872  for (unsigned int iRefl = 0; iRefl < 2; ++iRefl)
873  for (unsigned int iWave = 0; iWave < _nmbWavesRefl[iRefl]; ++iWave) {
874  // get normalization
875  const complexT normInt = _normMatrix[iRefl][iWave][iRefl][iWave];
876  vector<complexT> amps;
877  if (!firstWave) // number of events is known except for first wave that is read in
878  amps.reserve(nmbEvents);
879  // read decay amplitudes
880  string ampFilePath = ampDirName + "/" + _waveNames[iRefl][iWave];
881 #ifdef USE_STD_COMPLEX_TREE_LEAFS
882  if (useRootAmps) {
883  ampFilePath = changeFileExtension(ampFilePath, ".root");
884  printInfo << "loading amplitude data from '" << ampFilePath << "'" << endl;
885  // open amplitude file
886  TFile* ampFile = TFile::Open(ampFilePath.c_str(), "READ");
887  if (not ampFile or ampFile->IsZombie()) {
888  printWarn << "cannot open amplitude file '" << ampFilePath << "'. skipping." << endl;
889  continue;
890  }
891  // find amplitude tree
892  TTree* ampTree = 0;
893  const string ampTreeName = changeFileExtension(_waveNames[iRefl][iWave], ".amp");
894  ampFile->GetObject(ampTreeName.c_str(), ampTree);
895  if (not ampTree) {
896  printWarn << "cannot find tree '" << ampTreeName << "' in file "
897  << "'" << ampFilePath << "'. skipping." << endl;
898  continue;
899  }
900  // connect tree leaf
901  amplitudeTreeLeaf* ampTreeLeaf = 0;
902  ampTree->SetBranchAddress(ampLeafName.c_str(), &ampTreeLeaf);
903  for (long int eventIndex = 0; eventIndex < ampTree->GetEntriesFast(); ++eventIndex) {
904  ampTree->GetEntry(eventIndex);
905  if (!ampTreeLeaf) {
906  printWarn << "null pointer to amplitude leaf for event " << eventIndex << ". "
907  << "skipping." << endl;
908  continue;
909  }
910  assert(ampTreeLeaf->nmbIncohSubAmps() == 1);
911  complexT amp(ampTreeLeaf->incohSubAmp(0).real(), ampTreeLeaf->incohSubAmp(0).imag());
912  if (_useNormalizedAmps) // normalize data, if option is switched on
913  amp /= sqrt(normInt.real()); // rescale decay amplitude
914  amps.push_back(amp);
915  }
916  } else
917 #endif // USE_STD_COMPLEX_TREE_LEAFS
918  {
919  printInfo << "loading amplitude data from '" << ampFilePath << "'" << endl;
920  ifstream ampFile(ampFilePath.c_str());
921  if (not ampFile) {
922  printErr << "cannot open amplitude file '" << ampFilePath << "'. aborting." << endl;
923  throw;
924  }
925  complexT amp;
926  while (ampFile.read((char*)&amp, sizeof(complexT))) {
927  if (_useNormalizedAmps) { // normalize data, if option is switched on
928  value_type mynorm=normInt.real();
929  if(mynorm==0)mynorm=1;
930  amp /= sqrt(mynorm); // rescale decay amplitude
931  }
932  amps.push_back(amp);
933  }
934  }
935  if (firstWave)
936  nmbEvents = _nmbEvents = amps.size();
937  else {
938  nmbEvents = amps.size();
939  if (nmbEvents != _nmbEvents)
940  printWarn << "size mismatch in amplitude files: this file contains " << nmbEvents
941  << " events, previous file had " << _nmbEvents << " events." << endl;
942  nmbEvents = _nmbEvents;
943  }
944 
945  // copy decay amplitudes into array that is indexed [event index][reflectivity][wave index]
946  // this index scheme ensures a more linear memory access pattern in the likelihood function
947  _decayAmps.resize(extents[_nmbEvents][2][_nmbWavesReflMax]);
948  for (unsigned int iEvt = 0; iEvt < _nmbEvents; ++iEvt)
949  _decayAmps[iEvt][iRefl][iWave] = amps[iEvt];
950  if (_debug)
951  printDebug << "read " << _nmbEvents << " events from file "
952  << "'" << _waveNames[iRefl][iWave] << "'" << endl;
953  }
954  printInfo << "loaded decay amplitudes for " << _nmbEvents << " events into memory" << endl;
955 
956  // save phase space integrals
957  _phaseSpaceIntegral.resize(extents[2][_nmbWavesReflMax]);
958  for (unsigned int iRefl = 0; iRefl < 2; ++iRefl)
959  for (unsigned int iWave = 0; iWave < _nmbWavesRefl[iRefl]; ++iWave)
960  _phaseSpaceIntegral[iRefl][iWave] = sqrt(_normMatrix[iRefl][iWave][iRefl][iWave].real());
961 
962  // rescale integrals, if necessary
963  if (_useNormalizedAmps) {
964  // matrices _normMatrix and _accMatrix are already normalized to number of Monte Carlo events
965  printInfo << "rescaling integrals" << endl;
966  // rescale normalization and acceptance integrals
967  for (unsigned int iRefl = 0; iRefl < 2; ++iRefl)
968  for (unsigned int iWave = 0; iWave < _nmbWavesRefl[iRefl]; ++iWave) {
969  value_type norm_i = sqrt(_normMatrix[iRefl][iWave][iRefl][iWave].real());
970  if(norm_i==0)norm_i=1;
971  for (unsigned int jRefl = 0; jRefl < 2; ++jRefl)
972  for (unsigned int jWave = 0; jWave < _nmbWavesRefl[jRefl]; ++jWave) {
973  value_type norm_j = sqrt(_normMatrix[jRefl][jWave][jRefl][jWave].real());
974  // protect against empty amplitudes (which will not contribute anyway)
975  if(norm_j==0)norm_j=1;
976  if ((iRefl != jRefl) or (iWave != jWave))
977  // set diagonal terms later so that norm_i,j stay unaffected
978  _normMatrix[iRefl][iWave][jRefl][jWave] /= norm_i * norm_j;
979  _accMatrix[iRefl][iWave][jRefl][jWave] /= norm_i * norm_j;
980  }
981  }
982  // set diagonal elements of normalization matrix
983  for (unsigned int iRefl = 0; iRefl < 2; ++iRefl)
984  for (unsigned int iWave = 0; iWave < _nmbWavesRefl[iRefl]; ++iWave)
985  _normMatrix[iRefl][iWave][iRefl][iWave] = 1; // diagonal term
986  if (_debug) {
987  printDebug << "normalized integral matrices" << endl;
988  for (unsigned int iRefl = 0; iRefl < 2; ++iRefl)
989  for (unsigned int iWave = 0; iWave < _nmbWavesRefl[iRefl]; ++iWave)
990  for (unsigned int jRefl = 0; jRefl < 2; ++jRefl)
991  for (unsigned int jWave = 0; jWave < _nmbWavesRefl[jRefl]; ++jWave) {
992  cout << " normalization matrix [" << sign((int)iRefl * 2 - 1) << ", "
993  << setw(3) << iWave << ", " << sign((int)jRefl * 2 - 1) << ", "
994  << setw(3) << jWave << "] = "
995  << "(" << maxPrecisionAlign(_normMatrix[iRefl][iWave][jRefl][jWave].real())
996  << ", " << maxPrecisionAlign(_normMatrix[iRefl][iWave][jRefl][jWave].imag())
997  << "), acceptance matrix [" << setw(3) << iWave << ", "
998  << setw(3) << jWave << "] = "
999  << "(" << maxPrecisionAlign(_accMatrix[iRefl][iWave][jRefl][jWave].real())
1000  << ", " << maxPrecisionAlign(_accMatrix[iRefl][iWave][jRefl][jWave].imag()) << ")"
1001  << endl;
1002  }
1003  }
1004  } // _useNormalizedAmps
1005 }
1006 
1007 
1008 template<typename complexT>
1009 void
1010 TPWALikelihood<complexT>::getIntegralMatrices(TCMatrix& normMatrix,
1011  TCMatrix& accMatrix,
1012  vector<double>& phaseSpaceIntegral) const
1013 {
1014  phaseSpaceIntegral.clear();
1015  phaseSpaceIntegral.resize(_nmbWaves + 1, 0);
1016  unsigned int iIndex = 0;
1017  for (unsigned int iRefl = 0; iRefl < 2; ++iRefl) {
1018  for (unsigned int iWave = 0; iWave < _nmbWavesRefl[iRefl]; ++iWave) {
1019  phaseSpaceIntegral[iIndex] = _phaseSpaceIntegral[iRefl][iWave];
1020  unsigned int jIndex = 0;
1021  for (unsigned int jRefl = 0; jRefl < 2; ++jRefl) {
1022  for (unsigned int jWave = 0; jWave < _nmbWavesRefl[jRefl]; ++jWave) {
1023  const complexT normVal = _normMatrix[iRefl][iWave][jRefl][jWave];
1024  const complexT accVal = _accMatrix [iRefl][iWave][jRefl][jWave];
1025  normMatrix.set(iIndex, jIndex, complex<double>(normVal.real(), normVal.imag()));
1026  accMatrix.set (iIndex, jIndex, complex<double>(accVal.real(), accVal.imag() ));
1027  ++jIndex;
1028  }
1029  }
1030  ++iIndex;
1031  }
1032  }
1033  // add flat
1034  normMatrix.set(_nmbWaves, _nmbWaves, complexT(1,0));
1035  accMatrix.set (_nmbWaves, _nmbWaves, complexT(1,0));
1036  phaseSpaceIntegral[_nmbWaves] = 1;
1037 }
1038 
1039 
1040 // builds complex numbers from parameters
1041 // maps real and imaginary part of amplitudes to error matrix
1042 // for both rank restrictions are taken into account
1043 template<typename complexT>
1044 void
1045 TPWALikelihood<complexT>::buildProdAmpArrays(const double* inPar,
1046  vector<complex<double> >& prodAmps,
1047  vector<pair<int, int> >& parIndices,
1048  vector<string>& prodAmpNames,
1049  const bool withFlat) const
1050 {
1051  prodAmps.clear();
1052  parIndices.clear();
1053  prodAmpNames.clear();
1054  unsigned int parIndex = 0;
1055  for (unsigned int iRank = 0; iRank < _rank; ++iRank)
1056  for (unsigned int iRefl = 0; iRefl < 2; ++iRefl)
1057  for (unsigned int iWave = 0; iWave < _nmbWavesRefl[iRefl]; ++iWave) {
1058  double re, im;
1059  if (iWave < iRank) // zero production amplitude
1060  continue;
1061  else if (iWave == iRank) { // real production amplitude
1062  parIndices.push_back(make_pair(parIndex, -1));
1063  re = inPar[parIndex];
1064  im = 0;
1065  ++parIndex;
1066  } else { // complex production amplitude
1067  parIndices.push_back(make_pair(parIndex, parIndex + 1));
1068  re = inPar[parIndex];
1069  im = inPar[parIndex + 1];
1070  parIndex += 2;
1071  }
1072  prodAmps.push_back(complex<double>(re, im));
1073  stringstream prodAmpName;
1074  prodAmpName << "V" << iRank << "_" << _waveNames[iRefl][iWave];
1075  prodAmpNames.push_back(prodAmpName.str());
1076  }
1077  if (withFlat) {
1078  prodAmps.push_back(complex<double>(inPar[parIndex], 0));
1079  parIndices.push_back(make_pair(parIndex, -1));
1080  prodAmpNames.push_back("V_flat");
1081  }
1082 }
1083 
1084 
1085 template<typename complexT>
1086 void
1087 TPWALikelihood<complexT>::clear()
1088 {
1089  _decayAmps.resize(extents[0][0][0]);
1090 }
1091 
1092 
1093 // depends on naming convention for waves!!!
1094 // VR_IGJPCMEIso....
1095 template<typename complexT>
1096 int
1097 TPWALikelihood<complexT>::getReflectivity(const TString& waveName)
1098 {
1099  int refl = 0;
1100  unsigned int reflIndex = 6; // position of reflectivity in wave
1101  // check whether it is parameter or wave name
1102  if (waveName[0] == 'V')
1103  reflIndex = 9;
1104  if (waveName[reflIndex] == '-')
1105  refl= -1;
1106  else if (waveName[reflIndex] == '+')
1107  refl= +1;
1108  else {
1109  printErr << "cannot parse parameter/wave name '" << waveName << "'. "
1110  << "cannot not determine reflectivity. aborting." << endl;
1111  throw;
1112  }
1113  if (_debug)
1114  printDebug << "extracted reflectivity = " << refl << " from parameter name "
1115  << "'" << waveName << "' (char position " << reflIndex << ")" << endl;
1116  return refl;
1117 }
1118 
1119 
1120 // copy values from array that corresponds to the function parameters
1121 // to structure that corresponds to the complex production amplitudes
1122 // taking into account rank restrictions
1123 template<typename complexT>
1124 void
1125 TPWALikelihood<complexT>::copyFromParArray
1126 (const double* inPar, // input parameter array
1127  ampsArrayType& outVal, // array of complex output values [rank][reflectivity][wave index]
1128  value_type& outFlatVal) const // output value corresponding to flat wave
1129 {
1130  // group parameters by rank and wave index only; keeps the order defined in wave list
1131  outVal.resize(extents[_rank][2][max(_nmbWavesRefl[0], _nmbWavesRefl[1])]);
1132  unsigned int parIndex = 0;
1133  for (unsigned int iRank = 0; iRank < _rank; ++iRank)
1134  for (unsigned int iRefl = 0; iRefl < 2; ++iRefl)
1135  for (unsigned int iWave = 0; iWave < _nmbWavesRefl[iRefl]; ++iWave) {
1136  value_type re, im;
1137  if (iWave < iRank) // production amplitude is zero
1138  re = im = 0;
1139  else if (iWave == iRank) { // production amplitude is real
1140  re = inPar[parIndex];
1141  im = 0;
1142  ++parIndex;
1143  } else { // production amplitude is complex
1144  re = inPar[parIndex];
1145  im = inPar[parIndex + 1];
1146  parIndex += 2;
1147  }
1148  outVal[iRank][iRefl][iWave] = complexT(re, im);
1149  }
1150  outFlatVal = inPar[parIndex];
1151 }
1152 
1153 
1154 // copy values from structure that corresponds to complex
1155 // production amplitudes to array that corresponds to function
1156 // parameters taking into account rank restrictions
1157 template<typename complexT>
1158 void
1159 TPWALikelihood<complexT>::copyToParArray
1160 (const ampsArrayType& inVal, // values corresponding to production amplitudes
1161  const value_type inFlatVal, // value corresponding to flat wave
1162  double* outPar) const // output parameter array
1163 {
1164  for (unsigned int iRank = 0; iRank < _rank; ++iRank)
1165  for (unsigned int iRefl = 0; iRefl < 2; ++iRefl)
1166  for (unsigned int iWave = 0; iWave < _nmbWavesRefl[iRefl]; ++iWave) {
1167  tuple<int, int> parIndices = _prodAmpToFuncParMap[iRank][iRefl][iWave];
1168  if (get<0>(parIndices) >= 0) // real part
1169  outPar[get<0>(parIndices)] = inVal[iRank][iRefl][iWave].real();
1170  if (get<1>(parIndices) >= 0) // imaginary part
1171  outPar[get<1>(parIndices)] = inVal[iRank][iRefl][iWave].imag();
1172  }
1173  outPar[_nmbPars - 1] = inFlatVal;
1174 }
1175 
1176 
1177 template<typename complexT>
1178 ostream&
1179 TPWALikelihood<complexT>::print(ostream& out) const
1180 {
1181  out << "TPWALikelihood parameters:" << endl
1182  << "number of events ........................ " << _nmbEvents << endl
1183  << "rank .................................... " << _rank << endl
1184  << "number of waves ......................... " << _nmbWaves << endl
1185  << "number of positive reflectivity waves ... " << _nmbWavesRefl[1] << endl
1186  << "number of negative reflectivity waves ... " << _nmbWavesRefl[0] << endl
1187  << "number of function parameters ........... " << _nmbPars << endl
1188  << "print debug messages .................... " << _debug << endl
1189 #ifdef USE_CUDA
1190  << "use CUDA kernels ........................ " << _cudaEnabled << endl
1191 #endif
1192  << "use normalized amplitudes ............... " << _useNormalizedAmps << endl
1193  << "list of waves: " << endl;
1194  for (unsigned int iRefl = 0; iRefl < 2; ++iRefl)
1195  for (unsigned int iWave = 0; iWave < _nmbWavesRefl[iRefl]; ++iWave)
1196  out << " [" << setw(2) << sign((int)iRefl * 2 - 1) << " " << setw(3) << iWave << "] "
1197  << _waveNames[iRefl][iWave] << " threshold = "
1198  << _waveThresholds[iRefl][iWave] << " MeV/c^2" << endl;
1199  out << "list of function parameters: " << endl;
1200  for (unsigned int iPar = 0; iPar < _nmbPars; ++iPar)
1201  out << " [" << setw(3) << iPar << "] " << _parNames[iPar] << " "
1202  << "threshold = " << _parThresholds[iPar] << " MeV/c^2" << endl;
1203  return out;
1204 }
1205 
1206 
1207 template<typename complexT>
1208 ostream&
1209 TPWALikelihood<complexT>::printFuncInfo(ostream& out) const
1210 {
1211  const string funcNames[NMB_FUNCTIONCALLENUM] = {"FdF", "Gradient", "DoEval", "DoDerivative"};
1212  for (unsigned int i = 0; i < NMB_FUNCTIONCALLENUM; ++i)
1213  if (_funcCallInfo[i].nmbCalls > 0)
1214  out << " " << _funcCallInfo[i].nmbCalls
1215  << " calls to TPWALikelihood<complexT>::" << funcNames[i] << "()" << endl
1216  << " time spent in TPWALikelihood<complexT>::" << funcNames[i] << "(): " << endl
1217  << " total time ................. " << sum(_funcCallInfo[i].totalTime) << " sec" << endl
1218  << " return value calculation ... " << sum(_funcCallInfo[i].funcTime ) << " sec" << endl
1219  << " normalization .............. " << sum(_funcCallInfo[i].normTime ) << " sec" << endl;
1220  return out;
1221 }
1222 
1223 
1224 template<typename complexT>
1225 vector<unsigned int>
1226 TPWALikelihood<complexT>::orderedParIndices() const
1227 {
1228  vector<unsigned int> orderedIndices;
1229  for (unsigned int iRank = 0; iRank < _rank; ++iRank)
1230  for (unsigned int waveIndex = 0; waveIndex < _nmbWaves; ++waveIndex) {
1231  unsigned int iRefl, iWave;
1232  for (iRefl = 0; iRefl < 2; ++iRefl)
1233  for (iWave = 0; iWave < _nmbWavesRefl[iRefl]; ++iWave)
1234  if (_waveToWaveIndex[iRefl][iWave] == waveIndex)
1235  goto found;
1236  printWarn << "indices are inconsistent. cannot find wave with index " << waveIndex
1237  << " in wave list" << endl;
1238  continue;
1239 
1240  found:
1241 
1242  tuple<int, int> parIndices = _prodAmpToFuncParMap[iRank][iRefl][iWave];
1243  int parIndex = get<0>(parIndices);
1244  if (parIndex >= 0)
1245  orderedIndices.push_back(parIndex);
1246  parIndex = get<1>(parIndices);
1247  if (parIndex >= 0)
1248  orderedIndices.push_back(parIndex);
1249  }
1250  orderedIndices.push_back(_nmbPars - 1); // flat wave
1251  if (orderedIndices.size() != _nmbPars)
1252  printWarn << "ordered list of parameter indices has inconsistent size "
1253  << "(" << orderedIndices.size() << " vs. " << _nmbPars << ")" << endl;
1254  return orderedIndices;
1255 }
1256 
1257 
1258 template<typename complexT>
1259 vector<string>
1260 TPWALikelihood<complexT>::waveNames() const
1261 {
1262  vector<string> names(_nmbWaves, "");
1263  for (unsigned int iRefl = 0; iRefl < 2; ++iRefl)
1264  for (unsigned int iWave = 0; iWave < _nmbWavesRefl[iRefl]; ++iWave) {
1265  const unsigned int index = _waveToWaveIndex[iRefl][iWave];
1266  names[index] = _waveNames[iRefl][iWave];
1267  }
1268  return names;
1269 }
1270 
1271 
1272 template<typename complexT>
1273 void
1274 TPWALikelihood<complexT>::resetFuncCallInfo() const
1275 {
1276  for (unsigned int i = 0; i < NMB_FUNCTIONCALLENUM; ++i) {
1277  _funcCallInfo[i].nmbCalls = 0;
1278  // boost accumulators do not have a reset function
1279  _funcCallInfo[i].funcTime = typename functionCallInfo::timeAccType();
1280  _funcCallInfo[i].normTime = typename functionCallInfo::timeAccType();
1281  _funcCallInfo[i].totalTime = typename functionCallInfo::timeAccType();
1282  }
1283 }
1284 
1285 
1286 // explicit specializations
1287 template class TPWALikelihood<complex<float > >;
1288 template class TPWALikelihood<complex<double> >;
Generated by   doxygen 1.8.1.2