From b106ac61a110d7fb0ab2f498a3e144615c998ec5 Mon Sep 17 00:00:00 2001 From: Andreas Suter Date: Wed, 7 Dec 2016 17:02:18 +0100 Subject: [PATCH] added phase optimized Re FT for musrview. Still not perfect. --- src/classes/PMsrHandler.cpp | 12 +- src/classes/PMusrCanvas.cpp | 366 ++++++++++++++++++++++++++++-------- src/include/PMusrCanvas.h | 32 ++-- 3 files changed, 318 insertions(+), 92 deletions(-) diff --git a/src/classes/PMsrHandler.cpp b/src/classes/PMsrHandler.cpp index f5b75804..d126ccd9 100644 --- a/src/classes/PMsrHandler.cpp +++ b/src/classes/PMsrHandler.cpp @@ -1113,8 +1113,10 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t messages) fout << "POWER"; } else if (fFourier.fPlotTag == FOURIER_PLOT_PHASE) { fout << "PHASE"; + } else if (fFourier.fPlotTag == FOURIER_PLOT_PHASE_OPT_REAL) { + fout << "PHASE_OPT_REAL"; } - fout << " # REAL, IMAG, REAL_AND_IMAG, POWER, PHASE"; + fout << " # REAL, IMAG, REAL_AND_IMAG, POWER, PHASE, PHASE_OPT_REAL"; fout << endl; } else if (sstr.BeginsWith("phase")) { if (fFourier.fPhaseParamNo > 0) { @@ -2146,8 +2148,10 @@ Int_t PMsrHandler::WriteMsrFile(const Char_t *filename, map *co fout << "POWER"; } else if (fFourier.fPlotTag == FOURIER_PLOT_PHASE) { fout << "PHASE"; + } else if (fFourier.fPlotTag == FOURIER_PLOT_PHASE_OPT_REAL) { + fout << "PHASE_OPT_REAL"; } - fout << " # REAL, IMAG, REAL_AND_IMAG, POWER, PHASE"; + fout << " # REAL, IMAG, REAL_AND_IMAG, POWER, PHASE, PHASE_OPT_REAL"; fout << endl; } @@ -3951,6 +3955,8 @@ Bool_t PMsrHandler::HandleFourierEntry(PMsrLines &lines) fourier.fPlotTag = FOURIER_PLOT_POWER; } else if (!str.CompareTo("phase", TString::kIgnoreCase)) { fourier.fPlotTag = FOURIER_PLOT_PHASE; + } else if (!str.CompareTo("phase_opt_real", TString::kIgnoreCase)) { + fourier.fPlotTag = FOURIER_PLOT_PHASE_OPT_REAL; } else { // unrecognized plot tag error = true; continue; @@ -4087,7 +4093,7 @@ Bool_t PMsrHandler::HandleFourierEntry(PMsrLines &lines) cerr << endl << ">> 0 <= n <= 20 are allowed values"; cerr << endl << ">> [dc-corrected true | false]"; cerr << endl << ">> [apodization none | weak | medium | strong]"; - cerr << endl << ">> [plot real | imag | real_and_imag | power | phase]"; + cerr << endl << ">> [plot real | imag | real_and_imag | power | phase | phase_opt_real]"; cerr << endl << ">> [phase value]"; cerr << endl << ">> [range_for_phase_correction min max | all]"; cerr << endl << ">> [range min max]"; diff --git a/src/classes/PMusrCanvas.cpp b/src/classes/PMusrCanvas.cpp index 44d31939..2cae4d9e 100644 --- a/src/classes/PMusrCanvas.cpp +++ b/src/classes/PMusrCanvas.cpp @@ -815,6 +815,12 @@ void PMusrCanvas::UpdateDataTheoryPad() fPopupFourier->CheckEntry(P_MENU_ID_FOURIER+P_MENU_PLOT_OFFSET*fPlotNumber+P_MENU_ID_FOURIER_PHASE); } break; + case FOURIER_PLOT_PHASE_OPT_REAL: + fCurrentPlotView = PV_FOURIER_PHASE_OPT_REAL; + if (!fBatchMode) { + fPopupFourier->CheckEntry(P_MENU_ID_FOURIER+P_MENU_PLOT_OFFSET*fPlotNumber+P_MENU_ID_FOURIER_PHASE_OPT_REAL); + } + break; default: fCurrentPlotView = PV_FOURIER_PWR; if (!fBatchMode) { @@ -1096,6 +1102,11 @@ void PMusrCanvas::HandleCmdKey(Int_t event, Int_t x, Int_t y, TObject *selected) fCurrentPlotView = PV_FOURIER_PHASE; fPopupFourier->CheckEntry(P_MENU_ID_FOURIER+P_MENU_PLOT_OFFSET*fPlotNumber+P_MENU_ID_FOURIER_PHASE); break; + case FOURIER_PLOT_PHASE_OPT_REAL: + fPreviousPlotView = fCurrentPlotView; + fCurrentPlotView = PV_FOURIER_PHASE_OPT_REAL; + fPopupFourier->CheckEntry(P_MENU_ID_FOURIER+P_MENU_PLOT_OFFSET*fPlotNumber+P_MENU_ID_FOURIER_PHASE_OPT_REAL); + break; default: break; } @@ -1342,6 +1353,29 @@ void PMusrCanvas::HandleMenuPopup(Int_t id) HandleFourierDifference(); PlotFourierDifference(); } + } else if (id == P_MENU_ID_FOURIER+P_MENU_PLOT_OFFSET*fPlotNumber+P_MENU_ID_FOURIER_PHASE_OPT_REAL) { + // set appropriate plot view + fPreviousPlotView = fCurrentPlotView; + fCurrentPlotView = PV_FOURIER_PHASE_OPT_REAL; + // uncheck data + fPopupMain->UnCheckEntry(P_MENU_ID_DATA+P_MENU_PLOT_OFFSET*fPlotNumber); + // check appropriate fourier popup item + fPopupFourier->UnCheckEntries(); + fPopupFourier->CheckEntry(id); + // enable phase increment/decrement + fPopupFourier->EnableEntry(P_MENU_ID_FOURIER+P_MENU_PLOT_OFFSET*fPlotNumber+P_MENU_ID_FOURIER_PHASE_PLUS); + fPopupFourier->EnableEntry(P_MENU_ID_FOURIER+P_MENU_PLOT_OFFSET*fPlotNumber+P_MENU_ID_FOURIER_PHASE_MINUS); + // handle fourier phase + if (!fDifferenceView) { + HandleFourier(); + PlotFourier(); + } else { + if (previousPlotView == PV_DATA) + HandleDifferenceFourier(); + else + HandleFourierDifference(); + PlotFourierDifference(); + } } else if (id == P_MENU_ID_FOURIER+P_MENU_PLOT_OFFSET*fPlotNumber+P_MENU_ID_FOURIER_PHASE_PLUS) { IncrementFourierPhase(); } else if (id == P_MENU_ID_FOURIER+P_MENU_PLOT_OFFSET*fPlotNumber+P_MENU_ID_FOURIER_PHASE_MINUS) { @@ -1749,6 +1783,24 @@ void PMusrCanvas::ExportData(const Char_t *fileName) } } break; + case PV_FOURIER_PHASE_OPT_REAL: + // get current x-range + xminBin = fData[0].dataFourierPhaseOptReal->GetXaxis()->GetFirst(); // first bin of the zoomed range + xmaxBin = fData[0].dataFourierPhaseOptReal->GetXaxis()->GetLast(); // last bin of the zoomed range + xmin = fData[0].dataFourierPhaseOptReal->GetXaxis()->GetBinCenter(xminBin); + xmax = fData[0].dataFourierPhaseOptReal->GetXaxis()->GetBinCenter(xmaxBin); + + // fill ascii dump data + if (fAveragedView) { + GetExportDataSet(fDataAvg.dataFourierPhaseOptReal, xmin, xmax, dumpVector, false); + GetExportDataSet(fDataAvg.theoryFourierPhaseOptReal, xmin, xmax, dumpVector, false); + } else { // go through all the histogramms + for (UInt_t i=0; iAddEntry("Show Real+Imag", P_MENU_ID_FOURIER+P_MENU_PLOT_OFFSET*fPlotNumber+P_MENU_ID_FOURIER_REAL_AND_IMAG); fPopupFourier->AddEntry("Show Power", P_MENU_ID_FOURIER+P_MENU_PLOT_OFFSET*fPlotNumber+P_MENU_ID_FOURIER_PWR); fPopupFourier->AddEntry("Show Phase", P_MENU_ID_FOURIER+P_MENU_PLOT_OFFSET*fPlotNumber+P_MENU_ID_FOURIER_PHASE); + fPopupFourier->AddEntry("Show PhaseOptReal", P_MENU_ID_FOURIER+P_MENU_PLOT_OFFSET*fPlotNumber+P_MENU_ID_FOURIER_PHASE_OPT_REAL); fPopupFourier->AddSeparator(); fPopupFourier->AddEntry("Phase +", P_MENU_ID_FOURIER+P_MENU_PLOT_OFFSET*fPlotNumber+P_MENU_ID_FOURIER_PHASE_PLUS); fPopupFourier->AddEntry("Phase -", P_MENU_ID_FOURIER+P_MENU_PLOT_OFFSET*fPlotNumber+P_MENU_ID_FOURIER_PHASE_MINUS); @@ -2417,16 +2473,19 @@ void PMusrCanvas::InitDataSet(PMusrCanvasDataSet &dataSet) dataSet.dataFourierIm = 0; dataSet.dataFourierPwr = 0; dataSet.dataFourierPhase = 0; + dataSet.dataFourierPhaseOptReal = 0; dataSet.theory = 0; dataSet.theoryFourierRe = 0; dataSet.theoryFourierIm = 0; dataSet.theoryFourierPwr = 0; dataSet.theoryFourierPhase = 0; + dataSet.theoryFourierPhaseOptReal = 0; dataSet.diff = 0; dataSet.diffFourierRe = 0; dataSet.diffFourierIm = 0; dataSet.diffFourierPwr = 0; dataSet.diffFourierPhase = 0; + dataSet.diffFourierPhaseOptReal = 0; dataSet.dataRange = 0; } @@ -2488,6 +2547,10 @@ void PMusrCanvas::CleanupDataSet(PMusrCanvasDataSet &dataSet) delete dataSet.dataFourierPhase; dataSet.dataFourierPhase = 0; } + if (dataSet.dataFourierPhaseOptReal) { + delete dataSet.dataFourierPhaseOptReal; + dataSet.dataFourierPhaseOptReal = 0; + } if (dataSet.theory) { delete dataSet.theory; dataSet.theory = 0; @@ -2508,6 +2571,10 @@ void PMusrCanvas::CleanupDataSet(PMusrCanvasDataSet &dataSet) delete dataSet.theoryFourierPhase; dataSet.theoryFourierPhase = 0; } + if (dataSet.theoryFourierPhaseOptReal) { + delete dataSet.theoryFourierPhaseOptReal; + dataSet.theoryFourierPhaseOptReal = 0; + } if (dataSet.diff) { delete dataSet.diff; dataSet.diff = 0; @@ -2528,6 +2595,10 @@ void PMusrCanvas::CleanupDataSet(PMusrCanvasDataSet &dataSet) delete dataSet.diffFourierPhase; dataSet.diffFourierPhase = 0; } + if (dataSet.diffFourierPhaseOptReal) { + delete dataSet.diffFourierPhaseOptReal; + dataSet.diffFourierPhaseOptReal = 0; + } if (dataSet.dataRange) { delete dataSet.dataRange; dataSet.dataRange = 0; @@ -3260,6 +3331,11 @@ void PMusrCanvas::HandleDifference() */ void PMusrCanvas::HandleFourier() { + PDoubleVector phaseParam; + Double_t min=-1.0, max=-1.0; + Double_t re, im, ph; + char hName[1024]; + // check if plot type is appropriate for fourier if (fPlotType == MSR_PLOT_NON_MUSR) return; @@ -3269,7 +3345,7 @@ void PMusrCanvas::HandleFourier() Int_t bin; double startTime = fXmin; double endTime = fXmax; - if (!fStartWithFourier) { // fHistoFrame presen, hence get start/end from it + if (!fStartWithFourier) { // fHistoFrame present, hence get start/end from it bin = fHistoFrame->GetXaxis()->GetFirst(); startTime = fHistoFrame->GetBinCenter(bin); bin = fHistoFrame->GetXaxis()->GetLast(); @@ -3294,6 +3370,12 @@ void PMusrCanvas::HandleFourier() // get phase part of the data fData[i].dataFourierPhase = fourierData.GetPhaseFourier(); + // get phase optimized real fourier from data + min = fMsrHandler->GetMsrFourierList()->fRangeForPhaseCorrection[0]; + max = fMsrHandler->GetMsrFourierList()->fRangeForPhaseCorrection[1]; + // eventually min/max needs to be extracted from the 'range_for_phase_correction' parameter of the msr-file + fData[i].dataFourierPhaseOptReal = fourierData.GetPhaseOptRealFourier(phaseParam, scale, min, max); + // set marker and line color fData[i].dataFourierRe->SetMarkerColor(fData[i].data->GetMarkerColor()); fData[i].dataFourierRe->SetLineColor(fData[i].data->GetLineColor()); @@ -3303,17 +3385,22 @@ void PMusrCanvas::HandleFourier() fData[i].dataFourierPwr->SetLineColor(fData[i].data->GetLineColor()); fData[i].dataFourierPhase->SetMarkerColor(fData[i].data->GetMarkerColor()); fData[i].dataFourierPhase->SetLineColor(fData[i].data->GetLineColor()); + fData[i].dataFourierPhaseOptReal->SetMarkerColor(fData[i].data->GetMarkerColor()); + fData[i].dataFourierPhaseOptReal->SetLineColor(fData[i].data->GetLineColor()); // set marker size fData[i].dataFourierRe->SetMarkerSize(1); fData[i].dataFourierIm->SetMarkerSize(1); fData[i].dataFourierPwr->SetMarkerSize(1); fData[i].dataFourierPhase->SetMarkerSize(1); + fData[i].dataFourierPhaseOptReal->SetMarkerSize(1); + // set marker type fData[i].dataFourierRe->SetMarkerStyle(fData[i].data->GetMarkerStyle()); fData[i].dataFourierIm->SetMarkerStyle(fData[i].data->GetMarkerStyle()); fData[i].dataFourierPwr->SetMarkerStyle(fData[i].data->GetMarkerStyle()); fData[i].dataFourierPhase->SetMarkerStyle(fData[i].data->GetMarkerStyle()); + fData[i].dataFourierPhaseOptReal->SetMarkerStyle(fData[i].data->GetMarkerStyle()); // calculate fourier transform of the theory Int_t powerPad = (Int_t)round(log((endTime-startTime)/fData[i].theory->GetBinWidth(1))/log(2))+3; @@ -3333,16 +3420,33 @@ void PMusrCanvas::HandleFourier() // get phase part of the data fData[i].theoryFourierPhase = fourierTheory.GetPhaseFourier(); + // get phase optimized real fourier from data + + // clone theory Re FT + strcpy(hName, fData[i].theoryFourierPhase->GetName()); + strcat(hName, "_Opt_Real"); + fData[i].theoryFourierPhaseOptReal = (TH1F*) fData[i].theoryFourierRe->Clone(hName); + + // rotate the theory according to the optimized phase parameters + // first find minBin for min of the phase correction + Int_t minBin = fData[i].theoryFourierPhaseOptReal->GetXaxis()->FindBin(min); + Int_t maxBin = fData[i].theoryFourierPhaseOptReal->GetXaxis()->FindBin(max); + for (Int_t j=1; jGetNbinsX(); j++) { + ph = phaseParam[0] + phaseParam[1] * (Double_t)(j-minBin+1) / (Double_t)(maxBin-minBin); + re = fData[i].theoryFourierRe->GetBinContent(j) * cos(ph) - fData[i].theoryFourierIm->GetBinContent(j) * sin(ph); + fData[i].theoryFourierPhaseOptReal->SetBinContent(j, re); + } + // set line colors for the theory fData[i].theoryFourierRe->SetLineColor(fData[i].theory->GetLineColor()); fData[i].theoryFourierIm->SetLineColor(fData[i].theory->GetLineColor()); fData[i].theoryFourierPwr->SetLineColor(fData[i].theory->GetLineColor()); fData[i].theoryFourierPhase->SetLineColor(fData[i].theory->GetLineColor()); + fData[i].theoryFourierPhaseOptReal->SetLineColor(fData[i].theory->GetLineColor()); } // apply global phase if present if (fFourier.fPhase != 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()); @@ -3372,6 +3476,7 @@ void PMusrCanvas::HandleFourier() } } +/* as: will be obsolate ... // find optimal Fourier phase if range is given if ((fFourier.fRangeForPhaseCorrection[0] != -1.0) && (fFourier.fRangeForPhaseCorrection[1] != -1.0)) { @@ -3405,6 +3510,7 @@ void PMusrCanvas::HandleFourier() } } } +*/ } } @@ -3546,6 +3652,12 @@ void PMusrCanvas::HandleFourierDifference() fData[i].dataFourierPhase->GetXaxis()->GetXmin(), fData[i].dataFourierPhase->GetXaxis()->GetXmax()); + // phase optimized real part + name = TString(fData[i].dataFourierPhaseOptReal->GetTitle()) + "_diff"; + fData[i].diffFourierPhaseOptReal = new TH1F(name, name, fData[i].dataFourierPhaseOptReal->GetNbinsX(), + fData[i].dataFourierPhaseOptReal->GetXaxis()->GetXmin(), + fData[i].dataFourierPhaseOptReal->GetXaxis()->GetXmax()); + // calculate difference for (UInt_t j=1; jGetEntries(); j++) { dvalx = fData[i].dataFourierRe->GetXaxis()->GetBinCenter(j); // get x-axis value of bin j @@ -3564,6 +3676,10 @@ void PMusrCanvas::HandleFourierDifference() theoBin = fData[i].theoryFourierPhase->FindBin(dvalx); // get the theory x-axis bin dval = fData[i].dataFourierPhase->GetBinContent(j) - fData[i].theoryFourierPhase->GetBinContent(theoBin); fData[i].diffFourierPhase->SetBinContent(j, dval); + dvalx = fData[i].dataFourierPhaseOptReal->GetXaxis()->GetBinCenter(j); // get x-axis value of bin j + theoBin = fData[i].theoryFourierPhaseOptReal->FindBin(dvalx); // get the theory x-axis bin + dval = fData[i].dataFourierPhaseOptReal->GetBinContent(j) - fData[i].theoryFourierPhaseOptReal->GetBinContent(theoBin); + fData[i].diffFourierPhaseOptReal->SetBinContent(j, dval); } } @@ -3577,17 +3693,21 @@ void PMusrCanvas::HandleFourierDifference() fData[i].diffFourierPwr->SetLineColor(fData[i].dataFourierPwr->GetLineColor()); fData[i].diffFourierPhase->SetMarkerColor(fData[i].dataFourierPhase->GetMarkerColor()); fData[i].diffFourierPhase->SetLineColor(fData[i].dataFourierPhase->GetLineColor()); + fData[i].diffFourierPhaseOptReal->SetMarkerColor(fData[i].dataFourierPhaseOptReal->GetMarkerColor()); + fData[i].diffFourierPhaseOptReal->SetLineColor(fData[i].dataFourierPhaseOptReal->GetLineColor()); // set marker size fData[i].diffFourierRe->SetMarkerSize(1); fData[i].diffFourierIm->SetMarkerSize(1); fData[i].diffFourierPwr->SetMarkerSize(1); fData[i].diffFourierPhase->SetMarkerSize(1); + fData[i].diffFourierPhaseOptReal->SetMarkerSize(1); // set marker type fData[i].diffFourierRe->SetMarkerStyle(fData[i].dataFourierRe->GetMarkerStyle()); fData[i].diffFourierIm->SetMarkerStyle(fData[i].dataFourierIm->GetMarkerStyle()); fData[i].diffFourierPwr->SetMarkerStyle(fData[i].dataFourierPwr->GetMarkerStyle()); fData[i].diffFourierPhase->SetMarkerStyle(fData[i].dataFourierPhase->GetMarkerStyle()); + fData[i].diffFourierPhaseOptReal->SetMarkerStyle(fData[i].dataFourierPhaseOptReal->GetMarkerStyle()); // set diffFourierTag fData[i].diffFourierTag = 2; // f-d @@ -3646,6 +3766,12 @@ void PMusrCanvas::HandleAverage() fData[0].dataFourierPhase->GetXaxis()->GetXmin(), fData[0].dataFourierPhase->GetXaxis()->GetXmax()); } + if (fData[0].dataFourierPhaseOptReal != 0) { + name = TString(fData[0].dataFourierPhaseOptReal->GetTitle()) + "_avg"; + fDataAvg.dataFourierPhaseOptReal = new TH1F(name, name, fData[0].dataFourierPhaseOptReal->GetNbinsX(), + fData[0].dataFourierPhaseOptReal->GetXaxis()->GetXmin(), + fData[0].dataFourierPhaseOptReal->GetXaxis()->GetXmax()); + } if (fData[0].theory != 0) { name = TString(fData[0].theory->GetTitle()) + "_avg"; fDataAvg.theory = new TH1F(name, name, fData[0].theory->GetNbinsX(), @@ -3676,6 +3802,12 @@ void PMusrCanvas::HandleAverage() fData[0].theoryFourierPhase->GetXaxis()->GetXmin(), fData[0].theoryFourierPhase->GetXaxis()->GetXmax()); } + if (fData[0].theoryFourierPhaseOptReal != 0) { + name = TString(fData[0].theoryFourierPhaseOptReal->GetTitle()) + "_avg"; + fDataAvg.theoryFourierPhaseOptReal = new TH1F(name, name, fData[0].theoryFourierPhaseOptReal->GetNbinsX(), + fData[0].theoryFourierPhaseOptReal->GetXaxis()->GetXmin(), + fData[0].theoryFourierPhaseOptReal->GetXaxis()->GetXmax()); + } if (fData[0].diff != 0) { name = TString(fData[0].diff->GetTitle()) + "_avg"; fDataAvg.diff = new TH1F(name, name, fData[0].diff->GetNbinsX(), @@ -3706,6 +3838,12 @@ void PMusrCanvas::HandleAverage() fData[0].diffFourierPhase->GetXaxis()->GetXmin(), fData[0].diffFourierPhase->GetXaxis()->GetXmax()); } + if (fData[0].diffFourierPhaseOptReal != 0) { + name = TString(fData[0].diffFourierPhaseOptReal->GetTitle()) + "_avg"; + fDataAvg.diffFourierPhaseOptReal = new TH1F(name, name, fData[0].diffFourierPhaseOptReal->GetNbinsX(), + fData[0].diffFourierPhaseOptReal->GetXaxis()->GetXmin(), + fData[0].diffFourierPhaseOptReal->GetXaxis()->GetXmax()); + } // calculate all the average data sets double dval; @@ -3779,6 +3917,20 @@ void PMusrCanvas::HandleAverage() fDataAvg.dataFourierPhase->SetMarkerSize(fData[0].dataFourierPhase->GetMarkerSize()); fDataAvg.dataFourierPhase->SetMarkerStyle(fData[0].dataFourierPhase->GetMarkerStyle()); } + if (fDataAvg.dataFourierPhaseOptReal != 0) { + for (Int_t i=0; iGetNbinsX(); i++) { + dval = 0.0; + for (UInt_t j=0; jGetBinCenter(i)); + } + fDataAvg.dataFourierPhaseOptReal->SetBinContent(i, dval/fData.size()); + } + // set marker color, line color, maker size, marker type + fDataAvg.dataFourierPhaseOptReal->SetMarkerColor(fData[0].dataFourierPhaseOptReal->GetMarkerColor()); + fDataAvg.dataFourierPhaseOptReal->SetLineColor(fData[0].dataFourierPhaseOptReal->GetLineColor()); + fDataAvg.dataFourierPhaseOptReal->SetMarkerSize(fData[0].dataFourierPhaseOptReal->GetMarkerSize()); + fDataAvg.dataFourierPhaseOptReal->SetMarkerStyle(fData[0].dataFourierPhaseOptReal->GetMarkerStyle()); + } if (fDataAvg.theory != 0) { for (Int_t i=0; iGetNbinsX(); i++) { dval = 0.0; @@ -3845,6 +3997,20 @@ void PMusrCanvas::HandleAverage() fDataAvg.theoryFourierPhase->SetMarkerSize(fData[0].theoryFourierPhase->GetMarkerSize()); fDataAvg.theoryFourierPhase->SetMarkerStyle(fData[0].theoryFourierPhase->GetMarkerStyle()); } + if (fDataAvg.theoryFourierPhaseOptReal != 0) { + for (Int_t i=0; iGetNbinsX(); i++) { + dval = 0.0; + for (UInt_t j=0; jGetBinCenter(i)); + } + fDataAvg.theoryFourierPhaseOptReal->SetBinContent(i, dval/fData.size()); + } + // set marker color, line color, maker size, marker type + fDataAvg.theoryFourierPhaseOptReal->SetMarkerColor(fData[0].theoryFourierPhaseOptReal->GetMarkerColor()); + fDataAvg.theoryFourierPhaseOptReal->SetLineColor(fData[0].theoryFourierPhaseOptReal->GetLineColor()); + fDataAvg.theoryFourierPhaseOptReal->SetMarkerSize(fData[0].theoryFourierPhaseOptReal->GetMarkerSize()); + fDataAvg.theoryFourierPhaseOptReal->SetMarkerStyle(fData[0].theoryFourierPhaseOptReal->GetMarkerStyle()); + } if (fDataAvg.diff != 0) { for (Int_t i=0; iGetNbinsX(); i++) { dval = 0.0; @@ -3910,83 +4076,25 @@ void PMusrCanvas::HandleAverage() fDataAvg.diffFourierPhase->SetBinContent(i, dval/fData.size()); } // set marker color, line color, maker size, marker type - fDataAvg.diffFourierPhase->SetMarkerColor(fData[0].dataFourierRe->GetMarkerColor()); - fDataAvg.diffFourierPhase->SetLineColor(fData[0].dataFourierRe->GetLineColor()); - fDataAvg.diffFourierPhase->SetMarkerSize(fData[0].dataFourierRe->GetMarkerSize()); - fDataAvg.diffFourierPhase->SetMarkerStyle(fData[0].dataFourierRe->GetMarkerStyle()); + fDataAvg.diffFourierPhase->SetMarkerColor(fData[0].dataFourierPhase->GetMarkerColor()); + fDataAvg.diffFourierPhase->SetLineColor(fData[0].dataFourierPhase->GetLineColor()); + fDataAvg.diffFourierPhase->SetMarkerSize(fData[0].dataFourierPhase->GetMarkerSize()); + fDataAvg.diffFourierPhase->SetMarkerStyle(fData[0].dataFourierPhase->GetMarkerStyle()); } -} - -//-------------------------------------------------------------------------- -// FindOptimalFourierPhase (private) -//-------------------------------------------------------------------------- -/** - *

The idea to estimate the optimal phase is that the imaginary part of the fourier should be - * an antisymmetric function around the resonance, hence the asymmetry defined as asymmetry = max+min, - * where max/min is the maximum and minimum of the imaginary part, should be a minimum for the correct phase. - */ -double PMusrCanvas::FindOptimalFourierPhase() -{ - // check that Fourier is really present - if ((fData[0].dataFourierRe == 0) || (fData[0].dataFourierIm == 0)) - return 0.0; - - Double_t minPhase, x, valIm, val_xMin = 0.0, val_xMax = 0.0; - Double_t minIm = 0.0, maxIm = 0.0, asymmetry; - // get min/max of the imaginary part for phase = 0.0 as a starting point - minPhase = 0.0; - Bool_t first = true; - for (Int_t i=0; iGetNbinsX(); i++) { - x = fData[0].dataFourierIm->GetBinCenter(i); - if ((x > fFourier.fRangeForPhaseCorrection[0]) && (x < fFourier.fRangeForPhaseCorrection[1])) { - valIm = fData[0].dataFourierIm->GetBinContent(i); - if (first) { - minIm = valIm; - maxIm = valIm; - val_xMin = valIm; - first = false; - } else { - if (valIm < minIm) - minIm = valIm; - if (valIm > maxIm) - maxIm = valIm; - val_xMax = valIm; + if (fDataAvg.diffFourierPhaseOptReal != 0) { + for (Int_t i=0; iGetNbinsX(); i++) { + dval = 0.0; + for (UInt_t j=0; jGetBinCenter(i)); } + fDataAvg.diffFourierPhaseOptReal->SetBinContent(i, dval/fData.size()); } + // set marker color, line color, maker size, marker type + fDataAvg.diffFourierPhaseOptReal->SetMarkerColor(fData[0].dataFourierPhaseOptReal->GetMarkerColor()); + fDataAvg.diffFourierPhaseOptReal->SetLineColor(fData[0].dataFourierPhaseOptReal->GetLineColor()); + fDataAvg.diffFourierPhaseOptReal->SetMarkerSize(fData[0].dataFourierPhaseOptReal->GetMarkerSize()); + fDataAvg.diffFourierPhaseOptReal->SetMarkerStyle(fData[0].dataFourierPhaseOptReal->GetMarkerStyle()); } - asymmetry = (maxIm+minIm)*(val_xMin-val_xMax); - - // go through all phases an check if there is a larger max-min value of the imaginary part - double cp, sp; - for (double phase=0.1; phase < 180.0; phase += 0.1) { - cp = TMath::Cos(phase / 180.0 * TMath::Pi()); - sp = TMath::Sin(phase / 180.0 * TMath::Pi()); - first = true; - for (Int_t i=0; iGetNbinsX(); i++) { - x = fData[0].dataFourierIm->GetBinCenter(i); - if ((x > fFourier.fRangeForPhaseCorrection[0]) && (x < fFourier.fRangeForPhaseCorrection[1])) { - valIm = -sp * fData[0].dataFourierRe->GetBinContent(i) + cp * fData[0].dataFourierIm->GetBinContent(i); - if (first) { - minIm = valIm; - maxIm = valIm; - val_xMin = valIm; - first = false; - } else { - if (valIm < minIm) - minIm = valIm; - if (valIm > maxIm) - maxIm = valIm; - val_xMax = valIm; - } - } - } - if (fabs(asymmetry) > fabs((maxIm+minIm)*(val_xMin-val_xMax))) { - minPhase = phase; - asymmetry = (maxIm+minIm)*(val_xMin-val_xMax); - } - } - - return minPhase; } //-------------------------------------------------------------------------- @@ -4030,6 +4138,10 @@ void PMusrCanvas::CleanupFourier() delete fData[i].dataFourierPhase; fData[i].dataFourierPhase = 0; } + if (fData[i].dataFourierPhaseOptReal != 0) { + delete fData[i].dataFourierPhaseOptReal; + fData[i].dataFourierPhaseOptReal = 0; + } if (fData[i].theoryFourierRe != 0) { delete fData[i].theoryFourierRe; fData[i].theoryFourierRe = 0; @@ -4046,6 +4158,10 @@ void PMusrCanvas::CleanupFourier() delete fData[i].theoryFourierPhase; fData[i].theoryFourierPhase = 0; } + if (fData[i].theoryFourierPhaseOptReal != 0) { + delete fData[i].theoryFourierPhaseOptReal; + fData[i].theoryFourierPhaseOptReal = 0; + } } } @@ -4105,6 +4221,10 @@ void PMusrCanvas::CleanupAverage() delete fDataAvg.dataFourierPhase; fDataAvg.dataFourierPhase = 0; } + if (fDataAvg.dataFourierPhaseOptReal != 0) { + delete fDataAvg.dataFourierPhaseOptReal; + fDataAvg.dataFourierPhaseOptReal = 0; + } if (fDataAvg.theory != 0) { delete fDataAvg.theory; fDataAvg.theory = 0; @@ -4125,6 +4245,10 @@ void PMusrCanvas::CleanupAverage() delete fDataAvg.theoryFourierPhase; fDataAvg.theoryFourierPhase = 0; } + if (fDataAvg.theoryFourierPhaseOptReal != 0) { + delete fDataAvg.theoryFourierPhaseOptReal; + fDataAvg.theoryFourierPhaseOptReal = 0; + } if (fDataAvg.diff != 0) { delete fDataAvg.diff; fDataAvg.diff = 0; @@ -4145,6 +4269,10 @@ void PMusrCanvas::CleanupAverage() delete fDataAvg.diffFourierPhase; fDataAvg.diffFourierPhase = 0; } + if (fDataAvg.diffFourierPhaseOptReal != 0) { + delete fDataAvg.diffFourierPhaseOptReal; + fDataAvg.diffFourierPhaseOptReal = 0; + } } //-------------------------------------------------------------------------- @@ -5380,6 +5508,83 @@ void PMusrCanvas::PlotFourier(Bool_t unzoom) } break; + case PV_FOURIER_PHASE_OPT_REAL: + // set x-range + if ((fFourier.fPlotRange[0] != -1) && (fFourier.fPlotRange[1] != -1)) { + xmin = fFourier.fPlotRange[0]; + xmax = fFourier.fPlotRange[1]; + } else { + xmin = fData[0].dataFourierPhaseOptReal->GetBinLowEdge(1); + xmax = fData[0].dataFourierPhaseOptReal->GetBinLowEdge(fData[0].dataFourierPhaseOptReal->GetNbinsX())+fData[0].dataFourierPhaseOptReal->GetBinWidth(1); + } + + // set y-range + // first find minimum/maximum of all histos + ymin = GetMinimum(fData[0].dataFourierPhaseOptReal); + ymax = GetMaximum(fData[0].dataFourierPhaseOptReal); + binContent = GetMinimum(fData[0].theoryFourierPhaseOptReal); + if (binContent < ymin) + ymin = binContent; + binContent = GetMaximum(fData[0].theoryFourierPhaseOptReal); + if (binContent > ymax) + ymax = binContent; + for (UInt_t i=1; i ymax) + ymax = binContent; + binContent = GetMinimum(fData[i].theoryFourierPhaseOptReal); + if (binContent < ymin) + ymin = binContent; + binContent = GetMaximum(fData[i].theoryFourierPhaseOptReal); + if (binContent > ymax) + ymax = binContent; + } + + // delete old fHistoFrame if present + if (fHistoFrame) { + delete fHistoFrame; + fHistoFrame = 0; + } + + fHistoFrame = fDataTheoryPad->DrawFrame(xmin, 1.05*ymin, xmax, 1.05*ymax); + + // find the maximal number of points present in the histograms and increase the default number of points of fHistoFrame (1000) to the needed one + noOfPoints = 1000; + for (UInt_t i=0; iGetNbinsX() > (Int_t)noOfPoints) + noOfPoints = fData[i].dataFourierPhaseOptReal->GetNbinsX(); + } + noOfPoints *= 2; // make sure that there are enough points + fHistoFrame->SetBins(noOfPoints, xmin, xmax); + + for (UInt_t i=0; iGetXaxis()->SetRangeUser(xmin, xmax); + fData[i].dataFourierPhaseOptReal->GetYaxis()->SetRangeUser(1.05*ymin, 1.05*ymax); + fData[i].theoryFourierPhaseOptReal->GetXaxis()->SetRangeUser(xmin, xmax); + fData[i].theoryFourierPhaseOptReal->GetYaxis()->SetRangeUser(1.05*ymin, 1.05*ymax); + } + + // set x-axis title + fHistoFrame->GetXaxis()->SetTitle(xAxisTitle.Data()); + + // set y-axis title + fHistoFrame->GetYaxis()->SetTitleOffset(1.3); + fHistoFrame->GetYaxis()->SetTitle("Phase Opt. Real Fourier"); + + // plot data + for (UInt_t i=0; iDraw("psame"); + } + + // plot theories + for (UInt_t i=0; iDraw("same"); + } + + break; default: break; } @@ -5848,6 +6053,9 @@ void PMusrCanvas::PlotAverage(Bool_t unzoom) case PV_FOURIER_PHASE: yAxisTitle = ""; break; + case PV_FOURIER_PHASE_OPT_REAL: + yAxisTitle = ""; + break; default: yAxisTitle = "??"; break; @@ -5967,6 +6175,14 @@ void PMusrCanvas::PlotAverage(Bool_t unzoom) fDataAvg.diffFourierPhase->Draw("psame"); } break; + case PV_FOURIER_PHASE_OPT_REAL: + if (!fDifferenceView) { // averaged Fourier Phase Opt Real view + fDataAvg.dataFourierPhaseOptReal->Draw("psame"); + fDataAvg.theoryFourierPhaseOptReal->Draw("same"); + } else { // averaged diff Fourier Phase view + fDataAvg.diffFourierPhaseOptReal->Draw("psame"); + } + break; default: break; } diff --git a/src/include/PMusrCanvas.h b/src/include/PMusrCanvas.h index 42e73dcc..cdc6791a 100644 --- a/src/include/PMusrCanvas.h +++ b/src/include/PMusrCanvas.h @@ -55,12 +55,13 @@ #define XTHEO 0.75 // Current Plot Views -#define PV_DATA 1 -#define PV_FOURIER_REAL 2 -#define PV_FOURIER_IMAG 3 -#define PV_FOURIER_REAL_AND_IMAG 4 -#define PV_FOURIER_PWR 5 -#define PV_FOURIER_PHASE 6 +#define PV_DATA 1 +#define PV_FOURIER_REAL 2 +#define PV_FOURIER_IMAG 3 +#define PV_FOURIER_REAL_AND_IMAG 4 +#define PV_FOURIER_PWR 5 +#define PV_FOURIER_PHASE 6 +#define PV_FOURIER_PHASE_OPT_REAL 7 // Canvas menu id's #define P_MENU_ID_DATA 10001 @@ -71,13 +72,14 @@ #define P_MENU_PLOT_OFFSET 1000 -#define P_MENU_ID_FOURIER_REAL 100 -#define P_MENU_ID_FOURIER_IMAG 101 -#define P_MENU_ID_FOURIER_REAL_AND_IMAG 102 -#define P_MENU_ID_FOURIER_PWR 103 -#define P_MENU_ID_FOURIER_PHASE 104 -#define P_MENU_ID_FOURIER_PHASE_PLUS 105 -#define P_MENU_ID_FOURIER_PHASE_MINUS 106 +#define P_MENU_ID_FOURIER_REAL 100 +#define P_MENU_ID_FOURIER_IMAG 101 +#define P_MENU_ID_FOURIER_REAL_AND_IMAG 102 +#define P_MENU_ID_FOURIER_PWR 103 +#define P_MENU_ID_FOURIER_PHASE 104 +#define P_MENU_ID_FOURIER_PHASE_OPT_REAL 105 +#define P_MENU_ID_FOURIER_PHASE_PLUS 106 +#define P_MENU_ID_FOURIER_PHASE_MINUS 107 //------------------------------------------------------------------------ /** @@ -122,16 +124,19 @@ typedef struct { TH1F *dataFourierIm; ///< imaginary part of the Fourier transform of the data histogram TH1F *dataFourierPwr; ///< power spectrum of the Fourier transform of the data histogram TH1F *dataFourierPhase; ///< phase spectrum of the Fourier transform of the data histogram + TH1F *dataFourierPhaseOptReal; ///< phase optimized real part spectrum Fourier transform of the data histogram TH1F *theory; ///< theory histogram belonging to the data histogram TH1F *theoryFourierRe; ///< real part of the Fourier transform of the theory histogram TH1F *theoryFourierIm; ///< imaginary part of the Fourier transform of the theory histogram TH1F *theoryFourierPwr; ///< power spectrum of the Fourier transform of the theory histogram TH1F *theoryFourierPhase; ///< phase spectrum of the Fourier transform of the theory histogram + TH1F *theoryFourierPhaseOptReal; ///< phase optimized real part spectrum Fourier transform of the theory histogram TH1F *diff; ///< difference histogram, i.e. data-theory TH1F *diffFourierRe; ///< real part of the Fourier transform of the diff histogram TH1F *diffFourierIm; ///< imaginary part of the Fourier transform of the diff histogram TH1F *diffFourierPwr; ///< power spectrum of the Fourier transform of the diff histogram TH1F *diffFourierPhase; ///< phase spectrum of the Fourier transform of the diff histogram + TH1F *diffFourierPhaseOptReal; ///< phase optimized real part spectrum Fourier transform of the diff histogram PMusrCanvasPlotRange *dataRange; ///< keep the msr-file plot data range UInt_t diffFourierTag; ///< 0=not relevant, 1=d-f (Fourier of difference time spectra), 2=f-d (difference of Fourier spectra) } PMusrCanvasDataSet; @@ -302,7 +307,6 @@ class PMusrCanvas : public TObject, public TQObject virtual void HandleDifferenceFourier(); virtual void HandleFourierDifference(); virtual void HandleAverage(); - virtual Double_t FindOptimalFourierPhase(); virtual void CleanupDifference(); virtual void CleanupFourier(); virtual void CleanupFourierDifference();