PRunSingleHistoRRF.cpp

Go to the documentation of this file.

/***************************************************************************
 
  PRunSingleHistoRRF.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.             *
 ***************************************************************************/
 
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
 
#ifdef HAVE_GOMP
#include <omp.h>
#endif
 
#include <cmath>
#include <iostream>
#include <fstream>
#include <memory>
#include <string>
#include <vector>
 
#include <TString.h>
#include <TH1F.h>
 
#include "PMusr.h"
#include "PStringUtils.h"
#include "PFourier.h"
#include "PRunSingleHistoRRF.h"
 
//--------------------------------------------------------------------------
// Constructor
//--------------------------------------------------------------------------
PRunSingleHistoRRF::PRunSingleHistoRRF() : PRunBase()
{
  fNoOfFitBins  = 0;
  fBackground = 0.0;
  fBkgErr = 1.0;
  fRRFPacking = -1;
  fTheoAsData = false;
 
  // the 2 following variables are need in case fit range is given in bins, and since
  // the fit range can be changed in the command block, these variables need to be accessible
  fGoodBins[0] = -1;
  fGoodBins[1] = -1;
  
  fN0EstimateEndTime = 1.0; // end time in (us) over which N0 is estimated.
}
 
//--------------------------------------------------------------------------
// Constructor
//--------------------------------------------------------------------------
PRunSingleHistoRRF::PRunSingleHistoRRF(PMsrHandler *msrInfo, PRunDataHandler *rawData, UInt_t runNo, EPMusrHandleTag tag, Bool_t theoAsData) :
  PRunBase(msrInfo, rawData, runNo, tag), fTheoAsData(theoAsData)
{
  fNoOfFitBins  = 0;
 
  PMsrGlobalBlock *global = msrInfo->GetMsrGlobal();
 
  if (!global->IsPresent()) {
    std::cerr << std::endl << ">> PRunSingleHistoRRF::PRunSingleHistoRRF(): **SEVERE ERROR**: no GLOBAL-block present!";
    std::cerr << std::endl << ">> For Single Histo RRF the GLOBAL-block is mandatory! Please fix this first.";
    std::cerr << std::endl;
    fValid = false;
    return;
  }
 
  if (!global->GetRRFUnit().CompareTo("??")) {
    std::cerr << std::endl << ">> PRunSingleHistoRRF::PRunSingleHistoRRF(): **SEVERE ERROR**: no RRF-Frequency found!";
    std::cerr << std::endl;
    fValid = false;
    return;
  }
 
  fRRFPacking = global->GetRRFPacking();
  if (fRRFPacking == -1) {
    std::cerr << std::endl << ">> PRunSingleHistoRRF::PRunSingleHistoRRF(): **SEVERE ERROR**: no RRF-Packing found!";
    std::cerr << std::endl;
    fValid = false;
    return;
  }
 
  // the 2 following variables are need in case fit range is given in bins, and since
  // the fit range can be changed in the command block, these variables need to be accessible
  fGoodBins[0] = -1;
  fGoodBins[1] = -1;
 
  fN0EstimateEndTime = 1.0; // end time in (us) over which N0 is estimated.
 
  if (!PrepareData()) {
    std::cerr << std::endl << ">> PRunSingleHistoRRF::PRunSingleHistoRRF(): **SEVERE ERROR**: Couldn't prepare data for fitting!";
    std::cerr << std::endl << ">> This is very bad :-(, will quit ...";
    std::cerr << std::endl;
    fValid = false;
  }
}
 
//--------------------------------------------------------------------------
// Destructor
//--------------------------------------------------------------------------
PRunSingleHistoRRF::~PRunSingleHistoRRF()
{
  fForward.clear();
}
 
//--------------------------------------------------------------------------
// CalcChiSquare (public)
//--------------------------------------------------------------------------
Double_t PRunSingleHistoRRF::CalcChiSquare(const std::vector<Double_t>& par)
{
  Double_t chisq     = 0.0;
  Double_t diff      = 0.0;
 
  // calculate functions
  for (Int_t i=0; i<fMsrInfo->GetNoOfFuncs(); i++) {
    UInt_t funcNo = fMsrInfo->GetFuncNo(i);
    fFuncValues[i] = fMsrInfo->EvalFunc(funcNo, *fRunInfo->GetMap(), par, fMetaData);
  }
 
  // calculate chi square
  Double_t time(1.0);
  Int_t i;
 
  // Calculate the theory function once to ensure one function evaluation for the current set of parameters.
  // This is needed for the LF and user functions where some non-thread-save calculations only need to be calculated once
  // for a given set of parameters---which should be done outside of the parallelized loop.
  // For all other functions it means a tiny and acceptable overhead.
  time = fTheory->Func(time, par, fFuncValues);
 
  #ifdef HAVE_GOMP
  Int_t chunk = (fEndTimeBin - fStartTimeBin)/omp_get_num_procs();
  if (chunk < 10)
    chunk = 10;
  #pragma omp parallel for default(shared) private(i,time,diff) schedule(dynamic,chunk) reduction(+:chisq)
  #endif
  for (i=fStartTimeBin; i<fEndTimeBin; ++i) {
    time = fData.GetDataTimeStart() + static_cast<Double_t>(i)*fData.GetDataTimeStep();
    diff = fData.GetValue()->at(i) - fTheory->Func(time, par, fFuncValues);
    chisq += diff*diff / (fData.GetError()->at(i)*fData.GetError()->at(i));
  }
 
  return chisq;
}
 
//--------------------------------------------------------------------------
// CalcChiSquareExpected (public)
//--------------------------------------------------------------------------
Double_t PRunSingleHistoRRF::CalcChiSquareExpected(const std::vector<Double_t>& par)
{
  Double_t chisq = 0.0;
  Double_t diff  = 0.0;
  Double_t theo  = 0.0;
 
  // calculate functions
  for (Int_t i=0; i<fMsrInfo->GetNoOfFuncs(); i++) {
    UInt_t funcNo = fMsrInfo->GetFuncNo(i);
    fFuncValues[i] = fMsrInfo->EvalFunc(funcNo, *fRunInfo->GetMap(), par, fMetaData);
  }
 
  // calculate chi square
  Double_t time(1.0);
  Int_t i;
 
  // Calculate the theory function once to ensure one function evaluation for the current set of parameters.
  // This is needed for the LF and user functions where some non-thread-save calculations only need to be calculated once
  // for a given set of parameters---which should be done outside of the parallelized loop.
  // For all other functions it means a tiny and acceptable overhead.
  time = fTheory->Func(time, par, fFuncValues);
 
  #ifdef HAVE_GOMP
  Int_t chunk = (fEndTimeBin - fStartTimeBin)/omp_get_num_procs();
  if (chunk < 10)
    chunk = 10;
  #pragma omp parallel for default(shared) private(i,time,diff) schedule(dynamic,chunk) reduction(+:chisq)
  #endif
  for (i=fStartTimeBin; i < fEndTimeBin; ++i) {
    time = fData.GetDataTimeStart() + static_cast<Double_t>(i)*fData.GetDataTimeStep();
    theo = fTheory->Func(time, par, fFuncValues);
    diff = fData.GetValue()->at(i) - theo;
    chisq += diff*diff / theo;
  }
 
  return chisq;
}
 
//--------------------------------------------------------------------------
// CalcMaxLikelihood (public)
//--------------------------------------------------------------------------
Double_t PRunSingleHistoRRF::CalcMaxLikelihood(const std::vector<Double_t>& par)
{
  // not yet implemented
 
  return 0.0;
}
 
//--------------------------------------------------------------------------
// CalcTheory (public)
//--------------------------------------------------------------------------
void PRunSingleHistoRRF::CalcTheory()
{
  // feed the parameter vector
  std::vector<Double_t> par;
  PMsrParamList *paramList = fMsrInfo->GetMsrParamList();
  for (UInt_t i=0; i<paramList->size(); i++)
    par.push_back((*paramList)[i].fValue);
 
  // calculate functions
  for (Int_t i=0; i<fMsrInfo->GetNoOfFuncs(); i++) {
    fFuncValues[i] = fMsrInfo->EvalFunc(fMsrInfo->GetFuncNo(i), *fRunInfo->GetMap(), par, fMetaData);
  }
 
  // calculate theory
  UInt_t size = fData.GetValue()->size();
  Double_t start = fData.GetDataTimeStart();
  Double_t resolution = fData.GetDataTimeStep();
  Double_t time;
  for (UInt_t i=0; i<size; i++) {
    time = start + static_cast<Double_t>(i)*resolution;
    fData.AppendTheoryValue(fTheory->Func(time, par, fFuncValues));
  }
 
  // clean up
  par.clear();
}
 
//--------------------------------------------------------------------------
// GetNoOfFitBins (public)
//--------------------------------------------------------------------------
UInt_t PRunSingleHistoRRF::GetNoOfFitBins()
{
  CalcNoOfFitBins();
 
  return fNoOfFitBins;
}
 
//--------------------------------------------------------------------------
// SetFitRangeBin (public)
//--------------------------------------------------------------------------
void PRunSingleHistoRRF::SetFitRangeBin(const TString fitRange)
{
  TString str;
  Ssiz_t idx = -1;
  Int_t offset = 0;
 
  std::vector<std::string> tok = PStringUtils::Split(fitRange.Data(), " \t");
 
  if (tok.size() == 3) { // structure FIT_RANGE fgb+n0 lgb-n1
    // handle fgb+n0 entry
    str = tok[1];
    // check if there is an offset present
    idx = str.First("+");
    if (idx != -1) { // offset present
      str.Remove(0, idx+1);
      if (str.IsFloat()) // if str is a valid number, convert is to an integer
        offset = str.Atoi();
    }
    fFitStartTime = (fGoodBins[0] + offset - fT0s[0]) * fTimeResolution;
 
    // handle lgb-n1 entry
    str = tok[2];
    // check if there is an offset present
    idx = str.First("-");
    if (idx != -1) { // offset present
      str.Remove(0, idx+1);
      if (str.IsFloat()) // if str is a valid number, convert is to an integer
        offset = str.Atoi();
    }
    fFitEndTime = (fGoodBins[1] - offset - fT0s[0]) * fTimeResolution;
  } else if ((tok.size() > 3) && (tok.size() % 2 == 1)) { // structure FIT_RANGE fgb[+n00] lgb[-n01] [fgb[+n10] lgb[-n11] ... fgb[+nN0] lgb[-nN1]]
    UInt_t pos = 2*(fRunNo+1)-1;
 
    if (pos + 1 >= tok.size()) {
      std::cerr << std::endl << ">> PRunSingleHistoRRF::SetFitRangeBin(): **ERROR** invalid FIT_RANGE command found: '" << fitRange << "'";
      std::cerr << std::endl << ">> will ignore it. Sorry ..." << std::endl;
    } else {
      // handle fgb+n0 entry
      str = tok[pos];
      // check if there is an offset present
      idx = str.First("+");
      if (idx != -1) { // offset present
        str.Remove(0, idx+1);
        if (str.IsFloat()) // if str is a valid number, convert is to an integer
          offset = str.Atoi();
      }
      fFitStartTime = (fGoodBins[0] + offset - fT0s[0]) * fTimeResolution;
 
      // handle lgb-n1 entry
      str = tok[pos+1];
      // check if there is an offset present
      idx = str.First("-");
      if (idx != -1) { // offset present
        str.Remove(0, idx+1);
        if (str.IsFloat()) // if str is a valid number, convert is to an integer
          offset = str.Atoi();
      }
      fFitEndTime = (fGoodBins[1] - offset - fT0s[0]) * fTimeResolution;
    }
  } else { // error
    std::cerr << std::endl << ">> PRunSingleHistoRRF::SetFitRangeBin(): **ERROR** invalid FIT_RANGE command found: '" << fitRange << "'";
    std::cerr << std::endl << ">> will ignore it. Sorry ..." << std::endl;
  }
}
 
//--------------------------------------------------------------------------
// CalcNoOfFitBins (public)
//--------------------------------------------------------------------------
void PRunSingleHistoRRF::CalcNoOfFitBins()
{
  // In order not having to loop over all bins and to stay consistent with the chisq method, calculate the start and end bins explicitly
  fStartTimeBin = static_cast<Int_t>(ceil((fFitStartTime - fData.GetDataTimeStart())/fData.GetDataTimeStep()));
  if (fStartTimeBin < 0)
    fStartTimeBin = 0;
  fEndTimeBin = static_cast<Int_t>(floor((fFitEndTime - fData.GetDataTimeStart())/fData.GetDataTimeStep())) + 1;
  if (fEndTimeBin > static_cast<Int_t>(fData.GetValue()->size()))
    fEndTimeBin = fData.GetValue()->size();
 
  if (fEndTimeBin > fStartTimeBin)
    fNoOfFitBins = fEndTimeBin - fStartTimeBin;
  else
    fNoOfFitBins = 0;
}
 
//--------------------------------------------------------------------------
// PrepareData (protected)
//--------------------------------------------------------------------------
Bool_t PRunSingleHistoRRF::PrepareData()
{
  Bool_t success = true;
 
  if (!fValid)
    return false;
 
  // keep the Global block info
  PMsrGlobalBlock *globalBlock = fMsrInfo->GetMsrGlobal();
 
  // get the proper run
  PRawRunData* runData = fRawData->GetRunData(*fRunInfo->GetRunName());
  if (!runData) { // couldn't get run
    std::cerr << std::endl << ">> PRunSingleHistoRRF::PrepareData(): **ERROR** Couldn't get run " << fRunInfo->GetRunName()->Data() << "!";
    std::cerr << std::endl;
    return false;
  }
 
  // keep the field from the meta-data from the data-file
  fMetaData.fField = runData->GetField();
 
  // keep the energy from the meta-data from the data-file
  fMetaData.fEnergy = runData->GetEnergy();
 
  // keep the temperature(s) from the meta-data from the data-file
  for (unsigned int i=0; i<runData->GetNoOfTemperatures(); i++)
    fMetaData.fTemp.push_back(runData->GetTemperature(i));
 
  // collect histogram numbers
  PUIntVector histoNo; // histoNo = msr-file forward + redGreen_offset - 1
  for (UInt_t i=0; i<fRunInfo->GetForwardHistoNoSize(); i++) {
    histoNo.push_back(fRunInfo->GetForwardHistoNo(i));
 
    if (!runData->IsPresent(histoNo[i])) {
      std::cerr << std::endl << ">> PRunSingleHistoRRF::PrepareData(): **PANIC ERROR**:";
      std::cerr << std::endl << ">> histoNo found = " << histoNo[i] << ", which is NOT present in the data file!?!?";
      std::cerr << std::endl << ">> Will quit :-(";
      std::cerr << std::endl;
      histoNo.clear();
      return false;
    }
  }
 
  // keep the time resolution in (us)
  fTimeResolution = runData->GetTimeResolution()/1.0e3;
  std::cout.precision(10);
  std::cout << std::endl << ">> PRunSingleHisto::PrepareData(): time resolution=" << std::fixed << runData->GetTimeResolution() << "(ns)" << std::endl;
 
  // get all the proper t0's and addt0's for the current RUN block
  if (!GetProperT0(runData, globalBlock, histoNo)) {
    return false;
  }
 
  // keep the histo of each group at this point (addruns handled below)
  std::vector<PDoubleVector> forward;
  forward.resize(histoNo.size());   // resize to number of groups
  for (UInt_t i=0; i<histoNo.size(); i++) {
    forward[i].resize(runData->GetDataBin(histoNo[i])->size());
    forward[i] = *runData->GetDataBin(histoNo[i]);
  }
 
  // check if there are runs to be added to the current one
  if (fRunInfo->GetRunNameSize() > 1) { // runs to be added present
    PRawRunData *addRunData;
    std::vector<PDoubleVector> addForward;
    for (UInt_t i=1; i<fRunInfo->GetRunNameSize(); i++) {
 
      // get run to be added to the main one
      addRunData = fRawData->GetRunData(*fRunInfo->GetRunName(i));
      if (addRunData == nullptr) { // couldn't get run
        std::cerr << std::endl << ">> PRunSingleHistoRRF::PrepareData(): **ERROR** Couldn't get addrun " << fRunInfo->GetRunName(i)->Data() << "!";
        std::cerr << std::endl;
        return false;
      }
 
      addForward.clear();
      addForward.resize(histoNo.size());   // resize to number of groups
      for (UInt_t j=0; j<histoNo.size(); j++) {
        addForward[j].resize(addRunData->GetDataBin(histoNo[j])->size());
        addForward[j] = *addRunData->GetDataBin(histoNo[j]);
      }
      DeadTimeCorrection(addForward, histoNo);
 
      // add forward run
      UInt_t addRunSize;
      for (UInt_t k=0; k<histoNo.size(); k++) { // fill each group
        addRunSize = addForward[k].size();
        for (UInt_t j=0; j<addRunSize; j++) { // loop over the bin indices
          // make sure that the index stays in the proper range
          if ((static_cast<Int_t>(j)+static_cast<Int_t>(fAddT0s[i-1][k])-static_cast<Int_t>(fT0s[k]) >= 0) &&
              (j+static_cast<Int_t>(fAddT0s[i-1][k])-static_cast<Int_t>(fT0s[k]) < addRunSize)) {
            forward[k][j] += addForward[k][j+static_cast<Int_t>(fAddT0s[i-1][k])-static_cast<Int_t>(fT0s[k])];
          }
        }
      }
    }
  }
 
  // set forward histo data of the first group
  fForward.resize(forward[0].size());
  for (UInt_t i=0; i<fForward.size(); i++) {
    fForward[i] = forward[0][i];
  }
 
  // group histograms, add all the remaining forward histograms of the group
  for (UInt_t i=1; i<histoNo.size(); i++) { // loop over the groupings
    for (UInt_t j=0; j<runData->GetDataBin(histoNo[i])->size(); j++) { // loop over the bin indices
      // make sure that the index stays within proper range
      if ((static_cast<Int_t>(j)+static_cast<Int_t>(fT0s[i])-static_cast<Int_t>(fT0s[0]) >= 0) &&
          (j+static_cast<Int_t>(fT0s[i])-static_cast<Int_t>(fT0s[0]) < runData->GetDataBin(histoNo[i])->size())) {
        fForward[j] += forward[i][j+static_cast<Int_t>(fT0s[i])-static_cast<Int_t>(fT0s[0])];
      }
    }
  }
 
  // get the data range (fgb/lgb) for the current RUN block
  if (!GetProperDataRange()) {
    return false;
  }
 
  // get the fit range for the current RUN block
  GetProperFitRange(globalBlock);
 
  // do the more fit/view specific stuff
  if (fHandleTag == kFit)
    success = PrepareFitData(runData, histoNo[0]);
  else if (fHandleTag == kView)
    success = PrepareViewData(runData, histoNo[0]);
  else
    success = false;
 
  // cleanup
  histoNo.clear();
 
  return success;
}
 
//--------------------------------------------------------------------------
// PrepareFitData (protected)
//--------------------------------------------------------------------------
Bool_t PRunSingleHistoRRF::PrepareFitData(PRawRunData* runData, const UInt_t histoNo)
{
  // keep the raw data for the RRF asymmetry error estimate for later
  PDoubleVector rawNt;
  for (UInt_t i=0; i<fForward.size(); i++) {
    rawNt.push_back(fForward[i]);  // N(t) without any corrections
  }
  Double_t freqMax = GetMainFrequency(rawNt);
  std::cout << "info> freqMax=" << freqMax << " (MHz)" << std::endl;
 
  // "optimal packing"
  Double_t optNoPoints = 8;
  if (freqMax < 271.0) // < 271 MHz, i.e ~ 2T
    optNoPoints = 5;
  std::cout << "info> optimal packing: " << static_cast<Int_t>(1.0 / (fTimeResolution*(freqMax - fMsrInfo->GetMsrGlobal()->GetRRFFreq("MHz"))) / optNoPoints);
 
  // initially fForward is the "raw data set" (i.e. grouped histo and raw runs already added) to be fitted. This means fForward = N(t) at this point.
 
  // 1) check how the background shall be handled
  // subtract background from histogramms ------------------------------------------
  if (fRunInfo->GetBkgFix(0) == PMUSR_UNDEFINED) { // no fixed background given
    if (fRunInfo->GetBkgRange(0) >= 0) {
      if (!EstimateBkg(histoNo))
        return false;
    } else { // no background given to do the job, try estimate
      fRunInfo->SetBkgRange(static_cast<Int_t>(fT0s[0]*0.1), 0);
      fRunInfo->SetBkgRange(static_cast<Int_t>(fT0s[0]*0.6), 1);
      std::cerr << std::endl << ">> PRunSingleHistoRRF::PrepareFitData(): **WARNING** Neither fix background nor background bins are given!";
      std::cerr << std::endl << ">> Will try the following: bkg start = " << fRunInfo->GetBkgRange(0) << ", bkg end = " << fRunInfo->GetBkgRange(1);
      std::cerr << std::endl << ">> NO WARRANTY THAT THIS MAKES ANY SENSE! Better check ...";
      std::cerr << std::endl;
      if (!EstimateBkg(histoNo))
        return false;
    }
    // subtract background from fForward
    for (UInt_t i=0; i<fForward.size(); i++)
      fForward[i] -= fBackground;
  } else { // fixed background given
    for (UInt_t i=0; i<fForward.size(); i++) {
      fForward[i] -= fRunInfo->GetBkgFix(0);
    }
    fBackground = fRunInfo->GetBkgFix(0);
  }
  // here fForward = N(t) - Nbkg
 
  Int_t t0 = static_cast<Int_t>(fT0s[0]);
 
  // 2) N(t) - Nbkg -> exp(+t/tau) [N(t)-Nbkg]
  Double_t startTime = fTimeResolution * (static_cast<Double_t>(fGoodBins[0]) - static_cast<Double_t>(t0));
 
  Double_t time_tau=0.0;
  Double_t exp_t_tau=0.0;
  for (Int_t i=fGoodBins[0]; i<fGoodBins[1]; i++) {
    time_tau = (startTime + fTimeResolution * (i - fGoodBins[0])) / PMUON_LIFETIME;
    exp_t_tau = exp(time_tau);
    fForward[i] *= exp_t_tau;
    fM.push_back(fForward[i]);   // i.e. M(t) = [N(t)-Nbkg] exp(+t/tau); needed to estimate N0 later on
    fMerr.push_back(exp_t_tau*sqrt(rawNt[i]+fBkgErr*fBkgErr));
  }
 
  // calculate weights
  for (UInt_t i=0; i<fMerr.size(); i++) {
    if (fMerr[i] > 0.0)
      fW.push_back(1.0/(fMerr[i]*fMerr[i]));
    else
      fW.push_back(1.0);
  }
  // now fForward = exp(+t/tau) [N(t)-Nbkg] = M(t)
 
  // 3) estimate N0
  Double_t errN0 = 0.0;
  Double_t n0 = EstimateN0(errN0, freqMax);
 
  // 4a) A(t) = exp(+t/tau) [N(t)-Nbkg] / N0 - 1.0
  for (Int_t i=fGoodBins[0]; i<=fGoodBins[1]; i++) {
    fForward[i] = fForward[i] / n0 - 1.0;
  }
 
  // 4b) error estimate of A(t): errA(t) = exp(+t/tau)/N0 sqrt( N(t) + ([N(t)-N_bkg]/N0)^2 errN0^2 )
  for (Int_t i=fGoodBins[0]; i<=fGoodBins[1]; i++) {
    time_tau = (startTime + fTimeResolution * (i - fGoodBins[0])) / PMUON_LIFETIME;
    exp_t_tau = exp(time_tau);
    fAerr.push_back(exp_t_tau/n0*sqrt(rawNt[i]+pow(((rawNt[i]-fBackground)/n0)*errN0,2.0)));
  }
 
  // 5) rotate A(t): A(t) -> 2* A(t) * cos(wRRF t + phiRRF), the factor 2.0 is needed since the high frequency part is suppressed.
  PMsrGlobalBlock *globalBlock = fMsrInfo->GetMsrGlobal();
  Double_t wRRF = globalBlock->GetRRFFreq("Mc");
  Double_t phaseRRF = globalBlock->GetRRFPhase()*TMath::TwoPi()/180.0;
  Double_t time = 0.0;
  for (Int_t i=fGoodBins[0]; i<=fGoodBins[1]; i++) {
    time = startTime + fTimeResolution * (static_cast<Double_t>(i) - static_cast<Double_t>(fGoodBins[0]));
    fForward[i] *= 2.0*cos(wRRF * time + phaseRRF);
  }
 
  // 6) RRF packing
  Double_t dval=0.0;
  for (Int_t i=fGoodBins[0]; i<=fGoodBins[1]; i++) {
    if (fRRFPacking == 1) {
      fData.AppendValue(fForward[i]);
    } else { // RRF packing > 1
      if (((i-fGoodBins[0]) % fRRFPacking == 0) && (i != fGoodBins[0])) { // fill data
        dval /= fRRFPacking;
        fData.AppendValue(dval);
        // reset dval
        dval = 0.0;
      }
      dval += fForward[i];
    }
  }
 
  // 7) estimate packed RRF errors (see log-book p.204)
  //    the error estimate of the unpacked RRF asymmetry is: errA_RRF(t) \simeq exp(t/tau)/N0 sqrt( [N(t) + ((N(t)-N_bkg)/N0)^2 errN0^2] )
  dval = 0.0;
  // the packed RRF asymmetry error
  for (Int_t i=fGoodBins[0]; i<=fGoodBins[1]; i++) {
    if (((i-fGoodBins[0]) % fRRFPacking == 0) && (i != fGoodBins[0])) { // fill data
      fData.AppendErrorValue(sqrt(2.0*dval)/fRRFPacking); // the factor 2.0 is needed since the high frequency part is suppressed.
      dval = 0.0;
    }
    dval += fAerr[i-fGoodBins[0]]*fAerr[i-fGoodBins[0]];
  }
 
  // set start time and time step
  fData.SetDataTimeStart(fTimeResolution*(static_cast<Double_t>(fGoodBins[0])-static_cast<Double_t>(t0)+static_cast<Double_t>(fRRFPacking-1)/2.0));
  fData.SetDataTimeStep(fTimeResolution*fRRFPacking);
 
  CalcNoOfFitBins();
 
  return true;
}
 
//--------------------------------------------------------------------------
// PrepareViewData (protected)
//--------------------------------------------------------------------------
Bool_t PRunSingleHistoRRF::PrepareViewData(PRawRunData* runData, const UInt_t histoNo)
{
  // --------------
  // prepare data
  // --------------
 
  // prepare RRF single histo
  PrepareFitData(runData, histoNo);
 
  // check for view packing
  Int_t viewPacking = fMsrInfo->GetMsrPlotList()->at(0).fViewPacking;
  if (viewPacking > 0) {
    if (viewPacking < fRRFPacking) {
      std::cerr << ">> PRunSingleHistoRRF::PrepareViewData(): **WARNING** Found View Packing (" << viewPacking << ") < RRF Packing (" << fRRFPacking << ").";
      std::cerr << ">> Will ignore View Packing." << std::endl;
    } else {
      // STILL MISSING
    }
  }
 
  // --------------
  // prepare theory
  // --------------
 
  // feed the parameter vector
  std::vector<Double_t> par;
  PMsrParamList *paramList = fMsrInfo->GetMsrParamList();
  for (UInt_t i=0; i<paramList->size(); i++)
    par.push_back((*paramList)[i].fValue);
 
  // calculate functions
  for (Int_t i=0; i<fMsrInfo->GetNoOfFuncs(); i++) {
    fFuncValues[i] = fMsrInfo->EvalFunc(fMsrInfo->GetFuncNo(i), *fRunInfo->GetMap(), par, fMetaData);
  }
 
  // check if a finer binning for the theory is needed
  UInt_t size = fForward.size();
  Int_t factor = 8; // 8 times more points for the theory (if fTheoAsData == false)
  fData.SetTheoryTimeStart(fData.GetDataTimeStart());
  if (fTheoAsData) { // calculate theory only at the data points
    fData.SetTheoryTimeStep(fData.GetDataTimeStep());
  } else {
    // finer binning for the theory (8 times as many points = factor)
    size *= factor;
    fData.SetTheoryTimeStep(fData.GetDataTimeStep()/(Double_t)factor);
  }
 
  // calculate theory
  Double_t time = 0.0;
  Double_t theoryValue = 0.0;
  for (UInt_t i=0; i<size; i++) {
    time = fData.GetTheoryTimeStart() + static_cast<Double_t>(i)*fData.GetTheoryTimeStep();
    theoryValue = fTheory->Func(time, par, fFuncValues);
    if (fabs(theoryValue) > 10.0) { // dirty hack needs to be fixed!!
      theoryValue = 0.0;
    }
    fData.AppendTheoryValue(theoryValue);
  }
 
  return true;
}
 
//--------------------------------------------------------------------------
// GetProperT0 (private)
//--------------------------------------------------------------------------
Bool_t PRunSingleHistoRRF::GetProperT0(PRawRunData* runData, PMsrGlobalBlock *globalBlock, PUIntVector &histoNo)
{
  // feed all T0's
  // first init T0's, T0's are stored as (forward T0, backward T0, etc.)
  fT0s.clear();
  fT0s.resize(histoNo.size());
  for (UInt_t i=0; i<fT0s.size(); i++) {
    fT0s[i] = -1.0;
  }
 
  // fill in the T0's from the msr-file (if present)
  for (UInt_t i=0; i<fRunInfo->GetT0BinSize(); i++) {
    fT0s[i] = fRunInfo->GetT0Bin(i);
  }
 
  // fill in the T0's from the GLOBAL block section (if present)
  for (UInt_t i=0; i<globalBlock->GetT0BinSize(); i++) {
    if (fT0s[i] == -1.0) { // i.e. not given in the RUN block section
      fT0s[i] = globalBlock->GetT0Bin(i);
    }
  }
 
  // fill in the T0's from the data file, if not already present in the msr-file
  for (UInt_t i=0; i<histoNo.size(); i++) {
    if (fT0s[i] == -1.0) { // i.e. not present in the msr-file, try the data file
      if (runData->GetT0Bin(histoNo[i]) > 0.0) {
        fT0s[i] = runData->GetT0Bin(histoNo[i]);
        fRunInfo->SetT0Bin(fT0s[i], i); // keep value for the msr-file
      }
    }
  }
 
  // fill in the T0's gaps, i.e. in case the T0's are NOT in the msr-file and NOT in the data file
  for (UInt_t i=0; i<histoNo.size(); i++) {
    if (fT0s[i] == -1.0) { // i.e. not present in the msr-file and data file, use the estimated T0
      fT0s[i] = runData->GetT0BinEstimated(histoNo[i]);
      fRunInfo->SetT0Bin(fT0s[i], i); // keep value for the msr-file
 
      std::cerr << std::endl << ">> PRunSingleHistoRRF::GetProperT0(): **WARRNING** NO t0's found, neither in the run data nor in the msr-file!";
      std::cerr << std::endl << ">> run: " << fRunInfo->GetRunName()->Data();
      std::cerr << std::endl << ">> will try the estimated one: forward t0 = " << runData->GetT0BinEstimated(histoNo[i]);
      std::cerr << std::endl << ">> NO WARRANTY THAT THIS OK!! For instance for LEM this is almost for sure rubbish!";
      std::cerr << std::endl;
    }
  }
 
  // check if t0 is within proper bounds
  for (UInt_t i=0; i<fRunInfo->GetForwardHistoNoSize(); i++) {
    if ((fT0s[i] < 0.0) || (fT0s[i] > static_cast<Int_t>(runData->GetDataBin(histoNo[i])->size()))) {
      std::cerr << std::endl << ">> PRunSingleHistoRRF::GetProperT0(): **ERROR** t0 data bin (" << fT0s[i] << ") doesn't make any sense!";
      std::cerr << std::endl;
      return false;
    }
  }
 
  // check if there are runs to be added to the current one. If yes keep the needed t0's
  if (fRunInfo->GetRunNameSize() > 1) { // runs to be added present
    PRawRunData *addRunData;
    fAddT0s.resize(fRunInfo->GetRunNameSize()-1); // resize to the number of addruns
    for (UInt_t i=1; i<fRunInfo->GetRunNameSize(); i++) {
 
      // get run to be added to the main one
      addRunData = fRawData->GetRunData(*fRunInfo->GetRunName(i));
      if (addRunData == nullptr) { // couldn't get run
        std::cerr << std::endl << ">> PRunSingleHistoRRF::GetProperT0(): **ERROR** Couldn't get addrun " << fRunInfo->GetRunName(i)->Data() << "!";
        std::cerr << std::endl;
        return false;
      }
 
      // feed all T0's
      // first init T0's, T0's are stored as (forward T0, backward T0, etc.)
      fAddT0s[i-1].resize(histoNo.size());
      for (UInt_t j=0; j<fAddT0s[i-1].size(); j++) {
        fAddT0s[i-1][j] = -1.0;
      }
 
      // fill in the T0's from the msr-file (if present)
      for (UInt_t j=0; j<fRunInfo->GetT0BinSize(); j++) {
        fAddT0s[i-1][j] = fRunInfo->GetAddT0Bin(i-1,j); // addRunIdx starts at 0
      }
 
      // fill in the T0's from the data file, if not already present in the msr-file
      for (UInt_t j=0; j<histoNo.size(); j++) {
        if (fAddT0s[i-1][j] == -1.0) // i.e. not present in the msr-file, try the data file
          if (addRunData->GetT0Bin(histoNo[j]) > 0.0) {
            fAddT0s[i-1][j] = addRunData->GetT0Bin(histoNo[j]);
            fRunInfo->SetAddT0Bin(fAddT0s[i-1][j], i-1, j); // keep value for the msr-file
          }
      }
 
      // fill in the T0's gaps, i.e. in case the T0's are NOT in the msr-file and NOT in the data file
      for (UInt_t j=0; j<histoNo.size(); j++) {
        if (fAddT0s[i-1][j] == -1.0) { // i.e. not present in the msr-file and data file, use the estimated T0
          fAddT0s[i-1][j] = addRunData->GetT0BinEstimated(histoNo[j]);
          fRunInfo->SetAddT0Bin(fAddT0s[i-1][j], i-1, j); // keep value for the msr-file
 
          std::cerr << std::endl << ">> PRunSingleHistoRRF::GetProperT0(): **WARRNING** NO t0's found, neither in the run data nor in the msr-file!";
          std::cerr << std::endl << ">> run: " << fRunInfo->GetRunName(i)->Data();
          std::cerr << std::endl << ">> will try the estimated one: forward t0 = " << addRunData->GetT0BinEstimated(histoNo[j]);
          std::cerr << std::endl << ">> NO WARRANTY THAT THIS OK!! For instance for LEM this is almost for sure rubbish!";
          std::cerr << std::endl;
        }
      }
 
      // check if t0 is within proper bounds
      for (UInt_t j=0; j<fRunInfo->GetForwardHistoNoSize(); j++) {
        if ((fAddT0s[i-1][j] < 0) || (fAddT0s[i-1][j] > static_cast<Int_t>(addRunData->GetDataBin(histoNo[j])->size()))) {
          std::cerr << std::endl << ">> PRunSingleHistoRRF::GetProperT0(): **ERROR** addt0 data bin (" << fAddT0s[i-1][j] << ") doesn't make any sense!";
          std::cerr << std::endl;
          return false;
        }
      }
    }
  }
 
  return true;
}
 
//--------------------------------------------------------------------------
// GetProperDataRange (private)
//--------------------------------------------------------------------------
Bool_t PRunSingleHistoRRF::GetProperDataRange()
{
  // get start/end data
  Int_t start;
  Int_t end;
  start = fRunInfo->GetDataRange(0);
  end   = fRunInfo->GetDataRange(1);
 
  // check if data range has been given in the RUN block, if not try to get it from the GLOBAL block
  if (start < 0) {
    start = fMsrInfo->GetMsrGlobal()->GetDataRange(0);
  }
  if (end < 0) {
    end = fMsrInfo->GetMsrGlobal()->GetDataRange(1);
  }
 
  // check if data range has been provided, and if not try to estimate them
  if (start < 0) {
    Int_t offset = static_cast<Int_t>(10.0e-3/fTimeResolution);
    start = static_cast<Int_t>(fT0s[0])+offset;
    fRunInfo->SetDataRange(start, 0);
    std::cerr << std::endl << ">> PRunSingleHistoRRF::GetProperDataRange(): **WARNING** data range was not provided, will try data range start = t0+" << offset << "(=10ns) = " << start << ".";
    std::cerr << std::endl << ">> NO WARRANTY THAT THIS DOES MAKE ANY SENSE.";
    std::cerr << std::endl;
  }
  if (end < 0) {
    end = fForward.size();
    fRunInfo->SetDataRange(end, 1);
    std::cerr << std::endl << ">> PRunSingleHistoRRF::GetProperDataRange(): **WARNING** data range was not provided, will try data range end = " << end << ".";
    std::cerr << std::endl << ">> NO WARRANTY THAT THIS DOES MAKE ANY SENSE.";
    std::cerr << std::endl;
  }
 
  // check if start and end make any sense
  // 1st check if start and end are in proper order
  if (end < start) { // need to swap them
    Int_t keep = end;
    end = start;
    start = keep;
  }
  // 2nd check if start is within proper bounds
  if ((start < 0) || (start > static_cast<Int_t>(fForward.size()))) {
    std::cerr << std::endl << ">> PRunSingleHistoRRF::GetProperDataRange(): **ERROR** start data bin (" << start << ") doesn't make any sense!";
    std::cerr << std::endl;
    return false;
  }
  // 3rd check if end is within proper bounds
  if (end < 0) {
    std::cerr << std::endl << ">> PRunSingleHistoRRF::GetProperDataRange(): **ERROR** end data bin (" << end << ") doesn't make any sense!";
    std::cerr << std::endl;
    return false;
  }
  if (end > static_cast<Int_t>(fForward.size())) {
    std::cerr << std::endl << ">> PRunSingleHistoRRF::GetProperDataRange(): **WARNING** end data bin (" << end << ") > histo length (" << fForward.size() << ").";
    std::cerr << std::endl << ">>    Will set end = (histo length - 1). Consider to change it in the msr-file." << std::endl;
    std::cerr << std::endl;
    end = static_cast<Int_t>(fForward.size()-1);
  }
 
  // keep good bins for potential later use
  fGoodBins[0] = start;
  fGoodBins[1] = end;
 
  // make sure that fGoodBins are in proper range for fForward
  if (fGoodBins[0] < 0)
    fGoodBins[0]=0;
  if (fGoodBins[1] > fForward.size()) {
    std::cerr << std::endl << ">> PRunSingleHisto::GetProperDataRange **WARNING** needed to shift forward lgb,";
    std::cerr << std::endl << ">> from " << fGoodBins[1] << " to " << fForward.size()-1 << std::endl;
    fGoodBins[1]=fForward.size()-1;
  }
 
  return true;
}
 
//--------------------------------------------------------------------------
// GetProperFitRange (private)
//--------------------------------------------------------------------------
void PRunSingleHistoRRF::GetProperFitRange(PMsrGlobalBlock *globalBlock)
{
  // set fit start/end time; first check RUN Block
  fFitStartTime = fRunInfo->GetFitRange(0);
  fFitEndTime   = fRunInfo->GetFitRange(1);
  // if fit range is given in bins (and not time), the fit start/end time can be calculated at this point now
  if (fRunInfo->IsFitRangeInBin()) {
    fFitStartTime = (fGoodBins[0] + fRunInfo->GetFitRangeOffset(0) - fT0s[0]) * fTimeResolution; // (fgb+n0-t0)*dt
    fFitEndTime = (fGoodBins[1] - fRunInfo->GetFitRangeOffset(1) - fT0s[0]) * fTimeResolution;   // (lgb-n1-t0)*dt
    // write these times back into the data structure. This way it is available when writting the log-file
    fRunInfo->SetFitRange(fFitStartTime, 0);
    fRunInfo->SetFitRange(fFitEndTime, 1);
  }
  if (fFitStartTime == PMUSR_UNDEFINED) { // fit start/end NOT found in the RUN block, check GLOBAL block
    fFitStartTime = globalBlock->GetFitRange(0);
    fFitEndTime   = globalBlock->GetFitRange(1);
    // if fit range is given in bins (and not time), the fit start/end time can be calculated at this point now
    if (globalBlock->IsFitRangeInBin()) {
      fFitStartTime = (fGoodBins[0] + globalBlock->GetFitRangeOffset(0) - fT0s[0]) * fTimeResolution; // (fgb+n0-t0)*dt
      fFitEndTime = (fGoodBins[1] - globalBlock->GetFitRangeOffset(1) - fT0s[0]) * fTimeResolution;   // (lgb-n1-t0)*dt
      // write these times back into the data structure. This way it is available when writting the log-file
      globalBlock->SetFitRange(fFitStartTime, 0);
      globalBlock->SetFitRange(fFitEndTime, 1);
    }
  }
  if ((fFitStartTime == PMUSR_UNDEFINED) || (fFitEndTime == PMUSR_UNDEFINED)) {
    fFitStartTime = (fGoodBins[0] - fT0s[0]) * fTimeResolution; // (fgb-t0)*dt
    fFitEndTime = (fGoodBins[1] - fT0s[0]) * fTimeResolution;   // (lgb-t0)*dt
    std::cerr << ">> PRunSingleHistoRRF::GetProperFitRange(): **WARNING** Couldn't get fit start/end time!" << std::endl;
    std::cerr << ">>    Will set it to fgb/lgb which given in time is: " << fFitStartTime << "..." << fFitEndTime << " (usec)" << std::endl;
  }
}
 
//--------------------------------------------------------------------------
// GetMainFrequency (private)
//--------------------------------------------------------------------------
Double_t PRunSingleHistoRRF::GetMainFrequency(PDoubleVector &data)
{
  Double_t freqMax = 0.0;
 
  // create histo
  Double_t startTime = (fGoodBins[0]-fT0s[0]) * fTimeResolution;
  Int_t noOfBins = fGoodBins[1]-fGoodBins[0]+1;
  std::unique_ptr<TH1F> histo = std::make_unique<TH1F>("data", "data", noOfBins, startTime-fTimeResolution/2.0, startTime+fTimeResolution/2.0+noOfBins*fTimeResolution);
  for (Int_t i=fGoodBins[0]; i<=fGoodBins[1]; i++) {
    histo->SetBinContent(i-fGoodBins[0]+1, data[i]);
  }
 
  // Fourier transform
  std::unique_ptr<PFourier> ft = std::make_unique<PFourier>(histo.get(), FOURIER_UNIT_FREQ);
  ft->Transform(F_APODIZATION_STRONG);
 
  // find frequency maximum
  TH1F *power = ft->GetPowerFourier();
  Double_t maxFreqVal = 0.0;
  for (Int_t i=1; i<power->GetNbinsX(); i++) {
    // ignore dc part at 0 frequency
    if (i<power->GetNbinsX()-1) {
      if (power->GetBinContent(i)>power->GetBinContent(i+1))
        continue;
    }
    // ignore everything below 10 MHz
    if (power->GetBinCenter(i) < 10.0)
      continue;
    // check for maximum
    if (power->GetBinContent(i) > maxFreqVal) {
      maxFreqVal = power->GetBinContent(i);
      freqMax = power->GetBinCenter(i);
    }
  }
 
  // clean up
  if (power)
    delete power;
 
  return freqMax;
}
 
//--------------------------------------------------------------------------
// EstimateN0 (private)
//--------------------------------------------------------------------------
Double_t PRunSingleHistoRRF::EstimateN0(Double_t &errN0, Double_t freqMax)
{
  // endBin is estimated such that the number of full cycles (according to the maximum frequency of the data)
  // is approximately the time fN0EstimateEndTime.
  Int_t endBin = static_cast<Int_t>(round(ceil(fN0EstimateEndTime*freqMax/TMath::TwoPi()) * (TMath::TwoPi()/freqMax) / fTimeResolution));
 
  Double_t n0 = 0.0;
  Double_t wN = 0.0;
  for (Int_t i=0; i<endBin; i++) {
//    n0 += fW[i]*fM[i];
    n0 += fM[i];
    wN += fW[i];
  }
//  n0 /= wN;
  n0 /= endBin;
 
  errN0 = 0.0;
  for (Int_t i=0; i<endBin; i++) {
    errN0 += fW[i]*fW[i]*fMerr[i]*fMerr[i];
  }
  errN0 = sqrt(errN0)/wN;
 
  std::cout << "info> PRunSingleHistoRRF::EstimateN0(): N0=" << n0 << "(" << errN0 << ")" << std::endl;
 
  return n0;
}
 
//--------------------------------------------------------------------------
// EstimatBkg (private)
//--------------------------------------------------------------------------
Bool_t PRunSingleHistoRRF::EstimateBkg(UInt_t histoNo)
{
  Double_t beamPeriod = 0.0;
 
  // check if data are from PSI, RAL, or TRIUMF
  if (fRunInfo->GetInstitute()->Contains("psi"))
    beamPeriod = ACCEL_PERIOD_PSI;
  else if (fRunInfo->GetInstitute()->Contains("ral"))
    beamPeriod = ACCEL_PERIOD_RAL;
  else if (fRunInfo->GetInstitute()->Contains("triumf"))
    beamPeriod = ACCEL_PERIOD_TRIUMF;
  else
    beamPeriod = 0.0;
 
  // check if start and end are in proper order
  UInt_t start = fRunInfo->GetBkgRange(0);
  UInt_t end   = fRunInfo->GetBkgRange(1);
  if (end < start) {
    std::cout << std::endl << "PRunSingleHistoRRF::EstimatBkg(): end = " << end << " > start = " << start << "! Will swap them!";
    UInt_t keep = end;
    end = start;
    start = keep;
  }
 
  // calculate proper background range
  if (beamPeriod != 0.0) {
    Double_t timeBkg = static_cast<Double_t>(end-start)*fTimeResolution; // length of the background intervall in time
    UInt_t fullCycles = static_cast<UInt_t>(timeBkg/beamPeriod); // how many proton beam cylces can be placed within the proposed background intervall
    // correct the end of the background intervall such that the background is as close as possible to a multiple of the proton cylce
    end = start + static_cast<UInt_t>((fullCycles*beamPeriod)/fTimeResolution);
    std::cout << std::endl << "PRunSingleHistoRRF::EstimatBkg(): Background " << start << ", " << end;
    if (end == start)
      end = fRunInfo->GetBkgRange(1);
  }
 
  // check if start is within histogram bounds
  if (start >= fForward.size()) {
    std::cerr << std::endl << ">> PRunSingleHistoRRF::EstimatBkg(): **ERROR** background bin values out of bound!";
    std::cerr << std::endl << ">> histo lengths    = " << fForward.size();
    std::cerr << std::endl << ">> background start = " << start;
    std::cerr << std::endl;
    return false;
  }
 
  // check if end is within histogram bounds
  if (end >= fForward.size()) {
    std::cerr << std::endl << ">> PRunSingleHistoRRF::EstimatBkg(): **ERROR** background bin values out of bound!";
    std::cerr << std::endl << ">> histo lengths  = " << fForward.size();
    std::cerr << std::endl << ">> background end = " << end;
    std::cerr << std::endl;
    return false;
  }
 
  // calculate background
  Double_t bkg    = 0.0;
 
  // forward
  for (UInt_t i=start; i<end; i++)
    bkg += fForward[i];
  bkg /= static_cast<Double_t>(end - start + 1);
 
  fBackground = bkg;  // keep background (per bin)
 
  bkg = 0.0;
  for (UInt_t i=start; i<end; i++)
    bkg += pow(fForward[i]-fBackground, 2.0);
  fBkgErr = sqrt(bkg/(static_cast<Double_t>(end - start)));
 
  std::cout << std::endl << "info> fBackground=" << fBackground << "(" << fBkgErr << ")" << std::endl;
 
  fRunInfo->SetBkgEstimated(fBackground, 0);
 
  return true;
}