From 61c495f218fa75eda995375f0dd679e9ca4a17d2 Mon Sep 17 00:00:00 2001
From: Anna Bergamaschi <anna.bergamaschi@psi.ch>
Date: Wed, 18 Sep 2019 13:21:47 +0200
Subject: [PATCH] multithreading fixed for gain map, common mode and ghosting

---
 slsDetectorCalibration/MovingStat.h           |   8 +-
 slsDetectorCalibration/analogDetector.h       | 122 +++++++++----
 .../commonModeSubtractionNew.h                |   3 +
 slsDetectorCalibration/ghostSummation.h       |   1 +
 slsDetectorCalibration/moench03CommonMode.h   |  11 +-
 .../moench03GhostSummation.h                  |  27 ++-
 .../moenchExecutables/moenchPhotonCounter.cpp | 172 +++++++++++++-----
 .../moenchExecutables/moenchZmqProcess.cpp    |  38 +++-
 .../multiThreadedAnalogDetector.h             |   5 +-
 slsDetectorCalibration/singlePhotonDetector.h |   2 +-
 10 files changed, 283 insertions(+), 106 deletions(-)

diff --git a/slsDetectorCalibration/MovingStat.h b/slsDetectorCalibration/MovingStat.h
index 35249b1c9..46917c98e 100755
--- a/slsDetectorCalibration/MovingStat.h
+++ b/slsDetectorCalibration/MovingStat.h
@@ -31,18 +31,19 @@ class MovingStat
 	*/
         void Set(double val, double rms=0, int m=-1)
         {
-	  if (m>=0) m_n = m;  else m_n = n; 
+	  if (m>0) m_n = m;  else m_n = n; 
 	  m_newM=val*m_n; 
+	  // cout << "set " << val << " " << m << " " << m_n << " " << m_newM << endl;
 	  SetRMS(rms);
         }
-		/**
+	/**
 	   clears the moving average number of samples parameter, mean and standard deviation
 	*/
         void SetRMS(double rms)
         {
 	  if (rms<=0) {
 	    m_newM2=m_newM*m_newM/n;
-	    m_n=0;
+	    //m_n=0;
 	  } else {
 	    if (m_n>0)
 	      m_newM2=(m_n*rms*rms+m_newM*m_newM/m_n);
@@ -120,6 +121,7 @@ class MovingStat
 	 */
         inline double Mean() const
         {
+	  //  cout << "get " <<  m_n << " " << m_newM << " " << m_newM/m_n << endl;
             return (m_n > 0) ? m_newM/m_n : 0.0;
         }
 
diff --git a/slsDetectorCalibration/analogDetector.h b/slsDetectorCalibration/analogDetector.h
index 66fec8ceb..d8b0e141f 100644
--- a/slsDetectorCalibration/analogDetector.h
+++ b/slsDetectorCalibration/analogDetector.h
@@ -154,8 +154,15 @@ 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();
+     if (orig->cmSub) {
+       cmSub=(orig->cmSub)->Clone();
+       cout <<"cloning cm" << endl;
+     }
+     else cmSub=NULL;
+     if (orig->ghSum) {
+       ghSum=(orig->ghSum)->Clone();
+       cout <<"cloning gs" << endl;
+     }
      else ghSum=NULL;
   }
 
@@ -224,9 +231,10 @@ template <class dataType> class analogDetector {
     if (gm) {
       if (gmap) delete [] gmap; 
       gmap=new double[nnx*nny];
-	for (int iy=0; iy<nny; iy++) {
-      for (int ix=0; ix<nnx; ix++) {
+      for (int iy=0; iy<nny; iy++) {
+	for (int ix=0; ix<nnx; ix++) {
 	  gmap[iy*nnx+ix]=gm[iy*nnx+ix];
+	  // cout << gmap[iy*nnx+ix] << " " ;
 	}
       }
       return gmap;
@@ -280,7 +288,12 @@ template <class dataType> class analogDetector {
   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);};
+  virtual void newFrame(char *data){
+    iframe++; 
+    if (cmSub) cmSub->newFrame(); 
+    calcGhost(data); 
+    //    cout << getId() << " Calc ghost " <<  getGhost(15,15) << endl;
+      };
   
 
     /** sets the commonModeSubtraction algorithm to be used 
@@ -315,7 +328,7 @@ template <class dataType> class analogDetector {
        \param iy pixel y coordinate
        \param cm 1 adds the value to common mod, 0 skips it. Defaults to 0. - not properly implemented
     */
-    virtual void addToPedestal(double val, int ix, int iy=0, int cm=0){ 
+    virtual void addToPedestal(double val, int ix, int iy, int cm=0){ 
       if (ix>=0 && ix<nx && iy>=0 && iy<ny) { 
 	
 	//	cout << val << " " ;
@@ -335,15 +348,17 @@ template <class dataType> class analogDetector {
     }
     
     double getCommonMode(int ix, int iy) {
-      if (cmSub) return cmSub->getCommonMode(ix, iy);
-      else return 0;
+      if (cmSub) {
+	return cmSub->getCommonMode(ix, iy);
+      }
+      return 0;
     }
 
 
     virtual void addToCommonMode(char *data){ 
-      // cout << "+"<< endl;
+      // cout << "+"<< getId() << endl;
       if (cmSub) {
-	//	cout << "*" << endl;
+	//cout << "*" << endl;
 	  for (int iy=ymin; iy<ymax; iy++) { 
 	for (int ix=xmin; ix<xmax; ix++) {
 	    // if (getNumpedestals(ix,iy)>0) 
@@ -358,12 +373,22 @@ template <class dataType> class analogDetector {
     }      
 
     virtual void addToCommonMode(char *data, int ix, int iy=0){ 
+      //  cout <<".";
       if (cmSub) {
+	//cout <<":";
 	if (det) if (det->isGood(ix, iy)==0) return;
 	if (getNumpedestals(ix,iy)>0){
-	    cmSub->addToCommonMode(subtractPedestal(data,ix,iy,0), ix, iy);
-	    //  cout << ":";
-	    // cout << ix << " " <<" " << iy << subtractPedestal(data,ix,iy,0)  << endl;
+	  double val;
+	  if (det) {
+	  /* if (det->getChannel(data, ix, iy)>=0x3fff) */
+	  /*   cout << ix << " " << iy << " " << det->getChannel(data, ix, iy) <<endl; */
+	    val= (dataSign*det->getValue(data, ix, iy)-getPedestal(ix,iy,0));
+	  }	else
+	    val= (((double*)data)[iy*nx+ix]-getPedestal(ix,iy));
+	  	val+=getGhost(ix,iy);
+	
+	  cmSub->addToCommonMode(val, ix, iy);
+	  //cout << ":";
 	}
       }
     }      
@@ -375,11 +400,12 @@ template <class dataType> class analogDetector {
        \returns pedestal value
     */
     virtual double getPedestal (int ix, int iy, int cm=0){
-      if (ix>=0 && ix<nx && iy>=0 && iy<ny) 
-	if (cmSub && cm>0) 
+      if (ix>=0 && ix<nx && iy>=0 && iy<ny) {
+	if (cmSub && cm>0) {
 	  return stat[iy][ix].getPedestal()+getCommonMode(ix,iy); 
-	else return stat[iy][ix].getPedestal(); 
-      else return -1;
+	  
+	}     	else return stat[iy][ix].getPedestal(); 
+      } else return -1;
     };
 
 
@@ -396,6 +422,7 @@ template <class dataType> class analogDetector {
 	  g=gmap[iy*nx+ix];
 	  if (g==0) g=-1.;
 	}
+      
 	return stat[iy][ix].getPedestalRMS()/g;//divide by gain?
       }
       return -1;
@@ -474,13 +501,14 @@ template <class dataType> class analogDetector {
     */
     virtual void setPedestal(double *ped, double *rms=NULL, int m=-1){
       double rr=0;
-	for (int iy=ymin; iy<ymax; iy++) {
-      for (int ix=xmin; ix<xmax; ix++) {
+      for (int iy=ymin; iy<ymax; iy++) {
+	for (int ix=xmin; ix<xmax; ix++) {
 	  if (rms) rr=rms[iy*nx+ix];
 	  stat[iy][ix].setPedestal(ped[iy*nx+ix],rr, m);
+	  //  cout << ix << " " << iy << " " << ped[iy*nx+ix] << " " << stat[iy][ix].getPedestal() << endl;
 	};
       };
-
+      
     }
     
 
@@ -511,10 +539,24 @@ 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, 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;};
+ virtual void calcGhost(char *data){if (ghSum) {
+     ghSum->calcGhost(data);
+     
+     //  cout << getId() << " @ " << ghSum->getXTalk() << " " <<  ghSum->getGhost(15,15) << endl;
+   }
+     
+
+};
+
+
+ virtual double getGhost(int ix, int iy) {if (ghSum) return ghSum->getGhost(ix, iy); return 0;};
 
      /**
        write 32bit tiff file with detector image data
@@ -722,20 +764,21 @@ template <class dataType> class analogDetector {
     
     virtual void addToPedestal(char *data, int cm=0) {
       
-      //    cout << "add to pedestal " << endl;
+      // cout << "add to pedestal " << endl;
       newFrame(data);
       
       //calcGhost(data);
 
 
       if (cmSub && cm) {
+	//	cout <<",";
 	addToCommonMode(data);
-      }
+      } 
             
       //cout << xmin << " " << xmax << endl;
       // cout << ymin << " " << ymax << endl;
-	for (int iy=ymin; iy<ymax; iy++) {
-      for (int ix=xmin; ix<xmax; ix++) {
+      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,cm);
@@ -818,7 +861,7 @@ template <class dataType> class analogDetector {
     */
 
 
-    virtual void addToPedestal(char *data, int ix, int iy=0, int cm=0) {
+    virtual void addToPedestal(char *data, int ix, int iy, int cm=0) {
       
       
       double val;
@@ -899,7 +942,10 @@ template <class dataType> class analogDetector {
 	  val= (dataSign*det->getValue(data, ix, iy)-getPedestal(ix,iy,cm))/g;
 	}	else
 	  val= (((double*)data)[iy*nx+ix]-getPedestal(ix,iy))/g;
-    
+	
+	//if (val>=0.5*thr)
+	// cout << val << " " << (dataSign*det->getValue(data, ix, iy)-getPedestal(ix,iy,cm)) << " " << g << " ";
+
 	val+=getGhost(ix,iy)/g;
 
 #ifdef ROOTSPECTRUM
@@ -961,8 +1007,7 @@ template <class dataType> class analogDetector {
      int nph=0;
      double v;
      if (ix>=0 && ix<nx && iy>=0 && iy<ny) {
-       v=subtractPedestal(data,ix,iy);
-
+       v=subtractPedestal(data,ix,iy,1);
 /*        // cout << v << " " ; */
 /* #ifdef ROOTSPECTRUM */
 /*        // cout << (iy-ymin)*(xmax-xmin)+(ix-xmin) << endl; */
@@ -971,8 +1016,14 @@ template <class dataType> class analogDetector {
        if (thr>0) {
 	 v+=0.5*thr;
 	 nph=v/thr;
-	 if (nph>0)
+	 
+	 /* if (ix ==10 && iy<20) */
+	 /*   cout << det->getValue(data,ix,iy) << " " << stat[iy][ix].getPedestal() << " " << getCommonMode(ix,iy) << " " << getPedestal(ix,iy,0)<< " "  << getPedestal(ix,iy,1) << endl; */
+	 // cout << subtractPedestal(data,ix,iy,0) << " " << subtractPedestal(data,ix,iy,1) << " " << nph << endl;
+	 if (nph>0) {
+	   //cout << " " << nph << endl;
 	   return nph;
+	 }
 	 return 0;
        } 
        return v;
@@ -991,12 +1042,13 @@ template <class dataType> class analogDetector {
        //double val;
        if (nph==NULL)
 	 nph=image;
+       
        newFrame(data);
 
        //calcGhost(data);
        addToCommonMode(data);
-	 for (int iy=ymin; iy<ymax; iy++) {
-       for (int ix=xmin; ix<xmax; ix++) {
+       for (int iy=ymin; iy<ymax; iy++) {
+	 for (int ix=xmin; ix<xmax; ix++) {
 	   if (det->isGood(ix,iy))
 	     nph[iy*nx+ix]+=getNPhotons(data, ix, iy);
 	 }
@@ -1101,7 +1153,7 @@ template <class dataType> class analogDetector {
       switch(fMode) {
       case ePedestal:
 	//cout << "analog ped " << endl;
-	addToPedestal(data);
+	addToPedestal(data,1);
 	break;
       default:
 	//	cout << "analog " << endl;
diff --git a/slsDetectorCalibration/commonModeSubtractionNew.h b/slsDetectorCalibration/commonModeSubtractionNew.h
index 791e599a3..1add20cc2 100644
--- a/slsDetectorCalibration/commonModeSubtractionNew.h
+++ b/slsDetectorCalibration/commonModeSubtractionNew.h
@@ -41,6 +41,7 @@ class commonModeSubtraction {
     
     /** adds the average of pedestals to the moving average and reinitializes the calculation of the sum of pedestals for all ROIs. - virtual func*/
     virtual void  newFrame(){
+      //cout << "Reset CM" << endl;
       for (int i=0; i<nROI; i++) {
 	//	if (nCm[i]>0) cmStat[i].Calc(cmPed[i]/nCm[i]); 
 	nCm[i]=0; 
@@ -60,6 +61,7 @@ class commonModeSubtraction {
       // else val=-100;
       // if (isc>=0 && isc<nROI) {
       //	cout << ix << " " << iy << " " << iroi << endl;
+      //if (ix==15 && iy==15) cout << "=" << val << endl;
       if (iroi>=0 && iroi<nROI) {
 	//	cout << ix << " " << iy << " " << iroi << endl;
 	mean[iroi]+=val; 
@@ -80,6 +82,7 @@ class commonModeSubtraction {
       /* else */
       /* 	return -100; */
       // cout << "*" << ix << " " << iy << " " << iroi << " " << mean[iroi] << " " << nCm[iroi]<< endl;
+      // if (ix==15 && iy==15) cout << "-" << mean[iroi]/nCm[iroi] << endl;
       if (iroi>=0 && iroi<nROI) {
 	if (nCm[iroi]>0)  
 	  return mean[iroi]/nCm[iroi]; 
diff --git a/slsDetectorCalibration/ghostSummation.h b/slsDetectorCalibration/ghostSummation.h
index 5546c75ee..a104855a7 100644
--- a/slsDetectorCalibration/ghostSummation.h
+++ b/slsDetectorCalibration/ghostSummation.h
@@ -41,6 +41,7 @@ template <class dataType> class ghostSummation {
       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];
diff --git a/slsDetectorCalibration/moench03CommonMode.h b/slsDetectorCalibration/moench03CommonMode.h
index f4e3b1f43..0d21c884f 100644
--- a/slsDetectorCalibration/moench03CommonMode.h
+++ b/slsDetectorCalibration/moench03CommonMode.h
@@ -9,16 +9,24 @@ public:
   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) {
+     if (iy<rows || iy>399-rows) {
        int iroi=getROI(ix,iy);
+       // cout << iy << " " << ix << " " << iroi ;
        if (iroi>=0 && iroi<nROI) {
 	 mean[iroi]+=val; 
 	 mean2[iroi]+=val*val; 
 	 nCm[iroi]++;
+	 if (nCm[iroi]>rows) cout << "Too many pixels added " << nCm[iroi] << endl;
+	 /* if (ix==10 && iy<20) */
+	 /*   cout << " ** "<<val << " " << mean[iroi] << " " << nCm[iroi] << " " << getCommonMode(ix, iy) << endl; */
        }
      }
     };
 
+  virtual  commonModeSubtractionColumn *Clone() {
+      return new  commonModeSubtractionColumn(this->rows);
+    }
+
  private:
   int rows;
 };
@@ -47,7 +55,6 @@ class moench03CommonMode : public commonModeSubtractionColumn {
   moench03CommonMode(int nr=20) : commonModeSubtractionColumn(nr){} ;
     
 
-
 };
 
 
diff --git a/slsDetectorCalibration/moench03GhostSummation.h b/slsDetectorCalibration/moench03GhostSummation.h
index 88377dbe8..165e422d5 100644
--- a/slsDetectorCalibration/moench03GhostSummation.h
+++ b/slsDetectorCalibration/moench03GhostSummation.h
@@ -22,22 +22,29 @@ class moench03GhostSummation : public ghostSummation<uint16_t> {
     }
   };
 
- 
+ moench03GhostSummation(moench03GhostSummation *orig) : ghostSummation(orig) {
+  }
+
+  virtual  moench03GhostSummation *Clone() {
+    return new  moench03GhostSummation(this);
+  }
+
   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];
+    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;
+    // if (ix==15 && iy==15) cout << ":" << ghost[iy*nx+ix] << " " << val << endl; 
+    return ghost[iy*nx+ix];
   };
 
 
diff --git a/slsDetectorCalibration/moenchExecutables/moenchPhotonCounter.cpp b/slsDetectorCalibration/moenchExecutables/moenchPhotonCounter.cpp
index b371f8be3..97737d9bb 100644
--- a/slsDetectorCalibration/moenchExecutables/moenchPhotonCounter.cpp
+++ b/slsDetectorCalibration/moenchExecutables/moenchPhotonCounter.cpp
@@ -2,6 +2,10 @@
 #include <iostream>
 #define CORR
 
+#define C_GHOST 0.0004
+
+#define CM_ROWS 50
+
 //#define VERSION_V1
 
 //#include "moench03T1ZmqData.h"
@@ -48,7 +52,7 @@ int main(int argc, char *argv[]) {
 
 
   if (argc<4) {
-    cout << "Usage is " << argv[0] << "indir outdir fname [runmin] [runmax] [pedfile] [threshold] [nframes] [xmin xmax ymin ymax]" << endl;
+    cout << "Usage is " << argv[0] << "indir outdir fname [runmin] [runmax] [pedfile] [threshold] [nframes] [xmin xmax ymin ymax] [gainmap]" << endl;
     cout << "threshold <0 means analog; threshold=0 means cluster finder; threshold>0 means photon counting" << endl;
     cout << "nframes <0 means sum everything; nframes=0 means one file per run; nframes>0 means one file every nframes" << endl;
     return 1;
@@ -56,7 +60,7 @@ int main(int argc, char *argv[]) {
 
   int p=10000;
   int fifosize=1000;
-  int nthreads=1;
+  int nthreads=10;
   int nsubpix=25;
   int etabins=nsubpix*10;
   double etamin=-1, etamax=2;
@@ -92,32 +96,24 @@ int main(int argc, char *argv[]) {
 
   decoder->getDetectorSize(nx,ny);
 
-  int ncol_cm=20;
-  double xt_ghost=0.0004;
+  int ncol_cm=CM_ROWS;
+  double xt_ghost=C_GHOST;
   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);
+  float *gm;
+
+
 
   int size = 327680;////atoi(argv[3]);
   
   int* image;
 	//int* image =new int[327680/sizeof(int)];
-  filter->newDataSet();
-
 
   int ff, np;
-  int dsize=decoder->getDataSize();
   //cout << " data size is " << dsize;
   
 
-  char data[dsize];
 
   ifstream filebin;
   char *indir=argv[1];
@@ -144,7 +140,7 @@ int main(int argc, char *argv[]) {
   double thr1=1;
 
   if (argc>=8) {
-    thr=atoi(argv[7]);
+    thr=atof(argv[7]);
   }
   
 
@@ -162,6 +158,12 @@ int main(int argc, char *argv[]) {
   }
 
   
+  char *gainfname=NULL;
+  if (argc>13) {
+    gainfname=argv[13];
+    cout << "Gain map file name is: " << gainfname << endl;
+  }
+
   
 
 
@@ -181,13 +183,59 @@ int main(int argc, char *argv[]) {
   cout << "runmax is " << runmax << endl;
   if (pedfile)
     cout << "pedestal file is " << pedfile << endl;
+ if (thr>0) 
+    cout << "threshold is " << thr << endl;
+ 
+ uint32 nnx, nny;
+ double *gmap;
+
+  // if (gainfname) {
+  //   gm=ReadFromTiff(gainfname, nny, nnx);
+  //   if (gm && nnx==nx && nny==ny) {
+  //     gmap=new double[nx*ny];
+  //     for (int i=0; i<nx*ny; i++) {
+  // 	gmap[i]=gm[i];
+  //     }
+  //     delete gm;
+  //   } else
+  //     cout << "Could not open gain map " << gainfname << endl;
+  // }
+
+#ifdef CORR
+ cout << "Applying common mode  " << ncol_cm << endl;
+  cm=new moench03CommonMode(ncol_cm);
+  
+
+  cout << "Applying ghost corrections " << xt_ghost << endl;
+  gs=new moench03GhostSummation(decoder, xt_ghost);
+#endif
+
+  singlePhotonDetector *filter=new singlePhotonDetector(decoder,csize, nsigma, 1, cm, nped, 200, -1, -1, gainmap, gs);
+
+  if (gainfname) {
+
+    if (filter->readGainMap(gainfname))
+      cout << "using gain map " << gainfname << endl;
+    else
+      cout << "Could not open gain map " << gainfname << endl;
+  } else
+    thr=0.15*thr;
+  filter->newDataSet();
+  int dsize=decoder->getDataSize();
+
+
+  char data[dsize];
+
+
+
+
+
   //#ifndef ANALOG
   if (thr>0) {
     cout << "threshold is " << thr << endl;
-    
-#ifndef ANALOG
+    //#ifndef ANALOG
     filter->setThreshold(thr);
-#endif
+    //#endif
       cf=0;
 
   } else
@@ -215,7 +263,7 @@ int main(int argc, char *argv[]) {
       mt->setDetectorMode(eAnalog);
       cout << "Analog!" << endl;
       cf=0;
-      thr1=thr;
+      //thr1=thr;
 #endif  
       //  }
   
@@ -226,45 +274,67 @@ int main(int argc, char *argv[]) {
   //  cout << "mt " << endl;
 
   int ifr=0;
- 
 
+  double ped[nx*ny], *ped1;
 
 
   if (pedfile) {
-    cout << "PEDESTAL " ;
-    sprintf(fname,"%s.raw",pedfile);
-    cout << fname << endl ;
-    sprintf(imgfname,"%s/pedestals.tiff",outdir,fformat);
-    std::time(&end_time);
-    cout << "aaa" << std::ctime(&end_time) <<   endl;
 
+    cout << "PEDESTAL " << endl;
+    sprintf(imgfname,"%s/pedestals.tiff",outdir);
 
-    mt->setFrameMode(ePedestal);
-      // sprintf(fn,fformat,irun);
-    filebin.open((const char *)(fname), ios::in | ios::binary);
-    //      //open file
-    if (filebin.is_open()){
-      ff=-1;
-      while (decoder->readNextFrame(filebin, ff, np,buff)) {
-	if (np==40) {
-	  mt->pushData(buff);
-	  mt->nextThread();
-	  mt->popFree(buff);
-	  ifr++;
-	  if (ifr%100==0) 
-	    cout << ifr << " " << ff << " " << np << endl;
-	} else
-	  cout << ifr << " " << ff << " " << np << endl;
-	ff=-1;
-      }
-      filebin.close();	 
-      while (mt->isBusy()) {;}
-      mt->writePedestal(imgfname);
+    if (string(pedfile).find(".tif")==std::string::npos){
+      sprintf(fname,"%s.raw",pedfile);
+      cout << fname << endl ;
       std::time(&end_time);
-      cout << std::ctime(&end_time) <<   endl;
+      cout << "aaa" << std::ctime(&end_time) <<   endl;
       
+      
+      mt->setFrameMode(ePedestal);
+      // sprintf(fn,fformat,irun);
+      filebin.open((const char *)(fname), ios::in | ios::binary);
+      //      //open file
+      if (filebin.is_open()){
+	ff=-1;
+	while (decoder->readNextFrame(filebin, ff, np,buff)) {
+	  if (np==40) {
+	    mt->pushData(buff);
+	    mt->nextThread();
+	    mt->popFree(buff);
+	    ifr++;
+	    if (ifr%100==0) 
+	      cout << ifr << " " << ff << " " << np << endl;
+	  } else
+	    cout << ifr << " " << ff << " " << np << endl;
+	  ff=-1;
+	}
+	filebin.close();	 
+	while (mt->isBusy()) {;}
+	
       } else 
 	cout << "Could not open pedestal file "<< fname << " for reading " << endl;
+    } else {
+      float *pp=ReadFromTiff(pedfile, nny, nnx);
+      if (pp && nnx==nx && nny==ny) {
+	for (int i=0; i<nx*ny; i++) {
+	  ped[i]=pp[i];
+	}
+	delete [] pp;
+	mt->setPedestal(ped);
+	// ped1=mt->getPedestal();
+	
+	// for (int i=0; i<nx*ny; i++) {
+
+	//   cout << ped[i]<<"/"<<ped1[i] << " " ;
+	// }
+	cout << "Pedestal set from tiff file " << pedfile << endl;
+      } else {
+	cout << "Could not open pedestal tiff file "<< pedfile << " for reading " << endl;
+      }
+    }
+    mt->writePedestal(imgfname);
+    std::time(&end_time);
+    cout << std::ctime(&end_time) <<   endl;
   }
   
 
@@ -318,6 +388,10 @@ int main(int argc, char *argv[]) {
 	  mt->nextThread();
 	  // // 		//	cout << " " << (void*)buff;
 	  mt->popFree(buff);
+
+	  while (mt->isBusy()) {;}
+
+
 	  ifr++;
 	  if (ifr%100==0) cout << ifr << " " << ff << endl;
 	  if (nframes>0) {
diff --git a/slsDetectorCalibration/moenchExecutables/moenchZmqProcess.cpp b/slsDetectorCalibration/moenchExecutables/moenchZmqProcess.cpp
index c1517b970..2cd8db82e 100644
--- a/slsDetectorCalibration/moenchExecutables/moenchZmqProcess.cpp
+++ b/slsDetectorCalibration/moenchExecutables/moenchZmqProcess.cpp
@@ -48,7 +48,7 @@ int main(int argc, char *argv[]) {
   double etamin=-1, etamax=2;
 	// help
   if (argc < 3 ) {
-    cprintf(RED, "Help: ./trial [receive socket ip] [receive starting port number] [send_socket ip] [send starting port number] [nthreads] [nsubpix] [etafile]\n");
+    cprintf(RED, "Help: ./trial [receive socket ip] [receive starting port number] [send_socket ip] [send starting port number] [nthreads] [nsubpix] [gainmap]  [etafile]\n");
     return EXIT_FAILURE;  
   }
   
@@ -90,12 +90,18 @@ int main(int argc, char *argv[]) {
     nSubPixels=atoi(argv[6]);
   cout << "Number of subpixels is: " << nSubPixels << endl;
   
-  char *etafname=NULL;
+
+  char *gainfname=NULL;
   if (argc>7) {
-    etafname=argv[7];
-    cout << "Eta file name is: " << etafname << endl;
+    gainfname=argv[7];
+    cout << "Gain map file name is: " << gainfname << endl;
   }
 
+  char *etafname=NULL;
+  if (argc>8) {
+    etafname=argv[8];
+    cout << "Eta file name is: " << etafname << endl;
+  }
 
   //slsDetectorData *det=new moench03T1ZmqDataNew(); 
   moench03T1ZmqDataNew *det=new moench03T1ZmqDataNew(); 
@@ -118,7 +124,29 @@ int main(int argc, char *argv[]) {
   moench03CommonMode *cm=new moench03CommonMode(ncol_cm);
   moench03GhostSummation *gs=new moench03GhostSummation(det, xt_ghost);
   double *gainmap=NULL;
-   
+  float *gm;
+
+
+  if (gainfname) {
+    gm=ReadFromTiff(gainfname, npy, npx);
+    if (gm) {
+      gmap=new double[npx*npy];
+      for (int i=0; i<npx*npy; i++) {
+	gmap[i]=gm[i];
+      }
+      delete [] gm;
+    } else
+      cout << "Could not open gain map " << gainfname << endl;
+  }
+
+
+
+
+
+
+
+
+
   //analogDetector<uint16_t> *filter=new analogDetector<uint16_t>(det,1,NULL,1000);
 #ifndef INTERP
   singlePhotonDetector *filter=new singlePhotonDetector(det,3, 5, 1, cm, 1000, 10, -1, -1, gainmap, gs);
diff --git a/slsDetectorCalibration/multiThreadedAnalogDetector.h b/slsDetectorCalibration/multiThreadedAnalogDetector.h
index d69fe033c..3f4d4982e 100644
--- a/slsDetectorCalibration/multiThreadedAnalogDetector.h
+++ b/slsDetectorCalibration/multiThreadedAnalogDetector.h
@@ -100,7 +100,10 @@ public:
 
    /** Returns true if the thread was successfully started, false if there was an error starting the thread */
    virtual bool StartThread()
-   { stop=0; 
+   { stop=0;  
+     cout << "Detector number " << det->getId() << endl; 
+     cout << "common mode is " << det->getCommonModeSubtraction()<< endl;
+     cout << "ghos summation is " << det->getGhostSummation()<< endl;
      return (pthread_create(&_thread, NULL, processData, this) == 0);
    }
 
diff --git a/slsDetectorCalibration/singlePhotonDetector.h b/slsDetectorCalibration/singlePhotonDetector.h
index 24e75fed6..99554ab54 100644
--- a/slsDetectorCalibration/singlePhotonDetector.h
+++ b/slsDetectorCalibration/singlePhotonDetector.h
@@ -592,7 +592,7 @@ int *getClusters(char *data,  int *ph=NULL) {
       switch(fMode) {
       case ePedestal:
 	//cout <<"spc add to ped " << endl;
-	addToPedestal(data);
+	addToPedestal(data,1);
 	break;
       default:
 	switch (dMode) {