StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
dFitter3d.cxx
1 /***************************************************************************
2  *
3  *
4  * Author: Dominik Flierl, flierl@bnl.gov
5  ***************************************************************************
6  *
7  * Description: part of STAR HBT Framework: StHbtMaker package
8  * a object which uses TMinuit to fit correlationfunctions
9  *
10  *
11  **************************************************************************/
12 
13 #include "dFitter3d.h"
14 #include "Stiostream.h"
15 #include "TMath.h"
16 #include "TH3.h"
17 #include "TH1.h"
18 #include "TMinuit.h"
19 #include <TVector3.h>
20 #include "TString.h"
21 #include "TObject.h"
22 #include "TMath.h"
23 //#ifdef __ROOT__
24 ClassImp(dFitter3d)
25  //#endif
26 
27 // this is the stupid way ROOT wants it: the fcn which is going to be minimzed must
28 // either be static or a global function. here it's global so make sur that we got
29 // the connection to this class via TMinuit->SetObjectFit(this) !
30  void gfcn(Int_t &nParamters, Double_t *gin, Double_t &finalChi2, Double_t *parameterSet, Int_t iflag)
31 {
32  dFitter3d* dummy = dynamic_cast<dFitter3d*>(gMinuit->GetObjectFit()) ;
33  dummy->mfcn(nParamters, gin, finalChi2 , parameterSet, iflag) ;
34 }
35 
36 
37 //____________________________
38 dFitter3d::dFitter3d(TMinuit* gMinuit,TH3D* numerator, TH3D* denominator, TString opt , TString opt2)
39 {
40  // get TMinuit in
41  mMinuit = gMinuit ;
42 
43  // get numerator and denominator in
44  mNumerator = numerator ;
45  mDenominator = denominator ;
46 
47  // set type of correlationfunction
48  if ( opt == "ykp" )
49  {
50  mCorrFctnType = "ykp" ;
51  cout << "Fitting function set to Yano Koonin Podgoretskii. \n" ;
52  }
53  else if ( opt == "bp" )
54  {
55  mCorrFctnType = "bp" ;
56  cout << "Fitting function set to Bertsch Pratt. \n" ;
57  }
58  else
59  {
60  mCorrFctnType = "ykp" ;
61  cout << "Bad option : fitting function set to Yano Koonin Podgoretskii. \n" ;
62  }
63 
64  // set typ of fitting MML or Chi2
65  if ( opt2 == "chi2" )
66  {
67  mFitMethod = "chi2" ;
68  cout << "Using chi2.\n" ;
69  }
70  else if ( opt2 == "mml" )
71  {
72  mFitMethod = "mml" ;
73  cout << "Using mml.\n" ;
74  }
75  else
76  {
77  mFitMethod = "chi2" ;
78  cout << "Bad option : Using chi2 method.\n" ;
79  } ;
80 
81 
82  // set constants
83  mhc = 0.197326960277 ; // GeV * fm
84  mhc2 = 0.038937929230 ; // GeV^2 * fm^2 ;
85 
86  // set defaults for thresholds = 1 means accept any entry
87  mThresholdNumerator = 1.0 ;
88  mThresholdDenominator = 1.0 ;
89 }
90 //____________________________
91 dFitter3d::~dFitter3d()
92 {
93  // nothing to do here
94 }
95 //____________________________
96 void dFitter3d::FillInternalArrays()
97 {
98  // some output
99  cout << " Filling internal arrays with norm factor : " << mNorm
100  << " numerator threshold : " << mThresholdNumerator
101  << " denominator threshold : " << mThresholdDenominator << endl ;
102 
103 
104  // get the maximum number of bins
105  int maximumInternalArraySize = (mNumerator->GetNbinsX()+2) * (mNumerator->GetNbinsY()+2) * (mNumerator->GetNbinsZ()+2) ;
106  if ( maximumInternalArraySize != (mDenominator->GetNbinsX()+2) * (mDenominator->GetNbinsY()+2) * (mDenominator->GetNbinsZ()+2) )
107  std::cerr << "Warning : different histogram sizes in numerator and denominator\n " ;
108 
109  // now create arrays
110  mNumeratorInternalArray = new double[maximumInternalArraySize] ;
111  mDenominatorInternalArray = new double[maximumInternalArraySize] ;
112  mErrorInternalArray = new double[maximumInternalArraySize] ;
113  mVariablePositionArray = new TVector3[maximumInternalArraySize] ;
114 
115  // count number of actually used bins
116  mInternalArraySize = 0 ;
117 
118  // loop over histogram and fill values into vectors
119  for (int index = 0 ; index < maximumInternalArraySize ; index++)
120  {
121  if ( mNumerator->GetBinContent(index) > mThresholdNumerator && mDenominator->GetBinContent(index) > mThresholdDenominator )
122  {
123  // take center of the BIN ( as according position in qlong,qside,qnull )
124  int binX = 0 ;
125  int binY = 0 ;
126  int binZ = 0 ;
127  Bin1ToBin3(mNumerator, index, binX, binY, binZ) ;
128  double qx = mNumerator->GetXaxis()->GetBinCenter(binX) ;
129  double qy = mNumerator->GetYaxis()->GetBinCenter(binY) ;
130  double qz = mNumerator->GetZaxis()->GetBinCenter(binZ) ;
131 
132  if ( qx*qx+qy*qy+qz*qz < mSphereLimit )
133  {
134 
136  // fill the numerator and denominator histogram into vectors
138  double num = mNumerator->GetBinContent(index) ;
139  double denom = mDenominator->GetBinContent(index) ;
140  // fill numerator
141  mNumeratorInternalArray[mInternalArraySize] = num ;
142  // normalize and fill background
143  mDenominatorInternalArray[mInternalArraySize] = denom * mNorm ;
144 
145  // the error is cacluated using Gauss' error propagation : e = num/(denom*norm) * ::sqrt(1/num + 1/(denom))
146  double error = num/(denom*mNorm) * ::sqrt( (1/num) + (1/denom) ) ;
147  mErrorInternalArray[mInternalArraySize] = error ;
148  // the simple case : take center of the BIN ( as according position in qlong,qside,qnull )
149  mVariablePositionArray[mInternalArraySize].SetX(qx) ;
150  mVariablePositionArray[mInternalArraySize].SetY(qy) ;
151  mVariablePositionArray[mInternalArraySize].SetZ(qz) ;
152 
153  // some silly output to see where we are
154  if (0)
155  {
156  cout << (double) num << "\t" << (double) (denom*mNorm) ; //<<index << binX << binY << binZ << endl;
157  cout << "\t" << (double) (num/(denom*mNorm)) ;
158  double xx = mNumerator->GetXaxis()->GetBinCenter(binX) ;
159  double yy = mNumerator->GetYaxis()->GetBinCenter(binY) ;
160  double zz = mNumerator->GetZaxis()->GetBinCenter(binZ) ;
161  cout << "xa " << xx << "\t" << binX << "\t" ;
162  cout << "ya " << yy << "\t" << binY << "\t" ;
163  cout << "za " << zz << "\t" << binZ << endl ;
164  }
165 
166  // count entries
167  mInternalArraySize++ ;
168  }
169  }
170  }
171  // done
172  cout << "Internal array size :" << mInternalArraySize << endl ;
173 }
174 //____________________________
175 void dFitter3d::mfcn(Int_t &nParameters, Double_t *gin, Double_t &finalChi2, Double_t *parameterSet, Int_t iflag)
176 {
177  if ( mFitMethod == "chi2")
178  {
179  fcnChi2(nParameters,gin,finalChi2,parameterSet,iflag) ;
180  }
181  else if ( mFitMethod == "mml")
182  {
183  fcnMml(nParameters,gin,finalChi2,parameterSet,iflag) ;
184  } ;
185 }
186 //____________________________
187 void dFitter3d::fcnChi2(Int_t &nParameters, Double_t *gin, Double_t &finalChi2, Double_t *parameterSet, Int_t iflag)
188 {
190  // calculate chi2 of the given paramset 'parameterSet'
192  double chi2 = 0.0 ;
193 
194  // loop over Internal Array and caluclate chi2
195  for (Int_t index = 0 ; index < mInternalArraySize ; index++)
196  {
197  // calculate the expected value at this point (variablePosition) for this set of paramters
198  double theoreticalValue = mCorrelationFunction( mVariablePositionArray[index] , parameterSet) ;
199  // get the difference between measured and predicted value
200  double delta = (mNumeratorInternalArray[index]/mDenominatorInternalArray[index]-theoreticalValue)/mErrorInternalArray[index];
201  // sum the chi2
202  chi2 += delta*delta;
203 
204 
205  // some silly output
206  if(countMinuitCalls == 100 && index%10000==0)
207  {
208  cout << "ql : " << mVariablePositionArray[index].x() << "\t"
209  << "qp : " << mVariablePositionArray[index].y() << "\t"
210  << "q0 : " << mVariablePositionArray[index].z()<< "\t"
211  << "thval: " << theoreticalValue << "\t"
212  << "mes " << mNumeratorInternalArray[index]/mDenominatorInternalArray[index] << "\t"
213  << "num: " << mNumeratorInternalArray[index] << "\t"
214  << "th-mea" << fabs( mNumeratorInternalArray[index]/mDenominatorInternalArray[index] - theoreticalValue) << "\t"
215  << "err " << mErrorInternalArray[index] << endl ;
216  }
217 
218  }
219  // chi2 is calculated
220  finalChi2 = chi2 ;
221 
222  // count minuit calls
223  countMinuitCalls++ ;
224  //cout << "Minuit call : " << countMinuitCalls << " chi2 : " << finalChi2 << endl ;
225 };
226 //____________________________
227 void dFitter3d::fcnMml(Int_t &nParameters, Double_t *gin, Double_t &finalMLH, Double_t *parameterSet, Int_t iflag)
228 {
229  // calculate maximum likelihood asuming an underlying poisson distribution, for details see for example
230  // http://nedwww.ipac.caltech.edu/level5/Leo/Stats_contents.html
231  /*
232  poison distribution : p(r) = y^r * exp(-y) / r! y = mean
233  likelihood distribution to be maximized L, n number of measurements (x1....xn) :
234  L = -n*y + Sum(xi * ln(y)) - Sum(ln(xi!)) = Sum( -y + xi * ln(y) - ln(xi!) )
235  */
236 
237  // some definitions
238  double MLH = 0.0 ;
239 
240  // loop over Internal Array and caluclate chi2
241  for (Int_t index = 0 ; index < mInternalArraySize ; index++)
242  {
243  // calculate the expected value at this point (variablePosition) for this set of paramters
244  double expectedValue = mCorrelationFunction( mVariablePositionArray[index] , parameterSet) * mDenominatorInternalArray[index] ;
245  // measured value
246  double measueredValue = mNumeratorInternalArray[index] ;
247  // sum up
248  MLH += -expectedValue + measueredValue * ::log(expectedValue) - lnfactorial(measueredValue) ;
249  }
250  // likelihood for this set of paramters
251  finalMLH = - MLH ;
252 
253  // count number of calls from Minuit
254  countMinuitCalls++ ;
255 }
256 //________________________________
257 double dFitter3d::lnfactorial(double arg)
258 {
259  // return the factorial of arg
260  double fac = 1.0 ;
261  int iarg = (int) arg ;
262  if(iarg<50)
263  {
264  for (int index = 1; index < iarg+1 ; index++)
265  { fac *=(double)index; } ;
266  fac = ::log(fac) ;
267  }
268  else
269  {
270  //if the argument is too large : use Stirlings formula
271  fac = 0.91 + (iarg+0.5) * ::log(double(iarg)) - iarg ;
272  }
273 
274  return fac;
275 }
276 //____________________________
277 double dFitter3d::mCorrelationFunction(TVector3& position, double* parameterSet )
278 {
279  if ( mCorrFctnType == "ykp")
280  {
281  return ykpCorrelationFunction( position, parameterSet ) ;
282  }
283  else if ( mCorrFctnType == "bp")
284  {
285  return bpCorrelationFunction( position, parameterSet ) ;
286  }
287  else
288  {
289  // this should never happen
290  cout << "this should never happen.\n " ;
291  return 0.0 ;
292  }
293 }
294 //____________________________
295 double dFitter3d::ykpCorrelationFunction(TVector3& position, double* parameterSet )
296 {
297  // return C2 at 'position' with Parameterset 'parameterSet' for YKP parametrization
298 
299  // position.x = qlong
300  // position.y = qperp
301  // position.z = qnull
302  // parameterSet[0] = lambda
303  // parameterSet[1] = Rlong
304  // parameterSet[2] = Rperp
305  // parameterSet[3] = Rnull
306  // parameterSet[4] = beta
307  double gamma2 = 1.0/fabs(1.0-::pow(parameterSet[4],2)) ;
308  double mhcc = 0.038937929230 ;
309  double c2 = 1.0 + parameterSet[0] * TMath::Exp(
310  (
311  (-1.0) * ::pow(position.y(),2) * ::pow(parameterSet[2],2)
312  - gamma2 * ::pow((position.x() - parameterSet[4]*position.z()),2) * ::pow(parameterSet[1],2)
313  - gamma2 * ::pow((position.z() - parameterSet[4]*position.x()),2) * ::pow(parameterSet[3],2)
314  )
315  / mhcc );
316  return c2;
317 };
318 //____________________________
319 double dFitter3d::bpCorrelationFunction(TVector3& position, double* parameterSet )
320 {
321  // return C2 at 'position' with 'parameterSet' for Bertsch-Pratt parametrization
322 
323  // position.x = qlong
324  // position.y = qside
325  // position.z = qout
326  // parameterSet[0] = lambda
327  // parameterSet[1] = Rside
328  // parameterSet[2] = Rout
329  // parameterSet[3] = Rlong
330  // parameterSet[4] = Routlong
331  double c2 = 1.0 + parameterSet[0] * TMath::Exp(
332  (
333  (-1.0) * ::pow(position.y(),2) * ::pow(parameterSet[1],2)
334  - ::pow(position.z(),2) * ::pow(parameterSet[2],2)
335  - ::pow(position.x(),2) * ::pow(parameterSet[3],2)
336  - 2* position.x() * position.z() * ::pow(parameterSet[4],2)
337  )
338  / mhc2 );
339  return c2 ;
340 };
341 //____________________________
342 void dFitter3d::doFit()
343 {
344  // // Errorflag
345  int ierflg = 0 ;
346  double arglist[4];
347 
348  if ( mCorrFctnType == "ykp")
349  {
350  // set start values
351  double startValues[5] = { 0.2 , 5.0 , 7.0 , 6.0 , 0.5} ;
352  // set step size
353  double stepSize[5] = {0.1 ,0.1, 0.1 ,0.1 ,0.1 };
354  // set fit parameter names startvalue and stepsize
355  mMinuit->mnparm(0, "lambda", startValues[0], stepSize[0], 0,0,ierflg);
356  mMinuit->mnparm(1, "Rlong", startValues[1], stepSize[1], 0,0,ierflg);
357  mMinuit->mnparm(2, "Rperp", startValues[2], stepSize[2], 0,0,ierflg);
358  mMinuit->mnparm(3, "Rnull", startValues[3], stepSize[3], 0,0,ierflg);
359  mMinuit->mnparm(4, "beta", startValues[4], stepSize[4], 0,0,ierflg);
360  }
361  else if ( mCorrFctnType == "bp")
362  {
363  // set start values
364  double startValues[5] = { 0.5 , 6.0 , 6.0 , 6.0 , 6.0} ;
365  // set step size
366  double stepSize[5] = {0.1 ,0.1, 0.1 ,0.1 ,0.1 };
367  // set fit parameter names startvalue and stepsize
368  mMinuit->mnparm(0, "lambda", startValues[0], stepSize[0], 0,0,ierflg);
369  mMinuit->mnparm(1, "Rside", startValues[1], stepSize[1], 0,0,ierflg);
370  mMinuit->mnparm(2, "Rout", startValues[2], stepSize[2], 0,0,ierflg);
371  mMinuit->mnparm(3, "Rlong", startValues[3], stepSize[3], 0,0,ierflg);
372  mMinuit->mnparm(4, "Routlong", startValues[4], stepSize[4], 0,0,ierflg);
373  }
374  else
375  {
376  // this should never happen
377  cout << "this should never happen.\n " ;
378  return ;
379  }
380 
381 
382  // set function to call
383  mMinuit->SetObjectFit(this) ;
384  mMinuit->SetFCN(gfcn) ;
385 
386  // Call minuit with some instructions
387  //arglist[0] = 1; mMinuit->mnexcm("SET ERR", arglist ,1,ierflg); // set error ?
388  //arglist[0] = 1.0; mMinuit->mnexcm("SET GRAD", arglist ,1,ierflg); // set grad ?
389  //arglist[0] = 2; mMinuit->mnexcm("SET STRA", arglist ,1,ierflg); // set strategie ?
390  //arglist[0] = 12; mMinuit->mnexcm("CALL FCN", arglist ,1,ierflg); // call fcn
391 
392  // Fix some fit parameter to get the fit to converge
393  /*arglist[0] = 6; mMinuit->mnexcm("FIX", arglist ,1,ierflg);
394  arglist[0] = 7; mMinuit->mnexcm("FIX", arglist ,1,ierflg);
395  arglist[0] = 8; mMinuit->mnexcm("FIX", arglist ,1,ierflg);
396  arglist[0] = 9; mMinuit->mnexcm("FIX", arglist ,1,ierflg);
397  arglist[0] = 10; mMinuit->mnexcm("FIX", arglist ,1,ierflg);
398  arglist[0] = 11; mMinuit->mnexcm("FIX", arglist ,1,ierflg);*/
399 
400  // Call minuit again ... ???
401  //arglist[0] = 4; mMinuit->mnexcm("CALL FCN", arglist ,1,ierflg) ;
402  //arglist[0] = 1000 ;
403  //arglist[1] = 4; mMinuit->mnexcm("SEEK", arglist ,2,ierflg);
404  arglist[0] = 0; mMinuit->mnexcm("MIGRATE", arglist ,0,ierflg) ;
405  //arglist[0] = 0; mMinuit->mnexcm("MINOS", arglist ,0,ierflg) ;
406 
408  // get chi2pdf
410  // get parameter values
411  //double doubledummy1 = 0 ;
412  //double doubledummy2 = 0 ;
413  //double pars[5] ;
414 
415  //for(int i = 0 ; i < 5 ; i++)
416  // {
417  // mMinuit->GetParameter(i, doubledummy1, doubledummy2) ;
418  // pars[i] = doubledummy1;
419  // }
420  // call fcn
421  //int intdummy = 0 ;
422  //double* doubledummy3 ;
423  //double finalchi2 = 0.0;
424  //mfcn(intdummy, doubledummy3, finalchi2, pars, intdummy) ;
425 
426  //cout << " finalchi2 " << finalchi2 << endl ;
427  //cout << " chi2pdf = finalchi2 / InternalArraySize-5parameters = " << finalchi2/(double)(mInternalArraySize-5) << endl ;
428 
429  cout << "Number of to Minuit calls : " << countMinuitCalls << endl ;
430 }
431 
432 //________________________________
433 void dFitter3d::SetHistos(TH3D* numerator, TH3D* denominator)
434 {
435  // set the numerator and denominator
436  mNumerator = numerator ;
437  mDenominator = denominator ;
438 
439  // calculate the ratio
440  mRatio = (TH3D*) numerator->Clone() ;
441  mRatio->Reset() ;
442  mRatio->Divide(numerator,denominator,1.0,mNorm) ;
443 }
444 //____________________________
445 void dFitter3d::Bin1ToBin3(TH3D* histo, int bin, int& binx, int& biny, int& binz)
446 {
447  // // actually this should be provided by ROOT ... grrrr
448  // return bin number of every axis (binx,biny,binz) if the general bin number (bin) is known
449  // general bin = binX + binY * (nbinX+2) + binZ * ((nbinx+2) * (nbiny+2) ) (+2 acounts for over and underflow!)
450  int nbbinX = (histo->GetNbinsX()) + 2 ;
451  int nbbinY = (histo->GetNbinsY()) + 2 ;
452 
453  binx = bin%nbbinX ;
454  biny = (bin/nbbinX) % nbbinY ;
455  binz = (bin/nbbinX) / nbbinY ;
456 
457 }
458 //____________________________
459 void dFitter3d::setCorrFctn(TString opt)
460 {
461 
462 };
463 //____________________________
464 void dFitter3d::setFitMethod(TString opt2)
465 {
466  // set typ of fitting MML or Chi2
467  if ( opt2 == "chi2" )
468  {
469  mFitMethod = "chi2" ;
470  cout << "Using chi2.\n" ;
471  }
472  else if ( opt2 == "mml" )
473  {
474  mFitMethod = "mml" ;
475  cout << "Using mml.\n" ;
476  }
477  else
478  {
479  mFitMethod = "chi2" ;
480  cout << "Bad option : Using chi2 method.\n" ;
481  } ;
482 
483 };
Generated by   doxygen 1.8.5