GSI Object Oriented Online Offline (Go4)  GO4-6.1.4
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
TGo4FitAmplEstimation.cxx
Go to the documentation of this file.
1 // $Id: TGo4FitAmplEstimation.cxx 2769 2020-04-16 14:54:23Z linev $
2 //-----------------------------------------------------------------------
3 // The GSI Online Offline Object Oriented (Go4) Project
4 // Experiment Data Processing at EE department, GSI
5 //-----------------------------------------------------------------------
6 // Copyright (C) 2000- GSI Helmholtzzentrum fuer Schwerionenforschung GmbH
7 // Planckstr. 1, 64291 Darmstadt, Germany
8 // Contact: http://go4.gsi.de
9 //-----------------------------------------------------------------------
10 // This software can be used under the license agreements as stated
11 // in Go4License.txt file which is part of the distribution.
12 //-----------------------------------------------------------------------
13 
14 #include "TGo4FitAmplEstimation.h"
15 
16 #include <iostream>
17 
18 #include "TMath.h"
19 #include "TMatrixD.h"
20 #include "TVectorD.h"
21 #include "TArrayC.h"
22 #include "TArrayD.h"
23 
24 #include "TGo4Fitter.h"
25 #include "TGo4FitModel.h"
26 #include "TGo4FitData.h"
27 #include "TGo4FitParameter.h"
28 
30 }
31 
32 TGo4FitAmplEstimation::TGo4FitAmplEstimation(const char* Name, Int_t NumIters) :
33  TGo4FitterAction(Name,"Estimate amplitude action"), fiNumIters(NumIters) {
34 }
35 
37 }
38 
40  TGo4Fitter* f = dynamic_cast<TGo4Fitter*> (Fitter);
41  if(f==0) return;
42  if (!CalculateWithBuffers(f))
44 }
45 
46 Double_t TGo4FitAmplEstimation::PointWeight(Int_t niter, Int_t FFtype, Double_t value, Double_t modelvalue, Double_t standdev) {
47  switch (FFtype) {
48  case TGo4Fitter::ff_least_squares: return 1.;
49  case TGo4Fitter::ff_chi_square: if (standdev<=0.) return 0.; else return 1./standdev;
52  if (niter==0) if (value<0.) return 0.; else return 1./(value+1.);
53  else if (modelvalue<=0.) return 0.; else return 1./modelvalue;
54  case TGo4Fitter::ff_chi_Neyman: if (value<1.) return 1.; else return 1./value;
55  case TGo4Fitter::ff_chi_gamma: if (value<0.) return 0.; else return 1./(value+1.);
56  case TGo4Fitter::ff_user: return 1.;
57  default: return 1.;
58  }
59 }
60 
61 
63  if (!fitter->IsInitialized()) return kFALSE;
64 
65  for(Int_t n=0;n<fitter->GetNumData();n++)
66  if (!fitter->DataBuffersAllocated(fitter->GetData(n))) return kFALSE;
67 
68  for(Int_t n=0;n<fitter->GetNumModel();n++)
69  if (!fitter->ModelBuffersAllocated(fitter->GetModel(n))) return kFALSE;
70 
71  TObjArray comps, fixedcomps;
72 
73  for(Int_t n=0;n<fitter->GetNumModel();n++) {
74  TGo4FitModel *fModel = fitter->GetModel(n);
75 
76  Bool_t assigned = kFALSE;
77  for (Int_t na=0;na<fModel->NumAssigments();na++)
78  if (fitter->FindData(fModel->AssignmentName(na))) {
79  assigned = kTRUE;
80  break;
81  }
82  if (!assigned) continue;
83 
84  if ( (fModel->GetAmplPar()!=0) &&
85  (!fModel->GetAmplPar()->GetFixed()) ) comps.Add(fModel);
86  else fixedcomps.Add(fModel);
87  }
88 
89  Int_t ncomp = comps.GetLast()+1;
90  if (ncomp<=0) return kFALSE;
91 
92  TMatrixD matr(ncomp,ncomp);
93  TVectorD vector(ncomp);
94 
95  Double_t** DopData = new Double_t* [fitter->GetNumData()];
96  for(Int_t n=0;n<fitter->GetNumData();n++)
97  DopData[n] = new Double_t[fitter->GetDataBinsSize(fitter->GetData(n))];
98 
99  Double_t** Weights = new Double_t* [fitter->GetNumData()];
100  for(Int_t n=0;n<fitter->GetNumData();n++)
101  Weights[n] = new Double_t[fitter->GetDataBinsSize(fitter->GetData(n))];
102 
103 
104  Int_t niter = 0;
105  Int_t FFtype = fitter->GetFitFunctionType();
106  Bool_t changed = kFALSE;
107 
108  do {
109 
110  fitter->RebuildAll();
111 
112  for(Int_t ndata=0;ndata<fitter->GetNumData();ndata++) {
113  TGo4FitData* data = fitter->GetData(ndata);
114 
115  Int_t size = fitter->GetDataBinsSize(data);
116  Double_t ampl = data->GetAmplValue();
117 
118  Double_t* values = fitter->GetDataBinsValues(data);
119  Double_t* result = fitter->GetDataBinsResult(data);
120  Double_t* devs = fitter->GetDataBinsDevs(data);
121  Double_t* weight = Weights[ndata];
122 
123  for (Int_t nbin=0;nbin<size;nbin++)
124  weight[nbin] = PointWeight(niter,FFtype, ampl*values[nbin], result[nbin], ampl*ampl*devs[nbin]);
125  }
126 
127  for(Int_t n1=0;n1<ncomp;n1++)
128  for(Int_t n2=n1;n2<ncomp;n2++) {
129  Double_t fSum = 0.;
130  for(Int_t ndata=0;ndata<fitter->GetNumData();ndata++) {
131  Double_t* bins1 = fitter->GetModelBinsValues((TGo4FitModel*) comps[n1], fitter->GetDataName(ndata));
132  Double_t* bins2 = fitter->GetModelBinsValues((TGo4FitModel*) comps[n2], fitter->GetDataName(ndata));
133  if ((bins1!=0) && (bins2!=0)) {
134  Int_t size = fitter->GetDataBinsSize(fitter->GetData(ndata));
135  Double_t* weight = Weights[ndata];
136  for(Int_t nbin=0;nbin<size;nbin++)
137  fSum+=bins1[nbin]*bins2[nbin]*weight[nbin];
138  }
139  }
140  matr(n1,n2) = fSum;
141  matr(n2,n1) = fSum;
142  }
143 
144  // calculating difference between data and fixed components
145  for(Int_t ndata=0;ndata<fitter->GetNumData();ndata++) {
146  TGo4FitData* data = fitter->GetData(ndata);
147  Int_t size = fitter->GetDataBinsSize(data);
148  Double_t dampl = data->GetAmplValue();
149  Double_t* values = fitter->GetDataBinsValues(data);
150  Double_t* bins = DopData[ndata];
151 
152  for(Int_t nbin=0;nbin<size;nbin++)
153  bins[nbin] = dampl*values[nbin];
154 
155  for(Int_t n=0;n<=fixedcomps.GetLast();n++) {
156  TGo4FitModel *model = (TGo4FitModel*) fixedcomps[n];
157  Double_t* mbins = fitter->GetModelBinsValues(model,data->GetName());
158  if (mbins!=0) {
159  Double_t mampl = model->GetAmplValue();
160  for(Int_t nbin=0;nbin<size;nbin++)
161  bins[nbin] -= mampl*mbins[nbin];
162  }
163  }
164  }
165 
166  // filling vector
167 
168  for(Int_t n=0;n<ncomp;n++) {
169  Double_t fSum = 0.;
170  for(Int_t ndata=0;ndata<fitter->GetNumData();ndata++) {
171  Double_t* bins1 = fitter->GetModelBinsValues((TGo4FitModel*) comps[n], fitter->GetDataName(ndata));
172  if (bins1==0) continue;
173  Int_t size = fitter->GetDataBinsSize(fitter->GetData(ndata));
174  Double_t* bins = DopData[ndata];
175  Double_t* weight = Weights[ndata];
176  for(Int_t nbin=0;nbin<size;nbin++)
177  fSum+=bins[nbin]*bins1[nbin]*weight[nbin];
178 
179  }
180  vector(n) = fSum;
181  }
182 
183  if (matr.Determinant()==0.) {
184  std::cerr << " Amplitude estimation algorithm failed. " << std::endl;
185  std::cerr << " Determinant of matrix == 0.; This may be due to equivalent model components or zero model at all" << std::endl;
186  break;
187  }
188 
189  matr.Invert();
190  vector *= matr;
191 
192  for(Int_t n=0;n<ncomp;n++) {
193  TGo4FitModel *model = (TGo4FitModel*) comps[n];
194  model->SetAmplValue(vector(n));
195  if (matr(n,n)>0.)
196  model->SetAmplError(TMath::Sqrt(matr(n,n)));
197  }
198 
199  changed = kFALSE;
200  niter++;
201  if (niter<fiNumIters)
202  if ((FFtype==TGo4Fitter::ff_chi_Pearson) ||
203  (FFtype==TGo4Fitter::ff_ML_Poisson)) changed=kTRUE;
204 
205  } while (changed);
206 
207  for(Int_t n=0;n<fitter->GetNumData();n++) delete[] Weights[n];
208  delete[] Weights;
209 
210  for(Int_t n=0;n<fitter->GetNumData();n++) delete[] DopData[n];
211  delete[] DopData;
212 
213  return kTRUE;
214 }
215 
216 
218  if (fitter==0) return kFALSE;
219 
220  Int_t nmodel = fitter->GetNumModel();
221  TArrayI refindex(fitter->GetNumModel()); refindex.Reset(-1);
222  TArrayI fixedindex(fitter->GetNumModel()); fixedindex.Reset(-1);
223 
224  Int_t ncomp = 0;
225  Int_t nfixed = 0;
226 
227  for(Int_t n=0;n<nmodel;n++) {
228  TGo4FitModel *model = fitter->GetModel(n);
229  Bool_t assigned = kFALSE;
230  for (Int_t na=0;na<model->NumAssigments();na++)
231  if (fitter->FindData(model->AssignmentName(na))) {
232  assigned = kTRUE;
233  break;
234  }
235  if (!assigned) continue;
236 
237  if ( (model->GetAmplPar()!=0) &&
238  (!model->GetAmplPar()->GetFixed()) ) refindex[ncomp++] = n;
239  else fixedindex[nfixed++] = n;
240  }
241 
242  if (ncomp<=0) return kFALSE;
243 
244  TMatrixD matr(ncomp,ncomp);
245  TVectorD vector(ncomp);
246 
247  TArrayC Usage(nmodel);
248  TArrayD ModelValues(nmodel);
249  TArrayD ModelAmpls(nmodel);
250 
251  Int_t niter = 0;
252  Int_t FFtype = fitter->GetFitFunctionType();
253  Bool_t changed = kFALSE;
254 
255  do {
256  matr = 0.; vector = 0.;
257 
258  for(Int_t ndata=0;ndata<fitter->GetNumData();ndata++) {
259  TGo4FitData* data = fitter->GetData(ndata);
260  TGo4FitDataIter* iter = data->MakeIter();
261  if (iter==0) return kFALSE;
262 
263  Double_t dampl = data->GetAmplValue();
264 
265  if (iter->Reset()) {
266  Usage.Reset(0); ModelAmpls.Reset(0.); ModelValues.Reset(0.);
267  for(Int_t n=0;n<nmodel;n++) {
268  TGo4FitModel* model = fitter->GetModel(n);
269  if (model->IsAssignTo(data->GetName())) {
270  Usage[n] = 1;
271  ModelAmpls[n] = model->GetAmplValue() * model->GetRatioValueFor(data->GetName());
272  model->BeforeEval(iter->IndexesSize());
273  }
274  }
275 
276  do {
277  Double_t result = 0.;
278  for (Int_t n=0; n<nmodel;n++)
279  if (Usage[n]) {
280  ModelValues[n] = fitter->GetModel(n)->EvaluateAtPoint(iter);
281  result += ModelAmpls[n] * ModelValues[n];
282  }
283  Double_t weight = PointWeight(niter, FFtype, dampl*iter->Value(), result, dampl*dampl*iter->StandardDeviation());
284  if (weight<=0.) continue;
285 
286  for (Int_t i1=0;i1<ncomp;i1++)
287  for(Int_t i2=i1;i2<ncomp;i2++) {
288  Double_t zn = matr(i1,i2);
289  zn += ModelValues[refindex[i1]] * ModelValues[refindex[i2]] * weight;
290  matr(i1,i2) = zn;
291  if (i1!=i2) matr(i2,i1) = zn;
292  }
293 
294  Double_t dvalue = dampl*iter->Value();
295  for(Int_t nf=0;nf<nfixed;nf++)
296  dvalue -= ModelAmpls[fixedindex[nf]] * ModelValues[fixedindex[nf]];
297 
298  for (Int_t i1=0;i1<ncomp;i1++) {
299  Double_t zn = vector(i1);
300  zn += ModelValues[refindex[i1]] * dvalue * weight;
301  vector(i1) = zn;
302  }
303 
304  } while (iter->Next());
305 
306  for(Int_t n=0;n<nmodel;n++)
307  if (Usage[n]) fitter->GetModel(n)->AfterEval();
308  }
309  delete iter;
310  }
311 
312  if (matr.Determinant()==0.) {
313  std::cerr << " Amplitude estimation algorithm failed. " << std::endl;
314  std::cerr << " Determinant of matrix == 0.; This may be due to equivalent model components or zero model at all" << std::endl;
315  return kFALSE;
316  }
317 
318  matr.Invert();
319  vector *= matr;
320 
321  for(Int_t n=0;n<ncomp;n++) {
322  TGo4FitModel *model = fitter->GetModel(refindex[n]);
323  model->SetAmplValue(vector(n));
324  if (matr(n,n)>0.)
325  model->SetAmplError(TMath::Sqrt(matr(n,n)));
326  }
327 
328  changed = kFALSE;
329  niter++;
330  if (niter<fiNumIters)
331  if ((FFtype==TGo4Fitter::ff_chi_Pearson) ||
332  (FFtype==TGo4Fitter::ff_ML_Poisson)) changed = kTRUE;
333 
334  } while (changed);
335 
336  return kTRUE;
337 }
338 
339 void TGo4FitAmplEstimation::Print(Option_t* option) const
340 {
341  TGo4FitterAction::Print(option);
342  std::cout << " number of iterations " << fiNumIters << std::endl;
343 }
virtual void DoAction(TGo4FitterAbstract *Fitter)
Bool_t GetFixed() const
void RebuildAll(Bool_t ForceBuild=kFALSE)
Definition: TGo4Fitter.cxx:625
Double_t PointWeight(Int_t niter, Int_t FFtype, Double_t value, Double_t modelvalue, Double_t standdev)
virtual TGo4FitDataIter * MakeIter()
Definition: TGo4FitData.h:165
Int_t GetDataBinsSize(TGo4FitData *data)
Definition: TGo4Fitter.cxx:680
Double_t * GetDataBinsResult(TGo4FitData *data)
Definition: TGo4Fitter.cxx:695
void SetAmplError(Double_t iError)
virtual void Print(Option_t *option) const
void SetAmplValue(Double_t iAmpl)
const char * AssignmentName(Int_t n)
Definition: TGo4FitModel.h:129
Bool_t CalculateWithIterators(TGo4Fitter *fitter)
Double_t * GetDataBinsDevs(TGo4FitData *data)
Definition: TGo4Fitter.cxx:690
virtual Bool_t Reset(Bool_t UseRanges=kTRUE)
virtual Double_t EvaluateAtPoint(TGo4FitData *data, Int_t nbin, Bool_t UseRanges=kTRUE)
void Print(Option_t *option) const
TGo4FitModel * GetModel(Int_t n)
Definition: TGo4Fitter.cxx:191
virtual Bool_t Next(Bool_t UseRanges=kTRUE)
virtual Bool_t BeforeEval(Int_t ndim)
virtual void AfterEval()
Definition: TGo4FitModel.h:253
TGo4FitData * GetData(Int_t n)
Definition: TGo4Fitter.cxx:101
Int_t NumAssigments() const
Definition: TGo4FitModel.h:124
Int_t GetFitFunctionType()
Definition: TGo4Fitter.h:85
Double_t * GetModelBinsValues(TGo4FitModel *model, const char *DataName)
Definition: TGo4Fitter.cxx:700
Bool_t DataBuffersAllocated(TGo4FitData *data)
Definition: TGo4Fitter.cxx:675
Bool_t CalculateWithBuffers(TGo4Fitter *fitter)
const char * GetDataName(Int_t n)
Definition: TGo4Fitter.cxx:106
Double_t Value() const
Definition: TGo4FitData.h:514
Int_t GetNumModel() const
Definition: TGo4Fitter.h:178
Bool_t IsAssignTo(const char *DataName) const
Definition: TGo4FitModel.h:141
Double_t * GetDataBinsValues(TGo4FitData *data)
Definition: TGo4Fitter.cxx:685
TGo4FitData * FindData(const char *DataName)
Definition: TGo4Fitter.cxx:111
Int_t IndexesSize() const
Definition: TGo4FitData.h:459
Double_t StandardDeviation() const
Definition: TGo4FitData.h:519
Double_t GetRatioValueFor(const char *DataName)
Bool_t ModelBuffersAllocated(TGo4FitModel *model)
Definition: TGo4Fitter.cxx:670
TGo4FitParameter * GetAmplPar()
Int_t GetNumData() const
Definition: TGo4Fitter.h:123