more work towards a standalone Fourier (mainly for HAL-9500). Add msr-file support. Still some severe flaws.
This commit is contained in:
parent
f56c2b5fda
commit
d4ba2f6d81
@ -269,7 +269,7 @@ Int_t PMsrHandler::ReadMsrFile()
|
||||
}
|
||||
|
||||
// fill parameter-in-use vector
|
||||
if (result == PMUSR_SUCCESS)
|
||||
if ((result == PMUSR_SUCCESS) && !fFourierOnly)
|
||||
FillParameterInUse(theory, functions, run);
|
||||
|
||||
// check that each run fulfills the minimum requirements
|
||||
@ -278,7 +278,7 @@ Int_t PMsrHandler::ReadMsrFile()
|
||||
result = PMUSR_MSR_SYNTAX_ERROR;
|
||||
|
||||
// check that parameter names are unique
|
||||
if (result == PMUSR_SUCCESS) {
|
||||
if ((result == PMUSR_SUCCESS) && !fFourierOnly) {
|
||||
UInt_t parX, parY;
|
||||
if (!CheckUniquenessOfParamNames(parX, parY)) {
|
||||
cerr << endl << ">> PMsrHandler::ReadMsrFile: **SEVERE ERROR** parameter name " << fParam[parX].fName.Data() << " is identical for parameter no " << fParam[parX].fNo << " and " << fParam[parY].fNo << "!";
|
||||
@ -290,7 +290,7 @@ Int_t PMsrHandler::ReadMsrFile()
|
||||
|
||||
// check that if maps are present in the theory- and/or function-block,
|
||||
// that there are really present in the run block
|
||||
if (result == PMUSR_SUCCESS)
|
||||
if ((result == PMUSR_SUCCESS) && !fFourierOnly)
|
||||
if (!CheckMaps())
|
||||
result = PMUSR_MSR_SYNTAX_ERROR;
|
||||
|
||||
@ -2387,6 +2387,10 @@ Int_t PMsrHandler::ParameterInUse(UInt_t paramNo)
|
||||
*/
|
||||
Bool_t PMsrHandler::HandleFitParameterEntry(PMsrLines &lines)
|
||||
{
|
||||
// If msr-file is used for musrFT only, nothing needs to be done here.
|
||||
if (fFourierOnly)
|
||||
return true;
|
||||
|
||||
PMsrParamStructure param;
|
||||
Bool_t error = false;
|
||||
|
||||
@ -2599,6 +2603,10 @@ Bool_t PMsrHandler::HandleFitParameterEntry(PMsrLines &lines)
|
||||
*/
|
||||
Bool_t PMsrHandler::HandleTheoryEntry(PMsrLines &lines)
|
||||
{
|
||||
// If msr-file is used for musrFT only, nothing needs to be done here.
|
||||
if (fFourierOnly)
|
||||
return true;
|
||||
|
||||
// store the theory lines
|
||||
fTheory = lines;
|
||||
|
||||
@ -2619,6 +2627,10 @@ Bool_t PMsrHandler::HandleTheoryEntry(PMsrLines &lines)
|
||||
*/
|
||||
Bool_t PMsrHandler::HandleFunctionsEntry(PMsrLines &lines)
|
||||
{
|
||||
// If msr-file is used for musrFT only, nothing needs to be done here.
|
||||
if (fFourierOnly)
|
||||
return true;
|
||||
|
||||
// store the functions lines
|
||||
fFunctions = lines;
|
||||
|
||||
@ -3157,11 +3169,13 @@ Bool_t PMsrHandler::HandleRunEntry(PMsrLines &lines)
|
||||
}
|
||||
}
|
||||
// check map entries, i.e. if the map values are within parameter bounds
|
||||
for (UInt_t i=0; i<param.GetMap()->size(); i++) {
|
||||
if ((param.GetMap(i) < 0) || (param.GetMap(i) > (Int_t) fParam.size())) {
|
||||
cerr << endl << ">> PMsrHandler::HandleRunEntry: **SEVERE ERROR** map value " << param.GetMap(i) << " in line " << iter->fLineNo << " is out of range!";
|
||||
error = true;
|
||||
break;
|
||||
if (!fFourierOnly) {
|
||||
for (UInt_t i=0; i<param.GetMap()->size(); i++) {
|
||||
if ((param.GetMap(i) < 0) || (param.GetMap(i) > (Int_t) fParam.size())) {
|
||||
cerr << endl << ">> PMsrHandler::HandleRunEntry: **SEVERE ERROR** map value " << param.GetMap(i) << " in line " << iter->fLineNo << " is out of range!";
|
||||
error = true;
|
||||
break;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -3515,9 +3529,13 @@ Bool_t PMsrHandler::FilterNumber(TString str, const Char_t *filter, Int_t offset
|
||||
*/
|
||||
Bool_t PMsrHandler::HandleCommandsEntry(PMsrLines &lines)
|
||||
{
|
||||
// If msr-file is used for musrFT only, nothing needs to be done here.
|
||||
if (fFourierOnly)
|
||||
return true;
|
||||
|
||||
PMsrLines::iterator iter;
|
||||
|
||||
if (lines.empty() && !fFourierOnly) {
|
||||
if (lines.empty()) {
|
||||
cerr << endl << ">> PMsrHandler::HandleCommandsEntry(): **WARNING**: There is no COMMAND block! Do you really want this?";
|
||||
cerr << endl;
|
||||
}
|
||||
@ -4352,7 +4370,11 @@ Bool_t PMsrHandler::HandlePlotEntry(PMsrLines &lines)
|
||||
*/
|
||||
Bool_t PMsrHandler::HandleStatisticEntry(PMsrLines &lines)
|
||||
{
|
||||
if (lines.empty() && !fFourierOnly) {
|
||||
// If msr-file is used for musrFT only, nothing needs to be done here.
|
||||
if (fFourierOnly)
|
||||
return true;
|
||||
|
||||
if (lines.empty()) {
|
||||
cerr << endl << ">> PMsrHandler::HandleStatisticEntry: **WARNING** There is no STATISTIC block! Do you really want this?";
|
||||
cerr << endl;
|
||||
fStatistic.fValid = false;
|
||||
@ -5042,7 +5064,7 @@ Bool_t PMsrHandler::CheckRunBlockIntegrity()
|
||||
}
|
||||
}
|
||||
// check fit range
|
||||
if (!fRuns[i].IsFitRangeInBin()) { // fit range given as times in usec (RUN block)
|
||||
if (!fRuns[i].IsFitRangeInBin() && !fFourierOnly) { // fit range given as times in usec (RUN block)
|
||||
if ((fRuns[i].GetFitRange(0) == PMUSR_UNDEFINED) || (fRuns[i].GetFitRange(1) == PMUSR_UNDEFINED)) { // check fit range in RUN block
|
||||
if (!fGlobal.IsFitRangeInBin()) { // fit range given as times in usec (GLOBAL block)
|
||||
if ((fGlobal.GetFitRange(0) == PMUSR_UNDEFINED) || (fGlobal.GetFitRange(1) == PMUSR_UNDEFINED)) { // check fit range in GLOBAL block
|
||||
|
||||
@ -155,6 +155,11 @@ void PPrepFourier::DoBkgCorrection()
|
||||
{
|
||||
cout << endl << "debug> DoBkgCorrection ...";
|
||||
|
||||
// make sure fData are already present, and if not create the necessary data sets
|
||||
if (fData.size() != fRawData.size()) {
|
||||
InitData();
|
||||
}
|
||||
|
||||
// if no bkg-range is given, nothing needs to be done
|
||||
if ((fBkgRange[0] == -1) && (fBkgRange[1] == -1)) {
|
||||
cout << endl << "debug> no background range given ...";
|
||||
@ -169,11 +174,6 @@ void PPrepFourier::DoBkgCorrection()
|
||||
}
|
||||
}
|
||||
|
||||
// make sure fData are already present, and if not create the necessary data sets
|
||||
if (fData.size() != fRawData.size()) {
|
||||
InitData();
|
||||
}
|
||||
|
||||
Double_t bkg=0.0;
|
||||
for (unsigned int i=0; i<fRawData.size(); i++) {
|
||||
// calculate the bkg for the given range
|
||||
@ -198,7 +198,12 @@ void PPrepFourier::DoBkgCorrection()
|
||||
*/
|
||||
void PPrepFourier::DoPacking()
|
||||
{
|
||||
cout << endl << "debug> DoPacking ...";
|
||||
cout << endl << "debug> DoPacking : pack=" << fPacking << " ...";
|
||||
|
||||
// make sure fData are already present, and if not create the necessary data sets
|
||||
if (fData.size() != fRawData.size()) {
|
||||
InitData();
|
||||
}
|
||||
|
||||
if (fPacking == 1) // nothing to be done
|
||||
return;
|
||||
@ -236,6 +241,53 @@ void PPrepFourier::DoPacking()
|
||||
void PPrepFourier::DoFiltering()
|
||||
{
|
||||
cout << endl << "debug> DoFiltering not yet implemented ...";
|
||||
|
||||
// make sure fData are already present, and if not create the necessary data sets
|
||||
if (fData.size() != fRawData.size()) {
|
||||
InitData();
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------
|
||||
// DoLifeTimeCorrection
|
||||
//--------------------------------------------------------------------------
|
||||
/**
|
||||
* <p>
|
||||
*
|
||||
* \param fudge
|
||||
*/
|
||||
void PPrepFourier::DoLifeTimeCorrection(Double_t fudge)
|
||||
{
|
||||
// make sure fData are already present, and if not create the necessary data sets
|
||||
if (fData.size() != fRawData.size()) {
|
||||
InitData();
|
||||
}
|
||||
|
||||
// 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++) {
|
||||
scale = fRawData[i].timeResolution / PMUON_LIFETIME;
|
||||
for (unsigned int j=0; j<fData[i].size(); j++) {
|
||||
fData[i][j] *= exp(j*scale);
|
||||
}
|
||||
}
|
||||
|
||||
// calc N0
|
||||
Double_t dval;
|
||||
Double_t N0;
|
||||
for (unsigned int i=0; i<fData.size(); i++) {
|
||||
dval = 0.0;
|
||||
for (unsigned int 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++) {
|
||||
fData[i][j] -= N0;
|
||||
fData[i][j] /= N0;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------
|
||||
@ -268,6 +320,10 @@ vector<TH1F*> PPrepFourier::GetData()
|
||||
vector<TH1F*> data;
|
||||
data.resize(fData.size());
|
||||
|
||||
// if not data are present, just return an empty vector
|
||||
if (fData.size() == 0)
|
||||
return data;
|
||||
|
||||
TString name("");
|
||||
Double_t dt=0.0;
|
||||
Double_t start=0.0;
|
||||
@ -330,7 +386,10 @@ vector<TH1F*> PPrepFourier::GetData()
|
||||
*/
|
||||
TH1F *PPrepFourier::GetData(const UInt_t idx)
|
||||
{
|
||||
if (idx > fData.size())
|
||||
if (fData.size() == 0) // no data present
|
||||
return 0;
|
||||
|
||||
if (idx > fData.size()) // requested index out of range
|
||||
return 0;
|
||||
|
||||
TString name = TString::Format("histo%2d", idx);
|
||||
@ -381,8 +440,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++) {
|
||||
for (unsigned int j=fRawData[i].t0; j<fRawData[i].rawData.size(); j++) {
|
||||
if (fRawData[i].t0 >= 0)
|
||||
t0 = fRawData[i].t0;
|
||||
else
|
||||
t0 = 0;
|
||||
for (unsigned int j=t0; j<fRawData[i].rawData.size(); j++) {
|
||||
fData[i].push_back(fRawData[i].rawData[j]);
|
||||
}
|
||||
}
|
||||
|
||||
@ -704,7 +704,7 @@ typedef vector<PMsrRunBlock> PMsrRunList;
|
||||
*/
|
||||
typedef struct {
|
||||
Bool_t fFourierBlockPresent; ///< flag indicating if a Fourier block is present in the msr-file
|
||||
Int_t fUnits; ///< flag used to indicate the units. 0=field units (G); 1=frequency units (MHz); 2=Mc/s
|
||||
Int_t fUnits; ///< flag used to indicate the units. 1=field units (G); 2=field units (T); 3=frequency units (MHz); 4=Mc/s
|
||||
Bool_t fDCCorrected; ///< if set true, the dc offset of the signal/theory will be removed before the FFT is made.
|
||||
Int_t fFourierPower; ///< i.e. zero padding up to 2^fFourierPower, default = 0 which means NO zero padding
|
||||
Int_t fApodization; ///< tag indicating the kind of apodization wished, 0=no appodization (default), 1=weak, 2=medium, 3=strong (for details see the docu)
|
||||
|
||||
@ -67,6 +67,7 @@ class PPrepFourier {
|
||||
void DoBkgCorrection();
|
||||
void DoPacking();
|
||||
void DoFiltering();
|
||||
void DoLifeTimeCorrection(Double_t fudge);
|
||||
|
||||
TString GetInfo(const UInt_t idx);
|
||||
UInt_t GetNoOfData() { return fRawData.size(); }
|
||||
|
||||
364
src/musrFT.cpp
364
src/musrFT.cpp
@ -74,6 +74,7 @@ typedef struct {
|
||||
vector<int> t0;
|
||||
int packing;
|
||||
TString title;
|
||||
double lifetimecorrection; //< is == 0.0 for NO life time correction, otherwise it holds the fudge factor
|
||||
} musrFT_startup_param;
|
||||
|
||||
//-------------------------------------------------------------------------
|
||||
@ -100,7 +101,7 @@ void musrFT_syntax()
|
||||
cout << endl << " Supported graphic-format-extension: eps, pdf, gif, jpg, png, svg, xpm, root";
|
||||
cout << endl << " --dump <fln> : rather than starting a root session and showing Fourier graphs of the data,";
|
||||
cout << endl << " it will output the Fourier data in an ascii file <fln>.";
|
||||
cout << endl << " --filter : filter and filter-specific-information -- TO BE WRITTEN YET.";
|
||||
cout << endl << " --filter : filter and filter-specific-information -- ***TO BE WRITTEN YET***.";
|
||||
cout << endl << " -b, --background <start> <end>: background interval used to estimate the backround to be";
|
||||
cout << endl << " subtracted before the Fourier transform. <start>, <end> to be given in bins.";
|
||||
cout << endl << " -fo, --fourier-option <fopt>: <fopt> can be 'real', 'power', 'imag', 'real+imag', of 'phase'.";
|
||||
@ -128,12 +129,17 @@ void musrFT_syntax()
|
||||
cout << endl << " to all data-files given. If --histo is not given, all histos of a data file will be used.";
|
||||
cout << endl << " -a, --average : show the average of all Fourier transformed data.";
|
||||
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 inital 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";
|
||||
cout << endl << " to all histos.";
|
||||
cout << endl << " Example: musrFT -df lem15_his_01234.root -fo real --t0 2750 --histo 1 3";
|
||||
cout << endl << " --fourier-option real allows to get the proper inital 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";
|
||||
cout << endl << " to all histos.";
|
||||
cout << endl << " Example: musrFT -df lem15_his_01234.root -fo real --t0 2750 --histo 1 3";
|
||||
cout << endl << " -pa, --packing <N> : if <N> (an integer), the time domain data will first be packed/rebinned by <N>.";
|
||||
cout << endl << " -t, --title <title> : give a global title for the plot.";
|
||||
cout << endl << " --create-msr-file <fln> : creates a msr-file based on the command line options";
|
||||
cout << endl << " provided. This will help on the way to a full fitting model.";
|
||||
cout << endl << " ***TO BE WRITTEN YET.***";
|
||||
cout << endl << " -lc, --lifetimecorrection <fudge>: try to eliminate muon life time decay. Only makes sense for low";
|
||||
cout << endl << " transverse fields. <fudge> is a tweaking factor and should be kept around 1.0.";
|
||||
cout << endl << endl;
|
||||
}
|
||||
|
||||
@ -152,7 +158,7 @@ void musrFT_init(musrFT_startup_param &startupParam)
|
||||
startupParam.fourierOpt = TString("real");
|
||||
startupParam.apotization = TString("none");
|
||||
startupParam.fourierPower = -1;
|
||||
startupParam.fourierUnits = TString("MHz");
|
||||
startupParam.fourierUnits = TString("??");
|
||||
startupParam.initialPhase = 0.0;
|
||||
startupParam.fourierRange[0] = -1.0;
|
||||
startupParam.fourierRange[1] = -1.0;
|
||||
@ -161,6 +167,7 @@ void musrFT_init(musrFT_startup_param &startupParam)
|
||||
startupParam.showAverage = false;
|
||||
startupParam.packing = 1;
|
||||
startupParam.title = TString("");
|
||||
startupParam.lifetimecorrection = 0.0;
|
||||
}
|
||||
|
||||
//-------------------------------------------------------------------------
|
||||
@ -461,6 +468,18 @@ int musrFT_parse_options(int argc, char *argv[], musrFT_startup_param &startupPa
|
||||
return 2;
|
||||
}
|
||||
startupParam.packing = pack.Atoi();
|
||||
} else if (tstr.BeginsWith("-lc") || tstr.BeginsWith("--lifetimecorrection")) {
|
||||
if (i+1 >= argc) { // something is wrong since there needs to be an argument here
|
||||
cerr << endl << ">> musrFT **ERROR** found option --lifetimecorrection without argument!" << endl;
|
||||
return 2;
|
||||
}
|
||||
++i;
|
||||
TString fudge(argv[i]);
|
||||
if (!fudge.IsFloat()) {
|
||||
cerr << endl << ">> musrFT **ERROR** found option --lifetimecorrection with a fudge which is not an double '" << fudge << "'." << endl;
|
||||
return 2;
|
||||
}
|
||||
startupParam.lifetimecorrection = fudge.Atof();
|
||||
} else if (tstr.BeginsWith("-df") || tstr.BeginsWith("--data-file")) {
|
||||
while (++i < argc) {
|
||||
if (argv[i][0] == '-') {
|
||||
@ -520,6 +539,10 @@ int musrFT_parse_options(int argc, char *argv[], musrFT_startup_param &startupPa
|
||||
return 0;
|
||||
}
|
||||
|
||||
//----------------------------------------------------------------------------------------
|
||||
/**
|
||||
*
|
||||
*/
|
||||
void musrFT_getMetaInfo(const TString fln, PRawRunData *rawRunData, TString &metaInfo)
|
||||
{
|
||||
double dval;
|
||||
@ -572,7 +595,7 @@ void musrFT_estimateT0(musrFT_data &rd)
|
||||
idx = (int)i;
|
||||
}
|
||||
}
|
||||
cout << endl << "debug> estimated t0=" << idx << endl;
|
||||
cout << endl << ">> musrFT_estimateT0: estimated t0=" << idx << endl;
|
||||
rd.t0 = idx;
|
||||
}
|
||||
|
||||
@ -668,6 +691,94 @@ int musrFT_dumpData(TString fln, vector<PFourier*> &fourierData, double start, d
|
||||
return 0;
|
||||
}
|
||||
|
||||
//---------------------------------------------------
|
||||
/**
|
||||
* <p>
|
||||
*
|
||||
* \param runDataHandler
|
||||
* \param global
|
||||
* \param run
|
||||
* \param rd
|
||||
*/
|
||||
int musrFT_groupHistos(PRunDataHandler *runDataHandler, PMsrGlobalBlock *global, PMsrRunBlock &run, musrFT_data &rd)
|
||||
{
|
||||
// get proper raw run data set
|
||||
TString runName = *(run.GetRunName());
|
||||
PRawRunData *rawRunData = runDataHandler->GetRunData(runName);
|
||||
if (rawRunData == 0) {
|
||||
cerr << endl << ">> musrFT_groupHistos **ERROR** Couldn't get raw run data for run '" << runName << "'." << endl;
|
||||
return 1;
|
||||
}
|
||||
|
||||
// keep histo list
|
||||
PIntVector histoList;
|
||||
for (unsigned int i=0; i<run.GetForwardHistoNoSize(); i++) {
|
||||
histoList.push_back(run.GetForwardHistoNo(i));
|
||||
}
|
||||
|
||||
// check if t0's are found and that #t0 == #histos
|
||||
PDoubleVector t0;
|
||||
t0.resize(histoList.size());
|
||||
// init t0 vector
|
||||
for (unsigned int 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++) {
|
||||
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.";
|
||||
cerr << endl << ">> Please check your msr-file!!" << endl;
|
||||
} else {
|
||||
t0[i] = global->GetT0Bin(i);
|
||||
}
|
||||
}
|
||||
// 2nd: check in the run block
|
||||
for (unsigned int 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.";
|
||||
cerr << endl << ">> Please check your msr-file!!" << endl;
|
||||
} else {
|
||||
t0[i] = run.GetT0Bin(i);
|
||||
}
|
||||
}
|
||||
// if still some t0's are == -1, estimate t0
|
||||
unsigned int idx;
|
||||
double max;
|
||||
for (unsigned int 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++) {
|
||||
if (rawRunData->GetDataBin(histoList[i])->at(j) > max) {
|
||||
max = rawRunData->GetDataBin(histoList[i])->at(j);
|
||||
idx = j;
|
||||
}
|
||||
}
|
||||
cout << endl << ">> estimated t0=" << idx << endl;
|
||||
t0[i] = idx;
|
||||
}
|
||||
}
|
||||
|
||||
// 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++) {
|
||||
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);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
rd.rawData.clear();
|
||||
rd.rawData = data;
|
||||
rd.t0 = (int)t0[0];
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
||||
//---------------------------------------------------
|
||||
double millitime()
|
||||
{
|
||||
@ -686,6 +797,9 @@ double millitime()
|
||||
*/
|
||||
int main(int argc, char *argv[])
|
||||
{
|
||||
Int_t unitTag = FOURIER_UNIT_NOT_GIVEN;
|
||||
Int_t apodTag = F_APODIZATION_NONE;
|
||||
|
||||
// only program name alone
|
||||
if (argc == 1) {
|
||||
musrFT_syntax();
|
||||
@ -696,7 +810,7 @@ int main(int argc, char *argv[])
|
||||
// init startupParam
|
||||
musrFT_init(startupParam);
|
||||
|
||||
// parse options
|
||||
// parse command line options
|
||||
int status = musrFT_parse_options(argc, argv, startupParam);
|
||||
if (status != 0) {
|
||||
int retVal = PMUSR_SUCCESS;
|
||||
@ -755,7 +869,6 @@ int main(int argc, char *argv[])
|
||||
PPrepFourier data(startupParam.bkg, startupParam.packing);
|
||||
|
||||
// load msr-file(s)
|
||||
cout << endl << "debug> before reading msr-files ...";
|
||||
vector<PMsrHandler*> msrHandler;
|
||||
msrHandler.resize(startupParam.msrFln.size());
|
||||
for (unsigned int i=0; i<startupParam.msrFln.size(); i++) {
|
||||
@ -777,7 +890,6 @@ int main(int argc, char *argv[])
|
||||
}
|
||||
}
|
||||
|
||||
cout << endl << "debug> before reading data-files ...";
|
||||
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
|
||||
@ -791,110 +903,209 @@ 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++) {
|
||||
cout << endl << "debug> dataFln[" << i-msrHandler.size() << "]=" << startupParam.dataFln[i-msrHandler.size()];
|
||||
// cout << endl << "debug> dataFln[" << i-msrHandler.size() << "]=" << startupParam.dataFln[i-msrHandler.size()];
|
||||
// create run data handler
|
||||
if (startupHandler)
|
||||
runDataHandler[i] = new PRunDataHandler(startupParam.dataFln[i-msrHandler.size()], startupParam.dataFileFormat[i-msrHandler.size()], startupHandler->GetDataPathList());
|
||||
else
|
||||
runDataHandler[i] = new PRunDataHandler(startupParam.dataFln[i-msrHandler.size()], startupParam.dataFileFormat[i-msrHandler.size()]);
|
||||
}
|
||||
|
||||
// read the data files
|
||||
// read all the data files
|
||||
for (unsigned int i=0; i<runDataHandler.size(); i++) {
|
||||
runDataHandler[i]->ReadData();
|
||||
|
||||
if (!runDataHandler[i]->IsAllDataAvailable()) {
|
||||
cerr << endl << ">> musrFT **ERROR** couldn't read data-file '" << startupParam.dataFln[i-msrHandler.size()] << "'.";
|
||||
if (i < msrHandler.size()) {
|
||||
cerr << endl << ">> musrFT **ERROR** couldn't read data from msr-file '" << startupParam.msrFln[i] << "'." << endl;
|
||||
} else {
|
||||
cerr << endl << ">> musrFT **ERROR** couldn't read data-file '" << startupParam.dataFln[i] << "'." << endl;
|
||||
}
|
||||
return PMUSR_DATA_FILE_READ_ERROR;
|
||||
}
|
||||
|
||||
// dig out all the necessary time domain data
|
||||
PRawRunData *rawRunData = runDataHandler[i]->GetRunData();
|
||||
if (rawRunData == 0) {
|
||||
cerr << endl << ">> musrFT **ERROR** couldn't obtain the raw run data set for " << startupParam.dataFln[i-msrHandler.size()] << endl;
|
||||
if (i < msrHandler.size()) {
|
||||
cerr << endl << ">> musrFT **ERROR** couldn't obtain the raw run data set from msr-file " << startupParam.msrFln[i] << endl;
|
||||
} else {
|
||||
cerr << endl << ">> musrFT **ERROR** couldn't obtain the raw run data set for " << startupParam.dataFln[i-msrHandler.size()] << endl;
|
||||
}
|
||||
return PMUSR_DATA_FILE_READ_ERROR;
|
||||
}
|
||||
|
||||
// first check of histo list makes sense
|
||||
for (unsigned int j=0; j<startupParam.histo.size(); j++) {
|
||||
if ((unsigned int)startupParam.histo[j] > rawRunData->GetNoOfHistos()) {
|
||||
cerr << endl << ">> musrFT **ERROR** found histo no " << startupParam.histo[j] << " > # of histo in the file (";
|
||||
cerr << startupParam.dataFln[i-msrHandler.size()] << " // # histo: " << rawRunData->GetNoOfHistos() << ")." << endl;
|
||||
return PMUSR_DATA_FILE_READ_ERROR;
|
||||
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()) {
|
||||
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;
|
||||
}
|
||||
}
|
||||
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++)
|
||||
startupParam.histo.push_back(j+1);
|
||||
}
|
||||
}
|
||||
|
||||
musrFT_data rd;
|
||||
TString str("");
|
||||
TString str(""), fln("");
|
||||
unsigned int idx=0;
|
||||
musrFT_getMetaInfo(startupParam.dataFln[i-msrHandler.size()], rawRunData, str);
|
||||
for (unsigned int j=0; j<startupParam.histo.size(); j++) {
|
||||
idx = startupParam.histo[j];
|
||||
// meta information
|
||||
rd.info = TString::Format("h%d:", idx);
|
||||
rd.info += str;
|
||||
/*
|
||||
cout << endl << "debug> ++++++++";
|
||||
cout << endl << "debug> info: " << rd.info;
|
||||
*/
|
||||
// time resolution
|
||||
rd.timeResolution = rawRunData->GetTimeResolution() / 1.0e3; // time resolution in (us)
|
||||
// cout << endl << "debug> time resolution: " << rd.timeResolution << " (us)";
|
||||
// time range
|
||||
rd.timeRange[0] = startupParam.timeRange[0]; // in (us)
|
||||
rd.timeRange[1] = startupParam.timeRange[1]; // in (us)
|
||||
// cout << endl << "debug> time range: " << rd.timeRange[0] << ".." << rd.timeRange[1] << " (us)";
|
||||
// data set
|
||||
rd.rawData.clear();
|
||||
rd.rawData = *(rawRunData->GetDataBin(idx));
|
||||
// check for potentially given t0's
|
||||
rd.t0 = -1;
|
||||
if (startupParam.t0.size() == 1)
|
||||
rd.t0 = startupParam.t0[0];
|
||||
else if (j < startupParam.t0.size())
|
||||
rd.t0 = startupParam.t0[j];
|
||||
// cout << endl << "debug> t0 = " << rd.t0;
|
||||
if (rd.t0 == -1) { // no t0 given, try to estimate it
|
||||
musrFT_estimateT0(rd);
|
||||
// get meta info, time resolution, time range, raw data sets
|
||||
if (i < msrHandler.size()) { // obtain info from msr-files
|
||||
// keep PLOT block info
|
||||
PMsrPlotList *plot = msrHandler[i]->GetMsrPlotList();
|
||||
if (plot == 0) {
|
||||
cerr << endl << ">> musrFT **ERROR** couldn't obtain PLOT block from msr-handler." << endl;
|
||||
return PMUSR_DATA_FILE_READ_ERROR;
|
||||
}
|
||||
// keep RUN block(s) info
|
||||
PMsrRunList *runs = msrHandler[i]->GetMsrRunList();
|
||||
if (runs == 0) {
|
||||
cerr << endl << ">> musrFT **ERROR** couldn't obtain RUN block(s) from msr-handler." << endl;
|
||||
return PMUSR_DATA_FILE_READ_ERROR;
|
||||
}
|
||||
// keep GLOBAL block info
|
||||
PMsrGlobalBlock *global = msrHandler[i]->GetMsrGlobal();
|
||||
if (global == 0) {
|
||||
cerr << endl << ">> musrFT **ERROR** couldn't obtain GLOBAL block from msr-handler." << endl;
|
||||
return PMUSR_DATA_FILE_READ_ERROR;
|
||||
}
|
||||
// keep FOURIER block info
|
||||
PMsrFourierStructure *fourierBlock = msrHandler[i]->GetMsrFourierList();
|
||||
if (fourierBlock == 0) {
|
||||
cerr << endl << ">> msrFT **WARNING** couldn't obtain FOURIER block from msr-handler." << endl;
|
||||
return PMUSR_DATA_FILE_READ_ERROR;
|
||||
} else { // filter out all necessary info
|
||||
if (fourierBlock->fFourierBlockPresent) {
|
||||
// get units
|
||||
unitTag = fourierBlock->fUnits;
|
||||
// get fourier power
|
||||
if (startupParam.fourierPower == -1) { // no Fourier power given from the command line, hence check FOURIER block
|
||||
if (fourierBlock->fFourierPower > 1)
|
||||
startupParam.fourierPower = fourierBlock->fFourierPower;
|
||||
}
|
||||
// get apodization tag
|
||||
apodTag = fourierBlock->fApodization;
|
||||
// get range
|
||||
if ((startupParam.fourierRange[0] == -1) && (startupParam.fourierRange[1] == -1)) { // no Fourier range given from the command line
|
||||
startupParam.fourierRange[0] = fourierBlock->fPlotRange[0];
|
||||
startupParam.fourierRange[1] = fourierBlock->fPlotRange[1];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// get the run information from the msr-file PLOT block 'runs'
|
||||
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++) {
|
||||
|
||||
// keep forward histo list
|
||||
PIntVector histoList;
|
||||
for (unsigned int k=0; k<runs->at(runList[j]-1).GetForwardHistoNoSize(); k++) {
|
||||
histoList.push_back(runs->at(runList[j]-1).GetForwardHistoNo(k));
|
||||
}
|
||||
|
||||
// handle meta information
|
||||
fln = runDataHandler[i]->GetRunPathName();
|
||||
musrFT_getMetaInfo(fln, rawRunData, str);
|
||||
TString hh("");
|
||||
hh = TString::Format("h%d", histoList[0]);
|
||||
for (unsigned int k=1; k<histoList.size(); k++)
|
||||
hh += TString::Format("/%d", histoList[k]);
|
||||
hh += ":";
|
||||
rd.info = hh;
|
||||
rd.info += str;
|
||||
|
||||
// handle time resolution
|
||||
rd.timeResolution = rawRunData->GetTimeResolution() / 1.0e3; // time resolution in (us)
|
||||
|
||||
// handle time range
|
||||
// take it from msr-file PLOT block 'range' if not overwritten from the command line
|
||||
if ((startupParam.timeRange[0] != -1) && (startupParam.timeRange[1] != -1)) {
|
||||
rd.timeRange[0] = startupParam.timeRange[0];
|
||||
rd.timeRange[1] = startupParam.timeRange[1];
|
||||
} else {
|
||||
rd.timeRange[0] = plot->at(0).fTmin[0];
|
||||
rd.timeRange[1] = plot->at(0).fTmax[0];
|
||||
}
|
||||
|
||||
// handle data set(s)
|
||||
// group forward histos
|
||||
if (musrFT_groupHistos(runDataHandler[i], global, runs->at(runList[j]-1), rd)) {
|
||||
return PMUSR_DATA_FILE_READ_ERROR;
|
||||
}
|
||||
// keep data set
|
||||
data.AddData(rd);
|
||||
|
||||
// get packing
|
||||
Int_t pack = 1;
|
||||
if (global->GetPacking() != -1) {
|
||||
pack = global->GetPacking();
|
||||
}
|
||||
if (runs->at(runList[j]-1).GetPacking() != -1) {
|
||||
pack = runs->at(runList[j]-1).GetPacking();
|
||||
}
|
||||
if (startupParam.packing > 1)
|
||||
pack = startupParam.packing;
|
||||
data.SetPacking(pack);
|
||||
}
|
||||
} 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++) {
|
||||
idx = startupParam.histo[j];
|
||||
|
||||
// handle meta information
|
||||
rd.info = TString::Format("h%d:", idx);
|
||||
rd.info += str;
|
||||
|
||||
// handle time resolution
|
||||
rd.timeResolution = rawRunData->GetTimeResolution() / 1.0e3; // time resolution in (us)
|
||||
|
||||
// handle time range
|
||||
rd.timeRange[0] = startupParam.timeRange[0]; // in (us)
|
||||
rd.timeRange[1] = startupParam.timeRange[1]; // in (us)
|
||||
|
||||
// handle data set
|
||||
rd.rawData.clear();
|
||||
rd.rawData = *(rawRunData->GetDataBin(idx));
|
||||
|
||||
// handle t0's
|
||||
rd.t0 = -1;
|
||||
if (startupParam.t0.size() == 1)
|
||||
rd.t0 = startupParam.t0[0];
|
||||
else if (j < startupParam.t0.size())
|
||||
rd.t0 = startupParam.t0[j];
|
||||
if (rd.t0 == -1) { // no t0 given, try to estimate it
|
||||
musrFT_estimateT0(rd);
|
||||
}
|
||||
|
||||
data.AddData(rd);
|
||||
}
|
||||
/*
|
||||
cout << endl << "debug> rd.rawData.size()=" << rd.rawData.size();
|
||||
cout << endl << "debug> rd.rawData: h" << idx << ": ";
|
||||
*/
|
||||
for (unsigned int k=rd.t0; k<rd.t0+15; k++)
|
||||
cout << rd.rawData[k] << ", ";
|
||||
data.AddData(rd);
|
||||
}
|
||||
}
|
||||
|
||||
/*
|
||||
cout << endl << "debug> +++++++++";
|
||||
cout << endl << "debug> data.GetNoOfData()=" << data.GetNoOfData();
|
||||
*/
|
||||
|
||||
// calculate background levels and subtract them from the data
|
||||
data.DoBkgCorrection();
|
||||
|
||||
// do the time domain filtering now
|
||||
data.DoFiltering();
|
||||
|
||||
// do lifetime correction
|
||||
if (startupParam.lifetimecorrection != 0.0)
|
||||
data.DoLifeTimeCorrection(startupParam.lifetimecorrection);
|
||||
|
||||
// do packing
|
||||
data.DoPacking();
|
||||
|
||||
// get all the corrected data
|
||||
vector<TH1F*> histo = data.GetData();
|
||||
/*
|
||||
if (histo.size() == 0)
|
||||
cout << endl << "debug> histo.size()==0!! :-(" << endl;
|
||||
for (unsigned int i=0; i<histo.size(); i++) {
|
||||
cout << endl << "debug> histo: time resolution " << histo[i]->GetBinWidth(1) << " (ns).";
|
||||
cout << endl << "debug> histo[" << i << "]->GetNbinsX()=" << histo[i]->GetNbinsX();
|
||||
cout << endl << "debug> histo data: ";
|
||||
for (unsigned int j=1; j<=15; j++)
|
||||
cout << histo[i]->GetBinContent(j) << ", ";
|
||||
cout << endl;
|
||||
}
|
||||
*/
|
||||
|
||||
// prepare Fourier
|
||||
Int_t unitTag = FOURIER_UNIT_NOT_GIVEN;
|
||||
if (startupParam.fourierUnits.BeginsWith("gauss", TString::kIgnoreCase))
|
||||
unitTag = FOURIER_UNIT_GAUSS;
|
||||
else if (startupParam.fourierUnits.BeginsWith("tesla", TString::kIgnoreCase))
|
||||
@ -903,6 +1114,8 @@ int main(int argc, char *argv[])
|
||||
unitTag = FOURIER_UNIT_FREQ;
|
||||
else if (startupParam.fourierUnits.BeginsWith("mc/s", TString::kIgnoreCase))
|
||||
unitTag = FOURIER_UNIT_CYCLES;
|
||||
else if (startupParam.fourierUnits.BeginsWith("??", TString::kIgnoreCase) && (unitTag == FOURIER_UNIT_NOT_GIVEN))
|
||||
unitTag = FOURIER_UNIT_FREQ;
|
||||
|
||||
vector<PFourier*> fourier;
|
||||
fourier.resize(histo.size());
|
||||
@ -911,7 +1124,6 @@ int main(int argc, char *argv[])
|
||||
}
|
||||
|
||||
// Fourier transform data
|
||||
Int_t apodTag = F_APODIZATION_NONE;
|
||||
if (startupParam.apotization.BeginsWith("weak", TString::kIgnoreCase))
|
||||
apodTag = F_APODIZATION_WEAK;
|
||||
else if (startupParam.apotization.BeginsWith("medium", TString::kIgnoreCase))
|
||||
@ -919,6 +1131,8 @@ int main(int argc, char *argv[])
|
||||
else if (startupParam.apotization.BeginsWith("strong", TString::kIgnoreCase))
|
||||
apodTag = F_APODIZATION_STRONG;
|
||||
|
||||
cout << endl << "debug> apodTag = " << apodTag << endl;
|
||||
|
||||
double start = millitime();
|
||||
for (unsigned int i=0; i<fourier.size(); i++) {
|
||||
fourier[i]->Transform(apodTag);
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user