TGo4FitData.cxx

Go to the documentation of this file.

// $Id$
//-----------------------------------------------------------------------
//       The GSI Online Offline Object Oriented (Go4) Project
//         Experiment Data Processing at EE department, GSI
//-----------------------------------------------------------------------
// Copyright (C) 2000- GSI Helmholtzzentrum fuer Schwerionenforschung GmbH
//                     Planckstr. 1, 64291 Darmstadt, Germany
// Contact:            http://go4.gsi.de
//-----------------------------------------------------------------------
// This software can be used under the license agreements as stated
// in Go4License.txt file which is part of the distribution.
//-----------------------------------------------------------------------

#include "TGo4FitData.h"

#include <iostream>

#include "TBuffer.h"
#include "TMath.h"
#include "TH1.h"
#include "TH2.h"
#include "TH3.h"
#include "TGraph.h"
#include "TClass.h"

#include "TGo4FitAxisTrans.h"

TGo4FitData::TGo4FitData()
   : TGo4FitComponent(), fiDataType(0), fbUseBinScale(kFALSE), fiTakeSigmasFrom(1), fdSigmaValue(1.),
     fdExcludeLessThen(0.), fxAxisTrans()
{
   ResetAllPoinetrs();
}

TGo4FitData::TGo4FitData(const char *iName, const char *iTitle, Int_t iDataType, Bool_t AddAmpl)
   : TGo4FitComponent(iName, iTitle), fiDataType(iDataType), fbUseBinScale(kFALSE), fiTakeSigmasFrom(1),
     fdSigmaValue(1.), fdExcludeLessThen(0.), fxAxisTrans()
{

   ResetAllPoinetrs();

   if (AddAmpl)
      NewAmplitude("Ampl", 1.0, kTRUE);

   fxAxisTrans.SetOwner(kTRUE);
}

TGo4FitData::~TGo4FitData()
{
   ReleaseAllPointers();
}

Bool_t TGo4FitData::SetNumberOfTransSlots(Int_t nslots)
{
   Int_t oldnum = GetNumberOfTransSlots();
   if ((nslots < 0) || (nslots == oldnum))
      return kFALSE;

   if (oldnum < nslots)
      for (Int_t n = oldnum; n < nslots; n++) {
         TString name("Trans");
         name += n;
         fxAxisTrans.Add(new TGo4FitSlot(name.Data(), "Axis transformation", this, TGo4FitAxisTrans::Class(), kFALSE));
      }
   else
      for (Int_t n = oldnum; n > nslots; n--) {
         TObject *slot = fxAxisTrans.Last();
         fxAxisTrans.Remove(slot);
         fxAxisTrans.Compress();
         delete slot;
      }
   SetUpdateSlotList();
   return kTRUE;
}

TGo4FitSlot *TGo4FitData::GetAxisTransSlot(Int_t nslot) const
{
   return (nslot >= 0) && (nslot <= fxAxisTrans.GetLast()) ? dynamic_cast<TGo4FitSlot *>(fxAxisTrans[nslot]) : nullptr;
}

TGo4FitAxisTrans *TGo4FitData::GetAxisTrans(Int_t nslot) const
{
   TGo4FitSlot *slot = GetAxisTransSlot(nslot);
   return !slot ? nullptr : dynamic_cast<TGo4FitAxisTrans *>(slot->GetObject());
}

void TGo4FitData::SetAxisTrans(Int_t nslot, TGo4FitAxisTrans *Trans, Bool_t TransOwned)
{
   if (nslot < 0)
      return;
   if (nslot >= GetNumberOfTransSlots())
      SetNumberOfTransSlots(nslot + 1);
   ((TGo4FitSlot *)(fxAxisTrans[nslot]))->SetObject(Trans, TransOwned);
}

void TGo4FitData::AddAxisTrans(TGo4FitAxisTrans *Trans, Bool_t TransOwned)
{
   Int_t nslot = GetNumberOfTransSlots();
   SetNumberOfTransSlots(nslot + 1);
   ((TGo4FitSlot *)(fxAxisTrans[nslot]))->SetObject(Trans, TransOwned);
}

void TGo4FitData::SetAxisTransNeeded(Int_t nslot, Bool_t iNeeded)
{
   if (iNeeded && (nslot >= GetNumberOfTransSlots()))
      SetNumberOfTransSlots(nslot + 1);
   if ((nslot >= 0) && (nslot < GetNumberOfTransSlots()))
      ((TGo4FitSlot *)(fxAxisTrans[nslot]))->SetNeeded(iNeeded);
}

Bool_t TGo4FitData::IsAnyDataTransform() const
{
   if (GetUseBinScale() || (GetExcludeLessThen() > 0))
      return kTRUE;
   for (Int_t n = 0; n < GetNumberOfTransSlots(); n++)
      if (GetAxisTrans(n))
         return kTRUE;
   return kFALSE;
}

TObject *TGo4FitData::CreateDrawObject(const char *ObjName)
{
   auto iter = MakeIter();
   if (!iter)
      return nullptr;
   TObject *obj = iter->CreateDrawObject(ObjName);
   return obj;
}

Bool_t TGo4FitData::Initialize(Int_t UseBuffers)
{
   auto iter = MakeIter();
   if (!iter)
      return kFALSE;

   fiBinsSize = iter->CountPoints(kTRUE);

   fiIndexesSize = iter->IndexesSize();
   fiScalesSize = iter->ScalesSize();

   Bool_t use = ((UseBuffers < 0) && GetUseBuffers()) || (UseBuffers > 0);

   if (use)
      for (Int_t n = 0; n < GetNumberOfTransSlots(); n++) {
         TGo4FitAxisTrans *trans = GetAxisTrans(n);
         if (trans && !trans->IsAllParsFixed()) {
            use = kFALSE;
            break;
         }
      }

   if (use) {

      fxValues = new Double_t[fiBinsSize];
      fxStandDev = new Double_t[fiBinsSize];
      fxBinsResult = new Double_t[fiBinsSize];

      if (iter->HasIndexes())
         fxFullIndex = new Int_t[fiBinsSize * fiIndexesSize];
      fxFullScale = new Double_t[fiBinsSize * fiScalesSize];
      if (iter->HasWidths())
         fxFullWidth = new Double_t[fiBinsSize * fiScalesSize];

      Int_t nbin = 0;
      if (iter->Reset())
         do {

            fxValues[nbin] = iter->Value();
            fxStandDev[nbin] = iter->StandardDeviation();

            if (fxFullIndex)
               for (Int_t n = 0; n < fiIndexesSize; n++)
                  fxFullIndex[nbin * fiIndexesSize + n] = iter->Indexes()[n];

            if (fxFullScale)
               for (Int_t naxis = 0; naxis < fiScalesSize; naxis++)
                  fxFullScale[nbin * fiScalesSize + naxis] = iter->Scales()[naxis];

            if (fxFullWidth && iter->HasWidths())
               for (Int_t naxis = 0; naxis < fiScalesSize; naxis++)
                  fxFullWidth[nbin * fiScalesSize + naxis] = iter->Widths()[naxis];

            nbin++;
         } while (iter->Next());
   }

   return kTRUE;
}

void TGo4FitData::Finalize()
{
   ReleaseAllPointers();
}

void TGo4FitData::ResetAllPoinetrs()
{
   fiBinsSize = 0;
   fiIndexesSize = 0;
   fiScalesSize = 0;

   fxValues = nullptr;
   fxStandDev = nullptr;
   fxBinsResult = nullptr;

   fxFullScale = nullptr;
   fxFullWidth = nullptr;

   fxFullIndex = nullptr;
}

void TGo4FitData::ReleaseAllPointers()
{
   if (fxValues)
      delete[] fxValues;

   if (fxStandDev)
      delete[] fxStandDev;

   if (fxBinsResult)
      delete[] fxBinsResult;

   if (fxFullIndex)
      delete[] fxFullIndex;

   if (fxFullScale)
      delete[] fxFullScale;

   if (fxFullWidth)
      delete[] fxFullWidth;

   ResetAllPoinetrs();
}

Bool_t TGo4FitData::DefineScaleMinMax(Int_t naxis, Double_t &min, Double_t &max)
{
   auto iter = MakeIter();
   if (!iter)
      return kFALSE;
   Bool_t res = kFALSE;
   if (iter->Reset(kFALSE) && (iter->ScalesSize() <= naxis)) {
      min = iter->Scales()[naxis];
      max = min;
      do {
         Double_t value = iter->Scales()[naxis];
         if (value < min)
            min = value;
         else if (value > max)
            max = value;
      } while (iter->Next(kFALSE));
      res = kTRUE;
   }

   return res;
}

Int_t TGo4FitData::DefineDimensions()
{
   auto iter = MakeIter();
   if (!iter)
      return 0;
   Int_t res = 0;
   if (iter->Reset(kFALSE))
      res = iter->IndexesSize();
   return res;
}

Int_t TGo4FitData::DefineBinsSize()
{
   auto iter = MakeIter();

   return iter ? iter->CountPoints(kTRUE) : 0;
}

const Double_t *TGo4FitData::GetScaleValues(Int_t nbin) const
{
   if (fxFullScale)
      return &(fxFullScale[nbin * GetScalesSize()]);

   return nullptr;
}

const Double_t *TGo4FitData::GetWidthValues(Int_t nbin) const
{
   if (fxFullWidth)
      return &(fxFullWidth[nbin * GetScalesSize()]);
   return nullptr;
}

const Int_t *TGo4FitData::GetFullIndex(Int_t nbin) const
{
   if (fxFullIndex)
      return &(fxFullIndex[nbin * GetIndexesSize()]);
   return nullptr;
}

Bool_t TGo4FitData::IsCompatibleData(TGo4FitData *data)
{
   if (!data)
      return kFALSE;
   auto iter = data->MakeIter();
   if (!iter)
      return kFALSE;

   Bool_t res = kFALSE;
   if (iter->Reset(kFALSE))
      res = (iter->IndexesSize() == GetIndexesSize()) && (GetIndexesSize() > 0);

   return res;
}

void TGo4FitData::ApplyRangesForModelMask(TGo4FitComponent *model, Char_t *ModelMask)
{
   if (!ModelMask)
      return;

   if (BuffersAllocated()) {
      for (Int_t nbin = 0; nbin < GetBinsSize(); nbin++) {
         const Double_t *values = GetScaleValues(nbin);

         Bool_t res = model->CheckRangeConditions(values, GetScalesSize());

         ModelMask[nbin] = res ? 1 : 0;
      }
   } else {
      auto iter = MakeIter();
      Int_t nbin = 0;
      if (iter->Reset())
         do {
            Bool_t res = model->CheckRangeConditions(iter->Scales(), iter->ScalesSize());
            ModelMask[nbin] = res ? 1 : 0;
            nbin++;
         } while (iter->Next());
   }
}

void TGo4FitData::FillSlotList(TSeqCollection *list)
{
   TGo4FitComponent::FillSlotList(list);
   for (Int_t n = 0; n <= fxAxisTrans.GetLast(); n++)
      list->Add(fxAxisTrans[n]);
}

void TGo4FitData::Print(Option_t *option) const
{
   TGo4FitComponent::Print(option);
   std::cout << "   Data type: ";
   switch (fiDataType) {
      case 1: std::cout << "histogram" << std::endl; break;
      case 2: std::cout << "graph" << std::endl; break;
      default: std::cout << fiDataType << std::endl;
   }
   std::cout << "   Use bin scale:  " << fbUseBinScale << std::endl;
   std::cout << "   Take sigmas from:  ";
   switch (GetSigmaSource()) {
      case 0: std::cout << "none" << std::endl; break;
      case 1: std::cout << "data" << std::endl; break;
      case 2: std::cout << "const value " << GetSigmaValue() << std::endl; break;
   }
   std::cout << "   Exclude bins less then: " << GetExcludeLessThen() << std::endl;
   std::cout << "   Axis transformation data: " << std::endl;
   fxAxisTrans.Print(option);
}

void TGo4FitData::Streamer(TBuffer &b)
{
   if (b.IsReading()) {

      TGo4FitData::Class()->ReadBuffer(b, this);

      for (Int_t n = 0; n <= fxAxisTrans.GetLast(); n++) {
         TGo4FitSlot *dc = (TGo4FitSlot *)fxAxisTrans[n];
         dc->SetDefaults(this, TGo4FitAxisTrans::Class());
      }

   } else {
      TGo4FitData::Class()->WriteBuffer(b, this);
   }
}

// *******************************************************************************************

TGo4FitDataIter::TGo4FitDataIter()
   : TObject(), fxIndexes(), fxScales(), fxWidths(), fdValue(0.), fdStandardDeviation(1.), fiNumPoint(0),
     fbReachEnd(kTRUE)
{
}

TGo4FitDataIter::~TGo4FitDataIter() {}

Bool_t TGo4FitDataIter::ReserveArrays(Int_t NumDimen, Int_t NumOwnAxis, Bool_t HasWidth)
{
   TGo4FitData *data = GetData();
   if (!data)
      return kFALSE;

   fxIndexes.Set(NumDimen);
   fxIndexes.Reset(0);

   Int_t size = 0;
   if (data->GetUseBinScale())
      size = NumDimen;
   else
      size = NumOwnAxis;

   if (size <= 0)
      return kFALSE;

   fxScales.Set(size);
   fxScales.Reset(0.);
   if (HasWidth) {
      fxWidths.Set(size);
      fxWidths.Reset(1.);
   } else
      fxWidths.Set(0);

   return kTRUE;
}

void TGo4FitDataIter::TransformScales(Double_t *scales)
{
   TGo4FitData *data = GetData();
   for (Int_t nslot = 0; nslot < data->GetNumberOfTransSlots(); nslot++) {
      TGo4FitAxisTrans *trans = data->GetAxisTrans(nslot);
      if (trans)
         trans->Transformation(scales, ScalesSize());
   }
}

Bool_t TGo4FitDataIter::ProduceScales(const Int_t *index, const Double_t *ownscales, const Double_t *ownwidths)
{
   TGo4FitData *data = GetData();
   if (!data)
      return kFALSE;

   if (data->GetUseBinScale() || !ownscales) {
      if (!index)
         return kFALSE;
      Double_t add = (data->GetDataType() == TGo4FitData::dtHistogram) ? .5 : 0.;
      for (Int_t n = 0; n < fxScales.GetSize(); n++)
         fxScales[n] = index[n] + add;
      fxWidths.Reset(1.);
   } else {
      for (Int_t n = 0; n < fxScales.GetSize(); n++)
         fxScales[n] = ownscales[n];
      if (ownwidths)
         for (Int_t n = 0; n < fxWidths.GetSize(); n++)
            fxWidths[n] = ownwidths[n];
   }

   if (data->GetNumberOfTransSlots() > 0) {
      if (fxWidths.GetSize() == ScalesSize()) {
         TArrayD arr1(ScalesSize()), arr2(ScalesSize());
         for (Int_t n = 0; n < ScalesSize(); n++) {
            arr1[n] = fxScales[n] - fxWidths[n] / 2.;
            arr2[n] = fxScales[n] + fxWidths[n] / 2.;
         }
         TransformScales(arr1.GetArray());
         TransformScales(arr2.GetArray());
         for (Int_t n = 0; n < ScalesSize(); n++)
            fxWidths[n] = TMath::Abs(arr2[n] - arr1[n]);
      }

      TransformScales(fxScales.GetArray());
   }

   return kTRUE;
}

Bool_t TGo4FitDataIter::NextIndex(TArrayI &Index, TArrayI &Limits)
{
   Int_t n = 0;
   while (n < Index.GetSize()) {
      Index[n]++;
      if (Index[n] < Limits[n])
         return kTRUE;
      Index[n] = 0;
      n++;
   }
   return kFALSE;
}

Bool_t TGo4FitDataIter::GetDeviation()
{
   TGo4FitData *data = GetData();
   if (!data)
      return kFALSE;
   if (data->GetSigmaSource() == 2) {
      fdStandardDeviation = data->GetSigmaValue() * data->GetSigmaValue();
      return kTRUE;
   } else
      return kFALSE;
}

Bool_t TGo4FitDataIter::CheckPointForRange()
{
   TGo4FitData *data = GetData();
   if (!data)
      return kFALSE;
   if (Value() < data->GetExcludeLessThen())
      return kFALSE;
   return data->CheckRangeConditions(Scales(), ScalesSize());
}

Bool_t TGo4FitDataIter::Reset(Bool_t UseRanges)
{
   fbReachEnd = kTRUE;

   if (!StartReset())
      return kFALSE;

   fiNumPoint = 0;

   if (!ReadCurrentPoint())
      return kFALSE;
   if (!UseRanges) {
      fbReachEnd = kFALSE;
      return kTRUE;
   }

   while (!CheckPointForRange()) {
      if (!ShiftToNextPoint())
         return kFALSE;
      if (!ReadCurrentPoint())
         return kFALSE;
   }

   fbReachEnd = kFALSE;
   return kTRUE;
}

Bool_t TGo4FitDataIter::Next(Bool_t UseRanges)
{
   fiNumPoint++;

   if (fbReachEnd || !GetData()) {
      fbReachEnd = kTRUE;
      return kFALSE;
   }

   do {
      if (!ShiftToNextPoint()) {
         fbReachEnd = kTRUE;
         return kFALSE;
      }

      if (!ReadCurrentPoint()) {
         fbReachEnd = kTRUE;
         return kFALSE;
      }

      if (!UseRanges)
         return kTRUE;

   } while (!CheckPointForRange());

   return kTRUE;
}

Double_t TGo4FitDataIter::xWidths() const
{
   double res = 1.;
   if (HasWidths())
      for (int n = 0; n < fxWidths.GetSize(); n++)
         res = res * fxWidths[n];
   return res;
}

Int_t TGo4FitDataIter::CountPoints(Bool_t UseRanges)
{
   if (!Reset(UseRanges))
      return 0;
   Int_t cnt = 0;
   do {
      cnt += 1;
   } while (Next(UseRanges));
   return cnt;
}

Bool_t TGo4FitDataIter::DefineIndexesLimits(TArrayI &Limits)
{
   if (!Reset(kFALSE))
      return kFALSE;
   if (IndexesSize() <= 0)
      return kFALSE;
   Limits.Set(IndexesSize());
   Limits.Reset(0);
   do {
      for (Int_t n = 0; n < IndexesSize(); n++)
         if (Indexes()[n] > Limits[n])
            Limits[n] = Indexes()[n];
   } while (Next(kFALSE));
   return kTRUE;
}

TH1 *TGo4FitDataIter::CreateHistogram(const char *HistoName, Bool_t UseRanges, Bool_t SetBins)
{
   TArrayI Limits;
   if (!DefineIndexesLimits(Limits))
      return nullptr;
   if (!HasIndexes() || (IndexesSize() != ScalesSize()) || !HasWidths())
      return nullptr;

   Int_t NumDim = IndexesSize();
   if (NumDim > 3)
      NumDim = 3;

   Double_t *dummy = nullptr;
   TH1 *histo = nullptr;
   switch (NumDim) {
   case 1: histo = new TH1D(HistoName, "result", Limits[0] + 1, dummy); break;
   case 2: histo = new TH2D(HistoName, "result", Limits[0] + 1, dummy, Limits[1] + 1, dummy); break;
   case 3:
      histo = new TH3D(HistoName, "result", Limits[0] + 1, dummy, Limits[1] + 1, dummy, Limits[2] + 1, dummy);
      break;
   default: return nullptr;
   }

   histo->SetDirectory(nullptr);

   Double_t *Axises[3];
   for (Int_t n = 0; n < NumDim; n++)
      Axises[n] = new Double_t[Limits[n] + 2];

   Double_t ampl = GetData()->GetAmplValue();

   if (Reset(UseRanges))
      do {
         if (SetBins)
            switch (NumDim) {
            case 1: histo->SetBinContent(Indexes()[0] + 1, ampl * Value()); break;
            case 2: histo->SetBinContent(Indexes()[0] + 1, Indexes()[1] + 1, ampl * Value()); break;
            case 3: histo->SetBinContent(Indexes()[0] + 1, Indexes()[1] + 1, Indexes()[2] + 1, ampl * Value()); break;
            }
         for (Int_t n = 0; n < NumDim; n++) {
            Int_t indx = Indexes()[n];
            Axises[n][indx] = Scales()[n] - Widths()[n] / 2.;
            Axises[n][indx + 1] = Scales()[n] + Widths()[n] / 2.;
         }
      } while (Next(UseRanges));

   histo->GetXaxis()->Set(Limits[0] + 1, Axises[0]);
   if (NumDim > 1)
      histo->GetYaxis()->Set(Limits[1] + 1, Axises[1]);
   if (NumDim > 2)
      histo->GetZaxis()->Set(Limits[2] + 1, Axises[2]);

   for (Int_t n = 0; n < NumDim; n++)
      delete[] Axises[n];

   return histo;
}

TGraph *TGo4FitDataIter::CreateGraph(const char *GraphName, Bool_t UseRanges, Bool_t SetBins)
{
   Int_t NumPoints = CountPoints(UseRanges);
   if ((NumPoints <= 0) || (ScalesSize() < 1))
      return nullptr;

   TGraph *gr = new TGraph(NumPoints);
   gr->SetName(GraphName);
   if (Reset(UseRanges))
      do {
         (gr->GetX())[Point()] = x();
         if (SetBins)
            (gr->GetY())[Point()] = GetData()->GetAmplValue() * Value();
      } while (Next(UseRanges));

   return gr;
}

TObject *TGo4FitDataIter::CreateDrawObject(const char *ObjName)
{
   if (!Reset(kFALSE))
      return nullptr;

   if (HasIndexes() && (IndexesSize() == ScalesSize()) && HasWidths())
      return CreateHistogram(ObjName, kFALSE, kTRUE);

   return CreateGraph(ObjName, kFALSE, kTRUE);
}