LMU/musrfit
5
0
Fork 0
Code Issues 6 Pull Requests Actions Packages Projects Releases 4 Activity

implemented an average-per-data-set option for musrFT.

Browse Source
This commit is contained in:
suter_a
2015-02-23 16:30:14 +01:00
parent 80b3952649
commit c938589286
22 changed files with 941 additions and 515 deletions
Download Patch File Download Diff File Expand all files Collapse all files

967
src/classes/PFourierCanvas.cpp
View File

File diff suppressed because it is too large Load Diff

62
src/classes/PPrepFourier.cpp
View File

@ -79,7 +79,7 @@ PPrepFourier::~PPrepFourier()
*/
void PPrepFourier::SetBkgRange(const Int_t *bkgRange)
{
int err=0;
Int_t err=0;
if (bkgRange[0] >= -1) {
fBkgRange[0] = bkgRange[0];
} else {
@ -125,7 +125,7 @@ void PPrepFourier::SetBkgRange(const Int_t *bkgRange)
*/
void PPrepFourier::SetBkg(PDoubleVector bkg)
{
for (unsigned int i=0; i<bkg.size(); i++)
for (UInt_t i=0; i<bkg.size(); i++)
fBkg.push_back(bkg[i]);
}
@ -180,7 +180,7 @@ void PPrepFourier::DoBkgCorrection()
if ((fBkgRange[0] != -1) && (fBkgRange[1] != -1)) { // background range is given
// make sure that the bkg range is ok
for (unsigned int i=0; i<fRawData.size(); i++) {
for (UInt_t i=0; i<fRawData.size(); i++) {
if ((fBkgRange[0] >= fRawData[i].rawData.size()) || (fBkgRange[1] >= fRawData[i].rawData.size())) {
cerr << endl << "PPrepFourier::DoBkgCorrection() **ERROR** bkg-range out of data-range!";
return;
@ -188,16 +188,16 @@ void PPrepFourier::DoBkgCorrection()
}
Double_t bkg=0.0;
for (unsigned int i=0; i<fRawData.size(); i++) {
for (UInt_t i=0; i<fRawData.size(); i++) {
// calculate the bkg for the given range
for (int j=fBkgRange[0]; j<=fBkgRange[1]; j++) {
for (Int_t j=fBkgRange[0]; j<=fBkgRange[1]; j++) {
bkg += fRawData[i].rawData[j];
}
bkg /= (fBkgRange[1]-fBkgRange[0]+1);
cout << "debug> background " << i << ": " << bkg << endl;
cout << "info> background " << i << ": " << bkg << endl;
// correct data
for (unsigned int j=0; j<fData[i].size(); j++)
for (UInt_t j=0; j<fData[i].size(); j++)
fData[i][j] -= bkg;
}
} else { // there might be an explicit background list
@ -211,8 +211,8 @@ void PPrepFourier::DoBkgCorrection()
return;
}
for (unsigned int i=0; i<fData.size(); i++)
for (unsigned int j=0; j<fData[i].size(); j++)
for (UInt_t i=0; i<fData.size(); i++)
for (UInt_t j=0; j<fData[i].size(); j++)
fData[i][j] -= fBkg[i];
}
}
@ -235,10 +235,10 @@ void PPrepFourier::DoPacking()
PDoubleVector tmpData;
Double_t dval = 0.0;
for (unsigned int i=0; i<fData.size(); i++) {
for (UInt_t i=0; i<fData.size(); i++) {
tmpData.clear();
dval = 0.0;
for (unsigned int j=0; j<fData[i].size(); j++) {
for (UInt_t j=0; j<fData[i].size(); j++) {
if ((j % fPacking == 0) && (j != 0)) {
tmpData.push_back(dval);
dval = 0.0;
@ -271,9 +271,9 @@ void PPrepFourier::DoLifeTimeCorrection(Double_t fudge)
// calc exp(+t/tau)*N(t), where N(t) is already background corrected
Double_t scale;
for (unsigned int i=0; i<fData.size(); i++) {
for (UInt_t i=0; i<fData.size(); i++) {
scale = fRawData[i].timeResolution / PMUON_LIFETIME;
for (unsigned int j=0; j<fData[i].size(); j++) {
for (UInt_t j=0; j<fData[i].size(); j++) {
fData[i][j] *= exp(j*scale);
}
}
@ -281,14 +281,14 @@ void PPrepFourier::DoLifeTimeCorrection(Double_t fudge)
// calc N0
Double_t dval;
Double_t N0;
for (unsigned int i=0; i<fData.size(); i++) {
for (UInt_t i=0; i<fData.size(); i++) {
dval = 0.0;
for (unsigned int j=0; j<fData[i].size(); j++) {
for (UInt_t j=0; j<fData[i].size(); j++) {
dval += fData[i][j];
}
N0 = dval/fData[i].size();
N0 *= fudge;
for (unsigned int j=0; j<fData[i].size(); j++) {
for (UInt_t j=0; j<fData[i].size(); j++) {
fData[i][j] -= N0;
fData[i][j] /= N0;
}
@ -313,6 +313,24 @@ TString PPrepFourier::GetInfo(const UInt_t idx)
return info;
}
//--------------------------------------------------------------------------
// GetDataSetTag
//--------------------------------------------------------------------------
/**
* <p>Returns the data set tag of the object
*
* \param idx index of the object
*/
Int_t PPrepFourier::GetDataSetTag(const UInt_t idx)
{
Int_t result = -1;
if (idx < fRawData.size())
result = fRawData[idx].dataSetTag;
return result;
}
//--------------------------------------------------------------------------
// GetData
//--------------------------------------------------------------------------
@ -337,7 +355,7 @@ vector<TH1F*> PPrepFourier::GetData()
UInt_t startIdx;
UInt_t endIdx;
for (unsigned int i=0; i<fData.size(); i++) {
for (UInt_t i=0; i<fData.size(); i++) {
name = TString::Format("histo%2d", i);
dt = fRawData[i].timeResolution*fPacking;
start = fRawData[i].timeRange[0];
@ -373,7 +391,7 @@ vector<TH1F*> PPrepFourier::GetData()
}
data[i] = new TH1F(name.Data(), fRawData[i].info.Data(), size, start, end);
for (unsigned int j=startIdx; j<endIdx; j++)
for (UInt_t j=startIdx; j<endIdx; j++)
data[i]->SetBinContent(j-startIdx+1, fData[i][j]);
}
@ -430,7 +448,7 @@ TH1F *PPrepFourier::GetData(const UInt_t idx)
}
TH1F *data = new TH1F(name.Data(), fRawData[idx].info.Data(), size, start, end);
for (unsigned int i=startIdx; i<endIdx; i++)
for (UInt_t i=startIdx; i<endIdx; i++)
data->SetBinContent(i-startIdx+1, fData[idx][i]);
return data;
@ -445,13 +463,13 @@ TH1F *PPrepFourier::GetData(const UInt_t idx)
void PPrepFourier::InitData()
{
fData.resize(fRawData.size());
unsigned int t0;
for (unsigned int i=0; i<fRawData.size(); i++) {
UInt_t t0;
for (UInt_t i=0; i<fRawData.size(); i++) {
if (fRawData[i].t0 >= 0)
t0 = fRawData[i].t0;
else
t0 = 0;
for (unsigned int j=t0; j<fRawData[i].rawData.size(); j++) {
for (UInt_t j=t0; j<fRawData[i].rawData.size(); j++) {
fData[i].push_back(fRawData[i].rawData[j]);
}
}

21
src/include/PFourierCanvas.h
View File

@ -38,9 +38,10 @@
#include "PFourier.h"
// Canvas menu id's
#define P_MENU_ID_FOURIER 10001
#define P_MENU_ID_AVERAGE 10002
#define P_MENU_ID_EXPORT_DATA 10003
#define P_MENU_ID_FOURIER 10001
#define P_MENU_ID_AVERAGE 10002
#define P_MENU_ID_AVERAGE_PER_DATA_SET 10003
#define P_MENU_ID_EXPORT_DATA 10004
#define P_MENU_ID_FOURIER_REAL 100
#define P_MENU_ID_FOURIER_IMAG 101
@ -75,10 +76,12 @@ class PFourierCanvas : public TObject, public TQObject
{
public:
PFourierCanvas();
PFourierCanvas(vector<PFourier*> &fourier, const Char_t* title, const Bool_t showAverage,
PFourierCanvas(vector<PFourier*> &fourier, PIntVector dataSetTag, const Char_t* title,
const Bool_t showAverage, const Bool_t showAveragePerDataSet,
const Int_t fourierPlotOpt, Double_t fourierXrange[2], Double_t phase,
Int_t wtopx, Int_t wtopy, Int_t ww, Int_t wh, const Bool_t batch);
PFourierCanvas(vector<PFourier*> &fourier, const Char_t* title, const Bool_t showAverage,
PFourierCanvas(vector<PFourier*> &fourier, PIntVector dataSetTag, const Char_t* title,
const Bool_t showAverage, const Bool_t showAveragePerDataSet,
const Int_t fourierPlotOpt, Double_t fourierXrange[2], Double_t phase,
Int_t wtopx, Int_t wtopy, Int_t ww, Int_t wh,
const PIntVector markerList, const PIntVector colorList, const Bool_t batch);
@ -103,7 +106,9 @@ class PFourierCanvas : public TObject, public TQObject
Int_t fTimeout; ///< timeout after which the Done signal should be emited. If timeout <= 0, no timeout is taking place
Bool_t fBatchMode; ///< musrview in ROOT batch mode
Bool_t fValid; ///< if true, everything looks OK
Bool_t fAveragedView; ///< tag showing that the averaged view or normal view should be presented.
Bool_t fAveragedView; ///< tag showing that the averaged view for ALL data or normal view should be presented.
Bool_t fAveragedViewPerDataSet; ///< tag showing that the averaged view for individual data sets or normal view should be presented.
PIntVector fDataSetTag; ///< vector holding the data set tags
Int_t fCurrentPlotView; ///< tag showing what the current plot view is: real, imag, power, phase, ...
Double_t fInitialXRange[2]; ///< keeps the initial x-range
Double_t fInitialYRange[2]; ///< keeps the initial y-range
@ -112,7 +117,7 @@ class PFourierCanvas : public TObject, public TQObject
TString fXaxisTitle;
vector<PFourier*> fFourier; ///< keeps all the Fourier data, ownership is with the caller
PFourierCanvasDataList fFourierHistos; ///< keeps all the Fourier histos
PFourierCanvasDataSet fFourierAverage; ///< keeps the average of the Fourier histos
PFourierCanvasDataList fFourierAverage; ///< keeps the average of the Fourier histos
Double_t fCurrentFourierPhase; ///< keeps the current Fourier phase (real/imag)
TLatex *fCurrentFourierPhaseText; ///< used in Re/Im Fourier to show the current phase in the pad
@ -134,6 +139,7 @@ class PFourierCanvas : public TObject, public TQObject
TPaveText *fTitlePad; ///< title pad used to display a title
TPad *fFourierPad; ///< fourier pad used to display the fourier
TLegend *fInfoPad; ///< info pad used to display a legend of the data plotted
TLegend *fLegAvgPerDataSet; ///< legend used for averaged per data set view
virtual void CreateXaxisTitle();
virtual void CreateStyle();
@ -151,6 +157,7 @@ class PFourierCanvas : public TObject, public TQObject
virtual Double_t GetMaximum(TH1F* histo, Double_t xmin=-1.0, Double_t xmax=-1.0);
virtual Double_t GetMinimum(TH1F* histo, Double_t xmin=-1.0, Double_t xmax=-1.0);
virtual Double_t GetInterpolatedValue(TH1F* histo, Double_t xVal);
virtual TString GetDataSetName(TString title);
ClassDef(PFourierCanvas, 1)
};

12
src/include/PPrepFourier.h
View File

@ -45,11 +45,12 @@ using namespace std;
* necessary meta information.
*/
typedef struct {
TString info; ///< keeps all the meta information
double timeResolution; ///< time resolution in (usec)
int t0; ///< keep the t0 bin
Double_t timeRange[2]; ///< time range to be used, given in (usec).
PDoubleVector rawData; ///< a single time domain data vector
Int_t dataSetTag; ///< tag to label the data set. Needed for average-per-data-set
TString info; ///< keeps all the meta information
Double_t timeResolution; ///< time resolution in (usec)
Int_t t0; ///< keep the t0 bin
Double_t timeRange[2]; ///< time range to be used, given in (usec).
PDoubleVector rawData; ///< a single time domain data vector
} musrFT_data;
//----------------------------------------------------------------------------
@ -72,6 +73,7 @@ class PPrepFourier {
virtual void DoLifeTimeCorrection(Double_t fudge);
TString GetInfo(const UInt_t idx);
Int_t GetDataSetTag(const UInt_t idx);
UInt_t GetNoOfData() { return fRawData.size(); }
vector<TH1F*> GetData();
TH1F *GetData(const UInt_t idx);

201
src/musrFT.cpp
View File

@ -59,27 +59,28 @@ using namespace std;
* <p>Structure keeping the command line options.
*/
typedef struct {
vector<TString> msrFln; ///< msr-file names to be used.
vector<TString> dataFln; ///< raw-data-file names to be used.
vector<TString> dataFileFormat; ///< file format guess
PStringVector msrFln; ///< msr-file names to be used.
PStringVector dataFln; ///< raw-data-file names to be used.
PStringVector dataFileFormat; ///< file format guess
TString graphicFormat; ///< format for the graphical output dump
TString dumpFln; ///< dump file name for Fourier data output
TString msrFlnOut; ///< dump file name for msr-file generation
int bkg_range[2]; ///< background range
vector<double> bkg; ///< background value
Int_t bkg_range[2]; ///< background range
PDoubleVector bkg; ///< background value
TString fourierOpt; ///< Fourier options, i.e. real, imag, power, phase
TString apodization; ///< apodization setting: none, weak, medium, strong
int fourierPower; ///< Fourier power for zero padding, i.e. 2^fourierPower points
Int_t fourierPower; ///< Fourier power for zero padding, i.e. 2^fourierPower points
TString fourierUnits; ///< wished Fourier units: Gauss, Tesla, MHz, Mc/s
double initialPhase; ///< inital Fourier phase for Real/Imag
double fourierRange[2]; ///< Fourier range to be plotted. Given in the choosen units.
double timeRange[2]; ///< time range used for the Fourier
vector<int> histo; ///< selection of the histos used from at data file for Fourier
bool showAverage; ///< flag indicating if initially the Fourier average over the given histos shall be plotted.
vector<int> t0; ///< t0 vector for the histos. If not given t0's will be estimated.
int packing; ///< packing for rebinning the time histograms before Fourier transform.
Double_t initialPhase; ///< inital Fourier phase for Real/Imag
Double_t fourierRange[2]; ///< Fourier range to be plotted. Given in the choosen units.
Double_t timeRange[2]; ///< time range used for the Fourier
PIntVector histo; ///< selection of the histos used from at data file for Fourier
Bool_t showAverage; ///< flag indicating if initially the Fourier average over the given histos shall be plotted, this over ALL data sets.
Bool_t showAveragePerDataSet; ///< flag indicating if initially the Fourier average over the given histos shall be plotted, this per data set.
PIntVector t0; ///< t0 vector for the histos. If not given t0's will be estimated.
Int_t packing; ///< packing for rebinning the time histograms before Fourier transform.
TString title; ///< title to be shown for the Fourier plot.
double lifetimecorrection; ///< is == 0.0 for NO life time correction, otherwise it holds the fudge factor
Double_t lifetimecorrection; ///< is == 0.0 for NO life time correction, otherwise it holds the fudge factor
Int_t timeout; ///< timeout in (sec) after which musrFT will terminate. if <= 0, no automatic termination will take place.
} musrFT_startup_param;
@ -134,7 +135,8 @@ void musrFT_syntax()
cout << endl << " the option --data-file. If multiple data file are given, <list> will apply";
cout << endl << " to all data-files given. If --histo is not given, all histos of a data file will be used.";
cout << endl << " <list> can be anything like: 2 3 6, or 2-17, or 1-6 9, etc.";
cout << endl << " -a, --average : show the average of all Fourier transformed data.";
cout << endl << " -a, --average : show the average of ALL Fourier transformed data.";
cout << endl << " -ad, --average-per-data-set : show the average of the Fourier transformed data per data set.";
cout << endl << " --t0 <list> : A list of t0's can be provided. This in conjunction with --data-file and";
cout << endl << " --fourier-option real allows to get the proper initial phase if t0's are known.";
cout << endl << " If a single t0 for multiple histos is given, it is assume, that this t0 is common";
@ -174,6 +176,7 @@ void musrFT_init(musrFT_startup_param &startupParam)
startupParam.timeRange[0] = -1.0;
startupParam.timeRange[1] = -1.0;
startupParam.showAverage = false;
startupParam.showAveragePerDataSet = false;
startupParam.packing = 1;
startupParam.title = TString("");
startupParam.lifetimecorrection = 0.0;
@ -196,9 +199,9 @@ void musrFT_init(musrFT_startup_param &startupParam)
* \param argv list of command line tokens
* \param startupParam startup parameter structure
*/
bool musrFT_filter_histo(int &i, int argc, char *argv[], musrFT_startup_param &startupParam)
Bool_t musrFT_filter_histo(Int_t &i, Int_t argc, Char_t *argv[], musrFT_startup_param &startupParam)
{
int start = i+1, end = 0;
Int_t start = i+1, end = 0;
// find last element of histo option
while (++i < argc) {
@ -217,7 +220,7 @@ bool musrFT_filter_histo(int &i, int argc, char *argv[], musrFT_startup_param &s
// handle histo arguments
TString tstr("");
for (int j=start; j<end; j++) {
for (Int_t j=start; j<end; j++) {
tstr = argv[j];
if (!tstr.Contains("-")) { // likely to be a single number
if (tstr.IsDigit()) {
@ -236,7 +239,7 @@ bool musrFT_filter_histo(int &i, int argc, char *argv[], musrFT_startup_param &s
}
TObjString *ostr;
TString sstr("");
int first=0, last=0;
Int_t first=0, last=0;
ostr = dynamic_cast<TObjString*>(tok->At(0));
sstr = ostr->GetString();
if (sstr.IsDigit()) {
@ -264,7 +267,7 @@ bool musrFT_filter_histo(int &i, int argc, char *argv[], musrFT_startup_param &s
return false;
}
for (int k=first; k<=last; k++) {
for (Int_t k=first; k<=last; k++) {
startupParam.histo.push_back(k);
}
@ -287,11 +290,11 @@ bool musrFT_filter_histo(int &i, int argc, char *argv[], musrFT_startup_param &s
* \param argv command line argument array
* \param startupParam command line data structure
*/
int musrFT_parse_options(int argc, char *argv[], musrFT_startup_param &startupParam)
Int_t musrFT_parse_options(Int_t argc, Char_t *argv[], musrFT_startup_param &startupParam)
{
TString tstr("");
for (int i=1; i<argc; i++) {
for (Int_t i=1; i<argc; i++) {
tstr = argv[i];
if (tstr.BeginsWith("--version")) {
#ifdef HAVE_CONFIG_H
@ -452,8 +455,10 @@ int musrFT_parse_options(int argc, char *argv[], musrFT_startup_param &startupPa
startupParam.timeRange[0] = timeRange[0].Atof();
startupParam.timeRange[1] = timeRange[1].Atof();
i += 2;
} else if (tstr.BeginsWith("-a") || tstr.BeginsWith("--average")) {
} else if (!tstr.CompareTo("-a") || !tstr.CompareTo("--average")) {
startupParam.showAverage = true;
} else if (!tstr.CompareTo("-ad") || !tstr.CompareTo("--average-per-data-set")) {
startupParam.showAveragePerDataSet = true;
} else if (tstr.BeginsWith("--histo")) {
if (!musrFT_filter_histo(i, argc, argv, startupParam))
return 2;
@ -571,16 +576,20 @@ int musrFT_parse_options(int argc, char *argv[], musrFT_startup_param &startupPa
}
if (startupParam.bkg_range[0] > startupParam.bkg_range[1]) {
cerr << endl << ">> musrFT **WARNING** in --background-range, start=" << startupParam.bkg_range[0] << " > end=" << startupParam.bkg_range[1] << ", will swap them." << endl;
double swap = startupParam.bkg_range[0];
Double_t swap = startupParam.bkg_range[0];
startupParam.bkg_range[0] = startupParam.bkg_range[1];
startupParam.bkg_range[1] = swap;
}
if (startupParam.fourierRange[0] > startupParam.fourierRange[1]) {
cerr << endl << ">> musrFT **WARNING** in --fourier-range, start=" << startupParam.fourierRange[0] << " > end=" << startupParam.fourierRange[1] << ", will swap them." << endl;
double swap = startupParam.fourierRange[0];
Double_t swap = startupParam.fourierRange[0];
startupParam.fourierRange[0] = startupParam.fourierRange[1];
startupParam.fourierRange[1] = swap;
}
if (startupParam.showAverage && startupParam.showAveragePerDataSet) {
cerr << endl << ">> musrFT **WARNING** Options: --average and --average-per-data-set exclude each other, will choose the latter." << endl;
startupParam.showAverage = false;
}
return 0;
}
@ -595,7 +604,7 @@ int musrFT_parse_options(int argc, char *argv[], musrFT_startup_param &startupPa
*/
void musrFT_getMetaInfo(const TString fln, PRawRunData *rawRunData, TString &metaInfo)
{
double dval;
Double_t dval;
TString str = fln;
// file name
// trunc it in case a path-name is given
@ -644,12 +653,12 @@ void musrFT_estimateT0(musrFT_data &rd)
cout << endl << " '" << rd.info << "'";
cout << endl << " NO warranty this is sensible!" << endl;
unsigned int idx = 0;
double max = rd.rawData[0];
for (unsigned int i=1; i<rd.rawData.size(); i++) {
UInt_t idx = 0;
Double_t max = rd.rawData[0];
for (UInt_t i=1; i<rd.rawData.size(); i++) {
if (rd.rawData[i] > max) {
max = rd.rawData[i];
idx = (int)i;
idx = (Int_t)i;
}
}
cout << endl << ">> musrFT_estimateT0: estimated t0=" << idx << endl;
@ -679,7 +688,7 @@ void musrFT_cleanup(TH1F *h)
* \param start starting point from where the data shall be written to file.
* \param end ending point up to where the data shall be written to file.
*/
int musrFT_dumpData(TString fln, vector<PFourier*> &fourierData, double start, double end)
Int_t musrFT_dumpData(TString fln, vector<PFourier*> &fourierData, Double_t start, Double_t end)
{
vector<PDoubleVector> data;
PDoubleVector freq;
@ -697,19 +706,19 @@ int musrFT_dumpData(TString fln, vector<PFourier*> &fourierData, double start, d
end = hRe->GetBinCenter(hRe->GetNbinsX());
}
unsigned int minSize = hRe->GetNbinsX()-1;
UInt_t minSize = hRe->GetNbinsX()-1;
musrFT_cleanup(hRe);
for (unsigned int i=1; i<fourierData.size(); i++) {
for (UInt_t i=1; i<fourierData.size(); i++) {
hRe = fourierData[i]->GetRealFourier();
if (hRe->GetNbinsX()-1 < minSize)
minSize = hRe->GetNbinsX()-1;
musrFT_cleanup(hRe);
}
for (unsigned int i=0; i<fourierData.size(); i++) {
for (UInt_t i=0; i<fourierData.size(); i++) {
hRe = fourierData[i]->GetRealFourier();
hIm = fourierData[i]->GetImaginaryFourier();
for (int j=1; j<minSize; j++) {
for (Int_t j=1; j<minSize; j++) {
dval = hRe->GetBinCenter(j);
if ((dval >= start) && (dval <= end)) {
freq.push_back(dval);
@ -735,13 +744,13 @@ int musrFT_dumpData(TString fln, vector<PFourier*> &fourierData, double start, d
// write header
fout << "% ";
for (unsigned int i=0; i<fourierData.size()-1; i++)
for (UInt_t i=0; i<fourierData.size()-1; i++)
fout << "freq" << i << ", Re[d" << i << "], Im[d" << i << "], Pwr[d" << i << "], ";
fout << "freq" << fourierData.size()-1 << ", Re[d" << fourierData.size()-1 << "], Im[d" << fourierData.size()-1 << "], Pwr[d" << fourierData.size()-1 << "]" << endl;
// write data
for (unsigned int j=0; j<data[0].size(); j++) {
for (unsigned int i=0; i<data.size()-1; i++) {
for (UInt_t j=0; j<data[0].size(); j++) {
for (UInt_t i=0; i<data.size()-1; i++) {
fout << data[i][j] << ", ";
}
fout << data[data.size()-1][j] << endl;
@ -761,7 +770,7 @@ int musrFT_dumpData(TString fln, vector<PFourier*> &fourierData, double start, d
* \param run reference to the relevant RUN block of the msr-file
* \param rd data collection which will hold the grouped histograms.
*/
int musrFT_groupHistos(PRunDataHandler *runDataHandler, PMsrGlobalBlock *global, PMsrRunBlock &run, musrFT_data &rd)
Int_t musrFT_groupHistos(PRunDataHandler *runDataHandler, PMsrGlobalBlock *global, PMsrRunBlock &run, musrFT_data &rd)
{
// get proper raw run data set
TString runName = *(run.GetRunName());
@ -773,7 +782,7 @@ int musrFT_groupHistos(PRunDataHandler *runDataHandler, PMsrGlobalBlock *global,
// keep histo list
PIntVector histoList;
for (unsigned int i=0; i<run.GetForwardHistoNoSize(); i++) {
for (UInt_t i=0; i<run.GetForwardHistoNoSize(); i++) {
histoList.push_back(run.GetForwardHistoNo(i));
}
@ -781,10 +790,10 @@ int musrFT_groupHistos(PRunDataHandler *runDataHandler, PMsrGlobalBlock *global,
PDoubleVector t0;
t0.resize(histoList.size());
// init t0 vector
for (unsigned int i=0; i<t0.size(); i++)
for (UInt_t i=0; i<t0.size(); i++)
t0[i] = -1.0;
// 1st: check in the global block
for (unsigned int i=0; i<global->GetT0BinSize(); i++) {
for (UInt_t i=0; i<global->GetT0BinSize(); i++) {
if (i >= t0.size()) { // something is VERY strange
cerr << endl << ">> musrFT_groupHistos **WARNING** found #t0's in GLOBAL block > #histos!";
cerr << endl << ">> This should NEVER happen. Will ignore these entries.";
@ -794,7 +803,7 @@ int musrFT_groupHistos(PRunDataHandler *runDataHandler, PMsrGlobalBlock *global,
}
}
// 2nd: check in the run block
for (unsigned int i=0; i<run.GetT0BinSize(); i++) {
for (UInt_t i=0; i<run.GetT0BinSize(); i++) {
if (i >= t0.size()) { // something is VERY strange
cerr << endl << ">> musrFT_groupHistos **WARNING** found #t0's in RUN block > #histos!";
cerr << endl << ">> This should NEVER happen. Will ignore these entries.";
@ -804,15 +813,15 @@ int musrFT_groupHistos(PRunDataHandler *runDataHandler, PMsrGlobalBlock *global,
}
}
// if still some t0's are == -1, estimate t0
unsigned int idx;
double max;
for (unsigned int i=0; i<t0.size(); i++) {
UInt_t idx;
Double_t max;
for (UInt_t i=0; i<t0.size(); i++) {
if (t0[i] == -1.0) {
cout << endl << ">> musrFT_groupHistos **WARNING** try to estimate t0 from maximum in the data set";
cout << endl << ">> '" << runName << "', histo " << histoList[i] << ". NO warranty this is sensible!";
idx = 0;
max = rawRunData->GetDataBin(histoList[i])->at(0);
for (unsigned int j=1; j<rawRunData->GetDataBin(histoList[i])->size(); j++) {
for (UInt_t j=1; j<rawRunData->GetDataBin(histoList[i])->size(); j++) {
if (rawRunData->GetDataBin(histoList[i])->at(j) > max) {
max = rawRunData->GetDataBin(histoList[i])->at(j);
idx = j;
@ -825,8 +834,8 @@ int musrFT_groupHistos(PRunDataHandler *runDataHandler, PMsrGlobalBlock *global,
// group histos
PDoubleVector data = *(rawRunData->GetDataBin(histoList[0]));
for (unsigned int i=1; i<histoList.size(); i++) {
for (unsigned int j=0; j<data.size(); j++) {
for (UInt_t i=1; i<histoList.size(); i++) {
for (UInt_t j=0; j<data.size(); j++) {
if ((j+t0[i]-t0[0] >= 0) && (j+t0[i]-t0[0] < rawRunData->GetDataBin(histoList[i])->size())) {
data[j] += rawRunData->GetDataBin(histoList[i])->at(j);
}
@ -835,7 +844,7 @@ int musrFT_groupHistos(PRunDataHandler *runDataHandler, PMsrGlobalBlock *global,
rd.rawData.clear();
rd.rawData = data;
rd.t0 = (int)t0[0];
rd.t0 = (Int_t)t0[0];
return 0;
}
@ -881,14 +890,14 @@ void musrFT_dumpMsrFile(musrFT_startup_param &param)
// write RUN block
// get extension of the data file
TString fileFormat("MUSR-ROOT");
for (unsigned int i=0; i<param.dataFln.size(); i++) {
for (UInt_t i=0; i<param.dataFln.size(); i++) {
if (param.dataFileFormat[i].BeginsWith("PsiBin"))
fileFormat = TString("PSI-MDU");
else if (param.dataFileFormat[i].BeginsWith("NeXus"))
fileFormat = TString("NEXUS");
else if (param.dataFileFormat[i].BeginsWith("Mud"))
fileFormat = TString("MUD");
for (unsigned int j=0; j<param.histo.size(); j++) {
for (UInt_t j=0; j<param.histo.size(); j++) {
fout << "RUN " << param.dataFln[i] << " BXXX IXX " << fileFormat << " (name beamline institute data-file-format)" << endl;
fout << "forward " << param.histo[j] << endl;
if ((param.t0.size() > 1) && (j < param.t0.size())) {
@ -908,7 +917,7 @@ void musrFT_dumpMsrFile(musrFT_startup_param &param)
fout << "runs 1" << endl;
} else {
fout << "runs ";
for (unsigned int i=0; i<param.histo.size(); i++)
for (UInt_t i=0; i<param.histo.size(); i++)
fout << i+1 << " ";
fout << endl;
}
@ -949,12 +958,12 @@ void musrFT_dumpMsrFile(musrFT_startup_param &param)
*
* <b>return:</b> time stamp with msec resolution.
*/
double millitime()
Double_t millitime()
{
struct timeval now;
gettimeofday(&now, 0);
return ((double)now.tv_sec * 1.0e6 + (double)now.tv_usec)/1.0e3;
return ((Double_t)now.tv_sec * 1.0e6 + (Double_t)now.tv_usec)/1.0e3;
}
//-------------------------------------------------------------------------
@ -968,7 +977,7 @@ double millitime()
* \param argc number of command line arguments
* \param argv command line argument array
*/
int main(int argc, char *argv[])
Int_t main(Int_t argc, Char_t *argv[])
{
Int_t unitTag = FOURIER_UNIT_NOT_GIVEN;
Int_t apodTag = F_APODIZATION_NONE;
@ -985,9 +994,9 @@ int main(int argc, char *argv[])
musrFT_init(startupParam);
// parse command line options
int status = musrFT_parse_options(argc, argv, startupParam);
Int_t status = musrFT_parse_options(argc, argv, startupParam);
if (status != 0) {
int retVal = PMUSR_SUCCESS;
Int_t retVal = PMUSR_SUCCESS;
if (status == 2) {
musrFT_syntax();
retVal = PMUSR_WRONG_STARTUP_SYNTAX;
@ -1002,7 +1011,7 @@ int main(int argc, char *argv[])
}
// read startup file
char startup_path_name[128];
Char_t startup_path_name[128];
PStartupOptions startup_options;
startup_options.writeExpectedChisq = false;
startup_options.estimateN0 = true;
@ -1051,7 +1060,7 @@ int main(int argc, char *argv[])
// load msr-file(s)
vector<PMsrHandler*> msrHandler;
msrHandler.resize(startupParam.msrFln.size());
for (unsigned int i=0; i<startupParam.msrFln.size(); i++) {
for (UInt_t i=0; i<startupParam.msrFln.size(); i++) {
msrHandler[i] = new PMsrHandler(startupParam.msrFln[i].Data(), startupHandler->GetStartupOptions(), true);
status = msrHandler[i]->ReadMsrFile();
if (status != PMUSR_SUCCESS) {
@ -1073,7 +1082,7 @@ int main(int argc, char *argv[])
vector<PRunDataHandler*> runDataHandler;
runDataHandler.resize(startupParam.msrFln.size()+startupParam.dataFln.size()); // resize to the total number of run data provided
// load data-file(s) related to msr-file
for (unsigned int i=0; i<msrHandler.size(); i++) {
for (UInt_t i=0; i<msrHandler.size(); i++) {
// create run data handler
if (startupHandler)
runDataHandler[i] = new PRunDataHandler(msrHandler[i], startupHandler->GetDataPathList());
@ -1082,7 +1091,7 @@ int main(int argc, char *argv[])
}
// load data-file(s) provided directly
for (unsigned int i=msrHandler.size(); i<msrHandler.size()+startupParam.dataFln.size(); i++) {
for (UInt_t i=msrHandler.size(); i<msrHandler.size()+startupParam.dataFln.size(); i++) {
// create run data handler
if (startupHandler)
runDataHandler[i] = new PRunDataHandler(startupParam.dataFln[i-msrHandler.size()], startupParam.dataFileFormat[i-msrHandler.size()], startupHandler->GetDataPathList());
@ -1091,7 +1100,14 @@ int main(int argc, char *argv[])
}
// read all the data files
for (unsigned int i=0; i<runDataHandler.size(); i++) {
musrFT_data rd;
rd.dataSetTag = -1;
Int_t dataSetTagCounter = 0;
TString prevDataSetPathName("");
TString str(""), fln("");
UInt_t idx=0;
for (UInt_t i=0; i<runDataHandler.size(); i++) {
runDataHandler[i]->ReadData();
if (!runDataHandler[i]->IsAllDataAvailable()) {
@ -1116,8 +1132,8 @@ int main(int argc, char *argv[])
// first check of histo list makes sense
if (i >= msrHandler.size()) { // only check if originating from data-files (not msr-files)
for (unsigned int j=0; j<startupParam.histo.size(); j++) {
if ((unsigned int)startupParam.histo[j] > rawRunData->GetNoOfHistos()) {
for (UInt_t j=0; j<startupParam.histo.size(); j++) {
if ((UInt_t)startupParam.histo[j] > rawRunData->GetNoOfHistos()) {
cerr << endl << ">> musrFT **ERROR** found histo no " << startupParam.histo[j] << " > # of histo in the file (";
cerr << startupParam.dataFln[i] << " // # histo: " << rawRunData->GetNoOfHistos() << ")." << endl;
return PMUSR_DATA_FILE_READ_ERROR;
@ -1125,15 +1141,18 @@ int main(int argc, char *argv[])
}
if (startupParam.histo.size() == 0) { // no histo list given
// set histo list to ALL available histos for the data file
for (unsigned int j=0; j<rawRunData->GetNoOfHistos(); j++)
for (UInt_t j=0; j<rawRunData->GetNoOfHistos(); j++)
startupParam.histo.push_back(j+1);
}
}
musrFT_data rd;
TString str(""), fln("");
unsigned int idx=0;
// get meta info, time resolution, time range, raw data sets
// check if the data set path-name has changed
if (prevDataSetPathName.CompareTo(runDataHandler[i]->GetRunPathName())) { // i.e. data set path-name changed
rd.dataSetTag = dataSetTagCounter++;
prevDataSetPathName = runDataHandler[i]->GetRunPathName();
}
if (i < msrHandler.size()) { // obtain info from msr-files
// keep title if not overwritten by the command line
if (startupParam.title.Length() == 0)
@ -1206,11 +1225,11 @@ int main(int argc, char *argv[])
PIntVector runList = plot->at(0).fRuns;
// loop over all runs listed in the msr-file PLOT block
for (unsigned int j=0; j<runList.size(); j++) {
for (UInt_t j=0; j<runList.size(); j++) {
// keep forward histo list
PIntVector histoList;
for (unsigned int k=0; k<runs->at(runList[j]-1).GetForwardHistoNoSize(); k++) {
for (UInt_t k=0; k<runs->at(runList[j]-1).GetForwardHistoNoSize(); k++) {
histoList.push_back(runs->at(runList[j]-1).GetForwardHistoNo(k));
}
@ -1219,7 +1238,7 @@ int main(int argc, char *argv[])
musrFT_getMetaInfo(fln, rawRunData, str);
TString hh("");
hh = TString::Format("h%d", histoList[0]);
for (unsigned int k=1; k<histoList.size(); k++)
for (UInt_t k=1; k<histoList.size(); k++)
hh += TString::Format("/%d", histoList[k]);
hh += ":";
rd.info = hh;
@ -1262,7 +1281,7 @@ int main(int argc, char *argv[])
}
} else { // obtain info from command line options for direct data-file read
musrFT_getMetaInfo(startupParam.dataFln[i-msrHandler.size()], rawRunData, str);
for (unsigned int j=0; j<startupParam.histo.size(); j++) {
for (UInt_t j=0; j<startupParam.histo.size(); j++) {
idx = startupParam.histo[j];
// handle meta information
@ -1295,6 +1314,12 @@ int main(int argc, char *argv[])
}
}
// generate data set label vector
PIntVector dataSetTag;
for (UInt_t i=0; i<data.GetNoOfData(); i++) {
dataSetTag.push_back(data.GetDataSetTag(i));
}
// make sure Fourier plot tag is set
if (fourierPlotTag == FOURIER_PLOT_NOT_GIVEN) {
if (!startupParam.fourierOpt.CompareTo("real", TString::kIgnoreCase))
@ -1338,7 +1363,7 @@ int main(int argc, char *argv[])
vector<PFourier*> fourier;
fourier.resize(histo.size());
for (unsigned int i=0; i<fourier.size(); i++) {
for (UInt_t i=0; i<fourier.size(); i++) {
fourier[i] = new PFourier(histo[i], unitTag, 0.0, 0.0, true, startupParam.fourierPower);
}
@ -1350,11 +1375,11 @@ int main(int argc, char *argv[])
else if (startupParam.apodization.BeginsWith("strong", TString::kIgnoreCase))
apodTag = F_APODIZATION_STRONG;
double start = millitime();
for (unsigned int i=0; i<fourier.size(); i++) {
Double_t start = millitime();
for (UInt_t i=0; i<fourier.size(); i++) {
fourier[i]->Transform(apodTag);
}
double end = millitime();
Double_t end = millitime();
cout << endl << "info> after FFT. calculation time: " << (end-start)/1.0e3 << " (sec)." << endl;
// make sure that a Fourier range is provided, if not calculate one
@ -1374,7 +1399,7 @@ int main(int argc, char *argv[])
Bool_t batch = false;
if (startupParam.graphicFormat.Length() != 0) {
batch = true;
argv[argc] = (char*)malloc(16*sizeof(char));
argv[argc] = (Char_t*)malloc(16*sizeof(Char_t));
strcpy(argv[argc], "-b");
argc++;
}
@ -1383,17 +1408,17 @@ int main(int argc, char *argv[])
TApplication app("App", &argc, argv);
if (startupHandler) {
fourierCanvas = new PFourierCanvas(fourier, startupParam.title.Data(),
startupParam.showAverage, fourierPlotTag,
startupParam.fourierRange, startupParam.initialPhase,
fourierCanvas = new PFourierCanvas(fourier, dataSetTag, startupParam.title.Data(),
startupParam.showAverage, startupParam.showAveragePerDataSet,
fourierPlotTag, startupParam.fourierRange, startupParam.initialPhase,
10, 10, 800, 800,
startupHandler->GetMarkerList(),
startupHandler->GetColorList(),
batch);
} else {
fourierCanvas = new PFourierCanvas(fourier, startupParam.title.Data(),
startupParam.showAverage, fourierPlotTag,
startupParam.fourierRange, startupParam.initialPhase,
fourierCanvas = new PFourierCanvas(fourier, dataSetTag, startupParam.title.Data(),
startupParam.showAverage, startupParam.showAveragePerDataSet,
fourierPlotTag, startupParam.fourierRange, startupParam.initialPhase,
10, 10, 800, 800,
batch);
}
@ -1443,23 +1468,23 @@ int main(int argc, char *argv[])
if (startupHandler)
delete startupHandler;
for (unsigned int i=0; i<msrHandler.size(); i++)
for (UInt_t i=0; i<msrHandler.size(); i++)
if (msrHandler[i])
delete msrHandler[i];
msrHandler.clear();
for (unsigned int i=0; i<runDataHandler.size(); i++)
for (UInt_t i=0; i<runDataHandler.size(); i++)
if (runDataHandler[i])
delete runDataHandler[i];
runDataHandler.clear();
if (histo.size() > 0) {
for (unsigned int i=0; i<histo.size(); i++)
for (UInt_t i=0; i<histo.size(); i++)
delete histo[i];
histo.clear();
}
if (fourier.size() > 0) {
for (unsigned int i=0; i<fourier.size(); i++)
for (UInt_t i=0; i<fourier.size(); i++)
delete fourier[i];
fourier.clear();
}

34
src/musredit/PGetMusrFTOptionsDialog.cpp
View File

@ -197,10 +197,14 @@ PGetMusrFTOptionsDialog::PGetMusrFTOptionsDialog(QString currentMsrFile, QString
line.remove(line.length()-1, 1);
fHistoList_lineEdit->setText(line);
}
// average tag
// average ALL tag
if (prevCmd[i] == "-a") {
fAveragedView_checkBox->setCheckState(Qt::Checked);
}
// average per data set tag
if (prevCmd[i] == "-ad") {
fAveragePerDataSet_checkBox->setCheckState(Qt::Checked);
}
// t0 list
if (prevCmd[i] == "--t0") {
i++;
@ -239,6 +243,8 @@ PGetMusrFTOptionsDialog::PGetMusrFTOptionsDialog(QString currentMsrFile, QString
connect(fMsrFileNameClear_pushButton, SIGNAL (clicked() ), this, SLOT( clearMsrFileNames() ) );
connect(fDataFileNameClear_pushButton, SIGNAL (clicked() ), this, SLOT( clearDataFileNames() ) );
connect(fResetAll_pushButton, SIGNAL( clicked() ), this, SLOT( resetAll() ) );
connect(fAveragedView_checkBox, SIGNAL ( stateChanged(int) ), this, SLOT( averagedAll(int) ) );
connect(fAveragePerDataSet_checkBox, SIGNAL ( stateChanged(int) ), this, SLOT( averagedPerDataSet(int) ) );
}
//----------------------------------------------------------------------------------------------------
@ -329,10 +335,14 @@ QStringList PGetMusrFTOptionsDialog::getMusrFTOptions()
cmd << strList[i];
}
// averaged view
// averaged view ALL
if (fAveragedView_checkBox->checkState() == Qt::Checked)
cmd << "-a";
// averaged view per data set
if (fAveragePerDataSet_checkBox->checkState() == Qt::Checked)
cmd << "-ad";
// t0 list
if (fT0_lineEdit->text().length() > 0) {
cmd << "--t0";
@ -546,6 +556,26 @@ void PGetMusrFTOptionsDialog::resetAll()
fFourierTitle_lineEdit->setText("");
}
//----------------------------------------------------------------------------------------------------
/**
* <p>SLOT called when averaged view for ALL data is checked.
*/
void PGetMusrFTOptionsDialog::averagedAll(int state)
{
if ((state == Qt::Checked) && fAveragePerDataSet_checkBox->isChecked())
fAveragePerDataSet_checkBox->setCheckState(Qt::Unchecked);
}
//----------------------------------------------------------------------------------------------------
/**
* <p>SLOT called when averaged view per data set is checked.
*/
void PGetMusrFTOptionsDialog::averagedPerDataSet(int state)
{
if ((state == Qt::Checked) && fAveragedView_checkBox->isChecked())
fAveragedView_checkBox->setCheckState(Qt::Unchecked);
}
//----------------------------------------------------------------------------------------------------
/**
* <p>Generates a help content window showing the description for musrFT.

2
src/musredit/PGetMusrFTOptionsDialog.h
View File

@ -55,6 +55,8 @@ class PGetMusrFTOptionsDialog : public QDialog, private Ui::PGetMusrFTOptionsDia
void clearDataFileNames();
void createMsrFileChanged(int state);
void resetAll();
void averagedAll(int state);
void averagedPerDataSet(int state);
private:
QStringList fMsrFilePaths; ///< list keeping all the paths from the msr-file path-name list

37
src/musredit/forms/PGetMusrFTOptionsDialog.ui
View File

@ -34,11 +34,11 @@
<property name="title">
<string> Fourier </string>
</property>
<widget class="QWidget" name="layoutWidget">
<widget class="QWidget" name="">
<property name="geometry">
<rect>
<x>20</x>
<y>30</y>
<x>21</x>
<y>31</y>
<width>641</width>
<height>156</height>
</rect>
@ -264,7 +264,14 @@ p, li { white-space: pre-wrap; }
<item>
<widget class="QCheckBox" name="fAveragedView_checkBox">
<property name="text">
<string>Averaged View</string>
<string>Average All</string>
</property>
</widget>
</item>
<item>
<widget class="QCheckBox" name="fAveragePerDataSet_checkBox">
<property name="text">
<string>Average per Data Set</string>
</property>
</widget>
</item>
@ -275,7 +282,7 @@ p, li { white-space: pre-wrap; }
</property>
<property name="sizeHint" stdset="0">
<size>
<width>278</width>
<width>88</width>
<height>20</height>
</size>
</property>
@ -644,7 +651,14 @@ p, li { white-space: pre-wrap; }
</widget>
</item>
<item>
<widget class="QLineEdit" name="fMsrFileSelector_lineEdit"/>
<widget class="QLineEdit" name="fMsrFileSelector_lineEdit">
<property name="minimumSize">
<size>
<width>0</width>
<height>23</height>
</size>
</property>
</widget>
</item>
<item>
<widget class="QPushButton" name="fMsrFileNameClear_pushButton">
@ -668,7 +682,14 @@ p, li { white-space: pre-wrap; }
</widget>
</item>
<item>
<widget class="QLineEdit" name="fDataFileSelector_lineEdit"/>
<widget class="QLineEdit" name="fDataFileSelector_lineEdit">
<property name="minimumSize">
<size>
<width>0</width>
<height>23</height>
</size>
</property>
</widget>
</item>
<item>
<widget class="QPushButton" name="fDataFileNameClear_pushButton">
@ -685,7 +706,7 @@ p, li { white-space: pre-wrap; }
</layout>
</widget>
</widget>
<widget class="QWidget" name="">
<widget class="QWidget" name="layoutWidget">
<property name="geometry">
<rect>
<x>20</x>