dont add detector id for .trim

slsDetectorCalibration/circularFifo.h

@@ -0,0 +1,156 @@
+#pragma once
+/* CircularFifo.h 
+* Code & platform dependent issues with it was originally
+* published at http://www.kjellkod.cc/threadsafecircularqueue
+* 2009-11-02
+* @author Kjell Hedstr<74>m, hedstrom@kjellkod.cc
+* modified by the sls detetor group
+* */
+
+//#include "sls_receiver_defs.h"
+#include <semaphore.h>
+#include <vector>
+#include <iostream>
+
+typedef  double double32_t;
+typedef  float float32_t;
+typedef  int int32_t;
+
+
+
+/** Circular Fifo (a.k.a. Circular Buffer)
+* Thread safe for one reader, and one writer */
+template<typename Element>
+class CircularFifo {
+public:
+
+   CircularFifo(unsigned int Size) : tail(0), head(0){
+	   Capacity = Size + 1;
+	   array.resize(Capacity);
+	   sem_init(&data_mutex,0,0);
+	   sem_init(&free_mutex,0,Size);
+   }
+   virtual ~CircularFifo() {
+       sem_destroy(&data_mutex);
+	   sem_destroy(&free_mutex);
+   }
+
+   bool push(Element*& item_, bool no_block=false);
+   bool pop(Element*& item_, bool no_block=false);
+
+   bool isEmpty() const;
+   bool isFull() const;
+
+   int getDataValue() const;
+   int getFreeValue() const;
+
+private:
+   std::vector <Element*> array;
+   unsigned int tail; // input index
+   unsigned int head; // output index
+   unsigned int Capacity;
+   mutable sem_t data_mutex;
+   mutable sem_t free_mutex;
+   unsigned int increment(unsigned int idx_) const;
+};
+
+template<typename Element>
+int CircularFifo<Element>::getDataValue() const
+{
+   int value;
+   sem_getvalue(&data_mutex, &value);
+   return value;
+}
+
+template<typename Element>
+int CircularFifo<Element>::getFreeValue() const
+{
+    int value;
+    sem_getvalue(&free_mutex, &value);
+    return value;
+}
+
+
+/** Producer only: Adds item to the circular queue.
+* If queue is full at 'push' operation no update/overwrite
+* will happen, it is up to the caller to handle this case
+*
+* \param item_ copy by reference the input item
+* \param no_block if true, return immediately if fifo is full
+* \return whether operation was successful or not */
+template<typename Element>
+bool CircularFifo<Element>::push(Element*& item_, bool no_block)
+{
+    // check for fifo full
+    if (no_block && isFull())
+        return false;
+
+    sem_wait(&free_mutex);
+    array[tail] = item_;
+    tail = increment(tail);
+    sem_post(&data_mutex);
+    return true;
+}
+
+/** Consumer only: Removes and returns item from the queue
+* If queue is empty at 'pop' operation no retrieve will happen
+* It is up to the caller to handle this case
+*
+* \param item_ return by reference the wanted item
+* \param no_block if true, return immediately if fifo is full
+* \return whether operation was successful or not */
+template<typename Element>
+bool CircularFifo<Element>::pop(Element*& item_, bool no_block)
+{
+    // check for fifo empty
+    if (no_block && isEmpty())
+        return false;
+
+   sem_wait(&data_mutex);
+   item_ = array[head];
+   head = increment(head);
+   sem_post(&free_mutex);
+   return true;
+}
+
+/** Useful for testinng and Consumer check of status
+  * Remember that the 'empty' status can change quickly
+  * as the Procuder adds more items.
+  *
+  * \return true if circular buffer is empty */
+template<typename Element>
+bool CircularFifo<Element>::isEmpty() const
+{
+    return (getDataValue() == 0);
+}
+
+/** Useful for testing and Producer check of status
+  * Remember that the 'full' status can change quickly
+  * as the Consumer catches up.
+  *
+  * \return true if circular buffer is full.  */
+template<typename Element>
+bool CircularFifo<Element>::isFull() const
+{
+    return (getFreeValue() == 0);
+}
+
+/** Increment helper function for index of the circular queue
+* index is inremented or wrapped
+*
+*  \param idx_ the index to the incremented/wrapped
+*  \return new value for the index */
+template<typename Element>
+unsigned int CircularFifo<Element>::increment(unsigned int idx_) const
+{
+   // increment or wrap
+   // =================
+   //    index++;
+   //    if(index == array.lenght) -> index = 0;
+   //
+   //or as written below:
+   //    index = (index+1) % array.length
+   idx_ = (idx_+1) % Capacity;
+   return idx_;
+}
+

slsDetectorCalibration/dataStructures/jungfrauHighZSingleChipData.h

@@ -0,0 +1,238 @@
+#ifndef JUNGFRAUHIGHZSINGLECHIPDATA_H
+#define  JUNGFRAUHIGHZSINGLECHIPDATA_H
+#include "slsDetectorData.h"
+
+//#define VERSION_V2
+    /**
+        @short  structure for a Detector Packet or Image Header
+        @li frameNumber is the frame number
+        @li expLength is the subframe number (32 bit eiger) or real time exposure time in 100ns (others)
+        @li packetNumber is the packet number
+        @li bunchId is the bunch id from beamline
+        @li timestamp is the time stamp with 10 MHz clock
+        @li modId is the unique module id (unique even for left, right, top, bottom)
+        @li xCoord is the x coordinate in the complete detector system
+        @li yCoord is the y coordinate in the complete detector system
+        @li zCoord is the z coordinate in the complete detector system
+        @li debug is for debugging purposes
+        @li roundRNumber is the round robin set number
+        @li detType is the detector type see :: detectorType
+        @li version is the version number of this structure format
+    */
+    typedef struct {
+        uint64_t bunchNumber;    /**< is the frame number */
+        uint64_t pre;    /**< something */
+
+    } jf_header;
+
+
+
+
+class jungfrauHighZSingleChipData : public slsDetectorData<uint16_t> {
+ 
+ private:
+  
+  int iframe;
+ public:
+  /**
+     Implements the slsReceiverData structure for the moench02 prototype read out by a module i.e. using the slsReceiver
+     (160x160 pixels, 40 packets 1286 large etc.)
+     \param c crosstalk parameter for the output buffer
+
+  */
+ jungfrauHighZSingleChipData(): slsDetectorData<uint16_t>(256, 256, 256*256*2+sizeof(jf_header)) {
+
+
+     for (int ix=0; ix<256; ix++) {
+       for (int iy=0; iy<256; iy++) {
+	 dataMap[iy][ix]=sizeof(jf_header)+(256*iy+ix)*2;
+#ifdef HIGHZ
+	 dataMask[iy][ix]=0x3fff;
+#endif
+	 }
+     }
+    
+    iframe=0;
+    //  cout << "data struct created" << endl;
+  };
+    
+
+ 
+  /**
+     Returns the value of the selected channel for the given dataset as double.
+     \param data pointer to the dataset (including headers etc)
+     \param ix pixel number in the x direction
+     \param iy pixel number in the y direction
+     \returns data for the selected channel, with inversion if required as double
+
+  */
+  virtual double getValue(char *data, int ix, int iy=0) {
+   
+    uint16_t val=getChannel(data, ix, iy)&0x3fff;
+    return val;
+  };
+
+   
+
+  /* virtual void calcGhost(char *data, int ix, int iy) { */
+  /*   double val=0; */
+  /*   ghost[iy][ix]=0; */
+    
+  /* } */
+
+
+
+  /* virtual void calcGhost(char *data) { */
+  /*   for (int ix=0; ix<25; ix++){ */
+  /*     for (int iy=0; iy<200; iy++) { */
+  /* 	calcGhost(data, ix,iy); */
+  /*     }  */
+  /*   } */
+  /*   // cout << "*" << endl; */
+  /* } */
+
+
+  /* double getGhost(int ix, int iy) { */
+  /*   return 0; */
+  /* }; */
+
+     /**
+
+     Returns the frame number for the given dataset. Purely virtual func.
+     \param buff pointer to the dataset
+     \returns frame number
+
+  */
+
+/* class jfrau_packet_header_t { */
+/*  public: */
+/* 	unsigned char reserved[4]; */
+/* 	unsigned char packetNumber[1]; */
+/* 	unsigned char frameNumber[3]; */
+/* 	unsigned char bunchid[8]; */
+/* }; */
+
+
+
+  int getFrameNumber(char *buff){return ((jf_header*)buff)->bunchNumber;};//*((int*)(buff+5))&0xffffff;};   
+
+  /**
+
+     Returns the packet number for the given dataset. purely virtual func
+     \param buff pointer to the dataset
+     \returns packet number number
+
+
+
+  */
+   int getPacketNumber(char *buff){return 0;}//((*(((int*)(buff+4))))&0xff)+1;};   
+
+/*    /\** */
+
+/*      Loops over a memory slot until a complete frame is found (i.e. all packets 0 to nPackets, same frame number). purely virtual func */
+/*      \param data pointer to the memory to be analyzed */
+/*      \param ndata reference to the amount of data found for the frame, in case the frame is incomplete at the end of the memory slot */
+/*      \param dsize size of the memory slot to be analyzed */
+/*      \returns pointer to the beginning of the last good frame (might be incomplete if ndata smaller than dataSize), or NULL if no frame is found  */
+
+/*   *\/ */
+/*     virtual  char *findNextFrame(char *data, int &ndata, int dsize){ndata=dsize; setDataSize(dsize);  return data;}; */
+
+
+/*    /\** */
+
+/*      Loops over a file stream until a complete frame is found (i.e. all packets 0 to nPackets, same frame number). Can be overloaded for different kind of detectors!  */
+/*      \param filebin input file stream (binary) */
+/*      \returns pointer to the begin of the last good frame, NULL if no frame is found or last frame is incomplete */
+
+/*   *\/ */
+/*     virtual char *readNextFrame(ifstream &filebin){ */
+/*       //	int afifo_length=0;   */
+/*       uint16_t *afifo_cont;  */
+/*       int ib=0; */
+/*       if (filebin.is_open()) { */
+/* 	afifo_cont=new uint16_t[dataSize/2]; */
+/*  	while (filebin.read(((char*)afifo_cont)+ib,2)) { */
+/* 	  ib+=2; */
+/* 	  if (ib==dataSize) break; */
+/* 	} */
+/* 	if (ib>0) { */
+/* 	  iframe++; */
+/* 	  // cout << ib << "-" << endl; */
+/* 	  return (char*)afifo_cont; */
+/* 	} else { */
+/* 	  delete [] afifo_cont; */
+/* 	  return NULL; */
+/* 	} */
+/*       }      */
+/*       return NULL; */
+/*     }; */
+
+
+  virtual char *readNextFrame(ifstream &filebin) {
+    int ff=-1, np=-1;
+    return readNextFrame(filebin, ff, np);
+  };
+
+  virtual char *readNextFrame(ifstream &filebin, int &ff) {
+    int np=-1;
+    return readNextFrame(filebin, ff, np);
+  };
+
+  virtual char *readNextFrame(ifstream &filebin, int& ff, int &np) {
+	  char *data=new char[dataSize];
+	  char *d=readNextFrame(filebin, ff, np, data);
+	  if (d==NULL) {delete [] data; data=NULL;}
+	  return data;
+  }
+
+
+
+
+  virtual char *readNextFrame(ifstream &filebin, int& ff, int &np, char *data) {
+	  char *retval=0;
+	  int  nd;
+	  int fnum = -1;
+	  np=0;
+	  int  pn;
+	  
+	  //  cout << dataSize << endl;
+	  if (ff>=0)
+	    fnum=ff;
+
+	  if (filebin.is_open()) {
+	    if (filebin.read(data, dataSize) ){
+	      ff=getFrameNumber(data);
+	      np=getPacketNumber(data);
+	      return data;
+	    }
+	  }
+	  return NULL;
+  };
+
+
+
+  /**
+
+     Loops over a memory slot until a complete frame is found (i.e. all packets 0 to nPackets, same frame number). purely virtual func
+     \param data pointer to the memory to be analyzed
+     \param ndata reference to the amount of data found for the frame, in case the frame is incomplete at the end of the memory slot
+     \param dsize size of the memory slot to be analyzed
+     \returns pointer to the beginning of the last good frame (might be incomplete if ndata smaller than dataSize), or NULL if no frame is found 
+     
+  */
+  virtual  char *findNextFrame(char *data, int &ndata, int dsize){
+    if (dsize<dataSize) ndata=dsize;
+    else ndata=dataSize;
+    return data;
+
+  }
+
+  //int getPacketNumber(int x, int y) {return dataMap[y][x]/packetSize;};
+
+};
+
+
+
+
+#endif

slsDetectorCalibration/dataStructures/moench03T1ZmqDataNew.h

@@ -14,6 +14,8 @@ class moench03T1ZmqDataNew : public slsDetectorData<uint16_t> {
   const int nSamples;
   const int offset;
 
+
+
   double ghost[200][25];
   double xtalk;
 
@@ -29,17 +31,22 @@ class moench03T1ZmqDataNew : public slsDetectorData<uint16_t> {
 
   */
   // moench03T1ZmqDataNew(int ns=5000): slsDetectorData<uint16_t>(400, 400, ns*32*2+sizeof(int)), nSamples(ns), offset(sizeof(int)), xtalk(0.00021) {
- moench03T1ZmqDataNew(int ns=5000, int oo=0): slsDetectorData<uint16_t>(400, 400, ns*32*2+oo), nSamples(ns), offset(oo), xtalk(0.00021) {
-
+ moench03T1ZmqDataNew(int ns=5000, int oo=2*2): slsDetectorData<uint16_t>(400, 400, ns*32*2+oo), nSamples(ns), offset(oo), xtalk(0.00021) {
+    cout << "M0.3" << endl;
     int nadc=32;
     int sc_width=25;
     int sc_height=200;
-
-    int adc_nr[32]={300,325,350,375,300,325,350,375,		\
+    
+     int adc_nr[32]={300,325,350,375,300,325,350,375,	\
     		    200,225,250,275,200,225,250,275,\
     		    100,125,150,175,100,125,150,175,\
     		    0,25,50,75,0,25,50,75};
-
+   
+     /*   int adc_nr[32]={350,375,150,175,350,375,150,175,	\
+    		    300,325,100,125,300,325,100,125,\
+    		    250,275,50,75,250,275,50,75,\
+    		    200,225,0,25,200,225,0,25};
+     */
     int row, col;
     
     int isample;
@@ -50,11 +57,11 @@ class moench03T1ZmqDataNew : public slsDetectorData<uint16_t> {
     int i;
     int adc4(0);
   
-    for (int ip=0; ip<npackets; ip++) {
-      for (int is=0; is<128; is++) {
-
+    //for (int ip=0; ip<npackets; ip++) {
+    // for (int is=0; is<128; is++) {
+    for (i=0; i<nSamples; i++) {
 	for (iadc=0; iadc<nadc; iadc++) {
-	  i=128*ip+is;
+	  //i=128*ip+is;
 	  adc4=(int)iadc/4;
 	  if (i<sc_width*sc_height) {
 	    //  for (int i=0; i<sc_width*sc_height; i++) {
@@ -70,7 +77,7 @@ class moench03T1ZmqDataNew : public slsDetectorData<uint16_t> {
 	  }
 	}
       }
-    }
+  // }
 
     
     int ii=0;

slsDetectorCalibration/energyCalibration.cpp

@@ -0,0 +1,798 @@
+#include "energyCalibration.h"
+
+#ifdef __CINT
+#define MYROOT
+#endif
+
+
+#ifdef MYROOT
+#include <TMath.h>
+#include <TH1F.h>
+#include <TH2F.h>
+#include <TGraphErrors.h>
+#endif
+
+#include <iostream>
+
+#define max(a,b) ((a) > (b) ? (a) : (b))
+#define min(a,b) ((a) < (b) ? (a) : (b))
+#define ELEM_SWAP(a,b) { register int t=(a);(a)=(b);(b)=t; }
+
+
+using namespace std;
+
+#ifdef MYROOT
+
+Double_t energyCalibrationFunctions::pedestal(Double_t *x, Double_t *par) { 
+    return par[0]-par[1]*sign*x[0];                                      
+}
+
+
+Double_t energyCalibrationFunctions::gaussChargeSharing(Double_t *x, Double_t *par) { 
+  Double_t f, arg=0;
+  // Gaussian exponent
+  if (par[3]!=0) {
+    arg=sign*(x[0]-par[2])/par[3];
+  }
+  // the Gaussian
+  f=TMath::Exp(-1*arg*arg/2.); 
+  // Gaussian + error function
+  f=f+par[5]/2.*(TMath::Erfc(arg/(TMath::Sqrt(2.)))); 
+  // Gaussian + error function + pedestal
+  return par[4]*f+pedestal(x,par);                                       
+}
+
+Double_t energyCalibrationFunctions::gaussChargeSharingKb(Double_t *x, Double_t *par) { 
+  Double_t f, arg=0,argb=0;
+  // Gaussian exponent
+  if (par[3]!=0) {
+    arg=sign*(x[0]-par[2])/par[3];
+    argb=sign*(x[0]-(par[6]*par[2]))/par[3]; // using absolute kb mean might seem better but like this the ratio can be fixed
+  }
+  // the Gaussian
+  f=TMath::Exp(-1*arg*arg/2.); 
+  f=f+par[7]*(TMath::Exp(-1*argb*argb/2.));
+  // Gaussian + error function
+  f=f+par[5]/2.*(TMath::Erfc(arg/(TMath::Sqrt(2.)))); 
+  f=f+par[7]*par[5]/2.*(TMath::Erfc(argb/(TMath::Sqrt(2.)))); 
+  // Gaussian + error function + pedestal
+  return par[4]*f+pedestal(x,par);                                       
+}
+
+Double_t energyCalibrationFunctions::gaussChargeSharingKaDoublet(Double_t *x, Double_t *par) { 
+  Double_t f, f2, arg=0, arg2=0;
+  // Gaussian exponent
+  if (par[3]!=0) {
+    arg=sign*(x[0]-par[2])/par[3];
+    arg2=sign*(x[0]-par[6])/par[3];
+  }
+  // the Gaussian
+  f=TMath::Exp(-1*arg*arg/2.); 
+  f2=TMath::Exp(-1*arg2*arg2/2.); 
+  // Gaussian + error function
+  f=f+par[5]/2.*(TMath::Erfc(arg/(TMath::Sqrt(2.)))); 
+  f2=f2+par[5]/2.*(TMath::Erfc(arg/(TMath::Sqrt(2.)))); // shouldn't this be arg2?
+  // Gaussian + error function + pedestal
+  return par[4]*f+par[7]*f2+pedestal(x,par);                                       
+}
+
+Double_t energyCalibrationFunctions::gaussChargeSharingPixel(Double_t *x, Double_t *par) { 
+  Double_t f;
+  if (par[3]<=0 || par[2]*(*x)<=0 || par[5]<0 || par[4]<=0) return 0;
+
+  Double_t pp[3];
+
+  pp[0]=0;
+  pp[1]=par[2];
+  pp[2]=par[3];
+  
+
+  f=(par[5]-par[6]*(TMath::Log(*x/par[2])))*erfBox(x,pp);
+  f+=par[4]*TMath::Gaus(*x, par[2], par[3], kTRUE);
+  return f+pedestal(x,par);                                       
+}
+
+Double_t energyCalibrationFunctions::erfBox(Double_t *z, Double_t *par) {
+
+  
+  
+  Double_t m=par[0];
+  Double_t M=par[1];
+  
+  if (par[0]>par[1]) {
+    m=par[1];
+    M=par[0];
+  }
+  
+  if (m==M)
+    return 0;
+
+
+  if (par[2]<=0) {
+    if (*z>=m && *z<=M)
+      return 1./(M-m);
+    else
+      return 0;
+
+  }
+
+  return (TMath::Erfc((z[0]-M)/par[2])-TMath::Erfc((z[0]-m)/par[2]))*0.5/(M-m);
+
+}
+
+
+// basic erf function
+Double_t energyCalibrationFunctions::erfFunction(Double_t *x, Double_t *par) {	
+  double arg=0;
+  if (par[1]!=0) arg=(par[0]-x[0])/par[1];
+  return ((par[2]/2.*(1+TMath::Erf(sign*arg/(TMath::Sqrt(2))))));  
+};
+  
+
+Double_t energyCalibrationFunctions::erfFunctionChargeSharing(Double_t *x, Double_t *par) {	       
+  Double_t f;		
+  
+  f=erfFunction(x, par+2)*(1+par[5]*(par[2]-x[0]))+par[0]-par[1]*x[0]*sign;		
+  return f;								
+};
+  
+  
+Double_t energyCalibrationFunctions::erfFuncFluo(Double_t *x, Double_t *par) {
+  Double_t f;							       
+  f=erfFunctionChargeSharing(x, par)+erfFunction(x, par+6)*(1+par[9]*(par[6]-x[0])); 
+  return f;							
+};
+#endif
+
+double energyCalibrationFunctions::median(double *x, int n){
+  // sorts x into xmed array and returns median 
+  // n is number of values already in the xmed array
+  double xmed[n];
+  int k,i,j;
+
+  for (i=0; i<n; i++) {
+    k=0;
+    for (j=0; j<n; j++) {
+      if(*(x+i)>*(x+j))
+	k++;
+      if (*(x+i)==*(x+j)) {
+	if (i>j) 
+	  k++;
+      }
+    }
+    xmed[k]=*(x+i);
+  }
+  k=n/2;
+  return xmed[k];
+}
+
+
+int energyCalibrationFunctions::quick_select(int arr[], int n){
+    int low, high ;
+    int median;
+    int middle, ll, hh;
+
+    low = 0 ; high = n-1 ; median = (low + high) / 2;
+    for (;;) {
+        if (high <= low) /* One element only */
+            return arr[median] ;
+
+        if (high == low + 1) {  /* Two elements only */
+            if (arr[low] > arr[high])
+                ELEM_SWAP(arr[low], arr[high]) ;
+            return arr[median] ;
+        }
+
+    /* Find median of low, middle and high items; swap into position low */
+    middle = (low + high) / 2;
+    if (arr[middle] > arr[high])    ELEM_SWAP(arr[middle], arr[high]) ;
+    if (arr[low] > arr[high])       ELEM_SWAP(arr[low], arr[high]) ;
+    if (arr[middle] > arr[low])     ELEM_SWAP(arr[middle], arr[low]) ;
+
+    /* Swap low item (now in position middle) into position (low+1) */
+    ELEM_SWAP(arr[middle], arr[low+1]) ;
+
+    /* Nibble from each end towards middle, swapping items when stuck */
+    ll = low + 1;
+    hh = high;
+    for (;;) {
+        do ll++; while (arr[low] > arr[ll]) ;
+        do hh--; while (arr[hh]  > arr[low]) ;
+
+        if (hh < ll)
+        break;
+
+        ELEM_SWAP(arr[ll], arr[hh]) ;
+    }
+
+    /* Swap middle item (in position low) back into correct position */
+    ELEM_SWAP(arr[low], arr[hh]) ;
+
+    /* Re-set active partition */
+    if (hh <= median)
+        low = ll;
+        if (hh >= median)
+        high = hh - 1;
+    }
+}
+
+int energyCalibrationFunctions::kth_smallest(int *a, int n, int k){
+    register int i,j,l,m ;
+    register double x ;
+
+    l=0 ; m=n-1 ;
+    while (l<m) {
+        x=a[k] ;
+        i=l ;
+        j=m ;
+        do {
+            while (a[i]<x) i++ ;
+            while (x<a[j]) j-- ;
+            if (i<=j) {
+                ELEM_SWAP(a[i],a[j]) ;
+                i++ ; j-- ;
+            }
+        } while (i<=j) ;
+        if (j<k) l=i ;
+        if (k<i) m=j ;
+    }
+    return a[k] ;
+}
+
+
+
+#ifdef MYROOT
+Double_t energyCalibrationFunctions::spectrum(Double_t *x, Double_t *par) {
+  return gaussChargeSharing(x,par);	
+}
+
+Double_t energyCalibrationFunctions::spectrumkb(Double_t *x, Double_t *par) {
+  return gaussChargeSharingKb(x,par);	
+}
+
+Double_t energyCalibrationFunctions::spectrumkadoublet(Double_t *x, Double_t *par) {
+  return gaussChargeSharingKaDoublet(x,par);	
+}
+
+Double_t energyCalibrationFunctions::spectrumPixel(Double_t *x, Double_t *par) {
+  return gaussChargeSharingPixel(x,par);	
+}
+
+
+Double_t energyCalibrationFunctions::scurve(Double_t *x, Double_t *par) {
+  return erfFunctionChargeSharing(x,par);
+}
+
+
+Double_t energyCalibrationFunctions::scurveFluo(Double_t *x, Double_t *par) {
+  return erfFuncFluo(x,par);						
+}
+#endif
+
+energyCalibration::energyCalibration() : 
+#ifdef MYROOT	
+					 fit_min(-1),	
+					 fit_max(-1),	
+					 bg_offset(-1),	
+					 bg_slope(-1),	
+					 flex(-1),
+					 noise(-1),
+					 ampl(-1),
+					 cs_slope(-1),
+					 kb_mean(-1),
+					 kb_frac(-1),
+					 mean2(-1),
+					 ampl2(-1),
+					 fscurve(NULL),
+					 fspectrum(NULL),
+					 fspectrumkb(NULL),
+					 fspectrumkadoublet(NULL),
+#endif
+					 funcs(NULL),
+					 plot_flag(1), // fit parameters output to screen
+					 cs_flag(1)
+{
+
+#ifdef MYROOT	
+  funcs=new energyCalibrationFunctions();
+  
+  fscurve=new TF1("fscurve",funcs,&energyCalibrationFunctions::scurve,0,1000,6,"energyCalibrationFunctions","scurve");	
+  fscurve->SetParNames("Background Offset","Background Slope","Inflection Point","Noise RMS", "Number of Photons","Charge Sharing Slope");
+
+  fspectrum=new TF1("fspectrum",funcs,&energyCalibrationFunctions::spectrum,0,1000,6,"energyCalibrationFunctions","spectrum");	
+  fspectrum->SetParNames("Background Pedestal","Background slope", "Peak position","Noise RMS", "Number of Photons","Charge Sharing Pedestal");
+  fspectrumkb=new TF1("fspectrumkb",funcs,&energyCalibrationFunctions::spectrumkb,0,1000,8,"energyCalibrationFunctions","spectrumkb");	
+  fspectrumkb->SetParNames("Background Pedestal","Background slope", "Peak position","Noise RMS", "Number of Photons","Charge Sharing Pedestal","kb mean","kb frac");
+
+  fspectrumkadoublet=new TF1("fspectrumkadoublet",funcs,&energyCalibrationFunctions::spectrumkadoublet,0,1000,8,"energyCalibrationFunctions","spectrumkadoublet");	
+  fspectrumkadoublet->SetParNames("Background Pedestal","Background slope", "Peak position","Noise RMS", "Number of Photons","Charge Sharing Pedestal","ka2 mean","n2");
+
+  fspixel=new TF1("fspixel",funcs,&energyCalibrationFunctions::spectrumPixel,0,1000,7,"energyCalibrationFunctions","spectrumPixel");	
+  fspixel->SetParNames("Background Pedestal","Background slope", "Peak position","Noise RMS", "Number of Photons","Charge Sharing Pedestal","Corner");
+
+#endif
+
+
+}
+
+
+
+void energyCalibration::fixParameter(int ip, Double_t val){
+  
+  fscurve->FixParameter(ip, val);
+  fspectrum->FixParameter(ip, val);
+  fspectrumkb->FixParameter(ip, val);
+  fspectrumkadoublet->FixParameter(ip, val);
+}
+
+
+void energyCalibration::releaseParameter(int ip){
+
+  fscurve->ReleaseParameter(ip);
+  fspectrum->ReleaseParameter(ip);
+  fspectrumkb->ReleaseParameter(ip);
+  fspectrumkadoublet->ReleaseParameter(ip);
+}
+
+
+
+
+
+
+
+energyCalibration::~energyCalibration(){ 
+#ifdef MYROOT
+  delete fscurve;
+  delete fspectrum;
+  delete fspectrumkb;
+  delete fspectrumkadoublet;
+#endif
+  
+}
+
+#ifdef MYROOT
+
+
+
+
+
+TH1F* energyCalibration::createMedianHistogram(TH2F* h2, int ch0, int nch, int direction) {
+
+  if (h2==NULL || nch==0)
+    return NULL;
+
+  double *x=new double[nch];
+  TH1F *h1=NULL;
+  
+  double val=-1;
+
+  if (direction==0) {
+    h1=new TH1F("median","Median",h2->GetYaxis()->GetNbins(),h2->GetYaxis()->GetXmin(),h2->GetYaxis()->GetXmax());
+    for (int ib=0; ib<h1->GetXaxis()->GetNbins(); ib++) {
+      for (int ich=0; ich<nch; ich++) {
+	x[ich]=h2->GetBinContent(ch0+ich+1,ib+1);
+      }
+      val=energyCalibrationFunctions::median(x, nch);
+      h1->SetBinContent(ib+1,val);
+    }
+  } else if (direction==1) {   
+    h1=new TH1F("median","Median",h2->GetXaxis()->GetNbins(),h2->GetXaxis()->GetXmin(),h2->GetXaxis()->GetXmax());
+    for (int ib=0; ib<h1->GetYaxis()->GetNbins(); ib++) {
+      for (int ich=0; ich<nch; ich++) {
+	x[ich]=h2->GetBinContent(ib+1,ch0+ich+1);
+      }
+      val=energyCalibrationFunctions::median(x, nch);
+      h1->SetBinContent(ib+1,val);
+    }
+  }
+  delete [] x;
+
+  return h1;
+
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+void energyCalibration::setStartParameters(Double_t *par){	
+  bg_offset=par[0];						
+  bg_slope=par[1];						
+  flex=par[2];							
+  noise=par[3];							
+  ampl=par[4];							
+  cs_slope=par[5];						
+}
+
+void energyCalibration::setStartParametersKb(Double_t *par){	
+  bg_offset=par[0];						
+  bg_slope=par[1];						
+  flex=par[2];							
+  noise=par[3];							
+  ampl=par[4];							
+  cs_slope=par[5];						
+  kb_mean=par[6];
+  kb_frac=par[7];
+  //fit_min = 400; // used for soleil flat field
+  //fit_max = 800;
+}
+
+void energyCalibration::setStartParametersKaDoublet(Double_t *par){	
+  bg_offset=par[0];						
+  bg_slope=par[1];						
+  flex=par[2];							
+  noise=par[3];							
+  ampl=par[4];							
+  cs_slope=par[5];						
+  mean2=par[6];
+  ampl2=par[7];
+  //fit_min = 400; // used for soleil flat field
+  //fit_max = 800;
+}
+
+
+void energyCalibration::getStartParameters(Double_t *par){	
+  par[0]=bg_offset;
+  par[1]=bg_slope;
+  par[2]=flex;
+  par[3]=noise;
+  par[4]=ampl;				
+  par[5]=cs_slope;				
+}
+
+#endif
+int energyCalibration::setChargeSharing(int p) {
+  if (p>=0) {
+    cs_flag=p; 
+#ifdef MYROOT  
+    if (p) {
+      fscurve->ReleaseParameter(5);
+      fspectrum->ReleaseParameter(1);
+      fspectrumkb->ReleaseParameter(1);
+      fspectrumkadoublet->ReleaseParameter(1);
+    } else {
+      fscurve->FixParameter(5,0);
+      fspectrum->FixParameter(1,0);
+      fspectrumkb->FixParameter(1,0);
+      fspectrumkadoublet->FixParameter(1,0);
+    }
+#endif
+  } 
+  
+  return cs_flag;
+}
+
+ 
+#ifdef MYROOT  
+void energyCalibration::initFitFunction(TF1 *fun, TH1 *h1) {
+
+ Double_t min=fit_min, max=fit_max; 
+
+ Double_t mypar[6];
+  
+ if (max==-1) 
+    max=h1->GetXaxis()->GetXmax();
+ 
+  if (min==-1) 
+    min=h1->GetXaxis()->GetXmin();
+
+  
+  if (bg_offset==-1)
+    mypar[0]=0;
+  else
+    mypar[0]=bg_offset;
+
+
+  if (bg_slope==-1)
+    mypar[1]=0;
+  else
+    mypar[1]=bg_slope;
+
+
+  if (flex==-1)
+    mypar[2]=(min+max)/2.;
+  else
+    mypar[2]=flex;
+
+
+  if (noise==-1)
+    mypar[3]=0.1;
+  else
+    mypar[3]=noise;
+
+  if (ampl==-1)
+    mypar[4]=h1->GetBinContent(h1->GetXaxis()->FindBin(0.5*(max+min)));
+  else
+    mypar[4]=ampl;
+
+  if (cs_slope==-1)
+    mypar[5]=0;
+  else
+    mypar[5]=cs_slope;
+ 
+  fun->SetParameters(mypar);
+
+  fun->SetRange(min,max);
+
+}
+
+void energyCalibration::initFitFunctionKb(TF1 *fun, TH1 *h1) {
+
+ Double_t min=fit_min, max=fit_max; 
+
+ Double_t mypar[8];
+  
+ if (max==-1) 
+    max=h1->GetXaxis()->GetXmax();
+ 
+  if (min==-1) 
+    min=h1->GetXaxis()->GetXmin();
+
+  
+  if (bg_offset==-1)
+    mypar[0]=0;
+  else
+    mypar[0]=bg_offset;
+
+
+  if (bg_slope==-1)
+    mypar[1]=0;
+  else
+    mypar[1]=bg_slope;
+
+
+  if (flex==-1)
+    mypar[2]=(min+max)/2.;
+  else
+    mypar[2]=flex;
+
+
+  if (noise==-1)
+    mypar[3]=0.1;
+  else
+    mypar[3]=noise;
+
+  if (ampl==-1)
+    mypar[4]=h1->GetBinContent(h1->GetXaxis()->FindBin(0.5*(max+min)));
+  else
+    mypar[4]=ampl;
+
+  if (cs_slope==-1)
+    mypar[5]=0;
+  else
+    mypar[5]=cs_slope;
+
+  if (kb_mean==-1)
+    mypar[6]=0;
+  else
+    mypar[6]=kb_mean;
+
+  if (kb_frac==-1)
+    mypar[7]=0;
+  else
+    mypar[7]=kb_frac;
+ 
+  fun->SetParameters(mypar);
+
+  fun->SetRange(min,max);
+
+}
+
+void energyCalibration::initFitFunctionKaDoublet(TF1 *fun, TH1 *h1) {
+
+ Double_t min=fit_min, max=fit_max; 
+
+ Double_t mypar[8];
+  
+ if (max==-1) 
+    max=h1->GetXaxis()->GetXmax();
+ 
+  if (min==-1) 
+    min=h1->GetXaxis()->GetXmin();
+
+  
+  if (bg_offset==-1)
+    mypar[0]=0;
+  else
+    mypar[0]=bg_offset;
+
+
+  if (bg_slope==-1)
+    mypar[1]=0;
+  else
+    mypar[1]=bg_slope;
+
+
+  if (flex==-1)
+    mypar[2]=(min+max)/2.;
+  else
+    mypar[2]=flex;
+
+
+  if (noise==-1)
+    mypar[3]=0.1;
+  else
+    mypar[3]=noise;
+
+  if (ampl==-1)
+    mypar[4]=h1->GetBinContent(h1->GetXaxis()->FindBin(0.5*(max+min)));
+  else
+    mypar[4]=ampl;
+
+  if (cs_slope==-1)
+    mypar[5]=0;
+  else
+    mypar[5]=cs_slope;
+
+  if (mean2==-1)
+    mypar[6]=0;
+  else
+    mypar[6]=mean2;
+
+  if (ampl2==-1)
+    mypar[7]=0;
+  else
+    mypar[7]=ampl2;
+ 
+  fun->SetParameters(mypar);
+
+  fun->SetRange(min,max);
+
+}
+
+TF1* energyCalibration::fitFunction(TF1 *fun, TH1 *h1, Double_t *mypar, Double_t *emypar) {
+
+
+  TF1* fitfun;
+
+  char fname[100];
+
+  strcpy(fname, fun->GetName());
+
+  if (plot_flag) {
+    h1->Fit(fname,"R0Q");
+  } else
+    h1->Fit(fname,"R0Q");
+
+  fitfun= h1->GetFunction(fname);
+  fitfun->GetParameters(mypar);
+  for (int ip=0; ip<6; ip++) {
+    emypar[ip]=fitfun->GetParError(ip);
+  }
+  return fitfun;
+}
+
+TF1* energyCalibration::fitFunctionKb(TF1 *fun, TH1 *h1, Double_t *mypar, Double_t *emypar) {
+
+
+  TF1* fitfun;
+
+  char fname[100];
+
+  strcpy(fname, fun->GetName());
+
+  if (plot_flag) {
+    h1->Fit(fname,"R0Q");
+  } else
+    h1->Fit(fname,"R0Q");
+
+  fitfun= h1->GetFunction(fname);
+  fitfun->GetParameters(mypar);
+  for (int ip=0; ip<8; ip++) {
+    emypar[ip]=fitfun->GetParError(ip);
+  }
+  return fitfun;
+}
+
+TF1* energyCalibration::fitFunctionKaDoublet(TF1 *fun, TH1 *h1, Double_t *mypar, Double_t *emypar) {
+
+
+  TF1* fitfun;
+
+  char fname[100];
+
+  strcpy(fname, fun->GetName());
+
+  if (plot_flag) {
+    h1->Fit(fname,"R0Q");
+  } else
+    h1->Fit(fname,"R0Q");
+
+
+  fitfun= h1->GetFunction(fname);
+  fitfun->GetParameters(mypar);
+  for (int ip=0; ip<8; ip++) {
+    emypar[ip]=fitfun->GetParError(ip);
+  }
+  return fitfun;
+}
+
+TF1* energyCalibration::fitSCurve(TH1 *h1, Double_t *mypar, Double_t *emypar) {
+  initFitFunction(fscurve,h1);
+  return fitFunction(fscurve, h1, mypar, emypar);
+}
+
+
+
+
+
+TF1* energyCalibration::fitSpectrum(TH1 *h1, Double_t *mypar, Double_t *emypar) {
+  initFitFunction(fspectrum,h1);
+  return fitFunction(fspectrum, h1, mypar, emypar);
+}
+
+TF1* energyCalibration::fitSpectrumKb(TH1 *h1, Double_t *mypar, Double_t *emypar) {
+  initFitFunctionKb(fspectrumkb,h1);
+  return fitFunctionKb(fspectrumkb, h1, mypar, emypar);
+}
+
+TF1* energyCalibration::fitSpectrumKaDoublet(TH1 *h1, Double_t *mypar, Double_t *emypar) {
+  initFitFunctionKaDoublet(fspectrumkadoublet,h1);
+  return fitFunctionKaDoublet(fspectrumkadoublet, h1, mypar, emypar);
+}
+
+
+TGraphErrors* energyCalibration::linearCalibration(int nscan, Double_t *en, Double_t *een, Double_t *fl, Double_t *efl, Double_t &gain, Double_t &off, Double_t &egain, Double_t &eoff) { 
+ 
+  TGraphErrors *gr;
+ 
+  Double_t  mypar[2];
+ 
+  gr = new TGraphErrors(nscan,en,fl,een,efl);
+
+  if (plot_flag) {
+    gr->Fit("pol1");
+    gr->SetMarkerStyle(20);
+  } else
+    gr->Fit("pol1","0Q");
+        
+  TF1 *fitfun= gr->GetFunction("pol1"); 
+  fitfun->GetParameters(mypar);
+
+  egain=fitfun->GetParError(1);
+  eoff=fitfun->GetParError(0);
+
+  gain=funcs->setScanSign()*mypar[1];
+  
+  off=mypar[0]; 
+  
+  return gr;
+}
+
+
+TGraphErrors* energyCalibration::calibrate(int nscan, Double_t *en, Double_t *een, TH1F **h1, Double_t &gain, Double_t &off, Double_t &egain, Double_t &eoff, int integral) {
+
+  TH1F *h;
+
+  Double_t mypar[6], emypar[6];
+  Double_t fl[nscan], efl[nscan];
+  
+
+  for (int ien=0; ien<nscan; ien++) {
+    h=h1[ien];
+    if (integral)
+      fitSCurve(h,mypar,emypar);
+    else
+      fitSpectrum(h,mypar,emypar);
+
+    fl[ien]=mypar[2];
+    efl[ien]=emypar[2];
+  }
+  return linearCalibration(nscan,en,een,fl,efl,gain,off, egain, eoff);
+
+}
+
+#endif
+
+
+
+

slsDetectorCalibration/energyCalibration.h

@@ -98,6 +98,8 @@ class energyCalibrationFunctions {
       par[5] is the fractional height of the charge sharing pedestal (scales with par[3])
   */
   Double_t gaussChargeSharing(Double_t *x, Double_t *par);
+  Double_t gaussChargeSharingKb(Double_t *x, Double_t *par);
+  Double_t gaussChargeSharingKaDoublet(Double_t *x, Double_t *par);
   /** 
       Gaussian Function with charge sharing pedestal
       par[0] is the absolute height of the background pedestal
@@ -153,6 +155,8 @@ Double_t erfFuncFluo(Double_t *x, Double_t *par);
       par[5] is the fractional height of the charge sharing pedestal (scales with par[4]
   */
   Double_t spectrum(Double_t *x, Double_t *par);
+  Double_t spectrumkb(Double_t *x, Double_t *par);
+  Double_t spectrumkadoublet(Double_t *x, Double_t *par);
 
   /** 
       static function Gaussian with charge sharing pedestal with the correct scan sign
@@ -285,7 +289,9 @@ class energyCalibration  {
       par[5] is the angual coefficient of the charge sharing slope (scales with par[3]) -- always positive
   */
   void setStartParameters(Double_t *par);
-  
+  void setStartParametersKb(Double_t *par);
+  void setStartParametersKaDoublet(Double_t *par);
+
 /** get start parameters for the s-curve function
     \param par parameters, -1 means auto-calculated
       par[0] is the  pedestal 
@@ -315,16 +321,8 @@ class energyCalibration  {
      \returns the fitted function - can be used e.g. to get the Chi2 or similar
   */
   TF1 *fitSpectrum(TH1 *h1, Double_t *mypar, Double_t *emypar);
-
-
-  /**
-     fits histogram with the spectrum
-     \param h1 1d-histogram to be fitted
-     \param mypar pointer to fit parameters array 
-     \param emypar pointer to fit parameter errors
-     \returns the fitted function - can be used e.g. to get the Chi2 or similar
-  */
-  TF1 *fitSpectrumPixel(TH1 *h1, Double_t *mypar, Double_t *emypar);
+  TF1 *fitSpectrumKb(TH1 *h1, Double_t *mypar, Double_t *emypar);
+  TF1 *fitSpectrumKaDoublet(TH1 *h1, Double_t *mypar, Double_t *emypar);
 
 
   /**
@@ -399,6 +397,8 @@ class energyCalibration  {
 */
   
   void initFitFunction(TF1 *fun, TH1 *h1);
+  void initFitFunctionKb(TF1 *fun, TH1 *h1);
+  void initFitFunctionKaDoublet(TF1 *fun, TH1 *h1);
 
 
   /**
@@ -410,6 +410,8 @@ class energyCalibration  {
      \returns the fitted function - can be used e.g. to get the Chi2 or similar
   */
   TF1 *fitFunction(TF1 *fun, TH1 *h1, Double_t *mypar, Double_t *emypar);  
+  TF1 *fitFunctionKb(TF1 *fun, TH1 *h1, Double_t *mypar, Double_t *emypar);  
+  TF1 *fitFunctionKaDoublet(TF1 *fun, TH1 *h1, Double_t *mypar, Double_t *emypar);  
 
 #endif
 
@@ -423,11 +425,16 @@ class energyCalibration  {
   Double_t noise; /**< start value for the noise */
   Double_t ampl; /**< start value for the number of photons */
   Double_t cs_slope; /**< start value for the charge sharing slope */
-
+  Double_t kb_mean;
+  Double_t kb_frac;
+  Double_t mean2;
+  Double_t ampl2;
 
   TF1 *fscurve; /**< function with which the s-curve will be fitted */
 
   TF1 *fspectrum; /**< function with which the spectrum will be fitted */
+  TF1 *fspectrumkb; /**< function with which the spectrum will be fitted */
+  TF1 *fspectrumkadoublet; /**< function with which the spectrum will be fitted */
 
   TF1 *fspixel; /**< function with which the spectrum will be fitted */
 

slsDetectorCalibration/jungfrauExecutables/Makefile.cluster_finder

@@ -0,0 +1,47 @@
+#module add CBFlib/0.9.5
+INCDIR=-I. -I../ -I../interpolations  -I../interpolations/etaVEL  -I../dataStructures  -I../../slsSupportLib/include/ -I../../slsReceiverSoftware/include/ 
+
+LDFLAG= ../tiffIO.cpp -L/usr/lib64/ -lpthread -lm  -lstdc++   -pthread -lrt    -ltiff -O3  -std=c++11
+
+MAIN=jungfrauClusterFinder.cpp
+
+
+all: jungfrauClusterFinder jungfrauMakeEta jungfrauInterpolation jungfrauNoInterpolation jungfrauPhotonCounter jungfrauAnalog
+
+
+
+jungfrauClusterFinder:  jungfrauClusterFinder.cpp  $(INCS) clean
+			 g++ -o jungfrauClusterFinder  jungfrauClusterFinder.cpp $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF) -DSAVE_ALL  
+
+
+jungfrauClusterFinderHighZ:  jungfrauClusterFinder.cpp  $(INCS) clean
+			 g++ -o jungfrauClusterFinderHighZ  jungfrauClusterFinder.cpp $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF) -DSAVE_ALL  -DHIGHZ 
+
+
+
+
+jungfrauMakeEta:  jungfrauInterpolation.cpp  $(INCS) clean
+			 g++ -o jungfrauMakeEta jungfrauInterpolation.cpp    $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF) -DFF
+
+jungfrauInterpolation:  jungfrauInterpolation.cpp  $(INCS) clean
+			 g++ -o jungfrauInterpolation jungfrauInterpolation.cpp    $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF)  
+
+jungfrauNoInterpolation:  jungfrauNoInterpolation.cpp  $(INCS) clean
+			 g++ -o jungfrauNoInterpolation jungfrauNoInterpolation.cpp    $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF) 
+
+jungfrauPhotonCounter: jungfrauPhotonCounter.cpp   $(INCS) clean
+			 g++ -o jungfrauPhotonCounter  jungfrauPhotonCounter.cpp  $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF)  -DNEWRECEIVER
+
+jungfrauAnalog: jungfrauPhotonCounter.cpp  $(INCS) clean
+			 g++ -o jungfrauAnalog  jungfrauPhotonCounter.cpp  $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF)  -DNEWRECEIVER -DANALOG
+
+jungfrauPhotonCounterHighZ: jungfrauPhotonCounter.cpp   $(INCS) clean
+			 g++ -o jungfrauPhotonCounterHighZ  jungfrauPhotonCounter.cpp  $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF)  -DNEWRECEIVER -DHIGHZ
+
+jungfrauAnalogHighZ: jungfrauPhotonCounter.cpp  $(INCS) clean
+			 g++ -o jungfrauAnalogHighZ  jungfrauPhotonCounter.cpp  $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF)  -DNEWRECEIVER -DANALOG -DHIGHZ
+
+clean: 	
+	rm -f  jungfrauClusterFinder jungfrauMakeEta jungfrauInterpolation jungfrauNoInterpolation jungfrauPhotonCounter jungfrauAnalog
+
+

slsDetectorCalibration/jungfrauExecutables/Makefile.cluster_finder~

@@ -0,0 +1,47 @@
+#module add CBFlib/0.9.5
+INCDIR=-I. -I../ -I../interpolations  -I../interpolations/etaVEL  -I../dataStructures  -I../../slsSupportLib/include/ -I../../slsReceiverSoftware/include/
+
+LDFLAG= ../tiffIO.cpp -L/usr/lib64/ -lpthread -lm  -lstdc++   -pthread -lrt    -ltiff -O3  -std=c++11
+
+MAIN=moench03ClusterFinder.cpp
+
+
+all: moenchClusterFinder moenchMakeEta moenchInterpolation moenchNoInterpolation moenchPhotonCounter moenchAnalog
+
+
+
+moenchClusterFinder:  moench03ClusterFinder.cpp  $(INCS) clean
+			 g++ -o moenchClusterFinder  moench03ClusterFinder.cpp $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF) -DSAVE_ALL  -DNEWRECEIVER 
+
+
+moenchClusterFinderHighZ:  moench03ClusterFinder.cpp  $(INCS) clean
+			 g++ -o moenchClusterFinderHighZ  moench03ClusterFinder.cpp $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF) -DSAVE_ALL  -DNEWRECEIVER -DHIGHZ 
+
+
+
+
+moenchMakeEta:  moench03Interpolation.cpp  $(INCS) clean
+			 g++ -o moenchMakeEta moench03Interpolation.cpp    $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF) -DFF
+
+moenchInterpolation:  moench03Interpolation.cpp  $(INCS) clean
+			 g++ -o moenchInterpolation moench03Interpolation.cpp    $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF)  
+
+moenchNoInterpolation:  moench03NoInterpolation.cpp  $(INCS) clean
+			 g++ -o moenchNoInterpolation moench03NoInterpolation.cpp    $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF) 
+
+moenchPhotonCounter: moenchPhotonCounter.cpp   $(INCS) clean
+			 g++ -o moenchPhotonCounter  moenchPhotonCounter.cpp  $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF)  -DNEWRECEIVER
+
+moenchAnalog: moenchPhotonCounter.cpp  $(INCS) clean
+			 g++ -o moenchAnalog  moenchPhotonCounter.cpp  $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF)  -DNEWRECEIVER -DANALOG
+
+moenchPhotonCounterHighZ: moenchPhotonCounter.cpp   $(INCS) clean
+			 g++ -o moenchPhotonCounterHighZ  moenchPhotonCounter.cpp  $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF)  -DNEWRECEIVER -DHIGHZ
+
+moenchAnalogHighZ: moenchPhotonCounter.cpp  $(INCS) clean
+			 g++ -o moenchAnalogHighZ  moenchPhotonCounter.cpp  $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF)  -DNEWRECEIVER -DANALOG -DHIGHZ
+
+clean: 	
+	rm -f  moenchClusterFinder moenchMakeEta moenchInterpolation moenchNoInterpolation moenchPhotonCounter moenchAnalog
+
+

slsDetectorCalibration/jungfrauExecutables/Makefile.zmq

@@ -0,0 +1,23 @@
+
+INCDIR=   -I. -I../dataStructures ../tiffIO.cpp -I../ -I../interpolations/ -I../../slsSupportLib/include/ -I../../slsReceiverSoftware/include/ -I../../libs/rapidjson/
+LDFLAG= -L/usr/lib64/ -lpthread -lm  -lstdc++  -lzmq -pthread -lrt -ltiff   -O3  -std=c++11 -Wall -L../../build/bin/ -lSlsSupport
+#-L../../bin -lhdf5  -L.
+
+#DESTDIR?=../bin
+
+all:  moenchZmqProcess moenchZmq04Process
+ #moenchZmqProcessCtbGui
+
+moenchZmqProcess:  moenchZmqProcess.cpp   clean
+		      g++ -o moenchZmqProcess  moenchZmqProcess.cpp  $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF)  -DNEWZMQ -DINTERP 
+
+moenchZmq04Process:  moenchZmqProcess.cpp   clean
+		      g++ -o moench04ZmqProcess  moenchZmqProcess.cpp  $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF)  -DNEWZMQ -DINTERP -DMOENCH04 
+
+#moenchZmqProcessCtbGui:  moenchZmqProcess.cpp   clean
+#		      g++ -o moenchZmqProcessCtbGui  moenchZmqProcess.cpp  $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF)  -DNEWZMQ -DINTERP -DCTBGUI
+
+clean: 	
+	rm -f  moenchZmqProcess  
+
+

slsDetectorCalibration/jungfrauExecutables/jungfrauClusterFinder.cpp

@@ -0,0 +1,161 @@
+//#include "sls/ansi.h"
+#include <iostream>
+
+
+//#include "moench03T1ZmqData.h"
+#include "jungfrauHighZSingleChipData.h"
+
+
+#include "multiThreadedAnalogDetector.h"
+#include "singlePhotonDetector.h"
+
+#include <stdio.h>
+#include <map>
+#include <fstream>
+#include <sys/stat.h>
+
+#include <ctime>
+using namespace std;
+
+
+int main(int argc, char *argv[]) {
+
+
+  if (argc<6) {
+    cout << "Usage is " << argv[0] << "indir outdir fname runmin runmax " << endl;
+    return 1;
+  }
+  int p=10000;
+  int fifosize=1000;
+  int nthreads=1;
+  int nsubpix=25;
+  int etabins=nsubpix*10;
+  double etamin=-1, etamax=2;
+  int csize=3;
+  int nx=400, ny=400;
+  int save=1;
+  int nsigma=5;
+  int nped=1000;
+  int ndark=100;
+  int ok;
+  int iprog=0;
+
+
+
+
+  jungfrauHighZSingleChipData *decoder=new  jungfrauHighZSingleChipData();
+
+  decoder->getDetectorSize(nx,ny);
+  cout << "nx " << nx << " ny " << ny << endl;
+
+  //moench03T1ZmqData *decoder=new  moench03T1ZmqData();
+  singlePhotonDetector *filter=new singlePhotonDetector(decoder,csize, nsigma, 1, 0, nped, 200);
+  //  char tit[10000];
+  cout << "filter " << endl;
+
+
+
+
+  int* image;
+  filter->newDataSet();
+
+
+  int ff, np;
+  int dsize=decoder->getDataSize();
+  cout << " data size is " << dsize;
+  
+
+  char data[dsize];
+
+  ifstream filebin;
+  char *indir=argv[1];
+  char *outdir=argv[2];
+  char *fformat=argv[3];
+  int runmin=atoi(argv[4]);
+  int runmax=atoi(argv[5]);
+  
+  char fname[10000];
+  char outfname[10000];
+  char imgfname[10000];
+  char pedfname[10000];
+  char fn[10000];
+  
+  std::time_t end_time;
+
+  FILE *of=NULL;
+  cout << "input directory is " << indir << endl;
+  cout << "output directory is " << outdir << endl;
+  cout << "fileformat is " << fformat << endl;
+  
+
+  std::time(&end_time);
+  cout << std::ctime(&end_time) <<   endl;
+
+  char* buff;
+  multiThreadedAnalogDetector *mt=new multiThreadedAnalogDetector(filter,nthreads,fifosize);
+
+ 
+  mt->setDetectorMode(ePhotonCounting);
+  mt->setFrameMode(eFrame);
+  mt->StartThreads();
+  mt->popFree(buff);
+
+
+  cout << "mt " << endl;
+
+  int ifr=0;
+ 
+
+  for (int irun=runmin; irun<runmax; irun++) {
+    sprintf(fn,fformat,irun);
+    sprintf(fname,"%s/%s.raw",indir,fn);
+    sprintf(outfname,"%s/%s.clust",outdir,fn);
+    sprintf(imgfname,"%s/%s.tiff",outdir,fn);
+    std::time(&end_time);
+    cout << std::ctime(&end_time) <<    endl;
+    cout <<  fname << " " << outfname << " " << imgfname <<  endl;
+    filebin.open((const char *)(fname), ios::in | ios::binary);
+    //      //open file
+    if (filebin.is_open()){
+      of=fopen(outfname,"w");
+      if (of) {
+  	mt->setFilePointer(of);
+	//	cout << "file pointer set " << endl;
+      } else {
+  	cout << "Could not open "<< outfname << " for writing " << endl;
+  	mt->setFilePointer(NULL);
+  	return 1;
+      }
+      //     //while read frame 
+      ff=-1;
+      while (decoder->readNextFrame(filebin, ff, np,buff)) {
+
+  	mt->pushData(buff);
+  	mt->nextThread();
+  	mt->popFree(buff);
+	ifr++;
+	if (ifr%10000==0) cout << ifr << " " << ff << endl;
+	ff=-1;
+      }
+        cout << "--" << endl;
+      filebin.close();	 
+      while (mt->isBusy()) {;}//wait until all data are processed from the queues
+      if (of)
+  	fclose(of);
+      
+      mt->writeImage(imgfname);
+      mt->clearImage();
+   
+      std::time(&end_time);
+      cout << std::ctime(&end_time) <<   endl;
+
+    } else 
+     cout << "Could not open "<< fname << " for reading " << endl;
+      
+    
+  }
+    
+
+  return 0;
+}
+

slsDetectorCalibration/jungfrauExecutables/jungfrauClusterFinder.cpp~

@@ -0,0 +1,237 @@
+//#include "sls/ansi.h"
+#include <iostream>
+
+
+//#include "moench03T1ZmqData.h"
+#ifdef NEWRECEIVER
+#ifndef RECT
+#include "moench03T1ReceiverDataNew.h"
+#endif
+
+#ifdef RECT
+#include "moench03T1ReceiverDataNewRect.h"
+#endif
+
+#endif
+
+
+#ifdef CSAXS_FP
+#include "moench03T1ReceiverData.h"
+#endif 
+#ifdef OLDDATA
+#include "moench03Ctb10GbT1Data.h"
+#endif 
+
+#ifdef REORDERED
+#include "moench03T1ReorderedData.h"
+#endif
+
+// #include "interpolatingDetector.h"
+//#include "etaInterpolationPosXY.h"
+// #include "linearInterpolation.h"
+// #include "noInterpolation.h"
+#include "multiThreadedAnalogDetector.h"
+#include "singlePhotonDetector.h"
+//#include "interpolatingDetector.h"
+
+#include <stdio.h>
+#include <map>
+#include <fstream>
+#include <sys/stat.h>
+
+#include <ctime>
+using namespace std;
+
+
+int main(int argc, char *argv[]) {
+
+
+  if (argc<6) {
+    cout << "Usage is " << argv[0] << "indir outdir fname runmin runmax " << endl;
+    return 1;
+  }
+  int p=10000;
+  int fifosize=1000;
+  int nthreads=1;
+  int nsubpix=25;
+  int etabins=nsubpix*10;
+  double etamin=-1, etamax=2;
+  int csize=3;
+  int nx=400, ny=400;
+  int save=1;
+  int nsigma=5;
+  int nped=1000;
+  int ndark=100;
+  int ok;
+  int iprog=0;
+
+
+
+
+
+#ifdef NEWRECEIVER
+#ifdef RECT
+  cout << "Should be rectangular!" <<endl;
+#endif
+  moench03T1ReceiverDataNew *decoder=new  moench03T1ReceiverDataNew();
+  cout << "RECEIVER DATA WITH ONE HEADER!"<<endl;
+#endif
+
+#ifdef CSAXS_FP
+  moench03T1ReceiverData *decoder=new  moench03T1ReceiverData();
+  cout << "RECEIVER DATA WITH ALL HEADERS!"<<endl;
+#endif
+
+#ifdef OLDDATA
+  moench03Ctb10GbT1Data *decoder=new  moench03Ctb10GbT1Data();
+  cout << "OLD RECEIVER DATA!"<<endl;
+#endif
+
+#ifdef REORDERED
+  moench03T1ReorderedData *decoder=new  moench03T1ReorderedData();
+  cout << "REORDERED DATA!"<<endl;
+#endif
+
+
+  decoder->getDetectorSize(nx,ny);
+  cout << "nx " << nx << " ny " << ny << endl;
+
+  //moench03T1ZmqData *decoder=new  moench03T1ZmqData();
+  singlePhotonDetector *filter=new singlePhotonDetector(decoder,csize, nsigma, 1, 0, nped, 200);
+  //  char tit[10000];
+  cout << "filter " << endl;
+
+
+
+  // filter->readPedestals("/scratch/ped_100.tiff");
+  // interp->readFlatField("/scratch/eta_100.tiff",etamin,etamax);
+  // cout << "filter "<< endl;
+  
+
+  int size = 327680;////atoi(argv[3]);
+  
+  int* image;
+	//int* image =new int[327680/sizeof(int)];
+  filter->newDataSet();
+
+
+  int ff, np;
+  int dsize=decoder->getDataSize();
+  cout << " data size is " << dsize;
+  
+
+  char data[dsize];
+
+  ifstream filebin;
+  char *indir=argv[1];
+  char *outdir=argv[2];
+  char *fformat=argv[3];
+  int runmin=atoi(argv[4]);
+  int runmax=atoi(argv[5]);
+  
+  char fname[10000];
+  char outfname[10000];
+  char imgfname[10000];
+  char pedfname[10000];
+  //  strcpy(pedfname,argv[6]);
+  char fn[10000];
+  
+  std::time_t end_time;
+
+  FILE *of=NULL;
+  cout << "input directory is " << indir << endl;
+  cout << "output directory is " << outdir << endl;
+  cout << "fileformat is " << fformat << endl;
+  
+
+  std::time(&end_time);
+  cout << std::ctime(&end_time) <<   endl;
+ 
+
+
+
+
+
+
+
+  char* buff;
+  multiThreadedAnalogDetector *mt=new multiThreadedAnalogDetector(filter,nthreads,fifosize);
+
+ 
+  mt->setDetectorMode(ePhotonCounting);
+  mt->setFrameMode(eFrame);
+  mt->StartThreads();
+  mt->popFree(buff);
+
+
+  cout << "mt " << endl;
+
+  int ifr=0;
+ 
+
+  for (int irun=runmin; irun<runmax; irun++) {
+    sprintf(fn,fformat,irun);
+    sprintf(fname,"%s/%s.raw",indir,fn);
+    sprintf(outfname,"%s/%s.clust",outdir,fn);
+    sprintf(imgfname,"%s/%s.tiff",outdir,fn);
+    std::time(&end_time);
+    cout << std::ctime(&end_time) <<    endl;
+    cout <<  fname << " " << outfname << " " << imgfname <<  endl;
+    filebin.open((const char *)(fname), ios::in | ios::binary);
+    //      //open file
+    if (filebin.is_open()){
+      of=fopen(outfname,"w");
+      if (of) {
+  	mt->setFilePointer(of);
+	//	cout << "file pointer set " << endl;
+      } else {
+  	cout << "Could not open "<< outfname << " for writing " << endl;
+  	mt->setFilePointer(NULL);
+  	return 1;
+      }
+      //     //while read frame 
+      ff=-1;
+      while (decoder->readNextFrame(filebin, ff, np,buff)) {
+	//	cout << "*"<<ifr++<<"*"<<ff<< endl;
+	//	cout << ff << " " << np << endl;
+  	//         //push
+	// for (int ix=0; ix<400; ix++)
+	//   for (int iy=0; iy<400; iy++) {
+	//     if (decoder->getChannel(buff, ix, iy)<3000 || decoder->getChannel(buff, ix, iy)>8000) {
+	//       cout <<  ifr << " " << ff << " " << ix << " " << iy << " " << decoder->getChannel(buff, ix, iy) << endl ;
+	//     }
+	//   }
+
+  	mt->pushData(buff);
+  // 	//         //pop
+  	mt->nextThread();
+  // // 		//	cout << " " << (void*)buff;
+  	mt->popFree(buff);
+	ifr++;
+	if (ifr%10000==0) cout << ifr << " " << ff << endl;
+	ff=-1;
+      }
+        cout << "--" << endl;
+      filebin.close();	 
+      //      //close file 
+      //     //join threads
+      while (mt->isBusy()) {;}//wait until all data are processed from the queues
+      if (of)
+  	fclose(of);
+      
+      mt->writeImage(imgfname);
+      mt->clearImage();
+   
+      std::time(&end_time);
+      cout << std::ctime(&end_time) <<   endl;
+
+    } else 
+     cout << "Could not open "<< fname << " for reading " << endl;
+      
+    
+  }
+    
+
+  return 0;
+}
+

slsDetectorCalibration/jungfrauExecutables/jungfrauInterpolation.cpp

@@ -0,0 +1,250 @@
+
+#include "sls/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 "etaInterpolationPosXY.h"
+//#include "etaInterpolationCleverAdaptiveBins.h"
+//#include "etaInterpolationRandomBins.h"
+using namespace std;
+#define NC 400
+#define NR 400
+#define MAX_ITERATIONS (nSubPixels*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 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;
+#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
+
+
+ 
+#ifndef FF
+#ifndef NOINTERPOLATION
+    cout << "read ff " << argv[2] << endl;
+    sprintf(fname,"%s",argv[2]);
+    interp->readFlatField(fname);
+    interp->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->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
+#ifndef FF
+	    // 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);
+	    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%10000000==0) { 
+#ifndef FF
+		interp->writeInterpolatedImage(outfname);
+#endif
+#ifdef FF 
+		interp->writeFlatField(outfname);
+#endif
+		
+		}
+	  }
+	    
+	} 
+	  
+	  
+	fclose(f);
+#ifdef FF
+	interp->writeFlatField(outfname);
+#endif
+	
+#ifndef FF
+	interp->writeInterpolatedImage(outfname); 
+
+	img=interp->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)];
+		}
+	    }
+	  }
+	}
+	cout << "Read " << nframes << " frames (first frame: " << f0 << " last frame: " << lastframe << " delta:" << lastframe-f0 << ") nph="<< nph <<endl;
+	interp->clearInterpolatedImage();
+#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->writeFlatField(outfname);
+#endif
+
+    cout << "Filled " << nph << " (/"<< totph <<") " << endl; 
+    return 0;
+}
+

slsDetectorCalibration/jungfrauExecutables/jungfrauPhotonCounter.cpp

@@ -0,0 +1,453 @@
+//#include "sls/ansi.h"
+#include <iostream>
+#define CORR
+
+#define C_GHOST 0.0004
+
+#define CM_ROWS 50
+
+//#define VERSION_V1
+
+//#include "moench03T1ZmqData.h"
+#ifdef NEWRECEIVER
+#ifndef RECT
+#include "moench03T1ReceiverDataNew.h"
+#endif
+
+#ifdef RECT
+#include "moench03T1ReceiverDataNewRect.h"
+#endif
+
+#endif
+
+
+#ifdef CSAXS_FP
+#include "moench03T1ReceiverData.h"
+#endif 
+#ifdef OLDDATA
+#include "moench03Ctb10GbT1Data.h"
+#endif 
+
+// #include "interpolatingDetector.h"
+//#include "etaInterpolationPosXY.h"
+// #include "linearInterpolation.h"
+// #include "noInterpolation.h"
+#include "multiThreadedCountingDetector.h"
+//#include "multiThreadedAnalogDetector.h"
+#include "singlePhotonDetector.h"
+#include "moench03GhostSummation.h"
+#include "moench03CommonMode.h"
+//#include "interpolatingDetector.h"
+
+#include <stdio.h>
+#include <map>
+#include <fstream>
+#include <sys/stat.h>
+
+#include <ctime>
+using namespace std;
+
+
+int main(int argc, char *argv[]) {
+
+
+  if (argc<4) {
+    cout << "Usage is " << argv[0] << "indir outdir fname [runmin] [runmax] [pedfile] [threshold] [nframes] [xmin xmax ymin ymax] [gainmap]" << endl;
+    cout << "threshold <0 means analog; threshold=0 means cluster finder; threshold>0 means photon counting" << endl;
+    cout << "nframes <0 means sum everything; nframes=0 means one file per run; nframes>0 means one file every nframes" << endl;
+    return 1;
+  }
+
+  int p=10000;
+  int fifosize=1000;
+  int nthreads=10;
+  int nsubpix=25;
+  int etabins=nsubpix*10;
+  double etamin=-1, etamax=2;
+  int csize=3;
+  int save=1;
+  int nsigma=5;
+  int nped=10000;
+  int ndark=100;
+  int ok;
+  int iprog=0;
+
+  int cf=0;
+
+#ifdef NEWRECEIVER
+#ifdef RECT
+  cout << "Should be rectangular!" <<endl;
+#endif
+  moench03T1ReceiverDataNew *decoder=new  moench03T1ReceiverDataNew();
+  cout << "RECEIVER DATA WITH ONE HEADER!"<<endl;
+#endif
+
+#ifdef CSAXS_FP
+  moench03T1ReceiverData *decoder=new  moench03T1ReceiverData();
+  cout << "RECEIVER DATA WITH ALL HEADERS!"<<endl;
+#endif
+
+#ifdef OLDDATA
+  moench03Ctb10GbT1Data *decoder=new  moench03Ctb10GbT1Data();
+  cout << "OLD RECEIVER DATA!"<<endl;
+#endif
+
+  int nx=400, ny=400;
+
+  decoder->getDetectorSize(nx,ny);
+
+  int ncol_cm=CM_ROWS;
+  double xt_ghost=C_GHOST;
+  moench03CommonMode *cm=NULL;
+  moench03GhostSummation *gs;
+  double *gainmap=NULL;
+  float *gm;
+
+
+
+  int size = 327680;////atoi(argv[3]);
+  
+  int* image;
+	//int* image =new int[327680/sizeof(int)];
+
+  int ff, np;
+  //cout << " data size is " << dsize;
+  
+
+
+  ifstream filebin;
+  char *indir=argv[1];
+  char *outdir=argv[2];
+  char *fformat=argv[3];
+  int runmin=0;
+  
+  // cout << "argc is " << argc << endl;
+ if (argc>=5) {
+   runmin=atoi(argv[4]);
+  }
+
+ int runmax=runmin;
+
+ if (argc>=6) {
+   runmax=atoi(argv[5]);
+  }
+
+  char *pedfile=NULL;
+  if (argc>=7) {
+    pedfile=argv[6];
+  }
+  double thr=0;
+  double thr1=1;
+
+  if (argc>=8) {
+    thr=atof(argv[7]);
+  }
+  
+
+  int nframes=0;
+  
+  if (argc>=9) {
+    nframes=atoi(argv[8]);
+  }
+
+  int xmin=0, xmax=nx, ymin=0, ymax=ny;
+  if (argc>=13) {
+    xmin=atoi(argv[9]);
+    xmax=atoi(argv[10]);
+    ymin=atoi(argv[11]);
+    ymax=atoi(argv[12]);
+  }
+
+  
+  char *gainfname=NULL;
+  if (argc>13) {
+    gainfname=argv[13];
+    cout << "Gain map file name is: " << gainfname << endl;
+  }
+
+  
+
+
+  char ffname[10000];
+  char fname[10000];
+  char imgfname[10000];
+  char cfname[10000];
+  char fn[10000];
+  
+  std::time_t end_time;
+
+  FILE *of=NULL;
+  cout << "input directory is " << indir << endl;
+  cout << "output directory is " << outdir << endl;
+  cout << "input file is " << fformat << endl;
+  cout << "runmin is " << runmin << endl;
+  cout << "runmax is " << runmax << endl;
+  if (pedfile)
+    cout << "pedestal file is " << pedfile << endl;
+ if (thr>0) 
+    cout << "threshold is " << thr << endl;
+ cout << "Nframes is " << nframes << endl;
+ 
+ uint32 nnx, nny;
+ double *gmap;
+
+  // if (gainfname) {
+  //   gm=ReadFromTiff(gainfname, nny, nnx);
+  //   if (gm && nnx==nx && nny==ny) {
+  //     gmap=new double[nx*ny];
+  //     for (int i=0; i<nx*ny; i++) {
+  // 	gmap[i]=gm[i];
+  //     }
+  //     delete gm;
+  //   } else
+  //     cout << "Could not open gain map " << gainfname << endl;
+  // }
+
+#ifdef CORR
+ cout << "Applying common mode  " << ncol_cm << endl;
+  cm=new moench03CommonMode(ncol_cm);
+  
+
+  cout << "Applying ghost corrections " << xt_ghost << endl;
+  gs=new moench03GhostSummation(decoder, xt_ghost);
+#endif
+
+  singlePhotonDetector *filter=new singlePhotonDetector(decoder,csize, nsigma, 1, cm, nped, 200, -1, -1, gainmap, gs);
+
+  if (gainfname) {
+
+    if (filter->readGainMap(gainfname))
+      cout << "using gain map " << gainfname << endl;
+    else
+      cout << "Could not open gain map " << gainfname << endl;
+  } else
+    thr=0.15*thr;
+  filter->newDataSet();
+  int dsize=decoder->getDataSize();
+
+
+  char data[dsize];
+
+
+
+
+
+  //#ifndef ANALOG
+  if (thr>0) {
+    cout << "threshold is " << thr << endl;
+    //#ifndef ANALOG
+    filter->setThreshold(thr);
+    //#endif
+      cf=0;
+
+  } else
+    cf=1;
+  //#endif
+
+
+  filter->setROI(xmin,xmax,ymin,ymax);
+  std::time(&end_time);
+  cout << std::ctime(&end_time) <<   endl;
+
+  char* buff;
+
+  // multiThreadedAnalogDetector *mt=new multiThreadedAnalogDetector(filter,nthreads,fifosize);
+ multiThreadedCountingDetector *mt=new multiThreadedCountingDetector(filter,nthreads,fifosize);
+#ifndef ANALOG
+    mt->setDetectorMode(ePhotonCounting);
+    cout << "Counting!" << endl;
+  if (thr>0) {
+    cf=0;
+  } 
+#endif
+//{
+#ifdef ANALOG
+      mt->setDetectorMode(eAnalog);
+      cout << "Analog!" << endl;
+      cf=0;
+      //thr1=thr;
+#endif  
+      //  }
+  
+  mt->StartThreads();
+  mt->popFree(buff);
+
+
+  //  cout << "mt " << endl;
+
+  int ifr=0;
+
+  double ped[nx*ny], *ped1;
+
+
+  if (pedfile) {
+
+    cout << "PEDESTAL " << endl;
+    sprintf(imgfname,"%s/pedestals.tiff",outdir);
+
+    if (string(pedfile).find(".tif")==std::string::npos){
+      sprintf(fname,"%s.raw",pedfile);
+      cout << fname << endl ;
+      std::time(&end_time);
+      cout << "aaa" << std::ctime(&end_time) <<   endl;
+      
+      
+      mt->setFrameMode(ePedestal);
+      // sprintf(fn,fformat,irun);
+      filebin.open((const char *)(fname), ios::in | ios::binary);
+      //      //open file
+      if (filebin.is_open()){
+	ff=-1;
+	while (decoder->readNextFrame(filebin, ff, np,buff)) {
+	  if (np==40) {
+	    mt->pushData(buff);
+	    mt->nextThread();
+	    mt->popFree(buff);
+	    ifr++;
+	    if (ifr%100==0) 
+	      cout << ifr << " " << ff << " " << np << endl;
+	  } else
+	    cout << ifr << " " << ff << " " << np << endl;
+	  ff=-1;
+	}
+	filebin.close();	 
+	while (mt->isBusy()) {;}
+	
+      } else 
+	cout << "Could not open pedestal file "<< fname << " for reading " << endl;
+    } else {
+      float *pp=ReadFromTiff(pedfile, nny, nnx);
+      if (pp && nnx==nx && nny==ny) {
+	for (int i=0; i<nx*ny; i++) {
+	  ped[i]=pp[i];
+	}
+	delete [] pp;
+	mt->setPedestal(ped);
+	// ped1=mt->getPedestal();
+	
+	// for (int i=0; i<nx*ny; i++) {
+
+	//   cout << ped[i]<<"/"<<ped1[i] << " " ;
+	// }
+	cout << "Pedestal set from tiff file " << pedfile << endl;
+      } else {
+	cout << "Could not open pedestal tiff file "<< pedfile << " for reading " << endl;
+      }
+    }
+    mt->writePedestal(imgfname);
+    std::time(&end_time);
+    cout << std::ctime(&end_time) <<   endl;
+  }
+  
+
+
+  ifr=0;
+  int ifile=0;
+  
+  mt->setFrameMode(eFrame);
+
+  for (int irun=runmin; irun<=runmax; irun++) {
+    cout << "DATA " ;
+    // sprintf(fn,fformat,irun);
+    sprintf(ffname,"%s/%s.raw",indir,fformat);
+    sprintf(fname,ffname,irun);
+    sprintf(ffname,"%s/%s.tiff",outdir,fformat);
+    sprintf(imgfname,ffname,irun);
+    sprintf(ffname,"%s/%s.clust",outdir,fformat);
+    sprintf(cfname,ffname,irun);
+    cout << fname << " " ;
+    cout << imgfname << endl;
+    std::time(&end_time);
+    cout << std::ctime(&end_time) <<    endl;
+    //  cout <<  fname << " " << outfname << " " << imgfname <<  endl;
+    filebin.open((const char *)(fname), ios::in | ios::binary);
+    //      //open file
+    ifile=0;
+    if (filebin.is_open()){
+      if (thr<=0 && cf!=0) { //cluster finder
+	  if (of==NULL) {
+	    of=fopen(cfname,"w");
+	    if (of) {
+	    mt->setFilePointer(of);
+	    cout << "file pointer set " << endl;
+	    } else {
+	      cout << "Could not open "<< cfname << " for writing " << endl;
+	      mt->setFilePointer(NULL);
+	      return 1;
+	    }
+	  }
+      }
+      //     //while read frame 
+      ff=-1;
+      ifr=0;
+      while (decoder->readNextFrame(filebin, ff, np,buff)) {
+	if (np==40) {
+	//	cout << "*"<<ifr++<<"*"<<ff<< endl;
+	//	cout << ff << " " << np << endl;
+  	//         //push
+	  mt->pushData(buff);
+	  // 	//         //pop
+	  mt->nextThread();
+	  // // 		//	cout << " " << (void*)buff;
+	  mt->popFree(buff);
+
+
+
+
+	  ifr++;
+	  if (ifr%100==0) cout << ifr << " " << ff << endl;
+	  if (nframes>0) {
+	    if (ifr%nframes==0) {
+	      //The name has an additional "_fXXXXX" at the end, where "XXXXX" is the initial frame number of the image (0,1000,2000...)
+	      
+	      sprintf(ffname,"%s/%s_f%05d.tiff",outdir,fformat,ifile);
+	      sprintf(imgfname,ffname,irun);
+	      //cout << "Writing tiff to " << imgfname << " " << thr1 << endl;
+	      mt->writeImage(imgfname, thr1);
+	      mt->clearImage();
+	      ifile++;
+	    }
+	  }
+	} else
+	  cout << ifr << " " << ff << " " << np << endl;
+	ff=-1;
+      }
+      cout << "--" << endl;
+      filebin.close();	 
+      //      //close file 
+      //     //join threads
+      while (mt->isBusy()) {;}
+      if (nframes>=0) {
+	if (nframes>0) {
+	    sprintf(ffname,"%s/%s_f%05d.tiff",outdir,fformat,ifile);
+	    sprintf(imgfname,ffname,irun);
+	} else {
+	  sprintf(ffname,"%s/%s.tiff",outdir,fformat);
+	  sprintf(imgfname,ffname,irun);
+	}
+	cout << "Writing tiff to " << imgfname << " " << thr1 <<endl;
+	mt->writeImage(imgfname, thr1);
+	mt->clearImage();
+	if (of) {
+	  fclose(of);
+	  of=NULL;
+	  mt->setFilePointer(NULL);
+	}
+      }	
+      std::time(&end_time);
+      cout << std::ctime(&end_time) <<   endl;
+    } else 
+      cout << "Could not open "<< fname << " for reading " << endl;
+  }
+  if (nframes<0){
+    sprintf(ffname,"%s/%s.tiff",outdir,fformat);
+    strcpy(imgfname,ffname);
+    cout << "Writing tiff to " << imgfname << " " << thr1 <<endl;
+    mt->writeImage(imgfname, thr1);
+  }
+    
+
+
+  return 0;
+}
+

slsDetectorCalibration/moenchExecutables/CMakeLists.txt

@@ -3,11 +3,11 @@ find_package(TIFF REQUIRED)
 
 #Moench ZMQ
 add_executable(moenchZmqProcess moenchZmqProcess.cpp ../tiffIO.cpp)
-add_compile_definitions(moenchZmqProcess NEWZMQ INTERP)
+target_compile_definitions(moenchZmqProcess PRIVATE NEWZMQ INTERP)
 
 #Moench04 ZMQ
 add_executable(moench04ZmqProcess moenchZmqProcess.cpp ../tiffIO.cpp)
-add_compile_definitions(moenchZmqProcess NEWZMQ INTERP MOENCH04)
+target_compile_definitions(moench04ZmqProcess PRIVATE NEWZMQ INTERP MOENCH04)
 
 
 #Both executables should have the same includes and output dirs

slsDetectorCalibration/moenchExecutables/moenchZmqProcess.cpp

@@ -1,6 +1,8 @@
 //#define WRITE_QUAD
 #define DEVELOPER
 #undef CORR
+#undef MOENCH04
+
 #define C_GHOST 0.0004
 
 #define CM_ROWS 20
@@ -28,6 +30,7 @@
 #include <fstream> 
 #include "tiffIO.h"
 
+
 #include <rapidjson/document.h> //json header in zmq stream
 
 #include<iostream>
@@ -135,9 +138,11 @@ int main(int argc, char *argv[]) {
 
   //slsDetectorData *det=new moench03T1ZmqDataNew(); 
 #ifndef MOENCH04
+  cout << "This is a Moench03" << endl;
   moench03T1ZmqDataNew *det=new moench03T1ZmqDataNew(); 
 #endif
 #ifdef MOENCH04
+  cout << "This is a Moench04" << endl;
   moench04CtbZmq10GbData *det=new moench04CtbZmq10GbData(); 
 #endif
   cout << endl << " det" <<endl;
@@ -614,61 +619,8 @@ int main(int argc, char *argv[]) {
 	      version=zHeader.version;//doc["version"].GetUint();
 	      /*document["bitmode"].GetUint(); zHeader.dynamicRange
 
-document["fileIndex"].GetUint64(); zHeader.fileIndex
 
-document["detshape"][0].GetUint();
-zHeader.ndetx
-
-document["detshape"][1].GetUint();
-zHeader.ndety
-
-document["shape"][0].GetUint();
-zHeader.npixelsx
-
-document["shape"][1].GetUint();
-zHeader.npixelsy
-
-document["size"].GetUint(); zHeader.imageSize
-
-document["acqIndex"].GetUint64(); zHeader.acqIndex
-
-document["frameIndex"].GetUint64(); zHeader.frameIndex
-
-document["fname"].GetString(); zHeader.fname
-
-document["frameNumber"].GetUint64(); zHeader.frameNumber
-
-document["expLength"].GetUint(); zHeader.expLength
-
-document["packetNumber"].GetUint(); zHeader.packetNumber
-
-document["bunchId"].GetUint64(); zHeader.bunchId
-
-document["timestamp"].GetUint64(); zHeader.timestamp
-
-document["modId"].GetUint(); zHeader.modId
-
-document["row"].GetUint(); zHeader.row
-
-document["column"].GetUint(); zHeader.column
-
-document["reserved"].GetUint(); zHeader.reserved
-
-document["debug"].GetUint(); zHeader.debug
-
-document["roundRNumber"].GetUint(); zHeader.roundRNumber
-
-document["detType"].GetUint(); zHeader.detType
-
-document["version"].GetUint(); zHeader.version
-
-document["flippedDataX"].GetUint(); zHeader.flippedDataX
-
-document["quad"].GetUint(); zHeader.quad
-
-document["completeImage"].GetUint(); zHeader.completeImage
 	      */
-	      //dataSize=size;
 
 	      //strcpy(fname,filename.c_str());
 	      fname=filename;

slsDetectorCalibration/multiThreadedAnalogDetector.h

@@ -17,7 +17,7 @@
 #include <pthread.h>
 
 #include "analogDetector.h"
-#include "sls/CircularFifo.h"
+#include "circularFifo.h"
 #include <unistd.h>
 #include <stdio.h>
 #include <stdlib.h>
@@ -123,7 +123,7 @@ public:
      return fifoFree->pop(ptr);
    }
 
-   virtual int isBusy() {if (fifoData->isEmpty() && busy==0) return 0; else return 1;}
+   virtual int isBusy() {if (fifoData->isEmpty() && busy==0) return 0; return 1;}
    
    //protected:
    /** Implement this method in your subclass with the code you want your thread to run. */

slsDetectorSoftware/include/sls/Detector.h

@@ -498,7 +498,7 @@ class Detector {
 
     /** Non blocking: Abort detector acquisition. Status changes to IDLE or
      * STOPPED. Goes to stop server. */
-    void stopDetector();
+    void stopDetector(Positions pos = {});
 
     /** IDLE, ERROR, WAITING, RUN_FINISHED, TRANSMITTING, RUNNING, STOPPED \n
      * Goes to stop server */
@@ -1737,4 +1737,4 @@ class Detector {
     void updateRxRateCorrections();
 };
 
-} // namespace sls
+} // namespace sls

slsDetectorSoftware/src/CmdProxy.h

@@ -1444,7 +1444,7 @@ class CmdProxy {
         "\n\t[Mythen3] Starts detector readout. Status changes to TRANSMITTING "
         "and automatically returns to idle at the end of readout.");
 
-    EXECUTE_SET_COMMAND_NOID(stop, stopDetector,
+    EXECUTE_SET_COMMAND(stop, stopDetector,
                              "\n\tAbort detector acquisition. Status changes "
                              "to IDLE or STOPPED. Goes to stop server.");
 

slsDetectorSoftware/src/Detector.cpp

@@ -697,7 +697,7 @@ void Detector::startDetectorReadout() {
     pimpl->Parallel(&Module::startReadout, {});
 }
 
-void Detector::stopDetector() { pimpl->Parallel(&Module::stopAcquisition, {}); }
+void Detector::stopDetector(Positions pos) { pimpl->Parallel(&Module::stopAcquisition, pos); }
 
 Result<defs::runStatus> Detector::getDetectorStatus(Positions pos) const {
     return pimpl->Parallel(&Module::getRunStatus, pos);
@@ -2104,4 +2104,4 @@ std::vector<int> Detector::getPortNumbers(int start_port) {
     return res;
 }
 
-} // namespace sls
+} // namespace sls

slsDetectorSoftware/src/Module.cpp

@@ -412,7 +412,7 @@ void Module::loadSettingsFile(const std::string &fname) {
         if (shm()->myDetectorType == MYTHEN3) {
             serialNumberWidth = 4;
         }
-        if (fname.find(".sn") == std::string::npos) {
+        if ((fname.find(".sn") == std::string::npos) && (fname.find(".trim") == std::string::npos)) {
             ostfn << ".sn" << std::setfill('0') << std::setw(serialNumberWidth)
                   << std::dec << getSerialNumber();
         }