Merged muonspin/musrfit into master

ChangeLog

@@ -4,6 +4,10 @@
 
 changes since 0.13.0
 ===================================
+NEW 2015-02-13 first implementation of a standalone Fourier transform/plotter: 
+               musrFT. Initially it is meant to be used for HAL-9500, 
+               i.e. Fourier transform WITHOUT lifetime correction.
+               A first simple minded lifetime correction is implemented as well.                
 NEW 2014-12-18 first implementation of a GLOBAL block which allows to shorten
                a typical msr-file. Duplicate entries from the RUN blocks can be
                added here. Furthermore, the 'lifetimecorrection' flag is 

src/Makefile.am

@@ -34,13 +34,14 @@ SUBDIRS += $(EDITORDIR)
 EXTRA_DIST = $(EDITORDIR)/Makefile
 endif
 
-bin_PROGRAMS = musrfit musrview musrt0 msr2msr msr2data any2many
+bin_PROGRAMS = musrfit musrview musrt0 musrFT msr2msr msr2data any2many
 bin_PROGRAMS += write_musrRoot_runHeader musrRootValidation
 bin_PROGRAMS += dump_header
 
 musrfit_SOURCES = musrfit.cpp
 musrview_SOURCES = musrview.cpp
 musrt0_SOURCES = musrt0.cpp
+musrFT_SOURCES = musrFT.cpp
 msr2msr_SOURCES = msr2msr.cpp
 msr2data_SOURCES = msr2data.cpp
 any2many_SOURCES = any2many.cpp

src/classes/Makefile.am

@@ -4,6 +4,7 @@ h_sources = \
   ../include/PFitterFcn.h \
   ../include/PFitter.h \
   ../include/PFourier.h \
+  ../include/PFourierCanvas.h \
   ../include/PFunctionGrammar.h \
   ../include/PFunction.h \
   ../include/PFunctionHandler.h \
@@ -12,6 +13,7 @@ h_sources = \
   ../include/PMusrCanvas.h \
   ../include/PMusr.h \
   ../include/PMusrT0.h \
+  ../include/PPrepFourier.h \
   ../include/PRunAsymmetry.h \
   ../include/PRunBase.h \
   ../include/PRunDataHandler.h \
@@ -26,6 +28,7 @@ h_sources_userFcn = \
   ../include/PUserFcnBase.h
 
 h_linkdef = \
+  ../include/PFourierCanvasLinkDef.h \
   ../include/PMusrCanvasLinkDef.h \
   ../include/PMusrT0LinkDef.h \
   ../include/PStartupHandlerLinkDef.h
@@ -34,6 +37,7 @@ h_linkdef_userFcn = \
   ../include/PUserFcnBaseLinkDef.h
 
 dict_h_sources = \
+  PFourierCanvasDict.h \
   PMusrCanvasDict.h \
   PMusrT0Dict.h \
   PStartupHandlerDict.h
@@ -45,6 +49,7 @@ cpp_sources = \
   PFitter.cpp \
   PFitterFcn.cpp \
   PFourier.cpp \
+  PFourierCanvas.cpp \
   PFunction.cpp \
   PFunctionHandler.cpp \
   PMsr2Data.cpp \
@@ -52,6 +57,7 @@ cpp_sources = \
   PMusrCanvas.cpp \
   PMusr.cpp \
   PMusrT0.cpp \
+  PPrepFourier.cpp \
   PRunAsymmetry.cpp \
   PRunBase.cpp \
   PRunDataHandler.cpp \
@@ -66,6 +72,7 @@ cpp_sources_userFcn = \
   PUserFcnBase.cpp
 
 dict_cpp_sources = \
+  PFourierCanvasDict.cpp \
   PMusrCanvasDict.cpp \
   PMusrT0Dict.cpp \
   PStartupHandlerDict.cpp

src/classes/PMsrHandler.cpp

@@ -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,6 +3169,7 @@ Bool_t PMsrHandler::HandleRunEntry(PMsrLines &lines)
         }
       }
       // check map entries, i.e. if the map values are within parameter bounds
+      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!";
@@ -3165,6 +3178,7 @@ Bool_t PMsrHandler::HandleRunEntry(PMsrLines &lines)
           }
         }
       }
+    }
 
     // forward ------------------------------------------------
     if (iter->fLine.BeginsWith("forward", TString::kIgnoreCase)) {
@@ -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

src/classes/PMusrCanvas.cpp

@@ -3797,15 +3797,10 @@ void PMusrCanvas::HandleAverage()
   // calculate all the average data sets
   double dval;
   if (fDataAvg.data != 0) {
-    if (!CalcAlignment(eTime)) {
-      cerr << endl << ">> PMusrCanvas::HandleAverage: data: **WARNING** only approx. alignment possible." << endl;
-    }
     for (Int_t i=0; i<fData[0].data->GetNbinsX(); i++) {
       dval = 0.0;
       for (UInt_t j=0; j<fData.size(); j++) {
-        if ((i-fAlignmentOffset[j]) > 0) {
+        dval += GetInterpolatedValue(fData[j].data, fData[0].data->GetBinContent(i));
-          dval += fData[j].data->GetBinContent(i-fAlignmentOffset[j]);
-        }
       }
       fDataAvg.data->SetBinContent(i, dval/fData.size());
     }
@@ -3816,15 +3811,10 @@ void PMusrCanvas::HandleAverage()
     fDataAvg.data->SetMarkerStyle(fData[0].data->GetMarkerStyle());
   }
   if (fDataAvg.dataFourierRe != 0) {
-    if (!CalcAlignment(eFreq)) {
-      cerr << endl << ">> PMusrCanvas::HandleAverage: Fourier Re: **WARNING** only approx. alignment possible." << endl;
-    }
     for (Int_t i=0; i<fData[0].dataFourierRe->GetNbinsX(); i++) {
       dval = 0.0;
       for (UInt_t j=0; j<fData.size(); j++) {
-        if ((i-fAlignmentOffset[j]) > 0) {
+        dval += GetInterpolatedValue(fData[j].dataFourierRe, fData[0].dataFourierRe->GetBinContent(i));
-          dval += fData[j].dataFourierRe->GetBinContent(i-fAlignmentOffset[j]);
-        }
       }
       fDataAvg.dataFourierRe->SetBinContent(i, dval/fData.size());
     }
@@ -3835,15 +3825,10 @@ void PMusrCanvas::HandleAverage()
     fDataAvg.dataFourierRe->SetMarkerStyle(fData[0].dataFourierRe->GetMarkerStyle());
   }
   if (fDataAvg.dataFourierIm != 0) {
-    if (!CalcAlignment(eFreq)) {
-      cerr << endl << ">> PMusrCanvas::HandleAverage: Fourier Im: **WARNING** only approx. alignment possible." << endl;
-    }
     for (Int_t i=0; i<fData[0].dataFourierIm->GetNbinsX(); i++) {
       dval = 0.0;
       for (UInt_t j=0; j<fData.size(); j++) {
-        if ((i-fAlignmentOffset[j]) > 0) {
+        dval += GetInterpolatedValue(fData[j].dataFourierIm, fData[0].dataFourierIm->GetBinContent(i));
-          dval += fData[j].dataFourierIm->GetBinContent(i-fAlignmentOffset[j]);
-        }
       }
       fDataAvg.dataFourierIm->SetBinContent(i, dval/fData.size());
     }
@@ -3854,15 +3839,10 @@ void PMusrCanvas::HandleAverage()
     fDataAvg.dataFourierIm->SetMarkerStyle(fData[0].dataFourierIm->GetMarkerStyle());
   }
   if (fDataAvg.dataFourierPwr != 0) {
-    if (!CalcAlignment(eFreq)) {
-      cerr << endl << ">> PMusrCanvas::HandleAverage: Fourier Pwr: **WARNING** only approx. alignment possible." << endl;
-    }
     for (Int_t i=0; i<fData[0].dataFourierPwr->GetNbinsX(); i++) {
       dval = 0.0;
       for (UInt_t j=0; j<fData.size(); j++) {
-        if ((i-fAlignmentOffset[j]) > 0) {
+        dval += GetInterpolatedValue(fData[j].dataFourierPwr, fData[0].dataFourierPwr->GetBinContent(i));
-          dval += fData[j].dataFourierPwr->GetBinContent(i-fAlignmentOffset[j]);
-        }
       }
       fDataAvg.dataFourierPwr->SetBinContent(i, dval/fData.size());
     }
@@ -3873,15 +3853,10 @@ void PMusrCanvas::HandleAverage()
     fDataAvg.dataFourierPwr->SetMarkerStyle(fData[0].dataFourierPwr->GetMarkerStyle());
   }
   if (fDataAvg.dataFourierPhase != 0) {
-    if (!CalcAlignment(eFreq)) {
-      cerr << endl << ">> PMusrCanvas::HandleAverage: Fourier Phase: **WARNING** only approx. alignment possible." << endl;
-    }
     for (Int_t i=0; i<fData[0].dataFourierPhase->GetNbinsX(); i++) {
       dval = 0.0;
       for (UInt_t j=0; j<fData.size(); j++) {
-        if ((i-fAlignmentOffset[j]) > 0) {
+        dval += GetInterpolatedValue(fData[j].dataFourierPhase, fData[0].dataFourierPhase->GetBinContent(i));
-          dval += fData[j].dataFourierPhase->GetBinContent(i-fAlignmentOffset[j]);
-        }
       }
       fDataAvg.dataFourierPhase->SetBinContent(i, dval/fData.size());
     }
@@ -3892,15 +3867,10 @@ void PMusrCanvas::HandleAverage()
     fDataAvg.dataFourierPhase->SetMarkerStyle(fData[0].dataFourierPhase->GetMarkerStyle());
   }
   if (fDataAvg.theory != 0) {
-    if (!CalcAlignment(eTheoTime)) {
-      cerr << endl << ">> PMusrCanvas::HandleAverage: theory: **WARNING** only approx. alignment possible." << endl;
-    }
     for (Int_t i=0; i<fData[0].theory->GetNbinsX(); i++) {
       dval = 0.0;
       for (UInt_t j=0; j<fData.size(); j++) {
-        if ((i-fAlignmentOffset[j]) > 0) {
+        dval += GetInterpolatedValue(fData[j].theory, fData[0].theory->GetBinContent(i));
-          dval += fData[j].theory->GetBinContent(i-fAlignmentOffset[j]);
-        }
       }
       fDataAvg.theory->SetBinContent(i, dval/fData.size());
     }
@@ -3910,7 +3880,7 @@ void PMusrCanvas::HandleAverage()
     for (Int_t i=0; i<fData[0].theoryFourierRe->GetNbinsX(); i++) {
       dval = 0.0;
       for (UInt_t j=0; j<fData.size(); j++) {
-        dval += fData[j].theoryFourierRe->GetBinContent(i);
+        dval += GetInterpolatedValue(fData[j].theoryFourierRe, fData[0].theoryFourierRe->GetBinContent(i));
       }
       fDataAvg.theoryFourierRe->SetBinContent(i, dval/fData.size());
     }
@@ -3924,7 +3894,7 @@ void PMusrCanvas::HandleAverage()
     for (Int_t i=0; i<fData[0].theoryFourierIm->GetNbinsX(); i++) {
       dval = 0.0;
       for (UInt_t j=0; j<fData.size(); j++) {
-        dval += fData[j].theoryFourierIm->GetBinContent(i);
+        dval += GetInterpolatedValue(fData[j].theoryFourierIm, fData[0].theoryFourierIm->GetBinContent(i));
       }
       fDataAvg.theoryFourierIm->SetBinContent(i, dval/fData.size());
     }
@@ -3938,7 +3908,7 @@ void PMusrCanvas::HandleAverage()
     for (Int_t i=0; i<fData[0].theoryFourierPwr->GetNbinsX(); i++) {
       dval = 0.0;
       for (UInt_t j=0; j<fData.size(); j++) {
-        dval += fData[j].theoryFourierPwr->GetBinContent(i);
+        dval += GetInterpolatedValue(fData[j].theoryFourierPwr, fData[0].theoryFourierPwr->GetBinContent(i));
       }
       fDataAvg.theoryFourierPwr->SetBinContent(i, dval/fData.size());
     }
@@ -3952,7 +3922,7 @@ void PMusrCanvas::HandleAverage()
     for (Int_t i=0; i<fData[0].theoryFourierPhase->GetNbinsX(); i++) {
       dval = 0.0;
       for (UInt_t j=0; j<fData.size(); j++) {
-        dval += fData[j].theoryFourierPhase->GetBinContent(i);
+        dval += GetInterpolatedValue(fData[j].theoryFourierPhase, fData[0].theoryFourierPhase->GetBinContent(i));
       }
       fDataAvg.theoryFourierPhase->SetBinContent(i, dval/fData.size());
     }
@@ -3963,15 +3933,10 @@ void PMusrCanvas::HandleAverage()
     fDataAvg.theoryFourierPhase->SetMarkerStyle(fData[0].theoryFourierPhase->GetMarkerStyle());
   }
   if (fDataAvg.diff != 0) {
-    if (!CalcAlignment(eTime)) {
-      cerr << endl << ">> PMusrCanvas::HandleAverage: diff: **WARNING** only approx. alignment possible." << endl;
-    }
     for (Int_t i=0; i<fData[0].diff->GetNbinsX(); i++) {
       dval = 0.0;
       for (UInt_t j=0; j<fData.size(); j++) {
-        if ((i-fAlignmentOffset[j]) > 0) {
+        dval += GetInterpolatedValue(fData[j].diff, fData[0].diff->GetBinContent(i));
-          dval += fData[j].diff->GetBinContent(i-fAlignmentOffset[j]);
-        }
       }
       fDataAvg.diff->SetBinContent(i, dval/fData.size());
     }
@@ -3985,7 +3950,7 @@ void PMusrCanvas::HandleAverage()
     for (Int_t i=0; i<fData[0].diffFourierRe->GetNbinsX(); i++) {
       dval = 0.0;
       for (UInt_t j=0; j<fData.size(); j++) {
-        dval += fData[j].diffFourierRe->GetBinContent(i);
+        dval += GetInterpolatedValue(fData[j].diffFourierRe, fData[0].diffFourierRe->GetBinContent(i));
       }
       fDataAvg.diffFourierRe->SetBinContent(i, dval/fData.size());
     }
@@ -3999,7 +3964,7 @@ void PMusrCanvas::HandleAverage()
     for (Int_t i=0; i<fData[0].diffFourierIm->GetNbinsX(); i++) {
       dval = 0.0;
       for (UInt_t j=0; j<fData.size(); j++) {
-        dval += fData[j].diffFourierIm->GetBinContent(i);
+        dval += GetInterpolatedValue(fData[j].diffFourierIm, fData[0].diffFourierIm->GetBinContent(i));
       }
       fDataAvg.diffFourierIm->SetBinContent(i, dval/fData.size());
     }
@@ -4013,7 +3978,7 @@ void PMusrCanvas::HandleAverage()
     for (Int_t i=0; i<fData[0].diffFourierPwr->GetNbinsX(); i++) {
       dval = 0.0;
       for (UInt_t j=0; j<fData.size(); j++) {
-        dval += fData[j].diffFourierPwr->GetBinContent(i);
+        dval += GetInterpolatedValue(fData[j].diffFourierPwr, fData[0].diffFourierPwr->GetBinContent(i));
       }
       fDataAvg.diffFourierPwr->SetBinContent(i, dval/fData.size());
     }
@@ -4027,7 +3992,7 @@ void PMusrCanvas::HandleAverage()
     for (Int_t i=0; i<fData[0].diffFourierPhase->GetNbinsX(); i++) {
       dval = 0.0;
       for (UInt_t j=0; j<fData.size(); j++) {
-        dval += fData[j].diffFourierPhase->GetBinContent(i);
+        dval += GetInterpolatedValue(fData[j].diffFourierPhase, fData[0].diffFourierPhase->GetBinContent(i));
       }
       fDataAvg.diffFourierPhase->SetBinContent(i, dval/fData.size());
     }
@@ -6314,130 +6279,40 @@ UInt_t PMusrCanvas::GetNeededAccuracy(PMsrParamStructure param)
 
 
 //--------------------------------------------------------------------------
-// CalcAlignment (private)
+// GetInterpolatedValue (private)
 //--------------------------------------------------------------------------
 /**
- * <p>Calculates the alignment index for each data set needed to average the data.
+ * <p>search for xVal in histo. If xVal is not found exactly, interpolate and
+ * return the interpolated y-value.
  *
  * <b>return:</b>
- * - true for perfect alignment
+ * - interpolated value if xVal is within histo range, 0 otherwise.
- * - false for approximate alignment
  *
- * \param tag to distinguish time data sets from Fourier data sets.
+ * \param histo pointer of the histogram
+ * \param xVal x-value to be looked for
  */
-Bool_t PMusrCanvas::CalcAlignment(const EAlignTag tag)
+Double_t PMusrCanvas::GetInterpolatedValue(TH1F* histo, Double_t xVal)
 {
-  Bool_t result = true;
+  if (histo == 0)
+    return 0.0;
 
-  fAlignmentOffset.clear();
+  Int_t idx = histo->FindBin(xVal);
 
-  UInt_t idx=0;
+  // make sure idx is within range
-  Double_t dval;
+  if ((idx < 1) || (idx > histo->GetNbinsX()))
-  if (tag == eTime) {
+    return 0.0;
-    // first find the data vector with the lowest initial time
-    dval = fData[0].data->GetXaxis()->GetBinCenter(1);
-    for (UInt_t i=1; i<fData.size(); i++) {
-      if (fData[i].data->GetXaxis()->GetBinCenter(1) < dval) {
-        idx = i;
-        dval = fData[i].data->GetXaxis()->GetBinCenter(1);
-      }
-    }
 
-    // next setp: find all the alignment indices
+  // make sure the lower bound idx is found. This is needed since
-    fAlignmentOffset.resize(fData.size());
+  // FindBin rounds towards the closer idx.
-    for (UInt_t i=0; i<fData.size(); i++) {
+  if (histo->GetBinCenter(idx) > xVal)
-      dval = fData[i].data->GetXaxis()->GetBinCenter(1);
+    --idx;
-      Int_t j=1;
-      Bool_t found = false;
-      do {
-        if (fData[idx].data->GetXaxis()->GetBinCenter(j) >= dval) {
-          fAlignmentOffset[i] = (UInt_t)(j-1);
-          if (fData[idx].data->GetXaxis()->GetBinCenter(j) != dval) {
-            result = false;
-          }
-          found = true;
-        }
-      } while (!found && (++j<fData[idx].data->GetNbinsX()));
-    }
-  } else if (tag == eTheoTime) {
-    // first find the data vector with the lowest initial time
-    dval = fData[0].theory->GetXaxis()->GetBinCenter(1);
-    for (UInt_t i=1; i<fData.size(); i++) {
-      if (fData[i].theory->GetXaxis()->GetBinCenter(1) < dval) {
-        idx = i;
-        dval = fData[i].theory->GetXaxis()->GetBinCenter(1);
-      }
-    }
 
-    // next setp: find all the alignment indices
+  Double_t x0, x1, y0, y1;
-    fAlignmentOffset.resize(fData.size());
+  x0 = histo->GetBinCenter(idx);
-    for (UInt_t i=0; i<fData.size(); i++) {
+  x1 = histo->GetBinCenter(idx+1);
-      dval = fData[i].theory->GetXaxis()->GetBinCenter(1);
+  y0 = histo->GetBinContent(idx);
-      Int_t j=1;
+  y1 = histo->GetBinContent(idx+1);
-      Bool_t found = false;
-      do {
-        if (fData[idx].theory->GetXaxis()->GetBinCenter(j) >= dval) {
-          fAlignmentOffset[i] = (UInt_t)(j-1);
-          if (fData[idx].theory->GetXaxis()->GetBinCenter(j) != dval) {
-            result = false;
-          }
-          found = true;
-        }
-      } while (!found && (++j<fData[idx].theory->GetNbinsX()));
-    }
-  } else if (tag == eFreq) {
-    // first find the data vector with the lowest initial time
-    dval = fData[0].dataFourierRe->GetXaxis()->GetBinCenter(1);
-    for (UInt_t i=1; i<fData.size(); i++) {
-      if (fData[i].dataFourierRe->GetXaxis()->GetBinCenter(1) < dval) {
-        idx = i;
-        dval = fData[i].dataFourierRe->GetXaxis()->GetBinCenter(1);
-      }
-    }
 
-    // next setp: find all the alignment indices
+  return (y1-y0)*(xVal-x0)/(x1-x0)+y0;
-    fAlignmentOffset.resize(fData.size());
-    for (UInt_t i=0; i<fData.size(); i++) {
-      dval = fData[i].dataFourierRe->GetXaxis()->GetBinCenter(1);
-      Int_t j=1;
-      Bool_t found = false;
-      do {
-        if (fData[idx].dataFourierRe->GetXaxis()->GetBinCenter(j) >= dval) {
-          fAlignmentOffset[i] = (UInt_t)(j-1);
-          if (fData[idx].dataFourierRe->GetXaxis()->GetBinCenter(j) != dval) {
-            result = false;
-          }
-          found = true;
-        }
-      } while (!found && (++j<fData[idx].dataFourierRe->GetNbinsX()));
-    }
-  } else if (tag == eTheoFreq) {
-    // first find the data vector with the lowest initial time
-    dval = fData[0].theoryFourierRe->GetXaxis()->GetBinCenter(1);
-    for (UInt_t i=1; i<fData.size(); i++) {
-      if (fData[i].theoryFourierRe->GetXaxis()->GetBinCenter(1) < dval) {
-        idx = i;
-        dval = fData[i].theoryFourierRe->GetXaxis()->GetBinCenter(1);
-      }
-    }
-
-    // next setp: find all the alignment indices
-    fAlignmentOffset.resize(fData.size());
-    for (UInt_t i=0; i<fData.size(); i++) {
-      dval = fData[i].theoryFourierRe->GetXaxis()->GetBinCenter(1);
-      Int_t j=1;
-      Bool_t found = false;
-      do {
-        if (fData[idx].theoryFourierRe->GetXaxis()->GetBinCenter(j) >= dval) {
-          fAlignmentOffset[i] = (UInt_t)(j-1);
-          if (fData[idx].theoryFourierRe->GetXaxis()->GetBinCenter(j) != dval) {
-            result = false;
-          }
-          found = true;
-        }
-      } while (!found && (++j<fData[idx].theoryFourierRe->GetNbinsX()));
-    }
-  }
-
-  return result;
 }
+

src/classes/PPrepFourier.cpp

@@ -0,0 +1,438 @@
+/***************************************************************************
+
+  PPrepFourier.cpp
+
+  Author: Andreas Suter
+  e-mail: andreas.suter@psi.ch
+
+***************************************************************************/
+
+/***************************************************************************
+ *   Copyright (C) 2007-2015 by Andreas Suter                              *
+ *   andreas.suter@psi.ch                                                  *
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ *                                                                         *
+ *   This program is distributed in the hope that it will be useful,       *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
+ *   GNU General Public License for more details.                          *
+ *                                                                         *
+ *   You should have received a copy of the GNU General Public License     *
+ *   along with this program; if not, write to the                         *
+ *   Free Software Foundation, Inc.,                                       *
+ *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
+ ***************************************************************************/
+
+#include "PPrepFourier.h"
+
+//--------------------------------------------------------------------------
+// Constructor
+//--------------------------------------------------------------------------
+/**
+ * <p>Constructor.
+ */
+PPrepFourier::PPrepFourier()
+{
+  fBkgRange[0] = -1;
+  fBkgRange[1] = -1;
+  fPacking = 1;
+}
+
+//--------------------------------------------------------------------------
+// Constructor
+//--------------------------------------------------------------------------
+/**
+ * <p>Constructor.
+ */
+PPrepFourier::PPrepFourier(const Int_t *bkgRange, const Int_t packing) :
+  fPacking(packing)
+{
+  SetBkgRange(bkgRange);
+}
+
+//--------------------------------------------------------------------------
+// Destructor
+//--------------------------------------------------------------------------
+/**
+ * <p>Destructor.
+ */
+PPrepFourier::~PPrepFourier()
+{
+  fRawData.clear();
+  fData.clear();
+}
+
+//--------------------------------------------------------------------------
+// SetBkgRange
+//--------------------------------------------------------------------------
+/**
+ * <p>set the background range.
+ *
+ * \param bkgRange array with background range
+ */
+void PPrepFourier::SetBkgRange(const Int_t *bkgRange)
+{
+  int err=0;
+  if (bkgRange[0] >= -1) {
+    fBkgRange[0] = bkgRange[0];
+  } else {
+    err = 1;
+  }
+  if (bkgRange[1] >= -1) {
+    fBkgRange[1] = bkgRange[1];
+  } else {
+    if (err == 0)
+      err = 2;
+    else
+      err = 3;
+  }
+
+  TString errMsg("");
+  switch (err) {
+    case 1:
+    errMsg = TString::Format("start bkg range < 0 (given: %d), will ignore it.", bkgRange[0]);
+      break;
+    case 2:
+      errMsg = TString::Format("end bkg range < 0 (given: %d), will ignore it.", bkgRange[1]);
+      break;
+    case 3:
+      errMsg = TString::Format("start/end bkg range < 0 (given: %d/%d), will ignore it.", bkgRange[0], bkgRange[1]);
+      break;
+    default:
+      errMsg = TString("??");
+      break;
+  }
+
+  if (err != 0) {
+    cerr << endl << ">> PPrepFourier::SetBkgRange: **WARNING** " << errMsg << endl;
+  }
+}
+
+//--------------------------------------------------------------------------
+// SetPacking
+//--------------------------------------------------------------------------
+/**
+ * <p>set the packing for the histograms.
+ *
+ * \param packing number to be used.
+ */
+void PPrepFourier::SetPacking(const Int_t packing)
+{
+  if (packing > 0) {
+    fPacking = packing;
+  } else {
+    cerr << endl << ">> PPrepFourier::SetPacking: **WARNING** found packing=" << packing << " < 0, will ignore it." << endl;
+  }
+}
+
+//--------------------------------------------------------------------------
+// AddData
+//--------------------------------------------------------------------------
+/**
+ * <p>add a data-set (time domain data + meta information) to the internal
+ * data vector.
+ *
+ * \param data set to be added
+ */
+void PPrepFourier::AddData(musrFT_data &data)
+{
+  fRawData.push_back(data);
+}
+
+//--------------------------------------------------------------------------
+// DoBkgCorrection
+//--------------------------------------------------------------------------
+/**
+ * <p>Correct the internal data sets according to a background interval given.
+ */
+void PPrepFourier::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)) {
+    return;
+  }
+
+  // make sure that the bkg range is ok
+  for (unsigned int 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;
+    }
+  }
+
+  Double_t bkg=0.0;
+  for (unsigned int i=0; i<fRawData.size(); i++) {
+    // calculate the bkg for the given range
+    for (int j=fBkgRange[0]; j<=fBkgRange[1]; j++) {
+      bkg += fRawData[i].rawData[j];
+    }
+    bkg /= (fBkgRange[1]-fBkgRange[0]+1);
+
+    // correct data
+    for (unsigned int j=0; j<fData[i].size(); j++)
+      fData[i][j] -= bkg;
+  }
+}
+
+//--------------------------------------------------------------------------
+// DoPacking
+//--------------------------------------------------------------------------
+/**
+ * <p>Rebin (pack) the internal data.
+ */
+void PPrepFourier::DoPacking()
+{
+  // 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;
+
+  PDoubleVector tmpData;
+  Double_t dval = 0.0;
+  for (unsigned int i=0; i<fData.size(); i++) {
+    tmpData.clear();
+    dval = 0.0;
+    for (unsigned int j=0; j<fData[i].size(); j++) {
+      if ((j % fPacking == 0) && (j != 0)) {
+        tmpData.push_back(dval);
+        dval = 0.0;
+      } else {
+        dval += fData[i][j];
+      }
+    }
+    // change the original data set with the packed one
+    fData[i].clear();
+    fData[i] = tmpData;
+  }
+}
+
+//--------------------------------------------------------------------------
+// DoFiltering
+//--------------------------------------------------------------------------
+/**
+ * <p>Not implemented yet.
+ */
+void PPrepFourier::DoFiltering()
+{
+  // make sure fData are already present, and if not create the necessary data sets
+  if (fData.size() != fRawData.size()) {
+    InitData();
+  }
+}
+
+//--------------------------------------------------------------------------
+// DoLifeTimeCorrection
+//--------------------------------------------------------------------------
+/**
+ * <p>Try to do a muon life time correction. The idea is to estimate N0 without
+ * any theory. This will be OK for high fields (> couple kGauss) but not so good
+ * for low fields.
+ *
+ * \param fudge rescaling factor for the estimated N0. Should be around 1
+ */
+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;
+    }
+  }
+}
+
+//--------------------------------------------------------------------------
+// GetInfo
+//--------------------------------------------------------------------------
+/**
+ * <p>Returns the meta information of a data set.
+ *
+ * \param idx index of the object
+ */
+TString PPrepFourier::GetInfo(const UInt_t idx)
+{
+  TString info("");
+
+  if (idx < fRawData.size())
+    info = fRawData[idx].info;
+
+  return info;
+}
+
+//--------------------------------------------------------------------------
+// GetData
+//--------------------------------------------------------------------------
+/**
+ * <p>Creates the requested TH1F objects and returns them. The ownership is with
+ * the caller.
+ */
+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;
+  Double_t end=0.0;
+  UInt_t size;
+  UInt_t startIdx;
+  UInt_t endIdx;
+
+  for (unsigned int i=0; i<fData.size(); i++) {
+    name = TString::Format("histo%2d", i);
+    dt = fRawData[i].timeResolution*fPacking;
+    start = fRawData[i].timeRange[0];
+    end = fRawData[i].timeRange[1];
+
+    // init size and start/end indices
+    size = fData[i].size();
+    startIdx = 1;
+    endIdx = size;
+
+    // time range given, hence calculate the proper size
+    if (start != -1.0) {
+      size = (UInt_t)((end-start)/dt);
+      if (start >= 0.0) {
+        startIdx = (UInt_t)(start/dt)+1;
+        endIdx = (UInt_t)(end/dt)+1;
+      } else {
+        cerr << endl << ">> PPrepFourier::GetData **WARNING** found start time < 0.0, will set it to 0.0" << endl;
+        endIdx = static_cast<UInt_t>(end/dt)+1;
+      }
+    }
+
+    if (start == -1.0) { // no time range given, hence start == 0.0 - dt/2
+      start = -dt/2.0;
+    } else { // time range given
+      start -= dt/2.0;
+    }
+
+    if (end == -1.0) { // no time range given, hence end == (fData[idx].size()-1)*dt + dt/2
+      end = (fData[i].size()-1)*dt+dt/2.0;
+    } else { // time range given
+      end += dt/2.0;
+    }
+
+    data[i] = new TH1F(name.Data(), fRawData[i].info.Data(), size, start, end);
+    for (unsigned int j=startIdx; j<endIdx; j++)
+      data[i]->SetBinContent(j-startIdx+1, fData[i][j]);
+  }
+
+  return data;
+}
+
+//--------------------------------------------------------------------------
+// GetData
+//--------------------------------------------------------------------------
+/**
+ * <p>Creates the requested TH1F object and returns it. The ownership is with
+ * the caller.
+ *
+ * \param idx index of the requested histogram
+ */
+TH1F *PPrepFourier::GetData(const UInt_t idx)
+{
+  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);
+  Double_t dt = fRawData[idx].timeResolution*fPacking;
+  Double_t start = fRawData[idx].timeRange[0];
+  Double_t end = fRawData[idx].timeRange[1];
+  UInt_t size = fData[idx].size();
+  UInt_t startIdx = 1;
+  UInt_t endIdx = size;
+
+  // time range given, hence calculate the proper size
+  if (start != -1.0) {
+    size = (UInt_t)((end-start)/dt);
+    if (start >= 0.0) {
+      startIdx = (UInt_t)(start/dt)+1;
+      endIdx = (UInt_t)(end/dt)+1;
+    } else {
+      cerr << endl << ">> PPrepFourier::GetData **WARNING** found start time < 0.0, will set it to 0.0" << endl;
+      endIdx = (UInt_t)(end/dt)+1;
+    }
+  }
+
+  if (start == -1.0) { // no time range given, hence start == 0.0 - dt/2
+    start = -dt/2.0;
+  } else { // time range given
+    start -= dt/2.0;
+  }
+
+  if (end == -1.0) { // no time range given, hence end == (fData[idx].size()-1)*dt + dt/2
+    end = (fData[idx].size()-1)*dt+dt/2.0;
+  } else { // time range given
+    end += dt/2.0;
+  }
+
+  TH1F *data = new TH1F(name.Data(), fRawData[idx].info.Data(), size, start, end);
+  for (unsigned int i=startIdx; i<endIdx; i++)
+    data->SetBinContent(i-startIdx+1, fData[idx][i]);
+
+  return data;
+}
+
+//--------------------------------------------------------------------------
+// InitData
+//--------------------------------------------------------------------------
+/**
+ * <p>Copy raw-data to internal data from t0 to the size of raw-data.
+ */
+void PPrepFourier::InitData()
+{
+  fData.resize(fRawData.size());
+  unsigned int t0;
+  for (unsigned int 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++) {
+      fData[i].push_back(fRawData[i].rawData[j]);
+    }
+  }
+}

src/classes/PRunDataHandler.cpp

@@ -218,6 +218,7 @@ PRunDataHandler::~PRunDataHandler()
 //--------------------------------------------------------------------------
 /**
  * <p>Checks if runName is found, and if so return these data.
+ * runName is as given in the msr-file.
  *
  * <b>return:</b>
  * - if data are found: pointer to the data.
@@ -240,6 +241,26 @@ PRawRunData* PRunDataHandler::GetRunData(const TString &runName)
     return &fData[i];
 }
 
+//--------------------------------------------------------------------------
+// GetRunData
+//--------------------------------------------------------------------------
+/**
+ * <p>return data-set with index idx.
+ *
+ * <b>return:</b>
+ * - if data are found: pointer to the data.
+ * - otherwise the null pointer will be returned.
+ *
+ * \param idx index of the raw data set.
+ */
+PRawRunData* PRunDataHandler::GetRunData(const UInt_t idx)
+{
+  if (idx >= fData.size())
+    return 0;
+  else
+    return &fData[idx];
+}
+
 //--------------------------------------------------------------------------
 // ReadData
 //--------------------------------------------------------------------------

src/include/PFourier.h

@@ -8,7 +8,7 @@
 ***************************************************************************/
 
 /***************************************************************************
- *   Copyright (C) 2007-2014 by Andreas Suter                              *
+ *   Copyright (C) 2007-2015 by Andreas Suter                              *
  *   andreas.suter@psi.ch                                                  *
  *                                                                         *
  *   This program is free software; you can redistribute it and/or modify  *
@@ -52,6 +52,8 @@ class PFourier
 
     virtual void Transform(UInt_t apodizationTag = 0);
 
+    virtual const char* GetDataTitle() { return fData->GetTitle(); }
+    virtual const Int_t GetUnitTag() { return fUnitTag; }
     virtual Double_t GetResolution() { return fResolution; }
     virtual TH1F* GetRealFourier(const Double_t scale = 1.0);
     virtual TH1F* GetImaginaryFourier(const Double_t scale = 1.0);
@@ -68,7 +70,7 @@ class PFourier
 
     Int_t fApodization; ///< 0=none, 1=weak, 2=medium, 3=strong
 
-    Double_t fTimeResolution; ///< time resolution of the data histogram
+    Double_t fTimeResolution; ///< time resolution of the data histogram in (us)
     Double_t fStartTime; ///< start time of the data histogram
     Double_t fEndTime; ///< end time of the data histogram
     Bool_t fDCCorrected; ///< if true, removed DC offset from signal before Fourier transformation, otherwise not

src/include/PFourierCanvas.h

@@ -0,0 +1,159 @@
+/***************************************************************************
+ *   Copyright (C) 2007-2015 by Andreas Suter                              *
+ *   andreas.suter@psi.ch                                                  *
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ *                                                                         *
+ *   This program is distributed in the hope that it will be useful,       *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
+ *   GNU General Public License for more details.                          *
+ *                                                                         *
+ *   You should have received a copy of the GNU General Public License     *
+ *   along with this program; if not, write to the                         *
+ *   Free Software Foundation, Inc.,                                       *
+ *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
+ ***************************************************************************/
+
+#ifndef _PFOURIERCANVAS_H_
+#define _PFOURIERCANVAS_H_
+
+#include <TObject.h>
+#include <TQObject.h>
+#include <TTimer.h>
+#include <TStyle.h>
+#include <TRootCanvas.h>
+#include <TGMenu.h>
+#include <TCanvas.h>
+#include <TPaveText.h>
+#include <TLegend.h>
+#include <TPad.h>
+#include <TH1F.h>
+#include <TLatex.h>
+
+#include "PMusr.h"
+#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_REAL          100
+#define P_MENU_ID_FOURIER_IMAG          101
+#define P_MENU_ID_FOURIER_REAL_AND_IMAG 102
+#define P_MENU_ID_FOURIER_PWR           103
+#define P_MENU_ID_FOURIER_PHASE         104
+#define P_MENU_ID_FOURIER_PHASE_PLUS    105
+#define P_MENU_ID_FOURIER_PHASE_MINUS   106
+
+//------------------------------------------------------------------------
+/**
+ * <p>Structure holding all necessary Fourier histograms.
+ */
+typedef struct {
+  TH1F *dataFourierRe;       ///< real part of the Fourier transform of the data histogram
+  TH1F *dataFourierIm;       ///< imaginary part of the Fourier transform of the data histogram
+  TH1F *dataFourierPwr;      ///< power spectrum of the Fourier transform of the data histogram
+  TH1F *dataFourierPhase;    ///< phase spectrum of the Fourier transform of the data histogram
+} PFourierCanvasDataSet;
+
+//------------------------------------------------------------------------
+/**
+ * <p>typedef to make to code more readable: list of histogram data sets.
+ */
+typedef vector<PFourierCanvasDataSet> PFourierCanvasDataList;
+
+//--------------------------------------------------------------------------
+/**
+ * <p>
+ */
+class PFourierCanvas : public TObject, public TQObject
+{
+  public:
+    PFourierCanvas();
+    PFourierCanvas(vector<PFourier*> &fourier, const Char_t* title, const Bool_t showAverage,
+                   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,
+                   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);
+    virtual ~PFourierCanvas();
+
+    virtual void Done(Int_t status=0); // *SIGNAL*
+    virtual void HandleCmdKey(Int_t event, Int_t x, Int_t y, TObject *selected); // SLOT
+    virtual void HandleMenuPopup(Int_t id); // SLOT
+    virtual void LastCanvasClosed(); // SLOT
+
+    virtual void UpdateFourierPad();
+    virtual void UpdateInfoPad();
+
+    virtual Bool_t IsValid() { return fValid; }
+
+    virtual void SetTimeout(Int_t ival);
+
+    virtual void SaveGraphicsAndQuit(const Char_t *fileName);
+    virtual void ExportData(const Char_t *pathFileName);
+
+  private:
+    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.
+    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
+
+    TString fTitle;
+    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
+    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
+
+    TStyle *fStyle; ///< A collection of all graphics attributes
+
+    TTimer *fTimeoutTimer;  ///< timeout timer in order to terminate if no action is taking place for too long
+
+    PIntVector fMarkerList; ///< list of markers
+    PIntVector fColorList;  ///< list of colors
+
+    // canvas menu related variables
+    TRootCanvas *fImp;           ///< ROOT native GUI version of main window with menubar and drawing area
+    TGMenuBar   *fBar;           ///< menu bar
+    TGPopupMenu *fPopupMain;     ///< popup menu MusrFT in the main menu bar
+//    TGPopupMenu *fPopupSave;     ///< popup menu of the MusrFT/Save Data sub menu
+    TGPopupMenu *fPopupFourier;  ///< popup menu of the MusrFT/Fourier sub menu
+
+    // canvas related variables
+    TCanvas   *fMainCanvas;        ///< main canvas
+    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
+
+    virtual void CreateXaxisTitle();
+    virtual void CreateStyle();
+    virtual void InitFourierDataSets();
+    virtual void InitFourierCanvas(const Char_t* title, Int_t wtopx, Int_t wtopy, Int_t ww, Int_t wh);
+    virtual void CleanupAverage();
+    virtual void HandleAverage();
+
+    virtual void PlotFourier();
+    virtual void PlotFourierPhaseValue();
+    virtual void PlotAverage();
+    virtual void IncrementFourierPhase();
+    virtual void DecrementFourierPhase();
+
+    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);
+
+    ClassDef(PFourierCanvas, 1)
+};
+
+#endif // _PFOURIERCANVAS_H_

src/include/PFourierCanvasLinkDef.h

@@ -0,0 +1,39 @@
+/***************************************************************************
+
+  PMusrCanvasLinkDef.h
+
+  Author: Andreas Suter
+  e-mail: andreas.suter@psi.ch
+
+***************************************************************************/
+
+/***************************************************************************
+ *   Copyright (C) 2007-2015 by Andreas Suter                              *
+ *   andreas.suter@psi.ch                                                  *
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ *                                                                         *
+ *   This program is distributed in the hope that it will be useful,       *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
+ *   GNU General Public License for more details.                          *
+ *                                                                         *
+ *   You should have received a copy of the GNU General Public License     *
+ *   along with this program; if not, write to the                         *
+ *   Free Software Foundation, Inc.,                                       *
+ *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
+ ***************************************************************************/
+
+#ifdef __CINT__
+
+#pragma link off all globals;
+#pragma link off all classes;
+#pragma link off all functions;
+
+#pragma link C++ class PFourierCanvas+;
+
+#endif
+

src/include/PMusr.h

@@ -8,7 +8,7 @@
 ***************************************************************************/
 
 /***************************************************************************
- *   Copyright (C) 2007-2014 by Andreas Suter                              *
+ *   Copyright (C) 2007-2015 by Andreas Suter                              *
  *   andreas.suter@psi.ch                                                  *
  *                                                                         *
  *   This program is free software; you can redistribute it and/or modify  *
@@ -44,6 +44,7 @@ using namespace std;
 #define PMUSR_TOKENIZE_ERROR           -5
 #define PMUSR_MSR_LOG_FILE_WRITE_ERROR -6
 #define PMUSR_MSR_FILE_WRITE_ERROR     -7
+#define PMUSR_DATA_FILE_READ_ERROR     -8
 
 #define PRUN_SINGLE_HISTO 0
 #define PRUN_ASYMMETRY    2
@@ -703,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)

src/include/PMusrCanvas.h

@@ -278,9 +278,7 @@ class PMusrCanvas : public TObject, public TQObject
     PRunListCollection *fRunList;    ///< data handler
 #endif // __MAKECINT__
 
-    enum EAlignTag {eTime, eTheoTime, eFreq, eTheoFreq};
     PMusrCanvasDataSet fDataAvg; ///< set of all averaged data to be plotted (asymmetry/single histogram)
-    PIntVector fAlignmentOffset; ///< holds the vector with the time/freq alignment offsets
 
     PMusrCanvasDataList fData; ///< list of all histogram data to be plotted (asymmetry/single histogram)
     PMusrCanvasNonMusrDataList fNonMusrData; ///< list of all error graphs to be plotted (non-muSR)
@@ -333,7 +331,7 @@ class PMusrCanvas : public TObject, public TQObject
     virtual Bool_t IsScaleN0AndBkg();
     virtual UInt_t GetNeededAccuracy(PMsrParamStructure param);
 
-    virtual Bool_t CalcAlignment(const EAlignTag tag);
+    virtual Double_t GetInterpolatedValue(TH1F* histo, Double_t xVal);
 
   ClassDef(PMusrCanvas, 1)
 };

src/include/PPrepFourier.h

@@ -0,0 +1,88 @@
+/***************************************************************************
+
+  PPrepFourier.h
+
+  Author: Andreas Suter
+  e-mail: andreas.suter@psi.ch
+
+***************************************************************************/
+
+/***************************************************************************
+ *   Copyright (C) 2007-2015 by Andreas Suter                              *
+ *   andreas.suter@psi.ch                                                  *
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ *                                                                         *
+ *   This program is distributed in the hope that it will be useful,       *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
+ *   GNU General Public License for more details.                          *
+ *                                                                         *
+ *   You should have received a copy of the GNU General Public License     *
+ *   along with this program; if not, write to the                         *
+ *   Free Software Foundation, Inc.,                                       *
+ *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
+ ***************************************************************************/
+
+#ifndef _PPREPFOURIER_H_
+#define _PPREPFOURIER_H_
+
+#include <iostream>
+#include <vector>
+using namespace std;
+
+#include "TH1F.h"
+
+#include "PMusr.h"
+
+//----------------------------------------------------------------------------
+/**
+ * <p>Data structure holding raw time domain uSR data together with some
+ * 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
+} musrFT_data;
+
+//----------------------------------------------------------------------------
+/**
+ * <p>Little helper class to prepare time-domain data for Fourier transform, without
+ * theory, etc.
+ */
+class PPrepFourier {
+  public:
+    PPrepFourier();
+    PPrepFourier(const Int_t *bkgRange, const Int_t packing);
+    virtual ~PPrepFourier();
+
+    void SetBkgRange(const Int_t *bkgRange);
+    void SetPacking(const Int_t packing);
+    void AddData(musrFT_data &data);
+    void DoBkgCorrection();
+    void DoPacking();
+    void DoFiltering();
+    void DoLifeTimeCorrection(Double_t fudge);
+
+    TString GetInfo(const UInt_t idx);
+    UInt_t GetNoOfData() { return fRawData.size(); }
+    vector<TH1F*> GetData();
+    TH1F *GetData(const UInt_t idx);
+
+  private:
+    vector<musrFT_data> fRawData;
+    vector<PDoubleVector>fData;
+    Int_t fBkgRange[2];
+    Int_t fPacking;
+
+    void InitData();
+};
+
+#endif // _PPREPFOURIER_H_
+

src/include/PRunDataHandler.h

@@ -8,7 +8,7 @@
 ***************************************************************************/
 
 /***************************************************************************
- *   Copyright (C) 2007-2014 by Andreas Suter                              *
+ *   Copyright (C) 2007-2015 by Andreas Suter                              *
  *   andreas.suter@psi.ch                                                  *
  *                                                                         *
  *   This program is free software; you can redistribute it and/or modify  *
@@ -60,6 +60,7 @@ class PRunDataHandler
 
     virtual Bool_t IsAllDataAvailable() const { return fAllDataAvailable; }
     virtual PRawRunData* GetRunData(const TString &runName);
+    virtual PRawRunData* GetRunData(const UInt_t idx=0);
     virtual TString GetRunPathName() {return fRunPathName; }
 
   private:

src/musrview.cpp

@@ -360,7 +360,6 @@ int main(int argc, char *argv[])
       sprintf(canvasName, "fMainCanvas%d", i);
       if (gROOT->GetListOfCanvases()->FindObject(canvasName) != 0) {
         canvasVector[i]->~PMusrCanvas();
-      } else {
       }
     }
     canvasVector.empty();