Go4Fit/TGo4FitModelGaussN.cxx
Go to the documentation of this file.
00001 // $Id: TGo4FitModelGaussN.cxx 933 2013-01-29 15:27:58Z linev $ 00002 //----------------------------------------------------------------------- 00003 // The GSI Online Offline Object Oriented (Go4) Project 00004 // Experiment Data Processing at EE department, GSI 00005 //----------------------------------------------------------------------- 00006 // Copyright (C) 2000- GSI Helmholtzzentrum für Schwerionenforschung GmbH 00007 // Planckstr. 1, 64291 Darmstadt, Germany 00008 // Contact: http://go4.gsi.de 00009 //----------------------------------------------------------------------- 00010 // This software can be used under the license agreements as stated 00011 // in Go4License.txt file which is part of the distribution. 00012 //----------------------------------------------------------------------- 00013 00014 #include "TGo4FitModelGaussN.h" 00015 00016 #include "Riostream.h" 00017 #include "TMath.h" 00018 #include "TGo4FitParameter.h" 00019 00020 TGo4FitModelGaussN::TGo4FitModelGaussN() : TGo4FitModel(), fxIndexes(), 00021 Vect_mu(0),Matr_sig(0), Vect_x(0), Vect_dx(0) { 00022 } 00023 00024 TGo4FitModelGaussN::TGo4FitModelGaussN(const char* iName, Int_t iNDimension) : 00025 TGo4FitModel(iName,"N-dimensional Gaussian", kTRUE), fxIndexes(iNDimension), 00026 Vect_mu(0),Matr_sig(0), Vect_x(0), Vect_dx(0) { 00027 for(Int_t n=0;n<iNDimension;n++) fxIndexes[n] = n; 00028 00029 for(Int_t n=0;n<iNDimension;n++) 00030 NewParameter(GetPosParName(n),"position",1.); 00031 00032 for(Int_t n=0;n<iNDimension;n++) 00033 NewParameter(GetWidthParName(n),"width",.5); 00034 00035 for(Int_t n1=0;n1<iNDimension;n1++) 00036 for(Int_t n2=n1+1;n2<iNDimension;n2++) 00037 NewParameter(GetCovarParName(n1,n2),"covariation",0.); 00038 } 00039 00040 TGo4FitModelGaussN::~TGo4FitModelGaussN() { 00041 if(Vect_mu) { delete Vect_mu; Vect_mu = 0; } 00042 if(Matr_sig) { delete Matr_sig; Matr_sig = 0; } 00043 if(Vect_x) { delete Vect_x; Vect_x = 0; } 00044 if (Vect_dx) { delete Vect_dx; Vect_dx = 0; } 00045 } 00046 00047 TString TGo4FitModelGaussN::GetPosParName(Int_t naxis) 00048 { 00049 TString res; 00050 res.Form("Pos%d",naxis); 00051 return res; 00052 } 00053 00054 TString TGo4FitModelGaussN::GetWidthParName(Int_t naxis) 00055 { 00056 TString res; 00057 res.Form("Width%d",naxis); 00058 return res; 00059 } 00060 00061 TString TGo4FitModelGaussN::GetCovarParName(Int_t naxis1, Int_t naxis2) 00062 { 00063 TString res; 00064 res.Form("Cov%d_%d",naxis1,naxis2); 00065 return res; 00066 } 00067 00068 00069 00070 Bool_t TGo4FitModelGaussN::SetAxisNumbers(Int_t naxis) { 00071 Int_t oldnaxis = GetAxisNumbers(); 00072 if ((naxis<0) || (oldnaxis==naxis)) return kFALSE; 00073 if(naxis<oldnaxis) { 00074 for(Int_t n=naxis;n<oldnaxis;n++) { 00075 RemovePar(GetPosParName(n)); 00076 RemovePar(GetWidthParName(n)); 00077 } 00078 for(Int_t n1=0;n1<oldnaxis;n1++) 00079 for(Int_t n2=n1+1;n2<oldnaxis;n2++) 00080 if((n1>=naxis) || (n2>=naxis)) 00081 RemovePar(GetCovarParName(n1,n2)); 00082 fxIndexes.Set(naxis); 00083 } else { 00084 Int_t lastindx = 0; 00085 for(Int_t n=oldnaxis;n<naxis;n++) { 00086 lastindx = GetParIndex(FindPar(GetPosParName(n-1)))+1; 00087 NewParameter(GetPosParName(n),"position",1.,kFALSE,lastindx); 00088 lastindx = GetParIndex(FindPar(GetWidthParName(n-1)))+1; 00089 NewParameter(GetWidthParName(n),"width",5.,kFALSE,lastindx); 00090 } 00091 lastindx++; 00092 for(Int_t n1=0;n1<naxis;n1++) 00093 for(Int_t n2=n1+1;n2<naxis;n2++) { 00094 TGo4FitParameter* par = FindPar(GetCovarParName(n1,n2)); 00095 if (par) lastindx = GetParIndex(par)+1; 00096 else NewParameter(GetCovarParName(n1,n2),"covariation",0.,kFALSE,lastindx++); 00097 } 00098 00099 fxIndexes.Set(naxis); 00100 for (Int_t n=oldnaxis;n<naxis;n++) { 00101 Int_t test = 0; 00102 Bool_t find; 00103 do { 00104 find = kFALSE; 00105 for (Int_t n1=0;n1<n;n1++) 00106 if(test==fxIndexes[n1]) { find=kTRUE; test++; break; } 00107 } while(find); 00108 fxIndexes[n] = test; 00109 } 00110 } 00111 00112 return kTRUE; 00113 } 00114 00115 Bool_t TGo4FitModelGaussN::ResortIndexes(Int_t leaveaxis) { 00116 Bool_t res = kFALSE; 00117 for (Int_t n=fxIndexes.GetSize()-1;n>=0;n--) { 00118 if (n==leaveaxis) continue; 00119 Int_t test = fxIndexes[n]; 00120 Bool_t find, first = kTRUE; 00121 do { 00122 find = kFALSE; 00123 for (Int_t n1=0;n1<fxIndexes.GetSize();n1++) 00124 if((n!=n1) && (test==fxIndexes[n1])) { 00125 if (first) test=0; else test++; 00126 find = kTRUE; 00127 first = kFALSE; 00128 break; 00129 } 00130 } while(find); 00131 fxIndexes[n] = test; 00132 res = res || (!first); 00133 } 00134 return res; 00135 } 00136 00137 00138 Int_t TGo4FitModelGaussN::GetPosParIndex(Int_t naxis) { 00139 for(Int_t n=0;n<fxIndexes.GetSize();n++) 00140 if (naxis==fxIndexes[n]) return 1+n; 00141 return -1; 00142 } 00143 00144 Int_t TGo4FitModelGaussN::GetWidthParIndex(Int_t naxis) { 00145 for(Int_t n=0;n<fxIndexes.GetSize();n++) 00146 if (naxis==fxIndexes[n]) return 1+fxIndexes.GetSize()+n; 00147 return -1; 00148 } 00149 00150 void TGo4FitModelGaussN::FillMuVector(TVectorD& Mu) { 00151 Int_t ndim = fxIndexes.GetSize(); 00152 Mu.ResizeTo(ndim); 00153 for(Int_t n=0;n<ndim;n++) 00154 Mu(n) = GetPar(1+n)->GetValue(); 00155 } 00156 00157 void TGo4FitModelGaussN::FillSigmaMatrix(TMatrixD& Sigma) 00158 { 00159 Int_t ndim = fxIndexes.GetSize(); 00160 Sigma.ResizeTo(ndim,ndim); 00161 for(Int_t n=0;n<ndim;n++) 00162 Sigma(n,n) = GetPar(1+ndim+n)->GetValue() * GetPar(1+ndim+n)->GetValue(); 00163 Int_t indx = 1+2*ndim; 00164 for(Int_t n1=0;n1<ndim;n1++) 00165 for(Int_t n2=n1+1;n2<ndim;n2++) { 00166 Double_t zn = GetPar(indx++)->GetValue(); 00167 Sigma(n1,n2) = zn; 00168 Sigma(n2,n1) = zn; 00169 } 00170 } 00171 00172 Bool_t TGo4FitModelGaussN::BeforeEval(Int_t NDimension) 00173 { 00174 Par_ndim = fxIndexes.GetSize(); 00175 Par_indx = fxIndexes.GetArray(); 00176 for (Int_t i=0;i<Par_ndim;i++) 00177 if (Par_indx[i]>=NDimension) 00178 { std::cout << "TGo4FitModelGaussN:: invalid index " << std::endl; return kFALSE; } 00179 00180 Vect_mu = new TVectorD(Par_ndim); 00181 Matr_sig = new TMatrixD(Par_ndim,Par_ndim); 00182 Vect_x = new TVectorD(Par_ndim); 00183 Vect_dx = new TVectorD(Par_ndim); 00184 00185 FillMuVector(*Vect_mu); 00186 FillSigmaMatrix(*Matr_sig); 00187 Double_t determ; 00188 Matr_sig->Invert(&determ); 00189 if (determ==0.) 00190 { std::cout << "TGo4FitModelGaussN:: Invalid sigma matrice " << std::endl; AfterEval(); return kFALSE; } 00191 00192 return kTRUE; 00193 } 00194 00195 Double_t TGo4FitModelGaussN::EvalN(const Double_t* v) { 00196 for (Int_t n=0;n<Par_ndim;n++) 00197 (*Vect_x)(n) = v[Par_indx[n]]; 00198 *Vect_x -= *Vect_mu; 00199 *Vect_dx = *Vect_x; 00200 *Vect_dx *= *Matr_sig; 00201 Double_t z = *Vect_dx * *Vect_x; // see 11.1 00202 return TMath::Exp(-0.5*z); 00203 } 00204 00205 void TGo4FitModelGaussN::AfterEval() { 00206 if(Vect_mu) { delete Vect_mu; Vect_mu = 0; } 00207 if(Matr_sig) { delete Matr_sig; Matr_sig = 0; } 00208 if(Vect_x) { delete Vect_x; Vect_x = 0; } 00209 if (Vect_dx) { delete Vect_dx; Vect_dx = 0; } 00210 } 00211 00212 void TGo4FitModelGaussN::Print(Option_t* option) const 00213 { 00214 TGo4FitModel::Print(option); 00215 std::cout << " N-dimensional Gauss" << std::endl; 00216 std::cout << " axis indexes:"; 00217 for (Int_t i=0;i<fxIndexes.GetSize();i++) 00218 std::cout << " " << fxIndexes[i]; 00219 std::cout << std::endl; 00220 }
