fixed CM corrections

slsDetectorCalibration/analogDetector.h

@@ -866,7 +866,7 @@ template <class dataType> class analogDetector {
             for (ix = xmin; ix < xmax; ++ix) {
                 if (det->isGood(ix, iy)) {
                     // addToPedestal(data,ix,iy,1);
-                    addToPedestal(data, ix, iy, cm);
+                    addToPedestal(data, ix, iy, 0);
                     /* if (ix==50 && iy==50) */
                     /*   cout<< "*ped* " << id << " " << ix << " " << iy << " "
                      * << det->getChannel(data,ix,iy) << " " <<
@@ -1097,8 +1097,7 @@ template <class dataType> class analogDetector {
         return thr;
     };
 
-    virtual int setClusterSize(int n = -1) { 
-    };
+    virtual int setClusterSize(int n = -1) { return 0;     };
 
     /**
        gets threshold value for conversion into number of photons

slsDetectorCalibration/commonModeSubtractionNew.h

@@ -100,7 +100,10 @@ class commonModeSubtraction {
         if (iroi >= 0 && iroi < nROI) {
             if (nCm[iroi] > 0)
                 return mean[iroi] / nCm[iroi];
-        }
+	    else
+	      std::cout << ix << " " << iy << " no data in cm " << iroi << std::endl; 
+        } else
+	  std::cout << ix << " " << iy << " bad roi " << iroi << std::endl; 
         return 0;
     };
 
@@ -117,6 +120,8 @@ class commonModeSubtraction {
                 return sqrt(mean2[iroi] / nCm[iroi] -
                             (mean[iroi] / nCm[iroi]) *
                                 (mean[iroi] / nCm[iroi]));
+	    else
+	      std::cout << ix << " " << iy << " no data in cm " << iroi << std::endl; 
         }
         return 0;
     };

slsDetectorCalibration/dataStructures/jungfrauModuleData.h

@@ -3,6 +3,7 @@
 #ifndef JUNGFRAUMODULEDATA_H
 #define JUNGFRAUMODULEDATA_H
 #include "slsDetectorData.h"
+#include "sls/sls_detector_defs.h"
 
 //#define VERSION_V2
 /**
@@ -41,8 +42,15 @@ class jungfrauModuleData : public slsDetectorData<uint16_t> {
        1286 large etc.) \param c crosstalk parameter for the output buffer
 
    */
+   
+#ifdef ALDO //VH
+    using header = jf_header; //VH
+#else //VH
+    using header = sls::defs::sls_receiver_header;
+#endif //VH
+
 #ifndef ZMQ
-#define off sizeof(jf_header)
+#define off sizeof(header)
 #endif
 #ifdef ZMQ
 #define off 0
@@ -102,9 +110,13 @@ class jungfrauModuleData : public slsDetectorData<uint16_t> {
     /* 	unsigned char frameNumber[3]; */
     /* 	unsigned char bunchid[8]; */
     /* }; */
-
+    
     int getFrameNumber(char *buff) {
-        return ((jf_header *)buff)->bunchNumber;
+    #ifdef ALDO
+    return ((header *)buff)->bunchNumber; 
+    #else 
+    return ((header *)buff)->detHeader.frameNumber;
+    #endif 
     };
   
     /**

slsDetectorCalibration/interpolations/slsInterpolation.h

@@ -27,6 +27,8 @@ enum quadrant {
 //#ifdef MYROOT1
 //: public TObject
 //#endif
+using namespace std;
+
 class slsInterpolation {
 
   public:

slsDetectorCalibration/jungfrauCommonMode.h

@@ -0,0 +1,58 @@
+// SPDX-License-Identifier: LGPL-3.0-or-other
+// Copyright (C) 2021 Contributors to the SLS Detector Package
+#ifndef COMMONMODEJF_H
+#define COMMONMODEJF_H
+
+#include "commonModeSubtractionNew.h"
+
+class commonModeSubtractionSuperColumnJF : public commonModeSubtraction {
+  public:
+    commonModeSubtractionSuperColumnJF()
+        : commonModeSubtraction(32){};
+    virtual int getROI(int ix, int iy) { int top=iy/256; int sc=ix/64; return sc+top*16; };
+/*
+    virtual void addToCommonMode(double val, int ix = 0, int iy = 0) {
+            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] > 64*256)
+                    std::cout << "Too many pixels added " << nCm[iroi] << std::endl;
+           
+        }
+    };
+*/
+    virtual commonModeSubtractionSuperColumnJF *Clone() {
+        return new commonModeSubtractionSuperColumnJF();
+    };
+};
+
+
+class commonModeSubtractionChip : public commonModeSubtraction {
+  public:
+    commonModeSubtractionChip()
+        : commonModeSubtraction(8){};
+    virtual int getROI(int ix, int iy) { int top=iy/256; int sc=ix/256; return sc+top*4; };
+/*
+    virtual void addToCommonMode(double val, int ix = 0, int iy = 0) {
+            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] > 64*256)
+                    std::cout << "Too many pixels added " << nCm[iroi] << std::endl;
+           
+        }
+    };
+*/
+    virtual commonModeSubtractionChip *Clone() {
+        return new commonModeSubtractionChip();
+    };
+};
+
+
+#endif

slsDetectorCalibration/jungfrauExecutables/jungfrauRawDataProcess.cpp

@@ -2,8 +2,12 @@
 // Copyright (C) 2021 Contributors to the SLS Detector Package
 //#include "sls/ansi.h"
 #include <iostream>
-#undef CORR
 
+//enable common mode subtraction
+#define CMS
+//disable common mode subtraction
+//#undef CMS
+#undef CORR
 #define C_GHOST 0.0004
 
 #define CM_ROWS 50
@@ -36,15 +40,17 @@
 #include <map>
 #include <stdio.h>
 #include <sys/stat.h>
+//#include "commonModeSubtractionNew.h"
+#include "jungfrauCommonMode.h"
 
 #include <ctime>
 using namespace std;
 
 int main(int argc, char *argv[]) {
 
-    if (argc < 5) {
+    if (argc < 6) {
         cout << "Usage is " << argv[0]
-             << "indir outdir fname(no extension) fextension [runmin] [runmax] [pedfile (raw or tiff)] [threshold] "
+             << " indir outdir fname(no extension) fextension csize [runmin] [runmax] [pedfile (raw or tiff)] [threshold]"
                 "[nframes] [xmin xmax ymin ymax] [gainmap]"
              << endl;
         cout << "threshold <0 means analog; threshold=0 means cluster finder; "
@@ -53,11 +59,12 @@ int main(int argc, char *argv[]) {
         cout << "nframes <0 means sum everything; nframes=0 means one file per "
                 "run; nframes>0 means one file every nframes"
              << endl;
+
         return 1;
     }
 
     int fifosize = 1000;
-    int nthreads = 10;
+    int nthreads = 1;
     int csize = 3; //3
     int nsigma = 5;
     int nped = 10000;
@@ -67,10 +74,12 @@ int main(int argc, char *argv[]) {
 
 #if !defined JFSTRX && !defined JFSTRXOLD && !defined JFSTRXCHIP1 && !defined JFSTRXCHIP6
 #ifndef MODULE
+    cout << "This is a JF single chip!" <<endl;
     jungfrauHighZSingleChipData *decoder = new jungfrauHighZSingleChipData();
     int nx = 256, ny = 256;
 #endif
 #ifdef MODULE
+    cout << "This is a JF module!" <<endl;
     jungfrauModuleData *decoder = new jungfrauModuleData();
     int nx = 1024, ny = 512;
 #endif
@@ -80,7 +89,7 @@ int main(int argc, char *argv[]) {
     cout << "Jungfrau strixel full module readout" << endl;
     jungfrauLGADStrixelsData *decoder = new jungfrauLGADStrixelsData();
     int nx = 1024/5, ny = 512*5;
-#endif
+#endif 
 #ifdef JFSTRXCHIP1
     std::cout << "Jungfrau strixel LGAD single chip 1" << std::endl;
     jungfrauLGADStrixelsDataSingleChip *decoder = new jungfrauLGADStrixelsDataSingleChip(1);
@@ -99,20 +108,9 @@ int main(int argc, char *argv[]) {
 
 
     decoder->getDetectorSize(nx, ny);
-
-
-
-
-
+    
     cout << "Detector size is " << nx << " " << ny << endl;
 
-
-
-
-
-
-
-
     double *gainmap = NULL;
     //float *gm;
 
@@ -124,47 +122,48 @@ int main(int argc, char *argv[]) {
     char *outdir = argv[2];
     char *fformat = argv[3];
     char *fext = argv[4];
+    sscanf(argv[5], "%d", &csize);
     int runmin = 0;
 
     // cout << "argc is " << argc << endl;
-    if (argc >= 6) {
-        runmin = atoi(argv[5]);
+    if (argc >= 7) {
+        runmin = atoi(argv[6]);
     }
 
     int runmax = runmin;
 
-    if (argc >= 7) {
-        runmax = atoi(argv[6]);
+    if (argc >= 8) {
+        runmax = atoi(argv[7]);
     }
 
     char *pedfile = NULL;
-    if (argc >= 8) {
-        pedfile = argv[7];
+    if (argc >= 9) {
+        pedfile = argv[8];
     }
     double thr = 0;
     double thr1 = 1;
 
-    if (argc >= 9) {
-        thr = atof(argv[8]);
+    if (argc >= 10) {
+        thr = atof(argv[9]);
     }
 
     int nframes = 0;
 
-    if (argc >= 10) {
-        nframes = atoi(argv[9]);
+    if (argc >= 11) {
+        nframes = atoi(argv[10]);
     }
 
     int xmin = 0, xmax = nx, ymin = 0, ymax = ny;
-    if (argc >= 14) {
-        xmin = atoi(argv[10]);
-        xmax = atoi(argv[11]);
-        ymin = atoi(argv[12]);
-        ymax = atoi(argv[13]);
+    if (argc >= 15) {
+        xmin = atoi(argv[11]);
+        xmax = atoi(argv[12]);
+        ymin = atoi(argv[13]);
+        ymax = atoi(argv[14]);
     }
 
     char *gainfname = NULL;
-    if (argc > 14) {
-        gainfname = argv[14];
+    if (argc > 15) {
+        gainfname = argv[15];
         cout << "Gain map file name is: " << gainfname << endl;
     }
 
@@ -188,9 +187,13 @@ int main(int argc, char *argv[]) {
     cout << "Nframes is " << nframes << endl;
 
     uint32_t nnx, nny;
-
+    commonModeSubtraction *cm = NULL;
+#ifdef CMS
+    cm = new commonModeSubtractionChip();//commonModeSubtraction(1,5);//commonModeSubtractionSuperColumnJF();
+    std::cout << "Enabled common mode subtraction" << std::endl;
+#endif
     singlePhotonDetector *filter = new singlePhotonDetector(
-        decoder, 3, nsigma, 1, NULL, nped, 200, -1, -1, gainmap, NULL);
+        decoder, 3, nsigma, 1, cm, nped, 200, -1, -1, gainmap, NULL);
 
     if (gainfname) {
 
@@ -250,13 +253,14 @@ int main(int argc, char *argv[]) {
     double *ped=new double[nx * ny];//, *ped1;
 
     int pos,pos1;
+    int ipixX=351, ipixY=331;
 
     if (pedfile) {
 
-      if (string(pedfile).find(".dat") != std::string::npos) {
+      if (string(pedfile).find(".raw") != std::string::npos) {
 	pos1=string(pedfile).rfind("/");
 	strcpy(froot,pedfile+pos1);
-	pos=string(froot).find(".dat");
+	pos=string(froot).find(".raw");
 	froot[pos]='\0';
       }
 
@@ -277,19 +281,43 @@ int main(int argc, char *argv[]) {
                 ff = -1;
                 while (decoder->readNextFrame(filebin, ff, np, buff)) {
 		  // if (np == 40) {
-		  if ((ifr+1) % 100 == 0) {
-		    cout << " ****" << decoder->getValue(buff,20,20);// << endl;
-		  }
+		  // if ((ifr+1) % 100 == 0) {
+		  //   cout << " ****" << decoder->getValue(buff,20,20);// << endl;
+		  // }
                         mt->pushData(buff);
+
+
+			// while (mt->isBusy()) 
+			//   ;
+			//if ((ifr+1) % 100 == 0)
+
+			// double pix=0,ped=0,sub=0,cmp=0, sub0=0;
+
+			// for (int ix=-1; ix<2; ix++) {
+			//   for (int iy=-1; iy<2; iy++) {
+			//     pix+=decoder->getValue(buff, ipixX+ix, ipixY+iy);
+			//     sub+=filter->subtractPedestal(buff, ipixX+ix, ipixY+iy, 1);
+			//     sub0+=filter->subtractPedestal(buff, ipixX+ix, ipixY+iy, 0);
+			//     ped+=filter->getPedestal(ipixX+ix, ipixY+iy);
+			//     cmp+=filter->getCommonMode(ipixX+ix, ipixY+iy);
+			//   }
+			// }
+
+			// cout << pix << " " << sub << " " sub0 << " " << ped << " " << cmp  << endl;
+
+
                         mt->nextThread();
                         mt->popFree(buff);
                         ifr++;
-                        if (ifr % 100 == 0) {
-			  cout << " ****" << ifr << " " << ff << " " << np << endl;
-			} //else
+                        // if (ifr % 100 == 0) {
+			//   cout << " ****" << ifr << " " << ff << " " << np << endl;
+			// } 
+
+			//else
                         //cout << ifr << " " << ff << " " << np << endl;
-			if (ifr>=1000)
-			  break;
+			// if (ifr>=1000)
+			//   break;
+			
                     ff = -1;
                 }
                 filebin.close();
@@ -370,7 +398,25 @@ int main(int argc, char *argv[]) {
 		  if ((ifr+1) % 100 == 0) {
 		    cout << " ****" << decoder->getValue(buff,20,20);// << endl;
 		  }
-                    mt->pushData(buff);
+                    mt->pushData(buff);	
+
+
+		    // while (mt->isBusy()) 
+
+			// double pix=0,ped=0,sub=0,cmp=0, sub0=0;
+
+			// for (int ix=-1; ix<2; ix++) {
+			//   for (int iy=-1; iy<2; iy++) {
+			//     pix+=decoder->getValue(buff, ipixX+ix, ipixY+iy);
+			//     sub+=filter->subtractPedestal(buff, ipixX+ix, ipixY+iy, 1);
+			//     sub0+=filter->subtractPedestal(buff, ipixX+ix, ipixY+iy, 0);
+			//     ped+=filter->getPedestal(ipixX+ix, ipixY+iy);
+			//     cmp+=filter->getCommonMode(ipixX+ix, ipixY+iy);
+			//   }
+			// }
+
+			// cout << pix << " " << sub << " " sub0 << " " << ped << " " << cmp  << endl;
+
                     // 	//         //pop
                     mt->nextThread();
                     mt->popFree(buff);

slsDetectorCalibration/singlePhotonDetector.h

@@ -396,14 +396,14 @@ class singlePhotonDetector : public analogDetector<uint16_t> {
             return 0;
         }
         newFrame(data);
-
+	/*
         if (cmSub) {
             addToCommonMode(data);
             cm = 1;
         }
-
+	*/
         double *val = new double[ny * nx];
-
+	cm=0;
         for (iy = ymin; iy < ymax; ++iy) {
             for (ix = xmin; ix < xmax; ++ix) {
                 if (det->isGood(ix, iy) == 0)
@@ -547,12 +547,50 @@ class singlePhotonDetector : public analogDetector<uint16_t> {
 
                 } else if (ee == PEDESTAL) {
                     addToPedestal(data, ix, iy, cm);
+		    if (cmSub) {
+		      addToCommonMode(data, ix, iy);
+		    }
                 } /*else {
                     eventMask[iy][ix]=PHOTON;
                     }*/
                   // eventMask[iy][ix]=ee;
             }
         }
+	
+	if (cmSub) {
+	  
+	  double cmv=0;
+	  for (int iph=0; iph<nph; iph++) {
+	    ix=(clusters + iph)->x ;
+	    iy=(clusters + nph)->y;
+	    
+	    for (ir = -(clusterSizeY / 2); ir < (clusterSizeY / 2) + 1;
+		 ir++) {
+	      for (ic = -(clusterSize / 2);
+		   ic < (clusterSize / 2) + 1; ic++) {
+		if ((iy + ir) >= 0 && (iy + ir) < ny &&
+		    (ix + ic) >= 0 && (ix + ic) < nx) {
+		  
+		  cmv=getCommonMode(ix+ic, iy+ir);
+		  //if (ir==0 && ic==0)
+		    //cout << cmv << endl;
+		  (clusters + iph)
+		    ->set_data((clusters + iph)
+		    ->get_data(ic,ir)-cmv, ic,ir);
+		}
+		
+	      }
+	      
+	    }
+	    
+	  }
+	  
+	  
+	}
+	
+
+	
+
 
         nphFrame = nph;
         nphTot += nph;

slsDetectorCalibration/tiffio/src/tiffIO.cpp

@@ -3,7 +3,8 @@
 
 #include "sls/tiffIO.h"
 #include <iostream>
-#include <tiffio.h>
+#include "/usr/include/tiffio.h"
+//#include <tiffio.h>
 
 void *WriteToTiff(float *imgData, const char *imgname, int ncol, int nrow) {
     constexpr uint32_t sampleperpixel = 1;