diff --git a/BP_analysis.cpp b/BP_analysis.cpp
index a4d8159..753e93d 100644
--- a/BP_analysis.cpp
+++ b/BP_analysis.cpp
@@ -1,11 +1,11 @@
 // to analyse the backplane pulsing data per module
 // changes by VH 210906: to eliminate hardcoded absolute paths, uses location of the analysis root files as additional input argument (accordingly changed in filename_creator.sh)
 
-#include "../sls_detector_calibration/jungfrauCommonHeader.h"
-#include "../sls_detector_calibration/jungfrauCommonFunctions.h"
+#include "sls_detector_calibration/jungfrauCommonHeader.h"
+#include "sls_detector_calibration/jungfrauCommonFunctions.h"
 
-#include "../sls_detector_calibration/jungfrauFile.C"
-#include "../sls_detector_calibration/jungfrauPedestal.C"
+#include "sls_detector_calibration/jungfrauFile.C"
+#include "sls_detector_calibration/jungfrauPedestal.C"
 
 #include "TGraphErrors.h"
 #include "TF1.h"
diff --git a/CS_analysis.cpp b/CS_analysis.cpp
index fde431d..6f03497 100644
--- a/CS_analysis.cpp
+++ b/CS_analysis.cpp
@@ -2,12 +2,12 @@
 // changes by VH 210906: to eliminate hardcoded absolute paths, uses location of the analysis root files as additional input argument (accordingly changed in filename_creator.sh)
 
 
-#include "../sls_detector_calibration/jungfrauCommonHeader.h"
-#include "../sls_detector_calibration/jungfrauCommonFunctions.h"
+#include "sls_detector_calibration/jungfrauCommonHeader.h"
+#include "sls_detector_calibration/jungfrauCommonFunctions.h"
 
-#include "../sls_detector_calibration/jungfrauFile.C"
-#include "../sls_detector_calibration/jungfrauPedestal.C"
-#include "../sls_detector_calibration/jungfrauPixelMask.C"
+#include "sls_detector_calibration/jungfrauFile.C"
+#include "sls_detector_calibration/jungfrauPedestal.C"
+#include "sls_detector_calibration/jungfrauPixelMask.C"
 
 #include "TGraph.h"
 #include "TGraphErrors.h"
diff --git a/CuFluo_analysis.cpp b/CuFluo_analysis.cpp
index 2b72591..a3cce9a 100644
--- a/CuFluo_analysis.cpp
+++ b/CuFluo_analysis.cpp
@@ -4,15 +4,15 @@
 // and save peak position and uncertainty
 // changes by VH 210906: to eliminate hardcoded absolute paths, uses location of the analysis root files as additional input argument (accordingly changed in filename_creator.sh)
 
-#include "../sls_detector_calibration/jungfrauCommonHeader.h"
-#include "../sls_detector_calibration/jungfrauCommonFunctions.h"
+#include "sls_detector_calibration/jungfrauCommonHeader.h"
+#include "sls_detector_calibration/jungfrauCommonFunctions.h"
 
-#include "../sls_detector_calibration/jungfrauFile.C"
-#include "../sls_detector_calibration/jungfrauPixelMask.C"
-#include "../sls_detector_calibration/jungfrauPedestal.C"
+#include "sls_detector_calibration/jungfrauFile.C"
+#include "sls_detector_calibration/jungfrauPixelMask.C"
+#include "sls_detector_calibration/jungfrauPedestal.C"
 
-#include "../sls_detector_calibration/energyCalibration.h"
-#include "../sls_detector_calibration/energyCalibration.cpp"
+#include "sls_detector_calibration/energyCalibration.h"
+#include "sls_detector_calibration/energyCalibration.cpp"
 
 #include "TGraph.h"
 #include "TGraphErrors.h"
diff --git a/sls_detector_calibration/MovingStat.h b/sls_detector_calibration/MovingStat.h
new file mode 100755
index 0000000..a5f7086
--- /dev/null
+++ b/sls_detector_calibration/MovingStat.h
@@ -0,0 +1,145 @@
+#ifndef MOVINGSTAT_H
+#define MOVINGSTAT_H 
+
+#include <math.h>
+
+
+class MovingStat
+    {
+
+      /** @short approximated moving average structure */
+    public:
+
+      
+      /** constructor
+	  \param nn number of samples  parameter to be used
+      */
+    MovingStat(int nn=1000) : n(nn), m_n(0), u(0) {}
+
+	/**
+	   clears the moving average number of samples parameter, mean and standard deviation
+	*/
+        void Clear()
+        {
+	  m_n = 0; 
+	  m_newM=0; 
+	  m_newM2=0;
+	  u = 0;
+        }
+	
+	void ResetUpdates()
+        {
+	  u = 0;
+        }
+	
+	/** sets number of samples parameter
+	    \param i number of samples  parameter to be set 
+	*/
+	
+	void SetN(int i) {if (i>=1) n=i;};
+	
+	/** 
+	    gets number of samples  parameter
+	    \returns actual number of samples   parameter
+	*/
+	int GetN() {return n;};
+		
+	/** calculates the moving average i.e. adds if number of elements is lower than number of samples parameter, pushes otherwise
+	   \param x value to calculate the moving average
+	*/
+	inline void Calc(double x) {
+	  u++;
+	  if (m_n<n) Add(x);
+	  else Push(x);
+	}
+	/** adds the element to the accumulated average and standard deviation 
+	   \param x value to add
+	*/
+	inline void Add(double x) {
+	  m_n++;
+
+	  if (m_n == 1)
+            {
+	      m_newM = x;
+	      m_newM2 = x*x;
+            }             else             {
+                m_newM = m_newM + x;
+		m_newM2 = m_newM2 + x*x;
+	
+            }
+
+	}
+
+
+        inline void Push(double x)
+        {
+	  /** adds the element to the accumulated average and squared mean, while subtracting the current value of the average and squared average
+	      \param x value to push
+	  */
+	  if (m_n == 0)
+            {
+	      m_newM = x;
+	      m_newM2 = x*x;
+	      m_n++;
+            }             else             {
+                m_newM = m_newM + x - m_newM/m_n;
+		m_newM2 = m_newM2 + x*x - m_newM2/m_n;
+            }
+
+        }
+
+	/** returns the current number of elements of the moving average
+	    \returns  returns the current number of elements of the moving average
+	 */
+        int NumDataValues() const
+        {
+            return m_n;
+        }
+	/** returns the mean, 0 if no elements are inside
+	    \returns  returns the mean
+	 */
+        inline double Mean() const
+        {
+            return (m_n > 0) ? m_newM/m_n : 0.0;
+        }
+
+	/** returns the squared mean, 0 if no elements are inside
+	    \returns  returns the squared average
+	 */
+	double M2() const
+        {
+            return ( (m_n > 1) ? m_newM2/m_n : 0.0 );
+        }
+
+	/** returns the variance, 0 if no elements are inside
+	    \returns  returns the variance
+	 */
+        inline double Variance() const
+        {
+	  return ( (m_n > 1) ? (M2()-Mean()*Mean()) : 0.0 );
+        }
+
+	/** returns the standard deviation, 0 if no elements are inside
+	    \returns  returns the standard deviation
+	 */
+        inline double StandardDeviation() const
+        {
+	  return ( (Variance() > 0) ? sqrt( Variance() ) : -1 );
+        }
+
+	/** returns the number of times Calc was called
+	    \returns the number of times Calc was called
+	*/
+        inline double GetNUpdates() const
+        {
+	  return u;
+        }
+
+    private:
+	int n; /**< number of samples parameter */
+        int m_n; /**< current number of elements */
+	double m_newM; /**< accumulated average */
+	double m_newM2; /**< accumulated squared average */
+	int u; /**< number of times Calc was called (number of updates) */
+    };
+#endif
diff --git a/sls_detector_calibration/energyCalibration.cpp b/sls_detector_calibration/energyCalibration.cpp
new file mode 100644
index 0000000..7461eca
--- /dev/null
+++ b/sls_detector_calibration/energyCalibration.cpp
@@ -0,0 +1,798 @@
+#include "energyCalibration.h"
+
+#ifdef __CINT
+#define MYROOT
+#endif
+
+
+#ifdef MYROOT
+#include <TMath.h>
+#include <TH1F.h>
+#include <TH2F.h>
+#include <TGraphErrors.h>
+#endif
+
+#include <iostream>
+
+#define max(a,b) ((a) > (b) ? (a) : (b))
+#define min(a,b) ((a) < (b) ? (a) : (b))
+#define ELEM_SWAP(a,b) { register int t=(a);(a)=(b);(b)=t; }
+
+
+using namespace std;
+
+#ifdef MYROOT
+
+Double_t energyCalibrationFunctions::pedestal(Double_t *x, Double_t *par) { 
+    return par[0]-par[1]*sign*x[0];                                      
+}
+
+
+Double_t energyCalibrationFunctions::gaussChargeSharing(Double_t *x, Double_t *par) { 
+  Double_t f, arg=0;
+  // Gaussian exponent
+  if (par[3]!=0) {
+    arg=sign*(x[0]-par[2])/par[3];
+  }
+  // the Gaussian
+  f=TMath::Exp(-1*arg*arg/2.); 
+  // Gaussian + error function
+  f=f+par[5]/2.*(TMath::Erfc(arg/(TMath::Sqrt(2.)))); 
+  // Gaussian + error function + pedestal
+  return par[4]*f+pedestal(x,par);                                       
+}
+
+Double_t energyCalibrationFunctions::gaussChargeSharingKb(Double_t *x, Double_t *par) { 
+  Double_t f, arg=0,argb=0;
+  // Gaussian exponent
+  if (par[3]!=0) {
+    arg=sign*(x[0]-par[2])/par[3];
+    argb=sign*(x[0]-(par[6]*par[2]))/par[3]; // using absolute kb mean might seem better but like this the ratio can be fixed
+  }
+  // the Gaussian
+  f=TMath::Exp(-1*arg*arg/2.); 
+  f=f+par[7]*(TMath::Exp(-1*argb*argb/2.));
+  // Gaussian + error function
+  f=f+par[5]/2.*(TMath::Erfc(arg/(TMath::Sqrt(2.)))); 
+  f=f+par[7]*par[5]/2.*(TMath::Erfc(argb/(TMath::Sqrt(2.)))); 
+  // Gaussian + error function + pedestal
+  return par[4]*f+pedestal(x,par);                                       
+}
+
+Double_t energyCalibrationFunctions::gaussChargeSharingKaDoublet(Double_t *x, Double_t *par) { 
+  Double_t f, f2, arg=0, arg2=0;
+  // Gaussian exponent
+  if (par[3]!=0) {
+    arg=sign*(x[0]-par[2])/par[3];
+    arg2=sign*(x[0]-par[6])/par[3];
+  }
+  // the Gaussian
+  f=TMath::Exp(-1*arg*arg/2.); 
+  f2=TMath::Exp(-1*arg2*arg2/2.); 
+  // Gaussian + error function
+  f=f+par[5]/2.*(TMath::Erfc(arg/(TMath::Sqrt(2.)))); 
+  f2=f2+par[5]/2.*(TMath::Erfc(arg/(TMath::Sqrt(2.)))); // shouldn't this be arg2?
+  // Gaussian + error function + pedestal
+  return par[4]*f+par[7]*f2+pedestal(x,par);                                       
+}
+
+Double_t energyCalibrationFunctions::gaussChargeSharingPixel(Double_t *x, Double_t *par) { 
+  Double_t f;
+  if (par[3]<=0 || par[2]*(*x)<=0 || par[5]<0 || par[4]<=0) return 0;
+
+  Double_t pp[3];
+
+  pp[0]=0;
+  pp[1]=par[2];
+  pp[2]=par[3];
+  
+
+  f=(par[5]-par[6]*(TMath::Log(*x/par[2])))*erfBox(x,pp);
+  f+=par[4]*TMath::Gaus(*x, par[2], par[3], kTRUE);
+  return f+pedestal(x,par);                                       
+}
+
+Double_t energyCalibrationFunctions::erfBox(Double_t *z, Double_t *par) {
+
+  
+  
+  Double_t m=par[0];
+  Double_t M=par[1];
+  
+  if (par[0]>par[1]) {
+    m=par[1];
+    M=par[0];
+  }
+  
+  if (m==M)
+    return 0;
+
+
+  if (par[2]<=0) {
+    if (*z>=m && *z<=M)
+      return 1./(M-m);
+    else
+      return 0;
+
+  }
+
+  return (TMath::Erfc((z[0]-M)/par[2])-TMath::Erfc((z[0]-m)/par[2]))*0.5/(M-m);
+
+}
+
+
+// basic erf function
+Double_t energyCalibrationFunctions::erfFunction(Double_t *x, Double_t *par) {	
+  double arg=0;
+  if (par[1]!=0) arg=(par[0]-x[0])/par[1];
+  return ((par[2]/2.*(1+TMath::Erf(sign*arg/(TMath::Sqrt(2))))));  
+};
+  
+
+Double_t energyCalibrationFunctions::erfFunctionChargeSharing(Double_t *x, Double_t *par) {	       
+  Double_t f;		
+  
+  f=erfFunction(x, par+2)*(1+par[5]*(par[2]-x[0]))+par[0]-par[1]*x[0]*sign;		
+  return f;								
+};
+  
+  
+Double_t energyCalibrationFunctions::erfFuncFluo(Double_t *x, Double_t *par) {
+  Double_t f;							       
+  f=erfFunctionChargeSharing(x, par)+erfFunction(x, par+6)*(1+par[9]*(par[6]-x[0])); 
+  return f;							
+};
+#endif
+
+double energyCalibrationFunctions::median(double *x, int n){
+  // sorts x into xmed array and returns median 
+  // n is number of values already in the xmed array
+  double xmed[n];
+  int k,i,j;
+
+  for (i=0; i<n; i++) {
+    k=0;
+    for (j=0; j<n; j++) {
+      if(*(x+i)>*(x+j))
+	k++;
+      if (*(x+i)==*(x+j)) {
+	if (i>j) 
+	  k++;
+      }
+    }
+    xmed[k]=*(x+i);
+  }
+  k=n/2;
+  return xmed[k];
+}
+
+
+int energyCalibrationFunctions::quick_select(int arr[], int n){
+    int low, high ;
+    int median;
+    int middle, ll, hh;
+
+    low = 0 ; high = n-1 ; median = (low + high) / 2;
+    for (;;) {
+        if (high <= low) /* One element only */
+            return arr[median] ;
+
+        if (high == low + 1) {  /* Two elements only */
+            if (arr[low] > arr[high])
+                ELEM_SWAP(arr[low], arr[high]) ;
+            return arr[median] ;
+        }
+
+    /* Find median of low, middle and high items; swap into position low */
+    middle = (low + high) / 2;
+    if (arr[middle] > arr[high])    ELEM_SWAP(arr[middle], arr[high]) ;
+    if (arr[low] > arr[high])       ELEM_SWAP(arr[low], arr[high]) ;
+    if (arr[middle] > arr[low])     ELEM_SWAP(arr[middle], arr[low]) ;
+
+    /* Swap low item (now in position middle) into position (low+1) */
+    ELEM_SWAP(arr[middle], arr[low+1]) ;
+
+    /* Nibble from each end towards middle, swapping items when stuck */
+    ll = low + 1;
+    hh = high;
+    for (;;) {
+        do ll++; while (arr[low] > arr[ll]) ;
+        do hh--; while (arr[hh]  > arr[low]) ;
+
+        if (hh < ll)
+        break;
+
+        ELEM_SWAP(arr[ll], arr[hh]) ;
+    }
+
+    /* Swap middle item (in position low) back into correct position */
+    ELEM_SWAP(arr[low], arr[hh]) ;
+
+    /* Re-set active partition */
+    if (hh <= median)
+        low = ll;
+        if (hh >= median)
+        high = hh - 1;
+    }
+}
+
+int energyCalibrationFunctions::kth_smallest(int *a, int n, int k){
+    register int i,j,l,m ;
+    register double x ;
+
+    l=0 ; m=n-1 ;
+    while (l<m) {
+        x=a[k] ;
+        i=l ;
+        j=m ;
+        do {
+            while (a[i]<x) i++ ;
+            while (x<a[j]) j-- ;
+            if (i<=j) {
+                ELEM_SWAP(a[i],a[j]) ;
+                i++ ; j-- ;
+            }
+        } while (i<=j) ;
+        if (j<k) l=i ;
+        if (k<i) m=j ;
+    }
+    return a[k] ;
+}
+
+
+
+#ifdef MYROOT
+Double_t energyCalibrationFunctions::spectrum(Double_t *x, Double_t *par) {
+  return gaussChargeSharing(x,par);	
+}
+
+Double_t energyCalibrationFunctions::spectrumkb(Double_t *x, Double_t *par) {
+  return gaussChargeSharingKb(x,par);	
+}
+
+Double_t energyCalibrationFunctions::spectrumkadoublet(Double_t *x, Double_t *par) {
+  return gaussChargeSharingKaDoublet(x,par);	
+}
+
+Double_t energyCalibrationFunctions::spectrumPixel(Double_t *x, Double_t *par) {
+  return gaussChargeSharingPixel(x,par);	
+}
+
+
+Double_t energyCalibrationFunctions::scurve(Double_t *x, Double_t *par) {
+  return erfFunctionChargeSharing(x,par);
+}
+
+
+Double_t energyCalibrationFunctions::scurveFluo(Double_t *x, Double_t *par) {
+  return erfFuncFluo(x,par);						
+}
+#endif
+
+energyCalibration::energyCalibration() : 
+#ifdef MYROOT	
+					 fit_min(-1),	
+					 fit_max(-1),	
+					 bg_offset(-1),	
+					 bg_slope(-1),	
+					 flex(-1),
+					 noise(-1),
+					 ampl(-1),
+					 cs_slope(-1),
+					 kb_mean(-1),
+					 kb_frac(-1),
+					 mean2(-1),
+					 ampl2(-1),
+					 fscurve(NULL),
+					 fspectrum(NULL),
+					 fspectrumkb(NULL),
+					 fspectrumkadoublet(NULL),
+#endif
+					 funcs(NULL),
+					 plot_flag(1), // fit parameters output to screen
+					 cs_flag(1)
+{
+
+#ifdef MYROOT	
+  funcs=new energyCalibrationFunctions();
+  
+  fscurve=new TF1("fscurve",funcs,&energyCalibrationFunctions::scurve,0,1000,6,"energyCalibrationFunctions","scurve");	
+  fscurve->SetParNames("Background Offset","Background Slope","Inflection Point","Noise RMS", "Number of Photons","Charge Sharing Slope");
+
+  fspectrum=new TF1("fspectrum",funcs,&energyCalibrationFunctions::spectrum,0,1000,6,"energyCalibrationFunctions","spectrum");	
+  fspectrum->SetParNames("Background Pedestal","Background slope", "Peak position","Noise RMS", "Number of Photons","Charge Sharing Pedestal");
+  fspectrumkb=new TF1("fspectrumkb",funcs,&energyCalibrationFunctions::spectrumkb,0,1000,8,"energyCalibrationFunctions","spectrumkb");	
+  fspectrumkb->SetParNames("Background Pedestal","Background slope", "Peak position","Noise RMS", "Number of Photons","Charge Sharing Pedestal","kb mean","kb frac");
+
+  fspectrumkadoublet=new TF1("fspectrumkadoublet",funcs,&energyCalibrationFunctions::spectrumkadoublet,0,1000,8,"energyCalibrationFunctions","spectrumkadoublet");	
+  fspectrumkadoublet->SetParNames("Background Pedestal","Background slope", "Peak position","Noise RMS", "Number of Photons","Charge Sharing Pedestal","ka2 mean","n2");
+
+  fspixel=new TF1("fspixel",funcs,&energyCalibrationFunctions::spectrumPixel,0,1000,7,"energyCalibrationFunctions","spectrumPixel");	
+  fspixel->SetParNames("Background Pedestal","Background slope", "Peak position","Noise RMS", "Number of Photons","Charge Sharing Pedestal","Corner");
+
+#endif
+
+
+}
+
+
+
+void energyCalibration::fixParameter(int ip, Double_t val){
+  
+  fscurve->FixParameter(ip, val);
+  fspectrum->FixParameter(ip, val);
+  fspectrumkb->FixParameter(ip, val);
+  fspectrumkadoublet->FixParameter(ip, val);
+}
+
+
+void energyCalibration::releaseParameter(int ip){
+
+  fscurve->ReleaseParameter(ip);
+  fspectrum->ReleaseParameter(ip);
+  fspectrumkb->ReleaseParameter(ip);
+  fspectrumkadoublet->ReleaseParameter(ip);
+}
+
+
+
+
+
+
+
+energyCalibration::~energyCalibration(){ 
+#ifdef MYROOT
+  delete fscurve;
+  delete fspectrum;
+  delete fspectrumkb;
+  delete fspectrumkadoublet;
+#endif
+  
+}
+
+#ifdef MYROOT
+
+
+
+
+
+TH1F* energyCalibration::createMedianHistogram(TH2F* h2, int ch0, int nch, int direction) {
+
+  if (h2==NULL || nch==0)
+    return NULL;
+
+  double *x=new double[nch];
+  TH1F *h1=NULL;
+  
+  double val=-1;
+
+  if (direction==0) {
+    h1=new TH1F("median","Median",h2->GetYaxis()->GetNbins(),h2->GetYaxis()->GetXmin(),h2->GetYaxis()->GetXmax());
+    for (int ib=0; ib<h1->GetXaxis()->GetNbins(); ib++) {
+      for (int ich=0; ich<nch; ich++) {
+	x[ich]=h2->GetBinContent(ch0+ich+1,ib+1);
+      }
+      val=energyCalibrationFunctions::median(x, nch);
+      h1->SetBinContent(ib+1,val);
+    }
+  } else if (direction==1) {   
+    h1=new TH1F("median","Median",h2->GetXaxis()->GetNbins(),h2->GetXaxis()->GetXmin(),h2->GetXaxis()->GetXmax());
+    for (int ib=0; ib<h1->GetYaxis()->GetNbins(); ib++) {
+      for (int ich=0; ich<nch; ich++) {
+	x[ich]=h2->GetBinContent(ib+1,ch0+ich+1);
+      }
+      val=energyCalibrationFunctions::median(x, nch);
+      h1->SetBinContent(ib+1,val);
+    }
+  }
+  delete [] x;
+
+  return h1;
+
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+void energyCalibration::setStartParameters(Double_t *par){	
+  bg_offset=par[0];						
+  bg_slope=par[1];						
+  flex=par[2];							
+  noise=par[3];							
+  ampl=par[4];							
+  cs_slope=par[5];						
+}
+
+void energyCalibration::setStartParametersKb(Double_t *par){	
+  bg_offset=par[0];						
+  bg_slope=par[1];						
+  flex=par[2];							
+  noise=par[3];							
+  ampl=par[4];							
+  cs_slope=par[5];						
+  kb_mean=par[6];
+  kb_frac=par[7];
+  //fit_min = 400; // used for soleil flat field
+  //fit_max = 800;
+}
+
+void energyCalibration::setStartParametersKaDoublet(Double_t *par){	
+  bg_offset=par[0];						
+  bg_slope=par[1];						
+  flex=par[2];							
+  noise=par[3];							
+  ampl=par[4];							
+  cs_slope=par[5];						
+  mean2=par[6];
+  ampl2=par[7];
+  //fit_min = 400; // used for soleil flat field
+  //fit_max = 800;
+}
+
+
+void energyCalibration::getStartParameters(Double_t *par){	
+  par[0]=bg_offset;
+  par[1]=bg_slope;
+  par[2]=flex;
+  par[3]=noise;
+  par[4]=ampl;				
+  par[5]=cs_slope;				
+}
+
+#endif
+int energyCalibration::setChargeSharing(int p) {
+  if (p>=0) {
+    cs_flag=p; 
+#ifdef MYROOT  
+    if (p) {
+      fscurve->ReleaseParameter(5);
+      fspectrum->ReleaseParameter(1);
+      fspectrumkb->ReleaseParameter(1);
+      fspectrumkadoublet->ReleaseParameter(1);
+    } else {
+      fscurve->FixParameter(5,0);
+      fspectrum->FixParameter(1,0);
+      fspectrumkb->FixParameter(1,0);
+      fspectrumkadoublet->FixParameter(1,0);
+    }
+#endif
+  } 
+  
+  return cs_flag;
+}
+
+ 
+#ifdef MYROOT  
+void energyCalibration::initFitFunction(TF1 *fun, TH1 *h1) {
+
+ Double_t min=fit_min, max=fit_max; 
+
+ Double_t mypar[6];
+  
+ if (max==-1) 
+    max=h1->GetXaxis()->GetXmax();
+ 
+  if (min==-1) 
+    min=h1->GetXaxis()->GetXmin();
+
+  
+  if (bg_offset==-1)
+    mypar[0]=0;
+  else
+    mypar[0]=bg_offset;
+
+
+  if (bg_slope==-1)
+    mypar[1]=0;
+  else
+    mypar[1]=bg_slope;
+
+
+  if (flex==-1)
+    mypar[2]=(min+max)/2.;
+  else
+    mypar[2]=flex;
+
+
+  if (noise==-1)
+    mypar[3]=0.1;
+  else
+    mypar[3]=noise;
+
+  if (ampl==-1)
+    mypar[4]=h1->GetBinContent(h1->GetXaxis()->FindBin(0.5*(max+min)));
+  else
+    mypar[4]=ampl;
+
+  if (cs_slope==-1)
+    mypar[5]=0;
+  else
+    mypar[5]=cs_slope;
+ 
+  fun->SetParameters(mypar);
+
+  fun->SetRange(min,max);
+
+}
+
+void energyCalibration::initFitFunctionKb(TF1 *fun, TH1 *h1) {
+
+ Double_t min=fit_min, max=fit_max; 
+
+ Double_t mypar[8];
+  
+ if (max==-1) 
+    max=h1->GetXaxis()->GetXmax();
+ 
+  if (min==-1) 
+    min=h1->GetXaxis()->GetXmin();
+
+  
+  if (bg_offset==-1)
+    mypar[0]=0;
+  else
+    mypar[0]=bg_offset;
+
+
+  if (bg_slope==-1)
+    mypar[1]=0;
+  else
+    mypar[1]=bg_slope;
+
+
+  if (flex==-1)
+    mypar[2]=(min+max)/2.;
+  else
+    mypar[2]=flex;
+
+
+  if (noise==-1)
+    mypar[3]=0.1;
+  else
+    mypar[3]=noise;
+
+  if (ampl==-1)
+    mypar[4]=h1->GetBinContent(h1->GetXaxis()->FindBin(0.5*(max+min)));
+  else
+    mypar[4]=ampl;
+
+  if (cs_slope==-1)
+    mypar[5]=0;
+  else
+    mypar[5]=cs_slope;
+
+  if (kb_mean==-1)
+    mypar[6]=0;
+  else
+    mypar[6]=kb_mean;
+
+  if (kb_frac==-1)
+    mypar[7]=0;
+  else
+    mypar[7]=kb_frac;
+ 
+  fun->SetParameters(mypar);
+
+  fun->SetRange(min,max);
+
+}
+
+void energyCalibration::initFitFunctionKaDoublet(TF1 *fun, TH1 *h1) {
+
+ Double_t min=fit_min, max=fit_max; 
+
+ Double_t mypar[8];
+  
+ if (max==-1) 
+    max=h1->GetXaxis()->GetXmax();
+ 
+  if (min==-1) 
+    min=h1->GetXaxis()->GetXmin();
+
+  
+  if (bg_offset==-1)
+    mypar[0]=0;
+  else
+    mypar[0]=bg_offset;
+
+
+  if (bg_slope==-1)
+    mypar[1]=0;
+  else
+    mypar[1]=bg_slope;
+
+
+  if (flex==-1)
+    mypar[2]=(min+max)/2.;
+  else
+    mypar[2]=flex;
+
+
+  if (noise==-1)
+    mypar[3]=0.1;
+  else
+    mypar[3]=noise;
+
+  if (ampl==-1)
+    mypar[4]=h1->GetBinContent(h1->GetXaxis()->FindBin(0.5*(max+min)));
+  else
+    mypar[4]=ampl;
+
+  if (cs_slope==-1)
+    mypar[5]=0;
+  else
+    mypar[5]=cs_slope;
+
+  if (mean2==-1)
+    mypar[6]=0;
+  else
+    mypar[6]=mean2;
+
+  if (ampl2==-1)
+    mypar[7]=0;
+  else
+    mypar[7]=ampl2;
+ 
+  fun->SetParameters(mypar);
+
+  fun->SetRange(min,max);
+
+}
+
+TF1* energyCalibration::fitFunction(TF1 *fun, TH1 *h1, Double_t *mypar, Double_t *emypar) {
+
+
+  TF1* fitfun;
+
+  char fname[100];
+
+  strcpy(fname, fun->GetName());
+
+  if (plot_flag) {
+    h1->Fit(fname,"R0Q");
+  } else
+    h1->Fit(fname,"R0Q");
+
+  fitfun= h1->GetFunction(fname);
+  fitfun->GetParameters(mypar);
+  for (int ip=0; ip<6; ip++) {
+    emypar[ip]=fitfun->GetParError(ip);
+  }
+  return fitfun;
+}
+
+TF1* energyCalibration::fitFunctionKb(TF1 *fun, TH1 *h1, Double_t *mypar, Double_t *emypar) {
+
+
+  TF1* fitfun;
+
+  char fname[100];
+
+  strcpy(fname, fun->GetName());
+
+  if (plot_flag) {
+    h1->Fit(fname,"R0Q");
+  } else
+    h1->Fit(fname,"R0Q");
+
+  fitfun= h1->GetFunction(fname);
+  fitfun->GetParameters(mypar);
+  for (int ip=0; ip<8; ip++) {
+    emypar[ip]=fitfun->GetParError(ip);
+  }
+  return fitfun;
+}
+
+TF1* energyCalibration::fitFunctionKaDoublet(TF1 *fun, TH1 *h1, Double_t *mypar, Double_t *emypar) {
+
+
+  TF1* fitfun;
+
+  char fname[100];
+
+  strcpy(fname, fun->GetName());
+
+  if (plot_flag) {
+    h1->Fit(fname,"R0Q");
+  } else
+    h1->Fit(fname,"R0Q");
+
+
+  fitfun= h1->GetFunction(fname);
+  fitfun->GetParameters(mypar);
+  for (int ip=0; ip<8; ip++) {
+    emypar[ip]=fitfun->GetParError(ip);
+  }
+  return fitfun;
+}
+
+TF1* energyCalibration::fitSCurve(TH1 *h1, Double_t *mypar, Double_t *emypar) {
+  initFitFunction(fscurve,h1);
+  return fitFunction(fscurve, h1, mypar, emypar);
+}
+
+
+
+
+
+TF1* energyCalibration::fitSpectrum(TH1 *h1, Double_t *mypar, Double_t *emypar) {
+  initFitFunction(fspectrum,h1);
+  return fitFunction(fspectrum, h1, mypar, emypar);
+}
+
+TF1* energyCalibration::fitSpectrumKb(TH1 *h1, Double_t *mypar, Double_t *emypar) {
+  initFitFunctionKb(fspectrumkb,h1);
+  return fitFunctionKb(fspectrumkb, h1, mypar, emypar);
+}
+
+TF1* energyCalibration::fitSpectrumKaDoublet(TH1 *h1, Double_t *mypar, Double_t *emypar) {
+  initFitFunctionKaDoublet(fspectrumkadoublet,h1);
+  return fitFunctionKaDoublet(fspectrumkadoublet, h1, mypar, emypar);
+}
+
+
+TGraphErrors* energyCalibration::linearCalibration(int nscan, Double_t *en, Double_t *een, Double_t *fl, Double_t *efl, Double_t &gain, Double_t &off, Double_t &egain, Double_t &eoff) { 
+ 
+  TGraphErrors *gr;
+ 
+  Double_t  mypar[2];
+ 
+  gr = new TGraphErrors(nscan,en,fl,een,efl);
+
+  if (plot_flag) {
+    gr->Fit("pol1");
+    gr->SetMarkerStyle(20);
+  } else
+    gr->Fit("pol1","0Q");
+        
+  TF1 *fitfun= gr->GetFunction("pol1"); 
+  fitfun->GetParameters(mypar);
+
+  egain=fitfun->GetParError(1);
+  eoff=fitfun->GetParError(0);
+
+  gain=funcs->setScanSign()*mypar[1];
+  
+  off=mypar[0]; 
+  
+  return gr;
+}
+
+
+TGraphErrors* energyCalibration::calibrate(int nscan, Double_t *en, Double_t *een, TH1F **h1, Double_t &gain, Double_t &off, Double_t &egain, Double_t &eoff, int integral) {
+
+  TH1F *h;
+
+  Double_t mypar[6], emypar[6];
+  Double_t fl[nscan], efl[nscan];
+  
+
+  for (int ien=0; ien<nscan; ien++) {
+    h=h1[ien];
+    if (integral)
+      fitSCurve(h,mypar,emypar);
+    else
+      fitSpectrum(h,mypar,emypar);
+
+    fl[ien]=mypar[2];
+    efl[ien]=emypar[2];
+  }
+  return linearCalibration(nscan,en,een,fl,efl,gain,off, egain, eoff);
+
+}
+
+#endif
+
+
+
+
diff --git a/sls_detector_calibration/energyCalibration.h b/sls_detector_calibration/energyCalibration.h
new file mode 100644
index 0000000..cef73d3
--- /dev/null
+++ b/sls_detector_calibration/energyCalibration.h
@@ -0,0 +1,469 @@
+#ifndef ENERGYCALIBRATION_H
+#define ENERGYCALIBRATION_H
+
+#ifdef __CINT__
+#define MYROOT
+#endif
+
+#ifdef G__ROOT
+#define MYROOT
+#endif
+
+
+#ifdef __MAKECINT__
+#define MYROOT
+#endif
+
+#ifdef ROOT_VERSION
+#define MYROOT
+#endif
+
+#define MYROOT
+
+#ifdef MYROOT
+#include <TROOT.h>
+#include <TF1.h>
+class TH1F;
+class TH2F;
+class TGraphErrors;
+#endif
+
+
+using namespace std;
+
+
+
+
+const double conven=1000./3.6; /**< electrons/keV */
+const double el=1.67E-4; /**< electron charge in fC */
+
+
+
+/**
+ \mainpage Common Root library for SLS detectors data analysis
+ *
+ * \section intro_sec Introduction
+ We know very well s-curves etc. but at the end everybody uses different functions ;-).
+
+ * \subsection mot_sec Motivation
+ It would be greate to use everybody the same functions...
+
+*/
+
+
+/**
+ * 
+ *
+@libdoc The energiCalibration class contains all the necessary functions for s-curve fitting and linear calibration of the threshold.
+ *
+ * @short Energy calibration functions
+ * @author Anna Bergamaschi
+ * @version 0.1alpha
+
+
+ */
+
+/** 
+  class containing all the possible energy calibration functions (scurves with and without charge sharing, gaussian spectrum with and without charge sharing, possibility of chosing the sign of the X-axis)
+
+*/
+class energyCalibrationFunctions {
+
+ public:
+  
+  energyCalibrationFunctions(int s=-1) {setScanSign(s);};
+  
+    /** sets scan sign
+	\param s can be 1 (energy and x-axis have the same direction) or -1 (energy and x-axis have opposite directions) otherwise gets
+	\returns current scan sign can be 1 (energy and x-axis have the same direction) or -1 (energy and x-axis have opposite directions)
+    */
+  int setScanSign(int s=0) {if (s==1 || s==-1) sign=s; return sign;};;
+  
+
+#ifdef MYROOT
+  /** 
+      Gaussian Function with charge sharing pedestal
+      par[0] is the absolute height of the background pedestal
+      par[1] is the slope of the background pedestal
+  */
+  Double_t pedestal(Double_t *x, Double_t *par);
+
+  /** 
+      Gaussian Function with charge sharing pedestal
+      par[0] is the absolute height of the background pedestal
+      par[1] is the slope of the background pedestal
+      par[2] is the gaussian peak position
+      par[3] is the RMS of the gaussian (and of the pedestal)
+      par[4] is the height of the function
+      par[5] is the fractional height of the charge sharing pedestal (scales with par[3])
+  */
+  Double_t gaussChargeSharing(Double_t *x, Double_t *par);
+  Double_t gaussChargeSharingKb(Double_t *x, Double_t *par);
+  Double_t gaussChargeSharingKaDoublet(Double_t *x, Double_t *par);
+  /** 
+      Gaussian Function with charge sharing pedestal
+      par[0] is the absolute height of the background pedestal
+      par[1] is the slope of the background pedestal
+      par[2] is the gaussian peak position
+      par[3] is the RMS of the gaussian (and of the pedestal)
+      par[4] is the height of the function
+      par[5] is the fractional height of the charge sharing pedestal (scales with par[3])
+  */
+  Double_t gaussChargeSharingPixel(Double_t *x, Double_t *par);
+
+  /** 
+      Basic erf function
+      par[0] is the inflection point
+      par[1] is the RMS 
+      par[2] is the amplitude
+  */
+Double_t erfFunction(Double_t *x, Double_t *par) ;
+ Double_t erfBox(Double_t *z, Double_t *par);
+  /** Erf function with charge sharing slope
+      par[0] is the  pedestal 
+      par[1] is the  slope of the pedestal
+      par[2] is the inflection point
+      par[3] is the RMS 
+      par[4] is the amplitude
+      par[5] is the angual coefficient of the charge sharing slope (scales with par[3])
+  */
+Double_t erfFunctionChargeSharing(Double_t *x, Double_t *par);
+  
+  /** Double Erf function with charge sharing slope
+      par[0] is the  pedestal 
+      par[1] is the  slope of the pedestal
+      par[2] is the inflection point of the first energy
+      par[3] is the RMS of the first energy
+      par[4] is the amplitude of the first energy 
+      par[5] is the angual coefficient of the charge sharing slope of the first energy (scales with par[3])
+      par[6] is the inflection point of the second energy 
+      par[7] is the RMS of the second energy 
+      par[8] is the amplitude of the second energy 
+      par[9] is the angual coefficient of the charge sharing slope  of the second energy (scales with par[8])
+  */
+
+Double_t erfFuncFluo(Double_t *x, Double_t *par);
+
+  
+  /** 
+      static function Gaussian with charge sharing pedestal with the correct scan sign
+      par[0] is the absolute height of the background pedestal
+      par[1]  is the  slope of the pedestal 
+      par[2] is the gaussian peak position
+      par[3] is the RMS of the gaussian (and of the pedestal)
+      par[4] is the height of the function
+      par[5] is the fractional height of the charge sharing pedestal (scales with par[4]
+  */
+  Double_t spectrum(Double_t *x, Double_t *par);
+  Double_t spectrumkb(Double_t *x, Double_t *par);
+  Double_t spectrumkadoublet(Double_t *x, Double_t *par);
+
+  /** 
+      static function Gaussian with charge sharing pedestal with the correct scan sign
+      par[0] is the absolute height of the background pedestal
+      par[1]  is the  slope of the pedestal 
+      par[2] is the gaussian peak position
+      par[3] is the RMS of the gaussian (and of the pedestal)
+      par[4] is the height of the function
+      par[5] is the fractional height of the charge sharing pedestal (scales with par[4]
+  */
+  Double_t spectrumPixel(Double_t *x, Double_t *par);
+
+
+ /** Erf function with charge sharing slope with the correct scan sign
+      par[0] is the  pedestal 
+      par[1] is the  slope of the pedestal
+      par[2] is the inflection point
+      par[3] is the RMS 
+      par[4] is the amplitude
+      par[5] is the angual coefficient of the charge sharing slope (scales with par[3])
+  */
+  Double_t scurve(Double_t *x, Double_t *par);
+
+
+
+   /** Double Erf function with charge sharing slope
+      par[0] is the  pedestal 
+      par[1] is the  slope of the pedestal 
+      par[2] is the inflection point of the first energy
+      par[3] is the RMS of the first energy
+      par[4] is the amplitude of the first energy 
+      par[5] is the angual coefficient of the charge sharing slope of the first energy (scales with par[3]) 
+      par[6] is the inflection point of the second energy 
+      par[7] is the RMS of the second energy 
+      par[8] is the amplitude of the second energy 
+      par[9] is the angual coefficient of the charge sharing slope  of the second energy (scales with par[8]) 
+  */
+  Double_t scurveFluo(Double_t *x, Double_t *par);
+
+#endif
+
+/** Calculates the median of an array of n elements */
+  static double median(double *x, int n);
+/**  Calculates the median of an array of n elements  (swaps the arrays!)*/
+  static int quick_select(int arr[], int n);
+/** Calculates the median of an array of n elements (swaps the arrays!)*/
+  static int kth_smallest(int *a, int n, int k);
+
+ private:
+  int sign;
+    
+
+};
+
+/**
+   class alowing the energy calibration of photon counting and anlogue detectors
+
+*/
+
+class energyCalibration  {
+
+
+ public:
+  /** 
+      default constructor  -  creates the function with which the s-curves will be fitted
+  */
+  energyCalibration();
+    
+  /** 
+      default destructor - deletes the function with which the s-curves will be fitted 
+  */
+  ~energyCalibration();
+    
+    /** sets plot flag
+	\param p plot flag (-1 gets, 0 unsets, >0 plot)
+	\returns current plot flag
+    */
+  int setPlotFlag(int p=-1) {if (p>=0) plot_flag=p; return plot_flag;};
+
+    /** sets scan sign
+	\param s can be 1 (energy and x-axis have the same direction) or -1 (energy and x-axis have opposite directions) otherwise gets
+	\returns current scan sign can be 1 (energy and x-axis have the same direction) or -1 (energy and x-axis have opposite directions)
+    */
+  int setScanSign(int s=0) {return funcs->setScanSign(s);};
+  
+    /** sets plot flag
+	\param p plot flag (-1 gets, 0 unsets, >0 plot)
+	\returns current plot flag
+    */
+  int setChargeSharing(int p=-1);
+
+
+  void fixParameter(int ip, Double_t val);
+
+  void releaseParameter(int ip);
+
+#ifdef MYROOT
+
+  /**
+     Creates an histogram with the median of nchannels starting from a specified one. the direction on which it is mediated can be selected (defaults to x=0)
+     \param h2 2D histogram on which the median will be calculated
+     \param ch0 starting channel
+     \param nch number of channels to be mediated
+     \param direction can be either 0 (x, default) or 1 (y)
+     \returns a TH1F histogram with the X-axis as a clone of the h2 Y (if direction=0) or X (if direction=0) axis, and on the Y axis the median of the counts of the mediated channels f h2
+   */
+  static TH1F* createMedianHistogram(TH2F* h2, int ch0, int nch, int direction=0);
+
+
+  /** sets the s-curve fit range 
+      \param mi  minimum of the fit range (-1 is histogram x-min)
+      \param ma  maximum of the fit range (-1 is histogram x-max)
+  */
+  void setFitRange(Double_t mi, Double_t ma){fit_min=mi; fit_max=ma;};
+
+  /** gets the s-curve fit range 
+      \param mi reference for minimum of the fit range (-1 is histogram x-min)
+      \param ma reference for maximum of the fit range (-1 is histogram x-max)
+  */
+  void getFitRange(Double_t &mi, Double_t &ma){mi=fit_min; ma=fit_max;};
+
+
+/** set start parameters for the s-curve function
+    \param par parameters, -1 sets to auto-calculation
+      par[0] is the  pedestal 
+      par[1] is the  slope of the pedestal 
+      par[2] is the inflection point
+      par[3] is the RMS 
+      par[4] is the amplitude
+      par[5] is the angual coefficient of the charge sharing slope (scales with par[3]) -- always positive
+  */
+  void setStartParameters(Double_t *par);
+  void setStartParametersKb(Double_t *par);
+  void setStartParametersKaDoublet(Double_t *par);
+
+/** get start parameters for the s-curve function
+    \param par parameters, -1 means auto-calculated
+      par[0] is the  pedestal 
+      par[1] is the  slope of the pedestal 
+      par[2] is the inflection point
+      par[3] is the RMS 
+      par[4] is the amplitude
+      par[5] is the angual coefficient of the charge sharing slope (scales with par[3]) -- always positive
+  */
+  void getStartParameters(Double_t *par);
+
+  /**
+     fits histogram with the s-curve function
+     \param h1 1d-histogram to be fitted
+     \param mypar pointer to fit parameters array 
+     \param emypar pointer to fit parameter errors
+     \returns the fitted function - can be used e.g. to get the Chi2 or similar
+  */
+  TF1 *fitSCurve(TH1 *h1, Double_t *mypar, Double_t *emypar);
+
+
+  /**
+     fits histogram with the spectrum
+     \param h1 1d-histogram to be fitted
+     \param mypar pointer to fit parameters array 
+     \param emypar pointer to fit parameter errors
+     \returns the fitted function - can be used e.g. to get the Chi2 or similar
+  */
+  TF1 *fitSpectrum(TH1 *h1, Double_t *mypar, Double_t *emypar);
+  TF1 *fitSpectrumKb(TH1 *h1, Double_t *mypar, Double_t *emypar);
+  TF1 *fitSpectrumKaDoublet(TH1 *h1, Double_t *mypar, Double_t *emypar);
+
+
+  /**
+     calculates gain and offset for the set of inflection points
+     \param nscan number of energy scans
+     \param en array of energies (nscan long)
+     \param een array of errors on energies (nscan long) - can be NULL!
+     \param fl array of inflection points (nscan long)
+     \param efl array of errors on the inflection points (nscan long)
+     \param gain reference to gain resulting from the fit
+     \param off reference to offset resulting from the fit
+     \param egain reference to error on the gain resulting from the fit
+     \param eoff reference to the error on the offset resulting from the fit
+     \returns graph  energy vs inflection point
+  */
+  TGraphErrors* linearCalibration(int nscan, Double_t *en, Double_t *een, Double_t *fl, Double_t *efl, Double_t &gain, Double_t &off, Double_t &egain, Double_t &eoff);  
+
+  /**
+     calculates gain and offset for the set of energy scans
+     \param nscan number of energy scans
+     \param en array of energies (nscan long)
+     \param een array of errors on energies (nscan long) - can be NULL!
+     \param h1 array of TH1
+     \param gain reference to gain resulting from the fit
+     \param off reference to offset resulting from the fit
+     \param egain reference to error on the gain resulting from the fit
+     \param eoff reference to the error on the offset resulting from the fit
+     \returns graph  energy vs inflection point 
+  */
+  TGraphErrors* calibrateScurves(int nscan, Double_t *en, Double_t *een, TH1F **h1, Double_t &gain, Double_t &off, Double_t &egain, Double_t &eoff){return calibrate(nscan, en, een, h1, gain, off, egain, eoff, 1);};    
+
+  /**
+     calculates gain and offset for the set of energy spectra
+     \param nscan number of energy scans
+     \param en array of energies (nscan long)
+     \param een array of errors on energies (nscan long) - can be NULL!
+     \param h1 array of TH1
+     \param gain reference to gain resulting from the fit
+     \param off reference to offset resulting from the fit
+     \param egain reference to error on the gain resulting from the fit
+     \param eoff reference to the error on the offset resulting from the fit
+     \returns graph  energy vs peak
+  */
+  TGraphErrors* calibrateSpectra(int nscan, Double_t *en, Double_t *een, TH1F **h1, Double_t &gain, Double_t &off, Double_t &egain, Double_t &eoff){return calibrate(nscan, en, een, h1, gain, off, egain, eoff, 0);};  
+
+
+#endif
+ private:
+
+#ifdef MYROOT
+    /**
+     calculates gain and offset for the set of energies
+     \param nscan number of energy scans
+     \param en array of energies (nscan long)
+     \param een array of errors on energies (nscan long) - can be NULL!
+     \param h1 array of TH1
+     \param gain reference to gain resulting from the fit
+     \param off reference to offset resulting from the fit
+     \param egain reference to error on the gain resulting from the fit
+     \param eoff reference to the error on the offset resulting from the fit
+     \param integral 1 is an s-curve set (default), 0 spectra
+     \returns graph  energy vs peak/inflection point 
+  */
+  TGraphErrors* calibrate(int nscan, Double_t *en, Double_t *een, TH1F **h1, Double_t &gain, Double_t &off, Double_t &egain, Double_t &eoff, int integral=1);  
+
+
+  /** 
+      Initializes the start parameters and the range of the fit depending on the histogram characteristics and/or on the start parameters specified by the user 
+  \param fun pointer to function to be initialized
+  \param h1 histogram from which to extract the range and start parameters, if not already specified by the user
+  
+*/
+  
+  void initFitFunction(TF1 *fun, TH1 *h1);
+  void initFitFunctionKb(TF1 *fun, TH1 *h1);
+  void initFitFunctionKaDoublet(TF1 *fun, TH1 *h1);
+
+
+  /**
+     Performs the fit according to the flags specified and returns the fitted function
+     \param fun function to fit
+     \param h1 histogram to fit
+     \param mypar pointer to fit parameters array 
+     \param emypar pointer to fit parameter errors
+     \returns the fitted function - can be used e.g. to get the Chi2 or similar
+  */
+  TF1 *fitFunction(TF1 *fun, TH1 *h1, Double_t *mypar, Double_t *emypar);  
+  TF1 *fitFunctionKb(TF1 *fun, TH1 *h1, Double_t *mypar, Double_t *emypar);  
+  TF1 *fitFunctionKaDoublet(TF1 *fun, TH1 *h1, Double_t *mypar, Double_t *emypar);  
+
+#endif
+
+#ifdef MYROOT
+  Double_t fit_min; /**< minimum of the s-curve fitting range, -1 is histogram x-min */
+  Double_t fit_max; /**< maximum of the s-curve fitting range, -1 is histogram x-max */
+  
+  Double_t bg_offset; /**< start value for the background pedestal */
+  Double_t bg_slope; /**< start value for the background slope */
+  Double_t flex; /**< start value for the inflection point */
+  Double_t noise; /**< start value for the noise */
+  Double_t ampl; /**< start value for the number of photons */
+  Double_t cs_slope; /**< start value for the charge sharing slope */
+  Double_t kb_mean;
+  Double_t kb_frac;
+  Double_t mean2;
+  Double_t ampl2;
+
+  TF1 *fscurve; /**< function with which the s-curve will be fitted */
+
+  TF1 *fspectrum; /**< function with which the spectrum will be fitted */
+  TF1 *fspectrumkb; /**< function with which the spectrum will be fitted */
+  TF1 *fspectrumkadoublet; /**< function with which the spectrum will be fitted */
+
+  TF1 *fspixel; /**< function with which the spectrum will be fitted */
+
+#endif
+
+  energyCalibrationFunctions *funcs;
+  int plot_flag; /**< 0 does not plot, >0 plots (flags?) */
+  
+  int cs_flag; /**< 0 functions without charge sharing contribution, >0 with charge sharing contribution */
+
+};
+
+#endif
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/sls_detector_calibration/jungfrauCommonFunctions.h b/sls_detector_calibration/jungfrauCommonFunctions.h
new file mode 100644
index 0000000..3c78e64
--- /dev/null
+++ b/sls_detector_calibration/jungfrauCommonFunctions.h
@@ -0,0 +1,344 @@
+#ifndef JUNGFRAUCOMMONFUNCTIONS_H
+#define JUNGFRAUCOMMONFUNCTIONS_H
+
+#include <algorithm>
+
+int chipOfPixel(int i) { 
+  int chip = 0;
+  if (i/1024 >= 256 && i%1024 >= 0 && i%1024 < 256) {
+    chip = 1;
+  } else if (i/1024 >= 256 && i%1024 >= 256 && i%1024 < 512) {
+    chip = 2;
+  } else if (i/1024 >= 256 && i%1024 >= 512 && i%1024 < 768) {
+    chip = 3;
+  } else if (i/1024 >= 256 && i%1024 >= 768 && i%1024 < 1024) {
+    chip = 4;
+  } else if (i/1024 < 256 && i%1024 >= 0 && i%1024 < 256) {
+    chip = 5;
+  } else if (i/1024 < 256 && i%1024 >= 256 && i%1024 < 512) {
+    chip = 6;
+  } else if (i/1024 < 256 && i%1024 >= 512 && i%1024 < 768) {
+    chip = 7;
+  } else if (i/1024 < 256 && i%1024 >= 768 && i%1024 < 1024) {
+    chip = 8;
+  } else {
+    cout << "problem" << endl;
+  }
+  return chip;
+}
+
+
+int supercolumnOfPixel(int i) {
+  int sc = ((i-((i/256)*256))/64) + 1;
+  return sc;
+}
+
+
+int ADCOfPixel(int i) {
+  int adc = -1;
+  if (i%1024 < 256) {
+    adc = 0;
+  } else if (i%1024 < 512) {
+    adc = 1;
+  } else if (i%1024 < 768) {
+    adc = 2;
+  } else {
+    adc = 3;
+  }
+  return adc;
+}
+
+
+void LoadPaletteFalse() {
+
+  // jet
+  Double_t stops[5] = { 0.00, 0.34, 0.61, 0.84, 1.00 };
+  Double_t red[5] = { 0.00, 0.00, 0.87, 1.00, 0.51 };
+  Double_t green[5] = { 0.00, 0.81, 1.00, 0.20, 0.00 };
+  Double_t blue[5] = { 0.51, 1.00, 0.12, 0.00, 0.00 };
+  TColor::CreateGradientColorTable(5, stops, red, green, blue, 255, 1.0);
+  gStyle->SetNumberContours(255);
+
+}
+
+
+void UsePaletteViridis() {
+
+  // viridis 
+  Double_t stops[9] = {0.0000, 0.1250, 0.2500, 0.3750, 0.5000, 0.6250, 0.7500, 0.8750, 1.0000};
+  Double_t red[9] = {26./255., 51./255., 43./255., 33./255., 28./255., 35./255., 74./255., 144./255., 246./255.};
+  Double_t green[9] = {9./255., 24./255., 55./255., 87./255., 118./255., 150./255., 180./255., 200./255., 222./255.};
+  Double_t blue[9] = {30./255., 96./255., 112./255., 114./255., 112./255., 101./255., 72./255., 35./255., 0./255.};
+  TColor::CreateGradientColorTable(9, stops, red, green, blue, 255, 1.0);
+  gStyle->SetNumberContours(255);
+
+}
+
+
+void UsePaletteDBR() {
+
+  // dark body radiator
+  Double_t stops[9] = {0.0000, 0.1250, 0.2500, 0.3750, 0.5000, 0.6250, 0.7500, 0.8750, 1.0000};
+  Double_t red[9] = {0./255., 45./255., 99./255., 156./255., 212./255., 230./255., 237./255., 234./255., 242./255.};
+  Double_t green[9] = {0./255., 0./255., 0./255., 45./255., 101./255., 168./255., 238./255., 238./255., 243./255.};
+  Double_t blue[9] = {0./255., 1./255., 1./255., 3./255., 9./255., 8./255., 11./255., 95./255., 230./255.};
+  TColor::CreateGradientColorTable(9, stops, red, green, blue, 255, 1.0);
+  gStyle->SetNumberContours(255);
+
+}
+
+
+void jungfrauStyle()
+{
+  gROOT->SetStyle("Plain"); /*Default white background for all plots*/
+  /*set bkg color of all to white*/
+  gStyle->SetCanvasColor(kWhite);
+  gStyle->SetFrameFillColor(kWhite);
+  gStyle->SetStatColor(kWhite);
+  gStyle->SetPadColor(kWhite);
+  gStyle->SetFillColor(10);
+  gStyle->SetTitleFillColor(kWhite);
+  
+  /* SetPaperSize wants width & height in cm: A4 is 20,26 & US is 20,24*/
+  gStyle->SetPaperSize(20, 26); 
+  /* No yellow border around histogram*/
+  gStyle->SetDrawBorder(0);
+  /* remove border of canvas*/
+  gStyle->SetCanvasBorderMode(0);
+  /* remove border of pads*/
+  gStyle->SetPadBorderMode(0);
+  gStyle->SetFrameBorderMode(0);
+  gStyle->SetLegendBorderSize(0);
+  
+  /* default text size*/
+  gStyle->SetTextSize(0.05);
+  gStyle->SetTitleSize(0.07,"xyz");
+  gStyle->SetLabelSize(0.06,"xyz");
+  /* title offset: distance between given text and axis, here x,y,z*/
+  gStyle->SetLabelOffset(0.015,"xyz");
+  gStyle->SetTitleOffset(1.1,"yz");
+  gStyle->SetTitleOffset(1.0,"x");
+  
+  /* Use visible font for all text*/
+  int font = 42; 
+  gStyle->SetTitleFont(font);
+  gStyle->SetTitleFontSize(0.06);
+  gStyle->SetStatFont(font);
+  gStyle->SetStatFontSize(0.07);
+  gStyle->SetTextFont(font);
+  gStyle->SetLabelFont(font,"xyz");
+  gStyle->SetTitleFont(font,"xyz");
+  gStyle->SetTitleBorderSize(0);
+  gStyle->SetStatBorderSize(1);
+  gStyle->SetLegendFont(font);
+
+  /* big marker points*/
+  gStyle->SetMarkerStyle(1);
+  gStyle->SetLineWidth(2);  
+  gStyle->SetMarkerSize(1.2);
+  /*set palette in 2d histogram to nice and colorful one*/
+  gStyle->SetPalette(1,0); 
+  
+  /*No title for histograms*/
+  gStyle->SetOptTitle(0);
+  /* show the errors on the stat box */
+  gStyle->SetOptStat(0); 
+  /* show errors on fitted parameters*/
+  gStyle->SetOptFit(0); 
+  /* number of decimals used for errors*/
+  gStyle->SetEndErrorSize(5);   
+  
+  /* set line width to 2 by default so that histograms are visible when printed small
+     idea: emphasize the data, not the frame around*/
+  gStyle->SetHistLineWidth(2);
+  gStyle->SetFrameLineWidth(2);
+  gStyle->SetFuncWidth(2);
+  gStyle->SetHistLineColor(kBlack);
+  gStyle->SetFuncColor(kRed);
+  gStyle->SetLabelColor(kBlack,"xyz");
+  
+  //set the margins
+  gStyle->SetPadBottomMargin(0.2);
+  gStyle->SetPadTopMargin(0.05);
+  gStyle->SetPadRightMargin(0.05);
+  gStyle->SetPadLeftMargin(0.2);
+  
+  //set the number of divisions to show
+  gStyle->SetNdivisions(506, "xy");
+  
+  //turn off xy grids
+  gStyle->SetPadGridX(0);
+  gStyle->SetPadGridY(0);
+  
+  //set the tick mark style
+  gStyle->SetPadTickX(1);
+  gStyle->SetPadTickY(1);
+  
+  gStyle->SetCanvasDefW(800);
+  gStyle->SetCanvasDefH(700);
+  
+  LoadPaletteFalse();
+
+  gROOT->ForceStyle();
+}
+
+// types of pixels
+bool isQuad(int i) {
+  if ( (i/NC == 255 || i/NC == 256) && (i%NC == 255 || i%NC == 256 || i%NC == 511 || i%NC == 512 || i%NC == 767 || i%NC == 768) ) {
+    return true;
+  } else {
+    return false;
+  }
+}
+
+bool isDouble(int i) {
+  if(i/NC == 255 || i/NC == 256 || i%NC == 255 || i%NC == 256 || i%NC == 511 || i%NC == 512 || i%NC == 767 || i%NC == 768) {
+    if (not isQuad(i)) {
+	return true;
+      } else {
+      return false;
+    }
+  } else {
+    return false;
+  }
+}
+
+bool isNextToDouble(int i) {
+  if(i/NC == 254 || i/NC == 257 || i%NC == 254 || i%NC == 257 || i%NC == 510 || i%NC == 513 || i%NC == 766 || i%NC == 769) {
+    if (not isDouble(i) && not isQuad(i)) {
+      return true;
+    } else {
+      return false;
+    }
+  } else {
+    return false;
+  }
+}
+
+bool isEdge(int i) {
+  if (i%NC == 0 || i%NC == 1023 || i/NC == 0 || i/NC == 511) {
+    if (not isDouble(i) && not isNextToDouble(i)) {
+      return true;
+    } else {
+      return false;
+    }
+  } else {
+    return false;
+  }
+}
+
+bool isInnerEdge(int i) {
+  if (i%NC == 1 || i%NC == 1022 || i/NC == 1 || i/NC == 510) {
+    if (not isEdge(i) && not isDouble(i) && not isNextToDouble(i)) {
+      return true;
+    } else {
+      return false;
+    }
+  } else {
+    return false;
+  }
+}
+
+bool isBulk(int i) {
+  if (not isEdge(i) && not isInnerEdge(i) && not isDouble(i) && not isQuad(i) && not isNextToDouble(i)) {
+    return true;
+  } else {
+    return false;
+  }
+}
+
+double highestPointBeforeSwitching(const vector<double> &lower_filter, const vector<double> &higher_filter) {
+
+  double highest_point_before_switching;
+
+  if (higher_filter.size() > 0) {
+    // find the highest value in lower_filter that is lower than all entries in higher_filter
+    highest_point_before_switching = *min_element(lower_filter.begin(),lower_filter.end())-1;
+    
+    double lowest_entry_in_higher_filter = *min_element(higher_filter.begin(),higher_filter.end());
+    for(vector<double>::const_iterator it = lower_filter.begin(); it != lower_filter.end(); ++it) {
+      if (*it < lowest_entry_in_higher_filter) {
+	if (*it > highest_point_before_switching) {
+	  highest_point_before_switching = *it;
+	}
+      }
+    }
+  } else {
+    highest_point_before_switching = *max_element(lower_filter.begin(),lower_filter.end());
+  }
+
+  return highest_point_before_switching;
+
+}
+
+double lowestPointAfterSwitching(const vector<double> &higher_filter, const vector<double> &lower_filter) {
+
+  double lowest_point_after_switching;
+
+  if (lower_filter.size() > 0) {
+    // find the lowest value in higher_filter that is higher than all entries in lower_filter
+    lowest_point_after_switching = *max_element(higher_filter.begin(),higher_filter.end())+1;
+
+    double highest_entry_in_lower_filter = *max_element(lower_filter.begin(),lower_filter.end());
+    for(vector<double>::const_iterator it = higher_filter.begin(); it != higher_filter.end(); ++it) {
+      if (*it > highest_entry_in_lower_filter) {
+	if (*it < lowest_point_after_switching) {
+	  lowest_point_after_switching = *it;
+	}
+      }
+    }
+  } else {
+    lowest_point_after_switching = *min_element(higher_filter.begin(),higher_filter.end());
+  }
+
+  return lowest_point_after_switching;
+
+}
+
+double highestPointBeforeSaturation(const vector<double> &filter, const vector<double> &adc_g2) {
+
+  // first look if saturation ever occurs
+  bool saturation = false;
+  for(vector<double>::const_iterator it_adc = adc_g2.begin(); it_adc != adc_g2.end(); ++it_adc) {
+    if (*it_adc == 0) {
+      saturation = true;
+    }
+  }
+
+  if (saturation == false) {
+
+    double highest_point_before_saturation = *max_element(filter.begin(),filter.end());
+    return highest_point_before_saturation;
+
+  } else {
+
+    // find the lowest filter value which saturates
+    double lowest_point_which_saturates = *max_element(filter.begin(),filter.end());
+
+    vector<double>::const_iterator it_adc = adc_g2.begin();
+
+    for(vector<double>::const_iterator it_filter = filter.begin(); it_filter != filter.end(); ++it_filter) {
+      if (*it_adc == 0) {
+	if (*it_filter < lowest_point_which_saturates) {
+	  lowest_point_which_saturates = *it_filter;
+	}
+      }
+      ++it_adc;
+    }
+  
+    // find the highest filter value which is below saturation
+    double highest_point_before_saturation = *min_element(filter.begin(),filter.end());
+    for(vector<double>::const_iterator it_filter = filter.begin(); it_filter != filter.end(); ++it_filter) {
+      if (*it_filter < lowest_point_which_saturates) {
+	if (*it_filter > highest_point_before_saturation) {
+	  highest_point_before_saturation = *it_filter;
+	}
+      }
+    }
+
+    return highest_point_before_saturation;
+  }
+
+}
+
+#endif /* JUNGFRAUCOMMONFUNCTIONS_H */
diff --git a/sls_detector_calibration/jungfrauCommonHeader.h b/sls_detector_calibration/jungfrauCommonHeader.h
new file mode 100644
index 0000000..6f41a49
--- /dev/null
+++ b/sls_detector_calibration/jungfrauCommonHeader.h
@@ -0,0 +1,14 @@
+#include<iostream>  // std::cout
+
+// LoadPalette
+#include "TColor.h"
+#include "TStyle.h"
+// jungfrauStyle
+#include "TROOT.h"
+
+#define NC 1024
+#define NR 512
+#define NCH NC*NR
+#define NSC 16
+
+using namespace std;
diff --git a/sls_detector_calibration/jungfrauFile.C b/sls_detector_calibration/jungfrauFile.C
new file mode 100644
index 0000000..63fb943
--- /dev/null
+++ b/sls_detector_calibration/jungfrauFile.C
@@ -0,0 +1,150 @@
+//#include "jungfrauCommonHeader.h"
+
+#include <fstream>  // std::ifstream
+#include <stdint.h>  // uint16_t
+#include <ctime>  // clock
+#include <iomanip> // std::setprecision
+
+#include <time.h>
+#include <sys/time.h>
+
+class oneframe {
+
+public:
+  uint64_t framenumber;
+  uint64_t bunchid;
+  uint16_t imagedata[NCH];
+
+};
+
+class jungfrauFile {
+
+ private:
+
+ public:
+  ifstream filebin; // file descriptor
+  oneframe thisframe;
+  long long int filebinSize; // need 64 bit
+  clock_t filebinStart;
+  clock_t filebinStop;
+  struct  timeval tss,tsss; //for timing
+
+  jungfrauFile() { // constructor
+
+    cout << "jungfrauFile constructed" << endl;
+    //cout << "sizeof(thisframe) " << sizeof(thisframe) << endl;
+    //cout << "address of data array " <<  &(thisframe.imagedata)<< endl;
+
+  }
+  
+  uint16_t * getFrameDataHandle() {
+
+    return (uint16_t*)&(thisframe.imagedata);
+
+  }
+
+  bool readNextFrame() {
+ 
+    if (filebin.is_open()) {
+      if (!filebin.eof()) {
+	filebin.read(((char*)&thisframe),sizeof(thisframe));
+	
+	std::cout << std::fixed << std::setprecision(2)
+		  << "\r   [" << std::string(int(filebin.tellg()*10 / float(filebinSize)), '#')
+		  << std::string(10 + 1 - int(filebin.tellg()*10 / float(filebinSize)), ' ') 
+		  << "] " << 100 * filebin.tellg() / float(filebinSize) << "%";
+
+        if (filebin.good()) {return true;} // necessary to stop the frame 'after' the last frame returning true
+	else {
+	  filebinStop = clock();
+	  double elapsed_secs = double(filebinStop - filebinStart) / CLOCKS_PER_SEC;
+	  cout << " end of file reached in " << elapsed_secs << "s" << endl;
+	  return false;
+	}
+
+      } else {
+	filebinStop = clock();
+	double elapsed_secs = double(filebinStop - filebinStart) / CLOCKS_PER_SEC;
+	cout << " end of file reached in " << elapsed_secs << "s" << endl;
+	
+	return false;
+      }
+    } else {
+      cout << "file not open" << endl;
+      return false;
+    }
+
+  }
+
+  int currentFrameNumber() {
+ 
+    return int(thisframe.framenumber);
+
+  }
+
+  uint64_t currentBunchID() {
+ 
+    return uint64_t(thisframe.bunchid);
+
+  }
+
+  /*****************/
+  /* VH 2022-05-05 */
+  /*****************/
+  uint64_t currentSCnumber() {
+
+    return ( ( thisframe.bunchid >> 8 )&0xf );
+
+  }
+  /*****************/
+  /* VH 2022-05-05 */
+  /*****************/
+
+  void rewind() {
+    if (filebin.is_open()) {
+      cout << "rewinding file" << endl;
+      filebin.clear();
+      filebin.seekg(0);
+    } else {
+      cout << "file was not open" << endl;
+    }
+  }
+
+  void open(char *fformat, int fileindex) {
+
+    char fname[1000];
+    sprintf(fname,fformat,fileindex);
+    cout << "file name " << fname << endl;
+    if (filebin.is_open()) {
+      cout << "WARNING: a file is already open" << endl;
+      cout << "WARNING: closing this file" << endl;
+      filebin.close();
+    }
+    cout << "opening file" << endl;
+    filebin.open((const char *)(fname), ios::in | ios::binary);
+    filebin.seekg(0,filebin.end);
+    filebinSize = filebin.tellg();
+    filebin.seekg(0,filebin.beg);
+    filebinStart = clock();
+  }
+
+  void close() {
+
+    if (filebin.is_open()) {
+      filebin.close();
+      cout << "closed file" << endl;
+    } else {
+      cout << "file was not open" << endl;
+    }
+
+  }
+
+  bool isOpen() {
+    if (filebin.is_open()) {
+      return true;
+    } else {
+      return false;
+    }
+  }
+
+};
diff --git a/sls_detector_calibration/jungfrauPedestal.C b/sls_detector_calibration/jungfrauPedestal.C
new file mode 100644
index 0000000..8ee987f
--- /dev/null
+++ b/sls_detector_calibration/jungfrauPedestal.C
@@ -0,0 +1,217 @@
+#include "MovingStat.h"
+
+#include "TCanvas.h"
+#include "TH1F.h"
+#include "TH2F.h"
+
+class jungfrauPedestal {
+
+private:
+  MovingStat stat [NCH];
+
+public:
+  TH1F* pedestals1d;
+  TH2F* pedestals2d;
+  TH1F* pedestalvars1d;
+  TH2F* pedestalvars2d;
+  TH2F* pedestalsvars;
+
+  jungfrauPedestal() {
+    cout << "jungfrauPedestal constructed" << endl;
+  }
+
+  bool addFrameToPedestalCalculation(uint16_t *imagedata) {
+    for (int i = 0; i<NCH; i++) {
+      stat[i].Calc(double(imagedata[i]&0x3fff));
+    }
+
+    return true;
+  }
+
+  bool addG0FrameToPedestalCalculation(uint16_t *imagedata) {
+    for (int i = 0; i<NCH; i++) {
+      if (((imagedata[i]&0xc000) >> 14) == 0) {
+	stat[i].Calc(double(imagedata[i]&0x3fff));
+      }
+    }
+
+    return true;
+  }
+
+  bool addG0PixelToPedestalCalculation(int i, double adc_double) {
+    stat[i].Calc(adc_double);
+
+    return true;
+  }
+
+  bool addG1FrameToPedestalCalculation(uint16_t *imagedata) {
+    for (int i = 0; i<NCH; i++) {
+      //if ((((imagedata[i]&0xc000) >> 14) == 1)) {     
+     if ((((imagedata[i]&0xc000) >> 14) == 1)||(((imagedata[i]&0xc000) >> 14) == 2)) { //to fix FW error during JF11 FW develop. needed for modules 243-379 
+	stat[i].Calc(double(imagedata[i]&0x3fff));
+      }
+    }
+
+    return true;
+  }
+
+  bool addG2FrameToPedestalCalculation(uint16_t *imagedata) {
+    for (int i = 0; i<NCH; i++) {
+      if (((imagedata[i]&0xc000) >> 14) == 3) {
+	stat[i].Calc(double(imagedata[i]&0x3fff));
+      }
+    }
+    
+    return true;
+  }
+  
+  // single threshold
+  bool addG0FrameToPedestalCalculationWThreshold(uint16_t *imagedata, jungfrauPedestal *pedestalObj, uint16_t threshold) {
+    for (int i = 0; i<NCH; i++) {
+      if (((imagedata[i]&0xc000) >> 14) == 0) {
+  	if (fabs((imagedata[i]&0x3fff) - pedestalObj->pedestalOfChannel(i)) < threshold) {
+  	  stat[i].Calc(double(imagedata[i]&0x3fff));
+  	}
+      }
+    }
+    
+    return true;
+  }
+  
+  // array of thresholds
+  bool addG0FrameToPedestalCalculationWThreshold(uint16_t *imagedata, jungfrauPedestal *pedestalObj, double *thresholds) {
+    for (int i = 0; i<NCH; i++) {
+      if (((imagedata[i]&0xc000) >> 14) == 0) {
+	if (fabs((imagedata[i]&0x3fff) - pedestalObj->pedestalOfChannel(i)) < thresholds[i]) { 
+	  stat[i].Calc(double(imagedata[i]&0x3fff));
+	}
+      }
+    }
+    
+    return true;
+  }
+
+  double pedestalOfChannel(int ichannel) {
+    return stat[ichannel].Mean();
+  }
+
+  double semOfChannel(int ichannel) {
+    return stat[ichannel].StandardDeviation() / sqrt(stat[ichannel].NumDataValues());
+  }
+
+  double rmsOfChannel(int ichannel) {
+    return stat[ichannel].StandardDeviation();
+  }
+
+  double varianceOfChannel(int ichannel) {
+    return stat[ichannel].Variance();
+  }
+
+  bool pedestalData(double* pedestals) {
+    for (int i = 0; i < NCH; i++) {
+      pedestals[i] = stat[i].Mean();
+    }
+    return true;
+  }
+
+  bool pedestalData(uint16_t* pedestals) {
+    for (int i = 0; i < NCH; i++) {
+      pedestals[i] = uint16_t(round(stat[i].Mean()));
+    }
+    return true;
+  }
+
+  bool pedestalRMSData(double* pedestalRMSs) {
+    for (int i = 0; i < NCH; i++) {
+      pedestalRMSs[i] = stat[i].StandardDeviation();
+    }
+    return true;
+  }
+
+  int pedestalNFrames() {
+    return stat[0].GetN();
+  }
+
+  int pedestalUpdates(int ichannel) {
+    return stat[ichannel].GetNUpdates();
+  }
+
+  void pedestalResetUpdates() {
+    for (int i=0; i < NCH; i++){
+      stat[i].ResetUpdates();
+    }
+  }
+
+  void pedestalClear() {
+    for (int i=0; i < NCH; i++){
+      stat[i].Clear();
+    }
+  }
+
+  bool pedestalSetNFrames(int nframes) {
+    for (int i=0; i < NCH; i++){
+      stat[i].SetN(nframes);
+    }
+    return true;
+  }
+
+  void savePedestals(uint16_t* pedestals, string fileName) {
+    ofstream pedefile (fileName.c_str());
+    if (pedefile.is_open()) {
+      cout << "saving pedestals" << endl;
+      for (int i=0; i<NCH; i++){
+        pedefile << pedestals[i] << " " ;
+      }
+      pedefile.close();
+    } else {
+      cout << "Unable to open file" << endl;
+    }
+  }
+
+  void savePedestals(double* pedestals, string fileName) {
+    ofstream pedefile (fileName.c_str());
+    if (pedefile.is_open()) {
+      cout << "saving pedestals" << endl;
+      for (int i=0; i<NCH; i++){
+        pedefile << pedestals[i] << " " ;
+      }
+      pedefile.close();
+    } else {
+      cout << "Unable to open file" << endl;
+    }
+  }
+
+  bool readPedestals(uint16_t* pedestals, string fileName) {
+    ifstream pedefile(fileName.c_str());
+    if (pedefile.is_open()) {
+      cout << "reading pedestal file" << endl;
+      for (int i = 0; i < NCH; i++) {
+	pedefile >> pedestals[i];
+      }
+      return true;
+    } else {
+      cout << "Unable to open file" << endl;
+      return false;
+    }
+  }
+
+  bool readPedestals(double* pedestals, string fileName) {
+    ifstream pedefile(fileName.c_str());
+    if (pedefile.is_open()) {
+      cout << "reading pedestal file" << endl;
+      for (int i = 0; i < NCH; i++) {
+	pedefile >> pedestals[i];
+      }
+      return true;
+    } else {
+      cout << "Unable to open file" << endl;
+      return false;
+    }
+  }
+
+  bool setPedestalOfChannel(int ichannel, uint16_t set_pede) {
+    stat[ichannel].Calc(double(set_pede));
+    return true;
+  }
+
+};
diff --git a/sls_detector_calibration/jungfrauPixelMask.C b/sls_detector_calibration/jungfrauPixelMask.C
new file mode 100755
index 0000000..a66ecc7
--- /dev/null
+++ b/sls_detector_calibration/jungfrauPixelMask.C
@@ -0,0 +1,245 @@
+
+#include <algorithm> // std::count
+#include "TCanvas.h"
+#include "TH2F.h"
+
+class jungfrauPixelMask {
+
+private:
+
+public:
+  TH2F* pixel_mask2d;
+  TCanvas* maskcanvas;
+
+  jungfrauPixelMask() {
+    cout << "jungfrauPixelMask constructed" << endl;
+  }
+
+  void initialisePixelMask(bool* mask) {
+    std::fill_n(mask, NCH, true); // initialise with true (not masked)
+    cout << "jungfrauPixelMask initialised" << endl;
+  }
+
+  int getNMasked(bool* mask) {
+    return count(mask, mask + NCH, false);
+  }
+
+  void maskPixel(int pixel, bool *mask) {
+    mask[pixel]= false;
+  }
+
+  void maskChip(int chip, bool* mask) {
+    if (chip < 1 || chip > 8){
+      cout << "expected chip number 1-8" << endl;
+    }
+    for (int i=0; i<NCH; i++) {
+      if (chipOfPixel(i) == chip) {
+	mask[i] = false;
+      }
+    }
+  }
+  
+  void maskSupercolumn(int chip, int sc, bool* mask) {
+    if (chip < 1 || chip > 8 || sc < 1 || sc > 4){
+      cout << "expected chip number 1-8 and supercolumn number 1-4" << endl;
+    }
+    for (int i=0; i<NCH; i++) {
+      if (chipOfPixel(i) == chip){
+	if (supercolumnOfPixel(i) == sc) {
+	  mask[i] = false;
+	}
+      }
+    }
+  }
+
+  void maskColumn(int chip, int c, bool* mask) {
+    if (chip < 1 || chip > 8 || c < 0 || c >= NC){
+      cout << "expected chip number 1-8 and column number 0--(NC-1)" << endl;
+    }
+    for (int i=0; i<NCH; i++) {
+      if (chipOfPixel(i) == chip){
+	if (i%NC == c){
+	  mask[i] = false;
+	}
+      }
+    }
+  }
+
+  void maskChipEdges(bool* mask) {
+    for (int i=0; i<NCH; i++) {
+      if (i%1024 == 0 ||
+	  i%1024 == 255 ||
+	  i%1024 == 256 ||
+	  i%1024 == 511 ||
+	  i%1024 == 512 ||
+	  i%1024 == 767 ||
+	  i%1024 == 768 ||
+	  i%1024 == 1023 ||
+	  i/1024 == 0 ||
+	  i/1024 == 255 ||
+	  i/1024 == 256 ||
+	  i/1024 == 511) {
+	mask[i] = false;
+      }
+    }
+  }
+
+  void maskNeighbours(bool* mask) {
+    bool tmp_mask [NCH];
+    std::fill_n(tmp_mask, NCH, true);
+
+    for (int i=0; i<NCH; i++) {
+      if (mask[i] == false) {
+	tmp_mask[i] = false;
+
+	if (i<1024) { // bottom
+	  if (i==0) { // bottom left corner
+	    tmp_mask[i+1] = false;
+	    tmp_mask[i+1024] = false;
+	    tmp_mask[i+1024+1] = false;
+	  } else if (i==1023) { // bottom right corner
+	    tmp_mask[i-1] = false;
+	    tmp_mask[i+1024] = false;
+	    tmp_mask[i+1024-1] = false;
+	  } else { // bottom centre
+	    tmp_mask[i-1] = false;
+	    tmp_mask[i+1]  = false;
+	    tmp_mask[i+1024-1] = false;
+	    tmp_mask[i+1024] = false;
+	    tmp_mask[i+1024+1] = false;
+	  }
+	} else if (i>=523264){ // top
+	  if (i==523264){ // top left corner
+	    tmp_mask[i+1] = false;
+	    tmp_mask[i-1024] = false;
+	    tmp_mask[i-1024+1] = false;
+	  } else if (i==524287){ // top right corner
+	    tmp_mask[i-1] = false;
+	    tmp_mask[i-1024] = false;
+	    tmp_mask[i-1024-1] = false;
+	  } else { // top center
+	    tmp_mask[i+1] = false;
+	    tmp_mask[i-1] = false;
+	    tmp_mask[i-1024+1] = false;
+	    tmp_mask[i-1024] = false;
+	    tmp_mask[i-1024-1] = false;
+	  }
+	} else if ((i%1024)==0) { // left hand side center
+	  tmp_mask[i+1] = false;
+	  tmp_mask[i+1024] = false;
+	  tmp_mask[i+1024+1] = false;
+	  tmp_mask[i-1024] = false;
+	  tmp_mask[i-1024+1] = false;
+	} else if ((i%1024)==1023) { // right hand side center
+	  tmp_mask[i-1] = false;
+	  tmp_mask[i+1024] = false;
+	  tmp_mask[i+1024-1] = false;
+	  tmp_mask[i-1024] = false;
+	  tmp_mask[i-1024-1] = false;
+	} else { // center
+	  tmp_mask[i-1] = false;
+	  tmp_mask[i+1] = false;
+	  tmp_mask[i+1024-1] = false;
+	  tmp_mask[i+1024] = false;
+	  tmp_mask[i+1024+1] = false;
+	  tmp_mask[i-1024-1] = false;
+	  tmp_mask[i-1024] = false;
+	  tmp_mask[i-1024+1] = false;
+	}
+      }
+    }
+    
+    for (int i = 0; i < NCH; i++) {
+      mask[i] = tmp_mask[i];
+    }
+
+  }
+
+  void maskIfGainNot(int expected, uint16_t *imagedata, bool* mask) {
+    // if gain of a pixel is not as expected, mask pixel
+    for (int i = 0; i<NCH; i++) {
+      if (((imagedata[i]&0xc000) >> 14) != expected) {
+	mask[i] = false;
+	//cout << "gain is " << ((imagedata[i]&0xc000) >> 14) << " not " << expected << " as expected so masking" << endl;
+      }
+    }
+  }
+
+  void maskFromPedestal(bool* mask, double* pedestals, double* pedestalRMSs) {
+    // if pedestal RMS is outside limits, mask pixel
+    for (int i = 0; i<NCH; i++) {
+      if (pedestalRMSs[i] > 100 || pedestals[i] < 1000 || pedestals[i] > 4000) {
+	  mask[i] = false;
+	}
+    }
+  }
+
+  void maskIfPedeRMSG0GreaterThan(bool* mask, double* pedestalG0RMSs, int val) {
+    for (int i = 0; i < NCH; i++) {
+      if (pedestalG0RMSs[i] > val) {
+	mask[i] = false;
+      }
+    }
+  }
+
+  void maskIfPedestalStep(uint16_t *imagedata, uint16_t* last_frame, bool* mask, int step) {
+    for (int i = 0; i < NCH; i++) {
+      if ( fabs(last_frame[i] - (imagedata[i]&0x3fff)) > step) {
+  	mask[i] = false;
+      }
+    }
+
+  }
+
+  void plotPixelMask(bool* mask, string fileName) {
+    pixel_mask2d = new TH2F("pixel_mask2d","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);
+    for (int i = 0; i < NCH; i++) {
+      if (mask[i] == false) {
+	pixel_mask2d->Fill(i%NC,i/NC,1);
+      }
+    }
+
+    maskcanvas = new TCanvas("maskcanvas","",150,10,800,400);
+    maskcanvas->SetLeftMargin(0.1);
+    maskcanvas->SetRightMargin(0.13);
+    maskcanvas->SetTopMargin(0.08);
+    maskcanvas->SetBottomMargin(0.15);
+
+    pixel_mask2d->GetXaxis()->SetTitle("Column");
+    pixel_mask2d->GetYaxis()->SetTitle("Row");
+    pixel_mask2d->GetYaxis()->SetTitleOffset(0.7);
+    pixel_mask2d->Draw("colz");
+    maskcanvas->SaveAs(fileName.c_str());
+    maskcanvas->Close();
+    pixel_mask2d->Delete();
+  }
+
+  void saveMask(bool* mask, string fileName) {
+    ofstream maskfile (fileName.c_str());
+    if (maskfile.is_open()) {
+      cout << "saving pixel mask" << endl;
+      for (int i=0; i<NCH; i++){
+        maskfile << mask[i] << " " ;
+      }
+      maskfile.close();
+    } else {
+      cout << "Unable to open file" << endl;
+    }
+  }
+
+  bool readMask(bool* mask, string fileName) {
+    ifstream maskfile(fileName.c_str());
+    if (maskfile.is_open()) {
+      cout << "reading pixel mask" << endl;
+      for (int i = 0; i < NCH; ++i) {
+	maskfile >> mask[i];
+      }
+      return true;
+    } else {
+      cout << "Unable to open pixel mask" << endl;
+      return false;
+    }
+  }
+
+
+};