Go4-v3.04-1: Go4Fit/TGo4FitAmplEstimation.cxx Source File

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