Classes, Libraries and Functions for Characterization and Calibration of SLS Detectors. First import.

git-svn-id: file:///afs/psi.ch/project/sls_det_software/svn/slsDetectorCalibration@1 113b152e-814d-439b-b186-022a431db7b5

slsDetectorCalibration/MovingStat.h

@@ -0,0 +1,88 @@
+  class MovingStat
+    {
+    public:
+    MovingStat() : m_n(0), n(0) {}
+
+        void Clear()
+        {
+            m_n = 0;
+        }
+	
+	void SetN(int i) {n=i;};
+	int GetN() {return n;};
+	void Calc(double x) {
+	  if (m_n<n) Add(x);
+	  else Push(x);
+	}
+
+	void Add(double x) {
+	  m_n++;
+
+            // See Knuth TAOCP vol 2, 3rd edition, page 232
+	  if (m_n == 1)
+            {
+                m_oldM = m_newM = x;
+                m_oldM2 = x*x;
+            }
+            else
+            {
+                m_newM = m_oldM + x;
+		m_newM2 = m_oldM2 + x*x;
+		//m_newM2 = m_oldM2 + (x*x - m_oldM*m_oldM)/m_n;
+
+                // set up for next iteration
+                m_oldM = m_newM; 
+                m_oldM2 = m_newM2;
+            }
+
+	}
+
+
+        void Push(double x)
+        {
+	  if (m_n == 1)
+            {
+	      m_oldM = m_newM = x;
+	      m_oldM2 = x*x;
+            }
+            else
+            {
+                m_newM = m_oldM + (x - m_oldM/m_n);
+		m_newM2 = m_oldM2 + (x*x - m_oldM2/m_n);
+		//m_newM2 = m_oldM2 + (x*x/m_n - m_oldM*m_oldM/(m_n*m_n));
+                // set up for next iteration
+                m_oldM = m_newM; 
+                m_oldM2 = m_newM2;
+            }
+
+        }
+
+        int NumDataValues() const
+        {
+            return m_n;
+        }
+
+        double Mean() const
+        {
+            return (m_n > 0) ? m_newM/m_n : 0.0;
+        }
+	double M2() const
+        {
+            return ( (m_n > 1) ? m_newM2/m_n : 0.0 );
+        }
+
+        double Variance() const
+        {
+	  return ( (m_n > 1) ? (M2()-Mean()*Mean()) : 0.0 );
+        }
+
+        double StandardDeviation() const
+        {
+	  return ( (Variance() > 0) ? sqrt( Variance() ) : -1 );
+        }
+
+    private:
+	int n;
+        int m_n;
+        double m_oldM, m_newM, m_oldM2, m_newM2;
+    };

slsDetectorCalibration/RunningStat.h

@@ -0,0 +1,55 @@
+  class RunningStat
+    {
+    public:
+        RunningStat() : m_n(0) {}
+
+        void Clear()
+        {
+            m_n = 0;
+        }
+
+        void Push(double x)
+        {
+            m_n++;
+
+            // See Knuth TAOCP vol 2, 3rd edition, page 232
+            if (m_n == 1)
+            {
+                m_oldM = m_newM = x;
+                m_oldS = 0.0;
+            }
+            else
+            {
+                m_newM = m_oldM + (x - m_oldM)/m_n;
+                m_newS = m_oldS + (x - m_oldM)*(x - m_newM);
+    
+                // set up for next iteration
+                m_oldM = m_newM; 
+                m_oldS = m_newS;
+            }
+        }
+
+        int NumDataValues() const
+        {
+            return m_n;
+        }
+
+        double Mean() const
+        {
+            return (m_n > 0) ? m_newM : 0.0;
+        }
+
+        double Variance() const
+        {
+            return ( (m_n > 1) ? m_newS/(m_n - 1) : 0.0 );
+        }
+
+        double StandardDeviation() const
+        {
+            return sqrt( Variance() );
+        }
+
+    private:
+        int m_n;
+        double m_oldM, m_newM, m_oldS, m_newS;
+    };

slsDetectorCalibration/energyCalibration.cpp

@@ -0,0 +1,457 @@
+#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::gaussChargeSharing(Double_t *x, Double_t *par) { 
+  Double_t f, arg=0;
+  if (par[3]!=0) arg=(x[0]-par[2])/par[3];
+    f=par[4]*TMath::Exp(-1*arg*arg/2.); 
+    f=f+par[5]*(par[4]/2*(TMath::Erfc(sign*arg/(TMath::Sqrt(2.)))))+par[0]-par[1]*sign*x[0]; 
+    return f;                                       
+}
+
+// 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::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),
+					 fscurve(NULL),
+					 fspectrum(NULL),
+#endif
+					 funcs(NULL),
+					 plot_flag(1),
+					 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");
+
+#endif
+
+
+}
+
+
+
+void energyCalibration::fixParameter(int ip, Double_t val){
+  
+  fscurve->FixParameter(ip, val);
+  fspectrum->FixParameter(ip, val);
+}
+
+
+void energyCalibration::releaseParameter(int ip){
+
+  fscurve->ReleaseParameter(ip);
+  fspectrum->ReleaseParameter(ip);
+}
+
+
+
+
+
+
+
+energyCalibration::~energyCalibration(){ 
+#ifdef MYROOT
+  delete fscurve;
+  delete fspectrum;
+#endif
+  
+}
+
+#ifdef MYROOT
+
+
+
+TH1F* energyCalibration::createMedianHistogram(TH2F* h2, int ch0, int nch) {
+
+  if (h2==NULL || nch==0)
+    return NULL;
+
+  double *x=new double[nch];
+  TH1F *h1=NULL;
+  
+  double val=-1;
+
+  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);
+  }
+  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::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);
+    } else {
+      fscurve->FixParameter(5,0);
+      fspectrum->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);
+
+}
+
+
+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,"R");
+  } 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::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);
+}
+
+
+
+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
+
+
+
+

slsDetectorCalibration/energyCalibration.h

@@ -0,0 +1,412 @@
+#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
+
+
+#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
+      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);
+
+  /** 
+      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) ;
+
+  /** 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 fractional height of the charge sharing pedestal (scales with par[3]
+      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
+  */
+  Double_t spectrum(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
+
+  static TH1F* createMedianHistogram(TH2F* h2, int ch0, int nch);
+
+
+  /** 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);
+  
+/** 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);
+
+
+  /**
+     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);
+
+
+  /**
+     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);  
+
+#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 */
+
+
+  TF1 *fscurve; /**< function with which the s-curve will be fitted */
+
+  TF1 *fspectrum; /**< 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
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+

slsDetectorCalibration/energyCalibration_cpp.d

@@ -0,0 +1,6 @@
+
+# DO NOT DELETE
+
+./energyCalibration_cpp.so: energyCalibration.h
+./energyCalibration_cpp.so: /afs/psi.ch/project/sls_det_software/root_v5.32.01_sl5_32bit/include/cintdictversion.h /afs/psi.ch/project/sls_det_software/root_v5.32.01_sl5_32bit/include/RVersion.h
+energyCalibration_cpp__ROOTBUILDVERSION= 5.32/01