some more work towards Nonlocal fitting

src/external/Nonlocal/Makefile

@@ -63,6 +63,7 @@ GLIBS         = $(ROOTGLIBS) -lXMLParser
 # some definitions: headers (used to generate *Dict* stuff), sources, objects,...
 OBJS =
 OBJS += PNL_StartupHandler.o PNL_StartupHandlerDict.o
+OBJS += PNL_RgeHandler.o
 OBJS += PNL_PippardFitter.o PNL_PippardFitterDict.o
 
 SHLIB = libPNL_PippardFitter.so

src/external/Nonlocal/PNL_PippardFitter.cpp

@@ -29,6 +29,13 @@
  *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
  ***************************************************************************/
 
+#include <cassert>
+
+#include <iostream>
+using namespace std;
+
+#include <TSAXParser.h>
+
 #include "PNL_PippardFitter.h"
 
 ClassImp(PNL_PippardFitter)
@@ -41,6 +48,39 @@ ClassImp(PNL_PippardFitter)
  */
 PNL_PippardFitter::PNL_PippardFitter()
 {
+  // read XML startup file
+  char startup_path_name[128];
+  TSAXParser *saxParser = new TSAXParser();
+  PNL_StartupHandler *fStartupHandler = new PNL_StartupHandler();
+  strcpy(startup_path_name, fStartupHandler->GetStartupFilePath().Data());
+  saxParser->ConnectToHandler("PNL_StartupHandler", fStartupHandler);
+  Int_t status = saxParser->ParseFile(startup_path_name);
+  // check for parse errors
+  if (status) { // error
+    cout << endl << "**WARNING** reading/parsing nonlocal_startup.xml.";
+    cout << endl;
+    // clean up
+    if (saxParser) {
+      delete saxParser;
+      saxParser = 0;
+    }
+    if (fStartupHandler) {
+      delete fStartupHandler;
+      fStartupHandler = 0;
+    }
+    assert(false);
+  }
+
+  // clean up
+  if (saxParser) {
+    delete saxParser;
+    saxParser = 0;
+  }
+
+  // load all the TRIM.SP rge-files
+  fRgeHandler = new PNL_RgeHandler(fStartupHandler->GetTrimSpDataPathList());
+  if (!fRgeHandler->IsValid())
+    assert(false);
 }
 
 //--------------------------------------------------------------------------
@@ -51,6 +91,10 @@ PNL_PippardFitter::PNL_PippardFitter()
  */
 PNL_PippardFitter::~PNL_PippardFitter()
 {
+  if (fStartupHandler) {
+    delete fStartupHandler;
+    fStartupHandler = 0;
+  }
 }
 
 //--------------------------------------------------------------------------
@@ -61,7 +105,8 @@ PNL_PippardFitter::~PNL_PippardFitter()
  */
 Double_t PNL_PippardFitter::operator()(Double_t t, const std::vector<Double_t> &param) const
 {
-//  assert(param.size() >= 4);
+  // expected parameters: energy, temp, thickness, meanFreePath, xi0, lambdaL
+  assert(param.size() != 6);
 
   return 0.0;
 }

src/external/Nonlocal/PNL_PippardFitter.h

@@ -33,6 +33,8 @@
 #define _PNL_PIPPARDFITTER_H_
 
 #include "PUserFcnBase.h"
+#include "PNL_StartupHandler.h"
+#include "PNL_RgeHandler.h"
 
 class PNL_PippardFitter : public PUserFcnBase
 {
@@ -42,6 +44,10 @@ class PNL_PippardFitter : public PUserFcnBase
 
     virtual Double_t operator()(Double_t t, const std::vector<Double_t> &param) const;
 
+  private:
+    PNL_StartupHandler *fStartupHandler;
+    PNL_RgeHandler *fRgeHandler;
+
   ClassDef(PNL_PippardFitter, 1)
 };
 

src/external/Nonlocal/PNL_RgeHandler.cpp

@@ -0,0 +1,208 @@
+/***************************************************************************
+
+  PNL_RgeHandler.cpp
+
+  Author: Andreas Suter
+  e-mail: andreas.suter@psi.ch
+
+  $Id$
+
+***************************************************************************/
+
+/***************************************************************************
+ *   Copyright (C) 2009 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 <iostream>
+#include <fstream>
+using namespace std;
+
+#include "PNL_RgeHandler.h"
+
+//--------------------------------------------------------------------------
+// Constructor
+//--------------------------------------------------------------------------
+/**
+ *
+ */
+PNL_RgeHandler::PNL_RgeHandler(const PStringVector &rgeDataPathList)
+{
+  fIsValid = false;
+
+  fIsValid = LoadRgeData(rgeDataPathList);
+}
+
+//--------------------------------------------------------------------------
+// Destructor
+//--------------------------------------------------------------------------
+/**
+ *
+ */
+PNL_RgeHandler::~PNL_RgeHandler()
+{
+  fRgeDataList.clear();
+}
+
+//--------------------------------------------------------------------------
+// GetRgeEnergyIndex
+//--------------------------------------------------------------------------
+/**
+ *
+ */
+Int_t PNL_RgeHandler::GetRgeEnergyIndex(const Double_t energy)
+{
+  Int_t idx = -1;
+
+  for (UInt_t i=0; i<fRgeDataList.size(); i++) {
+    if (energy == fRgeDataList[i].energy) {
+      idx = i;
+      break;
+    }
+  }
+
+  return idx;
+}
+
+//--------------------------------------------------------------------------
+// GetRgeValue
+//--------------------------------------------------------------------------
+/**
+ *
+ */
+Double_t PNL_RgeHandler::GetRgeValue(const Int_t index, const Double_t dist)
+{
+  Double_t rgeVal = 0.0;
+
+  UInt_t distIdx = (UInt_t)(dist/(fRgeDataList[index].stoppingDistance[1]-fRgeDataList[index].stoppingDistance[0]));
+
+  if ((distIdx >= fRgeDataList[index].stoppingDistance.size()) || (distIdx < 0)) {
+    rgeVal = 0.0;
+  } else {
+    rgeVal = fRgeDataList[index].stoppingAmplitude[distIdx] +
+             (fRgeDataList[index].stoppingAmplitude[distIdx] - fRgeDataList[index].stoppingAmplitude[distIdx+1]) *
+             (fRgeDataList[index].stoppingDistance[distIdx+1]-dist)/(fRgeDataList[index].stoppingDistance[distIdx+1]-fRgeDataList[index].stoppingDistance[distIdx]);
+  }
+
+  return rgeVal;
+}
+
+//--------------------------------------------------------------------------
+// GetRgeValue
+//--------------------------------------------------------------------------
+/**
+ *
+ */
+Double_t PNL_RgeHandler::GetRgeValue(const Double_t energy, const Double_t dist)
+{
+  // check if energy is present in rge data list
+  Int_t idx = -1;
+
+  for (UInt_t i=0; i<fRgeDataList.size(); i++) {
+    if (energy == fRgeDataList[i].energy) {
+      idx = i;
+      break;
+    }
+  }
+
+  // energy not found
+  if (idx == -1)
+    return -1.0;
+
+  return GetRgeValue(idx, dist);
+}
+
+//--------------------------------------------------------------------------
+// LoadRgeData
+//--------------------------------------------------------------------------
+/**
+ *
+ */
+Bool_t PNL_RgeHandler::LoadRgeData(const PStringVector &rgeDataPathList)
+{
+  ifstream fin;
+  PNL_RgeData data;
+  Int_t idx=0;
+  TString dataName;
+  char line[512];
+  int result;
+  double dist, val;
+
+  for (UInt_t i=0; i<rgeDataPathList.size(); i++) {
+    // open rge-file for reading
+    fin.open(rgeDataPathList[i].Data(), iostream::in);
+    if (!fin.is_open()) {
+      return false;
+    }
+
+    // extract energy from rgeDataPathList name
+    dataName = rgeDataPathList[i];
+    // remove rge extension
+    dataName.Remove(dataName.Length()-4, 4);
+    // get energy part
+    dataName.Replace(0, dataName.Length()-3, "");
+    if (!dataName.IsDigit()) {
+      fin.close();
+      return false;
+    }
+    data.energy = dataName.Atof()/10.0;
+
+    // read msr-file
+    idx = 0;
+    while (!fin.eof()) {
+      // read a line
+      fin.getline(line, sizeof(line));
+      idx++;
+
+      // ignore first line
+      if (idx == 1)
+        continue;
+
+      // get values
+      result = sscanf(line, "%lf %lf", &dist, &val);
+      // check if data are valid
+      if (result != 2) {
+        fin.close();
+        return false;
+      }
+
+      // feed fRgeDataList
+      data.stoppingDistance.push_back(dist/10.0);
+      data.stoppingAmplitude.push_back(val);
+    }
+
+    // normalize stopping distribution
+    Double_t norm = 0.0;
+    for (UInt_t j=0; j<data.stoppingAmplitude.size(); j++)
+      norm += data.stoppingAmplitude[i];
+    norm *= (data.stoppingDistance[1] - data.stoppingDistance[0]);
+    for (UInt_t j=0; j<data.stoppingAmplitude.size(); j++)
+      data.stoppingAmplitude[i] /= norm;
+
+    // keep data
+    fRgeDataList.push_back(data);
+
+    // clean up
+    data.stoppingAmplitude.clear();
+    data.stoppingDistance.clear();
+
+    fin.close();
+  }
+
+  return true;
+}

src/external/Nonlocal/PNL_RgeHandler.h

@@ -0,0 +1,55 @@
+/***************************************************************************
+
+  PNL_RgeHandler.h
+
+  Author: Andreas Suter
+  e-mail: andreas.suter@psi.ch
+
+  $Id$
+
+***************************************************************************/
+
+/***************************************************************************
+ *   Copyright (C) 2009 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 _PNL_RGEHANDLER_H_
+#define _PNL_RGEHANDLER_H_
+
+#include "PNonlocal.h"
+
+class PNL_RgeHandler
+{
+  public:
+    PNL_RgeHandler(const PStringVector &rgeDataPathList);
+    virtual ~PNL_RgeHandler();
+
+    virtual Bool_t IsValid() { return fIsValid; }
+    virtual Int_t GetRgeEnergyIndex(const Double_t energy);
+    virtual Double_t GetRgeValue(const Double_t energy, const Double_t dist);
+    virtual Double_t GetRgeValue(const Int_t index, const Double_t dist);
+
+  private:
+    Bool_t fIsValid;
+    PNL_RgeDataList fRgeDataList;
+
+    virtual Bool_t LoadRgeData(const PStringVector &rgeDataPathList);
+};
+
+#endif // _PNL_RGEHANDLER_H_

src/external/Nonlocal/PNL_StartupHandler.h

@@ -55,6 +55,7 @@ class PNL_StartupHandler : public TObject
     virtual void OnFatalError(const char*); // SLOT
     virtual void OnCdataBlock(const char*, Int_t); // SLOT
 
+    virtual TString GetStartupFilePath() { return fStartupFilePath; }
     virtual const Int_t GetFourierPoints() const { return fFourierPoints; }
     virtual const PStringVector GetTrimSpDataPathList() const { return fTrimSpDataPathList; }
 

src/external/Nonlocal/PNonlocal.h

@@ -42,5 +42,26 @@
  */
 typedef std::vector<TString> PStringVector;
 
+//-------------------------------------------------------------
+/**
+ * <p> typedef to make to code more readable.
+ */
+typedef std::vector<Double_t> PDoubleVector;
+
+//-------------------------------------------------------------
+/**
+ * <p>
+ */
+typedef struct {
+  Double_t energy;
+  PDoubleVector stoppingDistance;
+  PDoubleVector stoppingAmplitude;
+} PNL_RgeData;
+
+//-------------------------------------------------------------
+/**
+ * <p>
+ */
+typedef std::vector<PNL_RgeData> PNL_RgeDataList;
 
 #endif // _PNONLOCAL_H_