/Go4Fit/TGo4FitAmplEstimation.cxx Source File

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