From 5dbfbdb82b69c63838549f1de478c653957baf8b Mon Sep 17 00:00:00 2001
From: Anna Bergamaschi <anna.bergamaschi@psi.ch>
Date: Tue, 30 May 2017 12:21:48 +0200
Subject: [PATCH] improved readNextFrame for moench 10G streaming; implemented
 many photons finder

---
 slsDetectorCalibration/MovingStat.h           |   7 +-
 slsDetectorCalibration/Mythen3_01_jctbData.h  |  56 ++-
 slsDetectorCalibration/etaVEL/imageMacro.C    |  21 +-
 slsDetectorCalibration/etaVEL/imageMacroG2D.C |  48 ++-
 .../etaVEL/imageMacroSample.C                 |   2 +-
 slsDetectorCalibration/moench03Ctb10GbData.h  | 143 +++++---
 slsDetectorCalibration/moench03TCtb10GbData.h | 320 ++++++++++--------
 slsDetectorCalibration/pedestalSubtraction.h  |   2 +-
 slsDetectorCalibration/singlePhotonDetector.h | 160 +++++++--
 9 files changed, 524 insertions(+), 235 deletions(-)

diff --git a/slsDetectorCalibration/MovingStat.h b/slsDetectorCalibration/MovingStat.h
index 968a3eab1..dd8ad7e8a 100755
--- a/slsDetectorCalibration/MovingStat.h
+++ b/slsDetectorCalibration/MovingStat.h
@@ -29,11 +29,14 @@ class MovingStat
 	/**
 	   clears the moving average number of samples parameter, mean and standard deviation
 	*/
-        void Set(double val)
+        void Set(double val, double rms=0)
         {
 	  m_n = n; 
 	  m_newM=val*n; 
-	  m_newM2=val*val*n;
+	  if (rms<=0)
+	    m_newM2=val*val*n;
+	  else
+	    m_newM2=(n*rms*rms+m_newM*m_newM/n);
         }
 	
 	
diff --git a/slsDetectorCalibration/Mythen3_01_jctbData.h b/slsDetectorCalibration/Mythen3_01_jctbData.h
index 1fa016e2f..101074dea 100644
--- a/slsDetectorCalibration/Mythen3_01_jctbData.h
+++ b/slsDetectorCalibration/Mythen3_01_jctbData.h
@@ -44,8 +44,9 @@ class mythen3_01_jctbData : public slsDetectorData<short unsigned int> {
  
 
  static short unsigned int* mythen03_frame(char *ptr, int dr=24,  int nch=64*3, int off=5) {
+   // off=0;
     int iarg;
-    int64_t word;
+    int64_t word, *wp;
     short unsigned int* val=new short unsigned int[nch];
     int bit[64];
     int nb=2;
@@ -53,30 +54,53 @@ class mythen3_01_jctbData : public slsDetectorData<short unsigned int> {
     int idr=0;
     int ib=0;
     int iw=0;
+    int ii=0;
     bit[0]=19;
     bit[1]=8;
     idr=0;
     for (ib=0; ib<nch; ib++) {
       val[ib]=0;
     }
+    wp=(int64_t*)ptr;
+    
     for (iw=0; iw<nch/nb; iw) {
-      word=*((int64_t*)(ptr+8));
-	if (ioff<off) ioff++;
-	else {
-	  if (idr<16) {
-	    for (ib=0; ib<nb; ib++) {
-	      if (word&(1<<bit[ib]))    val[iw+nch*ib/nb]|=(1<<idr);
-	    }//end for()
-	  }
-	  idr++;
-	  if (idr==dr) {
-	    idr=0;
-	    iw++;
-	  }//end if()
-	}//end else()
+      word=*wp;;
+      if (ioff<off) {
+	ioff++;
+	 cout <<"*";
+      }  else {
+	
+	if (idr<16) {
+	  for (ib=0; ib<nb; ib++) {
+	    if (word&(1<<bit[ib]))  {
+	      cout << "+" ;
+	      val[iw+nch*(ib/nb)]|=(1<<idr);
+	    }  else { 
+	      	cout << "-" ; 
+	       } 
+	  }//end for()
+	}
+	
+	idr++;
+	
+	
+	if (idr==dr) {
+	  idr=0;
+	  //  cout << dec << " " << iw << " " << val[iw] << " " << val[iw+nch/2] << endl;
+	  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;
+  
     return val;
-  }
+ }
 
  virtual int setFrameNumber(int f=0) {if (f>=0) frameNumber=f; return frameNumber; };
  virtual int setDynamicRange(int d=-1) {if (d>0 && d<=24) dynamicRange=d; return dynamicRange;};
diff --git a/slsDetectorCalibration/etaVEL/imageMacro.C b/slsDetectorCalibration/etaVEL/imageMacro.C
index b037604af..16dc9117d 100644
--- a/slsDetectorCalibration/etaVEL/imageMacro.C
+++ b/slsDetectorCalibration/etaVEL/imageMacro.C
@@ -1,11 +1,12 @@
-void imageMacro(char *name){
+TH2D *imageMacro(char *name) {
+
   //TH2D *makeNorm(){
 
   //TFile ff("/local_zfs_raid/tomcat_20160528/trees/img_blank_eta_gmap.root");
   //TH2D *hff=(TH2D*)ff.Get("imgHR");
-    TFile ff("/local_zfs_raid/tomcat_20160528/trees/img_blank_eta_nb25.root");
+    TFile *ff=new TFile("/mnt/moench_data/tomcat_20160528_img/img_blank_eta_nb25.root");
   //  TFile ff("/local_zfs_raid/tomcat_20160528/trees/img_blank_eta_gcorr_nb25.root");
-  TH2D *hff=(TH2D*)ff.Get("blankHR");
+  TH2D *hff=(TH2D*)ff->Get("blankHR");
   hff->SetName("imgBlank");
   TH2D *hpixel=new TH2D("hpixel","hpixel",25,0,25,25,0,25);
   for (int ibx=10*25; ibx<hff->GetNbinsX()-10*25; ibx++) {
@@ -45,9 +46,9 @@ void imageMacro(char *name){
   }    else {
     
   
-  sprintf(nn,"/local_zfs_raid/tomcat_20160528/trees/img_%s_eta_nb25.root", name);
+  sprintf(nn,"/mnt/moench_data/tomcat_20160528_img/img_%s_eta_nb25.root", name);
   // if (strcmp(name,"blank"))
-    TFile fg(nn);
+    TFile *fg=new TFile(nn);
     // else
     // TFile fg=gDirectory;
 
@@ -57,13 +58,21 @@ void imageMacro(char *name){
   //  TFile fg("/local_zfs_raid/tomcat_20160528/trees/img_sample_eta_nb25.root");
   // TFile fg("/local_zfs_raid/tomcat_20160528/trees/img_grating_1d_eta_gcorr_nb25.root");
    sprintf(nn,"%sHR",name);
-   TH2D *hg=(TH2D*)fg.Get(nn);
+   TH2D *hg=(TH2D*)fg->Get(nn);
    // hg->SetName("imgGrating");
 
   //hg->Divide(hff);
   }
    if (hpix)
      hg->Divide(hpix);
+   new TCanvas();
+   hg->Draw("colz");
+
+   return hg;
+}
+
+
+void imageMacro(TH2D *hg){
 
   Double_t imageData[200*400*25*25];
   const int nsigma=5.;
diff --git a/slsDetectorCalibration/etaVEL/imageMacroG2D.C b/slsDetectorCalibration/etaVEL/imageMacroG2D.C
index 8f2a8a598..6caa8eb31 100644
--- a/slsDetectorCalibration/etaVEL/imageMacroG2D.C
+++ b/slsDetectorCalibration/etaVEL/imageMacroG2D.C
@@ -1,7 +1,31 @@
-{
+{ 
+  TColor::InitializeColors();
+  const Int_t NRGBs = 5;
+  const Int_t NCont = 255;//90;
+
+  Double_t stops[NRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 };
+  Double_t red[NRGBs]   = { 0.00, 0.00, 0.87, 1.00, 0.51 };
+  Double_t green[NRGBs] = { 0.00, 0.81, 1.00, 0.20, 0.00 };
+  Double_t blue[NRGBs]  = { 0.51, 1.00, 0.12, 0.00, 0.00 };
+  Double_t gray[NRGBs]  = { 1., 0.34, 0.61, 0.84, 1.00};
+  Double_t zero[NRGBs]  = { 0., 0.0,0.0, 0.0, 0.00};
+  //TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
+   TColor::CreateGradientColorTable(NRGBs, stops, gray, gray, gray, NCont);
+ 
+  gStyle->SetNumberContours(NCont);
+
+  gStyle->SetPadTopMargin(0);
+  gStyle->SetPadRightMargin(0);
+  gStyle->SetPadBottomMargin(0);
+  gStyle->SetPadLeftMargin(0);
+
+
+
+  gROOT->ForceStyle();
+
   // TFile ff("/local_zfs_raid/tomcat_20160528/trees/img_blank_eta.root");
   // TH2D *hff=(TH2D*)ff.Get("imgHR");
-  TFile ff("/local_zfs_raid/tomcat_20160528/trees/img_blank_eta_gmap_v2.root");
+  TFile ff("/mnt/moench_data/tomcat_20160528_img/img_blank_eta_nb25.root");
   TH2D *hff=(TH2D*)ff.Get("blankHR");
   hff->SetName("imgBlank");
   TH2D *hpixel=new TH2D("hpixel","hpixel",25,0,25,25,0,25);
@@ -32,7 +56,7 @@
 
   // TFile fg("/local_zfs_raid/tomcat_20160528/trees/img_grating_2d_eta.root");
   // TH2D *hg=(TH2D*)fg.Get("imgHR");
-  TFile fg("/local_zfs_raid/tomcat_20160528/trees/img_grating_2d_eta_gmap_v2.root");
+  TFile fg("/mnt/moench_data/tomcat_20160528_img/img_grating_2d_eta_nb25.root");
   TH2D *hg=(TH2D*)fg.Get("grating_2dHR");
   hg->SetName("imgGrating");
 
@@ -50,9 +74,21 @@
  
    hg->Scale(1/ng);
 
-   new TCanvas();
-   hg->SetMaximum(1.);
-   hg->Draw("colz");
+   new TCanvas("c1","c1",800,800);
+   hg->SetMaximum(2.);
+   Double_t xmin=hg->GetXaxis()->GetXmin();
+   Double_t xmax=hg->GetXaxis()->GetXmax();
+   Double_t ymin=hg->GetYaxis()->GetXmin();
+   Double_t ymax=hg->GetYaxis()->GetXmax();
+   hg->SetTitle(";mm;mm;Normalized intensity (a.u.)");
+   hg->GetXaxis()->Set(hg->GetXaxis()->GetNbins(),xmin*0.025-180.*0.025+0.1,xmax*0.025-180.*0.025+0.1);
+   hg->GetYaxis()->Set(hg->GetYaxis()->GetNbins(),ymin*0.025-180.*0.025-1.5,ymax*0.025-180*0.025-1.5);
+   
+   hg->GetXaxis()->SetRangeUser(0,0.5);
+   hg->GetYaxis()->SetRangeUser(0,0.5);
+   
+   hg->Draw("col");
+   hg->SetStats(kFALSE);
    // new TCanvas();
    // hsg->SetMaximum(1.);
    // hsg->Draw("colz");
diff --git a/slsDetectorCalibration/etaVEL/imageMacroSample.C b/slsDetectorCalibration/etaVEL/imageMacroSample.C
index 6eb9627f3..5fced606c 100644
--- a/slsDetectorCalibration/etaVEL/imageMacroSample.C
+++ b/slsDetectorCalibration/etaVEL/imageMacroSample.C
@@ -35,7 +35,7 @@
 
   // TFile fg("/local_zfs_raid/tomcat_20160528/trees/img_grating_2d_eta.root");
   // TH2D *hg=(TH2D*)fg.Get("imgHR");
-  TFile fg("/local_zfs_raid/tomcat_20160528/trees/img_sample_eta_gmap_v2.root");
+  TFile fg("/mnt/moench_data/tomcat_20160528_img/img_sample_eta_gmap_v2.root");
   TH2D *hg=(TH2D*)fg.Get("sampleHR");
   hg->SetName("imgSample");
 
diff --git a/slsDetectorCalibration/moench03Ctb10GbData.h b/slsDetectorCalibration/moench03Ctb10GbData.h
index 52b7575e4..fb0b90bf7 100644
--- a/slsDetectorCalibration/moench03Ctb10GbData.h
+++ b/slsDetectorCalibration/moench03Ctb10GbData.h
@@ -6,12 +6,13 @@
 
 class moench03Ctb10GbData : public slsReceiverData<uint16_t> {
  
- private:
+ protected:
   
   int iframe;
   int nadc;
   int sc_width;
   int sc_height;
+
  public:
 
 
@@ -165,58 +166,116 @@ class moench03Ctb10GbData : public slsReceiverData<uint16_t> {
 /*     }; */
 
 
-  	virtual char *readNextFrame(ifstream &filebin, int& fnum) {
-	  char *data=new char[packetSize*nPackets];
-	  char *retval=0;
-	  int np=0, nd;
-	  fnum = -1;
-	  int  pn;
-	  char aa[8224];
-	  char *packet=(char *)aa;
-
-	  if (filebin.is_open()) {
+  virtual char *readNextFrame(ifstream &filebin, int& fnum, char *data=NULL) {
+    int dd=0;
+    if (data==NULL) {
+      data=new char[packetSize*nPackets];
+      dd=1;
+    }
+      char *retval=0;
+      int np=0, nd;
+      fnum = -1;
+      int  pn;
+      char aa[8224];
+      char *packet=(char *)aa;
+      int place;
+      
+      if (filebin.is_open()) {
 
 
 
-
-	    while(filebin.read((char*)packet, 8208) && np<nPackets){
-	      pn=getPacketNumber(packet);
-	      
-	      if (pn==1 && fnum<0)
-		fnum= getFrameNumber(packet);
-	      
-	      // cout << "fn: " << fnum << "\t pn: " << pn << endl;
-	      if (fnum>=0) {
-		if (getFrameNumber(packet) !=fnum) { 
-		  
-		  if (np==0){
-		    delete [] data;
-		    return NULL;
-		  } else
-		    return data;
-		}
-		
-		memcpy(data+(pn-1)*packetSize, packet, packetSize);
-		np++;
-		
+        place=filebin.tellg();
+        while(filebin.read((char*)packet, 8208) && np<nPackets){
+          pn=getPacketNumber(packet);
+          if (pn==1 && fnum<0)
+	    fnum= getFrameNumber(packet);
+          
+	  //	  cout <<getFrameNumber(packet)<< " fn: " << fnum << "\t pn: " << pn << " np " << np << endl;
+          if (fnum>=0) {
+	    if (getFrameNumber(packet) !=fnum) {
+	      cout << "-"<<endl;
+	      filebin.seekg(place);
+	      if (np==0){
+		if (dd)
+		  delete [] data;
+		return NULL;
+	      } else {
+		return data;
 	      }
 	    }
+        
+	    memcpy(data+(pn-1)*packetSize, packet, packetSize);
+	    np++;
+	    if (np==nPackets)
+	      return data;
+        
+          }
+          place=filebin.tellg();
+        }
+        
+      }
+
+	      
+      if (np==0){
+	cout << "-"<<endl;
+        filebin.seekg(place);
+	if (dd)
+	  delete [] data;
+        return NULL;
+      }
+
+      //filebin.seekg(place);
+      return data;
+    };
+
+/*   	virtual char *readNextFrame(ifstream &filebin, int& fnum) { */
+/* 	  char *data=new char[packetSize*nPackets]; */
+/* 	  char *retval=0; */
+/* 	  int np=0, nd; */
+/* 	  fnum = -1; */
+/* 	  int  pn; */
+/* 	  char aa[8224]; */
+/* 	  char *packet=(char *)aa; */
+
+/* 	  if (filebin.is_open()) { */
+
+
+
+
+/* 	    while(filebin.read((char*)packet, 8208) && np<nPackets){ */
+/* 	      pn=getPacketNumber(packet); */
+	      
+/* 	      if (pn==1 && fnum<0) */
+/* 		fnum= getFrameNumber(packet); */
+	      
+/* 	      // cout << "fn: " << fnum << "\t pn: " << pn << endl; */
+/* 	      if (fnum>=0) { */
+/* 		if (getFrameNumber(packet) !=fnum) {  */
+		  
+/* 		  if (np==0){ */
+/* 		    delete [] data; */
+/* 		    return NULL; */
+/* 		  } else */
+/* 		    return data; */
+/* 		} */
+		
+/* 		memcpy(data+(pn-1)*packetSize, packet, packetSize); */
+/* 		np++; */
+		
+/* 	      } */
+/* 	    } */
 	    
-	  }
+/* 	  } */
 
-	  if (np==0){
-	    delete [] data;
-	    return NULL;
-	  }
+/* 	  if (np==0){ */
+/* 	    delete [] data; */
+/* 	    return NULL; */
+/* 	  } */
 
-	};
+/* 	}; */
 
 int getPacketNumber(int x, int y) {return dataMap[y][x]/8208;};
 
-  	virtual char *readNextFrame(ifstream &filebin) {
-	  int fnum;
-	  return readNextFrame(filebin, fnum);
-	};
 
 };
 
diff --git a/slsDetectorCalibration/moench03TCtb10GbData.h b/slsDetectorCalibration/moench03TCtb10GbData.h
index 2eaeb0019..989012991 100644
--- a/slsDetectorCalibration/moench03TCtb10GbData.h
+++ b/slsDetectorCalibration/moench03TCtb10GbData.h
@@ -1,17 +1,12 @@
 #ifndef MOENCH03TCTB10GBDATA_H
 #define  MOENCH03TCTB10GBDATA_H
-#include "slsReceiverData.h"
+#include "moench03Ctb10GbData.h"
 
 
 
-class moench03TCtb10GbData : public slsReceiverData<uint16_t> {
+class moench03TCtb10GbData : public moench03Ctb10GbData { //slsReceiverData<uint16_t> {
  
- private:
-  
-  int iframe;
-  int nadc;
-  int sc_width;
-  int sc_height;
+
  public:
 
 
@@ -26,8 +21,9 @@ class moench03TCtb10GbData : public slsReceiverData<uint16_t> {
   
 
   //  moench03TCtb10GbData(int ns=5000): slsDetectorData<uint16_t>(400, 400, 8208*40, NULL, NULL) , nadc(32), sc_width(25), sc_height(200) {
- moench03TCtb10GbData(int ns=5000): slsReceiverData<uint16_t>(400, 400, 40, 8208), nadc(32), sc_width(25), sc_height(200) {
-
+ moench03TCtb10GbData(int ns=5000): moench03Ctb10GbData(ns) {//slsReceiverData<uint16_t>(400, 400, 40, 8208), nadc(32), sc_width(25), sc_height(200) {
+    // cout <<"constructor"<< endl;
+      // nadc=32;
     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,\
@@ -47,11 +43,13 @@ class moench03TCtb10GbData : public slsReceiverData<uint16_t> {
 
 	for (iadc=0; iadc<nadc; iadc++) {
 	  i=128*ip+is;
-	  adc4=(int)iadc/4;
+	  //cout << i << endl;
+	   adc4=(int)iadc/4;
 	  if (i<sc_width*sc_height) {
 	    //  for (int i=0; i<sc_width*sc_height; i++) {
 	    col=adc_nr[iadc]+(i%sc_width);
 	    if (adc4%2==0) {
+	      // if (iadc<(nadc/2)) {
 	      row=199-i/sc_width; 	            
 	    } else {
 	      row=200+i/sc_width;
@@ -63,181 +61,221 @@ class moench03TCtb10GbData : public slsReceiverData<uint16_t> {
 	}
       }
     }
-
+    //  cout <<"map"<< endl;
     int ipacket;
     int ibyte;
     int ii=0;
-    
     for (int ipacket=0; ipacket<npackets; ipacket++) {
       for (int ibyte=0;  ibyte< 8208/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;
-	  iadc=ii%nadc;
-	  adc4 = (int)iadc/4;
+    	i=ipacket*8208/2+ibyte;
+    	//	cout << i << " ";
+    	if (ibyte<8) {
+    	//header!
+    	  xmap[i]=-1;
+    	  ymap[i]=-1;
+    	} else {
+    	  // ii=ibyte+128*32*ipacket;
+    	  // cout << nadc << endl;
+    	  isample=ii/nadc;
+    	  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;
-
-	  }
+    	  ix=isample%sc_width;
+    	  iy=isample/sc_width;
+    	  //  cout << isample << " " << iadc << " " << ix << " " << iy ;//<< endl;
+    	  if (adc4%2==0) {//iadc<(nadc/2)) {
+    	  //if (iadc<(nadc/2)) {
+    	    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++;
+    	}
       }
     }
     
 
     
+    /* int ipacket; */
+    /* int ibyte; */
+    /* int ii=0; */
+    
+    /* for (int ipacket=0; ipacket<npackets; ipacket++) { */
+    /*   for (int ibyte=0;  ibyte< 8208/2; ibyte++) { */
+    /* 	i=ipacket*8208/2+ibyte; */
+    /* 	cout << i << endl; */
+    /* 	if (ibyte<8) { */
+    /* 	//header! */
+    /* 	  xmap[i]=-1; */
+    /* 	  ymap[i]=-1; */
+    /* 	} else { */
+    /* 	  ii=ibyte+128*32*ipacket; */
+    /* 	  isample=ii/nadc; */
+    /* 	  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++; */
+    /* 	} */
+    /*   } */
+    /* } */
+    
+
+    
+    //  cout <<"done"<< endl;
     
     iframe=0;
     //  cout << "data struct created" << endl;
   };
     
 
-     /**
+/*      /\** */
 
-     Returns the frame number for the given dataset. Purely virtual func.
-     \param buff pointer to the dataset
-     \returns frame number
-
-  */
-
-
-  int getFrameNumber(char *buff){return *((int*)buff)&0xffffffff;};   
-
-  /**
-
-     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 ((*(((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  */
+/*      Returns the frame number for the given dataset. Purely virtual func. */
+/*      \param buff pointer to the dataset */
+/*      \returns frame number */
 
 /*   *\/ */
-/*     virtual  char *findNextFrame(char *data, int &ndata, int dsize){ndata=dsize; setDataSize(dsize);  return data;}; */
 
 
-/*    /\** */
+/*   int getFrameNumber(char *buff){return *((int*)buff)&0xffffffff;};    */
+
+/*   /\** */
+
+/*      Returns the packet number for the given dataset. purely virtual func */
+/*      \param buff pointer to the dataset */
+/*      \returns packet number number */
+
 
-/*      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; */
-/*     }; */
+/*   int getPacketNumber(char *buff){return ((*(((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;}; *\/ */
 
 
-  virtual char *readNextFrame(ifstream &filebin, int& ff, int &np) {
-	  char *data=new char[packetSize*nPackets];
-	  char *retval=0;
-	  int nd;
-	  np=0;
-	  int  pn;
-	  char aa[8224]={0};
-	  char *packet=(char *)aa;
-	  int  fnum = -1;
+/* /\*    /\\** *\/ */
+
+/* /\*      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, int &np) { */
+/* 	  char *data=new char[packetSize*nPackets]; */
+/* 	  char *retval=0; */
+/* 	  int nd; */
+/* 	  np=0; */
+/* 	  int  pn; */
+/* 	  char aa[8224]={0}; */
+/* 	  char *packet=(char *)aa; */
+/* 	  int  fnum = -1; */
 	  
-	  if (ff>=0)
-	    fnum=ff;
+/* 	  if (ff>=0) */
+/* 	    fnum=ff; */
 
-	  if (filebin.is_open()) {
+/* 	  if (filebin.is_open()) { */
 
 
-	    cout << "+";
+/* 	    cout << "+"; */
 
-	    while(filebin.read((char*)packet, 8208)){
+/* 	    while(filebin.read((char*)packet, 8208)){ */
 	      
-	      pn=getPacketNumber(packet);
+/* 	      pn=getPacketNumber(packet); */
 	      
-	      if (fnum<0)
-		fnum= getFrameNumber(packet);
+/* 	      if (fnum<0) */
+/* 		fnum= getFrameNumber(packet); */
 	      
-	      if (fnum>=0) {
-		if (getFrameNumber(packet) !=fnum) { 
+/* 	      if (fnum>=0) { */
+/* 		if (getFrameNumber(packet) !=fnum) {  */
 		  
-		  if (np==0){
-		    cout << "-";
-		    delete [] data;
-		    return NULL;
-		  } else
-		    filebin.seekg(-8208,ios_base::cur);
+/* 		  if (np==0){ */
+/* 		    cout << "-"; */
+/* 		    delete [] data; */
+/* 		    return NULL; */
+/* 		  } else */
+/* 		    filebin.seekg(-8208,ios_base::cur); */
 
-		    return data;
-		}
+/* 		    return data; */
+/* 		} */
 		
-		memcpy(data+(pn-1)*packetSize, packet, packetSize);
-		np++;
-		if (np==nPackets)
-		  break;
-		if (pn==nPackets)
-		  break;
-	      }
-	    }
+/* 		memcpy(data+(pn-1)*packetSize, packet, packetSize); */
+/* 		np++; */
+/* 		if (np==nPackets) */
+/* 		  break; */
+/* 		if (pn==nPackets) */
+/* 		  break; */
+/* 	      } */
+/* 	    } */
 	    
-	  }
+/* 	  } */
 	  
-	  if (np==0){
-	    cout << "?";
-	    delete [] data;
-	    return NULL;
-	  }// else if (np<nPackets)
-	    //  cout << "Frame " << fnum << " lost "  << nPackets-np << " packets " << endl;
+/* 	  if (np==0){ */
+/* 	    cout << "?"; */
+/* 	    delete [] data; */
+/* 	    return NULL; */
+/* 	  }// else if (np<nPackets) */
+/* 	    //  cout << "Frame " << fnum << " lost "  << nPackets-np << " packets " << endl; */
 
-	  return data;
+/* 	  return data; */
 	  
-	};
+/* 	}; */
 
-int getPacketNumber(int x, int y) {return dataMap[y][x]/8208;};
+/* int getPacketNumber(int x, int y) {return dataMap[y][x]/8208;}; */
 
-  	virtual char *readNextFrame(ifstream &filebin) {
-	  int fnum=-1, np;
-	  return readNextFrame(filebin, fnum, np);
-	};
+/*   	virtual char *readNextFrame(ifstream &filebin) { */
+/* 	  int fnum=-1, np; */
+/* 	  return readNextFrame(filebin, fnum, np); */
+/* 	}; */
 
 };
 
diff --git a/slsDetectorCalibration/pedestalSubtraction.h b/slsDetectorCalibration/pedestalSubtraction.h
index 6ec9343bd..f03f58fd5 100644
--- a/slsDetectorCalibration/pedestalSubtraction.h
+++ b/slsDetectorCalibration/pedestalSubtraction.h
@@ -40,7 +40,7 @@ class pedestalSubtraction  {
     virtual int SetNPedestals(int i=-1) {if (i>0) stat.SetN(i); return stat.GetN();};
     
     /** sets the moving average */
-    virtual void setPedestal(double val) {stat.Set(val);}
+    virtual void setPedestal(double val, double rms=0) {stat.Set(val, rms);}
 
   
   
diff --git a/slsDetectorCalibration/singlePhotonDetector.h b/slsDetectorCalibration/singlePhotonDetector.h
index 1e112d917..993462031 100644
--- a/slsDetectorCalibration/singlePhotonDetector.h
+++ b/slsDetectorCalibration/singlePhotonDetector.h
@@ -31,6 +31,8 @@ using namespace std;
     UNDEFINED_EVENT=-1
   };
 
+#ifndef DEF_QUAD
+#define DEF_QUAD
   enum quadrant {
     TOP_LEFT=0,
     TOP_RIGHT=1,
@@ -38,7 +40,7 @@ using namespace std;
     BOTTOM_RIGHT=3,
     UNDEFINED_QUADRANT=-1
   };
-
+#endif
 
 
 template <class dataType>
@@ -158,7 +160,7 @@ class singlePhotonDetector {
        \param iy pixel y coordinate
        \param val value to set
     */
-    virtual void setPedestal(int ix, int iy, double val){if (ix>=0 && ix<nx && iy>=0 && iy<ny) stat[iy][ix].setPedestal(val);};
+    virtual void setPedestal(int ix, int iy, double val, double rms=0){if (ix>=0 && ix<nx && iy>=0 && iy<ny) stat[iy][ix].setPedestal(val,rms);};
 
  
 
@@ -187,7 +189,122 @@ class singlePhotonDetector {
       }
       return clusterSize;
     }
-    
+
+
+
+    int *getNPhotons(char *data, double thr=-1) {
+
+      double val;
+      int *nph=new int[nx*ny];
+      double rest[ny][nx];
+      int cy=(clusterSizeY+1)/2;
+      int cs=(clusterSize+1)/2;
+
+      int ccs=clusterSize;
+      int ccy=clusterSizeY;
+
+
+      double tthr=thr;
+      int nn;
+      double max=0, tl=0, tr=0, bl=0,br=0, v;
+
+
+      if (thr>=0) {
+	cy=1;
+	cs=1;
+	ccs=1;
+	ccy=1;
+      }
+
+
+      for (int ix=0; ix<nx; ix++) {
+      for (int iy=0; iy<ny; iy++) {
+	
+	val=dataSign*(det->getValue(data, ix, iy)-getPedestal(ix,iy,0));
+	
+	if (thr<=0) tthr=nSigma*getPedestalRMS(ix,iy);
+	
+	nn=val/tthr;
+	
+	if (nn>0) {
+	  rest[iy][ix]=val-nn*tthr;
+	  nph[ix+nx*iy]=nn;
+	} else {
+	  rest[iy][ix]=val;
+	  nph[ix+nx*iy]=0;
+	}
+      }
+      
+      }
+      for (int ix=clusterSize/2; ix<clusterSize/2-1; ix++) {
+	for (int iy=clusterSizeY/2; iy<ny-clusterSizeY/2; iy++) {
+	  eventMask[iy][ix]=PEDESTAL;
+	
+	  if (thr<=0) tthr=nSigma*getPedestalRMS(ix,iy);
+	  
+	  max=0; 
+	  tl=0;
+	  tr=0; 
+	  bl=0;
+	  br=0;
+	  
+	  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) {
+		//cluster->set_data(rest[iy+ir][ix+ic], ic, ir);
+		v=rest[iy+ir][ix+ic];//cluster->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>tthr) {
+		    eventMask[iy][ix]=PHOTON;
+		  } 
+		}
+	      }
+	    }
+	  }
+	  
+
+	  //if (cluster->get_data(0,0)>=max) {
+	    if (rest[iy][ix]>=max) {
+	    
+	    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 (rest[iy][ix]>tthr || tot>sqrt(ccy*ccs)*tthr || quadTot>sqrt(cy*cs)*tthr) {
+	      nph[ix+nx*iy]++;
+	      rest[iy][ix]-=tthr;
+	    }
+	  }      
+	}
+	
+      }
+
+      
+      return nph;
+
+    }
 
 
 
@@ -207,6 +324,11 @@ class singlePhotonDetector {
       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;
+
+
+
       tot=0;
       quadTot=0;
       quad=UNDEFINED_QUADRANT;
@@ -250,31 +372,17 @@ class singlePhotonDetector {
 	      if (ir>=0 && ic>=0)
 		tr+=v;
 	      
-	      if (cluster->get_data(ic,ir)>max) {
+	      if (v>max) {
 		max=v;
 	      }
 	      if (ir==0 && ic==0) {
-		if (v>nSigma*cluster->rms) {
-		  eventMask[iy][ix]=PHOTON;
-		} else if (cluster->get_data(ic,ir)<-nSigma*cluster->rms)
+		if (v<-nSigma*cluster->rms)
 		  eventMask[iy][ix]=NEGATIVE_PEDESTAL;
 	      }
 	    }
 	  }
 	}
 	
-	if (eventMask[iy][ix]!=PHOTON && tot>sqrt(clusterSizeY*clusterSize)*nSigma*cluster->rms) {
-	  eventMask[iy][ix]=NEIGHBOUR;
-	} else if (eventMask[iy][ix]==PHOTON) {
-	  if (cluster->get_data(0,0)>=max) {
-	    eventMask[iy][ix]=PHOTON_MAX;
-	  }
-	} else if (eventMask[iy][ix]==PEDESTAL) {
-	  if (cm==0)
-	    addToPedestal(det->getValue(data, ix, iy),ix,iy);
-	}
-      
-	
 	if (bl>=br && bl>=tl && bl>=tr) {
 	  quad=BOTTOM_LEFT;
 	  quadTot=bl;
@@ -288,7 +396,19 @@ class singlePhotonDetector {
 	  quad=TOP_RIGHT;
 	  quadTot=tr;
 	} 
-
+	
+	if (max>nSigma*cluster->rms || tot>sqrt(clusterSizeY*clusterSize)*nSigma*cluster->rms || quadTot>cy*cs*nSigma*cluster->rms) {
+	  if (cluster->get_data(0,0)>=max) {
+	    eventMask[iy][ix]=PHOTON_MAX;
+	  } else {
+	    eventMask[iy][ix]=PHOTON;
+	  }
+	} else if (eventMask[iy][ix]==PEDESTAL) {
+	  if (cm==0)
+	    addToPedestal(det->getValue(data, ix, iy),ix,iy);
+	}
+      
+	
 
       return  eventMask[iy][ix];