diff --git a/src/classes/PFourier.cpp b/src/classes/PFourier.cpp index b6034b26..8e1ad423 100644 --- a/src/classes/PFourier.cpp +++ b/src/classes/PFourier.cpp @@ -94,8 +94,8 @@ void PFTPhaseCorrection::Minimize() // create Minuit2 parameters ROOT::Minuit2::MnUserParameters upar; - upar.Add("c0", fPh_c0, 0.5); - upar.Add("c1", fPh_c1, 0.5); + upar.Add("c0", fPh_c0, 2.0); + upar.Add("c1", fPh_c1, 2.0); // create minimizer ROOT::Minuit2::MnMinimize mn_min(*this, upar); @@ -296,7 +296,7 @@ PFourier::PFourier(TH1F *data, Int_t unitTag, Double_t startTime, Double_t endTi fValid = true; fIn = 0; fOut = 0; - fPhCorrectedReFT = 0; +//as fPhCorrectedReFT = 0; fApodization = F_APODIZATION_NONE; @@ -384,8 +384,8 @@ PFourier::~PFourier() fftw_free(fIn); if (fOut) fftw_free(fOut); - if (fPhCorrectedReFT) - delete fPhCorrectedReFT; +//as if (fPhCorrectedReFT) +//as delete fPhCorrectedReFT; } //-------------------------------------------------------------------------- @@ -485,73 +485,79 @@ TH1F* PFourier::GetRealFourier(const Double_t scale) } //-------------------------------------------------------------------------- -// GetPhaseOptRealFourier (public) +// GetPhaseOptRealFourier (public, static) //-------------------------------------------------------------------------- /** *

returns the phase corrected real Fourier transform. * * \return the TH1F histo of the phase 'optimzed' real Fourier transform. * + * \param re real part Fourier histogram + * \param im imaginary part Fourier histogram * \param phase return value of the optimal phase dispersion phase[0]+phase[1]*i/N * \param scale normalisation factor * \param min minimal freq / field from which to optimise. Given in the choosen unit. * \param max maximal freq / field up to which to optimise. Given in the choosen unit. */ -TH1F* PFourier::GetPhaseOptRealFourier(vector &phase, const Double_t scale, const Double_t min, const Double_t max) +TH1F* PFourier::GetPhaseOptRealFourier(const TH1F *re, const TH1F *im, vector &phase, + const Double_t scale, const Double_t min, const Double_t max) { + if ((re == 0) || (im == 0)) + return 0; + phase.resize(2); // c0, c1 - UInt_t noOfFourierBins = 0; - if (fNoOfBins % 2 == 0) - noOfFourierBins = fNoOfBins/2; - else - noOfFourierBins = (fNoOfBins+1)/2; + TAxis *axis = re->GetXaxis(); - UInt_t minBin = 0; - UInt_t maxBin = noOfFourierBins; + Int_t minBin = 1; + Int_t maxBin = axis->GetNbins(); + Int_t noOfBins = axis->GetNbins(); + Double_t res = axis->GetBinWidth(1); // check if minimum frequency is given. If yes, get the proper minBin if (min != -1.0) { - minBin = (UInt_t)(min/fResolution); + minBin = axis->FindFixBin(min); + if ((minBin == 0) || (minBin > maxBin)) { + minBin = 1; + cerr << "**WARNING** minimum frequency/field out of range. Will adopted it." << endl; + } } // check if maximum frequency is given. If yes, get the proper maxBin if (max != -1.0) { - maxBin = (UInt_t)(max/fResolution); - if (maxBin >= noOfFourierBins) { - maxBin = noOfFourierBins; + maxBin = axis->FindFixBin(max); + if ((maxBin == 0) || (maxBin > axis->GetNbins())) { + maxBin = axis->GetNbins(); cerr << "**WARNING** maximum frequency/field out of range. Will adopted it." << endl; } } // copy the real/imag Fourier from min to max vector realF, imagF; - for (UInt_t i=minBin; i<=maxBin; i++) { - realF.push_back(fOut[i][0]); - imagF.push_back(fOut[i][1]); + for (Int_t i=minBin; i<=maxBin; i++) { + realF.push_back(re->GetBinContent(i)); + imagF.push_back(im->GetBinContent(i)); } // optimize real FT phase - fPhCorrectedReFT = new PFTPhaseCorrection(realF, imagF); - if (fPhCorrectedReFT == 0) { - fValid = false; + PFTPhaseCorrection *phCorrectedReFT = new PFTPhaseCorrection(realF, imagF); + if (phCorrectedReFT == 0) { cerr << endl << "**SEVERE ERROR** couldn't invoke PFTPhaseCorrection object." << endl; return 0; } - fPhCorrectedReFT->Minimize(); - if (!fPhCorrectedReFT->IsValid()) { - fValid = false; - cerr << endl << "**ERROR** could not fina a valid phase correction minimum." << endl; + phCorrectedReFT->Minimize(); + if (!phCorrectedReFT->IsValid()) { + cerr << endl << "**ERROR** could not find a valid phase correction minimum." << endl; return 0; } - phase[0] = fPhCorrectedReFT->GetPhaseCorrectionParam(0); - phase[1] = fPhCorrectedReFT->GetPhaseCorrectionParam(1); + phase[0] = phCorrectedReFT->GetPhaseCorrectionParam(0); + phase[1] = phCorrectedReFT->GetPhaseCorrectionParam(1); // clean up - if (fPhCorrectedReFT) { - delete fPhCorrectedReFT; - fPhCorrectedReFT = 0; + if (phCorrectedReFT) { + delete phCorrectedReFT; + phCorrectedReFT = 0; } realF.clear(); imagF.clear(); @@ -559,21 +565,20 @@ TH1F* PFourier::GetPhaseOptRealFourier(vector &phase, const Double_t s // invoke the real phase optimised histo to be filled. Caller is the owner! Char_t name[256]; Char_t title[256]; - snprintf(name, sizeof(name), "%s_Fourier_PhOptRe", fData->GetName()); - snprintf(title, sizeof(title), "%s_Fourier_PhOptRe", fData->GetTitle()); + snprintf(name, sizeof(name), "%s_Fourier_PhOptRe", re->GetName()); + snprintf(title, sizeof(title), "%s_Fourier_PhOptRe", re->GetTitle()); - TH1F *realPhaseOptFourier = new TH1F(name, title, noOfFourierBins, -fResolution/2.0, (Double_t)(noOfFourierBins-1)*fResolution+fResolution/2.0); + TH1F *realPhaseOptFourier = new TH1F(name, title, noOfBins, -res/2.0, (Double_t)(noOfBins-1)*res+res/2.0); if (realPhaseOptFourier == 0) { - fValid = false; cerr << endl << "**SEVERE ERROR** couldn't allocate memory for the real part of the Fourier transform." << endl; return 0; } // fill realFourier vector Double_t ph; - for (UInt_t i=0; iSetBinContent(i+1, scale*(fOut[i][0]*cos(ph) - fOut[i][1]*sin(ph))); + realPhaseOptFourier->SetBinContent(i+1, scale*(re->GetBinContent(i+1)*cos(ph) - im->GetBinContent(i+1)*sin(ph))); realPhaseOptFourier->SetBinError(i+1, 0.0); } diff --git a/src/classes/PFourierCanvas.cpp b/src/classes/PFourierCanvas.cpp index 0b9cbc1f..17e06de2 100644 --- a/src/classes/PFourierCanvas.cpp +++ b/src/classes/PFourierCanvas.cpp @@ -341,6 +341,16 @@ void PFourierCanvas::HandleCmdKey(Int_t event, Int_t x, Int_t y, TObject *select CleanupAverage(); PlotFourier(); } + } else if (x == 'c') { + Int_t state = fFourierPad->GetCrosshair(); + if (state == 0) { + fMainCanvas->ToggleEventStatus(); + fFourierPad->SetCrosshair(2); + } else { + fMainCanvas->ToggleEventStatus(); + fFourierPad->SetCrosshair(0); + } + fMainCanvas->Update(); } else if (x == '+') { // increment phase (Fourier real/imag) IncrementFourierPhase(); } else if (x == '-') { // decrement phase (Fourier real/imag) @@ -366,32 +376,62 @@ void PFourierCanvas::HandleMenuPopup(Int_t id) fPopupFourier->UnCheckEntries(); fPopupFourier->CheckEntry(P_MENU_ID_FOURIER+P_MENU_ID_FOURIER_REAL); fCurrentPlotView = FOURIER_PLOT_REAL; - PlotFourier(); + if (!fAveragedView) { + PlotFourier(); + } else { + HandleAverage(); + PlotAverage(); + } } else if (id == P_MENU_ID_FOURIER+P_MENU_ID_FOURIER_IMAG) { fPopupFourier->UnCheckEntries(); fPopupFourier->CheckEntry(P_MENU_ID_FOURIER+P_MENU_ID_FOURIER_IMAG); fCurrentPlotView = FOURIER_PLOT_IMAG; - PlotFourier(); + if (!fAveragedView) { + PlotFourier(); + } else { + HandleAverage(); + PlotAverage(); + } } else if (id == P_MENU_ID_FOURIER+P_MENU_ID_FOURIER_REAL_AND_IMAG) { fPopupFourier->UnCheckEntries(); fPopupFourier->CheckEntry(P_MENU_ID_FOURIER+P_MENU_ID_FOURIER_REAL_AND_IMAG); fCurrentPlotView = P_MENU_ID_FOURIER_REAL_AND_IMAG; - PlotFourier(); + if (!fAveragedView) { + PlotFourier(); + } else { + HandleAverage(); + PlotAverage(); + } } else if (id == P_MENU_ID_FOURIER+P_MENU_ID_FOURIER_PWR) { fPopupFourier->UnCheckEntries(); fPopupFourier->CheckEntry(P_MENU_ID_FOURIER+P_MENU_ID_FOURIER_PWR); fCurrentPlotView = FOURIER_PLOT_POWER; - PlotFourier(); + if (!fAveragedView) { + PlotFourier(); + } else { + HandleAverage(); + PlotAverage(); + } } else if (id == P_MENU_ID_FOURIER+P_MENU_ID_FOURIER_PHASE) { fPopupFourier->UnCheckEntries(); fPopupFourier->CheckEntry(P_MENU_ID_FOURIER+P_MENU_ID_FOURIER_PHASE); fCurrentPlotView = FOURIER_PLOT_PHASE; - PlotFourier(); + if (!fAveragedView) { + PlotFourier(); + } else { + HandleAverage(); + PlotAverage(); + } } else if (id == P_MENU_ID_FOURIER+P_MENU_ID_FOURIER_PHASE_OPT_REAL) { fPopupFourier->UnCheckEntries(); fPopupFourier->CheckEntry(P_MENU_ID_FOURIER+P_MENU_ID_FOURIER_PHASE_OPT_REAL); fCurrentPlotView = FOURIER_PLOT_PHASE_OPT_REAL; - PlotFourier(); + if (!fAveragedView) { + PlotFourier(); + } else { + HandleAverage(); + PlotAverage(); + } } else if (id == P_MENU_ID_FOURIER+P_MENU_ID_FOURIER_PHASE_PLUS) { IncrementFourierPhase(); } else if (id == P_MENU_ID_FOURIER+P_MENU_ID_FOURIER_PHASE_MINUS) { @@ -814,7 +854,10 @@ void PFourierCanvas::InitFourierDataSets() fFourierHistos[i].dataFourierIm = fFourier[i]->GetImaginaryFourier(); fFourierHistos[i].dataFourierPwr = fFourier[i]->GetPowerFourier(); fFourierHistos[i].dataFourierPhase = fFourier[i]->GetPhaseFourier(); - fFourierHistos[i].dataFourierPhaseOptReal = fFourier[i]->GetPhaseOptRealFourier(fFourierHistos[i].optPhase, 1.0, fInitialXRange[0], fInitialXRange[1]); + if (fCurrentPlotView == FOURIER_PLOT_PHASE_OPT_REAL) { + fFourierHistos[i].dataFourierPhaseOptReal = fFourier[i]->GetPhaseOptRealFourier(fFourierHistos[i].dataFourierRe, + fFourierHistos[i].dataFourierIm, fFourierHistos[i].optPhase, 1.0, fInitialXRange[0], fInitialXRange[1]); + } } // rescale histo to abs(maximum) == 1 @@ -879,17 +922,19 @@ void PFourierCanvas::InitFourierDataSets() } } - // phase opt real - for (UInt_t i=0; iGetBinContent(j); - if (fabs(dval) > max) - max = dval; + if (fCurrentPlotView == FOURIER_PLOT_PHASE_OPT_REAL) { + // phase opt real + for (UInt_t i=0; iGetBinContent(j); + if (fabs(dval) > max) + max = dval; + } } - } - for (UInt_t i=0; iGetNbinsX(); j++) { - fFourierHistos[i].dataFourierPhaseOptReal->SetBinContent(j, fFourierHistos[i].dataFourierPhaseOptReal->GetBinContent(j)/fabs(max)); + for (UInt_t i=0; iGetNbinsX(); j++) { + fFourierHistos[i].dataFourierPhaseOptReal->SetBinContent(j, fFourierHistos[i].dataFourierPhaseOptReal->GetBinContent(j)/fabs(max)); + } } } @@ -903,8 +948,10 @@ void PFourierCanvas::InitFourierDataSets() fFourierHistos[i].dataFourierPwr->SetLineColor(fColorList[i]); fFourierHistos[i].dataFourierPhase->SetMarkerColor(fColorList[i]); fFourierHistos[i].dataFourierPhase->SetLineColor(fColorList[i]); - fFourierHistos[i].dataFourierPhaseOptReal->SetMarkerColor(fColorList[i]); - fFourierHistos[i].dataFourierPhaseOptReal->SetLineColor(fColorList[i]); + if (fCurrentPlotView == FOURIER_PLOT_PHASE_OPT_REAL) { + fFourierHistos[i].dataFourierPhaseOptReal->SetMarkerColor(fColorList[i]); + fFourierHistos[i].dataFourierPhaseOptReal->SetLineColor(fColorList[i]); + } } // set the marker symbol and size @@ -917,8 +964,10 @@ void PFourierCanvas::InitFourierDataSets() fFourierHistos[i].dataFourierPwr->SetMarkerSize(0.7); fFourierHistos[i].dataFourierPhase->SetMarkerStyle(fMarkerList[i]); fFourierHistos[i].dataFourierPhase->SetMarkerSize(0.7); - fFourierHistos[i].dataFourierPhaseOptReal->SetMarkerStyle(fMarkerList[i]); - fFourierHistos[i].dataFourierPhaseOptReal->SetMarkerSize(0.7); + if (fCurrentPlotView == FOURIER_PLOT_PHASE_OPT_REAL) { + fFourierHistos[i].dataFourierPhaseOptReal->SetMarkerStyle(fMarkerList[i]); + fFourierHistos[i].dataFourierPhaseOptReal->SetMarkerSize(0.7); + } } // initialize average histos @@ -1116,6 +1165,10 @@ void PFourierCanvas::HandleAverage() TString name(""); Double_t dval=0.0; + Bool_t phaseOptRealPresent = false; + if (fFourierHistos[0].dataFourierPhaseOptReal != 0) + phaseOptRealPresent = true; + // check if ALL data shall be averaged if (fAveragedView) { fFourierAverage.resize(1); @@ -1141,10 +1194,12 @@ void PFourierCanvas::HandleAverage() fFourierHistos[0].dataFourierPhase->GetXaxis()->GetXmin(), fFourierHistos[0].dataFourierPhase->GetXaxis()->GetXmax()); - name = TString(fFourierHistos[0].dataFourierPhaseOptReal->GetTitle()) + "_avg"; - fFourierAverage[0].dataFourierPhaseOptReal = new TH1F(name, name, fFourierHistos[0].dataFourierPhaseOptReal->GetNbinsX(), - fFourierHistos[0].dataFourierPhaseOptReal->GetXaxis()->GetXmin(), - fFourierHistos[0].dataFourierPhaseOptReal->GetXaxis()->GetXmax()); + if (phaseOptRealPresent) { + name = TString(fFourierHistos[0].dataFourierPhaseOptReal->GetTitle()) + "_avg"; + fFourierAverage[0].dataFourierPhaseOptReal = new TH1F(name, name, fFourierHistos[0].dataFourierPhaseOptReal->GetNbinsX(), + fFourierHistos[0].dataFourierPhaseOptReal->GetXaxis()->GetXmin(), + fFourierHistos[0].dataFourierPhaseOptReal->GetXaxis()->GetXmax()); + } // real average for (Int_t j=0; jGetNbinsX(); j++) { @@ -1206,20 +1261,22 @@ void PFourierCanvas::HandleAverage() fFourierAverage[0].dataFourierPhase->SetMarkerSize(0.8); fFourierAverage[0].dataFourierPhase->SetMarkerStyle(22); - // phase optimised real average - for (Int_t j=0; jGetNbinsX(); j++) { - dval = 0.0; - for (UInt_t i=0; iGetNbinsX()) - dval += GetInterpolatedValue(fFourierHistos[i].dataFourierPhaseOptReal, fFourierHistos[0].dataFourierPhaseOptReal->GetBinCenter(j)); + if (phaseOptRealPresent) { + // phase optimised real average + for (Int_t j=0; jGetNbinsX(); j++) { + dval = 0.0; + for (UInt_t i=0; iGetNbinsX()) + dval += GetInterpolatedValue(fFourierHistos[i].dataFourierPhaseOptReal, fFourierHistos[0].dataFourierPhaseOptReal->GetBinCenter(j)); + } + fFourierAverage[0].dataFourierPhaseOptReal->SetBinContent(j, dval/fFourierHistos.size()); } - fFourierAverage[0].dataFourierPhaseOptReal->SetBinContent(j, dval/fFourierHistos.size()); + // set marker color, line color, maker size, marker type + fFourierAverage[0].dataFourierPhaseOptReal->SetMarkerColor(kBlack); + fFourierAverage[0].dataFourierPhaseOptReal->SetLineColor(kBlack); + fFourierAverage[0].dataFourierPhaseOptReal->SetMarkerSize(0.8); + fFourierAverage[0].dataFourierPhaseOptReal->SetMarkerStyle(22); } - // set marker color, line color, maker size, marker type - fFourierAverage[0].dataFourierPhaseOptReal->SetMarkerColor(kBlack); - fFourierAverage[0].dataFourierPhaseOptReal->SetLineColor(kBlack); - fFourierAverage[0].dataFourierPhaseOptReal->SetMarkerSize(0.8); - fFourierAverage[0].dataFourierPhaseOptReal->SetMarkerStyle(22); } // check if per data set shall be averaged @@ -1283,10 +1340,12 @@ void PFourierCanvas::HandleAverage() fFourierHistos[0].dataFourierPhase->GetXaxis()->GetXmin(), fFourierHistos[0].dataFourierPhase->GetXaxis()->GetXmax()); - name = TString(fFourierHistos[start].dataFourierPhaseOptReal->GetTitle()) + "_avg"; - fFourierAverage[i].dataFourierPhaseOptReal = new TH1F(name, name, fFourierHistos[0].dataFourierPhaseOptReal->GetNbinsX(), - fFourierHistos[0].dataFourierPhaseOptReal->GetXaxis()->GetXmin(), - fFourierHistos[0].dataFourierPhaseOptReal->GetXaxis()->GetXmax()); + if (phaseOptRealPresent) { + name = TString(fFourierHistos[start].dataFourierPhaseOptReal->GetTitle()) + "_avg"; + fFourierAverage[i].dataFourierPhaseOptReal = new TH1F(name, name, fFourierHistos[0].dataFourierPhaseOptReal->GetNbinsX(), + fFourierHistos[0].dataFourierPhaseOptReal->GetXaxis()->GetXmin(), + fFourierHistos[0].dataFourierPhaseOptReal->GetXaxis()->GetXmax()); + } // real average for (Int_t j=0; jGetNbinsX(); j++) { @@ -1348,24 +1407,61 @@ void PFourierCanvas::HandleAverage() fFourierAverage[i].dataFourierPhase->SetMarkerSize(0.8); fFourierAverage[i].dataFourierPhase->SetMarkerStyle(22); // closed up triangle - // phase optimised real average - for (Int_t j=0; jGetNbinsX(); j++) { - dval = 0.0; - for (Int_t k=start; k<=end; k++) { - if (j < fFourierHistos[k].dataFourierPhaseOptReal->GetNbinsX()) - dval += GetInterpolatedValue(fFourierHistos[k].dataFourierPhaseOptReal, fFourierHistos[0].dataFourierPhaseOptReal->GetBinCenter(j)); + if (phaseOptRealPresent) { + // phase optimised real average + for (Int_t j=0; jGetNbinsX(); j++) { + dval = 0.0; + for (Int_t k=start; k<=end; k++) { + if (j < fFourierHistos[k].dataFourierPhaseOptReal->GetNbinsX()) + dval += GetInterpolatedValue(fFourierHistos[k].dataFourierPhaseOptReal, fFourierHistos[0].dataFourierPhaseOptReal->GetBinCenter(j)); + } + fFourierAverage[i].dataFourierPhaseOptReal->SetBinContent(j, dval/(end-start+1)); } - fFourierAverage[i].dataFourierPhaseOptReal->SetBinContent(j, dval/(end-start+1)); + // set marker color, line color, maker size, marker type + fFourierAverage[i].dataFourierPhaseOptReal->SetMarkerColor(fColorList[i]); + fFourierAverage[i].dataFourierPhaseOptReal->SetLineColor(fColorList[i]); + fFourierAverage[i].dataFourierPhaseOptReal->SetMarkerSize(0.8); + fFourierAverage[i].dataFourierPhaseOptReal->SetMarkerStyle(22); // closed up triangle } - // set marker color, line color, maker size, marker type - fFourierAverage[i].dataFourierPhaseOptReal->SetMarkerColor(fColorList[i]); - fFourierAverage[i].dataFourierPhaseOptReal->SetLineColor(fColorList[i]); - fFourierAverage[i].dataFourierPhaseOptReal->SetMarkerSize(0.8); - fFourierAverage[i].dataFourierPhaseOptReal->SetMarkerStyle(22); // closed up triangle } } } +//-------------------------------------------------------------------------- +// CalcPhaseOptReal (private) +//-------------------------------------------------------------------------- +/** + *

calculate the phase opt. real FT + */ +void PFourierCanvas::CalcPhaseOptReal() +{ + Int_t start = fFourierHistos[0].dataFourierRe->FindFixBin(fInitialXRange[0]); + Int_t end = fFourierHistos[0].dataFourierRe->FindFixBin(fInitialXRange[1]); + + Double_t dval=0.0, max=0.0; + for (UInt_t i=0; iGetPhaseOptRealFourier(fFourierHistos[i].dataFourierRe, + fFourierHistos[i].dataFourierIm, fFourierHistos[i].optPhase, 1.0, fInitialXRange[0], fInitialXRange[1]); + // normalize it + max = 0.0; + for (Int_t j=start; j<=end; j++) { + dval = fFourierHistos[i].dataFourierPhaseOptReal->GetBinContent(j); + if (fabs(dval) > max) + max = dval; + } + for (Int_t j=1; jGetNbinsX(); j++) { + fFourierHistos[i].dataFourierPhaseOptReal->SetBinContent(j, fFourierHistos[i].dataFourierPhaseOptReal->GetBinContent(j)/fabs(max)); + } + // set the marker and line color + fFourierHistos[i].dataFourierPhaseOptReal->SetMarkerColor(fColorList[i]); + fFourierHistos[i].dataFourierPhaseOptReal->SetLineColor(fColorList[i]); + + // set the marker symbol and size + fFourierHistos[i].dataFourierPhaseOptReal->SetMarkerStyle(fMarkerList[i]); + fFourierHistos[i].dataFourierPhaseOptReal->SetMarkerSize(0.7); + } +} + //-------------------------------------------------------------------------- // PlotFourier (private) //-------------------------------------------------------------------------- @@ -1374,6 +1470,14 @@ void PFourierCanvas::HandleAverage() */ void PFourierCanvas::PlotFourier() { + // check if phase opt real Fourier spectra already exists if requested, + // and if not calculate them first + if (fCurrentPlotView == FOURIER_PLOT_PHASE_OPT_REAL) { + if (fFourierHistos[0].dataFourierPhaseOptReal == 0) { // not yet calculated + CalcPhaseOptReal(); + } + } + fFourierPad->cd(); Double_t xmin=0, xmax=0; @@ -1669,6 +1773,11 @@ void PFourierCanvas::PlotAverage() fFourierAverage[0].dataFourierPhase->Draw("p"); break; case FOURIER_PLOT_PHASE_OPT_REAL: + if (fFourierHistos[0].dataFourierPhaseOptReal == 0) { + cout << "debug> need to calculate phase opt. average first ..." << endl; + CalcPhaseOptReal(); + HandleAverage(); + } fFourierAverage[0].dataFourierPhaseOptReal->GetXaxis()->SetRangeUser(xmin, xmax); ymin = GetMinimum(fFourierAverage[0].dataFourierPhaseOptReal, xmin, xmax); ymax = GetMaximum(fFourierAverage[0].dataFourierPhaseOptReal, xmin, xmax); diff --git a/src/classes/PMusrCanvas.cpp b/src/classes/PMusrCanvas.cpp index 2cae4d9e..e335fac5 100644 --- a/src/classes/PMusrCanvas.cpp +++ b/src/classes/PMusrCanvas.cpp @@ -182,11 +182,14 @@ PMusrCanvas::PMusrCanvas() * \param wh height (in pixels) of the canvas. * \param batch flag: if set true, the canvas will not be displayed. This is used when just dumping of a * graphical output file is wished. + * \param fourier flag: if set true, the canvas will present the Fourier view. + * \param avg flag: if set true, the canvas will present the averages data/Fourier view. */ PMusrCanvas::PMusrCanvas(const Int_t number, const Char_t* title, Int_t wtopx, Int_t wtopy, Int_t ww, Int_t wh, - const Bool_t batch, const Bool_t fourier) : - fStartWithFourier(fourier), fBatchMode(batch), fPlotNumber(number) + const Bool_t batch, const Bool_t fourier, const Bool_t avg) : + fStartWithFourier(fourier), fStartWithAvg(avg), + fBatchMode(batch), fPlotNumber(number) { fTimeout = 0; fTimeoutTimer = 0; @@ -235,15 +238,17 @@ PMusrCanvas::PMusrCanvas(const Int_t number, const Char_t* title, * \param colorList pre-defined list of colors * \param batch flag: if set true, the canvas will not be displayed. This is used when just dumping of a * graphical output file is wished. + * \param fourier flag: if set true, the canvas will present the Fourier view. + * \param avg flag: if set true, the canvas will present the averages data/Fourier view. */ PMusrCanvas::PMusrCanvas(const Int_t number, const Char_t* title, Int_t wtopx, Int_t wtopy, Int_t ww, Int_t wh, PMsrFourierStructure fourierDefault, const PIntVector markerList, const PIntVector colorList, - const Bool_t batch, const Bool_t fourier) : - fStartWithFourier(fourier), fBatchMode(batch), - fPlotNumber(number), fFourier(fourierDefault), - fMarkerList(markerList), fColorList(colorList) + const Bool_t batch, const Bool_t fourier, const Bool_t avg) : + fStartWithFourier(fourier), fStartWithAvg(avg), fBatchMode(batch), + fPlotNumber(number), fFourier(fourierDefault), + fMarkerList(markerList), fColorList(colorList) { fTimeout = 0; fTimeoutTimer = 0; @@ -832,6 +837,12 @@ void PMusrCanvas::UpdateDataTheoryPad() HandleFourier(); PlotFourier(); } + + // if fStartWithAvg=true, start with averaged data/Fourier representation + // fStartWithAvg is given at the command line level + if (fStartWithAvg) { + HandleCmdKey(kKeyPress, (Int_t)'a', 0, 0); + } } //-------------------------------------------------------------------------- @@ -1357,6 +1368,13 @@ void PMusrCanvas::HandleMenuPopup(Int_t id) // set appropriate plot view fPreviousPlotView = fCurrentPlotView; fCurrentPlotView = PV_FOURIER_PHASE_OPT_REAL; + // make sure that phase opt. real indeed exists + if (fData[0].dataFourierPhaseOptReal == 0) { + if (fData[0].dataFourierRe == 0) + HandleFourier(); + else + CalcPhaseOptReFT(); + } // uncheck data fPopupMain->UnCheckEntry(P_MENU_ID_DATA+P_MENU_PLOT_OFFSET*fPlotNumber); // check appropriate fourier popup item @@ -3331,10 +3349,7 @@ void PMusrCanvas::HandleDifference() */ void PMusrCanvas::HandleFourier() { - PDoubleVector phaseParam; - Double_t min=-1.0, max=-1.0; - Double_t re, im, ph; - char hName[1024]; + Double_t re, im; // check if plot type is appropriate for fourier if (fPlotType == MSR_PLOT_NON_MUSR) @@ -3347,9 +3362,9 @@ void PMusrCanvas::HandleFourier() double endTime = fXmax; if (!fStartWithFourier) { // fHistoFrame present, hence get start/end from it bin = fHistoFrame->GetXaxis()->GetFirst(); - startTime = fHistoFrame->GetBinCenter(bin); + startTime = fHistoFrame->GetBinLowEdge(bin); bin = fHistoFrame->GetXaxis()->GetLast(); - endTime = fHistoFrame->GetBinCenter(bin); + endTime = fHistoFrame->GetBinLowEdge(bin)+fHistoFrame->GetBinWidth(bin); } for (UInt_t i=0; iGetMsrFourierList()->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()); @@ -3385,22 +3394,18 @@ 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; @@ -3418,31 +3423,18 @@ void PMusrCanvas::HandleFourier() // get power part of the data fData[i].theoryFourierPwr = fourierTheory.GetPowerFourier(scale); // 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); - } + fData[i].theoryFourierPhase = fourierTheory.GetPhaseFourier(); // 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()); + } + + // phase opt. real FT requested initially in the msr-file, hence calculate it here + if (fCurrentPlotView == PV_FOURIER_PHASE_OPT_REAL) { + CalcPhaseOptReFT(); } // apply global phase if present @@ -3475,42 +3467,6 @@ 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)) { - - fCurrentFourierPhase = FindOptimalFourierPhase(); - - // apply optimal Fourier phase - double re, im; - const double cp = TMath::Cos(fCurrentFourierPhase/180.0*TMath::Pi()); - const double sp = TMath::Sin(fCurrentFourierPhase/180.0*TMath::Pi()); - - for (UInt_t i=0; iGetNbinsX(); j++) { // loop over a fourier data set - // calculate new fourier data set value - re = fData[i].dataFourierRe->GetBinContent(j) * cp + fData[i].dataFourierIm->GetBinContent(j) * sp; - im = fData[i].dataFourierIm->GetBinContent(j) * cp - fData[i].dataFourierRe->GetBinContent(j) * sp; - // overwrite fourier data set value - fData[i].dataFourierRe->SetBinContent(j, re); - fData[i].dataFourierIm->SetBinContent(j, im); - } - } - if ((fData[i].theoryFourierRe != 0) && (fData[i].theoryFourierIm != 0)) { - for (Int_t j=0; jGetNbinsX(); j++) { // loop over a fourier data set - // calculate new fourier data set value - re = fData[i].theoryFourierRe->GetBinContent(j) * cp + fData[i].theoryFourierIm->GetBinContent(j) * sp; - im = fData[i].theoryFourierIm->GetBinContent(j) * cp - fData[i].theoryFourierRe->GetBinContent(j) * sp; - // overwrite fourier data set value - fData[i].theoryFourierRe->SetBinContent(j, re); - fData[i].theoryFourierIm->SetBinContent(j, im); - } - } - } - } -*/ } } @@ -4190,6 +4146,10 @@ void PMusrCanvas::CleanupFourierDifference() delete fData[i].diffFourierPhase; fData[i].diffFourierPhase = 0; } + if (fData[i].diffFourierPhaseOptReal != 0) { + delete fData[i].diffFourierPhaseOptReal; + fData[i].diffFourierPhaseOptReal = 0; + } } } @@ -4275,6 +4235,64 @@ void PMusrCanvas::CleanupAverage() } } +//-------------------------------------------------------------------------- +// CalculateDiff (private) +//-------------------------------------------------------------------------- +/** + * @brief PMusrCanvas::CalcPhaseOptReFT + */ +void PMusrCanvas::CalcPhaseOptReFT() +{ + Double_t min = fMsrHandler->GetMsrFourierList()->fRangeForPhaseCorrection[0]; + Double_t max = fMsrHandler->GetMsrFourierList()->fRangeForPhaseCorrection[1]; + + if ((min == -1.0) && (max == -1.0)) { + if ((fFourier.fPlotRange[0] != -1) && (fFourier.fPlotRange[1] != -1)) { + min = fFourier.fPlotRange[0]; + max = fFourier.fPlotRange[1]; + } else { + min = fData[0].dataFourierRe->GetBinLowEdge(1); + max = fData[0].dataFourierRe->GetBinLowEdge(fData[0].dataFourierRe->GetNbinsX())+fData[0].dataFourierRe->GetBinWidth(1); + } + } + + PDoubleVector phaseParam; + Char_t hName[1024]; + Double_t ph, re; + + for (UInt_t i=0; iSetMarkerColor(fData[i].data->GetMarkerColor()); + fData[i].dataFourierPhaseOptReal->SetLineColor(fData[i].data->GetLineColor()); + // set marker size + fData[i].dataFourierPhaseOptReal->SetMarkerSize(1); + // set marker type + fData[i].dataFourierPhaseOptReal->SetMarkerStyle(fData[i].data->GetMarkerStyle()); + + // handle Fourier theory part + // 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()->FindFixBin(min); + Int_t maxBin = fData[i].theoryFourierPhaseOptReal->GetXaxis()->FindFixBin(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].theoryFourierPhaseOptReal->SetLineColor(fData[i].theory->GetLineColor()); + } +} + //-------------------------------------------------------------------------- // CalculateDiff (private) //-------------------------------------------------------------------------- @@ -5924,6 +5942,58 @@ void PMusrCanvas::PlotFourierDifference(Bool_t unzoom) PlotFourierPhaseValue(); + 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].diffFourierPhaseOptReal->GetBinLowEdge(1); + xmax = fData[0].diffFourierPhaseOptReal->GetBinLowEdge(fData[0].diffFourierPhaseOptReal->GetNbinsX())+fData[0].diffFourierPhaseOptReal->GetBinWidth(1); + } + + // set y-range + // first find minimum/maximum of all histos + ymin = GetMinimum(fData[0].diffFourierPhaseOptReal); + ymax = GetMaximum(fData[0].diffFourierPhaseOptReal); + for (UInt_t i=1; i 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); + + // set ranges for phase opt. real Fourier difference + for (UInt_t i=0; iGetXaxis()->SetRangeUser(xmin, xmax); + fData[i].diffFourierPhaseOptReal->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); + if (fData[0].diffFourierTag == 1) + fHistoFrame->GetYaxis()->SetTitle("Real Fourier (d-f: data-theory)"); + else + fHistoFrame->GetYaxis()->SetTitle("Real Fourier (f-d: [(F data)-(F theory)]"); + + // plot data + for (UInt_t i=0; iDraw("plsame"); + } break; default: break; diff --git a/src/include/PFourier.h b/src/include/PFourier.h index 38f7e7ee..5edaddb9 100644 --- a/src/include/PFourier.h +++ b/src/include/PFourier.h @@ -108,11 +108,14 @@ class PFourier virtual Double_t GetResolution() { return fResolution; } virtual Double_t GetMaxFreq(); virtual TH1F* GetRealFourier(const Double_t scale = 1.0); - virtual TH1F* GetPhaseOptRealFourier(vector &phase, const Double_t scale = 1.0, const Double_t min = -1.0, const Double_t max = -1.0); +//as virtual TH1F* GetPhaseOptRealFourier(vector &phase, const Double_t scale = 1.0, const Double_t min = -1.0, const Double_t max = -1.0); virtual TH1F* GetImaginaryFourier(const Double_t scale = 1.0); virtual TH1F* GetPowerFourier(const Double_t scale = 1.0); virtual TH1F* GetPhaseFourier(const Double_t scale = 1.0); + static TH1F* GetPhaseOptRealFourier(const TH1F *re, const TH1F *im, vector &phase, + const Double_t scale = 1.0, const Double_t min = -1.0, const Double_t max = -1.0); + virtual Bool_t IsValid() { return fValid; } private: @@ -136,7 +139,7 @@ class PFourier fftw_complex *fIn; ///< real part of the Fourier transform fftw_complex *fOut; ///< imaginary part of the Fourier transform - PFTPhaseCorrection *fPhCorrectedReFT; +//as PFTPhaseCorrection *fPhCorrectedReFT; virtual void PrepareFFTwInputData(UInt_t apodizationTag); virtual void ApodizeData(Int_t apodizationTag); diff --git a/src/include/PFourierCanvas.h b/src/include/PFourierCanvas.h index 3aff49e9..05dad1e9 100644 --- a/src/include/PFourierCanvas.h +++ b/src/include/PFourierCanvas.h @@ -159,6 +159,7 @@ class PFourierCanvas : public TObject, public TQObject virtual void InitFourierCanvas(const Char_t* title, Int_t wtopx, Int_t wtopy, Int_t ww, Int_t wh); virtual void CleanupAverage(); virtual void HandleAverage(); + virtual void CalcPhaseOptReal(); virtual void PlotFourier(); virtual void PlotFourierPhaseValue(); diff --git a/src/musrFT.cpp b/src/musrFT.cpp index 84bc5991..85ca609f 100644 --- a/src/musrFT.cpp +++ b/src/musrFT.cpp @@ -713,7 +713,7 @@ Int_t musrFT_dumpData(TString fln, vector &fourierData, Double_t star musrFT_cleanup(hRe); for (UInt_t i=1; iGetRealFourier(); - if (hRe->GetNbinsX()-1 < minSize) + if (hRe->GetNbinsX()-1 < (Int_t)minSize) minSize = hRe->GetNbinsX()-1; musrFT_cleanup(hRe); } @@ -721,7 +721,7 @@ Int_t musrFT_dumpData(TString fln, vector &fourierData, Double_t star for (UInt_t i=0; iGetRealFourier(); hIm = fourierData[i]->GetImaginaryFourier(); - for (Int_t j=1; jGetBinCenter(j); if ((dval >= start) && (dval <= end)) { freq.push_back(dval); @@ -1020,7 +1020,6 @@ Int_t main(Int_t argc, Char_t *argv[]) PStartupOptions startup_options; startup_options.writeExpectedChisq = false; startup_options.estimateN0 = true; - startup_options.alphaEstimateN0 = 0.0; TSAXParser *saxParser = new TSAXParser(); PStartupHandler *startupHandler = new PStartupHandler(); if (!startupHandler->StartupFileFound()) { @@ -1057,7 +1056,6 @@ Int_t main(Int_t argc, Char_t *argv[]) } } } - startupHandler->SetStartupOptions(startup_options); // defines the raw time-domain data vector PPrepFourier data(startupParam.packing, startupParam.bkg_range, startupParam.bkg); @@ -1066,7 +1064,7 @@ Int_t main(Int_t argc, Char_t *argv[]) vector msrHandler; msrHandler.resize(startupParam.msrFln.size()); for (UInt_t i=0; iGetStartupOptions(), true); + msrHandler[i] = new PMsrHandler(startupParam.msrFln[i].Data(), &startup_options, true); status = msrHandler[i]->ReadMsrFile(); if (status != PMUSR_SUCCESS) { switch (status) { diff --git a/src/musrview.cpp b/src/musrview.cpp index f91dd6fd..c416b71d 100644 --- a/src/musrview.cpp +++ b/src/musrview.cpp @@ -63,6 +63,7 @@ void musrview_syntax() cout << endl << " --help : display this help and exit."; cout << endl << " --version : output version information and exit."; cout << endl << " -f, --fourier: will directly present the Fourier transform of the ."; + cout << endl << " -a, --avg: will directly present the averaged data/Fourier of the ."; cout << endl << " --: "; cout << endl << " will produce a graphics-output-file without starting a root session."; cout << endl << " the name is based on the , e.g. 3310.msr -> 3310_0.png"; @@ -102,6 +103,7 @@ int main(int argc, char *argv[]) bool success = true; char fileName[128]; bool fourier = false; + bool avg = false; bool graphicsOutput = false; bool asciiOutput = false; char graphicsExtension[128]; @@ -135,6 +137,8 @@ int main(int argc, char *argv[]) break; } else if (!strcmp(argv[i], "-f") || !strcmp(argv[i], "--fourier")) { fourier = true; + } else if (!strcmp(argv[i], "-a") || !strcmp(argv[i], "--avg")) { + avg = true; } else if (!strcmp(argv[i], "--eps") || !strcmp(argv[i], "--pdf") || !strcmp(argv[i], "--gif") || !strcmp(argv[i], "--jpg") || !strcmp(argv[i], "--png") || !strcmp(argv[i], "--svg") || !strcmp(argv[i], "--xpm") || !strcmp(argv[i], "--root")) { @@ -308,12 +312,12 @@ int main(int argc, char *argv[]) startupHandler->GetMarkerList(), startupHandler->GetColorList(), graphicsOutput||asciiOutput, - fourier); + fourier, avg); else musrCanvas = new PMusrCanvas(i, msrHandler->GetMsrTitle()->Data(), 10+i*100, 10+i*100, 800, 600, graphicsOutput||asciiOutput, - fourier); + fourier, avg); if (!musrCanvas->IsValid()) { cerr << endl << ">> musrview **SEVERE ERROR** Couldn't invoke all necessary objects, will quit.";