Rosenblatt transform works

examples/rawdata.config

@@ -5,18 +5,18 @@ nsigma 5
 gainfile none
 detectorMode counting
 threshold 0
-pedestalfile none
+pedestalfile /mnt/sls_det_storage/moench_data/testNewFW20230714/cu_half_speed_d0_f0_4.raw
 nframes 0
 xMin 0
 xMax 400
 yMin 0
 yMax 400
-outdir ./
-indir ./
+outdir /mnt/sls_det_storage/moench_data/testNewFW20230714/clustNew/
+indir /mnt/sls_det_storage/moench_data/testNewFW20230714/
 flist none
-fformat none
+fformat cu_half_speed_d0_f0_4
 runmin 0
-runmax -1
+runmax 0
 readnrows 400
 eMin 0
 eMax 16000

slsDetectorCalibration/interpolations/etaInterpolationBase.h

@@ -200,6 +200,9 @@ class etaInterpolationBase : public slsInterpolation {
     // virtual void prepareInterpolation(int &ok)=0;
 
     void debugSaveAll(int ind = 0) {
+
+      std::cout << "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++SAVING ETAS"<< std::endl;
+      
         int ibx, iby;
         char tit[10000];
 
@@ -214,21 +217,21 @@ class etaInterpolationBase : public slsInterpolation {
             etah[ii] = heta[ii];
             tot_eta += heta[ii];
         }
-        sprintf(tit, "/scratch/eta_%d.tiff", ind);
+        sprintf(tit, "eta_%d.tiff", ind);
         WriteToTiff(etah, tit, nbetaX, nbetaY);
 
         for (int ii = 0; ii < nbetaX * nbetaY; ii++) {
             ibb = (hhx[ii] * nSubPixelsX);
             etah[ii] = ibb;
         }
-        sprintf(tit, "/scratch/eta_hhx_%d.tiff", ind);
+        sprintf(tit, "eta_hhx_%d.tiff", ind);
         WriteToTiff(etah, tit, nbetaX, nbetaY);
 
         for (int ii = 0; ii < nbetaX * nbetaY; ii++) {
             ibb = hhy[ii] * nSubPixelsY;
             etah[ii] = ibb;
         }
-        sprintf(tit, "/scratch/eta_hhy_%d.tiff", ind);
+        sprintf(tit, "eta_hhy_%d.tiff", ind);
         WriteToTiff(etah, tit, nbetaX, nbetaY);
 
         float *ftest = new float[nSubPixelsX * nSubPixelsY];
@@ -262,7 +265,7 @@ class etaInterpolationBase : public slsInterpolation {
                      << std::endl;
         }
 
-        sprintf(tit, "/scratch/ftest_%d.tiff", ind);
+        sprintf(tit, "./ftest_%d.tiff", ind);
         WriteToTiff(ftest, tit, nSubPixelsX, nSubPixelsY);
 
         // int ibx=0, iby=0;
@@ -280,7 +283,7 @@ class etaInterpolationBase : public slsInterpolation {
             } else
                 etah[ii] = 0;
         }
-        sprintf(tit, "/scratch/eta_bad_%d.tiff", ind);
+        sprintf(tit, "./eta_bad_%d.tiff", ind);
         WriteToTiff(etah, tit, nbetaX, nbetaY);
         // std::cout << "Index: " << ind << "\t Bad bins: "<< nbad << std::endl;
         // int ibx=0, iby=0;

slsDetectorCalibration/interpolations/etaInterpolationRosenblatt.h

@@ -0,0 +1,232 @@
+// SPDX-License-Identifier: LGPL-3.0-or-other
+// Copyright (C) 2021 Contributors to the SLS Detector Package
+#ifndef ETA_INTERPOLATION_POSXY_H
+#define ETA_INTERPOLATION_POSXY_H
+//#define SAVE_ALL
+//#include "sls/tiffIO.h"
+#include "eta2InterpolationBase.h"
+#include "eta3InterpolationBase.h"
+#include "etaInterpolationBase.h"
+
+
+class etaInterpolationRosenblatt : public virtual etaInterpolationBase {
+  public:
+    etaInterpolationRosenblatt(int nx = 400, int ny = 400, int ns = 25, int nsy = 25,
+                          int nb = -1, int nby = -1, double emin = 1,
+                          double emax = 0)
+        : etaInterpolationBase(nx, ny, ns, nsy, nb, nby, emin, emax){
+              //  std::cout << "epxy " << nb << " " << emin << " " << emax << std::endl;
+              //  std::cout << nbeta << " " << etamin << " " << etamax << std::endl;
+          };
+
+    etaInterpolationRosenblatt(etaInterpolationRosenblatt *orig)
+        : etaInterpolationBase(orig){};
+
+    virtual etaInterpolationRosenblatt *Clone() = 0; /* { */
+
+    /* return new etaInterpolationPosXY(this); */
+
+    /* }; */
+
+    virtual void prepareInterpolation(int &ok) {
+        ok = 1;
+
+        ///*Eta Distribution Rebinning*///
+        // double bsize=1./nSubPixels; //precision
+        // std::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 < nbetaX * nbetaY; ip++)
+            tot_eta += heta[ip];
+        std::cout << "total eta entries is :" << tot_eta << std::endl;
+        if (tot_eta <= 0) {
+            ok = 0;
+            return;
+        };
+
+        double *hx = new double[nbetaX];  // profile x
+        double *hy = new double[nbetaY];  // profile y
+        double *hix = new double[nbetaX]; // integral of projection x
+        double *hiy = new double[nbetaY]; // integral of projection y
+        //  int ii=0;
+        double etax, etay;
+        for (int ib = 0; ib < nbetaX; ib++) {
+
+            // tot_eta_y=0;
+
+            for (int iby = 0; iby < nbetaY; iby++) {
+                etay = etamin + iby * etastepY;
+                // std::cout << etax << std::endl;
+
+                // tot_eta_x+=hx[iby];
+                if (etay >= 0 && etay <= 1)
+                    hy[iby] = heta[ib + iby * nbetaX];
+                else
+                    hy[iby] = 0;
+                // tot_eta_y+=hy[iby];
+            }
+
+            hiy[0] = hy[0];
+
+            for (int iby = 1; iby < nbetaY; iby++) {
+                hiy[iby] = hiy[iby - 1] + hy[iby];
+            }
+
+            tot_eta_y = hiy[nbetaY - 1] + 1;
+
+            for (int iby = 0; iby < nbetaY; iby++) {
+                if (tot_eta_y <= 0) {
+                    hhy[ib + iby * nbetaX] = -1;
+                    // ii=(ibx*nSubPixels)/nbeta;
+                } else {
+                    // if (hiy[ibx]>tot_eta_y*(ii+1)/nSubPixels) ii++;
+                    hhy[ib + iby * nbetaX] = hiy[iby] / tot_eta_y;
+                }
+            }
+        }
+
+        for (int ib = 0; ib < nbetaY; ib++) {
+
+            for (int ibx = 0; ibx < nbetaX; ibx++) {
+                etax = etamin + ibx * etastepX;
+                // std::cout << etax << std::endl;
+                if (etax >= 0 && etax <= 1)
+                    hx[ibx] = heta[ibx + ib * nbetaX];
+                else {
+                    hx[ibx] = 0;
+                }
+            }
+            hix[0] = hx[0];
+
+            for (int ibx = 1; ibx < nbetaX; ibx++) {
+                hix[ibx] = hix[ibx - 1] + hx[ibx];
+            }
+
+	    // tot_eta_x = hix[nbetaX - 1] + 1;
+
+            for (int ibx = 0; ibx < nbetaX; ibx++) {
+	      //if (tot_eta_x <= 0) {
+              //      hhx[ibx + ib * nbetaX] = -1;
+              //  } else {
+	      hhx[ibx + ib * nbetaX] = hix[ibx];// / tot_eta_x;
+	      //}
+            }
+
+            
+        }
+
+
+	  for (int ibx = 0; ibx < nbetaX; ibx++) {
+	  int val=0;
+	for (int ib = 0; ib < nbetaY; ib++) {
+	    val+=hhx[ibx + ib * nbetaX];
+	  }
+	// for (int ibx = 0; ibx < nbetaX; ibx++) {
+	for (int ib = 0; ib < nbetaY; ib++) {
+	    hhx[ibx + ib * nbetaX]=val;
+          }    
+	}
+	tot_eta_x = hhx[nbetaX - 1] + 1;
+	for (int ib = 0; ib < nbetaY; ib++) {
+	    
+	  //tot_eta_x = hix[nbetaX - 1] + 1;
+
+            for (int ibx = 0; ibx < nbetaX; ibx++) {
+                if (tot_eta_x <= 0) {
+                    hhx[ibx + ib * nbetaX] = -1;
+                } else {
+                    hhx[ibx + ib * nbetaX] = hhx[ibx + ib * nbetaX] / tot_eta_x;
+                }
+            }
+	}
+	/*
+        int ibx, iby, ib;
+
+        iby = 0;
+        while (hhx[iby * nbetaY + nbetaY / 2] < 0)
+            iby++;
+        for (ib = 0; ib < iby; ib++) {
+            for (ibx = 0; ibx < nbetaX; ibx++)
+                hhx[ibx + nbetaX * ib] = hhx[ibx + nbetaX * iby];
+        }
+        iby = nbetaY - 1;
+
+        while (hhx[iby * nbetaY + nbetaY / 2] < 0)
+            iby--;
+        for (ib = iby + 1; ib < nbetaY; ib++) {
+            for (ibx = 0; ibx < nbetaX; ibx++)
+                hhx[ibx + nbetaX * ib] = hhx[ibx + nbetaX * iby];
+        }
+
+        iby = 0;
+        while (hhy[nbetaX / 2 * nbetaX + iby] < 0)
+            iby++;
+        for (ib = 0; ib < iby; ib++) {
+            for (ibx = 0; ibx < nbetaY; ibx++)
+                hhy[ib + nbetaX * ibx] = hhy[iby + nbetaX * ibx];
+        }
+        iby = nbetaX - 1;
+
+        while (hhy[nbetaX / 2 * nbetaX + iby] < 0)
+            iby--;
+        for (ib = iby + 1; ib < nbetaX; ib++) {
+            for (ibx = 0; ibx < nbetaY; ibx++)
+                hhy[ib + nbetaX * ibx] = hhy[iby + nbetaX * ibx];
+        }
+	*/
+	//#ifdef SAVE_ALL
+        debugSaveAll();
+	//#endif
+        delete[] hx;
+        delete[] hy;
+        delete[] hix;
+        delete[] hiy;
+
+        return;
+    }
+};
+
+class eta2InterpolationRosenblatt : public virtual eta2InterpolationBase,
+                               public virtual etaInterpolationRosenblatt {
+  public:
+    eta2InterpolationRosenblatt(int nx = 400, int ny = 400, int ns = 25,
+                           int nsy = 25, int nb = -1, int nby = -1,
+                           double emin = 1, double emax = 0)
+        : etaInterpolationBase(nx, ny, ns, nsy, nb, nby, emin, emax),
+          eta2InterpolationBase(nx, ny, ns, nsy, nb, nby, emin, emax),
+          etaInterpolationRosenblatt(nx, ny, ns, nsy, nb, nby, emin, emax){
+              //  std::cout << "e2pxy " << nb << " " << emin << " " << emax << std::endl;
+          };
+
+    eta2InterpolationRosenblatt(eta2InterpolationRosenblatt *orig)
+        : etaInterpolationBase(orig), etaInterpolationRosenblatt(orig){};
+
+    virtual eta2InterpolationRosenblatt *Clone() {
+        return new eta2InterpolationRosenblatt(this);
+    };
+};
+
+class eta3InterpolationRosenblatt : public virtual eta3InterpolationBase,
+                               public virtual etaInterpolationRosenblatt {
+  public:
+    eta3InterpolationRosenblatt(int nx = 400, int ny = 400, int ns = 25,
+                           int nsy = 25, int nb = -1, int nby = -1,
+                           double emin = 1, double emax = 0)
+        : etaInterpolationBase(nx, ny, ns, nsy, nb, nby, emin, emax),
+          eta3InterpolationBase(nx, ny, ns, nsy, nb, nby, emin, emax),
+          etaInterpolationRosenblatt(nx, ny, ns, nsy, nb, nby, emin, emax){
+              //   std::cout << "e3pxy " << nbeta << " " << etamin << " " << etamax
+              //   << " " << nSubPixels<< std::endl;
+          };
+
+    eta3InterpolationRosenblatt(eta3InterpolationRosenblatt *orig)
+        : etaInterpolationBase(orig), etaInterpolationRosenblatt(orig){};
+
+    virtual eta3InterpolationRosenblatt *Clone() {
+        return new eta3InterpolationRosenblatt(this);
+    };
+};
+
+#endif

slsDetectorCalibration/moenchExecutables/moench03Interpolation.cpp

@@ -17,7 +17,7 @@
 #endif
 
 //#include "etaInterpolationPosXY.h"
-#include "etaInterpolationPosXY.h"
+#include "etaInterpolationRosenblatt.h"
 #include "noInterpolation.h"
 //#include "etaInterpolationCleverAdaptiveBins.h"
 //#include "etaInterpolationRandomBins.h"
@@ -114,8 +114,8 @@ int main(int argc, char *argv[]) {
     //int f0 = -1;
     //  int nSubPixels = nsubpix;
 #ifndef NOINTERPOLATION
-    eta2InterpolationPosXY *interp =
-      new eta2InterpolationPosXY(NC, NR, nsubpix, nsubpix, etabins, etabins, etamin, etamax);
+    eta2InterpolationRosenblatt *interp =
+      new eta2InterpolationRosenblatt(NC, NR, nsubpix, nsubpix, etabins, etabins, etamin, etamax);
     // eta2InterpolationCleverAdaptiveBins *interp=new
     // eta2InterpolationCleverAdaptiveBins(NC, NR, nsubpix, etabins, etamin,
     // etamax);
@@ -134,6 +134,7 @@ int main(int argc, char *argv[]) {
     cout << "read ff " << argv[2] << endl;
     sprintf(fname, "%s", argv[2]);
     interp->readFlatField(fname);
+    
     interp->prepareInterpolation(ok); //, MAX_ITERATIONS);
 #endif
     // return 0;

slsDetectorCalibration/moenchExecutables/moenchRawDataProcess.cpp

@@ -174,6 +174,7 @@ int main(int argc, char *argv[]) {
     int nped = 1000;
     
     int cf = 0;
+
     int numberOfPackets=nrows/8;
 #ifdef RECT
     cout << "Should be rectangular but now it will crash! No data structure defined!" << endl;
@@ -181,7 +182,7 @@ int main(int argc, char *argv[]) {
     
 #ifndef MOENCH04
     moench03v2Data *decoder = new moench03v2Data(nrows/2);
-    cout << "MOENCH03!" << endl;
+    cout << "MOENCH03 V2!" << endl;
 #endif
 
 #ifdef MOENCH04
@@ -298,7 +299,7 @@ int main(int argc, char *argv[]) {
     double *ped=new double[nx * ny];//, *ped1;
     int pos,pos1;
     //return  0;
-  
+
     if (pedfile.find(".raw") != std::string::npos) {
       pos1=pedfile.rfind("/");
       strcpy(froot,pedfile.substr(pos1).c_str());
@@ -309,6 +310,7 @@ int main(int argc, char *argv[]) {
     cout << "PEDESTAL " << endl;
     if (pedfile.find(".tif") == std::string::npos) {
       sprintf(fname, "%s", pedfile.c_str());
+
             cout << fname << endl;
             std::time(&end_time);
             //cout << "aaa" << std::ctime(&end_time) << endl;
@@ -320,24 +322,28 @@ int main(int argc, char *argv[]) {
             if (filebin.is_open()) {
 	      ff = -1;
                 while (decoder->readNextFrame(filebin, ff, np, buff)) {
+
 		  if (np == numberOfPackets) {
 		    mt->pushData(buff);
                         mt->nextThread();
                         mt->popFree(buff);
                         ifr++;
 			if (ifr % 100 == 0)
+
 			  cout << ifr << " " << ff << " " << np << endl;
 			    //	break;
 		  } else {
-		    cout << ifr << " " << ff << " " << np << endl;
+		    cout << "bp "<< ifr << " " << ff << " " << np << endl;
 		    break;
 		  }
                     ff = -1;
                 }
                 filebin.close();
                 while (mt->isBusy()) {
-		          ;
-		          }
+
+		  ;
+
+		}
 		
 		sprintf(imgfname, "%s/%s_ped.tiff", outdir.c_str(),froot);
 		mt->writePedestal(imgfname);
@@ -382,7 +388,7 @@ int main(int argc, char *argv[]) {
        runmax--;
        flist.close();
        cout << "Found " << runmax << " files " << endl;
-       flist.open (fformat, std::ifstream::in);
+       flist.open (args["flist"].c_str(), std::ifstream::in);
     } 
 
     for (int irun = runmin; irun <= runmax; irun++) {
@@ -435,6 +441,7 @@ int main(int argc, char *argv[]) {
             ff = -1;
             ifr = 0;
             while (decoder->readNextFrame(filebin, ff, np, buff)) {
+
 	      if (np == numberOfPackets) {
                     //         //push
                     mt->pushData(buff);
@@ -462,8 +469,19 @@ int main(int argc, char *argv[]) {
                     }
 	      } else {
 		cout << "bp " << ifr << " " << ff << " " << np << endl;
+		
 	      }
                 ff = -1;
+
+
+		if (ifr % 100 == 0){
+		  sprintf(imgfname, "%s/%s_f%05d_ped.tiff", outdir.c_str(),froot, ifile++); 
+		  while (mt->isBusy()) {
+		    ;
+		  }
+		  mt->writePedestal(imgfname);
+		}
+
             }
             cout << "--" << endl;
             filebin.close();