diff --git a/src/classes/PMsrHandler.cpp b/src/classes/PMsrHandler.cpp index ab0d8255..a0123467 100644 --- a/src/classes/PMsrHandler.cpp +++ b/src/classes/PMsrHandler.cpp @@ -112,9 +112,13 @@ int PMsrHandler::ReadMsrFile() PMsrLines functions; PMsrLines run; PMsrLines commands; + PMsrLines fourier; PMsrLines plot; PMsrLines statistic; + // init stuff + InitFourierParameterStructure(fFourier); + // open msr-file f.open(fFileName.Data(), iostream::in); if (!f.is_open()) { @@ -157,6 +161,9 @@ int PMsrHandler::ReadMsrFile() } else if (line.Contains("COMMANDS")) { // COMMANDS block tag fMsrBlockCounter = MSR_TAG_COMMANDS; commands.push_back(current); + } else if (line.Contains("FOURIER")) { // FOURIER block tag + fMsrBlockCounter = MSR_TAG_FOURIER; + fourier.push_back(current); } else if (line.Contains("PLOT")) { // PLOT block tag fMsrBlockCounter = MSR_TAG_PLOT; plot.push_back(current); @@ -181,6 +188,9 @@ int PMsrHandler::ReadMsrFile() case MSR_TAG_COMMANDS: // COMMANDS block commands.push_back(current); break; + case MSR_TAG_FOURIER: // FOURIER block + fourier.push_back(current); + break; case MSR_TAG_PLOT: // PLOT block plot.push_back(current); break; @@ -212,6 +222,9 @@ int PMsrHandler::ReadMsrFile() if (result == PMUSR_SUCCESS) if (!HandleCommandsEntry(commands)) result = PMUSR_MSR_SYNTAX_ERROR; + if (result == PMUSR_SUCCESS) + if (!HandleFourierEntry(fourier)) + result = PMUSR_MSR_SYNTAX_ERROR; if (result == PMUSR_SUCCESS) if (!HandlePlotEntry(plot)) result = PMUSR_MSR_SYNTAX_ERROR; @@ -252,6 +265,7 @@ int PMsrHandler::ReadMsrFile() functions.clear(); run.clear(); commands.clear(); + fourier.clear(); plot.clear(); statistic.clear(); @@ -261,6 +275,17 @@ int PMsrHandler::ReadMsrFile() // } // cout << endl; +cout << endl << ">> FOURIER Block:"; +cout << endl << ">> Fourier Block Present : " << fFourier.fFourierBlockPresent; +cout << endl << ">> Fourier Units : " << fFourier.fUnits; +cout << endl << ">> Fourier Power : " << fFourier.fFourierPower; +cout << endl << ">> Apodization : " << fFourier.fApodization; +cout << endl << ">> Plot Tag : " << fFourier.fPlotTag; +cout << endl << ">> Phase : " << fFourier.fPhase; +cout << endl << ">> Range for Freq. Corrections : " << fFourier.fRangeForPhaseCorrection[0] << ", " << fFourier.fRangeForPhaseCorrection[1]; +cout << endl << ">> Plot Range : " << fFourier.fPlotRange[0] << ", " << fFourier.fPlotRange[1]; +cout << endl; + return result; } @@ -642,6 +667,82 @@ int PMsrHandler::WriteMsrLogFile() CheckAndWriteComment(f, ++lineNo); } + // write fourier block if present + if (fFourier.fFourierBlockPresent) { + f << endl << "FOURIER"; + CheckAndWriteComment(f, ++lineNo); + + // write 'units' parameter if present + if (fFourier.fUnits != FOURIER_UNIT_NOT_GIVEN) { + f << endl << "units "; + if (fFourier.fUnits == FOURIER_UNIT_FIELD) { + f << "Gauss"; + } else if (fFourier.fUnits == FOURIER_UNIT_FREQ) { + f << "MHz "; + } + f << " # units either 'Gauss' or 'MHz'"; + CheckAndWriteComment(f, ++lineNo); + } + + // write 'fourier_power' parameter if present + if (fFourier.fFourierPower != 0) { + f << endl << "fourier_power " << fFourier.fFourierPower; + CheckAndWriteComment(f, ++lineNo); + } + + // write 'appodization' if present + if (fFourier.fApodization != FOURIER_APOD_NOT_GIVEN) { + f << endl << "apodization "; + if (fFourier.fApodization == FOURIER_APOD_NONE) { + f << "NONE "; + } else if (fFourier.fApodization == FOURIER_APOD_WEAK) { + f << "WEAK "; + } else if (fFourier.fApodization == FOURIER_APOD_MEDIUM) { + f << "MEDIUM"; + } else if (fFourier.fApodization == FOURIER_APOD_STRONG) { + f << "STRONG"; + } + f << " # NONE, WEAK, MEDIUM, STRONG"; + CheckAndWriteComment(f, ++lineNo); + } + + // write 'plot' if present + if (fFourier.fPlotTag != FOURIER_PLOT_NOT_GIVEN) { + f << endl << "plot "; + if (fFourier.fPlotTag == FOURIER_PLOT_REAL) { + f << "REAL "; + } else if (fFourier.fPlotTag == FOURIER_PLOT_IMAG) { + f << "IMAG "; + } else if (fFourier.fPlotTag == FOURIER_PLOT_REAL_AND_IMAG) { + f << "REAL_AND_IMAG"; + } else if (fFourier.fPlotTag == FOURIER_PLOT_POWER) { + f << "POWER"; + } else if (fFourier.fPlotTag == FOURIER_PLOT_PHASE) { + f << "PHASE"; + } + f << " # REAL, IMAG, REAL_AND_IMAG, POWER, PHASE"; + CheckAndWriteComment(f, ++lineNo); + } + + // write 'phase' if present + if (fFourier.fPhase != -999.0) { + f << endl << "phase " << fFourier.fPhase; + CheckAndWriteComment(f, ++lineNo); + } + + // write 'range_for_phase_correction' if present + if (fFourier.fRangeForPhaseCorrection[0] != -1.0) { + f << endl << "range_for_phase_correction " << fFourier.fRangeForPhaseCorrection[0] << " " << fFourier.fRangeForPhaseCorrection[1]; + CheckAndWriteComment(f, ++lineNo); + } + + // write 'range' if present + if (fFourier.fPlotRange[0] != -1.0) { + f << endl << "range " << fFourier.fPlotRange[0] << " " << fFourier.fPlotRange[1]; + CheckAndWriteComment(f, ++lineNo); + } + } + // write plot block for (unsigned int i=0; i + * + * \param fourier fourier parameters + */ +void PMsrHandler::InitFourierParameterStructure(PMsrFourierStructure &fourier) +{ + fourier.fFourierBlockPresent = false; // fourier block present + fourier.fUnits = FOURIER_UNIT_NOT_GIVEN; // fourier untis in field, i.e. Gauss + fourier.fFourierPower = 0; // no zero padding + fourier.fApodization = FOURIER_APOD_NOT_GIVEN; // no apodization + fourier.fPlotTag = FOURIER_PLOT_NOT_GIVEN; // initial plot tag: show real and imaginary part at once + fourier.fPhase = -999.0; // fourier phase + for (unsigned int i=0; i<2; i++) { + fourier.fRangeForPhaseCorrection[i] = -1.0; // frequency range for phase correction + fourier.fPlotRange[i] = -1.0; // fourier plot range + } +} + +//-------------------------------------------------------------------------- +// HandleFourierEntry (private) +//-------------------------------------------------------------------------- +/** + *

+ * + * \param lines is a list of lines containing the fourier parameter block + */ +bool PMsrHandler::HandleFourierEntry(PMsrLines &lines) +{ +cout << endl << ">> in PMsrHandler::HandleFourierEntry ..."; + + bool error = false; + + if (lines.empty()) // no fourier block present + return true; + +cout << endl << ">> in PMsrHandler::HandleFourierEntry, Fourier block present ..."; + + PMsrFourierStructure fourier; + + InitFourierParameterStructure(fourier); + + fourier.fFourierBlockPresent = true; + + PMsrLines::iterator iter; + + TObjArray *tokens; + TObjString *ostr; + TString str; + + int ival; + + iter = lines.begin(); + while ((iter != lines.end()) && !error) { + // tokenize line + tokens = iter->fLine.Tokenize(" \t"); + if (!tokens) { + cout << endl << ">> PMsrHandler::HandleRunEntry: **SEVERE ERROR** Couldn't tokenize Parameters in line " << iter->fLineNo; + cout << endl << endl; + return false; + } + + // units ----------------------------------------------- + if (iter->fLine.BeginsWith("units", TString::kIgnoreCase)) { + if (tokens->GetEntries() < 2) { // units are missing + error = true; + continue; + } else { + ostr = dynamic_cast(tokens->At(1)); + str = ostr->GetString(); + if (!str.CompareTo("gauss", TString::kIgnoreCase)) { + fourier.fUnits = FOURIER_UNIT_FIELD; + } else if (!str.CompareTo("mhz", TString::kIgnoreCase)) { + fourier.fUnits = FOURIER_UNIT_FREQ; + } else { + error = true; + continue; + } + } + } + + // fourier power (zero padding) ------------------------ + if (iter->fLine.BeginsWith("fourier_power", TString::kIgnoreCase)) { + if (tokens->GetEntries() < 2) { // fourier power exponent is missing + error = true; + continue; + } else { + ostr = dynamic_cast(tokens->At(1)); + str = ostr->GetString(); + if (str.IsDigit()) { + ival = str.Atoi(); + if ((ival >= 0) && (ival <= 20)) { + fourier.fFourierPower = ival; + } else { // fourier power out of range + error = true; + continue; + } + } else { // fourier power not a number + error = true; + continue; + } + } + } + + // apodization ----------------------------------------- + if (iter->fLine.BeginsWith("apodization", TString::kIgnoreCase)) { + if (tokens->GetEntries() < 2) { // apodization tag is missing + error = true; + continue; + } else { + ostr = dynamic_cast(tokens->At(1)); + str = ostr->GetString(); + if (!str.CompareTo("none", TString::kIgnoreCase)) { + fourier.fApodization = FOURIER_APOD_NONE; + } else if (!str.CompareTo("weak", TString::kIgnoreCase)) { + fourier.fApodization = FOURIER_APOD_WEAK; + } else if (!str.CompareTo("medium", TString::kIgnoreCase)) { + fourier.fApodization = FOURIER_APOD_MEDIUM; + } else if (!str.CompareTo("strong", TString::kIgnoreCase)) { + fourier.fApodization = FOURIER_APOD_STRONG; + } else { // unrecognized apodization tag + error = true; + continue; + } + } + } + + // plot tag -------------------------------------------- + if (iter->fLine.BeginsWith("plot", TString::kIgnoreCase)) { + if (tokens->GetEntries() < 2) { // plot tag is missing + error = true; + continue; + } else { + ostr = dynamic_cast(tokens->At(1)); + str = ostr->GetString(); + if (!str.CompareTo("real", TString::kIgnoreCase)) { + fourier.fPlotTag = FOURIER_PLOT_REAL; + } else if (!str.CompareTo("imag", TString::kIgnoreCase)) { + fourier.fPlotTag = FOURIER_PLOT_IMAG; + } else if (!str.CompareTo("real_and_imag", TString::kIgnoreCase)) { + fourier.fPlotTag = FOURIER_PLOT_REAL_AND_IMAG; + } else if (!str.CompareTo("power", TString::kIgnoreCase)) { + fourier.fPlotTag = FOURIER_PLOT_POWER; + } else if (!str.CompareTo("phase", TString::kIgnoreCase)) { + fourier.fPlotTag = FOURIER_PLOT_PHASE; + } else { // unrecognized plot tag + error = true; + continue; + } + } + } + + // phase ----------------------------------------------- + if (iter->fLine.BeginsWith("phase", TString::kIgnoreCase)) { + if (tokens->GetEntries() < 2) { // phase value is missing + error = true; + continue; + } else { + ostr = dynamic_cast(tokens->At(1)); + str = ostr->GetString(); + if (str.IsFloat()) { + fourier.fPhase = str.Atof(); + } else { + error = true; + continue; + } + } + } + + // range for automatic phase correction ---------------- + if (iter->fLine.BeginsWith("range_for_phase_correction", TString::kIgnoreCase)) { + if (tokens->GetEntries() < 3) { // range values are missing + error = true; + continue; + } else { + for (unsigned int i=0; i<2; i++) { + ostr = dynamic_cast(tokens->At(i+1)); + str = ostr->GetString(); + if (str.IsFloat()) { + fourier.fRangeForPhaseCorrection[i] = str.Atof(); + } else { + error = true; + continue; + } + } + } + } + + // fourier plot range ---------------------------------- + if (iter->fLine.BeginsWith("range", TString::kIgnoreCase)) { + if (tokens->GetEntries() < 3) { // plot range values are missing + error = true; + continue; + } else { + for (unsigned int i=0; i<2; i++) { + ostr = dynamic_cast(tokens->At(i+1)); + str = ostr->GetString(); + if (str.IsFloat()) { + fourier.fPlotRange[i] = str.Atof(); + } else { + error = true; + continue; + } + } + } + } + + ++iter; + } + + if (error) { + cout << endl << ">> PMsrHandler::HandleFourierEntry: **ERROR** in line " << iter->fLineNo << ":"; + cout << endl; + cout << endl << iter->fLine.Data(); + cout << endl; + cout << endl << "FOURIER block syntax, parameters in [] are optinal:"; + cout << endl; + cout << endl << "FOURIER"; + cout << endl << "[units Gauss | MHz]"; + cout << endl << "[fourier_power n # n is a number such that zero padding up to 2^n will be used]"; + cout << endl << " n=0 means no zero padding"; + cout << endl << " 0 <= n <= 20 are allowed values"; + cout << endl << "[apodization none | weak | medium | strong]"; + cout << endl << "[plot real | imag | real_and_imag | power | phase]"; + cout << endl << "[phase value]"; + cout << endl << "[range_for_phase_correction min max]"; + cout << endl << "[range min max]"; + cout << endl; + } else { // save last run found + fFourier = fourier; + } + + return !error; +} + //-------------------------------------------------------------------------- // HandlePlotEntry (private) //-------------------------------------------------------------------------- diff --git a/src/include/PMsrHandler.h b/src/include/PMsrHandler.h index 98174905..168256b9 100644 --- a/src/include/PMsrHandler.h +++ b/src/include/PMsrHandler.h @@ -59,6 +59,7 @@ class PMsrHandler virtual PMsrLines* GetMsrFunctions() { return &fFunctions; } virtual PMsrRunList* GetMsrRunList() { return &fRuns; } virtual PMsrLines* GetMsrCommands() { return &fCommands; } + virtual PMsrFourierStructure* GetMsrFourierList() { return &fFourier; } virtual PMsrPlotList* GetMsrPlotList() { return &fPlots; } virtual PMsrStatisticStructure* GetMsrStatistic() { return &fStatistic; } @@ -95,6 +96,7 @@ class PMsrHandler PMsrLines fFunctions; ///< holds the user defined functions PMsrRunList fRuns; ///< holds a list of run information PMsrLines fCommands; ///< holds a list of the minuit commands + PMsrFourierStructure fFourier; ///< holds the parameters used for the Fourier transform PMsrPlotList fPlots; ///< holds a list of the plot input parameters PMsrStatisticStructure fStatistic; ///< holds the statistic info @@ -109,6 +111,7 @@ class PMsrHandler virtual bool HandleFunctionsEntry(PMsrLines &line); virtual bool HandleRunEntry(PMsrLines &line); virtual bool HandleCommandsEntry(PMsrLines &line); + virtual bool HandleFourierEntry(PMsrLines &line); virtual bool HandlePlotEntry(PMsrLines &line); virtual bool HandleStatisticEntry(PMsrLines &line); @@ -117,6 +120,8 @@ class PMsrHandler virtual void FillParameterInUse(PMsrLines &theory, PMsrLines &funcs, PMsrLines &run); virtual void InitRunParameterStructure(PMsrRunStructure ¶m); + virtual void InitFourierParameterStructure(PMsrFourierStructure &fourier); + virtual bool FilterFunMapNumber(TString str, const char *filter, int &no); }; diff --git a/src/include/PMusr.h b/src/include/PMusr.h index cb85cf11..fdbca2c7 100644 --- a/src/include/PMusr.h +++ b/src/include/PMusr.h @@ -70,8 +70,9 @@ using namespace std; #define MSR_TAG_FUNCTIONS 3 #define MSR_TAG_RUN 4 #define MSR_TAG_COMMANDS 5 -#define MSR_TAG_PLOT 6 -#define MSR_TAG_STATISTIC 7 +#define MSR_TAG_FOURIER 6 +#define MSR_TAG_PLOT 7 +#define MSR_TAG_STATISTIC 8 //------------------------------------------------------------- // msr fit type tags @@ -92,6 +93,25 @@ using namespace std; #define MSR_PARAM_MAP_OFFSET 10000 #define MSR_PARAM_FUN_OFFSET 20000 +//------------------------------------------------------------- +// fourier related tags +#define FOURIER_UNIT_NOT_GIVEN 0 +#define FOURIER_UNIT_FIELD 1 +#define FOURIER_UNIT_FREQ 2 + +#define FOURIER_APOD_NOT_GIVEN 0 +#define FOURIER_APOD_NONE 1 +#define FOURIER_APOD_WEAK 2 +#define FOURIER_APOD_MEDIUM 3 +#define FOURIER_APOD_STRONG 4 + +#define FOURIER_PLOT_NOT_GIVEN 0 +#define FOURIER_PLOT_REAL 1 +#define FOURIER_PLOT_IMAG 2 +#define FOURIER_PLOT_REAL_AND_IMAG 3 +#define FOURIER_PLOT_POWER 4 +#define FOURIER_PLOT_PHASE 5 + //------------------------------------------------------------- /** *

typedef to make to code more readable. @@ -259,6 +279,21 @@ typedef struct { */ typedef vector PMsrRunList; +//------------------------------------------------------------- +/** + *

Holds the information of the Fourier block + */ +typedef struct { + bool fFourierBlockPresent; ///< flag indicating if a Fourier block is present in the msr-file + bool fUnits; ///< flag used to indicate the units. 0=field units (G); 1=frequency units (MHz) + int fFourierPower; ///< i.e. zero padding up to 2^fFourierPower, default = 0 which means NO zero padding + int fApodization; ///< tag indicating the kind of apodization wished, 0=no appodization (default), 1=weak, 2=medium, 3=strong (for details see the docu) + int fPlotTag; ///< tag used for initial plot. 0=real, 1=imaginary, 2=real & imaginary (default), 3=power, 4=phase + double fPhase; ///< phase + double fRangeForPhaseCorrection[2]; ///< field/frequency range for automatic phase correction + double fPlotRange[2]; ///< field/frequency plot range +} PMsrFourierStructure; + //------------------------------------------------------------- /** *

Holds the information of a single plot block diff --git a/src/tests/fourier/PMusrFourier.cpp b/src/tests/fourier/PMusrFourier.cpp index 2c317cfb..a64a3bd1 100644 --- a/src/tests/fourier/PMusrFourier.cpp +++ b/src/tests/fourier/PMusrFourier.cpp @@ -202,12 +202,14 @@ void PMusrFourier::Transform(unsigned int apodizationTag, unsigned int filterTag // loop over the phase double sum; double corr_phase; -double min, min_phase;; +double min, min_phase; TH1F sumHist("sumHist", "sumHist", 181, -90.5, 90.5); double dB = 1.0/(2.0 * F_GAMMA_BAR_MUON * (fEndTime-fStartTime)); double Bmax = 1.0/(2.0 * F_GAMMA_BAR_MUON * fTimeResolution); TH1F re("re", "re", fNoOfBins/2+1, -dB/2.0, Bmax+dB/2.0); TH1F im("im", "im", fNoOfBins/2+1, -dB/2.0, Bmax+dB/2.0); +TH1F pwr("pwr", "pwr", fNoOfBins/2+1, -dB/2.0, Bmax+dB/2.0); +/* for (int p=-90; p<90; p++) { // calculate sum of the rotated imaginary part of fOut sum = 0.0; @@ -229,15 +231,18 @@ cout << endl << "-> min = " << min << ", min_phase = " << min_phase/PI*180.0; sumHist.SetBinContent(p+91, fabs(sum)); } cout << endl << ">> min = " << min << ", min_phase = " << min_phase/PI*180.0; - +*/ +min_phase = 0.0; for (unsigned int i=0; i> phase = " << phaseAbs << endl; + + Int_t noOfEntries = re->GetEntries(); + Double_t min = re->GetBinLowEdge(1); + Int_t maxBin = re->GetNbinsX()-1; + Double_t max = re->GetBinLowEdge(maxBin)+re->GetBinWidth(maxBin); + TH1F reShift("reShift", "reShift", noOfEntries, min, max); + TH1F imShift("imShift", "imShift", noOfEntries, min, max); + + Double_t phaseRad = phase/180.0*TMath::Pi(); + for (Int_t i=1; iGetBinContent(i)*TMath::Cos(phaseRad) - im->GetBinContent(i)*TMath::Sin(phaseRad)); + imShift.SetBinContent(i, im->GetBinContent(i)*TMath::Cos(phaseRad) + re->GetBinContent(i)*TMath::Sin(phaseRad)); + } + + for (Int_t i=1; iSetBinContent(i, reShift.GetBinContent(i)); + im->SetBinContent(i, imShift.GetBinContent(i)); + } + + re->SetMarkerStyle(20); + re->SetMarkerColor(TColor::GetColor(0,0,0)); + re->SetLineColor(TColor::GetColor(0,0,0)); + im->SetMarkerStyle(21); + im->SetMarkerColor(TColor::GetColor(255,0,0)); + im->SetLineColor(TColor::GetColor(255,0,0)); + pwr->SetMarkerStyle(22); + pwr->SetMarkerColor(TColor::GetColor(0,0,255)); + pwr->SetLineColor(TColor::GetColor(0,0,255)); + re->Draw("alp"); + im->Draw("plsame"); + pwr->Draw("plsame"); +}