From 097a6db6b0a9f0a0411b194e76cb9dc5a8d29530 Mon Sep 17 00:00:00 2001 From: Suter Andreas Date: Thu, 7 Jan 2016 13:19:00 +0100 Subject: [PATCH] more work towards fully functional RRF single histo fit --- src/classes/PMsrHandler.cpp | 23 +- src/classes/PMusr.cpp | 4 - src/classes/PMusrCanvas.cpp | 180 ++++++---- src/classes/PRunListCollection.cpp | 20 +- src/classes/PRunSingleHistoRRF.cpp | 540 +++++------------------------ src/include/PMusr.h | 9 +- src/include/PRunSingleHistoRRF.h | 3 +- 7 files changed, 227 insertions(+), 552 deletions(-) diff --git a/src/classes/PMsrHandler.cpp b/src/classes/PMsrHandler.cpp index 3a3148d0..e1a19754 100644 --- a/src/classes/PMsrHandler.cpp +++ b/src/classes/PMsrHandler.cpp @@ -636,12 +636,8 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t messages) fout.width(16); fout << left << "rrf_freq "; fout.width(8); - -cout << "debug> PMsrHandler::WriteMsrLogFile(): fGlobal.GetRRFFreq(fGlobal.GetRRFUnit().Data())=" << fGlobal.GetRRFFreq(fGlobal.GetRRFUnit().Data()) << endl; -neededPrec = LastSignificant(fGlobal.GetRRFFreq(fGlobal.GetRRFUnit().Data()),10); -cout << "debug> PMsrHandler::WriteMsrLogFile(): neededPrec=" << neededPrec << endl; -fout.precision(neededPrec); - + neededPrec = LastSignificant(fGlobal.GetRRFFreq(fGlobal.GetRRFUnit().Data()),10); + fout.precision(neededPrec); fout << left << fixed << fGlobal.GetRRFFreq(fGlobal.GetRRFUnit().Data()); fout << " " << fGlobal.GetRRFUnit(); fout << endl; @@ -1140,6 +1136,9 @@ fout.precision(neededPrec); case MSR_PLOT_SINGLE_HISTO: fout << "PLOT " << fPlots[plotNo].fPlotType << " (single histo plot)" << endl; break; + case MSR_PLOT_SINGLE_HISTO_RRF: + fout << "PLOT " << fPlots[plotNo].fPlotType << " (single histo RRF plot)" << endl; + break; case MSR_PLOT_ASYM: fout << "PLOT " << fPlots[plotNo].fPlotType << " (asymmetry plot)" << endl; break; @@ -1656,7 +1655,6 @@ Int_t PMsrHandler::WriteMsrFile(const Char_t *filename, map *co fout.width(16); fout << left << "rrf_freq "; fout.width(8); -cout << "debug> PMsrHandler::WriteMsrFile(): fGlobal.GetRRFFreq(fGlobal.GetRRFUnit().Data())=" << fGlobal.GetRRFFreq(fGlobal.GetRRFUnit().Data()) << endl; fout << left << fGlobal.GetRRFFreq(fGlobal.GetRRFUnit().Data()); fout << " " << fGlobal.GetRRFUnit(); fout << endl; @@ -2166,6 +2164,9 @@ cout << "debug> PMsrHandler::WriteMsrFile(): fGlobal.GetRRFFreq(fGlobal.GetRRFUn case MSR_PLOT_SINGLE_HISTO: fout << "PLOT " << fPlots[i].fPlotType << " (single histo plot)" << endl; break; + case MSR_PLOT_SINGLE_HISTO_RRF: + fout << "PLOT " << fPlots[i].fPlotType << " (single histo RRF plot)" << endl; + break; case MSR_PLOT_ASYM: fout << "PLOT " << fPlots[i].fPlotType << " (asymmetry plot)" << endl; break; @@ -4136,6 +4137,7 @@ Bool_t PMsrHandler::HandlePlotEntry(PMsrLines &lines) error = true; break; case MSR_PLOT_SINGLE_HISTO: // like: runs 1 5 13 + case MSR_PLOT_SINGLE_HISTO_RRF: case MSR_PLOT_ASYM: case MSR_PLOT_NON_MUSR: case MSR_PLOT_MU_MINUS: @@ -5789,6 +5791,11 @@ Bool_t PMsrHandler::CheckRRFSettings() } } + // if not a RRF fit, done at this point + if (fittype != MSR_FITTYPE_SINGLE_HISTO_RRF) { + return true; + } + // second set of tests: if RRF fit is chosen, do I find the necessary RRF parameters? fittype = fGlobal.GetFitType(); if (fittype == MSR_FITTYPE_SINGLE_HISTO_RRF) { // make sure RRF freq and RRF packing are set @@ -5802,7 +5809,7 @@ Bool_t PMsrHandler::CheckRRFSettings() cerr << endl << ">> no RRF packing found in the GLOBAL section! Fix it."; return false; } - } else { + } else { // check single runs for RRF UInt_t rrfFitCounter = 0; for (UInt_t i=0; i PMsrGlobalBlock::GetRRFFreq(" << unit << "): unitTag=" << unitTag << ", fRRFFreq=" << fRRFFreq; - // calc the conversion factor if (unitTag == fRRFUnitTag) freq = fRRFFreq; @@ -770,8 +768,6 @@ cout << endl << "debug> PMsrGlobalBlock::GetRRFFreq(" << unit << "): unitTag=" < else if ((unitTag == RRF_UNIT_T) && (fRRFUnitTag == RRF_UNIT_Mcs)) freq = fRRFFreq/(TMath::TwoPi()*GAMMA_BAR_MUON)*1e-4; // 1e-4 need for G -> T since GAMMA_BAR_MUON is given in MHz/G -cout << endl << "debug> PMsrGlobalBlock::GetRRFFreq(" << unit << "): freq=" << freq << endl; - return freq; } diff --git a/src/classes/PMusrCanvas.cpp b/src/classes/PMusrCanvas.cpp index e34486bf..1982aa4b 100644 --- a/src/classes/PMusrCanvas.cpp +++ b/src/classes/PMusrCanvas.cpp @@ -75,7 +75,7 @@ PMusrCanvasPlotRange::PMusrCanvasPlotRange() void PMusrCanvasPlotRange::SetXRange(Double_t xmin, Double_t xmax) { if (xmin > xmax) { - cerr << endl << "PMusrCanvasPlotRange::SetXRange: **WARNING** xmin > xmax, will swap them." << endl; + cerr << endl << ">> PMusrCanvasPlotRange::SetXRange(): **WARNING** xmin > xmax, will swap them." << endl; fXmin = xmax; fXmax = xmin; } else { @@ -97,7 +97,7 @@ void PMusrCanvasPlotRange::SetXRange(Double_t xmin, Double_t xmax) void PMusrCanvasPlotRange::SetYRange(Double_t ymin, Double_t ymax) { if (ymin > ymax) { - cerr << endl << "PMusrCanvasPlotRange::SetYRange: **WARNING** ymin > ymax, will swap them." << endl; + cerr << endl << ">> PMusrCanvasPlotRange::SetYRange(): **WARNING** ymin > ymax, will swap them." << endl; fYmin = ymax; fYmax = ymin; } else { @@ -405,37 +405,72 @@ void PMusrCanvas::SetMsrHandler(PMsrHandler *msrHandler) } // check if RRF data are present - if ((fMsrHandler->GetMsrPlotList()->at(0).fRRFPacking > 0) && - (fMsrHandler->GetMsrPlotList()->at(0).fRRFFreq != 0.0)) { + if (((fMsrHandler->GetMsrPlotList()->at(0).fRRFPacking > 0) && + (fMsrHandler->GetMsrPlotList()->at(0).fRRFFreq != 0.0)) || + (fMsrHandler->GetMsrGlobal()->GetRRFPacking() > 0 && + fMsrHandler->GetMsrGlobal()->GetRRFUnit().CompareTo("??"))) { 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)"); + Int_t rrfUnitTag = -1; + Double_t rrfFreq = 0.0; + if (fMsrHandler->GetMsrPlotList()->at(0).fRRFPacking > 0) { // RRF single histo PLOT + fRRFText = new TString("RRF: "); + rrfUnitTag = fMsrHandler->GetMsrPlotList()->at(0).fRRFUnit; + rrfFreq = fMsrHandler->GetMsrPlotList()->at(0).fRRFFreq; + if (rrfUnitTag == RRF_UNIT_kHz) { + *fRRFText += TString("#nu_{RRF} = "); + *fRRFText += rrfFreq; + *fRRFText += TString(" (kHz)"); + } else if (rrfUnitTag == RRF_UNIT_MHz) { + *fRRFText += TString("#nu_{RRF} = "); + *fRRFText += rrfFreq; + *fRRFText += TString(" (MHz)"); + } else if (rrfUnitTag == RRF_UNIT_Mcs) { + *fRRFText += TString("#omega_{RRF} = "); + *fRRFText += rrfFreq; + *fRRFText += TString(" (Mc/s)"); + } else if (rrfUnitTag == RRF_UNIT_G) { + *fRRFText += TString("B_{RRF} = "); + *fRRFText += rrfFreq; + *fRRFText += TString(" (G)"); + } else if (rrfUnitTag == RRF_UNIT_T) { + *fRRFText += TString("B_{RRF} = "); + *fRRFText += rrfFreq; + *fRRFText += TString(" (T)"); + } + *fRRFText += TString(", RRF packing = "); + *fRRFText += fMsrHandler->GetMsrPlotList()->at(0).fRRFPacking; + } else { // RRF single histo FIT + fRRFText = new TString("RRF: "); + rrfUnitTag = fMsrHandler->GetMsrGlobal()->GetRRFUnitTag(); + rrfFreq = fMsrHandler->GetMsrGlobal()->GetRRFFreq(fMsrHandler->GetMsrGlobal()->GetRRFUnit().Data()); + if (rrfUnitTag == RRF_UNIT_kHz) { + *fRRFText += TString("#nu_{RRF} = "); + *fRRFText += rrfFreq; + *fRRFText += TString(" (kHz)"); + } else if (rrfUnitTag == RRF_UNIT_MHz) { + *fRRFText += TString("#nu_{RRF} = "); + *fRRFText += rrfFreq; + *fRRFText += TString(" (MHz)"); + } else if (rrfUnitTag == RRF_UNIT_Mcs) { + *fRRFText += TString("#omega_{RRF} = "); + *fRRFText += rrfFreq; + *fRRFText += TString(" (Mc/s)"); + } else if (rrfUnitTag == RRF_UNIT_G) { + *fRRFText += TString("B_{RRF} = "); + *fRRFText += rrfFreq; + *fRRFText += TString(" (G)"); + } else if (rrfUnitTag == RRF_UNIT_T) { + *fRRFText += TString("B_{RRF} = "); + *fRRFText += rrfFreq; + *fRRFText += TString(" (T)"); + } + *fRRFText += TString(", RRF packing = "); + *fRRFText += fMsrHandler->GetMsrGlobal()->GetRRFPacking(); } - *fRRFText += TString(", RRF packing = "); - *fRRFText += fMsrHandler->GetMsrPlotList()->at(0).fRRFPacking; } } @@ -609,7 +644,7 @@ void PMusrCanvas::UpdateDataTheoryPad() // first check that plot number is smaller than the maximal number of runs if ((Int_t)plotInfo.fRuns[i] > (Int_t)runs.size()) { fValid = false; - cerr << endl << "PMusrCanvas::UpdateDataTheoryPad: **ERROR** run plot number " << (Int_t)plotInfo.fRuns[i] << " is larger than the number of runs " << runs.size(); + cerr << endl << ">> PMusrCanvas::UpdateDataTheoryPad(): **ERROR** run plot number " << (Int_t)plotInfo.fRuns[i] << " is larger than the number of runs " << runs.size(); cerr << endl; return; } @@ -620,7 +655,7 @@ void PMusrCanvas::UpdateDataTheoryPad() } if (fitType == -1) { fValid = false; - cerr << endl << "PMusrCanvas::UpdateDataTheoryPad: **ERROR** plottype = " << fPlotType; + cerr << endl << ">> PMusrCanvas::UpdateDataTheoryPad(): **ERROR** plottype = " << fPlotType; cerr << ", fittype = " << runs[runNo].GetFitType() << "(RUN block)/"; cerr << "fittype = " << globalBlock->GetFitType() << "(GLOBAL block). However, they have to correspond!"; cerr << endl; @@ -643,7 +678,19 @@ void PMusrCanvas::UpdateDataTheoryPad() if (!data) { // something wrong fValid = false; // error message - cerr << endl << "PMusrCanvas::UpdateDataTheoryPad: **ERROR** couldn't obtain run no " << runNo << " for a single histogram plot"; + cerr << endl << ">> PMusrCanvas::UpdateDataTheoryPad(): **ERROR** couldn't obtain run no " << runNo << " for a single histogram plot"; + cerr << endl; + return; + } + // handle data + HandleDataSet(i, runNo, data); + break; + case MSR_FITTYPE_SINGLE_HISTO_RRF: + data = fRunList->GetSingleHistoRRF(runNo, PRunListCollection::kRunNo); + if (!data) { // something wrong + fValid = false; + // error message + cerr << endl << ">> PMusrCanvas::UpdateDataTheoryPad(): **ERROR** couldn't obtain run no " << runNo << " for a single histogram RRF plot"; cerr << endl; return; } @@ -655,7 +702,7 @@ void PMusrCanvas::UpdateDataTheoryPad() if (!data) { // something wrong fValid = false; // error message - cerr << endl << "PMusrCanvas::UpdateDataTheoryPad: **ERROR** couldn't obtain run no " << runNo << " for a asymmetry plot"; + cerr << endl << ">> PMusrCanvas::UpdateDataTheoryPad(): **ERROR** couldn't obtain run no " << runNo << " for a asymmetry plot"; cerr << endl; return; } @@ -667,7 +714,7 @@ void PMusrCanvas::UpdateDataTheoryPad() if (!data) { // something wrong fValid = false; // error message - cerr << endl << "PMusrCanvas::UpdateDataTheoryPad: **ERROR** couldn't obtain run no " << runNo << " for a mu minus single histogram plot"; + cerr << endl << ">> PMusrCanvas::UpdateDataTheoryPad(): **ERROR** couldn't obtain run no " << runNo << " for a mu minus single histogram plot"; cerr << endl; return; } @@ -679,7 +726,7 @@ void PMusrCanvas::UpdateDataTheoryPad() if (!data) { // something wrong fValid = false; // error message - cerr << endl << "PMusrCanvas::UpdateDataTheoryPad: **ERROR** couldn't obtain run no " << runNo << " for a none musr data plot"; + cerr << endl << ">> PMusrCanvas::UpdateDataTheoryPad(): **ERROR** couldn't obtain run no " << runNo << " for a none musr data plot"; cerr << endl; return; } @@ -699,7 +746,7 @@ void PMusrCanvas::UpdateDataTheoryPad() default: fValid = false; // error message - cerr << endl << "PMusrCanvas::UpdateDataTheoryPad: **ERROR** wrong plottype tag?!"; + cerr << endl << ">> PMusrCanvas::UpdateDataTheoryPad(): **ERROR** wrong plottype tag?!"; cerr << endl; return; break; @@ -1409,7 +1456,7 @@ void PMusrCanvas::SaveGraphicsAndQuit(Char_t *fileName, Char_t *graphicsFormat) } if (idx == -1) { - cerr << endl << "PMusrCanvas::SaveGraphicsAndQuit **ERROR**: fileName (" << fileName << ") is invalid." << endl; + cerr << endl << ">> PMusrCanvas::SaveGraphicsAndQuit(): **ERROR** fileName (" << fileName << ") is invalid." << endl; return; } @@ -1439,7 +1486,7 @@ void PMusrCanvas::SaveGraphicsAndQuit(Char_t *fileName, Char_t *graphicsFormat) void PMusrCanvas::ExportData(const Char_t *fileName) { if (fileName == 0) { // path file name NOT provided, generate a default path file name - cerr << endl << ">> PMusrCanvas::ExportData **ERROR** NO path file name provided. Will do nothing." << endl; + cerr << endl << ">> PMusrCanvas::ExportData(): **ERROR** NO path file name provided. Will do nothing." << endl; return; } @@ -2094,7 +2141,7 @@ void PMusrCanvas::ExportData(const Char_t *fileName) // open output data-file fout.open(fileName, iostream::out); if (!fout.is_open()) { - cerr << endl << ">> PMusrCanvas::ExportData: **ERROR** couldn't open file " << fileName << " for writing." << endl; + cerr << endl << ">> PMusrCanvas::ExportData(): **ERROR** couldn't open file " << fileName << " for writing." << endl; return; } @@ -2433,7 +2480,7 @@ void PMusrCanvas::InitMusrCanvas(const Char_t* title, Int_t wtopx, Int_t wtopy, canvasName += fPlotNumber; fMainCanvas = new TCanvas(canvasName.Data(), title, wtopx, wtopy, ww, wh); if (fMainCanvas == 0) { - cerr << endl << "PMusrCanvas::PMusrCanvas: **PANIC ERROR**: Couldn't invoke " << canvasName.Data(); + cerr << endl << ">> PMusrCanvas::PMusrCanvas(): **PANIC ERROR** Couldn't invoke " << canvasName.Data(); cerr << endl; return; } @@ -2479,7 +2526,7 @@ void PMusrCanvas::InitMusrCanvas(const Char_t* title, Int_t wtopx, Int_t wtopy, // title pad fTitlePad = new TPaveText(0.0, YTITLE, 1.0, 1.0, "NDC"); if (fTitlePad == 0) { - cerr << endl << "PMusrCanvas::PMusrCanvas: **PANIC ERROR**: Couldn't invoke fTitlePad"; + cerr << endl << ">> PMusrCanvas::PMusrCanvas(): **PANIC ERROR** Couldn't invoke fTitlePad"; cerr << endl; return; } @@ -2491,7 +2538,7 @@ void PMusrCanvas::InitMusrCanvas(const Char_t* title, Int_t wtopx, Int_t wtopy, // data/theory pad fDataTheoryPad = new TPad("dataTheoryPad", "dataTheoryPad", 0.0, YINFO, XTHEO, YTITLE); if (fDataTheoryPad == 0) { - cerr << endl << "PMusrCanvas::PMusrCanvas: **PANIC ERROR**: Couldn't invoke fDataTheoryPad"; + cerr << endl << ">> PMusrCanvas::PMusrCanvas(): **PANIC ERROR** Couldn't invoke fDataTheoryPad"; cerr << endl; return; } @@ -2501,7 +2548,7 @@ void PMusrCanvas::InitMusrCanvas(const Char_t* title, Int_t wtopx, Int_t wtopy, // parameter pad fParameterPad = new TPaveText(XTHEO, 0.5, 1.0, YTITLE, "NDC"); if (fParameterPad == 0) { - cerr << endl << "PMusrCanvas::PMusrCanvas: **PANIC ERROR**: Couldn't invoke fParameterPad"; + cerr << endl << ">> PMusrCanvas::PMusrCanvas(): **PANIC ERROR** Couldn't invoke fParameterPad"; cerr << endl; return; } @@ -2512,7 +2559,7 @@ void PMusrCanvas::InitMusrCanvas(const Char_t* title, Int_t wtopx, Int_t wtopy, // theory pad fTheoryPad = new TPaveText(XTHEO, 0.1, 1.0, 0.5, "NDC"); if (fTheoryPad == 0) { - cerr << endl << "PMusrCanvas::PMusrCanvas: **PANIC ERROR**: Couldn't invoke fTheoryPad"; + cerr << endl << ">> PMusrCanvas::PMusrCanvas(): **PANIC ERROR** Couldn't invoke fTheoryPad"; cerr << endl; return; } @@ -2524,7 +2571,7 @@ void PMusrCanvas::InitMusrCanvas(const Char_t* title, Int_t wtopx, Int_t wtopy, // info pad fInfoPad = new TLegend(0.0, 0.0, 1.0, YINFO, "NDC"); if (fInfoPad == 0) { - cerr << endl << "PMusrCanvas::PMusrCanvas: **PANIC ERROR**: Couldn't invoke fInfoPad"; + cerr << endl << ">> PMusrCanvas::PMusrCanvas(): **PANIC ERROR** Couldn't invoke fInfoPad"; cerr << endl; return; } @@ -2929,9 +2976,9 @@ void PMusrCanvas::HandleDataSet(UInt_t plotNo, UInt_t runNo, PRunData *data) Double_t dval = (startFitRange - data->GetDataTimeStart())/data->GetDataTimeStep(); if (dval < 0.0) { // make sure that startBin >= 0 startBin = 0; - cerr << endl << "PMusrCanvas::HandleDataSet() **WARNING** found startBin data < 0 for 'use_fit_range', will set it to 0" << endl << endl; + cerr << endl << ">> PMusrCanvas::HandleDataSet(): **WARNING** found startBin data < 0 for 'use_fit_range', will set it to 0" << endl << endl; } else if (dval >= (Double_t)data->GetValue()->size()) { // make sure that startBin <= length of data vector - cerr << endl << ">> PMusrCanvas::HandleDataSet() **WARNING** found startBin data=" << (UInt_t)dval << " >= data vector size=" << data->GetValue()->size() << " for 'use_fit_range',"; + cerr << endl << ">> PMusrCanvas::HandleDataSet(): **WARNING** found startBin data=" << (UInt_t)dval << " >= data vector size=" << data->GetValue()->size() << " for 'use_fit_range',"; cerr << endl << ">> will set it to data vector size" << endl << endl; startBin = data->GetValue()->size(); } else { @@ -2945,9 +2992,9 @@ void PMusrCanvas::HandleDataSet(UInt_t plotNo, UInt_t runNo, PRunData *data) dval = (endFitRange - data->GetDataTimeStart())/data->GetDataTimeStep(); if (dval < 0.0) { // make sure that endBin >= 0 endBin = 0; - cerr << endl << "PMusrCanvas::HandleDataSet() **WARNING** found endBin data < 0 for 'use_fit_range', will set it to 0" << endl << endl; + cerr << endl << ">> PMusrCanvas::HandleDataSet(): **WARNING** found endBin data < 0 for 'use_fit_range', will set it to 0" << endl << endl; } else if (dval >= (Double_t)data->GetValue()->size()) { // make sure that endBin <= length of data vector - cerr << endl << ">> PMusrCanvas::HandleDataSet() **WARNING** found endBin data=" << (UInt_t)dval << " >= data vector size=" << data->GetValue()->size() << " for 'use_fit_range',"; + cerr << endl << ">> PMusrCanvas::HandleDataSet(): **WARNING** found endBin data=" << (UInt_t)dval << " >= data vector size=" << data->GetValue()->size() << " for 'use_fit_range',"; cerr << endl << ">> will set it to data vector size" << endl << endl; endBin = data->GetValue()->size(); } else { @@ -2960,9 +3007,9 @@ void PMusrCanvas::HandleDataSet(UInt_t plotNo, UInt_t runNo, PRunData *data) Double_t dval = (fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fTmin[runNo] - data->GetDataTimeStart())/data->GetDataTimeStep(); if (dval < 0.0) { // make sure that startBin >= 0 startBin = 0; - cerr << endl << "PMusrCanvas::HandleDataSet() **WARNING** found startBin data < 0 for 'sub_ranges', will set it to 0" << endl << endl; + cerr << endl << ">> PMusrCanvas::HandleDataSet(): **WARNING** found startBin data < 0 for 'sub_ranges', will set it to 0" << endl << endl; } else if (dval >= (Double_t)data->GetValue()->size()) { // make sure that startBin <= length of data vector - cerr << endl << ">> PMusrCanvas::HandleDataSet() **WARNING** found startBin data=" << (UInt_t)dval << " >= data vector size=" << data->GetValue()->size() << " for 'sub_ranges',"; + cerr << endl << ">> PMusrCanvas::HandleDataSet(): **WARNING** found startBin data=" << (UInt_t)dval << " >= data vector size=" << data->GetValue()->size() << " for 'sub_ranges',"; cerr << endl << ">> will set it to data vector size" << endl << endl; startBin = data->GetValue()->size(); } else { @@ -2972,9 +3019,9 @@ void PMusrCanvas::HandleDataSet(UInt_t plotNo, UInt_t runNo, PRunData *data) dval = (fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fTmax[runNo] - data->GetDataTimeStart())/data->GetDataTimeStep(); if (dval < 0.0) { // make sure that endBin >= 0 endBin = 0; - cerr << endl << "PMusrCanvas::HandleDataSet() **WARNING** found endBin data < 0 for 'sub_ranges', will set it to 0" << endl << endl; + cerr << endl << ">> PMusrCanvas::HandleDataSet(): **WARNING** found endBin data < 0 for 'sub_ranges', will set it to 0" << endl << endl; } else if (dval >= (Double_t)data->GetValue()->size()) { // make sure that endtBin <= length of data vector - cerr << endl << ">> PMusrCanvas::HandleDataSet() **WARNING** found endBin data=" << (UInt_t)dval << " >= data vector size=" << data->GetValue()->size() << " for 'sub_ranges',"; + cerr << endl << ">> PMusrCanvas::HandleDataSet(): **WARNING** found endBin data=" << (UInt_t)dval << " >= data vector size=" << data->GetValue()->size() << " for 'sub_ranges',"; cerr << endl << ">> will set it to data vector size" << endl << endl; endBin = data->GetValue()->size(); } else { @@ -3061,9 +3108,9 @@ void PMusrCanvas::HandleDataSet(UInt_t plotNo, UInt_t runNo, PRunData *data) Double_t dval = (startFitRange - data->GetDataTimeStart())/data->GetTheoryTimeStep(); if (dval < 0.0) { // make sure that startBin >= 0 startBin = 0; - cerr << endl << "PMusrCanvas::HandleDataSet() **WARNING** found startBin theory < 0 for 'use_fit_range', will set it to 0" << endl << endl; + cerr << endl << ">> PMusrCanvas::HandleDataSet(): **WARNING** found startBin theory < 0 for 'use_fit_range', will set it to 0" << endl << endl; } else if (dval >= (Double_t)data->GetTheory()->size()) { // make sure that startBin <= length of theory vector - cerr << endl << ">> PMusrCanvas::HandleDataSet() **WARNING** found startBin theory=" << (UInt_t)dval << " >= theory vector size=" << data->GetTheory()->size() << " for 'use_fit_range',"; + cerr << endl << ">> PMusrCanvas::HandleDataSet(): **WARNING** found startBin theory=" << (UInt_t)dval << " >= theory vector size=" << data->GetTheory()->size() << " for 'use_fit_range',"; cerr << endl << ">> will set it to theory vector size" << endl << endl; startBin = data->GetTheory()->size(); } else { @@ -3077,9 +3124,9 @@ void PMusrCanvas::HandleDataSet(UInt_t plotNo, UInt_t runNo, PRunData *data) dval = (endFitRange - data->GetDataTimeStart())/data->GetTheoryTimeStep(); if (dval < 0.0) { // make sure that endBin >= 0 endBin = 0; - cerr << endl << "PMusrCanvas::HandleDataSet() **WARNING** found endBin theory < 0 for 'use_fit_range', will set it to 0" << endl << endl; + cerr << endl << ">> PMusrCanvas::HandleDataSet(): **WARNING** found endBin theory < 0 for 'use_fit_range', will set it to 0" << endl << endl; } else if (dval >= (Double_t)data->GetTheory()->size()) { // make sure that endBin <= length of theory vector - cerr << endl << ">> PMusrCanvas::HandleDataSet() **WARNING** found endBin theory=" << (UInt_t)dval << " >= theory vector size=" << data->GetTheory()->size() << " for 'use_fit_range',"; + cerr << endl << ">> PMusrCanvas::HandleDataSet(): **WARNING** found endBin theory=" << (UInt_t)dval << " >= theory vector size=" << data->GetTheory()->size() << " for 'use_fit_range',"; cerr << endl << ">> will set it to theory vector size" << endl << endl; endBin = data->GetTheory()->size(); } else { @@ -3095,9 +3142,9 @@ void PMusrCanvas::HandleDataSet(UInt_t plotNo, UInt_t runNo, PRunData *data) Double_t dval = (fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fTmin[runNo] -data->GetDataTimeStart())/data->GetTheoryTimeStep(); if (dval < 0.0) { // make sure that startBin >= 0 startBin = 0; - cerr << endl << "PMusrCanvas::HandleDataSet() **WARNING** found startBin theory < 0 for 'sub_ranges', will set it to 0" << endl << endl; + cerr << endl << ">> PMusrCanvas::HandleDataSet(): **WARNING** found startBin theory < 0 for 'sub_ranges', will set it to 0" << endl << endl; } else if (dval >= (Double_t)data->GetTheory()->size()) { // make sure that startBin <= length of theory vector - cerr << endl << ">> PMusrCanvas::HandleDataSet() **WARNING** found startBin theory=" << (UInt_t)dval << " >= theory vector size=" << data->GetTheory()->size() << " for 'sub_ranges',"; + cerr << endl << ">> PMusrCanvas::HandleDataSet(): **WARNING** found startBin theory=" << (UInt_t)dval << " >= theory vector size=" << data->GetTheory()->size() << " for 'sub_ranges',"; cerr << endl << ">> will set it to theory vector size" << endl << endl; startBin = data->GetTheory()->size(); } else { @@ -3107,9 +3154,9 @@ void PMusrCanvas::HandleDataSet(UInt_t plotNo, UInt_t runNo, PRunData *data) dval = (fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fTmax[runNo] -data->GetDataTimeStart())/data->GetTheoryTimeStep(); if (dval < 0.0) { // make sure that endBin >= 0 endBin = 0; - cerr << endl << "PMusrCanvas::HandleDataSet() **WARNING** found endBin theory < 0 for 'sub_ranges', will set it to 0" << endl << endl; + cerr << endl << ">> PMusrCanvas::HandleDataSet(): **WARNING** found endBin theory < 0 for 'sub_ranges', will set it to 0" << endl << endl; } else if (dval >= (Double_t)data->GetTheory()->size()) { // make sure that endtBin <= length of theory vector - cerr << endl << ">> PMusrCanvas::HandleDataSet() **WARNING** found endBin theory=" << (UInt_t)dval << " >= theory vector size=" << data->GetTheory()->size() << " for 'sub_ranges',"; + cerr << endl << ">> PMusrCanvas::HandleDataSet(): **WARNING** found endBin theory=" << (UInt_t)dval << " >= theory vector size=" << data->GetTheory()->size() << " for 'sub_ranges',"; cerr << endl << ">> will set it to theory vector size" << endl << endl; endBin = data->GetTheory()->size(); } else { @@ -3418,7 +3465,7 @@ void PMusrCanvas::HandleFourier() // calculate fourier transform of the data PFourier fourierData(fData[i].data, fFourier.fUnits, startTime, endTime, fFourier.fDCCorrected, fFourier.fFourierPower); if (!fourierData.IsValid()) { - cerr << endl << "**SEVERE ERROR** PMusrCanvas::HandleFourier: couldn't invoke PFourier to calculate the Fourier data ..." << endl; + cerr << endl << ">> PMusrCanvas::HandleFourier(): **SEVERE ERROR** couldn't invoke PFourier to calculate the Fourier data ..." << endl; return; } fourierData.Transform(fFourier.fApodization); @@ -3458,7 +3505,7 @@ void PMusrCanvas::HandleFourier() Int_t powerPad = (Int_t)round(log((endTime-startTime)/fData[i].theory->GetBinWidth(1))/log(2))+3; PFourier fourierTheory(fData[i].theory, fFourier.fUnits, startTime, endTime, fFourier.fDCCorrected, powerPad); if (!fourierTheory.IsValid()) { - cerr << endl << "**SEVERE ERROR** PMusrCanvas::HandleFourier: couldn't invoke PFourier to calculate the Fourier theory ..." << endl; + cerr << endl << ">> PMusrCanvas::HandleFourier(): **SEVERE ERROR** couldn't invoke PFourier to calculate the Fourier theory ..." << endl; return; } fourierTheory.Transform(fFourier.fApodization); @@ -3578,7 +3625,7 @@ void PMusrCanvas::HandleDifferenceFourier() // calculate fourier transform of the data PFourier fourierData(fData[i].diff, fFourier.fUnits, startTime, endTime, fFourier.fDCCorrected, fFourier.fFourierPower); if (!fourierData.IsValid()) { - cerr << endl << "**SEVERE ERROR** PMusrCanvas::HandleFourier: couldn't invoke PFourier to calculate the Fourier diff ..." << endl; + cerr << endl << ">> PMusrCanvas::HandleFourier(): **SEVERE ERROR** couldn't invoke PFourier to calculate the Fourier diff ..." << endl; return; } fourierData.Transform(fFourier.fApodization); @@ -4648,14 +4695,14 @@ void PMusrCanvas::PlotData(Bool_t unzoom) fDataTheoryPad->SetLogy(1); // set x-axis label - fHistoFrame->GetXaxis()->SetTitle("time (#mus)"); + fHistoFrame->GetXaxis()->SetTitle("Time (#mus)"); // set y-axis label TString yAxisTitle; PMsrRunList *runList = fMsrHandler->GetMsrRunList(); switch (fPlotType) { case MSR_PLOT_SINGLE_HISTO: if (runList->at(0).IsLifetimeCorrected()) { // lifetime correction - yAxisTitle = "asymmetry"; + yAxisTitle = "Asymmetry"; } else { // no liftime correction if (fScaleN0AndBkg) yAxisTitle = "N(t) per nsec"; @@ -4663,8 +4710,11 @@ void PMusrCanvas::PlotData(Bool_t unzoom) yAxisTitle = "N(t) per bin"; } break; + case MSR_PLOT_SINGLE_HISTO_RRF: + yAxisTitle = "RRF Asymmetry"; + break; case MSR_PLOT_ASYM: - yAxisTitle = "asymmetry"; + yAxisTitle = "Asymmetry"; break; case MSR_PLOT_MU_MINUS: yAxisTitle = "N(t) per bin"; diff --git a/src/classes/PRunListCollection.cpp b/src/classes/PRunListCollection.cpp index bc58df51..d5b1d94e 100644 --- a/src/classes/PRunListCollection.cpp +++ b/src/classes/PRunListCollection.cpp @@ -106,8 +106,6 @@ Bool_t PRunListCollection::Add(Int_t runNo, EPMusrHandleTag tag) fitType = (*fMsrInfo->GetMsrGlobal()).GetFitType(); } -cout << "debug> PRunListCollection::Add(): (runNo: " << runNo << "), fitType = " << fitType << endl; - switch (fitType) { case PRUN_SINGLE_HISTO: fRunSingleHistoList.push_back(new PRunSingleHisto(fMsrInfo, fData, runNo, tag)); @@ -115,7 +113,6 @@ cout << "debug> PRunListCollection::Add(): (runNo: " << runNo << "), fitType = " success = false; break; case PRUN_SINGLE_HISTO_RRF: -cout << "debug> PRunListCollection::Add(): add RRF single histo run to PRunListCollection (runNo: " << runNo << ")" << endl; fRunSingleHistoRRFList.push_back(new PRunSingleHistoRRF(fMsrInfo, fData, runNo, tag)); if (!fRunSingleHistoRRFList[fRunSingleHistoRRFList.size()-1]->IsValid()) success = false; @@ -378,16 +375,17 @@ Double_t PRunListCollection::GetSingleRunChisq(const std::vector& par, return chisq; } + Int_t subIdx = 0; Int_t type = fMsrInfo->GetMsrRunList()->at(idx).GetFitType(); - if (type == -1) { // i.e. not forun in the RUN block, try the GLOBAL block + if (type == -1) { // i.e. not found in the RUN block, try the GLOBAL block type = fMsrInfo->GetMsrGlobal()->GetFitType(); - } - - // count how many entries of this fit-type are present up to idx - UInt_t subIdx = 0; - for (UInt_t i=0; iGetMsrRunList()->at(i).GetFitType() == type) - subIdx++; + subIdx = idx; + } else { // found in the RUN block + // count how many entries of this fit-type are present up to idx + for (UInt_t i=0; iGetMsrRunList()->at(i).GetFitType() == type) + subIdx++; + } } // return the chisq of the single run diff --git a/src/classes/PRunSingleHistoRRF.cpp b/src/classes/PRunSingleHistoRRF.cpp index 60f69ca7..8f1dfb50 100644 --- a/src/classes/PRunSingleHistoRRF.cpp +++ b/src/classes/PRunSingleHistoRRF.cpp @@ -57,7 +57,7 @@ PRunSingleHistoRRF::PRunSingleHistoRRF() : PRunBase() { fNoOfFitBins = 0; fBackground = 0; - fPacking = -1; + fRRFPacking = -1; // the 2 following variables are need in case fit range is given in bins, and since // the fit range can be changed in the command block, these variables need to be accessible @@ -80,13 +80,26 @@ PRunSingleHistoRRF::PRunSingleHistoRRF(PMsrHandler *msrInfo, PRunDataHandler *ra { fNoOfFitBins = 0; - fPacking = fRunInfo->GetPacking(); - if (fPacking == -1) { // i.e. packing is NOT given in the RUN-block, it must be given in the GLOBAL-block - fPacking = fMsrInfo->GetMsrGlobal()->GetPacking(); + PMsrGlobalBlock *global = msrInfo->GetMsrGlobal(); + + if (!global->IsPresent()) { + cerr << endl << ">> PRunSingleHistoRRF::PRunSingleHistoRRF(): **SEVERE ERROR**: no GLOBAL-block present!"; + cerr << endl << ">> For Single Histo RRF the GLOBAL-block is mandatory! Please fix this first."; + cerr << endl; + fValid = false; + return; } - if (fPacking == -1) { // this should NOT happen, somethin is severely wrong - cerr << endl << ">> PRunSingleHistoRRF::PRunSingleHistoRRF(): **SEVERE ERROR**: Couldn't find any packing information!"; - cerr << endl << ">> This is very bad :-(, will quit ..."; + + if (!global->GetRRFUnit().CompareTo("??")) { + cerr << endl << ">> PRunSingleHistoRRF::PRunSingleHistoRRF(): **SEVERE ERROR**: no RRF-Frequency found!"; + cerr << endl; + fValid = false; + return; + } + + fRRFPacking = global->GetRRFPacking(); + if (fRRFPacking == -1) { + cerr << endl << ">> PRunSingleHistoRRF::PRunSingleHistoRRF(): **SEVERE ERROR**: no RRF-Packing found!"; cerr << endl; fValid = false; return; @@ -521,17 +534,10 @@ Bool_t PRunSingleHistoRRF::PrepareData() // get the fit range for the current RUN block GetProperFitRange(globalBlock); - // get the lifetimecorrection flag - Bool_t lifetimecorrection = false; - PMsrPlotList *plot = fMsrInfo->GetMsrPlotList(); - lifetimecorrection = plot->at(0).fLifeTimeCorrection; - // do the more fit/view specific stuff if (fHandleTag == kFit) success = PrepareFitData(runData, histoNo[0]); - else if ((fHandleTag == kView) && !lifetimecorrection) - success = PrepareRawViewData(runData, histoNo[0]); - else if ((fHandleTag == kView) && lifetimecorrection) + else if (fHandleTag == kView) success = PrepareViewData(runData, histoNo[0]); else success = false; @@ -546,11 +552,13 @@ Bool_t PRunSingleHistoRRF::PrepareData() // PrepareFitData (protected) //-------------------------------------------------------------------------- /** - *

Take the pre-processed data (i.e. grouping and addrun are preformed) and form the histogram for fitting. + *

Take the pre-processed data (i.e. grouping and addrun are preformed) and form the RRF histogram for fitting. * The following steps are preformed: * -# get fit start/stop time * -# check that 'first good data bin', 'last good data bin', and 't0' make any sense * -# check how the background shall be handled, i.e. fitted, subtracted from background estimate data range, or subtacted from a given fixed background. + * -# estimate N0 + * -# RRF transformation * -# packing (i.e rebinning) * * return: @@ -571,28 +579,26 @@ Bool_t PRunSingleHistoRRF::PrepareFitData(PRawRunData* runData, const UInt_t his // initially fForward is the "raw data set" (i.e. grouped histo and raw runs already added) to be fitted. This means fForward = N(t) at this point. // 1) check how the background shall be handled - if (fRunInfo->GetBkgFitParamNo() == -1) { // bkg shall **NOT** be fitted - // subtract background from histogramms ------------------------------------------ - if (fRunInfo->GetBkgFix(0) == PMUSR_UNDEFINED) { // no fixed background given - if (fRunInfo->GetBkgRange(0) >= 0) { - if (!EstimateBkg(histoNo)) - return false; - } else { // no background given to do the job, try estimate - fRunInfo->SetBkgRange(static_cast(fT0s[0]*0.1), 0); - fRunInfo->SetBkgRange(static_cast(fT0s[0]*0.6), 1); - cerr << endl << ">> PRunSingleHistoRRF::PrepareFitData(): **WARNING** Neither fix background nor background bins are given!"; - cerr << endl << ">> Will try the following: bkg start = " << fRunInfo->GetBkgRange(0) << ", bkg end = " << fRunInfo->GetBkgRange(1); - cerr << endl << ">> NO WARRANTY THAT THIS MAKES ANY SENSE! Better check ..."; - cerr << endl; - if (!EstimateBkg(histoNo)) - return false; - } - } else { // fixed background given - for (UInt_t i=0; iGetBkgFix(0); - } - fBackground = fRunInfo->GetBkgFix(0); + // subtract background from histogramms ------------------------------------------ + if (fRunInfo->GetBkgFix(0) == PMUSR_UNDEFINED) { // no fixed background given + if (fRunInfo->GetBkgRange(0) >= 0) { + if (!EstimateBkg(histoNo)) + return false; + } else { // no background given to do the job, try estimate + fRunInfo->SetBkgRange(static_cast(fT0s[0]*0.1), 0); + fRunInfo->SetBkgRange(static_cast(fT0s[0]*0.6), 1); + cerr << endl << ">> PRunSingleHistoRRF::PrepareFitData(): **WARNING** Neither fix background nor background bins are given!"; + cerr << endl << ">> Will try the following: bkg start = " << fRunInfo->GetBkgRange(0) << ", bkg end = " << fRunInfo->GetBkgRange(1); + cerr << endl << ">> NO WARRANTY THAT THIS MAKES ANY SENSE! Better check ..."; + cerr << endl; + if (!EstimateBkg(histoNo)) + return false; } + } else { // fixed background given + for (UInt_t i=0; iGetBkgFix(0); + } + fBackground = fRunInfo->GetBkgFix(0); } // here fForward = N(t) - Nbkg @@ -600,9 +606,6 @@ Bool_t PRunSingleHistoRRF::PrepareFitData(PRawRunData* runData, const UInt_t his // 2) N(t) - Nbkg -> exp(+t/tau) [N(t)-Nbkg] Double_t startTime = fTimeResolution * ((Double_t)fGoodBins[0] - (Double_t)t0); - cout << endl << "debug> PRunSingleHistoRRF::PrepareFitData(): t0=" << t0 << ", fGoodBins[0]=" << fGoodBins[0] << ", fGoodBins[1]=" << fGoodBins[1] << endl; - cout << endl << "debug> PRunSingleHistoRRF::PrepareFitData(): startTime=" << startTime << endl; - cout << endl << "debug> PRunSingleHistoRRF::PrepareFitData(): endTime =" << startTime+fTimeResolution*((Double_t)fGoodBins[1]-(Double_t)t0) << endl; Double_t time_tau=0.0; Double_t exp_t_tau=0.0; @@ -642,8 +645,6 @@ Bool_t PRunSingleHistoRRF::PrepareFitData(PRawRunData* runData, const UInt_t his PMsrGlobalBlock *globalBlock = fMsrInfo->GetMsrGlobal(); Double_t wRRF = globalBlock->GetRRFFreq("Mc"); Double_t phaseRRF = globalBlock->GetRRFPhase()*TMath::TwoPi()/180.0; - cout << "debug> PRunSingleHistoRRF::PrepareFitData(): wRRF =" << wRRF << endl; - cout << "debug> PRunSingleHistoRRF::PrepareFitData(): phaseRRF =" << phaseRRF << endl; Double_t time = 0.0; for (Int_t i=fGoodBins[0]; iGetRRFPacking(); Double_t dval=0.0; for (Int_t i=fGoodBins[0]; i 1 - if (((i-fGoodBins[0]) % packingRRF == 0) && (i != fGoodBins[0])) { // fill data - dval /= packingRRF; + if (((i-fGoodBins[0]) % fRRFPacking == 0) && (i != fGoodBins[0])) { // fill data + dval /= fRRFPacking; fData.AppendValue(dval); // reset dval dval = 0.0; @@ -672,200 +672,22 @@ Bool_t PRunSingleHistoRRF::PrepareFitData(PRawRunData* runData, const UInt_t his dval = 0.0; // the packed RRF asymmetry error for (Int_t i=fGoodBins[0]; iTake the pre-processed data (i.e. grouping and addrun are preformed) and form the histogram for viewing - * without any life time correction. - *

The following steps are preformed: - * -# check if view packing is whished. - * -# check that 'first good data bin', 'last good data bin', and 't0' makes any sense - * -# packing (i.e. rebinnig) - * -# calculate theory - * - * return: - * - true, if everything went smooth - * - false, otherwise. - * - * \param runData raw run data handler - * \param histoNo forward histogram number - */ -Bool_t PRunSingleHistoRRF::PrepareRawViewData(PRawRunData* runData, const UInt_t histoNo) -{ -/* - // check if view_packing is wished - Int_t packing = fPacking; - if (fMsrInfo->GetMsrPlotList()->at(0).fViewPacking > 0) { - packing = fMsrInfo->GetMsrPlotList()->at(0).fViewPacking; - } - - // calculate necessary norms - Double_t dataNorm = 1.0, theoryNorm = 1.0; - if (fScaleN0AndBkg) { - dataNorm = 1.0/ (packing * (fTimeResolution * 1.0e3)); // fTimeResolution us->ns - } else if (!fScaleN0AndBkg && (fMsrInfo->GetMsrPlotList()->at(0).fViewPacking > 0)) { - theoryNorm = (Double_t)fMsrInfo->GetMsrPlotList()->at(0).fViewPacking/(Double_t)fPacking; - } - - // raw data, since PMusrCanvas is doing ranging etc. - // start = the first bin which is a multiple of packing backward from first good data bin - Int_t start = fGoodBins[0] - (fGoodBins[0]/packing)*packing; - // end = last bin starting from start which is a multipl of packing and still within the data - Int_t end = start + ((fForward.size()-start)/packing)*packing; - // check if data range has been provided, and if not try to estimate them - if (start < 0) { - Int_t offset = (Int_t)(10.0e-3/fTimeResolution); - start = ((Int_t)fT0s[0]+offset) - (((Int_t)fT0s[0]+offset)/packing)*packing; - end = start + ((fForward.size()-start)/packing)*packing; - cerr << endl << ">> PRunSingleHistoRRF::PrepareRawViewData(): **WARNING** data range was not provided, will try data range start = " << start << "."; - cerr << endl << ">> NO WARRANTY THAT THIS DOES MAKE ANY SENSE."; - cerr << endl; - } - // check if start, end, and t0 make any sense - // 1st check if start and end are in proper order - if (end < start) { // need to swap them - Int_t keep = end; - end = start; - start = keep; - } - // 2nd check if start is within proper bounds - if ((start < 0) || (start > (Int_t)fForward.size())) { - cerr << endl << ">> PRunSingleHistoRRF::PrepareRawViewData(): **ERROR** start data bin doesn't make any sense!"; - cerr << endl; - return false; - } - // 3rd check if end is within proper bounds - if ((end < 0) || (end > (Int_t)fForward.size())) { - cerr << endl << ">> PRunSingleHistoRRF::PrepareRawViewData(): **ERROR** end data bin doesn't make any sense!"; - cerr << endl; - return false; - } - - // everything looks fine, hence fill data set - Int_t t0 = (Int_t)fT0s[0]; - Double_t value = 0.0; - // data start at data_start-t0 - // time shifted so that packing is included correctly, i.e. t0 == t0 after packing - fData.SetDataTimeStart(fTimeResolution*((Double_t)start-(Double_t)t0+(Double_t)(packing-1)/2.0)); - fData.SetDataTimeStep(fTimeResolution*packing); - - for (Int_t i=start; i par; - PMsrParamList *paramList = fMsrInfo->GetMsrParamList(); - for (UInt_t i=0; isize(); i++) - par.push_back((*paramList)[i].fValue); - - // calculate asymmetry - Double_t N0; - // check if norm is a parameter or a function - if (fRunInfo->GetNormParamNo() < MSR_PARAM_FUN_OFFSET) { // norm is a parameter - N0 = par[fRunInfo->GetNormParamNo()-1]; - } else { // norm is a function - // get function number - UInt_t funNo = fRunInfo->GetNormParamNo()-MSR_PARAM_FUN_OFFSET; - // evaluate function - N0 = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par); - } - N0 *= theoryNorm; - - // get tau - Double_t tau; - if (fRunInfo->GetLifetimeParamNo() != -1) - tau = par[fRunInfo->GetLifetimeParamNo()-1]; - else - tau = PMUON_LIFETIME; - - // get background - Double_t bkg; - if (fRunInfo->GetBkgFitParamNo() == -1) { // bkg not fitted - if (fRunInfo->GetBkgFix(0) == PMUSR_UNDEFINED) { // no fixed background given (background interval) - if (fRunInfo->GetBkgRange(0) >= 0) { // background range given - if (!EstimateBkg(histoNo)) - return false; - } else { // no background given to do the job, try estimate - fRunInfo->SetBkgRange(static_cast(fT0s[0]*0.1), 0); - fRunInfo->SetBkgRange(static_cast(fT0s[0]*0.6), 1); - cerr << endl << ">> PRunSingleHistoRRF::PrepareRawViewData(): **WARNING** Neither fix background nor background bins are given!"; - cerr << endl << ">> Will try the following: bkg start = " << fRunInfo->GetBkgRange(0) << ", bkg end = " << fRunInfo->GetBkgRange(1); - cerr << endl << ">> NO WARRANTY THAT THIS MAKES ANY SENSE! Better check ..."; - cerr << endl; - if (!EstimateBkg(histoNo)) - return false; - } - bkg = fBackground; - } else { // fixed bkg given - bkg = fRunInfo->GetBkgFix(0); - } - } else { // bkg fitted - bkg = par[fRunInfo->GetBkgFitParamNo()-1]; - } - bkg *= theoryNorm; - - // calculate functions - for (Int_t i=0; iGetNoOfFuncs(); i++) { - fFuncValues[i] = fMsrInfo->EvalFunc(fMsrInfo->GetFuncNo(i), *fRunInfo->GetMap(), par); - } - - // calculate theory - UInt_t size = fForward.size(); - Double_t factor = 1.0; - if (fData.GetValue()->size() * 10 > fForward.size()) { - size = fData.GetValue()->size() * 10; - factor = (Double_t)fForward.size() / (Double_t)size; - } - Double_t time; - Double_t theoryValue; - fData.SetTheoryTimeStart(fData.GetDataTimeStart()); - fData.SetTheoryTimeStep(fTimeResolution*factor); - for (UInt_t i=0; iFunc(time, par, fFuncValues); - if (fabs(theoryValue) > 1.0e10) { // dirty hack needs to be fixed!! - theoryValue = 0.0; - } - fData.AppendTheoryValue(N0*TMath::Exp(-time/tau)*(1.0+theoryValue)+bkg); - } - - // clean up - par.clear(); -*/ - return true; -} - //-------------------------------------------------------------------------- // PrepareViewData (protected) //-------------------------------------------------------------------------- @@ -878,16 +700,6 @@ Bool_t PRunSingleHistoRRF::PrepareRawViewData(PRawRunData* runData, const UInt_t * -# transform data sets (see below). * -# calculate theory * - *

Muon life time corrected data: Starting from - * \f[ N(t) = N_0 e^{-t/\tau} [ 1 + A(t) ] + \mathrm{Bkg} \f] - * it follows that - * \f[ A(t) = (-1) + e^{+t/\tau}\, \frac{N(t)-\mathrm{Bkg}}{N_0}. \f] - * For the error estimate only the statistical error of \f$ N(t) \f$ is used, and hence - * \f[ \Delta A(t) = \frac{e^{t/\tau}}{N_0}\,\sqrt{\frac{N(t)}{p}} \f] - * where \f$ p \f$ is the packing, and \f$ N(t) \f$ are the packed data, i.e. - * \f[ N(t_i) = \frac{1}{p}\, \sum_{j=i}^{i+p} n_j \f] - * with \f$ n_j \f$ the raw histogram data bins. - * * return: * - true, if everything went smooth * - false, otherwise @@ -897,66 +709,27 @@ Bool_t PRunSingleHistoRRF::PrepareRawViewData(PRawRunData* runData, const UInt_t */ Bool_t PRunSingleHistoRRF::PrepareViewData(PRawRunData* runData, const UInt_t histoNo) { -/* - // check if view_packing is wished. This is a global option for all PLOT blocks! - Int_t packing = fPacking; - 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; + // -------------- + // prepare data + // -------------- + + // prepare RRF single histo + PrepareFitData(runData, histoNo); + + // check for view packing + Int_t viewPacking = fMsrInfo->GetMsrPlotList()->at(0).fViewPacking; + if (viewPacking > 0) { + if (viewPacking < fRRFPacking) { + cerr << ">> PRunSingleHistoRRF::PrepareViewData(): **WARNING** Found View Packing (" << viewPacking << ") < RRF Packing (" << fRRFPacking << ")."; + cerr << ">> Will ignore View Packing." << endl; + } else { + // STILL MISSING + } } - // calculate necessary norms - Double_t dataNorm = 1.0, theoryNorm = 1.0; - if (fScaleN0AndBkg) { - dataNorm = 1.0/ (packing * (fTimeResolution * 1.0e3)); // fTimeResolution us->ns - } else if (!fScaleN0AndBkg && (fMsrInfo->GetMsrPlotList()->at(0).fViewPacking > 0)) { - theoryNorm = (Double_t)fMsrInfo->GetMsrPlotList()->at(0).fViewPacking/(Double_t)fPacking; - } - - // 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 - // first get start data, end data, and t0 - Int_t t0 = (Int_t)fT0s[0]; - - // start = the first bin which is a multiple of packing backward from first good data bin - Int_t start = fGoodBins[0] - (fGoodBins[0]/packing)*packing; - // end = last bin starting from start which is a multiple of packing and still within the data - Int_t end = start + ((fForward.size()-start)/packing)*packing; - - // check if data range has been provided, and if not try to estimate them - if (start < 0) { - Int_t offset = (Int_t)(10.0e-3/fTimeResolution); - start = ((Int_t)fT0s[0]+offset) - (((Int_t)fT0s[0]+offset)/packing)*packing; - end = start + ((fForward.size()-start)/packing)*packing; - cerr << endl << ">> PRunSingleHistoRRF::PrepareViewData(): **WARNING** data range was not provided, will try data range start = " << start << "."; - cerr << endl << ">> NO WARRANTY THAT THIS DOES MAKE ANY SENSE."; - cerr << endl; - } - - // check if start, end, and t0 make any sense - // 1st check if start and end are in proper order - if (end < start) { // need to swap them - Int_t keep = end; - end = start; - start = keep; - } - // 2nd check if start is within proper bounds - if ((start < 0) || (start > (Int_t)fForward.size())) { - cerr << endl << ">> PRunSingleHistoRRF::PrepareViewData(): **ERROR** start data bin doesn't make any sense!"; - cerr << endl; - return false; - } - // 3rd check if end is within proper bounds - if ((end < 0) || (end > (Int_t)fForward.size())) { - cerr << endl << ">> PRunSingleHistoRRF::PrepareViewData(): **ERROR** end data bin doesn't make any sense!"; - cerr << endl; - return false; - } - - // everything looks fine, hence fill data set + // -------------- + // prepare theory + // -------------- // feed the parameter vector std::vector par; @@ -964,182 +737,34 @@ Bool_t PRunSingleHistoRRF::PrepareViewData(PRawRunData* runData, const UInt_t hi for (UInt_t i=0; isize(); i++) par.push_back((*paramList)[i].fValue); - // calculate asymmetry - Double_t N0; - // check if norm is a parameter or a function - if (fRunInfo->GetNormParamNo() < MSR_PARAM_FUN_OFFSET) { // norm is a parameter - N0 = par[fRunInfo->GetNormParamNo()-1]; - } else { // norm is a function - // get function number - UInt_t funNo = fRunInfo->GetNormParamNo()-MSR_PARAM_FUN_OFFSET; - // evaluate function - N0 = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par); - } - N0 *= theoryNorm; - - // get tau - Double_t tau; - if (fRunInfo->GetLifetimeParamNo() != -1) - tau = par[fRunInfo->GetLifetimeParamNo()-1]; - else - tau = PMUON_LIFETIME; - - // get background - Double_t bkg; - if (fRunInfo->GetBkgFitParamNo() == -1) { // bkg not fitted - if (fRunInfo->GetBkgFix(0) == PMUSR_UNDEFINED) { // no fixed background given (background interval) - if (fRunInfo->GetBkgRange(0) >= 0) { // background range given - if (!EstimateBkg(histoNo)) - return false; - } else { // no background given to do the job, try estimate - fRunInfo->SetBkgRange(static_cast(fT0s[0]*0.1), 0); - fRunInfo->SetBkgRange(static_cast(fT0s[0]*0.6), 1); - cerr << endl << ">> PRunSingleHistoRRF::PrepareViewData(): **WARNING** Neither fix background nor background bins are given!"; - cerr << endl << ">> Will try the following: bkg start = " << fRunInfo->GetBkgRange(0) << ", bkg end = " << fRunInfo->GetBkgRange(1); - cerr << endl << ">> NO WARRANTY THAT THIS MAKES ANY SENSE! Better check ..."; - cerr << endl; - if (!EstimateBkg(histoNo)) - return false; - } - bkg = fBackground; - } else { // fixed bkg given - bkg = fRunInfo->GetBkgFix(0); - } - } else { // bkg fitted - bkg = par[fRunInfo->GetBkgFitParamNo()-1]; - } - bkg *= theoryNorm; - - 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)); - fData.SetDataTimeStep(fTimeResolution*packing); - - // data is always normalized to (per nsec!!) - 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; iGetMsrPlotList()->at(0).fRRFUnit) { - case RRF_UNIT_kHz: - gammaRRF = TMath::TwoPi()*1.0e-3; - break; - case RRF_UNIT_MHz: - gammaRRF = TMath::TwoPi(); - break; - case RRF_UNIT_Mcs: - gammaRRF = 1.0; - break; - case RRF_UNIT_G: - gammaRRF = GAMMA_BAR_MUON*TMath::TwoPi(); - break; - case RRF_UNIT_T: - gammaRRF = GAMMA_BAR_MUON*TMath::TwoPi()*1.0e4; - break; - default: - gammaRRF = TMath::TwoPi(); - break; - } - wRRF = gammaRRF * fMsrInfo->GetMsrPlotList()->at(0).fRRFFreq; - phaseRRF = fMsrInfo->GetMsrPlotList()->at(0).fRRFPhase / 180.0 * TMath::Pi(); - - Double_t error = 0.0; - for (Int_t i=start; iGetNoOfFuncs(); i++) { fFuncValues[i] = fMsrInfo->EvalFunc(fMsrInfo->GetFuncNo(i), *fRunInfo->GetMap(), par); } - // calculate theory - Double_t theoryValue; + // check if a finer binning for the theory is needed UInt_t size = fForward.size(); Double_t factor = 1.0; - UInt_t rebinRRF = 0; + Double_t time = 0.0; + Double_t theoryValue = 0.0; - if (wRRF == 0) { // no RRF - // check if a finer binning for the theory is needed - if (fData.GetValue()->size() * 10 > fForward.size()) { - size = fData.GetValue()->size() * 10; - factor = (Double_t)fForward.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 + if (fData.GetValue()->size() * 10 > fForward.size()) { + size = fData.GetValue()->size() * 10; + factor = (Double_t)fForward.size() / (Double_t)size; } + fData.SetTheoryTimeStart(fData.GetDataTimeStart()); + fData.SetTheoryTimeStep(fTimeResolution*factor); + // calculate theory 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 theory - 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) - FilterTheo(); - } - - // clean up - par.clear(); -*/ return true; } @@ -1416,7 +1041,6 @@ Double_t PRunSingleHistoRRF::EstimateN0(Double_t &errN0) Int_t endBin = N0EstimateEndTime / fTimeResolution - fGoodBins[0]; Double_t n0 = 0.0; Double_t wN = 0.0; - cout << "debug> PRunSingleHistoRRF::EstimateN0(): endBin=" << endBin << endl; for (Int_t i=0; i PRunSingleHistoRRF::EstimateN0(): N0=" << n0 << "(" << errN0 << ")" << endl; + cout << "info> PRunSingleHistoRRF::EstimateN0(): N0=" << n0 << "(" << errN0 << ")" << endl; return n0; } @@ -1472,10 +1096,10 @@ Bool_t PRunSingleHistoRRF::EstimateBkg(UInt_t histoNo) // calculate proper background range if (beamPeriod != 0.0) { - Double_t timeBkg = (Double_t)(end-start)*(fTimeResolution*fPacking); // length of the background intervall in time + Double_t timeBkg = (Double_t)(end-start)*fTimeResolution; // length of the background intervall in time UInt_t fullCycles = (UInt_t)(timeBkg/beamPeriod); // how many proton beam cylces can be placed within the proposed background intervall // correct the end of the background intervall such that the background is as close as possible to a multiple of the proton cylce - end = start + (UInt_t) ((fullCycles*beamPeriod)/(fTimeResolution*fPacking)); + end = start + (UInt_t) ((fullCycles*beamPeriod)/fTimeResolution); cout << endl << "PRunSingleHistoRRF::EstimatBkg(): Background " << start << ", " << end; if (end == start) end = fRunInfo->GetBkgRange(1); @@ -1509,7 +1133,7 @@ Bool_t PRunSingleHistoRRF::EstimateBkg(UInt_t histoNo) fBackground = bkg; // keep background (per bin) - cout << endl << "debug> fBackground=" << fBackground << endl; + cout << endl << "info> fBackground=" << fBackground << endl; fRunInfo->SetBkgEstimated(fBackground, 0); diff --git a/src/include/PMusr.h b/src/include/PMusr.h index 3374558f..efa0cf69 100644 --- a/src/include/PMusr.h +++ b/src/include/PMusr.h @@ -101,10 +101,11 @@ typedef struct { char a[7]; } __float128; // needed since cint doesn't know it //------------------------------------------------------------- // msr plot type tags -#define MSR_PLOT_SINGLE_HISTO 0 -#define MSR_PLOT_ASYM 2 -#define MSR_PLOT_MU_MINUS 4 -#define MSR_PLOT_NON_MUSR 8 +#define MSR_PLOT_SINGLE_HISTO 0 +#define MSR_PLOT_SINGLE_HISTO_RRF 1 +#define MSR_PLOT_ASYM 2 +#define MSR_PLOT_MU_MINUS 4 +#define MSR_PLOT_NON_MUSR 8 //------------------------------------------------------------- // map and fun offsets for parameter parsing diff --git a/src/include/PRunSingleHistoRRF.h b/src/include/PRunSingleHistoRRF.h index 0f60d6c7..1c09c314 100644 --- a/src/include/PRunSingleHistoRRF.h +++ b/src/include/PRunSingleHistoRRF.h @@ -55,7 +55,6 @@ class PRunSingleHistoRRF : public PRunBase virtual void CalcNoOfFitBins(); virtual Bool_t PrepareData(); virtual Bool_t PrepareFitData(PRawRunData* runData, const UInt_t histoNo); - virtual Bool_t PrepareRawViewData(PRawRunData* runData, const UInt_t histoNo); virtual Bool_t PrepareViewData(PRawRunData* runData, const UInt_t histoNo); private: @@ -63,7 +62,7 @@ class PRunSingleHistoRRF : public PRunBase UInt_t fNoOfFitBins; ///< number of bins to be fitted Double_t fBackground; ///< needed if background range is given (units: 1/bin) - Int_t fPacking; ///< packing for this particular run. Either given in the RUN- or GLOBAL-block. + Int_t fRRFPacking; ///< RRF packing for this particular run. Given in the GLOBAL-block. Int_t fGoodBins[2]; ///< keep first/last good bins. 0=fgb, 1=lgb