Implemenging common mode and ghost corrections for moench03 (and generalized)

slsDetectorCalibration/analogDetector.h

@@ -6,7 +6,8 @@
 #include <pthread.h>
 #include "slsDetectorData.h"
 #include "pedestalSubtraction.h"
-#include "commonModeSubtraction.h"
+#include "commonModeSubtractionNew.h"
+#include "ghostSummation.h"
 #include "tiffIO.h"
 #include "slsInterpolation.h"
 
@@ -61,7 +62,7 @@ template <class dataType> class analogDetector {
   
 
  analogDetector(slsDetectorData<dataType> *d, int sign=1, 
-		commonModeSubtraction *cm=NULL, int nped=1000, int nnx=-1, int nny=-1, double *gm=NULL) : det(d), nx(nnx), ny(nny), stat(NULL), cmSub(cm),  iframe(-1), dataSign(sign), gmap(gm), id(0) {
+		commonModeSubtraction *cm=NULL, int nped=1000, int nnx=-1, int nny=-1, double *gm=NULL, ghostSummation<dataType> *gs=NULL) : det(d), nx(nnx), ny(nny), stat(NULL), cmSub(cm),  iframe(-1), dataSign(sign), gmap(gm), ghSum(gs), id(0) {
     
     if (det)
       det->getDetectorSize(nx,ny);
@@ -118,7 +119,7 @@ template <class dataType> class analogDetector {
     dataSign=orig->dataSign;
     iframe=orig->iframe;
     gmap=orig->gmap;
-    cmSub=orig->cmSub;
+    //  cmSub=orig->cmSub;
     id=orig->id;
     xmin=orig->xmin;
     xmax=orig->xmax;
@@ -153,7 +154,9 @@ template <class dataType> class analogDetector {
      hs9=(TH2F*)(orig->hs9)->Clone();//new TH2F("hs","hs",(orig->hs)-getNbins,-500,9500,nx*ny,-0.5,nx*ny-0.5); 
 #endif
 #endif
-
+     if (orig->cmSub) cmSub=(orig->cmSub)->Clone();
+     if (orig->ghSum) ghSum=(orig->ghSum)->Clone();
+     else ghSum=NULL;
   }
 
 
@@ -276,6 +279,9 @@ template <class dataType> class analogDetector {
   /** resets the commonModeSubtraction and increases the frame index */
   virtual void newFrame(){iframe++; if (cmSub) cmSub->newFrame();};
 
+  /** resets the commonModeSubtraction and increases the frame index */
+  virtual void newFrame(char *data){iframe++; if (cmSub) cmSub->newFrame(); calcGhost(data);};
+  
 
     /** sets the commonModeSubtraction algorithm to be used 
 	\param cm commonModeSubtraction algorithm to be used (NULL unsets) 
@@ -289,6 +295,11 @@ template <class dataType> class analogDetector {
     commonModeSubtraction *getCommonModeSubtraction() {return cmSub;};
 
 
+    ghostSummation<dataType> *getGhostSummation(){return ghSum;};
+    ghostSummation<dataType> *setGhostSummation(ghostSummation<dataType> *gs){ghSum=gs; return ghSum;}; 
+
+
+
     /**
        sets the sign of the data
        \param sign 1 means positive values for photons, -1 negative, 0 gets
@@ -306,9 +317,15 @@ template <class dataType> class analogDetector {
     */
     virtual void addToPedestal(double val, int ix, int iy=0, int cm=0){ 
       if (ix>=0 && ix<nx && iy>=0 && iy<ny) { 
+	
+	//	cout << val << " " ;
 	if (cmSub && cm>0)  { 
 	  val-=	getCommonMode(ix, iy);
 	}
+	//	cout << val << " " ;
+	val+=getGhost(ix,iy);
+	//	cout << val ;
+	//	cout << endl;
 	stat[iy][ix].addToPedestal(val); 
 	/* if (cmSub && cm>0)  { */
 	/*   if (det) if (det->isGood(ix, iy)==0) return; */
@@ -324,23 +341,29 @@ template <class dataType> class analogDetector {
 
 
     virtual void addToCommonMode(char *data){ 
+      // cout << "+"<< endl;
       if (cmSub) {
+	//	cout << "*" << endl;
 	  for (int iy=ymin; iy<ymax; iy++) { 
 	for (int ix=xmin; ix<xmax; ix++) {
 	    // if (getNumpedestals(ix,iy)>0) 
-	    if (det->isGood(ix,iy))
+	  //  if (det->isGood(ix,iy)) {
 	      addToCommonMode(data, ix, iy);
+	      //  cout << ":";
+	      // }
 	}
 	  }
 	//cout << "cm " << getCommonMode(0,0) << " " << getCommonMode(1,0) << endl;
       }
     }      
+
     virtual void addToCommonMode(char *data, int ix, int iy=0){ 
       if (cmSub) {
 	if (det) if (det->isGood(ix, iy)==0) return;
 	if (getNumpedestals(ix,iy)>0){
 	    cmSub->addToCommonMode(subtractPedestal(data,ix,iy,0), ix, iy);
-	    // cout << ix << " " <<subtractPedestal(data,ix,iy,0)  << endl;
+	    //  cout << ":";
+	    // cout << ix << " " <<" " << iy << subtractPedestal(data,ix,iy,0)  << endl;
 	}
       }
     }      
@@ -351,7 +374,7 @@ template <class dataType> class analogDetector {
        \param cm 0 (default) without common mode subtraction, 1 with common mode subtraction (if defined)
        \returns pedestal value
     */
-    virtual double getPedestal(int ix, int iy, int cm=0){
+    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()+getCommonMode(ix,iy); 
@@ -482,7 +505,10 @@ template <class dataType> class analogDetector {
 
     }
     
+  virtual void calcGhost(char *data, int ix, int iy=1) {if (ghSum) ghSum->calcGhost(data, ix, iy);};
 
+  virtual void calcGhost(char *data){if (ghSum) ghSum->calcGhost(data);};
+  virtual double getGhost(int ix, int iy) {if (ghSum) return ghSum->getGhost(ix, iy); return 0;};
 
      /**
        write 32bit tiff file with detector image data
@@ -691,9 +717,12 @@ template <class dataType> class analogDetector {
     virtual void addToPedestal(char *data, int cm=0) {
       
       //    cout << "add to pedestal " << endl;
-      newFrame();
+      newFrame(data);
       
-      if (cmSub) {
+      //calcGhost(data);
+
+
+      if (cmSub && cm) {
 	addToCommonMode(data);
       }
             
@@ -702,7 +731,8 @@ template <class dataType> class analogDetector {
 	for (int iy=ymin; iy<ymax; iy++) {
       for (int ix=xmin; ix<xmax; ix++) {
 	  if (det->isGood(ix,iy)) {
-	    addToPedestal(data,ix,iy,1);
+	    // addToPedestal(data,ix,iy,1);
+	    addToPedestal(data,ix,iy,cm);
 	  //if (ix==10 && iy==10) 
 	  // cout <<ix << " " << iy << " " << getPedestal(ix,iy)<< endl;
 #ifdef ROOTSPECTRUM 
@@ -792,12 +822,14 @@ template <class dataType> class analogDetector {
 	  val=dataSign*det->getValue(data, ix, iy);
 	else
 	  val=((double*)data)[iy*nx+ix];
-	
+	//	cout << val << endl;
+	val+=getGhost(ix,iy);
 	 /* if (ix==10 && iy==10)  */
 	 /*     cout << ix << " " << iy << " " << val ;  */
 	  /* if (ix==100 && iy==100) */
 	  /*   cout << ix << " " << iy << " " << val; */
 	addToPedestal(val,ix,iy);
+	//	cout << val << endl;
 	  /* if (ix==10 && iy==10)  */
 	  /*   cout <<" " << getPedestal(ix,iy)<< endl; */
 	  /* if (ix==100 && iy==100) */
@@ -818,11 +850,13 @@ template <class dataType> class analogDetector {
       
     virtual  int *subtractPedestal(char *data, int *val=NULL, int cm=0) {
       
-      newFrame();
+      newFrame(data);
       
       if (val==NULL)
 	val=image;//new double[nx*ny];
       
+      //calcGhost(data);
+
       for (int iy=ymin; iy<ymax; iy++) {
 	for (int ix=xmin; ix<xmax; ix++) {
 	  if (det->isGood(ix,iy))
@@ -861,6 +895,8 @@ template <class dataType> class analogDetector {
 	}	else
 	  val= (((double*)data)[iy*nx+ix]-getPedestal(ix,iy))/g;
     
+	val+=getGhost(ix,iy)/g;
+
 #ifdef ROOTSPECTRUM
 	hs->Fill(val,(iy-ymin)*(xmax-xmin)+(ix-xmin));
 #ifdef ROOTCLUST
@@ -949,10 +985,10 @@ template <class dataType> class analogDetector {
        double val;
        if (nph==NULL)
 	 nph=image;
-       newFrame();
+       newFrame(data);
 
+       //calcGhost(data);
        addToCommonMode(data);
-
 	 for (int iy=ymin; iy<ymax; iy++) {
        for (int ix=xmin; ix<xmax; ix++) {
 	   if (det->isGood(ix,iy))
@@ -1032,6 +1068,9 @@ template <class dataType> class analogDetector {
       if (ymi<0) ymi=ymin;
       if (yma<0) yma=ymax;
       
+      newFrame(data);
+      //calcGhost(data);
+      addToCommonMode(data);
       for (int iy=ymi; iy<yma; iy++)
 	for (int ix=xmi; ix<xma; ix++)
 	  if (det->isGood(ix,iy)) {
@@ -1140,6 +1179,7 @@ FILE *getFilePointer(){return myFile;};
     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 */
     double *gmap;
+    ghostSummation<dataType> *ghSum;/**< ghostSummation class */
     int *image;
     int id;
     //int xmin, xmax, ymin, ymax;  int xmin; /**< minimum x of the region of interest */

slsDetectorCalibration/commonModeSubtractionNew.h

@@ -23,6 +23,13 @@ class commonModeSubtraction {
     /** destructor - deletes the moving average(s) and the sum of pedestals calculator(s) */
     virtual ~commonModeSubtraction() {delete [] mean; delete [] mean2; delete [] nCm;};
     
+    /* commonModeSubtraction(commonModeSubtraction *cs) { */
+    /*   if (cs) new commonModeSubtraction(cs->getNRoi(), cs->nsigma); */
+    /* } */
+    
+    virtual commonModeSubtraction *Clone() {
+      new commonModeSubtraction(this->nROI, this->nsigma);
+    }
 
     /** clears the moving average and the sum of pedestals calculation - virtual func*/
     virtual void  Clear(){
@@ -52,7 +59,9 @@ class commonModeSubtraction {
       // if (iroi==0) val=100;
       // else val=-100;
       // if (isc>=0 && isc<nROI) {
+      //	cout << ix << " " << iy << " " << iroi << endl;
       if (iroi>=0 && iroi<nROI) {
+	//	cout << ix << " " << iy << " " << iroi << endl;
 	mean[iroi]+=val; 
 	mean2[iroi]+=val*val; 
 	nCm[iroi]++;
@@ -70,7 +79,7 @@ class commonModeSubtraction {
       /* 	return 100; */
       /* else */
       /* 	return -100; */
-
+      // cout << "*" << ix << " " << iy << " " << iroi << " " << mean[iroi] << " " << nCm[iroi]<< endl;
       if (iroi>=0 && iroi<nROI) {
 	if (nCm[iroi]>0)  
 	  return mean[iroi]/nCm[iroi]; 
@@ -96,8 +105,7 @@ class commonModeSubtraction {
 	gets the common mode ROI for pixel ix, iy -should be overloaded!
      */
     virtual int getROI(int ix, int iy){ (void) ix; (void) iy; return 0;};
-
+    int getNRoi(){return nROI;};
-
 
  protected:
   double *mean; /**<array of moving average of the pedestal average per region of interest */

slsDetectorCalibration/dataStructures/moench03T1ReceiverDataNew.h

@@ -52,6 +52,7 @@ class moench03T1ReceiverDataNew : public slsDetectorData<uint16_t> {
   int sc_height;
   const int nSamples;
 
+  double ghost[200][25];
 
  public:
 
@@ -106,6 +107,11 @@ class moench03T1ReceiverDataNew : public slsDetectorData<uint16_t> {
 	}
       }
     }
+    // double ghost[200][25];
+
+    for (int ix=0; ix<25; ix++)
+      for (int iy=0; iy<200; iy++)
+	ghost[iy][ix]=0.;
 
     int ipacket;
     int ibyte;
@@ -143,6 +149,73 @@ class moench03T1ReceiverDataNew : public slsDetectorData<uint16_t> {
     
 
  
+  /**
+     Returns the value of the selected channel for the given dataset as double.
+     \param data pointer to the dataset (including headers etc)
+     \param ix pixel number in the x direction
+     \param iy pixel number in the y direction
+     \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;*/
+    /* double val=0, vout=getChannel(data, ix, iy); */
+    /* int x1=ix%25; */
+    /* for (int ix=0; ix<16; ix++) { */
+    /*   for (int ii=0; ii<2; ii++) { */
+    /* 	val+=getChannel(data,x1+25*ix,iy); */
+    /* 	val+=getChannel(data,x1+25*ix,399-iy); */
+    /*   } */
+    /* } */
+    /* vout+=0.0008*val-6224; */
+    /* return vout; //(double)getChannel(data, ix, iy);
+     */
+    uint16_t val=getChannel(data, ix, iy)&0x3fff;
+    return val;
+  };
+
+   
+
+  virtual void calcGhost(char *data, int ix, int iy) {
+    double val=0;
+    /* for (int ix=0; ix<25; ix++){ */
+    /*   for (int iy=0; iy<200; iy++) { */
+	val=0;
+	//	cout << "** "; 
+	for (int isc=0; isc<16; isc++) {
+	  //  for (int ii=0; ii<2; ii++) {
+	   val+=getChannel(data,ix+25*isc,iy);
+	   //  cout << "(" << isc << "," << val << " " ;
+	    val+=getChannel(data,ix+25*isc,399-iy);
+	    // cout << val << " " ;
+	    // }
+	}
+	ghost[iy][ix]=val;//-6224;
+	//	cout << " --"<< endl;
+    /*   }  */
+    /* } */
+    // cout << "*" << endl;
+    
+  }
+
+
+
+  virtual void calcGhost(char *data) {
+    for (int ix=0; ix<25; ix++){
+      for (int iy=0; iy<200; iy++) {
+	calcGhost(data, ix,iy);
+      } 
+    }
+    // cout << "*" << endl;
+  }
+
+
+  double getGhost(int ix, int iy) {
+    if (iy<200) return ghost[iy][ix%25]; 
+    if (iy<400) return ghost[399-iy][ix%25];
+    return 0;
+  };
+
      /**
 
      Returns the frame number for the given dataset. Purely virtual func.
@@ -278,8 +351,6 @@ class moench03T1ReceiverDataNew : public slsDetectorData<uint16_t> {
 
   }
 
-
-
   //int getPacketNumber(int x, int y) {return dataMap[y][x]/packetSize;};
 
 };

slsDetectorCalibration/dataStructures/moench03T1ZmqDataNew.h

@@ -14,6 +14,9 @@ class moench03T1ZmqDataNew : public slsDetectorData<uint16_t> {
   const int nSamples;
   const int offset;
 
+  double ghost[200][25];
+  double xtalk;
+
  public:
 
 
@@ -25,7 +28,7 @@ class moench03T1ZmqDataNew : public slsDetectorData<uint16_t> {
      \param c crosstalk parameter for the output buffer
 
   */
- moench03T1ZmqDataNew(int ns=5000): slsDetectorData<uint16_t>(400, 400, ns*32*2+sizeof(int)), nSamples(ns), offset(sizeof(int)) {
+ moench03T1ZmqDataNew(int ns=5000): slsDetectorData<uint16_t>(400, 400, ns*32*2+sizeof(int)), nSamples(ns), offset(sizeof(int)), xtalk(0.00021) {
 
     int nadc=32;
     int sc_width=25;
@@ -100,12 +103,84 @@ class moench03T1ZmqDataNew : public slsDetectorData<uint16_t> {
     
 
     
+    for (int ix=0; ix<25; ix++)
+      for (int iy=0; iy<200; iy++)
+	ghost[iy][ix]=0.;
+    
     
     // iframe=0;
     //  cout << "data struct created" << endl;
   };
     
 
+  double getXTalk(){return xtalk;};
+  void setXTalk(double g) {xtalk=g;};
+  
+  /**
+     Returns the value of the selected channel for the given dataset as double.
+     \param data pointer to the dataset (including headers etc)
+     \param ix pixel number in the x direction
+     \param iy pixel number in the y direction
+     \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;*/
+    /* double val=0, vout=getChannel(data, ix, iy); */
+    /* int x1=ix%25; */
+    /* for (int ix=0; ix<16; ix++) { */
+    /*   for (int ii=0; ii<2; ii++) { */
+    /* 	val+=getChannel(data,x1+25*ix,iy); */
+    /* 	val+=getChannel(data,x1+25*ix,399-iy); */
+    /*   } */
+    /* } */
+    /* vout+=0.0008*val-6224; */
+    /* return vout; //(double)getChannel(data, ix, iy);
+     */
+    return ((double)getChannel(data, ix, iy))+xtalk*getGhost(iy,iy);
+  };
+
+
+
+  virtual void calcGhost(char *data, int ix, int iy) {
+    double val=0;
+    /* for (int ix=0; ix<25; ix++){ */
+    /*   for (int iy=0; iy<200; iy++) { */
+	val=0;
+	for (int isc=0; isc<16; isc++) {
+	  for (int ii=0; ii<2; ii++) {
+	   val+=getChannel(data,ix+25*isc,iy);
+	   //  cout << val << " " ;
+	    val+=getChannel(data,ix+25*isc,399-iy);
+	    //  cout << val << " " ;
+	  }
+	}
+	ghost[iy][ix]=val;//-6224;
+	//	cout << endl;
+    /*   }  */
+    /* } */
+    // cout << "*" << endl;
+    
+  }
+
+
+
+  virtual void calcGhost(char *data) {
+    for (int ix=0; ix<25; ix++){
+      for (int iy=0; iy<200; iy++) {
+	calcGhost(data, ix,iy);
+      } 
+    }
+    // cout << "*" << endl;
+  }
+
+
+  double getGhost(int ix, int iy) {
+    if (iy<200) return ghost[iy][ix%25]; 
+    if (iy<400) return ghost[399-iy][ix%25];
+    return 0;
+  };
+
 
      /**
 

slsDetectorCalibration/ghostSummation.h

@@ -0,0 +1,61 @@
+#ifndef GHOSTSUMMATION_H
+#define GHOSTSUMMATION_H
+
+#include <cmath>
+#include "slsDetectorData.h"
+
+template <class dataType> class ghostSummation {
+
+  /** @short virtual calss to handle ghosting*/ 
+
+ public:
+  
+  /** constructor
+      \param xt crosstalk
+  */
+ ghostSummation(slsDetectorData<dataType> *d, double xt) : xtalk(xt),det(d)  {
+    if (det)
+      det->getDetectorSize(nx,ny);
+    ghost=new double[nx*ny];
+  };
+
+  ghostSummation(ghostSummation *orig) {
+    xtalk=orig->xtalk; 
+    det=orig->det;
+    ghost=new double[nx*ny];
+    
+  }
+  ~ghostSummation() {delete [] ghost;};
+
+  virtual ghostSummation *Clone() {
+    new ghostSummation(this);
+  }
+
+  double getXTalk(){return xtalk;};
+  void setXTalk(double g) {xtalk=g;};
+
+  virtual double calcGhost(char *data, int ix, int iy=1){ghost[iy*nx+ix]=0;};
+
+  virtual void calcGhost(char *data){
+    for (int iy=0; iy<ny; iy++)
+      for (int ix=0; ix<nx; ix++)
+	ghost[iy*nx+ix]=calcGhost(data, ix, iy);
+  }
+  virtual double getGhost(int ix, int iy) {
+    if (ix<0 || ix>=nx || iy<0 || iy>=ny) return 0;
+    return ghost[iy*nx+ix];
+  }
+
+ protected:
+    double xtalk;
+    slsDetectorData<dataType> *det;
+    double *ghost;
+    int nx, ny;
+
+};
+
+
+
+
+
+#endif

slsDetectorCalibration/interpolatingDetector.h

@@ -46,8 +46,8 @@ class interpolatingDetector : public singlePhotonDetector {
 		       int sign=1, 
 		       commonModeSubtraction *cm=NULL,
 		       int nped=1000, 
-		       int nd=100, int nnx=-1, int nny=-1) : 
+		       int nd=100, int nnx=-1, int nny=-1, double *gm=NULL, ghostSummation<uint16_t> *gs=NULL) : 
-  singlePhotonDetector(d, 3,nsigma,sign, cm, nped, nd, nnx, nny) , interp(inte), id(0)  {
+  singlePhotonDetector(d, 3,nsigma,sign, cm, nped, nd, nnx, nny, gm, gs) , interp(inte), id(0)  {
     //cout << "**"<< xmin << " " << xmax << " " << ymin << " " << ymax << endl;
     fi=new pthread_mutex_t ;
 

slsDetectorCalibration/moench03CommonMode.h

@@ -1,43 +1,52 @@
 #ifndef MOENCH03COMMONMODE_H
 #define MOENCH03COMMONMODE_H
 
-#include "commonModeSubtraction.h"
+#include "commonModeSubtractionNew.h"
 
-class moench03CommonMode : public commonModeSubtraction {
+class commonModeSubtractionColumn:  public commonModeSubtraction{
+public:
+ commonModeSubtractionColumn(int nr=200) : commonModeSubtraction(800), rows(nr) {};
+  virtual int getROI(int ix, int iy){return ix+(iy/200)*400;};
+  
+   virtual void addToCommonMode(double val, int ix=0, int iy=0) {
+     if (ix<rows || ix>399-rows) {
+       int iroi=getROI(ix,iy);
+       if (iroi>=0 && iroi<nROI) {
+	 mean[iroi]+=val; 
+	 mean2[iroi]+=val*val; 
+	 nCm[iroi]++;
+       }
+     }
+    };
+
+ private:
+  int rows;
+};
+
+
+class commonModeSubtractionSuperColumn:  public commonModeSubtraction{
+public:
+commonModeSubtractionSuperColumn() : commonModeSubtraction(32) {};
+  virtual int getROI(int ix, int iy){  return ix/25+(iy/200)*16;};
+};
+
+
+class commonModeSubtractionHalf:  public commonModeSubtraction{
+public:
+commonModeSubtractionHalf() : commonModeSubtraction(2) {};
+  virtual int getROI(int ix, int iy){  (void) ix; return iy/200;};
+};
+
+
+class moench03CommonMode : public commonModeSubtractionColumn {
   /** @short class to calculate the common mode noise for moench02 i.e. on 4 supercolumns separately */
  public:
   /** constructor -  initalizes a commonModeSubtraction with 4 different regions of interest 
       \param nn number of samples for the moving average
   */
-  moench03CommonMode(int nn=1000) : commonModeSubtraction(nn,32){} ;
+  moench03CommonMode(int nr=20) : commonModeSubtractionColumn(nr){} ;
     
 
-    /** add value to common mode as a function of the pixel value, subdividing the region of interest in the 4 supercolumns of 40 columns each;
-	\param val value to add to the common mode
-	\param ix pixel coordinate in the x direction
-	\param iy pixel coordinate in the y direction
-    */
-    virtual void addToCommonMode(double val, int ix=0, int iy=0) {
-      // (void) iy;  
-      int isc=ix/25+(iy/200)*16;
-      if (isc>=0 && isc<nROI) {
-	cmPed[isc]+=val; 
-	nCm[isc]++;
-      }
-    };
-     /**returns common mode value as a function of the pixel value, subdividing the region of interest in the 4 supercolumns of 40 columns each;
-	\param ix pixel coordinate in the x direction
-	\param iy pixel coordinate in the y direction
-	\returns common mode value
-    */
-    virtual double getCommonMode(int ix=0, int iy=0) {
-      (void) iy;
-      int isc=ix/25+(iy/200)*16;
-      if (isc>=0 && isc<nROI) { 
-	if (nCm[isc]>0) return cmPed[isc]/nCm[isc]-cmStat[isc].Mean(); 
-      }
-      return 0;
-    };
 
 };
 

slsDetectorCalibration/moench03GhostSummation.h

@@ -0,0 +1,57 @@
+#ifndef MOENCH03GHOSTSUMMATION_H
+#define MOENCH03GHOSTSUMMATION_H
+
+#include "ghostSummation.h"
+
+class moench03GhostSummation : public ghostSummation<uint16_t> {
+
+  /** @short virtual calss to handle ghosting*/ 
+
+ public:
+  
+  /** constructor
+      \param xt crosstalk
+  */
+ moench03GhostSummation(slsDetectorData<uint16_t> *d, double xt) : ghostSummation<uint16_t>(d, xt)  {}
+  
+  virtual void calcGhost(char *data){
+    for (int iy=0; iy<200; iy++){
+      for (int ix=0; ix<25; ix++){
+	calcGhost(data,ix,iy);
+      }
+    }
+  };
+
+ 
+  virtual double calcGhost(char *data, int x, int y=0){
+    int ix=x%25;
+    int iy=y;
+    if (y>=200) iy=399-y;
+    if (iy<0 || ix<0) return 0;
+    double val=0;
+     val=0;
+     for (int isc=0; isc<16; isc++) {
+	 val+=det->getChannel(data,ix+25*isc,iy);
+	 //  cout << val << " " ;
+	 val+=det->getChannel(data,ix+25*isc,399-iy);
+	    //  cout << val << " " ;
+	}
+	ghost[iy*nx+ix]=xtalk*val;
+	return ghost[iy*nx+ix];
+  };
+
+
+  virtual double getGhost(int ix, int iy) {
+    if (iy >=0 && iy<200) return ghost[iy*nx+(ix%25)]; 
+    if (iy<400) return ghost[(399-iy)*nx+(ix%25)];
+    return 0;
+  };
+
+
+};
+
+
+
+
+
+#endif

slsDetectorCalibration/moenchExecutables/moenchPhotonCounter.cpp

@@ -1,6 +1,6 @@
 //#include "ansi.h"
 #include <iostream>
-
+#define CORR
 
 //#define VERSION_V1
 
@@ -28,8 +28,11 @@
 //#include "etaInterpolationPosXY.h"
 // #include "linearInterpolation.h"
 // #include "noInterpolation.h"
-#include "multiThreadedAnalogDetector.h"
+#include "multiThreadedCountingDetector.h"
+//#include "multiThreadedAnalogDetector.h"
 #include "singlePhotonDetector.h"
+#include "moench03GhostSummation.h"
+#include "moench03CommonMode.h"
 //#include "interpolatingDetector.h"
 
 #include <stdio.h>
@@ -89,7 +92,18 @@ int main(int argc, char *argv[]) {
 
   decoder->getDetectorSize(nx,ny);
 
-  singlePhotonDetector *filter=new singlePhotonDetector(decoder,csize, nsigma, 1, 0, nped, 200);
+  int ncol_cm=20;
+  double xt_ghost=0.00045;
+  moench03CommonMode *cm=NULL;
+  moench03GhostSummation *gs;
+  double *gainmap=NULL;
+#ifdef CORR
+  cout << "Applying common mode and ghost correction " << endl;
+  cm=new moench03CommonMode(ncol_cm);
+ gs=new moench03GhostSummation(decoder, xt_ghost);
+#endif
+  
+  singlePhotonDetector *filter=new singlePhotonDetector(decoder,csize, nsigma, 1, cm, nped, 200, -1, -1, gainmap, gs);
 
   int size = 327680;////atoi(argv[3]);
   
@@ -171,9 +185,10 @@ int main(int argc, char *argv[]) {
   if (thr>0) {
     cout << "threshold is " << thr << endl;
     
-    //#ifndef ANALOG
+#ifndef ANALOG
     filter->setThreshold(thr);
-    //#endif
+#endif
+      cf=0;
 
   } else
     cf=1;
@@ -181,17 +196,13 @@ int main(int argc, char *argv[]) {
 
 
   filter->setROI(xmin,xmax,ymin,ymax);
-#ifdef SOLEIL
-  filter->setROI(150,210,170,230);
-  nframes=-1;
-#endif
   std::time(&end_time);
   cout << std::ctime(&end_time) <<   endl;
 
   char* buff;
 
-  multiThreadedAnalogDetector *mt=new multiThreadedAnalogDetector(filter,nthreads,fifosize);
+  // multiThreadedAnalogDetector *mt=new multiThreadedAnalogDetector(filter,nthreads,fifosize);
-
+ multiThreadedCountingDetector *mt=new multiThreadedCountingDetector(filter,nthreads,fifosize);
 #ifndef ANALOG
     mt->setDetectorMode(ePhotonCounting);
     cout << "Counting!" << endl;
@@ -204,7 +215,7 @@ int main(int argc, char *argv[]) {
       mt->setDetectorMode(eAnalog);
       cout << "Analog!" << endl;
       cf=0;
-      // thr1=thr;
+      thr1=thr;
 #endif  
       //  }
   
@@ -240,7 +251,7 @@ int main(int argc, char *argv[]) {
 	  mt->nextThread();
 	  mt->popFree(buff);
 	  ifr++;
-	  if (ifr%10000==0) 
+	  if (ifr%100==0) 
 	    cout << ifr << " " << ff << " " << np << endl;
 	} else
 	  cout << ifr << " " << ff << " " << np << endl;
@@ -308,7 +319,7 @@ int main(int argc, char *argv[]) {
 	  // // 		//	cout << " " << (void*)buff;
 	  mt->popFree(buff);
 	  ifr++;
-	  if (ifr%1000==0) cout << ifr << " " << ff << endl;
+	  if (ifr%100==0) cout << ifr << " " << ff << endl;
 	  if (nframes>0) {
 	    if (ifr%nframes==0) {
 	      //The name has an additional "_fXXXXX" at the end, where "XXXXX" is the initial frame number of the image (0,1000,2000...)
@@ -334,9 +345,10 @@ int main(int argc, char *argv[]) {
 	if (nframes>0) {
 	    sprintf(ffname,"%s/%s_f%05d.tiff",outdir,fformat,ifile);
 	    sprintf(imgfname,ffname,irun);
-	}
+	} else {
 	  sprintf(ffname,"%s/%s.tiff",outdir,fformat);
 	  sprintf(imgfname,ffname,irun);
+	}
 	cout << "Writing tiff to " << imgfname << " " << thr1 <<endl;
 	mt->writeImage(imgfname, thr1);
 	mt->clearImage();
@@ -351,6 +363,13 @@ int main(int argc, char *argv[]) {
     } else 
       cout << "Could not open "<< fname << " for reading " << endl;
   }
+  if (nframes<0){
+    sprintf(ffname,"%s/%s.tiff",outdir,fformat);
+    strcpy(imgfname,ffname);
+    cout << "Writing tiff to " << imgfname << " " << thr1 <<endl;
+    mt->writeImage(imgfname, thr1);
+  }
+    
 
 
   return 0;

slsDetectorCalibration/moenchExecutables/moenchZmqProcess.cpp

@@ -3,13 +3,15 @@
 #include "sls_detector_defs.h"
 #include "ZmqSocket.h"
 #include "moench03T1ZmqDataNew.h"
+#include "moench03GhostSummation.h"
+#include "moench03CommonMode.h"
 #include <vector> 
 #include <string> 
 #include <sstream> 
 #include <iomanip> 
 #include <fstream> 
 #include "tiffIO.h"
-//#include <zmq.h>
+
 #include <rapidjson/document.h> //json header in zmq stream
 
 #include<iostream>
@@ -111,11 +113,15 @@ int main(int argc, char *argv[]) {
 
   char dummybuff[size];
   
-
+  int ncol_cm=20;
+  double xt_ghost=0.00045;
+  moench03CommonMode *cm=new moench03CommonMode(ncol_cm);
+  moench03GhostSummation *gs=new moench03GhostSummation(det, xt_ghost);
+  double *gainmap=NULL;
    
   //analogDetector<uint16_t> *filter=new analogDetector<uint16_t>(det,1,NULL,1000);
 #ifndef INTERP
-  singlePhotonDetector *filter=new singlePhotonDetector(det,3, 5, 1, 0, 1000, 10);
+  singlePhotonDetector *filter=new singlePhotonDetector(det,3, 5, 1, cm, 1000, 10, -1, -1, gainmap, gs);
 
     multiThreadedCountingDetector *mt=new multiThreadedCountingDetector(filter,nthreads,fifosize);
 
@@ -126,7 +132,7 @@ int main(int argc, char *argv[]) {
 
   if (etafname) interp->readFlatField(etafname);
 
-  interpolatingDetector *filter=new interpolatingDetector(det,interp, 5, 1, 0, 1000, 10);
+  interpolatingDetector *filter=new interpolatingDetector(det,interp, 5, 1, cm,  1000, 10, -1, -1, gainmap, gs);
   multiThreadedInterpolatingDetector *mt=new multiThreadedInterpolatingDetector(filter,nthreads,fifosize);
 #endif
 
@@ -654,7 +660,12 @@ int main(int argc, char *argv[]) {
 	    // cout << "file" << endl;
 	    //  cout << "data " << endl;
 	  if (of==NULL) {
+#ifdef WRITE_QUAD
+	    sprintf(ofname,"%s_%d.clust2",filename.c_str(),fileindex);
+#endif
+#ifndef WRITE_QUAD
 	    sprintf(ofname,"%s_%d.clust",filename.c_str(),fileindex);
+#endif
 	    of=fopen(ofname,"w");
 	    if (of) {
 	      mt->setFilePointer(of);

slsDetectorCalibration/singlePhotonDetector.h

@@ -58,7 +58,7 @@ public analogDetector<uint16_t> {
 		      int sign=1,
 		      commonModeSubtraction *cm=NULL,
 		      int nped=1000,
-		      int nd=100, int nnx=-1, int nny=-1, double *gm=NULL) : analogDetector<uint16_t>(d, sign, cm, nped, nnx, nny, gm),   nDark(nd), eventMask(NULL),nSigma (nsigma), clusterSize(csize), clusterSizeY(csize), clusters(NULL),   quad(UNDEFINED_QUADRANT), tot(0), quadTot(0), eMin(-1), eMax(-1) {
+		      int nd=100, int nnx=-1, int nny=-1, double *gm=NULL, ghostSummation<uint16_t> *gs=NULL) : analogDetector<uint16_t>(d, sign, cm, nped, nnx, nny, gm, gs),   nDark(nd), eventMask(NULL),nSigma (nsigma), clusterSize(csize), clusterSizeY(csize), clusters(NULL),   quad(UNDEFINED_QUADRANT), tot(0), quadTot(0), eMin(-1), eMax(-1) {
     
     
     
@@ -213,7 +213,7 @@ public analogDetector<uint16_t> {
 	return nph;
       }	else {
 	if (thr>0) {
-	  newFrame();
+	  newFrame(data);
 	  if (cmSub) {
 	    cout << "add to common mode?"<< endl;
 	    addToCommonMode(data);
@@ -362,7 +362,7 @@ int *getClusters(char *data,  int *ph=NULL) {
     addToPedestal(data);
     return 0;
   }
-  newFrame();
+  newFrame(data);