GSI Object Oriented Online Offline (Go4): Go4Fit/TGo4FitAmplEstimation.cxx Source File

Go to the documentation of this file.
00001 // $Id: TGo4FitAmplEstimation.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 "TGo4FitAmplEstimation.h"
00015 
00016 #include "Riostream.h"
00017 
00018 #include "TMath.h"
00019 #include "TMatrixD.h"
00020 #include "TVectorD.h"
00021 #include "TArrayC.h"
00022 #include "TArrayD.h"
00023 
00024 #include "TGo4Fitter.h"
00025 #include "TGo4FitModel.h"
00026 #include "TGo4FitData.h"
00027 #include "TGo4FitParameter.h"
00028 
00029 TGo4FitAmplEstimation::TGo4FitAmplEstimation() : TGo4FitterAction(), fiNumIters(1) {
00030 }
00031 
00032 TGo4FitAmplEstimation::TGo4FitAmplEstimation(const char* Name, Int_t NumIters) :
00033   TGo4FitterAction(Name,"Estimate amplitude action"), fiNumIters(NumIters) {
00034 }
00035 
00036 TGo4FitAmplEstimation::~TGo4FitAmplEstimation() {
00037 }
00038 
00039 void TGo4FitAmplEstimation::DoAction(TGo4FitterAbstract* Fitter) {
00040    TGo4Fitter* f = dynamic_cast<TGo4Fitter*> (Fitter);
00041    if(f==0) return;
00042    if (!CalculateWithBuffers(f))
00043      CalculateWithIterators(f);
00044 }
00045 
00046 Double_t TGo4FitAmplEstimation::PointWeight(Int_t niter, Int_t FFtype, Double_t value, Double_t modelvalue, Double_t standdev) {
00047    switch (FFtype) {
00048       case TGo4Fitter::ff_least_squares: return 1.;
00049       case TGo4Fitter::ff_chi_square: if (standdev<=0.) return 0.; else return 1./standdev;
00050       case TGo4Fitter::ff_chi_Pearson:
00051       case TGo4Fitter::ff_ML_Poisson:
00052         if (niter==0) if (value<0.) return 0.; else return 1./(value+1.);
00053                  else if (modelvalue<=0.) return 0.; else return 1./modelvalue;
00054       case TGo4Fitter::ff_chi_Neyman: if (value<1.) return 1.; else return 1./value;
00055       case TGo4Fitter::ff_chi_gamma: if (value<0.) return 0.; else return 1./(value+1.);
00056       case TGo4Fitter::ff_user: return 1.;
00057       default: return 1.;
00058    }
00059 }
00060 
00061 
00062 Bool_t TGo4FitAmplEstimation::CalculateWithBuffers(TGo4Fitter* fitter) {
00063    if (!fitter->IsInitialized()) return kFALSE;
00064 
00065    for(Int_t n=0;n<fitter->GetNumData();n++)
00066      if (!fitter->DataBuffersAllocated(fitter->GetData(n))) return kFALSE;
00067 
00068    for(Int_t n=0;n<fitter->GetNumModel();n++)
00069      if (!fitter->ModelBuffersAllocated(fitter->GetModel(n))) return kFALSE;
00070 
00071    TObjArray comps, fixedcomps;
00072 
00073    for(Int_t n=0;n<fitter->GetNumModel();n++) {
00074      TGo4FitModel *fModel = fitter->GetModel(n);
00075 
00076      Bool_t assigned = kFALSE;
00077      for (Int_t na=0;na<fModel->NumAssigments();na++)
00078        if (fitter->FindData(fModel->AssignmentName(na))) {
00079           assigned = kTRUE;
00080           break;
00081        }
00082      if (!assigned) continue;
00083 
00084      if ( (fModel->GetAmplPar()!=0) &&
00085           (!fModel->GetAmplPar()->GetFixed()) ) comps.Add(fModel);
00086                                            else fixedcomps.Add(fModel);
00087    }
00088 
00089    Int_t ncomp = comps.GetLast()+1;
00090    if (ncomp<=0) return kFALSE;
00091 
00092    TMatrixD matr(ncomp,ncomp);
00093    TVectorD vector(ncomp);
00094 
00095    Double_t** DopData = new Double_t* [fitter->GetNumData()];
00096    for(Int_t n=0;n<fitter->GetNumData();n++)
00097       DopData[n] = new Double_t[fitter->GetDataBinsSize(fitter->GetData(n))];
00098 
00099    Double_t** Weights = new Double_t* [fitter->GetNumData()];
00100    for(Int_t n=0;n<fitter->GetNumData();n++)
00101       Weights[n] = new Double_t[fitter->GetDataBinsSize(fitter->GetData(n))];
00102 
00103 
00104    Int_t niter = 0;
00105    Int_t FFtype = fitter->GetFitFunctionType();
00106    Bool_t changed = kFALSE;
00107 
00108    do {
00109 
00110        fitter->RebuildAll();
00111 
00112        for(Int_t ndata=0;ndata<fitter->GetNumData();ndata++) {
00113            TGo4FitData* data = fitter->GetData(ndata);
00114 
00115            Int_t size = fitter->GetDataBinsSize(data);
00116            Double_t ampl = data->GetAmplValue();
00117 
00118            Double_t* values = fitter->GetDataBinsValues(data);
00119            Double_t* result = fitter->GetDataBinsResult(data);
00120            Double_t* devs = fitter->GetDataBinsDevs(data);
00121            Double_t* weight = Weights[ndata];
00122 
00123           for (Int_t nbin=0;nbin<size;nbin++)
00124              weight[nbin] = PointWeight(niter,FFtype, ampl*values[nbin], result[nbin], ampl*ampl*devs[nbin]);
00125        }
00126 
00127        for(Int_t n1=0;n1<ncomp;n1++)
00128           for(Int_t n2=n1;n2<ncomp;n2++) {
00129              Double_t fSum = 0.;
00130              for(Int_t ndata=0;ndata<fitter->GetNumData();ndata++) {
00131                 Double_t* bins1 = fitter->GetModelBinsValues((TGo4FitModel*) comps[n1], fitter->GetDataName(ndata));
00132                 Double_t* bins2 = fitter->GetModelBinsValues((TGo4FitModel*) comps[n2], fitter->GetDataName(ndata));
00133                 if ((bins1!=0) && (bins2!=0)) {
00134                    Int_t size = fitter->GetDataBinsSize(fitter->GetData(ndata));
00135                    Double_t* weight = Weights[ndata];
00136                    for(Int_t nbin=0;nbin<size;nbin++)
00137                      fSum+=bins1[nbin]*bins2[nbin]*weight[nbin];
00138                 }
00139              }
00140              matr(n1,n2) = fSum;
00141              matr(n2,n1) = fSum;
00142           }
00143 
00144         // caluclating difference between data and fixed components
00145         for(Int_t ndata=0;ndata<fitter->GetNumData();ndata++) {
00146            TGo4FitData* data = fitter->GetData(ndata);
00147            Int_t size = fitter->GetDataBinsSize(data);
00148            Double_t dampl = data->GetAmplValue();
00149            Double_t* values = fitter->GetDataBinsValues(data);
00150            Double_t* bins = DopData[ndata];
00151 
00152            for(Int_t nbin=0;nbin<size;nbin++)
00153              bins[nbin] = dampl*values[nbin];
00154 
00155            for(Int_t n=0;n<=fixedcomps.GetLast();n++) {
00156               TGo4FitModel *model = (TGo4FitModel*) fixedcomps[n];
00157               Double_t* mbins = fitter->GetModelBinsValues(model,data->GetName());
00158               if (mbins!=0) {
00159                  Double_t mampl = model->GetAmplValue();
00160                  for(Int_t nbin=0;nbin<size;nbin++)
00161                    bins[nbin] -= mampl*mbins[nbin];
00162               }
00163            }
00164         }
00165 
00166         // filling vector
00167 
00168         for(Int_t n=0;n<ncomp;n++) {
00169            Double_t fSum = 0.;
00170            for(Int_t ndata=0;ndata<fitter->GetNumData();ndata++) {
00171                Double_t* bins1 = fitter->GetModelBinsValues((TGo4FitModel*) comps[n], fitter->GetDataName(ndata));
00172                if (bins1==0) continue;
00173                Int_t size = fitter->GetDataBinsSize(fitter->GetData(ndata));
00174                Double_t* bins = DopData[ndata];
00175                Double_t* weight = Weights[ndata];
00176                for(Int_t nbin=0;nbin<size;nbin++)
00177                   fSum+=bins[nbin]*bins1[nbin]*weight[nbin];
00178 
00179            }
00180            vector(n) = fSum;
00181         }
00182 
00183         if (matr.Determinant()==0.) {
00184            std::cerr << "  Amplitude estimation algorithm failed. " << std::endl;
00185            std::cerr << "  Determinant of matrice == 0.; This may be due to equiavalent model components or zero model at all" << std::endl;
00186            break;
00187         }
00188 
00189         matr.Invert();
00190         vector*=matr;
00191 
00192         for(Int_t n=0;n<ncomp;n++) {
00193            TGo4FitModel *model = (TGo4FitModel*) comps[n];
00194            model->SetAmplValue(vector(n));
00195            if (matr(n,n)>0.)
00196              model->SetAmplError(TMath::Sqrt(matr(n,n)));
00197         }
00198 
00199         changed = kFALSE;
00200         niter++;
00201         if (niter<fiNumIters)
00202           if ((FFtype==TGo4Fitter::ff_chi_Pearson) ||
00203               (FFtype==TGo4Fitter::ff_ML_Poisson)) changed=kTRUE;
00204 
00205     } while (changed);
00206 
00207    for(Int_t n=0;n<fitter->GetNumData();n++) delete[] Weights[n];
00208    delete[] Weights;
00209 
00210    for(Int_t n=0;n<fitter->GetNumData();n++) delete[] DopData[n];
00211    delete[] DopData;
00212 
00213    return kTRUE;
00214 }
00215 
00216 
00217 Bool_t TGo4FitAmplEstimation::CalculateWithIterators(TGo4Fitter* fitter) {
00218    if (fitter==0) return kFALSE;
00219 
00220    Int_t nmodel = fitter->GetNumModel();
00221    TArrayI refindex(fitter->GetNumModel()); refindex.Reset(-1);
00222    TArrayI fixedindex(fitter->GetNumModel()); fixedindex.Reset(-1);
00223 
00224    Int_t ncomp = 0;
00225    Int_t nfixed = 0;
00226 
00227    for(Int_t n=0;n<nmodel;n++) {
00228      TGo4FitModel *model = fitter->GetModel(n);
00229      Bool_t assigned = kFALSE;
00230      for (Int_t na=0;na<model->NumAssigments();na++)
00231        if (fitter->FindData(model->AssignmentName(na))) {
00232           assigned = kTRUE;
00233           break;
00234        }
00235      if (!assigned) continue;
00236 
00237      if ( (model->GetAmplPar()!=0) &&
00238           (!model->GetAmplPar()->GetFixed()) ) refindex[ncomp++] = n;
00239                                           else fixedindex[nfixed++] = n;
00240    }
00241 
00242    if (ncomp<=0) return kFALSE;
00243 
00244    TMatrixD matr(ncomp,ncomp);
00245    TVectorD vector(ncomp);
00246 
00247    TArrayC Usage(nmodel);
00248    TArrayD ModelValues(nmodel);
00249    TArrayD ModelAmpls(nmodel);
00250 
00251    Int_t niter = 0;
00252    Int_t FFtype = fitter->GetFitFunctionType();
00253    Bool_t changed = kFALSE;
00254 
00255    do {
00256        matr = 0.; vector = 0.;
00257 
00258        for(Int_t ndata=0;ndata<fitter->GetNumData();ndata++) {
00259           TGo4FitData* data = fitter->GetData(ndata);
00260           TGo4FitDataIter* iter = data->MakeIter();
00261           if (iter==0) return kFALSE;
00262 
00263           Double_t dampl = data->GetAmplValue();
00264 
00265           if (iter->Reset()) {
00266              Usage.Reset(0); ModelAmpls.Reset(0.); ModelValues.Reset(0.);
00267              for(Int_t n=0;n<nmodel;n++) {
00268                  TGo4FitModel* model = fitter->GetModel(n);
00269                  if (model->IsAssignTo(data->GetName())) {
00270                    Usage[n] = 1;
00271                    ModelAmpls[n] = model->GetAmplValue() * model->GetRatioValueFor(data->GetName());
00272                    model->BeforeEval(iter->IndexesSize());
00273                  }
00274              }
00275 
00276              do {
00277                Double_t result = 0.;
00278                for (Int_t n=0; n<nmodel;n++)
00279                  if (Usage[n]) {
00280                    ModelValues[n] = fitter->GetModel(n)->EvaluateAtPoint(iter);
00281                    result += ModelAmpls[n] * ModelValues[n];
00282                  }
00283                Double_t weight = PointWeight(niter, FFtype, dampl*iter->Value(), result, dampl*dampl*iter->StandardDeviation());
00284                if (weight<=0.) continue;
00285 
00286                for (Int_t i1=0;i1<ncomp;i1++)
00287                  for(Int_t i2=i1;i2<ncomp;i2++) {
00288                     Double_t zn = matr(i1,i2);
00289                     zn += ModelValues[refindex[i1]] * ModelValues[refindex[i2]] * weight;
00290                     matr(i1,i2) = zn;
00291                     if (i1!=i2) matr(i2,i1) = zn;
00292                  }
00293 
00294                Double_t dvalue = dampl*iter->Value();
00295                for(Int_t nf=0;nf<nfixed;nf++)
00296                  dvalue -= ModelAmpls[fixedindex[nf]] * ModelValues[fixedindex[nf]];
00297 
00298                for (Int_t i1=0;i1<ncomp;i1++) {
00299                  Double_t zn = vector(i1);
00300                  zn += ModelValues[refindex[i1]] * dvalue * weight;
00301                  vector(i1) = zn;
00302                }
00303 
00304              } while (iter->Next());
00305 
00306              for(Int_t n=0;n<nmodel;n++)
00307                if (Usage[n]) fitter->GetModel(n)->AfterEval();
00308           }
00309           delete iter;
00310        }
00311 
00312        if (matr.Determinant()==0.) {
00313            std::cerr << "  Amplitude estimation algorithm failed. " << std::endl;
00314            std::cerr << "  Determinant of matrice == 0.; This may be due to equiavalent model components or zero model at all" << std::endl;
00315            return kFALSE;
00316         }
00317 
00318         matr.Invert();
00319         vector*=matr;
00320 
00321         for(Int_t n=0;n<ncomp;n++) {
00322            TGo4FitModel *model = fitter->GetModel(refindex[n]);
00323            model->SetAmplValue(vector(n));
00324            if (matr(n,n)>0.)
00325              model->SetAmplError(TMath::Sqrt(matr(n,n)));
00326         }
00327 
00328         changed = kFALSE;
00329         niter++;
00330         if (niter<fiNumIters)
00331           if ((FFtype==TGo4Fitter::ff_chi_Pearson) ||
00332               (FFtype==TGo4Fitter::ff_ML_Poisson)) changed=kTRUE;
00333 
00334     } while (changed);
00335 
00336    return kTRUE;
00337 }
00338 
00339 void TGo4FitAmplEstimation::Print(Option_t* option) const
00340 {
00341    TGo4FitterAction::Print(option);
00342    std::cout << "  number of iterations " << fiNumIters << std::endl;
00343 }