diff --git a/src/tests/nonlocal/Makefile b/src/tests/nonlocal/Makefile
index 71101ed1..fced9d42 100644
--- a/src/tests/nonlocal/Makefile
+++ b/src/tests/nonlocal/Makefile
@@ -87,6 +87,9 @@ nonlocal$(SUFFIX): $(OBJS)
 #
 clean:; 	@rm -f $(OBJS)
 	@echo "---> removing $(OBJS)"
+	
+install:
+	cp nonlocal $(HOME)/bin
 
 #
 $(OBJS): %.o: %.cpp
diff --git a/src/tests/nonlocal/PPippard.cpp b/src/tests/nonlocal/PPippard.cpp
index 79edc591..def47843 100644
--- a/src/tests/nonlocal/PPippard.cpp
+++ b/src/tests/nonlocal/PPippard.cpp
@@ -3,6 +3,7 @@
 // $Id$
 // -----------------------------------------------------------------------
 
+#include <cassert>
 #include <cstdio>
 #include <cmath>
 
@@ -25,9 +26,6 @@ PPippard::PPippard(const PippardParams &params) : fParams(params)
   fFieldq = 0;
   fFieldB = 0;
 
-  fSecondDerivativeMatrix = 0;
-  fKernelMatrix = 0;
-  fBoundaryCondition = 0;
   fFieldDiffuse = 0;
 
   f_dx = 0.02;
@@ -48,29 +46,48 @@ PPippard::~PPippard()
     fftw_free(fFieldq);
     fFieldq = 0;
   }
-
   if (fFieldB) {
     fftw_free(fFieldq);
     fFieldB = 0;
   }
-  if (fSecondDerivativeMatrix) {
-    delete fSecondDerivativeMatrix;
-    fSecondDerivativeMatrix = 0;
-  }
-  if (fKernelMatrix) {
-    delete fKernelMatrix;
-    fKernelMatrix = 0;
-  }
-  if (fBoundaryCondition) {
-    delete fBoundaryCondition;
-    fBoundaryCondition = 0;
-  }
   if (fFieldDiffuse) {
     delete fFieldDiffuse;
     fFieldDiffuse = 0;
   }
 }
 
+//-----------------------------------------------------------------------------------------------------------
+/**
+ *
+ */
+void PPippard::DumpParams()
+{
+  cout << endl << ">> Parameters: ";
+  cout << endl << ">> reduced temperature: " << fParams.t;
+  cout << endl << ">> lambdaL            : " << fParams.lambdaL << " (nm)";
+  cout << endl << ">> xi0                : " << fParams.xi0 << " (nm)";
+  cout << endl << ">> mean free path     : " << fParams.meanFreePath << " (nm)";
+  cout << endl << ">> film thickness     : " << fParams.filmThickness << " (nm)";
+  if (fParams.specular)
+    cout << endl << ">> specular scattering boundary conditions";
+  else
+    cout << endl << ">> diffuse scattering boundary conditions";
+  cout << endl << ">> Bext               : " << fParams.b_ext << " (G)";
+  cout << endl << ">> dead layer         : " << fParams.deadLayer << " (nm)";
+  for (UInt_t i=0; i<fParams.rgeFileName.size(); i++)
+    cout << endl << ">> rge-file-name      : " << fParams.rgeFileName[i].Data();
+  if (fParams.outputFileName.Length() > 0)
+    cout << endl << ">> output file name   : " << fParams.outputFileName.Data();
+  else
+    cout << endl << ">> output file name   : n/a";
+  if (fParams.outputFileNameBmean.Length() > 0)
+    cout << endl << ">> output file name Bmean : " << fParams.outputFileNameBmean.Data();
+  else
+    cout << endl << ">> output file name Bmean : n/a";
+  cout << endl << ">> input file name    : " << fParams.inputFileName.Data();
+  cout << endl << endl;
+}
+
 //-----------------------------------------------------------------------------------------------------------
 /**
  *
@@ -248,8 +265,6 @@ void PPippard::CalculateFieldSpecular()
     }
   }
 
-cout << endl << "debug> fShift = " << fShift;
-
   for (Int_t i=0; i<PippardFourierPoints; i++) {
     fFieldB[i][1] /= norm;
   }
@@ -282,69 +297,64 @@ cout << endl << "debug> fShift = " << fShift;
  */
 void PPippard::CalculateFieldDiffuse()
 {
+  MatrixXd fSecondDerivativeMatrix(PippardDiffusePoints+1, PippardDiffusePoints+1); // 2nd derivative matrix
+  MatrixXd fKernelMatrix(PippardDiffusePoints+1, PippardDiffusePoints+1);           // kernel matrix
+  VectorXd fBoundaryCondition(PippardDiffusePoints+1);                              // boundary condition vector
+
   f_dz = 5.0/XiP_T(fParams.t);
 
   Double_t invL = 1/f_dz;
   Double_t ampl = 1.0/pow(f_dz,2.0)/(3.0/4.0*pow(XiP_T(fParams.t),3.0)/(fParams.xi0*pow(LambdaL_T(fParams.t),2.0)));
 
-cout << endl << ">> 1/alpha = " << 1.0/(3.0/4.0*pow(XiP_T(fParams.t),3.0)/(fParams.xi0*pow(LambdaL_T(fParams.t),2.0)));
-cout << endl << ">> 1/l^2   = " << 1.0/pow(f_dz,2.0);
-cout << endl << ">> ampl    = " << ampl << endl;
+  cout << endl << ">> f_dz    = " << f_dz << ", invL = " << invL;
+  cout << endl << ">> 1/alpha = " << 1.0/(3.0/4.0*pow(XiP_T(fParams.t),3.0)/(fParams.xi0*pow(LambdaL_T(fParams.t),2.0)));
+  cout << endl << ">> 1/l^2   = " << 1.0/pow(f_dz,2.0);
+  cout << endl << ">> ampl    = " << ampl << endl;
 
   // 2nd derivative matrix
-  if (fSecondDerivativeMatrix == 0) { // first time call, hence generate the 2nd derivative matrix
-    fSecondDerivativeMatrix = new MatrixXd(PippardDiffusePoints+1, PippardDiffusePoints+1);
-    fSecondDerivativeMatrix->setZero(PippardDiffusePoints+1, PippardDiffusePoints+1);
-    for (Int_t i=1; i<PippardDiffusePoints; i++) {
-      (*fSecondDerivativeMatrix)(i,i-1) = ampl;
-      (*fSecondDerivativeMatrix)(i,i) = -2.0*ampl;
-      (*fSecondDerivativeMatrix)(i,i+1) = ampl;
-    }
+  fSecondDerivativeMatrix.setZero();
+  for (Int_t i=1; i<PippardDiffusePoints; i++) {
+    fSecondDerivativeMatrix(i,i-1) = ampl;
+    fSecondDerivativeMatrix(i,i) = -2.0*ampl;
+    fSecondDerivativeMatrix(i,i+1) = ampl;
   }
 
-//cout << endl << "fSecondDerivativeMatrix = \n" << *fSecondDerivativeMatrix << endl;
-
   // kernel matrix
-  if (fKernelMatrix == 0) { // first time call, hence generate the kernel matrix
-    fKernelMatrix = new MatrixXd(PippardDiffusePoints+1, PippardDiffusePoints+1);
-    fKernelMatrix->setZero(PippardDiffusePoints+1, PippardDiffusePoints+1);
-    // 1st line (dealing with boundary conditions)
-    (*fKernelMatrix)(0,0) = -1.5*invL;
-    (*fKernelMatrix)(0,1) =  2.0*invL;
-    (*fKernelMatrix)(0,2) = -0.5*invL;
-    // Nth line (dealing with boundary conditions)
-    (*fKernelMatrix)(PippardDiffusePoints,PippardDiffusePoints-2) =  0.5*invL;
-    (*fKernelMatrix)(PippardDiffusePoints,PippardDiffusePoints-1) = -2.0*invL;
-    (*fKernelMatrix)(PippardDiffusePoints,PippardDiffusePoints)   =  1.5*invL;
-    // the real kernel
-    for (Int_t i=1; i<PippardDiffusePoints; i++) {
-      (*fKernelMatrix)(i,0)=Calc_a(i);
-      (*fKernelMatrix)(i,PippardDiffusePoints)=Calc_b(i);
-      for (Int_t j=1; j<PippardDiffusePoints; j++) {
-        (*fKernelMatrix)(i,j) = Calc_c(i,j);
-      }
+  fKernelMatrix.setZero();
+
+  // 1st line (dealing with boundary conditions)
+  fKernelMatrix(0,0) = -1.5*invL;
+  fKernelMatrix(0,1) =  2.0*invL;
+  fKernelMatrix(0,2) = -0.5*invL;
+
+  // Nth line (dealing with boundary conditions)
+  fKernelMatrix(PippardDiffusePoints,PippardDiffusePoints-2) =  0.5*invL;
+  fKernelMatrix(PippardDiffusePoints,PippardDiffusePoints-1) = -2.0*invL;
+  fKernelMatrix(PippardDiffusePoints,PippardDiffusePoints)   =  1.5*invL;
+
+  // the real kernel
+  for (Int_t i=1; i<PippardDiffusePoints; i++) {
+    fKernelMatrix(i,0)=Calc_a(i);
+    fKernelMatrix(i,PippardDiffusePoints)=Calc_b(i);
+    for (Int_t j=1; j<PippardDiffusePoints; j++) {
+      fKernelMatrix(i,j) = Calc_c(i,j);
     }
   }
 
-//cout << endl << "fKernelMatrix = \n" << *fKernelMatrix << endl;
-
   // boundary condition vector
-  if (fBoundaryCondition == 0) {
-    fBoundaryCondition = new VectorXd(PippardDiffusePoints+1);
-    fBoundaryCondition->setZero(PippardDiffusePoints+1);
-    (*fBoundaryCondition)(0) = 1.0;
-  }
+  fBoundaryCondition.setZero();
+  fBoundaryCondition(0) = 1.0;
 
-//cout << endl << "fBoundaryCondition = " << *fBoundaryCondition << endl;
-
-  if (fFieldDiffuse == 0) {
-    fFieldDiffuse = new VectorXd(PippardDiffusePoints+1);
-    fFieldDiffuse->setZero(PippardDiffusePoints+1);
+  if (fFieldDiffuse != 0) {
+    delete fFieldDiffuse;
+    fFieldDiffuse = 0;
   }
+  fFieldDiffuse = new VectorXd(PippardDiffusePoints+1);
+  fFieldDiffuse->setZero();
 
   // solve equation
-  *fSecondDerivativeMatrix = (*fSecondDerivativeMatrix)-(*fKernelMatrix);
-  fSecondDerivativeMatrix->lu().solve(*fBoundaryCondition, fFieldDiffuse);
+  fSecondDerivativeMatrix = fSecondDerivativeMatrix-fKernelMatrix;
+  *fFieldDiffuse = fSecondDerivativeMatrix.colPivHouseholderQr().solve(fBoundaryCondition);
 
   // normalize field
   Double_t norm = 0.0;
@@ -352,6 +362,8 @@ cout << endl << ">> ampl    = " << ampl << endl;
     if (norm < (*fFieldDiffuse)(i))
       norm = (*fFieldDiffuse)(i);
 
+  assert(norm != 0);
+
   for (Int_t i=0; i<PippardDiffusePoints+1; i++)
     (*fFieldDiffuse)(i) /= norm;
 
@@ -392,12 +404,15 @@ void PPippard::SaveField()
     fprintf(fp, "%%   Boundary Conditions: Diffuse\n");
   fprintf(fp, "%%   Bext                    = %lf (G)\n", fParams.b_ext);
   fprintf(fp, "%%   deadLayer               = %lf (nm)\n", fParams.deadLayer);
-  if (fParams.rgeFileName.Length() > 0)
-    fprintf(fp, "%%   rge file name : %s\n", fParams.rgeFileName.Data());
-  if (fParams.meanB != 0.0) {
-    fprintf(fp, "%%   Mean Distance           = %lf (nm)\n", fParams.meanX);
-    fprintf(fp, "%%   Mean Field/Bext         = %lf, Mean Field = %lf (G)\n", fParams.meanB, fParams.meanB * fParams.b_ext);
-    fprintf(fp, "%%   Var Field/Bext^2        = %lf, Var Field  = %lf (G^2)\n", fParams.varB, fParams.varB * fParams.b_ext * fParams.b_ext);
+  for (UInt_t i=0; i<fParams.rgeFileName.size(); i++) {
+    fprintf(fp, "%%   rge file name : %s\n", fParams.rgeFileName[i].Data());
+  }
+  for (UInt_t i=0; i<fParams.energy.size(); i++) {
+    fprintf(fp, "%%------------------------------------------\n");
+    fprintf(fp, "%%   Implantation Energy     = %lf (eV)\n", fParams.energy[i]);
+    fprintf(fp, "%%   Mean Distance           = %lf (nm)\n", fParams.meanX[i]);
+    fprintf(fp, "%%   Mean Field/Bext         = %lf, Mean Field = %lf (G)\n", fParams.meanB[i], fParams.meanB[i] * fParams.b_ext);
+    fprintf(fp, "%%   Var Field/Bext^2        = %lf, Var Field  = %lf (G^2)\n", fParams.varB[i], fParams.varB[i] * fParams.b_ext * fParams.b_ext);
   }
   fprintf(fp, "%%\n");
 
@@ -417,6 +432,38 @@ void PPippard::SaveField()
   fclose(fp);
 }
 
+//-----------------------------------------------------------------------------------------------------------
+/**
+ *
+ */
+void PPippard::SaveBmean()
+{
+  FILE *fp;
+
+  fp = fopen(fParams.outputFileNameBmean.Data(), "w");
+  if (fp == NULL) {
+    cout << endl << "Coudln't open " << fParams.outputFileNameBmean.Data() << " for writting, sorry ...";
+    cout << endl << endl;
+    return;
+  }
+
+  // write header
+  fprintf(fp, "%% Header -----------------------\n");
+  fprintf(fp, "%% generated from input file: %s\n", fParams.inputFileName.Data());
+  fprintf(fp, "%%\n");
+  fprintf(fp, "%% Data -------------------------\n");
+  fprintf(fp, "%% energy (eV), meanZ (nm), meanB/Bext, meanB (G), varB/Bext^2, varB (G^2)\n");
+
+  // write data
+  for (UInt_t i=0; i<fParams.energy.size(); i++) {
+    fprintf(fp, "%lf, %lf, %lf, %lf, %lf, %lf\n", fParams.energy[i], fParams.meanX[i],
+            fParams.meanB[i], fParams.b_ext*fParams.meanB[i],
+            fParams.varB[i], (fParams.b_ext*fParams.b_ext)*fParams.varB[i]);
+  }
+
+  fclose(fp);
+}
+
 //-----------------------------------------------------------------------------------------------------------
 /**
  *
diff --git a/src/tests/nonlocal/PPippard.h b/src/tests/nonlocal/PPippard.h
index 6903acec..cc0ae556 100644
--- a/src/tests/nonlocal/PPippard.h
+++ b/src/tests/nonlocal/PPippard.h
@@ -6,13 +6,15 @@
 #ifndef _PPIPPARD_H_
 #define _PPIPPARD_H_
 
+#include <vector>
+using namespace std;
+
 #include <fftw3.h>
 
 #include <Rtypes.h>
 #include <TString.h>
 
-#include <Eigen/Core>
-#include <Eigen/LU>
+#include <Eigen/Dense>
 using namespace Eigen;
 
 const Int_t PippardFourierPoints = 200000;
@@ -29,11 +31,14 @@ typedef struct {
   Double_t specularIntegral; // int_x=0^infty B(x) dx / Bext
   Double_t b_ext;         // external field in (G)
   Double_t deadLayer;     // dead layer in (nm)
-  TString rgeFileName;
+  TString inputFileName;
+  vector<TString> rgeFileName;
   TString outputFileName;
-  Double_t meanB;         // int_x=0^infty B(x) n(x) dx / int_x=0^infty n(x) dx / Bext
-  Double_t varB;          // int_x=0^infty (B(x)-<B>)^2 n(x) dx / int_x=0^infty n(x) dx / Bext^2
-  Double_t meanX;         // int_x=0^infty x n(x) dx / int_x=0^infty n(x) dx
+  TString outputFileNameBmean;
+  vector<Double_t> energy;
+  vector<Double_t> meanB;         // int_x=0^infty B(x) n(x) dx / int_x=0^infty n(x) dx / Bext
+  vector<Double_t> varB;          // int_x=0^infty (B(x)-<B>)^2 n(x) dx / int_x=0^infty n(x) dx / Bext^2
+  vector<Double_t> meanX;         // int_x=0^infty x n(x) dx / int_x=0^infty n(x) dx
 } PippardParams;
 
 class PPippard
@@ -42,15 +47,18 @@ class PPippard
     PPippard(const PippardParams &params);
     virtual ~PPippard();
 
+    virtual void DumpParams();
+
     virtual void SetTemp(Double_t temp) { fParams.t = temp; }
     virtual void SetLambdaL(Double_t lambdaL) { fParams.lambdaL = lambdaL; }
     virtual void SetXi0(Double_t xi0) { fParams.xi0 = xi0; }
     virtual void SetMeanFreePath(Double_t meanFreePath) { fParams.meanFreePath = meanFreePath; }
     virtual void SetFilmThickness(Double_t thickness) { fParams.filmThickness = thickness; }
     virtual void SetSpecular(Bool_t specular) { fParams.specular = specular; }
-    virtual void SetMeanX(Double_t meanX) { fParams.meanX = meanX; }
-    virtual void SetMeanB(Double_t meanB) { fParams.meanB = meanB; }
-    virtual void SetVarB(Double_t BB) { fParams.varB = BB; }
+    virtual void SetEnergy(Double_t energy) { fParams.energy.push_back(energy); }
+    virtual void SetMeanX(Double_t meanX) { fParams.meanX.push_back(meanX); }
+    virtual void SetMeanB(Double_t meanB) { fParams.meanB.push_back(meanB); }
+    virtual void SetVarB(Double_t BB) { fParams.varB.push_back(BB); }
 
     virtual void CalculateField();
 
@@ -58,6 +66,7 @@ class PPippard
     virtual Double_t GetMagneticField(const Double_t x) const;
 
     virtual void SaveField();
+    virtual void SaveBmean();
 
   private:
     PippardParams fParams;
@@ -72,9 +81,6 @@ class PPippard
 
     Int_t fShift; // shift needed to pick up fFieldB at the maximum for B->0
 
-    MatrixXd *fSecondDerivativeMatrix; // 2nd derivative matrix
-    MatrixXd *fKernelMatrix;           // kernel matrix
-    VectorXd *fBoundaryCondition;      // boundary condition vector
     VectorXd *fFieldDiffuse;           // resulting B(z)/B(0) field
 
     virtual Double_t DeltaBCS(const Double_t t) const;
diff --git a/src/tests/nonlocal/nonlocal.cpp b/src/tests/nonlocal/nonlocal.cpp
index 9851d2d7..aabeefec 100644
--- a/src/tests/nonlocal/nonlocal.cpp
+++ b/src/tests/nonlocal/nonlocal.cpp
@@ -10,6 +10,12 @@ using namespace std;
 
 #include "PPippard.h"
 
+typedef struct {
+  Double_t energy;
+  vector<Double_t> x;
+  vector<Double_t> n;
+} Pnofx;
+
 //-----------------------------------------------------------------------------------------
 /**
  *
@@ -29,8 +35,9 @@ void syntax()
   cout << endl << "%   specular         : 1=specular, 0=diffuse";
   cout << endl << "%   [Bext]           : external magnetic field (G), i.e. this is optional";
   cout << endl << "%   [deadLayer]      : dead layer (nm), i.e. this is optional";
-  cout << endl << "%   [rgeFileName]    : rge input file name, i.e. this is optional";
+  cout << endl << "%   [rgeFileName]    : rge input file name, i.e. this is optional. Multiple rge-files possible";
   cout << endl << "%   [outputFileName] : output file name, i.e. this is optional";
+  cout << endl << "%   [outputFileNameBmean] : output file name for <B> data";
   cout << endl;
   cout << endl << "reducedTemp = 0.8287";
   cout << endl << "lambdaL = 30.0";
@@ -40,8 +47,10 @@ void syntax()
   cout << endl << "specular = 1";
   cout << endl << "Bext = 47.14";
   cout << endl << "deadLayer = 0.4";
-  cout << endl << "rgeFileName = /afs/psi.ch/project/nemu/analysis/2009/Nonlocal/trimsp/InSne060.rge";
+  cout << endl << "rgeFileName = /afs/psi.ch/project/nemu/analysis/2009/Nonlocal/trimsp/InSn_E6000.rge";
+  cout << endl << "rgeFileName = /afs/psi.ch/project/nemu/analysis/2009/Nonlocal/trimsp/InSn_E14100.rge";
   cout << endl << "outputFileName = In_37_T2P83.dat";
+  cout << endl << "outputFileNameBmean = In_37_T2P83_Bmean.dat";
   cout << endl << endl;
   cout << endl << "lines starting with a '%' or a '#' are comment lines.";
   cout << endl << "empty lines are ignored.";
@@ -183,7 +192,7 @@ int readInputFile(char *fln, PippardParams &param)
         param.valid = false;
         break;
       }
-      param.rgeFileName = TString(str);
+      param.rgeFileName.push_back(TString(str));
     }
 
     if (tstr.BeginsWith("outputFileName =")) {
@@ -197,6 +206,18 @@ int readInputFile(char *fln, PippardParams &param)
       }
       param.outputFileName = TString(str);
     }
+
+    if (tstr.BeginsWith("outputFileNameBmean =")) {
+      result = sscanf(line, "outputFileNameBmean = %s", str);
+      if (result != 1) {
+        cout << endl << "**ERROR** while trying to extracted the Bmean output file name";
+        cout << endl << "line = " << line;
+        cout << endl;
+        param.valid = false;
+        break;
+      }
+      param.outputFileNameBmean = TString(str);
+    }
   }
 
   fclose(fp);
@@ -208,10 +229,28 @@ int readInputFile(char *fln, PippardParams &param)
 /**
  *
  */
-int readRgeFile(const char *fln, vector<Double_t> &x, vector<Double_t> &n)
+int readRgeFile(const char *fln, vector<Double_t> &x, vector<Double_t> &n, Double_t &energy)
 {
   FILE *fp;
 
+  x.clear();
+  n.clear();
+
+  TString str = TString(fln);
+  str.Remove(0, str.Index("_E")+2);
+  str.Remove(str.Index(".rge"));
+  if (str.Length() > 0) {
+    if (str.IsFloat()) {
+      energy = str.Atof();
+    } else {
+      cerr << endl << "**ERROR** couldn't extract energy from the file name '" << fln << "', will quit." << endl << endl;
+      return 0;
+    }
+  } else {
+    cerr << endl << "**ERROR** couldn't extract energy from the file name '" << fln << "', will quit." << endl << endl;
+    return 0;
+  }
+
   fp = fopen(fln, "r");
   if (fp == NULL) {
     cout << endl << "Couldn't read rge file " << fln;
@@ -253,15 +292,6 @@ int readRgeFile(const char *fln, vector<Double_t> &x, vector<Double_t> &n)
     n[i] /= norm * (x[1]-x[0]);
   }
 
-/*
-Double_t xx=0.0;
-for (unsigned int i=0; i<x.size(); i++) {
-  xx += n[i];
-}
-xx *= (x[1]-x[0]);
-cout << "debug> xx=" << xx << endl;
-*/
-
   fclose(fp);
 
   return 1;
@@ -291,10 +321,12 @@ int main(int argc, char *argv[])
   params.specularIntegral = -1.0;
   params.b_ext = -1.0;
   params.deadLayer = -1.0;
-  params.rgeFileName = "";
+  params.rgeFileName.clear();
+  params.inputFileName = TString(argv[1]);
   params.outputFileName = "";
-  params.meanX = 0.0;
-  params.meanB = 0.0;
+  params.outputFileNameBmean = "";
+  params.meanX.clear();
+  params.meanB.clear();
 
   if (!readInputFile(argv[1], params)) {
     cout << endl << "**ERROR** Couldn't open the input file " << argv[1] << ". Will quit.";
@@ -335,12 +367,34 @@ int main(int argc, char *argv[])
   }
 
   int result;
+  Pnofx rgeSet;
+  vector<Pnofx> nOfx;
+  double energy;
   vector<Double_t> x, n;
-  if (params.rgeFileName.Length() > 0) {
-    result = readRgeFile(params.rgeFileName.Data(), x, n);
-    if (result == 0)
-      params.rgeFileName = "";
+  for (UInt_t i=0; i<params.rgeFileName.size(); i++) {
+    x.clear();
+    n.clear();
+    energy = 0.0;
+    rgeSet.energy = 0.0;
+    rgeSet.x.clear();
+    rgeSet.n.clear();
+
+    result = readRgeFile(params.rgeFileName[i].Data(), x, n, energy);
+    if (result == 0) {
+      cerr << endl << "**ERROR** Couldn't read rge-file " << params.rgeFileName[i].Data() << ", will quit." << endl << endl;
+      return 0;
+    } else {
+      cout << endl << ">> read rge-file '" << params.rgeFileName[i].Data() << "'";
+    }
+
+    cout << endl << ">> energy        = " << energy << " eV";
+    rgeSet.energy = energy;
+    rgeSet.x = x;
+    rgeSet.n = n;
+    nOfx.push_back(rgeSet);
   }
+  cout << endl << ">> found " << nOfx.size() << " rge-files.";
+
 
   cout << endl << ">> temp          = " << params.t;
   cout << endl << ">> lambdaL       = " << params.lambdaL;
@@ -348,60 +402,75 @@ int main(int argc, char *argv[])
   cout << endl << ">> meanFreePath  = " << params.meanFreePath;
   cout << endl << ">> filmThickness = " << params.filmThickness;
   if (params.specular)
-    cout << endl << ">> specular     = true";
+    cout << endl << ">> specular      = true";
   else
-    cout << endl << ">> specular     = false";
+    cout << endl << ">> specular      = false";
   cout << endl;
 
   PPippard *pippard = new PPippard(params);
 
+//  pippard->DumpParams();
+
   pippard->CalculateField();
 
   // check if it is necessary to calculate the <B(x)> and <[B(x)-<B(x)>]^2>
-  if (x.size() > 0) {
+  if (nOfx.size() > 0) {
     if ((params.b_ext == -1.0) || (params.deadLayer == -1.0)) {
       cout << endl << "**ERROR** Bext or deadLayer missing :-(" << endl;
       return 0;
     }
 
     Double_t meanX = 0.0;
-    for (unsigned int i=0; i<x.size(); i++) {
-      meanX += x[i]*n[i];
-    }
-    meanX *= (x[1]-x[0]);
-
     Double_t meanB = 0.0;
     Double_t varB = 0.0;
     Double_t BB = 0.0;
-    for (unsigned int i=0; i<x.size()-1; i++) {
-      if (x[i] <= params.deadLayer) {
-        meanB += 1.0 * n[i];
-      } else if (x[i] > params.deadLayer+params.filmThickness) {
-        meanB += 1.0 * n[i];
-      } else {
-        BB = pippard->GetMagneticField(x[i]-params.deadLayer);
-        meanB += BB * n[i];
+    for (UInt_t i=0; i<nOfx.size(); i++) {
+      meanX = 0.0;
+      for (unsigned int j=0; j<nOfx[i].x.size(); j++) {
+        meanX += nOfx[i].x[j]*nOfx[i].n[j];
       }
-    }
-    meanB *= (x[1]-x[0]);
-    for (unsigned int i=0; i<x.size()-1; i++) {
-      if (x[i] > params.deadLayer) {
-        BB = pippard->GetMagneticField(x[i]-params.deadLayer);
-        varB += (BB-meanB) * (BB-meanB) * n[i];
-      }
-    }
-    varB *= (x[1]-x[0]);
+      meanX *= (nOfx[i].x[1]-nOfx[i].x[0]);
 
-    cout << endl << ">> mean x = " << meanX << " (nm), mean field = " << params.b_ext * meanB << " (G)";
-    cout << " <(B-<B>)^2> = " << varB * params.b_ext * params.b_ext << " (G^2)";
-    pippard->SetMeanX(meanX);
-    pippard->SetVarB(varB);
-    pippard->SetMeanB(meanB);
+      meanB = 0.0;
+      varB = 0.0;
+      BB = 0.0;
+      for (unsigned int j=0; j<nOfx[i].x.size()-1; j++) {
+        if (nOfx[i].x[j] <= params.deadLayer) {
+          meanB += 1.0 * nOfx[i].n[j];
+        } else if (nOfx[i].x[j] > params.deadLayer+params.filmThickness) {
+          meanB += 1.0 * nOfx[i].n[j];
+        } else {
+          BB = pippard->GetMagneticField(nOfx[i].x[j]-params.deadLayer);
+          meanB += BB * nOfx[i].n[j];
+        }
+      }
+      meanB *= (nOfx[i].x[1]-nOfx[i].x[0]);
+      for (unsigned int j=0; j<nOfx[i].x.size()-1; j++) {
+        if (nOfx[i].x[j] > params.deadLayer) {
+          BB = pippard->GetMagneticField(nOfx[i].x[j]-params.deadLayer);
+          varB += (BB-meanB) * (BB-meanB) * nOfx[i].n[j];
+        }
+      }
+      varB *= (nOfx[i].x[1]-nOfx[i].x[0]);
+
+      cout << endl << "----------------------------------------";
+      cout << endl << ">> energy = " << nOfx[i].energy << "eV";
+      cout << endl << ">> mean x = " << meanX << " (nm), mean field = " << params.b_ext * meanB << " (G),";
+      cout << " <(B-<B>)^2> = " << varB * params.b_ext * params.b_ext << " (G^2)";
+
+      pippard->SetEnergy(nOfx[i].energy);
+      pippard->SetMeanX(meanX);
+      pippard->SetVarB(varB);
+      pippard->SetMeanB(meanB);
+    }
   }
 
   if (params.outputFileName.Length() > 0)
     pippard->SaveField();
 
+  if (params.outputFileNameBmean.Length() > 0)
+    pippard->SaveBmean();
+
   cout << endl << endl;
 
   if (pippard) {