From 6ee5d76b35614ba144123d7a992054d194762459 Mon Sep 17 00:00:00 2001 From: Suter Andreas Date: Mon, 21 Dec 2015 16:54:12 +0100 Subject: [PATCH] first implementation of the single histo RRF fit (no plotting yet). Seems to work but the error estimate of the RRF asymmetry seems to be too large. Needs to be checked. --- src/classes/Makefile.am | 2 + src/classes/PFitterFcn.cpp | 2 + src/classes/PMsrHandler.cpp | 9 +- src/classes/PMusr.cpp | 12 +- src/classes/PRunListCollection.cpp | 119 ++- src/classes/PRunSingleHistoRRF.cpp | 1483 ++++++++++++++++++++++++++++ src/include/PRunListCollection.h | 16 +- src/include/PRunSingleHistoRRF.h | 79 ++ src/musrfit.cpp | 30 +- 9 files changed, 1738 insertions(+), 14 deletions(-) create mode 100644 src/classes/PRunSingleHistoRRF.cpp create mode 100644 src/include/PRunSingleHistoRRF.h diff --git a/src/classes/Makefile.am b/src/classes/Makefile.am index f6c13d5d..7e70efb4 100644 --- a/src/classes/Makefile.am +++ b/src/classes/Makefile.am @@ -21,6 +21,7 @@ h_sources = \ ../include/PRunNonMusr.h \ ../include/PRunMuMinus.h \ ../include/PRunSingleHisto.h \ + ../include/PRunSingleHistoRRF.h \ ../include/PStartupHandler.h \ ../include/PTheory.h @@ -65,6 +66,7 @@ cpp_sources = \ PRunNonMusr.cpp \ PRunMuMinus.cpp \ PRunSingleHisto.cpp \ + PRunSingleHistoRRF.cpp \ PStartupHandler.cpp \ PTheory.cpp diff --git a/src/classes/PFitterFcn.cpp b/src/classes/PFitterFcn.cpp index e67ac505..ccae307c 100644 --- a/src/classes/PFitterFcn.cpp +++ b/src/classes/PFitterFcn.cpp @@ -77,11 +77,13 @@ Double_t PFitterFcn::operator()(const std::vector& par) const if (fUseChi2) { // chi square value += fRunListCollection->GetSingleHistoChisq(par); + value += fRunListCollection->GetSingleHistoRRFChisq(par); value += fRunListCollection->GetAsymmetryChisq(par); value += fRunListCollection->GetMuMinusChisq(par); value += fRunListCollection->GetNonMusrChisq(par); } else { // max likelihood value += fRunListCollection->GetSingleHistoMaximumLikelihood(par); + value += fRunListCollection->GetSingleHistoRRFMaximumLikelihood(par); value += fRunListCollection->GetAsymmetryMaximumLikelihood(par); value += fRunListCollection->GetMuMinusMaximumLikelihood(par); value += fRunListCollection->GetNonMusrMaximumLikelihood(par); diff --git a/src/classes/PMsrHandler.cpp b/src/classes/PMsrHandler.cpp index 284b950b..3a3148d0 100644 --- a/src/classes/PMsrHandler.cpp +++ b/src/classes/PMsrHandler.cpp @@ -636,7 +636,13 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t messages) fout.width(16); fout << left << "rrf_freq "; fout.width(8); - fout << left << fGlobal.GetRRFFreq(fGlobal.GetRRFUnit().Data()); + +cout << "debug> PMsrHandler::WriteMsrLogFile(): fGlobal.GetRRFFreq(fGlobal.GetRRFUnit().Data())=" << fGlobal.GetRRFFreq(fGlobal.GetRRFUnit().Data()) << endl; +neededPrec = LastSignificant(fGlobal.GetRRFFreq(fGlobal.GetRRFUnit().Data()),10); +cout << "debug> PMsrHandler::WriteMsrLogFile(): neededPrec=" << neededPrec << endl; +fout.precision(neededPrec); + + fout << left << fixed << fGlobal.GetRRFFreq(fGlobal.GetRRFUnit().Data()); fout << " " << fGlobal.GetRRFUnit(); fout << endl; } else if (sstr.BeginsWith("rrf_phase", TString::kIgnoreCase) && (fGlobal.GetFitType() == MSR_FITTYPE_SINGLE_HISTO_RRF)) { @@ -1650,6 +1656,7 @@ Int_t PMsrHandler::WriteMsrFile(const Char_t *filename, map *co fout.width(16); fout << left << "rrf_freq "; fout.width(8); +cout << "debug> PMsrHandler::WriteMsrFile(): fGlobal.GetRRFFreq(fGlobal.GetRRFUnit().Data())=" << fGlobal.GetRRFFreq(fGlobal.GetRRFUnit().Data()) << endl; fout << left << fGlobal.GetRRFFreq(fGlobal.GetRRFUnit().Data()); fout << " " << fGlobal.GetRRFUnit(); fout << endl; diff --git a/src/classes/PMusr.cpp b/src/classes/PMusr.cpp index f77015b1..e2491b1d 100644 --- a/src/classes/PMusr.cpp +++ b/src/classes/PMusr.cpp @@ -752,21 +752,25 @@ Double_t PMsrGlobalBlock::GetRRFFreq(const char *unit) return freq; } +cout << endl << "debug> PMsrGlobalBlock::GetRRFFreq(" << unit << "): unitTag=" << unitTag << ", fRRFFreq=" << fRRFFreq; + // calc the conversion factor if (unitTag == fRRFUnitTag) freq = fRRFFreq; else if ((unitTag == RRF_UNIT_MHz) && (fRRFUnitTag == RRF_UNIT_Mcs)) freq = fRRFFreq/TMath::TwoPi(); else if ((unitTag == RRF_UNIT_MHz) && (fRRFUnitTag == RRF_UNIT_T)) - freq = fRRFFreq*GAMMA_BAR_MUON; + freq = fRRFFreq*1e4*GAMMA_BAR_MUON; // 1e4 need for T -> G since GAMMA_BAR_MUON is given in MHz/G else if ((unitTag == RRF_UNIT_Mcs) && (fRRFUnitTag == RRF_UNIT_MHz)) freq = fRRFFreq*TMath::TwoPi(); else if ((unitTag == RRF_UNIT_Mcs) && (fRRFUnitTag == RRF_UNIT_T)) - freq = fRRFFreq*TMath::TwoPi()*GAMMA_BAR_MUON; + freq = fRRFFreq*1e4*TMath::TwoPi()*GAMMA_BAR_MUON; // 1e4 need for T -> G since GAMMA_BAR_MUON is given in MHz/G else if ((unitTag == RRF_UNIT_T) && (fRRFUnitTag == RRF_UNIT_MHz)) - freq = fRRFFreq/GAMMA_BAR_MUON; + freq = fRRFFreq/GAMMA_BAR_MUON*1e-4; // 1e-4 need for G -> T since GAMMA_BAR_MUON is given in MHz/G else if ((unitTag == RRF_UNIT_T) && (fRRFUnitTag == RRF_UNIT_Mcs)) - freq = fRRFFreq/(TMath::TwoPi()*GAMMA_BAR_MUON); + freq = fRRFFreq/(TMath::TwoPi()*GAMMA_BAR_MUON)*1e-4; // 1e-4 need for G -> T since GAMMA_BAR_MUON is given in MHz/G + +cout << endl << "debug> PMsrGlobalBlock::GetRRFFreq(" << unit << "): freq=" << freq << endl; return freq; } diff --git a/src/classes/PRunListCollection.cpp b/src/classes/PRunListCollection.cpp index 4082e91e..bc58df51 100644 --- a/src/classes/PRunListCollection.cpp +++ b/src/classes/PRunListCollection.cpp @@ -58,6 +58,12 @@ PRunListCollection::~PRunListCollection() } fRunSingleHistoList.clear(); + for (UInt_t i=0; iCleanUp(); + fRunSingleHistoRRFList[i]->~PRunSingleHistoRRF(); + } + fRunSingleHistoRRFList.clear(); + for (UInt_t i=0; iCleanUp(); fRunAsymmetryList[i]->~PRunAsymmetry(); @@ -100,12 +106,20 @@ Bool_t PRunListCollection::Add(Int_t runNo, EPMusrHandleTag tag) fitType = (*fMsrInfo->GetMsrGlobal()).GetFitType(); } +cout << "debug> PRunListCollection::Add(): (runNo: " << runNo << "), fitType = " << fitType << endl; + switch (fitType) { case PRUN_SINGLE_HISTO: fRunSingleHistoList.push_back(new PRunSingleHisto(fMsrInfo, fData, runNo, tag)); if (!fRunSingleHistoList[fRunSingleHistoList.size()-1]->IsValid()) success = false; break; + case PRUN_SINGLE_HISTO_RRF: +cout << "debug> PRunListCollection::Add(): add RRF single histo run to PRunListCollection (runNo: " << runNo << ")" << endl; + fRunSingleHistoRRFList.push_back(new PRunSingleHistoRRF(fMsrInfo, fData, runNo, tag)); + if (!fRunSingleHistoRRFList[fRunSingleHistoRRFList.size()-1]->IsValid()) + success = false; + break; case PRUN_ASYMMETRY: fRunAsymmetryList.push_back(new PRunAsymmetry(fMsrInfo, fData, runNo, tag)); if (!fRunAsymmetryList[fRunAsymmetryList.size()-1]->IsValid()) @@ -147,6 +161,8 @@ void PRunListCollection::SetFitRange(const TString fitRange) { for (UInt_t i=0; iSetFitRangeBin(fitRange); + for (UInt_t i=0; iSetFitRangeBin(fitRange); for (UInt_t i=0; iSetFitRangeBin(fitRange); for (UInt_t i=0; iSetFitRange(fitRange); + for (UInt_t i=0; iSetFitRange(fitRange); for (UInt_t i=0; iSetFitRange(fitRange); for (UInt_t i=0; i& pa return chisq; } +//-------------------------------------------------------------------------- +// GetSingleHistoRRFChisq (public) +//-------------------------------------------------------------------------- +/** + *

Calculates chi-square of all single histogram RRF runs of a msr-file. + * + * return: + * - chi-square of all single histogram RRF runs of the msr-file + * + * \param par fit parameter vector + */ +Double_t PRunListCollection::GetSingleHistoRRFChisq(const std::vector& par) const +{ + Double_t chisq = 0.0; + + for (UInt_t i=0; iCalcChiSquare(par); + + return chisq; +} + //-------------------------------------------------------------------------- // GetAsymmetryChisq (public) //-------------------------------------------------------------------------- @@ -299,6 +338,9 @@ Double_t PRunListCollection::GetSingleHistoChisqExpected(const std::vectorCalcChiSquareExpected(par); break; + case PRUN_SINGLE_HISTO_RRF: + expectedChisq = fRunSingleHistoRRFList[subIdx]->CalcChiSquareExpected(par); + break; case PRUN_ASYMMETRY: expectedChisq = fRunAsymmetryList[subIdx]->CalcChiSquareExpected(par); break; @@ -353,6 +395,9 @@ Double_t PRunListCollection::GetSingleRunChisq(const std::vector& par, case PRUN_SINGLE_HISTO: chisq = fRunSingleHistoList[subIdx]->CalcChiSquare(par); break; + case PRUN_SINGLE_HISTO_RRF: + chisq = fRunSingleHistoRRFList[subIdx]->CalcChiSquare(par); + break; case PRUN_ASYMMETRY: chisq = fRunAsymmetryList[subIdx]->CalcChiSquare(par); break; @@ -376,7 +421,7 @@ Double_t PRunListCollection::GetSingleRunChisq(const std::vector& par, *

Calculates log max-likelihood of all single histogram runs of a msr-file. * * return: - * - chi-square of all single histogram runs of the msr-file + * - log max-likelihood of all single histogram runs of the msr-file * * \param par fit parameter vector */ @@ -390,6 +435,27 @@ Double_t PRunListCollection::GetSingleHistoMaximumLikelihood(const std::vectorCalculates log max-likelihood of all single histogram RRF runs of a msr-file. + * + * return: + * - log max-likelihood of all single histogram runs of the msr-file + * + * \param par fit parameter vector + */ +Double_t PRunListCollection::GetSingleHistoRRFMaximumLikelihood(const std::vector& par) const +{ + Double_t mlh = 0.0; + + for (UInt_t i=0; iCalcMaxLikelihood(par); + + return mlh; +} + //-------------------------------------------------------------------------- // GetAsymmetryMaximumLikelihood (public) //-------------------------------------------------------------------------- @@ -419,7 +485,7 @@ Double_t PRunListCollection::GetAsymmetryMaximumLikelihood(const std::vectorCalculates log max-likelihood of all mu minus runs of a msr-file. * * return: - * - chi-square of all mu minus runs of the msr-file + * - log max-likelihood of all mu minus runs of the msr-file * * \param par fit parameter vector */ @@ -493,6 +559,9 @@ UInt_t PRunListCollection::GetNoOfBinsFitted(const UInt_t idx) const case PRUN_SINGLE_HISTO: result = fRunSingleHistoList[subIdx]->GetNoOfFitBins(); break; + case PRUN_SINGLE_HISTO_RRF: + result = fRunSingleHistoRRFList[subIdx]->GetNoOfFitBins(); + break; case PRUN_ASYMMETRY: result = fRunAsymmetryList[subIdx]->GetNoOfFitBins(); break; @@ -526,6 +595,9 @@ UInt_t PRunListCollection::GetTotalNoOfBinsFitted() const for (UInt_t i=0; iGetNoOfFitBins(); + for (UInt_t i=0; iGetNoOfFitBins(); + for (UInt_t i=0; iGetNoOfFitBins(); @@ -581,6 +653,49 @@ PRunData* PRunListCollection::GetSingleHisto(UInt_t index, EDataSwitch tag) return data; } +//-------------------------------------------------------------------------- +// GetSingleHistoRRF (public) +//-------------------------------------------------------------------------- +/** + *

Get a processed single histogram RRF data set. + * + * return: + * - pointer to the run data set (processed data) if data set is found + * - null pointer otherwise + * + * \param index msr-file run index + * \param tag kIndex -> data at index, kRunNo -> data of given run no + */ +PRunData* PRunListCollection::GetSingleHistoRRF(UInt_t index, EDataSwitch tag) +{ + PRunData *data = 0; + + switch (tag) { + case kIndex: + if ((index < 0) || (index >= fRunSingleHistoRRFList.size())) { + cerr << endl << "PRunListCollection::GetSingleHistoRRF: **ERROR** index = " << index << " out of bounds"; + cerr << endl; + return 0; + } + + fRunSingleHistoRRFList[index]->CalcTheory(); + data = fRunSingleHistoRRFList[index]->GetData(); + break; + case kRunNo: + for (UInt_t i=0; iGetRunNo() == index) { + data = fRunSingleHistoRRFList[i]->GetData(); + break; + } + } + break; + default: // error + break; + } + + return data; +} + //-------------------------------------------------------------------------- // GetAsymmetry (public) //-------------------------------------------------------------------------- diff --git a/src/classes/PRunSingleHistoRRF.cpp b/src/classes/PRunSingleHistoRRF.cpp new file mode 100644 index 00000000..856935af --- /dev/null +++ b/src/classes/PRunSingleHistoRRF.cpp @@ -0,0 +1,1483 @@ +/*************************************************************************** + + PRunSingleHistoRRF.cpp + + Author: Andreas Suter + e-mail: andreas.suter@psi.ch + +***************************************************************************/ + +/*************************************************************************** + * Copyright (C) 2007-2016 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 +#endif + +#include +#include +#include +using namespace std; + +#include +#include +#include + +#include "PMusr.h" +#include "PRunSingleHistoRRF.h" + +//-------------------------------------------------------------------------- +// Constructor +//-------------------------------------------------------------------------- +/** + *

Constructor + */ +PRunSingleHistoRRF::PRunSingleHistoRRF() : PRunBase() +{ + fNoOfFitBins = 0; + fBackground = 0; + fPacking = -1; + + // 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; +} + +//-------------------------------------------------------------------------- +// Constructor +//-------------------------------------------------------------------------- +/** + *

Constructor + * + * \param msrInfo pointer to the msr-file handler + * \param rawData raw run data + * \param runNo number of the run within the msr-file + * \param tag tag showing what shall be done: kFit == fitting, kView == viewing + */ +PRunSingleHistoRRF::PRunSingleHistoRRF(PMsrHandler *msrInfo, PRunDataHandler *rawData, UInt_t runNo, EPMusrHandleTag tag) : PRunBase(msrInfo, rawData, runNo, tag) +{ + fNoOfFitBins = 0; + + fPacking = fRunInfo->GetPacking(); + if (fPacking == -1) { // i.e. packing is NOT given in the RUN-block, it must be given in the GLOBAL-block + fPacking = fMsrInfo->GetMsrGlobal()->GetPacking(); + } + if (fPacking == -1) { // this should NOT happen, somethin is severely wrong + cerr << endl << ">> PRunSingleHistoRRF::PRunSingleHistoRRF(): **SEVERE ERROR**: Couldn't find any packing information!"; + cerr << endl << ">> This is very bad :-(, will quit ..."; + cerr << 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; + + if (!PrepareData()) { + cerr << endl << ">> PRunSingleHistoRRF::PRunSingleHistoRRF(): **SEVERE ERROR**: Couldn't prepare data for fitting!"; + cerr << endl << ">> This is very bad :-(, will quit ..."; + cerr << endl; + fValid = false; + } +} + +//-------------------------------------------------------------------------- +// Destructor +//-------------------------------------------------------------------------- +/** + *

Destructor + */ +PRunSingleHistoRRF::~PRunSingleHistoRRF() +{ + fForward.clear(); +} + +//-------------------------------------------------------------------------- +// CalcChiSquare (public) +//-------------------------------------------------------------------------- +/** + *

Calculate chi-square. + * + * return: + * - chisq value + * + * \param par parameter vector iterated by minuit2 + */ +Double_t PRunSingleHistoRRF::CalcChiSquare(const std::vector& par) +{ + Double_t chisq = 0.0; + Double_t diff = 0.0; + + // calculate functions + for (Int_t i=0; iGetNoOfFuncs(); i++) { + Int_t funcNo = fMsrInfo->GetFuncNo(i); + fFuncValues[i] = fMsrInfo->EvalFunc(funcNo, *fRunInfo->GetMap(), par); + } + + // calculate chi square + Double_t time(1.0); + Int_t i, N(static_cast(fData.GetValue()->size())); + + // In order not to have an IF in the next loop, determine the start and end bins for the fit range now + Int_t startTimeBin = static_cast(ceil((fFitStartTime - fData.GetDataTimeStart())/fData.GetDataTimeStep())); + if (startTimeBin < 0) + startTimeBin = 0; + Int_t endTimeBin = static_cast(floor((fFitEndTime - fData.GetDataTimeStart())/fData.GetDataTimeStep())) + 1; + if (endTimeBin > N) + endTimeBin = N; + + // 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 = (endTimeBin - startTimeBin)/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=startTimeBin; iat(i) - fTheory->Func(time, par, fFuncValues); + chisq += diff*diff / (fData.GetError()->at(i)*fData.GetError()->at(i)); + } + + return chisq; +} + +//-------------------------------------------------------------------------- +// CalcChiSquareExpected (public) +//-------------------------------------------------------------------------- +/** + *

Calculate expected chi-square. + * + * return: + * - chisq value + * + * \param par parameter vector iterated by minuit2 + */ +Double_t PRunSingleHistoRRF::CalcChiSquareExpected(const std::vector& par) +{ + Double_t chisq = 0.0; + Double_t diff = 0.0; + Double_t theo = 0.0; + + // calculate functions + for (Int_t i=0; iGetNoOfFuncs(); i++) { + Int_t funcNo = fMsrInfo->GetFuncNo(i); + fFuncValues[i] = fMsrInfo->EvalFunc(funcNo, *fRunInfo->GetMap(), par); + } + + // calculate chi square + Double_t time(1.0); + Int_t i, N(static_cast(fData.GetValue()->size())); + + // In order not to have an IF in the next loop, determine the start and end bins for the fit range now + Int_t startTimeBin = static_cast(ceil((fFitStartTime - fData.GetDataTimeStart())/fData.GetDataTimeStep())); + if (startTimeBin < 0) + startTimeBin = 0; + Int_t endTimeBin = static_cast(floor((fFitEndTime - fData.GetDataTimeStart())/fData.GetDataTimeStep())) + 1; + if (endTimeBin > N) + endTimeBin = N; + + // 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 = (endTimeBin - startTimeBin)/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=startTimeBin; i < endTimeBin; ++i) { + time = fData.GetDataTimeStart() + (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) +//-------------------------------------------------------------------------- +/** + *

Calculate log maximum-likelihood. See http://pdg.lbl.gov/index.html + * + * return: + * - log maximum-likelihood value + * + * \param par parameter vector iterated by minuit2 + */ +Double_t PRunSingleHistoRRF::CalcMaxLikelihood(const std::vector& par) +{ + // not yet implemented + + return 0.0; +} + +//-------------------------------------------------------------------------- +// CalcTheory (public) +//-------------------------------------------------------------------------- +/** + *

Calculate theory for a given set of fit-parameters. + */ +void PRunSingleHistoRRF::CalcTheory() +{ + // feed the parameter vector + std::vector par; + PMsrParamList *paramList = fMsrInfo->GetMsrParamList(); + for (UInt_t i=0; isize(); i++) + par.push_back((*paramList)[i].fValue); + + // calculate functions + for (Int_t i=0; iGetNoOfFuncs(); i++) { + fFuncValues[i] = fMsrInfo->EvalFunc(fMsrInfo->GetFuncNo(i), *fRunInfo->GetMap(), par); + } + + // 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; iFunc(time, par, fFuncValues)); + } + + // clean up + par.clear(); +} + +//-------------------------------------------------------------------------- +// GetNoOfFitBins (public) +//-------------------------------------------------------------------------- +/** + *

Calculate the number of fitted bins for the current fit range. + * + * return: number of fitted bins. + */ +UInt_t PRunSingleHistoRRF::GetNoOfFitBins() +{ + CalcNoOfFitBins(); + + return fNoOfFitBins; +} + +//-------------------------------------------------------------------------- +// SetFitRangeBin (public) +//-------------------------------------------------------------------------- +/** + *

Allows to change the fit range on the fly. Used in the COMMAND block. + * The syntax of the string is: FIT_RANGE fgb[+n00] lgb[-n01] [fgb[+n10] lgb[-n11] ... fgb[+nN0] lgb[-nN1]]. + * If only one pair of fgb/lgb is given, it is used for all runs in the RUN block section. + * If multiple fgb/lgb's are given, the number N has to be the number of RUN blocks in + * the msr-file. + * + *

nXY are offsets which can be used to shift, limit the fit range. + * + * \param fitRange string containing the necessary information. + */ +void PRunSingleHistoRRF::SetFitRangeBin(const TString fitRange) +{ + TObjArray *tok = 0; + TObjString *ostr = 0; + TString str; + Ssiz_t idx = -1; + Int_t offset = 0; + + tok = fitRange.Tokenize(" \t"); + + if (tok->GetEntries() == 3) { // structure FIT_RANGE fgb+n0 lgb-n1 + // handle fgb+n0 entry + ostr = (TObjString*) tok->At(1); + str = ostr->GetString(); + // 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 + ostr = (TObjString*) tok->At(2); + str = ostr->GetString(); + // 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->GetEntries() > 3) && (tok->GetEntries() % 2 == 1)) { // structure FIT_RANGE fgb[+n00] lgb[-n01] [fgb[+n10] lgb[-n11] ... fgb[+nN0] lgb[-nN1]] + Int_t pos = 2*(fRunNo+1)-1; + + if (pos + 1 >= tok->GetEntries()) { + cerr << endl << ">> PRunSingleHistoRRF::SetFitRangeBin(): **ERROR** invalid FIT_RANGE command found: '" << fitRange << "'"; + cerr << endl << ">> will ignore it. Sorry ..." << endl; + } else { + // handle fgb+n0 entry + ostr = (TObjString*) tok->At(pos); + str = ostr->GetString(); + // 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 + ostr = (TObjString*) tok->At(pos+1); + str = ostr->GetString(); + // 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 + cerr << endl << ">> PRunSingleHistoRRF::SetFitRangeBin(): **ERROR** invalid FIT_RANGE command found: '" << fitRange << "'"; + cerr << endl << ">> will ignore it. Sorry ..." << endl; + } + + // clean up + if (tok) { + delete tok; + } +} + +//-------------------------------------------------------------------------- +// CalcNoOfFitBins (protected) +//-------------------------------------------------------------------------- +/** + *

Calculate the number of fitted bins for the current fit range. + */ +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 + Int_t startTimeBin = static_cast(ceil((fFitStartTime - fData.GetDataTimeStart())/fData.GetDataTimeStep())); + if (startTimeBin < 0) + startTimeBin = 0; + Int_t endTimeBin = static_cast(floor((fFitEndTime - fData.GetDataTimeStart())/fData.GetDataTimeStep())) + 1; + if (endTimeBin > static_cast(fData.GetValue()->size())) + endTimeBin = fData.GetValue()->size(); + + if (endTimeBin > startTimeBin) + fNoOfFitBins = endTimeBin - startTimeBin; + else + fNoOfFitBins = 0; +} + +//-------------------------------------------------------------------------- +// PrepareData (protected) +//-------------------------------------------------------------------------- +/** + *

Prepare data for fitting or viewing. What is already processed at this stage: + * -# get proper raw run data + * -# get all needed forward histograms + * -# get time resolution + * -# get t0's and perform necessary cross checks (e.g. if t0 of msr-file (if present) are consistent with t0 of the data files, etc.) + * -# add runs (if addruns are present) + * -# group histograms (if grouping is present) + * + * return: + * - true if everthing went smooth + * - false, otherwise. + */ +Bool_t PRunSingleHistoRRF::PrepareData() +{ + Bool_t success = true; + + // 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 + cerr << endl << ">> PRunSingleHistoRRF::PrepareData(): **ERROR** Couldn't get run " << fRunInfo->GetRunName()->Data() << "!"; + cerr << endl; + return false; + } + + // collect histogram numbers + PUIntVector histoNo; // histoNo = msr-file forward + redGreen_offset - 1 + for (UInt_t i=0; iGetForwardHistoNoSize(); i++) { + histoNo.push_back(fRunInfo->GetForwardHistoNo(i)); + + if (!runData->IsPresent(histoNo[i])) { + cerr << endl << ">> PRunSingleHistoRRF::PrepareData(): **PANIC ERROR**:"; + cerr << endl << ">> histoNo found = " << histoNo[i] << ", which is NOT present in the data file!?!?"; + cerr << endl << ">> Will quit :-("; + cerr << endl; + histoNo.clear(); + return false; + } + } + + // keep the time resolution in (us) + fTimeResolution = runData->GetTimeResolution()/1.0e3; + cout.precision(10); + cout << endl << ">> PRunSingleHisto::PrepareData(): time resolution=" << fixed << runData->GetTimeResolution() << "(ns)" << 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) + vector forward; + forward.resize(histoNo.size()); // resize to number of groups + for (UInt_t i=0; iGetDataBin(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; + for (UInt_t i=1; iGetRunNameSize(); i++) { + + // get run to be added to the main one + addRunData = fRawData->GetRunData(*fRunInfo->GetRunName(i)); + if (addRunData == 0) { // couldn't get run + cerr << endl << ">> PRunSingleHistoRRF::PrepareData(): **ERROR** Couldn't get addrun " << fRunInfo->GetRunName(i)->Data() << "!"; + cerr << endl; + return false; + } + + // add forward run + UInt_t addRunSize; + for (UInt_t k=0; kGetDataBin(histoNo[k])->size(); + for (UInt_t j=0; jGetDataBin(histoNo[k])->size(); j++) { // loop over the bin indices + // make sure that the index stays in the proper range + if ((j+(Int_t)fAddT0s[i-1][k]-(Int_t)fT0s[k] >= 0) && (j+(Int_t)fAddT0s[i-1][k]-(Int_t)fT0s[k] < addRunSize)) { + forward[k][j] += addRunData->GetDataBin(histoNo[k])->at(j+(Int_t)fAddT0s[i-1][k]-(Int_t)fT0s[k]); + } + } + } + } + } + + // set forward histo data of the first group + fForward.resize(forward[0].size()); + for (UInt_t i=0; iGetDataBin(histoNo[i])->size(); j++) { // loop over the bin indices + // make sure that the index stays within proper range + if ((j+fT0s[i]-fT0s[0] >= 0) && (j+fT0s[i]-fT0s[0] < runData->GetDataBin(histoNo[i])->size())) { + fForward[j] += forward[i][j+(Int_t)fT0s[i]-(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); + + // get the lifetimecorrection flag + Bool_t lifetimecorrection = false; + PMsrPlotList *plot = fMsrInfo->GetMsrPlotList(); + lifetimecorrection = plot->at(0).fLifeTimeCorrection; + + // do the more fit/view specific stuff + if (fHandleTag == kFit) + success = PrepareFitData(runData, histoNo[0]); + else if ((fHandleTag == kView) && !lifetimecorrection) + success = PrepareRawViewData(runData, histoNo[0]); + else if ((fHandleTag == kView) && lifetimecorrection) + success = PrepareViewData(runData, histoNo[0]); + else + success = false; + + // cleanup + histoNo.clear(); + + return success; +} + +//-------------------------------------------------------------------------- +// PrepareFitData (protected) +//-------------------------------------------------------------------------- +/** + *

Take the pre-processed data (i.e. grouping and addrun are preformed) and form the histogram for fitting. + * The following steps are preformed: + * -# get fit start/stop time + * -# check that 'first good data bin', 'last good data bin', and 't0' make any sense + * -# check how the background shall be handled, i.e. fitted, subtracted from background estimate data range, or subtacted from a given fixed background. + * -# packing (i.e rebinning) + * + * return: + * - true, if everything went smooth + * - false, otherwise + * + * \param runData raw run data handler + * \param histoNo forward histogram number + */ +Bool_t PRunSingleHistoRRF::PrepareFitData(PRawRunData* runData, const UInt_t histoNo) +{ + // keep the raw data for the RRF asymmetry error estimate for later + PDoubleVector sqrtNt; + for (UInt_t i=0; iGetBkgFitParamNo() == -1) { // bkg shall **NOT** be fitted + // 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(fT0s[0]*0.1), 0); + fRunInfo->SetBkgRange(static_cast(fT0s[0]*0.6), 1); + cerr << endl << ">> PRunSingleHistoRRF::PrepareFitData(): **WARNING** Neither fix background nor background bins are given!"; + cerr << endl << ">> Will try the following: bkg start = " << fRunInfo->GetBkgRange(0) << ", bkg end = " << fRunInfo->GetBkgRange(1); + cerr << endl << ">> NO WARRANTY THAT THIS MAKES ANY SENSE! Better check ..."; + cerr << endl; + if (!EstimateBkg(histoNo)) + return false; + } + } else { // fixed background given + for (UInt_t i=0; iGetBkgFix(0); + } + } + } + + Int_t t0 = (Int_t)fT0s[0]; + + // 2) N(t) - Nbkg -> exp(+t/tau) [N(t)-Nbkg] + Double_t startTime = fTimeResolution * ((Double_t)fGoodBins[0] - (Double_t)t0); + cout << endl << "debug> PRunSingleHistoRRF::PrepareFitData(): t0=" << t0 << ", fGoodBins[0]=" << fGoodBins[0] << ", fGoodBins[1]=" << fGoodBins[1] << endl; + cout << endl << "debug> PRunSingleHistoRRF::PrepareFitData(): startTime=" << startTime << endl; + cout << endl << "debug> PRunSingleHistoRRF::PrepareFitData(): endTime =" << startTime+fTimeResolution*((Double_t)fGoodBins[1]-(Double_t)t0) << endl; + + Double_t time_tau=0.0; + for (Int_t i=fGoodBins[0]; i A(t) * cos(wRRF t + phiRRF) + PMsrGlobalBlock *globalBlock = fMsrInfo->GetMsrGlobal(); + Double_t wRRF = globalBlock->GetRRFFreq("Mc"); + Double_t phaseRRF = globalBlock->GetRRFPhase()*TMath::TwoPi()/180.0; + cout << "debug> PRunSingleHistoRRF::PrepareFitData(): wRRF =" << wRRF << endl; + cout << "debug> PRunSingleHistoRRF::PrepareFitData(): phaseRRF =" << phaseRRF << endl; + Double_t time = 0.0; + for (Int_t i=fGoodBins[0]; iGetRRFPacking(); + Double_t dval=0.0; + for (Int_t i=fGoodBins[0]; i 1 + if (((i-fGoodBins[0]) % packingRRF == 0) && (i != fGoodBins[0])) { // fill data + dval /= packingRRF; + fData.AppendValue(dval); + // reset dval + dval = 0.0; + } + dval += fForward[i]; + } + } + + // 7) estimate 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)) + for (Int_t i=fGoodBins[0]; iTake the pre-processed data (i.e. grouping and addrun are preformed) and form the histogram for viewing + * without any life time correction. + *

The following steps are preformed: + * -# check if view packing is whished. + * -# check that 'first good data bin', 'last good data bin', and 't0' makes any sense + * -# packing (i.e. rebinnig) + * -# calculate theory + * + * return: + * - true, if everything went smooth + * - false, otherwise. + * + * \param runData raw run data handler + * \param histoNo forward histogram number + */ +Bool_t PRunSingleHistoRRF::PrepareRawViewData(PRawRunData* runData, const UInt_t histoNo) +{ +/* + // check if view_packing is wished + Int_t packing = fPacking; + if (fMsrInfo->GetMsrPlotList()->at(0).fViewPacking > 0) { + packing = fMsrInfo->GetMsrPlotList()->at(0).fViewPacking; + } + + // calculate necessary norms + Double_t dataNorm = 1.0, theoryNorm = 1.0; + if (fScaleN0AndBkg) { + dataNorm = 1.0/ (packing * (fTimeResolution * 1.0e3)); // fTimeResolution us->ns + } else if (!fScaleN0AndBkg && (fMsrInfo->GetMsrPlotList()->at(0).fViewPacking > 0)) { + theoryNorm = (Double_t)fMsrInfo->GetMsrPlotList()->at(0).fViewPacking/(Double_t)fPacking; + } + + // raw data, since PMusrCanvas is doing ranging etc. + // start = the first bin which is a multiple of packing backward from first good data bin + Int_t start = fGoodBins[0] - (fGoodBins[0]/packing)*packing; + // end = last bin starting from start which is a multipl of packing and still within the data + Int_t end = start + ((fForward.size()-start)/packing)*packing; + // check if data range has been provided, and if not try to estimate them + if (start < 0) { + Int_t offset = (Int_t)(10.0e-3/fTimeResolution); + start = ((Int_t)fT0s[0]+offset) - (((Int_t)fT0s[0]+offset)/packing)*packing; + end = start + ((fForward.size()-start)/packing)*packing; + cerr << endl << ">> PRunSingleHistoRRF::PrepareRawViewData(): **WARNING** data range was not provided, will try data range start = " << start << "."; + cerr << endl << ">> NO WARRANTY THAT THIS DOES MAKE ANY SENSE."; + cerr << endl; + } + // check if start, end, and t0 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 > (Int_t)fForward.size())) { + cerr << endl << ">> PRunSingleHistoRRF::PrepareRawViewData(): **ERROR** start data bin doesn't make any sense!"; + cerr << endl; + return false; + } + // 3rd check if end is within proper bounds + if ((end < 0) || (end > (Int_t)fForward.size())) { + cerr << endl << ">> PRunSingleHistoRRF::PrepareRawViewData(): **ERROR** end data bin doesn't make any sense!"; + cerr << endl; + return false; + } + + // everything looks fine, hence fill data set + Int_t t0 = (Int_t)fT0s[0]; + Double_t value = 0.0; + // data start at data_start-t0 + // time shifted so that packing is included correctly, i.e. t0 == t0 after packing + fData.SetDataTimeStart(fTimeResolution*((Double_t)start-(Double_t)t0+(Double_t)(packing-1)/2.0)); + fData.SetDataTimeStep(fTimeResolution*packing); + + for (Int_t i=start; i par; + PMsrParamList *paramList = fMsrInfo->GetMsrParamList(); + for (UInt_t i=0; isize(); i++) + par.push_back((*paramList)[i].fValue); + + // calculate asymmetry + Double_t N0; + // check if norm is a parameter or a function + if (fRunInfo->GetNormParamNo() < MSR_PARAM_FUN_OFFSET) { // norm is a parameter + N0 = par[fRunInfo->GetNormParamNo()-1]; + } else { // norm is a function + // get function number + UInt_t funNo = fRunInfo->GetNormParamNo()-MSR_PARAM_FUN_OFFSET; + // evaluate function + N0 = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par); + } + N0 *= theoryNorm; + + // get tau + Double_t tau; + if (fRunInfo->GetLifetimeParamNo() != -1) + tau = par[fRunInfo->GetLifetimeParamNo()-1]; + else + tau = PMUON_LIFETIME; + + // get background + Double_t bkg; + if (fRunInfo->GetBkgFitParamNo() == -1) { // bkg not fitted + if (fRunInfo->GetBkgFix(0) == PMUSR_UNDEFINED) { // no fixed background given (background interval) + if (fRunInfo->GetBkgRange(0) >= 0) { // background range given + if (!EstimateBkg(histoNo)) + return false; + } else { // no background given to do the job, try estimate + fRunInfo->SetBkgRange(static_cast(fT0s[0]*0.1), 0); + fRunInfo->SetBkgRange(static_cast(fT0s[0]*0.6), 1); + cerr << endl << ">> PRunSingleHistoRRF::PrepareRawViewData(): **WARNING** Neither fix background nor background bins are given!"; + cerr << endl << ">> Will try the following: bkg start = " << fRunInfo->GetBkgRange(0) << ", bkg end = " << fRunInfo->GetBkgRange(1); + cerr << endl << ">> NO WARRANTY THAT THIS MAKES ANY SENSE! Better check ..."; + cerr << endl; + if (!EstimateBkg(histoNo)) + return false; + } + bkg = fBackground; + } else { // fixed bkg given + bkg = fRunInfo->GetBkgFix(0); + } + } else { // bkg fitted + bkg = par[fRunInfo->GetBkgFitParamNo()-1]; + } + bkg *= theoryNorm; + + // calculate functions + for (Int_t i=0; iGetNoOfFuncs(); i++) { + fFuncValues[i] = fMsrInfo->EvalFunc(fMsrInfo->GetFuncNo(i), *fRunInfo->GetMap(), par); + } + + // calculate theory + UInt_t size = fForward.size(); + Double_t factor = 1.0; + if (fData.GetValue()->size() * 10 > fForward.size()) { + size = fData.GetValue()->size() * 10; + factor = (Double_t)fForward.size() / (Double_t)size; + } + Double_t time; + Double_t theoryValue; + fData.SetTheoryTimeStart(fData.GetDataTimeStart()); + fData.SetTheoryTimeStep(fTimeResolution*factor); + for (UInt_t i=0; iFunc(time, par, fFuncValues); + if (fabs(theoryValue) > 1.0e10) { // dirty hack needs to be fixed!! + theoryValue = 0.0; + } + fData.AppendTheoryValue(N0*TMath::Exp(-time/tau)*(1.0+theoryValue)+bkg); + } + + // clean up + par.clear(); +*/ + return true; +} + +//-------------------------------------------------------------------------- +// PrepareViewData (protected) +//-------------------------------------------------------------------------- +/** + *

Take the pre-processed data (i.e. grouping and addrun are preformed) and form the histogram for viewing + * with life time correction, i.e. the exponential decay is removed. + *

The following steps are preformed: + * -# check if view packing is whished. + * -# check that 'first good data bin', 'last good data bin', and 't0' makes any sense + * -# transform data sets (see below). + * -# calculate theory + * + *

Muon life time corrected data: Starting from + * \f[ N(t) = N_0 e^{-t/\tau} [ 1 + A(t) ] + \mathrm{Bkg} \f] + * it follows that + * \f[ A(t) = (-1) + e^{+t/\tau}\, \frac{N(t)-\mathrm{Bkg}}{N_0}. \f] + * For the error estimate only the statistical error of \f$ N(t) \f$ is used, and hence + * \f[ \Delta A(t) = \frac{e^{t/\tau}}{N_0}\,\sqrt{\frac{N(t)}{p}} \f] + * where \f$ p \f$ is the packing, and \f$ N(t) \f$ are the packed data, i.e. + * \f[ N(t_i) = \frac{1}{p}\, \sum_{j=i}^{i+p} n_j \f] + * with \f$ n_j \f$ the raw histogram data bins. + * + * return: + * - true, if everything went smooth + * - false, otherwise + * + * \param runData raw run data handler + * \param histoNo forward histogram number + */ +Bool_t PRunSingleHistoRRF::PrepareViewData(PRawRunData* runData, const UInt_t histoNo) +{ +/* + // check if view_packing is wished. This is a global option for all PLOT blocks! + Int_t packing = fPacking; + if (fMsrInfo->GetMsrPlotList()->at(0).fViewPacking > 0) { + packing = fMsrInfo->GetMsrPlotList()->at(0).fViewPacking; + } + // check if rrf_packing is present. This is a global option for all PLOT blocks, since operated on a single set of data. + if (fMsrInfo->GetMsrPlotList()->at(0).fRRFPacking > 0) { + packing = fMsrInfo->GetMsrPlotList()->at(0).fRRFPacking; + } + + // calculate necessary norms + Double_t dataNorm = 1.0, theoryNorm = 1.0; + if (fScaleN0AndBkg) { + dataNorm = 1.0/ (packing * (fTimeResolution * 1.0e3)); // fTimeResolution us->ns + } else if (!fScaleN0AndBkg && (fMsrInfo->GetMsrPlotList()->at(0).fViewPacking > 0)) { + theoryNorm = (Double_t)fMsrInfo->GetMsrPlotList()->at(0).fViewPacking/(Double_t)fPacking; + } + + // transform raw histo data. This is done the following way (for details see the manual): + // for the single histo fit, just the rebinned raw data are copied + // first get start data, end data, and t0 + Int_t t0 = (Int_t)fT0s[0]; + + // start = the first bin which is a multiple of packing backward from first good data bin + Int_t start = fGoodBins[0] - (fGoodBins[0]/packing)*packing; + // end = last bin starting from start which is a multiple of packing and still within the data + Int_t end = start + ((fForward.size()-start)/packing)*packing; + + // check if data range has been provided, and if not try to estimate them + if (start < 0) { + Int_t offset = (Int_t)(10.0e-3/fTimeResolution); + start = ((Int_t)fT0s[0]+offset) - (((Int_t)fT0s[0]+offset)/packing)*packing; + end = start + ((fForward.size()-start)/packing)*packing; + cerr << endl << ">> PRunSingleHistoRRF::PrepareViewData(): **WARNING** data range was not provided, will try data range start = " << start << "."; + cerr << endl << ">> NO WARRANTY THAT THIS DOES MAKE ANY SENSE."; + cerr << endl; + } + + // check if start, end, and t0 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 > (Int_t)fForward.size())) { + cerr << endl << ">> PRunSingleHistoRRF::PrepareViewData(): **ERROR** start data bin doesn't make any sense!"; + cerr << endl; + return false; + } + // 3rd check if end is within proper bounds + if ((end < 0) || (end > (Int_t)fForward.size())) { + cerr << endl << ">> PRunSingleHistoRRF::PrepareViewData(): **ERROR** end data bin doesn't make any sense!"; + cerr << endl; + return false; + } + + // everything looks fine, hence fill data set + + // feed the parameter vector + std::vector par; + PMsrParamList *paramList = fMsrInfo->GetMsrParamList(); + for (UInt_t i=0; isize(); i++) + par.push_back((*paramList)[i].fValue); + + // calculate asymmetry + Double_t N0; + // check if norm is a parameter or a function + if (fRunInfo->GetNormParamNo() < MSR_PARAM_FUN_OFFSET) { // norm is a parameter + N0 = par[fRunInfo->GetNormParamNo()-1]; + } else { // norm is a function + // get function number + UInt_t funNo = fRunInfo->GetNormParamNo()-MSR_PARAM_FUN_OFFSET; + // evaluate function + N0 = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par); + } + N0 *= theoryNorm; + + // get tau + Double_t tau; + if (fRunInfo->GetLifetimeParamNo() != -1) + tau = par[fRunInfo->GetLifetimeParamNo()-1]; + else + tau = PMUON_LIFETIME; + + // get background + Double_t bkg; + if (fRunInfo->GetBkgFitParamNo() == -1) { // bkg not fitted + if (fRunInfo->GetBkgFix(0) == PMUSR_UNDEFINED) { // no fixed background given (background interval) + if (fRunInfo->GetBkgRange(0) >= 0) { // background range given + if (!EstimateBkg(histoNo)) + return false; + } else { // no background given to do the job, try estimate + fRunInfo->SetBkgRange(static_cast(fT0s[0]*0.1), 0); + fRunInfo->SetBkgRange(static_cast(fT0s[0]*0.6), 1); + cerr << endl << ">> PRunSingleHistoRRF::PrepareViewData(): **WARNING** Neither fix background nor background bins are given!"; + cerr << endl << ">> Will try the following: bkg start = " << fRunInfo->GetBkgRange(0) << ", bkg end = " << fRunInfo->GetBkgRange(1); + cerr << endl << ">> NO WARRANTY THAT THIS MAKES ANY SENSE! Better check ..."; + cerr << endl; + if (!EstimateBkg(histoNo)) + return false; + } + bkg = fBackground; + } else { // fixed bkg given + bkg = fRunInfo->GetBkgFix(0); + } + } else { // bkg fitted + bkg = par[fRunInfo->GetBkgFitParamNo()-1]; + } + bkg *= theoryNorm; + + Double_t value = 0.0; + Double_t expval = 0.0; + Double_t rrf_val = 0.0; + Double_t time = 0.0; + + // data start at data_start-t0 shifted by (pack-1)/2 + fData.SetDataTimeStart(fTimeResolution*((Double_t)start-(Double_t)t0+(Double_t)(packing-1)/2.0)); + fData.SetDataTimeStep(fTimeResolution*packing); + + // data is always normalized to (per nsec!!) + Double_t gammaRRF = 0.0, wRRF = 0.0, phaseRRF = 0.0; + if (fMsrInfo->GetMsrPlotList()->at(0).fRRFFreq == 0.0) { // normal Data representation + for (Int_t i=start; iGetMsrPlotList()->at(0).fRRFUnit) { + case RRF_UNIT_kHz: + gammaRRF = TMath::TwoPi()*1.0e-3; + break; + case RRF_UNIT_MHz: + gammaRRF = TMath::TwoPi(); + break; + case RRF_UNIT_Mcs: + gammaRRF = 1.0; + break; + case RRF_UNIT_G: + gammaRRF = GAMMA_BAR_MUON*TMath::TwoPi(); + break; + case RRF_UNIT_T: + gammaRRF = GAMMA_BAR_MUON*TMath::TwoPi()*1.0e4; + break; + default: + gammaRRF = TMath::TwoPi(); + break; + } + wRRF = gammaRRF * fMsrInfo->GetMsrPlotList()->at(0).fRRFFreq; + phaseRRF = fMsrInfo->GetMsrPlotList()->at(0).fRRFPhase / 180.0 * TMath::Pi(); + + Double_t error = 0.0; + for (Int_t i=start; iGetNoOfFuncs(); i++) { + fFuncValues[i] = fMsrInfo->EvalFunc(fMsrInfo->GetFuncNo(i), *fRunInfo->GetMap(), par); + } + + // calculate theory + Double_t theoryValue; + UInt_t size = fForward.size(); + Double_t factor = 1.0; + UInt_t rebinRRF = 0; + + if (wRRF == 0) { // no RRF + // check if a finer binning for the theory is needed + if (fData.GetValue()->size() * 10 > fForward.size()) { + size = fData.GetValue()->size() * 10; + factor = (Double_t)fForward.size() / (Double_t)size; + } + fData.SetTheoryTimeStart(fData.GetDataTimeStart()); + fData.SetTheoryTimeStep(fTimeResolution*factor); + } else { // RRF + rebinRRF = static_cast((TMath::Pi()/2.0/wRRF)/fTimeResolution); // RRF time resolution / data time resolution + fData.SetTheoryTimeStart(fData.GetDataTimeStart()); + fData.SetTheoryTimeStep(TMath::Pi()/2.0/wRRF/rebinRRF); // = theory time resolution as close as possible to the data time resolution compatible with wRRF + } + + for (UInt_t i=0; iFunc(time, par, fFuncValues); + if (wRRF != 0.0) { + theoryValue *= 2.0*TMath::Cos(wRRF * time + phaseRRF); + } + if (fabs(theoryValue) > 10.0) { // dirty hack needs to be fixed!! + theoryValue = 0.0; + } + fData.AppendTheoryValue(theoryValue); + } + + // if RRF filter the theory with a FIR Kaiser low pass filter + if (wRRF != 0.0) { + // rebin theory to the RRF frequency + if (rebinRRF != 0) { + Double_t dval = 0.0; + PDoubleVector theo; + for (UInt_t i=0; isize(); i++) { + if ((i % rebinRRF == 0) && (i != 0)) { + theo.push_back(dval/rebinRRF); + dval = 0.0; + } + dval += fData.GetTheory()->at(i); + } + fData.SetTheoryTimeStart(fData.GetTheoryTimeStart()+static_cast(rebinRRF-1)*fData.GetTheoryTimeStep()/2.0); + fData.SetTheoryTimeStep(rebinRRF*fData.GetTheoryTimeStep()); + fData.ReplaceTheory(theo); + theo.clear(); + } + + // filter theory + CalculateKaiserFilterCoeff(wRRF, 60.0, 0.2); // w_c = wRRF, A = -20 log_10(delta), Delta w / w_c = (w_s - w_p) / (2 w_c) + FilterTheo(); + } + + // clean up + par.clear(); +*/ + return true; +} + +//-------------------------------------------------------------------------- +// GetProperT0 (private) +//-------------------------------------------------------------------------- +/** + *

Get the proper t0 for the single histogram run. + * -# the t0 vector size = number of detectors (grouping) for forward. + * -# initialize t0's with -1 + * -# fill t0's from RUN block + * -# if t0's are missing (i.e. t0 == -1), try to fill from the GLOBAL block. + * -# if t0's are missing, try t0's from the data file + * -# if t0's are missing, try to estimate them + * + * \param runData pointer to the current RUN block entry from the msr-file + * \param globalBlock pointer to the GLOBLA block entry from the msr-file + * \param histoNo histogram number vector of forward; histoNo = msr-file forward + redGreen_offset - 1 + * + * return: + * - true if everthing went smooth + * - false, otherwise. + */ +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; iGetT0BinSize(); i++) { + fT0s[i] = fRunInfo->GetT0Bin(i); + } + + // fill in the T0's from the GLOBAL block section (if present) + for (UInt_t i=0; iGetT0BinSize(); i++) { + if (fT0s[i] == -1) { // 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; iGetT0Bin(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; iGetT0BinEstimated(histoNo[i]); + fRunInfo->SetT0Bin(fT0s[i], i); // keep value for the msr-file + + cerr << endl << ">> PRunSingleHistoRRF::GetProperT0(): **WARRNING** NO t0's found, neither in the run data nor in the msr-file!"; + cerr << endl << ">> run: " << fRunInfo->GetRunName()->Data(); + cerr << endl << ">> will try the estimated one: forward t0 = " << runData->GetT0BinEstimated(histoNo[i]); + cerr << endl << ">> NO WARRANTY THAT THIS OK!! For instance for LEM this is almost for sure rubbish!"; + cerr << endl; + } + } + + // check if t0 is within proper bounds + for (UInt_t i=0; iGetForwardHistoNoSize(); i++) { + if ((fT0s[i] < 0) || (fT0s[i] > (Int_t)runData->GetDataBin(histoNo[i])->size())) { + cerr << endl << ">> PRunSingleHistoRRF::GetProperT0(): **ERROR** t0 data bin (" << fT0s[i] << ") doesn't make any sense!"; + cerr << 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; iGetRunNameSize(); i++) { + + // get run to be added to the main one + addRunData = fRawData->GetRunData(*fRunInfo->GetRunName(i)); + if (addRunData == 0) { // couldn't get run + cerr << endl << ">> PRunSingleHistoRRF::GetProperT0(): **ERROR** Couldn't get addrun " << fRunInfo->GetRunName(i)->Data() << "!"; + cerr << 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; jGetT0BinSize(); 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; jGetT0Bin(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; jGetT0BinEstimated(histoNo[j]); + fRunInfo->SetAddT0Bin(fAddT0s[i-1][j], i-1, j); // keep value for the msr-file + + cerr << endl << ">> PRunSingleHistoRRF::GetProperT0(): **WARRNING** NO t0's found, neither in the run data nor in the msr-file!"; + cerr << endl << ">> run: " << fRunInfo->GetRunName(i)->Data(); + cerr << endl << ">> will try the estimated one: forward t0 = " << addRunData->GetT0BinEstimated(histoNo[j]); + cerr << endl << ">> NO WARRANTY THAT THIS OK!! For instance for LEM this is almost for sure rubbish!"; + cerr << endl; + } + } + + // check if t0 is within proper bounds + for (UInt_t j=0; jGetForwardHistoNoSize(); j++) { + if ((fAddT0s[i-1][j] < 0) || (fAddT0s[i-1][j] > (Int_t)addRunData->GetDataBin(histoNo[j])->size())) { + cerr << endl << ">> PRunSingleHistoRRF::GetProperT0(): **ERROR** addt0 data bin (" << fAddT0s[i-1][j] << ") doesn't make any sense!"; + cerr << endl; + return false; + } + } + } + } + + return true; +} + +//-------------------------------------------------------------------------- +// GetProperDataRange (private) +//-------------------------------------------------------------------------- +/** + *

Get the proper data range, i.e. first/last good bin (fgb/lgb). + * -# get fgb/lgb from the RUN block + * -# if fgb/lgb still undefined, try to get it from the GLOBAL block + * -# if fgb/lgb still undefined, try to estimate them. + * + * return: + * - true if everthing went smooth + * - false, otherwise. + */ +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 = (Int_t)(10.0e-3/fTimeResolution); + start = (Int_t)fT0s[0]+offset; + fRunInfo->SetDataRange(start, 0); + cerr << endl << ">> PRunSingleHistoRRF::GetProperDataRange(): **WARNING** data range was not provided, will try data range start = t0+" << offset << "(=10ns) = " << start << "."; + cerr << endl << ">> NO WARRANTY THAT THIS DOES MAKE ANY SENSE."; + cerr << endl; + } + if (end < 0) { + end = fForward.size(); + fRunInfo->SetDataRange(end, 1); + cerr << endl << ">> PRunSingleHistoRRF::GetProperDataRange(): **WARNING** data range was not provided, will try data range end = " << end << "."; + cerr << endl << ">> NO WARRANTY THAT THIS DOES MAKE ANY SENSE."; + cerr << 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 > (Int_t)fForward.size())) { + cerr << endl << ">> PRunSingleHistoRRF::GetProperDataRange(): **ERROR** start data bin (" << start << ") doesn't make any sense!"; + cerr << endl; + return false; + } + // 3rd check if end is within proper bounds + if ((end < 0) || (end > (Int_t)fForward.size())) { + cerr << endl << ">> PRunSingleHistoRRF::GetProperDataRange(): **ERROR** end data bin (" << end << ") doesn't make any sense!"; + cerr << endl; + return false; + } + + // keep good bins for potential later use + fGoodBins[0] = start; + fGoodBins[1] = end; + + return true; +} + +//-------------------------------------------------------------------------- +// GetProperFitRange (private) +//-------------------------------------------------------------------------- +/** + *

Get the proper fit range. There are two possible fit range commands: + * fit given in (usec), or + * fit fgb+offset_0 lgb-offset_1 given in (bins), therefore it works the following way: + * -# get fit range assuming given in time from RUN block + * -# if fit range in RUN block is given in bins, replace start/end + * -# if fit range is NOT given yet, try fit range assuming given in time from GLOBAL block + * -# if fit range in GLOBAL block is given in bins, replace start/end + * -# if still no fit range is given, use fgb/lgb. + * + * \param globalBlock pointer to the GLOBAL block information form the msr-file. + */ +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 + cerr << ">> PRunSingleHistoRRF::GetProperFitRange(): **WARNING** Couldn't get fit start/end time!" << endl; + cerr << ">> Will set it to fgb/lgb which given in time is: " << fFitStartTime << "..." << fFitEndTime << " (usec)" << endl; + } +} + +//-------------------------------------------------------------------------- +// EstimateN0 (private) +//-------------------------------------------------------------------------- +/** + *

Estimate the N0 for the given run. + */ +Double_t PRunSingleHistoRRF::EstimateN0() +{ + Int_t endBin = N0EstimateEndTime / fTimeResolution; + Double_t n0 = 0.0; + cout << "debug> PRunSingleHistoRRF::EstimateN0(): startBin=" << fGoodBins[0] << ", endBin=" << endBin << endl; + for (Int_t i=fGoodBins[0]; i PRunSingleHistoRRF::EstimateN0(): N0=" << n0 << endl; + + return n0; +} + +//-------------------------------------------------------------------------- +// EstimatBkg (private) +//-------------------------------------------------------------------------- +/** + *

Estimate the background for a given interval. + * + * return: + * - true, if everything went smooth + * - false, otherwise + * + * \param histoNo forward histogram number of the run + */ +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) { + cout << 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 = (Double_t)(end-start)*(fTimeResolution*fPacking); // length of the background intervall in time + UInt_t fullCycles = (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 + (UInt_t) ((fullCycles*beamPeriod)/(fTimeResolution*fPacking)); + cout << endl << "PRunSingleHistoRRF::EstimatBkg(): Background " << start << ", " << end; + if (end == start) + end = fRunInfo->GetBkgRange(1); + } + + // check if start is within histogram bounds + if ((start < 0) || (start >= fForward.size())) { + cerr << endl << ">> PRunSingleHistoRRF::EstimatBkg(): **ERROR** background bin values out of bound!"; + cerr << endl << ">> histo lengths = " << fForward.size(); + cerr << endl << ">> background start = " << start; + cerr << endl; + return false; + } + + // check if end is within histogram bounds + if ((end < 0) || (end >= fForward.size())) { + cerr << endl << ">> PRunSingleHistoRRF::EstimatBkg(): **ERROR** background bin values out of bound!"; + cerr << endl << ">> histo lengths = " << fForward.size(); + cerr << endl << ">> background end = " << end; + cerr << endl; + return false; + } + + // calculate background + Double_t bkg = 0.0; + + // forward + for (UInt_t i=start; i(end - start + 1); + + fBackground = bkg * fPacking; // keep background (per bin) + + fRunInfo->SetBkgEstimated(fBackground, 0); + + return true; +} diff --git a/src/include/PRunListCollection.h b/src/include/PRunListCollection.h index 8a741416..acc4018a 100644 --- a/src/include/PRunListCollection.h +++ b/src/include/PRunListCollection.h @@ -8,7 +8,7 @@ ***************************************************************************/ /*************************************************************************** - * Copyright (C) 2007-2014 by Andreas Suter * + * Copyright (C) 2007-2016 by Andreas Suter * * andreas.suter@psi.ch * * * * This program is free software; you can redistribute it and/or modify * @@ -37,6 +37,7 @@ using namespace std; #include "PMsrHandler.h" #include "PRunDataHandler.h" #include "PRunSingleHisto.h" +#include "PRunSingleHistoRRF.h" #include "PRunAsymmetry.h" #include "PRunMuMinus.h" #include "PRunNonMusr.h" @@ -58,6 +59,7 @@ class PRunListCollection virtual void SetFitRange(const TString fitRange); virtual Double_t GetSingleHistoChisq(const std::vector& par) const; + virtual Double_t GetSingleHistoRRFChisq(const std::vector& par) const; virtual Double_t GetAsymmetryChisq(const std::vector& par) const; virtual Double_t GetMuMinusChisq(const std::vector& par) const; virtual Double_t GetNonMusrChisq(const std::vector& par) const; @@ -66,6 +68,7 @@ class PRunListCollection virtual Double_t GetSingleRunChisq(const std::vector& par, const UInt_t idx) const; virtual Double_t GetSingleHistoMaximumLikelihood(const std::vector& par) const; + virtual Double_t GetSingleHistoRRFMaximumLikelihood(const std::vector& par) const; virtual Double_t GetAsymmetryMaximumLikelihood(const std::vector& par) const; virtual Double_t GetMuMinusMaximumLikelihood(const std::vector& par) const; virtual Double_t GetNonMusrMaximumLikelihood(const std::vector& par) const; @@ -74,11 +77,13 @@ class PRunListCollection virtual UInt_t GetTotalNoOfBinsFitted() const; virtual UInt_t GetNoOfSingleHisto() const { return fRunSingleHistoList.size(); } ///< returns the number of single histogram data sets present in the msr-file + virtual UInt_t GetNoOfSingleHistoRRF() const { return fRunSingleHistoRRFList.size(); } ///< returns the number of single histogram RRF data sets present in the msr-file virtual UInt_t GetNoOfAsymmetry() const { return fRunAsymmetryList.size(); } ///< returns the number of asymmetry data sets present in the msr-file virtual UInt_t GetNoOfMuMinus() const { return fRunMuMinusList.size(); } ///< returns the number of mu minus data sets present in the msr-file virtual UInt_t GetNoOfNonMusr() const { return fRunNonMusrList.size(); } ///< returns the number of non-muSR data sets present in the msr-file virtual PRunData* GetSingleHisto(UInt_t index, EDataSwitch tag=kIndex); + virtual PRunData* GetSingleHistoRRF(UInt_t index, EDataSwitch tag=kIndex); virtual PRunData* GetAsymmetry(UInt_t index, EDataSwitch tag=kIndex); virtual PRunData* GetMuMinus(UInt_t index, EDataSwitch tag=kIndex); virtual PRunData* GetNonMusr(UInt_t index, EDataSwitch tag=kIndex); @@ -94,10 +99,11 @@ class PRunListCollection PMsrHandler *fMsrInfo; ///< pointer to the msr-file handler PRunDataHandler *fData; ///< pointer to the run-data handler - vector fRunSingleHistoList; ///< stores all processed single histogram data - vector fRunAsymmetryList; ///< stores all processed asymmetry data - vector fRunMuMinusList; ///< stores all processed mu-minus data - vector fRunNonMusrList; ///< stores all processed non-muSR data + vector fRunSingleHistoList; ///< stores all processed single histogram data + vector fRunSingleHistoRRFList; ///< stores all processed single histogram RRF data + vector fRunAsymmetryList; ///< stores all processed asymmetry data + vector fRunMuMinusList; ///< stores all processed mu-minus data + vector fRunNonMusrList; ///< stores all processed non-muSR data }; #endif // _PRUNLISTCOLLECTION_H_ diff --git a/src/include/PRunSingleHistoRRF.h b/src/include/PRunSingleHistoRRF.h new file mode 100644 index 00000000..17980533 --- /dev/null +++ b/src/include/PRunSingleHistoRRF.h @@ -0,0 +1,79 @@ +/*************************************************************************** + + PRunSingleHistoRRF.h + + Author: Andreas Suter + e-mail: andreas.suter@psi.ch + +***************************************************************************/ + +/*************************************************************************** + * Copyright (C) 2007-2016 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. * + ***************************************************************************/ + +#ifndef _PRUNSINGLEHISTORRF_H_ +#define _PRUNSINGLEHISTORRF_H_ + +#include "PRunBase.h" + +/** + *

Class handling single histogram fit type. + */ +class PRunSingleHistoRRF : public PRunBase +{ + public: + PRunSingleHistoRRF(); + PRunSingleHistoRRF(PMsrHandler *msrInfo, PRunDataHandler *rawData, UInt_t runNo, EPMusrHandleTag tag); + virtual ~PRunSingleHistoRRF(); + + virtual Double_t CalcChiSquare(const std::vector& par); + virtual Double_t CalcChiSquareExpected(const std::vector& par); + virtual Double_t CalcMaxLikelihood(const std::vector& par); + virtual void CalcTheory(); + + virtual UInt_t GetNoOfFitBins(); + + virtual void SetFitRangeBin(const TString fitRange); + + protected: + virtual void CalcNoOfFitBins(); + virtual Bool_t PrepareData(); + virtual Bool_t PrepareFitData(PRawRunData* runData, const UInt_t histoNo); + virtual Bool_t PrepareRawViewData(PRawRunData* runData, const UInt_t histoNo); + virtual Bool_t PrepareViewData(PRawRunData* runData, const UInt_t histoNo); + + private: + static const Double_t N0EstimateEndTime = 2.0; ///< end time in (us) over which N0 is estimated. Should evetually go into the musrfit_startup.xml + + UInt_t fNoOfFitBins; ///< number of bins to be fitted + Double_t fBackground; ///< needed if background range is given (units: 1/bin) + Int_t fPacking; ///< packing for this particular run. Either given in the RUN- or GLOBAL-block. + + Int_t fGoodBins[2]; ///< keep first/last good bins. 0=fgb, 1=lgb + + PDoubleVector fForward; ///< forward histo data + + virtual Bool_t GetProperT0(PRawRunData* runData, PMsrGlobalBlock *globalBlock, PUIntVector &histoNo); + virtual Bool_t GetProperDataRange(); + virtual void GetProperFitRange(PMsrGlobalBlock *globalBlock); + virtual Double_t EstimateN0(); + virtual Bool_t EstimateBkg(UInt_t histoNo); +}; + +#endif // _PRUNSINGLEHISTORRF_H_ diff --git a/src/musrfit.cpp b/src/musrfit.cpp index 060534f0..13e72c18 100644 --- a/src/musrfit.cpp +++ b/src/musrfit.cpp @@ -175,10 +175,10 @@ void musrfit_dump_ascii(char *fileName, PRunListCollection *runList) // go through the run list, get the data and dump it in a file int runCounter = 0; + PRunData *data; // single histos unsigned int size = runList->GetNoOfSingleHisto(); - PRunData *data; if (size > 0) { for (unsigned int i=0; iGetSingleHisto(i); @@ -190,6 +190,19 @@ void musrfit_dump_ascii(char *fileName, PRunListCollection *runList) } } + // single histos + size = runList->GetNoOfSingleHistoRRF(); + if (size > 0) { + for (unsigned int i=0; iGetSingleHistoRRF(i); + if (data) { + // dump data + musrfit_write_ascii(fln, data, runCounter); + runCounter++; + } + } + } + // asymmetry size = runList->GetNoOfAsymmetry(); if (size > 0) { @@ -308,10 +321,10 @@ void musrfit_dump_root(char *fileName, PRunListCollection *runList) // go through the run list, get the data and dump it in a file int runCounter = 0; + PRunData *data; // single histos unsigned int size = runList->GetNoOfSingleHisto(); - PRunData *data; if (size > 0) { for (unsigned int i=0; iGetSingleHisto(i); @@ -323,6 +336,19 @@ void musrfit_dump_root(char *fileName, PRunListCollection *runList) } } + // single histo RRF + size = runList->GetNoOfSingleHistoRRF(); + if (size > 0) { + for (unsigned int i=0; iGetSingleHistoRRF(i); + if (data) { + // dump data + musrfit_write_root(f, fln, data, runCounter); + runCounter++; + } + } + } + // asymmetry size = runList->GetNoOfAsymmetry(); if (size > 0) {