improved library structure for analog detectors and interpolation t.b.t.

slsDetectorCalibration/Makefile.moench

@@ -0,0 +1,20 @@
+
+CBFLIBDIR=/afs/psi.ch/project/sls_det_software/CBFlib-0.9.5
+LIBRARYCBF=$(CBFLIBDIR)/lib/*.o
+INCDIR=-IslsDetectorCalibration -I../slsReceiverSoftware/include -I$(CBFLIBDIR)/include/ -I. -IetaVEL
+LIBHDF5=-L$(CBFLIBDIR)/lib/ -lhdf5
+LDFLAG= -L/usr/lib64/ -lpthread
+#-L../../bin 
+MAIN=moench03OnTheFlyAnalysis.C
+
+#DESTDIR?=../bin
+
+all: moench03OnTheFlyAnalysis
+
+
+
+moench03OnTheFlyAnalysis:  $(MAIN)  $(INCS)
+			   g++ -o moench03OnTheFlyAnalysis  $(MAIN)  -lm -ltiff -lstdc++ $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF)
+
+clean:
+	rm -f moench03OnTheFlyAnalysis

slsDetectorCalibration/analogDetector.h

@@ -0,0 +1,254 @@
+#ifndef ANALOGDETECTOR_H
+#define  ANALOGDETECTOR_H
+
+
+#include "slsDetectorData.h"
+#include "pedestalSubtraction.h"
+#include "commonModeSubtraction.h"
+
+
+template <class dataType> class analogDetector {
+
+  /** @short class to perform pedestal subtraction etc. for an analog detector */
+
+ public:
+
+
+  /**
+
+     Constructor (no error checking if datasize and offsets are compatible!)
+     \param d detector data structure to be used
+     \param csize cluster size (should be an odd number). Defaults to 3
+     \param nsigma number of rms to discriminate from the noise. Defaults to 5
+     \param sign 1 if photons are positive, -1 if  negative
+     \param cm common mode subtraction algorithm, if any. Defaults to NULL i.e. none
+     \param nped number of samples for pedestal averaging
+     \param nd number of dark frames to average as pedestals without photon discrimination at the beginning of the measurement
+
+
+  */
+  
+
+ analogDetector(slsDetectorData<dataType> *d, int sign=1, 
+		       commonModeSubtraction *cm=NULL, int nnx=-1, int nny=-1) : det(d), nx(nnx), ny(nny), stat(NULL), cmSub(cm),  iframe(-1), dataSign(sign){
+    
+    if (det)
+      det->getDetectorSize(nx,ny);
+    
+    stat=new pedestalSubtraction*[ny];
+    for (int i=0; i<ny; i++) {
+      stat[i]=new pedestalSubtraction[nx];
+    }
+   
+    
+  };
+  /**
+     destructor. Deletes the cluster structure and the pdestalSubtraction array
+  */
+  virtual ~analogDetector() {for (int i=0; i<ny; i++) delete [] stat[i]; delete [] stat;};
+    
+    
+  /** resets the pedestalSubtraction array and the commonModeSubtraction */
+  void newDataSet(){iframe=-1; for (int iy=0; iy<ny; iy++) for (int ix=0; ix<nx; ix++) stat[iy][ix].Clear(); if (cmSub) cmSub->Clear(); };  
+  
+  /** resets the eventMask to undefined and the commonModeSubtraction */
+  void newFrame(){iframe++; if (cmSub) cmSub->newFrame();};
+  
+
+    /** sets the commonModeSubtraction algorithm to be used 
+	\param cm commonModeSubtraction algorithm to be used (NULL unsets) 
+	\returns pointer to the actual common mode subtraction algorithm
+    */
+    commonModeSubtraction *setCommonModeSubtraction(commonModeSubtraction *cm) {cmSub=cm; return cmSub;};
+
+
+    /**
+       sets the sign of the data
+       \param sign 1 means positive values for photons, -1 negative, 0 gets
+       \returns current sign for the data
+    */
+    int setDataSign(int sign=0) {if (sign==1 || sign==-1) dataSign=sign; return dataSign;};
+
+    
+    /**
+       adds value to pedestal (and common mode) for the given pixel
+       \param val value to be added
+       \param ix pixel x coordinate
+       \param iy pixel y coordinate
+    */
+    virtual void addToPedestal(double val, int ix, int iy=0){ 
+      if (ix>=0 && ix<nx && iy>=0 && iy<ny) {
+	stat[iy][ix].addToPedestal(val); 
+	if (cmSub)  {
+	  if (det) if (det->isGood(ix, iy)==0) return;
+	  cmSub->addToCommonMode(val, ix, iy);
+	};
+      };
+    }
+    
+    /**
+       gets  pedestal (and common mode)
+       \param ix pixel x coordinate
+       \param iy pixel y coordinate
+       \param cm 0 (default) without common mode subtraction, 1 with common mode subtraction (if defined)
+    */
+    virtual double getPedestal(int ix, int iy, int cm=0){if (ix>=0 && ix<nx && iy>=0 && iy<ny) if (cmSub && cm>0) return stat[iy][ix].getPedestal()-cmSub->getCommonMode(); else return stat[iy][ix].getPedestal(); else return -1;};
+
+    /**
+       gets  pedestal rms (i.e. noise)
+       \param ix pixel x coordinate
+       \param iy pixel y coordinate
+    */
+    double getPedestalRMS(int ix, int iy){if (ix>=0 && ix<nx && iy>=0 && iy<ny) return stat[iy][ix].getPedestalRMS();else return -1;};
+
+    
+    /**
+       sets  pedestal
+       \param ix pixel x coordinate
+       \param iy pixel y coordinate
+       \param val value to set
+    */
+    virtual void setPedestal(int ix, int iy, double val, double rms=0){if (ix>=0 && ix<nx && iy>=0 && iy<ny) stat[iy][ix].setPedestal(val,rms);};
+    
+    
+    
+    
+    
+    virtual void addToPedestal(char *data) {
+      
+      
+      newFrame();
+      
+      
+      for (int ix=0; ix<nx; ix++) {
+	for (int iy=0; iy<ny; iy++) {
+	  
+	  
+	  addToPedestal(data,ix,iy);
+	  
+	  
+	}
+      }
+      return ;
+      
+    };
+    virtual void addToPedestal(char *data, int ix, int iy=0) {
+      
+      
+      double val;
+      if (ix>=0 && ix<nx && iy>=0 && iy<ny) {
+	if (det)
+	  val=dataSign*det->getValue(data, ix, iy);
+	else
+	  val=((double*)data)[iy*nx+ix];
+	
+	addToPedestal(val,ix,iy);
+      }
+      return  ;
+      
+    };
+      
+ 
+        
+   virtual  double *subtractPedestal(char *data, double *val=NULL) {
+      
+      newFrame();
+
+      if (val==NULL)
+	val=new double[nx*ny];
+
+      for (int ix=0; ix<nx; ix++) {
+	for (int iy=0; iy<ny; iy++) {
+	  val[iy*nx+ix]+=subtractPedestal(data, ix, iy);
+	}
+      }
+      return  val;
+    };
+      
+   
+
+
+
+   virtual  double subtractPedestal(char *data, int ix, int iy=0) {
+      if (ix>=0 && ix<nx && iy>=0 && iy<ny) {
+	if (det)
+	  return dataSign*det->getValue(data, ix, iy)-getPedestal(ix,iy);
+	else
+	  return ((double*)data)[iy*nx+ix]-getPedestal(ix,iy);
+      }
+   };
+      
+
+
+   virtual  int getNPhotons(char *data, int ix, int iy=0, int thr=-1) {
+      
+      if (ix>=0 && ix<nx && iy>=0 && iy<ny) {
+	if (thr<=0) thr=-1*thr*getPedestalRMS(ix,iy);
+	
+	if (det)
+	  return (dataSign*det->getValue(data, ix, iy)-getPedestal(ix,iy))/thr;
+	else
+	  return (((double*)data)[(iy)*nx+ix]-getPedestal(ix,iy))/thr;
+	
+      }
+      return 0;
+      
+   };
+
+   int *getNPhotons(char *data, double thr=-1, int *nph=NULL) {
+
+      double val;
+      if (nph==NULL)
+	nph=new int[nx*ny];   
+      double tthr=thr;
+      newFrame();
+      for (int ix=0; ix<nx; ix++) {
+	for (int iy=0; iy<ny; iy++) {
+	  nph[iy*nx+ix]+=getNPhotons(data, ix, iy, thr);
+	}
+      }
+      return nph;
+
+    }
+
+
+
+    /** sets/gets number of samples for moving average pedestal calculation
+	\param i number of samples to be set (0 or negative gets)
+	\returns actual number of samples
+    */
+    int SetNPedestals(int i=-1) {
+      int ix=0, iy=0; 
+      if (i>0) 
+	for (ix=0; ix<nx; ix++) 
+	  for (iy=0; iy<ny; iy++) 
+	    stat[iy][ix].SetNPedestals(i); 
+      return stat[0][0].SetNPedestals();
+    };
+
+ 
+    double getROI(char *data, int xmin, int xmax, int ymin=0, int ymax=1) {
+      double val=0;
+      for (int ix=xmin; ix<xmax; ix++)
+	for (int iy=ymin; iy<ymax; iy++)
+	  if (ix>=0 && ix<nx && iy>=0 && iy<ny) 
+	    val+=subtractPedestal(data, ix, iy);
+      return val;
+      
+    }
+
+
+
+ protected:
+  
+    slsDetectorData<dataType> *det; /**< slsDetectorData to be used */
+    int nx; /**< Size of the detector in x direction */
+    int ny; /**< Size of the detector in y direction */
+    pedestalSubtraction **stat; /**< pedestalSubtraction class */
+    commonModeSubtraction *cmSub;/**< commonModeSubtraction class */
+    int dataSign; /**< sign of the data i.e. 1 if photon is positive, -1 if negative */
+    int iframe;  /**< frame number (not from file but incremented within the dataset every time newFrame is called */
+    
+};
+
+#endif

slsDetectorCalibration/etaVEL/etaInterpolationBase.h

@@ -0,0 +1,272 @@
+#ifndef ETA_INTERPOLATION_BASE_H
+#define ETA_INTERPOLATION_BASE_H
+
+#ifdef MYROOT1
+#include <TObject.h>
+#include <TTree.h>
+#include <TH2D.h>
+#include <TH2F.h>
+#endif
+
+#include "slsInterpolation.h"
+
+class etaInterpolationBase : public slsInterpolation {
+  
+ public:
+ etaInterpolationBase(int nx=400, int ny=400, int ns=25, int nb=-1, double emin=1, double emax=0) : slsInterpolation(nx,ny,ns), hhx(NULL), hhy(NULL), heta(NULL),nbeta(nb),etamin(emin), etamax(emax) {
+    if (nb<=0) 
+      nbeta=nSubPixels*10;   
+    if (etamin>=etamax) {
+      etamin=-0.1;
+      etamax=1.1;
+    }
+    etastep=(etamax-etamin)/nbeta;
+#ifdef MYROOT1
+    heta=new TH2D("heta","heta",nbeta,etamin,etamax,nbeta,etamin,etamax);
+    hhx=new TH2D("hhx","hhx",nbeta,etamin,etamax,nbeta,etamin,etamax);
+    hhy=new TH2D("hhy","hhy",nbeta,etamin,etamax,nbeta,etamin,etamax);
+#endif
+#ifndef MYROOT1
+    heta=new int[nbeta*nbeta];
+    hhx=new int[nbeta*nbeta];
+    hhy=new int[nbeta*nbeta];
+    
+#endif
+    
+  };
+  
+#ifdef MYROOT1
+  TH2D *setEta(TH2D *h, int nb=-1, double emin=1, double emax=0)
+  {  
+     if (h) { heta=h;
+       nbeta=heta->GetNbinsX();
+       etamin=heta->GetXaxis()->GetXmin();
+       etamax=heta->GetXaxis()->GetXmax();
+       etastep=(etamax-etamin)/nbeta;
+     }
+    return heta;
+  };
+  TH2D *getFlatField(){return setEta(NULL);};
+#endif
+  
+#ifndef MYROOT1
+  int *setEta(int *h, int nb=-1, double emin=1, double emax=0)
+  {  
+    if (h) {heta=h;
+      nbeta=nb;
+      if (nb<=0) nbeta=nSubPixels*10;
+      etamin=emin;
+      etamax=emax;
+      if (etamin>=etamax) {
+	etamin=-0.1;
+	etamax=1.1;
+      }
+      etastep=(etamax-etamin)/nbeta;
+    }
+    return heta;
+  };
+  int *getFlatField(){return setEta(NULL);};
+#endif
+  
+
+  
+virtual void prepareInterpolation(int &ok)=0;
+ 
+  /* ////////////////////////////////////////////////////////////////////////////// */
+
+#ifdef MYROOT1
+TH2D *gethhx()
+  {
+    hhx->Scale((double)nSubPixels);
+    return hhx;
+  };
+  
+  TH2D *gethhy()
+  {
+    hhy->Scale((double)nSubPixels);
+    return hhy;
+    };
+#endif
+    
+#ifndef MYROOT1
+int *gethhx()
+  {
+    // hhx->Scale((double)nSubPixels);
+    return hhx;
+  };
+  
+  int *gethhy()
+  {
+    // hhy->Scale((double)nSubPixels);
+    return hhy;
+    };
+#endif
+    
+  //////////////////////////////////////////////////////////////////////////////
+  //////////// /*It return position hit for the event in input */ //////////////
+  virtual void getInterpolatedPosition(int x, int y, double *data, double &int_x, double &int_y)
+  {
+    double sDum[2][2];
+    double tot, totquad;
+    double etax,etay;
+    
+    int corner;
+    corner=calcQuad(data, tot, totquad, sDum); 
+    
+    calcEta(totquad, sDum, etax, etay); 
+    getInterpolatedPosition(x,y,etax,etay,corner,int_x,int_y);
+    return;
+  };
+  
+
+  virtual void getInterpolatedPosition(int x, int y, double etax, double etay, int corner, double &int_x, double &int_y)
+  {
+
+
+    double xpos_eta,ypos_eta;
+    double dX,dY;
+    double ex,ey;
+    switch (corner)
+      {
+      case TOP_LEFT:
+	dX=-1.; 
+	dY=+1.; 
+	break;
+      case TOP_RIGHT:
+	dX=+1.; 
+	dY=+1.; 
+	break;
+      case BOTTOM_LEFT:
+	dX=-1.; 
+	dY=-1.; 
+	break;
+      case BOTTOM_RIGHT:
+	dX=+1.; 
+	dY=-1.; 
+	break;
+      default:
+	dX=0.; 
+	dY=0.;
+      }
+    
+
+
+#ifdef MYROOT1
+    xpos_eta=(hhx->GetBinContent(hhx->GetXaxis()->FindBin(etax),hhy->GetYaxis()->FindBin(etay)))/((double)nSubPixels);
+    ypos_eta=(hhy->GetBinContent(hhx->GetXaxis()->FindBin(etax),hhy->GetYaxis()->FindBin(etay)))/((double)nSubPixels);
+#endif
+#ifndef MYROOT1
+    ex=(etax-etamin)/etastep;
+    ey=(etay-etamin)/etastep;
+    if (ex<0) ex=0;
+    if (ex>=nSubPixels) ex=nSubPixels-1;
+    if (ey<0) ey=0;
+    if (ey>=nSubPixels) ey=nSubPixels-1;
+    
+   
+    xpos_eta=(((double)hhx[(int)(ey*nbeta+ex)]))/((double)nSubPixels);
+    ypos_eta=(((double)hhy[(int)(ey*nbeta+ex)]))/((double)nSubPixels);
+      //else
+      //return 0;
+
+#endif
+    
+    int_x=((double)x) + 0.5*dX + xpos_eta;
+    int_y=((double)y) + 0.5*dY + ypos_eta;
+    //return 1;
+
+  }
+
+
+
+
+
+
+  /////////////////////////////////////////////////////////////////////////////////////////////////
+   virtual void getPositionETA3(int x, int y, double *data, double &int_x, double &int_y)
+   {
+    double sDum[2][2];
+    double tot, totquad;
+    double eta3x,eta3y;
+    double ex,ey;
+    
+    calcQuad(data, tot, totquad, sDum); 
+    calcEta3(data,eta3x, eta3y,tot); 
+    
+    double xpos_eta,ypos_eta;
+    
+#ifdef MYROOT1
+    xpos_eta=((hhx->GetBinContent(hhx->GetXaxis()->FindBin(eta3x),hhy->GetYaxis()->FindBin(eta3y))))/((double)nSubPixels);
+    ypos_eta=((hhy->GetBinContent(hhx->GetXaxis()->FindBin(eta3x),hhy->GetYaxis()->FindBin(eta3y))))/((double)nSubPixels);
+    
+#endif
+#ifndef MYROOT1
+    ex=(eta3x-etamin)/etastep;
+    ey=(eta3y-etamin)/etastep;
+    
+    if (ex<0) ex=0;
+    if (ex>=nSubPixels) ex=nSubPixels-1;
+    if (ey<0) ey=0;
+    if (ey>=nSubPixels) ey=nSubPixels-1;
+      xpos_eta=(((double)hhx[(int)(ey*nbeta+ex)]))/((double)nSubPixels);
+      ypos_eta=(((double)hhy[(int)(ey*nbeta+ex)]))/((double)nSubPixels);
+#endif
+    
+    int_x=((double)x) + xpos_eta;
+    int_y=((double)y) + ypos_eta;
+    
+    return;
+  };
+  
+  //////////////////////////////////////////////////////////////////////////////////////
+  virtual int addToFlatField(double *cluster, double &etax, double &etay){
+    double sDum[2][2];
+    double tot, totquad;
+    int corner;
+    corner=calcQuad(cluster, tot, totquad, sDum); 
+    
+    double xpos_eta,ypos_eta;
+    double dX,dY;
+  
+    
+    calcEta(totquad, sDum, etax, etay); 
+   
+    return addToFlatField(etax,etay);
+
+    };
+
+
+  virtual int addToFlatField(double etax, double etay){
+
+    int ex,ey; 
+
+#ifdef MYROOT1
+    heta->Fill(etax,etay);
+#endif
+#ifndef MYROOT1
+    ex=(etax-etamin)/etastep;
+    ey=(etay-etamin)/etastep;
+    if (ey<nbeta && ex<nbeta && ex>=0 && ey>=0)
+      heta[ey*nbeta+ex]++;    
+#endif
+    return 0;    
+};
+  
+ protected:
+  
+#ifdef MYROOT1
+  TH2D *heta;
+  TH2D *hhx;
+  TH2D *hhy;
+#endif
+#ifndef MYROOT1
+  int *heta;
+  int *hhx;
+  int *hhy;
+#endif
+  int nbeta;
+  double etamin, etamax, etastep;
+
+};
+
+#endif

slsDetectorCalibration/etaVEL/etaInterpolationGlobal.h

@@ -0,0 +1,85 @@
+#ifndef ETA_INTERPOLATION_GLOBAL_H
+#define ETA_INTERPOLATION_GLOBAL_H
+
+
+#include "etaInterpolationBase.h"
+
+class etaInterpolationGlobal : public etaInterpolationBase{
+ public:
+ globalEtaInterpolation(int nx=400, int ny=400, int ns=25, int nb=-1, double emin=1, double emax=0) : etaInterpolationBase(nx,ny,ns, nb, emin,emax){};
+  
+
+
+  virtual void prepareInterpolation(int &ok)
+  {
+   ok=1;  
+#ifdef MYROOT1   
+   if (hhx) delete hhx;
+   if (hhy) delete hhy;
+
+   hhx=new TH2D("hhx","hhx",heta->GetNbinsX(),heta->GetXaxis()->GetXmin(),heta->GetXaxis()->GetXmax(), heta->GetNbinsY(),heta->GetYaxis()->GetXmin(),heta->GetYaxis()->GetXmax()); 
+   hhy=new TH2D("hhy","hhy",heta->GetNbinsX(),heta->GetXaxis()->GetXmin(),heta->GetXaxis()->GetXmax(), heta->GetNbinsY(),heta->GetYaxis()->GetXmin(),heta->GetYaxis()->GetXmax()); 
+   
+#endif
+
+
+   ///*Eta Distribution Rebinning*///
+   double bsize=1./nSubPixels; //precision
+   // cout<<"nPixelsX = "<<nPixelsX<<" nPixelsY = "<<nPixelsY<<" nSubPixels = "<<nSubPixels<<endl;
+   double tot_eta=0;
+   for (int ip=0; ip<nbeta*nbeta; ip++)
+     tot_eta+=heta[ip];
+   cout << "total eta entries is :"<< tot_eta << endl;   
+   if (tot_eta<=0) {ok=0; return;};
+
+
+   double hx[nbeta]; //projection x
+   double hy[nbeta]; //projection y
+
+   for (int ibx=0; ibx<nbeta; ibx++) {
+     for (int iby=0; iby<nbeta; iby++) {
+       hx[ibx]=hx[ibx]+heta[ibx+iby*nbeta];
+       hy[iby]=hx[iby]+heta[ibx+iby*nbeta];
+     }
+   }
+  double hix[nbeta]; //integral of projection x
+  double hiy[nbeta]; //integral of projection y
+  hix[0]=hx[0];
+  hiy[0]=hy[0];
+  
+  for (int ib=1; ib<nbeta; ib++) {
+    hix[ib]=hix[ib-1]+hx[ib];
+    hiy[ib]=hiy[ib-1]+hx[ib];
+  }
+
+  int ib=0;
+   for (int ibx=0; ibx<nbeta; ibx++) {
+     if (hix[ibx]>(ib+1)*tot_eta*bsize) ib++;
+     for (int iby=0; iby<nbeta; iby++) {
+#ifdef MYROOT1  
+       hhx->SetBinContent(ibx+1,iby+1,ib);
+#endif
+#ifndef MYROOT1  
+       hhx[ibx+iby*nbeta]=ib; 
+#endif
+     }
+   }
+   ib=0;
+   for (int iby=0; iby<nbeta; iby++) {
+     if (hiy[iby]>(ib+1)*tot_eta*bsize) ib++;
+     for (int ibx=0; ibx<nbeta; ibx++) {
+#ifdef MYROOT1 
+       hhy->SetBinContent(ibx+1,iby+1,ib); 
+#endif
+#ifndef MYROOT1  
+       hhy[ibx+iby*nbeta]=ib; 
+#endif
+     }
+   }
+
+   return ;
+  };
+
+};
+
+#endif

slsDetectorCalibration/etaVEL/etaInterpolationPosXY.h

@@ -0,0 +1,99 @@
+#ifndef ETA_INTERPOLATION_POSXY_H
+#define ETA_INTERPOLATION_POSXY_H
+
+
+#include "etaInterpolationBase.h"
+
+class etaInterpolationPosXY : public etaInterpolationBase{
+ public:
+ etaInterpolationPosXY(int nx=400, int ny=400, int ns=25, int nb=-1, double emin=1, double emax=0) : etaInterpolationBase(nx,ny,ns, nb, emin,emax){};
+  
+
+
+  virtual void prepareInterpolation(int &ok)
+  {
+   ok=1;  
+#ifdef MYROOT1   
+   if (hhx) delete hhx;
+   if (hhy) delete hhy;
+
+   hhx=new TH2D("hhx","hhx",heta->GetNbinsX(),heta->GetXaxis()->GetXmin(),heta->GetXaxis()->GetXmax(), heta->GetNbinsY(),heta->GetYaxis()->GetXmin(),heta->GetYaxis()->GetXmax()); 
+   hhy=new TH2D("hhy","hhy",heta->GetNbinsX(),heta->GetXaxis()->GetXmin(),heta->GetXaxis()->GetXmax(), heta->GetNbinsY(),heta->GetYaxis()->GetXmin(),heta->GetYaxis()->GetXmax()); 
+   
+#endif
+
+
+   ///*Eta Distribution Rebinning*///
+   double bsize=1./nSubPixels; //precision
+   // cout<<"nPixelsX = "<<nPixelsX<<" nPixelsY = "<<nPixelsY<<" nSubPixels = "<<nSubPixels<<endl;
+   double tot_eta=0;
+   double tot_eta_x=0;
+   double tot_eta_y=0;
+   for (int ip=0; ip<nbeta*nbeta; ip++)
+     tot_eta+=heta[ip];
+   cout << "total eta entries is :"<< tot_eta << endl;   
+   if (tot_eta<=0) {ok=0; return;};
+
+
+   double hx[nbeta]; //profile x
+   double hy[nbeta]; //profile y
+   double hix[nbeta]; //integral of projection x
+   double hiy[nbeta]; //integral of projection y
+   int ii=0;
+  
+   for (int ib=0; ib<nbeta; ib++) {
+
+     tot_eta_x=0;
+     tot_eta_y=0;
+
+     for (int iby=0; iby<nbeta; iby++) {
+       hx[iby]=heta[iby+ib*nbeta];
+       tot_eta_x+=hx[iby];
+       hy[iby]=heta[ib+iby*nbeta];
+       tot_eta_y+=hy[iby];
+     }
+
+     hix[0]=hx[0];
+     hiy[0]=hy[0];
+
+     for (int iby=1; iby<nbeta; iby++) {
+       hix[iby]=hix[iby-1]+hx[iby];
+       hiy[iby]=hiy[iby-1]+hx[iby];
+     }
+
+     ii=0;
+     
+     for (int ibx=0; ibx<nbeta; ibx++) {
+       if (hix[ibx]>(ii+1)*tot_eta_x*bsize) ii++;
+#ifdef MYROOT1  
+       hhx->SetBinContent(ibx+1,ib+1,ii);
+#endif
+#ifndef MYROOT1  
+       hhx[ibx+ib*nbeta]=ii; 
+#endif
+     }
+     
+     ii=0;
+     
+     for (int ibx=0; ibx<nbeta; ibx++) {
+       if (hiy[ibx]>(ii+1)*tot_eta_y*bsize) ii++;
+#ifdef MYROOT1  
+       hhy->SetBinContent(ib+1,ibx+1,ii);
+#endif
+#ifndef MYROOT1  
+       hhy[ib+ibx*nbeta]=ii; 
+#endif
+     }
+     
+
+
+
+     
+   }
+  
+   return ;
+  }
+
+};
+
+#endif

slsDetectorCalibration/etaVEL/linearInterpolation.h

@@ -0,0 +1,111 @@
+#ifndef LINEAR_INTERPOLATION_H
+#define LINEAR_INTERPOLATION_H
+
+//#include <TObject.h>
+//#include <TTree.h>
+//#include <TH2F.h>
+
+#include "slsInterpolation.h"
+
+class linearInterpolation : public slsInterpolation{
+
+ public:
+ linearInterpolation(int nx=400, int ny=400, int ns=25) : slsInterpolation(nx,ny,ns) {};
+ 
+  virtual void prepareInterpolation(int &ok){ok=1;};
+
+   //////////////////////////////////////////////////////////////////////////////
+  //////////// /*It return position hit for the event in input */ //////////////
+  virtual void getInterpolatedPosition(Int_t x, Int_t y, double *data, double &int_x, double &int_y)
+  {
+    double sDum[2][2];
+    double tot, totquad;
+    double etax,etay;
+    
+    int corner;
+    corner=calcQuad(data, tot, totquad, sDum); 
+    
+    calcEta(totquad, sDum, etax, etay); 
+    getInterpolatedPosition(x, y, etax,etay, corner, int_x, int_y);
+    
+    return;
+  };
+
+
+  virtual void getInterpolatedPosition(Int_t x, Int_t y, double etax, double etay, int corner, double &int_x, double &int_y)
+  {
+   
+    double xpos_eta,ypos_eta;
+    double dX,dY;
+    switch (corner)
+      {
+      case TOP_LEFT:
+	dX=-1.; 
+	dY=+1.; 
+	break;
+      case TOP_RIGHT:
+	dX=+1.; 
+	dY=+1.; 
+	break;
+      case BOTTOM_LEFT:
+	dX=-1.; 
+	dY=-1.; 
+	break;
+      case BOTTOM_RIGHT:
+	dX=+1.; 
+	dY=-1.; 
+	break;
+      default:
+	dX=0.; 
+	dY=0.;
+      }
+    
+    
+    xpos_eta=(etax);
+    ypos_eta=(etay);
+    
+    int_x=((double)x) + 0.5*dX + xpos_eta;
+    int_y=((double)y) + 0.5*dY + ypos_eta;
+    
+    return;
+  };
+
+
+
+
+
+
+
+  ////////////////////////////////////////////////////////////////////////////////////////////////////////
+  virtual void getPositionETA3(Int_t x, Int_t y, double *data, double &int_x, double &int_y)
+   {
+    double sDum[2][2];
+    double tot, totquad;
+    double eta3x,eta3y;
+    
+    calcQuad(data, tot, totquad, sDum); 
+    calcEta3(data,eta3x, eta3y,tot); 
+    
+    double xpos_eta,ypos_eta;
+    
+    xpos_eta=eta3x;
+    ypos_eta=eta3y;
+    
+    int_x=((double)x) + xpos_eta;
+    int_y=((double)y) + ypos_eta;
+    
+    return;
+  };
+  ////////////////////////////////////////////////////////////////////////////////////////////////////////
+ 
+  virtual int addToFlatField(double *cluster, double &etax, double &etay){};
+  virtual int addToFlatField(double etax, double etay){};
+  
+
+ protected:
+  ;
+
+
+};
+
+#endif

slsDetectorCalibration/etaVEL/noInterpolation.h

@@ -0,0 +1,58 @@
+#ifndef NO_INTERPOLATION_H
+#define NO_INTERPOLATION_H
+
+/* #ifdef MYROOT1 */
+/* #include <TObject.h> */
+/* #include <TTree.h> */
+/* #include <TH2F.h> */
+/* #include <TROOT.h> */
+/* #include <TRandom.h> */
+/* #endif */
+
+#include  <cstdlib>
+#include "slsInterpolation.h"
+
+
+
+class noInterpolation : public slsInterpolation{
+ public:
+ noInterpolation(int nx=400, int ny=400, int ns=25) : slsInterpolation(nx,ny,ns) {};// {eventGenerator=new TRandom();};
+  virtual void prepareInterpolation(int &ok){ok=1;};
+  
+  //////////////////////////////////////////////////////////////////////////////
+  //////////// /*It return position hit for the event in input */ //////////////
+
+  virtual void getInterpolatedPosition(Int_t x, Int_t y, double *data, double &int_x, double &int_y)
+  {
+    //Random coordinate in the Pixel reference
+    int_x = x + ((double)rand())/((double)RAND_MAX) -0.5;//eventGenerator->Uniform(-0.5,0.5);
+    int_y = y + ((double)rand())/((double)RAND_MAX) -0.5;//eventGenerator->Uniform(-0.5,0.5);
+   
+    return ;
+  };
+  virtual void getInterpolatedPosition(Int_t x, Int_t y, double etax, double etay, int corner, double &int_x, double &int_y)
+  {
+    return getInterpolatedPosition(x, y, NULL, int_x, int_y);
+  };
+  
+  //////////////////////////////////////////////////////////////////////////////////////
+  virtual void getPositionETA3(Int_t x, Int_t y, double *data, double &int_x, double &int_y)
+  {
+    //Random coordinate in the Pixel reference
+    int_x = x + ((double)rand())/((double)RAND_MAX) -0.5;//eventGenerator->Uniform(-0.5,0.5);
+    int_y = y + ((double)rand())/((double)RAND_MAX) -0.5;//eventGenerator->Uniform(-0.5,0.5);
+   
+    return ;
+  };
+  //////////////////////////////////////////////////////////////////////////////////////
+  virtual int addToFlatField(double *cluster, double &etax, double &etay){return 0;};
+  virtual int addToFlatField(double etax, double etay){return 0;};
+  
+ protected:
+  ;
+  // TRandom *eventGenerator;
+  //  ClassDefNV(slsInterpolation,1);
+  // #pragma link C++ class slsInterpolation-;
+};
+
+#endif

slsDetectorCalibration/etaVEL/slsInterpolation.h

@@ -1,9 +1,11 @@
 #ifndef SLS_INTERPOLATION_H
 #define SLS_INTERPOLATION_H
 
+#ifdef MYROOT1
 #include <TObject.h>
 #include <TTree.h>
 #include <TH2F.h>
+#endif
 
 #ifndef DEF_QUAD
 #define DEF_QUAD
@@ -16,35 +18,76 @@
   };
 #endif
 
-class slsInterpolation : public TObject{
+//#ifdef MYROOT1
+//: public TObject
+//#endif
+class slsInterpolation 
+{
 
  public:
- slsInterpolation(int nx=40, int ny=160, int ns=25) :nPixelsX(nx), nPixelsY(ny),  nSubPixels(ns) {hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny);};
+ slsInterpolation(int nx=400, int ny=400, int ns=25) :nPixelsX(nx), nPixelsY(ny),  nSubPixels(ns) {
- 
+   
+#ifdef MYROOT1
+hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny);
+#endif
+   
+#ifndef MYROOT1
+ hint=new int[ns*nx*ns*ny];
+#endif
+
+};
+  
+  int getNSubPixels() {return nSubPixels;};
+  
+  
   //create eta distribution, eta rebinnining etc.
   //returns flat field image
   virtual void prepareInterpolation(int &ok)=0;
 
   //create interpolated image
   //returns interpolated image
+#ifdef MYROOT1
   virtual TH2F *getInterpolatedImage(){return hint;};
+#endif
+
+#ifndef MYROOT1
+  virtual int *getInterpolatedImage(){return hint;};
+#endif
   //return position inside the pixel for the given photon
-  virtual void getInterpolatedPosition(Int_t x, Int_t y, Double_t *data, Double_t &int_x, Double_t &int_y)=0;
+  virtual void getInterpolatedPosition(int x, int y, double *data, double &int_x, double &int_y)=0;
+  //return position inside the pixel for the given photon
+  virtual void getInterpolatedPosition(int x, int y, double etax, double etay, int quad, double &int_x, double &int_y)=0;
 
-  TH2F *addToImage(Double_t int_x, Double_t int_y){hint->Fill(int_x, int_y); return hint;};
+#ifdef MYROOT1
+  TH2F *addToImage(double int_x, double int_y){hint->Fill(int_x, int_y); return hint;};
+#endif
+
+#ifndef MYROOT1
+  virtual int *addToImage(double int_x, double int_y){ int iy=nSubPixels*int_y; int ix=nSubPixels*int_x; 
+    if (ix>=0 && ix<(nPixelsX*nSubPixels) && iy<(nSubPixels*nPixelsY) && iy>=0 )(*(hint+ix+iy*nPixelsX))+=1; 
+    return hint;
+  };
+#endif
 
 
-
+  virtual int addToFlatField(double *cluster, double &etax, double &etay)=0;
-  virtual int addToFlatField(Double_t *cluster, Double_t &etax, Double_t &etay)=0;
+  virtual int addToFlatField(double etax, double etay)=0;
-  virtual int addToFlatField(Double_t etax, Double_t etay)=0;
   
+#ifdef MYROOT1
+  virtual TH2D *getFlatField(){return NULL;};
+#endif
+
+#ifndef MYROOT1
+  virtual int *getFlatField(){return NULL;};
+#endif
+
   //virtual void Streamer(TBuffer &b);
 
   
-  static int calcQuad(Double_t *cl, Double_t &sum, Double_t &totquad, Double_t sDum[2][2]){
+  static int calcQuad(double *cl, double &sum, double &totquad, double sDum[2][2]){
     
     int corner = UNDEFINED_QUADRANT;
-    Double_t *cluster[3];
+    double *cluster[3];
     cluster[0]=cl;
     cluster[1]=cl+3;
     cluster[2]=cl+6;
@@ -95,8 +138,8 @@ class slsInterpolation : public TObject{
     
   } 
 
-  static int calcEta(Double_t totquad, Double_t sDum[2][2], Double_t &etax, Double_t &etay){
+  static int calcEta(double totquad, double sDum[2][2], double &etax, double &etay){
-    Double_t t,r;
+    double t,r;
     
     if (totquad>0) {
       t = sDum[1][0] + sDum[1][1];
@@ -109,7 +152,7 @@ class slsInterpolation : public TObject{
   } 
 
 
-  static int calcEta(Double_t *cl, Double_t &etax, Double_t &etay, Double_t &sum, Double_t &totquad, Double_t sDum[2][2]) {
+  static int calcEta(double *cl, double &etax, double &etay, double &sum, double &totquad, double sDum[2][2]) {
     int corner = calcQuad(cl,sum,totquad,sDum);
     calcEta(totquad, sDum, etax, etay);
     
@@ -117,8 +160,8 @@ class slsInterpolation : public TObject{
   }
 
 
-  static int calcEtaL(Double_t totquad, int corner, Double_t sDum[2][2], Double_t &etax, Double_t &etay){
+  static int calcEtaL(double totquad, int corner, double sDum[2][2], double &etax, double &etay){
-    Double_t t,r, toth, totv;
+    double t,r, toth, totv;
     if (totquad>0) {
       switch(corner) {
       case TOP_LEFT:
@@ -156,7 +199,7 @@ class slsInterpolation : public TObject{
     return 0;
   }
 
-  static int calcEtaL(Double_t *cl, Double_t &etax, Double_t &etay, Double_t &sum, Double_t &totquad, Double_t sDum[2][2]) {
+  static int calcEtaL(double *cl, double &etax, double &etay, double &sum, double &totquad, double sDum[2][2]) {
     int corner = calcQuad(cl,sum,totquad,sDum);
     calcEtaL(totquad, corner, sDum, etax, etay);
     
@@ -165,7 +208,7 @@ class slsInterpolation : public TObject{
 
 
 
-  static int calcEtaC3(Double_t *cl, Double_t &etax, Double_t &etay, Double_t &sum, Double_t &totquad, Double_t sDum[2][2]){
+  static int calcEtaC3(double *cl, double &etax, double &etay, double &sum, double &totquad, double sDum[2][2]){
     
     int corner = calcQuad(cl,sum,totquad,sDum);
     calcEta(sum, sDum, etax, etay);
@@ -175,8 +218,8 @@ class slsInterpolation : public TObject{
 
 
 
-  static int calcEta3(Double_t *cl, Double_t &etax, Double_t &etay, Double_t &sum) {
+  static int calcEta3(double *cl, double &etax, double &etay, double &sum) {
-    Double_t l,r,t,b;
+    double l,r,t,b;
     sum=cl[0]+cl[1]+cl[2]+cl[3]+cl[4]+cl[5]+cl[6]+cl[7]+cl[8];
     if (sum>0) {
       l=cl[0]+cl[3]+cl[6];
@@ -192,8 +235,8 @@ class slsInterpolation : public TObject{
 
 
 
-  static int calcEta3X(Double_t *cl, Double_t &etax, Double_t &etay, Double_t &sum) {
+  static int calcEta3X(double *cl, double &etax, double &etay, double &sum) {
-    Double_t l,r,t,b;
+    double l,r,t,b;
     sum=cl[0]+cl[1]+cl[2]+cl[3]+cl[4]+cl[5]+cl[6]+cl[7]+cl[8];
     if (sum>0) {
       l=cl[3];
@@ -213,12 +256,14 @@ class slsInterpolation : public TObject{
 
  protected:
   int nPixelsX, nPixelsY;
-  int nSubPixels;
+  int nSubPixels; 
+#ifdef MYROOT1
   TH2F *hint;
+#endif
+#ifndef MYROOT1
+  int *hint;
+#endif
 
-
-  //  ClassDefNV(slsInterpolation,1);
-  // #pragma link C++ class slsInterpolation-;
 };
 
 #endif

slsDetectorCalibration/interpolatingDetector.h

@@ -0,0 +1,274 @@
+#ifndef INTERPOLATINGDETECTOR_H
+#define  INTERPOLATINGDETECTOR_H
+
+
+#include "singlePhotonDetector.h"
+
+#include "slsInterpolation.h"
+
+
+#ifdef MYROOT1
+#include <TTree.h>
+
+#endif
+
+
+#include <iostream>
+
+using namespace std;
+
+
+class interpolatingDetector : public singlePhotonDetector {
+
+  /** @short class to perform pedestal subtraction etc. and find single photon clusters for an analog detector */
+
+ public:
+
+
+  /**
+
+     Constructor (no error checking if datasize and offsets are compatible!)
+     \param d detector data structure to be used
+     \param csize cluster size (should be an odd number). Defaults to 3
+     \param nsigma number of rms to discriminate from the noise. Defaults to 5
+     \param sign 1 if photons are positive, -1 if  negative
+     \param cm common mode subtraction algorithm, if any. Defaults to NULL i.e. none
+     \param nped number of samples for pedestal averaging
+     \param nd number of dark frames to average as pedestals without photon discrimination at the beginning of the measurement
+
+
+  */
+  
+  
+ interpolatingDetector(slsDetectorData<uint16_t> *d, slsInterpolation *inte,
+		       double nsigma=5,  
+		       int sign=1, 
+		       commonModeSubtraction *cm=NULL,
+		       int nped=1000, 
+		       int nd=100, int nnx=-1, int nny=-1) : 
+  singlePhotonDetector(d, 3,nsigma,sign, cm, nped, nd, nnx, nny) , interp(inte)  {};
+  
+  
+  
+#ifdef MYROOT1
+  virtual TH2F *addToInterpolatedImage(char *data, single_photon_hit *clusters, int &nph)
+#endif
+#ifndef MYROOT1
+  virtual int *addToInterpolatedImage(char *data, single_photon_hit *clusters, int &nph)
+#endif
+  {
+    nph=addFrame(data,clusters,0); 
+    if (interp)
+      return interp->getInterpolatedImage();
+    else
+      return NULL;
+  };
+  
+
+#ifdef MYROOT1
+  virtual TH2F *addToFlatField(char *data, single_photon_hit *clusters, int &nph)
+#endif
+#ifndef MYROOT1
+    virtual int *addToFlatField(char *data, single_photon_hit *clusters, int &nph)
+#endif
+  {
+    nph=addFrame(data,clusters,1); 
+    if (interp)
+      return interp->getFlatField();
+    else
+      return NULL;
+  };
+  
+  
+  
+  
+  int addFrame(char *data, single_photon_hit *clusters, int ff=0) {
+    
+    
+    
+    double int_x,int_y, eta_x, eta_y;
+    int nph=0;
+    double val[ny][nx];
+    int cy=(clusterSizeY+1)/2;
+    int cs=(clusterSize+1)/2;
+    int ir, ic;
+    double cc[2][2];
+    double max=0, tl=0, tr=0, bl=0,br=0, *v, vv;
+    int xoff,yoff;
+    int skip=0;
+    if (iframe<nDark) {
+      //cout << iframe << "+"<< nDark <<endl;
+      addToPedestal(data);
+      return 0;
+    }
+
+    newFrame();
+    for (int ix=1; ix<nx-1; ix++) {
+      for (int iy=1; iy<ny-1; iy++) {
+	skip=0;
+	max=0;
+	tl=0;
+	tr=0;
+	bl=0;
+	br=0;
+	tot=0;
+	quadTot=0;
+	quad=UNDEFINED_QUADRANT;
+      
+      
+      
+      
+      
+	(clusters+nph)->rms=getPedestalRMS(ix,iy);
+
+	//	cout << iframe << " " << nph << " " << ix << " " << iy << endl;
+	
+      for (int ir=-(clusterSizeY/2); ir<(clusterSizeY/2)+1; ir++) {
+	for (int ic=-(clusterSize/2); ic<(clusterSize/2)+1; ic++) {
+	  
+
+	  //  if ((iy+ir)>=0 && (iy+ir)<ny && (ix+ic)>=0 && (ix+ic)<nx) {
+	    
+	    if ((iy>1 && ir<=0)) {
+	      ;
+	    } else if (ix>1 && ic<=0) {
+	      ;
+	    } else {
+	      val[iy+ir][ix+ic]=subtractPedestal(data,ix+ic,iy+ir);
+	      eventMask[iy+ir][ix+ic]=PEDESTAL; 
+	    }
+
+	    /* if ((iy+ir)>=iy  && (ix+ic)>=ix) { */
+	    /*   val[iy+ir][ix+ic]=subtractPedestal(data,ix+ic,iy+ir); */
+	    /*   eventMask[iy+ir][ix+ic]=PEDESTAL;  */
+	    /* }  */
+	  
+	  //   cout << ir << " " << ic << " " << val[iy+ir][ix+ic] << endl;
+	    v=&(val[iy+ir][ix+ic]);
+	
+	    if (ir==0 && ic==0) {
+	      if (*v<-nSigma*(clusters+nph)->rms) {
+		eventMask[iy][ix]=NEGATIVE_PEDESTAL;
+		//	cout << "neg ped" << endl;
+	      }
+	    } 
+
+	    // if (skip==0) {
+	      tot+=*v;
+	      if (ir<=0 && ic<=0)
+		bl+=*v;
+	      if (ir<=0 && ic>=0)
+		br+=*v;
+	      if (ir>=0 && ic<=0)
+		tl+=*v;
+	      if (ir>=0 && ic>=0)
+		tr+=*v;
+	      if (*v>max) {
+		max=*v;
+	      }
+	      // }
+	    
+	    //  } else skip=1;
+	}
+      }
+	
+	if (bl>=br && bl>=tl && bl>=tr) {
+	  (clusters+nph)->quad=BOTTOM_LEFT;
+	  (clusters+nph)->quadTot=bl;
+	  xoff=0;
+	  yoff=0;
+	} else if (br>=bl && br>=tl && br>=tr) {
+	  (clusters+nph)->quad=BOTTOM_RIGHT;
+	  (clusters+nph)->quadTot=br;
+	  xoff=1;
+	  yoff=0;
+	} else if (tl>=br && tl>=bl && tl>=tr) {
+	  (clusters+nph)->quad=TOP_LEFT;
+	  (clusters+nph)->quadTot=tl;
+	  xoff=0;
+	  yoff=1;
+	} else if   (tr>=bl && tr>=tl && tr>=br) {
+	  (clusters+nph)->quad=TOP_RIGHT;
+	  (clusters+nph)->quadTot=tr;
+	  xoff=1;
+	  yoff=1;
+	} 
+	
+	if (max>nSigma*(clusters+nph)->rms || tot>sqrt(clusterSizeY*clusterSize)*nSigma*(clusters+nph)->rms || ((clusters+nph)->quadTot)>sqrt(cy*cs)*nSigma*(clusters+nph)->rms) {
+	  if (val[iy][ix]>=max) {
+	    // cout << "max" << endl;
+	    eventMask[iy][ix]=PHOTON_MAX;
+	    (clusters+nph)->tot=tot;
+	    (clusters+nph)->x=ix;
+	    (clusters+nph)->y=iy;
+	    (clusters+nph)->ped=getPedestal(ix,iy, 0);
+	    //	  cout << iframe << " " << ix  << " " << iy << " "<< (clusters+nph)->tot << " " << (clusters+nph)->quadTot << " " << (clusters+nph)->ped<< " " << (clusters+nph)->rms << endl; 
+	    for (int ir=-(clusterSizeY/2); ir<(clusterSizeY/2)+1; ir++) {
+	      for (int ic=-(clusterSize/2); ic<(clusterSize/2)+1; ic++) {
+		if ((iy+ir)>=0 && (iy+ir)<ny && (ix+ic)>=0 && (ix+ic)<nx) {
+		  (clusters+nph)->set_data(val[iy+ir][ix+ic],ic,ir);
+		  cout << val[iy+ir][ix+ic] << " " ;
+		}
+	      }
+	      cout << endl << " " ;
+	    }
+	    
+	     cout << endl << " " ;
+	    
+	  cc[0][0]=(clusters+nph)->get_data(-1+xoff,-1+yoff);
+	  cc[1][0]=(clusters+nph)->get_data(-1+xoff,0+yoff);
+	  cc[0][1]=(clusters+nph)->get_data(0+xoff,-1+yoff);
+	  cc[1][1]=(clusters+nph)->get_data(0+xoff,0+yoff);
+	  
+	  cout << cc[0][0] << " " << cc[0][1] << endl;
+	  cout << cc[1][0] << " " << cc[1][1] << endl;
+	  
+
+	  cout << endl << " " ;
+	  
+	if (interp) {
+	  interp->calcEta((clusters+nph)->quadTot,cc,eta_x, eta_y);
+	  cout << eta_x << " " << eta_y << endl;
+	  cout << "**************************************************************************"<< endl;
+	  // cout << "eta" << endl;
+	  if (ff) {
+	    interp->addToFlatField(eta_x,eta_y);
+	  } else {
+	    //   cout << "interp" << endl;
+	    interp->getInterpolatedPosition(ix,iy,eta_x,eta_y,(clusters+nph)->quad,int_x,int_y);
+	    // cout << "add" << endl;
+	    interp->addToImage(int_x, int_y);
+	  }
+	  // cout << "done" << endl;
+	}	      
+	nph++;
+	  } else {
+	    //   cout << "ph" << endl;
+	    eventMask[iy][ix]=PHOTON;
+	  }
+      } else if (eventMask[iy][ix]==PEDESTAL) {
+	  //  cout << "ped" << endl;
+	addToPedestal(data,ix,iy);
+      }
+	
+      }
+    }
+    
+    return  nph;
+    
+  };
+
+  
+  
+ protected:
+  
+  
+  slsInterpolation *interp;
+
+};
+
+
+
+
+
+#endif

slsDetectorCalibration/moench03OnTheFlyAnalysis.C

@@ -0,0 +1,130 @@
+#include <vector>
+#include <string>
+#include <sstream>
+#include <iomanip>
+#include <fstream>
+#include <stdio.h>
+//#include <deque>
+//#include <list>
+//#include <queue>
+#include <fstream>
+
+#include "moench03Ctb10GbT1Data.h"
+
+#include "interpolatingDetector.h"
+#include "etaInterpolationPosXY.h"
+
+
+#include <ctime>
+
+#define NC 400
+#define NR 400
+
+#include "tiffIO.h"
+
+
+void *moenchProcessFrame() {
+  char fname[10000];
+  strcpy(fname,"/mnt/moench_data/m03-15_mufocustube/plant_40kV_10uA/m03-15_100V_g4hg_300us_dtf_0.raw");
+  
+
+  int nph, nph1;
+  single_photon_hit clusters[NR*NC];
+  // cout << "hits "<< endl;
+  int etabins=250;
+  double etamin=-0.1, etamax=1.1;
+  int nsubpix=5;
+  float *etah=new float[etabins*etabins];
+  // cout << "etah "<< endl;
+  cout << "image size "<< nsubpix*nsubpix*NC*NR << endl;
+  float *image=new float[nsubpix*nsubpix*NC*NR];
+  int *heta, *himage;
+
+  moench03Ctb10GbT1Data *decoder=new  moench03Ctb10GbT1Data();
+  // cout << "decoder "<< endl;
+  etaInterpolationPosXY *interp=new etaInterpolationPosXY(NC, NR, nsubpix, etabins, etamin, etamax);
+  // cout << "interp "<< endl;
+
+  interpolatingDetector *filter=new interpolatingDetector(decoder,interp, 5, 1, 0, 10000, 10000);
+  cout << "filter "<< endl;
+  
+
+  char *buff;
+  int nf=0;
+  int ok=0;
+  ifstream filebin;
+  std::time_t end_time;
+ 
+  int iFrame=-1;
+  int np=-1;
+
+  
+  filter->newDataSet();
+
+
+  nph=0;
+  nph1=0;
+  int iph;
+
+  cout << "file name " << fname << endl;
+  filebin.open((const char *)(fname), ios::in | ios::binary);
+  if (filebin.is_open())
+     cout << "Opened file " << fname<< endl;
+  else
+    cout << "Could not open file " << fname<< endl;
+    while ((buff=decoder->readNextFrame(filebin, iFrame)) && nf<5E4) {
+      if (nf<1.1E4)
+	filter->addToFlatField(buff,clusters,nph1);
+      else {
+	if (ok==0) {
+	  cout << "**************************************************************************"<< endl;
+	  heta=interp->getFlatField();
+	  for (int ii=0; ii<etabins*etabins; ii++) {
+	    etah[ii]=(float)heta[ii];
+	  }
+
+
+	  std::time(&end_time);
+	  cout << std::ctime(&end_time) << " " << nf << " " << nph1 << " " << nph << endl;
+	  WriteToTiff(etah, "/scratch/eta.tiff", etabins, etabins);
+	  
+	  interp->prepareInterpolation(ok);
+	  cout << "**************************************************************************"<< endl;
+	  std::time(&end_time);
+	  cout << std::ctime(&end_time) << " " << nf << " " << nph1 << " " << nph << endl;
+	}
+	filter->addToInterpolatedImage(buff,clusters,nph1);
+      }
+      
+      nph+=nph1;
+      if (nf%1000==0) {
+	std::time(&end_time);
+	cout << std::ctime(&end_time) << " " << nf << " " << nph1 << " " << nph << endl;
+      }
+      nf++;
+      delete [] buff;
+      iFrame=-1;
+    } 
+ 
+    if (filebin.is_open())
+      filebin.close();	  
+    else
+      cout << "could not open file " << fname << endl;
+    
+    
+    himage=interp->getInterpolatedImage();
+    for (int ii=0; ii<nsubpix*nsubpix*NC*NR; ii++) {
+      image[ii]=(float)himage[ii];
+    }
+    WriteToTiff(image, "/scratch/int_image.tiff", nsubpix*NR, nsubpix*NC);
+      delete interp;
+    delete decoder;
+    cout << "Read " << nf << " frames" << endl;
+    return NULL;
+}
+  
+int main(int argc, char *argv[]){
+
+   moenchProcessFrame();
+
+}

slsDetectorCalibration/singlePhotonDetector.h

@@ -1,50 +1,35 @@
 #ifndef SINGLEPHOTONDETECTOR_H
-#define  SINGLEPHOTONDETECTOR_H
+#define SINGLEPHOTONDETECTOR_H
 
-
+#include "analogDetector.h" 
-#include "slsDetectorData.h"
 
 #include "single_photon_hit.h"
-#include "pedestalSubtraction.h"
+
-#include "commonModeSubtraction.h"
 
 
 //#define MYROOT1
 
 #ifdef MYROOT1
 #include <TTree.h>
-
 #endif
 
 
-#include <iostream>
+#ifndef EVTYPE_DEF
-
+#define EVTYPE_DEF
-using namespace std;
+enum eventType {
-
+  PEDESTAL=0,
-
+  NEIGHBOUR=1,
-  enum eventType {
+  PHOTON=2,
-    PEDESTAL=0,
+  PHOTON_MAX=3,
-    NEIGHBOUR=1,
+  NEGATIVE_PEDESTAL=4,
-    PHOTON=2,
+  UNDEFINED_EVENT=-1
-    PHOTON_MAX=3,
+};
-    NEGATIVE_PEDESTAL=4,
-    UNDEFINED_EVENT=-1
-  };
-
-#ifndef DEF_QUAD
-#define DEF_QUAD
-  enum quadrant {
-    TOP_LEFT=0,
-    TOP_RIGHT=1,
-    BOTTOM_LEFT=2,
-    BOTTOM_RIGHT=3,
-    UNDEFINED_QUADRANT=-1
-  };
 #endif
 
-
+//template <class dataType> class singlePhotonDetector : 
-template <class dataType>
+//public analogDetector<dataType> {
-class singlePhotonDetector {
+class singlePhotonDetector : 
+public analogDetector<uint16_t> {
 
   /** @short class to perform pedestal subtraction etc. and find single photon clusters for an analog detector */
 
@@ -66,27 +51,25 @@ class singlePhotonDetector {
   */
   
 
-  singlePhotonDetector(slsDetectorData<dataType> *d, 
+ singlePhotonDetector(slsDetectorData<uint16_t> *d,
-		       int csize=3, 
+		      int csize=3,
-		       double nsigma=5,  
+		      double nsigma=5,
-		       int sign=1, 
+		      int sign=1,
-		       commonModeSubtraction *cm=NULL,
+		      commonModeSubtraction *cm=NULL,
-		       int nped=1000, 
+		      int nped=1000,
-		       int nd=100) : det(d), nx(0), ny(0), stat(NULL), cmSub(cm),  nDark(nd), eventMask(NULL),nSigma (nsigma), clusterSize(csize), clusterSizeY(csize), cluster(NULL), iframe(-1), dataSign(sign), quad(UNDEFINED_QUADRANT), tot(0), quadTot(0) {
+		      int nd=100, int nnx=-1, int nny=-1) : analogDetector<uint16_t>(d, sign, cm, nnx, nny),   nDark(nd), eventMask(NULL),nSigma (nsigma), clusterSize(csize), clusterSizeY(csize), cluster(NULL),   quad(UNDEFINED_QUADRANT), tot(0), quadTot(0) {
-
+    
-   
+    
-    det->getDetectorSize(nx,ny);
+    
     
-
-
-    stat=new pedestalSubtraction*[ny];
     eventMask=new eventType*[ny];
     for (int i=0; i<ny; i++) {
-      stat[i]=new pedestalSubtraction[nx];
-      stat[i]->SetNPedestals(nped);
       eventMask[i]=new eventType[nx];
+      for (int ix=0; ix<nx; ix++) {
+	stat[i][ix].SetNPedestals(nped);
+      }
     }
-   
+    
     if (ny==1)
       clusterSizeY=1;
 
@@ -97,71 +80,14 @@ class singlePhotonDetector {
     /**
        destructor. Deletes the cluster structure and the pdestalSubtraction array
     */
-    virtual ~singlePhotonDetector() {delete cluster; for (int i=0; i<ny; i++) delete [] stat[i]; delete [] stat;};
+  virtual ~singlePhotonDetector() {delete cluster;};
 
     
-    /** resets the pedestalSubtraction array and the commonModeSubtraction */
+  
-    void newDataSet(){iframe=-1; for (int iy=0; iy<ny; iy++) for (int ix=0; ix<nx; ix++) stat[iy][ix].Clear(); if (cmSub) cmSub->Clear(); };  
+    /* /\** resets the eventMask to undefined and the commonModeSubtraction *\/ */
-
+    /* void newFrame(){analogDetector::newFrame(); for (int iy=0; iy<ny; iy++) for (int ix=0; ix<nx; ix++) eventMask[iy][ix]=UNDEFINED_EVENT; }; */
-    /** resets the eventMask to undefined and the commonModeSubtraction */
-    void newFrame(){iframe++; for (int iy=0; iy<ny; iy++) for (int ix=0; ix<nx; ix++) eventMask[iy][ix]=UNDEFINED_EVENT; if (cmSub) cmSub->newFrame();};
-
-
-    /** sets the commonModeSubtraction algorithm to be used 
-	\param cm commonModeSubtraction algorithm to be used (NULL unsets) 
-	\returns pointer to the actual common mode subtraction algorithm
-    */
-    commonModeSubtraction *setCommonModeSubtraction(commonModeSubtraction *cm) {cmSub=cm; return cmSub;};
-
-
-    /**
-       sets the sign of the data
-       \param sign 1 means positive values for photons, -1 negative, 0 gets
-       \returns current sign for the data
-    */
-    int setDataSign(int sign=0) {if (sign==1 || sign==-1) dataSign=sign; return dataSign;};
 
     
-    /**
-       adds value to pedestal (and common mode) for the given pixel
-       \param val value to be added
-       \param ix pixel x coordinate
-       \param iy pixel y coordinate
-    */
-    virtual void addToPedestal(double val, int ix, int iy){ 
-      //	cout << "*"<< ix << " " << iy << " " << val << endl;
-      if (ix>=0 && ix<nx && iy>=0 && iy<ny) {
-	//	cout << ix << " " << iy << " " << val << endl;
-	stat[iy][ix].addToPedestal(val); 
-	if (cmSub && det->isGood(ix, iy) ) 
-	  cmSub->addToCommonMode(val, ix, iy);
-      };
-    };
-
-  /**
-       gets  pedestal (and common mode)
-       \param ix pixel x coordinate
-       \param iy pixel y coordinate
-       \param cm 0 (default) without common mode subtraction, 1 with common mode subtraction (if defined)
-    */
-    virtual double getPedestal(int ix, int iy, int cm=0){if (ix>=0 && ix<nx && iy>=0 && iy<ny) if (cmSub && cm>0) return stat[iy][ix].getPedestal()-cmSub->getCommonMode(); else return stat[iy][ix].getPedestal(); else return -1;};
-
-  /**
-       gets  pedestal rms (i.e. noise)
-       \param ix pixel x coordinate
-       \param iy pixel y coordinate
-    */
-    double getPedestalRMS(int ix, int iy){if (ix>=0 && ix<nx && iy>=0 && iy<ny) return stat[iy][ix].getPedestalRMS();else return -1;};
-
-
-  /**
-       sets  pedestal
-       \param ix pixel x coordinate
-       \param iy pixel y coordinate
-       \param val value to set
-    */
-    virtual void setPedestal(int ix, int iy, double val, double rms=0){if (ix>=0 && ix<nx && iy>=0 && iy<ny) stat[iy][ix].setPedestal(val,rms);};
-
  
 
   
@@ -180,9 +106,9 @@ class singlePhotonDetector {
       if (n>0 && n!=clusterSize) {
 	if (n%2==0)
 	  n+=1;
-	clusterSize=n; 
+	clusterSize=n;
 	if (cluster)
-	  delete cluster;    
+	  delete cluster;
 	if (ny>1)
 	  clusterSizeY=clusterSize;
 	cluster=new single_photon_hit(clusterSize,clusterSizeY);
@@ -192,13 +118,15 @@ class singlePhotonDetector {
 
 
 
-    int *getNPhotons(char *data, double thr=-1) {
+    virtual int *getNPhotons(char *data, double thr=-1) {
 
       double val;
       int *nph=new int[nx*ny];
       double rest[ny][nx];
       int cy=(clusterSizeY+1)/2;
       int cs=(clusterSize+1)/2;
+      
+      
 
       int ccs=clusterSize;
       int ccy=clusterSizeY;
@@ -219,20 +147,16 @@ class singlePhotonDetector {
 
       for (int ix=0; ix<nx; ix++) {
       for (int iy=0; iy<ny; iy++) {
-	
+	 
-	val=dataSign*(det->getValue(data, ix, iy)-getPedestal(ix,iy,0));
-	
 	if (thr<=0) tthr=nSigma*getPedestalRMS(ix,iy);
+    
+	val=subtractPedestal(data,ix,iy);
+
+	nph[ix+nx*iy]=analogDetector<uint16_t>::getNPhotons(data,ix,iy,tthr);
 	
-	nn=val/tthr;
+	rest[iy][ix]=subtractPedestal(data,ix,iy)-nph[ix+nx*iy]*tthr;
-	
+	  
-	if (nn>0) {
+	  
-	  rest[iy][ix]=val-nn*tthr;
-	  nph[ix+nx*iy]=nn;
-	} else {
-	  rest[iy][ix]=val;
-	  nph[ix+nx*iy]=0;
-	}
       }
       
       }
@@ -242,9 +166,9 @@ class singlePhotonDetector {
 	
 	  if (thr<=0) tthr=nSigma*getPedestalRMS(ix,iy);
 	  
-	  max=0; 
+	  max=0;
 	  tl=0;
-	  tr=0; 
+	  tr=0;
 	  bl=0;
 	  br=0;
 	  
@@ -269,16 +193,13 @@ class singlePhotonDetector {
 		if (ir==0 && ic==0) {
 		  if (v>tthr) {
 		    eventMask[iy][ix]=PHOTON;
-		  } 
+		  }
 		}
 	      }
 	    }
 	  }
-	  
-
 	  //if (cluster->get_data(0,0)>=max) {
 	    if (rest[iy][ix]>=max) {
-	    
 	    if (bl>=br && bl>=tl && bl>=tr) {
 	      quad=BOTTOM_LEFT;
 	      quadTot=bl;
@@ -291,17 +212,15 @@ class singlePhotonDetector {
 	    } else if   (tr>=bl && tr>=tl && tr>=br) {
 	      quad=TOP_RIGHT;
 	      quadTot=tr;
-	    } 
+	    }
 	    if (rest[iy][ix]>tthr || tot>sqrt(ccy*ccs)*tthr || quadTot>sqrt(cy*cs)*tthr) {
 	      nph[ix+nx*iy]++;
 	      rest[iy][ix]-=tthr;
 	    }
-	  }      
+	  }
 	}
 	
       }
-
-      
       return nph;
 
     }
@@ -315,7 +234,7 @@ class singlePhotonDetector {
 	/param data pointer to the data
 	/param ix pixel x coordinate
 	/param iy pixel y coordinate
-	/param cm enable(1)/disable(0) common mode subtraction (if defined). 
+	/param cm enable(1)/disable(0) common mode subtraction (if defined).
 	/returns event type for the given pixel
     */
     eventType getEventType(char *data, int ix, int iy, int cm=0) {
@@ -326,42 +245,39 @@ class singlePhotonDetector {
     
       int cy=(clusterSizeY+1)/2;
       int cs=(clusterSize+1)/2;
-
+      double val;
-
-
       tot=0;
       quadTot=0;
       quad=UNDEFINED_QUADRANT;
 
-
       if (iframe<nDark) {
-	if (cm==0) {
+	addToPedestal(data, ix,iy);
-	  //  cout << "=" << endl;
-	  addToPedestal(det->getValue(data, ix, iy),ix,iy);
-	  cluster->set_data(dataSign*(det->getValue(data, ix, iy)-getPedestal(ix,iy,cm)), 0,0       );
-	  // cout << "=" << endl;
-	}
 	return UNDEFINED_EVENT;
       }
-      
+
+
+     
       
 
       //   if (eventMask[iy][ix]==UNDEFINED) {
 	
-	eventMask[iy][ix]=PEDESTAL;
+      eventMask[iy][ix]=PEDESTAL;
 	
 	
-	cluster->x=ix;
+      cluster->x=ix;
-	cluster->y=iy;
+      cluster->y=iy;
-	cluster->rms=getPedestalRMS(ix,iy);
+      cluster->rms=getPedestalRMS(ix,iy);
-	cluster->ped=getPedestal(ix,iy, cm);
+      cluster->ped=getPedestal(ix,iy, cm);
 	
 
-	for (int ir=-(clusterSizeY/2); ir<(clusterSizeY/2)+1; ir++) {
+      for (int ir=-(clusterSizeY/2); ir<(clusterSizeY/2)+1; ir++) {
-	  for (int ic=-(clusterSize/2); ic<(clusterSize/2)+1; ic++) {
+	for (int ic=-(clusterSize/2); ic<(clusterSize/2)+1; ic++) {
 	    if ((iy+ir)>=0 && (iy+ir)<ny && (ix+ic)>=0 && (ix+ic)<nx) {
-	      cluster->set_data(dataSign*(det->getValue(data, ix+ic, iy+ir)-getPedestal(ix+ic,iy+ir,cm)), ic, ir       );
+	      
-	      v=cluster->get_data(ic,ir);
+	      v=subtractPedestal(data, ix+ic, iy+ir);
+	      
+	      cluster->set_data(v, ic, ir);
+	      //  v=cluster->get_data(ic,ir);
 	      tot+=v;
 	      if (ir<=0 && ic<=0)
 		bl+=v;
@@ -395,7 +311,7 @@ class singlePhotonDetector {
      	} else if   (tr>=bl && tr>=tl && tr>=br) {
 	  quad=TOP_RIGHT;
 	  quadTot=tr;
-	} 
+	}
 	
 	if (max>nSigma*cluster->rms || tot>sqrt(clusterSizeY*clusterSize)*nSigma*cluster->rms || quadTot>cy*cs*nSigma*cluster->rms) {
 	  if (cluster->get_data(0,0)>=max) {
@@ -404,18 +320,136 @@ class singlePhotonDetector {
 	    eventMask[iy][ix]=PHOTON;
 	  }
 	} else if (eventMask[iy][ix]==PEDESTAL) {
-	  if (cm==0)
+	  if (cm==0) {
-	    addToPedestal(det->getValue(data, ix, iy),ix,iy);
+	    if (det)
+	      val=dataSign*det->getValue(data, ix, iy);
+	    else
+	      val=((double**)data)[iy][ix];
+	    addToPedestal(val,ix,iy);
+	  }
 	}
       
-	
+
 
       return  eventMask[iy][ix];
 
   };
 
+
+int getClusters(char *data, single_photon_hit *clusters) {
+
+ 
+  int nph=0;
+  double val[ny][nx];
+  int cy=(clusterSizeY+1)/2;
+  int cs=(clusterSize+1)/2;
+  int ir, ic;
+  
+  double max=0, tl=0, tr=0, bl=0,br=0, *v, vv;
+
+  if (iframe<nDark) {
+    addToPedestal(data);
+    return 0;
+  }
+  newFrame();
+  for (int ix=0; ix<nx; ix++) {
+    for (int iy=0; iy<ny; iy++) {
+      
+      max=0;
+      tl=0;
+      tr=0;
+      bl=0;
+      br=0;
+      tot=0;
+      quadTot=0;
+      quad=UNDEFINED_QUADRANT;
+
+     
+
+      eventMask[iy][ix]=PEDESTAL;
+      
+	
+      (clusters+nph)->rms=getPedestalRMS(ix,iy);
+      
+      
+      
+      for (int ir=-(clusterSizeY/2); ir<(clusterSizeY/2)+1; ir++) {
+	for (int ic=-(clusterSize/2); ic<(clusterSize/2)+1; ic++) {
+	  
+	  if ((iy+ir)>=iy && (iy+ir)<ny && (ix+ic)>=ix && (ix+ic)<nx) {
+	    val[iy+ir][ix+ic]=subtractPedestal(data,ix+ic,iy+ir);
+	  }
+	  
+	  v=&(val[iy+ir][ix+ic]);
+	  tot+=*v;
+	  if (ir<=0 && ic<=0)
+	    bl+=*v;
+	  if (ir<=0 && ic>=0)
+	    br+=*v;
+	  if (ir>=0 && ic<=0)
+	    tl+=*v;
+	  if (ir>=0 && ic>=0)
+	    tr+=*v;
+	  if (*v>max) {
+	    max=*v;
+	  }
+	  
+	  
+	  if (ir==0 && ic==0) {
+	    if (*v<-nSigma*cluster->rms)
+	      eventMask[iy][ix]=NEGATIVE_PEDESTAL;
+	  }
+	  
+	}
+      }
+      
+      if (bl>=br && bl>=tl && bl>=tr) {
+	(clusters+nph)->quad=BOTTOM_LEFT;
+	(clusters+nph)->quadTot=bl;
+      } else if (br>=bl && br>=tl && br>=tr) {
+	(clusters+nph)->quad=BOTTOM_RIGHT;
+	(clusters+nph)->quadTot=br;
+      } else if (tl>=br && tl>=bl && tl>=tr) {
+	(clusters+nph)->quad=TOP_LEFT;
+	(clusters+nph)->quadTot=tl;
+      } else if   (tr>=bl && tr>=tl && tr>=br) {
+	(clusters+nph)->quad=TOP_RIGHT;
+	(clusters+nph)->quadTot=tr;
+      }
+      
+      if (max>nSigma*cluster->rms || tot>sqrt(clusterSizeY*clusterSize)*nSigma*cluster->rms || ((clusters+nph)->quadTot)>sqrt(cy*cs)*nSigma*cluster->rms) {
+	if (val[iy][ix]>=max) {
+	  eventMask[iy][ix]=PHOTON_MAX;
+	  (clusters+nph)->tot=tot;
+	  (clusters+nph)->x=ix;
+	  (clusters+nph)->y=iy;
+	  (clusters+nph)->ped=getPedestal(ix,iy,0);
+	  for (int ir=-(clusterSizeY/2); ir<(clusterSizeY/2)+1; ir++) {
+	    for (int ic=-(clusterSize/2); ic<(clusterSize/2)+1; ic++) {
+	      (clusters+nph)->set_data(val[iy+ir][ix+ic],ic,ir);
+	    }
+	  }
+	  nph++;
+
+
+	  } else {
+	    eventMask[iy][ix]=PHOTON;
+	  }
+      } else if (eventMask[iy][ix]==PEDESTAL) {
+	addToPedestal(data,ix,iy);
+      }
+
+
+    }
+  }
+
+  return  nph;
+  
+};
+
+
     /**<
-       retrurns the total signal in a cluster
+       returns the total signal in a cluster
        \param size cluser size  should be 1,2 or 3
        \returns cluster center if size=1, sum of the maximum quadrant if size=2, total of the cluster if size=3 or anything else
     */
@@ -437,12 +471,7 @@ class singlePhotonDetector {
 
     quadrant getQuadrant() {return quad;};
 
-    /** sets/gets number of samples for moving average pedestal calculation
+  
-	\param i number of samples to be set (0 or negative gets)
-	\returns actual number of samples
-    */
-    int SetNPedestals(int i=-1) {int ix=0, iy=0; if (i>0) for (ix=0; ix<nx; ix++) for (iy=0; iy<ny; iy++) stat[iy][ix].SetNPedestals(i); return stat[0][0].SetNPedestals();};
-
     /** returns value for cluster element in relative coordinates
 	\param ic x coordinate (center is (0,0))
 	\param ir y coordinate (center is (0,0))
@@ -458,7 +487,7 @@ class singlePhotonDetector {
     eventType getEventMask(int ic, int ir=0){return eventMask[ir][ic];};
  
 
-#ifdef MYROOT1  
+#ifdef MYROOT1
     /** generates a tree and maps the branches
 	\param tname name for the tree
 	\param iFrame pointer to the frame number
@@ -479,6 +508,9 @@ class singlePhotonDetector {
       tall->Branch("data",cluster->data,tit);
       tall->Branch("pedestal",&(cluster->ped),"pedestal/D");
       tall->Branch("rms",&(cluster->rms),"rms/D");
+      tall->Branch("tot",&(cluster->tot),"tot/D");
+      tall->Branch("quadTot",&(cluster->quadTot),"quadTot/D");
+      tall->Branch("quad",&(cluster->quad),"quad/I");
       return tall;
     };
 #else
@@ -488,23 +520,15 @@ class singlePhotonDetector {
 #endif
 
 
- private:
+ protected:
   
-    slsDetectorData<dataType> *det; /**< slsDetectorData to be used */
-    int nx; /**< Size of the detector in x direction */
-    int ny; /**< Size of the detector in y direction */
 
-
-    pedestalSubtraction **stat; /**< pedestalSubtraction class */
-    commonModeSubtraction *cmSub;/**< commonModeSubtraction class */
     int nDark; /**< number of frames to be used at the beginning of the dataset to calculate pedestal without applying photon discrimination */
     eventType **eventMask; /**< matrix of event type or each pixel */
     double nSigma; /**< number of sigma parameter for photon discrimination */
     int clusterSize; /**< cluster size in the x direction */
     int clusterSizeY; /**< cluster size in the y direction i.e. 1 for strips, clusterSize for pixels */
     single_photon_hit *cluster; /**< single photon hit data structure */
-    int iframe;  /**< frame number (not from file but incremented within the dataset every time newFrame is called */
-    int dataSign; /**< sign of the data i.e. 1 if photon is positive, -1 if negative */
     quadrant quad; /**< quadrant where the photon is located */
     double tot; /**< sum of the 3x3 cluster */
     double quadTot; /**< sum of the maximum 2x2cluster */

slsDetectorCalibration/single_photon_hit.h

@@ -1,10 +1,21 @@
 #ifndef SINGLE_PHOTON_HIT_H
-#define SINGLE_PHOTON_HIT_h
+#define SINGLE_PHOTON_HIT_H
 
 typedef  double double32_t;
 typedef  float float32_t;
 typedef  int int32_t;
 
+#ifndef DEF_QUAD
+#define DEF_QUAD
+  enum quadrant {
+    TOP_LEFT=0,
+    TOP_RIGHT=1,
+    BOTTOM_LEFT=2,
+    BOTTOM_RIGHT=3,
+    UNDEFINED_QUADRANT=-1
+  };
+#endif
+
 
 class single_photon_hit {
 
@@ -15,20 +26,20 @@ class single_photon_hit {
       \param nx cluster size in x direction
       \param ny cluster size in y direction (defaults to 1 for 1D detectors)
   */
-  single_photon_hit(int nx, int ny=1): dx(nx), dy(ny) {data=new double[dx*dy];};
+  single_photon_hit(int nx=3, int ny=3): dx(nx), dy(ny) {data=new double[dx*dy];};
 
   ~single_photon_hit(){delete [] data;};  /**< destructor, deletes the data array */
 
   /** binary write to file of all elements of the structure, except size of the cluster
       \param myFile file descriptor
   */
-  void write(FILE *myFile) {fwrite((void*)this, 1, 3*sizeof(int)+2*sizeof(double), myFile); fwrite((void*)data, 1, dx*dy*sizeof(double), myFile);};   
+  void write(FILE *myFile) {fwrite((void*)this, 1, 3*sizeof(int)+4*sizeof(double)+sizeof(quad), myFile); fwrite((void*)data, 1, dx*dy*sizeof(double), myFile);};   
   
   /** 
       binary read from file of all elements of the structure, except size of the cluster. The structure is then filled with those args
       \param myFile file descriptor
   */
-  void read(FILE *myFile) {fread((void*)this, 1,  3*sizeof(int)+2*sizeof(double), myFile); fread((void*)data, 1, dx*dy*sizeof(double), myFile);};
+  void read(FILE *myFile) {fread((void*)this, 1,  3*sizeof(int)+4*sizeof(double)+sizeof(quad), myFile); fread((void*)data, 1, dx*dy*sizeof(double), myFile);};
 
     /** 
 	assign the value to the element of the cluster matrix, with relative coordinates where the center of the cluster is (0,0)
@@ -52,9 +63,12 @@ class single_photon_hit {
   double	rms; 		/**< noise of central pixel l -- at some point it can be removed*/
   double 	ped; 		/**< pedestal of the central pixel -- at some point it can be removed*/
   int	iframe; 	        /**< frame number */
-  double *data;		        /**< pointer to data */
   const int dx;	                /**< size of data cluster in x */
-  const int dy;                 /**< size of data cluster in y */
+  const int dy;                 /**< size of data cluster in y */   
+  quadrant quad; /**< quadrant where the photon is located */
+  double tot; /**< sum of the 3x3 cluster */
+  double quadTot; /**< sum of the maximum 2x2cluster */
+  double *data;		        /**< pointer to data */
 };
 
 

slsDetectorCalibration/slsDetectorData.h

@@ -18,60 +18,64 @@ class slsDetectorData {
   int **dataMap; /**< Array of size nx*ny storing the pointers to the data in the dataset (as offset)*/
   dataType **dataMask; /**< Array of size nx*ny storing the polarity of the data in the dataset (should be 0 if no inversion is required, 0xffffffff is inversion is required) */
   int **dataROIMask; /**< Array of size nx*ny 1 if channel is good (or in the ROI), 0 if bad channel (or out of ROI)   */
- int *xmap;
+  int *xmap;
- int *ymap;
+  int *ymap;
-
+  
  public:
-
+  
   /**
+     
+     
+     General slsDetectors data structure. Works for data acquired using the slsDetectorReceiver. Can be generalized to other detectors (many virtual funcs).
 
-
+     Constructor (no error checking if datasize and offsets are compatible!)
-  General slsDetectors data structure. Works for data acquired using the slsDetectorReceiver. Can be generalized to other detectors (many virtual funcs).
-
-  Constructor (no error checking if datasize and offsets are compatible!)
   \param npx number of pixels in the x direction
   \param npy number of pixels in the y direction (1 for strips)
   \param dsize size of the data
   \param dMap array of size nx*ny storing the pointers to the data in the dataset (as offset)
   \param dMask Array of size nx*ny storing the polarity of the data in the dataset (should be 0 if no inversion is required, 0xffffffff is inversion is required)
   \param dROI Array of size nx*ny. The elements are 1s if the channel is good or in the ROI, 0 is bad or out of the ROI. NULL (default) means all 1s. 
-
+  
   */
-  slsDetectorData(int npx, int npy, int dsize, int **dMap=NULL, dataType **dMask=NULL, int **dROI=NULL): nx(npx), ny(npy), dataSize(dsize) {
+ slsDetectorData(int npx, int npy, int dsize, int **dMap=NULL, dataType **dMask=NULL, int **dROI=NULL): nx(npx), ny(npy), dataSize(dsize) {
-
+    
-   
+    
     xmap=new int[dsize/sizeof(dataType)];
     ymap=new int[dsize/sizeof(dataType)];
-
+    
 
     //  if (dataMask==NULL) {
-      dataMask=new dataType*[ny];
+    dataMask=new dataType*[ny];
-      for(int i = 0; i < ny; i++) {
+    for(int i = 0; i < ny; i++) {
-	dataMask[i] = new dataType[nx];
+      dataMask[i] = new dataType[nx];
-      }
+    }
-      // }
+    // }
     
     // if (dataMap==NULL) {
-      dataMap=new int*[ny];
+    dataMap=new int*[ny];
-      for(int i = 0; i < ny; i++) {
+    for(int i = 0; i < ny; i++) {
-	dataMap[i] = new int[nx];
+      dataMap[i] = new int[nx];
       }
-      // }
+    // }
-      // if (dataROIMask==NULL) {
+    // if (dataROIMask==NULL) {
-      dataROIMask=new int*[ny];
+    dataROIMask=new int*[ny];
-      for(int i = 0; i < ny; i++) {
+    for(int i = 0; i < ny; i++) {
-	dataROIMask[i] = new int[nx];
+      dataROIMask[i] = new int[nx];
-	for (int j=0; j<nx; j++)
+      for (int j=0; j<nx; j++)
-	  dataROIMask[i][j]=1;
+	dataROIMask[i][j]=1;
+    }
+    // }
+    for (int ip=0; ip<dsize/sizeof(dataType); ip++){
+      xmap[ip]=-1;
+      ymap[ip]=-1;
       }
-      // }
+    
-
     setDataMap(dMap);
     setDataMask(dMask);
     setDataROIMask(dROI);
-
+    
   };
-
+  
   virtual ~slsDetectorData() {
     for(int i = 0; i < ny; i++) {
       delete [] dataMap[i];
@@ -83,38 +87,48 @@ class slsDetectorData {
     delete [] dataROIMask;
     delete [] xmap;
     delete [] ymap;
-  }
+  };
-
+  
-
+  
   /**
      defines the data map (as offset) - no error checking if datasize and offsets are compatible!
      \param dMap array of size nx*ny storing the pointers to the data in the dataset (as offset). If NULL (default),the data are arranged as if read out row by row (dataMap[iy][ix]=(iy*nx+ix)*sizeof(dataType);)
   */
   void setDataMap(int **dMap=NULL) {
-
+    
-   /* int ip;*/
+    int ip=0;
-
+    int ix, iy;
-	  if (dMap==NULL) {
+    if (dMap==NULL) {
-		  for (int iy=0; iy<ny; iy++)
+      for (iy=0; iy<ny; iy++) {
-			  for (int ix=0; ix<nx; ix++)
+	for (ix=0; ix<nx; ix++) {
-				  dataMap[iy][ix]=(iy*nx+ix)*sizeof(dataType);
+		dataMap[iy][ix]=(iy*nx+ix)*sizeof(dataType);
-	  } else {
+	}
-		  //cout << "set dmap "<< dataMap << " " << dMap << endl;
+      }
-		  for (int iy=0; iy<ny; iy++){
+    } else {
+      //cout << "set dmap "<< dataMap << " " << dMap << endl;
+      for (iy=0; iy<ny; iy++){
 			 // cout << iy << endl;
-			  for (int ix=0; ix<nx; ix++) {
+	for (ix=0; ix<nx; ix++) {
-				  dataMap[iy][ix]=dMap[iy][ix];
+	  dataMap[iy][ix]=dMap[iy][ix];
-				 // cout << ix << " " << iy << endl;
+	  // cout << ix << " " << iy << endl;
-				 /*ip=dataMap[ix][iy]/sizeof(dataType);
+	  /*ip=dataMap[ix][iy]/sizeof(dataType);
-				  xmap[ip]=ix;
+	    xmap[ip]=ix;
-				  ymap[ip]=iy;Annaa*/
+	    ymap[ip]=iy;Annaa*/
-			  }
+	}
-		  }
+      }
-	  }
+    }
-	 // cout << "nx:" <<nx << " ny:" << ny << endl;
+    for (iy=0; iy<ny; iy++){
-
+      for (ix=0; ix<nx; ix++) {
+	ip=dataMap[ix][iy]/sizeof(dataType);
+	xmap[ip]=ix;
+	ymap[ip]=iy;
+      }
+    }
+    
+	  // cout << "nx:" <<nx << " ny:" << ny << endl;
+    
   };
-
+  
 
   /**
      defines the data mask i.e. the polarity of the data
@@ -122,41 +136,41 @@ class slsDetectorData {
 
   */
   void setDataMask(dataType **dMask=NULL){
-
+    
     if (dMask!=NULL) {
-
+      
       for (int iy=0; iy<ny; iy++)
 	 for (int ix=0; ix<nx; ix++)
 	   dataMask[iy][ix]=dMask[iy][ix];
     } else {
       for (int iy=0; iy<ny; iy++)
-	 for (int ix=0; ix<nx; ix++)
+	for (int ix=0; ix<nx; ix++)
-	   dataMask[iy][ix]=0;
+	  dataMask[iy][ix]=0;
     }
   };
   /**
      defines the region of interest and/or the bad channels mask
      \param dROI Array of size nx*ny. The lements are 1s if the channel is good or in the ROI, 0 is bad or out of the ROI. NULL (default) means all 1s. 
-
+     
   */
   void setDataROIMask(int **dROI=NULL){
-
+    
     if (dROI!=NULL) {
-
+      
       for (int iy=0; iy<ny; iy++)
-	 for (int ix=0; ix<nx; ix++)
+	for (int ix=0; ix<nx; ix++)
 	   dataROIMask[iy][ix]=dROI[iy][ix];
     } else {
-
-      for (int iy=0; iy<ny; iy++)
-	 for (int ix=0; ix<nx; ix++)
-	   dataROIMask[iy][ix]=1;
-
-    }
       
+      for (int iy=0; iy<ny; iy++)
+	for (int ix=0; ix<nx; ix++)
+	   dataROIMask[iy][ix]=1;
+      
+    }
+    
     
   };
-
+  
   /**
      Define bad channel or roi mask for a single channel
      \param ix channel x coordinate
@@ -172,7 +186,7 @@ class slsDetectorData {
      \returns 1 if pixel is good, 0 if it's bad, -1 if pixel is out of range
   */
   int isGood(int ix, int iy) {  if (ix>=0 && ix<nx && iy>=0 && iy<ny) return dataROIMask[iy][ix]; else return -1;};
-
+  
   /**
      Returns detector size in x,y
      \param npx reference to number of channels in x
@@ -180,18 +194,18 @@ class slsDetectorData {
      \returns total number of channels
   */
   int getDetectorSize(int &npx, int &npy){npx=nx; npy=ny; return nx*ny;};
-
+  
   /** Returns the size of the data frame */
   int getDataSize() {return dataSize;};
-
+  
   /** changes the size of the data frame */
   int setDataSize(int d) {dataSize=d; return dataSize;};
 
-
+  
   virtual void getPixel(int ip, int &x, int &y) {x=xmap[ip]; y=ymap[ip];};
-
+  
   virtual dataType **getData(char *ptr, int dsize=-1) {
-   
+    
     dataType **data;
     int ix,iy;
     data=new dataType*[ny];
@@ -206,9 +220,28 @@ class slsDetectorData {
       }
     }
     return data;
-
+    
   };
-
+  
+  virtual double **getImage(char *ptr, int dsize=-1) {
+    
+    double **data;
+    int ix,iy;
+    data=new double*[ny];
+    for(int i = 0; i < ny; i++) {
+      data[i]=new double[nx];
+    }
+    if (dsize<=0 || dsize>dataSize) dsize=dataSize;
+    for (int ip=0; ip<(dsize/sizeof(dataType)); ip++) {
+      getPixel(ip,ix,iy);
+      if (ix>=0 && ix<nx && iy>=0 && iy<ny) {
+	data[iy][ix]=getValue(ptr,ix,iy);
+      }
+    }
+    return data;
+    
+  };
+  
   /** 
      Returns the value of the selected channel for the given dataset. Virtual function, can be overloaded.
      \param data pointer to the dataset (including headers etc)
@@ -237,11 +270,11 @@ class slsDetectorData {
      \returns data for the selected channel, with inversion if required or -1 if its a missing packet
 
   */
-
+  
   virtual int getChannelwithMissingPackets(char *data, int ix, int iy) {
-	  return 0;
+    return 0;
   };
-
+  
 
 
   /**
@@ -252,48 +285,56 @@ class slsDetectorData {
      \returns data for the selected channel, with inversion if required as double
 
   */
-  virtual double getValue(char *data, int ix, int iy=0) {/* cout << " x "<< ix << " y"<< iy << " val " << getChannel(data, ix, iy)<< endl;*/return (double)getChannel(data, ix, iy);};
+  virtual double getValue(char *data, int ix, int iy=0) {
+    /* cout << " x "<< ix << " y"<< iy << " val " << getChannel(data, ix, iy)<< endl;*/
+    return (double)getChannel(data, ix, iy);
+  };
 
-
+  
-   /**
+  /**
      Returns the frame number for the given dataset. Purely virtual func.
      \param buff pointer to the dataset
      \returns frame number
-
+     
   */
   virtual  int getFrameNumber(char *buff)=0;   
-
+  
   /**
-
+     
      Returns the packet number for the given dataset. purely virtual func
      \param buff pointer to the dataset
      \returns packet number number
 
-
+     
-  virtual int getPacketNumber(char *buff)=0;
+     virtual int getPacketNumber(char *buff)=0;
-
+     
   */
-
+  
-   /**
+  /**
 
      Loops over a memory slot until a complete frame is found (i.e. all packets 0 to nPackets, same frame number). purely virtual func
      \param data pointer to the memory to be analyzed
      \param ndata reference to the amount of data found for the frame, in case the frame is incomplete at the end of the memory slot
      \param dsize size of the memory slot to be analyzed
      \returns pointer to the beginning of the last good frame (might be incomplete if ndata smaller than dataSize), or NULL if no frame is found 
-
+     
   */
   virtual  char *findNextFrame(char *data, int &ndata, int dsize)=0;
-
+  
-
+  
    /**
-
+      
      Loops over a file stream until a complete frame is found (i.e. all packets 0 to nPackets, same frame number). Can be overloaded for different kind of detectors! 
      \param filebin input file stream (binary)
      \returns pointer to the begin of the last good frame, NULL if no frame is found or last frame is incomplete
 
-  */
+   */
   virtual char *readNextFrame(ifstream &filebin)=0;
+  
+  
+
+
+
 
 };
 

slsDetectorCalibration/tiffIO.h

@@ -0,0 +1,92 @@
+#ifndef TIFF_IO_H
+#define  TIFF_IO_H
+
+#include <vector>
+#include <string>
+#include <sstream>
+#include <iomanip>
+#include <fstream>
+#include <stdio.h>
+#include <fstream>
+
+
+/*****************************************************************************/
+//
+//CBFlib must be installed to use this program
+//
+/*****************************************************************************/
+#include "tiffio.h"
+
+#undef cbf_failnez 
+#define cbf_failnez(x)				\
+  {						\
+  int err;	        			\
+  err = (x);					\
+  if (err) {						\
+  fprintf(stderr,"\nCBFlib fatal error %x \n",err);	\
+  exit(-1);						\
+  }							\
+  }
+
+void *WriteToTiff(float * imgData, const char * imgname, int nrow, int ncol){
+  int sampleperpixel=1;
+  // unsigned char * buff=NULL;
+  tsize_t linebytes;
+
+  TIFF * tif = TIFFOpen(imgname,"w");
+  if (tif) {
+    TIFFSetField(tif,TIFFTAG_IMAGEWIDTH,ncol);
+    TIFFSetField(tif, TIFFTAG_IMAGELENGTH, nrow);
+    TIFFSetField(tif, TIFFTAG_SAMPLESPERPIXEL,sampleperpixel);
+    TIFFSetField(tif, TIFFTAG_BITSPERSAMPLE, sizeof(float)*8);
+    TIFFSetField(tif, TIFFTAG_ORIENTATION, ORIENTATION_BOTLEFT); 
+    TIFFSetField(tif, TIFFTAG_PLANARCONFIG, PLANARCONFIG_CONTIG);
+    TIFFSetField(tif, TIFFTAG_PHOTOMETRIC, PHOTOMETRIC_MINISBLACK);
+    TIFFSetField(tif, TIFFTAG_SAMPLEFORMAT, SAMPLEFORMAT_IEEEFP);
+    
+    linebytes = sampleperpixel*ncol;
+    TIFFSetField(tif, TIFFTAG_ROWSPERSTRIP, TIFFDefaultStripSize(tif, ncol*sampleperpixel));
+    for(int irow=0; irow<nrow; irow++){
+      TIFFWriteScanline(tif,&imgData[irow*ncol],irow,0);
+    }
+    
+    TIFFClose(tif);
+  } else
+    cout << "could not open file " << imgname << " for writing " << endl;
+
+  return NULL;
+};
+
+float *ReadFromTiff( const char * imgname, uint32 &nrow, uint32 &ncol){
+  // unsigned char * buff=NULL;
+
+  TIFF * tif = TIFFOpen(imgname,"r");
+  if (tif){
+        uint32 bps;
+	uint32 sampleperpixel=1;
+  tsize_t linebytes;
+
+ uint32 imagelength;
+
+    TIFFGetField(tif,TIFFTAG_IMAGEWIDTH,&ncol);
+    TIFFGetField(tif, TIFFTAG_IMAGELENGTH, &nrow);
+    TIFFSetField(tif, TIFFTAG_SAMPLESPERPIXEL,sampleperpixel);
+    TIFFSetField(tif, TIFFTAG_BITSPERSAMPLE, &bps);
+    TIFFGetField(tif, TIFFTAG_IMAGELENGTH, &imagelength);
+
+    float * imgData=new float[ncol*nrow];
+    linebytes = sampleperpixel*ncol;
+    //  TIFFSetField(tif, TIFFTAG_ROWSPERSTRIP, TIFFDefaultStripSize(tif, ncol*sampleperpixel));
+    for(int irow=0; irow<nrow; irow++){
+      //tiffreadscanline(tif, buf, row);
+      TIFFReadScanline(tif,&imgData[irow*ncol],irow);
+    }
+
+    TIFFClose(tif);
+    return imgData;
+  }  else
+    cout << "could not open file " << imgname << " for reading " << endl;
+  return NULL;
+};
+
+#endif