TXXXCalibPar.cxx

Go to the documentation of this file.

// $Id: TXXXCalibPar.cxx 478 2009-10-29 12:26:09Z linev $
//-----------------------------------------------------------------------
//       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 "TXXXCalibPar.h"

#include "TMath.h"
#include "TH1.h"

#include "TGo4Log.h"
#include "TGo4Fitter.h"
#include "TGo4FitModel.h"
#include "TGo4Analysis.h"

//***********************************************************
TXXXCalibPar::TXXXCalibPar() :
   TGo4Parameter(),
   fbRecalibrate(kFALSE),
   fbReadDatabase(kFALSE)
{
   fxDatabase = "calilines.txt";
   for (Int_t ord = 0; ord < __POLORDER__; ++ord)
      fdA[ord] = 0;
   for (Int_t ix = 0; ix < __LINESNUMBER__; ++ix) {
      fiLinesChannel[ix] = 0;
      ffLinesEnergy[ix] = 0;
      fxLinesNames[ix] = TString::Format("Defaultline-%d",ix);
   }
}
//***********************************************************
TXXXCalibPar::TXXXCalibPar(const char *name, TH1 *spectrum, TGraph *curve) :
   TGo4Parameter(name),
   fbRecalibrate(kFALSE),
   fbReadDatabase(kFALSE),
   fxCalibCurve(curve),
   fxCalibSpectrum(spectrum)
{
   // Set up fitters:
   fxLinesFinder = new TGo4Fitter("Linefinder", TGo4Fitter::ff_least_squares, kTRUE);
   fxCalibrator = new TGo4Fitter("Calibrator", TGo4Fitter::ff_least_squares, kTRUE);
   if (fxCalibSpectrum) {
      fxLinesFinder->AddH1(__DATANAME__, fxCalibSpectrum, kFALSE);
      fxSpectrumName = fxCalibSpectrum->GetName();
   } else {
      fxSpectrumName = "Please specify calibration spectrum";
   }
   if (fxCalibCurve) {
      fxCalibrator->AddGraph(__GRAPHNAME__, fxCalibCurve, kFALSE);
      fxGraphName = fxCalibCurve->GetName();
   } else {
      fxSpectrumName = "Please specify fit graph name";
   }
   fxCalibrator->AddPolynomX(__GRAPHNAME__, "A", __POLORDER__ - 1);
   // note that __POLORDER__ is number of polynom parameters here
   // i.e. true order of polynom +1
   for (Int_t i = 0; i < __POLORDER__; ++i) {
      fdA[i] = 1 / (i + 1);
      TString modname = TString::Format("A_%d",i);
      TGo4FitModel *mod = fxCalibrator->FindModel(modname.Data());
      if(mod) {
         // for the beginning, disable models beyond order 1:
         if(i>1) mod->ClearAssignmentTo(__GRAPHNAME__);
      } else
         TGo4Log::Error("could not find model %s", modname.Data());
   }

   for (Int_t ix = 0; ix < __LINESNUMBER__; ++ix) {
      fiLinesChannel[ix] = 0;
      ffLinesEnergy[ix] = 0;
   }
   fxDatabase = "calilines.txt";
   ReadDatabase();
}

//***********************************************************
TXXXCalibPar::~TXXXCalibPar()
{
  delete fxLinesFinder;
  delete fxCalibrator;
}

//***********************************************************
 void TXXXCalibPar::SetCalibSpectrum(TH1 *h1)
 {
   fxCalibSpectrum = h1;
   if(fxLinesFinder)
      fxLinesFinder->SetH1(__DATANAME__, fxCalibSpectrum, kFALSE);
 }

//***********************************************************
Bool_t TXXXCalibPar::UpdateFrom(TGo4Parameter *source)
{
   auto from = dynamic_cast<TXXXCalibPar *> (source);
   if (!from) {
      TGo4Log::Error("Wrong parameter class: %s", source->ClassName());
      return kFALSE;
   }

   for (Int_t ord = 0; ord < __POLORDER__; ++ord)
      fdA[ord] = from->fdA[ord];
   fbRecalibrate = from->fbRecalibrate;
   fbReadDatabase = from->fbReadDatabase;
   if(fxLinesFinder) delete fxLinesFinder;
   fxLinesFinder = from->fxLinesFinder;
   from->fxLinesFinder = nullptr; // adopt lines finder
   if(fxCalibrator) delete fxCalibrator;
   fxCalibrator = from->fxCalibrator;
   from->fxCalibrator = nullptr; // adopt calibration fitter

   // note: graph with calibration curve is not copied!

   for (Int_t ix = 0; ix < __LINESNUMBER__; ++ix) {
      fiLinesChannel[ix] = from->fiLinesChannel[ix];
      ffLinesEnergy[ix] = from->ffLinesEnergy[ix];
      fxLinesNames[ix] = from->fxLinesNames[ix];
   }
   std::cout << "Updated Parameter:" << std::endl;
   //Print();
   // get references to graph and histogram from analysis:
   // note that updatefrom is only used on analysis side here!
   fxCalibCurve = dynamic_cast<TGraph *>(TGo4Analysis::Instance()->GetObject(fxGraphName.Data()));

   if(!fxCalibCurve)
   {
      std::cout <<"Graph "<<fxGraphName.Data() << " not existing in analysis"<< std::endl;
   }
   else
   {
      std::cout <<"Updated graph pointer ref to "<<fxCalibCurve << std::endl;
   }

   // now reread database if desired:
   if(fbReadDatabase)
   {
      std::cout <<"Reread database" << std::endl;
      ReadDatabase();
   }

   if(fbRecalibrate)
   {
      std::cout <<"Recalibrating..." << std::endl;
      // first we get the channels from the linesfinder fitter:

      for (Int_t i = 0; i < __LINESNUMBER__; ++i) {
         const char *linename = fxLinesNames[i];
         TGo4FitModel *mod = fxLinesFinder->FindModel(linename);
         if(mod) {
            // check here if component is active or not
            if(mod->IsAssignTo(__DATANAME__))
               fiLinesChannel[i]=(Int_t) mod->GetParValue("Pos");
            else
               fiLinesChannel[i] = 0; // mark not active lines
         }
      }

      // setup calibration graph with the new channel coords:
      if(fxCalibCurve)
      {
         fxCalibCurve->Set(0);
         Int_t point = 0;
         for (Int_t ix = 0; ix < __LINESNUMBER__; ++ix) {
            if (fiLinesChannel[ix] != 0) {
               fxCalibCurve->SetPoint(point, fiLinesChannel[ix], ffLinesEnergy[ix]);
               // we only fit active lines
               ++point;
            }
         } // for
         // now perform fit of calibration graph:
         fxCalibrator->SetObject(__GRAPHNAME__, fxCalibCurve, kFALSE);
         fxCalibrator->DoActions();
         fxCalibrator->PrintLines();
         printf("fxCalibrator = %p\n", fxCalibrator);
         // finally, copy results of calibration to the parameter fields:
         for (Int_t i = 0; i < __POLORDER__; ++i) {
            TString modname = TString::Format("A_%d", i);
            TGo4FitModel *mod = fxCalibrator->FindModel(modname.Data());
            if (mod) {
               // check here if component is active or not
               if (mod->IsAssignTo(__GRAPHNAME__))
                  fdA[i] = mod->GetParValue("Ampl");
               else
                  fdA[i] = 0;
            }
         }
      }
      else
      {
         TGo4Log::Error("Calibration parameter %s has no TGraph!",
               GetName());
      }
   }

  return kTRUE;
}

//***********************************************************
void TXXXCalibPar::ReadDatabase()
{
   // read energies from file:
   char nextline[__TEXTMAX__];
   char buf[__TEXTMAX__];
   std::ifstream database(fxDatabase.Data());
   if(!database)
   {
      TGo4Log::Error("Open error of calibration energy file %s",
            fxDatabase.Data());
   }
   else
   {
      Int_t ix = 0;
      while(1){
         do{
            database.getline(nextline,__TEXTMAX__,'\n' ); // read whole line
            if(database.eof() || !database.good())
            {
               break;
            }
         }while(strstr(nextline,"#") || strstr(nextline,"!") ); // skip any comments
         if(database.eof() || !database.good()) break;
         sscanf(nextline,"%s %f %d",buf,
               &ffLinesEnergy[ix],
               &fiLinesChannel[ix]);
         fxLinesNames[ix]=buf;
         //      std::cout <<"\tname:"<<fxLinesNames[ix].Data() << std::endl;
         //      std::cout <<"\te:"<<ffLinesEnergy[ix] << std::endl;
         //      std::cout <<"\tch:"<<fiLinesChannel[ix] << std::endl;

         fxLinesFinder->AddGauss1(__DATANAME__,
               fxLinesNames[ix].Data(),
               fiLinesChannel[ix],
               TMath::Sqrt((Long_t) fiLinesChannel[ix]));
         ix++;
      } // while(1)
   }
}

//***********************************************************
Double_t TXXXCalibPar::Energy(Int_t channel)
{
   Double_t result = 0.;
   for (Int_t ord = 0; ord < __POLORDER__; ++ord)
      result += fdA[ord] * TMath::Power(channel, ord);

   return result;
}