diff --git a/slsDetectorCalibration/MovingStat.h b/slsDetectorCalibration/MovingStat.h
index 9c8584d33..d856228dd 100755
--- a/slsDetectorCalibration/MovingStat.h
+++ b/slsDetectorCalibration/MovingStat.h
@@ -7,11 +7,15 @@
 class MovingStat
     {
     public:
-      MovingStat() : n(1000), m_n(0) {}
+      MovingStat(int nn=1000) : n(nn), m_n(0) {}
 
         void Clear()
         {
-            m_n = 0;
+	  m_n = 0; 
+	  m_oldM=0;  
+	  m_newM=0; 
+	  m_oldM2=0; 
+	  m_newM2=0;
         }
 	
 	void SetN(int i) {n=i;};
diff --git a/slsDetectorCalibration/commonModeSubtraction.h b/slsDetectorCalibration/commonModeSubtraction.h
new file mode 100644
index 000000000..89984ed81
--- /dev/null
+++ b/slsDetectorCalibration/commonModeSubtraction.h
@@ -0,0 +1,53 @@
+#ifndef COMMONMODESUBTRACTION_H
+#define COMMONMODESUBTRACTION_H
+
+#include "MovingStat.h"
+
+class commonModeSubtraction {
+
+ public:
+  commonModeSubtraction(int nn=1000, int iroi=1) : cmStat(NULL), cmPed(NULL), nCm(NULL), nROI(iroi)  {cmStat=new MovingStat[nROI]; for (int i=0; i<nROI; i++) cmStat[i].SetN(nn); cmPed=new double[nROI]; nCm=new double[nROI];};
+    virtual ~commonModeSubtraction() {delete [] cmStat; delete [] cmPed; delete [] nCm;};
+    
+    virtual void  Clear(){
+      for (int i=0; i<nROI; i++) {
+	cmStat[i].Clear();  
+	nCm[i]=0; 
+	cmPed[i]=0;
+      }};
+    
+    virtual void  newFrame(){
+      for (int i=0; i<nROI; i++) {
+	if (nCm[i]>0) cmStat[i].Calc(cmPed[i]/nCm[i]);  
+	nCm[i]=0; 
+	cmPed[i]=0;
+      }};
+    
+    virtual void addToCommonMode(double val, int ix=0, int iy=0) {
+      (void)ix; (void)iy;
+      cmPed[0]+=val; 
+      nCm[0]++;};
+    
+    virtual double getCommonMode(int ix=0, int iy=0) {
+      (void)ix; (void)iy;
+      if (nCm[0]>0) return cmPed[0]/nCm[0]-cmStat[0].Mean(); 
+      else return 0;};
+    
+
+
+
+
+ protected:
+  MovingStat *cmStat; //stores the common mode average per supercolumn */
+  double *cmPed; // stores the common mode for this frame per supercolumn */
+  double *nCm; // stores the number of pedestals to calculate the cm per supercolumn */
+  const int nROI;
+
+
+};
+
+
+
+
+
+#endif
diff --git a/slsDetectorCalibration/moench02ModuleData.h b/slsDetectorCalibration/moench02ModuleData.h
index 8a9a2d579..1791c74e0 100644
--- a/slsDetectorCalibration/moench02ModuleData.h
+++ b/slsDetectorCalibration/moench02ModuleData.h
@@ -33,38 +33,29 @@ class moench02ModuleData : public slsDetectorData<uint16_t> {
   */
   
 
-  moench02ModuleData(int nbg=1000): slsDetectorData<uint16_t>(160, 160, 40, 1286), nSigma(5), buff(NULL), oldbuff(NULL), pedSub(1) {
-    int headerOff=0;
-    int footerOff=1284;
-    int dataOff=4;
+  moench02ModuleData(double c=0): slsDetectorData<uint16_t>(160, 160, 40, 1286), xtalk(c) {
 
-    int ix, iy, ir, ic, i, isc, ip;
 
-    int **dMask;
+
+
+    uint16_t **dMask;
     int **dMap;
-
-    for (ir=0; ir<160; ir++) {
-      for (ic=0; ic<160; ic++) {
-	stat[ir][ic].Clear();
-	stat[ir][ic].SetN(nbg);
-      }
-    }
+    int ix, iy;
 
 
-    dMask=new int*[160];
-    for(i = 0; i < 160; i++) {
-      dMask[i] = new int[160];
-    }
+
+    dMask=new uint16_t*[160];
     dMap=new int*[160];
-    for(i = 0; i < 160; i++) {
+    for (int i = 0; i < 160; i++) {
       dMap[i] = new int[160];
+      dMask[i] = new uint16_t[160];
     }
 
-    for (isc=0; isc<4; isc++) {
-      for (ip=0; ip<10; ip++) {
+    for (int isc=0; isc<4; isc++) {
+      for (int ip=0; ip<10; ip++) {
 
-	for (ir=0; ir<16; ir++) {
-	  for (ic=0; ic<40; ic++) {
+	for (int ir=0; ir<16; ir++) {
+	  for (int ic=0; ic<40; ic++) {
 
 	    ix=isc*40+ic;
 	    iy=ip*16+ir;
@@ -88,9 +79,13 @@ class moench02ModuleData : public slsDetectorData<uint16_t> {
 
     setDataMap(dMap);
     setDataMask(dMask);
+
+
+    
+
   };
     
-    ~moench02ModuleData() {if (buff) delete [] buff; if (oldbuff) delete [] oldbuff; };
+    //   ~moench02ModuleData() {if (buff) delete [] buff; if (oldbuff) delete [] oldbuff; };
   
 
   int getPacketNumber(char *buff){
@@ -100,176 +95,27 @@ class moench02ModuleData : public slsDetectorData<uint16_t> {
     return np;
   };
 
-
-  void calculateCommonMode(int xmin=0, int xmax=160, int ymin=0, int ymax=160, double hc=0, double tc=0) {
-
-    int scmin=xmin/40;
-    int scmax=(xmax-1)/40+1;
-    int isc=0, ix, iy;
-
-    for (isc=scmin; isc<scmax; isc++) {
-      nCm[isc]=0;
-      cmPed[isc]=0;
-    }
-    
-    
-
-
-    for (ix=xmin-1; ix<xmax+1; ix++) {
-      isc=ix/40;
-      for (iy=ymin-1; iy<ymax+1; iy++) {
-	if (isGood(ix,iy)) {
-	  if (getEventType(ix, iy, hc, tc, 1,0)==PEDESTAL) {
-	    cmPed[isc]+=getChannelShort(ix, iy, hc, tc);
-	    nCm[isc]++;
-	  }
-	}
-	
-      }
-
-    }
-
-    
-    for (isc=scmin; isc<scmax; isc++) {
-
-      if (nCm[isc]>0)
-	cmStat[isc].Calc(cmPed[isc]/nCm[isc]);
-
-      //   cout << "SC " << isc << nCm[isc] << " " << cmPed[isc]/nCm[isc] << " " << cmStat[isc].Mean() << endl;
-
-    }
-    
-
-  };
-
-  double getCommonMode(int ix, int iy) {
-    int isc=ix/40;
-    if (nCm[isc]>0) {
-     /*  if (ix==20 && iy==20) */
-/* 	cout << cmPed[isc] << " " << nCm[isc] << " " << cmStat[isc].Mean() << " " << cmPed[isc]/nCm[isc]-cmStat[isc].Mean() << endl; */
-      return cmPed[isc]/nCm[isc]-cmStat[isc].Mean();
-    } else {
-  /*     if (ix==20 && iy==20) */
-/* 	cout << "No common mode!" << endl; */
-      
-      return 0;
-    }
-  };
-
-
-
-
-  void addToPedestal(double val, int ix, int iy){stat[iy][ix].Calc(val);};
-  double getPedestal(int ix, int iy){return stat[iy][ix].Mean();};
-  double getPedestalRMS(int ix, int iy){return stat[iy][ix].StandardDeviation();};
-  double setNSigma(double n){if (n>0) nSigma=n; return nSigma;}
-  char *getBuff(){return buff;}
-  char *getOldBuff(){return oldbuff;}
-
-  int setPedestalSubstraction(int i=-1){if (i>=0) pedSub=i; return pedSub;};
-
-
-  double getChannelShort(int ix, int iy, double hc=0, double tc=0) {
-    
-    double val=slsDetectorData<uint16_t>::getChannel(buff,ix,iy);
-
-    if (ix>0 && hc!=0)
-      val-=hc*slsDetectorData<uint16_t>::getChannel(buff,ix-1,iy);
-    
-    if (oldbuff && tc!=0) 
-      val+=tc*slsDetectorData<uint16_t>::getChannel(oldbuff,ix,iy);
-
-
-    return val;
+  double getValue(char *data, int ix, int iy=0) {
+    if (xtalk==0 || ix%40==0)
+      return (double)getValue(data, ix, iy);
+    else
+      return slsDetectorData<uint16_t>::getValue(data, ix, iy)-xtalk*slsDetectorData<uint16_t>::getValue(data, ix-1, iy);
   };
   
+  double setXTalk(double c) {xtalk=c; return xtalk;}
+  double getXTalk() {return xtalk;}
+  
 
 
-  eventType getEventType(int ix, int iy, double hcorr=0, double tcorr=0, int sign=1, int cm=0) {
 
- /*    PEDESTAL, */
-/*     NEIGHBOUR, */
-/*     PHOTON, */
-/*     PHOTON_MAX, */
- /*    NEGATIVE_PEDESTAL */
-      eventType ret=PEDESTAL;
-      double tot=0, v, max=0, sig;
-      
-
-      sig=getPedestalRMS(ix,iy);
-
-      for (int ir=-1; ir<2; ir++) {
-	for (int ic=-1; ic<2; ic ++) {
-	  if ((iy+ir)>=0 && (iy+ir)<160 && (ix+ic)>=0 && (ix+ic)<160) {
-	    v=sign*getChannelShort(ix+ic, iy+ir, hcorr, tcorr);
-	    if (pedSub)
-	      v-=sign*getPedestal(ix+ic,iy+ir);
-	    if (cm)
-	      v-=sign*getCommonMode(ix+ic, iy+ic);
-	    cluster[ic+1][ir+1]=v;
-	    tot+=v;
-	    if (v>max) {
-	      max=v;
-	    }
-	    if (ir==0 && ic==0) {
-	      if (v>nSigma*getPedestalRMS(ix,iy)) {
-		ret=PHOTON;
-	      } else if (v<-nSigma*getPedestalRMS(ix,iy))
-	       ret=NEGATIVE_PEDESTAL;
-	    }
-	  }
-	}
-      }
-      
-      if (ret!=PHOTON && tot>3*nSigma*sig)
-	ret=NEIGHBOUR;
-      else if (ret==PHOTON) {
-	if (cluster[1][1]>=max)
-	  ret=PHOTON_MAX;
-      }
-      
-   
-      cx=ix;
-      cy=iy;
-
-      return ret;
-
-  };
-
-
-  double getClusterElement(int ic, int ir, int cmsub=0){return cluster[ic+1][ir+1];};
-  
-  double  *getCluster(){return &cluster[0][0];};
-
-
-  TTree *initMoenchTree(char *tname, Int_t *iframe, Int_t * x, Int_t* y, Double_t data[3][3], Double_t *ped, Double_t *sigma) {
-    TTree* tall=new TTree(tname,tname);
-    tall->Branch("iFrame",iframe,"iframe/I");
-    tall->Branch("x",x,"x/I");
-    tall->Branch("y",y,"y/I");
-    tall->Branch("data",data,"data[3][3]/D");
-    tall->Branch("pedestal",ped,"pedestal/D");
-    tall->Branch("rms",sigma,"rms/D");
-    return tall;
-
-  };
 
   
 
  private:
   
-  MovingStat stat[160][160];
-  double nSigma;
-  double cluster[3][3];
-  char *buff;
-  char *oldbuff;
-  int pedSub;
-  int cmSub;
+  double xtalk;
+
 
-  MovingStat cmStat[4]; //stores the common mode average per supercolumn
-  double cmPed[4]; // stores the common mode for this frame per supercolumn
-  double nCm[4]; // stores the number of pedestals to calculate the cm per supercolumn
-  int cx, cy;
 };
 
 
diff --git a/slsDetectorCalibration/moenchCommonMode.h b/slsDetectorCalibration/moenchCommonMode.h
new file mode 100644
index 000000000..9a3a45c85
--- /dev/null
+++ b/slsDetectorCalibration/moenchCommonMode.h
@@ -0,0 +1,34 @@
+#ifndef MOENCHCOMMONMODE_H
+#define MOENCHCOMMONMODE_H
+
+#include "commonModeSubtraction.h"
+
+class moenchCommonMode : public commonModeSubtraction {
+
+ public:
+  moenchCommonMode(int nn=1000) : commonModeSubtraction(nn,4){} ;
+   
+
+    virtual void addToCommonMode(double val, int ix=0, int iy=0) {
+      (void)iy;
+      int isc=ix/40;
+      if (isc>=0 && isc<4) {
+	cmPed[isc]+=val; 
+	nCm[isc]++;
+      }
+
+    };
+    
+    virtual double getCommonMode(int ix=0, int iy=0) {
+      (void)iy;
+      int isc=ix/40;
+      if (isc>=0 && isc<4) {
+	if (nCm[isc]>0) return cmPed[isc]/nCm[isc]-cmStat[isc].Mean(); 
+      }
+      return 0;};
+
+
+};
+
+
+#endif
diff --git a/slsDetectorCalibration/moenchReadData.C b/slsDetectorCalibration/moenchReadData.C
index 468817e45..09388fe4a 100644
--- a/slsDetectorCalibration/moenchReadData.C
+++ b/slsDetectorCalibration/moenchReadData.C
@@ -15,6 +15,8 @@
 //#include <queue>
 #include <fstream>
 #include "moench02ModuleData.h"
+#include "moenchCommonMode.h"
+#include "singlePhotonDetector.h"
 
 //#include "MovingStat.h"
 
@@ -51,13 +53,19 @@ using namespace std;
 */
 
 
-THStack *moenchReadData(char *fformat, char *tit, int runmin, int runmax, int nbins=1500, int hmin=-500, int hmax=1000, int sign=1, double hc=0, double tc=0, int xmin=0, int xmax=NC, int ymin=0, int ymax=NR, int pedsub=1, int cmsub=0) {
+THStack *moenchReadData(char *fformat, char *tit, int runmin, int runmax, int nbins=1500, int hmin=-500, int hmax=1000, int sign=1, double hc=0, int xmin=0, int xmax=NC, int ymin=0, int ymax=NR, int cmsub=0) {
   
-  moench02ModuleData *decoder=new moench02ModuleData();
+  moench02ModuleData *decoder=new moench02ModuleData(hc);
+  moenchCommonMode *cmSub=NULL;
+  if (cmsub)
+    cmSub=new moenchCommonMode();
+
+    
+
+  singlePhotonDetector<uint16_t> *filter=new singlePhotonDetector<uint16_t>(decoder, 3, 5, sign, cmSub);
+
   char *buff;
-  char *oldbuff=NULL;
   char fname[10000];
-  double oldval;
   int nf=0;
   
   eventType thisEvent=PEDESTAL;
@@ -80,50 +88,39 @@ THStack *moenchReadData(char *fformat, char *tit, int runmin, int runmax, int nb
   TH2F *h1=new TH2F("h1",tit,nbins,hmin-0.5,hmax-0.5,NC*NR,-0.5,NC*NR-0.5);
   hs->Add(h1);
 
-  if (pedsub) {
 
-    h2=new TH2F("h2",tit,nbins,hmin-0.5,hmax-0.5,NC*NR,-0.5,NC*NR-0.5);
-    h3=new TH2F("h3",tit,nbins,hmin-0.5,hmax-0.5,NC*NR,-0.5,NC*NR-0.5);
-    hetaX=new TH2F("hetaX",tit,nbins,-1,2,NC*NR,-0.5,NC*NR-0.5);
-    hetaY=new TH2F("hetaY",tit,nbins,-1,2,NC*NR,-0.5,NC*NR-0.5);
-    hs->Add(h2);
-    hs->Add(h3);
-    hs->Add(hetaX);
-    hs->Add(hetaY);
-  }
+  h2=new TH2F("h2",tit,nbins,hmin-0.5,hmax-0.5,NC*NR,-0.5,NC*NR-0.5);
+  h3=new TH2F("h3",tit,nbins,hmin-0.5,hmax-0.5,NC*NR,-0.5,NC*NR-0.5);
+  hetaX=new TH2F("hetaX",tit,nbins,-1,2,NC*NR,-0.5,NC*NR-0.5);
+  hetaY=new TH2F("hetaY",tit,nbins,-1,2,NC*NR,-0.5,NC*NR-0.5);
+  hs->Add(h2);
+  hs->Add(h3);
+  hs->Add(hetaX);
+  hs->Add(hetaY);
 
   
-
+  
   ifstream filebin;
 
 
   int ix=20, iy=20, ir, ic;
 
 
-  Int_t iFrame, x, y;
-  Double_t data[3][3],ped,sigma;
+  Int_t iFrame;
+  TTree *tall=filter->initEventTree(tit, &iFrame);
+
 
-  TTree *tall;
-  if (pedsub) {
-    cout << "Subtracting pedestal!!!! " << endl;
-    tall=decoder->initMoenchTree(tit, &iFrame, &x, &y, data, &ped, &sigma);
-  }
   double tot, tl, tr, bl, br, v;
   
 
-  int scmin=xmin/40, scmax=(xmax-1)/40+1;
-
-  // 6% x-talk from previous pixel
-  // 12% x-talk from previous frame
-
 #ifdef MY_DEBUG
   TCanvas *myC=new TCanvas();
-  TH2F *he=new TH2F("he","Event",3,-1.5,1.5,3,-1.5,1.5);
+  TH2F *he=new TH2F("he","Event Mask",xmax-xmin, xmin, xmax, ymax-ymin, ymin, ymax);
   he->SetStats(kFALSE);
   he->Draw("colz");
-  he->SetMinimum(0);
-  he->SetMaximum(0.5*hmax);
 #endif
+  filter->newDataSet();
+
 
   for (int irun=runmin; irun<runmax; irun++) {
     sprintf(fname,fformat,irun);
@@ -133,77 +130,45 @@ THStack *moenchReadData(char *fformat, char *tit, int runmin, int runmax, int nb
     while ((buff=decoder->readNextFrame(filebin))) {
 
 
+      filter->newFrame();
 
-      if (nf>100) {
-	if (cmsub) {
-	  for (int isc=scmin; isc<scmax; isc++) {
-	    decoder->calculateCommonMode(3+isc*40, 40*(isc+1)-3, 3, NR-3, hc, tc);
-#ifdef MY_DEBUG
-	    if (nf%1000==0) cout << "sc=" << isc << " CM="<< decoder->getCommonMode(3+isc*40, NR/2)<< endl;
-#endif
+	//calculate pedestals and common modes
+      if (cmsub) {
+	for (ix=xmin-1; ix<xmax+1; ix++)
+	  for (iy=ymin-1; iy<ymax+1; iy++) {
+	    thisEvent=filter->getEventType(buff, ix, iy,0);
 	  }
-	}
-	
       }
-      
+
 
 
 
       for (ix=xmin-1; ix<xmax+1; ix++)
 	for (iy=ymin-1; iy<ymax+1; iy++) {
-	  ////////////////////////////////////////////////////////
-	  if (pedsub) {
+	  thisEvent=filter->getEventType(buff, ix, iy, cmsub);
 
-
-	    if (nf>100) {
-	      thisEvent= decoder->getEventType(ix, iy, hc, tc, 1,1);
-	    } 
-	    
+#ifdef MY_DEBUG
+	  if (iev%10000==0)
+	    he->SetBinContent(ix+1-xmin, iy+1-ymin, (int)thisEvent);
+#endif	  
 	  
-
-
-	    if (thisEvent==PEDESTAL) {
-	      if (cmsub && nf>1000)
-		decoder->addToPedestal(decoder->getChannelShort(ix, iy, hc, tc)-decoder->getCommonMode(ix,iy), ix, iy); 
-	      else
-		decoder->addToPedestal( decoder->getChannelShort(ix, iy, hc, tc), ix, iy); 
-	     
-	     
-	    }
-	    
-	      
-	      /***********************************************************
-Add  here the function that you want to call: fill histos, make trees etc.
-	      ************************************************************/
-	      // getClusterElement to access quickly the photon and the neighbours
 	  if (nf>1000) {
 	    tot=0;
 	    tl=0;
 	    tr=0;
 	    bl=0;
 	    br=0;
-	    h1->Fill(decoder->getClusterElement(0,0), iy+NR*ix);
+	    h1->Fill(filter->getClusterElement(0,0), iy+NR*ix);
 	    
 	    //  if (nf%1000==0 && ix==20 && iy==20) cout <<  " val="<< decoder->getClusterElement(0,0)<< endl;
 
 	    if (thisEvent==PHOTON_MAX ) {
-	      #ifdef MY_DEBUG
-		  if (iev%100000==0) {
-		  cout << "Event " << iev << " Frame "<< nf << endl;
-		  }
-#endif		  
-	
-	      
+     
 	      for (ir=-1; ir<2; ir++) {
 		for (ic=-1; ic<2; ic++) {
-		  v=decoder->getClusterElement(ic,ir,cmsub);
-		  data[ic+1][ir+1]=v;
-#ifdef MY_DEBUG
-		  if (iev%100000==0) {
-		    he->SetBinContent(ic+2,ir+2,v);
-		    cout << "Histo("<< ix+ic << ","<< iy+ir <<")" << v << endl;
-		  }
-#endif		  
+		  v=filter->getClusterElement(ic,ir);
+		  //	  data[ic+1][ir+1]=v;
+  
 		  
 
 		  tot+=v;
@@ -233,34 +198,34 @@ Add  here the function that you want to call: fill histos, make trees etc.
 	      if (bl>br && bl>tl && bl>tr) {
 		  h2->Fill(bl, iy+NR*ix);
 		if (bl>0) {
-		  hetaX->Fill((decoder->getClusterElement(0,0,cmsub)+decoder->getClusterElement(0,-1,cmsub))/bl,iy+NR*ix);
-		  hetaY->Fill((decoder->getClusterElement(0,0,cmsub)+decoder->getClusterElement(-1,0,cmsub))/bl,iy+NR*ix);
-		  hetaX->Fill(1.-(decoder->getClusterElement(0,0,cmsub)+decoder->getClusterElement(0,-1,cmsub))/bl,iy+NR*(ix-1));
-		  hetaY->Fill(1.-(decoder->getClusterElement(0,0,cmsub)+decoder->getClusterElement(-1,0,cmsub))/bl,(iy-1)+NR*ix);
+		  hetaX->Fill((filter->getClusterElement(0,0)+filter->getClusterElement(0,-1))/bl,iy+NR*ix);
+		  hetaY->Fill((filter->getClusterElement(0,0)+filter->getClusterElement(-1,0))/bl,iy+NR*ix);
+		  hetaX->Fill(1.-(filter->getClusterElement(0,0)+filter->getClusterElement(0,-1))/bl,iy+NR*(ix-1));
+		  hetaY->Fill(1.-(filter->getClusterElement(0,0)+filter->getClusterElement(-1,0))/bl,(iy-1)+NR*ix);
 		}
 	      } else if (br>bl && br>tl && br>tr) {
 		  h2->Fill(br, iy+NR*ix);
 		if (br>0) {
-		  hetaX->Fill((decoder->getClusterElement(0,0,cmsub)+decoder->getClusterElement(0,-1,cmsub))/br,iy+NR*ix);
-		  hetaY->Fill((decoder->getClusterElement(0,0,cmsub)+decoder->getClusterElement(1,0,cmsub))/br,iy+NR*ix);
-		  hetaX->Fill(1.-(decoder->getClusterElement(0,0,cmsub)+decoder->getClusterElement(0,-1,cmsub))/br,iy+NR*(ix+1));
-		  hetaY->Fill(1.-(decoder->getClusterElement(0,0,cmsub)+decoder->getClusterElement(1,0,cmsub))/br,iy-1+NR*ix);
+		  hetaX->Fill((filter->getClusterElement(0,0)+filter->getClusterElement(0,-1))/br,iy+NR*ix);
+		  hetaY->Fill((filter->getClusterElement(0,0)+filter->getClusterElement(1,0))/br,iy+NR*ix);
+		  hetaX->Fill(1.-(filter->getClusterElement(0,0)+filter->getClusterElement(0,-1))/br,iy+NR*(ix+1));
+		  hetaY->Fill(1.-(filter->getClusterElement(0,0)+filter->getClusterElement(1,0))/br,iy-1+NR*ix);
 		}
 	      } else if (tl>br && tl>bl && tl>tr) {
 		  h2->Fill(tl, iy+NR*ix);
 		if (tl>0) {
-		  hetaX->Fill((decoder->getClusterElement(0,0,cmsub)+decoder->getClusterElement(0,1,cmsub))/tl,iy+NR*ix);
-		  hetaY->Fill((decoder->getClusterElement(0,0,cmsub)+decoder->getClusterElement(-1,0,cmsub))/tl,iy+NR*ix);
-		  hetaX->Fill(1.-(decoder->getClusterElement(0,0,cmsub)+decoder->getClusterElement(0,1,cmsub))/tl,iy+NR*(ix-1));
-		  hetaY->Fill(1.-(decoder->getClusterElement(0,0,cmsub)+decoder->getClusterElement(-1,0,cmsub))/tl,iy+1+NR*ix);
+		  hetaX->Fill((filter->getClusterElement(0,0)+filter->getClusterElement(0,1))/tl,iy+NR*ix);
+		  hetaY->Fill((filter->getClusterElement(0,0)+filter->getClusterElement(-1,0))/tl,iy+NR*ix);
+		  hetaX->Fill(1.-(filter->getClusterElement(0,0)+filter->getClusterElement(0,1))/tl,iy+NR*(ix-1));
+		  hetaY->Fill(1.-(filter->getClusterElement(0,0)+filter->getClusterElement(-1,0))/tl,iy+1+NR*ix);
 		}
 	      } else if (tr>bl && tr>tl && tr>br) {
 		  h2->Fill(tr, iy+NR*ix);
 		if (tr>0) {
-		  hetaX->Fill((decoder->getClusterElement(0,0,cmsub)+decoder->getClusterElement(0,1,cmsub))/tr,iy+NR*ix);
-		  hetaY->Fill((decoder->getClusterElement(0,0,cmsub)+decoder->getClusterElement(1,0,cmsub))/tr,iy+NR*ix);
-		  hetaX->Fill(1.-(decoder->getClusterElement(0,0,cmsub)+decoder->getClusterElement(0,1,cmsub))/tr,iy+NR*(ix+1));
-		  hetaY->Fill(1.-(decoder->getClusterElement(0,0,cmsub)+decoder->getClusterElement(1,0,cmsub))/tr,iy+1+NR*ix);
+		  hetaX->Fill((filter->getClusterElement(0,0)+filter->getClusterElement(0,1))/tr,iy+NR*ix);
+		  hetaY->Fill((filter->getClusterElement(0,0)+filter->getClusterElement(1,0))/tr,iy+NR*ix);
+		  hetaX->Fill(1.-(filter->getClusterElement(0,0)+filter->getClusterElement(0,1))/tr,iy+NR*(ix+1));
+		  hetaY->Fill(1.-(filter->getClusterElement(0,0)+filter->getClusterElement(1,0))/tr,iy+1+NR*ix);
 		}
 
 
@@ -269,10 +234,6 @@ Add  here the function that you want to call: fill histos, make trees etc.
 
 		h3->Fill(tot, iy+NR*ix);
 		iFrame=decoder->getFrameNumber(buff);
-		ped=decoder->getPedestal(ix,iy);
-		sigma=decoder->getPedestalRMS(ix,iy);
-		x=ix;
-		y=iy;
 
 		tall->Fill();
 	   
@@ -290,29 +251,25 @@ Add  here the function that you want to call: fill histos, make trees etc.
 	    
 	    
 	  }
-	  } else {
-	    //cout << ix << " " << iy << endl;
-	    h1->Fill(decoder->getChannelShort(ix, iy), iy+NR*ix);
-	  }
+	} 
 	  //////////////////////////////////////////////////////////
-	}
+    
 
       cout << "=" ;
       nf++;
-	}
+    }
     cout << endl;
     if (filebin.is_open())
       filebin.close();	  
     else
       cout << "could not open file " << fname << endl;
-    }
-  if (pedsub)
-    tall->Write(tall->GetName(),TObject::kOverwrite);
-
+  }
+  tall->Write(tall->GetName(),TObject::kOverwrite);
+  
   
   
   delete decoder;
   cout << "Read " << nf << " frames" << endl;
   return hs;
-  }
+}
 
diff --git a/slsDetectorCalibration/pedestalSubtraction.h b/slsDetectorCalibration/pedestalSubtraction.h
new file mode 100644
index 000000000..f80991ec0
--- /dev/null
+++ b/slsDetectorCalibration/pedestalSubtraction.h
@@ -0,0 +1,24 @@
+#ifndef PEDESTALSUBTRACTION_H
+#define PEDESTALSUBTRACTION_H
+
+#include "MovingStat.h"
+
+class pedestalSubtraction  {
+
+ public:
+  pedestalSubtraction (int nn=1000) : stat(nn) {};
+    virtual ~pedestalSubtraction() {};
+    
+    virtual void Clear() {stat.Clear();}
+    virtual void addToPedestal(double val){ stat.Calc(val);};
+    virtual double getPedestal(){return stat.Mean();};
+    virtual double getPedestalRMS(){return stat.StandardDeviation();};
+    virtual int SetNPedestals(int i=-1) {if (i>0) stat.SetN(i); return stat.GetN();};
+    
+
+  
+ private:
+  MovingStat stat;
+
+};
+#endif
diff --git a/slsDetectorCalibration/singlePhotonDetector.h b/slsDetectorCalibration/singlePhotonDetector.h
index 2c226a69e..886896971 100644
--- a/slsDetectorCalibration/singlePhotonDetector.h
+++ b/slsDetectorCalibration/singlePhotonDetector.h
@@ -1,9 +1,22 @@
 #ifndef SINGLEPHOTONDETECTOR_H
 #define  SINGLEPHOTONDETECTOR_H
+
+
+#include <inttypes.h>
 #include "slsDetectorData.h"
 
-#include "MovingStat.h"
 #include "single_photon_hit.h"
+#include "pedestalSubtraction.h"
+#include "commonModeSubtraction.h"
+
+
+#define MYROOT1
+
+#ifdef MYROOT1
+#include <TTree.h>
+
+#endif
+
 
 #include <iostream>
 
@@ -11,11 +24,12 @@ using namespace std;
 
 
   enum eventType {
-    PEDESTAL,
-    NEIGHBOUR,
-    PHOTON,
-    PHOTON_MAX,
-    NEGATIVE_PEDESTAL,
+    PEDESTAL=0,
+    NEIGHBOUR=1,
+    PHOTON=2,
+    PHOTON_MAX=3,
+    NEGATIVE_PEDESTAL=4,
+    UNDEFINED=-1
   };
 
 
@@ -40,7 +54,13 @@ class singlePhotonDetector {
   */
   
 
-  singlePhotonDetector(slsDetectorData<dataType> *d, int csize=3, double nsigma=5) : det(d), clusterSize(csize), clusterSizeY(csize),nSigma(nsigma) {
+  singlePhotonDetector(slsDetectorData<dataType> *d, 
+		       int csize=3, 
+		       double nsigma=5,  
+		       int sign=1, 
+		       commonModeSubtraction *cm=NULL,
+		       int nped=1000, 
+		       int nd=100) : det(d), nx(0), ny(0), stat(NULL), cmSub(cm),  nDark(nd), eventMask(NULL),nSigma (nsigma), clusterSize(csize), clusterSizeY(csize), cluster(NULL), iframe(-1), dataSign(sign) {
 
    
 
@@ -48,9 +68,13 @@ class singlePhotonDetector {
     
 
 
-    stat=new MovingStat*[ny];
-    for (int i=0; i<ny; i++)
-      stat[i]=new MovingStat;
+    stat=new pedestalSubtraction*[ny];
+    eventMask=new eventType*[ny];
+    for (int i=0; i<ny; i++) {
+      stat[i]=new pedestalSubtraction[nx];
+      stat[i]->SetNPedestals(nped);
+      eventMask[i]=new eventType[nx];
+    }
    
     if (ny==1)
       clusterSizeY=1;
@@ -59,13 +83,24 @@ class singlePhotonDetector {
     
   };
     
-    ~singlePhotonDetector() {delete cluster; delete [] stat;};
-  
+    virtual ~singlePhotonDetector() {delete cluster; for (int i=0; i<ny; i++) delete [] stat[i]; delete [] stat;};
 
 
-    void addToPedestal(double val, int ix, int iy){ if (ix>=0 && ix<nx && iy>=0 && iy<ny) stat[iy][ix].Calc(val);};
-    double getPedestal(int ix, int iy){if (ix>=0 && ix<nx && iy>=0 && iy<ny) return stat[iy][ix].Mean(); else return -1;};
-    double getPedestalRMS(int ix, int iy){if (ix>=0 && ix<nx && iy>=0 && iy<ny) return stat[iy][ix].StandardDeviation();else return -1;};
+    void newDataSet(){iframe=-1; for (int iy=0; iy<ny; iy++) for (int ix=0; ix<nx; ix++) stat[iy][ix].Clear(); if (cmSub) cmSub->Clear(); };  
+    void newFrame(){iframe++; for (int iy=0; iy<ny; iy++) for (int ix=0; ix<nx; ix++) eventMask[iy][ix]=UNDEFINED; if (cmSub) cmSub->newFrame();};
+
+
+    commonModeSubtraction setCommonModeSubtraction(commonModeSubtraction *cm) {cmSub=cm; return cmSub;};
+
+    int setDataSign(int sign=0) {if (sign==1 || sign==-1) dataSign=sign; return dataSign;};
+
+    virtual void addToPedestal(double val, int ix, int iy){ if (ix>=0 && ix<nx && iy>=0 && iy<ny) { stat[iy][ix].addToPedestal(val); if (cmSub && det->isGood(ix, iy) ) cmSub->addToCommonMode(val, ix, iy);}};
+
+
+    virtual double getPedestal(int ix, int iy, int cm=0){if (ix>=0 && ix<nx && iy>=0 && iy<ny) if (cmSub && cm>0) return stat[iy][ix].getPedestal()-cmSub->getCommonMode(); else return stat[iy][ix].getPedestal(); else return -1;};
+
+
+    double getPedestalRMS(int ix, int iy){if (ix>=0 && ix<nx && iy>=0 && iy<ny) return stat[iy][ix].getPedestalRMS();else return -1;};
   
     double setNSigma(double n=-1){if (n>0) nSigma=n; return nSigma;}
  
@@ -85,49 +120,74 @@ class singlePhotonDetector {
     
 
 
-    eventType getEventType(char *data, int ix, int iy, double hcorr=0, double tcorr=0, int sign=1, int cm=0) {
+    eventType getEventType(char *data, int ix, int iy, int cm=0) {
 
-      eventType ret=PEDESTAL;
-      double tot=0, v, max=0, sig;
+      // eventType ret=PEDESTAL;
+      double tot=0, max=0;
+  
+      if (iframe<nDark) {
+	if (cm==0)
+	  addToPedestal(det->getValue(data, ix, iy),ix,iy);
+	return UNDEFINED;
+      }
+      
       
 
-      sig=getPedestalRMS(ix,iy);
+      if (eventMask[iy][ix]==UNDEFINED) {
+	
+	eventMask[iy][ix]=PEDESTAL;
+	
+	
+	cluster->x=ix;
+	cluster->y=iy;
+	cluster->rms=getPedestalRMS(ix,iy);
+	cluster->ped=getPedestal(ix,iy, cm);
+	
 
-      cluster->x=ix;
-      cluster->y=iy;
-      cluster->rms=sig;
-      cluster->ped=getPedestal(ix,iy);
-      
-
-      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)<160 && (ix+ic)>=0 && (ix+ic)<160) {
-	    v=sign*(det->getValue(data, ix+ic, iy+ir)-getPedestal(ix+ic,iy+ir));
-	    // if (cm)
-	    //  v-=sign*getCommonMode(ix+ic, iy+ic);
-	    cluster->set_data(ic, ir, v);
-	    tot+=v;
-	    if (v>max) {
-	      max=v;
-	    }
-	    if (ir==0 && ic==0) {
-	      if (v>nSigma*getPedestalRMS(ix,iy)) {
-		ret=PHOTON;
-	      } else if (v<-nSigma*getPedestalRMS(ix,iy))
-		ret=NEGATIVE_PEDESTAL;
+	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(dataSign*(det->getValue(data, ix+ic, iy+ir)-getPedestal(ix+ic,iy+ir,cm)), ic, ir       );
+	      tot+=cluster->get_data(ic,ir);
+	      if (cluster->get_data(ic,ir)>max) {
+		max=cluster->get_data(ic,ir);
+	      }
+	      if (ir==0 && ic==0) {
+		if (cluster->get_data(ic,ir)>nSigma*cluster->rms) {
+		  eventMask[iy][ix]=PHOTON;
+		} else if (cluster->get_data(ic,ir)<-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) {
+	    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) {
+		  if (eventMask[iy+ir][ix+ic]==UNDEFINED) 
+		    if (ir==0 && ic==0)   eventMask[iy+ir][ix+ic]=PHOTON_MAX;
+		    else eventMask[iy+ir][ix+ic]=NEIGHBOUR;
+		}
+		
+	      }
+	    }
+	    
+
+	  }
+	} else if (eventMask[iy][ix]==PEDESTAL) {
+	  if (cm==0)
+	    addToPedestal(det->getValue(data, ix, iy),ix,iy);
+	}
       }
-      
-      if (ret!=PHOTON && tot>3*nSigma*sig)
-	ret=NEIGHBOUR;
-      else if (ret==PHOTON) {
-	if (cluster->get_data(0,0)>=max)
-	  ret=PHOTON_MAX;
-      }
-      
-      return ret;
+           
+      return  eventMask[iy][ix];
 
   };
 
@@ -135,10 +195,31 @@ class singlePhotonDetector {
 
 
 
+    int SetNPedestals(int i=-1) {int ix=0, iy=0; if (i>0) for (ix=0; ix<nx; ix++) for (iy=0; iy<ny; iy++) stat[iy][ix].SetNPedestals(i); return stat[0][0].SetNPedestals();};
 
+    double getClusterElement(int ic, int ir){return cluster->get_data(ic,ir);};
+    eventType getEventMask(int ic, int ir){return eventMask[ir][ic];};
+ 
+#ifdef MYROOT1
+    TTree *initEventTree(char *tname, int *iFrame=NULL) {
+      TTree* tall=new TTree(tname,tname);
+
+      if (iFrame)
+	tall->Branch("iFrame",iFrame,"iframe/I");
+      else
+	tall->Branch("iFrame",&(cluster->iframe),"iframe/I");
+
+      tall->Branch("x",&(cluster->x),"x/I");
+      tall->Branch("y",&(cluster->y),"y/I");
+      char tit[100];
+      sprintf(tit,"data[%d]/D",clusterSize*clusterSizeY);
+      tall->Branch("data",cluster->data,tit);
+      tall->Branch("pedestal",&(cluster->ped),"pedestal/D");
+      tall->Branch("rms",&(cluster->rms),"rms/D");
+      return tall;
+    };
+#endif
 
-  double getClusterElement(int ic, int ir, int cmsub=0){return cluster->get_data(ic,ir);};
-  
 
  private:
   
@@ -146,80 +227,22 @@ class singlePhotonDetector {
   int nx, ny;
 
 
-
-  MovingStat **stat;
+  pedestalSubtraction **stat;
+  commonModeSubtraction *cmSub;
+  //MovingStat **stat;
+  int nDark;
+  eventType **eventMask;
   double nSigma;
   int clusterSize, clusterSizeY;
   single_photon_hit *cluster;
+  int iframe;
+  int dataSign;
 
-  MovingStat cmStat[4]; //stores the common mode average per supercolumn
-  double cmPed[4]; // stores the common mode for this frame per supercolumn
-  double nCm[4]; // stores the number of pedestals to calculate the cm per supercolumn
-  int cx, cy;
+  /*   int cx, cy; */
 };
 
 
 
 
 
-
-
-    /////////////////////////////////////////////////////// DONE UNTIL HERE!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- /*  void calculateCommonMode(int xmin=0, int xmax=160, int ymin=0, int ymax=160, double hc=0, double tc=0) { */
-
-/*     int scmin=xmin/40; */
-/*     int scmax=(xmax-1)/40+1; */
-/*     int isc=0, ix, iy; */
-
-/*     for (isc=scmin; isc<scmax; isc++) { */
-/*       nCm[isc]=0; */
-/*       cmPed[isc]=0; */
-/*     } */
-    
-    
-
-
-/*     for (ix=xmin-1; ix<xmax+1; ix++) { */
-/*       isc=ix/40; */
-/*       for (iy=ymin-1; iy<ymax+1; iy++) { */
-/* 	if (isGood(ix,iy)) { */
-/* 	  if (getEventType(ix, iy, hc, tc, 1,0)==PEDESTAL) { */
-/* 	    cmPed[isc]+=getChannelShort(ix, iy, hc, tc); */
-/* 	    nCm[isc]++; */
-/* 	  } */
-/* 	} */
-	
-/*       } */
-
-/*     } */
-
-    
-/*     for (isc=scmin; isc<scmax; isc++) { */
-
-/*       if (nCm[isc]>0) */
-/* 	cmStat[isc].Calc(cmPed[isc]/nCm[isc]); */
-
-/*       //   cout << "SC " << isc << nCm[isc] << " " << cmPed[isc]/nCm[isc] << " " << cmStat[isc].Mean() << endl; */
-
-/*     } */
-    
-
-/*   }; */
-
-/*   double getCommonMode(int ix, int iy) { */
-/*     int isc=ix/40; */
-/*     if (nCm[isc]>0) { */
-/*      /\*  if (ix==20 && iy==20) *\/ */
-/* /\* 	cout << cmPed[isc] << " " << nCm[isc] << " " << cmStat[isc].Mean() << " " << cmPed[isc]/nCm[isc]-cmStat[isc].Mean() << endl; *\/ */
-/*       return cmPed[isc]/nCm[isc]-cmStat[isc].Mean(); */
-/*     } else { */
-/*   /\*     if (ix==20 && iy==20) *\/ */
-/* /\* 	cout << "No common mode!" << endl; *\/ */
-      
-/*       return 0; */
-/*     } */
-/*   }; */
-
-
-
 #endif
diff --git a/slsDetectorCalibration/slsDetectorData.h b/slsDetectorCalibration/slsDetectorData.h
index 527e4e9bb..59bd872f0 100644
--- a/slsDetectorCalibration/slsDetectorData.h
+++ b/slsDetectorCalibration/slsDetectorData.h
@@ -24,13 +24,13 @@ class slsDetectorData {
 
 
   */
-  slsDetectorData(int npx, int npy, int np, int psize, int **dMap=NULL, int **dMask=NULL, int **dROI=NULL): nx(npx), ny(npy), nPackets(np), packetSize(psize) {
+  slsDetectorData(int npx, int npy, int np, int psize, int **dMap=NULL, dataType **dMask=NULL, int **dROI=NULL): nx(npx), ny(npy), nPackets(np), packetSize(psize) {
 
    
 
-    dataMask=new int*[ny];
+    dataMask=new dataType*[ny];
     for(int i = 0; i < ny; i++) {
-      dataMask[i] = new int[nx];
+      dataMask[i] = new dataType[nx];
     }
     dataMap=new int*[ny];
     for(int i = 0; i < ny; i++) {
@@ -74,7 +74,7 @@ class slsDetectorData {
     
   };
 
-  void setDataMask(int **dMask=NULL){
+  void setDataMask(dataType **dMask=NULL){
 
     if (dMask!=NULL) {
 
@@ -102,9 +102,9 @@ class slsDetectorData {
     
   };
 
-  int setGood(int ix, int iy, int i=1) {dataROIMask[iy][ix]=i; return isGood(ix,iy);};
+  int setGood(int ix, int iy, int i=1) {  if (ix>=0 && ix<nx && iy>=0 && iy<ny) dataROIMask[iy][ix]=i; return isGood(ix,iy);};
 
-  int isGood(int ix, int iy) {return dataROIMask[iy][ix];};
+  int isGood(int ix, int iy) {  if (ix>=0 && ix<nx && iy>=0 && iy<ny) return dataROIMask[iy][ix]; else return -1;};
 
 
   void getDetectorSize(int &npx, int &npy){npx=nx; npy=ny;};
@@ -123,7 +123,11 @@ class slsDetectorData {
 
 
   virtual dataType getChannel(char *data, int ix, int iy=0) {
-    dataType m=dataMask[iy][ix], d=*((dataType*)(data+dataMap[iy][ix]));
+    dataType m=0, d=0;
+    if (ix>=0 && ix<nx && iy>=0 && iy<ny && dataMap[iy][ix]>=0 && dataMap[iy][ix]<nPackets*packetSize) {
+      m=dataMask[iy][ix];
+      d=*((dataType*)(data+dataMap[iy][ix]));
+    }  
     return d^m;
   };
 
@@ -234,7 +238,7 @@ class slsDetectorData {
   const int nPackets; /**<number of UDP packets constituting one frame */
   const int packetSize; /**< size of a udp packet */
   int **dataMap; /**< Array of size nx*ny storing the pointers to the data in the dataset (as offset)*/
-  int **dataMask; /**< Array of size nx*ny storing the polarity of the data in the dataset (should be 0 if no inversion is required, 0xffffffff is inversion is required) */
+  dataType **dataMask; /**< Array of size nx*ny storing the polarity of the data in the dataset (should be 0 if no inversion is required, 0xffffffff is inversion is required) */
   int **dataROIMask; /**< Array of size nx*ny 1 if channel is good (or in the ROI), 0 if bad channel (or out of ROI)   */
 
 };