From 02a69a0a6ca4495f9e527b185ad3d6e602fca29f Mon Sep 17 00:00:00 2001
From: Anna Bergamaschi <anna.bergamaschi@psi.ch>
Date: Tue, 27 Jan 2015 17:34:20 +0100
Subject: [PATCH] moench03Ctb data structure and readout funcs added

---
 slsDetectorCalibration/moench03CtbData.h  | 150 ++++++++++++
 slsDetectorCalibration/moench03ReadData.C | 267 ++++++++++++++++++++++
 slsDetectorCalibration/slsDetectorData.h  |  25 +-
 3 files changed, 433 insertions(+), 9 deletions(-)
 create mode 100644 slsDetectorCalibration/moench03CtbData.h
 create mode 100644 slsDetectorCalibration/moench03ReadData.C

diff --git a/slsDetectorCalibration/moench03CtbData.h b/slsDetectorCalibration/moench03CtbData.h
new file mode 100644
index 000000000..226e1f237
--- /dev/null
+++ b/slsDetectorCalibration/moench03CtbData.h
@@ -0,0 +1,150 @@
+#ifndef MOENCH03CTBDATA_H
+#define  MOENCH03CTBDATA_H
+#include "slsDetectorData.h"
+
+
+
+class moench03CtbData : public slsDetectorData<uint16_t> {
+
+ private:
+  
+  int iframe;
+  int *xmap, *ymap;
+  int nadc;
+  int sc_width;
+  int sc_height;
+ 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
+
+  */
+  
+
+  moench03CtbData(int ns=5000): slsDetectorData<uint16_t>(400, 400, ns*2*32, NULL, NULL) , nadc(32), sc_width(25), sc_height(200) {
+
+    
+    int adc_nr[32]={200,225,250,275,300,325,350,375,0,25,50,75,100,125,150,175,175,150,125,100,75,50,25,0,375,350,325,300,275,250,225,200};
+    int row, col;
+
+    int isample;
+    int iadc;
+    int ix, iy;
+
+    xmap=new int[nx*ny];
+    ymap=new int[nx*ny];
+    
+
+
+
+    for (iadc=0; iadc<nadc; iadc++) {
+      for (int i=0; i<sc_width*sc_height; i++) {
+	col=adc_nr[iadc]+(i%sc_width);
+	if (iadc<16) {
+	  row=199-i/sc_width;
+	} else {
+	  row=200+i/sc_width;
+	}
+	dataMap[row][col]=(nadc*i+iadc)*2;
+	if (dataMap[row][col]<0 || dataMap[row][col]>=2*400*400)
+	  cout << "Error: pointer " << dataMap[row][col] << " out of range "<< endl;
+	  
+      }
+    }
+    for (int i=0; i<nx*ny; i++) {
+      isample=i/nadc;
+      iadc=i%nadc;
+      ix=isample%sc_width;
+      iy=isample/sc_width;
+      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;
+      }
+
+
+    }
+
+
+
+    
+    
+    iframe=0;
+    //  cout << "data struct created" << endl;
+  };
+    
+    void getPixel(int ip, int &x, int &y) {if (ip>=0 && ip<nx*ny) {x=xmap[ip]; y=ymap[ip];}};
+  
+
+     /**
+
+     Returns the frame number for the given dataset. Purely virtual func.
+     \param buff pointer to the dataset
+     \returns frame number
+
+  */
+
+
+    virtual  int getFrameNumber(char *buff){(void)buff; return iframe;};   
+
+  /**
+
+     Returns the packet number for the given dataset. purely virtual func
+     \param buff pointer to the dataset
+     \returns packet number number
+
+
+  virtual int getPacketNumber(char *buff)=0;
+
+  */
+
+   /**
+
+     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; 
+    if (filebin.is_open()) {
+	afifo_cont=new uint16_t[dataSize/2];
+ 	if (filebin.read((char*)afifo_cont,dataSize)) {
+	  iframe++;
+	  return (char*)afifo_cont;
+	} else {
+	delete [] afifo_cont;
+	return NULL;
+	}
+    }     
+    return NULL;
+    };
+
+
+  
+
+};
+
+
+
+#endif
diff --git a/slsDetectorCalibration/moench03ReadData.C b/slsDetectorCalibration/moench03ReadData.C
new file mode 100644
index 000000000..886d9621f
--- /dev/null
+++ b/slsDetectorCalibration/moench03ReadData.C
@@ -0,0 +1,267 @@
+#include <TH1D.h>
+#include <TH2D.h>
+#include <TPad.h>
+#include <TDirectory.h>
+#include <TEntryList.h>
+#include <TFile.h>
+#include <TMath.h>
+#include <TTree.h>
+#include <TChain.h>
+#include <THStack.h>
+#include <TCanvas.h>
+#include <stdio.h>
+//#include <deque>
+//#include <list>
+//#include <queue>
+#include <fstream>
+#include "moench03CtbData.h"
+#include "moenchCommonMode.h"
+#define MYROOT1
+#include "singlePhotonDetector.h"
+
+//#include "MovingStat.h"
+
+using namespace std;
+
+#define NC 400
+#define NR 400
+
+
+#define MY_DEBUG 1
+#ifdef MY_DEBUG
+#include <TCanvas.h>
+#endif
+
+
+
+TH2F *readImage(ifstream &filebin, int dsize=5000*32) {
+  moench03CtbData *decoder=new moench03CtbData();
+  char *buff=decoder->readNextFrame(filebin);
+  TH2F *h2=NULL;
+
+//   TH1F *h1=new TH1F("h1","",400*400,0,400*400);
+//   int ip=0;
+  if (buff) {
+    h2=new TH2F("h2","",400,0,400,400,0,400);
+    h2->SetStats(kFALSE);
+    for (int ix=0; ix<400; ix++) {
+      for (int iy=0; iy<400; iy++) {
+	//	cout <<  decoder->getDataSize() << " " << decoder->getValue(buff,ix,iy)<< endl;
+	h2->SetBinContent(ix+1,iy+1,decoder->getValue(buff,ix,iy));
+	//	h1->SetBinContent(++ip,decoder->getValue(buff,ix,iy));
+      }
+    }
+  }
+  return h2;
+}
+
+
+TH2F *readImage(char *fname) {
+  ifstream filebin;
+  filebin.open((const char *)(fname), ios::in | ios::binary);
+  TH2F *h2=readImage(filebin);
+  filebin.close();
+  return h2;
+}
+/**
+
+Loops over data file to find single photons, fills the tree (and writes it to file, althoug the root file should be opened before) and creates 1x1, 2x2, 3x3 cluster histograms with ADCu on the x axis, channel number (160*x+y) on the y axis.
+
+  \param fformat file name format
+  \param tit title of the tree etc.
+  \param runmin minimum run number
+  \param runmax max run number
+  \param nbins  number of bins for spectrum hists
+  \param hmin histo minimum for spectrum hists
+  \param hmax histo maximum for spectrum hists
+  \param sign sign of the spectrum to find hits
+  \param hc readout correlation coefficient with previous pixel
+  \param xmin minimum x coordinate
+  \param xmax maximum x coordinate
+  \param ymin minimum y coordinate
+  \param ymax maximum y coordinate
+  \param cmsub  enable commonmode subtraction
+  \returns pointer to histo stack with cluster spectra
+*/
+
+THStack *moench03ReadData(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=1, int xmax=NC-1, int ymin=1, int ymax=NR-1, int cmsub=0, int hitfinder=1) {
+  
+  moench03CtbData *decoder=new moench03CtbData();
+   moenchCommonMode *cmSub=NULL;
+//   if (cmsub)
+//     cmSub=new moenchCommonMode();
+   int iev=0;
+  int nph=0;
+
+  singlePhotonDetector<uint16_t> *filter=new singlePhotonDetector<uint16_t>(decoder, 3, 5, sign, cmSub);
+
+  char *buff;
+  char fname[10000];
+  int nf=0;
+  
+  eventType thisEvent=PEDESTAL;
+
+  // int iframe;
+  // double *data, ped, sigma;
+
+  // data=decoder->getCluster();
+
+  TH2F *h2;
+  TH2F *h3;
+  TH2F *hetaX;
+  TH2F *hetaY;
+
+  THStack *hs=new THStack("hs",fformat);
+
+
+
+  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 (hitfinder) {
+    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;
+  TTree *tall;
+  if (hitfinder)
+    tall=filter->initEventTree(tit, &iFrame);
+
+
+  
+
+#ifdef MY_DEBUG
+  
+  TCanvas *myC;
+  TH2F *he;
+  TCanvas *cH1;
+  TCanvas *cH2;
+  TCanvas *cH3;
+
+  if (hitfinder) {
+    myC=new TCanvas();
+    he=new TH2F("he","Event Mask",xmax-xmin, xmin, xmax, ymax-ymin, ymin, ymax);
+    he->SetStats(kFALSE);
+    he->Draw("colz");
+    cH1=new TCanvas();
+    cH1->SetLogz();
+    h1->Draw("colz");
+    cH2=new TCanvas();
+    cH2->SetLogz();
+    h2->Draw("colz");
+    cH3=new TCanvas();
+    cH3->SetLogz();
+    h3->Draw("colz");
+  }
+#endif
+  filter->newDataSet();
+
+
+  for (int irun=runmin; irun<runmax; irun++) {
+    sprintf(fname,fformat,irun);
+    cout << "file name " << fname << endl;
+    filebin.open((const char *)(fname), ios::in | ios::binary);
+    nph=0;
+    while ((buff=decoder->readNextFrame(filebin))) {
+      
+
+      if (hitfinder) {
+	filter->newFrame();
+
+	//calculate pedestals and common modes
+	if (cmsub) {
+	//	cout << "cm" << endl;
+	  for (ix=xmin-1; ix<xmax+1; ix++)
+	    for (iy=ymin-1; iy<ymax+1; iy++) {
+	      thisEvent=filter->getEventType(buff, ix, iy,0);
+	    }
+	}
+      }
+
+      //   cout << "new frame " << endl;
+
+      for (ix=xmin-1; ix<xmax+1; ix++)
+	for (iy=ymin-1; iy<ymax+1; iy++) {
+	  //	  cout << ix << " " << iy << endl;
+	  thisEvent=filter->getEventType(buff, ix, iy, cmsub);
+
+#ifdef MY_DEBUG
+	  if (hitfinder) {
+	    if (iev%1000==0)
+	      he->SetBinContent(ix+1-xmin, iy+1-ymin, (int)thisEvent);
+	  }
+#endif	  
+	  
+	  if (nf>1000) {
+	    h1->Fill(filter->getClusterTotal(1), iy+NR*ix);
+	    if (hitfinder) {
+	      
+	      if (thisEvent==PHOTON_MAX ) {
+		nph++;
+
+		h2->Fill(filter->getClusterTotal(2), iy+NR*ix);
+		h3->Fill(filter->getClusterTotal(3), iy+NR*ix);
+		iFrame=decoder->getFrameNumber(buff);
+
+		tall->Fill();
+	   
+       
+	      }
+	  
+	    
+	    }
+	    
+
+	  }
+	}
+	  //////////////////////////////////////////////////////////
+    
+#ifdef MY_DEBUG
+      if (iev%1000==0) {
+	myC->Modified();
+	myC->Update();
+	cH1->Modified();
+	cH1->Update();
+	cH2->Modified();
+	cH2->Update();
+	cH3->Modified();
+	cH3->Update();
+      }
+      iev++;
+#endif		    
+      nf++;
+      
+      cout << "=" ;
+      delete [] buff;
+    }
+    cout << nph << endl;
+    if (filebin.is_open())
+    filebin.close();	  
+    else
+      cout << "could not open file " << fname << endl;
+  } 
+  if (hitfinder)
+    tall->Write(tall->GetName(),TObject::kOverwrite);
+  
+    
+  
+  delete decoder;
+  cout << "Read " << nf << " frames" << endl;
+  return hs;
+}
+  
diff --git a/slsDetectorCalibration/slsDetectorData.h b/slsDetectorCalibration/slsDetectorData.h
index 294f7c871..8240a7484 100644
--- a/slsDetectorCalibration/slsDetectorData.h
+++ b/slsDetectorCalibration/slsDetectorData.h
@@ -11,6 +11,14 @@ using namespace std;
 template <class dataType>
 class slsDetectorData {
 
+ protected:
+  const int nx; /**< Number of pixels in the x direction */
+  const int ny; /**< Number of pixels in the y direction */
+  int dataSize; /**<size of the data constituting one frame */
+  int **dataMap; /**< Array of size nx*ny storing the pointers to the data in the dataset (as offset)*/
+  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)   */
+
 
  public:
 
@@ -65,6 +73,10 @@ class slsDetectorData {
     delete [] dataROIMask;
   }
 
+  virtual void getPixel(int ip, int &x, int &y) {x=ip; y=0;};
+
+
+
   /**
      defines the data map (as offset) - no error checking if datasize and offsets are compatible!
      \param dMap array of size nx*ny storing the pointers to the data in the dataset (as offset). If NULL (default),the data are arranged as if read out row by row (dataMap[iy][ix]=(iy*nx+ix)*sizeof(dataType);)
@@ -80,9 +92,12 @@ class slsDetectorData {
 	 for (int ix=0; ix<nx; ix++)
 	   dataMap[iy][ix]=(iy*nx+ix)*sizeof(dataType);
     } else {
+      cout << "set dmap "<< dataMap << " " << dMap << endl;
       for (int iy=0; iy<ny; iy++)
-	 for (int ix=0; ix<nx; ix++)
+	for (int ix=0; ix<nx; ix++) {
 	   dataMap[iy][ix]=dMap[iy][ix];
+	   cout << ix << " " << iy << endl;
+	}
     }
     
   };
@@ -236,14 +251,6 @@ class slsDetectorData {
   */
   virtual char *readNextFrame(ifstream &filebin)=0;
 
- protected:
-  const int nx; /**< Number of pixels in the x direction */
-  const int ny; /**< Number of pixels in the y direction */
-  int dataSize; /**<size of the data constituting one frame */
-  int **dataMap; /**< Array of size nx*ny storing the pointers to the data in the dataset (as offset)*/
-  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)   */
-
 };