started with the implementation of a Fourier transform

This commit is contained in:
nemu 2008-12-15 12:27:49 +00:00
parent f353993859
commit 3af125eac7
5 changed files with 425 additions and 4 deletions

View File

@ -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<fPlots.size(); i++) {
// plot header
@ -1684,6 +1785,244 @@ bool PMsrHandler::HandleCommandsEntry(PMsrLines &lines)
return true;
}
//--------------------------------------------------------------------------
// InitFourierParameterStructure (private)
//--------------------------------------------------------------------------
/**
* <p>
*
* \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)
//--------------------------------------------------------------------------
/**
* <p>
*
* \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<TObjString*>(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<TObjString*>(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<TObjString*>(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<TObjString*>(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<TObjString*>(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<TObjString*>(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<TObjString*>(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)
//--------------------------------------------------------------------------

View File

@ -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 &param);
virtual void InitFourierParameterStructure(PMsrFourierStructure &fourier);
virtual bool FilterFunMapNumber(TString str, const char *filter, int &no);
};

View File

@ -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
//-------------------------------------------------------------
/**
* <p> typedef to make to code more readable.
@ -259,6 +279,21 @@ typedef struct {
*/
typedef vector<PMsrRunStructure> PMsrRunList;
//-------------------------------------------------------------
/**
* <p> 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;
//-------------------------------------------------------------
/**
* <p> Holds the information of a single plot block

View File

@ -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<fNoOfBins/2; i++) {
re.SetBinContent(i+1, fOut[i][0]*cos(min_phase) - fOut[i][1]*sin(min_phase));
im.SetBinContent(i+1, fOut[i][0]*sin(min_phase) + fOut[i][1]*cos(min_phase));
pwr.SetBinContent(i+1, sqrt(fOut[i][0]*fOut[i][0] + fOut[i][1]*fOut[i][1]));
}
TFile f("test_out.root", "RECREATE");
re.Write();
im.Write();
pwr.Write();
sumHist.Write();
f.Close();
}

View File

@ -0,0 +1,37 @@
void PhaseShifter(Double_t phase, TH1F *re, TH1F *im, TH1F *pwr)
{
static Double_t phaseAbs = 0.0;
phaseAbs += phase;
cout << ">> 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; i<noOfEntries-1; i++) {
reShift.SetBinContent(i, re->GetBinContent(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; i<noOfEntries-1; i++) {
re->SetBinContent(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");
}