diff --git a/src/classes/CMakeLists.txt b/src/classes/CMakeLists.txt index 358f9c9a..25e99c8a 100644 --- a/src/classes/CMakeLists.txt +++ b/src/classes/CMakeLists.txt @@ -71,6 +71,7 @@ add_library(PMusr SHARED PMusrT0Dict.cxx PPrepFourier.cpp PRunAsymmetry.cpp + PRunAsymmetryBNMR.cpp PRunAsymmetryRRF.cpp PRunBase.cpp PRunDataHandler.cpp diff --git a/src/classes/Makefile.am b/src/classes/Makefile.am new file mode 100644 index 00000000..fa3a0fba --- /dev/null +++ b/src/classes/Makefile.am @@ -0,0 +1,122 @@ +## Process this file with automake to create Makefile.in + +h_sources = \ + ../include/PFitterFcn.h \ + ../include/PFitter.h \ + ../include/PFourier.h \ + ../include/PFourierCanvas.h \ + ../include/PFunctionGrammar.h \ + ../include/PFunction.h \ + ../include/PFunctionHandler.h \ + ../include/PMsr2Data.h \ + ../include/PMsrHandler.h \ + ../include/PMusrCanvas.h \ + ../include/PMusr.h \ + ../include/PMusrT0.h \ + ../include/PPrepFourier.h \ + ../include/PRunAsymmetry.h \ + ../include/PRunAsymmetryBNMR.h \ + ../include/PRunAsymmetryRRF.h \ + ../include/PRunBase.h \ + ../include/PRunDataHandler.h \ + ../include/PRunListCollection.h \ + ../include/PRunNonMusr.h \ + ../include/PRunMuMinus.h \ + ../include/PRunSingleHisto.h \ + ../include/PRunSingleHistoRRF.h \ + ../include/PStartupHandler.h \ + ../include/PTheory.h + +h_sources_userFcn = \ + ../include/PUserFcnBase.h + +h_linkdef = \ + ../include/PFourierCanvasLinkDef.h \ + ../include/PMusrCanvasLinkDef.h \ + ../include/PMusrT0LinkDef.h \ + ../include/PStartupHandlerLinkDef.h + +h_linkdef_userFcn = \ + ../include/PUserFcnBaseLinkDef.h + +dict_h_sources = \ + PFourierCanvasDict.h \ + PMusrCanvasDict.h \ + PMusrT0Dict.h \ + PStartupHandlerDict.h + +dict_h_sources_userFcn = \ + PUserFcnBaseDict.h + +cpp_sources = \ + PFitter.cpp \ + PFitterFcn.cpp \ + PFourier.cpp \ + PFourierCanvas.cpp \ + PFunction.cpp \ + PFunctionHandler.cpp \ + PMsr2Data.cpp \ + PMsrHandler.cpp \ + PMusrCanvas.cpp \ + PMusr.cpp \ + PMusrT0.cpp \ + PPrepFourier.cpp \ + PRunAsymmetry.cpp \ + PRunAsymmetryBNMR.cpp \ + PRunAsymmetryRRF.cpp \ + PRunBase.cpp \ + PRunDataHandler.cpp \ + PRunListCollection.cpp \ + PRunNonMusr.cpp \ + PRunMuMinus.cpp \ + PRunSingleHisto.cpp \ + PRunSingleHistoRRF.cpp \ + PStartupHandler.cpp \ + PTheory.cpp + +cpp_sources_userFcn = \ + PUserFcnBase.cpp + +dict_cpp_sources = \ + PFourierCanvasDict.cpp \ + PMusrCanvasDict.cpp \ + PMusrT0Dict.cpp \ + PStartupHandlerDict.cpp + +dict_cpp_sources_userFcn = \ + PUserFcnBaseDict.cpp + +pcmdir = $(libdir) +pcm_DATA = \ + PFourierCanvasDict_rdict.pcm \ + PMusrCanvasDict_rdict.pcm \ + PMusrT0Dict_rdict.pcm \ + PStartupHandlerDict_rdict.pcm \ + PUserFcnBaseDict_rdict.pcm + +include_HEADERS = $(h_sources) $(h_sources_userFcn) +noinst_HEADERS = $(h_linkdef) $(dict_h_sources) $(h_linkdef_userFcn) $(dict_h_sources_userFcn) + +AM_CPPFLAGS = -I$(top_srcdir)/src/include $(MUSR_ROOT_CFLAGS) $(PSIBIN_CFLAGS) $(MUD_CFLAGS) $(LEM_CFLAGS) $(FFTW3_CFLAGS) $(GSL_CFLAGS) $(BOOST_CFLAGS) -I$(ROOTINCDIR) $(PNEXUS_CXXFLAGS) $(NEXUS_CFLAGS) +AM_CXXFLAGS = $(LOCAL_LIB_CXXFLAGS) + +BUILT_SOURCES = $(dict_cpp_sources) $(dict_h_sources) $(dict_cpp_sources_userFcn) $(dict_h_sources_userFcn) +AM_LDFLAGS = $(LOCAL_LIB_LDFLAGS) -L@ROOTLIBDIR@ +CLEANFILES = *Dict.cpp *Dict.h *Dict* *~ core + +%Dict.cpp %Dict.h: ../include/%.h ../include/%LinkDef.h + @ROOTCLING@ -v -f $*Dict.cpp -c -p $(AM_CPPFLAGS) $^ + +lib_LTLIBRARIES = libPUserFcnBase.la libPMusr.la + +libPUserFcnBase_la_SOURCES = $(h_sources_userFcn) $(cpp_sources_userFcn) $(dict_h_sources_userFcn) $(dict_cpp_sources_userFcn) +libPUserFcnBase_la_LIBADD = $(ROOT_LIBS) +libPUserFcnBase_la_LDFLAGS = -version-info $(MUSR_LIBRARY_VERSION) -release $(MUSR_RELEASE) $(AM_LDFLAGS) + +libPMusr_la_SOURCES = $(h_sources) $(cpp_sources) $(dict_h_sources) $(dict_cpp_sources) +libPMusr_la_LIBADD = libPUserFcnBase.la $(MUSR_ROOT_LIBS) $(LEM_LIBS) $(PSIBIN_LIBS) $(MUD_LIBS) $(PNEXUS_LIBS) $(FFTW3_LIBS) $(GSL_LIBS) $(ROOT_LIBS) +libPMusr_la_LDFLAGS = -version-info $(MUSR_LIBRARY_VERSION) -release $(MUSR_RELEASE) $(AM_LDFLAGS) + +pkgconfigdir = $(libdir)/pkgconfig +pkgconfig_DATA = PUserFcnBase.pc PMusr.pc + diff --git a/src/classes/PRunAsymmetryBNMR.cpp b/src/classes/PRunAsymmetryBNMR.cpp new file mode 100644 index 00000000..b55052a0 --- /dev/null +++ b/src/classes/PRunAsymmetryBNMR.cpp @@ -0,0 +1,1869 @@ +/*************************************************************************** + + PRunAsymmetryBNMR.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 +using namespace std; + +#include +#include +#include + +#include "PMusr.h" +#include "PRunAsymmetryBNMR.h" + +//-------------------------------------------------------------------------- +// Constructor +//-------------------------------------------------------------------------- +/** + *

Constructor + */ +PRunAsymmetryBNMR::PRunAsymmetryBNMR() : PRunBase() +{ + fNoOfFitBins = 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; + + fStartTimeBin = -1; + fEndTimeBin = -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 + */ +PRunAsymmetryBNMR::PRunAsymmetryBNMR(PMsrHandler *msrInfo, PRunDataHandler *rawData, UInt_t runNo, EPMusrHandleTag tag) : PRunBase(msrInfo, rawData, runNo, tag) +{ + // 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; + + fStartTimeBin = -1; + fEndTimeBin = -1; + + 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 << ">> PRunAsymmetryBNMR::PRunAsymmetryBNMR(): **SEVERE ERROR**: Couldn't find any packing information!"; + cerr << endl << ">> This is very bad :-(, will quit ..."; + cerr << endl; + fValid = false; + return; + } + + // check if alpha and/or beta is fixed -------------------- + + PMsrParamList *param = msrInfo->GetMsrParamList(); + + // check if alpha is given + if (fRunInfo->GetAlphaParamNo() == -1) { // no alpha given + cerr << endl << ">> PRunAsymmetryBNMR::PRunAsymmetryBNMR(): **ERROR** no alpha parameter given! This is needed for an asymmetry fit!"; + cerr << endl; + fValid = false; + return; + } + // check if alpha parameter is within proper bounds + if ((fRunInfo->GetAlphaParamNo() < 0) || (fRunInfo->GetAlphaParamNo() > (Int_t)param->size())) { + cerr << endl << ">> PRunAsymmetryBNMR::PRunAsymmetryBNMR(): **ERROR** alpha parameter no = " << fRunInfo->GetAlphaParamNo(); + cerr << endl << ">> This is out of bound, since there are only " << param->size() << " parameters."; + cerr << endl; + fValid = false; + return; + } + // check if alpha is fixed + Bool_t alphaFixedToOne = false; + if (((*param)[fRunInfo->GetAlphaParamNo()-1].fStep == 0.0) && + ((*param)[fRunInfo->GetAlphaParamNo()-1].fValue == 1.0)) + alphaFixedToOne = true; + + // check if beta is given + Bool_t betaFixedToOne = false; + if (fRunInfo->GetBetaParamNo() == -1) { // no beta given hence assuming beta == 1 + betaFixedToOne = true; + } else if ((fRunInfo->GetBetaParamNo() < 0) || (fRunInfo->GetBetaParamNo() > (Int_t)param->size())) { // check if beta parameter is within proper bounds + cerr << endl << ">> PRunAsymmetryBNMR::PRunAsymmetryBNMR(): **ERROR** beta parameter no = " << fRunInfo->GetBetaParamNo(); + cerr << endl << ">> This is out of bound, since there are only " << param->size() << " parameters."; + cerr << endl; + fValid = false; + return; + } else { // check if beta is fixed + if (((*param)[fRunInfo->GetBetaParamNo()-1].fStep == 0.0) && + ((*param)[fRunInfo->GetBetaParamNo()-1].fValue == 1.0)) + betaFixedToOne = true; + } + + // set fAlphaBetaTag + if (alphaFixedToOne && betaFixedToOne) // alpha == 1, beta == 1 + fAlphaBetaTag = 1; + else if (!alphaFixedToOne && betaFixedToOne) // alpha != 1, beta == 1 + fAlphaBetaTag = 2; + else if (alphaFixedToOne && !betaFixedToOne) // alpha == 1, beta != 1 + fAlphaBetaTag = 3; + else + fAlphaBetaTag = 4; + + // calculate fData + if (!PrepareData()) + fValid = false; +} + +//-------------------------------------------------------------------------- +// Destructor +//-------------------------------------------------------------------------- +/** + *

Destructor. + */ +PRunAsymmetryBNMR::~PRunAsymmetryBNMR() +{ + fForward.clear(); + fForwardErr.clear(); + fBackward.clear(); + fBackwardErr.clear(); +} + +//-------------------------------------------------------------------------- +// CalcChiSquare (public) +//-------------------------------------------------------------------------- +/** + *

Calculate chi-square. + * + * return: + * - chisq value + * + * \param par parameter vector iterated by minuit2 + */ +Double_t PRunAsymmetryBNMR::CalcChiSquare(const std::vector& par) +{ + Double_t chisq = 0.0; + Double_t diff = 0.0; + Double_t asymFcnValue = 0.0; + Double_t a, b, f; + + // calculate functions + for (Int_t i=0; iGetNoOfFuncs(); i++) { + fFuncValues[i] = fMsrInfo->EvalFunc(fMsrInfo->GetFuncNo(i), *fRunInfo->GetMap(), par); + } + + // 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. + asymFcnValue = 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,asymFcnValue,a,b,f) schedule(dynamic,chunk) reduction(+:chisq) + #endif + for (i=fStartTimeBin; iFunc(time, par, fFuncValues); + break; + case 2: // alpha != 1, beta == 1 + a = par[fRunInfo->GetAlphaParamNo()-1]; + f = fTheory->Func(time, par, fFuncValues); + asymFcnValue = (f*(a+1.0)-(a-1.0))/((a+1.0)-f*(a-1.0)); + break; + case 3: // alpha == 1, beta != 1 + b = par[fRunInfo->GetBetaParamNo()-1]; + f = fTheory->Func(time, par, fFuncValues); + asymFcnValue = f*(b+1.0)/(2.0-f*(b-1.0)); + break; + case 4: // alpha != 1, beta != 1 + a = par[fRunInfo->GetAlphaParamNo()-1]; + b = par[fRunInfo->GetBetaParamNo()-1]; + f = fTheory->Func(time, par, fFuncValues); + asymFcnValue = (f*(a*b+1.0)-(a-1.0))/((a+1.0)-f*(a*b-1.0)); + break; + default: + asymFcnValue = 0.0; + break; + } + diff = fData.GetValue()->at(i) - asymFcnValue; + chisq += diff*diff / (fData.GetError()->at(i)*fData.GetError()->at(i)); + } + + return chisq; +} + +//-------------------------------------------------------------------------- +// CalcChiSquareExpected (public) +//-------------------------------------------------------------------------- +/** + *

Calculate expected chi-square. Currently not implemented since not clear what to be done. + * + * return: + * - chisq value == 0.0 + * + * \param par parameter vector iterated by minuit2 + */ +Double_t PRunAsymmetryBNMR::CalcChiSquareExpected(const std::vector& par) +{ + return 0.0; +} + +//-------------------------------------------------------------------------- +// CalcMaxLikelihood (public) +//-------------------------------------------------------------------------- +/** + *

NOT IMPLEMENTED!! + * + * \param par parameter vector iterated by minuit2 + */ +Double_t PRunAsymmetryBNMR::CalcMaxLikelihood(const std::vector& par) +{ + cout << endl << "PRunAsymmetryBNMR::CalcMaxLikelihood(): not implemented yet ..." << endl; + + return 1.0; +} + +//-------------------------------------------------------------------------- +// GetNoOfFitBins (public) +//-------------------------------------------------------------------------- +/** + *

Calculate the number of fitted bins for the current fit range. + * + * return: number of fitted bins. + */ +UInt_t PRunAsymmetryBNMR::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 PRunAsymmetryBNMR::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 << ">> PRunAsymmetryBNMR::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 << ">> PRunAsymmetryBNMR::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 PRunAsymmetryBNMR::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(ceil((fFitStartTime - fData.GetDataTimeStart())/fData.GetDataTimeStep())); + if (fStartTimeBin < 0) + fStartTimeBin = 0; + fEndTimeBin = static_cast(floor((fFitEndTime - fData.GetDataTimeStart())/fData.GetDataTimeStep())) + 1; + if (fEndTimeBin > static_cast(fData.GetValue()->size())) + fEndTimeBin = fData.GetValue()->size(); + + if (fEndTimeBin > fStartTimeBin) + fNoOfFitBins = fEndTimeBin - fStartTimeBin; + else + fNoOfFitBins = 0; +} + +//-------------------------------------------------------------------------- +// CalcTheory (protected) +//-------------------------------------------------------------------------- +/** + *

Calculate theory for a given set of fit-parameters. + */ +void PRunAsymmetryBNMR::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 asymmetry + Double_t asymFcnValue = 0.0; + Double_t a, b, f; + Double_t time; + for (UInt_t i=0; isize(); i++) { + time = fData.GetDataTimeStart() + (Double_t)i*fData.GetDataTimeStep(); + switch (fAlphaBetaTag) { + case 1: // alpha == 1, beta == 1 + asymFcnValue = fTheory->Func(time, par, fFuncValues); + break; + case 2: // alpha != 1, beta == 1 + a = par[fRunInfo->GetAlphaParamNo()-1]; + f = fTheory->Func(time, par, fFuncValues); + asymFcnValue = (f*(a+1.0)-(a-1.0))/((a+1.0)-f*(a-1.0)); + break; + case 3: // alpha == 1, beta != 1 + b = par[fRunInfo->GetBetaParamNo()-1]; + f = fTheory->Func(time, par, fFuncValues); + asymFcnValue = f*(b+1.0)/(2.0-f*(b-1.0)); + break; + case 4: // alpha != 1, beta != 1 + a = par[fRunInfo->GetAlphaParamNo()-1]; + b = par[fRunInfo->GetBetaParamNo()-1]; + f = fTheory->Func(time, par, fFuncValues); + asymFcnValue = (f*(a*b+1.0)-(a-1.0))/((a+1.0)-f*(a*b-1.0)); + break; + default: + asymFcnValue = 0.0; + break; + } + fData.AppendTheoryValue(asymFcnValue); + } + + // clean up + par.clear(); +} + +//-------------------------------------------------------------------------- +// PrepareData (protected) +//-------------------------------------------------------------------------- +/** + *

Prepare data for fitting or viewing. What is already processed at this stage: + * - get all needed forward/backward histograms + * - get time resolution + * - get start/stop fit time + * - 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) + * - subtract background + * + * Error propagation for \f$ A_i = (f_i^{\rm c}-b_i^{\rm c})/(f_i^{\rm c}+b_i^{\rm c})\f$: + * \f[ \Delta A_i = \pm\frac{2}{(f_i^{\rm c}+b_i^{\rm c})^2}\left[ + * (b_i^{\rm c})^2 (\Delta f_i^{\rm c})^2 + + * (\Delta b_i^{\rm c})^2 (f_i^{\rm c})^2\right]^{1/2}\f] + * + * return: + * - true if everthing went smooth + * - false, otherwise. + */ +Bool_t PRunAsymmetryBNMR::PrepareData() +{ + // keep the Global block info + PMsrGlobalBlock *globalBlock = fMsrInfo->GetMsrGlobal(); + + // get forward/backward histo from PRunDataHandler object ------------------------ + // get the correct run + PRawRunData *runData = fRawData->GetRunData(*(fRunInfo->GetRunName())); + if (!runData) { // run not found + cerr << endl << ">> PRunAsymmetryBNMR::PrepareData(): **ERROR** Couldn't get run " << fRunInfo->GetRunName()->Data() << "!"; + cerr << endl; + return false; + } + + // collect histogram numbers + PUIntVector forwardHistoNo; + PUIntVector backwardHistoNo; + for (UInt_t i=0; iGetForwardHistoNoSize(); i++) { + forwardHistoNo.push_back(fRunInfo->GetForwardHistoNo(i)); + + if (!runData->IsPresent(forwardHistoNo[i])) { + cerr << endl << ">> PRunAsymmetryBNMR::PrepareData(): **PANIC ERROR**:"; + cerr << endl << ">> forwardHistoNo found = " << forwardHistoNo[i] << ", which is NOT present in the data file!?!?"; + cerr << endl << ">> Will quit :-("; + cerr << endl; + // clean up + forwardHistoNo.clear(); + backwardHistoNo.clear(); + return false; + } + } + for (UInt_t i=0; iGetBackwardHistoNoSize(); i++) { + backwardHistoNo.push_back(fRunInfo->GetBackwardHistoNo(i)); + + if (!runData->IsPresent(backwardHistoNo[i])) { + cerr << endl << ">> PRunAsymmetryBNMR::PrepareData(): **PANIC ERROR**:"; + cerr << endl << ">> backwardHistoNo found = " << backwardHistoNo[i] << ", which is NOT present in the data file!?!?"; + cerr << endl << ">> Will quit :-("; + cerr << endl; + // clean up + forwardHistoNo.clear(); + backwardHistoNo.clear(); + return false; + } + } + if (forwardHistoNo.size() != backwardHistoNo.size()) { + cerr << endl << ">> PRunAsymmetryBNMR::PrepareData(): **PANIC ERROR**:"; + cerr << endl << ">> # of forward histograms different from # of backward histograms."; + cerr << endl << ">> Will quit :-("; + cerr << endl; + // clean up + forwardHistoNo.clear(); + backwardHistoNo.clear(); + return false; + } + + // keep the time resolution in (us) + fTimeResolution = runData->GetTimeResolution()/1.0e3; + cout.precision(10); + cout << endl << ">> PRunAsymmetryBNMR::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, forwardHistoNo, backwardHistoNo)) { + return false; + } + + // keep the histo of each group at this point (addruns handled below) + vector forward, backward; + forward.resize(forwardHistoNo.size()); // resize to number of groups + backward.resize(backwardHistoNo.size()); // resize to numer of groups + for (UInt_t i=0; iGetDataBin(forwardHistoNo[i])->size()); + backward[i].resize(runData->GetDataBin(backwardHistoNo[i])->size()); + forward[i] = *runData->GetDataBin(forwardHistoNo[i]); + backward[i] = *runData->GetDataBin(backwardHistoNo[i]); + } + + // check if addrun's are present, and if yes add data + // 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 << ">> PRunAsymmetryBNMR::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(forwardHistoNo[k])->size(); + for (UInt_t j=0; jGetDataBin(forwardHistoNo[k])->size(); j++) { // loop over the bin indices + // make sure that the index stays in the proper range + if (((Int_t)j+(Int_t)fAddT0s[i-1][2*k]-(Int_t)fT0s[2*k] >= 0) && (j+(Int_t)fAddT0s[i-1][2*k]-(Int_t)fT0s[2*k] < addRunSize)) { + forward[k][j] += addRunData->GetDataBin(forwardHistoNo[k])->at(j+(Int_t)fAddT0s[i-1][2*k]-(Int_t)fT0s[2*k]); + } + } + } + + // add backward run + for (UInt_t k=0; kGetDataBin(backwardHistoNo[k])->size(); + for (UInt_t j=0; jGetDataBin(backwardHistoNo[k])->size(); j++) { // loop over the bin indices + // make sure that the index stays in the proper range + if (((Int_t)j+(Int_t)fAddT0s[i-1][2*k+1]-(Int_t)fT0s[2*k+1] >= 0) && (j+(Int_t)fAddT0s[i-1][2*k+1]-(Int_t)fT0s[2*k+1] < addRunSize)) { + backward[k][j] += addRunData->GetDataBin(backwardHistoNo[k])->at(j+(Int_t)fAddT0s[i-1][2*k+1]-(Int_t)fT0s[2*k+1]); + } + } + } + } + } + + // set forward/backward histo data of the first group + fForward.resize(forward[0].size()); + fBackward.resize(backward[0].size()); + for (UInt_t i=0; iGetDataBin(forwardHistoNo[i])->size(); j++) { // loop over the bin indices + // make sure that the index stays within proper range + if ((j+fT0s[2*i]-fT0s[0] >= 0) && (j+fT0s[2*i]-fT0s[0] < runData->GetDataBin(forwardHistoNo[i])->size())) { + fForward[j] += forward[i][j+(Int_t)fT0s[2*i]-(Int_t)fT0s[0]]; + } + } + } + + // group histograms, add all the remaining backward histograms of the group + for (UInt_t i=1; iGetDataBin(backwardHistoNo[i])->size(); j++) { // loop over the bin indices + // make sure that the index stays within proper range + if ((j+fT0s[2*i+1]-fT0s[1] >= 0) && (j+fT0s[2*i+1]-fT0s[1] < runData->GetDataBin(backwardHistoNo[i])->size())) { + fBackward[j] += backward[i][j+(Int_t)fT0s[2*i+1]-(Int_t)fT0s[1]]; + } + } + } + + // subtract background from histogramms ------------------------------------------ + if (fRunInfo->GetBkgFix(0) == PMUSR_UNDEFINED) { // no fixed background given + if (fRunInfo->GetBkgRange(0) >= 0) { // background range given + if (!SubtractEstimatedBkg()) + return false; + } else { // no background given to do the job, try to estimate it + fRunInfo->SetBkgRange(static_cast(fT0s[0]*0.1), 0); + fRunInfo->SetBkgRange(static_cast(fT0s[0]*0.6), 1); + fRunInfo->SetBkgRange(static_cast(fT0s[1]*0.1), 2); + fRunInfo->SetBkgRange(static_cast(fT0s[1]*0.6), 3); + cerr << endl << ">> PRunAsymmetryBNMR::PrepareData(): **WARNING** Neither fix background nor background bins are given!"; + cerr << endl << ">> Will try the following:"; + cerr << endl << ">> forward: bkg start = " << fRunInfo->GetBkgRange(0) << ", bkg end = " << fRunInfo->GetBkgRange(1); + cerr << endl << ">> backward: bkg start = " << fRunInfo->GetBkgRange(2) << ", bkg end = " << fRunInfo->GetBkgRange(3); + cerr << endl << ">> NO GUARANTEE THAT THIS MAKES ANY SENSE! Better check ..."; + cerr << endl; + if (!SubtractEstimatedBkg()) + return false; + } + } else { // fixed background given + if (!SubtractFixBkg()) + return false; + } + + UInt_t histoNo[2] = {forwardHistoNo[0], backwardHistoNo[0]}; + + // get the data range (fgb/lgb) for the current RUN block + if (!GetProperDataRange(runData, histoNo)) { + return false; + } + + // get the fit range for the current RUN block + GetProperFitRange(globalBlock); + + // everything looks fine, hence fill data set + Bool_t status; + switch(fHandleTag) { + case kFit: + status = PrepareFitData(); + break; + case kView: + if (fMsrInfo->GetMsrPlotList()->at(0).fRRFPacking == 0) + status = PrepareViewData(runData, histoNo); + else + status = PrepareRRFViewData(runData, histoNo); + break; + default: + status = false; + break; + } + + // clean up + forwardHistoNo.clear(); + backwardHistoNo.clear(); + + return status; +} + +//-------------------------------------------------------------------------- +// SubtractFixBkg (private) +//-------------------------------------------------------------------------- +/** + *

Subtracts a fixed background from the raw data. The background is given + * in units of (1/bin); for the Asymmetry representation (1/ns) doesn't make too much sense. + * The error propagation is done the following way: it is assumed that the error of the background + * is Poisson like, i.e. \f$\Delta\mathrm{bkg} = \sqrt{\mathrm{bkg}}\f$. + * + * Error propagation: + * \f[ \Delta f_i^{\rm c} = \pm\left[ (\Delta f_i)^2 + (\Delta \mathrm{bkg})^2 \right]^{1/2} = + * \pm\left[ f_i + \mathrm{bkg} \right]^{1/2}, \f] + * where \f$ f_i^{\rm c} \f$ is the background corrected histogram, \f$ f_i \f$ the raw histogram + * and \f$ \mathrm{bkg} \f$ the fix given background. + * + * return: + * - true + * + */ +Bool_t PRunAsymmetryBNMR::SubtractFixBkg() +{ + Double_t dval; + for (UInt_t i=0; iGetBkgFix(0); + + // keep the error, and make sure that the bin is NOT empty + if (fBackward[i] != 0.0) + dval = TMath::Sqrt(fBackward[i]); + else + dval = 1.0; + fBackwardErr.push_back(dval); + fBackward[i] -= fRunInfo->GetBkgFix(1); + } + + return true; +} + +//-------------------------------------------------------------------------- +// SubtractEstimatedBkg (private) +//-------------------------------------------------------------------------- +/** + *

Subtracts the background which is estimated from a given interval (typically before t0). + * + * The background corrected histogramms are: + * \f$ f_i^{\rm c} = f_i - \mathrm{bkg} \f$, where \f$ f_i \f$ is the raw data histogram, + * \f$ \mathrm{bkg} \f$ the background estimate, and \f$ f_i^{\rm c} \f$ background corrected + * histogram. The error on \f$ f_i^{\rm c} \f$ is + * \f[ \Delta f_i^{\rm c} = \pm \sqrt{ (\Delta f_i)^2 + (\Delta \mathrm{bkg})^2 } = + * \pm \sqrt{f_i + (\Delta \mathrm{bkg})^2} \f] + * The background error \f$ \Delta \mathrm{bkg} \f$ is + * \f[ \Delta \mathrm{bkg} = \pm\frac{1}{N}\left[\sum_{i=0}^N (\Delta f_i)^2\right]^{1/2} = + * \pm\frac{1}{N}\left[\sum_{i=0}^N f_i \right]^{1/2},\f] + * where \f$N\f$ is the number of bins over which the background is formed. + * + * return: + * - true + */ +Bool_t PRunAsymmetryBNMR::SubtractEstimatedBkg() +{ + 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 + Int_t start[2] = {fRunInfo->GetBkgRange(0), fRunInfo->GetBkgRange(2)}; + Int_t end[2] = {fRunInfo->GetBkgRange(1), fRunInfo->GetBkgRange(3)}; + for (UInt_t i=0; i<2; i++) { + if (end[i] < start[i]) { + cout << endl << "PRunAsymmetryBNMR::SubtractEstimatedBkg(): end = " << end[i] << " > start = " << start[i] << "! Will swap them!"; + UInt_t keep = end[i]; + end[i] = start[i]; + start[i] = keep; + } + } + + // calculate proper background range + for (UInt_t i=0; i<2; i++) { + if (beamPeriod != 0.0) { + Double_t timeBkg = (Double_t)(end[i]-start[i])*(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[i] = start[i] + (UInt_t) ((fullCycles*beamPeriod)/(fTimeResolution*fPacking)); + cout << "PRunAsymmetryBNMR::SubtractEstimatedBkg(): Background " << start[i] << ", " << end[i] << endl; + if (end[i] == start[i]) + end[i] = fRunInfo->GetBkgRange(2*i+1); + } + } + + // check if start is within histogram bounds + if ((start[0] < 0) || (start[0] >= fForward.size()) || + (start[1] < 0) || (start[1] >= fBackward.size())) { + cerr << endl << ">> PRunAsymmetryBNMR::SubtractEstimatedBkg(): **ERROR** background bin values out of bound!"; + cerr << endl << ">> histo lengths (f/b) = (" << fForward.size() << "/" << fBackward.size() << ")."; + cerr << endl << ">> background start (f/b) = (" << start[0] << "/" << start[1] << ")."; + return false; + } + + // check if end is within histogram bounds + if ((end[0] < 0) || (end[0] >= fForward.size()) || + (end[1] < 0) || (end[1] >= fBackward.size())) { + cerr << endl << ">> PRunAsymmetryBNMR::SubtractEstimatedBkg(): **ERROR** background bin values out of bound!"; + cerr << endl << ">> histo lengths (f/b) = (" << fForward.size() << "/" << fBackward.size() << ")."; + cerr << endl << ">> background end (f/b) = (" << end[0] << "/" << end[1] << ")."; + return false; + } + + // calculate background + Double_t bkg[2] = {0.0, 0.0}; + Double_t errBkg[2] = {0.0, 0.0}; + + // forward + for (UInt_t i=start[0]; i<=end[0]; i++) + bkg[0] += fForward[i]; + errBkg[0] = TMath::Sqrt(bkg[0])/(end[0] - start[0] + 1); + bkg[0] /= static_cast(end[0] - start[0] + 1); + cout << endl << ">> estimated forward histo background: " << bkg[0]; + + // backward + for (UInt_t i=start[1]; i<=end[1]; i++) + bkg[1] += fBackward[i]; + errBkg[1] = TMath::Sqrt(bkg[1])/(end[1] - start[1] + 1); + bkg[1] /= static_cast(end[1] - start[1] + 1); + cout << endl << ">> estimated backward histo background: " << bkg[1] << endl; + + // correct error for forward, backward + Double_t errVal = 0.0; + for (UInt_t i=0; i 0.0) + errVal = TMath::Sqrt(fForward[i]+errBkg[0]*errBkg[0]); + else + errVal = 1.0; + fForwardErr.push_back(errVal); + if (fBackward[i] > 0.0) + errVal = TMath::Sqrt(fBackward[i]+errBkg[1]*errBkg[1]); + else + errVal = 1.0; + fBackwardErr.push_back(errVal); + } + + // subtract background from data + for (UInt_t i=0; iSetBkgEstimated(bkg[0], 0); + fRunInfo->SetBkgEstimated(bkg[1], 1); + + return true; +} + +//-------------------------------------------------------------------------- +// PrepareFitData (protected) +//-------------------------------------------------------------------------- +/** + *

Take the pre-processed data (i.e. grouping and addrun are preformed) and form the asymmetry for fitting. + * Before forming the asymmetry, the following checks will be performed: + * -# check if data range is given, if not try to estimate one. + * -# check that if a data range is present, that it makes any sense. + * -# check that 'first good bin'-'t0' is the same for forward and backward histogram. If not adjust it. + * -# pack data (rebin). + * -# if packed forward size != backward size, truncate the longer one such that an asymmetry can be formed. + * -# calculate the asymmetry: \f$ A_i = (f_i^c-b_i^c)/(f_i^c+b_i^c) \f$ + * -# calculate the asymmetry errors: \f$ \delta A_i = 2 \sqrt{(b_i^c)^2 (\delta f_i^c)^2 + (\delta b_i^c)^2 (f_i^c)^2}/(f_i^c+b_i^c)^2\f$ + */ +Bool_t PRunAsymmetryBNMR::PrepareFitData() +{ + // transform raw histo data. This is done the following way (for details see the manual): + // first rebin the data, than calculate the asymmetry + + // everything looks fine, hence fill packed forward and backward histo + PRunData forwardPacked; + PRunData backwardPacked; + Double_t value = 0.0; + Double_t error = 0.0; + // forward + for (Int_t i=fGoodBins[0]; i 1 + if (((i-fGoodBins[0]) % fPacking == 0) && (i != fGoodBins[0])) { // fill data + // in order that after rebinning the fit does not need to be redone (important for plots) + // the value is normalize to per bin + value /= fPacking; + forwardPacked.AppendValue(value); + if (value == 0.0) + forwardPacked.AppendErrorValue(1.0); + else + forwardPacked.AppendErrorValue(TMath::Sqrt(error)/fPacking); + value = 0.0; + error = 0.0; + } + value += fForward[i]; + error += fForwardErr[i]*fForwardErr[i]; + } + } + // backward + for (Int_t i=fGoodBins[2]; i 1 + if (((i-fGoodBins[2]) % fPacking == 0) && (i != fGoodBins[2])) { // fill data + // in order that after rebinning the fit does not need to be redone (important for plots) + // the value is normalize to per bin + value /= fPacking; + backwardPacked.AppendValue(value); + if (value == 0.0) + backwardPacked.AppendErrorValue(1.0); + else + backwardPacked.AppendErrorValue(TMath::Sqrt(error)/fPacking); + value = 0.0; + error = 0.0; + } + value += fBackward[i]; + error += fBackwardErr[i]*fBackwardErr[i]; + } + } + + // check if packed forward and backward hist have the same size, otherwise take the minimum size + UInt_t noOfBins = forwardPacked.GetValue()->size(); + if (forwardPacked.GetValue()->size() != backwardPacked.GetValue()->size()) { + if (forwardPacked.GetValue()->size() > backwardPacked.GetValue()->size()) + noOfBins = backwardPacked.GetValue()->size(); + } + + // form asymmetry including error propagation + Double_t asym; + Double_t f, b, ef, eb; + // fill data time start, and step + // data start at data_start-t0 shifted by (pack-1)/2 + fData.SetDataTimeStart(fTimeResolution*((Double_t)fGoodBins[0]-fT0s[0]+(Double_t)(fPacking-1)/2.0)); + fData.SetDataTimeStep(fTimeResolution*(Double_t)fPacking); + for (UInt_t i=0; iat(i); + b = backwardPacked.GetValue()->at(i); + ef = forwardPacked.GetError()->at(i); + eb = backwardPacked.GetError()->at(i); + // check that there are indeed bins + if (f+b != 0.0) + asym = (f-b) / (f+b); + else + asym = 0.0; + fData.AppendValue(asym); + // calculate the error + if (f+b != 0.0) + error = 2.0/((f+b)*(f+b))*TMath::Sqrt(b*b*ef*ef+eb*eb*f*f); + else + error = 1.0; + fData.AppendErrorValue(error); + } + + CalcNoOfFitBins(); + + // clean up + fForward.clear(); + fForwardErr.clear(); + fBackward.clear(); + fBackwardErr.clear(); + + return true; +} + +//-------------------------------------------------------------------------- +// PrepareViewData (protected) +//-------------------------------------------------------------------------- +/** + *

Take the pre-processed data (i.e. grouping and addrun are preformed) and form the asymmetry for view representation. + * Before forming the asymmetry, the following checks will be performed: + * -# check if view packing is whished. + * -# check if data range is given, if not try to estimate one. + * -# check that data range is present, that it makes any sense. + * -# check that 'first good bin'-'t0' is the same for forward and backward histogram. If not adjust it. + * -# pack data (rebin). + * -# if packed forward size != backward size, truncate the longer one such that an asymmetry can be formed. + * -# calculate the asymmetry: \f$ A_i = (\alpha f_i^c-b_i^c)/(\alpha \beta f_i^c+b_i^c) \f$ + * -# calculate the asymmetry errors: \f$ \delta A_i = 2 \sqrt{(b_i^c)^2 (\delta f_i^c)^2 + (\delta b_i^c)^2 (f_i^c)^2}/(f_i^c+b_i^c)^2\f$ + * -# calculate the theory vector. + * + * \param runData raw run data needed to perform some crosschecks + * \param histoNo histogram number (within a run). histoNo[0]: forward histogram number, histNo[1]: backward histogram number + */ +Bool_t PRunAsymmetryBNMR::PrepareViewData(PRawRunData* runData, UInt_t histoNo[2]) +{ + // check if view_packing is wished + Int_t packing = fPacking; + if (fMsrInfo->GetMsrPlotList()->at(0).fViewPacking > 0) { + packing = fMsrInfo->GetMsrPlotList()->at(0).fViewPacking; + } + + // feed the parameter vector + std::vector par; + PMsrParamList *paramList = fMsrInfo->GetMsrParamList(); + for (UInt_t i=0; isize(); i++) + par.push_back((*paramList)[i].fValue); + + // transform raw histo data. This is done the following way (for details see the manual): + // first rebin the data, than calculate the asymmetry + + // first get start data, end data, and t0 + Int_t start[2] = {fGoodBins[0], fGoodBins[2]}; + Int_t end[2] = {fGoodBins[1], fGoodBins[3]}; + Int_t t0[2] = {(Int_t)fT0s[0], (Int_t)fT0s[1]}; + + // check if the data ranges and t0's between forward/backward are compatible + Int_t fgb[2]; + if (start[0]-t0[0] != start[1]-t0[1]) { // wrong fgb aligning + if (abs(start[0]-t0[0]) > abs(start[1]-t0[1])) { + fgb[0] = start[0]; + fgb[1] = t0[1] + start[0]-t0[0]; + cerr << endl << ">> PRunAsymmetryBNMR::PrepareViewData(): **WARNING** needed to shift backward fgb from "; + cerr << start[1] << " to " << fgb[1] << endl; + } else { + fgb[0] = t0[0] + start[1]-t0[1]; + fgb[1] = start[1]; + cerr << endl << ">> PRunAsymmetryBNMR::PrepareViewData(): **WARNING** needed to shift forward fgb from "; + cerr << start[0] << " to " << fgb[0] << endl; + } + } else { // fgb aligning is correct + fgb[0] = start[0]; + fgb[1] = start[1]; + } + + Int_t val = fgb[0]-packing*(fgb[0]/packing); + do { + if (fgb[1] - fgb[0] < 0) + val += packing; + } while (val + fgb[1] - fgb[0] < 0); + + start[0] = val; + start[1] = val + fgb[1] - fgb[0]; + + // make sure that there are equal number of rebinned bins in forward and backward + UInt_t noOfBins0 = (runData->GetDataBin(histoNo[0])->size()-start[0])/packing; + UInt_t noOfBins1 = (runData->GetDataBin(histoNo[1])->size()-start[1])/packing; + if (noOfBins0 > noOfBins1) + noOfBins0 = noOfBins1; + end[0] = start[0] + noOfBins0 * packing; + end[1] = start[1] + noOfBins0 * packing; + + // check if start, end, and t0 make any sense + // 1st check if start and end are in proper order + for (UInt_t i=0; i<2; i++) { + if (end[i] < start[i]) { // need to swap them + Int_t keep = end[i]; + end[i] = start[i]; + start[i] = keep; + } + // 2nd check if start is within proper bounds + if ((start[i] < 0) || (start[i] > (Int_t)runData->GetDataBin(histoNo[i])->size())) { + cerr << endl << ">> PRunAsymmetryBNMR::PrepareViewData(): **ERROR** start data bin doesn't make any sense!"; + cerr << endl; + return false; + } + // 3rd check if end is within proper bounds + if ((end[i] < 0) || (end[i] > (Int_t)runData->GetDataBin(histoNo[i])->size())) { + cerr << endl << ">> PRunAsymmetryBNMR::PrepareViewData(): **ERROR** end data bin doesn't make any sense!"; + cerr << endl; + return false; + } + // 4th check if t0 is within proper bounds + if ((t0[i] < 0) || (t0[i] > (Int_t)runData->GetDataBin(histoNo[i])->size())) { + cerr << endl << ">> PRunAsymmetryBNMR::PrepareViewData(): **ERROR** t0 data bin doesn't make any sense!"; + cerr << endl; + return false; + } + } + + // everything looks fine, hence fill packed forward and backward histo + PRunData forwardPacked; + PRunData backwardPacked; + Double_t value = 0.0; + Double_t error = 0.0; + + // forward + for (Int_t i=start[0]; i 1 + if (((i-start[0]) % packing == 0) && (i != start[0])) { // fill data + // in order that after rebinning the fit does not need to be redone (important for plots) + // the value is normalize to per bin + value /= packing; + forwardPacked.AppendValue(value); + if (value == 0.0) + forwardPacked.AppendErrorValue(1.0); + else + forwardPacked.AppendErrorValue(TMath::Sqrt(error)/packing); + value = 0.0; + error = 0.0; + } + value += fForward[i]; + error += fForwardErr[i]*fForwardErr[i]; + } + } + + // backward + for (Int_t i=start[1]; i 1 + if (((i-start[1]) % packing == 0) && (i != start[1])) { // fill data + // in order that after rebinning the fit does not need to be redone (important for plots) + // the value is normalize to per bin + value /= packing; + backwardPacked.AppendValue(value); + if (value == 0.0) + backwardPacked.AppendErrorValue(1.0); + else + backwardPacked.AppendErrorValue(TMath::Sqrt(error)/packing); + value = 0.0; + error = 0.0; + } + value += fBackward[i]; + error += fBackwardErr[i]*fBackwardErr[i]; + } + } + + // check if packed forward and backward hist have the same size, otherwise take the minimum size + UInt_t noOfBins = forwardPacked.GetValue()->size(); + if (forwardPacked.GetValue()->size() != backwardPacked.GetValue()->size()) { + if (forwardPacked.GetValue()->size() > backwardPacked.GetValue()->size()) + noOfBins = backwardPacked.GetValue()->size(); + } + + // form asymmetry including error propagation + Double_t asym; + Double_t f, b, ef, eb, alpha = 1.0, beta = 1.0; + // set data time start, and step + // data start at data_start-t0 + fData.SetDataTimeStart(fTimeResolution*((Double_t)start[0]-t0[0]+(Double_t)(packing-1)/2.0)); + fData.SetDataTimeStep(fTimeResolution*(Double_t)packing); + + // get the proper alpha and beta + switch (fAlphaBetaTag) { + case 1: // alpha == 1, beta == 1 + alpha = 1.0; + beta = 1.0; + break; + case 2: // alpha != 1, beta == 1 + alpha = par[fRunInfo->GetAlphaParamNo()-1]; + beta = 1.0; + break; + case 3: // alpha == 1, beta != 1 + alpha = 1.0; + beta = par[fRunInfo->GetBetaParamNo()-1]; + break; + case 4: // alpha != 1, beta != 1 + alpha = par[fRunInfo->GetAlphaParamNo()-1]; + beta = par[fRunInfo->GetBetaParamNo()-1]; + break; + default: + break; + } + + for (UInt_t i=0; isize(); i++) { + // to make the formulae more readable + f = forwardPacked.GetValue()->at(i); + b = backwardPacked.GetValue()->at(i); + ef = forwardPacked.GetError()->at(i); + eb = backwardPacked.GetError()->at(i); + // check that there are indeed bins + if (f+b != 0.0) + asym = (alpha*f-b) / (alpha*beta*f+b); + else + asym = 0.0; + fData.AppendValue(asym); + // calculate the error + if (f+b != 0.0) + error = 2.0/((f+b)*(f+b))*TMath::Sqrt(b*b*ef*ef+eb*eb*f*f); + else + error = 1.0; + fData.AppendErrorValue(error); + } + + CalcNoOfFitBins(); + + // clean up + fForward.clear(); + fForwardErr.clear(); + fBackward.clear(); + fBackwardErr.clear(); + + // fill theory vector for kView + // calculate functions + for (Int_t i=0; iGetNoOfFuncs(); i++) { + fFuncValues[i] = fMsrInfo->EvalFunc(fMsrInfo->GetFuncNo(i), *fRunInfo->GetMap(), par); + } + + // calculate theory + Double_t time; + UInt_t size = runData->GetDataBin(histoNo[0])->size(); + Double_t factor = 1.0; + if (fData.GetValue()->size() * 10 > runData->GetDataBin(histoNo[0])->size()) { + size = fData.GetValue()->size() * 10; + factor = (Double_t)runData->GetDataBin(histoNo[0])->size() / (Double_t)size; + } + fData.SetTheoryTimeStart(fData.GetDataTimeStart()); + fData.SetTheoryTimeStep(fTimeResolution*factor); + for (UInt_t i=0; iFunc(time, par, fFuncValues); + if (fabs(value) > 10.0) { // dirty hack needs to be fixed!! + value = 0.0; + } + fData.AppendTheoryValue(value); + } + + // clean up + par.clear(); + + return true; +} + +//-------------------------------------------------------------------------- +// PrepareRRFViewData (protected) +//-------------------------------------------------------------------------- +/** + *

Prepares the RRF data set for visual representation. This is done the following way: + * -# make all necessary checks + * -# build the asymmetry, \f$ A(t) \f$, WITHOUT packing. + * -# \f$ A_R(t) = A(t) \cdot 2 \cos(\omega_R t + \phi_R) \f$ + * -# do the packing of \f$ A_R(t) \f$ + * -# calculate theory, \f$ T(t) \f$, as close as possible to the time resolution [compatible with the RRF frequency] + * -# \f$ T_R(t) = T(t) \cdot 2 \cos(\omega_R t + \phi_R) \f$ + * -# do the packing of \f$ T_R(t) \f$ + * -# calculate the Kaiser FIR filter coefficients + * -# filter \f$ T_R(t) \f$. + * + * \param runData raw run data needed to perform some crosschecks + * \param histoNo array of the histo numbers form which to build the asymmetry + */ +Bool_t PRunAsymmetryBNMR::PrepareRRFViewData(PRawRunData* runData, UInt_t histoNo[2]) +{ + // feed the parameter vector + std::vector par; + PMsrParamList *paramList = fMsrInfo->GetMsrParamList(); + for (UInt_t i=0; isize(); i++) + par.push_back((*paramList)[i].fValue); + + // ------------------------------------------------------------ + // 1. make all necessary checks + // ------------------------------------------------------------ + + // first get start data, end data, and t0 + Int_t start[2] = {fGoodBins[0], fGoodBins[2]}; + Int_t end[2] = {fGoodBins[1], fGoodBins[3]}; + Int_t t0[2] = {(Int_t)fT0s[0], (Int_t)fT0s[1]}; + UInt_t packing = fMsrInfo->GetMsrPlotList()->at(0).fRRFPacking; + + // check if the data ranges and t0's between forward/backward are compatible + Int_t fgb[2]; + if (start[0]-t0[0] != start[1]-t0[1]) { // wrong fgb aligning + if (abs(start[0]-t0[0]) > abs(start[1]-t0[1])) { + fgb[0] = start[0]; + fgb[1] = t0[1] + start[0]-t0[0]; + cerr << endl << ">> PRunAsymmetryBNMR::PrepareRRFViewData(): **WARNING** needed to shift backward fgb from "; + cerr << start[1] << " to " << fgb[1] << endl; + } else { + fgb[0] = t0[0] + start[1]-t0[1]; + fgb[1] = start[1]; + cerr << endl << ">> PRunAsymmetryBNMR::PrepareRRFViewData(): **WARNING** needed to shift forward fgb from "; + cerr << start[1] << " to " << fgb[0] << endl; + } + } else { // fgb aligning is correct + fgb[0] = start[0]; + fgb[1] = start[1]; + } + + Int_t val = fgb[0]-packing*(fgb[0]/packing); + do { + if (fgb[1] - fgb[0] < 0) + val += packing; + } while (val + fgb[1] - fgb[0] < 0); + + start[0] = val; + start[1] = val + fgb[1] - fgb[0]; + + // make sure that there are equal number of rebinned bins in forward and backward + UInt_t noOfBins0 = runData->GetDataBin(histoNo[0])->size()-start[0]; + UInt_t noOfBins1 = runData->GetDataBin(histoNo[1])->size()-start[1]; + if (noOfBins0 > noOfBins1) + noOfBins0 = noOfBins1; + end[0] = start[0] + noOfBins0; + end[1] = start[1] + noOfBins0; + + // check if start, end, and t0 make any sense + // 1st check if start and end are in proper order + for (UInt_t i=0; i<2; i++) { + if (end[i] < start[i]) { // need to swap them + Int_t keep = end[i]; + end[i] = start[i]; + start[i] = keep; + } + // 2nd check if start is within proper bounds + if ((start[i] < 0) || (start[i] > (Int_t)runData->GetDataBin(histoNo[i])->size())) { + cerr << endl << ">> PRunAsymmetryBNMR::PrepareRRFViewData(): **ERROR** start data bin doesn't make any sense!"; + cerr << endl; + return false; + } + // 3rd check if end is within proper bounds + if ((end[i] < 0) || (end[i] > (Int_t)runData->GetDataBin(histoNo[i])->size())) { + cerr << endl << ">> PRunAsymmetryBNMR::PrepareRRFViewData(): **ERROR** end data bin doesn't make any sense!"; + cerr << endl; + return false; + } + // 4th check if t0 is within proper bounds + if ((t0[i] < 0) || (t0[i] > (Int_t)runData->GetDataBin(histoNo[i])->size())) { + cerr << endl << ">> PRunAsymmetryBNMR::PrepareRRFViewData(): **ERROR** t0 data bin doesn't make any sense!"; + cerr << endl; + return false; + } + } + + // ------------------------------------------------------------ + // 2. build the asymmetry [A(t)] WITHOUT packing. + // ------------------------------------------------------------ + + PDoubleVector forward, forwardErr; + PDoubleVector backward, backwardErr; + Double_t error = 0.0; + // forward + for (Int_t i=start[0]; i backward.size()) + noOfBins = backward.size(); + } + + // form asymmetry including error propagation + PDoubleVector asymmetry, asymmetryErr; + Double_t asym; + Double_t f, b, ef, eb, alpha = 1.0, beta = 1.0; + + // get the proper alpha and beta + switch (fAlphaBetaTag) { + case 1: // alpha == 1, beta == 1 + alpha = 1.0; + beta = 1.0; + break; + case 2: // alpha != 1, beta == 1 + alpha = par[fRunInfo->GetAlphaParamNo()-1]; + beta = 1.0; + break; + case 3: // alpha == 1, beta != 1 + alpha = 1.0; + beta = par[fRunInfo->GetBetaParamNo()-1]; + break; + case 4: // alpha != 1, beta != 1 + alpha = par[fRunInfo->GetAlphaParamNo()-1]; + beta = par[fRunInfo->GetBetaParamNo()-1]; + break; + default: + break; + } + + for (UInt_t i=0; 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(); + + for (UInt_t i=0; i(start[0])-t0[0]+static_cast(i)); + asymmetry[i] *= 2.0*TMath::Cos(wRRF*time+phaseRRF); + } + + // ------------------------------------------------------------ + // 4. do the packing of A_R(t) + // ------------------------------------------------------------ + Double_t value = 0.0; + error = 0.0; + for (UInt_t i=0; i(packing-1)/2.0)); + fData.SetDataTimeStep(fTimeResolution*static_cast(packing)); + + // ------------------------------------------------------------ + // 5. calculate theory [T(t)] as close as possible to the time resolution [compatible with the RRF frequency] + // 6. T_R(t) = T(t) * 2 cos(w_R t + phi_R) + // ------------------------------------------------------------ + UInt_t 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 + + // calculate functions + for (Int_t i=0; iGetNoOfFuncs(); i++) { + fFuncValues[i] = fMsrInfo->EvalFunc(fMsrInfo->GetFuncNo(i), *fRunInfo->GetMap(), par); + } + + Double_t theoryValue; + for (UInt_t i=0; i(i)*fData.GetTheoryTimeStep(); + theoryValue = fTheory->Func(time, par, fFuncValues); + 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); + } + + // ------------------------------------------------------------ + // 7. do the packing of T_R(t) + // ------------------------------------------------------------ + + PDoubleVector theo; + Double_t dval = 0.0; + 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.ReplaceTheory(theo); + theo.clear(); + + // set the theory time start and step size + fData.SetTheoryTimeStart(fData.GetTheoryTimeStart()+static_cast(rebinRRF-1)*fData.GetTheoryTimeStep()/2.0); + fData.SetTheoryTimeStep(rebinRRF*fData.GetTheoryTimeStep()); + + // ------------------------------------------------------------ + // 8. calculate the Kaiser FIR filter coefficients + // ------------------------------------------------------------ + 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) + + // ------------------------------------------------------------ + // 9. filter T_R(t) + // ------------------------------------------------------------ + FilterTheo(); + + // clean up + par.clear(); + forward.clear(); + forwardErr.clear(); + backward.clear(); + backwardErr.clear(); + asymmetry.clear(); + asymmetryErr.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 forwardHistoNo histogram number vector of forward; forwardHistoNo = msr-file forward + redGreen_offset - 1 + * \param backwardHistoNo histogram number vector of backwardward; backwardHistoNo = msr-file backward + redGreen_offset - 1 + * + * return: + * - true if everthing went smooth + * - false, otherwise. + */ +Bool_t PRunAsymmetryBNMR::GetProperT0(PRawRunData* runData, PMsrGlobalBlock *globalBlock, PUIntVector &forwardHistoNo, PUIntVector &backwardHistoNo) +{ + // feed all T0's + // first init T0's, T0's are stored as (forward T0, backward T0, etc.) + fT0s.clear(); + fT0s.resize(2*forwardHistoNo.size()); + for (UInt_t i=0; iGetT0BinSize(); i++) { + fT0s[i] = fRunInfo->GetT0Bin(i); + } + + // fill in the missing 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 missing T0's from the data file, if not already present in the msr-file + for (UInt_t i=0; iGetT0Bin(forwardHistoNo[i]) > 0.0) { + fT0s[2*i] = runData->GetT0Bin(forwardHistoNo[i]); + fRunInfo->SetT0Bin(fT0s[2*i], 2*i); + } + } + for (UInt_t i=0; iGetT0Bin(backwardHistoNo[i]) > 0.0) { + fT0s[2*i+1] = runData->GetT0Bin(backwardHistoNo[i]); + fRunInfo->SetT0Bin(fT0s[2*i+1], 2*i+1); + } + } + + // fill in the T0's gaps, i.e. in case the T0's are NEITHER in the msr-file and NOR in the data file + for (UInt_t i=0; iGetT0BinEstimated(forwardHistoNo[i]); + fRunInfo->SetT0Bin(fT0s[2*i], 2*i); + + cerr << endl << ">> PRunAsymmetryBNMR::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(forwardHistoNo[i]); + cerr << endl << ">> NO GUARANTEE THAT THIS OK!! For instance for LEM this is almost for sure rubbish!"; + cerr << endl; + } + } + for (UInt_t i=0; iGetT0BinEstimated(backwardHistoNo[i]); + fRunInfo->SetT0Bin(fT0s[2*i+1], 2*i+1); + + cerr << endl << ">> PRunAsymmetryBNMR::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: backward t0 = " << runData->GetT0BinEstimated(backwardHistoNo[i]); + cerr << endl << ">> NO GUARANTEE 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; i (Int_t)runData->GetDataBin(forwardHistoNo[i])->size())) { + cerr << endl << ">> PRunAsymmetryBNMR::GetProperT0(): **ERROR** t0 data bin (" << fT0s[2*i] << ") doesn't make any sense!"; + cerr << endl << ">> forwardHistoNo " << forwardHistoNo[i]; + cerr << endl; + return false; + } + if ((fT0s[2*i+1] < 0) || (fT0s[2*i+1] > (Int_t)runData->GetDataBin(backwardHistoNo[i])->size())) { + cerr << endl << ">> PRunAsymmetryBNMR::PrepareData(): **ERROR** t0 data bin (" << fT0s[2*i+1] << ") doesn't make any sense!"; + cerr << endl << ">> backwardHistoNo " << backwardHistoNo[i]; + cerr << endl; + return false; + } + } + + // check if addrun's are present, and if yes add the necessary 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 << ">> PRunAsymmetryBNMR::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].clear(); + fAddT0s[i-1].resize(2*forwardHistoNo.size()); + for (UInt_t j=0; jGetAddT0BinSize(i); j++) { + fAddT0s[i-1][j] = fRunInfo->GetAddT0Bin(i, j); + } + + // fill in the T0's from the data file, if not already present in the msr-file + for (UInt_t j=0; jGetT0Bin(forwardHistoNo[j]) > 0.0) { + fAddT0s[i-1][2*j] = addRunData->GetT0Bin(forwardHistoNo[j]); + fRunInfo->SetAddT0Bin(fAddT0s[i-1][2*j], i-1, 2*j); + } + } + for (UInt_t j=0; jGetT0Bin(backwardHistoNo[j]) > 0.0) { + fAddT0s[i-1][2*j+1] = addRunData->GetT0Bin(backwardHistoNo[j]); + fRunInfo->SetAddT0Bin(fAddT0s[i-1][2*j+1], i-1, 2*j+1); + } + } + + // 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(forwardHistoNo[j]); + fRunInfo->SetAddT0Bin(fAddT0s[i-1][2*j], i-1, 2*j); + + cerr << endl << ">> PRunAsymmetryBNMR::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(forwardHistoNo[j]); + cerr << endl << ">> NO GUARANTEE THAT THIS OK!! For instance for LEM this is almost for sure rubbish!"; + cerr << endl; + } + } + for (UInt_t j=0; jGetT0BinEstimated(backwardHistoNo[j]); + fRunInfo->SetAddT0Bin(fAddT0s[i-1][2*j+1], i-1, 2*j+1); + + cerr << endl << ">> PRunAsymmetryBNMR::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: backward t0 = " << runData->GetT0BinEstimated(backwardHistoNo[j]); + cerr << endl << ">> NO GUARANTEE THAT THIS OK!! For instance for LEM this is almost for sure rubbish!"; + cerr << endl; + } + } + } + } + + 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. + * + * \param runData raw run data needed to perform some crosschecks + * \param histoNo histogram number (within a run). histoNo[0]: forward histogram number, histNo[1]: backward histogram number + * + * return: + * - true if everthing went smooth + * - false, otherwise. + */ +Bool_t PRunAsymmetryBNMR::GetProperDataRange(PRawRunData* runData, UInt_t histoNo[2]) +{ + // first get start/end data + Int_t start[2] = {fRunInfo->GetDataRange(0), fRunInfo->GetDataRange(2)}; + Int_t end[2] = {fRunInfo->GetDataRange(1), fRunInfo->GetDataRange(3)}; + // check if data range has been provided in the RUN block. If not, try the GLOBAL block + if (start[0] == -1) { + start[0] = fMsrInfo->GetMsrGlobal()->GetDataRange(0); + } + if (start[1] == -1) { + start[1] = fMsrInfo->GetMsrGlobal()->GetDataRange(2); + } + if (end[0] == -1) { + end[0] = fMsrInfo->GetMsrGlobal()->GetDataRange(1); + } + if (end[1] == -1) { + end[1] = fMsrInfo->GetMsrGlobal()->GetDataRange(3); + } + + Double_t t0[2] = {fT0s[0], fT0s[1]}; + Int_t offset = (Int_t)(10.0e-3/fTimeResolution); // needed in case first good bin is not given, default = 10ns + + // check if data range has been provided, and if not try to estimate them + if (start[0] < 0) { + start[0] = (Int_t)t0[0]+offset; + fRunInfo->SetDataRange(start[0], 0); + cerr << endl << ">> PRunAsymmetryBNMR::GetProperDataRange(): **WARNING** data range (forward) was not provided, will try data range start = t0+" << offset << "(=10ns) = " << start[0] << "."; + cerr << endl << ">> NO GUARANTEE THAT THIS DOES MAKE ANY SENSE."; + cerr << endl; + } + if (start[1] < 0) { + start[1] = (Int_t)t0[1]+offset; + fRunInfo->SetDataRange(start[1], 2); + cerr << endl << ">> PRunAsymmetryBNMR::GetProperDataRange(): **WARNING** data range (backward) was not provided, will try data range start = t0+" << offset << "(=10ns) = " << start[1] << "."; + cerr << endl << ">> NO GUARANTEE THAT THIS DOES MAKE ANY SENSE."; + cerr << endl; + } + if (end[0] < 0) { + end[0] = runData->GetDataBin(histoNo[0])->size(); + fRunInfo->SetDataRange(end[0], 1); + cerr << endl << ">> PRunAsymmetryBNMR::GetProperDataRange(): **WARNING** data range (forward) was not provided, will try data range end = " << end[0] << "."; + cerr << endl << ">> NO GUARANTEE THAT THIS DOES MAKE ANY SENSE."; + cerr << endl; + } + if (end[1] < 0) { + end[1] = runData->GetDataBin(histoNo[1])->size(); + fRunInfo->SetDataRange(end[1], 3); + cerr << endl << ">> PRunAsymmetryBNMR::GetProperDataRange(): **WARNING** data range (backward) was not provided, will try data range end = " << end[1] << "."; + cerr << endl << ">> NO GUARANTEE 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 + for (UInt_t i=0; i<2; i++) { + if (end[i] < start[i]) { // need to swap them + Int_t keep = end[i]; + end[i] = start[i]; + start[i] = keep; + } + // 2nd check if start is within proper bounds + if ((start[i] < 0) || (start[i] > (Int_t)runData->GetDataBin(histoNo[i])->size())) { + cerr << endl << ">> PRunAsymmetryBNMR::GetProperDataRange(): **ERROR** start data bin doesn't make any sense!"; + cerr << endl; + return false; + } + // 3rd check if end is within proper bounds + if (end[i] < 0) { + cerr << endl << ">> PRunAsymmetryBNMR::GetProperDataRange(): **ERROR** end data bin (" << end[i] << ") doesn't make any sense!"; + cerr << endl; + return false; + } + if (end[i] > (Int_t)runData->GetDataBin(histoNo[i])->size()) { + cerr << endl << ">> PRunAsymmetryBNMR::GetProperDataRange(): **WARNING** end data bin (" << end[i] << ") > histo length (" << (Int_t)runData->GetDataBin(histoNo[i])->size() << ")."; + cerr << endl << ">> Will set end = (histo length - 1). Consider to change it in the msr-file." << endl; + cerr << endl; + end[i] = (Int_t)runData->GetDataBin(histoNo[i])->size()-1; + } + // 4th check if t0 is within proper bounds + if ((t0[i] < 0) || (t0[i] > (Int_t)runData->GetDataBin(histoNo[i])->size())) { + cerr << endl << ">> PRunAsymmetryBNMR::GetProperDataRange(): **ERROR** t0 data bin doesn't make any sense!"; + cerr << endl; + return false; + } + } + + // check that start-t0 is the same for forward as for backward, otherwise take max(start[i]-t0[i]) + if (fabs(static_cast(start[0])-t0[0]) > fabs(static_cast(start[1])-t0[1])){ + start[1] = static_cast(t0[1] + static_cast(start[0]) - t0[0]); + end[1] = static_cast(t0[1] + static_cast(end[0]) - t0[0]); + cerr << endl << ">> PRunAsymmetryBNMR::GetProperDataRange **WARNING** needed to shift backward data range."; + cerr << endl << ">> given: " << fRunInfo->GetDataRange(2) << ", " << fRunInfo->GetDataRange(3); + cerr << endl << ">> used : " << start[1] << ", " << end[1]; + cerr << endl; + } + if (fabs(static_cast(start[0])-t0[0]) < fabs(static_cast(start[1])-t0[1])){ + start[0] = static_cast(t0[0] + static_cast(start[1]) - t0[1]); + end[0] = static_cast(t0[0] + static_cast(end[1]) - t0[1]); + cerr << endl << ">> PRunAsymmetryBNMR::GetProperDataRange **WARNING** needed to shift forward data range."; + cerr << endl << ">> given: " << fRunInfo->GetDataRange(0) << ", " << fRunInfo->GetDataRange(1); + cerr << endl << ">> used : " << start[0] << ", " << end[0]; + cerr << endl; + } + + // keep good bins for potential latter use + fGoodBins[0] = start[0]; + fGoodBins[1] = end[0]; + fGoodBins[2] = start[1]; + fGoodBins[3] = end[1]; + + 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 PRunAsymmetryBNMR::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 << ">> PRunSingleHisto::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; + } +} diff --git a/src/classes/PRunListCollection.cpp b/src/classes/PRunListCollection.cpp index 7013ee94..dd0ecca1 100644 --- a/src/classes/PRunListCollection.cpp +++ b/src/classes/PRunListCollection.cpp @@ -1054,49 +1054,6 @@ PRunData* PRunListCollection::GetAsymmetryRRF(UInt_t index, EDataSwitch tag) return data; } -//-------------------------------------------------------------------------- -// GetAsymmetryBNMR (public) -//-------------------------------------------------------------------------- -/** - *

Get a processed asymmetry BNMR 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::GetAsymmetryBNMR(UInt_t index, EDataSwitch tag) -{ - PRunData *data = 0; - - switch (tag) { - case kIndex: // called from musrfit when dumping the data - if (index > fRunAsymmetryBNMRList.size()) { - cerr << endl << ">> PRunListCollection::GetAsymmetryBNMR(): **ERROR** index = " << index << " out of bounds"; - cerr << endl; - return 0; - } - - fRunAsymmetryBNMRList[index]->CalcTheory(); - data = fRunAsymmetryBNMRList[index]->GetData(); - break; - case kRunNo: // called from PMusrCanvas - for (UInt_t i=0; iGetRunNo() == index) { - data = fRunAsymmetryBNMRList[i]->GetData(); - break; - } - } - break; - default: // error - break; - } - - return data; -} - //-------------------------------------------------------------------------- // GetMuMinus (public) //-------------------------------------------------------------------------- diff --git a/src/include/PRunAsymmetryBNMR.h b/src/include/PRunAsymmetryBNMR.h new file mode 100644 index 00000000..ec2f80a7 --- /dev/null +++ b/src/include/PRunAsymmetryBNMR.h @@ -0,0 +1,89 @@ +/*************************************************************************** + + PRunAsymmetry.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 _PRUNASYMMETRYBNMR_H_ +#define _PRUNASYMMETRYBNMR_H_ + +#include "PRunBase.h" + +//--------------------------------------------------------------------------- +/** + *

Class handling the asymmetry fit. + */ +class PRunAsymmetryBNMR : public PRunBase +{ + public: + PRunAsymmetryBNMR(); + PRunAsymmetryBNMR(PMsrHandler *msrInfo, PRunDataHandler *rawData, UInt_t runNo, EPMusrHandleTag tag); + virtual ~PRunAsymmetryBNMR(); + + 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); + + virtual Int_t GetStartTimeBin() { return fStartTimeBin; } + virtual Int_t GetEndTimeBin() { return fEndTimeBin; } + virtual Int_t GetPacking() { return fPacking; } + + protected: + virtual void CalcNoOfFitBins(); + virtual Bool_t PrepareData(); + virtual Bool_t PrepareFitData(); + virtual Bool_t PrepareViewData(PRawRunData* runData, UInt_t histoNo[2]); + virtual Bool_t PrepareRRFViewData(PRawRunData* runData, UInt_t histoNo[2]); + + private: + UInt_t fAlphaBetaTag; ///< \f$ 1 \to \alpha = \beta = 1\f$; \f$ 2 \to \alpha \neq 1, \beta = 1\f$; \f$ 3 \to \alpha = 1, \beta \neq 1\f$; \f$ 4 \to \alpha \neq 1, \beta \neq 1\f$. + UInt_t fNoOfFitBins; ///< number of bins to be be fitted + Int_t fPacking; ///< packing for this particular run. Either given in the RUN- or GLOBAL-block. + + PDoubleVector fForward; ///< forward histo data + PDoubleVector fForwardErr; ///< forward histo errors + PDoubleVector fBackward; ///< backward histo data + PDoubleVector fBackwardErr; ///< backward histo errors + + Int_t fGoodBins[4]; ///< keep first/last good bins. 0=fgb, 1=lgb (forward); 2=fgb, 3=lgb (backward) + + Int_t fStartTimeBin; ///< bin at which the fit starts + Int_t fEndTimeBin; ///< bin at which the fit ends + + Bool_t SubtractFixBkg(); + Bool_t SubtractEstimatedBkg(); + + virtual Bool_t GetProperT0(PRawRunData* runData, PMsrGlobalBlock *globalBlock, PUIntVector &forwardHisto, PUIntVector &backwardHistoNo); + virtual Bool_t GetProperDataRange(PRawRunData* runData, UInt_t histoNo[2]); + virtual void GetProperFitRange(PMsrGlobalBlock *globalBlock); +}; + +#endif // _PRUNASYMMETRYBNMR_H_