PRunBase.cpp

Go to the documentation of this file.

/***************************************************************************
 
  PRunBase.cpp
 
  Author: Andreas Suter
  e-mail: andreas.suter@psi.ch
 
***************************************************************************/
 
/***************************************************************************
 *   Copyright (C) 2007-2026 by Andreas Suter                              *
 *   andreas.suter@psi.ch                                                  *
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 *   This program is distributed in the hope that it will be useful,       *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
 *   GNU General Public License for more details.                          *
 *                                                                         *
 *   You should have received a copy of the GNU General Public License     *
 *   along with this program; if not, write to the                         *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/
 
#include <iostream>
 
#include <TROOT.h>
#include <TSystem.h>
#include <TString.h>
#include <TObjArray.h>
#include <TObjString.h>
#include <TFolder.h>
 
#include "PRunBase.h"
 
//--------------------------------------------------------------------------
// Constructor
//--------------------------------------------------------------------------
PRunBase::PRunBase()
{
  fRunNo = -1;
  fRunInfo = nullptr;
  fRawData = nullptr;
  fMetaData.fField  = PMUSR_UNDEFINED;
  fMetaData.fEnergy = PMUSR_UNDEFINED;
  fTimeResolution = -1.0;
 
  fFitStartTime = PMUSR_UNDEFINED;
  fFitEndTime   = PMUSR_UNDEFINED;
 
  fValid = true;
  fHandleTag = kEmpty;
}
 
//--------------------------------------------------------------------------
// Constructor
//--------------------------------------------------------------------------
PRunBase::PRunBase(PMsrHandler *msrInfo, PRunDataHandler *rawData, UInt_t runNo, EPMusrHandleTag tag) :
    fHandleTag(tag), fMsrInfo(msrInfo), fRawData(rawData)
{
  fValid = true;
 
  fRunNo = static_cast<Int_t>(runNo);
  if (runNo > fMsrInfo->GetMsrRunList()->size()) {
    fRunInfo = nullptr;
    return;
  }
 
  // keep the run header info for this run
  fRunInfo = &(*fMsrInfo->GetMsrRunList())[runNo];
 
  // check the parameter and map range of the functions
  if (!fMsrInfo->CheckMapAndParamRange(static_cast<UInt_t>(fRunInfo->GetMap()->size()), fMsrInfo->GetNoOfParams())) {
    std::cerr << std::endl << "**SEVERE ERROR** PRunBase::PRunBase: map and/or parameter out of range in FUNCTIONS." << std::endl;
    exit(0);
  }
 
  // init private variables
  fMetaData.fField  = PMUSR_UNDEFINED;
  fMetaData.fEnergy = PMUSR_UNDEFINED;
  fTimeResolution = -1.0;
  for (Int_t i=0; i<fMsrInfo->GetNoOfFuncs(); i++)
    fFuncValues.push_back(0.0);
 
  // generate theory
  fTheory = std::make_unique<PTheory>(fMsrInfo, runNo);
  if (fTheory == nullptr) {
    std::cerr << std::endl << "**SEVERE ERROR** PRunBase::PRunBase: Couldn't create an instance of PTheory :-(, will quit" << std::endl;
    exit(0);
  }
  if (!fTheory->IsValid()) {
    std::cerr << std::endl << "**SEVERE ERROR** PRunBase::PRunBase: Theory is not valid :-(, will quit" << std::endl;
    exit(0);
  }
 
  // set fit time ranges
  fFitStartTime = PMUSR_UNDEFINED;
  fFitEndTime   = PMUSR_UNDEFINED;
}
 
//--------------------------------------------------------------------------
// Destructor
//--------------------------------------------------------------------------
PRunBase::~PRunBase()
{
  fT0s.clear();
  for (UInt_t i=0; i<fAddT0s.size(); i++)
    fAddT0s[i].clear();
  fAddT0s.clear();
 
  fFuncValues.clear();
}
 
//--------------------------------------------------------------------------
// DeadTimeCorrection (protected)
//--------------------------------------------------------------------------
void PRunBase::DeadTimeCorrection(std::vector<PDoubleVector> &histos, PUIntVector &histoNo)
{
  PRawRunData* runData = fRawData->GetRunData(*fRunInfo->GetRunName());
 
  if (runData == nullptr) {
    std::cerr << std::endl << "**WARNING** PRunBase::DeadTimeCorrection: couldn't get data from '" << *fRunInfo->GetRunName() << "'" << std::endl;
    return;
  }
 
  bool checkDeadTime{false};
  PMsrRunList *rl = fMsrInfo->GetMsrRunList();
  for (unsigned int i=0; i<rl->size(); i++) {
    if (!rl->at(i).GetFileFormat()->CompareTo("nexus", TString::kIgnoreCase)) {
      checkDeadTime = true;
      break;
    }
  }
 
  if (!runData->DeadTimeCorrectionReady() && checkDeadTime) {
    std::cerr << std::endl << "**WARNING** PRunBase::DeadTimeCorrection: missing input for dead time correction" << std::endl;
    return;
  }
 
 
  // check if a dead time correction has to be done
  // first check the global block
  TString dtcg = fMsrInfo->GetMsrGlobal()->GetDeadTimeCorrection();
  Int_t dtcg_tag = 0; // 0=no, 1=file, 2=estimate
  if (dtcg.Contains("file", TString::kIgnoreCase)) {
    dtcg_tag = 1;
  } else if (dtcg.Contains("estimate", TString::kIgnoreCase)) {
    dtcg_tag = 2;
  }
  // now check each run
  TString dtcr{"no"};
  Bool_t needToCheck{true};
  for (UInt_t i=0; i<histoNo.size(); i++) {
    dtcr = fRunInfo->GetDeadTimeCorrection();
    if (dtcr.Contains("file", TString::kIgnoreCase) || (dtcg_tag != 0)) {
      if (runData->DeadTimeCorrectionReady()) {
        needToCheck = false;
        // Dead time correction: n_true = n_obs / (1 - n_obs * t_dt / (good_frames * dt))
        Double_t n_true;
        Int_t gf = runData->GetNumberOfGoodFrames();
        std::vector<float> t_dt = runData->GetDeadTimeParam();
        for (UInt_t j=0; j<histos[i].size(); j++) {
          n_true = histos[i][j] / (1.0 - histos[i][j] * t_dt[histoNo[i]-1] / (gf * fTimeResolution));
          histos[i][j] = n_true;
        }
      }
    }
    if ((dtcr.Contains("estimate", TString::kIgnoreCase) || (dtcg_tag != 0)) && needToCheck) {
      std::cerr << std::endl << "**INFO** deadtime correction estimate not yet implemented." << std::endl;
    }
  }
}
 
//--------------------------------------------------------------------------
// SetFitRange (public)
//--------------------------------------------------------------------------
void PRunBase::SetFitRange(PDoublePairVector fitRange)
{
  Double_t start=0.0, end=0.0;
 
  assert(fitRange.size()); // make sure fitRange is not empty
 
  if (fitRange.size()==1) { // one fit range for all
    start = fitRange[0].first;
    end   = fitRange[0].second;
  } else {
    // check that fRunNo is found within fitRange
    UInt_t idx=static_cast<UInt_t>(fRunNo);
    if (idx < fitRange.size()) { // fRunNo present
      start = fitRange[idx].first;
      end   = fitRange[idx].second;
    } else { // fRunNo NOT present
      std::cerr << std::endl << ">> PRunBase::SetFitRange(): **ERROR** msr-file run entry " << fRunNo << " not present in fit range vector.";
      std::cerr << std::endl << ">>    Will not do anything! Please check, this shouldn't happen." << std::endl;
      return;
    }
  }
 
  // check that start is before end
  if (start > end) {
    std::cerr << std::endl << ">> PRunBase::SetFitRange(): **WARNING** start=" << start << " is > as end=" << end;
    std::cerr << std::endl << ">>    Will swap them, since otherwise chisq/logLH == 0.0" << std::endl;
    fFitStartTime = end;
    fFitEndTime   = start;
  } else {
    fFitStartTime = start;
    fFitEndTime   = end;
  }
}
 
//--------------------------------------------------------------------------
// CleanUp (public)
//--------------------------------------------------------------------------
void PRunBase::CleanUp()
{
  fTheory.reset();
}
 
//--------------------------------------------------------------------------
// CalculateKaiserFilterCoeff (protected)
//--------------------------------------------------------------------------
void PRunBase::CalculateKaiserFilterCoeff(Double_t wc, Double_t A, Double_t dw)
{
  Double_t beta;
  Double_t dt = fData.GetTheoryTimeStep();
  Int_t m;
 
  // estimate beta (see reference above, p.574ff)
  if (A > 50.0) {
    beta = 0.1102*(A-8.7);
  } else if ((A >= 21.0) && (A <= 50.0)) {
    beta = 0.5842*TMath::Power(A-21.0, 0.4) + 0.07886*(A-21.0);
  } else {
    beta = 0.0;
  }
 
  m = TMath::FloorNint((A-8.0)/(2.285*dw*TMath::Pi()));
  // make sure m is odd
  if (m % 2 == 0)
    m++;
 
  Double_t alpha = static_cast<Double_t>(m)/2.0;
  Double_t dval;
  Double_t dsum = 0.0;
  for (Int_t i=0; i<=m; i++) {
    dval  = TMath::Sin(wc*(i-alpha)*dt)/(TMath::Pi()*(i-alpha)*dt);
    dval *= TMath::BesselI0(beta*TMath::Sqrt(1.0-TMath::Power((i-alpha)*dt/alpha, 2.0)))/TMath::BesselI0(beta);
    dsum += dval;
    fKaiserFilter.push_back(dval);
  }  
  for (UInt_t i=0; i<=static_cast<UInt_t>(m); i++) {
    fKaiserFilter[i] /= dsum;
  }
}
 
//--------------------------------------------------------------------------
// FilterTheo (protected)
//--------------------------------------------------------------------------
void PRunBase::FilterTheo()
{
  PDoubleVector theoFiltered;
  Double_t dval = 0.0;
  const PDoubleVector *theo = fData.GetTheory();
 
  for (UInt_t i=0; i<theo->size(); i++) {
    for (UInt_t j=0; j<fKaiserFilter.size(); j++) {
      if (i<j)
        dval = 0.0;
      else
        dval += fKaiserFilter[j]*theo->at(i-j);
    }
    theoFiltered.push_back(dval);
    dval = 0.0;
  }
 
  fData.ReplaceTheory(theoFiltered);
 
  // shift time start by half the filter length
  dval = fData.GetTheoryTimeStart() - 0.5*static_cast<Double_t>(fKaiserFilter.size())*fData.GetTheoryTimeStep();
  fData.SetTheoryTimeStart(dval);
 
  theoFiltered.clear();
}