GSI Object Oriented Online Offline (Go4)  GO4-6.1.4
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
TGo4FitModelGaussN.cxx
Go to the documentation of this file.
1 // $Id: TGo4FitModelGaussN.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 "TGo4FitModelGaussN.h"
15 
16 #include <iostream>
17 
18 #include "TMath.h"
19 #include "TGo4FitParameter.h"
20 
22  Vect_mu(0),Matr_sig(0), Vect_x(0), Vect_dx(0) {
23 }
24 
25 TGo4FitModelGaussN::TGo4FitModelGaussN(const char* iName, Int_t iNDimension) :
26  TGo4FitModel(iName,"N-dimensional Gaussian", kTRUE), fxIndexes(iNDimension),
27  Vect_mu(0),Matr_sig(0), Vect_x(0), Vect_dx(0) {
28  for(Int_t n=0;n<iNDimension;n++) fxIndexes[n] = n;
29 
30  for(Int_t n=0;n<iNDimension;n++)
31  NewParameter(GetPosParName(n),"position",1.);
32 
33  for(Int_t n=0;n<iNDimension;n++)
34  NewParameter(GetWidthParName(n),"width",.5);
35 
36  for(Int_t n1=0;n1<iNDimension;n1++)
37  for(Int_t n2=n1+1;n2<iNDimension;n2++)
38  NewParameter(GetCovarParName(n1,n2),"covariation",0.);
39 }
40 
42  if(Vect_mu) { delete Vect_mu; Vect_mu = 0; }
43  if(Matr_sig) { delete Matr_sig; Matr_sig = 0; }
44  if(Vect_x) { delete Vect_x; Vect_x = 0; }
45  if (Vect_dx) { delete Vect_dx; Vect_dx = 0; }
46 }
47 
49 {
50  TString res;
51  res.Form("Pos%d",naxis);
52  return res;
53 }
54 
56 {
57  TString res;
58  res.Form("Width%d",naxis);
59  return res;
60 }
61 
62 TString TGo4FitModelGaussN::GetCovarParName(Int_t naxis1, Int_t naxis2)
63 {
64  TString res;
65  res.Form("Cov%d_%d",naxis1,naxis2);
66  return res;
67 }
68 
69 
70 
72  Int_t oldnaxis = GetAxisNumbers();
73  if ((naxis<0) || (oldnaxis==naxis)) return kFALSE;
74  if(naxis<oldnaxis) {
75  for(Int_t n=naxis;n<oldnaxis;n++) {
78  }
79  for(Int_t n1=0;n1<oldnaxis;n1++)
80  for(Int_t n2=n1+1;n2<oldnaxis;n2++)
81  if((n1>=naxis) || (n2>=naxis))
82  RemovePar(GetCovarParName(n1,n2));
83  fxIndexes.Set(naxis);
84  } else {
85  Int_t lastindx = 0;
86  for(Int_t n=oldnaxis;n<naxis;n++) {
87  lastindx = GetParIndex(FindPar(GetPosParName(n-1)))+1;
88  NewParameter(GetPosParName(n),"position",1.,kFALSE,lastindx);
89  lastindx = GetParIndex(FindPar(GetWidthParName(n-1)))+1;
90  NewParameter(GetWidthParName(n),"width",5.,kFALSE,lastindx);
91  }
92  lastindx++;
93  for(Int_t n1=0;n1<naxis;n1++)
94  for(Int_t n2=n1+1;n2<naxis;n2++) {
96  if (par) lastindx = GetParIndex(par)+1;
97  else NewParameter(GetCovarParName(n1,n2),"covariation",0.,kFALSE,lastindx++);
98  }
99 
100  fxIndexes.Set(naxis);
101  for (Int_t n=oldnaxis;n<naxis;n++) {
102  Int_t test = 0;
103  Bool_t find;
104  do {
105  find = kFALSE;
106  for (Int_t n1=0;n1<n;n1++)
107  if(test==fxIndexes[n1]) { find=kTRUE; test++; break; }
108  } while(find);
109  fxIndexes[n] = test;
110  }
111  }
112 
113  return kTRUE;
114 }
115 
116 Bool_t TGo4FitModelGaussN::ResortIndexes(Int_t leaveaxis) {
117  Bool_t res = kFALSE;
118  for (Int_t n=fxIndexes.GetSize()-1;n>=0;n--) {
119  if (n==leaveaxis) continue;
120  Int_t test = fxIndexes[n];
121  Bool_t find, first = kTRUE;
122  do {
123  find = kFALSE;
124  for (Int_t n1=0;n1<fxIndexes.GetSize();n1++)
125  if((n!=n1) && (test==fxIndexes[n1])) {
126  if (first) test=0; else test++;
127  find = kTRUE;
128  first = kFALSE;
129  break;
130  }
131  } while(find);
132  fxIndexes[n] = test;
133  res = res || (!first);
134  }
135  return res;
136 }
137 
138 
140  for(Int_t n=0;n<fxIndexes.GetSize();n++)
141  if (naxis==fxIndexes[n]) return 1+n;
142  return -1;
143 }
144 
146  for(Int_t n=0;n<fxIndexes.GetSize();n++)
147  if (naxis==fxIndexes[n]) return 1+fxIndexes.GetSize()+n;
148  return -1;
149 }
150 
152  Int_t ndim = fxIndexes.GetSize();
153  Mu.ResizeTo(ndim);
154  for(Int_t n=0;n<ndim;n++)
155  Mu(n) = GetPar(1+n)->GetValue();
156 }
157 
159 {
160  Int_t ndim = fxIndexes.GetSize();
161  Sigma.ResizeTo(ndim,ndim);
162  for(Int_t n=0;n<ndim;n++)
163  Sigma(n,n) = GetPar(1+ndim+n)->GetValue() * GetPar(1+ndim+n)->GetValue();
164  Int_t indx = 1+2*ndim;
165  for(Int_t n1=0;n1<ndim;n1++)
166  for(Int_t n2=n1+1;n2<ndim;n2++) {
167  Double_t zn = GetPar(indx++)->GetValue();
168  Sigma(n1,n2) = zn;
169  Sigma(n2,n1) = zn;
170  }
171 }
172 
173 Bool_t TGo4FitModelGaussN::BeforeEval(Int_t NDimension)
174 {
175  Par_ndim = fxIndexes.GetSize();
176  Par_indx = fxIndexes.GetArray();
177  for (Int_t i=0;i<Par_ndim;i++)
178  if (Par_indx[i]>=NDimension)
179  { std::cout << "TGo4FitModelGaussN:: invalid index " << std::endl; return kFALSE; }
180 
181  Vect_mu = new TVectorD(Par_ndim);
182  Matr_sig = new TMatrixD(Par_ndim,Par_ndim);
183  Vect_x = new TVectorD(Par_ndim);
184  Vect_dx = new TVectorD(Par_ndim);
185 
188  Double_t determ;
189  Matr_sig->Invert(&determ);
190  if (determ==0.)
191  { std::cout << "TGo4FitModelGaussN:: Invalid sigma matrice " << std::endl; AfterEval(); return kFALSE; }
192 
193  return kTRUE;
194 }
195 
196 Double_t TGo4FitModelGaussN::EvalN(const Double_t* v) {
197  for (Int_t n=0;n<Par_ndim;n++)
198  (*Vect_x)(n) = v[Par_indx[n]];
199  *Vect_x -= *Vect_mu;
200  *Vect_dx = *Vect_x;
201  *Vect_dx *= *Matr_sig;
202  Double_t z = *Vect_dx * *Vect_x; // see 11.1
203  return TMath::Exp(-0.5*z);
204 }
205 
207  if(Vect_mu) { delete Vect_mu; Vect_mu = 0; }
208  if(Matr_sig) { delete Matr_sig; Matr_sig = 0; }
209  if(Vect_x) { delete Vect_x; Vect_x = 0; }
210  if (Vect_dx) { delete Vect_dx; Vect_dx = 0; }
211 }
212 
213 void TGo4FitModelGaussN::Print(Option_t* option) const
214 {
215  TGo4FitModel::Print(option);
216  std::cout << " N-dimensional Gauss" << std::endl;
217  std::cout << " axis indexes:";
218  for (Int_t i=0;i<fxIndexes.GetSize();i++)
219  std::cout << " " << fxIndexes[i];
220  std::cout << std::endl;
221 }
Bool_t ResortIndexes(Int_t leaveaxis=-1)
TString GetPosParName(Int_t naxis)
virtual Double_t EvalN(const Double_t *v)
TGo4FitParameter * GetPar(Int_t n)
void FillMuVector(TVectorD &Mu)
TGo4FitParameter * NewParameter(const char *Name, const char *Title, Double_t iValue=0., Bool_t Fixed=kFALSE, Int_t AtIndx=-1)
void FillSigmaMatrix(TMatrixD &Sigma)
Bool_t SetAxisNumbers(Int_t naxis)
virtual Int_t GetWidthParIndex(Int_t naxis)
Bool_t RemovePar(const char *name)
virtual Bool_t BeforeEval(Int_t)
TString GetCovarParName(Int_t naxis1, Int_t naxis2)
TGo4FitParameter * FindPar(const char *ParName)
virtual void Print(Option_t *option) const
TString GetWidthParName(Int_t naxis)
Int_t GetParIndex(const TGo4FitParameter *par)
Double_t GetValue() const
virtual Int_t GetPosParIndex(Int_t naxis)
virtual void Print(Option_t *option) const