first implementation of NonMusr stuff

src/ToDo.txt

@@ -58,7 +58,8 @@ short term:
 
 * implement access to user-function, i.e. functions not
   defined within musrfit, using the ROOT dictionary feature
-  to look for them and use them. **UNDER WAY, NOT TESTED YET** 08-06-05
+  to look for them and use them. 
+  **DONE**
 
 * tshift implementation under way. Wrong and needs to be fixed: 
   **FIXED** 08-06-16
@@ -66,6 +67,10 @@ short term:
 * faulty dependency check in the Makefile when handling the PSI-Bin-Class stuff.
   Needs to be fixed! **FIXED** 08-06-16
 
+* implement NonMuSR stuff **DONE** 08-06-20, some more testing needed
+
+* implement FFT with msr-interface
+
 ---------------------
 intermediate term:
 ---------------------
@@ -74,7 +79,6 @@ intermediate term:
   describe new features with examples!
 * introduce error numbers with corresponding docu/explanation in the latex-docu.  
 * implement RRF stuff
-* implement NonMuSR stuff
 * think about if it is worth to modify wkm in order to read
   the new msr-files, so that wkmview can be used for an intermediate
   time.

src/classes/PMusrCanvas.cpp

@@ -129,6 +129,11 @@ cout << "~PMusrCanvas() called" << endl;
       CleanupDataSet(fData[i]);
     fData.clear();
   }
+  if (fNonMusrData.size() > 0) {
+    for (unsigned int i=0; i<fNonMusrData.size(); i++)
+      CleanupDataSet(fNonMusrData[i]);
+    fNonMusrData.clear();
+  }
 }
 
 //--------------------------------------------------------------------------
@@ -466,7 +471,7 @@ cout << endl;
           return;
         }
         // handle data
-        HandleDataSet(runNo, data);
+        HandleNonMusrDataSet(runNo, data);
         break;
       default:
         fValid = false;
@@ -480,25 +485,49 @@ cout << endl;
 
   // generate the histo plot
   fDataTheoryPad->cd();
-  if (fData.size() > 0) {
+  if (plotInfo.fPlotType != MSR_PLOT_NO_MUSR) {
-    fData[0].data->Draw("pe");
+    if (fData.size() > 0) {
-    // set time range
+      fData[0].data->Draw("pe");
-    Double_t xmin = fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fTmin;
+      // set time range
-    Double_t xmax = fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fTmax;
+      Double_t xmin = fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fTmin;
-    fData[0].data->GetXaxis()->SetRangeUser(xmin, xmax);
+      Double_t xmax = fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fTmax;
-    // check if it is necessary to set the y-axis range
+      fData[0].data->GetXaxis()->SetRangeUser(xmin, xmax);
-    Double_t ymin = fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fYmin;
+      // check if it is necessary to set the y-axis range
-    Double_t ymax = fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fYmax;
+      Double_t ymin = fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fYmin;
-    if ((ymin != -999.0) && (ymax != -999.0)) {
+      Double_t ymax = fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fYmax;
-      fData[0].data->GetYaxis()->SetRangeUser(ymin, ymax);
+      if ((ymin != -999.0) && (ymax != -999.0)) {
+        fData[0].data->GetYaxis()->SetRangeUser(ymin, ymax);
+      }
+      // plot all remaining data
+      for (unsigned int i=1; i<fData.size(); i++) {
+        fData[i].data->Draw("pesame");
+      }
+      // plot all the theory
+      for (unsigned int i=0; i<fData.size(); i++) {
+        fData[i].theory->Draw("csame");
+      }
     }
-    // plot all remaining data
+  } else { // plotInfo.fPlotType == MSR_PLOT_NO_MUSR
-    for (unsigned int i=1; i<fData.size(); i++) {
+    if (fNonMusrData.size() > 0) {
-      fData[i].data->Draw("pesame");
+      fNonMusrData[0].data->Draw("ap");
-    }
+      // set x-range
-    // plot all the theory
+      Double_t xmin = fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fTmin;
-    for (unsigned int i=0; i<fData.size(); i++) {
+      Double_t xmax = fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fTmax;
-      fData[i].theory->Draw("csame");
+      fNonMusrData[0].data->GetXaxis()->SetRangeUser(xmin, xmax);
+      // check if it is necessary to set the y-axis range
+      Double_t ymin = fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fYmin;
+      Double_t ymax = fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fYmax;
+      if ((ymin != -999.0) && (ymax != -999.0)) {
+        fNonMusrData[0].data->GetYaxis()->SetRangeUser(ymin, ymax);
+      }
+      // plot all remaining data
+      for (unsigned int i=1; i<fNonMusrData.size(); i++) {
+        fNonMusrData[i].data->Draw("psame");
+      }
+      // plot all the theory
+      for (unsigned int i=0; i<fNonMusrData.size(); i++) {
+        fNonMusrData[i].theory->Draw("csame");
+      }
     }
   }
   fMainCanvas->cd();
@@ -633,6 +662,33 @@ void PMusrCanvas::InitDataSet(PMusrCanvasDataSet &dataSet)
   dataSet.diffFourierPhase = 0;
 }
 
+//--------------------------------------------------------------------------
+// InitDataSet
+//--------------------------------------------------------------------------
+/**
+ * <p>
+ *
+ * \param dataSet
+ */
+void PMusrCanvas::InitDataSet(PMusrCanvasNonMusrDataSet &dataSet)
+{
+  dataSet.data = 0;
+  dataSet.dataFourierRe = 0;
+  dataSet.dataFourierIm = 0;
+  dataSet.dataFourierPwr = 0;
+  dataSet.dataFourierPhase = 0;
+  dataSet.theory = 0;
+  dataSet.theoryFourierRe = 0;
+  dataSet.theoryFourierIm = 0;
+  dataSet.theoryFourierPwr = 0;
+  dataSet.theoryFourierPhase = 0;
+  dataSet.diff = 0;
+  dataSet.diffFourierRe = 0;
+  dataSet.diffFourierIm = 0;
+  dataSet.diffFourierPwr = 0;
+  dataSet.diffFourierPhase = 0;
+}
+
 //--------------------------------------------------------------------------
 // CleanupDataSet
 //--------------------------------------------------------------------------
@@ -705,6 +761,78 @@ void PMusrCanvas::CleanupDataSet(PMusrCanvasDataSet &dataSet)
   }
 }
 
+//--------------------------------------------------------------------------
+// CleanupDataSet
+//--------------------------------------------------------------------------
+/**
+ * <p>
+ *
+ * \param dataSet
+ */
+void PMusrCanvas::CleanupDataSet(PMusrCanvasNonMusrDataSet &dataSet)
+{
+  if (dataSet.data) {
+    delete dataSet.data;
+    dataSet.data = 0;
+  }
+  if (dataSet.dataFourierRe) {
+    delete dataSet.dataFourierRe;
+    dataSet.dataFourierRe = 0;
+  }
+  if (dataSet.dataFourierIm) {
+    delete dataSet.dataFourierIm;
+    dataSet.dataFourierIm = 0;
+  }
+  if (dataSet.dataFourierPwr) {
+    delete dataSet.dataFourierPwr;
+    dataSet.dataFourierPwr = 0;
+  }
+  if (dataSet.dataFourierPhase) {
+    delete dataSet.dataFourierPhase;
+    dataSet.dataFourierPhase = 0;
+  }
+  if (dataSet.theory) {
+    delete dataSet.theory;
+    dataSet.theory = 0;
+  }
+  if (dataSet.theoryFourierRe) {
+    delete dataSet.theoryFourierRe;
+    dataSet.theoryFourierRe = 0;
+  }
+  if (dataSet.theoryFourierIm) {
+    delete dataSet.theoryFourierIm;
+    dataSet.theoryFourierIm = 0;
+  }
+  if (dataSet.theoryFourierPwr) {
+    delete dataSet.theoryFourierPwr;
+    dataSet.theoryFourierPwr = 0;
+  }
+  if (dataSet.theoryFourierPhase) {
+    delete dataSet.theoryFourierPhase;
+    dataSet.theoryFourierPhase = 0;
+  }
+  if (dataSet.diff) {
+    delete dataSet.diff;
+    dataSet.diff = 0;
+  }
+  if (dataSet.diffFourierRe) {
+    delete dataSet.diffFourierRe;
+    dataSet.diffFourierRe = 0;
+  }
+  if (dataSet.diffFourierIm) {
+    delete dataSet.diffFourierIm;
+    dataSet.diffFourierIm = 0;
+  }
+  if (dataSet.diffFourierPwr) {
+    delete dataSet.diffFourierPwr;
+    dataSet.diffFourierPwr = 0;
+  }
+  if (dataSet.diffFourierPhase) {
+    delete dataSet.diffFourierPhase;
+    dataSet.diffFourierPhase = 0;
+  }
+}
+
 //--------------------------------------------------------------------------
 // HandleDataSet
 //--------------------------------------------------------------------------
@@ -794,3 +922,78 @@ void PMusrCanvas::HandleDataSet(unsigned int runNo, PRunData *data)
 
   fData.push_back(dataSet);
 }
+
+//--------------------------------------------------------------------------
+// HandleNonMusrDataSet
+//--------------------------------------------------------------------------
+/**
+ * <p>
+ *
+ * \param runNo
+ * \param data
+ */
+void PMusrCanvas::HandleNonMusrDataSet(unsigned int runNo, PRunData *data)
+{
+  PMusrCanvasNonMusrDataSet dataSet;
+  TGraphErrors *dataHisto;
+  TGraphErrors *theoHisto;
+
+  InitDataSet(dataSet);
+
+  // dataHisto -------------------------------------------------------------
+
+  // invoke graph
+  dataHisto = new TGraphErrors(data->fX.size());
+
+  // fill graph
+  for (unsigned int i=0; i<data->fX.size(); i++) {
+    dataHisto->SetPoint(i, data->fX[i], data->fValue[i]);
+    dataHisto->SetPointError(i, 0.0, data->fError[i]);
+  }
+
+  // set marker and line color
+  if (runNo < fColorList.size()) {
+    dataHisto->SetMarkerColor(fColorList[runNo]);
+    dataHisto->SetLineColor(fColorList[runNo]);
+  } else {
+    TRandom rand(runNo);
+    Int_t color = TColor::GetColor((Int_t)rand.Integer(255), (Int_t)rand.Integer(255), (Int_t)rand.Integer(255));
+    dataHisto->SetMarkerColor(color);
+    dataHisto->SetLineColor(color);
+  }
+  // set marker size
+  dataHisto->SetMarkerSize(1);
+  // set marker type
+  if (runNo < fMarkerList.size()) {
+    dataHisto->SetMarkerStyle(fMarkerList[runNo]);
+  } else {
+    TRandom rand(runNo);
+    dataHisto->SetMarkerStyle(20+(Int_t)rand.Integer(10));
+  }
+
+  // theoHisto -------------------------------------------------------------
+
+  // invoke graph
+  theoHisto = new TGraphErrors(data->fXTheory.size());
+
+  // fill graph
+  for (unsigned int i=0; i<data->fXTheory.size(); i++) {
+    theoHisto->SetPoint(i, data->fXTheory[i], data->fTheory[i]);
+    theoHisto->SetPointError(i, 0.0, 0.0);
+  }
+
+  // set the line color
+  if (runNo < fColorList.size()) {
+    theoHisto->SetLineColor(fColorList[runNo]);
+  } else {
+    TRandom rand(runNo);
+    Int_t color = TColor::GetColor((Int_t)rand.Integer(255), (Int_t)rand.Integer(255), (Int_t)rand.Integer(255));
+    theoHisto->SetLineColor(color);
+  }
+
+  // fill handler list -----------------------------------------------------
+  dataSet.data   = dataHisto;
+  dataSet.theory = theoHisto;
+
+  fNonMusrData.push_back(dataSet);
+}

src/classes/PRunDataHandler.cpp

@@ -837,8 +837,6 @@ bool PRunDataHandler::ReadAsciiFile()
 {
   bool success = true;
 
-  cout << endl << "PRunDataHandler::ReadAsciiFile(): Sorry, not yet implemented ...";
-
   // open file
   ifstream f;
 

src/classes/PRunListCollection.cpp

@@ -318,7 +318,6 @@ PRunData* PRunListCollection::GetSingleHisto(unsigned int index, EDataSwitch tag
       }
       break;
     default: // error
-      return 0;
       break;
   }
 
@@ -357,7 +356,6 @@ PRunData* PRunListCollection::GetAsymmetry(unsigned int index, EDataSwitch tag)
       }
       break;
     default: // error
-      return 0;
       break;
   }
 
@@ -387,11 +385,10 @@ PRunData* PRunListCollection::GetRRF(unsigned int index, EDataSwitch tag)
     case kRunNo:
       break;
     default: // error
-      return 0;
       break;
   }
 
-  return 0;
+  return data;
 }
 
 //--------------------------------------------------------------------------
@@ -415,13 +412,18 @@ PRunData* PRunListCollection::GetNonMusr(unsigned int index, EDataSwitch tag)
       }
       break;
     case kRunNo:
+      for (unsigned int i=0; i<fRunNonMusrList.size(); i++) {
+        if (fRunNonMusrList[i]->GetRunNo() == index) {
+          data = fRunNonMusrList[i]->GetData();
+          break;
+        }
+      }
       break;
     default: // error
-      return 0;
       break;
   }
 
-  return 0;
+  return data;
 }
 
 //--------------------------------------------------------------------------

src/classes/PRunNonMusr.cpp

@@ -224,5 +224,111 @@ bool PRunNonMusr::PrepareViewData()
 {
   bool success = true;
 
+  // get the proper run
+  PRawRunData* runData = fRawData->GetRunData(fRunInfo->fRunName);
+  if (!runData) { // couldn't get run
+    cout << endl << "PRunNonMusr::PrepareViewData(): **ERROR** Couldn't get run " << fRunInfo->fRunName.Data() << "!";
+    return false;
+  }
+
+  // fill data histo
+  // pack the raw data
+  double value  = 0.0;
+  double err = 0.0;
+cout << endl << ">> runData->fXData.size()=" << runData->fXData.size();
+  for (unsigned int i=0; i<runData->fXData.size(); i++) {
+cout << endl << ">> i=" << i << ", packing=" << fRunInfo->fPacking;
+    if ((i % fRunInfo->fPacking == 0) && (i != 0)) { // fill data
+cout << endl << "-> i=" << i;
+      fData.fX.push_back(runData->fXData[i]-(runData->fXData[i]-runData->fXData[i-fRunInfo->fPacking])/2.0);
+      fData.fValue.push_back(value);
+      fData.fError.push_back(TMath::Sqrt(err));
+      value = 0.0;
+      err = 0.0;
+    }
+    // sum raw data values
+    value += runData->fYData[i];
+    err += runData->fErrYData[i]*runData->fErrYData[i];
+  }
+cout << endl << ">> fData.fValue.size()=" << fData.fValue.size();
+
+  // count the number of bins to be fitted
+  fNoOfFitBins = fData.fValue.size();
+cout << endl << ">> fNoOfFitBins=" << fNoOfFitBins;
+
+  // fill theory histo
+  // feed the parameter vector
+  std::vector<double> par;
+  PMsrParamList *paramList = fMsrInfo->GetMsrParamList();
+  for (unsigned int i=0; i<paramList->size(); i++)
+    par.push_back((*paramList)[i].fValue);
+  // calculate functions
+  for (int i=0; i<fMsrInfo->GetNoOfFuncs(); i++) {
+    fFuncValues[i] = fMsrInfo->EvalFunc(fMsrInfo->GetFuncNo(i), fRunInfo->fMap, par);
+  }
+
+cout << endl << ">> after parameter fill" << endl;
+  // get plot range
+  PMsrPlotList *plotList;
+  PMsrPlotStructure plotBlock;
+  plotList = fMsrInfo->GetMsrPlotList();
+  // find the proper plot block
+  // Here a small complication has to be handled: there are potentially multiple
+  // run blocks and the run might be present in various of these run blocks. In
+  // order to get a nice resolution on the theory the following procedure will be
+  // followed: the smallest x-interval found will be used to for the fXTheory resolution
+  // which is 1000 function points. The function will be calculated from the smallest
+  // xmin found up to the largest xmax found.
+  double xMin, xMax;
+  double xAbsMin, xAbsMax;
+  bool first = true;
+cout << endl << ">> plotList->size()=" << plotList->size();
+  for (unsigned int i=0; i<plotList->size(); i++) {
+    plotBlock = plotList->at(i);
+cout << endl << ">> plotBlock.fRuns.size()=" << plotBlock.fRuns.size() << endl;
+    for (unsigned int j=0; j<plotBlock.fRuns.size(); j++) {
+cout << endl << ">> j=" << j;
+cout << endl << ">> fRunNo=" << fRunNo;
+cout << endl << ">> plotBlock.fRuns[j].Re()=" << plotBlock.fRuns[j].Re();
+cout << endl;
+      if (fRunNo == plotBlock.fRuns[j].Re()-1) { // run found
+        if (first) {
+          first = false;
+          xMin = plotBlock.fTmin;
+          xMax = plotBlock.fTmax;
+          xAbsMin = xMin;
+          xAbsMax = xMax;
+cout << endl << ">> first: xMin=" << xMin << ", xMax=" << xMax << endl;
+        } else {
+          if (fabs(xMax-xMin) > fabs(plotBlock.fTmax-plotBlock.fTmin)) {
+            xMin = plotBlock.fTmin;
+            xMax = plotBlock.fTmax;
+          }
+          if (xMin < xAbsMin)
+            xAbsMin = xMin;
+          if (xMax > xAbsMax)
+            xAbsMax = xMax;
+cout << endl << ">> !first: xMin=" << xMin << ", xMax=" << xMax << endl;
+        }
+cout << endl << ">> xMin=" << xMin << ", xMax=" << xMax << endl;
+      }
+    }
+  }
+cout << endl << ">> after the xmin/xmax loop." << endl;
+
+  double xStep = (xMax-xMin)/1000.0;
+  double xx = xAbsMin;
+  do {
+    // fill x-vector
+    fData.fXTheory.push_back(xx);
+    // fill y-vector
+    fData.fTheory.push_back(fTheory->Func(xx, par, fFuncValues));
+    // calculate next xx
+    xx += xStep;
+  } while (xx < xAbsMax);
+
+  // clean up
+  par.clear();
+
   return success;
 }

src/include/PMusr.h

@@ -130,12 +130,13 @@ typedef struct {
   // data related info
   double fDataTimeStart;
   double fDataTimeStep;
-  PDoubleVector fX; // only used for non-muSR
+  PDoubleVector fX;       // only used for non-muSR
   PDoubleVector fValue;
   PDoubleVector fError;
   // theory related info
   double fTheoryTimeStart;
   double fTheoryTimeStep;
+  PDoubleVector fXTheory; // only used for non-muSR
   PDoubleVector fTheory;
 } PRunData;
 

src/include/PMusrCanvas.h

@@ -40,6 +40,7 @@
 #include <TLegend.h>
 #include <TPad.h>
 #include <TH1F.h>
+#include <TGraphErrors.h>
 
 #include "PMusr.h"
 #ifndef __MAKECINT__
@@ -79,6 +80,34 @@ typedef struct {
  */
 typedef vector<PMusrCanvasDataSet> PMusrCanvasDataList;
 
+//------------------------------------------------------------------------
+/**
+ * <p>
+ */
+typedef struct {
+  TGraphErrors *data;
+  TGraphErrors *dataFourierRe;
+  TGraphErrors *dataFourierIm;
+  TGraphErrors *dataFourierPwr;
+  TGraphErrors *dataFourierPhase;
+  TGraphErrors *theory;
+  TGraphErrors *theoryFourierRe;
+  TGraphErrors *theoryFourierIm;
+  TGraphErrors *theoryFourierPwr;
+  TGraphErrors *theoryFourierPhase;
+  TGraphErrors *diff;
+  TGraphErrors *diffFourierRe;
+  TGraphErrors *diffFourierIm;
+  TGraphErrors *diffFourierPwr;
+  TGraphErrors *diffFourierPhase;
+} PMusrCanvasNonMusrDataSet;
+
+//------------------------------------------------------------------------
+/**
+ * <p>
+ */
+typedef vector<PMusrCanvasNonMusrDataSet> PMusrCanvasNonMusrDataList;
+
 //--------------------------------------------------------------------------
 /**
  * <p>The preprocessor tag __MAKECINT__ is used to hide away from rootcint
@@ -127,6 +156,7 @@ class PMusrCanvas : public TObject, public TQObject
 #endif // __MAKECINT__
 
     PMusrCanvasDataList fData;
+    PMusrCanvasNonMusrDataList fNonMusrData;
 
     PIntVector fMarkerList;
     PIntVector fColorList;
@@ -134,8 +164,11 @@ class PMusrCanvas : public TObject, public TQObject
     virtual void CreateStyle();
     virtual void InitMusrCanvas(const char* title, Int_t wtopx, Int_t wtopy, Int_t ww, Int_t wh);
     virtual void InitDataSet(PMusrCanvasDataSet &dataSet);
+    virtual void InitDataSet(PMusrCanvasNonMusrDataSet &dataSet);
     virtual void CleanupDataSet(PMusrCanvasDataSet &dataSet);
+    virtual void CleanupDataSet(PMusrCanvasNonMusrDataSet &dataSet);
     virtual void HandleDataSet(unsigned int runNo, PRunData *data);
+    virtual void HandleNonMusrDataSet(unsigned int runNo, PRunData *data);
 
   ClassDef(PMusrCanvas, 1)
 };