From 357d46aac4e9a5d30e69cc41fb65aa56cc5d9bdd Mon Sep 17 00:00:00 2001 From: Andreas Suter Date: Mon, 15 Oct 2018 16:09:17 +0200 Subject: [PATCH] allow multiple Fourier phase parameters for phase shifted real Fourier. Autophasing still missing. --- src/classes/PMsr2Data.cpp | 5 +- src/classes/PMsrHandler.cpp | 433 +++++++++++++++++++++++++++++--- src/classes/PMusrCanvas.cpp | 54 ++-- src/classes/PStartupHandler.cpp | 3 +- src/include/PMsrHandler.h | 5 + src/include/PMusr.h | 4 +- src/include/PMusrCanvas.h | 2 +- 7 files changed, 436 insertions(+), 70 deletions(-) diff --git a/src/classes/PMsr2Data.cpp b/src/classes/PMsr2Data.cpp index 39752f9d..4dd082dd 100644 --- a/src/classes/PMsr2Data.cpp +++ b/src/classes/PMsr2Data.cpp @@ -1025,12 +1025,11 @@ bool PMsr2Data::PrepareGlobalInputFile(unsigned int tempRun, const string &msrOu } // FOURIER block - in case a parameter is used for the phase - tempPar = fMsrHandler->GetMsrFourierList()->fPhaseParamNo; - if (tempPar > 0) { + if (fMsrHandler->GetMsrFourierList()->fPhaseParamNo.size() > 0) { // go through the whole parameter list ... for (unsigned int k(0); k < msrParamList->size(); ++k) { if (tempPar == msrParamList->at(k).fNo) { - fMsrHandler->GetMsrFourierList()->fPhaseParamNo = k + 1; + fMsrHandler->GetMsrFourierList()->fPhaseParamNo.push_back(k + 1); break; } } diff --git a/src/classes/PMsrHandler.cpp b/src/classes/PMsrHandler.cpp index 0f7f2008..d460dabc 100644 --- a/src/classes/PMsrHandler.cpp +++ b/src/classes/PMsrHandler.cpp @@ -326,6 +326,15 @@ Int_t PMsrHandler::ReadMsrFile() CheckLegacyLifetimecorrection(); // check if lifetimecorrection is found in RUN blocks, if yes transfer it to PLOT blocks } + // check if the given phases in the Fourier block are in agreement with the Plot block settings + if ((fFourier.fPhase.size() > 1) && (fPlots.size() > 0)) { + if (fFourier.fPhase.size() != fPlots[0].fRuns.size()) { + cerr << endl << ">> PMsrHandler::ReadMsrFile: **ERROR** if more than one phase is given in the Fourier block,"; + cerr << endl << ">> it needs to correspond to the number of runs in the Plot block!" << endl; + result = PMUSR_MSR_SYNTAX_ERROR; + } + } + // clean up fit_parameter.clear(); theory.clear(); @@ -1118,12 +1127,15 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t messages) fout << " # REAL, IMAG, REAL_AND_IMAG, POWER, PHASE, PHASE_OPT_REAL"; fout << endl; } else if (sstr.BeginsWith("phase")) { - if (fFourier.fPhaseParamNo > 0) { - fout << "phase par" << fFourier.fPhaseParamNo << endl; - } else { - if (fFourier.fPhase != -999.0) { - fout << "phase " << fFourier.fPhase << endl; + if (fFourier.fPhaseParamNo.size() > 0) { + TString phaseParamStr = BeautifyFourierPhaseParameterString(); + fout << "phase " << phaseParamStr << endl; + } else if (fFourier.fPhase.size() > 0) { + fout << "phase "; + for (UInt_t i=0; i *co fout << endl; } - // phase - if (fFourier.fPhaseParamNo > 0) { - fout << "phase par" << fFourier.fPhaseParamNo << endl; - } else if (fFourier.fPhase != -999.0) { - fout << "phase " << fFourier.fPhase << endl; + // phase + if (fFourier.fPhaseParamNo.size() > 0) { + TString phaseParamStr = BeautifyFourierPhaseParameterString(); + fout << "phase " << phaseParamStr << endl; + } else if (fFourier.fPhase.size() > 0) { + fout << "phase "; + for (UInt_t i=0; iRemoves a potentially present comment from str and returns the truncated + * string in truncStr. A comment starts with '#' + * + * @param str original string which might contain a comment + * @param truncStr string from which the comment has been removed + */ +void PMsrHandler::RemoveComment(const TString &str, TString &truncStr) +{ + truncStr = str; + Ssiz_t idx = str.First('#'); // find the index of the comment character + + // truncate string if comment is found + if (idx > 0) { + truncStr.Resize(idx-1); + } +} + +//-------------------------------------------------------------------------- +// ParseFourierPhaseValueVector (private) +//-------------------------------------------------------------------------- +/** + *

examines if str has the form 'phase val0 [sep val1 ... sep valN]'. + * If this form is found, fill in val0 ... valN to fFourier.fPhase + * vector. + * + * @param fourier msr-file Fourier structure + * @param str string to be analyzed + * @param error flag needed to propagate a fatal error + * + * @return true if a phase value form is found, otherwise return false + */ +Bool_t PMsrHandler::ParseFourierPhaseValueVector(PMsrFourierStructure &fourier, const TString &str, Bool_t &error) +{ + Bool_t result = true; + + TObjArray *tok = str.Tokenize(" ,;\t"); + if (tok == 0) { + cerr << endl << ">> PMsrHandler::ParseFourierPhaseValueVector: **ERROR** couldn't tokenize Fourier phase line." << endl << endl; + return false; + } + + // make sure there are enough tokens + if (tok->GetEntries() < 2) { + error = true; + return false; + } + + // convert all acceptable tokens + TObjString *ostr=0; + TString sstr(""); + for (Int_t i=1; iGetEntries(); i++) { + ostr = dynamic_cast(tok->At(i)); + sstr = ostr->GetString(); + if (sstr.IsFloat()) { + fourier.fPhase.push_back(sstr.Atof()); + } else { + result = false; + if (i>1) { // make sure that no 'phase val, parX' mixture is present + cerr << endl << ">> PMsrHandler::ParseFourierPhaseValueVector: **ERROR** in Fourier phase line."; + cerr << endl << ">> Attempt to mix val, parX? This is currently not supported." << endl << endl; + error = true; + } + break; + } + } + + // clean up + if (tok) { + delete tok; + } + + return result; +} + +//-------------------------------------------------------------------------- +// ParseFourierPhaseParVector (private) +//-------------------------------------------------------------------------- +/** + *

examines if str has the form 'phase parX0 [sep parX1 ... sep parXN]'. + * If this form is found, fill in parX0 ... parXN to fFourier.fPhaseParamNo + * and furthermore fill fFourier.fPhase accordingly. + * + * @param fourier msr-file Fourier structure + * @param str string to be analyzed + * @param error flag needed to propagate a fatal error + * + * @return true if a phase parameter form is found, otherwise return false + */ +Bool_t PMsrHandler::ParseFourierPhaseParVector(PMsrFourierStructure &fourier, const TString &str, Bool_t &error) +{ + Bool_t result = true; + + TObjArray *tok = str.Tokenize(" ,;\t"); + if (tok == 0) { + cerr << endl << ">> PMsrHandler::ParseFourierPhaseParVector: **ERROR** couldn't tokenize Fourier phase line." << endl << endl; + return false; + } + + // make sure there are enough tokens + if (tok->GetEntries() < 2) { + error = true; + return false; + } + + // check that all tokens start with par + TString sstr; + for (Int_t i=1; iGetEntries(); i++) { + TObjString *ostr = dynamic_cast(tok->At(i)); + sstr = ostr->GetString(); + if (!sstr.BeginsWith("par")) { + cerr << ">> PMsrHandler::ParseFourierPhaseParVector: **ERROR** found unhandable token '" << sstr << "'" << endl; + error = true; + result = false; + break; + } + + // rule out par(X, offset, #Param) syntax + if (sstr.BeginsWith("par(")) { + result = false; + break; + } + } + + // check that token has the form parX, where X is an int + if (result != false) { + for (Int_t i=1; iGetEntries(); i++) { + TObjString *ostr = dynamic_cast(tok->At(i)); + sstr = ostr->GetString(); + sstr.Remove(0, 3); // remove 'par' part. Rest should be an integer + if (sstr.IsDigit()) { + fourier.fPhaseParamNo.push_back(sstr.Atoi()); + } else { + cerr << ">> PMsrHandler::ParseFourierPhaseParVector: **ERROR** found token '" << ostr->GetString() << "' which is not parX with X an integer." << endl; + fourier.fPhaseParamNo.clear(); + error = true; + break; + } + } + } + + if (fourier.fPhaseParamNo.size() == tok->GetEntries()-1) { // everything as expected + result = true; + } else { + result = false; + } + + // clean up + if (tok) { + delete tok; + } + + return result; +} + +//-------------------------------------------------------------------------- +// ParseFourierPhaseParIterVector (private) +//-------------------------------------------------------------------------- +/** + *

examines if str has the form 'phase par(X0, offset, #params)'. + * If this form is found, fill in parX0 ... parXN to fFourier.fPhaseParamNo + * and furthermore fill fFourier.fPhase accordingly. + * + * @param fourier msr-file Fourier structure + * @param str string to be analyzed + * @param error flag needed to propagate a fatal error + * + * @return true if a phase parameter iterator form is found, otherwise return false + */ +Bool_t PMsrHandler::ParseFourierPhaseParIterVector(PMsrFourierStructure &fourier, const TString &str, Bool_t &error) +{ + TString wstr = str; + + // remove 'phase' from string + wstr.Remove(0, 5); + wstr = wstr.Strip(TString::kLeading, ' '); + + // remove 'par(' from string if present, otherwise and error is issued + if (!wstr.BeginsWith("par(")) { + cout << ">> PMsrHandler::ParseFourierPhaseParIterVector: **ERROR** token should start with 'par(', found: '" << wstr << "' -> ERROR" << endl; + error = true; + return false; + } + wstr.Remove(0, 4); + + // remove trailing white spaces + wstr = wstr.Strip(TString::kTrailing, ' '); + + // remove last ')' + Ssiz_t idx=wstr.Last(')'); + wstr.Remove(idx, wstr.Length()-idx); + + // tokenize rest which should have the form 'X0, offset, #Param' + TObjArray *tok = wstr.Tokenize(",;"); + if (tok == 0) { + cerr << ">> PMsrHandler::ParseFourierPhaseParIterVector: **ERROR** tokenize failed." << endl; + error = true; + return false; + } + + // check for proper number of expected elements + if (tok->GetEntries() != 3) { + cerr << ">> PMsrHandler::ParseFourierPhaseParIterVector: **ERROR** wrong syntax for the expected par(X0, offset, #param)." << endl; + error = true; + delete tok; + return false; + } + + Int_t x0, offset, noParam; + + // get X0 + TObjString *ostr = dynamic_cast(tok->At(0)); + wstr = ostr->GetString(); + if (wstr.IsDigit()) { + x0 = wstr.Atoi(); + } else { + cerr << ">> PMsrHandler::ParseFourierPhaseParIterVector: **ERROR** X0='" << wstr << "' is not an integer." << endl; + error = true; + delete tok; + return false; + } + + // get offset + ostr = dynamic_cast(tok->At(1)); + wstr = ostr->GetString(); + if (wstr.IsDigit()) { + offset = wstr.Atoi(); + } else { + cerr << ">> PMsrHandler::ParseFourierPhaseParIterVector: **ERROR** offset='" << wstr << "' is not an integer." << endl; + error = true; + delete tok; + return false; + } + + // get noParam + ostr = dynamic_cast(tok->At(2)); + wstr = ostr->GetString(); + if (wstr.IsDigit()) { + noParam = wstr.Atoi(); + } else { + cerr << ">> PMsrHandler::ParseFourierPhaseParIterVector: **ERROR** #Param='" << wstr << "' is not an integer." << endl; + error = true; + delete tok; + return false; + } + + for (Int_t i=0; ifLine.BeginsWith("phase", TString::kIgnoreCase)) { // phase - if (tokens->GetEntries() < 2) { // phase value is missing + if (tokens->GetEntries() < 2) { // phase value(s)/par(s) is(are) missing error = true; continue; } else { - ostr = dynamic_cast(tokens->At(1)); - str = ostr->GetString(); - if (str.BeginsWith("par", TString::kIgnoreCase)) { // parameter value - if (fFourierOnly) { - cerr << endl << ">> PMsrHandler::HandleFourierEntry: **WARNING** Found phase parameter for Fourier only."; - cerr << endl << ">> This is currently not supported. Will set the phase to 0." << endl; - fourier.fPhase = 0.0; - } else { - Int_t no = 0; - if (FilterNumber(str, "par", 0, no)) { - // check that the parameter is in range - if ((Int_t)fParam.size() < no) { - error = true; - continue; - } - // keep the parameter number - fourier.fPhaseParamNo = no; - // get parameter value - fourier.fPhase = fParam[no-1].fValue; - } else { + // allowed phase parameter patterns: + // (i) phase val [sep val sep val ...] [# comment], val=double, sep=' ,;\t' + // (ii) phase parX0 [sep parX1 sep parX2 ...] [# comment], val=double, sep=' ,;\t' + // (iii) phase par(X0 sep1 offset sep1 #param) [# comment], sep1= ',;' + + // remove potential comment + TString wstr(""); + RemoveComment(iter->fLine, wstr); + + // check for 'phase val ...' + Bool_t result = ParseFourierPhaseValueVector(fourier, wstr, error); + if (error) + continue; + + // check for 'phase parX0 ...' if not already val are found + if (!result) { + result = ParseFourierPhaseParVector(fourier, wstr, error); + if (error) + continue; + } + + // check for 'phase par(X0, offset, #param)' if not already covered by the previous ones + if (!result) { + result = ParseFourierPhaseParIterVector(fourier, wstr, error); + } + + if (!result || error) { + continue; + } + + // if parameter vector is given: check that all parameters are within range + if (fourier.fPhaseParamNo.size() > 0) { + for (UInt_t i=0; i fParam.size()) { + cerr << ">> PMsrHandler::HandleFourierEntry: found Fourier parameter entry par" << fourier.fPhaseParamNo[i] << " > #Param = " << fParam.size() << endl; error = true; + --iter; continue; } } - } else if (str.IsFloat()) { // phase value - fourier.fPhase = str.Atof(); - } else { - error = true; - continue; + } + + // if parameter vector is given -> fill corresponding phase values + if (fourier.fPhaseParamNo.size() > 0) { + fourier.fPhase.clear(); + UInt_t idx; + for (UInt_t i=0; ifLine.BeginsWith("range_for_phase_correction", TString::kIgnoreCase)) { @@ -4144,7 +4441,11 @@ Bool_t PMsrHandler::HandleFourierEntry(PMsrLines &lines) cerr << endl << ">> [dc-corrected true | false]"; cerr << endl << ">> [apodization none | weak | medium | strong]"; cerr << endl << ">> [plot real | imag | real_and_imag | power | phase | phase_opt_real]"; - cerr << endl << ">> [phase value]"; + cerr << endl << ">> [phase valList | parList | parIterList [# comment]]"; + cerr << endl << ">> valList : val [sep val ... sep val]. sep=' ,;\\t'"; + cerr << endl << ">> parList : parX0 [sep parX1 ... sep parX1]"; + cerr << endl << ">> parIterList : par(X0,offset,#param), with X0=first parameter number"; + cerr << endl << ">> offset=parameter offset, #param=number of phase parameters."; cerr << endl << ">> [range_for_phase_correction min max | all]"; cerr << endl << ">> [range min max]"; cerr << endl; @@ -6561,6 +6862,56 @@ void PMsrHandler::MakeDetectorGroupingString(TString str, PIntVector &group, TSt } while (iReturns the Fourier phase string if the phase is either of type + * phase parX0 sep ... sep parXn where sep = ',' or + * phase par(X0, offset, #param) + * + * @return Fourier phase parameter string if phase parameter(s) is(are) given, "??" otherwise + */ +TString PMsrHandler::BeautifyFourierPhaseParameterString() +{ + TString str("??"); + + if (fFourier.fPhaseParamNo.size() == 0) + return str; + + if (fFourier.fPhaseParamNo.size() == 1) { + str = TString::Format("par%d", fFourier.fPhaseParamNo[0]); + } else if (fFourier.fPhaseParamNo.size() == 2) { + str = TString::Format("par%d, par%d", fFourier.fPhaseParamNo[0], fFourier.fPhaseParamNo[1]); + } else { + Bool_t phaseIter = true; + + // first check if fPhaseParamNo vector can be compacted into par(X0, offset, #param) form + Int_t offset = fFourier.fPhaseParamNo[1] - fFourier.fPhaseParamNo[0]; + for (Int_t i=2; iClear(); delete fMultiGraphLegend; @@ -396,9 +389,8 @@ void PMusrCanvas::SetMsrHandler(PMsrHandler *msrHandler) if (fMsrHandler->GetMsrFourierList()->fPlotTag != FOURIER_PLOT_NOT_GIVEN) { fFourier.fPlotTag = fMsrHandler->GetMsrFourierList()->fPlotTag; } - if (fMsrHandler->GetMsrFourierList()->fPhase != -999.0) { - fFourier.fPhase = fMsrHandler->GetMsrFourierList()->fPhase; - } + fFourier.fPhase = fMsrHandler->GetMsrFourierList()->fPhase; + if ((fMsrHandler->GetMsrFourierList()->fRangeForPhaseCorrection[0] != -1.0) && (fMsrHandler->GetMsrFourierList()->fRangeForPhaseCorrection[1] != -1.0)) { fFourier.fRangeForPhaseCorrection[0] = fMsrHandler->GetMsrFourierList()->fRangeForPhaseCorrection[0]; @@ -2314,7 +2306,8 @@ void PMusrCanvas::InitFourier() fFourier.fFourierPower = 0; // no zero padding fFourier.fApodization = FOURIER_APOD_NONE; // no apodization fFourier.fPlotTag = FOURIER_PLOT_REAL_AND_IMAG; // initial plot tag, plot real and imaginary part - fFourier.fPhase = 0.0; // fourier phase 0° + fFourier.fPhaseParamNo.clear(); + fFourier.fPhase.clear(); for (UInt_t i=0; i<2; i++) { fFourier.fRangeForPhaseCorrection[i] = -1.0; // frequency range for phase correction, default: {-1, -1} = NOT GIVEN fFourier.fPlotRange[i] = -1.0; // fourier plot range, default: {-1, -1} = NOT GIVEN @@ -3471,13 +3464,20 @@ void PMusrCanvas::HandleFourier() } // apply global phase if present - if (fFourier.fPhase != 0.0) { - const double cp = TMath::Cos(fFourier.fPhase/180.0*TMath::Pi()); - const double sp = TMath::Sin(fFourier.fPhase/180.0*TMath::Pi()); + if (fFourier.fPhase.size() != 0.0) { + double cp; + double sp; fCurrentFourierPhase = fFourier.fPhase; for (UInt_t i=0; iGetNbinsX(); j++) { // loop over a fourier data set // calculate new fourier data set value @@ -3574,16 +3574,23 @@ void PMusrCanvas::HandleDifferenceFourier() fData[i].diffFourierTag = 1; // d-f } - // apply global phase - if (fFourier.fPhase != 0.0) { + // apply phase + if (fFourier.fPhase.size() != 0.0) { double re, im; - const double cp = TMath::Cos(fFourier.fPhase/180.0*TMath::Pi()); - const double sp = TMath::Sin(fFourier.fPhase/180.0*TMath::Pi()); + double cp; + double sp; fCurrentFourierPhase = fFourier.fPhase; for (UInt_t i=0; iGetNbinsX(); j++) { // loop over a fourier data set // calculate new fourier data set value re = fData[i].diffFourierRe->GetBinContent(j) * cp + fData[i].diffFourierIm->GetBinContent(j) * sp; @@ -6046,7 +6053,10 @@ void PMusrCanvas::PlotFourierPhaseValue(Bool_t unzoom) // plot Fourier phase str = TString("phase = "); - str += fCurrentFourierPhase; + str += fCurrentFourierPhase[0]; + if (fFourier.fPhase.size() > 1) { // if more than one phase is present, do NOT plot phase info + str = TString(""); + } x = 0.7; y = 0.85; fCurrentFourierPhaseText = new TLatex(); @@ -6298,7 +6308,8 @@ void PMusrCanvas::IncrementFourierPhase() const double cp = TMath::Cos(fFourier.fPhaseIncrement/180.0*TMath::Pi()); const double sp = TMath::Sin(fFourier.fPhaseIncrement/180.0*TMath::Pi()); - fCurrentFourierPhase += fFourier.fPhaseIncrement; + for (UInt_t i=0; i