Merge branch 'developer' of github.com:slsdetectorgroup/slsDetectorPackage into developer

.travis.yml

@@ -0,0 +1,45 @@
+sudo: false
+
+language: cpp
+
+matrix:
+  include:
+    - os: linux
+      env: CONDA_PY=3.6
+
+dist: trusty
+
+install:
+  - sudo apt-get update
+  - ldd --version
+  - wget https://repo.continuum.io/miniconda/Miniconda3-latest-Linux-x86_64.sh -O miniconda.sh;
+  - bash miniconda.sh -b -p $HOME/miniconda
+  - export PATH="$HOME/miniconda/bin:$PATH"
+  - rm -f miniconda.sh
+  - hash -r
+  - conda config --set always_yes yes --set changeps1 no
+  - conda config --add channels conda-forge
+  - conda config --add channels slsdetectorgroup
+  - conda update conda
+  - conda update --all
+  - conda install conda-build anaconda-client
+
+  # Useful for debugging any issues with conda
+  - conda info -a
+
+
+  # Replace dep1 dep2 ... with your dependencies
+  - conda create -q -n test-environment python=$CONDA_PY
+  - source activate test-environment
+  - conda-build .
+
+script:
+ - echo "No test scripts to be run!"
+
+
+deploy:
+  provider: script
+  script: find $HOME/miniconda/conda-bld/${TRAVIS_OS_NAME}-64 -name "*.tar.bz2" -exec anaconda -t $CONDA_TOKEN upload --force {} \;
+  on:
+    branch: developer
+

examples/moench03_T1_lab.config

@@ -1,4 +1,4 @@
-hostname bchip011
+hostname bchip085+
 
 patword 0000 0000000000000000
 patword 0001 0000000000000000
@@ -415,27 +415,45 @@ patwait2 0400
 patwaittime2 0
 
 
+####mcp2011
+#0:rx_udpip 10.1.1.102
+#0:detectorip 10.1.1.19
+#0:rx_udpport 32410
+####gui listening to 
+#zmqip 129.129.202.131
+#zmqport 30001
+####data streaming out of
+#rx_zmqip 10.1.2.103
+#rx_zmqport 30003
+#0:rx_hostname mpc2011
 
-0:rx_udpip 10.1.1.102
+####mx-test-1
+0:rx_udpip 10.1.1.100
 0:detectorip 10.1.1.19
 0:rx_udpport 32410
-#0:detectormac 00:ab:bc:cd:de:ef
-#0:rx_udpmac 70:10:6f:a0:b5:b1
-
-#gui listening to
-zmqip 129.129.202.131
+####gui listening to (on receiver pc)
+zmqip 129.129.202.92
 zmqport 30001
-
-#data streaming out of
-rx_zmqip 10.1.2.103
+####data streaming out of
+rx_zmqip 10.1.1.100
 rx_zmqport 30003
+0:rx_hostname pcmoench01
+
+
+
+
+
+
+
+
+
+
 
 #turn on datastream from commandline
 rx_datastream 1
 
 r_readfreq 1
 
-0:rx_hostname mpc2011
 #0:configuremac -1
 rx_datastream 1
 

manual/manual-client/Eiger_short.tex

@@ -64,7 +64,8 @@ The directory contains some executables that are needed to make your detector to
  \begin{verbatim}
 ./on #to switch modules on
 ./off #to switch modules off
-./hvget #gets the current HV value
+./state #tells you if is ON or OFF
+cat /var/log/pcu.log #displays the log if there are problem
 ./waterflow #returns the current waterflow returned by the flowmeter
 ./temp #returns the water temperature returned by the flowmeter
  \end{verbatim}
@@ -560,12 +561,23 @@ Here are the implemented options so far:
 \item {\tt{auto}} is the software controlled acquisition (does not use triggers), where {\tt{exptime}} and {\tt{period}} have to be set. Set number of cycles (i.e. triggers) to 1 using {\tt{cycles}}. Set number of frames using {\tt{frames}}.
 \item {\tt{trigger}} 1 frame taken for 1 trigger. Your {\tt{frames}} needs to be  1 always, {\tt{cycles}} can be changed and defines how many triggers are considered.  {\tt{exptime}} needs to be set. In the GUI this is called trigger exposure series.  
 \item {\tt{burst\_trigger}} gets only 1 trigger, but allows to take many frames. With {\tt{frames}} one can change the number of frames. {\tt{cycles}} needs to be 1.  {\tt{exptime}} and {\tt{period}} have to be set. In the gui it is called trigger readout.
-\item{\tt{gating}} allows to get a frame only when the trigger pulse is gating. Note that in this case the exp time and period only depend on the gating signal. {\tt{cycles}} allows to select how many gates to consider.  Set number of frames to 1 using {\tt{frames}}.
+\item{\tt{gating}} allows to get a frame only when the trigger pulse is gating. Note that in this case the exp time and period only depend on the gating signal. {\tt{cycles}} allows to select how many gates to consider.  Set number of frames to 1 using {\tt{frames}}. ATTENTION: if you are in 16 bit mode and you are applying online rate corrections, as now the exptime is generated by the trigger, you might not have correct rate corrections. If you know what the exposure time is in the gating signal, then you can set the {\tt{exptime}} once and the rate corrections will be correct. If the exposure time is unknow, it is recommended that you switch off the rate corrections. In 32 bit mode, it does not matter as the rate corrections depends on the {\tt{subexptime}} which is software set independently from the gate exptime.
+
 \end{itemize}
 
 Hardware-wise, the ENABLE OUT signal outputs when the chips are really acquiring. This means that the single subframes will be output in 32 bit mode. The TRIGGER OUT outputs the sum-up-signal at the moment (which is useless). This will be changed in the future to output the envelop of the enable signal. 
 
-We are planning to change some functionality, i.e. unify the {\tt{trigger}} and {\tt{burst}} trigger modes and make both {\tt{frames}} and {\tt{cycles}} configurable at the same time.
+We are planning to change some functionality, i.e. unify the {\tt{trigger}} and {\tt{burst\_trigger}} trigger modes and make both {\tt{frames}} and {\tt{cycles}} configurable at the same time.
+
+There is the possibility to use {\tt{timing trigger/burst\_trigger}} and send software single commands to fake the trigger. This is done with:
+\begin{verbatim}
+sls_detector_put 0-timing [trigger/burst_trigger]
+sls_detector_put 0-frames x
+sls_detector_put 0-cycles y
+sls_detector_status trigger
+\end{verbatim}
+ Note that this functionality is very (!) useful if you need to do something between and acquisition and the next. This can be used to do a fast threshold scan for example. See section~\ref{Sec:fastthresholdscan}.
+
 
 \section{Autosumming and rate corrections} \label{advanced}
 
@@ -1046,6 +1058,12 @@ To load the special noise file look at {\tt{settingsdir/eiger/standard/eigernois
 \begin{verbatim}
 sls_detector_put trimbits ../settingsdir/eiger/standard/eigernoise
 \end{verbatim}
+To exit from this pattern noise, just set the theshold to something known. 
+\begin{verbatim}
+\item sls_detector_put threshold 50000 standard
+\end{verbatim}
+where 5000 would be a value in eV and {/tt{standard}} is important in this case.
+
 
 \section{Troubleshooting}
 \subsection{Cannot successfully finish an acquisition}
@@ -1167,6 +1185,13 @@ If you see strange lines in vertical occurring at period patterns, it is a memor
 \subsection{ssh to the boards takes long}
 Depending on your network setup,  to speed up the ssh to the boards from a pc with internal dhcp server running: \textbf{iptables -t nat -A POSTROUTING -o eth1 -j MASQUERADE; echo "1" > /proc/sys/net/ipv4/ip\_forward}, where eth1 has to be the 1Gb network device on the pc
 
+\subsection{Generate keys on the boards not to have to type the password}
+\begin{verbatim}
+export AFSDIRS64=/afs/psi.ch/intranet/Controls/Software/Trolltech/SL6-x86_64
+ssh-copy-id -i /afs/psi.ch/user/t/tinti_g/.ssh/id_rsa.pub  root@beb100
+ssh-keygen
+\end{verbatim}
+
 \subsection{Check firmware version installed on BEB}
 You can either ask in the client as described in section~\ref{api}, or login to the boards directly. Follow some steps described in Section~\ref{server}. 
 \begin{verbatim}
@@ -1190,6 +1215,13 @@ Scroll up in the terminal till you find:\\
 *************** MASTER/SLAVE ***************\\
 *************** NORMAL/SPECIAL ***************\\
 
+There is also an easier way, that is that only the master module will reaturn the real value of the HV. If you have more than 1 detector system, then you will have more than 1 physical master, as the HV needs to be applied to all the systems.
+
+\begin{verbatim}
+for i in $(seq 0 36); do sls_detector_put $i:vhighvoltage; done
+\end{verbatim}
+ Only the master will return to you a sensible number (150 normally). the others will return -999.
+
 \subsection{'Cannot connect to socket'}
 This error is typically due to the detector server not running. For why, see section~\ref{servernot}.
 

recipe/build.sh

@@ -0,0 +1,15 @@
+mkdir build
+mkdir install
+cd build
+cmake .. \
+      -DCMAKE_PREFIX_PATH=$CONDA_PREFIX \
+      -DCMAKE_INSTALL_PREFIX=install \
+      -DUSE_TEXTCLIENT=ON \
+      -DUSE_RECEIVER=ON \
+      -DUSE_GUI=ON \
+      -DCMAKE_BUILD_TYPE=Release \
+      -DUSE_HDF5=OFF\
+     
+
+cmake --build . -- -j10
+cmake --build . --target install

recipe/copy_gui.sh

@@ -0,0 +1,15 @@
+
+mkdir $PREFIX/lib
+mkdir $PREFIX/bin
+mkdir $PREFIX/include
+
+#No libs for gui?
+
+#Binaries
+cp build/bin/gui_client $PREFIX/bin/.
+cp build/bin/slsDetectorGui $PREFIX/bin/.
+
+
+#Which headers do we need for development??
+
+# cp include/some_lib.h $PREFIX/include/.

recipe/copy_lib.sh

@@ -0,0 +1,23 @@
+
+mkdir $PREFIX/lib
+mkdir $PREFIX/bin
+mkdir $PREFIX/include
+mkdir $PREFIX/include/slsDetectorPackage
+
+#Shared and static libraries
+cp build/bin/libSlsDetector.so $PREFIX/lib/.
+cp build/bin/libSlsDetector.a $PREFIX/lib/.
+cp build/bin/libSlsReceiver.so $PREFIX/lib/.
+cp build/bin/libSlsReceiver.a $PREFIX/lib/.
+
+#Binaries
+cp build/bin/sls_detector_acquire $PREFIX/bin/.
+cp build/bin/sls_detector_get $PREFIX/bin/.
+cp build/bin/sls_detector_put $PREFIX/bin/.
+cp build/bin/sls_detector_help $PREFIX/bin/.
+cp build/bin/slsReceiver $PREFIX/bin/.
+cp build/bin/slsMultiReceiver $PREFIX/bin/.
+
+#Which headers do we need for development??
+cp build/install/include/* $PREFIX/include/slsDetectorPackage/
+# cp include/some_lib.h $PREFIX/include/.

recipe/meta.yaml

@@ -0,0 +1,89 @@
+
+package:
+  name: sls_detector_software
+  version: "developer"
+
+source:
+  - path: ..
+
+build:
+  number: 0
+  rpaths:
+    - lib/
+
+requirements:
+  build:
+    - {{ compiler('c') }}
+    - {{compiler('cxx')}}
+    - cmake
+    - qwt 6.*
+    - qt=4.8.7=7
+    - zeromq=4.2.5=hfc679d8_5
+    - pyzmq
+    - xorg-libx11
+    - xorg-libice
+    - xorg-libxext
+    - xorg-libsm
+    - xorg-libxau
+    - xorg-libxrender
+    - xorg-libxfixes
+    - {{ cdt('mesa-libgl-devel') }}      # [linux]
+    - {{ cdt('mesa-libegl-devel') }}     # [linux]
+    - {{ cdt('mesa-dri-drivers') }}      # [linux]
+    - {{ cdt('libselinux') }}      # [linux]
+    - {{ cdt('libxdamage') }}            # [linux]
+    - {{ cdt('libxxf86vm') }}            # [linux]
+
+  host:
+    - libstdcxx-ng
+    - libgcc-ng
+    - libpng >=1.6.32,<1.6.35
+    - xorg-libx11
+    - xorg-libice
+    - xorg-libxext
+    - xorg-libsm
+    - xorg-libxau
+    - xorg-libxrender
+    - xorg-libxfixes
+
+  run:
+    - libstdcxx-ng
+    - libgcc-ng 
+
+outputs:
+  - name: sls_detector_lib
+    version: "developer"
+    script: copy_lib.sh
+    requirements:
+      build:
+        - {{ compiler('c') }}
+        - {{compiler('cxx')}}
+  - name: sls_detector_gui
+    version: "developer"
+    script: copy_gui.sh
+    requirements:
+      build:
+      - {{ compiler('c') }}
+      - {{compiler('cxx')}}
+      - cmake
+      - qwt 6.*
+      - qt=4.8.7=7
+      - zeromq=4.2.5=hfc679d8_5
+      - pyzmq
+      - xorg-libx11
+      - xorg-libice
+      - xorg-libxext
+      - xorg-libsm
+      - xorg-libxau
+      - xorg-libxrender
+      - xorg-libxfixes
+      - {{ cdt('mesa-libgl-devel') }}      # [linux]
+      - {{ cdt('mesa-libegl-devel') }}     # [linux]
+      - {{ cdt('mesa-dri-drivers') }}      # [linux]
+      - {{ cdt('libselinux') }}      # [linux]
+      - {{ cdt('libxdamage') }}            # [linux]
+      - {{ cdt('libxxf86vm') }}            # [linux]
+      run:
+        - sls_detector_lib=developer
+        - qwt 6.*
+        - qt=4.8.7=7

slsDetectorCalibration/Stat.h

@@ -0,0 +1,45 @@
+class Stat
+  {
+  public:
+    
+  Stat() : n(0), m(0.), m2(0.) {}
+    
+    void Clear()
+    {
+      n = 0;
+      m=0;
+      m2=0;
+    }
+    
+        void Push(double x)
+        {
+	  
+	  m+=x;
+	  m2+=x*x;
+            n++;
+        }
+
+        int NumDataValues() const
+        {
+            return n;
+        }
+
+        double Mean() const
+        {
+            return (n > 0) ? m/n : 0.0;
+        }
+
+        double Variance() const
+        {
+	  return ( (n >0 ) ? (m2/n-m*m/(n*n)) : 0.0 );
+        }
+
+        double StandardDeviation() const
+        {
+            return sqrt( Variance() );
+        }
+
+    private:
+        int n;
+        double m, m2;
+    };

slsDetectorCalibration/analogDetector.h

@@ -221,8 +221,8 @@ template <class dataType> class analogDetector {
     if (gm) {
       if (gmap) delete [] gmap; 
       gmap=new double[nnx*nny];
-      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];
 	}
       }
@@ -240,8 +240,8 @@ template <class dataType> class analogDetector {
     void *ret;
     if (gmap)  {
       gm=new float[nx*ny];
-      for (int ix=0; ix<nx; ix++) {
 	for (int iy=0; iy<ny; iy++) {
+      for (int ix=0; ix<nx; ix++) {
 	  gm[iy*nx+ix]=gmap[iy*nx+ix];
 	}
       }
@@ -325,8 +325,8 @@ template <class dataType> class analogDetector {
 
     virtual void addToCommonMode(char *data){ 
       if (cmSub) {
-	for (int ix=xmin; ix<xmax; ix++) {
 	  for (int iy=ymin; iy<ymax; iy++) { 
+	for (int ix=xmin; ix<xmax; ix++) {
 	    // if (getNumpedestals(ix,iy)>0) 
 	    if (det->isGood(ix,iy))
 	      addToCommonMode(data, ix, iy);
@@ -387,8 +387,8 @@ template <class dataType> class analogDetector {
     virtual double* getPedestal(double *ped){
       if (ped==NULL)
 	ped=new double[nx*ny];
-      for (int ix=0; ix<nx; ix++) {
 	for (int iy=0; iy<ny; iy++) {
+      for (int ix=0; ix<nx; ix++) {
 	  ped[iy*nx+ix]=stat[iy][ix].getPedestal();
 	  //cout << ped[iy*nx+ix] << " " ;
 	}
@@ -405,8 +405,8 @@ template <class dataType> class analogDetector {
     virtual double* getPedestalRMS(double *ped=NULL){
       if (ped==NULL)
 	ped=new double[nx*ny];
-      for (int ix=0; ix<nx; ix++) {
 	for (int iy=0; iy<ny; iy++) {
+      for (int ix=0; ix<nx; ix++) {
 	  ped[iy*nx+ix]=stat[iy][ix].getPedestalRMS();
 	}
       }
@@ -445,8 +445,8 @@ template <class dataType> class analogDetector {
     */
     virtual void setPedestal(double *ped, double *rms=NULL, int m=-1){
       double rr=0;
-      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);
 	};
@@ -474,8 +474,8 @@ template <class dataType> class analogDetector {
        \param rms pointer to array of pedestal rms
     */
  virtual void setPedestalRMS(double *rms){
-      for (int ix=xmin; ix<xmax; ix++) {
 	for (int iy=ymin; iy<ymax; iy++) {
+      for (int ix=xmin; ix<xmax; ix++) {
 	  stat[iy][ix].setPedestalRMS(rms[iy*nx+ix]);
 	};
       };
@@ -498,8 +498,8 @@ template <class dataType> class analogDetector {
 
 #endif
       gm=new float[nx*ny];
-      for (int ix=0; ix<nx; ix++) {
 	for (int iy=0; iy<ny; iy++) {
+      for (int ix=0; ix<nx; ix++) {
 	    gm[iy*nx+ix]=image[iy*nx+ix];
 #ifdef ROOTSPECTRUM
 	    hmap->SetBinContent(ix+1, iy+1,image[iy*nx+ix]);
@@ -546,8 +546,8 @@ template <class dataType> class analogDetector {
       TH2F *hmap=new TH2F("hmap","hmap",nx, -0.5,nx-0.5, ny, -0.5, ny-0.5);
 
 #endif
-    for (int ix=0; ix<nx; ix++) {
       for (int iy=0; iy<ny; iy++) {
+    for (int ix=0; ix<nx; ix++) {
 	/* if (cmSub)  */
 	/*     gm[iy*nx+ix]=stat[iy][ix].getPedestal()-cmSub->getCommonMode(); */
 	/* else */
@@ -593,8 +593,8 @@ template <class dataType> class analogDetector {
 
 
     if (gm) {
-      for (int ix=0; ix<nnx; ix++) {
 	for (int iy=0; iy<nny; iy++) {
+      for (int ix=0; ix<nnx; ix++) {
 	  stat[iy][ix].setPedestal(gm[iy*nx+ix],-1,-1);
 	}
       }
@@ -618,8 +618,8 @@ template <class dataType> class analogDetector {
 
 
     if (gm) {
-      for (int ix=0; ix<nnx; ix++) {
 	for (int iy=0; iy<nny; iy++) {
+      for (int ix=0; ix<nnx; ix++) {
 	  image[iy*nx+ix]=gm[iy*nx+ix];
 	}
       }
@@ -644,8 +644,8 @@ template <class dataType> class analogDetector {
     float *gm=NULL;
     void *ret;
       gm=new float[nx*ny];
-      for (int ix=0; ix<nx; ix++) {
 	for (int iy=0; iy<ny; iy++) {
+      for (int ix=0; ix<nx; ix++) {
 	    gm[iy*nx+ix]=stat[iy][ix].getPedestalRMS();
 	}
       }
@@ -666,8 +666,8 @@ template <class dataType> class analogDetector {
     if (nnx>nx) nnx=nx;
     if (nny>ny) nny=ny;
     if (gm) {
-      for (int ix=0; ix<nnx; ix++) {
 	for (int iy=0; iy<nny; iy++) {
+      for (int ix=0; ix<nnx; ix++) {
 	  stat[iy][ix].setPedestalRMS(gm[iy*nx+ix]);
 	}
       }
@@ -699,8 +699,8 @@ template <class dataType> class analogDetector {
             
       //cout << xmin << " " << xmax << endl;
       // cout << ymin << " " << ymax << endl;
-      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);
 	  //if (ix==10 && iy==10) 
@@ -823,8 +823,8 @@ template <class dataType> class analogDetector {
       if (val==NULL)
 	val=image;//new double[nx*ny];
       
-      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))
 	      val[iy*nx+ix]+=subtractPedestal(data, ix, iy,cm);
 	}
@@ -951,8 +951,8 @@ template <class dataType> class analogDetector {
 
        addToCommonMode(data);
 
-       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);
 	 }
@@ -966,8 +966,8 @@ template <class dataType> class analogDetector {
      
     */
    virtual void clearImage(){  
-     for (int ix=0; ix<nx; ix++) {
        for (int iy=0; iy<ny; iy++) { 
+     for (int ix=0; ix<nx; ix++) {
 	 image[iy*nx+ix]=0;
        }
      }
@@ -997,8 +997,8 @@ template <class dataType> class analogDetector {
     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++) 
+	for (ix=0; ix<nx; ix++) 
 	    stat[iy][ix].SetNPedestals(i); 
       return stat[0][0].SetNPedestals();
     };
@@ -1030,8 +1030,8 @@ template <class dataType> class analogDetector {
       if (ymi<0) ymi=ymin;
       if (yma<0) yma=ymax;
       
-      for (int ix=xmi; ix<xma; ix++)
 	for (int iy=ymi; iy<yma; iy++)
+      for (int ix=xmi; ix<xma; ix++)
 	  if (det->isGood(ix,iy)) {
 	    if (ix>=0 && ix<nx && iy>=0 && iy<ny) 
 	      val+=getNPhotons(data, ix, iy);

slsDetectorCalibration/commonModeSubtractionNew.h

@@ -0,0 +1,116 @@
+#ifndef COMMONMODESUBTRACTION_H
+#define COMMONMODESUBTRACTION_H
+
+#include <cmath>
+
+class commonModeSubtraction {
+
+  /** @short class to calculate the common mode of the pedestals based on an approximated moving average*/ 
+
+ public:
+  
+  /** constructor
+      \param nn number of samples for the moving average to calculate the average common mode
+      \param iroi number of regions on which one can calculate the common mode separately. Defaults to 1 i.e. whole detector
+
+  */
+ commonModeSubtraction(int iroi=1, int ns=3) : nROI(iroi), nsigma(ns)  {
+    mean=new double[nROI];  
+    mean2=new double[nROI]; 
+    nCm=new double[nROI];
+  };
+
+    /** destructor - deletes the moving average(s) and the sum of pedestals calculator(s) */
+    virtual ~commonModeSubtraction() {delete [] mean; delete [] mean2; delete [] nCm;};
+    
+
+    /** clears the moving average and the sum of pedestals calculation - virtual func*/
+    virtual void  Clear(){
+      for (int i=0; i<nROI; i++) {
+	mean[i]=0;  
+	nCm[i]=0; 
+	mean2[i]=0;
+      }};
+    
+    /** 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(){
+      for (int i=0; i<nROI; i++) {
+	//	if (nCm[i]>0) cmStat[i].Calc(cmPed[i]/nCm[i]); 
+	nCm[i]=0; 
+	mean[i]=0;
+	mean2[i]=0;
+      }};
+    
+    /** adds the pixel to the sum of pedestals -- virtual func must be overloaded to define the regions of interest
+	\param val value to add
+	\param ix pixel x coordinate
+	\param iy pixel y coordinate
+    */
+    virtual void addToCommonMode(double val, int ix=0, int iy=0) {
+      
+      int iroi=getROI(ix,iy);
+      // if (iroi==0) val=100;
+      // else val=-100;
+      // if (isc>=0 && isc<nROI) {
+      if (iroi>=0 && iroi<nROI) {
+	mean[iroi]+=val; 
+	mean2[iroi]+=val*val; 
+	nCm[iroi]++;
+      }
+    };
+
+    /** gets the common mode i.e. the difference between the current average sum of pedestals mode and the average pedestal 
+	\param ix pixel x coordinate
+	\param iy pixel y coordinate
+	\return the difference  between the current average sum of pedestals and the average pedestal
+    */
+    virtual double getCommonMode(int ix=0, int iy=0) {
+      int iroi=getROI(ix,iy); 
+      /* if (iroi==0) */
+      /* 	return 100; */
+      /* else */
+      /* 	return -100; */
+
+      if (iroi>=0 && iroi<nROI) {
+	if (nCm[iroi]>0) 
+	  return mean[iroi]/nCm[iroi]; 
+      }
+      return 0;
+    };
+    
+    /** gets the common mode i.e. the difference between the current average sum of pedestals mode and the average pedestal 
+	\param ix pixel x coordinate
+	\param iy pixel y coordinate
+	\return the difference  between the current average sum of pedestals and the average pedestal
+    */
+    virtual double getCommonModeRMS(int ix=0, int iy=0) {
+      int iroi=getROI(ix,iy);
+      if (iroi>=0 && iroi<nROI) {
+	if (nCm[iroi]>0) 
+	  return sqrt(mean2[iroi]/nCm[iroi]-(mean[iroi]/nCm[iroi])*(mean[iroi]/nCm[iroi])); 
+      }
+      return 0;
+    };
+
+    /** 
+	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;};
+
+
+
+ protected:
+  double *mean; /**<array of moving average of the pedestal average per region of interest */
+  double *mean2; /**< array storing the sum of pedestals per region of interest */
+  double *nCm; /**< array storing the number of pixels currently contributing to the pedestals */
+  int nsigma; /** number of rms above which the pedestal should be considered as a photon */
+  const int nROI; /**< constant parameter for number of regions on which the common mode should be calculated separately e.g. supercolumns */
+
+
+};
+
+
+
+
+
+#endif

slsDetectorCalibration/dataStructures/Mythen3_01_jctbData.h

@@ -74,7 +74,7 @@ class mythen3_01_jctbData : public slsDetectorData<short unsigned int> {
 	  for (ib=0; ib<nb; ib++) {
 	    if (word&(1<<bit[ib]))  {
 	      cout << "+" ;
-	      val[iw+nch*(ib/nb)]|=(1<<idr);
+	      val[iw+nch/nb*(ib)]|=(1<<idr);
 	    }  else { 
 	      	cout << "-" ; 
 	       } 
@@ -96,8 +96,8 @@ class mythen3_01_jctbData : public slsDetectorData<short unsigned int> {
       ii++;
     }//end for
     
-    cout << "Decoded "<<ii << " samples"<< endl;
-    cout << "Should be "<< nch/nb*dr+off << " samples"<< endl;
+    cout << "M3.01 Decoded "<<ii << " samples"<< endl;
+    cout << "M3.01 Should be "<< nch/nb*dr+off << " samples"<< endl;
   
     return val;
  }

slsDetectorCalibration/dataStructures/Mythen3_02_jctbData.h

@@ -57,8 +57,9 @@ class mythen3_02_jctbData : public mythen3_01_jctbData {
     int ioff=0;
     int idr=0;
     int ib=0;
-    int iw=0;
+    int ich=0;
     int ii=0;
+    int iw=0;
     bit[0]=17;//19;
     bit[1]=6;//8;
     idr=0;
@@ -68,7 +69,8 @@ class mythen3_02_jctbData : public mythen3_01_jctbData {
     wp=(int64_t*)ptr;
     
     for (iw=0; iw<nch/nb; iw) {
-      word=*wp;;
+      word=*wp;
+
       if (ioff<off) {
 	ioff++;
 	 cout <<"*";
@@ -78,10 +80,11 @@ class mythen3_02_jctbData : public mythen3_01_jctbData {
 	  for (ib=0; ib<nb; ib++) {
 	    if (word&(1<<bit[ib]))  {
 	      cout << "+" ;
-	      val[iw+nch*(ib/nb)]|=(1<<idr);
+	      val[iw+nch/nb*(ib)]|=(1<<idr);
 	    }  else { 
 	      	cout << "-" ; 
-	       } 
+	    } 
+	    //    cout <<  iw+nch/nb*(ib)<< " " ;
 	  }//end for()
 	}
 	
@@ -94,14 +97,14 @@ class mythen3_02_jctbData : public mythen3_01_jctbData {
 	  cout <<dec << iw<<endl;
 	  iw++;
 	}//end if()
-
+	
       }//end else()
       wp+=1;
       ii++;
     }//end for
     
-    cout << "Decoded "<<ii << " samples"<< endl;
-    cout << "Should be "<< nch/nb*dr+off << " samples"<< endl;
+    cout << "M3.02 Decoded "<<ii << " samples"<< endl;
+    cout << "M3.02 Should be "<< nch/nb*dr+off << " samples"<< endl;
   
     return val;
  }

slsDetectorCalibration/dataStructures/moench02CtbData.h

@@ -35,7 +35,7 @@ class moench02CtbData : public slsDetectorData<uint16_t> {
 
     
     int adc_off[4]={40,0,120,80};
-    int adc_nr[4]={8,10,20,23};
+    int adc_nr[4]={8,10,20,22};
     int row, col;
 
     int isample;

slsDetectorCalibration/dataStructures/moench03T1ReceiverDataNew.h

@@ -2,6 +2,7 @@
 #define  MOENCH03T1RECDATANEW_H
 #include "slsDetectorData.h"
 
+#define VERSION_V2
     /**
         @short  structure for a Detector Packet or Image Header
         @li frameNumber is the frame number
@@ -32,6 +33,10 @@
         uint16_t roundRNumber;    /**< is the round robin set number */
         uint8_t detType;        /**< is the detector type see :: detectorType */
         uint8_t version;        /**< is the version number of this structure format */
+#ifdef VERSION_V2
+      uint64_t packetCaught[8];        /**< is the version number of this structure format */
+#endif
+
     } sls_detector_header;
 
 

slsDetectorCalibration/dataStructures/moench03T1ZmqDataNew.h

@@ -3,16 +3,16 @@
 #include "slsDetectorData.h"
 
 
-
 class moench03T1ZmqDataNew : public slsDetectorData<uint16_t> {
  
  private:
   
-  int iframe;
+  // int iframe;
   int nadc;
   int sc_width;
   int sc_height;
   const int nSamples;
+  const int offset;
 
  public:
 
@@ -25,7 +25,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), nSamples(ns) {
+ moench03T1ZmqDataNew(int ns=5000): slsDetectorData<uint16_t>(400, 400, ns*32*2+sizeof(int)), nSamples(ns), offset(sizeof(int)) {
 
     int nadc=32;
     int sc_width=25;
@@ -37,7 +37,7 @@ class moench03T1ZmqDataNew : public slsDetectorData<uint16_t> {
     		    0,25,50,75,0,25,50,75};
 
     int row, col;
-
+    
     int isample;
     int iadc;
     int ix, iy;
@@ -60,50 +60,48 @@ class moench03T1ZmqDataNew : public slsDetectorData<uint16_t> {
 	    } else {
 	      row=200+i/sc_width;
 	    }
-	    dataMap[row][col]=(nadc*i+iadc)*2;//+16*(ip+1);
-	    if (dataMap[row][col]<0 || dataMap[row][col]>=nSamples*2*32)
+	    dataMap[row][col]=(nadc*i+iadc)*2+offset;//+16*(ip+1);
+	    if (dataMap[row][col]<0 || dataMap[row][col]>=dataSize)
 	      cout << "Error: pointer " << dataMap[row][col] << " out of range "<< endl;
 	  }
 	}
       }
     }
 
-    int ipacket;
-    int ibyte;
+    
     int ii=0;
-    for (int ipacket=0; ipacket<npackets; ipacket++) {
-      for (int ibyte=0;  ibyte< 8192/2; ibyte++) {
-	i=ipacket*8208/2+ibyte;
-	/* if (ibyte<8) { */
-	/* //header! */
-	/*   xmap[i]=-1; */
-	/*   ymap[i]=-1; */
-	/* } else { */
-	  // ii=ibyte+128*32*ipacket;
-	  isample=ii/nadc;
-	  if (isample<nSamples) {
-	  iadc=ii%nadc;
-	  adc4 = (int)iadc/4;
-	  ix=isample%sc_width;
-	  iy=isample/sc_width;
-	  if (adc4%2==0) {
-	    xmap[i]=adc_nr[iadc]+ix;
-	    ymap[i]=ny/2-1-iy;
-	  } else {
-	    xmap[i]=adc_nr[iadc]+ix;
-	    ymap[i]=ny/2+iy;
-	  }
-	  }
+    
+    for (i=0; i<  dataSize; i++) {
+      if (i<offset) {
+    	//header! */
+    	xmap[i]=-1;
+    	ymap[i]=-1;
+      } else {
+    	  // ii=ibyte+128*32*ipacket;
+    	isample=ii/nadc;
+    	if (isample<nSamples) {
+    	  iadc=ii%nadc;
+    	  adc4 = (int)iadc/4;
+    	  ix=isample%sc_width;
+    	  iy=isample/sc_width;
+    	  if (adc4%2==0) {
+    	    xmap[i]=adc_nr[iadc]+ix;
+    	    ymap[i]=ny/2-1-iy;
+    	  } else {
+    	    xmap[i]=adc_nr[iadc]+ix;
+    	    ymap[i]=ny/2+iy;
+    	  }
+    	}
 	  
-	ii++;
-	//	}
+    	ii++;
       }
+  
     }
     
 
     
     
-    iframe=0;
+    // iframe=0;
     //  cout << "data struct created" << endl;
   };
     
@@ -127,7 +125,7 @@ class moench03T1ZmqDataNew : public slsDetectorData<uint16_t> {
 
 
 
-  int getFrameNumber(char *buff){return iframe;};//*((int*)(buff+5))&0xffffff;};   
+  int getFrameNumber(char *buff){return *((int*)buff);};//*((int*)(buff+5))&0xffffff;};   
 
   /**
 
@@ -215,8 +213,8 @@ class moench03T1ZmqDataNew : public slsDetectorData<uint16_t> {
 
 	  if (filebin.is_open()) {
 	    if (filebin.read(data, 32*2*nSamples) ){
-	      iframe++;
-	      ff=iframe;
+	      // iframe++;
+	      //ff=iframe;
 	      return data;
 	    }
 	  }
@@ -249,6 +247,8 @@ class moench03T1ZmqDataNew : public slsDetectorData<uint16_t> {
 
 
 
+  // virtual int setFrameNumber(int ff){iframe=ff};
+
 
 
 

slsDetectorCalibration/interpolatingDetector.h

@@ -178,7 +178,7 @@ int addFrame(char *data,  int *ph=NULL, int ff=0) {
   double int_x, int_y;
   double eta_x, eta_y;
   if (interp) {
-  cout << "int" << endl;
+    // cout << "int" << endl;
   pthread_mutex_lock(fi); 
     for (nph=0; nph<nphFrame; nph++) {
       if (ff) {
@@ -239,7 +239,7 @@ int addFrame(char *data,  int *ph=NULL, int ff=0) {
       if (interp)
 	interp->prepareInterpolation(ok);
 	pthread_mutex_unlock(fi);  */
-      cout << "det" << endl;
+      // cout << "det" << endl;
       return interp;
     };
   

slsDetectorCalibration/interpolations/eta2InterpolationBase.h

@@ -14,47 +14,18 @@ class eta2InterpolationBase : public virtual etaInterpolationBase {
   
  public:
  eta2InterpolationBase(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) {
-    // cout << "e2ib " << nb << " " << emin << " " << emax << endl; 
-
-    /* if (nbeta<=0) { */
-    /*   nbeta=nSubPixels*10;  */
+    
+    /* if (etamin>=etamax) { */
+    /*   etamin=-1; */
+    /*   etamax=2;  */
+    /*   //  cout << ":" <<endl; */
     /* } */
-    if (etamin>=etamax) {
-      etamin=-1;
-      etamax=2; 
-      //  cout << ":" <<endl;
-    }
-    etastep=(etamax-etamin)/nbeta;
-#ifdef MYROOT1
-    delete heta;
-    delete hhx;
-    delete hhy;
-    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
-    /* delete [] heta; */
-    /* delete [] hhx; */
-    /* delete [] hhy; */
-    /* heta=new int[nbeta*nbeta]; */
-    /* hhx=new float[nbeta*nbeta]; */
-    /* hhy=new float[nbeta*nbeta]; */
+    /* etastep=(etamax-etamin)/nbeta; */
     
-#endif
-    
-    //   cout << nbeta << " " << etamin << " " << etamax << endl;
   };
   
  eta2InterpolationBase(eta2InterpolationBase *orig): etaInterpolationBase(orig){ };
 
-/*  virtual eta2InterpolationBase* Clone()=0; {
-
-    return new eta2InterpolationBase(this);
-
-  };
-
-					   */
  
   //////////////////////////////////////////////////////////////////////////////
   //////////// /*It return position hit for the event in input */ //////////////
@@ -100,12 +71,7 @@ class eta2InterpolationBase : public virtual etaInterpolationBase {
   virtual void getInterpolatedPosition(int x, int y, double totquad,int quad,double *cl,double &int_x, double &int_y) {
     
      double cc[2][2];
-     double *cluster[3];
      int xoff, yoff;
-     cluster[0]=cl;
-     cluster[1]=cl+3;
-     cluster[2]=cl+6;
-     
      switch (quad) {
      case BOTTOM_LEFT:
        xoff=0;
@@ -128,10 +94,10 @@ class eta2InterpolationBase : public virtual etaInterpolationBase {
      } 
      double etax, etay;
      if (nSubPixels>2) { 
-       cc[0][0]=cluster[yoff][xoff];
-       cc[1][0]=cluster[yoff+1][xoff];
-       cc[0][1]=cluster[yoff][xoff+1];
-       cc[1][1]=cluster[yoff+1][xoff+1];
+       cc[0][0]=cl[xoff+3*yoff];
+       cc[1][0]=cl[xoff+3*(yoff+1)];
+       cc[0][1]=cl[xoff+1+3*yoff];
+       cc[1][1]=cl[xoff+1+3*(yoff+1)];
        calcEta(totquad,cc,etax,etay);
      }
      return getInterpolatedPosition(x,y,etax, etay,quad,int_x,int_y);
@@ -143,11 +109,7 @@ class eta2InterpolationBase : public virtual etaInterpolationBase {
   virtual void getInterpolatedPosition(int x, int y, double totquad,int quad,int *cl,double &int_x, double &int_y) {
     
      double cc[2][2];
-     int *cluster[3];
      int xoff, yoff;
-     cluster[0]=cl;
-     cluster[1]=cl+3;
-     cluster[2]=cl+6;
      
      switch (quad) {
      case BOTTOM_LEFT:
@@ -171,10 +133,10 @@ class eta2InterpolationBase : public virtual etaInterpolationBase {
      } 
      double etax, etay;
      if (nSubPixels>2) { 
-       cc[0][0]=cluster[yoff][xoff];
-       cc[1][0]=cluster[yoff+1][xoff];
-       cc[0][1]=cluster[yoff][xoff+1];
-       cc[1][1]=cluster[yoff+1][xoff+1];
+       cc[0][0]=cl[xoff+3*yoff];
+       cc[1][0]=cl[xoff+3*(yoff+1)];
+       cc[0][1]=cl[xoff+1+3*yoff];
+       cc[1][1]=cl[xoff+1+3*(xoff+1)];
        calcEta(totquad,cc,etax,etay);
      }
      return getInterpolatedPosition(x,y,etax, etay,quad,int_x,int_y);
@@ -221,60 +183,46 @@ class eta2InterpolationBase : public virtual etaInterpolationBase {
 
      if (nSubPixels>2) { 
 
-#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) {
-      cout << "x*"<< ex << endl;
-      ex=0;
-	} 
-    if (ex>=nbeta) {
-      cout << "x?"<< ex << endl;
-      ex=nbeta-1;
-      
-    }
-    if (ey<0) {
-      cout << "y*"<< ey << endl;
-      ey=0;
-	} 
-    if (ey>=nbeta) {
-      cout << "y?"<< ey << endl;
-      ey=nbeta-1;
-      
-    }
+       ex=(etax-etamin)/etastep;
+       ey=(etay-etamin)/etastep;
+       if (ex<0) {
+	 cout << "x*"<< ex << endl;
+	 ex=0;
+       } 
+       if (ex>=nbeta) {
+	 cout << "x?"<< ex << endl;
+	 ex=nbeta-1;
+       }
+       if (ey<0) {
+	 cout << "y*"<< ey << endl;
+	 ey=0;
+       } 
+       if (ey>=nbeta) {
+	 cout << "y?"<< ey << endl;
+	 ey=nbeta-1;
+       }
     
     
    
-    xpos_eta=(((double)hhx[(ey*nbeta+ex)]))+dX ;///((double)nSubPixels);
-    ypos_eta=(((double)hhy[(ey*nbeta+ex)]))+dY ;///((double)nSubPixels);
-    //else
-    //return 0;
-    
-#endif
+       xpos_eta=(((double)hhx[(ey*nbeta+ex)]))+dX ;///((double)nSubPixels);
+       ypos_eta=(((double)hhy[(ey*nbeta+ex)]))+dY ;///((double)nSubPixels);
+       
      } else {
        xpos_eta=0.5*dX+0.25;
        ypos_eta=0.5*dY+0.25;
      }
-       
-       int_x=((double)x) + xpos_eta+0.5;
-       int_y=((double)y) +  ypos_eta+0.5;
-    
+     
+     int_x=((double)x) + xpos_eta+0.5;
+     int_y=((double)y) +  ypos_eta+0.5;
+     
 
   }
-
+  
 
   
-   virtual int addToFlatField(double totquad,int quad,int *cl,double &etax, double &etay) {
+  virtual int addToFlatField(double totquad,int quad,int *cl,double &etax, double &etay) {
      double cc[2][2];
-     int *cluster[3];
      int xoff, yoff;
-     cluster[0]=cl;
-     cluster[1]=cl+3;
-     cluster[2]=cl+6;
      
      switch (quad) {
      case BOTTOM_LEFT:
@@ -296,17 +244,11 @@ class eta2InterpolationBase : public virtual etaInterpolationBase {
      default:
        ;
      } 
-     cc[0][0]=cluster[yoff][xoff];
-     cc[1][0]=cluster[yoff+1][xoff];
-     cc[0][1]=cluster[yoff][xoff+1];
-     cc[1][1]=cluster[yoff+1][xoff+1];
+     cc[0][0]=cl[xoff+3*yoff];
+     cc[1][0]=cl[xoff+3*(yoff+1)];
+     cc[0][1]=cl[xoff+1+3*yoff];
+     cc[1][1]=cl[xoff+1+3*(yoff+1)];
      
-      /* cout << cl[0] << " " << cl[1] << " " << cl[2] << endl;   */
-      /* cout << cl[3] << " " << cl[4] << " " << cl[5] << endl;   */
-      /* cout << cl[6] << " " << cl[7] << " " << cl[8] << endl;   */
-      /* cout <<"******"<<totquad << " " << quad << endl;  */
-      /* cout << cc[0][0]<< " " << cc[0][1] << endl;  */
-      /* cout << cc[1][0]<< " " << cc[1][1] << endl;  */
      //calcMyEta(totquad,quad,cl,etax, etay);
      calcEta(totquad, cc,etax, etay);
 
@@ -318,11 +260,7 @@ class eta2InterpolationBase : public virtual etaInterpolationBase {
 
    virtual int addToFlatField(double totquad,int quad,double *cl,double &etax, double &etay) {
      double cc[2][2];
-     double *cluster[3];
      int xoff, yoff;
-     cluster[0]=cl;
-     cluster[1]=cl+3;
-     cluster[2]=cl+6;
      
      switch (quad) {
      case BOTTOM_LEFT:
@@ -344,10 +282,10 @@ class eta2InterpolationBase : public virtual etaInterpolationBase {
      default:
        ;
      } 
-     cc[0][0]=cluster[yoff][xoff];
-     cc[1][0]=cluster[yoff+1][xoff];
-     cc[0][1]=cluster[yoff][xoff+1];
-     cc[1][1]=cluster[yoff+1][xoff+1];
+     cc[0][0]=cl[xoff+3*yoff];
+     cc[1][0]=cl[(yoff+1)*3+xoff];
+     cc[0][1]=cl[yoff*3+xoff+1];
+     cc[1][1]=cl[(yoff+1)*3+xoff+1];
      
       /* cout << cl[0] << " " << cl[1] << " " << cl[2] << endl;   */
       /* cout << cl[3] << " " << cl[4] << " " << cl[5] << endl;   */
@@ -413,7 +351,38 @@ class eta2InterpolationBase : public virtual etaInterpolationBase {
 #endif
     return 0;    
   };
-  
+
+ virtual int *getInterpolatedImage(){
+   int ipx, ipy;
+   // cout << "ff" << endl;
+   calcDiff(1, hhx, hhy); //get flat
+   double avg=0;
+   for (ipx=0; ipx<nSubPixels; ipx++)
+     for (ipy=0; ipy<nSubPixels; ipy++)
+       avg+=flat[ipx+ipy*nSubPixels];
+   avg/=nSubPixels*nSubPixels;
+
+   for (int ibx=0 ; ibx<nSubPixels*nPixelsX; ibx++) {
+     ipx=ibx%nSubPixels-nSubPixels/2;
+       if (ipx<0) ipx=nSubPixels+ipx;
+     for (int iby=0 ; iby<nSubPixels*nPixelsY; iby++) {
+       ipy=iby%nSubPixels-nSubPixels/2;
+     if (ipy<0) ipy=nSubPixels+ipy;
+     // cout << ipx << " " << ipy << " " << ibx << " " << iby << endl;
+     if (flat[ipx+ipy*nSubPixels]>0)
+       hintcorr[ibx+iby*nSubPixels*nPixelsX]=hint[ibx+iby*nSubPixels*nPixelsX]*(avg/flat[ipx+ipy*nSubPixels]);
+     else
+       hintcorr[ibx+iby*nSubPixels*nPixelsX]=hint[ibx+iby*nSubPixels*nPixelsX];
+
+
+     }
+   }
+
+
+
+    return hintcorr;
+  };
+
 /*  protected: */
   
 /* #ifdef MYROOT1 */

slsDetectorCalibration/interpolations/etaInterpolationAdaptiveBins.h

@@ -1,131 +1,137 @@
 #ifndef ETA_INTERPOLATION_ADAPTIVEBINS_H
 #define ETA_INTERPOLATION_ADAPTIVEBINS_H
 
-
+#include <cmath>
 #include "tiffIO.h"
 //#include "etaInterpolationBase.h"
 #include "etaInterpolationPosXY.h"
 
+
+
 class etaInterpolationAdaptiveBins : public etaInterpolationPosXY {
 
- private:
-  double calcDiff(double avg, float *hx, float *hy) {
-    double p_tot=0;
-    double diff=0;
-    double bsize=1./nSubPixels;
-    
-    for (int ipx=0; ipx<nSubPixels; ipx++) {
-      for (int ipy=0; ipy<nSubPixels; ipy++) {
-	p_tot=0;
-	for (int ibx=0; ibx<nbeta; ibx++) {
-	  for (int iby=0; iby<nbeta; iby++) {
-	    if ( hx[ibx+iby*nbeta]>=((ipx)*bsize) && hx[ibx+iby*nbeta]<((ipx+1)*bsize) &&  hy[ibx+iby*nbeta]>=((ipy)*bsize) && hy[ibx+iby*nbeta]<((ipy+1)*bsize)) {
-	      p_tot+=heta[ibx+iby*nbeta];
-	    } 
-	  }
-	}
-	
-	cout <<  p_tot << " \t ";
+  // protected:
 
 
+ private: 
 
-	diff+=(p_tot-avg)*(p_tot-avg);
-	
-      }
-      cout << "\n";
-    }
-    return diff;
-  }
-  
-  void iterate(float *newhhx, float *newhhy) {
+  virtual void iterate(float *newhhx, float *newhhy) {
 
     double bsize=1./nSubPixels;
     
-    double hy[nbeta]; //profile y
-    double hx[nbeta]; //profile x
-    double hix[nbeta]; //integral of projection x
-    double hiy[nbeta]; //integral of projection y
-    
-   double tot_eta_x=0;
-   double tot_eta_y=0;
-    for (int ipy=0; ipy<nSubPixels; ipy++) {
+    double hy[nSubPixels][nbeta]; //profile y
+    double hx[nSubPixels][nbeta]; //profile x
+    double hix[nSubPixels][nbeta]; //integral of projection x
+    double hiy[nSubPixels][nbeta]; //integral of projection y
+    int ipy, ipx;
+   double tot_eta_x[nSubPixels];
+   double tot_eta_y[nSubPixels];
+   //for (int ipy=0; ipy<nSubPixels; ipy++) {
      
-      for (int ibx=0; ibx<nbeta; ibx++) {
-	hx[ibx]=0;
-	hy[ibx]=0;
+       for (ipy=0; ipy<nSubPixels; ipy++) {
+     for (int ibx=0; ibx<nbeta; ibx++) { 
+      	hx[ipy][ibx]=0; 
+      	hy[ipy][ibx]=0;
       }
+       }
 
-      tot_eta_x=0;
-      tot_eta_y=0;
       // cout << ipy << " " << ((ipy)*bsize) << " " << ((ipy+1)*bsize) << endl;
        for (int ibx=0; ibx<nbeta; ibx++) {
 	 for (int iby=0; iby<nbeta; iby++) {
-	   
-	  if (hhy[ibx+iby*nbeta]>=((ipy)*bsize) && hhy[ibx+iby*nbeta]<=((ipy+1)*bsize)) {
-	    hx[ibx]+=heta[ibx+iby*nbeta];
-	    tot_eta_x+=heta[ibx+iby*nbeta];
-	  } 
-	  
-	  
-	  if (hhx[ibx+iby*nbeta]>=((ipy)*bsize) && hhx[ibx+iby*nbeta]<=((ipy+1)*bsize)) {
-	    hy[iby]+=heta[ibx+iby*nbeta];
-	    tot_eta_y+=heta[ibx+iby*nbeta];
-	  } 
-	 }
+	   ipy=hhy[ibx+iby*nbeta]*nSubPixels;
+	   if (ipy<0) ipy=0;
+	   if (ipy>=nSubPixels) ipy=nSubPixels-1;
+	   hx[ipy][ibx]+=heta[ibx+iby*nbeta];
+
+
+	   ipx=hhx[ibx+iby*nbeta]*nSubPixels;
+	   if (ipx<0) ipx=0;
+	   if (ipx>=nSubPixels) ipx=nSubPixels-1;
+	   hy[ipx][iby]+=heta[ibx+iby*nbeta];
+	 } 
        }
       
  
-       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]+hy[ib];
+       
+       for (ipy=0; ipy<nSubPixels; ipy++) {
+	 hix[ipy][0]=hx[ipy][0];
+	 hiy[ipy][0]=hy[ipy][0];
+	 for (int ib=1; ib<nbeta; ib++) {
+	   hix[ipy][ib]=hix[ipy][ib-1]+hx[ipy][ib];
+	   hiy[ipy][ib]=hiy[ipy][ib-1]+hy[ipy][ib];
+	 }
+	 tot_eta_x[ipy]=hix[ipy][nbeta-1]+1;
+	 tot_eta_y[ipy]=hiy[ipy][nbeta-1]+1;
+	 // cout << ipy << " " << tot_eta_x[ipy] << " " << tot_eta_y[ipy] << endl;
+	 
        }
-       // tot_eta_x=hix[nbeta-1];
-       // tot_eta_y=hiy[nbeta-1];
-       /* cout << "ipx " << ipy << " x: " << tot_eta_x << " " << hix[10]<< " " << hix[nbeta-1] << endl; */
-       /* cout << "ipy " << ipy << " y: " << tot_eta_y << " " << hiy[10]<< " " << hiy[nbeta-1] << endl; */
-  
-  // for (int ipy=0; ipy<nSubPixels; ipy++) {
+       // for (int ipy=0; ipy<nSubPixels; ipy++) {
        
        for (int ibx=0; ibx<nbeta; ibx++) {
 	 for (int iby=0; iby<nbeta; iby++) {
 	   
-	   if ( hhy[ibx+iby*nbeta]>=((ipy)*bsize) && hhy[ibx+iby*nbeta]<=((ipy+1)*bsize)) {
-	     newhhx[ibx+iby*nbeta]=hix[ibx]/((double)tot_eta_x);
-	     if (newhhx[ibx+iby*nbeta]>1) cout << "***"<< ibx << " " << iby << newhhx[ibx+iby*nbeta] << endl;
+	   //   if ( hhy[ibx+iby*nbeta]>=((ipy)*bsize) && hhy[ibx+iby*nbeta]<=((ipy+1)*bsize)) {
+	   ipy=hhy[ibx+iby*nbeta]*nSubPixels;
+
+	   if (ipy<0) ipy=0;
+	   if (ipy>=nSubPixels) ipy=nSubPixels-1;
+
+	   if (ipy>=0 && ipy<nSubPixels)
+	     if (tot_eta_x[ipy]>0)
+	       newhhx[ibx+iby*nbeta]=hix[ipy][ibx]/(tot_eta_x[ipy]);
+	     else
+	       cout << "Bad tot_etax " << ipy << " " << tot_eta_x[ipy] << endl;
+	   else
+	     cout << "** Bad value ipy " << ibx << " " << iby << " "<< ipy << " " << hhy[ibx+iby*nbeta]*nSubPixels << endl;
+	   // if (newhhx[ibx+iby*nbeta]>=1 || newhhx[ibx+iby*nbeta]<0 ) cout << "***"<< ibx << " " << iby << newhhx[ibx+iby*nbeta] << endl;
 	     // if (ipy==3 && ibx==10) cout << newhhx[ibx+iby*nbeta] << " " << hix[ibx] << " " << ibx+iby*nbeta << endl;
-	   }
-	   if (hhx[ibx+iby*nbeta]>=((ipy)*bsize) && hhx[ibx+iby*nbeta]<=((ipy+1)*bsize)) {
-	     newhhy[ibx+iby*nbeta]=hiy[iby]/((double)tot_eta_y);
-	     if (newhhy[ibx+iby*nbeta]>1) cout << "***"<< ibx << " " << iby << newhhy[ibx+iby*nbeta] << endl;
+	   // }
+	   ipy=hhx[ibx+iby*nbeta]*nSubPixels;
+	   //if (hhx[ibx+iby*nbeta]>=((ipy)*bsize) && hhx[ibx+iby*nbeta]<=((ipy+1)*bsize)) {
+	   if (ipy<0) ipy=0;
+	   if (ipy>=nSubPixels) ipy=nSubPixels-1;
+
+	   if (ipy>=0 && ipy<nSubPixels)
+	     if (tot_eta_y[ipy]>0)
+	       newhhy[ibx+iby*nbeta]=hiy[ipy][iby]/(tot_eta_y[ipy]);
+	     else
+	       cout << "Bad tot_etay " << ipy << " " << tot_eta_y[ipy] << endl;
+	   else
+	     cout << "** Bad value ipx " << ibx << " " << iby << " "<< ipy << " " << hhx[ibx+iby*nbeta]*nSubPixels << endl;
+	   // if (newhhy[ibx+iby*nbeta]>=1 || newhhy[ibx+iby*nbeta]<0 ) cout << "***"<< ibx << " " << iby << newhhy[ibx+iby*nbeta] << endl;
 	     //  if (ipy==3 && iby==10) cout << newhhy[ibx+iby*nbeta] << " " << hiy[iby] << " " << ibx+iby*nbeta << endl;
-	   }
+	   //   }
 	 }
        }
-  }
+       // }
     
     
   }
   
 
  public:
- etaInterpolationAdaptiveBins(int nx=400, int ny=400, int ns=25, int nb=-1, double emin=1, double emax=0) : etaInterpolationPosXY(nx,ny,ns, nb, emin,emax){};
+ etaInterpolationAdaptiveBins(int nx=400, int ny=400, int ns=25, int nb=-1, double emin=1, double emax=0) : etaInterpolationPosXY(nx,ny,ns, nb, emin,emax){
+    //   flat=new double[nSubPixels*nSubPixels]; flat_x=new double[nSubPixels]; flat_y=new double[nSubPixels];
+    //    flat=new double[nSubPixels*nSubPixels];
+};
 
- etaInterpolationAdaptiveBins(etaInterpolationAdaptiveBins *orig): etaInterpolationPosXY(orig){};
+ etaInterpolationAdaptiveBins(etaInterpolationAdaptiveBins *orig): etaInterpolationPosXY(orig){hintcorr=new int[nPixelsX*nPixelsY*nSubPixels];}; 
 
-  virtual etaInterpolationAdaptiveBins* Clone() {
+  virtual etaInterpolationAdaptiveBins* Clone()=0;
 
-    return new etaInterpolationAdaptiveBins(this);
+  /*   return new etaInterpolationAdaptiveBins(this); */
 
-  };
+  /* }; */
 
 
 
-  virtual void prepareInterpolation(int &ok)
+  virtual void prepareInterpolation(int &ok) {
+    prepareInterpolation(ok, 1000);
+  }
+
+  virtual void prepareInterpolation(int &ok, int nint)
   {
    ok=1;  
-   cout << "Adaptive bins" << endl;
 
    ///*Eta Distribution Rebinning*///
    double bsize=1./nSubPixels; //precision
@@ -153,65 +159,36 @@ class etaInterpolationAdaptiveBins : public etaInterpolationPosXY {
    //  }
    // }
    
-       etaInterpolationPosXY::prepareInterpolation(ok);
+   etaInterpolationPosXY::prepareInterpolation(ok);
 
-#ifdef SAVE_ALL
-   
-   char tit[10000];
-   float *etah=new float[nbeta*nbeta];
-  int etabins=nbeta;
-
-  for (int ii=0; ii<etabins*etabins; ii++) {
-    
-    etah[ii]=hhx[ii];
-  }
-  sprintf(tit,"/scratch/start_hhx.tiff");
-  WriteToTiff(etah, tit, etabins, etabins);
-	  
-  for (int ii=0; ii<etabins*etabins; ii++) {
-    etah[ii]=hhy[ii];
-  }
-  sprintf(tit,"/scratch/start_hhy.tiff");
-  WriteToTiff(etah, tit, etabins, etabins);
-	  
-#endif
-   int nint=1000;
 
    double thr=1./((double)nSubPixels);
    double avg=tot_eta/((double)(nSubPixels*nSubPixels));
-   cout << "total eta entries is :"<< tot_eta << " avg: "<< avg << endl;   
-   cout << "Start " << endl;
+   double rms=sqrt(tot_eta);
+   cout << "total eta entries is :"<< tot_eta << " avg: "<< avg << " rms: " << sqrt(tot_eta) << endl;   
    double old_diff=calcDiff(avg, hhx, hhy), new_diff=old_diff+1, best_diff=old_diff;
-   cout << " diff= " << old_diff << endl;
-   
-
+   // cout << " chi2= " << old_diff << " (rms= " << sqrt(tot_eta) << ")" << endl;
+   cout << endl;
+   cout << endl;
+   debugSaveAll(0);
    int iint=0;
    float *newhhx=new float[nbeta*nbeta]; //profile x
    float *newhhy=new float[nbeta*nbeta]; //profile y
    float *besthhx=hhx; //profile x
    float *besthhy=hhy; //profile y
-   while (iint<nint) {
-     
-     cout << "Iteration " << iint << endl;
+
+   cout << "Iteration "<< iint << " Chi2: " << old_diff <<  endl; //" Best: "<< best_diff << " RMS: "<< rms<< endl;
+   while (iint<nint && best_diff > rms) {
+
+/* #ifdef SAVE_ALL */
+/*      if (iint%10==0) */
+/*        debugSaveAll(iint); */
+/* #endif */
+     //  cout << "Iteration " << iint << endl;
      iterate(newhhx,newhhy);
      new_diff=calcDiff(avg, newhhx, newhhy);
-     cout << " diff= " << new_diff << endl;
-/* #ifdef SAVE_ALL */
-/*   for (int ii=0; ii<etabins*etabins; ii++) { */
-/*     etah[ii]=newhhx[ii]; */
-/*     if (etah[ii]>1 || etah[ii]<0  ) cout << "***"<< ii << etah[ii] << endl; */
-   
-/*   } */
-/*   sprintf(tit,"/scratch/neweta_hhx_%d.tiff",iint); */
-/*   WriteToTiff(etah, tit, etabins, etabins); */
-	  
-/*   for (int ii=0; ii<etabins*etabins; ii++) { */
-/*     etah[ii]=newhhy[ii]; */
-/*     if (etah[ii]>1 || etah[ii]<0  ) cout << "***"<< ii << etah[ii] << endl; */
-/*   } */
-/*   sprintf(tit,"/scratch/neweta_hhy_%d.tiff",iint); */
-/*   WriteToTiff(etah, tit, etabins, etabins); */
-/* #endif */
+     //   cout << " chi2= " << new_diff << " (rms= " << sqrt(tot_eta) << ")"<<endl;
+
   if (new_diff<best_diff) {
     best_diff=new_diff;
     besthhx=newhhx;
@@ -225,19 +202,33 @@ class etaInterpolationAdaptiveBins : public etaInterpolationPosXY {
   
   hhx=newhhx;
   hhy=newhhy;
+
+  
+#ifdef SAVE_ALL
+  if (new_diff<=best_diff) {
+       debugSaveAll(iint);
+  }
+#endif
+
+
   newhhx=new float[nbeta*nbeta]; //profile x
   newhhy=new float[nbeta*nbeta]; //profile y
 
-  old_diff=new_diff;
-  //} /* else { */
-     /*   cout << "Difference not decreasing after "<< iint << " iterations (" << old_diff << " < " << new_diff << ")"<< endl; */
-     /*   break; */
-     /* } */
   
-  iint++;
+    /* if (new_diff<old_diff){ */
+    /*   cout << "best difference at iteration "<< iint << " (" << new_diff << " < " << old_diff << ")"<< "Best: "<< best_diff << " RMS: "<< sqrt(tot_eta) << endl; */
+    /*   ; */
+    /* } else {  */
+     // break;
+     // }
+  
+    old_diff=new_diff;
+
+    iint++;
+    cout << "Iteration "<< iint << " Chi2: " << new_diff <<  endl; //" Best: "<< best_diff << " RMS: "<< rms<< endl;
    }
-     delete [] newhhx;
-     delete [] newhhy;
+   delete [] newhhx;
+   delete [] newhhy;
      
   if (hhx!=besthhx)
     delete [] hhx;
@@ -247,35 +238,48 @@ class etaInterpolationAdaptiveBins : public etaInterpolationPosXY {
   hhx=besthhx;
   hhy=besthhy;
 
-   
-
-
-
+  
+    cout << "Iteration "<< iint << " Chi2: " << best_diff <<  endl; //" Best: "<< best_diff << " RMS: "<< rms<< endl;
 #ifdef SAVE_ALL
-   
-
-  for (int ii=0; ii<etabins*etabins; ii++) {
-    etah[ii]=hhx[ii];
-  }
-  sprintf(tit,"/scratch/eta_hhx_%d.tiff",id);
-  WriteToTiff(etah, tit, etabins, etabins);
-	  
-  for (int ii=0; ii<etabins*etabins; ii++) {
-    etah[ii]=hhy[ii];
-  }
-  sprintf(tit,"/scratch/eta_hhy_%d.tiff",id);
-  WriteToTiff(etah, tit, etabins, etabins);
-	  
-  for (int ii=0; ii<etabins*etabins; ii++) {
-    etah[ii]=heta[ii];
-  }
-  sprintf(tit,"/scratch/eta_%d.tiff",id);
-  WriteToTiff(etah, tit, etabins, etabins);
-  delete [] etah;
+  debugSaveAll(iint);
 #endif
+
   return ;
   }
 
+
+
+
+
+};
+
+class eta2InterpolationAdaptiveBins : public virtual eta2InterpolationBase, public virtual etaInterpolationAdaptiveBins {
+ public:
+ eta2InterpolationAdaptiveBins(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),eta2InterpolationBase(nx,ny,ns, nb, emin,emax),etaInterpolationAdaptiveBins(nx,ny,ns, nb, emin,emax){
+    cout << "NSUBPIX is " <<  ns << " " << nSubPixels << endl;
+    //  cout << "e2pxy " << nb << " " << emin << " " << emax << endl; 
+  };
+ 
+ eta2InterpolationAdaptiveBins(eta2InterpolationAdaptiveBins *orig): etaInterpolationBase(orig), etaInterpolationAdaptiveBins(orig)  {};
+
+  virtual eta2InterpolationAdaptiveBins* Clone() { return new eta2InterpolationAdaptiveBins(this);};
+
+
+
+
+};
+
+
+
+class eta3InterpolationAdaptiveBins : public virtual eta3InterpolationBase, public virtual etaInterpolationAdaptiveBins {
+ public:
+ eta3InterpolationAdaptiveBins(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),eta3InterpolationBase(nx,ny,ns, nb, emin,emax), etaInterpolationAdaptiveBins(nx,ny,ns, nb, emin,emax){
+      cout << "e3pxy " << nbeta << " " << etamin << " " << etamax << " " << nSubPixels<< endl;  
+  };
+
+ eta3InterpolationAdaptiveBins(eta3InterpolationAdaptiveBins *orig): etaInterpolationBase(orig), etaInterpolationAdaptiveBins(orig)  {};
+
+  virtual eta3InterpolationAdaptiveBins* Clone() { return new eta3InterpolationAdaptiveBins(this);};
 };
 
 #endif

slsDetectorCalibration/interpolations/etaInterpolationBase.h

@@ -23,97 +23,53 @@ class etaInterpolationBase : public slsInterpolation {
       nbeta=nSubPixels*10;
     }
     if (etamin>=etamax) {
-      // cout << "aaa:" <<endl;
       etamin=-1;
       etamax=2;
     }
     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 float[nbeta*nbeta];
     hhy=new float[nbeta*nbeta];
+    rangeMin=etamin;
+    rangeMax=etamax;
+    flat= new double[nSubPixels*nSubPixels];
+    hintcorr=new int [nSubPixels*nSubPixels*nPixelsX*nPixelsY];
     
-#endif
-    
-    //cout << nbeta << " " << etamin << " " << etamax << endl;
   };
   
  etaInterpolationBase(etaInterpolationBase *orig): slsInterpolation(orig){
    nbeta=orig->nbeta;
    etamin=orig->etamin;
    etamax=orig->etamax;
+   rangeMin=orig->rangeMin;
+   rangeMax=orig->rangeMax;
+
 
     etastep=(etamax-etamin)/nbeta;
-#ifdef MYROOT1
-    heta=(TH2D*)(orig->heta)->Clone("heta");
-    hhx=(TH2D*)(orig->hhx)->Clone("hhx");
-    hhy=(TH2D*)(orig->hhy)->Clone("hhy");
-#endif
-
-#ifndef MYROOT1
     heta=new int[nbeta*nbeta];
     memcpy(heta,orig->heta,nbeta*nbeta*sizeof(int));
     hhx=new float[nbeta*nbeta];
     memcpy(hhx,orig->hhx,nbeta*nbeta*sizeof(float));
     hhy=new float[nbeta*nbeta];
     memcpy(hhy,orig->hhy,nbeta*nbeta*sizeof(float));
-    
-#endif
-    
+    hintcorr=new int [nSubPixels*nSubPixels*nPixelsX*nPixelsY];
 
  };
 
   
- /* virtual etaInterpolationBase* Clone()=0; {
-    return new etaInterpolationBase(this);
-    };*/
- 
 
 
 
   virtual void resetFlatField() {
-#ifndef MYROOT1
     for (int ibx=0; ibx<nbeta*nbeta; ibx++) {
       heta[ibx]=0;
       hhx[ibx]=0;
       hhy[ibx]=0;
     }
-#endif
 
-#ifdef MYROOT1
-    heta->Reset();
-    hhx->Reset();
-    hhy->Reset();
-#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 *setFlatField(TH2D *h, int nb=-1, double emin=1, double emax=0)
-  {  
-    return setEta(h, nb, emin, emax);
-  };
-  
-  TH2D *getFlatField(){return setEta(NULL);};
-#endif
-  
-#ifndef MYROOT1
   int *setEta(int *h, int nb=-1, double emin=1, double emax=0)
   {  
     if (h) {
@@ -127,6 +83,8 @@ class etaInterpolationBase : public slsInterpolation {
 	etamin=-1;
 	etamax=2;
       }
+      rangeMin=etamin;
+      rangeMax=etamax;
       etastep=(etamax-etamin)/nbeta;
     }
     return heta;
@@ -143,7 +101,6 @@ class etaInterpolationBase : public slsInterpolation {
   
   int *getFlatField(int &nb, double &emin, double &emax){
     nb=nbeta; 
-    //cout << "igff* ff has " << nb << " bins " << endl; 
     emin=etamin; 
     emax=etamax; 
     return getFlatField();
@@ -208,28 +165,6 @@ class etaInterpolationBase : public slsInterpolation {
   
     
 
-
-#endif
-  
-
- 
-  /* ////////////////////////////////////////////////////////////////////////////// */
-
-#ifdef MYROOT1
-TH2D *gethhx()
-  {
-    hhx->Scale((double)nSubPixels);
-    return hhx;
-  };
-  
-  TH2D *gethhy()
-  {
-    hhy->Scale((double)nSubPixels);
-    return hhy;
-    };
-#endif
-    
-#ifndef MYROOT1
 float *gethhx()
   {
     // hhx->Scale((double)nSubPixels);
@@ -241,53 +176,193 @@ float *gethhx()
     // hhy->Scale((double)nSubPixels);
     return hhy;
     };
-#endif
-    
-  //////////////////////////////////////////////////////////////////////////////
-  //////////// /*It return position hit for the event in input */ //////////////
-  /* virtual void getInterpolatedPosition(int x, int y, int *data, double &int_x, double &int_y)=0; */
-  /* virtual void getInterpolatedPosition(int x, int y, double *data, double &int_x, double &int_y)=0; */
-  /* virtual void getInterpolatedPosition(int x, int y, double totquad,int quad,double *cl,double &int_x, double &int_y)=0; */
-  /* virtual void getInterpolatedPosition(int x, int y, double totquad,int quad,int *cl,double &int_x, double &int_y)=0;  */
-  /* virtual void getInterpolatedPosition(int x, int y, double etax, double etay, int corner, double &int_x, double &int_y)=0; */
-
-
-  /* virtual int addToFlatField(double totquad,int quad,int *cl,double &etax, double &etay)=0; */
-  /* virtual int addToFlatField(double totquad,int quad,double *cl,double &etax, double &etay)=0; */
-  /* virtual int addToFlatField(double *cluster, double &etax, double &etay)=0; */
-  /* virtual int addToFlatField(int *cluster, double &etax, double &etay)=0; */
-  
   virtual int addToFlatField(double etax, double etay){
-#ifdef MYROOT1
-    heta->Fill(etax,etay);
-#endif
-#ifndef MYROOT1
     int ex,ey; 
     ex=(etax-etamin)/etastep;
     ey=(etay-etamin)/etastep;
     if (ey<nbeta && ex<nbeta && ex>=0 && ey>=0)
       heta[ey*nbeta+ex]++; 
-#endif
     return 0;    
   };
   
   
   // virtual void prepareInterpolation(int &ok)=0;
 
+
+  void debugSaveAll(int ind=0) {
+    int ib, ibx, iby;
+   char tit[10000];
+   
+   float tot_eta=0;
+
+  float *etah=new float[nbeta*nbeta];
+  int etabins=nbeta;
+  int ibb=0;
+
+  for (int ii=0; ii<etabins*etabins; ii++) {
+   
+      etah[ii]=heta[ii];
+      tot_eta+=heta[ii];
+  }
+  sprintf(tit,"/scratch/eta.tiff",ind);
+  WriteToTiff(etah, tit, etabins, etabins);
+
+
+  for (int ii=0; ii<etabins*etabins; ii++) {
+    ibb=(hhx[ii]*nSubPixels);
+    etah[ii]=ibb;
+  }
+  sprintf(tit,"/scratch/eta_hhx_%d.tiff",ind);
+  WriteToTiff(etah, tit, etabins, etabins);
+	  
+  for (int ii=0; ii<etabins*etabins; ii++) {
+    ibb=hhy[ii]*nSubPixels;
+    etah[ii]=ibb;
+  }
+  sprintf(tit,"/scratch/eta_hhy_%d.tiff",ind);
+  WriteToTiff(etah, tit, etabins, etabins);
+	  
+  
+  float *ftest=new float[nSubPixels*nSubPixels];
+
+  for (int ib=0; ib<nSubPixels*nSubPixels; ib++) ftest[ib]=0;
+  
+
+  //int ibx=0, iby=0;
+  
+  for (int ii=0; ii<nbeta*nbeta; ii++) {
+    
+    ibx=nSubPixels*hhx[ii];
+    iby=nSubPixels*hhy[ii];
+    if (ibx<0) ibx=0;
+    if (iby<0) iby=0;
+    if (ibx>=nSubPixels) ibx=nSubPixels-1;
+    if (iby>=nSubPixels) iby=nSubPixels-1;
+    
+
+    if (ibx>=0 && ibx<nSubPixels && iby>=0 && iby<nSubPixels) {
+      //
+      // if (ibx>0 && iby>0) cout << ibx << " " << iby << " " << ii << endl;
+      ftest[ibx+iby*nSubPixels]+=heta[ii];
+    } else
+      cout << "Bad interpolation "<< ii << " " << ibx << " " << iby<< endl; 
+
+  }
+
+  sprintf(tit,"/scratch/ftest_%d.tiff",ind);
+  WriteToTiff(ftest, tit, nSubPixels, nSubPixels);
+
+  //int ibx=0, iby=0;
+  tot_eta/=nSubPixels*nSubPixels;
+  int nbad=0;
+  for (int ii=0; ii<etabins*etabins; ii++) {
+    ibx=nSubPixels*hhx[ii];
+    iby=nSubPixels*hhy[ii];
+    if (ftest[ibx+iby*nSubPixels]<tot_eta*0.5) {
+      etah[ii]=1;
+      nbad++;
+    } else if(ftest[ibx+iby*nSubPixels]>tot_eta*2.){
+      etah[ii]=2;
+      nbad++;
+    } else
+      etah[ii]=0;
+  }
+  sprintf(tit,"/scratch/eta_bad_%d.tiff",ind);
+  WriteToTiff(etah, tit, etabins, etabins);
+  // cout << "Index: " << ind << "\t Bad bins: "<< nbad << endl;
+  //int ibx=0, iby=0;
+
+  delete [] ftest;
+  delete [] etah;
+
+  }
+
+
  protected:
   
-#ifdef MYROOT1
-  TH2D *heta;
-  TH2D *hhx;
-  TH2D *hhy;
-#endif
-#ifndef MYROOT1
+
+  double calcDiff(double avg, float *hx, float *hy) {
+    //double p_tot=0;
+    double diff=0, d;
+    double bsize=1./nSubPixels;
+    int nbad=0;
+    double p_tot_x[nSubPixels], p_tot_y[nSubPixels], p_tot[nSubPixels*nSubPixels];
+    double maxdiff=0, mindiff=avg*nSubPixels*nSubPixels;
+   
+    int ipx, ipy;
+    for (ipy=0; ipy<nSubPixels; ipy++) {
+      for (ipx=0; ipx<nSubPixels; ipx++) {
+	p_tot[ipx+ipy*nSubPixels]=0;
+      }
+      p_tot_y[ipy]=0;
+      p_tot_x[ipy]=0;
+    }
+   
+    for (int ibx=0; ibx<nbeta; ibx++) {
+      for (int iby=0; iby<nbeta; iby++) {
+	ipx=hx[ibx+iby*nbeta]*nSubPixels;
+	if (ipx<0) ipx=0;
+	if (ipx>=nSubPixels) ipx=nSubPixels-1;
+	
+	ipy=hy[ibx+iby*nbeta]*nSubPixels;
+	if (ipy<0) ipy=0;
+	if (ipy>=nSubPixels) ipy=nSubPixels-1;
+	
+	p_tot[ipx+ipy*nSubPixels]+=heta[ibx+iby*nbeta];
+	p_tot_y[ipy]+=heta[ibx+iby*nbeta];
+	p_tot_x[ipx]+=heta[ibx+iby*nbeta];
+	
+
+      }
+    }
+    
+    
+    //  cout << endl << endl;
+    for (ipy=0; ipy<nSubPixels; ipy++) { 
+     cout.width(5);
+     //flat_y[ipy]=p_tot_y[ipy];//avg/nSubPixels;
+      for (ipx=0; ipx<nSubPixels; ipx++) {
+
+	//	flat_x[ipx]=p_tot_x[ipx];///avg/nSubPixels;
+	flat[ipx+nSubPixels*ipy]=p_tot[ipx+nSubPixels*ipy];///avg;
+	d=p_tot[ipx+nSubPixels*ipy]-avg;
+	if (d<0) d*=-1.;
+	if (d>5*sqrt(avg) )
+	  nbad++;
+	diff+=d*d;
+	if (d<mindiff) mindiff=d;
+	if (d>maxdiff) maxdiff=d;
+	//	cout << setprecision(4) << p_tot[ipx+nSubPixels*ipy] << "   ";
+      }
+
+      /* cout << "** "  << setprecision(4) <<  flat_y[ipy]; */
+      //cout << "\n";
+    }
+    /* cout << "**" << endl; cout.width(5); */
+    /* for (ipx=0; ipx<nSubPixels; ipx++) { */
+    /*   cout  << setprecision(4) <<  flat_x[ipx] << " "; */
+    /* } */
+    //cout << "**" << endl; cout.width(5);
+    //cout << "Min diff: " << mindiff/sqrt(avg) << " Max diff: " << maxdiff/sqrt(avg) << " Nbad: " << nbad << endl;
+    
+
+
+    //   cout << "Bad pixels: " << 100.*(float)nbad/((float)(nSubPixels*nSubPixels)) << " %" << endl;
+    return sqrt(diff);
+  }
+  
   int *heta;
   float *hhx;
   float *hhy;
-#endif
   int nbeta;
   double etamin, etamax, etastep;
+  double rangeMin, rangeMax;
+
+
+
+  double *flat;
+  int *hintcorr;
+  
 
 };
 

slsDetectorCalibration/interpolations/etaInterpolationCleverAdaptiveBins.h

@@ -0,0 +1,263 @@
+#ifndef ETA_INTERPOLATION_CLEVER_ADAPTIVEBINS_H
+#define ETA_INTERPOLATION_CLEVER_ADAPTIVEBINS_H
+
+#include <cmath>
+#include "tiffIO.h"
+//#include "etaInterpolationBase.h"
+#include "etaInterpolationAdaptiveBins.h"
+
+
+//#define HSIZE 1
+
+
+class etaInterpolationCleverAdaptiveBins : public etaInterpolationAdaptiveBins {
+
+ private:
+
+
+  //  double *gradientX, *gradientY, *gradientXY;
+
+
+  virtual void iterate(float *newhhx, float *newhhy) {
+
+    double bsize=1./nSubPixels;
+    
+    /* double hy[nSubPixels*HSIZE][nbeta]; //profile y */
+    /* double hx[nSubPixels*HSIZE][nbeta]; //profile x */
+    //  double hix[nSubPixels*HSIZE][nbeta]; //integral of projection x
+    // double hiy[nSubPixels*HSIZE][nbeta]; //integral of projection y
+    int ipy, ipx, ippx, ippy;
+    // double tot_eta_x[nSubPixels*HSIZE];
+    //double tot_eta_y[nSubPixels*HSIZE];
+
+    double mean=0;
+    double maxflat=0, minflat=0, maxgradX=0, mingradX=0, maxgradY=0, mingradY=0, maxgr=0, mingr=0;
+
+    int ix_maxflat, iy_maxflat, ix_minflat, iy_minflat, ix_maxgrX, iy_maxgrX, ix_mingrX, iy_mingrX,ix_maxgrY, iy_maxgrY, ix_mingrY, iy_mingrY, ix_mingr, iy_mingr, ix_maxgr, iy_maxgr; 
+    int maskMin[nSubPixels*nSubPixels], maskMax[nSubPixels*nSubPixels];
+
+
+    //for (int ipy=0; ipy<nSubPixels; ipy++) {
+     
+    for (ipy=0; ipy<nSubPixels; ipy++) {
+      for (ipx=0; ipx<nSubPixels; ipx++) {
+	//	cout << ipx << " " << ipy << endl;
+	mean+=flat[ipx+nSubPixels*ipy]/((double)(nSubPixels*nSubPixels));
+      }
+    }
+
+    //  cout << "Mean is " << mean << endl;
+
+    /*** Find local minima and maxima within the staistical uncertainty **/
+    for (ipy=0; ipy<nSubPixels; ipy++) {
+      for (ipx=0; ipx<nSubPixels; ipx++) {
+	if (flat[ipx+nSubPixels*ipy]<mean-3.*sqrt(mean))maskMin[ipx+nSubPixels*ipy]=1; else maskMin[ipx+nSubPixels*ipy]=0;
+	if (flat[ipx+nSubPixels*ipy]>mean+3.*sqrt(mean)) maskMax[ipx+nSubPixels*ipy]=1; else maskMax[ipx+nSubPixels*ipy]=0;
+	if (ipx>0 && ipy>0) {
+	  if (flat[ipx+nSubPixels*ipy]<flat[ipx-1+nSubPixels*(ipy-1)]) 	maskMax[ipx+nSubPixels*ipy]=0;
+	  if (flat[ipx+nSubPixels*ipy]>flat[ipx-1+nSubPixels*(ipy-1)]) 	maskMin[ipx+nSubPixels*ipy]=0;
+	}
+	if (ipx>0 && ipy<nSubPixels-1) {
+	  if (flat[ipx+nSubPixels*ipy]<flat[ipx-1+nSubPixels*(ipy+1)]) 	maskMax[ipx+nSubPixels*ipy]=0;
+	  if (flat[ipx+nSubPixels*ipy]>flat[ipx-1+nSubPixels*(ipy+1)]) 	maskMin[ipx+nSubPixels*ipy]=0;
+	}
+	if (ipy>0 && ipx<nSubPixels-1) {
+	  if (flat[ipx+nSubPixels*ipy]<flat[ipx+1+nSubPixels*(ipy-1)]) 	maskMax[ipx+nSubPixels*ipy]=0;
+	  if (flat[ipx+nSubPixels*ipy]>flat[ipx+1+nSubPixels*(ipy-1)]) 	maskMin[ipx+nSubPixels*ipy]=0;
+	}
+	if (ipy<nSubPixels-1 && ipx<nSubPixels-1) {
+	  if (flat[ipx+nSubPixels*ipy]<flat[ipx+1+nSubPixels*(ipy+1)]) 	maskMax[ipx+nSubPixels*ipy]=0;
+	  if (flat[ipx+nSubPixels*ipy]>flat[ipx+1+nSubPixels*(ipy+1)]) 	maskMin[ipx+nSubPixels*ipy]=0;
+	}
+	if (ipy<nSubPixels-1 ) {
+	  if (flat[ipx+nSubPixels*ipy]<flat[ipx+nSubPixels*(ipy+1)]) 	maskMax[ipx+nSubPixels*ipy]=0;
+	  if (flat[ipx+nSubPixels*ipy]>flat[ipx+nSubPixels*(ipy+1)]) 	maskMin[ipx+nSubPixels*ipy]=0;
+	}
+	if (ipx<nSubPixels-1) {
+	  if (flat[ipx+nSubPixels*ipy]<flat[ipx+1+nSubPixels*(ipy)]) 	maskMax[ipx+nSubPixels*ipy]=0;
+	  if (flat[ipx+nSubPixels*ipy]>flat[ipx+1+nSubPixels*(ipy)]) 	maskMin[ipx+nSubPixels*ipy]=0;
+	}
+	
+	if (ipy>0 ) {
+	  if (flat[ipx+nSubPixels*ipy]<flat[ipx+nSubPixels*(ipy-1)]) 	maskMax[ipx+nSubPixels*ipy]=0;
+	  if (flat[ipx+nSubPixels*ipy]>flat[ipx+nSubPixels*(ipy-1)]) 	maskMin[ipx+nSubPixels*ipy]=0;
+	}
+
+	if (ipx>0 ) {
+	  if (flat[ipx+nSubPixels*ipy]<flat[ipx-1+nSubPixels*(ipy)]) 	maskMax[ipx+nSubPixels*ipy]=0;
+	  if (flat[ipx+nSubPixels*ipy]>flat[ipx-1+nSubPixels*(ipy)]) 	maskMin[ipx+nSubPixels*ipy]=0;
+	}
+
+	//	if (maskMin[ipx+nSubPixels*ipy]) cout << ipx << " " << ipy << " is a local minimum " << flat[ipx+nSubPixels*ipy] << endl;
+	//	if (maskMax[ipx+nSubPixels*ipy]) cout << ipx << " " << ipy << " is a local maximum "<< flat[ipx+nSubPixels*ipy] << endl;
+
+      }
+    }
+    int is_a_border=0;
+
+    //initialize the new partition to the previous one
+    // int ibx_p, iby_p, ibx_n, iby_n;
+    int ibbx, ibby;
+
+  
+    memcpy(newhhx,hhx,nbeta*nbeta*sizeof(float));
+    memcpy(newhhy,hhy,nbeta*nbeta*sizeof(float));
+
+    for (int ibx=0; ibx<nbeta; ibx++) {
+      for (int iby=0; iby<nbeta; iby++) { 
+ 
+	ippy=hhy[ibx+iby*nbeta]*nSubPixels;
+	ippx=hhx[ibx+iby*nbeta]*nSubPixels;
+	
+	
+	is_a_border=0;
+	
+	if (maskMin[ippx+nSubPixels*ippy] || maskMax[ippx+nSubPixels*ippy]) {
+	  
+
+	  for (int ix=-1; ix<2; ix++) {
+	    ibbx=ibx+ix;
+	    if (ibbx<0) ibbx=0;
+	    if (ibbx>nbeta-1) ibbx=nbeta-1;
+	    for (int iy=-1; iy<2; iy++) {
+	      ibby=iby+iy;
+	      if (ibby<0) ibby=0;
+	      if (ibby>nbeta-1) ibby=nbeta-1;
+	      
+	      
+	      ipy=hhy[ibbx+ibby*nbeta]*nSubPixels;
+	      ipx=hhx[ibbx+ibby*nbeta]*nSubPixels;
+
+	      
+	      if (ipx!=ippx || ipy!=ippy) {
+		is_a_border=1;
+		if (maskMin[ippx+nSubPixels*ippy]) {
+		  //increase the region
+		  newhhx[ibbx+ibby*nbeta]=((double)ippx+0.5)/((double)nSubPixels);
+		  newhhy[ibbx+ibby*nbeta]=((double)ippy+0.5)/((double)nSubPixels);
+		}
+		if (maskMax[ippx+nSubPixels*ippy]) {
+		  //reduce the region
+		  newhhx[ibx+iby*nbeta]=((double)ipx+0.5)/((double)nSubPixels);
+		  newhhy[ibx+iby*nbeta]=((double)ipy+0.5)/((double)nSubPixels);
+		}
+
+		//	cout << ippx << " " << ippy << " " << ibx << " " << iby << " * " << ipx << " " << ipy << " " << ibbx << " " << ibby << endl;
+	      }
+
+
+	    }
+	  }
+	} 	
+
+
+
+      }
+
+    }
+  
+
+    //Check that the resulting histograms are monotonic and they don't have holes!
+  
+    for (int ibx=0; ibx<nbeta-1; ibx++) {
+      for (int iby=0; iby<nbeta-1; iby++) {
+	
+    	ippy=newhhy[ibx+iby*nbeta]*nSubPixels;
+    	ippx=newhhx[ibx+iby*nbeta]*nSubPixels;
+	
+    	ipy=newhhy[ibx+(iby+1)*nbeta]*nSubPixels;
+    	ipx=newhhx[ibx+1+iby*nbeta]*nSubPixels;
+
+    	if ( ippx>ipx)
+    	  newhhx[ibx+1+iby*nbeta]=newhhx[ibx+iby*nbeta];
+    	else if (ipx >ippx+1)
+    	  newhhx[ibx+1+iby*nbeta]=((double)(ippx+1+0.5))/((double)nSubPixels);
+	  
+    	if ( ippy>ipy)
+    	  newhhy[ibx+(iby+1)*nbeta]=newhhy[ibx+iby*nbeta];
+    	else if (ipy >ippy+1)
+    	  newhhy[ibx+(iby+1)*nbeta]=((double)(ippy+1+0.5))/((double)nSubPixels);
+	  
+      }
+    }
+
+    
+  }
+  
+
+ public:
+ etaInterpolationCleverAdaptiveBins(int nx=400, int ny=400, int ns=25, int nb=-1, double emin=1, double emax=0) : etaInterpolationAdaptiveBins(nx,ny,ns, nb, emin,emax){
+  
+    
+  };
+
+ etaInterpolationCleverAdaptiveBins(etaInterpolationCleverAdaptiveBins *orig): etaInterpolationAdaptiveBins(orig){}; 
+
+  virtual etaInterpolationCleverAdaptiveBins* Clone()=0;
+
+  /*   return new etaInterpolationCleverAdaptiveBins(this); */
+
+  /* }; */
+
+
+};
+
+class eta2InterpolationCleverAdaptiveBins : public virtual eta2InterpolationBase, public virtual etaInterpolationCleverAdaptiveBins {
+ public:
+ eta2InterpolationCleverAdaptiveBins(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),eta2InterpolationBase(nx,ny,ns, nb, emin,emax),etaInterpolationCleverAdaptiveBins(nx,ny,ns, nb, emin,emax){
+  };
+ 
+ eta2InterpolationCleverAdaptiveBins(eta2InterpolationCleverAdaptiveBins *orig): etaInterpolationBase(orig), etaInterpolationCleverAdaptiveBins(orig)  {};
+
+  virtual eta2InterpolationCleverAdaptiveBins* Clone() { return new eta2InterpolationCleverAdaptiveBins(this);};
+
+  // virtual int *getInterpolatedImage(){return eta2InterpolationBase::getInterpolatedImage();};
+
+  /* virtual int *getInterpolatedImage(){ */
+  /*   int ipx, ipy; */
+  /*   cout << "ff" << endl; */
+  /*   calcDiff(1, hhx, hhy); //get flat */
+  /*   double avg=0; */
+  /*  for (ipx=0; ipx<nSubPixels; ipx++) */
+  /*    for (ipy=0; ipy<nSubPixels; ipy++) */
+  /*      avg+=flat[ipx+ipy*nSubPixels]; */
+  /*  avg/=nSubPixels*nSubPixels; */
+   
+  /*  for (int ibx=0 ; ibx<nSubPixels*nPixelsX; ibx++) { */
+  /*    ipx=ibx%nSubPixels-nSubPixels; */
+  /*    if (ipx<0) ipx=nSubPixels+ipx; */
+  /*    for (int iby=0 ; iby<nSubPixels*nPixelsY; iby++) { */
+  /*    ipy=iby%nSubPixels-nSubPixels; */
+  /*    if (ipy<0) ipy=nSubPixels+ipy; */
+     
+  /*    if (flat[ipx+ipy*nSubPixels]>0) */
+  /*      hintcorr[ibx+iby*nSubPixels*nPixelsX]=hint[ibx+iby*nSubPixels*nPixelsX]*(avg/flat[ipx+ipy*nSubPixels]); */
+  /*    else */
+  /*      hintcorr[ibx+iby*nSubPixels*nPixelsX]=hint[ibx+iby*nSubPixels*nPixelsX]; */
+
+
+  /*    } */
+  /*  } */
+
+
+
+  /*   return hintcorr; */
+  /* }; */
+
+
+};
+
+
+
+class eta3InterpolationCleverAdaptiveBins : public virtual eta3InterpolationBase, public virtual etaInterpolationCleverAdaptiveBins {
+ public:
+ eta3InterpolationCleverAdaptiveBins(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),eta3InterpolationBase(nx,ny,ns, nb, emin,emax), etaInterpolationCleverAdaptiveBins(nx,ny,ns, nb, emin,emax){
+     
+  };
+
+ eta3InterpolationCleverAdaptiveBins(eta3InterpolationCleverAdaptiveBins *orig): etaInterpolationBase(orig), etaInterpolationCleverAdaptiveBins(orig)  {};
+
+  virtual eta3InterpolationCleverAdaptiveBins* Clone() { return new eta3InterpolationCleverAdaptiveBins(this);};
+};
+
+#endif

slsDetectorCalibration/interpolations/etaInterpolationPosXY.h

@@ -51,16 +51,26 @@ class etaInterpolationPosXY : public virtual etaInterpolationBase{
    double hix[nbeta]; //integral of projection x
    double hiy[nbeta]; //integral of projection y
    int ii=0;
+   double etax, etay;
    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];
+       etax=etamin+iby*etastep;
+       //cout << etax << endl;
+       if (etax>=0 && etax<=1)
+	 hx[iby]=heta[iby+ib*nbeta];
+       else {
+	 hx[iby]=0;
+       }
+       // tot_eta_x+=hx[iby];
+       if (etax>=0 && etax<=1)
+	 hy[iby]=heta[ib+iby*nbeta];
+       else
+	 hy[iby]=0;
+       // tot_eta_y+=hy[iby];
      }
 
      hix[0]=hx[0];
@@ -72,22 +82,17 @@ class etaInterpolationPosXY : public virtual etaInterpolationBase{
      }
 
      ii=0;
+     tot_eta_x=hix[nbeta-1]+1;
+     tot_eta_y=hiy[nbeta-1]+1;
      
      for (int ibx=0; ibx<nbeta; ibx++) {
-       if (tot_eta_x==0) {
-	 hhx[ibx+ib*nbeta]=((float)ibx)/((float)nbeta);
-	 ii=(ibx)/nbeta;
+       if (tot_eta_x<=0) {
+	 hhx[ibx+ib*nbeta]=-1;
+	 //ii=(ibx)/nbeta;
        } else //if (hix[ibx]>(ii+1)*tot_eta_x*bsize) 
-	 {
-	 //ii++;
-	 // cout << ib << " x " << ibx << " " << tot_eta_x << " " << (ii)*tot_eta_x*bsize << " " << ii << endl;
-	 // }
-#ifdef MYROOT1  
-       hhx->SetBinContent(ibx+1,ib+1,ii);
-#endif
-#ifndef MYROOT1  
-       hhx[ibx+ib*nbeta]=hix[ibx]/((float)tot_eta_x);//ii; 
-#endif
+	 { 
+	   //if (hix[ibx]>tot_eta_x*(ii+1)/nSubPixels) ii++;
+	   hhx[ibx+ib*nbeta]=hix[ibx]/tot_eta_x;
 	 }
      }
      /* if (ii!=(nSubPixels-1)) */
@@ -96,53 +101,55 @@ class etaInterpolationPosXY : public virtual etaInterpolationBase{
      ii=0;
      
      for (int ibx=0; ibx<nbeta; ibx++) {
-       if (tot_eta_y==0) {
-	 hhx[ibx+ib*nbeta]=((float)ibx)/((float)nbeta);
-	 ii=(ibx*nSubPixels)/nbeta;
-       } else //if (hiy[ibx]>(ii+1)*tot_eta_y*bsize) 
-	 {
-	   //ii++;
-	 //cout << ib << " y " << ibx << " " << tot_eta_y << " "<< (ii)*tot_eta_y*bsize << " " << ii << endl;
-	   //}
-#ifdef MYROOT1  
-       hhy->SetBinContent(ib+1,ibx+1,ii);
-#endif
-#ifndef MYROOT1  
-       hhy[ib+ibx*nbeta]=hiy[ibx]/((float)tot_eta_y);//ii; 
-#endif
-	 }
+       if (tot_eta_y<=0) {
+	 hhy[ib+ibx*nbeta]=-1;
+	 //ii=(ibx*nSubPixels)/nbeta;
+       } else {   
+	 //if (hiy[ibx]>tot_eta_y*(ii+1)/nSubPixels) ii++;
+	 hhy[ib+ibx*nbeta]=hiy[ibx]/tot_eta_y;
+       }
      }
-     /* if (ii!=(nSubPixels-1)) */
-     /*   cout << ib << " y " <<  tot_eta_y << " " << (ii+1)*tot_eta_y*bsize << " " << ii << " " << hiy[nbeta-1]<< endl; */
-     
-     //	 cout << "y " << nbeta << " " << (ii+1)*tot_eta_x*bsize << " " << ii << endl;
-     
    }
 
-#ifdef SAVE_ALL
-   char tit[10000];
-   
-  float *etah=new float[nbeta*nbeta];
-  int etabins=nbeta;
 
-  for (int ii=0; ii<etabins*etabins; ii++) {
-    etah[ii]=hhx[ii];
-  }
-  sprintf(tit,"/scratch/eta_hhx_%d.tiff",id);
-  WriteToTiff(etah, tit, etabins, etabins);
-	  
-  for (int ii=0; ii<etabins*etabins; ii++) {
-    etah[ii]=hhy[ii];
-  }
-  sprintf(tit,"/scratch/eta_hhy_%d.tiff",id);
-  WriteToTiff(etah, tit, etabins, etabins);
-	  
-  for (int ii=0; ii<etabins*etabins; ii++) {
-    etah[ii]=heta[ii];
-  }
-  sprintf(tit,"/scratch/eta_%d.tiff",id);
-  WriteToTiff(etah, tit, etabins, etabins);
-  delete [] etah;
+   int ibx, iby, ib; 
+   
+   iby=0;
+   while (hhx[iby*nbeta+nbeta/2]<0) iby++;
+   for (ib=0; ib<iby;ib++) {
+     for (ibx=0; ibx<nbeta;ibx++)
+       hhx[ibx+nbeta*ib]=hhx[ibx+nbeta*iby];
+   }
+   iby=nbeta-1;
+   
+   while (hhx[iby*nbeta+nbeta/2]<0) iby--;
+   for (ib=iby+1; ib<nbeta;ib++) {
+     for (ibx=0; ibx<nbeta;ibx++)
+       hhx[ibx+nbeta*ib]=hhx[ibx+nbeta*iby];
+   }
+   
+   iby=0;
+   while (hhy[nbeta/2*nbeta+iby]<0) iby++;
+   for (ib=0; ib<iby;ib++) {
+     for (ibx=0; ibx<nbeta;ibx++)
+       hhy[ib+nbeta*ibx]=hhy[iby+nbeta*ibx];
+   }
+   iby=nbeta-1;
+   
+   while (hhy[nbeta/2*nbeta+iby]<0) iby--;
+   for (ib=iby+1; ib<nbeta;ib++) {
+     for (ibx=0; ibx<nbeta;ibx++)
+       hhy[ib+nbeta*ibx]=hhy[iby+nbeta*ibx];
+   }
+
+
+
+     
+   
+
+
+#ifdef SAVE_ALL
+   debugSaveAll();
 #endif	  
   return ;
   }

slsDetectorCalibration/interpolations/slsInterpolation.h

@@ -1,12 +1,6 @@
 #ifndef SLS_INTERPOLATION_H
 #define SLS_INTERPOLATION_H
 
-#ifdef MYROOT1
-#include <TObject.h>
-#include <TTree.h>
-#include <TH2F.h>
-#endif
-
 #include <cstdlib> 
 #ifndef MY_TIFF_IO_H
 #include "tiffIO.h"
@@ -37,13 +31,7 @@ class slsInterpolation
  public:
  slsInterpolation(int nx=400, int ny=400, int ns=25) :nPixelsX(nx), nPixelsY(ny),  nSubPixels(ns), id(0) {
    
-#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
 
 };
   
@@ -51,14 +39,9 @@ hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny);
    nPixelsX=orig->nPixelsX;
    nPixelsY=orig->nPixelsY;
    nSubPixels=orig->nSubPixels;
-#ifdef MYROOT1
-   hint=(TH2F*)(orig->hint)->Clone("hint");
-#endif
    
-#ifndef MYROOT1
    hint=new int[nSubPixels*nPixelsX*nSubPixels*nPixelsY];
    memcpy(hint, orig->hint,nSubPixels*nPixelsX*nSubPixels*nPixelsY*sizeof(int));
-#endif
  
  };
 
@@ -94,11 +77,6 @@ hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny);
 
   //create interpolated image
   //returns interpolated image
-#ifdef MYROOT1
-  virtual TH2F *getInterpolatedImage(){return hint;};
-#endif
-
-#ifndef MYROOT1
   virtual int *getInterpolatedImage(){
     //  cout << "return interpolated image " << endl;
     /* for (int i=0; i<nSubPixels*  nSubPixels* nPixelsX*nPixelsY; i++) { */
@@ -106,7 +84,6 @@ hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny);
     /* } */
     return hint;
   };
-#endif
 
 
 
@@ -114,11 +91,12 @@ hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny);
   void *writeInterpolatedImage(const char * imgname) {
     //cout << "!" <<endl;
     float *gm=NULL;
+    int *dummy=getInterpolatedImage();
     gm=new float[ nSubPixels*  nSubPixels* nPixelsX*nPixelsY];
     if (gm) {
       for (int ix=0; ix<nPixelsX*nSubPixels; ix++) {
 	for (int iy=0; iy<nPixelsY*nSubPixels; iy++) {
-	  gm[iy*nPixelsX*nSubPixels+ix]=hint[iy*nPixelsX*nSubPixels+ix];
+	  gm[iy*nPixelsX*nSubPixels+ix]=dummy[iy*nPixelsX*nSubPixels+ix];
 	}
       }
       WriteToTiff(gm, imgname,nSubPixels* nPixelsX ,nSubPixels* nPixelsY); 
@@ -142,16 +120,11 @@ hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny);
   virtual void clearInterpolatedImage() {
 
 
-#ifdef MYROOT1
-    hint->Reset();
-#endif
-#ifndef MYROOT1
       for (int ix=0; ix<nPixelsX*nSubPixels; ix++) {
 	for (int iy=0; iy<nPixelsY*nSubPixels; iy++) {
 	  hint[iy*nPixelsX*nSubPixels+ix]=0;
 	}
       }
-#endif
 
 
   };
@@ -159,11 +132,6 @@ hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny);
 
 
 
-#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=((double)nSubPixels)*int_y; 
     int ix=((double)nSubPixels)*int_x; 
@@ -176,7 +144,6 @@ hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny);
        
     return hint;
   };
-#endif
 
 
   virtual int addToFlatField(double *cluster, double &etax, double &etay)=0;
@@ -185,19 +152,11 @@ hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny);
   virtual int  addToFlatField(double totquad,int quad,double *cluster,double &etax, double &etay)=0;
   virtual int addToFlatField(double etax, double etay)=0;
 
-#ifdef MYROOT1
-  virtual TH2D *getFlatField(){return NULL;};
-  virtual TH2D *setFlatField(TH2D *h, int nb=-1, double emin=-1, double emax=-1){return NULL;};
- virtual TH2D *getFlatField(int &nb, double &emin, double &emax){nb=0; emin=0; emax=0; return getFlatField();}; 
-#endif
-  
-#ifndef MYROOT1
   virtual int *getFlatField(){return NULL;};
   virtual int *setFlatField(int *h, int nb=-1, double emin=-1, double emax=-1){return NULL;}; 
   virtual void *writeFlatField(const char * imgname){return NULL;};
   virtual void *readFlatField(const char * imgname, int nb=-1, double emin=1, double emax=0){return NULL;};
  virtual int *getFlatField(int &nb, double &emin, double &emax){nb=0; emin=0; emax=0; return getFlatField();}; 
-#endif
 
  virtual void resetFlatField()=0;
 
@@ -215,10 +174,10 @@ hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny);
   static int calcQuad(double *cl, double &sum, double &totquad, double sDum[2][2]){
     
     int corner = UNDEFINED_QUADRANT;
-    double *cluster[3];
-    cluster[0]=cl;
-    cluster[1]=cl+3;
-    cluster[2]=cl+6;
+    /* double *cluster[3]; */
+    /* cluster[0]=cl; */
+    /* cluster[1]=cl+3; */
+    /* cluster[2]=cl+6; */
     
     sum=0;
     double sumBL=0;
@@ -228,11 +187,11 @@ hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny);
     int xoff=0, yoff=0;
     for (int ix=0; ix<3; ix++) {
       for (int iy=0; iy<3; iy++) {
-	sum+=cluster[iy][ix];
-	 if (ix<=1 && iy<=1) sumBL+=cluster[iy][ix];
-	 if (ix<=1 && iy>=1) sumTL+=cluster[iy][ix];
-	 if (ix>=1 && iy<=1) sumBR+=cluster[iy][ix];
-	 if (ix>=1 && iy>=1) sumTR+=cluster[iy][ix];
+	sum+=cl[ix+3*iy];
+	 if (ix<=1 && iy<=1) sumBL+=cl[ix+iy*3];
+	 if (ix<=1 && iy>=1) sumTL+=cl[ix+iy*3];
+	 if (ix>=1 && iy<=1) sumBR+=cl[ix+iy*3];
+	 if (ix>=1 && iy>=1) sumTR+=cl[ix+iy*3];
       }
     }
 
@@ -240,6 +199,8 @@ hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny);
     /* sDum[0][1] = cluster[0][1]; sDum[1][1] = cluster[1][1]; */
     corner = BOTTOM_LEFT;
     totquad=sumBL;
+    xoff=0;
+    yoff=0;
 
     
     if(sumTL  >= totquad){
@@ -274,7 +235,7 @@ hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny);
 
     for (int ix=0; ix<2; ix++) {
       for (int iy=0; iy<2; iy++) {
-	sDum[iy][ix] = cluster[iy+yoff][ix+xoff]; 
+	sDum[iy][ix] = cl[ix+xoff+(iy+yoff)*3]; 
       }
     }
   
@@ -396,7 +357,7 @@ hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny);
     //  int quad;
     for (int ix=0; ix<3; ix++) {
       for (int iy=0; iy<3; iy++) {
-	val=cl[iy+3*ix];
+	val=cl[ix+3*iy];
 	sum+=val;
 	if (iy==0) l+=val;
 	if (iy==2) r+=val;
@@ -440,92 +401,92 @@ hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny);
   }
 
 
-  static int calcMyEta(double totquad, int quad, double *cl, double &etax, double &etay) {
-    double l,r,t,b, sum;
-    int yoff;
-    switch (quad) {
-     case BOTTOM_LEFT:
-     case BOTTOM_RIGHT:
-       yoff=0;
-       break;
-     case TOP_LEFT:
-     case TOP_RIGHT:
-       yoff=1;
-       break;
-     default:
-       ;
-     } 
-      l=cl[0+yoff*3]+cl[0+yoff*3+3];
-      r=cl[2+yoff*3]+cl[2+yoff*3+3];
-      b=cl[0+yoff*3]+cl[1+yoff*3]*cl[2+yoff*3];
-      t=cl[0+yoff*3+3]+cl[1+yoff*3+3]*cl[0+yoff*3+3];
-      sum=t+b;
-    if (sum>0) {
-      etax=(-l+r)/sum;
-      etay=(+t)/sum;
-    }
+  /* static int calcMyEta(double totquad, int quad, double *cl, double &etax, double &etay) { */
+  /*   double l,r,t,b, sum; */
+  /*   int yoff; */
+  /*   switch (quad) { */
+  /*    case BOTTOM_LEFT: */
+  /*    case BOTTOM_RIGHT: */
+  /*      yoff=0; */
+  /*      break; */
+  /*    case TOP_LEFT: */
+  /*    case TOP_RIGHT: */
+  /*      yoff=1; */
+  /*      break; */
+  /*    default: */
+  /*      ; */
+  /*    }  */
+  /*     l=cl[0+yoff*3]+cl[0+yoff*3+3]; */
+  /*     r=cl[2+yoff*3]+cl[2+yoff*3+3]; */
+  /*     b=cl[0+yoff*3]+cl[1+yoff*3]*cl[2+yoff*3]; */
+  /*     t=cl[0+yoff*3+3]+cl[1+yoff*3+3]*cl[0+yoff*3+3]; */
+  /*     sum=t+b; */
+  /*   if (sum>0) { */
+  /*     etax=(-l+r)/sum; */
+  /*     etay=(+t)/sum; */
+  /*   } */
     
-    return -1;
-  }
+  /*   return -1; */
+  /* } */
 
-  static int calcMyEta(double totquad, int quad, int *cl, double &etax, double &etay) {
-    double l,r,t,b, sum;
-    int yoff;
-    switch (quad) {
-     case BOTTOM_LEFT:
-     case BOTTOM_RIGHT:
-       yoff=0;
-       break;
-     case TOP_LEFT:
-     case TOP_RIGHT:
-       yoff=1;
-       break;
-     default:
-       ;
-     } 
-      l=cl[0+yoff*3]+cl[0+yoff*3+3];
-      r=cl[2+yoff*3]+cl[2+yoff*3+3];
-      b=cl[0+yoff*3]+cl[1+yoff*3]*cl[2+yoff*3];
-      t=cl[0+yoff*3+3]+cl[1+yoff*3+3]*cl[0+yoff*3+3];
-      sum=t+b;
-    if (sum>0) {
-      etax=(-l+r)/sum;
-      etay=(+t)/sum;
-    }
+  /* static int calcMyEta(double totquad, int quad, int *cl, double &etax, double &etay) { */
+  /*   double l,r,t,b, sum; */
+  /*   int yoff; */
+  /*   switch (quad) { */
+  /*    case BOTTOM_LEFT: */
+  /*    case BOTTOM_RIGHT: */
+  /*      yoff=0; */
+  /*      break; */
+  /*    case TOP_LEFT: */
+  /*    case TOP_RIGHT: */
+  /*      yoff=1; */
+  /*      break; */
+  /*    default: */
+  /*      ; */
+  /*    }  */
+  /*     l=cl[0+yoff*3]+cl[0+yoff*3+3]; */
+  /*     r=cl[2+yoff*3]+cl[2+yoff*3+3]; */
+  /*     b=cl[0+yoff*3]+cl[1+yoff*3]*cl[2+yoff*3]; */
+  /*     t=cl[0+yoff*3+3]+cl[1+yoff*3+3]*cl[0+yoff*3+3]; */
+  /*     sum=t+b; */
+  /*   if (sum>0) { */
+  /*     etax=(-l+r)/sum; */
+  /*     etay=(+t)/sum; */
+  /*   } */
     
-    return -1;
-  }
+  /*   return -1; */
+  /* } */
 
 
 
-  static int calcEta3X(double *cl, double &etax, double &etay, double &sum) {
-    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];
-      r=cl[5];
-      b=cl[1];
-      t=cl[7];
-      etax=(-l+r)/sum;
-      etay=(-b+t)/sum;
-    }
-    return -1;
-  }
+  /* static int calcEta3X(double *cl, double &etax, double &etay, double &sum) { */
+  /*   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]; */
+  /*     r=cl[5]; */
+  /*     b=cl[1]; */
+  /*     t=cl[7]; */
+  /*     etax=(-l+r)/sum; */
+  /*     etay=(-b+t)/sum; */
+  /*   } */
+  /*   return -1; */
+  /* } */
 
 
-  static int calcEta3X(int *cl, double &etax, double &etay, double &sum) {
-    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];
-      r=cl[5];
-      b=cl[1];
-      t=cl[7];
-      etax=(-l+r)/sum;
-      etay=(-b+t)/sum;
-    }
-    return -1;
-  }
+  /* static int calcEta3X(int *cl, double &etax, double &etay, double &sum) { */
+  /*   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]; */
+  /*     r=cl[5]; */
+  /*     b=cl[1]; */
+  /*     t=cl[7]; */
+  /*     etax=(-l+r)/sum; */
+  /*     etay=(-b+t)/sum; */
+  /*   } */
+  /*   return -1; */
+  /* } */
 
 
 
@@ -534,15 +495,9 @@ hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny);
 
  protected:
   int nPixelsX, nPixelsY;
-  int nSubPixels; 
-#ifdef MYROOT1
-  TH2F *hint;
-#endif
-#ifndef MYROOT1
-  int *hint;
-#endif
+  int nSubPixels;
   int id;
-
+  int *hint;
 };
 
 #endif

slsDetectorCalibration/moenchExecutables/Makefile.moench_zmq_rh7

@@ -1,7 +1,7 @@
 
 ZMQLIB=../../slsReceiverSoftware/include
 INCDIR= -I$(ZMQLIB)  -I. -I../dataStructures ../tiffIO.cpp -I../ -I../interpolations/
-LDFLAG= -L/usr/lib64/ -lpthread -lm  -lstdc++  -lzmq -pthread -lrt -ltiff -L$(ZMQLIB)  -O3 
+LDFLAG= -L/usr/lib64/ -lpthread -lm  -lstdc++  -lzmq -pthread -lrt -ltiff -L$(ZMQLIB)  -O3 -g
 #-L../../bin -lhdf5  -L.
 
 #DESTDIR?=../bin

slsDetectorCalibration/moenchExecutables/moench03ClusterFinder.cpp

@@ -48,7 +48,7 @@ int main(int argc, char *argv[]) {
   }
   int p=10000;
   int fifosize=1000;
-  int nthreads=1;
+  int nthreads=24;
   int nsubpix=25;
   int etabins=nsubpix*10;
   double etamin=-1, etamax=2;

slsDetectorCalibration/moenchExecutables/moench03Interpolation.cpp

@@ -14,14 +14,14 @@
 #include "single_photon_hit.h"
 #endif
 
-#include "etaInterpolationPosXY.h"
+//#include "etaInterpolationPosXY.h"
 #include "noInterpolation.h"
-//#include "etaInterpolationAdaptiveBins.h"
+#include "etaInterpolationCleverAdaptiveBins.h"
 //#include "etaInterpolationRandomBins.h"
 using namespace std;
 #define NC 400
 #define NR 400
-
+#define MAX_ITERATIONS (nSubPixels*100)
 
 #define XTALK
 
@@ -65,9 +65,10 @@ int main(int argc, char *argv[]) {
   float cmax=atof(argv[iarg++]);
   cout << "Energy min: " << cmin << endl;
   cout << "Energy max: " << cmax << endl;
-  
+  //int etabins=500;
   int etabins=1000;//nsubpix*2*100;
   double etamin=-1, etamax=2;
+  //double etamin=-0.1, etamax=1.1;
   double eta3min=-2, eta3max=2;
   int quad;
   double sum, totquad;
@@ -90,15 +91,25 @@ int main(int argc, char *argv[]) {
     single_photon_hit cl(3,3);
 #endif
 
+    int nSubPixels=nsubpix;
+#ifndef NOINTERPOLATION
+   eta2InterpolationPosXY *interp=new eta2InterpolationPosXY(NC, NR, nsubpix, etabins, etamin, etamax);
+   //eta2InterpolationCleverAdaptiveBins *interp=new eta2InterpolationCleverAdaptiveBins(NC, NR, nsubpix, etabins, etamin, etamax);  
+#endif
+#ifdef NOINTERPOLATION
+    noInterpolation *interp=new noInterpolation(NC, NR, nsubpix);  
+#endif
 
-    eta2InterpolationPosXY *interp=new eta2InterpolationPosXY(NC, NR, nsubpix, etabins, etamin, etamax);
-   
+
+ 
 #ifndef FF
+#ifndef NOINTERPOLATION
     cout << "read ff " << argv[2] << endl;
     sprintf(fname,"%s",argv[2]);
     interp->readFlatField(fname);
-    interp->prepareInterpolation(ok);
-
+    interp->prepareInterpolation(ok);//, MAX_ITERATIONS);
+#endif
+    // return 0;
 #endif
 #ifdef FF
     cout << "Will write eta file " << argv[2] << endl;
@@ -142,27 +153,49 @@ int main(int argc, char *argv[]) {
 	    if (nframes==0) f0=lastframe;
 	    nframes++;
 	  }
-	  quad=interp->calcQuad(cl.get_cluster(), sum, totquad, sDum);
-	  if (sum>cmin && totquad/sum>0.8 && totquad/sum<1.2 && sum<cmax ) {
+	  //quad=interp->calcQuad(cl.get_cluster(), sum, totquad, sDum);
+	  quad=interp->calcEta(cl.get_cluster(), etax, etay, sum, totquad, sDum);
+	    if (sum>cmin && totquad/sum>0.8 && totquad/sum<1.2 && sum<cmax ) {
 	    nph++;
 	    //  if (sum>200 && sum<580) {
 	    //  interp->getInterpolatedPosition(cl.x,cl.y, totquad,quad,cl.get_cluster(),int_x, int_y);
-#ifdef SOLEIL
-	    if (cl.x>210 && cl.x<240 && cl.y>210 && cl.y<240) {
-#endif
+// #ifdef SOLEIL
+// 	    if (cl.x>210 && cl.x<240 && cl.y>210 && cl.y<240) {
+// #endif
 #ifndef FF
-	    interp->getInterpolatedPosition(cl.x,cl.y, cl.get_cluster(),int_x, int_y);
+	    // interp->getInterpolatedPosition(cl.x,cl.y, cl.get_cluster(),int_x, int_y);  
+	    interp->getInterpolatedPosition(cl.x,cl.y, etax, etay, quad,int_x, int_y);
+	      // cout <<"**************"<< endl;
+	      // cout << cl.x << " " << cl.y << " " << sum << endl;
+	      // cl.print();
+	      // cout << int_x << " " << int_y << endl;
+	      // cout <<"**************"<< endl;
+	    //  if (etax!=0 && etay!=0 && etax!=1 && etay!=1)
 	    interp->addToImage(int_x, int_y);
-#endif
-#ifdef FF
-	    interp->addToFlatField(cl.get_cluster(), etax, etay);
-#endif
-#ifdef SOLEIL
+	    if (int_x<0 || int_y<0 || int_x>400 || int_y>400) {
+	      cout <<"**************"<< endl;
+	      cout << cl.x << " " << cl.y << " " << sum << endl;
+	      cl.print();
+	      cout << int_x << " " << int_y << endl;
+	      cout <<"**************"<< endl;
 	    }
 #endif
+#ifdef FF
+	    //	interp->addToFlatField(cl.get_cluster(), etax, etay);
+// #ifdef UCL
+// 	    if (cl.x>50)
+// #endif
+// 	    if (etax!=0 && etay!=0 && etax!=1 && etay!=1)
+	      interp->addToFlatField(etax, etay);
+	    //  if (etax==0 || etay==0) cout << cl.x << " " << cl.y << endl;
+	     
+#endif
+// #ifdef SOLEIL
+// 	    }
+// #endif
 	    
 	    if (nph%1000000==0) cout << nph << endl;
-	    if (nph%100000000==0) { 
+	    if (nph%10000000==0) { 
 #ifndef FF
 		interp->writeInterpolatedImage(outfname);
 #endif
@@ -170,12 +203,12 @@ int main(int argc, char *argv[]) {
 		interp->writeFlatField(outfname);
 #endif
 		
-	    }
+		}
+	  }
 	    
-	  } 
+	} 
 	  
 	  
-	}
 	fclose(f);
 #ifdef FF
 	interp->writeFlatField(outfname);
@@ -194,7 +227,7 @@ int main(int argc, char *argv[]) {
 	    }
 	  }
 	}
-	cout << "Read " << nframes << " frames (first frame: " << f0 << " last frame: " << lastframe << " delta:" << lastframe-f0 << ")"<<endl;
+	cout << "Read " << nframes << " frames (first frame: " << f0 << " last frame: " << lastframe << " delta:" << lastframe-f0 << ") nph="<< nph <<endl;
 	interp->clearInterpolatedImage();
 #endif
 	

slsDetectorCalibration/moenchExecutables/moench03InterpolationColor.cpp

@@ -0,0 +1,269 @@
+
+#include "ansi.h"
+#include <iostream>
+
+//#include "moench03T1ZmqData.h"
+//#define DOUBLE_SPH
+//#define MANYFILES
+
+#ifdef DOUBLE_SPH
+#include "single_photon_hit_double.h"
+#endif
+
+#ifndef DOUBLE_SPH
+#include "single_photon_hit.h"
+#endif
+
+//#include "etaInterpolationPosXY.h"
+#include "noInterpolation.h"
+#include "etaInterpolationCleverAdaptiveBins.h"
+//#include "etaInterpolationRandomBins.h"
+using namespace std;
+#define NC 400
+#define NR 400
+#define MAX_ITERATIONS (nSubPixels*100)
+#define MAX_EBINS 100
+#define XTALK
+
+int main(int argc, char *argv[]) {
+
+#ifndef FF
+  if (argc<9) {
+    cout << "Wrong usage! Should be: "<< argv[0] << " infile  etafile outfile runmin runmax ns cmin cmax" << endl;
+    return 1;
+  }
+#endif
+  
+#ifdef FF
+  if (argc<7) {
+    cout << "Wrong usage! Should be: "<< argv[0] << " infile  etafile runmin runmax cmin cmax" << endl;
+    return 1;
+  }
+#endif
+  int iarg=4;
+  char infname[10000];
+  char fname[10000];
+  char outfname[10000];
+#ifndef FF
+  iarg=4;
+#endif
+
+#ifdef FF
+  iarg=3;
+#endif
+  int runmin=atoi(argv[iarg++]);
+  int runmax=atoi(argv[iarg++]); 
+  cout << "Run min: " << runmin << endl;
+  cout << "Run max: " << runmax << endl;
+  
+  int nsubpix=4;
+#ifndef FF
+  nsubpix=atoi(argv[iarg++]);
+  cout << "Subpix: " << nsubpix << endl;
+#endif
+  float cmin=atof(argv[iarg++]);
+  float cmax=atof(argv[iarg++]);
+  cout << "Energy min: " << cmin << endl;
+  cout << "Energy max: " << cmax << endl;
+  int n_ebins=1;
+  if (argc>iarg)
+    n_ebins=atoi(argv[iarg++]);
+  //int etabins=500;
+  int etabins=1000;//nsubpix*2*100;
+  double etamin=-1, etamax=2;
+  //double etamin=-0.1, etamax=1.1;
+  double eta3min=-2, eta3max=2;
+  int quad;
+  double sum, totquad;
+  double sDum[2][2];
+  double etax, etay, int_x, int_y;
+  double eta3x, eta3y, int3_x, int3_y, noint_x, noint_y;
+  int ok;
+  int f0=-1;
+  int ix, iy, isx, isy;
+  int nframes=0, lastframe=-1;
+  double d_x, d_y, res=5, xx, yy;
+  int nph=0, badph=0, totph=0;
+    FILE *f=NULL;
+
+#ifdef DOUBLE_SPH
+    single_photon_hit_double cl(3,3);
+#endif
+
+#ifndef DOUBLE_SPH
+    single_photon_hit cl(3,3);
+#endif
+
+    int nSubPixels=nsubpix;
+
+    int iebin=0;
+    double eb_size=(cmax-cmin)/n_ebins;
+#ifndef NOINTERPOLATION
+    //  eta2InterpolationPosXY *interp[MAX_EBINS];
+    eta2InterpolationCleverAdaptiveBins *interp[MAX_EBINS];
+    for (int i=0; i< n_ebins; i++) {
+      //interp[i]=new eta2InterpolationPosXY(NC, NR, nsubpix, etabins, etamin, etamax);
+      interp[i]=new eta2InterpolationCleverAdaptiveBins(NC, NR, nsubpix, etabins, etamin, etamax);  
+    }
+#endif
+#ifdef NOINTERPOLATION
+    noInterpolation *interp=new noInterpolation(NC, NR, nsubpix);  
+#endif
+
+
+ 
+#ifndef FF
+#ifndef NOINTERPOLATION
+    cout << "read ff " << argv[2] << endl; 
+    for (int i=0; i< n_ebins; i++) {
+      sprintf(fname,argv[2],i);
+      interp[i]->readFlatField(fname);
+      interp[i]->prepareInterpolation(ok);//, MAX_ITERATIONS);
+    }
+#endif
+    // return 0;
+#endif
+#ifdef FF
+    cout << "Will write eta file " << argv[2] << endl;
+#endif
+
+    int *img;
+    float *totimg=new float[NC*NR*nsubpix*nsubpix];
+    for (ix=0; ix<NC; ix++) {
+      for (iy=0; iy<NR; iy++) {
+	for (isx=0; isx<nsubpix; isx++) {
+	  for (isy=0; isy<nsubpix; isy++) {
+	    totimg[ix*nsubpix+isx+(iy*nsubpix+isy)*(NC*nsubpix)]=0;
+	      }
+	}
+      }
+    }
+
+#ifdef FF
+      sprintf(outfname,argv[2]);
+#endif     
+
+    int irun;
+    for (irun=runmin; irun<runmax; irun++) {
+      sprintf(infname,argv[1],irun);
+#ifndef FF
+      sprintf(outfname,argv[3],irun);
+#endif
+ 
+      f=fopen(infname,"r");
+      if (f) {
+	cout << infname << endl;
+	nframes=0;
+	f0=-1;
+
+	while (cl.read(f)) {
+	  totph++;
+	  if (lastframe!=cl.iframe) {
+	    lastframe=cl.iframe;
+	    // cout << cl.iframe << endl;
+	    // f0=cl.iframe;
+	    if (nframes==0) f0=lastframe;
+	    nframes++;
+	  }
+	  //quad=interp->calcQuad(cl.get_cluster(), sum, totquad, sDum);
+	  quad=interp[0]->calcEta(cl.get_cluster(), etax, etay, sum, totquad, sDum);
+	  if (sum>cmin && totquad/sum>0.8 && totquad/sum<1.2 && sum<cmax ) {
+	    nph++;
+	    iebin=(sum-cmin)/eb_size;
+	    if (iebin>=0 && iebin<n_ebins) {
+	    //  if (sum>200 && sum<580) {
+	    //  interp->getInterpolatedPosition(cl.x,cl.y, totquad,quad,cl.get_cluster(),int_x, int_y);
+// #ifdef SOLEIL
+// 	    if (cl.x>210 && cl.x<240 && cl.y>210 && cl.y<240) {
+// #endif
+#ifndef FF
+	    // interp->getInterpolatedPosition(cl.x,cl.y, cl.get_cluster(),int_x, int_y);  
+	    interp[iebin]->getInterpolatedPosition(cl.x,cl.y, etax, etay, quad,int_x, int_y);
+	      // cout <<"**************"<< endl;
+	      // cout << cl.x << " " << cl.y << " " << sum << endl;
+	      // cl.print();
+	      // cout << int_x << " " << int_y << endl;
+	      // cout <<"**************"<< endl;
+	    if (etax!=0 && etay!=0 && etax!=1 && etay!=1)
+	      interp[iebin]->addToImage(int_x, int_y);
+#endif
+#ifdef FF
+	    //	interp->addToFlatField(cl.get_cluster(), etax, etay);
+#ifdef UCL
+	    if (cl.x>50)
+#endif
+	      if (etax!=0 && etay!=0 && etax!=1 && etay!=1)
+		interp[iebin]->addToFlatField(etax, etay);
+	    //  if (etax==0 || etay==0) cout << cl.x << " " << cl.y << endl;
+	     
+#endif
+// #ifdef SOLEIL
+// 	    }
+// #endif
+	    
+	    if (nph%1000000==0) cout << nph << endl;
+	    if (nph%10000000==0) { 
+#ifndef FF
+	      for (int i=0; i<n_ebins; i++) {
+		sprintf(outfname,argv[3],i);
+		interp[i]->writeInterpolatedImage(outfname);
+	      }
+#endif
+#ifdef FF 
+	      for (int i=0; i<n_ebins; i++) {
+		sprintf(outfname,argv[2],i);
+		cout << outfname << " " << argv[2] << " " << i << endl;
+		interp[i]->writeFlatField(outfname);
+	      }
+#endif
+	      
+	    }
+	    }
+	    
+	  } 
+	  
+	}
+	fclose(f);
+#ifdef FF
+	for (int i=0; i<n_ebins; i++) {
+	  sprintf(outfname,argv[2],i);
+	  cout << outfname << " " << argv[2] << " " << i << endl;
+	  interp[i]->writeFlatField(outfname);
+	}
+#endif
+	
+#ifndef FF 
+	      for (int i=0; i<n_ebins; i++) {
+		sprintf(outfname,argv[3],i,irun);
+		interp[i]->writeInterpolatedImage(outfname);
+		img=interp[i]->getInterpolatedImage();
+		for (ix=0; ix<NC; ix++) {
+		  for (iy=0; iy<NR; iy++) {
+		    for (isx=0; isx<nsubpix; isx++) {
+		      for (isy=0; isy<nsubpix; isy++) {
+			totimg[ix*nsubpix+isx+(iy*nsubpix+isy)*(NC*nsubpix)]+=img[ix*nsubpix+isx+(iy*nsubpix+isy)*(NC*nsubpix)];
+		      }
+		    }
+		  }
+		}
+		//interp[i]->clearInterpolatedImage();
+	      }
+	      cout << "Read " << nframes << " frames (first frame: " << f0 << " last frame: " << lastframe << " delta:" << lastframe-f0 << ")"<<endl;
+#endif
+	
+      } else
+	cout << "could not open file " << infname << endl;
+    }
+#ifndef FF
+    sprintf(outfname,argv[3], 11111);
+    WriteToTiff(totimg, outfname,NC*nsubpix,NR*nsubpix);
+#endif 
+    
+#ifdef FF
+	interp[iebin]->writeFlatField(outfname);
+#endif
+
+    cout << "Filled " << nph << " (/"<< totph <<") " << endl; 
+    return 0;
+}
+

slsDetectorCalibration/moenchExecutables/moenchZmqProcess.cpp

@@ -1,3 +1,5 @@
+#define WRITE_QUAD
+
 #include "sls_receiver_defs.h"
 #include "ZmqSocket.h"
 #include "moench03T1ZmqDataNew.h"
@@ -8,7 +10,6 @@
 #include <fstream> 
 #include "tiffIO.h"
 
-
 //#define NEWZMQ
 #ifdef NEWZMQ
 #include <rapidjson/document.h> //json header in zmq stream
@@ -25,10 +26,15 @@
 #include "etaInterpolationPosXY.h"
 #include "ansi.h"
 #include <iostream>
+
+//#include <chrono>
+#include <ctime> // time_t
+#include <cstdio>
+
 using namespace std;
+//using namespace std::chrono;
 
-
-#define SLS_DETECTOR_JSON_HEADER_VERSION 0x2
+//#define SLS_DETECTOR_JSON_HEADER_VERSION 0x2
 
 	//	myDet->setNetworkParameter(ADDITIONAL_JSON_HEADER, " \"what\":\"nothing\" ");
 
@@ -38,7 +44,7 @@ int main(int argc, char *argv[]) {
  *
  */
   FILE *of=NULL;
-  int fifosize=1000;
+  int fifosize=5000;
   int etabins=1000;//nsubpix*2*100;
   double etamin=-1, etamax=2;
 	// help
@@ -58,8 +64,9 @@ int main(int argc, char *argv[]) {
 
   int ok;
 
-
-
+  // high_resolution_clock::time_point t1;
+  // high_resolution_clock::time_point t2 ;
+  time_t begin,end,finished;
 
 
   if (argc > 4) {
@@ -105,6 +112,7 @@ int main(int argc, char *argv[]) {
   int multisize=size;
   int dataSize=size;
 
+  char dummybuff[size];
 
 
 
@@ -283,13 +291,13 @@ int main(int argc, char *argv[]) {
 #ifdef NEWZMQ	
 	    rapidjson::Document doc;
 	    if (!zmqsocket->ReceiveHeader(0, doc, SLS_DETECTOR_JSON_HEADER_VERSION)) {
-	      zmqsocket->CloseHeaderMessage();
+	      /* zmqsocket->CloseHeaderMessage();*/
 
 #endif
 	    //	  if (!zmqsocket->ReceiveHeader(0, acqIndex, frameIndex, subframeIndex, filename, fileindex)) {
-	    cprintf(RED, "Got Dummy\n");
-
-
+	      //  cprintf(RED, "Got Dummy\n");
+	      //   t1=high_resolution_clock::now();
+	      time(&end);
 
 
 	    while (mt->isBusy()) {;}//wait until all data are processed from the queues
@@ -413,6 +421,14 @@ int main(int argc, char *argv[]) {
 	    mt->clearImage();
 	    
 	    newFrame=1;
+	    //t2 = high_resolution_clock::now();
+	    
+	      time(&finished);
+	      //   auto meas_duration = duration_cast<microseconds>( t2 - t0 ).count();
+	      //   auto real_duration = duration_cast<microseconds>( t2 - t1 ).count();
+	    
+	    cout << "Measurement lasted " << difftime(end,begin) << endl;
+	    cout << "Processing lasted " << difftime(finished,begin) << endl;
 	    continue; //continue to not get out
 
 
@@ -420,31 +436,29 @@ int main(int argc, char *argv[]) {
 
 #ifdef NEWZMQ
 	    if (newFrame) {
-
+	      time(&begin);
+	      // t0 = high_resolution_clock::now();
+	      //cout <<"new frame" << endl;
 
 	      //  acqIndex, frameIndex, subframeIndex, filename, fileindex
 	      size       = doc["size"].GetUint();
-	      multisize  = size;// * zmqsocket->size();
+	      // multisize  = size;// * zmqsocket->size();
 	      dynamicRange  = doc["bitmode"].GetUint();
-	      nPixelsX = doc["shape"][0].GetUint();
-	      nPixelsY = doc["shape"][1].GetUint();
+	      //  nPixelsX = doc["shape"][0].GetUint();
+	      // nPixelsY = doc["shape"][1].GetUint();
 	      filename        = doc["fname"].GetString();
-	      acqIndex = doc["acqIndex"].GetUint64();
-	      frameIndex       = doc["fIndex"].GetUint64();
+	      //acqIndex = doc["acqIndex"].GetUint64();
+	      //frameIndex       = doc["fIndex"].GetUint64();
 	      fileindex        = doc["fileIndex"].GetUint64();
-	      subFrameIndex    = doc["expLength"].GetUint();
-	      xCoord                  = doc["xCoord"].GetUint();
-	      yCoord                  = doc["yCoord"].GetUint();
-	      zCoord                  = doc["zCoord"].GetUint();
-	      flippedDataX=doc["flippedDataX"].GetUint();
-	      packetNumber=doc["packetNumber"].GetUint();
-	      bunchId=doc["bunchId"].GetUint();
-	      timestamp=doc["timestamp"].GetUint();
-	      modId=doc["modId"].GetUint();
-	      debug=doc["debug"].GetUint();
-	      roundRNumber=doc["roundRNumber"].GetUint();
-	      detType=doc["detType"].GetUint();
-	      version=doc["version"].GetUint();
+	      //subFrameIndex    = doc["expLength"].GetUint();
+	      //packetNumber=doc["packetNumber"].GetUint();
+	      //bunchId=doc["bunchId"].GetUint();
+	      //timestamp=doc["timestamp"].GetUint();
+	      //modId=doc["modId"].GetUint();
+	      //debug=doc["debug"].GetUint();
+	      //roundRNumber=doc["roundRNumber"].GetUint();
+	      //detType=doc["detType"].GetUint();
+	      //version=doc["version"].GetUint();
 
 	      dataSize=size;
 
@@ -636,7 +650,7 @@ int main(int argc, char *argv[]) {
 	      
 	      
 	      newFrame=0;
-	      zmqsocket->CloseHeaderMessage();
+	      /* zmqsocket->CloseHeaderMessage();*/
 	    }
 #endif
 
@@ -656,11 +670,26 @@ int main(int argc, char *argv[]) {
 
 	  // cout << "data" << endl;
 	  // get data
-	  length = zmqsocket->ReceiveData(0, buff, size);
-	  mt->pushData(buff);
-	  mt->nextThread();
-	  mt->popFree(buff);
-	  
+	  // acqIndex = doc["acqIndex"].GetUint64();
+	  frameIndex       = doc["fIndex"].GetUint64();
+	  // subFrameIndex    = doc["expLength"].GetUint();
+
+	  //  bunchId=doc["bunchId"].GetUint();
+	  //  timestamp=doc["timestamp"].GetUint();
+	      packetNumber=doc["packetNumber"].GetUint();
+	      // cout << acqIndex << " " << frameIndex << " " << subFrameIndex << " "<< bunchId << " " << timestamp << " " << packetNumber << endl;
+	      if (packetNumber>=40) {
+		//*((int*)buff)=frameIndex;
+		memcpy(buff,&frameIndex,sizeof(int));
+		length = zmqsocket->ReceiveData(0, buff+sizeof(int), size);
+		mt->pushData(buff);
+		mt->nextThread();
+		mt->popFree(buff);
+	      } else {
+		cprintf(RED, "Incomplete frame: received only %d packet\n", packetNumber);
+		length = zmqsocket->ReceiveData(0, dummybuff, size);
+
+	      }
 	
 	  
 	  

slsDetectorCalibration/multiThreadedAnalogDetector.h

@@ -27,6 +27,7 @@
 using namespace std;
 
 
+
 class threadedAnalogDetector
 {
 public:
@@ -35,12 +36,25 @@ public:
     det=d;
     fifoFree=new CircularFifo<char>(fs);
     fifoData=new CircularFifo<char>(fs);
-    mem=(char*)malloc(fs*det->getDataSize());
-    // cout << "data size is " << det->getDataSize()*fs << endl;
-    for (int i=0; i<fs; i++) {
-      mm=mem+i*det->getDataSize();
-      fifoFree->push(mm);
+   
+    /* mem=(char*)calloc(fs, det->getDataSize()); */
+    /* if (mem) */
+    /*   memset(mem,0, fs*det->getDataSize()); */
+    int i;
+    for (i=0; i<fs; i++) { 
+      //
+      //  mm=mem+i*det->getDataSize();
+      // cout << i << endl;
+      mm=(char*)calloc(1, det->getDataSize());
+      if (mm) { 
+	//memset(mm,0, det->getDataSize());
+	fifoFree->push(mm);
+      } else
+	break;
     }  
+    if (i<fs) cout << "Could allocate only "<< i <<" frames";
+    
+
     busy=0;
     stop=1;
     fMode=eFrame;
@@ -264,6 +278,7 @@ public:
     }
     
     for (int i=0; i<nThreads; i++) {
+      cout << "**" << i << endl;
       dets[i]=new threadedAnalogDetector(dd[i], fs);
     }
     

slsDetectorCalibration/singlePhotonDetector.h

@@ -218,8 +218,8 @@ public analogDetector<uint16_t> {
 	    cout << "add to common mode?"<< endl;
 	    addToCommonMode(data);
 	  }
-	  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)) {
 		val=subtractPedestal(data,ix,iy, cm);
 		
@@ -236,8 +236,8 @@ public analogDetector<uint16_t> {
 	    }
 	  }
 	  
-	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)) {
 	      eventMask[iy][ix]=PEDESTAL;
@@ -276,7 +276,6 @@ public analogDetector<uint16_t> {
 		      //}
 		    }
 		  }
-		}
 		
 		if (rest[iy][ix]>=max) { 
 		  if (bl>=br && bl>=tl && bl>=tr) {
@@ -324,6 +323,7 @@ public analogDetector<uint16_t> {
 	      }
 	    }
 	  }
+		}
 	}
 	} else return getClusters(data, nph);
       }
@@ -331,113 +331,6 @@ public analogDetector<uint16_t> {
     };
 
 
-  /*   /\** finds event type for pixel and fills cluster structure. The algorithm loops only if the evenMask for this pixel is still undefined. */
-  /* 	if pixel or cluster around it are above threshold (nsigma*pedestalRMS) cluster is filled and pixel mask is PHOTON_MAX (if maximum in cluster) or NEIGHBOUR; If PHOTON_MAX, the elements of the cluster are also set as NEIGHBOURs in order to speed up the looping */
-  /* 	if below threshold the pixel is either marked as PEDESTAL (and added to the pedestal calculator) or NEGATIVE_PEDESTAL is case it's lower than -threshold, otherwise the pedestal average would drift to negative values while it should be 0. */
-
-  /* 	/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). */
-  /* 	/returns event type for the given pixel */
-  /*   *\/ */
-  /*   eventType getEventType(char *data, int ix, int iy, int cm=0) { */
-
-  /*     // eventType ret=PEDESTAL; */
-  /*     double max=0, tl=0, tr=0, bl=0,br=0, v; */
-  /*     //  cout << iframe << endl; */
-    
-  /*     int cy=(clusterSizeY+1)/2; */
-  /*     int cs=(clusterSize+1)/2; */
-  /*     double val; */
-  /*     tot=0; */
-  /*     quadTot=0; */
-  /*     quad=UNDEFINED_QUADRANT; */
-
-  /*     if (iframe<nDark) { */
-  /* 	addToPedestal(data, ix,iy); */
-  /* 	return UNDEFINED_EVENT; */
-  /*     } */
-
-
-     
-      
-
-  /*     //   if (eventMask[iy][ix]==UNDEFINED) { */
-	
-  /*     eventMask[iy][ix]=PEDESTAL; */
-	
-	
-  /*     clusters->x=ix; */
-  /*     clusters->y=iy; */
-  /*     clusters->rms=getPedestalRMS(ix,iy); */
-  /*     clusters->ped=getPedestal(ix,iy, cm); */
-	
-
-  /*     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) { */
-	      
-  /* 	      v=subtractPedestal(data, ix+ic, iy+ir); */
-	      
-  /* 	      clusters->set_data(v, ic, ir); */
-  /* 	      //  v=clusters->get_data(ic,ir); */
-  /* 	      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*clusters->rms) */
-  /* 		  eventMask[iy][ix]=NEGATIVE_PEDESTAL; */
-  /* 	      } */
-  /* 	    } */
-  /* 	  } */
-  /* 	} */
-	
-  /* 	if (bl>=br && bl>=tl && bl>=tr) { */
-  /* 	  quad=BOTTOM_LEFT; */
-  /* 	  quadTot=bl; */
-  /* 	} else if (br>=bl && br>=tl && br>=tr) { */
-  /* 	  quad=BOTTOM_RIGHT; */
-  /* 	  quadTot=br; */
-  /* 	} else if (tl>=br && tl>=bl && tl>=tr) { */
-  /* 	  quad=TOP_LEFT; */
-  /* 	  quadTot=tl; */
-  /*    	} else if   (tr>=bl && tr>=tl && tr>=br) { */
-  /* 	  quad=TOP_RIGHT; */
-  /* 	  quadTot=tr; */
-  /* 	} */
-	
-  /* 	if (max>nSigma*clusters->rms || tot>sqrt(clusterSizeY*clusterSize)*nSigma*clusters->rms || quadTot>cy*cs*nSigma*clusters->rms) { */
-  /* 	  if (clusters->get_data(0,0)>=max) { */
-  /* 	    eventMask[iy][ix]=PHOTON_MAX; */
-  /* 	  } else { */
-  /* 	    eventMask[iy][ix]=PHOTON; */
-  /* 	  } */
-  /* 	} else if (eventMask[iy][ix]==PEDESTAL) { */
-  /* 	  if (cm==0) { */
-  /* 	    if (det) */
-  /* 	      val=dataSign*det->getValue(data, ix, iy); */
-  /* 	    else */
-  /* 	      val=((double**)data)[iy][ix]; */
-  /* 	    addToPedestal(val,ix,iy); */
-  /* 	  } */
-  /* 	} */
-      
-
-
-  /*     return  eventMask[iy][ix]; */
-
-  /* }; */
 
 
 
@@ -477,8 +370,8 @@ int *getClusters(char *data,  int *ph=NULL) {
     addToCommonMode(data);
 
 
-  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)) {
 	max=0;
 	tl=0;
@@ -523,10 +416,14 @@ int *getClusters(char *data,  int *ph=NULL) {
 	    if (ir==0 && ic==0) {
 	      if (*v<-nSigma*(clusters+nph)->rms)
 		eventMask[iy][ix]=NEGATIVE_PEDESTAL;
+	      else if (*v>nSigma*(clusters+nph)->rms)
+		eventMask[iy][ix]=PHOTON;
 	    }
 	    
+	  }
 	}
-	}
+	if (eventMask[iy][ix]==PHOTON && val[iy][ix]<max) 
+	  continue;
 	
 	if (bl>=br && bl>=tl && bl>=tr) {
 	  (clusters+nph)->quad=BOTTOM_LEFT;
@@ -548,7 +445,8 @@ int *getClusters(char *data,  int *ph=NULL) {
 	  (clusters+nph)->tot=tot;
 	  (clusters+nph)->x=ix;
 	  (clusters+nph)->y=iy;
-	  (clusters+nph)->iframe=det->getFrameNumber(data);
+	  //	  (clusters+nph)->iframe=det->getFrameNumber(data);
+	  // cout << det->getFrameNumber(data) << " " << (clusters+nph)->iframe << endl;
 	  (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++) {
@@ -576,8 +474,8 @@ int *getClusters(char *data,  int *ph=NULL) {
   nphFrame=nph;
   nphTot+=nph;
   //cout << nphFrame << endl;
-  // cout <<"**********************************"<< endl;
-  writeClusters();
+  //  cout <<"**********************************"<< det->getFrameNumber(data) << " " << nphFrame << endl;
+  writeClusters(det->getFrameNumber(data));
   return  image;
   
 };
@@ -662,13 +560,27 @@ int *getClusters(char *data,  int *ph=NULL) {
 
 */
 
-static void writeClusters(FILE *f, single_photon_hit *clusters, int nph){for (int i=0; i<nph; i++) (clusters+i)->write(f);};
-void writeClusters(FILE *f){for (int i=0; i<nphFrame; i++) (clusters+i)->write(f);};
- void writeClusters(){if (myFile) {  
+    static void writeClusters(FILE *f, single_photon_hit *cl, int nph, int fn=0){
+  
+/* #ifndef OLDFORMAT   */
+/*       if (fwrite((void*)&fn, 1,  sizeof(int), f))  */
+/* 	if (fwrite((void*)&nph, 1,  sizeof(int), f))  */
+/* #endif    */
+	  for (int i=0; i<nph; i++) (cl+i)->write(f);
+};
+ void writeClusters(FILE *f, int fn=0){
+   writeClusters(f,clusters,nphFrame, fn);
+  //for (int i=0; i<nphFrame; i++) 
+  //(clusters+i)->write(f);
+};
+ void writeClusters(int fn){
+   if (myFile) {  
      //cout << "++" << endl;  
      pthread_mutex_lock(fm);
-     for (int i=0; i<nphFrame; i++) 
-       (clusters+i)->write(myFile);
+     //   cout <<"**********************************"<< fn << " " << nphFrame << endl;
+     writeClusters(myFile,clusters,nphFrame, fn);
+     // for (int i=0; i<nphFrame; i++) 
+     //  (clusters+i)->write(myFile);
      pthread_mutex_unlock(fm); 
      //cout << "--" << endl; 
    } 

slsDetectorCalibration/single_photon_hit.h

@@ -38,8 +38,62 @@ class single_photon_hit {
   size_t write(FILE *myFile) {
     //fwrite((void*)this, 1, 3*sizeof(int)+4*sizeof(double)+sizeof(quad), myFile); 
     
-    if (fwrite((void*)this, 1, sizeof(int)+2*sizeof(int16_t), myFile)) 
+    // if (fwrite((void*)this, 1, sizeof(int)+2*sizeof(int16_t), myFile))
+    //#ifdef OLDFORMAT 
+    if (fwrite((void*)&iframe, 1,  sizeof(int), myFile)) {};
+      //#endif  
+#ifndef WRITE_QUAD
+      //printf("no quad ");
+    if (fwrite((void*)&x, 2,  sizeof(int16_t), myFile)) 
       return fwrite((void*)data, 1, dx*dy*sizeof(int), myFile);
+#endif 
+#ifdef WRITE_QUAD
+    // printf("quad ");
+    int qq[4];
+    switch(quad) {
+    case TOP_LEFT:
+      qq[0]=data[3];
+      qq[1]=data[4];
+      qq[2]=data[6];
+      qq[3]=data[7];
+      x=x-1;
+      y=y;
+      break;
+      
+    case TOP_RIGHT:
+      qq[0]=data[4];
+      qq[1]=data[5];
+      qq[2]=data[7];
+      qq[3]=data[8];
+      x=x;
+      y=y;
+      break;
+
+
+    case BOTTOM_LEFT:
+      qq[0]=data[0];
+      qq[1]=data[1];
+      qq[2]=data[3];
+      qq[3]=data[4];
+      x=x-1;
+      y=y-1;
+      break;
+    case BOTTOM_RIGHT:
+      qq[0]=data[1];
+      qq[1]=data[2];
+      qq[2]=data[4];
+      qq[3]=data[5];
+      x=x;
+      y=y-1;
+      break;
+
+
+    default:
+      ;
+    }
+    if (fwrite((void*)&x, 2,  sizeof(int16_t), myFile)) 
+      return fwrite((void*)qq, 1, 4*sizeof(int), myFile);
+#endif
     return 0;
   };   
   
@@ -50,11 +104,126 @@ class single_photon_hit {
   size_t read(FILE *myFile) {
     //fread((void*)this, 1,  3*sizeof(int)+4*sizeof(double)+sizeof(quad), myFile); 
     
-    if (fread((void*)this, 1,  sizeof(int)+2*sizeof(int16_t), myFile))  
+    //#ifdef OLDFORMAT 
+    if (fread((void*)&iframe, 1,  sizeof(int), myFile))  {}
+      //#endif   
+#ifndef WRITE_QUAD
+      // printf( "no quad \n");
+    if (fread((void*)&x, 2,  sizeof(int16_t), myFile))  
       return fread((void*)data, 1, dx*dy*sizeof(int), myFile);
+#endif 
+#ifdef WRITE_QUAD
+    int qq[4];
+    // printf( "quad \n");
+    if (fread((void*)&x, 2,  sizeof(int16_t), myFile))  
+      if (fread((void*)qq, 1, 4*sizeof(int), myFile)) {
+
+	quad=TOP_RIGHT;
+	int mm=qq[0];
+	for (int i=1; i<4; i++) {
+	  if (qq[i]<mm) {
+	    switch (i) {
+	    case 1: 
+	      quad=TOP_LEFT;
+	      break;
+	    case 2:
+	      quad=BOTTOM_RIGHT;
+	      break;
+	    case 3:
+	      quad=BOTTOM_LEFT;
+	      break;
+	    default:
+	      ;
+	    }
+
+	  }
+
+	}
+
+
+    switch(quad) {
+    case TOP_LEFT:
+      data[0]=0;
+      data[1]=0;
+      data[2]=0;
+      data[3]=qq[0];
+      data[4]=qq[1];
+      data[5]=0;
+      data[6]=qq[2];
+      data[7]=qq[3];
+      data[8]=0;
+      x=x+1;
+      y=y;
+      break;
+      
+    case TOP_RIGHT:
+      data[0]=0;
+      data[1]=0;
+      data[2]=0;
+      data[3]=0;
+      data[4]=qq[0];
+      data[5]=qq[1];
+      data[6]=0;
+      data[7]=qq[2];
+      data[8]=qq[3];
+      x=x;
+      y=y;
+      break;
+
+
+    case BOTTOM_LEFT:
+      data[0]=qq[0];
+      data[1]=qq[1];
+      data[2]=0;
+      data[3]=qq[2];
+      data[4]=qq[3];
+      data[5]=0;
+      data[6]=0;
+      data[7]=0;
+      data[8]=0;
+      x=x+1;
+      y=y+1;
+      break;
+    case BOTTOM_RIGHT:
+      data[0]=0;
+      data[1]=qq[0];
+      data[2]=qq[1];
+      data[3]=0;
+      data[4]=qq[2];
+      data[5]=qq[3];
+      data[6]=0;
+      data[7]=0;
+      data[8]=0;
+      x=x;
+      y=y+1;
+      break;
+
+
+    default:
+      ;
+    }
+    return 1;
+      }
+#endif
+
     return 0;
   };
 
+  void print() {
+
+    int ix, iy;
+    
+    for (int iy=0; iy<dy; iy++) {
+      for (int ix=0; ix<dx; ix++) {
+	printf("%d \t",data[ix+iy*dx]);
+
+      }
+    printf("\n");
+    }
+
+  }
+
+
     /** 
 	assign the value to the element of the cluster matrix, with relative coordinates where the center of the cluster is (0,0)
 	\param v value to be set

slsDetectorCalibration/slsDetectorData.h

@@ -331,8 +331,6 @@ class slsDetectorData {
    */
   virtual char *readNextFrame(ifstream &filebin)=0;
   
-  
-
 };
 
 

slsDetectorCalibration/tiffIO.cpp

@@ -18,7 +18,7 @@ void *WriteToTiff(float * imgData, const char * imgname, int nrow, int ncol){
   int sampleperpixel=1;
   // unsigned char * buff=NULL;
   tsize_t linebytes;
- cout << "--" <<endl;
+  // cout << "--" <<endl;
   TIFF * tif = TIFFOpen(imgname,"w");
   if (tif) {
     TIFFSetField(tif,TIFFTAG_IMAGEWIDTH,ncol);

slsDetectorSoftware/commonFiles/communication_funcs.c

@@ -31,6 +31,9 @@ int myport=-1;
 //struct sockaddr_in address;
 //#define VERBOSE
 
+#define BLACKFIN_DRVR_SND_LMT	(30000) // rough limit
+#define BLACKFIN_RSND_PCKT_LOOP	(10)
+#define BLACKFIN_RSND_WAIT_US	(1)
 
 int bindSocket(unsigned short int port_number) {
   int i;
@@ -314,13 +317,46 @@ int receiveData(int file_des, void* buf,int length, intType itype){
 	return ret;
 }
 
-
- int sendDataOnly(int file_des, void* buf,int length) {
+int sendDataOnly(int file_des, void* buf,int length) {
      if (!length)
          return 0;
-	 int ret =  write(file_des, buf, length); //value of -1 is other end socket crash as sigpipe is ignored
-	 if (ret < 0) cprintf(BG_RED, "Error writing to socket. Possible socket crash\n");
-	 return ret;
+
+
+	  int bytesSent = 0;
+	  int retry = 0; // retry index when buffer is blocked (write returns 0)
+	  while (bytesSent < length) {
+
+		  // setting a max packet size for blackfin driver (and network driver does not do a check if packets sent)
+		  int bytesToSend = length - bytesSent;
+		  if (bytesToSend > BLACKFIN_DRVR_SND_LMT)
+			  bytesToSend = BLACKFIN_DRVR_SND_LMT;
+
+		  // send
+		  int rc = write(file_des, (char*)((char*)buf + bytesSent), bytesToSend);
+		  // error
+		  if (rc < 0) {
+		    cprintf(BG_RED, "Error writing to socket. Possible socket crash: left=%d rc=%d length=%d sent=%d\n", bytesToSend, rc, length, bytesSent);
+			  return bytesSent;
+		  }
+		  // also error, wrote nothing, buffer blocked up, too fast sending for client
+		  if (rc == 0) {
+			  cprintf(RED, "Error writing to socket. Buffer full. Retry: %d\n", retry);
+			  ++retry;
+			  // wrote nothing for many loops
+			  if (retry >= BLACKFIN_RSND_PCKT_LOOP) {
+				  cprintf(BG_RED, "Error writing to socket. Buffer full! Too fast! No more.\n");
+				  return bytesSent;
+			  }
+			  usleep(BLACKFIN_RSND_WAIT_US);
+		  }
+		  // wrote something, reset retry
+		  else  {
+			  retry = 0;
+		  }
+		  bytesSent += rc;
+	  }
+
+	  return bytesSent;
 }
 
 

slsDetectorSoftware/jctbDetectorServer/firmware_funcs.c

@@ -75,7 +75,8 @@ typedef struct ip_header_struct {
 
 struct timeval tss,tse,tsss; //for timing
 
-
+#define MAXDAC 32
+u_int32_t dac[MAXDAC];
 
 FILE *debugfp, *datafp;
 
@@ -890,107 +891,33 @@ int setTiming(int ti) {
 
 
   int ret=GET_EXTERNAL_COMMUNICATION_MODE;
-
+  int val;
   int g=-1, t=-1, rot=-1;
 
   int i;
 
+  val=bus_r(EXT_SIGNAL_REG); 
   switch (ti) {
   case AUTO_TIMING:
-    timingMode=ti;
-    // disable all gates/triggers in except if used for master/slave synchronization
-    for (i=0; i<4; i++) {
-      if (getFPGASignal(i)>0 && getFPGASignal(i)<GATE_OUT_ACTIVE_HIGH && signals[i]!=MASTER_SLAVE_SYNCHRONIZATION)
-	setFPGASignal(i,SIGNAL_OFF);
-    }
+    bus_w(EXT_SIGNAL_REG,val&~(0x1));
     break;
 
   case   TRIGGER_EXPOSURE:
-    timingMode=ti;
-    // if one of the signals is configured to be trigger, set it and unset possible gates
-    for (i=0; i<4; i++) {
-      if (signals[i]==TRIGGER_IN_RISING_EDGE ||  signals[i]==TRIGGER_IN_FALLING_EDGE)
-	setFPGASignal(i,signals[i]);
-      else if (signals[i]==GATE_IN_ACTIVE_HIGH || signals[i]==GATE_IN_ACTIVE_LOW)
-	setFPGASignal(i,SIGNAL_OFF);
-      else if (signals[i]==RO_TRIGGER_IN_RISING_EDGE ||  signals[i]==RO_TRIGGER_IN_FALLING_EDGE)
-	setFPGASignal(i,SIGNAL_OFF);
-
-    }
+    bus_w(EXT_SIGNAL_REG,val|(0x1));
     break;
 
 
-
-  case  TRIGGER_READOUT:
-    timingMode=ti;
-    // if one of the signals is configured to be trigger, set it and unset possible gates
-    for (i=0; i<4; i++) {
-      if (signals[i]==RO_TRIGGER_IN_RISING_EDGE ||  signals[i]==RO_TRIGGER_IN_FALLING_EDGE)
-	setFPGASignal(i,signals[i]);
-      else if (signals[i]==GATE_IN_ACTIVE_HIGH || signals[i]==GATE_IN_ACTIVE_LOW)
-	setFPGASignal(i,SIGNAL_OFF);
-      else if (signals[i]==TRIGGER_IN_RISING_EDGE ||  signals[i]==TRIGGER_IN_FALLING_EDGE)
-	setFPGASignal(i,SIGNAL_OFF);
-    }
-    break;
-
-  case GATE_FIX_NUMBER:
-    timingMode=ti;
-    // if one of the signals is configured to be trigger, set it and unset possible gates
-    for (i=0; i<4; i++) {
-      if (signals[i]==RO_TRIGGER_IN_RISING_EDGE ||  signals[i]==RO_TRIGGER_IN_FALLING_EDGE)
-	setFPGASignal(i,SIGNAL_OFF);
-      else if (signals[i]==GATE_IN_ACTIVE_HIGH || signals[i]==GATE_IN_ACTIVE_LOW)
-	setFPGASignal(i,signals[i]);
-      else if (signals[i]==TRIGGER_IN_RISING_EDGE ||  signals[i]==TRIGGER_IN_FALLING_EDGE)
-	setFPGASignal(i,SIGNAL_OFF);
-    }
-    break;
-
-
-
-  case GATE_WITH_START_TRIGGER:
-    timingMode=ti;
-    for (i=0; i<4; i++) {
-      if (signals[i]==RO_TRIGGER_IN_RISING_EDGE ||  signals[i]==RO_TRIGGER_IN_FALLING_EDGE)
-	setFPGASignal(i,SIGNAL_OFF);
-      else if (signals[i]==GATE_IN_ACTIVE_HIGH || signals[i]==GATE_IN_ACTIVE_LOW)
-	setFPGASignal(i,signals[i]);
-      else if (signals[i]==TRIGGER_IN_RISING_EDGE ||  signals[i]==TRIGGER_IN_FALLING_EDGE)
-	setFPGASignal(i,signals[i]);
-    }
-    break;
-
   default:
 	  break;
 
   }
 
-
-
-  for (i=0; i<4; i++) {
-    if (signals[i]!=MASTER_SLAVE_SYNCHRONIZATION) {
-      if (getFPGASignal(i)==RO_TRIGGER_IN_RISING_EDGE ||  getFPGASignal(i)==RO_TRIGGER_IN_FALLING_EDGE)
-	rot=i;
-      else if (getFPGASignal(i)==GATE_IN_ACTIVE_HIGH || getFPGASignal(i)==GATE_IN_ACTIVE_LOW)
-	g=i;
-      else if (getFPGASignal(i)==TRIGGER_IN_RISING_EDGE ||  getFPGASignal(i)==TRIGGER_IN_FALLING_EDGE)
-	t=i;
-    }
-  }
-
-
-  if (g>=0 && t>=0 && rot<0) {
-    ret=GATE_WITH_START_TRIGGER;
-  } else if (g<0 && t>=0 && rot<0) {
+  if (bus_r(EXT_SIGNAL_REG)&0x1)
     ret=TRIGGER_EXPOSURE;
-  } else if (g>=0 && t<0 && rot<0) {
-    ret=GATE_FIX_NUMBER;
-  } else if (g<0 && t<0 && rot>0) {
-    ret=TRIGGER_READOUT;
-  } else if (g<0 && t<0 && rot<0) {
+  else
     ret=AUTO_TIMING;
-  }
+
+
 
   // timingMode=ret;
 
@@ -1671,22 +1598,28 @@ int initHighVoltage(int val, int imod){
        codata=((dacvalue)&0xff);
      
     
-       valw=bus_r(offw)&0x7fff; //switch off HV 
+       valw=(bus_r(offw) & 0x7fff) | 0x00ff;
+
        bus_w(offw,(valw)); // start point
        valw=((valw&(~(0x1<<csdx))));bus_w(offw,valw); //chip sel bar down
+       
        for (i=0;i<8;i++) {
-	valw=(valw&(~(0x1<<cdx)));bus_w(offw,valw); //cldwn
-	valw=((valw&(~(0x1<<ddx)))+(((codata>>(7-i))&0x1)<<ddx));bus_w(offw,valw);//write data (i)
-	valw=((valw&(~(0x1<<cdx)))+(0x1<<cdx));bus_w(offw,valw);//clkup
+	 valw=(valw&(~(0x1<<cdx)));bus_w(offw,valw); //cldwn
+	 valw=((valw&(~(0x1<<ddx)))+(((codata>>(7-i))&0x1)<<ddx));bus_w(offw,valw);//write data (i)
+	 valw=((valw&(~(0x1<<cdx)))+(0x1<<cdx));bus_w(offw,valw);//clkup
        }
        valw=((valw&(~(0x1<<csdx)))+(0x1<<csdx));bus_w(offw,valw); //csup
        
        valw=(valw&(~(0x1<<cdx)));bus_w(offw,valw); //cldwn
+
+
        
-       if (dacvalue>=0) {
-	 valw=bus_r(offw)|0xff00;; //switch on HV
+       if (val>=0) {
+	 valw=bus_r(offw) | 0x8000;; //switch on HV
+	 printf("switch on HV in CTB4 %x\n");
        } else {
-	 valw=bus_r(offw)&0x7fff;//switch off HV
+	 valw=bus_r(offw) & 0x7fff;//switch off HV
+	 printf("switch off HV in CTB4\n");
        }
 
 
@@ -1695,6 +1628,22 @@ int initHighVoltage(int val, int imod){
        
        bus_w(HV_REG,val);
        
+
+       //added for CTB5!
+
+       // valw=bus_r(POWER_ON_REG);
+       if (val>=0) {
+	 i=bus_r(POWER_ON_REG) | 0x80000000;; //switch on HV
+	 printf("switch on HV in CTB5 %08x\n", i);
+       } else {
+	 i=bus_r(POWER_ON_REG)  & 0x7fffffff;//switch off HV
+	 printf("switch off HV in CTB5 %08x\n", i);
+       }
+       
+       bus_w(POWER_ON_REG,i);
+
+
+
 	 // } 
   }
 
@@ -3212,12 +3161,22 @@ int setDacRegister(int dacnum,int dacvalue) {
 /*   printf("Dac register %x wrote %08x\n",(DAC_REG_OFF+dacnum/2)<<11,val); */
 /*     bus_w((DAC_REG_OFF+dacnum/2)<<11, val); */
     
-
+#ifdef CTB4
   bus_w(DAC_NUM_REG, dacnum);
   bus_w(DAC_VAL_REG, dacvalue);
   bus_w(DAC_NUM_REG, dacnum | (1<<16));
   bus_w(DAC_NUM_REG, dacnum);
   printf("Wrote dac register value %d address %d\n",bus_r(DAC_VAL_REG),bus_r(DAC_NUM_REG)) ;
+#endif
+  
+#ifndef CTB4
+  
+  if (dacnum<MAXDAC) {
+    printf("dac variable %d set to  value %d\n",dacnum, dacvalue) ;
+    dac[dacnum]=dacvalue;
+  } else
+    printf("Bad dac index %d\n", dacnum);
+#endif
   return getDacRegister(dacnum);
   
 
@@ -3225,9 +3184,22 @@ int setDacRegister(int dacnum,int dacvalue) {
 int getDacRegister(int dacnum) {
 
   
+#ifdef CTB4
   bus_w(DAC_NUM_REG, dacnum);
   printf("READ dac register value %d address %d\n",(int16_t)bus_r(DAC_VAL_OUT_REG),bus_r(DAC_NUM_REG)) ;
   return (int16_t)bus_r(DAC_VAL_OUT_REG);
+#endif
+
+#ifndef CTB4
+  if (dacnum<MAXDAC) {
+    printf("dac variable %d is %d\n",dacnum,dac[dacnum]) ;
+    return dac[dacnum];
+  }else
+    printf("Bad dac index %d\n", dacnum);
+#endif
+
+  return -1;
+
 /* #define DAC_VAL_REG 121<<11 */
 /* #define DAC_NUM_REG 122<<11 */
 /* #define DAC_VAL_OUT_REG 42<<11 */
@@ -3533,6 +3505,7 @@ int powerChip(int arg) {
   
   u_int32_t preg=bus_r(POWER_ON_REG); 
   if (myDetectorType!=JUNGFRAUCTB) {
+    printf("This is not a JCTB!\n");
     if (arg>=0) {
       if (arg)
 	bus_w(POWER_ON_REG,preg|0xffff0000);

slsDetectorSoftware/jctbDetectorServer/gitInfo.txt

@@ -1,9 +1,9 @@
 Path: slsDetectorsPackage/slsDetectorSoftware/jctbDetectorServer 
 URL: origin git@github.com:slsdetectorgroup/slsDetectorPackage.git 
 Repository Root: origin git@github.com:slsdetectorgroup/slsDetectorPackage.git 
-Repsitory UUID: 83600fcb15c8261173ab15a8ba8d1009693f2d23 
-Revision: 27 
-Branch: anna 
-Last Changed Author: Anna_Bergamaschi 
-Last Changed Rev: 3962 
-Last Changed Date: 2018-09-04 11:59:37.000000002 +0200 ./firmware_funcs.c
+Repsitory UUID: 9e5ec6a57bf206fe9384260fc4040867a9bbd71c 
+Revision: 33 
+Branch: developer 
+Last Changed Author: Erik_Frojdh 
+Last Changed Rev: 4065 
+Last Changed Date: 2019-01-30 17:38:35.000000002 +0100 ./firmware_funcs.c

slsDetectorSoftware/jctbDetectorServer/gitInfoMoench.h

@@ -1,6 +1,6 @@
 #define GITURL "git@github.com:slsdetectorgroup/slsDetectorPackage.git"
-#define GITREPUUID "83600fcb15c8261173ab15a8ba8d1009693f2d23"
-#define GITAUTH "Anna_Bergamaschi"
-#define GITREV 0x3962
-#define GITDATE 0x20180904
-#define GITBRANCH "anna"
+#define GITREPUUID "9e5ec6a57bf206fe9384260fc4040867a9bbd71c"
+#define GITAUTH "Erik_Frojdh"
+#define GITREV 0x4065
+#define GITDATE 0x20190130
+#define GITBRANCH "developer"

slsDetectorSoftware/jctbDetectorServer/server_funcs.c

@@ -1877,6 +1877,7 @@ int start_acquisition(int file_des) {
     ret=FAIL;
     sprintf(mess,"Detector locked by %s\n",lastClientIP);  
   } else {
+	  nframes = 0;
     ret=startStateMachine();
   }
   if (ret==FAIL)
@@ -2057,7 +2058,7 @@ int read_frame(int file_des) {
 	  printf("sending pointer %x of size %d\n",(unsigned int)(dataretval),dataBytes);
 	 #endif
 	  n=sendDataOnly(file_des,dataretval,dataBytes);
-	  printf("Sent %d bytes\n",n);
+	  printf("Frame %d, Sent %d bytes\n", nframes, n);
   } else  {
       if (getFrames()>-1) {
 	dataret=FAIL;
@@ -2122,7 +2123,7 @@ int start_and_read_all(int file_des) {
     return dataret;
 
   }
-
+  nframes = 0;
   startStateMachine();
 
   /*  ret=startStateMachine();  

slsDetectorSoftware/multiSlsDetector/multiSlsDetector.cpp

@@ -2380,7 +2380,7 @@ int* multiSlsDetector::startAndReadAll() {
 		while ((retval = getDataFromDetector())) {
 			++i;
 #ifdef VERBOSE
-			std::cout << i << std::endl;
+			  std::cout << i << " " retval << std::endl;
 #endif
 			dataQueue.push(retval);
 		}
@@ -2458,7 +2458,7 @@ int* multiSlsDetector::getDataFromDetector() {
 	int nodatadet          = -1;
 	int nodatadetectortype = false;
 	detectorType types     = getDetectorsType();
-	if (types == EIGER || types == JUNGFRAU || GOTTHARD || PROPIX) {
+	if (types == EIGER || types == JUNGFRAU || types == GOTTHARD || types == PROPIX) {
 		nodatadetectortype = true;
 	}
 

slsDetectorSoftware/slsDetector/slsDetector.cpp

@@ -145,7 +145,7 @@ double* slsDetector::decodeData(int *datain, int &nn, double *fdata) {
 				dataout[ichan]=*((u_int16_t*)ptr);
 				ptr+=2;
 			}
-			std::cout<< "decoded "<< ichan << " channels" << std::endl;
+			//std::cout<< "decoded "<< ichan << " channels" << std::endl;
 		} else {
 			switch (nbits) {
 			case 1:
@@ -3821,8 +3821,8 @@ int* slsDetector::startAndReadAll() {
 	//#ifdef VERBOSE
 #ifdef VERBOSE
 	int i=0;
-#endif
 	//#endif
+#endif
 	if(thisDetector->myDetectorType == EIGER) {
 		if (prepareAcquisition() == FAIL)
 			return NULL;
@@ -3832,15 +3832,15 @@ int* slsDetector::startAndReadAll() {
 	// std::cout<< "started" << std::endl;
 	//#endif
 	while ((retval=getDataFromDetector())){
-#ifdef VERBOSE
+	#ifdef VERBOSE
 		++i;
-		std::cout<< i << std::endl;
+		//	std::cout<< i << std::endl;
 		//#else
 		//std::cout<< "-" << flush;
-#endif
+	#endif
 		dataQueue.push(retval);
 
-		//std::cout<< "pushed" << std::endl;
+		//	std::cout<< "pushed" << retval << std::endl;
 	}
 	disconnectControl();
 
@@ -3848,7 +3848,7 @@ int* slsDetector::startAndReadAll() {
 	std::cout<< "received "<< i<< " frames" << std::endl;
 	//#else
 	// std::cout << std::endl;
-#endif
+	#endif
 	return dataQueue.front(); // check what we return!
 	/* while ((retval=getDataFromDetectorNoWait()))
      ++i;
@@ -3892,10 +3892,13 @@ int* slsDetector::getDataFromDetector(int *retval) {
 
 	int nodatadetectortype = false;
 	detectorType types = getDetectorsType();
-	if(types == EIGER || types == JUNGFRAU || GOTTHARD || PROPIX){
+	//	cout << types << endl;
+	
+	if(types == EIGER || types == JUNGFRAU || types == GOTTHARD || types == PROPIX){
 		nodatadetectortype = true;
 	}
 
+	//cout << "nodata det" << nodatadetectortype << endl;
 
 	if (!nodatadetectortype && retval==NULL)
 		retval=new int[nel];
@@ -3917,22 +3920,22 @@ int* slsDetector::getDataFromDetector(int *retval) {
 			std::cout<< "Detector returned: " << mess << " " << n << std::endl;
 		} else {
 			;
-#ifdef VERBOSE
+			#ifdef VERBOSE
 			std::cout<< "Detector successfully returned: " << mess << " " << n
 					<< std::endl;
-#endif
+		#endif
 		}
 		if ((!nodatadetectortype) && (r==NULL)){
 			delete [] retval;
 		}
 		return NULL;
 	} else 	if (!nodatadetectortype){
+	  // cout <<"??" << endl;
+	  n=controlSocket->ReceiveDataOnly(retval,thisDetector->dataBytes);
 
-		n=controlSocket->ReceiveDataOnly(retval,thisDetector->dataBytes);
-
-#ifdef VERBOSE
+	#ifdef VERBOSE
 		std::cout<< "Received "<< n << " data bytes" << std::endl;
-#endif
+	#endif
 		if (n!=thisDetector->dataBytes) {
 			std::cout<< "wrong data size received from detector: received " <<
 					n << " but expected " << thisDetector->dataBytes << std::endl;
@@ -3948,7 +3951,8 @@ int* slsDetector::getDataFromDetector(int *retval) {
 
 
 	}
-	//	cout << "get data returning " << endl;
+	
+	//	cout << "get data returning " << retval << endl;
 	//	cout << endl;
 	return retval;
 
@@ -4962,9 +4966,9 @@ int slsDetector::setReadOutFlags(readOutFlags flag) {
 	char mess[MAX_STR_LENGTH]="";
 	int ret=OK;
 
-#ifdef VERBOSE
-	std::cout<< "Setting readout flags to "<< flag << std::endl;
-#endif
+	//#ifdef VERBOSE
+	std::cout<< "Setting readout flags to "<< hex << flag << dec << std::endl;
+	//#endif
 
 	if (thisDetector->onlineFlag==ONLINE_FLAG) {
 		if (connectControl() == OK){
@@ -4980,7 +4984,8 @@ int slsDetector::setReadOutFlags(readOutFlags flag) {
 				thisDetector->roFlags=(readOutFlags)retval;
 				if (thisDetector->myDetectorType==JUNGFRAUCTB) {
 
-					getTotalNumberOfChannels();
+					int nn=getTotalNumberOfChannels();
+					cout << "Total number of channels is " << nn << endl;
 					//thisDetector->dataBytes=getTotalNumberOfChannels()*
 					//thisDetector->dynamicRange/8*thisDetector->timerValue[SAMPLES_JCTB];
 				}
@@ -5033,27 +5038,27 @@ int slsDetector::setReadOutFlags(readOutFlags flag) {
 	  char h[1000];
 	  switch (flag) {
 	  case PEDESTAL:
-	    retval=PEDESTAL;
+	    retval|=PEDESTAL;
 	    strcpy(h,"\"frameMode\":\"pedestal\"");
 	    break;
 	  case NEWPEDESTAL:
-	    retval=NEWPEDESTAL;
+	    retval|=NEWPEDESTAL;
 	    strcpy(h,"\"frameMode\":\"newPedestal\"");
 	    break;
 	  case FLAT:
-	    retval=FLAT;
+	    retval|=FLAT;
 	    strcpy(h,"\"frameMode\":\"flatfield\"");
 	    break;
 	  case NEWFLAT:
-	    retval=NEWFLAT;
+	    retval|=NEWFLAT;
 	    strcpy(h,"\"frameMode\":\"newFlatfield\"");
 	    break;
 	  default:
-	    retval=FRAME;
+	    retval|=FRAME;
 	    strcpy(h,"\"frameMode\":\"frame\"");
 	  }
 	  if (flag!=GET_READOUT_FLAGS) {
-	  if (header.length()>0) {
+	    if (header.length()>0) {
 	    if (header.at(0)==',')
 	      header.erase(0,1);
 	    if (header.length()>0)
@@ -5103,14 +5108,14 @@ int slsDetector::setReadOutFlags(readOutFlags flag) {
 	  switch (flag) {
 	  case COUNTING:
 	    strcpy(h,"\"detectorMode\":\"counting\"");
-	    retval=COUNTING;
+	    retval|=COUNTING;
 	    break;
 	  case INTERPOLATING:
-	    retval=INTERPOLATING;
+	    retval|=INTERPOLATING;
 	    strcpy(h,"\"detectorMode\":\"interpolating\"");
 	    break;
 	  default:
-	    retval=ANALOG;
+	    retval|=ANALOG;
 	    strcpy(h,"\"detectorMode\":\"analog\"");
 	  }
 	  if (flag!=GET_READOUT_FLAGS) {

slsDetectorSoftware/slsDetector/slsDetectorActions.cpp

@@ -176,31 +176,33 @@ int slsDetectorActions::setScan(int iscan, string script, int nvalues, double *v
     if (par!="")
       strcpy(scanParameter[iscan],par.c_str());
       
-      if (nvalues>=0) {    
-	if (nvalues==0)
-	  scanMode[iscan]=noScan;
-	else {
-	  nScanSteps[iscan]=nvalues;
+    if (nvalues>=0) {
+      //  cout << "nvalues " << nvalues << endl;
+      if (nvalues==0)
+	scanMode[iscan]=noScan;
+      else {
+	nScanSteps[iscan]=nvalues;
 	  if (nvalues>MAX_SCAN_STEPS)
 	    nScanSteps[iscan]=MAX_SCAN_STEPS;
-	}
+      }
       }
       
-      if (values && 	scanMode[iscan]>0 ) {
-	for (int iv=0; iv<nScanSteps[iscan]; iv++) {
-	  scanSteps[iscan][iv]=values[iv];
-	}
+    if (values && 	scanMode[iscan]>0 ) {
+      for (int iv=0; iv<nScanSteps[iscan]; iv++) {
+	scanSteps[iscan][iv]=values[iv];
+	//	cout << values[iv] << endl;
       }
+    }
 
-      if (precision>=0)
-	scanPrecision[iscan]=precision;
-      
-      if (scanMode[iscan]>0){
+    if (precision>=0)
+      scanPrecision[iscan]=precision;
+    
+    if (scanMode[iscan]>0){
 	*actionMask |= 1<< (iscan+MAX_ACTIONS);
-      } else {
-	*actionMask &= ~(1 <<  (iscan+MAX_ACTIONS));
-      }
-
+    } else {
+      *actionMask &= ~(1 <<  (iscan+MAX_ACTIONS));
+    }
+    
 
 
       setTotalProgress();

slsDetectorSoftware/slsDetectorAnalysis/fileIO.cpp

@@ -125,7 +125,7 @@ int fileIO::writeDataFile(void *data, int iframe) {
 
 
 int fileIO::closeDataFile() {
-  cout << "close file...." << endl;
+  // cout << "close file...." << endl;
   if (filefd)
       fclose(filefd);
     filefd=NULL;

slsReceiverSoftware/include/genericSocket.h

@@ -666,7 +666,9 @@ public:
 			if (tcpfd<0) return -1;
 			while(length>0){
 				nsending = (length>packet_size) ? packet_size:length;
+				//	std::cout << "*"<<nsending << std::endl;
 				nsent = read(tcpfd,(char*)buf+total_sent,nsending);
+				//	std::cout << "+"<<nsent << std::endl;
 				if(!nsent) {
 					if(!total_sent) {
 						return -1; //to handle it
@@ -675,6 +677,7 @@ public:
 				}
 				length-=nsent;
 				total_sent+=nsent;
+				// std::cout << "+"<< length << " " << total_sent << std::endl;
 			}
 
 			if (total_sent>0)

slsReceiverSoftware/src/Listener.cpp

@@ -557,14 +557,20 @@ uint32_t Listener::ListenToAnImage(char* buf) {
 		lastCaughtFrameIndex = fnum;
 
 
-#ifdef VERBOSE
+		//#ifdef VERBOSE
 		//if (!index)
 		cprintf(GREEN,"Listening %d: currentfindex:%lu, fnum:%lu,   pnum:%u numpackets:%u\n",
 				index,currentFrameIndex, fnum, pnum, numpackets);
-#endif
+		//#endif
 		if (!measurementStartedFlag)
 			RecordFirstIndices(fnum);
 
+		  if (pnum >= pperFrame ) {
+			cprintf(RED,"bad packet, throwing away. packets caught so far: %d\n", numpackets);
+
+			return 0;	// bad packet
+		  }
+
 		//future packet	by looking at image number  (all other detectors)
 		if (fnum != currentFrameIndex) {
 			//cprintf(RED,"setting carry over flag to true num:%llu nump:%u\n",fnum, numpackets );
@@ -600,6 +606,7 @@ uint32_t Listener::ListenToAnImage(char* buf) {
 				memcpy(buf + fifohsize + (pnum * dsize) - 2, listeningPacket + hsize, dsize+2);
 			break;
 		case JUNGFRAUCTB:
+
 			if (pnum == (pperFrame-1))
 				memcpy(buf + fifohsize + (pnum * dsize), listeningPacket + hsize, corrected_dsize);
 			else