ROOTPWA
TPWALikelihood.h
Go to the documentation of this file.
1 
2 //
3 // Copyright 2009 Sebastian Neubert
4 //
5 // This file is part of rootpwa
6 //
7 // rootpwa is free software: you can redistribute it and/or modify
8 // it under the terms of the GNU General Public License as published by
9 // the Free Software Foundation, either version 3 of the License, or
10 // (at your option) any later version.
11 //
12 // rootpwa is distributed in the hope that it will be useful,
13 // but WITHOUT ANY WARRANTY; without even the implied warranty of
14 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 // GNU General Public License for more details.
16 //
17 // You should have received a copy of the GNU General Public License
18 // along with rootpwa. If not, see <http://www.gnu.org/licenses/>.
19 //
21 //-----------------------------------------------------------
22 // File and Version Information:
23 // $Id$
24 //
25 // Description:
26 // Likelihood function Object to use with ROOT minimizers
27 //
28 //
29 // Author List:
30 // Sebastian Neubert TUM (original author)
31 //
32 //
33 //-----------------------------------------------------------
34 
35 
36 #ifndef TPWALIKELIHOOD_H
37 #define TPWALIKELIHOOD_H
38 
39 
40 #include <vector>
41 #include <string>
42 #include <iostream>
43 
44 #define BOOST_DISABLE_ASSERTS
45 #include "boost/multi_array.hpp"
46 #include "boost/tuple/tuple.hpp"
47 
48 #include "Math/IFunction.h"
49 #include "TFile.h"
50 #include "TH1.h"
51 
52 // PWA2000 classes
53 #include "integral.h"
54 #include "matrix.h"
55 
56 #include "sumAccumulators.hpp"
57 #include "ampIntegralMatrix.h"
58 
59 
60 class TString;
61 class TCMatrix;
62 
63 
64 template<typename complexT> // type of internal variables used for intermediate results
65 class TPWALikelihood : public ROOT::Math::IGradientFunctionMultiDim {
66 
67 public:
68 
69  typedef typename complexT::value_type value_type;
70 
71  // define array types
72  typedef boost::multi_array<std::string, 2> waveNameArrayType; // array for wave names
73  typedef boost::multi_array<double, 2> waveThrArrayType; // array for wave thresholds
74  typedef boost::multi_array<unsigned int, 2> waveToIntMapType; // array for mapping of waves to integral indices
75  typedef boost::multi_array<unsigned int, 2> waveToListMapType; // array for mapping of waves to position in wave list
76  typedef boost::multi_array<boost::tuple<int, int>, 3> ampToParMapType; // array for mapping of amplitudes to parameters
77  typedef boost::multi_array<complexT, 3> ampsArrayType; // array for production and decay amplitudes
78  typedef boost::multi_array<complexT, 4> normMatrixArrayType; // array for normalization matrices
79  typedef boost::multi_array<value_type, 2> phaseSpaceIntType; // array for phase space integrals
80 
81 
82 public:
83 
84  // enum for function call counters
85  enum functionCallEnum {
86  FDF = 0,
87  GRADIENT = 1,
88  DOEVAL = 2,
89  DODERIVATIVE = 3,
90  NMB_FUNCTIONCALLENUM = 4
91  };
92 
93  struct functionCallInfo {
94  typedef boost::accumulators::accumulator_set
95  <double, boost::accumulators::stats
96  <boost::accumulators::tag::sum(boost::accumulators::compensated)> > timeAccType;
97  unsigned int nmbCalls; // number of times function was called
98  timeAccType funcTime; // time needed to calculate function value(s) (w/o normalization)
99  timeAccType normTime; // time needed to normalize function value(s)
100  timeAccType totalTime; // total execution time of function
101  };
102 
103  TPWALikelihood();
104  ~TPWALikelihood();
105 
106  // overload public IGradientFunctionMultiDim member functions:
108  virtual TPWALikelihood* Clone() const { return new TPWALikelihood(*this); }
110  virtual unsigned int NDim() const { return nmbPars(); }
112  virtual void FdF(const double* par,
113  double& funcVal,
114  double* gradient) const;
116  virtual void Gradient(const double* par,
117  double* gradient) const;
118 
119  // overload private IGradientFunctionMultiDim member functions
120  virtual double DoEval (const double* par) const;
121  virtual double DoDerivative(const double* par,
122  unsigned int derivativeIndex) const;
123 
124  unsigned int nmbEvents () const { return _nmbEvents; }
125  unsigned int rank () const { return _rank; }
126  unsigned int nmbWaves (const int reflectivity = 0) const;
127  unsigned int nmbPars () const { return _nmbPars; }
128  // std::string waveName (const unsigned int waveIndex) const { return _waveNames[waveIndex]; } ///< returns name of wave at waveIndex
129  std::vector<std::string> waveNames () const;
130  std::string parName (const unsigned int parIndex) const { return _parNames[parIndex]; }
131  double parThreshold(const unsigned int parIndex) const { return _parThresholds[parIndex]; }
132 
133  double dLcache(const unsigned int i) const { return _derivCache[i]; }
134  unsigned int ncalls(const functionCallEnum func = FDF) const
135  { return _funcCallInfo[func].nmbCalls; }
136  double Ltime(const functionCallEnum func = FDF) const
137  { return boost::accumulators::sum(_funcCallInfo[func].funcTime); }
138  double Ntime(const functionCallEnum func = FDF) const
139  { return boost::accumulators::sum(_funcCallInfo[func].normTime); }
140  //const integral& normInt() const { return _normInt; }
141 
142  // modifiers
143  void enableCuda (const bool enableCuda = true);
144  bool cudaEnabled () const;
145  void useNormalizedAmps(const bool useNorm = true) { _useNormalizedAmps = useNorm; }
146  static void setQuiet (const bool flag = true) { _debug = !flag; }
147 
148  // operations
149  void init(const unsigned int rank,
150  const std::string& waveListFileName,
151  const std::string& normIntFileName,
152  const std::string& accIntFileName,
153  const std::string& ampDirName = ".",
154  const unsigned int numbAccEvents = 0,
155  const bool useRootAmps = false);
156 
157  void getIntegralMatrices(TCMatrix& normMatrix,
158  TCMatrix& accMatrix,
159  std::vector<double>& phaseSpaceIntegral) const;
160 
161  // note: amplitudes which do not exist in higher ranks are NOT built!
162  void buildProdAmpArrays(const double* inPar,
163  std::vector<std::complex<double> >& prodAmps,
164  std::vector<std::pair<int,int> >& parIndices,
165  std::vector<std::string>& prodAmpNames,
166  const bool withFlat = false) const;
167 
168  std::ostream& print(std::ostream& out = std::cout) const;
169  std::ostream& printFuncInfo(std::ostream& out = std::cout) const;
170  friend std::ostream& operator << (std::ostream& out,
171  const TPWALikelihood& func) { return func.print(out); }
172 
173  std::vector<unsigned int> orderedParIndices() const; // helper function for backwards-compatibility
174 
175 
176 private:
177 
178  // helper functions
179  void readWaveList (const std::string& waveListFileName);
180  void buildParDataStruct (const unsigned int rank);
181  void readIntegrals (const std::string& normIntFileName,
182  const std::string& accIntFileName,
183  const std::string& integralTKeyName = "integral");
184 
185  void readDecayAmplitudes(const std::string& ampDirName = ".",
186  const bool useRootAmps = false,
187  const std::string& ampLeafName = "amplitude");
188 
189 
190  void clear();
191  static int getReflectivity(const TString& waveName);
192 
193  void reorderIntegralMatrix(integral& integral,
194  normMatrixArrayType& reorderedMatrix) const;
195  void reorderIntegralMatrix(const rpwa::ampIntegralMatrix& integral,
196  normMatrixArrayType& reorderedMatrix) const;
197 
198 public:
199 
200  void copyFromParArray(const double* inPar, // input parameter array
201  ampsArrayType& outVal, // output values organized as 3D array of complex numbers with [rank][reflectivity][wave index]
202  value_type& outFlatVal) const; // output value corresponding to flat wave
203  void copyToParArray(const ampsArrayType& inVal, // values corresponding to production amplitudes [rank][reflectivity][wave index]
204  const value_type inFlatVal, // value corresponding to flat wave
205  double* outPar) const; // output parameter array
206 
207 private:
208 
209  void resetFuncCallInfo() const;
210 
211  unsigned int _nmbEvents; // number of events
212  unsigned int _rank; // rank of spin density matrix
213  unsigned int _nmbWaves; // number of waves
214  unsigned int _nmbWavesRefl[2]; // number of negative (= 0) and positive (= 1) reflectivity waves
215  unsigned int _nmbWavesReflMax; // maximum of number of negative and positive reflectivity waves
216  unsigned int _nmbPars; // number of function parameters
217 
218 #ifdef USE_CUDA
219  bool _cudaEnabled; // if true CUDA kernels are used for some calculations
220 #endif
221  bool _useNormalizedAmps; // if true normalized amplitudes are used
222  static bool _debug; // if true debug messages are printed
223 
224  unsigned int _numbAccEvents; // number of input events used for acceptance integrals (accepted + rejected!)
225  double _totAcc; // total acceptance in this bin
226 
227  waveNameArrayType _waveNames; // wave names [reflectivity][wave index]
228  waveThrArrayType _waveThresholds; // mass thresholds of waves
229  waveToListMapType _waveToWaveIndex; // maps wave to its index in wave list
230  std::vector<std::string> _parNames; // function parameter names
231  std::vector<double> _parThresholds; // mass thresholds of parameters
232  ampToParMapType _prodAmpToFuncParMap; // maps each production amplitude to the indices
233  // of its real and imginary part in the parameter
234  // array; negative indices mean that the parameter
235  // is not existing due to rank restrictions
236 
237  ampsArrayType _decayAmps; // precalculated decay amplitudes [event index][reflectivity][wave index]
238 
239  mutable std::vector<double> _parCache; // parameter cache for derivative calc.
240  mutable std::vector<double> _derivCache; // cache for derivatives
241 
242  // normalization integrals
243  normMatrixArrayType _normMatrix; // normalization matrix w/o acceptance [reflectivity 1][wave index 1][reflectivity 2][wave index 2]
244  normMatrixArrayType _accMatrix; // normalization matrix with acceptance [reflectivity 1][wave index 1][reflectivity 2][wave index 2]
245  phaseSpaceIntType _phaseSpaceIntegral; // phase space integrals
246 
247  mutable functionCallInfo _funcCallInfo[NMB_FUNCTIONCALLENUM]; // collects function call statistics
248 
249 };
250 
251 
252 #endif // TPWALIKELIHOOD_H
Generated by   doxygen 1.8.1.2