From e0c893503f00fbb881a137ec46c2b0a2d5219ddf Mon Sep 17 00:00:00 2001 From: nemu Date: Wed, 23 Dec 2009 13:43:23 +0000 Subject: [PATCH] first implementation of a RRF. So far, for single histo only. --- src/classes/PMsrHandler.cpp | 283 +++++++----------- src/classes/PMusr.cpp | 41 ++- src/classes/PMusrCanvas.cpp | 75 +++++ src/classes/PRunBase.cpp | 85 ++++++ src/classes/PRunSingleHisto.cpp | 136 +++++++-- src/include/PMusr.h | 33 +- src/include/PMusrCanvas.h | 2 + src/include/PRunBase.h | 9 +- src/musrt0.cpp | 2 + .../t0NotEqFirstGoodData.C | 18 +- 10 files changed, 444 insertions(+), 240 deletions(-) diff --git a/src/classes/PMsrHandler.cpp b/src/classes/PMsrHandler.cpp index d9fa8db5..17ab8782 100644 --- a/src/classes/PMsrHandler.cpp +++ b/src/classes/PMsrHandler.cpp @@ -599,16 +599,6 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t messages) default: break; } - } else if (sstr.BeginsWith("rrffrequency")) { - fout.width(16); - fout << left << "rrffrequency"; - fout.precision(prec); - fout << fRuns[runNo].GetRRFFreq() << endl; - } else if (sstr.BeginsWith("rrfpacking")) { - fout.width(16); - fout << left << "rrfpacking"; - fout.precision(prec); - fout << fRuns[runNo].GetRRFPacking() << endl; } else if (sstr.BeginsWith("alpha ")) { fout.width(16); fout << left << "alpha"; @@ -617,14 +607,6 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t messages) fout.width(16); fout << left << "beta"; fout << fRuns[runNo].GetBetaParamNo() << endl; - } else if (sstr.BeginsWith("alpha2")) { - fout.width(16); - fout << left << "alpha2"; - fout << fRuns[runNo].GetAlpha2ParamNo() << endl; - } else if (sstr.BeginsWith("beta2")) { - fout.width(16); - fout << left << "beta2"; - fout << fRuns[runNo].GetBeta2ParamNo() << endl; } else if (sstr.BeginsWith("norm")) { fout.width(16); fout << left << "norm"; @@ -670,14 +652,6 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t messages) fout.width(16); fout << left << "backward"; fout << fRuns[runNo].GetBackwardHistoNo() << endl; - } else if (sstr.BeginsWith("right")) { - fout.width(16); - fout << left << "right"; - fout << fRuns[runNo].GetRightHistoNo() << endl; - } else if (sstr.BeginsWith("left")) { - fout.width(16); - fout << left << "left"; - fout << fRuns[runNo].GetLeftHistoNo() << endl; } else if (sstr.BeginsWith("backgr.fix")) { fout.width(15); fout << left << "backgr.fix"; @@ -1731,28 +1705,6 @@ Bool_t PMsrHandler::HandleRunEntry(PMsrLines &lines) } } - // rphase ------------------------------------------------ - if (iter->fLine.BeginsWith("rphase", TString::kIgnoreCase)) { - - runLinePresent = false; // this is needed to make sure that a run line is present before and ADDRUN is following - - if (tokens->GetEntries() < 2) { - error = true; - } else { - ostr = dynamic_cast(tokens->At(1)); - str = ostr->GetString(); - if (str.IsDigit()) { - dval = str.Atoi(); - if (dval > 0) - param.SetPhaseParamNo(dval); - else - error = true; - } else { - error = true; - } - } - } - // lifetime ------------------------------------------------ if (iter->fLine.BeginsWith("lifetime ", TString::kIgnoreCase)) { @@ -1987,133 +1939,6 @@ Bool_t PMsrHandler::HandleRunEntry(PMsrLines &lines) } } - // rrffrequency -------------------------------------------------- - if (iter->fLine.BeginsWith("rrffrequency", TString::kIgnoreCase)) { - - runLinePresent = false; // this is needed to make sure that a run line is present before and ADDRUN is following - - if (tokens->GetEntries() != 2) { - error = true; - } else { - ostr = dynamic_cast(tokens->At(1)); - str = ostr->GetString(); - if (str.IsFloat()) - param.SetRRFFreq(str.Atof()); - else - error = true; - } - } - - // rrfpacking -------------------------------------------------- - if (iter->fLine.BeginsWith("rrfpacking", TString::kIgnoreCase)) { - - runLinePresent = false; // this is needed to make sure that a run line is present before and ADDRUN is following - - if (tokens->GetEntries() != 2) { - error = true; - } else { - ostr = dynamic_cast(tokens->At(1)); - str = ostr->GetString(); - if (str.IsDigit()) { - dval = str.Atoi(); - if (dval > 0) - param.SetRRFPacking(dval); - else - error = true; - } else { - error = true; - } - } - } - - // alpha2 -------------------------------------------------- - if (iter->fLine.BeginsWith("alpha2", TString::kIgnoreCase)) { - - runLinePresent = false; // this is needed to make sure that a run line is present before and ADDRUN is following - - if (tokens->GetEntries() != 2) { - error = true; - } else { - ostr = dynamic_cast(tokens->At(1)); - str = ostr->GetString(); - if (str.IsDigit()) { - dval = str.Atoi(); - if (dval > 0) - param.SetAlpha2ParamNo(dval); - else - error = true; - } else { - error = true; - } - } - } - - // beta2 -------------------------------------------------- - if (iter->fLine.BeginsWith("beta2", TString::kIgnoreCase)) { - - runLinePresent = false; // this is needed to make sure that a run line is present before and ADDRUN is following - - if (tokens->GetEntries() != 2) { - error = true; - } else { - ostr = dynamic_cast(tokens->At(1)); - str = ostr->GetString(); - if (str.IsDigit()) { - dval = str.Atoi(); - if (dval > 0) - param.SetBeta2ParamNo(dval); - else - error = true; - } else { - error = true; - } - } - } - - // right -------------------------------------------------- - if (iter->fLine.BeginsWith("right", TString::kIgnoreCase)) { - - runLinePresent = false; // this is needed to make sure that a run line is present before and ADDRUN is following - - if (tokens->GetEntries() != 2) { - error = true; - } else { - ostr = dynamic_cast(tokens->At(1)); - str = ostr->GetString(); - if (str.IsDigit()) { - dval = str.Atoi(); - if (dval > 0) - param.SetRightHistoNo(dval); - else - error = true; - } else { - error = true; - } - } - } - - // left -------------------------------------------------- - if (iter->fLine.BeginsWith("left", TString::kIgnoreCase)) { - - runLinePresent = false; // this is needed to make sure that a run line is present before and ADDRUN is following - - if (tokens->GetEntries() != 2) { - error = true; - } else { - ostr = dynamic_cast(tokens->At(1)); - str = ostr->GetString(); - if (str.IsDigit()) { - dval = str.Atoi(); - if (dval > 0) - param.SetLeftHistoNo(str.Atoi()); - else - error = true; - } else { - error = true; - } - } - } - // xy-data ----------------------------------------------- if (iter->fLine.BeginsWith("xy-data", TString::kIgnoreCase)) { @@ -2536,6 +2361,10 @@ Bool_t PMsrHandler::HandlePlotEntry(PMsrLines &lines) param.fLogX = false; // i.e. if not overwritten use linear x-axis param.fLogY = false; // i.e. if not overwritten use linear y-axis param.fViewPacking = -1; // i.e. if not overwritten use the packing of the run blocks + param.fRRFPacking = 0; // i.e. if not overwritten it will not be a valid RRF + param.fRRFFreq = 0.0; // i.e. no RRF whished + param.fRRFUnit = RRF_UNIT_MHz; + param.fRRFPhase = 0.0; // find next plot if any is present iter2 = iter1; @@ -2777,7 +2606,7 @@ Bool_t PMsrHandler::HandlePlotEntry(PMsrLines &lines) } else if (iter1->fLine.Contains("view_packing", TString::kIgnoreCase)) { tokens = iter1->fLine.Tokenize(" \t"); if (!tokens) { - cerr << endl << ">> PMsrHandler::HandlePlotEntry: **SEVERE ERROR** Couldn't tokenize PLOT in line " << iter1->fLineNo; + cerr << endl << ">> PMsrHandler::HandlePlotEntry: **SEVERE ERROR** Couldn't tokenize view_packing in line " << iter1->fLineNo; cerr << endl << endl; return false; } @@ -2797,6 +2626,97 @@ Bool_t PMsrHandler::HandlePlotEntry(PMsrLines &lines) } } + // clean up + if (tokens) { + delete tokens; + tokens = 0; + } + } else if (iter1->fLine.Contains("rrf_freq", TString::kIgnoreCase)) { + // expected entry: rrf_freq value unit + // allowed units: kHz, MHz, Mc/s, G, T + tokens = iter1->fLine.Tokenize(" \t"); + if (!tokens) { + cerr << endl << ">> PMsrHandler::HandlePlotEntry: **SEVERE ERROR** Couldn't tokenize rrf_freq in line " << iter1->fLineNo; + cerr << endl << endl; + return false; + } + if (tokens->GetEntries() != 3) { + error = true; + } else { + // get rrf frequency + ostr = dynamic_cast(tokens->At(1)); + str = ostr->GetString(); + if (str.IsFloat()) { + param.fRRFFreq = str.Atof(); + } else { + error = true; + } + // get unit + ostr = dynamic_cast(tokens->At(2)); + str = ostr->GetString(); + if (str.Contains("kHz", TString::kIgnoreCase)) + param.fRRFUnit = RRF_UNIT_kHz; + else if (str.Contains("MHz", TString::kIgnoreCase)) + param.fRRFUnit = RRF_UNIT_MHz; + else if (str.Contains("Mc/s", TString::kIgnoreCase)) + param.fRRFUnit = RRF_UNIT_Mcs; + else if (str.Contains("G", TString::kIgnoreCase)) + param.fRRFUnit = RRF_UNIT_G; + else if (str.Contains("T", TString::kIgnoreCase)) + param.fRRFUnit = RRF_UNIT_T; + else + error = true; + } + // clean up + if (tokens) { + delete tokens; + tokens = 0; + } + } else if (iter1->fLine.Contains("rrf_phase", TString::kIgnoreCase)) { + // expected entry: rrf_phase value. value given in units of degree + tokens = iter1->fLine.Tokenize(" \t"); + if (!tokens) { + cerr << endl << ">> PMsrHandler::HandlePlotEntry: **SEVERE ERROR** Couldn't tokenize rrf_phase in line " << iter1->fLineNo; + cerr << endl << endl; + return false; + } + if (tokens->GetEntries() != 2) { + error = true; + } else { + // get rrf phase + ostr = dynamic_cast(tokens->At(1)); + str = ostr->GetString(); + if (str.IsFloat()) { + param.fRRFPhase = str.Atof(); + } else { + error = true; + } + } + // clean up + if (tokens) { + delete tokens; + tokens = 0; + } + } else if (iter1->fLine.Contains("rrf_packing", TString::kIgnoreCase)) { + // expected entry: rrf_phase value. value given in units of degree + tokens = iter1->fLine.Tokenize(" \t"); + if (!tokens) { + cerr << endl << ">> PMsrHandler::HandlePlotEntry: **SEVERE ERROR** Couldn't tokenize rrf_packing in line " << iter1->fLineNo; + cerr << endl << endl; + return false; + } + if (tokens->GetEntries() != 2) { + error = true; + } else { + // get rrf packing + ostr = dynamic_cast(tokens->At(1)); + str = ostr->GetString(); + if (str.IsDigit()) { + param.fRRFPacking = str.Atoi(); + } else { + error = true; + } + } // clean up if (tokens) { delete tokens; @@ -2835,6 +2755,17 @@ Bool_t PMsrHandler::HandlePlotEntry(PMsrLines &lines) } } } + + // check RRF entries + if (param.fRRFFreq != 0.0) { + if (param.fRRFPacking == 0) { + cerr << endl << ">> PMsrHandler::HandlePlotEntry(): **ERROR** found RRF frequency but no required RRF packing."; + cerr << endl << ">> Will ignore the RRF option."; + cerr << endl; + param.fRRFFreq = 0.0; + } + } + fPlots.push_back(param); } } diff --git a/src/classes/PMusr.cpp b/src/classes/PMusr.cpp index f3dd02de..813cbd15 100644 --- a/src/classes/PMusr.cpp +++ b/src/classes/PMusr.cpp @@ -74,6 +74,35 @@ PRunData::~PRunData() fTheory.clear(); } +//-------------------------------------------------------------------------- +// SetTheoryValue +//-------------------------------------------------------------------------- +/** + *

Sets a value of the theory vector + * + * \param i index of the theory vector + */ +void PRunData::SetTheoryValue(UInt_t i, Double_t dval) +{ + if (i > fTheory.size()) + fTheory.resize(i+1); + + fTheory[i] = dval; +} + +//-------------------------------------------------------------------------- +// ReplaceTheory +//-------------------------------------------------------------------------- +/** + *

Replaces the theory vector. + * + * \param theo vector which is replacing the current theory vector + */ +void PRunData::ReplaceTheory(const PDoubleVector &theo) +{ + fTheory = theo; +} + //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ // implementation PNonMusrRawRunData //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ @@ -504,12 +533,6 @@ PMsrRunBlock::PMsrRunBlock() fFitRange[0] = PMUSR_UNDEFINED; // undefined start fit range fFitRange[1] = PMUSR_UNDEFINED; // undefined end fit range fPacking = -1; // undefined packing - fRRFFreq = PMUSR_UNDEFINED; // undefined RRF frequency - fRRFPacking = -1; // undefined RRF packing - fAlpha2ParamNo = -1; // undefined alpha2 parameter number - fBeta2ParamNo = -1; // undefined beta2 parameter number - fRightHistoNo = -1; // undefined right histogram number - fLeftHistoNo = -1; // undefined left histogram number fXYDataIndex[0] = -1; // undefined x data index (NonMusr) fXYDataIndex[1] = -1; // undefined y data index (NonMusr) fXYDataLabel[0] = TString(""); // undefined x data label (NonMusr) @@ -556,12 +579,6 @@ void PMsrRunBlock::CleanUp() fFitRange[0] = PMUSR_UNDEFINED; // undefined start fit range fFitRange[1] = PMUSR_UNDEFINED; // undefined end fit range fPacking = -1; // undefined packing - fRRFFreq = PMUSR_UNDEFINED; // undefined RRF frequency - fRRFPacking = -1; // undefined RRF packing - fAlpha2ParamNo = -1; // undefined alpha2 parameter number - fBeta2ParamNo = -1; // undefined beta2 parameter number - fRightHistoNo = -1; // undefined right histogram number - fLeftHistoNo = -1; // undefined left histogram number fXYDataIndex[0] = -1; // undefined x data index (NonMusr) fXYDataIndex[1] = -1; // undefined y data index (NonMusr) fXYDataLabel[0] = TString("??"); // undefined x data label (NonMusr) diff --git a/src/classes/PMusrCanvas.cpp b/src/classes/PMusrCanvas.cpp index 5ac07a7b..9304589a 100644 --- a/src/classes/PMusrCanvas.cpp +++ b/src/classes/PMusrCanvas.cpp @@ -81,6 +81,9 @@ PMusrCanvas::PMusrCanvas() fCurrentFourierPhase = fFourier.fPhaseIncrement; fCurrentFourierPhaseText = 0; + + fRRFText = 0; + fRRFLatexText = 0; } //-------------------------------------------------------------------------- @@ -105,6 +108,9 @@ PMusrCanvas::PMusrCanvas(const Int_t number, const Char_t* title, fCurrentFourierPhase = 0.0; fCurrentFourierPhaseText = 0; + + fRRFText = 0; + fRRFLatexText = 0; } //-------------------------------------------------------------------------- @@ -132,6 +138,9 @@ PMusrCanvas::PMusrCanvas(const Int_t number, const Char_t* title, fCurrentFourierPhase = 0.0; fCurrentFourierPhaseText = 0; + + fRRFText = 0; + fRRFLatexText = 0; } //-------------------------------------------------------------------------- @@ -144,6 +153,18 @@ PMusrCanvas::~PMusrCanvas() { //cout << "~PMusrCanvas() called. fMainCanvas name=" << fMainCanvas->GetName() << endl; // cleanup + if (fCurrentFourierPhaseText) { + delete fCurrentFourierPhaseText; + fCurrentFourierPhaseText = 0; + } + if (fRRFLatexText) { + delete fRRFLatexText; + fRRFLatexText = 0; + } + if (fRRFText) { + delete fRRFText; + fRRFText = 0; + } if (fStyle) { delete fStyle; fStyle = 0; @@ -244,6 +265,40 @@ void PMusrCanvas::SetMsrHandler(PMsrHandler *msrHandler) fFourier.fPlotRange[1] = fMsrHandler->GetMsrFourierList().fPlotRange[1]; } } + + // check if RRF data are present + if ((fMsrHandler->GetMsrPlotList()->at(0).fRRFPacking > 0) && + (fMsrHandler->GetMsrPlotList()->at(0).fRRFFreq != 0.0)) { + fRRFLatexText = new TLatex(); + fRRFLatexText->SetNDC(kTRUE); + fRRFLatexText->SetTextFont(62); + fRRFLatexText->SetTextSize(0.03); + + fRRFText = new TString("RRF: "); + if (fMsrHandler->GetMsrPlotList()->at(0).fRRFUnit == RRF_UNIT_kHz) { + *fRRFText += TString("#nu_{RRF} = "); + *fRRFText += fMsrHandler->GetMsrPlotList()->at(0).fRRFFreq; + *fRRFText += TString(" (kHz)"); + } else if (fMsrHandler->GetMsrPlotList()->at(0).fRRFUnit == RRF_UNIT_MHz) { + *fRRFText += TString("#nu_{RRF} = "); + *fRRFText += fMsrHandler->GetMsrPlotList()->at(0).fRRFFreq; + *fRRFText += TString(" (MHz)"); + } else if (fMsrHandler->GetMsrPlotList()->at(0).fRRFUnit == RRF_UNIT_Mcs) { + *fRRFText += TString("#omega_{RRF} = "); + *fRRFText += fMsrHandler->GetMsrPlotList()->at(0).fRRFFreq; + *fRRFText += TString(" (Mc/s)"); + } else if (fMsrHandler->GetMsrPlotList()->at(0).fRRFUnit == RRF_UNIT_G) { + *fRRFText += TString("B_{RRF} = "); + *fRRFText += fMsrHandler->GetMsrPlotList()->at(0).fRRFFreq; + *fRRFText += TString(" (G)"); + } else if (fMsrHandler->GetMsrPlotList()->at(0).fRRFUnit == RRF_UNIT_T) { + *fRRFText += TString("B_{RRF} = "); + *fRRFText += fMsrHandler->GetMsrPlotList()->at(0).fRRFFreq; + *fRRFText += TString(" (T)"); + } + *fRRFText += TString(", RRF packing = "); + *fRRFText += fMsrHandler->GetMsrPlotList()->at(0).fRRFPacking; + } } //-------------------------------------------------------------------------- @@ -2531,6 +2586,11 @@ void PMusrCanvas::PlotData() fData[i].theory->Draw("lsame"); } } + + // check if RRF and if yes show a label + if ((fRRFText != 0) && (fRRFLatexText != 0)) { + fRRFLatexText->DrawLatex(0.1, 0.92, fRRFText->Data()); + } } else { // fPlotType == MSR_PLOT_NO_MUSR // ugly workaround since multigraphs axis are not going away when switching TMultiGraphs delete fDataTheoryPad; @@ -2649,6 +2709,11 @@ void PMusrCanvas::PlotDifference() for (UInt_t i=0; iDraw("pesame"); } + + // check if RRF and if yes show a label + if ((fRRFText != 0) && (fRRFLatexText != 0)) { + fRRFLatexText->DrawLatex(0.1, 0.92, fRRFText->Data()); + } } else { // fPlotType == MSR_PLOT_NON_MUSR // ugly workaround since multigraphs axis are not going away when switching TMultiGraphs delete fDataTheoryPad; @@ -2980,6 +3045,11 @@ void PMusrCanvas::PlotFourier() break; } + // check if RRF and if yes show a label + if ((fRRFText != 0) && (fRRFLatexText != 0)) { + fRRFLatexText->DrawLatex(0.1, 0.92, fRRFText->Data()); + } + fDataTheoryPad->Update(); fMainCanvas->cd(); @@ -3249,6 +3319,11 @@ void PMusrCanvas::PlotFourierDifference() break; } + // check if RRF and if yes show a label + if ((fRRFText != 0) && (fRRFLatexText != 0)) { + fRRFLatexText->DrawLatex(0.1, 0.92, fRRFText->Data()); + } + fDataTheoryPad->Update(); fMainCanvas->cd(); diff --git a/src/classes/PRunBase.cpp b/src/classes/PRunBase.cpp index 89a7ed21..e8e9ac3a 100644 --- a/src/classes/PRunBase.cpp +++ b/src/classes/PRunBase.cpp @@ -145,3 +145,88 @@ void PRunBase::CleanUp() fTheory = 0; } } + +//-------------------------------------------------------------------------- +// CalculateKaiserFilterCoeff (protected) +//-------------------------------------------------------------------------- +/** + *

Calculates the Kaiser filter coefficients for a low pass filter with + * a cut off frequency wc. + * For details see "Zeitdiskrete Signalverarbeitung", A.V. Oppenheim, R.W. Schafer, J.R. Buck. Pearson 2004. + * + * \param wc cut off frequency + * \param A defined as \f[ A = -\log_{10}(\delta) \f], where \f[\delta\f] is the tolerance band. + * \param dw defined as \f[ \Delta\omega = \omega_{\rm S} - \omega_{\rm P} \f], where \f[ \omega_{\rm S} \f] is the + * stop band frequency, and \f[ \omega_{\rm P} \f] is the pass band frequency. + */ +void PRunBase::CalculateKaiserFilterCoeff(Double_t wc, Double_t A, Double_t dw) +{ + Double_t beta; + Double_t dt = fData.GetTheoryTimeStep(); + UInt_t m; + + // estimate beta (see reference above, p.574ff) + if (A > 50.0) { + beta = 0.1102*(A-8.7); + } else if ((A >= 21.0) && (A <= 50.0)) { + beta = 0.5842*TMath::Power(A-21.0, 0.4) + 0.07886*(A-21.0); + } else { + beta = 0.0; + } + + m = TMath::FloorNint((A-8.0)/(2.285*dw*TMath::Pi())); + // make sure m is odd + if (m % 2 == 0) + m++; + + Double_t alpha = static_cast(m)/2.0; + Double_t dval; + Double_t dsum = 0.0; + for (UInt_t i=0; i<=m; i++) { + dval = TMath::Sin(wc*(i-alpha)*dt)/(TMath::Pi()*(i-alpha)*dt); + dval *= TMath::BesselI0(beta*TMath::Sqrt(1.0-TMath::Power((i-alpha)*dt/alpha, 2.0)))/TMath::BesselI0(beta); + dsum += dval; + fKaiserFilter.push_back(dval); + } + for (UInt_t i=0; i<=m; i++) { + fKaiserFilter[i] /= dsum; + } + +/* +for (UInt_t i=0; i<=m; i++) +cout << endl << ">> " << i << ", fKaiserFilter[" << i << "]=" << fKaiserFilter[i]; +cout << endl; +*/ +} + +//-------------------------------------------------------------------------- +// FilterData (protected) +//-------------------------------------------------------------------------- +/** + *

Filters the data with a Kaiser FIR filter. + */ +void PRunBase::FilterData() +{ + PDoubleVector theoFiltered; + Double_t dval = 0.0; + const PDoubleVector *theo = fData.GetTheory(); + + for (UInt_t i=0; isize(); i++) { + for (UInt_t j=0; jat(i-j); + } + theoFiltered.push_back(dval); + dval = 0.0; + } + + fData.ReplaceTheory(theoFiltered); + + // shift time start by half the filter length + dval = fData.GetTheoryTimeStart() - 0.5*static_cast(fKaiserFilter.size())*fData.GetTheoryTimeStep(); + fData.SetTheoryTimeStart(dval); + + theoFiltered.clear(); +} diff --git a/src/classes/PRunSingleHisto.cpp b/src/classes/PRunSingleHisto.cpp index 38f89a02..5634f294 100644 --- a/src/classes/PRunSingleHisto.cpp +++ b/src/classes/PRunSingleHisto.cpp @@ -757,11 +757,15 @@ cout << endl << ">> data start time = " << fData.GetDataTimeStart(); */ Bool_t PRunSingleHisto::PrepareViewData(PRawRunData* runData, const UInt_t histoNo) { - // check if view_packing is wished + // check if view_packing is wished. This is a global option for all PLOT blocks! Int_t packing = fRunInfo->GetPacking(); if (fMsrInfo->GetMsrPlotList()->at(0).fViewPacking > 0) { packing = fMsrInfo->GetMsrPlotList()->at(0).fViewPacking; } + // check if rrf_packing is present. This is a global option for all PLOT blocks, since operated on a single set of data. + if (fMsrInfo->GetMsrPlotList()->at(0).fRRFPacking > 0) { + packing = fMsrInfo->GetMsrPlotList()->at(0).fRRFPacking; + } // transform raw histo data. This is done the following way (for details see the manual): // for the single histo fit, just the rebinned raw data are copied @@ -855,9 +859,10 @@ Bool_t PRunSingleHisto::PrepareViewData(PRawRunData* runData, const UInt_t histo bkg = par[fRunInfo->GetBkgFitParamNo()-1]; } - Double_t value = 0.0; - Double_t expval; - Double_t time; + Double_t value = 0.0; + Double_t expval = 0.0; + Double_t rrf_val = 0.0; + Double_t time = 0.0; // data start at data_start-t0 shifted by (pack-1)/2 fData.SetDataTimeStart(fTimeResolution*((Double_t)start-(Double_t)t0+(Double_t)(packing-1)/2.0)); @@ -870,21 +875,67 @@ cout << endl << "--------------------------------" << endl; */ Double_t normalizer = 1.0; - for (Int_t i=start; ins - value /= normalizer; - time = (((Double_t)i-(Double_t)(packing-1)/2.0)-t0)*fTimeResolution; - expval = TMath::Exp(+time/tau)/N0; - fData.AppendValue(-1.0+expval*(value-bkg)); -//cout << endl << ">> i=" << i << ",t0=" << t0 << ",time=" << time << ",expval=" << expval << ",value=" << value << ",bkg=" << bkg << ",expval*(value-bkg)-1=" << expval*(value-bkg)-1.0; - fData.AppendErrorValue(expval*TMath::Sqrt(value/normalizer)); -//cout << endl << ">> " << time << ", " << expval << ", " << -1.0+expval*(value-bkg) << ", " << expval*TMath::Sqrt(value/packing); - value = 0.0; + Double_t gammaRRF = 0.0, wRRF = 0.0, phaseRRF = 0.0; + if (fMsrInfo->GetMsrPlotList()->at(0).fRRFFreq == 0.0) { // normal Data representation + for (Int_t i=start; ins + value /= normalizer; + time = (((Double_t)i-(Double_t)(packing-1)/2.0)-t0)*fTimeResolution; + expval = TMath::Exp(+time/tau)/N0; + fData.AppendValue(-1.0+expval*(value-bkg)); + //cout << endl << ">> i=" << i << ",t0=" << t0 << ",time=" << time << ",expval=" << expval << ",value=" << value << ",bkg=" << bkg << ",expval*(value-bkg)-1=" << expval*(value-bkg)-1.0; + fData.AppendErrorValue(expval*TMath::Sqrt(value/normalizer)); + //cout << endl << ">> " << time << ", " << expval << ", " << -1.0+expval*(value-bkg) << ", " << expval*TMath::Sqrt(value/packing); + value = 0.0; + } + value += runData->GetDataBin(histoNo)->at(i); + } + } else { // RRF representation + // check which units shall be used + switch (fMsrInfo->GetMsrPlotList()->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 = 0.0135538817*TMath::TwoPi(); + break; + case RRF_UNIT_T: + gammaRRF = 0.0135538817*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(); + + // in order that after rebinning the fit does not need to be redone (important for plots) + // the value is normalize to per 1 nsec + normalizer = (fTimeResolution * 1.0e3); // fTimeResolution us->ns + + Double_t error = 0.0; + for (Int_t i=start; iGetDataBin(histoNo)->at(i)/normalizer-bkg))*TMath::Cos(wRRF * time + phaseRRF); + value += rrf_val; + error += runData->GetDataBin(histoNo)->at(i); } - value += runData->GetDataBin(histoNo)->at(i); } // count the number of bins to be fitted @@ -904,23 +955,58 @@ cout << endl << "--------------------------------" << endl; Double_t theoryValue; UInt_t size = runData->GetDataBin(histoNo)->size(); Double_t factor = 1.0; - if (fData.GetValue()->size() * 10 > runData->GetDataBin(histoNo)->size()) { - size = fData.GetValue()->size() * 10; - factor = (Double_t)runData->GetDataBin(histoNo)->size() / (Double_t)size; + UInt_t rebinRRF = 0; + + if (wRRF == 0) { // no RRF + // check if a finer binning for the theory is needed + if (fData.GetValue()->size() * 10 > runData->GetDataBin(histoNo)->size()) { + size = fData.GetValue()->size() * 10; + factor = (Double_t)runData->GetDataBin(histoNo)->size() / (Double_t)size; + } + fData.SetTheoryTimeStart(fData.GetDataTimeStart()); + fData.SetTheoryTimeStep(fTimeResolution*factor); + } else { // RRF + rebinRRF = static_cast((TMath::Pi()/2.0/wRRF)/fTimeResolution); // RRF time resolution / data time resolution + fData.SetTheoryTimeStart(fData.GetDataTimeStart()); + fData.SetTheoryTimeStep(TMath::Pi()/2.0/wRRF/rebinRRF); // = theory time resolution as close as possible to the data time resolution compatible with wRRF } -//cout << endl << ">> runData->GetDataBin(histoNo).size() = " << runData->GetDataBin(histoNo).size() << ", fData.GetValue()->size() * 10 = " << fData.GetValue()->size() * 10 << ", size = " << size << ", factor = " << factor << endl; - fData.SetTheoryTimeStart(fData.GetDataTimeStart()); - fData.SetTheoryTimeStep(fTimeResolution*factor); -//cout << endl << ">> size=" << size << ", startTime=" << startTime << ", fTimeResolution=" << fTimeResolution; + for (UInt_t i=0; iFunc(time, par, fFuncValues); + if (wRRF != 0.0) { + theoryValue *= 2.0*TMath::Cos(wRRF * time + phaseRRF); + } if (fabs(theoryValue) > 10.0) { // dirty hack needs to be fixed!! theoryValue = 0.0; } fData.AppendTheoryValue(theoryValue); } + // if RRF filter the theory with a FIR Kaiser low pass filter + if (wRRF != 0.0) { + // rebin theory to the RRF frequency + if (rebinRRF != 0) { + Double_t dval = 0.0; + PDoubleVector theo; + for (UInt_t i=0; isize(); i++) { + if ((i % rebinRRF == 0) && (i != 0)) { + theo.push_back(dval/rebinRRF); + dval = 0.0; + } + dval += fData.GetTheory()->at(i); + } + fData.SetTheoryTimeStart(fData.GetTheoryTimeStart()+static_cast(rebinRRF-1)*fData.GetTheoryTimeStep()/2.0); + fData.SetTheoryTimeStep(rebinRRF*fData.GetTheoryTimeStep()); + fData.ReplaceTheory(theo); + theo.clear(); + } + + // filter data + 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) + FilterData(); + } + // clean up par.clear(); diff --git a/src/include/PMusr.h b/src/include/PMusr.h index da90d725..32a42ae7 100644 --- a/src/include/PMusr.h +++ b/src/include/PMusr.h @@ -116,6 +116,14 @@ using namespace std; #define FOURIER_PLOT_POWER 4 #define FOURIER_PLOT_PHASE 5 +//------------------------------------------------------------- +// RRF related tags +#define RRF_UNIT_kHz 0 +#define RRF_UNIT_MHz 1 +#define RRF_UNIT_Mcs 2 +#define RRF_UNIT_G 3 +#define RRF_UNIT_T 4 + //------------------------------------------------------------- /** *

typedef to make to code more readable. @@ -201,6 +209,9 @@ class PRunData { virtual void AppendXTheoryValue(Double_t dval) { fXTheory.push_back(dval); } virtual void AppendTheoryValue(Double_t dval) { fTheory.push_back(dval); } + virtual void SetTheoryValue(UInt_t i, Double_t dval); + virtual void ReplaceTheory(const PDoubleVector &theo); + private: // data related info Double_t fDataTimeStart; @@ -403,12 +414,6 @@ class PMsrRunBlock { virtual Int_t GetT0(UInt_t i=0); virtual Double_t GetFitRange(UInt_t i); virtual Int_t GetPacking() { return fPacking; } - virtual Double_t GetRRFFreq() { return fRRFFreq; } - virtual Int_t GetRRFPacking() { return fRRFPacking; } - virtual Int_t GetAlpha2ParamNo() { return fAlpha2ParamNo; } - virtual Int_t GetBeta2ParamNo() { return fBeta2ParamNo; } - virtual Int_t GetRightHistoNo() { return fRightHistoNo; } - virtual Int_t GetLeftHistoNo() { return fLeftHistoNo; } virtual Int_t GetXDataIndex() { return fXYDataIndex[0]; } virtual Int_t GetYDataIndex() { return fXYDataIndex[1]; } virtual TString* GetXDataLabel() { return &fXYDataLabel[0]; } @@ -444,12 +449,6 @@ class PMsrRunBlock { virtual void SetT0(Int_t ival, UInt_t idx); virtual void SetFitRange(Double_t dval, UInt_t idx); virtual void SetPacking(Int_t ival) { fPacking = ival; } - virtual void SetRRFFreq(Double_t dval) { fRRFFreq = dval; } - virtual void SetRRFPacking(Int_t ival) { fRRFPacking = ival; } - virtual void SetAlpha2ParamNo(Int_t ival) { fAlpha2ParamNo = ival; } - virtual void SetBeta2ParamNo(Int_t ival) { fBeta2ParamNo = ival; } - virtual void SetRightHistoNo(Int_t ival) { fRightHistoNo = ival; } - virtual void SetLeftHistoNo(Int_t ival) { fLeftHistoNo = ival; } virtual void SetXDataIndex(Int_t ival) { fXYDataIndex[0] = ival; } virtual void SetYDataIndex(Int_t ival) { fXYDataIndex[1] = ival; } virtual void SetXDataLabel(TString& str) { fXYDataLabel[0] = str; } @@ -477,12 +476,6 @@ class PMsrRunBlock { PIntVector fT0; ///< t0 bins (fit type 0, 2, 4). if fit type 0 -> f0, f1, f2, ...; if fit type 2, 4 -> f0, b0, f1, b1, ... Double_t fFitRange[2]; ///< fit range in (us) Int_t fPacking; ///< packing/rebinning - Double_t fRRFFreq; ///< rotating reference frequency (fit type 4) - Int_t fRRFPacking; ///< rotating reference packing (fit type 4) - Int_t fAlpha2ParamNo; ///< rotating reference alpha2 (fit type 4) - Int_t fBeta2ParamNo; ///< rotating reference beta2 (fit type 4) - Int_t fRightHistoNo; ///< rotating reference right histogram number (fit type 4) - Int_t fLeftHistoNo; ///< rotating reference left histogram number (fit type 4) 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) }; @@ -525,6 +518,10 @@ typedef struct { PDoubleVector fTmax; ///< time maximum PDoubleVector fYmin; ///< asymmetry/counts minimum PDoubleVector fYmax; ///< asymmetry/counts maximum + UInt_t fRRFPacking; ///< rotating reference frame (RRF) packing + Double_t fRRFFreq; ///< RRF frequency + UInt_t fRRFUnit; ///< RRF frequency unit. 0=kHz, 1=MHz, 2=Mc/s, 3=Gauss, 4=Tesla + Double_t fRRFPhase; ///< RRF phase } PMsrPlotStructure; //------------------------------------------------------------- diff --git a/src/include/PMusrCanvas.h b/src/include/PMusrCanvas.h index fd2a0b8c..382108e1 100644 --- a/src/include/PMusrCanvas.h +++ b/src/include/PMusrCanvas.h @@ -204,6 +204,8 @@ class PMusrCanvas : public TObject, public TQObject Double_t fCurrentFourierPhase; /// holds the current Fourier phase TLatex *fCurrentFourierPhaseText; /// used in Re/Im Fourier to show the current phase in the pad + TString *fRRFText; + TLatex *fRRFLatexText; /// used to display RRF info TStyle *fStyle; diff --git a/src/include/PRunBase.h b/src/include/PRunBase.h index 78780964..87113656 100644 --- a/src/include/PRunBase.h +++ b/src/include/PRunBase.h @@ -79,10 +79,15 @@ class PRunBase Double_t fTimeResolution; ///< time resolution in (us) PIntVector fT0s; ///< all t0's of a run! The derived classes will handle it - virtual Bool_t PrepareData() = 0; // pure virtual, i.e. needs to be implemented by the deriving class!! - PDoubleVector fFuncValues; ///< is keeping the values of the functions from the FUNCTIONS block PTheory *fTheory; ///< theory needed to calculate chi-square + + PDoubleVector fKaiserFilter; + + virtual Bool_t PrepareData() = 0; // pure virtual, i.e. needs to be implemented by the deriving class!! + + virtual void CalculateKaiserFilterCoeff(Double_t wc, Double_t A, Double_t dw); + virtual void FilterData(); }; #endif // _PRUNBASE_H_ diff --git a/src/musrt0.cpp b/src/musrt0.cpp index e14166ca..259d66b4 100644 --- a/src/musrt0.cpp +++ b/src/musrt0.cpp @@ -265,6 +265,7 @@ int main(int argc, char *argv[]) } break; case MSR_FITTYPE_ASYM_RRF: +/* for (unsigned int j=0; jat(i).GetRunNameSize(); j++) { // necessary in case of ADDRUN if (!musrt0_item(app, msrHandler, dataHandler->GetRunData(*(runList->at(i).GetRunName(j))), i, runList->at(i).GetForwardHistoNo(), 0, j)) { musrt0_cleanup(saxParser, startupHandler, msrHandler, dataHandler); @@ -283,6 +284,7 @@ int main(int argc, char *argv[]) exit(0); } } +*/ break; default: break; diff --git a/src/tests/t0NotEqFirstGoodData/t0NotEqFirstGoodData.C b/src/tests/t0NotEqFirstGoodData/t0NotEqFirstGoodData.C index 7bb59b75..f596fa9c 100644 --- a/src/tests/t0NotEqFirstGoodData/t0NotEqFirstGoodData.C +++ b/src/tests/t0NotEqFirstGoodData/t0NotEqFirstGoodData.C @@ -68,17 +68,20 @@ void t0NotEqFirstGoodData() // create histos UInt_t t0[4] = {3419, 3419, 3419, 3419}; - UInt_t n0[4] = {200.0, 205.0, 203.0, 198.0}; + UInt_t n0[4] = {2000.0, 2005.0, 2003.0, 1998.0}; UInt_t bkg[4] = {11.0, 11.0, 5.0, 8.0}; const Double_t tau = 2197.019; // ns // asymmetry related stuff const Double_t timeResolution = 0.1953125; // ns - Double_t a[4] = {0.16, 0.161, 0.16, 0.159}; - Double_t phase[4] = {13.0*TMath::Pi()/180.0, 103.0*TMath::Pi()/180.0, 193.0*TMath::Pi()/180.0, 283.0*TMath::Pi()/180.0}; + Double_t a0[4] = {0.16, 0.16, 0.16, 0.16}; + Double_t a1[4] = {0.0, 0.0, 0.0, 0.0}; + Double_t phase[4] = {0.0*TMath::Pi()/180.0, 90.0*TMath::Pi()/180.0, 180.0*TMath::Pi()/180.0, 270.0*TMath::Pi()/180.0}; const Double_t gamma = 0.00001355; // gamma/(2pi) - Double_t bb = 15000.0; // field in Gauss + Double_t bb0 = 15000.0; // field in Gauss + Double_t bb1 = 15702.0; // field in Gauss // Double_t bb = 100.0; // field in Gauss - Double_t sigma = 1.05/1000.0; // Gaussian sigma in 1/ns + Double_t sigma0 = 0.05/1000.0; // Gaussian sigma in 1/ns + Double_t sigma1 = 0.10/1000.0; // Gaussian sigma in 1/ns TH1F *histo[8]; char str[128]; @@ -96,7 +99,8 @@ void t0NotEqFirstGoodData() histo[i]->SetBinContent(j, bkg[i]); } else { time = (Double_t)(j-t0[i])*timeResolution; - dval = (Double_t)n0[i]*TMath::Exp(-time/tau)*(1.0+a[i]*TMath::Exp(-0.5*TMath::Power(time*sigma,2.0))*TMath::Cos(TMath::TwoPi()*gamma*bb*time+phase[i]))+(Double_t)bkg[i]; + dval = (Double_t)n0[i]*TMath::Exp(-time/tau)*(1.0+a0[i]*TMath::Exp(-0.5*TMath::Power(time*sigma0,2.0))*TMath::Cos(TMath::TwoPi()*gamma*bb0*time+phase[i]) + +a1[i]*TMath::Exp(-0.5*TMath::Power(time*sigma1,2.0))*TMath::Cos(TMath::TwoPi()*gamma*bb1*time+phase[i]))+(Double_t)bkg[i]; histo[i]->SetBinContent(j, (UInt_t)dval); } } @@ -134,7 +138,7 @@ void t0NotEqFirstGoodData() decayAnaModule->Add(histo[i]); // write file - TFile *fout = new TFile("010105.root", "RECREATE", "Midas Fake Histograms"); + TFile *fout = new TFile("010108.root", "RECREATE", "Midas Fake Histograms"); if (fout == 0) { cout << endl << "**ERROR** Couldn't create ROOT file"; cout << endl << endl;