From e2a59af6e125b46ef7062a4149fc7b38b6fca101 Mon Sep 17 00:00:00 2001 From: Zaher Salman Date: Wed, 15 Aug 2018 20:35:06 +0200 Subject: [PATCH 01/20] Start branch to develop full featured beta-NMR support --- src/classes/PMsrHandler.cpp | 16 ++- src/classes/PMusrCanvas.cpp | 15 ++- src/classes/PRunListCollection.cpp | 156 +++++++++++++++++++++++++++++ src/include/PMusr.h | 3 + src/include/PRunListCollection.h | 6 ++ 5 files changed, 194 insertions(+), 2 deletions(-) diff --git a/src/classes/PMsrHandler.cpp b/src/classes/PMsrHandler.cpp index a1b20238..8685c8ba 100644 --- a/src/classes/PMsrHandler.cpp +++ b/src/classes/PMsrHandler.cpp @@ -632,6 +632,9 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t messages) case MSR_FITTYPE_MU_MINUS: fout << left << "fittype" << MSR_FITTYPE_MU_MINUS << " (mu minus fit)" << endl ; break; + case MSR_FITTYPE_BNMR: + fout << left << "fittype" << MSR_FITTYPE_BNMR << " (beta-NMR fit)" << endl ; + break; case MSR_FITTYPE_NON_MUSR: fout << left << "fittype" << MSR_FITTYPE_NON_MUSR << " (non muSR fit)" << endl ; break; @@ -796,6 +799,9 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t messages) case MSR_FITTYPE_MU_MINUS: fout << left << "fittype" << MSR_FITTYPE_MU_MINUS << " (mu minus fit)" << endl ; break; + case MSR_FITTYPE_BNMR: + fout << left << "fittype" << MSR_FITTYPE_BNMR << " (beta-NMR fit)" << endl ; + break; case MSR_FITTYPE_NON_MUSR: fout << left << "fittype" << MSR_FITTYPE_NON_MUSR << " (non muSR fit)" << endl ; break; @@ -1705,6 +1711,9 @@ Int_t PMsrHandler::WriteMsrFile(const Char_t *filename, map *co case MSR_FITTYPE_MU_MINUS: fout << left << "fittype" << MSR_FITTYPE_MU_MINUS << " (mu minus fit)" << endl ; break; + case MSR_FITTYPE_BNMR: + fout << left << "fittype" << MSR_FITTYPE_BNMR << " (beta-NMR fit)" << endl ; + break; case MSR_FITTYPE_NON_MUSR: fout << left << "fittype" << MSR_FITTYPE_NON_MUSR << " (non muSR fit)" << endl ; break; @@ -1896,6 +1905,9 @@ Int_t PMsrHandler::WriteMsrFile(const Char_t *filename, map *co case MSR_FITTYPE_MU_MINUS: fout << left << "fittype" << MSR_FITTYPE_MU_MINUS << " (mu minus fit)" << endl ; break; + case MSR_FITTYPE_BNMR: + fout << left << "fittype" << MSR_FITTYPE_BNMR << " (beta-NMR fit)" << endl ; + break; case MSR_FITTYPE_NON_MUSR: fout << left << "fittype" << MSR_FITTYPE_NON_MUSR << " (non muSR fit)" << endl ; break; @@ -2965,6 +2977,7 @@ Bool_t PMsrHandler::HandleGlobalEntry(PMsrLines &lines) (fittype == MSR_FITTYPE_ASYM) || (fittype == MSR_FITTYPE_ASYM_RRF) || (fittype == MSR_FITTYPE_MU_MINUS) || + (fittype == MSR_FITTYPE_BNMR) || (fittype == MSR_FITTYPE_NON_MUSR)) { global.SetFitType(fittype); } else { @@ -3306,6 +3319,7 @@ Bool_t PMsrHandler::HandleRunEntry(PMsrLines &lines) (fittype == MSR_FITTYPE_ASYM) || (fittype == MSR_FITTYPE_ASYM_RRF) || (fittype == MSR_FITTYPE_MU_MINUS) || + (fittype == MSR_FITTYPE_BNMR) || (fittype == MSR_FITTYPE_NON_MUSR)) { param.SetFitType(fittype); } else { @@ -5865,7 +5879,7 @@ Bool_t PMsrHandler::CheckHistoGrouping() Bool_t result = true; for (UInt_t i=0; i> PMsrHandler::CheckHistoGrouping: **ERROR** # of forward histos != # of backward histos."; cerr << endl << ">> Run #" << i+1; diff --git a/src/classes/PMusrCanvas.cpp b/src/classes/PMusrCanvas.cpp index c6240b7c..7e632b0f 100644 --- a/src/classes/PMusrCanvas.cpp +++ b/src/classes/PMusrCanvas.cpp @@ -742,6 +742,18 @@ void PMusrCanvas::UpdateDataTheoryPad() // handle data HandleDataSet(i, runNo, data); break; + case MSR_FITTYPE_BNMR: + data = fRunList->GetAsymmetryBNMR(runNo, PRunListCollection::kRunNo); + if (!data) { // something wrong + fValid = false; + // error message + cerr << endl << ">> PMusrCanvas::UpdateDataTheoryPad(): **ERROR** couldn't obtain run no " << runNo << " for a beta-NMR asymmetry plot"; + cerr << endl; + return; + } + // handle data + HandleDataSet(i, runNo, data); + break; case MSR_FITTYPE_NON_MUSR: data = fRunList->GetNonMusr(runNo, PRunListCollection::kRunNo); if (!data) { // something wrong @@ -910,7 +922,8 @@ void PMusrCanvas::UpdateInfoPad() tstr += grouping; tstr += TString(","); } else if ((runs[runNo].GetFitType() == MSR_FITTYPE_ASYM) || - (runs[runNo].GetFitType() == MSR_FITTYPE_ASYM_RRF)) { + (runs[runNo].GetFitType() == MSR_FITTYPE_ASYM_RRF) || + (runs[runNo].GetFitType() == MSR_FITTYPE_BNMR)) { tstr += TString("h:"); TString grouping; fMsrHandler->GetGroupingString(runNo, "forward", grouping); diff --git a/src/classes/PRunListCollection.cpp b/src/classes/PRunListCollection.cpp index e11bb0c6..dc95aac2 100644 --- a/src/classes/PRunListCollection.cpp +++ b/src/classes/PRunListCollection.cpp @@ -76,6 +76,12 @@ PRunListCollection::~PRunListCollection() } fRunAsymmetryRRFList.clear(); + for (UInt_t i=0; iCleanUp(); + fRunAsymmetryBNMRList[i]->~PRunAsymmetryBNMR(); + } + fRunAsymmetryBNMRList.clear(); + for (UInt_t i=0; iCleanUp(); fRunMuMinusList[i]->~PRunMuMinus(); @@ -133,6 +139,11 @@ Bool_t PRunListCollection::Add(Int_t runNo, EPMusrHandleTag tag) if (!fRunAsymmetryRRFList[fRunAsymmetryRRFList.size()-1]->IsValid()) success = false; break; + case PRUN_ASYMMETRY_BNMR: + fRunAsymmetryBNMRList.push_back(new PRunAsymmetryBNMR(fMsrInfo, fData, runNo, tag)); + if (!fRunAsymmetryBNMRList[fRunAsymmetryBNMRList.size()-1]->IsValid()) + success = false; + break; case PRUN_MU_MINUS: fRunMuMinusList.push_back(new PRunMuMinus(fMsrInfo, fData, runNo, tag)); if (!fRunMuMinusList[fRunMuMinusList.size()-1]->IsValid()) @@ -175,6 +186,8 @@ void PRunListCollection::SetFitRange(const TString fitRange) fRunAsymmetryList[i]->SetFitRangeBin(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; iSetFitRange(fitRange); for (UInt_t i=0; i& p return chisq; } +//-------------------------------------------------------------------------- +// GetAsymmetryBNMRChisq (public) +//-------------------------------------------------------------------------- +/** + *

Calculates chi-square of all asymmetry BNMR runs of a msr-file. + * + * return: + * - chi-square of all asymmetry BNMR runs of the msr-file + * + * \param par fit parameter vector + */ +Double_t PRunListCollection::GetAsymmetryBNMRChisq(const std::vector& par) const +{ + Double_t chisq = 0.0; + + for (UInt_t i=0; iCalcChiSquare(par); + + return chisq; +} + //-------------------------------------------------------------------------- // GetMuMinusChisq (public) //-------------------------------------------------------------------------- @@ -380,6 +416,9 @@ Double_t PRunListCollection::GetSingleHistoChisqExpected(const std::vectorCalcChiSquareExpected(par); break; + case PRUN_ASYMMETRY_BNMR: + expectedChisq = fRunAsymmetryBNMRList[subIdx]->CalcChiSquareExpected(par); + break; case PRUN_MU_MINUS: expectedChisq = fRunMuMinusList[subIdx]->CalcChiSquareExpected(par); break; @@ -441,6 +480,9 @@ Double_t PRunListCollection::GetSingleRunChisq(const std::vector& par, case PRUN_ASYMMETRY_RRF: chisq = fRunAsymmetryRRFList[subIdx]->CalcChiSquare(par); break; + case PRUN_ASYMMETRY_BNMR: + chisq = fRunAsymmetryBNMRList[subIdx]->CalcChiSquare(par); + break; case PRUN_MU_MINUS: chisq = fRunMuMinusList[subIdx]->CalcChiSquare(par); break; @@ -540,6 +582,28 @@ Double_t PRunListCollection::GetAsymmetryRRFMaximumLikelihood(const std::vector< return mlh; } +//-------------------------------------------------------------------------- +// GetAsymmetryBNMRMaximumLikelihood (public) +//-------------------------------------------------------------------------- +/** + *

Since it is not clear yet how to handle asymmetry fits with max likelihood + * the chi square will be used! + * + * return: + * - chi-square of all asymmetry BNMR runs of the msr-file + * + * \param par fit parameter vector + */ +Double_t PRunListCollection::GetAsymmetryBNMRMaximumLikelihood(const std::vector& par) const +{ + Double_t mlh = 0.0; + + for (UInt_t i=0; iCalcChiSquare(par); + + return mlh; +} + //-------------------------------------------------------------------------- // GetMuMinusMaximumLikelihood (public) //-------------------------------------------------------------------------- @@ -721,6 +785,9 @@ UInt_t PRunListCollection::GetNoOfBinsFitted(const UInt_t idx) const case PRUN_ASYMMETRY_RRF: result = fRunAsymmetryRRFList[subIdx]->GetNoOfFitBins(); break; + case PRUN_ASYMMETRY_BNMR: + result = fRunAsymmetryBNMRList[subIdx]->GetNoOfFitBins(); + break; case PRUN_MU_MINUS: result = fRunMuMinusList[subIdx]->GetNoOfFitBins(); break; @@ -760,6 +827,9 @@ UInt_t PRunListCollection::GetTotalNoOfBinsFitted() const for (UInt_t i=0; iGetNoOfFitBins(); + for (UInt_t i=0; iGetNoOfFitBins(); + for (UInt_t i=0; iGetNoOfFitBins(); @@ -898,6 +968,49 @@ PRunData* PRunListCollection::GetAsymmetry(UInt_t index, EDataSwitch tag) return data; } +//-------------------------------------------------------------------------- +// GetAsymmetryBNMR (public) +//-------------------------------------------------------------------------- +/** + *

Get a processed asymmetry 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 > fRunAsymmetryList.size()) { + cerr << endl << ">> PRunListCollection::GetAsymmetryBNMR(): **ERROR** index = " << index << " out of bounds"; + cerr << endl; + return 0; + } + + fRunAsymmetryList[index]->CalcTheory(); + data = fRunAsymmetryList[index]->GetData(); + break; + case kRunNo: // called from PMusrCanvas + for (UInt_t i=0; iGetRunNo() == index) { + data = fRunAsymmetryList[i]->GetData(); + break; + } + } + break; + default: // error + break; + } + + return data; +} + //-------------------------------------------------------------------------- // GetAsymmetryRRF (public) //-------------------------------------------------------------------------- @@ -941,6 +1054,49 @@ 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/PMusr.h b/src/include/PMusr.h index 11884323..cd9d6375 100644 --- a/src/include/PMusr.h +++ b/src/include/PMusr.h @@ -57,6 +57,7 @@ using namespace std; #define PRUN_ASYMMETRY 2 #define PRUN_ASYMMETRY_RRF 3 #define PRUN_MU_MINUS 4 +#define PRUN_ASYMMETRY_BNMR 5 #define PRUN_NON_MUSR 8 // muon life time in (us), see PRL99, 032001 (2007) @@ -99,6 +100,7 @@ using namespace std; #define MSR_FITTYPE_ASYM 2 #define MSR_FITTYPE_ASYM_RRF 3 #define MSR_FITTYPE_MU_MINUS 4 +#define MSR_FITTYPE_BNMR 5 #define MSR_FITTYPE_NON_MUSR 8 //------------------------------------------------------------- @@ -108,6 +110,7 @@ using namespace std; #define MSR_PLOT_ASYM 2 #define MSR_PLOT_ASYM_RRF 3 #define MSR_PLOT_MU_MINUS 4 +#define MSR_PLOT_BNMR 5 #define MSR_PLOT_NON_MUSR 8 //------------------------------------------------------------- diff --git a/src/include/PRunListCollection.h b/src/include/PRunListCollection.h index 1c22ee51..4d3a83ef 100644 --- a/src/include/PRunListCollection.h +++ b/src/include/PRunListCollection.h @@ -40,6 +40,7 @@ using namespace std; #include "PRunSingleHistoRRF.h" #include "PRunAsymmetry.h" #include "PRunAsymmetryRRF.h" +#include "PRunAsymmetryBNMR.h" #include "PRunMuMinus.h" #include "PRunNonMusr.h" @@ -63,6 +64,7 @@ class PRunListCollection virtual Double_t GetSingleHistoRRFChisq(const std::vector& par) const; virtual Double_t GetAsymmetryChisq(const std::vector& par) const; virtual Double_t GetAsymmetryRRFChisq(const std::vector& par) const; + virtual Double_t GetAsymmetryBNMRChisq(const std::vector& par) const; virtual Double_t GetMuMinusChisq(const std::vector& par) const; virtual Double_t GetNonMusrChisq(const std::vector& par) const; @@ -73,6 +75,7 @@ class PRunListCollection virtual Double_t GetSingleHistoRRFMaximumLikelihood(const std::vector& par) const; virtual Double_t GetAsymmetryMaximumLikelihood(const std::vector& par) const; virtual Double_t GetAsymmetryRRFMaximumLikelihood(const std::vector& par) const; + virtual Double_t GetAsymmetryBNMRMaximumLikelihood(const std::vector& par) const; virtual Double_t GetMuMinusMaximumLikelihood(const std::vector& par) const; virtual Double_t GetNonMusrMaximumLikelihood(const std::vector& par) const; @@ -86,6 +89,7 @@ class PRunListCollection 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 GetNoOfAsymmetryRRF() const { return fRunAsymmetryRRFList.size(); } ///< returns the number of asymmetry RRF data sets present in the msr-file + virtual UInt_t GetNoOfAsymmetryBNMR() const { return fRunAsymmetryBNMRList.size(); } ///< returns the number of asymmetry BNMR 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 @@ -93,6 +97,7 @@ class PRunListCollection virtual PRunData* GetSingleHistoRRF(UInt_t index, EDataSwitch tag=kIndex); virtual PRunData* GetAsymmetry(UInt_t index, EDataSwitch tag=kIndex); virtual PRunData* GetAsymmetryRRF(UInt_t index, EDataSwitch tag=kIndex); + virtual PRunData* GetAsymmetryBNMR(UInt_t index, EDataSwitch tag=kIndex); virtual PRunData* GetMuMinus(UInt_t index, EDataSwitch tag=kIndex); virtual PRunData* GetNonMusr(UInt_t index, EDataSwitch tag=kIndex); @@ -111,6 +116,7 @@ class PRunListCollection vector fRunSingleHistoRRFList; ///< stores all processed single histogram RRF data vector fRunAsymmetryList; ///< stores all processed asymmetry data vector fRunAsymmetryRRFList; ///< stores all processed asymmetry RRF data + vector fRunAsymmetryBNMRList; ///< stores all processed asymmetry BNMR data vector fRunMuMinusList; ///< stores all processed mu-minus data vector fRunNonMusrList; ///< stores all processed non-muSR data }; From 4e54d9c0b966e5160d88d87cc061c86933fd4de6 Mon Sep 17 00:00:00 2001 From: Zaher Salman Date: Wed, 15 Aug 2018 21:35:56 +0200 Subject: [PATCH 02/20] First successful compilation of AsymmetryBNMR --- src/classes/CMakeLists.txt | 1 + src/classes/Makefile.am | 2 + src/classes/PRunAsymmetryBNMR.cpp | 1869 ++++++++++++++++++++++++++++ src/classes/PRunListCollection.cpp | 43 - src/include/PRunAsymmetryBNMR.h | 89 ++ 5 files changed, 1961 insertions(+), 43 deletions(-) create mode 100644 src/classes/PRunAsymmetryBNMR.cpp create mode 100644 src/include/PRunAsymmetryBNMR.h diff --git a/src/classes/CMakeLists.txt b/src/classes/CMakeLists.txt index f64feace..40eb0378 100644 --- a/src/classes/CMakeLists.txt +++ b/src/classes/CMakeLists.txt @@ -65,6 +65,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 index 77cb1d8f..fa3a0fba 100644 --- a/src/classes/Makefile.am +++ b/src/classes/Makefile.am @@ -15,6 +15,7 @@ h_sources = \ ../include/PMusrT0.h \ ../include/PPrepFourier.h \ ../include/PRunAsymmetry.h \ + ../include/PRunAsymmetryBNMR.h \ ../include/PRunAsymmetryRRF.h \ ../include/PRunBase.h \ ../include/PRunDataHandler.h \ @@ -61,6 +62,7 @@ cpp_sources = \ PMusrT0.cpp \ PPrepFourier.cpp \ PRunAsymmetry.cpp \ + PRunAsymmetryBNMR.cpp \ PRunAsymmetryRRF.cpp \ PRunBase.cpp \ PRunDataHandler.cpp \ 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 dc95aac2..2419521e 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_ From b044d97a55db70074867bbe1b77b4612de836da9 Mon Sep 17 00:00:00 2001 From: Zaher Salman Date: Thu, 16 Aug 2018 23:17:39 +0200 Subject: [PATCH 03/20] More work towards implementation of beta-NMR asymmetry calculation --- src/classes/PMsrHandler.cpp | 60 +++++++++++++++++++++++++++++ src/classes/PMusrCanvas.cpp | 11 ++++++ src/classes/PRunAsymmetryBNMR.cpp | 1 + src/external/libBNMR/ExpRlx-MUD.msr | 6 +-- src/include/PMusrCanvas.h | 2 +- 5 files changed, 76 insertions(+), 4 deletions(-) diff --git a/src/classes/PMsrHandler.cpp b/src/classes/PMsrHandler.cpp index 8685c8ba..031c0088 100644 --- a/src/classes/PMsrHandler.cpp +++ b/src/classes/PMsrHandler.cpp @@ -1162,6 +1162,9 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t messages) case MSR_PLOT_MU_MINUS: fout << "PLOT " << fPlots[plotNo].fPlotType << " (mu minus plot)" << endl; break; + case MSR_PLOT_BNMR: + fout << "PLOT " << fPlots[plotNo].fPlotType << " (beta-NMR asymmetry plot)" << endl; + break; case MSR_PLOT_NON_MUSR: fout << "PLOT " << fPlots[plotNo].fPlotType << " (non muSR plot)" << endl; break; @@ -2259,6 +2262,9 @@ Int_t PMsrHandler::WriteMsrFile(const Char_t *filename, map *co case MSR_PLOT_MU_MINUS: fout << "PLOT " << fPlots[i].fPlotType << " (mu minus plot)" << endl; break; + case MSR_PLOT_BNMR: + fout << "PLOT " << fPlots[i].fPlotType << " (beta-NMR asymmetry plot)" << endl; + break; case MSR_PLOT_NON_MUSR: fout << "PLOT " << fPlots[i].fPlotType << " (non muSR plot)" << endl; break; @@ -4264,6 +4270,7 @@ Bool_t PMsrHandler::HandlePlotEntry(PMsrLines &lines) case MSR_PLOT_SINGLE_HISTO: // like: runs 1 5 13 case MSR_PLOT_SINGLE_HISTO_RRF: case MSR_PLOT_ASYM: + case MSR_PLOT_BNMR: case MSR_PLOT_ASYM_RRF: case MSR_PLOT_NON_MUSR: case MSR_PLOT_MU_MINUS: @@ -4681,6 +4688,7 @@ Bool_t PMsrHandler::HandlePlotEntry(PMsrLines &lines) cerr << endl << ">> 2=forward-backward asym,"; cerr << endl << ">> 3=forward-backward RRF asym,"; cerr << endl << ">> 4=mu minus single histo,"; + cerr << endl << ">> 5=forward-backward beta-NMR asym,"; cerr << endl << ">> 8=non muSR."; cerr << endl << ">> is the list of runs, e.g. runs 1 3"; cerr << endl << ">> range is optional"; @@ -5532,6 +5540,58 @@ Bool_t PMsrHandler::CheckRunBlockIntegrity() fRuns[i].SetPacking(1); } break; + case PRUN_ASYMMETRY_BNMR: + // check alpha + if ((fRuns[i].GetAlphaParamNo() == -1) && !fFourierOnly) { + cerr << endl << ">> PMsrHandler::CheckRunBlockIntegrity(): **ERROR** in RUN block number " << i+1; + cerr << endl << ">> alpha parameter number missing which is needed for an asymmetry fit."; + cerr << endl << ">> Consider to check the manual ;-)" << endl; + return false; + } + // check that there is a forward parameter number + if (fRuns[i].GetForwardHistoNo() == -1) { + cerr << endl << ">> PMsrHandler::CheckRunBlockIntegrity(): **ERROR** in RUN block number " << i+1; + cerr << endl << ">> forward histogram number not defined. Necessary for asymmetry fits." << endl; + return false; + } + // check that there is a backward parameter number + if (fRuns[i].GetBackwardHistoNo() == -1) { + cerr << endl << ">> PMsrHandler::CheckRunBlockIntegrity(): **ERROR** in RUN block number " << i+1; + cerr << endl << ">> backward histogram number not defined. Necessary for asymmetry fits." << endl; + return false; + } + // check fit range + if (!fRuns[i].IsFitRangeInBin()) { // fit range given as times in usec + if ((fRuns[i].GetFitRange(0) == PMUSR_UNDEFINED) || (fRuns[i].GetFitRange(1) == PMUSR_UNDEFINED)) { + if ((fGlobal.GetFitRange(0) == PMUSR_UNDEFINED) || (fGlobal.GetFitRange(1) == PMUSR_UNDEFINED)) { + cerr << endl << ">> PMsrHandler::CheckRunBlockIntegrity(): **ERROR** in RUN block number " << i+1; + cerr << endl << ">> Fit range is not defined, also NOT present in the GLOBAL block. Necessary for asymmetry fits." << endl; + return false; + } + } + } + // check number of T0's provided + if ((fRuns[i].GetT0BinSize() > 2*fRuns[i].GetForwardHistoNoSize()) && + (fGlobal.GetT0BinSize() > 2*fRuns[i].GetForwardHistoNoSize())) { + cerr << endl << ">> PMsrHandler::CheckRunBlockIntegrity(): **ERROR** in RUN block number " << i+1; + cerr << endl << ">> Found " << fRuns[i].GetT0BinSize() << " T0 entries. Expecting only " << 2*fRuns[i].GetForwardHistoNoSize() << " in forward. Needs to be fixed." << endl; + cerr << endl << ">> In GLOBAL block: " << fGlobal.GetT0BinSize() << " T0 entries. Expecting only " << 2*fRuns[i].GetForwardHistoNoSize() << ". Needs to be fixed." << endl; + return false; + } + if ((fRuns[i].GetT0BinSize() > 2*fRuns[i].GetBackwardHistoNoSize()) && + (fGlobal.GetT0BinSize() > 2*fRuns[i].GetBackwardHistoNoSize())) { + cerr << endl << ">> PMsrHandler::CheckRunBlockIntegrity(): **ERROR** in RUN block number " << i+1; + cerr << endl << ">> Found " << fRuns[i].GetT0BinSize() << " T0 entries. Expecting only " << 2*fRuns[i].GetBackwardHistoNoSize() << " in backward. Needs to be fixed." << endl; + cerr << endl << ">> In GLOBAL block: " << fGlobal.GetT0BinSize() << " T0 entries. Expecting only " << 2*fRuns[i].GetBackwardHistoNoSize() << ". Needs to be fixed." << endl; + return false; + } + // check packing + if ((fRuns[i].GetPacking() == -1) && (fGlobal.GetPacking() == -1)) { + cerr << endl << ">> PMsrHandler::CheckRunBlockIntegrity(): **WARNING** in RUN block number " << i+1; + cerr << endl << ">> Packing is neither defined here, nor in the GLOBAL block, will set it to 1." << endl; + fRuns[i].SetPacking(1); + } + break; case PRUN_ASYMMETRY_RRF: // check alpha if ((fRuns[i].GetAlphaParamNo() == -1) && !fFourierOnly) { diff --git a/src/classes/PMusrCanvas.cpp b/src/classes/PMusrCanvas.cpp index 7e632b0f..f94af908 100644 --- a/src/classes/PMusrCanvas.cpp +++ b/src/classes/PMusrCanvas.cpp @@ -1618,6 +1618,7 @@ void PMusrCanvas::ExportData(const Char_t *fileName) case MSR_PLOT_SINGLE_HISTO: case MSR_PLOT_SINGLE_HISTO_RRF: case MSR_PLOT_ASYM: + case MSR_PLOT_BNMR: case MSR_PLOT_ASYM_RRF: case MSR_PLOT_MU_MINUS: if (fDifferenceView) { // difference view plot @@ -2778,6 +2779,8 @@ void PMusrCanvas::HandleDataSet(UInt_t plotNo, UInt_t runNo, PRunData *data) // make sure that for asymmetry the y-range is initialized reasonably if (fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fPlotType == MSR_PLOT_ASYM) dataSet.dataRange->SetYRange(-0.4, 0.4); + if (fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fPlotType == MSR_PLOT_BNMR) + dataSet.dataRange->SetYRange(-0.4, 0.4); // extract necessary range information if ((fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fTmin.size() == 0) && !fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fUseFitRanges) { // no range information at all @@ -2792,6 +2795,7 @@ void PMusrCanvas::HandleDataSet(UInt_t plotNo, UInt_t runNo, PRunData *data) fXmax = end; } if ((fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fPlotType == MSR_PLOT_ASYM) || + (fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fPlotType == MSR_PLOT_BNMR) || (fMsrHandler->GetMsrRunList()->at(runNo).IsLifetimeCorrected())) { fYRangePresent = true; fYmin = -0.4; @@ -2840,6 +2844,7 @@ void PMusrCanvas::HandleDataSet(UInt_t plotNo, UInt_t runNo, PRunData *data) // make sure that for asymmetry the y-range is initialized reasonably if ((fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fPlotType == MSR_PLOT_ASYM) || + (fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fPlotType == MSR_PLOT_BNMR) || (fMsrHandler->GetMsrRunList()->at(runNo).IsLifetimeCorrected())) { dataSet.dataRange->SetYRange(-0.4, 0.4); } @@ -4714,6 +4719,9 @@ void PMusrCanvas::PlotData(Bool_t unzoom) case MSR_PLOT_ASYM: yAxisTitle = "Asymmetry"; break; + case MSR_PLOT_BNMR: + yAxisTitle = "Asymmetry"; + break; case MSR_PLOT_MU_MINUS: yAxisTitle = "N(t) per bin"; break; @@ -6136,6 +6144,9 @@ void PMusrCanvas::PlotAverage(Bool_t unzoom) case MSR_PLOT_ASYM: yAxisTitle = ""; break; + case MSR_PLOT_BNMR: + yAxisTitle = ""; + break; case MSR_PLOT_MU_MINUS: yAxisTitle = " per bin"; break; diff --git a/src/classes/PRunAsymmetryBNMR.cpp b/src/classes/PRunAsymmetryBNMR.cpp index b55052a0..a3ff5a29 100644 --- a/src/classes/PRunAsymmetryBNMR.cpp +++ b/src/classes/PRunAsymmetryBNMR.cpp @@ -535,6 +535,7 @@ Bool_t PRunAsymmetryBNMR::PrepareData() } // keep the time resolution in (us) + // possibility to rescale for betaNMR fTimeResolution = runData->GetTimeResolution()/1.0e3; cout.precision(10); cout << endl << ">> PRunAsymmetryBNMR::PrepareData(): time resolution=" << fixed << runData->GetTimeResolution() << "(ns)" << endl; diff --git a/src/external/libBNMR/ExpRlx-MUD.msr b/src/external/libBNMR/ExpRlx-MUD.msr index ef423491..404a46ef 100644 --- a/src/external/libBNMR/ExpRlx-MUD.msr +++ b/src/external/libBNMR/ExpRlx-MUD.msr @@ -24,7 +24,7 @@ fun1 = 0.5 * map1 * map2 ############################################################### RUN 045674 BNMR TRIUMF MUD (name beamline institute data-file-format) -fittype 2 (asymmetry fit) +fittype 5 (asymmetry fit) alpha 1 forward 3 backward 4 @@ -37,7 +37,7 @@ fit 0.5 8 packing 5 RUN 045674 BNMR TRIUMF MUD (name beamline institute data-file-format) -fittype 2 (asymmetry fit) +fittype 5 (asymmetry fit) alpha 1 forward 5 backward 6 @@ -57,7 +57,7 @@ HESSE SAVE ############################################################### -PLOT 2 (asymmetry plot) +PLOT 5 (asymmetry plot) runs 1 2 use_fit_ranges view_packing 10 diff --git a/src/include/PMusrCanvas.h b/src/include/PMusrCanvas.h index a3b1cffa..f50b8b57 100644 --- a/src/include/PMusrCanvas.h +++ b/src/include/PMusrCanvas.h @@ -246,7 +246,7 @@ class PMusrCanvas : public TObject, public TQObject Bool_t fDifferenceView; ///< tag showing that the shown data, fourier, are the difference between data and theory Int_t fCurrentPlotView; ///< tag showing what the current plot view is: data, fourier, ... Int_t fPreviousPlotView; ///< tag showing the previous plot view - Int_t fPlotType; ///< plot type tag: -1 == undefined, MSR_PLOT_SINGLE_HISTO == single histogram, MSR_PLOT_ASYM == asymmetry, MSR_PLOT_MU_MINUS == mu minus (not yet implemented), MSR_PLOT_NON_MUSR == non-muSR + Int_t fPlotType; ///< plot type tag: -1 == undefined, MSR_PLOT_SINGLE_HISTO == single histogram, MSR_PLOT_ASYM == asymmetry, MSR_PLOT_BNMR == beta-NMR asymmetry, MSR_PLOT_MU_MINUS == mu minus (not yet implemented), MSR_PLOT_NON_MUSR == non-muSR Int_t fPlotNumber; ///< plot number Bool_t fXRangePresent, fYRangePresent; ///< flag indicating if x-/y-range is present From 758f338972b412296d7cfb503255ca37cdada75d Mon Sep 17 00:00:00 2001 From: Zaher Salman Date: Fri, 17 Aug 2018 16:08:04 +0200 Subject: [PATCH 04/20] First working duplicate asymmetry class for beta-NMR. --- src/classes/PMusrCanvas.cpp | 4 +- src/classes/PRunListCollection.cpp | 12 +-- src/musrt0.cpp | 164 +++++++++++++++++++++++++++++ 3 files changed, 171 insertions(+), 9 deletions(-) diff --git a/src/classes/PMusrCanvas.cpp b/src/classes/PMusrCanvas.cpp index f94af908..60d067a2 100644 --- a/src/classes/PMusrCanvas.cpp +++ b/src/classes/PMusrCanvas.cpp @@ -2777,9 +2777,7 @@ void PMusrCanvas::HandleDataSet(UInt_t plotNo, UInt_t runNo, PRunData *data) size = data->GetValue()->size(); dataSet.dataRange->SetXRange(start, end); // full possible range // make sure that for asymmetry the y-range is initialized reasonably - if (fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fPlotType == MSR_PLOT_ASYM) - dataSet.dataRange->SetYRange(-0.4, 0.4); - if (fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fPlotType == MSR_PLOT_BNMR) + if ((fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fPlotType == MSR_PLOT_ASYM) || (fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fPlotType == MSR_PLOT_BNMR)) dataSet.dataRange->SetYRange(-0.4, 0.4); // extract necessary range information if ((fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fTmin.size() == 0) && diff --git a/src/classes/PRunListCollection.cpp b/src/classes/PRunListCollection.cpp index 2419521e..6733062c 100644 --- a/src/classes/PRunListCollection.cpp +++ b/src/classes/PRunListCollection.cpp @@ -987,19 +987,19 @@ PRunData* PRunListCollection::GetAsymmetryBNMR(UInt_t index, EDataSwitch tag) switch (tag) { case kIndex: // called from musrfit when dumping the data - if (index > fRunAsymmetryList.size()) { + if (index > fRunAsymmetryBNMRList.size()) { cerr << endl << ">> PRunListCollection::GetAsymmetryBNMR(): **ERROR** index = " << index << " out of bounds"; cerr << endl; return 0; } - fRunAsymmetryList[index]->CalcTheory(); - data = fRunAsymmetryList[index]->GetData(); + fRunAsymmetryBNMRList[index]->CalcTheory(); + data = fRunAsymmetryBNMRList[index]->GetData(); break; case kRunNo: // called from PMusrCanvas - for (UInt_t i=0; iGetRunNo() == index) { - data = fRunAsymmetryList[i]->GetData(); + for (UInt_t i=0; iGetRunNo() == index) { + data = fRunAsymmetryBNMRList[i]->GetData(); break; } } diff --git a/src/musrt0.cpp b/src/musrt0.cpp index dab51fa7..d300ef7a 100644 --- a/src/musrt0.cpp +++ b/src/musrt0.cpp @@ -636,6 +636,169 @@ Int_t main(Int_t argc, Char_t *argv[]) } } } + break; + case MSR_FITTYPE_BNMR: + if ((runList->at(i).GetRunNameSize() == 1) && (runList->at(i).GetForwardHistoNoSize() == 1)) { // no addruns / no grouping + // handle forward histo + // get histo number + histoNo = runList->at(i).GetForwardHistoNo(); + runName = runList->at(i).GetRunName(); + // get bin position of maximal data + t0Bin = musrt0_getMaxBin(dataHandler->GetRunData(*runName)->GetDataBin(histoNo)); + // set t0 to maximum data position + runList->at(i).SetT0Bin(t0Bin, 0); + // set data range as well if firstGoodBinOffset is given + if (firstGoodBinOffsetPresent) { + start = t0Bin + firstGoodBinOffset; + end = dataHandler->GetRunData(*runName)->GetDataBin(histoNo)->size(); + runList->at(i).SetDataRange(start, 0); + runList->at(i).SetDataRange(end, 1); + } + // handle backward histo + // get histo number + histoNo = runList->at(i).GetBackwardHistoNo(); + runName = runList->at(i).GetRunName(); + // get bin position of maximal data + t0Bin = musrt0_getMaxBin(dataHandler->GetRunData(*runName)->GetDataBin(histoNo)); + // set t0 to maximum data position + runList->at(i).SetT0Bin(t0Bin, 1); + // set data range as well if firstGoodBinOffset is given + if (firstGoodBinOffsetPresent) { + start = t0Bin + firstGoodBinOffset; + end = dataHandler->GetRunData(*runName)->GetDataBin(histoNo)->size(); + runList->at(i).SetDataRange(start, 2); + runList->at(i).SetDataRange(end, 3); + } + } else if ((runList->at(i).GetRunNameSize() > 1) && (runList->at(i).GetForwardHistoNoSize() == 1)) { // addruns / no grouping + // handle forward histo + // get histo number + histoNo = runList->at(i).GetForwardHistoNo(); + runName = runList->at(i).GetRunName(); + // get bin position of maximal data + t0Bin = musrt0_getMaxBin(dataHandler->GetRunData(*runName)->GetDataBin(histoNo)); + // set t0 to maximum data position + runList->at(i).SetT0Bin(t0Bin, 0); + // set data range as well if firstGoodBinOffset is given + if (firstGoodBinOffsetPresent) { + start = t0Bin + firstGoodBinOffset; + end = dataHandler->GetRunData(*runName)->GetDataBin(histoNo)->size(); + runList->at(i).SetDataRange(start, 0); + runList->at(i).SetDataRange(end, 1); + } + // handle addruns + for (UInt_t j=1; jat(i).GetRunNameSize(); j++) { + runName = runList->at(i).GetRunName(j); + // get bin position of maximal data + t0Bin = musrt0_getMaxBin(dataHandler->GetRunData(*runName)->GetDataBin(histoNo)); + // set t0 to maximum data position + runList->at(i).SetAddT0Bin(t0Bin, j-1, 0); + } + // handle backward histo + // get histo number + histoNo = runList->at(i).GetBackwardHistoNo(); + runName = runList->at(i).GetRunName(); + // get bin position of maximal data + t0Bin = musrt0_getMaxBin(dataHandler->GetRunData(*runName)->GetDataBin(histoNo)); + // set t0 to maximum data position + runList->at(i).SetT0Bin(t0Bin, 1); + // set data range as well if firstGoodBinOffset is given + if (firstGoodBinOffsetPresent) { + start = t0Bin + firstGoodBinOffset; + end = dataHandler->GetRunData(*runName)->GetDataBin(histoNo)->size(); + runList->at(i).SetDataRange(start, 2); + runList->at(i).SetDataRange(end, 3); + } + // handle addruns + for (UInt_t j=1; jat(i).GetRunNameSize(); j++) { + runName = runList->at(i).GetRunName(j); + // get bin position of maximal data + t0Bin = musrt0_getMaxBin(dataHandler->GetRunData(*runName)->GetDataBin(histoNo)); + // set t0 to maximum data position + runList->at(i).SetAddT0Bin(t0Bin, j-1, 1); + } + } else if ((runList->at(i).GetRunNameSize() == 1) && (runList->at(i).GetForwardHistoNoSize() > 1)) { // no addruns / grouping + // handle forward histo + for (UInt_t j=0; jat(i).GetForwardHistoNoSize(); j++) { + // get histo number + histoNo = runList->at(i).GetForwardHistoNo(j); + runName = runList->at(i).GetRunName(); + // get bin position of maximal data + t0Bin = musrt0_getMaxBin(dataHandler->GetRunData(*runName)->GetDataBin(histoNo)); + // set t0 to maximum data position + runList->at(i).SetT0Bin(t0Bin, 2*j); + if (firstGoodBinOffsetPresent && (j==0)) { + start = t0Bin + firstGoodBinOffset; + end = dataHandler->GetRunData(*runName)->GetDataBin(histoNo)->size(); + runList->at(i).SetDataRange(start, 0); + runList->at(i).SetDataRange(end, 1); + } + } + // handle backward histo + for (UInt_t j=0; jat(i).GetBackwardHistoNoSize(); j++) { + // get histo number + histoNo = runList->at(i).GetBackwardHistoNo(j); + runName = runList->at(i).GetRunName(); + // get bin position of maximal data + t0Bin = musrt0_getMaxBin(dataHandler->GetRunData(*runName)->GetDataBin(histoNo)); + // set t0 to maximum data position + runList->at(i).SetT0Bin(t0Bin, 2*j+1); + if (firstGoodBinOffsetPresent && (j==0)) { + start = t0Bin + firstGoodBinOffset; + end = dataHandler->GetRunData(*runName)->GetDataBin(histoNo)->size(); + runList->at(i).SetDataRange(start, 2); + runList->at(i).SetDataRange(end, 3); + } + } + } else { // addruns / grouping + // handle forward histo + for (UInt_t j=0; jat(i).GetForwardHistoNoSize(); j++) { + // get histo number + histoNo = runList->at(i).GetForwardHistoNo(j); + runName = runList->at(i).GetRunName(); + // get bin position of maximal data + t0Bin = musrt0_getMaxBin(dataHandler->GetRunData(*runName)->GetDataBin(histoNo)); + // set t0 to maximum data position + runList->at(i).SetT0Bin(t0Bin, 2*j); + if (firstGoodBinOffsetPresent && (j==0)) { + start = t0Bin + firstGoodBinOffset; + end = dataHandler->GetRunData(*runName)->GetDataBin(histoNo)->size(); + runList->at(i).SetDataRange(start, 0); + runList->at(i).SetDataRange(end, 1); + } + // handle addruns + for (UInt_t k=1; kat(i).GetRunNameSize(); k++) { + runName = runList->at(i).GetRunName(k); + // get bin position of maximal data + t0Bin = musrt0_getMaxBin(dataHandler->GetRunData(*runName)->GetDataBin(histoNo)); + // set t0 to maximum data position + runList->at(i).SetAddT0Bin(t0Bin, k-1, 2*j); + } + } + // handle backward histo + for (UInt_t j=0; jat(i).GetBackwardHistoNoSize(); j++) { + // get histo number + histoNo = runList->at(i).GetBackwardHistoNo(j); + runName = runList->at(i).GetRunName(); + // get bin position of maximal data + t0Bin = musrt0_getMaxBin(dataHandler->GetRunData(*runName)->GetDataBin(histoNo)); + // set t0 to maximum data position + runList->at(i).SetT0Bin(t0Bin, 2*j+1); + if (firstGoodBinOffsetPresent && (j==0)) { + start = t0Bin + firstGoodBinOffset; + end = dataHandler->GetRunData(*runName)->GetDataBin(histoNo)->size(); + runList->at(i).SetDataRange(start, 2); + runList->at(i).SetDataRange(end, 3); + } + // handle addruns + for (UInt_t k=1; kat(i).GetRunNameSize(); k++) { + runName = runList->at(i).GetRunName(k); + // get bin position of maximal data + t0Bin = musrt0_getMaxBin(dataHandler->GetRunData(*runName)->GetDataBin(histoNo)); + // set t0 to maximum data position + runList->at(i).SetAddT0Bin(t0Bin, k-1, 2*j+1); + } + } + } break; default: break; @@ -797,6 +960,7 @@ Int_t main(Int_t argc, Char_t *argv[]) } break; case MSR_FITTYPE_ASYM: + case MSR_FITTYPE_BNMR: case MSR_FITTYPE_ASYM_RRF: if ((runList->at(i).GetRunNameSize() == 1) && (runList->at(i).GetForwardHistoNoSize() == 1)) { // no addruns / no grouping // feed necessary data forward From e5c1450fc87d0cadc953d93740700bfe1824bd87 Mon Sep 17 00:00:00 2001 From: Zaher Salman Date: Sat, 18 Aug 2018 23:45:28 +0200 Subject: [PATCH 05/20] Fixed chisqr calculation bug in BNMR --- src/classes/PFitterFcn.cpp | 2 ++ src/external/libBNMR/ExpRlx-MUD.msr | 10 +++++----- 2 files changed, 7 insertions(+), 5 deletions(-) diff --git a/src/classes/PFitterFcn.cpp b/src/classes/PFitterFcn.cpp index 07c4831a..378fd55d 100644 --- a/src/classes/PFitterFcn.cpp +++ b/src/classes/PFitterFcn.cpp @@ -80,6 +80,7 @@ Double_t PFitterFcn::operator()(const std::vector& par) const value += fRunListCollection->GetSingleHistoRRFChisq(par); value += fRunListCollection->GetAsymmetryChisq(par); value += fRunListCollection->GetAsymmetryRRFChisq(par); + value += fRunListCollection->GetAsymmetryBNMRChisq(par); value += fRunListCollection->GetMuMinusChisq(par); value += fRunListCollection->GetNonMusrChisq(par); } else { // max likelihood @@ -87,6 +88,7 @@ Double_t PFitterFcn::operator()(const std::vector& par) const value += fRunListCollection->GetSingleHistoRRFMaximumLikelihood(par); value += fRunListCollection->GetAsymmetryMaximumLikelihood(par); value += fRunListCollection->GetAsymmetryRRFMaximumLikelihood(par); + value += fRunListCollection->GetAsymmetryBNMRMaximumLikelihood(par); value += fRunListCollection->GetMuMinusMaximumLikelihood(par); value += fRunListCollection->GetNonMusrMaximumLikelihood(par); } diff --git a/src/external/libBNMR/ExpRlx-MUD.msr b/src/external/libBNMR/ExpRlx-MUD.msr index 404a46ef..da4e9a70 100644 --- a/src/external/libBNMR/ExpRlx-MUD.msr +++ b/src/external/libBNMR/ExpRlx-MUD.msr @@ -24,7 +24,7 @@ fun1 = 0.5 * map1 * map2 ############################################################### RUN 045674 BNMR TRIUMF MUD (name beamline institute data-file-format) -fittype 5 (asymmetry fit) +fittype 5 (beta-NMR fit) alpha 1 forward 3 backward 4 @@ -37,7 +37,7 @@ fit 0.5 8 packing 5 RUN 045674 BNMR TRIUMF MUD (name beamline institute data-file-format) -fittype 5 (asymmetry fit) +fittype 5 (beta-NMR fit) alpha 1 forward 5 backward 6 @@ -57,7 +57,7 @@ HESSE SAVE ############################################################### -PLOT 5 (asymmetry plot) +PLOT 5 (beta-NMR asymmetry plot) runs 1 2 use_fit_ranges view_packing 10 @@ -68,9 +68,9 @@ FOURIER units MHz # units either 'Gauss', 'Tesla', 'MHz', or 'Mc/s' fourier_power 12 apodization STRONG # NONE, WEAK, MEDIUM, STRONG -plot POWER # REAL, IMAG, REAL_AND_IMAG, POWER, PHASE +plot POWER # REAL, IMAG, REAL_AND_IMAG, POWER, PHASE, PHASE_OPT_REAL phase 8 #range FRQMIN FRQMAX ############################################################### -STATISTIC --- 2015-04-14 11:00:19 +STATISTIC --- 2018-08-18 23:43:35 chisq = 399.5, NDF = 290, chisq/NDF = 1.377736 From 782451e825694cec49b4c3144979a45bc84d1649 Mon Sep 17 00:00:00 2001 From: Zaher Salman Date: Sun, 19 Aug 2018 23:37:48 +0200 Subject: [PATCH 06/20] Start substract asymmetry implementation. --- src/classes/PRunAsymmetryBNMR.cpp | 738 +++++++++++------------------- src/include/PRunAsymmetryBNMR.h | 12 +- 2 files changed, 286 insertions(+), 464 deletions(-) diff --git a/src/classes/PRunAsymmetryBNMR.cpp b/src/classes/PRunAsymmetryBNMR.cpp index a3ff5a29..07bb2c3c 100644 --- a/src/classes/PRunAsymmetryBNMR.cpp +++ b/src/classes/PRunAsymmetryBNMR.cpp @@ -164,10 +164,15 @@ PRunAsymmetryBNMR::PRunAsymmetryBNMR(PMsrHandler *msrInfo, PRunDataHandler *rawD */ PRunAsymmetryBNMR::~PRunAsymmetryBNMR() { - fForward.clear(); - fForwardErr.clear(); - fBackward.clear(); - fBackwardErr.clear(); + fForwardp.clear(); + fForwardpErr.clear(); + fBackwardp.clear(); + fBackwardpErr.clear(); + + fForwardm.clear(); + fForwardmErr.clear(); + fBackwardm.clear(); + fBackwardmErr.clear(); } //-------------------------------------------------------------------------- @@ -595,32 +600,34 @@ Bool_t PRunAsymmetryBNMR::PrepareData() } // 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 forward histograms of the group + // for (UInt_t i=1; 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]]; - } - } - } + // // 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 @@ -663,10 +670,7 @@ Bool_t PRunAsymmetryBNMR::PrepareData() status = PrepareFitData(); break; case kView: - if (fMsrInfo->GetMsrPlotList()->at(0).fRRFPacking == 0) - status = PrepareViewData(runData, histoNo); - else - status = PrepareRRFViewData(runData, histoNo); + status = PrepareViewData(runData, histoNo); break; default: status = false; @@ -783,64 +787,90 @@ Bool_t PRunAsymmetryBNMR::SubtractEstimatedBkg() } // check if start is within histogram bounds - if ((start[0] < 0) || (start[0] >= fForward.size()) || - (start[1] < 0) || (start[1] >= fBackward.size())) { + if ((start[0] < 0) || (start[0] >= fForwardp.size()) || + (start[1] < 0) || (start[1] >= fBackwardp.size())) { cerr << endl << ">> PRunAsymmetryBNMR::SubtractEstimatedBkg(): **ERROR** background bin values out of bound!"; - cerr << endl << ">> histo lengths (f/b) = (" << fForward.size() << "/" << fBackward.size() << ")."; + cerr << endl << ">> histo lengths (f/b) = (" << fForwardp.size() << "/" << fBackwardp.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())) { + if ((end[0] < 0) || (end[0] >= fForwardp.size()) || + (end[1] < 0) || (end[1] >= fBackwardp.size())) { cerr << endl << ">> PRunAsymmetryBNMR::SubtractEstimatedBkg(): **ERROR** background bin values out of bound!"; - cerr << endl << ">> histo lengths (f/b) = (" << fForward.size() << "/" << fBackward.size() << ")."; + cerr << endl << ">> histo lengths (f/b) = (" << fForwardp.size() << "/" << fBackwardp.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}; + Double_t bkgp[2] = {0.0, 0.0}; + Double_t errBkgp[2] = {0.0, 0.0}; + Double_t bkgn[2] = {0.0, 0.0}; + Double_t errBkgn[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]; + for (UInt_t i=start[0]; i<=end[0]; i++) { + bkgp[0] += fForwardp[i]; + bkgm[0] += fForwardm[i]; + } + errBkgp[0] = TMath::Sqrt(bkgp[0])/(end[0] - start[0] + 1); + bkgp[0] /= static_cast(end[0] - start[0] + 1); + cout << endl << ">> estimated pos hel forward histo background: " << bkgp[0]; + errBkgm[0] = TMath::Sqrt(bkgp[0])/(end[0] - start[0] + 1); + bkgm[0] /= static_cast(end[0] - start[0] + 1); + cout << endl << ">> estimated neg hel forward histo background: " << bkgm[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; + for (UInt_t i=start[1]; i<=end[1]; i++) { + bkgp[1] += fBackwardp[i]; + bkgm[1] += fBackwardm[i]; + } + errBkgp[1] = TMath::Sqrt(bkgp[1])/(end[1] - start[1] + 1); + bkgp[1] /= static_cast(end[1] - start[1] + 1); + cout << endl << ">> estimated pos hel backward histo background: " << bkgp[1] << endl; + errBkgm[1] = TMath::Sqrt(bkgm[1])/(end[1] - start[1] + 1); + bkgm[1] /= static_cast(end[1] - start[1] + 1); + cout << endl << ">> estimated neg hel backward histo background: " << bkgm[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]); + for (UInt_t i=0; i 0.0) + errVal = TMath::Sqrt(fForwardp[i]+errBkgp[0]*errBkgp[0]); else errVal = 1.0; - fForwardErr.push_back(errVal); - if (fBackward[i] > 0.0) - errVal = TMath::Sqrt(fBackward[i]+errBkg[1]*errBkg[1]); + fForwardpErr.push_back(errVal); + if (fBackwardp[i] > 0.0) + errVal = TMath::Sqrt(fBackwardp[i]+errBkgp[1]*errBkgp[1]); else errVal = 1.0; - fBackwardErr.push_back(errVal); + fBackwardpErr.push_back(errVal); + if (fForwardm[i] > 0.0) + errVal = TMath::Sqrt(fForwardm[i]+errBkgm[0]*errBkgm[0]); + else + errVal = 1.0; + fForwardmErr.push_back(errVal); + if (fBackwardm[i] > 0.0) + errVal = TMath::Sqrt(fBackwardm[i]+errBkgm[1]*errBkgm[1]); + else + errVal = 1.0; + fBackwardmErr.push_back(errVal); } // subtract background from data - for (UInt_t i=0; iSetBkgEstimated(bkg[0], 0); - fRunInfo->SetBkgEstimated(bkg[1], 1); + fRunInfo->SetBkgEstimated(bkgp[0], 0); + fRunInfo->SetBkgEstimated(bkgp[1], 1); + fRunInfo->SetBkgEstimated(bkgm[0], 3); + fRunInfo->SetBkgEstimated(bkgm[1], 4); return true; } @@ -865,96 +895,143 @@ Bool_t PRunAsymmetryBNMR::PrepareFitData() // 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; + PRunData forwardpPacked; + PRunData backwardpPacked; + PRunData forwardmPacked; + PRunData backwarmpPacked; + Double_t valuep = 0.0; + Double_t errorp = 0.0; + Double_t valuem = 0.0; + Double_t errorm = 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; + valuep /= fPacking; + valuem /= fPacking; + forwardpPacked.AppendValue(valuep); + forwardmPacked.AppendValue(valuem); + if (valuep == 0.0) { + forwardpPacked.AppendErrorValue(1.0); + } else { + forwardpPacked.AppendErrorValue(TMath::Sqrt(errorp)/fPacking); + } + if (valuem == 0.0) { + forwardmPacked.AppendErrorValue(1.0); + } else { + forwardmPacked.AppendErrorValue(TMath::Sqrt(errorm)/fPacking); + } + valuep = 0.0; + errorp = 0.0; + valuem = 0.0; + errorm = 0.0; } - value += fForward[i]; - error += fForwardErr[i]*fForwardErr[i]; + valuep += fForwardp[i]; + errorp += fForwardpErr[i]*fForwardpErr[i]; + valuem += fForwardm[i]; + errorm += fForwardmErr[i]*fForwardmErr[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; + valuep /= fPacking; + valuem /= fPacking; + backwardpPacked.AppendValue(valuep); + backwardmPacked.AppendValue(valuem); + if (valuep == 0.0) { + backwardpPacked.AppendErrorValue(1.0); + } else { + backwardpPacked.AppendErrorValue(TMath::Sqrt(errorp)/fPacking); + } + if (valuem == 0.0) { + backwardmPacked.AppendErrorValue(1.0); + } else { + backwardmPacked.AppendErrorValue(TMath::Sqrt(errorm)/fPacking); + } + valuep = 0.0; + errorp = 0.0; + valuem = 0.0; + errorm = 0.0; } - value += fBackward[i]; - error += fBackwardErr[i]*fBackwardErr[i]; + valuep += fBackwardp[i]; + errorp += fBackwardpErr[i]*fBackwardpErr[i]; + valuem += fBackwardm[i]; + errorm += fBackwardmErr[i]*fBackwardmErr[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(); + UInt_t noOfBins = forwardpPacked.GetValue()->size(); + if (forwardpPacked.GetValue()->size() != backwardpPacked.GetValue()->size()) { + if (forwardpPacked.GetValue()->size() > backwardpPacked.GetValue()->size()) + noOfBins = backwardpPacked.GetValue()->size(); } // form asymmetry including error propagation - Double_t asym; - Double_t f, b, ef, eb; + Double_t asym,error; + Double_t fp, bp, efp, ebp; + Double_t fm, bm, efm, ebm; // 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); + fp = forwardpPacked.GetValue()->at(i); + bp = backwardpPacked.GetValue()->at(i); + efp = forwardpPacked.GetError()->at(i); + ebp = backwardpPacked.GetError()->at(i); + fm = forwardmPacked.GetValue()->at(i); + bm = backwardmPacked.GetValue()->at(i); + efm = forwardmPacked.GetError()->at(i); + ebm = backwardmPacked.GetError()->at(i); // check that there are indeed bins - if (f+b != 0.0) - asym = (f-b) / (f+b); + if (fp+bp != 0.0) + asym = (fp-bp) / (fp+bp) - (fm-bm) / (fm+bm); 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); + if (fp+bp != 0.0) + errorp = 2.0/((fp+bp)*(fp+bp))*TMath::Sqrt(bp*bp*efp*efp+ebp*ebp*fp*fp); else - error = 1.0; + errorp = 1.0; + if (fm+bm != 0.0) + errorm = 2.0/((fm+bm)*(fm+bm))*TMath::Sqrt(bm*bm*efm*efm+ebm*ebm*fm*fm); + else + errorp = 1.0; + + error = TMath::Sqrt(errorp*errorp+errorm*errorm); fData.AppendErrorValue(error); } CalcNoOfFitBins(); // clean up - fForward.clear(); - fForwardErr.clear(); - fBackward.clear(); - fBackwardErr.clear(); + fForwardp.clear(); + fForwardpErr.clear(); + fBackwardp.clear(); + fBackwardpErr.clear(); + fForwardm.clear(); + fForwardmErr.clear(); + fBackwardm.clear(); + fBackwardmErr.clear(); return true; } @@ -1065,67 +1142,100 @@ Bool_t PRunAsymmetryBNMR::PrepareViewData(PRawRunData* runData, UInt_t histoNo[2 } // everything looks fine, hence fill packed forward and backward histo - PRunData forwardPacked; - PRunData backwardPacked; - Double_t value = 0.0; - Double_t error = 0.0; + PRunData forwardpPacked; + PRunData backwardpPacked; + PRunData forwardmPacked; + PRunData backwardmPacked; + Double_t valuep = 0.0; + Double_t errorp = 0.0; + Double_t valuem = 0.0; + Double_t errorm = 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; + valuep /= packing; + forwardpPacked.AppendValue(valuep); + valuem /= packing; + forwardmPacked.AppendValue(valuem); + if (valuep == 0.0) { + forwardpPacked.AppendErrorValue(1.0); + } else { + forwardpPacked.AppendErrorValue(TMath::Sqrt(errorp)/packing); + } + if (valuem == 0.0) { + forwardmPacked.AppendErrorValue(1.0); + } else { + forwardmPacked.AppendErrorValue(TMath::Sqrt(errorm)/packing); + } + valuep = 0.0; + errorp = 0.0; + valuem = 0.0; + errorm = 0.0; } - value += fForward[i]; - error += fForwardErr[i]*fForwardErr[i]; + valuep += fForwardp[i]; + errorp += fForwardpErr[i]*fForwardpErr[i]; + valuem += fForwardm[i]; + errorm += fForwardmErr[i]*fForwardmErr[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; + valuep /= packing; + valuem /= packing; + backwardpPacked.AppendValue(valuep); + backwardmPacked.AppendValue(valuem); + if (valuep == 0.0) { + backwardpPacked.AppendErrorValue(1.0); + } else { + backwardpPacked.AppendErrorValue(TMath::Sqrt(errorp)/packing); + } + if (valuem == 0.0) { + backwardmPacked.AppendErrorValue(1.0); + } else { + backwardmPacked.AppendErrorValue(TMath::Sqrt(errorm)/packing); + } + valuep = 0.0; + errorp = 0.0; + valuem = 0.0; + errorm = 0.0; } - value += fBackward[i]; - error += fBackwardErr[i]*fBackwardErr[i]; + valuep += fBackwardp[i]; + errorp += fBackwardpErr[i]*fBackwardpErr[i]; + valuem += fBackwardm[i]; + errorm += fBackwardmErr[i]*fBackwardmErr[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(); + UInt_t noOfBins = forwardpPacked.GetValue()->size(); + if (forwardpPacked.GetValue()->size() != backwardpPacked.GetValue()->size()) { + if (forwardpPacked.GetValue()->size() > backwardpPacked.GetValue()->size()) + noOfBins = backwardpPacked.GetValue()->size(); } // form asymmetry including error propagation Double_t asym; - Double_t f, b, ef, eb, alpha = 1.0, beta = 1.0; + Double_t fp, bp, efp, ebp, alpha = 1.0, beta = 1.0; + Double_t fm, bm, efm, ebm; // 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)); @@ -1153,33 +1263,46 @@ Bool_t PRunAsymmetryBNMR::PrepareViewData(PRawRunData* runData, UInt_t histoNo[2 break; } - for (UInt_t i=0; isize(); i++) { + 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); + fp = forwardpPacked.GetValue()->at(i); + bp = backwardpPacked.GetValue()->at(i); + efp = forwardpPacked.GetError()->at(i); + ebp = backwardpPacked.GetError()->at(i); + fm = forwardmPacked.GetValue()->at(i); + bm = backwardmPacked.GetValue()->at(i); + efm = forwardmPacked.GetError()->at(i); + ebm = backwardmPacked.GetError()->at(i); // check that there are indeed bins - if (f+b != 0.0) - asym = (alpha*f-b) / (alpha*beta*f+b); + if (fp+bp != 0.0) + asym = (alpha*fp-bp) / (alpha*beta*fp+bp) - (alpha*fm-bm) / (alpha*beta*fm+bm); 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); + if (fp+bp != 0.0) + errorp = 2.0/((fp+bp)*(fp+bp))*TMath::Sqrt(bp*bp*efp*efp+ebp*ebp*fp*fp); else - error = 1.0; + errorp = 1.0; + if (fm+bm != 0.0) + errorm = 2.0/((fm+bm)*(fm+bm))*TMath::Sqrt(bm*bm*efm*efm+ebm*ebm*fm*fm); + else + errorm = 1.0; + error = TMath::Sqrt(errorp*errorp+errorm*errorm); fData.AppendErrorValue(error); } CalcNoOfFitBins(); // clean up - fForward.clear(); - fForwardErr.clear(); - fBackward.clear(); - fBackwardErr.clear(); + fForwardp.clear(); + fForwardpErr.clear(); + fBackwardp.clear(); + fBackwardpErr.clear(); + fForwardm.clear(); + fForwardmErr.clear(); + fBackwardm.clear(); + fBackwardmErr.clear(); // fill theory vector for kView // calculate functions @@ -1212,311 +1335,6 @@ Bool_t PRunAsymmetryBNMR::PrepareViewData(PRawRunData* runData, UInt_t histoNo[2 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) diff --git a/src/include/PRunAsymmetryBNMR.h b/src/include/PRunAsymmetryBNMR.h index ec2f80a7..3e3201d9 100644 --- a/src/include/PRunAsymmetryBNMR.h +++ b/src/include/PRunAsymmetryBNMR.h @@ -68,10 +68,14 @@ class PRunAsymmetryBNMR : public PRunBase 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 + PDoubleVector fForwardp; ///< pos hel forward histo data + PDoubleVector fForwardpErr; ///< pos hel forward histo errors + PDoubleVector fBackwardp; ///< pos hel backward histo data + PDoubleVector fBackwardpErr; ///< pos hel backward histo errors + PDoubleVector fForwardm; ///< neg hel forward histo data + PDoubleVector fForwardmErr; ///< neg hel forward histo errors + PDoubleVector fBackwardm; ///< neg hel backward histo data + PDoubleVector fBackwardmErr; ///< neg hel backward histo errors Int_t fGoodBins[4]; ///< keep first/last good bins. 0=fgb, 1=lgb (forward); 2=fgb, 3=lgb (backward) From 52a5bc25b14e79833ac01cbf880eee676ae42236 Mon Sep 17 00:00:00 2001 From: Zaher Salman Date: Sun, 19 Aug 2018 23:40:55 +0200 Subject: [PATCH 07/20] Start substract asymmetry implementation. --- src/include/PRunAsymmetryBNMR.h | 1 - 1 file changed, 1 deletion(-) diff --git a/src/include/PRunAsymmetryBNMR.h b/src/include/PRunAsymmetryBNMR.h index 3e3201d9..a1fc5377 100644 --- a/src/include/PRunAsymmetryBNMR.h +++ b/src/include/PRunAsymmetryBNMR.h @@ -61,7 +61,6 @@ class PRunAsymmetryBNMR : public PRunBase 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$. From e9ec126d9e3d26e8af69ba96e51655f36b497c92 Mon Sep 17 00:00:00 2001 From: Zaher Salman Date: Mon, 20 Aug 2018 17:51:39 +0200 Subject: [PATCH 08/20] Beta-NMR helicity subtraction seems to work in the new class. Still needs debugging and tests. Some of the changes from Jonas were reverted since SExp does not work otherwise. --- src/classes/PRunAsymmetryBNMR.cpp | 70 ++++++++++++++-------------- src/external/libBNMR/ExpRlx-MUD.msr | 42 ++++++----------- src/external/libBNMR/SExpRlx-MUD.msr | 61 ++++++++++++++++++++++++ src/external/libBNMR/TBNMR.cpp | 42 +++++++---------- src/external/libBNMR/TBNMR.h | 9 ++-- 5 files changed, 131 insertions(+), 93 deletions(-) create mode 100644 src/external/libBNMR/SExpRlx-MUD.msr diff --git a/src/classes/PRunAsymmetryBNMR.cpp b/src/classes/PRunAsymmetryBNMR.cpp index 07bb2c3c..f9aa4cfe 100644 --- a/src/classes/PRunAsymmetryBNMR.cpp +++ b/src/classes/PRunAsymmetryBNMR.cpp @@ -539,11 +539,11 @@ Bool_t PRunAsymmetryBNMR::PrepareData() return false; } - // keep the time resolution in (us) + // keep the time resolution in (ms) // possibility to rescale for betaNMR fTimeResolution = runData->GetTimeResolution()/1.0e3; cout.precision(10); - cout << endl << ">> PRunAsymmetryBNMR::PrepareData(): time resolution=" << fixed << runData->GetTimeResolution() << "(ns)" << endl; + cout << endl << ">> PRunAsymmetryBNMR::PrepareData(): time resolution=" << fixed << runData->GetTimeResolution() << "(ms)" << endl; // get all the proper t0's and addt0's for the current RUN block if (!GetProperT0(runData, globalBlock, forwardHistoNo, backwardHistoNo)) { @@ -602,6 +602,8 @@ Bool_t PRunAsymmetryBNMR::PrepareData() // set forward/backward histo data of the first group fForwardp.resize(forward[0].size()); fBackwardp.resize(backward[0].size()); + fForwardm.resize(forward[0].size()); + fBackwardm.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 @@ -706,22 +688,40 @@ Bool_t PRunAsymmetryBNMR::PrepareData() Bool_t PRunAsymmetryBNMR::SubtractFixBkg() { Double_t dval; - for (UInt_t i=0; iGetBkgFix(0); + fForwardpErr.push_back(dval); + fForwardp[i] -= fRunInfo->GetBkgFix(0); // keep the error, and make sure that the bin is NOT empty - if (fBackward[i] != 0.0) - dval = TMath::Sqrt(fBackward[i]); + if (fForwardm[i] != 0.0) + dval = TMath::Sqrt(fForwardm[i]); else dval = 1.0; - fBackwardErr.push_back(dval); - fBackward[i] -= fRunInfo->GetBkgFix(1); + fForwardmErr.push_back(dval); + fForwardm[i] -= fRunInfo->GetBkgFix(1); + + // keep the error, and make sure that the bin is NOT empty + if (fBackwardp[i] != 0.0) + dval = TMath::Sqrt(fBackwardp[i]); + else + dval = 1.0; + fBackwardpErr.push_back(dval); + fBackwardp[i] -= fRunInfo->GetBkgFix(2); + + // keep the error, and make sure that the bin is NOT empty + if (fBackwardm[i] != 0.0) + dval = TMath::Sqrt(fBackwardm[i]); + else + dval = 1.0; + fBackwardmErr.push_back(dval); + fBackwardm[i] -= fRunInfo->GetBkgFix(3); } return true; @@ -807,8 +807,8 @@ Bool_t PRunAsymmetryBNMR::SubtractEstimatedBkg() // calculate background Double_t bkgp[2] = {0.0, 0.0}; Double_t errBkgp[2] = {0.0, 0.0}; - Double_t bkgn[2] = {0.0, 0.0}; - Double_t errBkgn[2] = {0.0, 0.0}; + Double_t bkgm[2] = {0.0, 0.0}; + Double_t errBkgm[2] = {0.0, 0.0}; // forward for (UInt_t i=start[0]; i<=end[0]; i++) { @@ -898,7 +898,7 @@ Bool_t PRunAsymmetryBNMR::PrepareFitData() PRunData forwardpPacked; PRunData backwardpPacked; PRunData forwardmPacked; - PRunData backwarmpPacked; + PRunData backwardmPacked; Double_t valuep = 0.0; Double_t errorp = 0.0; Double_t valuem = 0.0; @@ -1150,6 +1150,8 @@ Bool_t PRunAsymmetryBNMR::PrepareViewData(PRawRunData* runData, UInt_t histoNo[2 Double_t errorp = 0.0; Double_t valuem = 0.0; Double_t errorm = 0.0; + Double_t value = 0.0; + Double_t error = 0.0; // forward for (Int_t i=start[0]; i &par) const { assert(par.size()==2); // make sure the number of parameters handed to the function is correct - // par[0] time of beam off - // par[1] is the relaxation rate + // par[0] time of beam on (pulse length) in seconds + // par[1] is the relaxation rate in 1/s double tau_p; - double y; - tau_p = (tau_Li/(1.+par[1]*tau_Li)); - + + // x should be in seconds, otherwise it should be rescaled here if ( x <= par[0] && x >= 0) { - y=(tau_p/tau_Li)*(1-exp(-x/tau_p))/(1-exp(-x/tau_Li)); + return (tau_p/tau_Li)*(1-exp(-x/tau_p))/(1-exp(-x/tau_Li)); } else if ( x > par[0] ){ - y=(tau_p/tau_Li)*(1-exp(-par[0]/tau_p))/(1-exp(-par[0]/tau_Li))*exp(-par[1]*(x-par[0])); - } else { - y = 0; - } - - return y; + return (tau_p/tau_Li)*(1-exp(-par[0]/tau_p))/(1-exp(-par[0]/tau_Li))*exp(-par[1]*(x-par[0])); + } + return 0; } - -//initialize Integrators -TF1 SExpRlx::sexp1=TF1("sexp", "exp(-([0]-x)/[3])*exp(-pow(([1]*([0]-x)),[2]))", 0.0, 20000.0); -TF1 SExpRlx::sexp2=TF1("sexp", "exp(-([3]-x)/[4])*exp(-pow(([1]*([0]-x)),[2]))", 0.0, 20000.0); - double SExpRlx::operator()(double x, const vector &par) const { assert(par.size()==3); // make sure the number of parameters handed to the function is correct - // par[0] time of beam off - // par[1] is the relaxation rate + // par[0] beam of beam on (pulse length) in seconds + // par[1] is the relaxation rate in 1/s // par[2] is the exponent + // x should be in seconds, otherwise it should be rescaled here if ( x >= 0 && x <= par[0] ) { + TF1 sexp1("sexp1", "exp(-([0]-x)/[3])*exp(-pow(([1]*([0]-x)),[2]))", 0.0, 20.0); sexp1.SetParameters(x, par[1], par[2],tau_Li); + sexp1.SetNpx(1000); return sexp1.Integral(0.0,x)/(1-exp(-x/tau_Li))/tau_Li; } else if ( x > par[0] ) { + TF1 sexp2("sexp2", "exp(-([3]-x)/[4])*exp(-pow(([1]*([0]-x)),[2]))", 0.0, 20.0); sexp2.SetParameters(x, par[1], par[2], par[0],tau_Li); + sexp2.SetNpx(1000); return sexp2.Integral(0.0,par[0])/(1-exp(-par[0]/tau_Li))/tau_Li; - } + } return 0; } diff --git a/src/external/libBNMR/TBNMR.h b/src/external/libBNMR/TBNMR.h index 24d95172..042feb03 100644 --- a/src/external/libBNMR/TBNMR.h +++ b/src/external/libBNMR/TBNMR.h @@ -41,6 +41,10 @@ #ifndef LIBBNMRH #define LIBBNMRH +#define tau_Li 1.210 // In seconds +#define PI 3.14159265358979323846 +#define TWOPI 6.28318530717958647692 + using namespace std; @@ -66,7 +70,7 @@ class SExpRlx : public PUserFcnBase { public: // default constructor and destructor - SExpRlx(){sexp1.SetNpx(1000); sexp2.SetNpx(1000);} + SExpRlx(){} ~SExpRlx(){} Bool_t NeedGlobalPart() const { return false; } @@ -75,9 +79,6 @@ public: // function operator double operator()(double, const vector&) const; -private: - static TF1 sexp1; - static TF1 sexp2; // definition of the class for the ROOT-dictionary ClassDef(SExpRlx,1) From fa7e74e56096ad354f4c79095e50fee3d922b6ee Mon Sep 17 00:00:00 2001 From: Zaher Salman Date: Mon, 20 Aug 2018 20:04:43 +0200 Subject: [PATCH 09/20] Remove beta-NMR ascci based templates to avoid confusion. --- src/external/libBNMR/45377.dat | 96 ---------------------------- src/external/libBNMR/ExpRlx.msr | 45 ------------- src/external/libBNMR/SExpRlx-MUD.msr | 14 ++-- src/external/libBNMR/SExpRlx.msr | 46 ------------- 4 files changed, 7 insertions(+), 194 deletions(-) delete mode 100644 src/external/libBNMR/45377.dat delete mode 100644 src/external/libBNMR/ExpRlx.msr delete mode 100644 src/external/libBNMR/SExpRlx.msr diff --git a/src/external/libBNMR/45377.dat b/src/external/libBNMR/45377.dat deleted file mode 100644 index 36f5f791..00000000 --- a/src/external/libBNMR/45377.dat +++ /dev/null @@ -1,96 +0,0 @@ -DATA - 592.75 0.0466642055 0.00075601151 - 693.25 0.042699595 0.0007136273 - 793.75 0.03943051 0.000679583922 - 894.25 0.0366291886 0.000652378784 - 994.75 0.0355316716 0.000664640277 - 1095.25 0.033304272 0.000610910492 - 1195.75 0.0322397741 0.000595454518 - 1296.25 0.0308209315 0.000581648459 - 1396.75 0.0301227967 0.000569861788 - 1497.25 0.0289808195 0.000559611024 - 1597.75 0.0284020831 0.000550685264 - 1698.25 0.0272996592 0.000542785334 - 1798.75 0.0259033059 0.000536112979 - 1899.25 0.0257421462 0.000529928146 - 1999.75 0.0252476607 0.000524276201 - 2100.25 0.0252439822 0.00049508141 - 2200.75 0.0244494569 0.000514789454 - 2301.25 0.0248826473 0.000510792766 - 2401.75 0.0252364689 0.000507093213 - 2502.25 0.024077753 0.000503959631 - 2602.75 0.0239705744 0.000501170092 - 2703.25 0.023648613 0.000498351823 - 2803.75 0.0236132062 0.000496113677 - 2904.25 0.0232728359 0.000494038794 - 3004.75 0.0229806315 0.00049216698 - 3105.25 0.0226408705 0.000490122174 - 3205.75 0.022949639 0.000488517677 - 3306.25 0.0236152274 0.000487206227 - 3406.75 0.022514855 0.000485776619 - 3507.25 0.0232919509 0.000484617162 - 3607.75 0.0227762445 0.000483435077 - 3708.25 0.0225740091 0.000482352785 - 3808.75 0.0217810425 0.000481423643 - 3909.25 0.0210459719 0.000480488676 - 4009.75 0.0200829485 0.0004835672 - 4110.25 0.0144273479 0.000480165881 - 4210.75 0.0104726937 0.000526210282 - 4311.25 0.00977670745 0.000548717033 - 4411.75 0.00734786823 0.000572069626 - 4512.25 0.0071749696 0.000596458669 - 4612.75 0.0061713664 0.000621795033 - 4713.25 0.00457656715 0.000648552578 - 4813.75 0.00253947053 0.000676064828 - 4914.25 0.00298436429 0.000705158693 - 5014.75 0.00520775648 0.000735097152 - 5115.25 0.00372267868 0.000766213043 - 5215.75 0.00284020846 0.000798853062 - 5316.25 0.00158717036 0.000832862333 - 5416.75 0.00200819085 0.000868342561 - 5517.25 0.0044706107 0.000905327348 - 5617.75 0.00135523051 0.000943421323 - 5718.25 0.000261760966 0.000983714257 - 5818.75 0.000210009562 0.00102515713 - 5919.25 -0.000157024642 0.00106890614 - 6019.75 0.000611697915 0.00111378026 - 6120.25 0.000353532075 0.00111004546 - 6220.75 0.00120570511 0.00121565301 - 6321.25 -0.00126499964 0.00126730997 - 6421.75 0.0029594966 0.00132005144 - 6522.25 -0.00180278144 0.00137777581 - 6622.75 0.00154819637 0.00143478041 - 6723.25 0.00137296473 0.00149757271 - 6823.75 -0.000271202669 0.00155948172 - 6924.25 -0.00306837271 0.00162568622 - 7024.75 0.00346602562 0.00169188976 - 7125.25 -8.40576554E-05 0.00176525436 - 7225.75 0.00141055893 0.00183909095 - 7326.25 -0.00170315946 0.00191744238 - 7426.75 0.00399357243 0.00199844555 - 7527.25 0.00378603394 0.00208493839 - 7627.75 -0.0015303158 0.00216992033 - 7728.25 -0.000931658073 0.00226567509 - 7828.75 0.0011671311 0.00235924947 - 7929.25 0.00147725609 0.00245769641 - 8029.75 -0.00100983925 0.00256942322 - 8130.25 -0.00329542124 0.00255286858 - 8230.75 0.000636890228 0.00279597795 - 8331.25 -0.000814018633 0.00291290155 - 8431.75 0.000101912547 0.00303673191 - 8532.25 0.00311038784 0.0031603991 - 8632.75 0.000473049572 0.00329964113 - 8733.25 8.36732E-05 0.00343974064 - 8833.75 0.000275585165 0.00358184548 - 8934.25 0.0059381551 0.00373021006 - 9034.75 -0.0014588358 0.0038926573 - 9135.25 -0.00158156038 0.00405870998 - 9235.75 -0.00125287442 0.00422594705 - 9336.25 -0.00212622165 0.00439704849 - 9436.75 -0.00477056006 0.0045911195 - 9537.25 -0.00229735693 0.00479087708 - 9637.75 0.00575469607 0.00498561109 - 9738.25 -7.52893075E-05 0.00519483122 - 9838.75 -0.00368605583 0.0054198485 - 9939.25 -0.00862419517 0.00564201231 - 10039.75 -0.0323936833 0.0181186246 diff --git a/src/external/libBNMR/ExpRlx.msr b/src/external/libBNMR/ExpRlx.msr deleted file mode 100644 index 2b3f52ac..00000000 --- a/src/external/libBNMR/ExpRlx.msr +++ /dev/null @@ -1,45 +0,0 @@ -Title -############################################################### -FITPARAMETER -# Nr. Name Value Step Pos_Error Boundaries - 1 Asy1 0.0812706 0.00149848 none 0 none - 2 T 4000 0 none 0 none - 3 Lam1 0.00239816 6.05947e-05 none 0 none - -############################################################### -THEORY -asymmetry 1 -userFcn /usr/local/lib/libBNMR.so ExpRlx 2 3 - -############################################################### -RUN 45377 MUE4 PSI ASCII (name beamline institute data-file-format) -fittype 8 (non muSR fit) -map 0 0 0 0 0 0 0 0 0 0 -xy-data 1 2 -fit 0.00 8000.00 -packing 1 - -############################################################### -COMMANDS -STRATEGY 1 -MINIMIZE -#MINOS -SAVE -END RETURN - -############################################################### -PLOT 8 (non muSR plot) -runs 1 -range 0.00 8000.00 - -############################################################### -FOURIER -units MHz # units either 'Gauss', 'MHz', or 'Mc/s' -fourier_power 12 -apodization STRONG # NONE, WEAK, MEDIUM, STRONG -plot POWER # REAL, IMAG, REAL_AND_IMAG, POWER, PHASE -phase 8.50 -#range FRQMIN FRQMAX -############################################################### -STATISTIC --- 2010-09-02 17:44:20 - chisq = 186.322493286053, NDF = 72, chisq/NDF = 2.5878124067507362 diff --git a/src/external/libBNMR/SExpRlx-MUD.msr b/src/external/libBNMR/SExpRlx-MUD.msr index fd76633d..f689af83 100644 --- a/src/external/libBNMR/SExpRlx-MUD.msr +++ b/src/external/libBNMR/SExpRlx-MUD.msr @@ -5,10 +5,10 @@ FITPARAMETER ############################################################### # No Name Value Err Min Max 1 Alpha 1 0 none - 2 Asy 0.108 0.011 none 0 0.2 + 2 Asy 0.1152 0.0099 none 0 0.2 3 T 1 0 none - 4 Rlx 3.63 0.89 none 0 15000 - 5 Beta 0.452 0.039 none 0.3 2 + 4 Rlx 4.19 0.91 none 0 15000 + 5 Beta 0.432 0.032 none 0.3 2 ############################################################### THEORY @@ -27,12 +27,12 @@ fittype 5 (beta-NMR fit) alpha 1 forward 3 5 backward 4 6 -data 11 800 11 800 +data 11 1000 11 1000 #backgr.fix 0 background 1 9 1 9 # estimated bkg: 8.0000 / 9.6250 t0 10.0 10.0 10.0 10.0 map 0 0 0 0 0 0 0 0 0 0 -fit 0.5 8 +fit 0.2 10 packing 5 ############################################################### @@ -57,5 +57,5 @@ plot POWER # REAL, IMAG, REAL_AND_IMAG, POWER, PHASE, PHASE_OPT_RE phase 8 #range FRQMIN FRQMAX ############################################################### -STATISTIC --- 2018-08-20 17:48:41 - chisq = 177.2, NDF = 144, chisq/NDF = 1.230711 +STATISTIC --- 2018-08-20 19:54:17 + chisq = 181.6, NDF = 150, chisq/NDF = 1.210991 diff --git a/src/external/libBNMR/SExpRlx.msr b/src/external/libBNMR/SExpRlx.msr deleted file mode 100644 index 636583f8..00000000 --- a/src/external/libBNMR/SExpRlx.msr +++ /dev/null @@ -1,46 +0,0 @@ -Title -############################################################### -FITPARAMETER -# Nr. Name Value Step Pos_Error Boundaries - 1 Asy1 0.262789 0.0465582 none 0 none - 2 T 4000 0 none 0 none - 3 Lam1 0.0233569 0.0094885 none 0 none - 4 Bet1 0.399259 0.0353915 none 0 none - -############################################################### -THEORY -asymmetry 1 -userFcn /usr/local/lib/libBNMR.so SExpRlx 2 3 4 - -############################################################### -RUN 45377 MUE4 PSI ASCII (name beamline institute data-file-format) -fittype 8 (non muSR fit) -map 0 0 0 0 0 0 0 0 0 0 -xy-data 1 2 -fit 0.00 8000.00 -packing 1 - -############################################################### -COMMANDS -STRATEGY 1 -MINIMIZE -#MINOS -SAVE -END RETURN - -############################################################### -PLOT 8 (non muSR plot) -runs 1 -range 0.00 8000.00 - -############################################################### -FOURIER -units MHz # units either 'Gauss', 'MHz', or 'Mc/s' -fourier_power 12 -apodization STRONG # NONE, WEAK, MEDIUM, STRONG -plot POWER # REAL, IMAG, REAL_AND_IMAG, POWER, PHASE -phase 8.50 -#range FRQMIN FRQMAX -############################################################### -STATISTIC --- 2010-09-02 17:51:10 - chisq = 93.71265965991455, NDF = 71, chisq/NDF = 1.319896614928374 From af4040228b3cb4cdf2b0d4a0bc8865de8ed4aba9 Mon Sep 17 00:00:00 2001 From: Zaher Salman Date: Mon, 20 Aug 2018 22:04:50 +0200 Subject: [PATCH 10/20] Template files for two helicities and subtracted helicities updated. --- src/external/libBNMR/ExpRlx-hel.msr | 72 ++++++++++++++++++ .../libBNMR/{ExpRlx-MUD.msr => ExpRlx.msr} | 21 ++---- src/external/libBNMR/SExpRlx-hel.msr | 73 +++++++++++++++++++ .../libBNMR/{SExpRlx-MUD.msr => SExpRlx.msr} | 23 ++---- 4 files changed, 160 insertions(+), 29 deletions(-) create mode 100644 src/external/libBNMR/ExpRlx-hel.msr rename src/external/libBNMR/{ExpRlx-MUD.msr => ExpRlx.msr} (77%) create mode 100644 src/external/libBNMR/SExpRlx-hel.msr rename src/external/libBNMR/{SExpRlx-MUD.msr => SExpRlx.msr} (74%) diff --git a/src/external/libBNMR/ExpRlx-hel.msr b/src/external/libBNMR/ExpRlx-hel.msr new file mode 100644 index 00000000..048996bc --- /dev/null +++ b/src/external/libBNMR/ExpRlx-hel.msr @@ -0,0 +1,72 @@ +LaAlO3, SLR, 30G, 150K, Bias=23 kV, Ix=0A, Iy=+3A, aligned, beam on=1s +############################################################### +FITPARAMETER +############################################################### +# No Name Value Err Min Max + 1 Alpha 1.11592 0.00035 none + 2 Asy 0.05869 0.00091 none 0 0.2 + 3 T 1 0 none + 4 Rlx 1.142 0.023 none 0 15000 + 5 One 1 0 none + 6 FlHel -1.047 0.022 none -2 0 + +############################################################### +THEORY +############################################################### +asymmetry fun1 +userFcn .libs/libBNMR.so ExpRlx 3 4 + +############################################################### +FUNCTIONS +############################################################### +fun1 = 0.5 * map1 * map2 + +############################################################### +RUN 045674 BNMR TRIUMF MUD (name beamline institute data-file-format) +fittype 2 (asymmetry fit) +alpha 1 +forward 3 +backward 4 +data 11 1000 11 1000 +background 1 9 1 9 # estimated bkg: 8.0000 / 9.6250 +t0 10.0 10.0 10.0 10.0 +map 2 5 0 0 0 0 0 0 0 0 +fit 0.2 9 +packing 5 + +RUN 045674 BNMR TRIUMF MUD (name beamline institute data-file-format) +fittype 2 (asymmetry fit) +alpha 1 +forward 5 +backward 6 +data 11 1000 11 1000 +background 1 9 1 9 # estimated bkg: 11.6250 / 15.6250 +t0 10.0 10.0 10.0 10.0 +map 2 6 0 0 0 0 0 0 0 0 +fit 0.2 9 +packing 5 + +############################################################### +COMMANDS +MINIMIZE +HESSE +SAVE + +############################################################### +PLOT 2 (asymmetry plot) +runs 1 2 +use_fit_ranges +view_packing 10 + + +############################################################### +FOURIER +units MHz # units either 'Gauss', 'Tesla', 'MHz', or 'Mc/s' +fourier_power 12 +apodization STRONG # NONE, WEAK, MEDIUM, STRONG +plot POWER # REAL, IMAG, REAL_AND_IMAG, POWER, PHASE, PHASE_OPT_REAL +phase 8 +#range FRQMIN FRQMAX +############################################################### +STATISTIC --- 2018-08-20 22:01:15 + chisq = 566.6, NDF = 346, chisq/NDF = 1.637584 diff --git a/src/external/libBNMR/ExpRlx-MUD.msr b/src/external/libBNMR/ExpRlx.msr similarity index 77% rename from src/external/libBNMR/ExpRlx-MUD.msr rename to src/external/libBNMR/ExpRlx.msr index 909af362..eb356304 100644 --- a/src/external/libBNMR/ExpRlx-MUD.msr +++ b/src/external/libBNMR/ExpRlx.msr @@ -1,13 +1,12 @@ - -# Run Numbers: 1111 +LaAlO3, SLR, 30G, 150K, Bias=23 kV, Ix=0A, Iy=+3A, aligned, beam on=1s ############################################################### FITPARAMETER ############################################################### # No Name Value Err Min Max 1 Alpha 1 0 none - 2 Asy 0.05615 0.00073 none + 2 Asy 0.05992 0.00063 none 0 0.2 3 T 1 0 none - 4 Rlx 1.046 0.024 none + 4 Rlx 1.143 0.023 none 0 100 ############################################################### THEORY @@ -15,23 +14,17 @@ THEORY asymmetry 2 userFcn .libs/libBNMR.so ExpRlx 3 4 -############################################################### -#FUNCTIONS -############################################################### -#fun1 = 0.5 * map1 * map2 - ############################################################### RUN 045674 BNMR TRIUMF MUD (name beamline institute data-file-format) fittype 5 (beta-NMR fit) alpha 1 forward 3 5 backward 4 6 -data 11 800 11 800 -#backgr.fix 0 +data 11 1000 11 1000 background 1 9 1 9 # estimated bkg: 8.0000 / 9.6250 t0 10.0 10.0 10.0 10.0 map 0 0 0 0 0 0 0 0 0 0 -fit 0.5 8 +fit 0.2 9 packing 5 ############################################################### @@ -56,5 +49,5 @@ plot POWER # REAL, IMAG, REAL_AND_IMAG, POWER, PHASE, PHASE_OPT_RE phase 8 #range FRQMIN FRQMAX ############################################################### -STATISTIC --- 2018-08-20 17:31:40 - chisq = 309.1, NDF = 145, chisq/NDF = 2.131395 +STATISTIC --- 2018-08-20 22:02:02 + chisq = 419.1, NDF = 173, chisq/NDF = 2.422357 diff --git a/src/external/libBNMR/SExpRlx-hel.msr b/src/external/libBNMR/SExpRlx-hel.msr new file mode 100644 index 00000000..7165dd2c --- /dev/null +++ b/src/external/libBNMR/SExpRlx-hel.msr @@ -0,0 +1,73 @@ +LaAlO3, SLR, 30G, 150K, Bias=23 kV, Ix=0A, Iy=+3A, aligned, beam on=1s +############################################################### +FITPARAMETER +############################################################### +# No Name Value Err Min Max + 1 Alpha 1.11589 0.00037 none + 2 Asy 0.1121 0.0094 none 0 0.2 + 3 T 1 0 none + 4 Rlx 4.12 0.87 none 0 15000 + 5 Beta 0.434 0.031 none 0.3 2 + 6 One 1 0 none + 7 FlHel -1.048 0.023 none -2 0 + +############################################################### +THEORY +############################################################### +asymmetry fun1 +userFcn .libs/libBNMR.so SExpRlx 3 4 5 + +############################################################### +FUNCTIONS +############################################################### +fun1 = 0.5 * map1 * map2 + +############################################################### +RUN 045674 BNMR TRIUMF MUD (name beamline institute data-file-format) +fittype 2 (asymmetry fit) +alpha 1 +forward 3 +backward 4 +data 11 1000 11 1000 +background 1 9 1 9 # estimated bkg: 8.0000 / 9.6250 +t0 10.0 10.0 10.0 10.0 +map 2 6 0 0 0 0 0 0 0 0 +fit 0.2 9 +packing 5 + +RUN 045674 BNMR TRIUMF MUD (name beamline institute data-file-format) +fittype 2 (asymmetry fit) +alpha 1 +forward 5 +backward 6 +data 11 1000 11 1000 +background 1 9 1 9 # estimated bkg: 11.6250 / 15.6250 +t0 10.0 10.0 10.0 10.0 +map 2 7 0 0 0 0 0 0 0 0 +fit 0.2 9 +packing 5 + +############################################################### +COMMANDS +MINIMIZE +HESSE +SAVE + +############################################################### +PLOT 2 (asymmetry plot) +runs 1 2 +use_fit_ranges +view_packing 10 + + +############################################################### +FOURIER +units MHz # units either 'Gauss', 'Tesla', 'MHz', or 'Mc/s' +fourier_power 12 +apodization STRONG # NONE, WEAK, MEDIUM, STRONG +plot POWER # REAL, IMAG, REAL_AND_IMAG, POWER, PHASE, PHASE_OPT_REAL +phase 8 +#range FRQMIN FRQMAX +############################################################### +STATISTIC --- 2018-08-20 21:59:04 + chisq = 358.6, NDF = 345, chisq/NDF = 1.039300 diff --git a/src/external/libBNMR/SExpRlx-MUD.msr b/src/external/libBNMR/SExpRlx.msr similarity index 74% rename from src/external/libBNMR/SExpRlx-MUD.msr rename to src/external/libBNMR/SExpRlx.msr index f689af83..a530e625 100644 --- a/src/external/libBNMR/SExpRlx-MUD.msr +++ b/src/external/libBNMR/SExpRlx.msr @@ -1,14 +1,13 @@ - -# Run Numbers: 1111 +LaAlO3, SLR, 30G, 150K, Bias=23 kV, Ix=0A, Iy=+3A, aligned, beam on=1s ############################################################### FITPARAMETER ############################################################### # No Name Value Err Min Max 1 Alpha 1 0 none - 2 Asy 0.1152 0.0099 none 0 0.2 + 2 Asy 0.1148 0.0096 none 0 0.2 3 T 1 0 none - 4 Rlx 4.19 0.91 none 0 15000 - 5 Beta 0.432 0.032 none 0.3 2 + 4 Rlx 4.16 0.88 none 0 100 + 5 Beta 0.433 0.031 none 0.3 2 ############################################################### THEORY @@ -16,23 +15,17 @@ THEORY asymmetry 2 userFcn .libs/libBNMR.so SExpRlx 3 4 5 -############################################################### -#FUNCTIONS -############################################################### -#fun1 = 0.5 * map1 * map2 - ############################################################### RUN 045674 BNMR TRIUMF MUD (name beamline institute data-file-format) fittype 5 (beta-NMR fit) alpha 1 forward 3 5 backward 4 6 -data 11 1000 11 1000 -#backgr.fix 0 +data 11 1000 11 1000 background 1 9 1 9 # estimated bkg: 8.0000 / 9.6250 t0 10.0 10.0 10.0 10.0 map 0 0 0 0 0 0 0 0 0 0 -fit 0.2 10 +fit 0.2 9 packing 5 ############################################################### @@ -57,5 +50,5 @@ plot POWER # REAL, IMAG, REAL_AND_IMAG, POWER, PHASE, PHASE_OPT_RE phase 8 #range FRQMIN FRQMAX ############################################################### -STATISTIC --- 2018-08-20 19:54:17 - chisq = 181.6, NDF = 150, chisq/NDF = 1.210991 +STATISTIC --- 2018-08-20 21:59:14 + chisq = 210.2, NDF = 172, chisq/NDF = 1.222196 From b8dc67fa1adc97c18bc7c19f2a488d5a8df21e20 Mon Sep 17 00:00:00 2001 From: Zaher Salman Date: Tue, 21 Aug 2018 10:32:57 +0200 Subject: [PATCH 11/20] Fix x-axis label bug (reported by Gerland Morris). --- src/classes/PMusrCanvas.cpp | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/src/classes/PMusrCanvas.cpp b/src/classes/PMusrCanvas.cpp index 60d067a2..a5a3ec2d 100644 --- a/src/classes/PMusrCanvas.cpp +++ b/src/classes/PMusrCanvas.cpp @@ -4690,7 +4690,8 @@ void PMusrCanvas::PlotData(Bool_t unzoom) // set x-axis label PMsrRunList runs = *fMsrHandler->GetMsrRunList(); TString setup = fRunList->GetSetup(*runs[0].GetRunName()); - if (strcmp(setup, "TRIUMF/BNQR") || strcmp(setup, "TRIUMF/BNMR")) { + // For BNMR/BNQR runs use seconds + if (fPlotType == MSR_PLOT_BNMR) { fHistoFrame->GetXaxis()->SetTitle("time (s)"); } else { fHistoFrame->GetXaxis()->SetTitle("time (#mus)"); @@ -4984,7 +4985,8 @@ void PMusrCanvas::PlotDifference(Bool_t unzoom) // set x-axis label PMsrRunList runs = *fMsrHandler->GetMsrRunList(); TString setup = fRunList->GetSetup(*runs[0].GetRunName()); - if (strcmp(setup, "TRIUMF/BNQR") || strcmp(setup, "TRIUMF/BNMR")) { + // For BNMR/BNQR runs use seconds + if (fPlotType == MSR_PLOT_BNMR) { fHistoFrame->GetXaxis()->SetTitle("time (s)"); } else { fHistoFrame->GetXaxis()->SetTitle("time (#mus)"); @@ -6105,7 +6107,8 @@ void PMusrCanvas::PlotAverage(Bool_t unzoom) if (fCurrentPlotView == PV_DATA) { PMsrRunList runs = *fMsrHandler->GetMsrRunList(); TString setup = fRunList->GetSetup(*runs[0].GetRunName()); - if (strcmp(setup, "TRIUMF/BNQR") || strcmp(setup, "TRIUMF/BNMR")) { + // For BNMR/BNQR runs use seconds + if (fPlotType == MSR_PLOT_BNMR) { xAxisTitle = TString("time (s)"); } else { xAxisTitle = TString("time (#mus)"); From dfea8334b78612d5bb900f966a6d5a06150fd871 Mon Sep 17 00:00:00 2001 From: Zaher Salman Date: Fri, 24 Aug 2018 13:03:28 +0200 Subject: [PATCH 12/20] Possible fix for alpha problem for beta-NMR asymmetry calculation. To be tested. --- src/classes/PMusrCanvas.cpp | 24 ++++++++++++------------ src/classes/PRunAsymmetryBNMR.cpp | 14 +++++++------- src/classes/PRunListCollection.cpp | 2 +- 3 files changed, 20 insertions(+), 20 deletions(-) diff --git a/src/classes/PMusrCanvas.cpp b/src/classes/PMusrCanvas.cpp index a5a3ec2d..9e607797 100644 --- a/src/classes/PMusrCanvas.cpp +++ b/src/classes/PMusrCanvas.cpp @@ -718,6 +718,18 @@ void PMusrCanvas::UpdateDataTheoryPad() // handle data HandleDataSet(i, runNo, data); break; + case MSR_FITTYPE_BNMR: + data = fRunList->GetAsymmetryBNMR(runNo, PRunListCollection::kRunNo); + if (!data) { // something wrong + fValid = false; + // error message + cerr << endl << ">> PMusrCanvas::UpdateDataTheoryPad(): **ERROR** couldn't obtain run no " << runNo << " for a beta-NMR asymmetry plot"; + cerr << endl; + return; + } + // handle data + HandleDataSet(i, runNo, data); + break; case MSR_FITTYPE_ASYM_RRF: data = fRunList->GetAsymmetryRRF(runNo, PRunListCollection::kRunNo); if (!data) { // something wrong @@ -742,18 +754,6 @@ void PMusrCanvas::UpdateDataTheoryPad() // handle data HandleDataSet(i, runNo, data); break; - case MSR_FITTYPE_BNMR: - data = fRunList->GetAsymmetryBNMR(runNo, PRunListCollection::kRunNo); - if (!data) { // something wrong - fValid = false; - // error message - cerr << endl << ">> PMusrCanvas::UpdateDataTheoryPad(): **ERROR** couldn't obtain run no " << runNo << " for a beta-NMR asymmetry plot"; - cerr << endl; - return; - } - // handle data - HandleDataSet(i, runNo, data); - break; case MSR_FITTYPE_NON_MUSR: data = fRunList->GetNonMusr(runNo, PRunListCollection::kRunNo); if (!data) { // something wrong diff --git a/src/classes/PRunAsymmetryBNMR.cpp b/src/classes/PRunAsymmetryBNMR.cpp index f9aa4cfe..0de439aa 100644 --- a/src/classes/PRunAsymmetryBNMR.cpp +++ b/src/classes/PRunAsymmetryBNMR.cpp @@ -222,19 +222,19 @@ Double_t PRunAsymmetryBNMR::CalcChiSquare(const std::vector& par) 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)); + f = fTheory->Func(time, par, fFuncValues)/2; + asymFcnValue = (f*(a+1.0)-(a-1.0))/((a+1.0)-f*(a-1.0)) - (-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)); + f = fTheory->Func(time, par, fFuncValues)/2; + asymFcnValue = f*(b+1.0)/(2.0-f*(b-1.0))-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)); + f = fTheory->Func(time, par, fFuncValues)/2; + asymFcnValue = (f*(a*b+1.0)-(a-1.0))/((a+1.0)-f*(a*b-1.0))-(-f*(a*b+1.0)-(a-1.0))/((a+1.0)+f*(a*b-1.0)); break; default: asymFcnValue = 0.0; @@ -1276,7 +1276,7 @@ Bool_t PRunAsymmetryBNMR::PrepareViewData(PRawRunData* runData, UInt_t histoNo[2 efm = forwardmPacked.GetError()->at(i); ebm = backwardmPacked.GetError()->at(i); // check that there are indeed bins - if (fp+bp != 0.0) + if (fp+bp != 0.0) asym = (alpha*fp-bp) / (alpha*beta*fp+bp) - (alpha*fm-bm) / (alpha*beta*fm+bm); else asym = 0.0; diff --git a/src/classes/PRunListCollection.cpp b/src/classes/PRunListCollection.cpp index 6733062c..4fc9e032 100644 --- a/src/classes/PRunListCollection.cpp +++ b/src/classes/PRunListCollection.cpp @@ -972,7 +972,7 @@ PRunData* PRunListCollection::GetAsymmetry(UInt_t index, EDataSwitch tag) // GetAsymmetryBNMR (public) //-------------------------------------------------------------------------- /** - *

Get a processed asymmetry data set. + *

Get a processed asymmetry from beta-NMR data set. * * return: * - pointer to the run data set (processed data) if data set is found From b5f5d26e9c571832ff9308e0d3bb842a833d3189 Mon Sep 17 00:00:00 2001 From: Zaher Salman Date: Fri, 24 Aug 2018 14:08:28 +0200 Subject: [PATCH 13/20] Added estimation for the value of alpha from integrated counts. --- src/classes/PRunAsymmetryBNMR.cpp | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/src/classes/PRunAsymmetryBNMR.cpp b/src/classes/PRunAsymmetryBNMR.cpp index 0de439aa..5fa90517 100644 --- a/src/classes/PRunAsymmetryBNMR.cpp +++ b/src/classes/PRunAsymmetryBNMR.cpp @@ -903,8 +903,15 @@ Bool_t PRunAsymmetryBNMR::PrepareFitData() Double_t errorp = 0.0; Double_t valuem = 0.0; Double_t errorm = 0.0; + + Double_t SumF[2] = {0.0, 0.0}; + Double_t SumB[2] = {0.0, 0.0}; + Double_t alphaest = 1; + // forward for (Int_t i=fGoodBins[0]; i> PRunAsymmetryBNMR::PrepareFitData(): alpha estimate=" << alphaest << endl; + // check if packed forward and backward hist have the same size, otherwise take the minimum size UInt_t noOfBins = forwardpPacked.GetValue()->size(); if (forwardpPacked.GetValue()->size() != backwardpPacked.GetValue()->size()) { From f81c6aa2bf76647c9c543e6e5f5f919c64054717 Mon Sep 17 00:00:00 2001 From: Zaher Salman Date: Fri, 24 Aug 2018 17:42:23 +0200 Subject: [PATCH 14/20] Deal with missing alpha line for beta-NMR. --- src/classes/PMsrHandler.cpp | 12 ++++++------ src/classes/PRunAsymmetryBNMR.cpp | 31 +++++++++++++++---------------- 2 files changed, 21 insertions(+), 22 deletions(-) diff --git a/src/classes/PMsrHandler.cpp b/src/classes/PMsrHandler.cpp index 031c0088..08faad67 100644 --- a/src/classes/PMsrHandler.cpp +++ b/src/classes/PMsrHandler.cpp @@ -5542,12 +5542,12 @@ Bool_t PMsrHandler::CheckRunBlockIntegrity() break; case PRUN_ASYMMETRY_BNMR: // check alpha - if ((fRuns[i].GetAlphaParamNo() == -1) && !fFourierOnly) { - cerr << endl << ">> PMsrHandler::CheckRunBlockIntegrity(): **ERROR** in RUN block number " << i+1; - cerr << endl << ">> alpha parameter number missing which is needed for an asymmetry fit."; - cerr << endl << ">> Consider to check the manual ;-)" << endl; - return false; - } + // if ((fRuns[i].GetAlphaParamNo() == -1) && !fFourierOnly) { + // cerr << endl << ">> PMsrHandler::CheckRunBlockIntegrity(): **ERROR** in RUN block number " << i+1; + // cerr << endl << ">> alpha parameter number missing which is needed for an asymmetry fit."; + // cerr << endl << ">> Consider to check the manual ;-)" << endl; + // return false; + // } // check that there is a forward parameter number if (fRuns[i].GetForwardHistoNo() == -1) { cerr << endl << ">> PMsrHandler::CheckRunBlockIntegrity(): **ERROR** in RUN block number " << i+1; diff --git a/src/classes/PRunAsymmetryBNMR.cpp b/src/classes/PRunAsymmetryBNMR.cpp index 5fa90517..83d27d2a 100644 --- a/src/classes/PRunAsymmetryBNMR.cpp +++ b/src/classes/PRunAsymmetryBNMR.cpp @@ -105,26 +105,25 @@ PRunAsymmetryBNMR::PRunAsymmetryBNMR(PMsrHandler *msrInfo, PRunDataHandler *rawD PMsrParamList *param = msrInfo->GetMsrParamList(); // check if alpha is given + Bool_t alphaFixedToOne = false; 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** no alpha parameter given! This is needed for an asymmetry fit!"; + // cerr << endl; + // fValid = false; + // return; + alphaFixedToOne = true; + } else if ((fRunInfo->GetAlphaParamNo() < 0) || (fRunInfo->GetAlphaParamNo() > (Int_t)param->size())) { // check if alpha parameter is within proper bounds 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; + } else { // check if alpha is fixed + if (((*param)[fRunInfo->GetAlphaParamNo()-1].fStep == 0.0) && + ((*param)[fRunInfo->GetAlphaParamNo()-1].fValue == 1.0)) + alphaFixedToOne = true; } - // 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 @@ -222,18 +221,18 @@ Double_t PRunAsymmetryBNMR::CalcChiSquare(const std::vector& par) break; case 2: // alpha != 1, beta == 1 a = par[fRunInfo->GetAlphaParamNo()-1]; - f = fTheory->Func(time, par, fFuncValues)/2; + f = fTheory->Func(time, par, fFuncValues)/2.0; asymFcnValue = (f*(a+1.0)-(a-1.0))/((a+1.0)-f*(a-1.0)) - (-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)/2; + f = fTheory->Func(time, par, fFuncValues)/2.0; asymFcnValue = f*(b+1.0)/(2.0-f*(b-1.0))-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)/2; + f = fTheory->Func(time, par, fFuncValues)/2.0; asymFcnValue = (f*(a*b+1.0)-(a-1.0))/((a+1.0)-f*(a*b-1.0))-(-f*(a*b+1.0)-(a-1.0))/((a+1.0)+f*(a*b-1.0)); break; default: From 8c2537efae09826cd2c6f9aa0f755a3d282dd072 Mon Sep 17 00:00:00 2001 From: Zaher Salman Date: Sat, 25 Aug 2018 23:06:57 +0200 Subject: [PATCH 15/20] Estimate alpha for beta-NMR asymmetry calculation from the ratio of B/F counts if no alpha line is given in the run block. --- src/classes/PMusr.cpp | 13 ++++ src/classes/PRunAsymmetryBNMR.cpp | 123 ++++++++++++++++++++++++------ src/include/PMusr.h | 5 ++ src/include/PRunAsymmetryBNMR.h | 1 + 4 files changed, 120 insertions(+), 22 deletions(-) diff --git a/src/classes/PMusr.cpp b/src/classes/PMusr.cpp index a84f4765..6a2142df 100644 --- a/src/classes/PMusr.cpp +++ b/src/classes/PMusr.cpp @@ -1879,6 +1879,19 @@ void PMsrRunBlock::SetMapGlobal(UInt_t idx, Int_t ival) return; } +//-------------------------------------------------------------------------- +// SetEstimatedAlpha (public) +//-------------------------------------------------------------------------- +/** + *

set the value of estimated alpha at position idx + * + * \param alpha is the estimated value + */ +void PMsrRunBlock::SetEstimatedAlpha(Double_t dval) +{ + fAlpha = dval; +} + //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ // implementation PStringNumberList //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ diff --git a/src/classes/PRunAsymmetryBNMR.cpp b/src/classes/PRunAsymmetryBNMR.cpp index 83d27d2a..8afa4937 100644 --- a/src/classes/PRunAsymmetryBNMR.cpp +++ b/src/classes/PRunAsymmetryBNMR.cpp @@ -106,12 +106,13 @@ PRunAsymmetryBNMR::PRunAsymmetryBNMR(PMsrHandler *msrInfo, PRunDataHandler *rawD // check if alpha is given Bool_t alphaFixedToOne = false; + Bool_t alphaSetDefault = false; 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; - alphaFixedToOne = true; + alphaSetDefault = true; } else if ((fRunInfo->GetAlphaParamNo() < 0) || (fRunInfo->GetAlphaParamNo() > (Int_t)param->size())) { // check if alpha parameter is within proper bounds cerr << endl << ">> PRunAsymmetryBNMR::PRunAsymmetryBNMR(): **ERROR** alpha parameter no = " << fRunInfo->GetAlphaParamNo(); cerr << endl << ">> This is out of bound, since there are only " << param->size() << " parameters."; @@ -141,12 +142,16 @@ PRunAsymmetryBNMR::PRunAsymmetryBNMR(PMsrHandler *msrInfo, PRunDataHandler *rawD } // set fAlphaBetaTag - if (alphaFixedToOne && betaFixedToOne) // alpha == 1, beta == 1 + if (alphaFixedToOne && betaFixedToOne) // alpha == 1, beta == 1 fAlphaBetaTag = 1; - else if (!alphaFixedToOne && betaFixedToOne) // alpha != 1, beta == 1 + else if (!alphaFixedToOne && betaFixedToOne & !alphaSetDefault) // alpha != 1, beta == 1 fAlphaBetaTag = 2; - else if (alphaFixedToOne && !betaFixedToOne) // alpha == 1, beta != 1 + else if (alphaFixedToOne && !betaFixedToOne) // alpha == 1, beta != 1 fAlphaBetaTag = 3; + else if (!alphaFixedToOne && betaFixedToOne & alphaSetDefault) // alpha ??, beta == 1 + fAlphaBetaTag = 5; + else if (!alphaFixedToOne && !betaFixedToOne & alphaSetDefault) // alpha ??, beta != 1 + fAlphaBetaTag = 6; else fAlphaBetaTag = 4; @@ -198,7 +203,7 @@ Double_t PRunAsymmetryBNMR::CalcChiSquare(const std::vector& par) } // calculate chi square - Double_t time(1.0); + Double_t time(1.0),alphaest; Int_t i; // Calculate the theory function once to ensure one function evaluation for the current set of parameters. @@ -206,6 +211,7 @@ Double_t PRunAsymmetryBNMR::CalcChiSquare(const std::vector& par) // 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); + alphaest = fRunInfo->GetEstimatedAlpha(); #ifdef HAVE_GOMP Int_t chunk = (fEndTimeBin - fStartTimeBin)/omp_get_num_procs(); @@ -235,6 +241,17 @@ Double_t PRunAsymmetryBNMR::CalcChiSquare(const std::vector& par) f = fTheory->Func(time, par, fFuncValues)/2.0; asymFcnValue = (f*(a*b+1.0)-(a-1.0))/((a+1.0)-f*(a*b-1.0))-(-f*(a*b+1.0)-(a-1.0))/((a+1.0)+f*(a*b-1.0)); break; + case 5: // alpha ?? , beta == 1 + a = alphaest; + f = fTheory->Func(time, par, fFuncValues)/2.0; + asymFcnValue = (f*(a+1.0)-(a-1.0))/((a+1.0)-f*(a-1.0)) - (-f*(a+1.0)-(a-1.0))/((a+1.0)+f*(a-1.0)); + break; + case 6: // alpha ??, beta != 1 + a = alphaest; + b = par[fRunInfo->GetBetaParamNo()-1]; + f = fTheory->Func(time, par, fFuncValues)/2.0; + asymFcnValue = (f*(a*b+1.0)-(a-1.0))/((a+1.0)-f*(a*b-1.0))-(-f*(a*b+1.0)-(a-1.0))/((a+1.0)+f*(a*b-1.0)); + break; default: asymFcnValue = 0.0; break; @@ -425,8 +442,12 @@ void PRunAsymmetryBNMR::CalcTheory() // calculate asymmetry Double_t asymFcnValue = 0.0; - Double_t a, b, f; + Double_t a, b, f, alphaest; Double_t time; + + // Get estimated alpha + alphaest = fRunInfo->GetEstimatedAlpha(); + for (UInt_t i=0; isize(); i++) { time = fData.GetDataTimeStart() + (Double_t)i*fData.GetDataTimeStep(); switch (fAlphaBetaTag) { @@ -436,18 +457,29 @@ void PRunAsymmetryBNMR::CalcTheory() 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)); + asymFcnValue = (f*(a+1.0)-(a-1.0))/((a+1.0)-f*(a-1.0)) - (-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)); + asymFcnValue = f*(b+1.0)/(2.0-f*(b-1.0))-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)); + asymFcnValue = (f*(a*b+1.0)-(a-1.0))/((a+1.0)-f*(a*b-1.0))-(-f*(a*b+1.0)-(a-1.0))/((a+1.0)+f*(a*b-1.0)); + break; + case 5: // alpha ?? , beta == 1 + a = alphaest; + f = fTheory->Func(time, par, fFuncValues); + asymFcnValue = (f*(a+1.0)-(a-1.0))/((a+1.0)-f*(a-1.0)) - (-f*(a+1.0)-(a-1.0))/((a+1.0)+f*(a-1.0)); + break; + case 6: // alpha ??, beta != 1 + a = alphaest; + 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))-(-f*(a*b+1.0)-(a-1.0))/((a+1.0)+f*(a*b-1.0)); break; default: asymFcnValue = 0.0; @@ -871,6 +903,14 @@ Bool_t PRunAsymmetryBNMR::SubtractEstimatedBkg() fRunInfo->SetBkgEstimated(bkgm[0], 3); fRunInfo->SetBkgEstimated(bkgm[1], 4); + cout << "I am here 1 " << endl; + // Get estimate for alpha once + Double_t alpha; + alpha = EstimateAlpha(); + cout << "I am here 2 " << alpha << endl; + fRunInfo->SetEstimatedAlpha(alpha); + + return true; } @@ -903,14 +943,8 @@ Bool_t PRunAsymmetryBNMR::PrepareFitData() Double_t valuem = 0.0; Double_t errorm = 0.0; - Double_t SumF[2] = {0.0, 0.0}; - Double_t SumB[2] = {0.0, 0.0}; - Double_t alphaest = 1; - // forward for (Int_t i=fGoodBins[0]; i> PRunAsymmetryBNMR::PrepareFitData(): alpha estimate=" << alphaest << endl; - // check if packed forward and backward hist have the same size, otherwise take the minimum size UInt_t noOfBins = forwardpPacked.GetValue()->size(); if (forwardpPacked.GetValue()->size() != backwardpPacked.GetValue()->size()) { @@ -1249,13 +1277,16 @@ Bool_t PRunAsymmetryBNMR::PrepareViewData(PRawRunData* runData, UInt_t histoNo[2 // form asymmetry including error propagation Double_t asym; - Double_t fp, bp, efp, ebp, alpha = 1.0, beta = 1.0; + Double_t fp, bp, efp, ebp, alpha, beta = 1.0; Double_t fm, bm, efm, ebm; // 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 estimate for alpha once + alpha = EstimateAlpha(); + // get the proper alpha and beta switch (fAlphaBetaTag) { case 1: // alpha == 1, beta == 1 @@ -1274,6 +1305,14 @@ Bool_t PRunAsymmetryBNMR::PrepareViewData(PRawRunData* runData, UInt_t histoNo[2 alpha = par[fRunInfo->GetAlphaParamNo()-1]; beta = par[fRunInfo->GetBetaParamNo()-1]; break; + case 5: // alpha ?? , beta == 1 + // use estimated value + beta = 1.0; + break; + case 6: // alpha ??, beta != 1 + // use estimated value + beta = par[fRunInfo->GetBetaParamNo()-1]; + break; default: break; } @@ -1701,3 +1740,43 @@ void PRunAsymmetryBNMR::GetProperFitRange(PMsrGlobalBlock *globalBlock) cerr << ">> Will set it to fgb/lgb which given in time is: " << fFitStartTime << "..." << fFitEndTime << " (usec)" << endl; } } + + +//-------------------------------------------------------------------------- +// EstimateAlpha (private) +//-------------------------------------------------------------------------- +/** + *

Get an estimate for alpha from the forward and backward histograms + * + * \param globalBlock pointer to the GLOBAL block information form the msr-file. + */ +Double_t PRunAsymmetryBNMR::EstimateAlpha() +{ + + Double_t SumF[2] = {0.0, 0.0}; + Double_t SumB[2] = {0.0, 0.0}; + Double_t alpha = 1; + + + // forward + for (Int_t i=0; i> PRunAsymmetryBNMR::EstimateAlpha(): alpha estimate=" << alpha << endl; + + return alpha; +} diff --git a/src/include/PMusr.h b/src/include/PMusr.h index cd9d6375..357cf6b4 100644 --- a/src/include/PMusr.h +++ b/src/include/PMusr.h @@ -570,6 +570,7 @@ class PMsrGlobalBlock { virtual Double_t GetFitRange(UInt_t idx); virtual Int_t GetFitRangeOffset(UInt_t idx); virtual Int_t GetPacking() { return fPacking; } + virtual Double_t GetEstimatedAlpha() { return fAlpha; } virtual void SetGlobalPresent(Bool_t bval) { fGlobalPresent = bval; } virtual void SetRRFFreq(Double_t freq, const char *unit); @@ -598,6 +599,7 @@ class PMsrGlobalBlock { Double_t fFitRange[2]; ///< fit range in (us) Int_t fFitRangeOffset[2]; ///< if fit range is given in bins it can have the form fit fgb+n0 lgb-n1. This variable holds the n0 and n1. Int_t fPacking; ///< packing/rebinning + Double_t fAlpha; ///< estimated alpha value from F/B counts }; //------------------------------------------------------------- @@ -645,6 +647,7 @@ class PMsrRunBlock { virtual Double_t GetFitRange(UInt_t idx); virtual Int_t GetFitRangeOffset(UInt_t idx); virtual Int_t GetPacking() { return fPacking; } + virtual Double_t GetEstimatedAlpha() { return fAlpha; } virtual Int_t GetXDataIndex() { return fXYDataIndex[0]; } virtual Int_t GetYDataIndex() { return fXYDataIndex[1]; } virtual TString* GetXDataLabel() { return &fXYDataLabel[0]; } @@ -667,6 +670,7 @@ class PMsrRunBlock { virtual void SetForwardHistoNo(Int_t histoNo, Int_t idx=-1); virtual void SetBackwardHistoNo(Int_t histoNo, Int_t idx=-1); virtual void SetBkgEstimated(Double_t dval, Int_t idx); + virtual void SetEstimatedAlpha(Double_t dval); virtual void SetBkgFix(Double_t dval, Int_t idx); virtual void SetBkgRange(Int_t ival, Int_t idx); virtual void SetDataRange(Int_t ival, Int_t idx); @@ -707,6 +711,7 @@ class PMsrRunBlock { Bool_t fFitRangeInBins; ///< flag telling if fit range is given in time or in bins Double_t fFitRange[2]; ///< fit range in (us) Int_t fFitRangeOffset[2]; ///< if fit range is given in bins it can have the form fit fgb+n0 lgb-n1. This variable holds the n0 and n1. + Double_t fAlpha; ///< estimated alpha value from F/B counts Int_t fPacking; ///< packing/rebinning Int_t fXYDataIndex[2]; ///< used to get the data indices when using db-files (fit type 8) TString fXYDataLabel[2]; ///< used to get the indices via labels when using db-files (fit type 8) diff --git a/src/include/PRunAsymmetryBNMR.h b/src/include/PRunAsymmetryBNMR.h index a1fc5377..5ed85c9f 100644 --- a/src/include/PRunAsymmetryBNMR.h +++ b/src/include/PRunAsymmetryBNMR.h @@ -87,6 +87,7 @@ class PRunAsymmetryBNMR : public PRunBase 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); + virtual Double_t EstimateAlpha(); }; #endif // _PRUNASYMMETRYBNMR_H_ From e19289d0c38d09787172339c0e3669bf20dd5e2c Mon Sep 17 00:00:00 2001 From: Zaher Salman Date: Sat, 25 Aug 2018 23:51:08 +0200 Subject: [PATCH 16/20] Cleanup --- src/classes/PRunAsymmetryBNMR.cpp | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/src/classes/PRunAsymmetryBNMR.cpp b/src/classes/PRunAsymmetryBNMR.cpp index 8afa4937..63b3c271 100644 --- a/src/classes/PRunAsymmetryBNMR.cpp +++ b/src/classes/PRunAsymmetryBNMR.cpp @@ -860,10 +860,10 @@ Bool_t PRunAsymmetryBNMR::SubtractEstimatedBkg() } errBkgp[1] = TMath::Sqrt(bkgp[1])/(end[1] - start[1] + 1); bkgp[1] /= static_cast(end[1] - start[1] + 1); - cout << endl << ">> estimated pos hel backward histo background: " << bkgp[1] << endl; + cout << endl << ">> estimated pos hel backward histo background: " << bkgp[1]; errBkgm[1] = TMath::Sqrt(bkgm[1])/(end[1] - start[1] + 1); bkgm[1] /= static_cast(end[1] - start[1] + 1); - cout << endl << ">> estimated neg hel backward histo background: " << bkgm[1] << endl; + cout << endl << ">> estimated neg hel backward histo background: " << bkgm[1]; // correct error for forward, backward Double_t errVal = 0.0; @@ -903,11 +903,9 @@ Bool_t PRunAsymmetryBNMR::SubtractEstimatedBkg() fRunInfo->SetBkgEstimated(bkgm[0], 3); fRunInfo->SetBkgEstimated(bkgm[1], 4); - cout << "I am here 1 " << endl; // Get estimate for alpha once Double_t alpha; alpha = EstimateAlpha(); - cout << "I am here 2 " << alpha << endl; fRunInfo->SetEstimatedAlpha(alpha); From e428554f7e9c24d1921fa9491b6b8b218e311100 Mon Sep 17 00:00:00 2001 From: Zaher Salman Date: Sun, 26 Aug 2018 11:33:13 +0200 Subject: [PATCH 17/20] Added documentation relevant to beta-NMR asymmetry fits, fittype 5. --- doc/html/_sources/user-manual.txt | 78 +++++++++++++++++++++++++++---- 1 file changed, 69 insertions(+), 9 deletions(-) diff --git a/doc/html/_sources/user-manual.txt b/doc/html/_sources/user-manual.txt index 4fd2b4d2..fb1f67f3 100644 --- a/doc/html/_sources/user-manual.txt +++ b/doc/html/_sources/user-manual.txt @@ -845,7 +845,7 @@ Currently the supported GLOBAL block entries are: * ``rrf_freq`` for fittype 1, 3 * ``rrf_packing`` for fittype 1, 3 * ``rrf_phase`` for fittype 1, 3 -* ``packing`` for fittype 0, 2, 4 +* ``packing`` for fittype 0, 2, 4, 5 For a detailed discussion of these entries see the section :ref:`RUN block `. @@ -1039,6 +1039,8 @@ In order to describe the operations needed for fitting and plotting, quite some Asymmetry RRF Fit (only for online analysis) **4** MuMinus Fit. This is a single histogram fit especially for negative muon |mgr|\SR + **5** + beta-NMR Asymmetry Fit **8** Non-|mgr|\SR Fit @@ -1051,8 +1053,8 @@ In order to describe the operations needed for fitting and plotting, quite some .. index:: alpha-beta .. _msr-alpha-beta: -**alpha, beta** (fit type 2, 3) - These parameters are used to correct the asymmetry for different detector efficiencies, solid angles and initial asymmetries. They are defined as :math:`\alpha = N_{0,b}/N_{0,f}` and :math:`\beta = A_{0,b}/A_{0,f}`. If the parameters are not specified in the :ref:`RUN block `, for each one the value of 1 is assumed. Example for alpha with fit parameter number 1: +**alpha, beta** (fit type 2, 3, 5) + These parameters are used to correct the asymmetry for different detector efficiencies, solid angles and initial asymmetries. They are defined as :math:`\alpha = N_{0,b}/N_{0,f}` and :math:`\beta = A_{0,b}/A_{0,f}`. If the parameters are not specified in the :ref:`RUN block `, for each one the value of 1 is assumed (for fittype 5 alpha is estimated from the ration of sum of Bp+Bm and Fp+Fm). Example for alpha with fit parameter number 1: :: @@ -1131,10 +1133,19 @@ In order to describe the operations needed for fitting and plotting, quite some forward 1-3 backward 7-9 +**forward, backward** (fit type 5) + Numbers of the histograms in the data file that should be taken to calculate the asymmetry. Two forward and backward histograms should be given indicationg positive and negative helicities. The asymmetry from opposite helicities will be subtracted. Examples: + + :: + + # build forward/backward asymmetry with histogram 1 and 3 then subtract asymmetry built with histograms 2 and 4 + forward 1 2 + backward 3 4 + .. index:: backgr.fix .. _msr-backgr.fix: -**backgr.fix** (fit types 0, 1, 2, 3) +**backgr.fix** (fit types 0, 1, 2, 3, 5) A fixed constant background in counts per nanosecond or per bin (see :ref:`below `) may be given at this point. The background is specified for all histograms in the order :math:`B_f B_b [B_r B_l]`. If this keyword is present, *any* information on a ``background`` line is ignored. @@ -1152,7 +1163,7 @@ In order to describe the operations needed for fitting and plotting, quite some .. index:: background-asymmetry .. _msr-background-asymmetry: -**background** (fit types 2, 3) +**background** (fit types 2, 3, 5) The numbers of the first and the last channel of an interval from which the constant background should be calculated are specified here. For all the histograms this is done together in the following order: :math:`k_{f,\rm first} k_{f,\rm last} k_{b,\rm first} k_{b, \rm last} [k_{r,\rm first} k_{r,\rm last} k_{l,\rm first} k_{l,\rm last}]`. In case histograms are being grouped, the specified channels are interpreted with respect to the first histograms. Example: @@ -1176,7 +1187,7 @@ In order to describe the operations needed for fitting and plotting, quite some .. index:: data-asymmetry .. _msr-data-asymmetry: -**data** (fit type 2, 3) +**data** (fit type 2, 3, 5) The numbers of the first and the last channel of an interval from which the data is taken are specified here. Typically these channels are referred to as first good bin / last good bin (fgb/lgb). For all the histograms this is done together in the following order: :math:`k_{f,\rm first} k_{f,\rm last} k_{b,\rm first} k_{b, \rm last} [k_{r,\rm first} k_{r,\rm last} k_{l,\rm first} k_{l,\rm last}]`. @@ -1200,8 +1211,8 @@ In order to describe the operations needed for fitting and plotting, quite some .. index:: t0-asymmetry .. _msr-t0-asymmetry: -**t0** (fit type 2, 3) - The numbers of time-zero channels of the histograms in the order :math:`t_{0,f} t_{0,b}`. Example: +**t0** (fit type 2, 3, 5) + The numbers of time-zero channels of the histograms in the order :math:`t_{0,f} t_{0,b}`. For fit type 5, the time-zero is the channel of the start of beam pulse. Example: :: @@ -1219,7 +1230,7 @@ In order to describe the operations needed for fitting and plotting, quite some .. index:: addt0-asymmetry .. _msr-addt0-asymmetry: -**addt0** (fit type 2, 3) +**addt0** (fit type 2, 3, 5) The numbers of time-zero channels of the histograms in the order :math:`t_{0,f} t_{0,b} [t_{0,r} t_{0,l}]`. If grouping of histograms is present (see :ref:`forward `) the same syntax as for :ref:`t0 ` applies. If one addt0 is given, the total number of addt0's needs to be equal to the total number of :ref:`ADDRUN `\'s! @@ -1782,6 +1793,55 @@ where :math:`i` runs over the different lifetime channels of :math:`\mu^{-}`, an Since MuMinus is quite generic, the full functional depends has to be written in the :ref:`THEORY Block `. +.. index:: bnmr-asymmetry-fit +.. _bnmr-asymmetry-fit: + +beta-NMR Asymmetry Fit (fit type 5) ++++++++++++++++++++++++++++++++++++ + +For a beta-NMR asymmetry fit (fit type 5) four histograms are needed, two for positive and two for negative helicities. These are given by the :ref:`forward ` and :ref:`backward ` keywords +in the :ref:`RUN block `. Additionally, the parameters :ref:`alpha ` and :ref:`beta ` which relate the detector +efficiencies, solid angles and initial asymmetries of the two detectors can be supplied. The constant background for the two histograms is either given by +:ref:`background-determined intervals ` or specified through :ref:`backgr.fix ` in the :ref:`RUN-block `. + +The experimental asymmetry :math:`a(k)` then is calculated from the four histograms: + +.. math:: + + a(k)=\frac{\left[N_{\mathrm{fp}}(k)-B_{\mathrm{fp}}\right]-\left[N_{\mathrm{bp}}(k)-B_{\mathrm{bp}}\right]}{\left[N_{\mathrm{fp}}(k)-B_{\mathrm{fp}}\right]+\left[N_{\mathrm{bp}}(k)-B_{\mathrm{bp}}\right]} + - \frac{\left[N_{\mathrm{fm}}(k)-B_{\mathrm{fm}}\right]-\left[N_{\mathrm{bm}}(k)-B_{\mathrm{bm}}\right]}{\left[N_{\mathrm{fm}}(k)-B_{\mathrm{fm}}\right]+\left[N_{\mathrm{bm}}(k)-B_{\mathrm{bm}}\right]}, + +with + + * :math:`N_{\mathrm{fp}}(k)`: counts in the **forward** histogram channel with positive helicity :math:`k` + * :math:`N_{\mathrm{bp}}(k)`: counts in the **backward** histogram channel with positive helicity :math:`k` + * :math:` B_{\mathrm{fp}}`: constant background in the **forward** histogram with positive helicity (RUN block: :ref:`backgr.fix ` or :ref:`background `) + * :math:` B_{\mathrm{bp}}`: constant background in the **backward** histogram with positive helicity (RUN block: :ref:`backgr.fix ` or :ref:`background `) + * :math:`N_{\mathrm{fm}}(k)`: counts in the **forward** histogram channel with negative helicity :math:`k` + * :math:`N_{\mathrm{bm}}(k)`: counts in the **backward** histogram channel with negative helicity :math:`k` + * :math:` B_{\mathrm{fm}}`: constant background in the **forward** histogram with negative helicity (RUN block: :ref:`backgr.fix ` or :ref:`background `) + * :math:` B_{\mathrm{bm}}`: constant background in the **backward** histogram with negative helicity (RUN block: :ref:`backgr.fix ` or :ref:`background `) + +This theoretical asymmetry :math:`a(t)` is used to fit the function + +.. math:: + + a(t)=\frac{(\alpha\beta +1)A(t)-(\alpha -1)}{(\alpha +1)-(\alpha\beta -1)A(t)} - \frac{(\alpha -1)-(\alpha\beta 1)A(t)}{(\alpha +1)+(\alpha\beta -1)Am(t)}, + +where + + * :math:`\alpha`: accounts for the different detector efficiencies and solid angles (RUN block: :ref:`alpha `). + * :math:`\beta`: accounts for the different detector asymmetries (RUN block: :ref:`beta `). + * :math:`A(t)`: is the depolarization function as given in the :ref:`THEORY block `. + +For the graphical representation in plot type 5 the equation above is rearranged to get :math:`A(t)`: + +.. math:: + + A(t)=\frac{(\alpha -1)+(\alpha +1)a(t)}{(\alpha\beta +1)+(\alpha\beta -1)a(t)}-\frac{(\alpha +1)a(t)-(\alpha -1)}{(\alpha\beta +1)+(1-\alpha\beta)a(t)}=\frac{\alpha\left[N_{\mathrm{fp}}(t)-B_{\mathrm{fp}}\right]-\left[N_{\mathrm{bp}}(t)-B_{\mathrm{bp}}\right]}{\alpha\beta\left[N_{\mathrm{fp}}(t)-B_{\mathrm{fp}}\right]+\left[N_{\mathrm{bp}}(t)-B_{\mathrm{bp}}\right]} -\frac{\alpha\left[N_{\mathrm{fm}}(t)-B_{\mathrm{fm}}\right]-\left[N_{\mathrm{bm}}(t)-B_{\mathrm{bm}}\right]}{\alpha\beta\left[N_{\mathrm{fm}}(t)-B_{\mathrm{fm}}\right]+\left[N_{\mathrm{bm}}(t)-B_{\mathrm{bm}}\right]} + +and plotted together with the function given in the THEORY block. + .. index:: non-musr-fit .. _non-musr-fit: From 2cd81cf2354b7f3cff89a62c5f99a7822614bcc9 Mon Sep 17 00:00:00 2001 From: Zaher Salman Date: Mon, 15 Oct 2018 10:16:41 +0200 Subject: [PATCH 18/20] More recent working version --- src/external/MuSRFitGUI/MSR.pm | 50 +++++++++++++++++++++++----------- 1 file changed, 34 insertions(+), 16 deletions(-) diff --git a/src/external/MuSRFitGUI/MSR.pm b/src/external/MuSRFitGUI/MSR.pm index 87b4d120..cb8c6086 100644 --- a/src/external/MuSRFitGUI/MSR.pm +++ b/src/external/MuSRFitGUI/MSR.pm @@ -27,10 +27,11 @@ my $SUMM_DIR="/afs/psi.ch/project/nemu/data/summ/"; my %DBDIR=("LEM","/afs/psi.ch/project/nemu/data/log/", "GPS","/afs/psi.ch/project/bulkmusr/olddata/list/", "Dolly","/afs/psi.ch/project/bulkmusr/olddata/list/", - "GPD","/afs/psi.ch/project/bulkmusr/olddata/list/", - "ALC","/afs/psi.ch/project/bulkmusr/olddata/list/", - "HAL","/afs/psi.ch/project/bulkmusr/olddata/list/", - "LTF","/afs/psi.ch/project/bulkmusr/olddata/list/"); + "GPD","/afs/psi.ch/project/bulkmusr/olddata/list/", + "ALC","/afs/psi.ch/project/bulkmusr/olddata/list/", + "HAL","/afs/psi.ch/project/bulkmusr/olddata/list/", + "LTF","/afs/psi.ch/project/bulkmusr/olddata/list/", + "ALL","/afs/psi.ch/project/bulkmusr/olddata/list/"); # Information available since my %MinYears=("LEM","2001", @@ -662,12 +663,13 @@ sub CreateTheory { # Start from this theory line for the different fitting functions my %THEORY = ( - "asymmetry", "Asy", - "simplExpo", "Lam", - "generExpo", "Lam Bet", - "simpleGss", "Sgm", - "statGssKT", "Sgm", - "statGssKTLF", "Frqg Sgm", + "asymmetry", "Asy", + "simplExpo", "Lam", + "generExpo", "Lam Bet", + "abragam", "Del Lam", + "simpleGss", "Sgm", + "statGssKT", "Sgm", + "statGssKTLF", "Frqg Sgm", "dynGssKTLF", "Frql Sgm Lam", "statExpKT", "Lam", "statExpKTLF", "Frq Aa", @@ -680,10 +682,9 @@ sub CreateTheory { "internFld", "Alp Phi Frq LamT LamL", "Bessel", "Phi Frq", "internBsl", "Alp Phi Frq LamT LamL", - "abragam", "Sgm gam", - "Meissner", "Phi Energy Field DeadLayer Lambda", - "skewedGss", "Phi Frq Sgmm Sgmp" - ); + "Meissner", "Phi Energy Field DeadLayer Lambda", + "skewedGss", "Phi Frq Sgmm Sgmp" + ); ####################################################################### # Generate the full THEORY block @@ -752,6 +753,14 @@ sub CreateTheory { $Parameters = join( $SPACE, $Parameters, $THEORY{'TFieldCos'} ); } + # Oscillationg Abragam function + elsif ( $FitType eq "AbragamCos" ) { + $T_Block = $T_Block . "\n" . "abragam " . $THEORY{'abragam'}; + $Parameters = join( $SPACE, $Parameters, $THEORY{'abragam'} ); + $T_Block = $T_Block . "\n" . "TFieldCos " . $THEORY{'TFieldCos'}; + $Parameters = join( $SPACE, $Parameters, $THEORY{'TFieldCos'} ); + } + # Oscillationg Bessel Gaussian elsif ( $FitType eq "GaussianBessel" ) { $T_Block = $T_Block . "\n" . "simpleGss " . $THEORY{'simpleGss'}; @@ -1015,8 +1024,16 @@ sub T0BgData { ($RV{"t0"},$RV{"Bg1"},$RV{"Bg2"},$RV{"Data1"},$RV{"Data2"})=split(/,/,$HistParams); } elsif ($BeamLine eq "GPS") { - my $HistParams=$GPS{$Hists[0]}; - ($RV{"t0"},$RV{"Bg1"},$RV{"Bg2"},$RV{"Data1"},$RV{"Data2"})=split(/,/,$HistParams); + my $HistParams=$GPS{$Hists[0]}; + ($RV{"t0"},$RV{"Bg1"},$RV{"Bg2"},$RV{"Data1"},$RV{"Data2"})=split(/,/,$HistParams); + } + elsif ($BeamLine eq "GPD") { + my $HistParams=$GPD{$Hists[0]}; + ($RV{"t0"},$RV{"Bg1"},$RV{"Bg2"},$RV{"Data1"},$RV{"Data2"})=split(/,/,$HistParams); + } + else { + my $HistParams=$GPS{$Hists[0]}; + ($RV{"t0"},$RV{"Bg1"},$RV{"Bg2"},$RV{"Data1"},$RV{"Data2"})=split(/,/,$HistParams); } return $RV{$Name}; @@ -1262,6 +1279,7 @@ sub ExportParams { push( @FitTypes, $All{"FitType$i"} ); } } + # Get theory block to determine the size of the table my ($Full_T_Block,$Paramcomp_ref)= MSR::CreateTheory(@FitTypes); # For now the line below does not work. Why? From 8207aea183703c59eec398946535a3baafaa3661 Mon Sep 17 00:00:00 2001 From: Zaher Salman Date: Mon, 15 Oct 2018 10:19:02 +0200 Subject: [PATCH 19/20] Changes to template files --- src/external/libBNMR/ExpRlx.msr | 2 +- src/external/libBNMR/SExpRlx.msr | 10 +++++----- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/external/libBNMR/ExpRlx.msr b/src/external/libBNMR/ExpRlx.msr index eb356304..59c94244 100644 --- a/src/external/libBNMR/ExpRlx.msr +++ b/src/external/libBNMR/ExpRlx.msr @@ -17,7 +17,7 @@ userFcn .libs/libBNMR.so ExpRlx 3 4 ############################################################### RUN 045674 BNMR TRIUMF MUD (name beamline institute data-file-format) fittype 5 (beta-NMR fit) -alpha 1 +#alpha 1 forward 3 5 backward 4 6 data 11 1000 11 1000 diff --git a/src/external/libBNMR/SExpRlx.msr b/src/external/libBNMR/SExpRlx.msr index a530e625..5ed46968 100644 --- a/src/external/libBNMR/SExpRlx.msr +++ b/src/external/libBNMR/SExpRlx.msr @@ -3,8 +3,8 @@ LaAlO3, SLR, 30G, 150K, Bias=23 kV, Ix=0A, Iy=+3A, aligned, beam on=1s FITPARAMETER ############################################################### # No Name Value Err Min Max - 1 Alpha 1 0 none - 2 Asy 0.1148 0.0096 none 0 0.2 + 1 Alpha 0.9 1.0 none + 2 Asy 0.115 0.011 none 0 0.2 3 T 1 0 none 4 Rlx 4.16 0.88 none 0 100 5 Beta 0.433 0.031 none 0.3 2 @@ -13,7 +13,7 @@ FITPARAMETER THEORY ############################################################### asymmetry 2 -userFcn .libs/libBNMR.so SExpRlx 3 4 5 +userFcn libBNMR SExpRlx 3 4 5 ############################################################### RUN 045674 BNMR TRIUMF MUD (name beamline institute data-file-format) @@ -50,5 +50,5 @@ plot POWER # REAL, IMAG, REAL_AND_IMAG, POWER, PHASE, PHASE_OPT_RE phase 8 #range FRQMIN FRQMAX ############################################################### -STATISTIC --- 2018-08-20 21:59:14 - chisq = 210.2, NDF = 172, chisq/NDF = 1.222196 +STATISTIC --- 2018-08-24 16:52:40 + chisq = 210.2, NDF = 171, chisq/NDF = 1.229343 From 6c88bce19c71e7562f1f7cec2da3d4e368d2756a Mon Sep 17 00:00:00 2001 From: Zaher Salman Date: Fri, 5 Apr 2019 10:14:16 +0200 Subject: [PATCH 20/20] Removed automake file to resolve conflict --- src/classes/Makefile.am | 122 ---------------------------------------- 1 file changed, 122 deletions(-) delete mode 100644 src/classes/Makefile.am diff --git a/src/classes/Makefile.am b/src/classes/Makefile.am deleted file mode 100644 index fa3a0fba..00000000 --- a/src/classes/Makefile.am +++ /dev/null @@ -1,122 +0,0 @@ -## 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 -