Go4Fit/TGo4FitModel.cxx
Go to the documentation of this file.
00001 // $Id: TGo4FitModel.cxx 999 2013-07-25 11:58:59Z 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 "TGo4FitModel.h" 00015 00016 #include "Riostream.h" 00017 00018 #include "TArrayD.h" 00019 #include "TArrayI.h" 00020 #include "TClass.h" 00021 00022 #include "TGo4FitData.h" 00023 #include "TGo4FitParameter.h" 00024 00025 TGo4FitAssignment::TGo4FitAssignment() : 00026 TNamed(), fxRatio(0), fxData(0), fxModelMask(0), fxModelBins(0) { 00027 } 00028 00029 TGo4FitAssignment::TGo4FitAssignment(const char* DataName) : 00030 TNamed(DataName,""), 00031 fxRatio(0), fxData(0), fxModelMask(0), fxModelBins(0) { 00032 } 00033 00034 TGo4FitAssignment::~TGo4FitAssignment() { 00035 if (fxRatio) delete fxRatio; 00036 if (fxModelMask) delete[] fxModelMask; 00037 if (fxModelBins) delete[] fxModelBins; 00038 } 00039 00040 Double_t TGo4FitAssignment::RatioValue() { 00041 return (fxRatio==0) ? 1. : fxRatio->GetValue(); 00042 } 00043 00044 void TGo4FitAssignment::Print(Option_t* option) const 00045 { 00046 std::cout << " " << GetName(); 00047 } 00048 00049 // ************************************************************************************************** 00050 00051 TGo4FitModel::TGo4FitModel() : TGo4FitComponent(), 00052 fiMinIntegrDepth(0), fiMaxIntegrDepth(0), fdIntegrEps(0.), fbAbsoluteEps(kFALSE), fbIntegrScaling(kFALSE), 00053 fxAssigments(), fiGroupIndex(-1), 00054 fxCurrentPars(), fxCurrentParsArray(0), 00055 fbNeedToRebuild(kFALSE), fxAllPars(0), fxAllParsValues(0) { 00056 } 00057 00058 TGo4FitModel::TGo4FitModel(const char* iName, const char* iTitle, Bool_t MakeAmplitude) : 00059 TGo4FitComponent(iName,iTitle), 00060 fiMinIntegrDepth(0), fiMaxIntegrDepth(0), fdIntegrEps(0.), fbAbsoluteEps(kFALSE), fbIntegrScaling(kFALSE), 00061 fxAssigments(), fiGroupIndex(-1), 00062 fxCurrentPars(), fxCurrentParsArray(0), 00063 fbNeedToRebuild(kFALSE), fxAllPars(0), fxAllParsValues(0) { 00064 fxAssigments.SetOwner(kTRUE); 00065 if (MakeAmplitude) NewAmplitude(); 00066 } 00067 00068 TGo4FitModel::~TGo4FitModel() { 00069 RemoveAllPars(); 00070 } 00071 00072 Int_t TGo4FitModel::NumPars() { 00073 Int_t res = TGo4FitParsList::NumPars(); 00074 if (NumAssigments()>0) res += NumAssigments()-1; 00075 return res; 00076 } 00077 00078 TGo4FitParameter* TGo4FitModel::Get(Int_t n) { 00079 Int_t numpars = TGo4FitParsList::NumPars(); 00080 if (n<numpars) return TGo4FitParsList::Get(n); 00081 n-=numpars-1; 00082 return (n>=0) && (n<NumAssigments()) ? GetAssigment(n)->fxRatio : 0; 00083 } 00084 00085 void TGo4FitModel::AssignToData(const char* DataName, Double_t RatioValue, Bool_t FixRatio) { 00086 if (FindAssigment(DataName)==0) { 00087 TGo4FitAssignment* ass = new TGo4FitAssignment(DataName); 00088 if (NumAssigments()>0) { 00089 ass->fxRatio = new TGo4FitParameter(GetRatioName(NumAssigments()), 00090 "Amplitude ratio to first data", 00091 RatioValue); 00092 ass->fxRatio->SetOwner(this); 00093 if (FixRatio) ass->fxRatio->SetFixed(kTRUE); 00094 } 00095 fxAssigments.Add(ass); 00096 } 00097 } 00098 00099 void TGo4FitModel::ChangeDataNameInAssignments(const char* oldname, const char* newname) { 00100 TGo4FitAssignment* ass = FindAssigment(oldname); 00101 if (ass) ass->SetName(newname); 00102 } 00103 00104 void TGo4FitModel::ClearAssignmentTo(const char* DataName) { 00105 TGo4FitAssignment* ass = FindAssigment(DataName); 00106 if (ass==0) return; 00107 00108 fxAssigments.Remove(ass); 00109 delete ass; 00110 fxAssigments.Compress(); 00111 00112 for (Int_t n=0;n<NumAssigments();n++) { 00113 ass = GetAssigment(n); 00114 if (ass->fxRatio==0) continue; 00115 if (n==0) { 00116 delete ass->fxRatio; 00117 ass->fxRatio = 0; 00118 } else 00119 ass->fxRatio->SetName(GetRatioName(n+1)); 00120 } 00121 } 00122 00123 void TGo4FitModel::ClearAssignments() { 00124 fxAssigments.Clear(); 00125 fxAssigments.Compress(); 00126 } 00127 00128 void TGo4FitModel::ConnectToDataIfAssigned(TGo4FitData* data) { 00129 if (data==0) return; 00130 TGo4FitAssignment* ass = FindAssigment(data->GetName()); 00131 if ( ass != 0 ) ass->fxData = data; 00132 } 00133 00134 TGo4FitData* TGo4FitModel::GetAssignedConnection(Int_t n) { 00135 TGo4FitAssignment* ass = GetAssigment(n); 00136 if(ass) return ass->fxData; 00137 else return 0; 00138 } 00139 00140 Bool_t TGo4FitModel::GetPosition(Int_t naxis, Double_t& pos) { 00141 if(GetPosPar(naxis)) { 00142 pos = GetPosPar(naxis)->GetValue(); 00143 return kTRUE; 00144 } else return kFALSE; 00145 } 00146 00147 Bool_t TGo4FitModel::SetPosition(Int_t naxis, Double_t pos) { 00148 if(GetPosPar(naxis)) { 00149 GetPosPar(naxis)->SetValue(pos); 00150 return kTRUE; 00151 } else return kFALSE; 00152 } 00153 00154 Bool_t TGo4FitModel::GetWidth(Int_t naxis, Double_t& width) { 00155 if(GetWidthPar(naxis)) { 00156 width = GetWidthPar(naxis)->GetValue(); 00157 return kTRUE; 00158 } else return kFALSE; 00159 } 00160 00161 Bool_t TGo4FitModel::SetWidth(Int_t naxis, Double_t width) { 00162 if(GetWidthPar(naxis)) { 00163 GetWidthPar(naxis)->SetValue(width); 00164 return kTRUE; 00165 } else return kFALSE; 00166 } 00167 00168 void TGo4FitModel::SetIntegrationsProperty(Int_t iMinIntegrDepth, Int_t iMaxIntegrDepth, Double_t iIntegrEps, Bool_t iAbsoluteEps, Bool_t iIntegrScaling) { 00169 fiMinIntegrDepth = iMinIntegrDepth; 00170 fiMaxIntegrDepth = iMaxIntegrDepth; 00171 fdIntegrEps = iIntegrEps; 00172 fbAbsoluteEps = iAbsoluteEps; 00173 fbIntegrScaling = iIntegrScaling; 00174 if (fdIntegrEps<=0.) 00175 fiMaxIntegrDepth = fiMinIntegrDepth; 00176 } 00177 00178 void TGo4FitModel::RemoveAllPars() { 00179 if (fxAllPars) { delete fxAllPars; fxAllPars=0; } 00180 if (fxAllParsValues) { delete fxAllParsValues; fxAllParsValues=0; } 00181 } 00182 00183 Bool_t TGo4FitModel::Initialize(Int_t UseBuffers) { 00184 Bool_t use = ((UseBuffers<0) && GetUseBuffers()) || (UseBuffers>0); 00185 for(Int_t n=0;n<NumAssigments();n++) { 00186 TGo4FitAssignment* ass = GetAssigment(n); 00187 if (ass->fxData==0) { 00188 std::cout << "Data " << ass->GetName() << " not assigned to model " << GetName() << std::endl; 00189 continue; 00190 } 00191 00192 if (use && ass->fxData->BuffersAllocated()) { 00193 if (IsAnyRangeLimits()) { 00194 ass->fxModelMask = new Char_t[ass->fxData->GetBinsSize()]; 00195 ass->fxData->ApplyRangesForModelMask(this,ass->fxModelMask); 00196 } 00197 00198 ass->fxModelBins = new Double_t[ass->fxData->GetBinsSize()]; 00199 } 00200 } 00201 00202 SetNeedToRebuild(); 00203 00204 if (use) { 00205 RemoveAllPars(); 00206 00207 fxAllPars = new TGo4FitParsList(kFALSE); 00208 CollectParsTo(*fxAllPars); 00209 00210 fxAllParsValues = new TArrayD(fxAllPars->NumPars()); 00211 fxAllPars->GetParsValues(fxAllParsValues->GetArray()); 00212 } 00213 00214 return kTRUE; 00215 } 00216 00217 void TGo4FitModel::Finalize() { 00218 RemoveAllPars(); 00219 fxCurrentPars.Set(0); 00220 for(Int_t n=0;n<NumAssigments();n++) { 00221 TGo4FitAssignment* ass = GetAssigment(n); 00222 if (ass->fxModelMask) { delete[] ass->fxModelMask; ass->fxModelMask = 0; } 00223 if (ass->fxModelBins) { delete[] ass->fxModelBins; ass->fxModelBins = 0; } 00224 } 00225 } 00226 00227 Bool_t TGo4FitModel::BuffersAllocated() const { 00228 Bool_t res = (NumAssigments()>0); 00229 for(Int_t n=0;n<NumAssigments();n++) 00230 res = res && (GetAssigment(n)->fxModelBins!=0); 00231 return res; 00232 } 00233 00234 Double_t* TGo4FitModel::GetModelBins(const char* DataName) const { 00235 TGo4FitAssignment* ass = FindAssigment(DataName); 00236 return ass != 0 ? ass->fxModelBins : 0; 00237 } 00238 00239 Double_t TGo4FitModel::EvaluateAndIntegrate(Int_t NumScales, const Double_t* Scales, const Double_t* Widths) { 00240 if ((NumScales<1) || (Scales==0)) return 0.; 00241 00242 if ( (fiMinIntegrDepth>0) && (fiMaxIntegrDepth>0) && (Widths!=0)) { 00243 TArrayI IntegrIndexes(NumScales); 00244 TArrayD ScaleValues(NumScales, Scales); 00245 TArrayD WidthValues(NumScales, Widths); 00246 TArrayD dScaleValues(NumScales); 00247 00248 Int_t* dindx = IntegrIndexes.GetArray(); 00249 Double_t* vector = ScaleValues.GetArray(); 00250 Double_t* width = WidthValues.GetArray(); 00251 Double_t* dvector = dScaleValues.GetArray(); 00252 00253 Int_t n=0; 00254 00255 // set value and range for current point of model 00256 for(n=0;n<NumScales;n++) 00257 vector[n] -= 0.5*width[n]; 00258 00259 // include left bound of interval 00260 Double_t TotalSum = EvalN(vector); 00261 Int_t TotalNumPnt = 1; 00262 00263 Int_t IntegrDepth = 0; 00264 Int_t NumStep = 1; 00265 Double_t Step = 1.0; 00266 Bool_t stopcondition = kFALSE; 00267 00268 // cycle over integration depths 00269 do { 00270 00271 IntegrIndexes.Reset(0); 00272 Double_t Sum = 0.; 00273 Int_t NumPnt = 0; 00274 00275 // run over all integration points 00276 do { 00277 for(n=0;n<NumScales;n++) 00278 dvector[n] = vector[n]+Step*width[n]*(dindx[n]+0.5); 00279 00280 Sum+=EvalN(dvector); 00281 NumPnt++; 00282 00283 n=0; 00284 do { 00285 dindx[n]++; 00286 if (dindx[n]<NumStep) break; 00287 dindx[n]=0; n++; 00288 } while (n<NumScales); 00289 } while (n<NumScales); 00290 00291 // check stop condition for integrations depth 00292 stopcondition = kFALSE; 00293 if (IntegrDepth>=fiMinIntegrDepth) { 00294 if (IntegrDepth>=fiMaxIntegrDepth) stopcondition = kTRUE; else { 00295 Double_t v1 = TotalSum/TotalNumPnt; 00296 Double_t v2 = Sum/NumPnt; 00297 Double_t v = TMath::Abs(v1-v2); 00298 if (!fbAbsoluteEps) { 00299 if ((v1!=0) || (v2!=0)) v=v/(TMath::Abs(v1)+TMath::Abs(v2)); 00300 else v=0.; 00301 } 00302 if (v<=fdIntegrEps) stopcondition = kTRUE; 00303 } 00304 } 00305 TotalSum+=Sum; 00306 TotalNumPnt+=NumPnt; 00307 IntegrDepth++; 00308 NumStep*=2; 00309 Step = Step/2.; 00310 } while (!stopcondition); 00311 00312 Double_t value = 0; 00313 if (TotalNumPnt>0) value = TotalSum / TotalNumPnt; 00314 if (fbIntegrScaling) 00315 for(n=0;n<NumScales;n++) value*=width[n]; 00316 return value; 00317 } else return EvalN(Scales); 00318 } 00319 00320 Double_t TGo4FitModel::EvaluateAtPoint(TGo4FitData* data, Int_t nbin, Bool_t UseRanges) { 00321 if (data==0) return 0.; 00322 const Double_t* scales = data->GetScaleValues(nbin); 00323 Int_t scalessize = data->GetScalesSize(); 00324 if (UseRanges && !CheckRangeConditions(scales, scalessize)) return 0.; 00325 return EvaluateAndIntegrate(scalessize, scales, data->GetWidthValues(nbin)); 00326 } 00327 00328 Double_t TGo4FitModel::EvaluateAtPoint(TGo4FitDataIter* iter, Bool_t UseRanges) { 00329 if (iter==0) return 0; 00330 const Double_t* scales = iter->Scales(); 00331 Int_t scalessize = iter->ScalesSize(); 00332 if (UseRanges && !CheckRangeConditions(scales, scalessize)) return 0.; 00333 return EvaluateAndIntegrate(scalessize, scales, iter->Widths()); 00334 } 00335 00336 void TGo4FitModel::RebuildShape(Bool_t ForceBuild) { 00337 if ((fxAllPars==0) || (fxAllParsValues==0)) return; 00338 00339 Bool_t fill = ForceBuild || fbNeedToRebuild; 00340 if (!fill) 00341 for(Int_t n=0;n<fxAllPars->NumPars();n++) { 00342 if (n==GetAmplIndex()) continue; 00343 if ((*fxAllParsValues)[n] != fxAllPars->GetPar(n)->GetValue()) { 00344 fill = kTRUE; 00345 break; 00346 } 00347 } 00348 00349 if(fill) 00350 for(Int_t n=0;n<NumAssigments();n++) { 00351 TGo4FitAssignment* ass = GetAssigment(n); 00352 00353 if ((ass->fxData==0) || (ass->fxModelBins==0)) continue; 00354 00355 TGo4FitData* data = ass->fxData; 00356 00357 if (!BeforeEval(data->GetScalesSize())) continue; 00358 00359 Int_t size = data->GetBinsSize(); 00360 Double_t* bins = ass->fxModelBins; 00361 Char_t* mask = ass->fxModelMask; 00362 Double_t ratio = ass->RatioValue(); 00363 00364 for(Int_t nbin=0;nbin<size;nbin++) { 00365 if ((mask!=0) && (mask[nbin]==0)) continue; 00366 bins[nbin] = ratio * EvaluateAtPoint(data, nbin, kFALSE); 00367 } 00368 00369 AfterEval(); 00370 } 00371 00372 00373 for(Int_t n=0;n<fxAllPars->NumPars();n++) 00374 (*fxAllParsValues)[n] = fxAllPars->GetPar(n)->GetValue(); 00375 00376 fbNeedToRebuild = kFALSE; 00377 } 00378 00379 Bool_t TGo4FitModel::AddModelToDataResult(TGo4FitData* data) { 00380 if ((data==0) || (!data->BuffersAllocated())) return kFALSE; 00381 00382 Double_t* result = data->GetBinsResult(); 00383 if (result==0) return kFALSE; 00384 Int_t size = data->GetBinsSize(); 00385 00386 Double_t* modelbins = GetModelBins(data->GetName()); 00387 if (modelbins) { 00388 Double_t ampl = GetAmplValue(); 00389 for (Int_t nbin=0;nbin<size;nbin++) 00390 result[nbin] += ampl * modelbins[nbin]; 00391 return kTRUE; 00392 } 00393 00394 if (!BeforeEval(data->GetScalesSize())) return kFALSE; 00395 00396 Double_t ampl = GetAmplValue() * GetRatioValueFor(data->GetName()); 00397 00398 for(Int_t nbin=0;nbin<size;nbin++) 00399 result[nbin] += ampl * EvaluateAtPoint(data, nbin); 00400 00401 AfterEval(); 00402 00403 return kTRUE; 00404 } 00405 00406 Bool_t TGo4FitModel::BeforeEval(Int_t ndim) { 00407 fxCurrentPars.Set(NumPars()); 00408 GetParsValues(fxCurrentPars.GetArray()); 00409 fxCurrentParsArray = fxCurrentPars.GetArray(); 00410 if (GetAmplPar()!=0) fxCurrentParsArray++; 00411 return kTRUE; 00412 } 00413 00414 Double_t TGo4FitModel::EvalN(const Double_t* v) { 00415 return UserFunction((Double_t*)v, fxCurrentParsArray); 00416 } 00417 00418 Double_t TGo4FitModel::Evaluate(Double_t x) { 00419 Double_t zn = 0.; 00420 if (BeforeEval(1)) { 00421 zn = GetAmplValue() * EvalN(&x); 00422 AfterEval(); 00423 } 00424 return zn; 00425 } 00426 00427 Double_t TGo4FitModel::Evaluate(Double_t x, Double_t y) { 00428 Double_t zn = 0.; 00429 if (BeforeEval(2)) { 00430 Double_t vector[2] = {x,y}; 00431 zn = GetAmplValue() * EvalN(vector); 00432 AfterEval(); 00433 } 00434 return zn; 00435 } 00436 00437 Double_t TGo4FitModel::Evaluate(Double_t x, Double_t y, Double_t z) { 00438 Double_t zn = 0.; 00439 if (BeforeEval(3)) { 00440 Double_t vector[3] = {x,y,z}; 00441 zn = GetAmplValue() * EvalN(vector); 00442 AfterEval(); 00443 } 00444 return zn;} 00445 00446 Double_t TGo4FitModel::Evaluate(Double_t* v, Int_t ndim) { 00447 Double_t zn = 0.; 00448 if (BeforeEval(ndim)) { 00449 zn = GetAmplValue() * EvalN(v); 00450 AfterEval(); 00451 } 00452 return zn; 00453 } 00454 00455 Double_t TGo4FitModel::Integral() { 00456 return 0; 00457 } 00458 00459 00460 const Int_t* TGo4FitModel::GetDataFullIndex(TGo4FitData* data, Int_t nbin) { 00461 return data==0 ? 0 : data->GetFullIndex(nbin); 00462 } 00463 00464 Int_t TGo4FitModel::GetDataIndexesSize(TGo4FitData* data) { 00465 return data==0 ? 0 : data->GetIndexesSize(); 00466 } 00467 00468 TGo4FitAssignment* TGo4FitModel::FindAssigment(const char* DataName) const { 00469 for (Int_t n=0;n<NumAssigments();n++) 00470 if (strcmp(DataName, GetAssigment(n)->GetName()) == 0) 00471 return GetAssigment(n); 00472 return 0; 00473 } 00474 00475 TString TGo4FitModel::GetRatioName(Int_t n) 00476 { 00477 TString res; 00478 res.Form("Ratio%d",n); 00479 return res; 00480 } 00481 00482 Double_t TGo4FitModel::GetRatioValueFor(const char* DataName) { 00483 TGo4FitAssignment* ass = FindAssigment(DataName); 00484 return (ass==0) ? 1. : ass->RatioValue(); 00485 } 00486 00487 void TGo4FitModel::Print(Option_t* option) const { 00488 TGo4FitComponent::Print(option); 00489 std::cout << " Assigned to: "; 00490 fxAssigments.Print(option); 00491 std::cout << std::endl; 00492 if ( (fiMinIntegrDepth>0) && (fiMaxIntegrDepth>0) ) { 00493 std::cout << " Integration property: depths from " << fiMinIntegrDepth << " to " << fiMaxIntegrDepth; 00494 if (fbAbsoluteEps) std::cout << " absolute"; else std::cout << " relative"; 00495 std::cout << " error " << fdIntegrEps << std::endl; 00496 } 00497 } 00498 00499 void TGo4FitModel::Streamer(TBuffer& b) { 00500 if (b.IsReading()) { 00501 00502 TGo4FitModel::Class()->ReadBuffer(b, this); 00503 00504 for (Int_t n=0;n<NumAssigments();n++) 00505 if (GetAssigment(n)->fxRatio) 00506 GetAssigment(n)->fxRatio->SetOwner(this); 00507 00508 } else { 00509 TGo4FitModel::Class()->WriteBuffer(b, this); 00510 } 00511 }
