detectors/slsDetectorPackage
25
0
Fork 0
mirror of https://github.com/slsdetectorgroup/slsDetectorPackage.git synced 2025-06-23 10:07:59 +02:00
Code Issues Actions Packages Releases Activity

Merge branch 'master' of gitorious.psi.ch:sls_det_software/sls_detector_calibration

Browse Source
This commit is contained in:
Dhanya Maliakal
2015-06-26 11:53:01 +02:00
parent 9041c48bbf b820ae0cad
commit 1f9f474670
9 changed files with 2006 additions and 134 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
slsDetectorCalibration/energyCalibration.cpp
View File

@ -524,4 +524,4 @@ TGraphErrors* energyCalibration::calibrate(int nscan, Double_t *en, Double_t *ee
Fit data is empty

8
slsDetectorCalibration/makefile Normal file
View File

@ -0,0 +1,8 @@
CC=gcc
CFLAGS= -g -c -W -lstdc++
ROOTINCLUDE=$(ROOTSYS)/include
energyCalibration.o: energyCalibration.cpp
$(CC) $(CFLAGS) energyCalibration.cpp -I $(ROOTINCLUDE)

225
slsDetectorCalibration/moench03Ctb10GbData.h Normal file
View File

@ -0,0 +1,225 @@
#ifndef MOENCH03CTB10GBDATA_H
#define MOENCH03CTB10GBDATA_H
#include "slsReceiverData.h"
class moench03Ctb10GbData : public slsReceiverData<uint16_t> {
private:
int iframe;
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
*/
// moench03Ctb10GbData(int ns=5000): slsDetectorData<uint16_t>(400, 400, 8208*40, NULL, NULL) , nadc(32), sc_width(25), sc_height(200) {
moench03Ctb10GbData(int ns=5000): slsReceiverData<uint16_t>(400, 400, 40, 8208), 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;
int npackets=40;
int i;
for (int ip=0; ip<npackets; ip++) {
for (int is=0; is<128; is++) {
for (iadc=0; iadc<nadc; iadc++) {
i=128*ip+is;
if (i<sc_width*sc_height) {
// 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+16*(ip+1);
if (dataMap[row][col]<0 || dataMap[row][col]>=8208*40)
cout << "Error: pointer " << dataMap[row][col] << " out of range "<< endl;
}
}
}
}
int ipacket;
int ibyte;
int ii=0;
for (int ipacket=0; ipacket<npackets; ipacket++) {
for (int ibyte=0; ibyte< 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;
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;
}
ii++;
}
}
}
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 */
/* *\/ */
/* 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& 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;
}
};
virtual char *readNextFrame(ifstream &filebin) {
int fnum;
return readNextFrame(filebin, fnum);
};
};
#endif

854
slsDetectorCalibration/moench03ReadData.C
View File

@ -9,11 +9,14 @@
#include <TChain.h>
#include <THStack.h>
#include <TCanvas.h>
#include <TF1.h>
#include <stdio.h>
//#include <deque>
//#include <list>
//#include <queue>
#include <map>
#include <fstream>
#include <sys/stat.h>
#include "moench03CtbData.h"
#include "moench03CommonMode.h"
#define MYROOT1
@ -27,7 +30,7 @@ using namespace std;
#define NR 400
//#define MY_DEBUG 1
#define MY_DEBUG 1
#ifdef MY_DEBUG
#include <TCanvas.h>
@ -47,7 +50,7 @@ TH2F *readImage(ifstream &filebin, TH2F *h2=NULL, TH2F *hped=NULL) {
h2=new TH2F("h2","",400,0,400,400,0,400);
h2->SetStats(kFALSE);
}
// cout << "." << endl;
cout << "." << endl;
for (int ix=0; ix<400; ix++) {
for (int iy=0; iy<400; iy++) {
// cout << decoder->getDataSize() << " " << decoder->getValue(buff,ix,iy)<< endl;
@ -77,7 +80,7 @@ TH2F *readImage(char *fname, int iframe=0, TH2F *hped=NULL) {
if (hh==NULL) break;
hh=readImage(filebin, h2, hped );
if (hh)
;// cout << "="<< endl;
cout << "="<< endl;
else {
delete h2;
return NULL;
@ -114,7 +117,7 @@ TH2F *calcPedestal(char *fformat, int runmin, int runmax){
}
}
delete [] buff;
//cout << "="<< endl;
cout << "="<< endl;
ii++;
}
if (filebin.is_open())
@ -134,7 +137,400 @@ TH2F *calcPedestal(char *fformat, int runmin, int runmax){
return h2;
}
/******************************************************************************/
TH2F * calcNoiseRMS(char *fformat, int runmin, int runmax){
ifstream filebin;
char fname[10000];
moench03CtbData *decoder=new moench03CtbData();
singlePhotonDetector<uint16_t> *filter=new singlePhotonDetector<uint16_t>(decoder, 3, 5, 1, NULL);
char *buff;
int ix,iy;
int ii=0;
TH2F* h2=NULL;
for (int irun=runmin; irun<=runmax; irun++) {
sprintf(fname,fformat,irun);
filebin.open((const char *)(fname), ios::in | ios::binary);
while ((buff=decoder->readNextFrame(filebin))) {
for (ix=0; ix<400; ix++) {
for (iy=0; iy<400; iy++) {
filter->addToPedestal(decoder->getValue(buff,ix,iy), ix, iy);
}
}
delete [] buff;
cout << "="<< endl;
ii++;
}
if (filebin.is_open())
filebin.close();
}
if (ii>0) {
h2=new TH2F("hped","",400,0,400,400,0,400);
for (ix=0; ix<400; ix++) {
for (iy=0; iy<400; iy++) {
h2->SetBinContent(ix+1, iy+1,filter->getPedestalRMS(ix,iy));
}
}
}
return h2;
}
/**********************************************************************************/
THStack * calcNoise(char *fformat,std::string flag, int runmin, int runmax, int nfiles){
std::map<int,TH1F *> pixelMap;
cerr<<"Creating pedestal map for frames: "<<flag<<endl;
cerr<<"/***************************************/"<<endl;
THStack * hs = new THStack();
int n(2);
int r(0);
if(flag=="even"){
n=2;
r=0;
}else if(flag=="odd"){
n=2;
r=1;
}else if(flag=="all"){
n=1;
r=0;
}
ifstream filebin;
char fname[10000];
moench03CtbData *decoder=new moench03CtbData();
singlePhotonDetector<uint16_t> *filter=new singlePhotonDetector<uint16_t>(decoder, 3, 5, 1, NULL);
char *buff;
int ix,iy;
int ii=0;
TH2F* hped=NULL;
TH2F* hnoise=NULL;
int ibin(0);
int iframe(0);
Float_t sigma(0.);
int filefound(0);
pixelMap.clear();
for (int irun=runmin; irun<=runmax; irun++) {
sprintf(fname,fformat,irun);
if( access(fname,F_OK) != -1){
if(filefound<nfiles){
filebin.open((const char *)(fname), ios::in | ios::binary);
while ((buff=decoder->readNextFrame(filebin))) {
iframe=decoder->getFrameNumber(buff);
if(iframe%n==r){
for (ix=0; ix<400; ix++) {
for (iy=0; iy<400; iy++) {
ibin = ix*NR+iy;
if(pixelMap.find(ibin)!=pixelMap.end()){
pixelMap[ibin]->Fill(decoder->getValue(buff,ix,iy));
}else{
pixelMap[ibin] = new TH1F(Form("h_%d",ibin),Form("h_%d",ibin),8000,-0.5,15999.5);
}
}
}
}
delete [] buff;
if(ii%1000==0)cout << "="<<ii<< endl;
ii++;
}
if (filebin.is_open())
filebin.close();
filefound++;
}else{
break;
}
}else{//end chek if file exists
cerr<<fname<<" does not exist!"<<endl;
}
}
if(ii>0){
hped=new TH2F("hped","",400,-0.5,399.5,400,-0.5,399.5);
hnoise=new TH2F("hnoise","",400,-0.5,399.5,400,-0.5,399.5);
hs->Add(hped);
hs->Add(hnoise);
Float_t p(0.);
Float_t perr(0.);
Float_t amp(0.);
Float_t sigmaerr(0.);
for(ix=0; ix<400; ix++){
for (iy=0; iy<400; iy++){
sigma=0.;
sigmaerr=0.;
p=0.;
perr=0.;
ibin = ix*NR+iy;
if(ibin%1000==0)cerr<<ibin<<endl;
if(pixelMap.find(ibin)!=pixelMap.end()){
// cerr<<"/****************/"<<endl;
// cerr<<"("<<ix<<","<<iy<<")"<<endl;
p=pixelMap[ibin]->GetBinCenter(pixelMap[ibin]->GetMaximumBin());
// cerr<<"peak: "<<p<<endl;
amp=pixelMap[ibin]->GetBinContent(pixelMap[ibin]->GetMaximumBin());
sigma=pixelMap[ibin]->GetRMS();
TF1 * g = new TF1("g","gaus(0)",p-sigma,p+sigma);
g->SetParameters(amp,p,sigma);
g->SetParLimits(1,p-sigma,p+sigma);
g->SetParLimits(2,0.1*sigma,2.*sigma);
pixelMap[ibin]->Fit(g,"RQN0");
p=g->GetParameter(1);
perr=g->GetParError(1);
sigma=g->GetParameter(2);
sigmaerr=g->GetParError(2);
// cerr<<"sigma: "<<sigma<<endl;
delete g;
delete pixelMap[ibin];
}
hped->SetBinContent(ix+1,iy+1,p);
hped->SetBinError(ix+1,iy+1,perr);
hnoise->SetBinContent(ix+1,iy+1,sigma);
hnoise->SetBinError(ix+1,iy+1,sigmaerr);
}
}
}
return hs;
}
/****************************************************/
THStack * DrawFrames(char *fformat,int framemin, int framemax){
int nbins(400);
int xmin(0);
int xmax(400);
THStack * hs = new THStack();
ifstream filebin;
char fname[10000];
moench03CtbData *decoder=new moench03CtbData();
char *buff;
int ix,iy;
int iframe(0);
if(access(fformat,F_OK) != -1){
filebin.open((const char *)(fformat), ios::in | ios::binary);
while ((buff=decoder->readNextFrame(filebin))) {
iframe=decoder->getFrameNumber(buff);
if(iframe>=framemin && iframe<framemax){//check on frames
ix=0;
iy=0;
TH2F * h = new TH2F(Form("h_%d",iframe),Form("h_%d",iframe),nbins,xmin,xmax,nbins,xmin,xmax);
for(ix=xmin; ix<xmax; ix++){//loop on pixels
for(iy=xmin; iy<xmax; iy++){
h->SetBinContent(ix+1,iy+1,decoder->getValue(buff,ix,iy));
}
}//end loop on pixels
hs->Add(h);
}//end check on frames
delete [] buff;
}
}else{//check if file exists
cerr<<fformat<<" does not exists!"<<endl;
}
return hs;
}
/****************************************************/
TH1F * SinglePixelHisto(char *fformat,int nbins, float xmin, float xmax, int runmin, int runmax, int x0, int y0){
moench03CtbData *decoder=new moench03CtbData();
cout << "decoder allocated " << endl;
TH1F * _hout = new TH1F("hpix","",nbins,xmin,xmax);
char *buff;
char fname[10000];
ifstream filebin;
int iframe(0);
for (int irun=runmin; irun<=runmax; irun++){//loop on files
sprintf(fname,fformat,irun);
iframe=0;
if( access(fname,F_OK) != -1){//check if file exists
filebin.open((const char *)(fname), ios::in | ios::binary);
while ((buff=decoder->readNextFrame(filebin))){
iframe=decoder->getFrameNumber(buff);
if(iframe%1000==0) cerr<<iframe<<" frames"<<endl;
_hout->Fill(decoder->getValue(buff,x0,y0));
delete [] buff;
}//end readout of frame
}else{//end check on file
cerr<<fname<<" does not exists!"<<endl;
}
if(filebin.is_open()) filebin.close();
}//end loop on files
return _hout;
}
/****************************************************/
void DrawAllPixelCorrelation(char *fformat, std::string frame,int runmin, int runmax, int xmin, int xmax, int ymin, int ymax, int sc_width,TH2F* hcorr){
int cmsub(0);
moench03CtbData *decoder=new moench03CtbData();
cout << "decoder allocated " << endl;
moench03CommonMode *cmSub=NULL;
if (cmsub) {
cmSub=new moench03CommonMode(100);
cout << "common mode allocated " << endl;
} else {
cout << "non allocating common mode " << endl;
}
int n(2);
int r(0);
if(frame=="even"){
n=2;
r=0;
}else if(frame=="odd"){
n=2;
r=1;
}else if(frame=="all"){
n=1;
r=0;
}
int sign(1);
singlePhotonDetector<uint16_t> *filter=new singlePhotonDetector<uint16_t>(decoder, 3, 5, sign, cmSub, 1000, 100);
char *buff;
char fname[10000];
float pix1(0.);
float pix2(0.);
ifstream filebin;
int iframe(0);
for (int irun=runmin; irun<=runmax; irun++){//loop on files
sprintf(fname,fformat,irun);
iframe=0;
pix1=0.;
pix2=0.;
if( access(fname,F_OK) != -1){//check if file exists
filebin.open((const char *)(fname), ios::in | ios::binary);
while ((buff=decoder->readNextFrame(filebin))){//loop on frames
iframe=decoder->getFrameNumber(buff);
if(iframe%n==r){//check if frame is even or odd
for(int ix=xmin; ix<xmax; ix++){//loop on pixels
for(int iy=ymin; iy<ymax; iy++){
if((ix+1)%sc_width != 0){//check that pixel is not at the edge of supercolumn
pix1=decoder->getValue(buff,ix,iy);
pix2=decoder->getValue(buff,ix+1,iy);
hcorr->Fill(pix1,pix2);
}//end check on s_c edge
}
}//end loop on pixels
}//edn check on frame number
delete [] buff;
}//end loop on frames
}else{//end check if file exists
cerr<<fname<<" does not exist!"<<endl;
}
if(filebin.is_open())filebin.close();
}//end loop on files
return;
}
/***************************************************/
void DrawPixelCorrelation(char *fformat, std::string frame,int runmin, int runmax, int px, int py,TH1F * h1,TH1F * h2, TH2F * hcorr){
int cmsub(0);
moench03CtbData *decoder=new moench03CtbData();
cout << "decoder allocated " << endl;
moench03CommonMode *cmSub=NULL;
if (cmsub) {
cmSub=new moench03CommonMode(100);
cout << "common mode allocated " << endl;
} else {
cout << "non allocating common mode " << endl;
}
int n(2);
int r(0);
if(frame=="even"){
n=2;
r=0;
}else if(frame=="odd"){
n=2;
r=1;
}else if(frame=="all"){
n=1;
r=0;
}
int sign(1);
singlePhotonDetector<uint16_t> *filter=new singlePhotonDetector<uint16_t>(decoder, 3, 5, sign, cmSub, 1000, 100);
char *buff;
char fname[10000];
float pix1(0.);
float pix2(0.);
ifstream filebin;
int ii=0;
int iframe(0);
for (int irun=runmin; irun<=runmax; irun++){//loop on files
sprintf(fname,fformat,irun);
pix1=0.;
pix2=0.;
iframe=0;
if( access(fname,F_OK) != -1){//check if file exists
filebin.open((const char *)(fname), ios::in | ios::binary);
while ((buff=decoder->readNextFrame(filebin))){//loop on frames
iframe=decoder->getFrameNumber(buff);
if(iframe%n==r){
pix1=decoder->getValue(buff,px,py);
pix2=decoder->getValue(buff,px+1,py);
h1->Fill(pix1);
h2->Fill(pix2);
hcorr->Fill(pix1,pix2);
delete [] buff;
}
if(ii%1000==0)cerr<<"reading frame "<<ii<<endl;
ii++;
}//end loop on frames
if(filebin.is_open())filebin.close();
}else{//end check on file
cerr<<fname<<" does not exists!"<<endl;
}
if(filebin.is_open())filebin.close();
}//end loop on files
return;
}
/**
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.
@ -146,6 +542,8 @@ Loops over data file to find single photons, fills the tree (and writes it to fi
\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
@ -155,9 +553,9 @@ Loops over data file to find single photons, fills the tree (and writes it to fi
*/
THStack *moench03ReadData(char *fformat, char *tit, int runmin, int runmax, int nbins=1500, int hmin=-500, int hmax=1000, int xmin=1, int xmax=NC-1, int ymin=1, int ymax=NR-1, int cmsub=0, int hitfinder=1) {
double hc=0;
int sign=1;
double hc(0);
int sign(1);
moench03CtbData *decoder=new moench03CtbData();
cout << "decoder allocated " << endl;
@ -173,7 +571,7 @@ THStack *moench03ReadData(char *fformat, char *tit, int runmin, int runmax, int
int iev=0;
int nph=0;
singlePhotonDetector<uint16_t> *filter=new singlePhotonDetector<uint16_t>(decoder, 3, 5, sign, cmSub, 100, 10);
singlePhotonDetector<uint16_t> *filter=new singlePhotonDetector<uint16_t>(decoder, 3, 5, sign, cmSub, 1000, 100);
cout << "filter allocated " << endl;
char *buff;
@ -197,8 +595,230 @@ THStack *moench03ReadData(char *fformat, char *tit, int runmin, int runmax, int
hs->Add(h1);
cout << "h1 allocated " << endl;
TH2F *h2;
TH2F *h3;
// TH2F *h2=NULL;
// TH2F *h3=NULL;
if (hitfinder) {
// h2=new TH2F("h2",tit,nbins,hmin-0.5,hmax-0.5,NC*NR,-0.5,NC*NR-0.5);
// cout << "h2 allocated " << endl;
// h3=new TH2F("h3",tit,nbins,hmin-0.5,hmax-0.5,NC*NR,-0.5,NC*NR-0.5);
// cout << "h3 allocated " << endl;
// // 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);
}
if (hs->GetHists()) {
for (int i=0; i<1; i++)
if (hs->GetHists()->At(i)) cout << i << " " ;
cout << " histos allocated " << endl;
} else
cout << "no hists in stack " << endl;
int ix=20, iy=20, ir, ic;
int adc(0);
ifstream filebin;
Int_t iFrame;
TTree *tall;
if (hitfinder){
tall=filter->initEventTree(tit, &iFrame);
}
#ifdef MY_DEBUG
cout << "debug mode " << endl;
TCanvas *myC;
TH2F *he;
TCanvas *cH1;
// TCanvas *cH2;
// TCanvas *cH3;
cH1=new TCanvas("ch1");
cH1->SetLogz();
h1->Draw("colz");
if (hitfinder) {
myC=new TCanvas("myc");
he=new TH2F("he","Event Mask",xmax-xmin, xmin, xmax, ymax-ymin, ymin, ymax);
he->SetStats(kFALSE);
he->Draw("colz");
// cH2=new TCanvas("ch2");
// cH2->SetLogz();
// h2->Draw("colz");
// cH3=new TCanvas("ch3");
// cH3->SetLogz();
// h3->Draw("colz");
}
#endif
filter->newDataSet();
for (int irun=runmin; irun<runmax; irun++) {
sprintf(fname,fformat,irun);
cout << "file name " << fname << endl;
if( access(fname,F_OK) != -1){//check if file exists
filebin.open((const char *)(fname), ios::in | ios::binary);
nph=0;
while ( (buff=decoder->readNextFrame(filebin))) {//getting frames on the buffer
if (hitfinder) {
filter->newFrame();
//calculate pedestals and common modes
if (cmsub) {
// cout << "cm" << endl;
for (ix=xmin; ix<xmax+1; ix++)
for (iy=ymin; iy<ymax+1; iy++) {
thisEvent=filter->getEventType(buff, ix, iy,0);
}//loop on pixels
}//end if cmsub
}//end if hitfinder
// cout << "new frame " << endl;
for (ix=xmin; ix<xmax+1; ix++)
for (iy=ymin; iy<ymax+1; iy++) {
// cout << ix << " " << iy << endl;
thisEvent=filter->getEventType(buff, ix, iy, cmsub);
iFrame=decoder->getFrameNumber(buff);
#ifdef MY_DEBUG
if (hitfinder) {
// if (iev%10==0)
he->SetBinContent(ix+1-xmin, iy+1-ymin, (int)thisEvent);
}
#endif
// if (nf>10) {
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);
tall->Fill();
}//end if PHOTON_MAX
}//end if hitfinder
}//end loop on pixels
//////////////////////////////////////////////////////////
#ifdef MY_DEBUG
// cout << iev << " " << h1->GetEntries() << " " << h2->GetEntries() << endl;
// if (iev%10==0) {
// myC->Modified();
// myC->Update();
// cH1->Modified();
// cH1->Update();
// cH2->Modified();
// cH2->Update();
// cH3->Modified();
// cH3->Update();
// }
iev++;
#endif
nf++;
cout << "=" ;
delete [] buff;
}//end frame
cout << nph <<" Photons found"<< endl;
if (filebin.is_open())
filebin.close();
else
cout << "could not open file " << fname << endl;
}else{
cerr<<"File does not exist"<<endl;
}//end check on file path
}//end loop on files
if (hitfinder){
tall->Write(tall->GetName(),TObject::kOverwrite);
}
//////////////////////////////////////////////////////////
#ifdef MY_DEBUG
if(hitfinder){
myC->Modified();
myC->Update();
}
cH1->Modified();
cH1->Update();
// cH2->Modified();
// cH2->Update();
// cH3->Modified();
// cH3->Update();
#endif
delete decoder;
cout << "Read " << nf << " frames" << endl;
return hs;
}
/*********************************************************************************************/
THStack *moench03ReadDataEvenOdd(char *fformat, char *tit, int runmin, int runmax, int nbins=1500, int hmin=-500, int hmax=1000, int xmin=1, int xmax=NC-1, int ymin=1, int ymax=NR-1, int cmsub=0, int hitfinder=1){
double hc(0);
int sign(1);
moench03CtbData *decoder=new moench03CtbData();
cout << "decoder allocated " << endl;
moench03CommonMode *cmSub=NULL;
if (cmsub) {
cmSub=new moench03CommonMode(100);
cout << "common mode allocated " << endl;
} else {
cout << "non allocating common mode " << endl;
}
int iev=0;
int nph=0;
singlePhotonDetector<uint16_t> *filterEven=new singlePhotonDetector<uint16_t>(decoder, 3, 5, sign, cmSub, 1000, 100);
cout << "filter for even frames allocated " << endl;
singlePhotonDetector<uint16_t> *filterOdd=new singlePhotonDetector<uint16_t>(decoder, 3, 5, sign, cmSub, 1000, 100);
cout << "filter for odd frames allocated " << endl;
char *buff;
char fname[10000];
int nf=0;
eventType thisEvent=PEDESTAL;
// int iframe;
// double *data, ped, sigma;
// data=decoder->getCluster();
THStack *hs=new THStack("hs",fformat);
cout << "hstack allocated " << endl;
TH2F *h1=new TH2F("h1",tit,nbins,hmin-0.5,hmax-0.5,NC*NR,-0.5,NC*NR-0.5);
hs->Add(h1);
cout << "h1 allocated " << endl;
TH2F *h2=NULL;
TH2F *h3=NULL;
if (hitfinder) {
h2=new TH2F("h2",tit,nbins,hmin-0.5,hmax-0.5,NC*NR,-0.5,NC*NR-0.5);
cout << "h2 allocated " << endl;
@ -211,9 +831,12 @@ THStack *moench03ReadData(char *fformat, char *tit, int runmin, int runmax, int
// hs->Add(hetaX);
// hs->Add(hetaY);
}
if (hs->GetHists()) {
for (int i=0; i<3; i++)
if (hs->GetHists()->At(1)) cout << i << " " ;
if (hs->GetHists()->At(i)) cout << i << " " ;
cout << " histos allocated " << endl;
} else
cout << "no hists in stack " << endl;
@ -225,12 +848,16 @@ THStack *moench03ReadData(char *fformat, char *tit, int runmin, int runmax, int
int ix=20, iy=20, ir, ic;
Int_t iFrame;
TTree *tall;
if (hitfinder)
tall=filter->initEventTree(tit, &iFrame);
TTree *tEven;
TTree *tOdd;
if (hitfinder){
tEven=filterEven->initEventTree(Form("%s_Even",tit), &iFrame);
cout<< "TTree for even frames initialized"<<endl;
tOdd=filterOdd->initEventTree(Form("%s_Odd",tit), &iFrame);
cout<< "TTree for odd frames initialized"<<endl;
}
#ifdef MY_DEBUG
@ -257,10 +884,15 @@ THStack *moench03ReadData(char *fformat, char *tit, int runmin, int runmax, int
cH3=new TCanvas("ch3");
cH3->SetLogz();
h3->Draw("colz");
cout<<"Canvas allocated"<<endl;
}
#endif
filter->newDataSet();
filterEven->newDataSet();
filterOdd->newDataSet();
for (int irun=runmin; irun<runmax; irun++) {
@ -269,54 +901,70 @@ THStack *moench03ReadData(char *fformat, char *tit, int runmin, int runmax, int
filebin.open((const char *)(fname), ios::in | ios::binary);
nph=0;
while ((buff=decoder->readNextFrame(filebin))) {
iFrame=decoder->getFrameNumber(buff);
filter->newFrame();
if (hitfinder) {
if(iFrame%2==0){
filterEven->newFrame();
}else{
filterOdd->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);
if(iFrame%2==0){
thisEvent=filterEven->getEventType(buff, ix, iy,0);
}else{
thisEvent=filterOdd->getEventType(buff, ix, iy,0);
}
}
}
}
// if (hitfinder) {
// //calculate pedestals and common modes
// }
// 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);
// if (nf>10) {
h1->Fill(filter->getClusterTotal(1), iy+NR*ix);
if(iFrame%2==0){
thisEvent=filterEven->getEventType(buff, ix, iy, cmsub);
}else{
thisEvent=filterOdd->getEventType(buff, ix, iy,cmsub);
}
#ifdef MY_DEBUG
// if (iev%10==0)
if (hitfinder) {
he->SetBinContent(ix+1-xmin, iy+1-ymin, (int)thisEvent);
}
#endif
if(iFrame%2==0){
h1->Fill(filterEven->getClusterTotal(1), iy+NR*ix);
}else{
h1->Fill(filterOdd->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();
if(iFrame%2==0){
h2->Fill(filterEven->getClusterTotal(2), iy+NR*ix);
h3->Fill(filterEven->getClusterTotal(3), iy+NR*ix);
tEven->Fill();
}else{
h2->Fill(filterOdd->getClusterTotal(2), iy+NR*ix);
h3->Fill(filterOdd->getClusterTotal(3), iy+NR*ix);
tOdd->Fill();
}
} // else {
// filter->addToPedestal(decoder->getValue(buff,ix,iy, cmsub));
// }
}
}
// }
@ -325,32 +973,34 @@ THStack *moench03ReadData(char *fformat, char *tit, int runmin, int runmax, int
#ifdef MY_DEBUG
// cout << iev << " " << h1->GetEntries() << " " << h2->GetEntries() << endl;
// if (iev%10==0) {
// myC->Modified();
// myC->Update();
// cH1->Modified();
// cH1->Update();
// cH2->Modified();
// cH2->Update();
// cH3->Modified();
// cH3->Update();
// }
// if (iev%10==0) {
// myC->Modified();
// myC->Update();
// cH1->Modified();
// cH1->Update();
// cH2->Modified();
// cH2->Update();
// cH3->Modified();
// cH3->Update();
// }
iev++;
#endif
nf++;
// cout << "=" ;
cout << "=" ;
delete [] buff;
}
// cout << nph << endl;
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);
if (hitfinder){
tEven->Write(tEven->GetName(),TObject::kOverwrite);
tOdd->Write(tOdd->GetName(),TObject::kOverwrite);
}
//////////////////////////////////////////////////////////
#ifdef MY_DEBUG
@ -369,3 +1019,99 @@ THStack *moench03ReadData(char *fformat, char *tit, int runmin, int runmax, int
return hs;
}
/***********************************************************************************/
Int_t moench03DrawPedestals(char *fformat, char *tit, int runmin, int runmax, int nbins=1500, int hmin=-500, int hmax=1000, int xmin=1, int xmax=NC-1, int ymin=1, int ymax=NR-1, int cmsub=0, int hitfinder=1){
double hc(0);
int sign(1);
moench03CtbData *decoder=new moench03CtbData();
cout << "decoder allocated " << endl;
moench03CommonMode *cmSub=NULL;
if (cmsub) {
cmSub=new moench03CommonMode(100);
cout << "common mode allocated " << endl;
} else {
cout << "non allocating common mode " << endl;
}
int iev=0;
int nph=0;
singlePhotonDetector<uint16_t> *filter=new singlePhotonDetector<uint16_t>(decoder, 3, 5, sign, cmSub, 1000, 100);
cout << "filter allocated " << endl;
char *buff;
char fname[10000];
int nf=0;
eventType thisEvent=PEDESTAL;
int iFrame(0);
// int iframe;
// double *data, ped, sigma;
// data=decoder->getCluster();
int nChx(xmax-xmin+2);
int nChy(ymax-ymin+2);
TH1F * HistoMap_All[nChx][nChy];
TH1F * HistoMap_Even[nChx][nChy];
TH1F * HistoMap_Odd[nChx][nChy];
for(Int_t ix=0; ix<nChx; ix++){
for(Int_t iy=0; iy<nChy; iy++){
HistoMap_All[ix][iy] = new TH1F(Form("h_%d_%d_All",ix+xmin,iy+ymin),"",8000,0,16000);
HistoMap_Even[ix][iy] = new TH1F(Form("h_%d_%d_Even",ix+xmin,iy+ymin),"",8000,0,16000);
HistoMap_Odd[ix][iy] = new TH1F(Form("h_%d_%d_Odd",ix+xmin,iy+ymin),"",8000,0,16000);
}
}
filter->newDataSet();
float adc(0.);
ifstream filebin;
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))) {
iFrame=decoder->getFrameNumber(buff);
filter->newFrame();
for(Int_t ix=0; ix<nChx; ix++){//loop on pixels
for(Int_t iy=0; iy<nChy; iy++){
thisEvent=filter->getEventType(buff, ix+xmin, iy+ymin,cmsub);
adc=decoder->getValue(buff,ix+xmin,iy+ymin);
HistoMap_All[ix][iy]->Fill(adc);
if(iFrame%2==0){
HistoMap_Even[ix][iy]->Fill(adc);
}else{
HistoMap_Odd[ix][iy]->Fill(adc);
}
}
}//end loop on pixels
}
if (filebin.is_open())
filebin.close();
else
cout << "could not open file " << fname << endl;
}//end loop on files
for(Int_t ix=0; ix<nChx; ix++){
for(Int_t iy=0; iy<nChy; iy++){
HistoMap_All[ix][iy]->Write();
HistoMap_Even[ix][iy]->Write();
HistoMap_Odd[ix][iy]->Write();
}
}
delete decoder;
return 1;
}

490
slsDetectorCalibration/moench03ReadData10Gb.C Normal file
View File

@ -0,0 +1,490 @@
#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 <TGraph.h>
#include <stdio.h>
//#include <deque>
//#include <list>
//#include <queue>
#include <fstream>
#include "moench03Ctb10GbData.h"
#include "moench03CommonMode.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, TH2F *h2=NULL, TH2F *hped=NULL) {
moench03Ctb10GbData *decoder=new moench03Ctb10GbData();
char *buff=decoder->readNextFrame(filebin);
// TH1F *h1=new TH1F("h1","",400*400,0,400*400);
// int ip=0;
if (buff) {
if (h2==NULL) {
h2=new TH2F("h2","",400,0,400,400,0,400);
h2->SetStats(kFALSE);
}
// cout << "." << endl;
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));
}
}
if (hped) h2->Add(hped,-1);
return h2;
}
return NULL;
}
TH2F *readImage(char *fname, int iframe=0, TH2F *hped=NULL) {
ifstream filebin;
filebin.open((const char *)(fname), ios::in | ios::binary);
TH2F *h2=new TH2F("h2","",400,0,400,400,0,400);
TH2F *hh;
hh=readImage(filebin, h2, hped );
if (hh==NULL) {
delete h2;
return NULL;
}
for (int i=0; i<iframe; i++) {
if (hh==NULL) break;
hh=readImage(filebin, h2, hped );
if (hh)
;// cout << "="<< endl;
else {
delete h2;
return NULL;
}
}
if (filebin.is_open())
filebin.close();
if (hped!=NULL)
h2->Add(hped,-1);
return h2;
}
TH2F *calcPedestal(char *fformat, int runmin, int runmax){
ifstream filebin;
char fname[10000];
moench03Ctb10GbData *decoder=new moench03Ctb10GbData();
singlePhotonDetector<uint16_t> *filter=new singlePhotonDetector<uint16_t>(decoder, 3, 5, 1, NULL);
char *buff;
int ix,iy;
int ii=0;
TH2F* h2=NULL;
for (int irun=runmin; irun<=runmax; irun++) {
sprintf(fname,fformat,irun);
cout << fname << endl;
filebin.open((const char *)(fname), ios::in | ios::binary);
while ((buff=decoder->readNextFrame(filebin))) {
for (ix=0; ix<400; ix++) {
for (iy=0; iy<400; iy++) {
if (decoder->getValue(buff,ix,iy)>1000)
filter->addToPedestal(decoder->getValue(buff,ix,iy), ix, iy);
}
}
delete [] buff;
//cout << "="<< endl;
ii++;
}
if (filebin.is_open())
filebin.close();
}
if (ii>0) {
h2=new TH2F("hped","",400,0,400,400,0,400);
for (ix=0; ix<400; ix++) {
for (iy=0; iy<400; iy++) {
h2->SetBinContent(ix+1, iy+1,filter->getPedestal(ix,iy));
}
}
}
return h2;
}
TH1D *calcSpectrum(char *fformat, int runmin, int runmax, TH2F *hped=NULL){
ifstream filebin;
char fname[10000];
moench03Ctb10GbData *decoder=new moench03Ctb10GbData();
TH1D *hspectrum=new TH1D("hsp","hsp",2500,-500,10000);
char *buff;
int ix,iy;
int ii=0;
TH2F* h2=NULL;
int ich=0;
Double_t ped=0;
for (int irun=runmin; irun<=runmax; irun++) {
sprintf(fname,fformat,irun);
cout << fname << endl;
filebin.open((const char *)(fname), ios::in | ios::binary);
while ((buff=decoder->readNextFrame(filebin))) {
for (ix=0; ix<200; ix++) {
for (iy=200; iy<400; iy++) {
if (decoder->getValue(buff,ix,iy)>1000) {
if(hped) ped=hped->GetBinContent(ix+1,iy+1);
hspectrum->Fill(decoder->getValue(buff,ix,iy)-ped);
}
ich++;
}
}
delete [] buff;
//cout << "="<< endl;
ii++;
}
if (filebin.is_open())
filebin.close();
}
return hspectrum;
}
TH2F *drawImage(char *fformat, int runmin, int runmax, TH2F *hped=NULL){
ifstream filebin;
char fname[10000];
moench03Ctb10GbData *decoder=new moench03Ctb10GbData();
TH2F *hspectrum=new TH2F("hsp","hsp",400,0,400,400,0,400);
char *buff;
int ix,iy;
int ii=0;
TH2F* h2=NULL;
int ich=0;
Double_t ped=0;
for (int irun=runmin; irun<=runmax; irun++) {
sprintf(fname,fformat,irun);
cout << fname << endl;
filebin.open((const char *)(fname), ios::in | ios::binary);
while ((buff=decoder->readNextFrame(filebin))) {
for (ix=0; ix<400; ix++) {
for (iy=0; iy<400; iy++) {
if (decoder->getValue(buff,ix,iy)>1000) {
if(hped) ped=hped->GetBinContent(ix+1,iy+1);
hspectrum->Fill(ix, iy, decoder->getValue(buff,ix,iy)-ped);
}
ich++;
}
}
delete [] buff;
//cout << "="<< endl;
ii++;
}
if (filebin.is_open())
filebin.close();
}
return hspectrum;
}
/**
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 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 xmin=1, int xmax=NC-1, int ymin=1, int ymax=NR-1, int cmsub=0, int hitfinder=1) {
double hc=0;
int sign=1;
moench03Ctb10GbData *decoder=new moench03Ctb10GbData();
cout << "decoder allocated " << endl;
moench03CommonMode *cmSub=NULL;
if (cmsub) {
cmSub=new moench03CommonMode(100);
cout << "common mode allocated " << endl;
} else {
cout << "non allocating common mode " << endl;
}
int iev=0;
int nph=0;
singlePhotonDetector<uint16_t> *filter=new singlePhotonDetector<uint16_t>(decoder, 3, 5, sign, cmSub, 100, 10);
cout << "filter allocated " << endl;
char *buff;
char fname[10000];
int nf=0;
eventType thisEvent=PEDESTAL;
// int iframe;
// double *data, ped, sigma;
// data=decoder->getCluster();
THStack *hs=new THStack("hs",fformat);
cout << "hstack allocated " << endl;
TH2F *h1=new TH2F("h1",tit,nbins,hmin-0.5,hmax-0.5,NC*NR,-0.5,NC*NR-0.5);
hs->Add(h1);
cout << "h1 allocated " << endl;
TH2F *h2;
TH2F *h3;
if (hitfinder) {
h2=new TH2F("h2",tit,nbins,hmin-0.5,hmax-0.5,NC*NR,-0.5,NC*NR-0.5);
cout << "h2 allocated " << endl;
h3=new TH2F("h3",tit,nbins,hmin-0.5,hmax-0.5,NC*NR,-0.5,NC*NR-0.5);
cout << "h3 allocated " << endl;
// 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);
}
if (hs->GetHists()) {
for (int i=0; i<3; i++)
if (hs->GetHists()->At(1)) cout << i << " " ;
cout << " histos allocated " << endl;
} else
cout << "no hists in stack " << endl;
ifstream filebin;
int ix=20, iy=20, ir, ic;
Int_t iFrame;
TTree *tall;
if (hitfinder)
tall=filter->initEventTree(tit, &iFrame);
#ifdef MY_DEBUG
cout << "debug mode " << endl;
TCanvas *myC;
TH2F *he;
TCanvas *cH1;
TCanvas *cH2;
TCanvas *cH3;
if (hitfinder) {
myC=new TCanvas("myc");
he=new TH2F("he","Event Mask",xmax-xmin, xmin, xmax, ymax-ymin, ymin, ymax);
he->SetStats(kFALSE);
he->Draw("colz");
cH1=new TCanvas("ch1");
cH1->SetLogz();
h1->Draw("colz");
cH2=new TCanvas("ch2");
cH2->SetLogz();
h2->Draw("colz");
cH3=new TCanvas("ch3");
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))) {
filter->newFrame();
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);
}
}
// if (hitfinder) {
// //calculate pedestals and common modes
// }
// 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);
// if (nf>10) {
h1->Fill(filter->getClusterTotal(1), iy+NR*ix);
#ifdef MY_DEBUG
// if (iev%10==0)
he->SetBinContent(ix+1-xmin, iy+1-ymin, (int)thisEvent);
#endif
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();
}
} // else {
// filter->addToPedestal(decoder->getValue(buff,ix,iy, cmsub));
// }
// }
}
//////////////////////////////////////////////////////////
#ifdef MY_DEBUG
// cout << iev << " " << h1->GetEntries() << " " << h2->GetEntries() << endl;
// if (iev%10==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);
//////////////////////////////////////////////////////////
#ifdef MY_DEBUG
myC->Modified();
myC->Update();
cH1->Modified();
cH1->Update();
cH2->Modified();
cH2->Update();
cH3->Modified();
cH3->Update();
#endif
delete decoder;
cout << "Read " << nf << " frames" << endl;
return hs;
}
TGraph* checkFrameNumber(char *fformat, int runmin, int runmax, int ix, int iy){
ifstream filebin;
char fname[10000];
moench03Ctb10GbData *decoder=new moench03Ctb10GbData();
char *buff;
int ii=0;
TGraph *g=new TGraph();
for (int irun=runmin; irun<=runmax; irun++) {
sprintf(fname,fformat,irun);
cout << fname << endl;
filebin.open((const char *)(fname), ios::in | ios::binary);
if (filebin.is_open()) {
while ((buff=decoder->readNextFrame(filebin))) {
g->SetPoint(ii,decoder->getFrameNumber(buff),decoder->getValue(buff,ix,iy));
ii++;
delete [] buff;
}
//cout << "="<< endl;
filebin.close();
} else
cout << "Could not open file " << fname << endl;
}
return g;
}

243
slsDetectorCalibration/readMoench03Data.C
View File

@ -2,48 +2,271 @@
void readMoench03Data(char *tit, int sign=1){
void readMoench03Data(std::string path,char* tit, std::string phase, std::string wtime,int sign=1,int runmin=0, int runmax=100, int sc_num=0,int hitfinder=1){
// int runmin(150);
// int runmax(200);
int nbins(2000);
int hmin(-1000);
int hmax(3000);
int xmin(1);
int xmax(399);
int ymin(1);
int ymax(399);
int cmsub(0);
char fname[1000];
TFile *fout;
// TFile *fouth1;
THStack *hs;
sprintf(fname,"/mnt/moenchnas/big_moench_xbox_20150223/Mo.root");
sprintf(fname,"%s/%s_%s_%s_%d_%d_Csub%d_HF%d.root",path.c_str(),tit,phase.c_str(),wtime.c_str(),runmin,runmax-1,cmsub,hitfinder);
fout=new TFile(fname,"RECREATE");
cerr<<"creating output file:"<<endl;
cerr<<fname<<endl;
cerr<<"/****************************/"<<endl;
sprintf(fname,"%s/%s_%s_%s_f0_%%d.raw",path.c_str(),tit,phase.c_str(),wtime.c_str()); // sprintf(fname,"%s/%s_phase%s_wtime%s_period0.075_OD%1.1f_f0_%%d.raw",path.c_str(),tit,phase.c_str(),wtime.c_str(),OD);
cerr<<fname<<endl;
// Int_t out = moench03DrawPedestals(fname,"Mo",runmin,runmax,nbins,hmin,hmax,xmin,xmax,ymin,ymax, cmsub,hitfinder);
fout->cd();
hs=moench03ReadData(fname,tit,runmin,runmax,nbins,hmin,hmax,xmin,xmax,ymin,ymax, cmsub,hitfinder);
sprintf(fname,"/mnt/moenchnas/big_moench_xbox_20150223/Mo_f0_%%d.raw");
hs=moench03ReadData(fname,"Mo",25133,25187,1500,-500,2500,1,399,1,399, 0,1);
cout << "returned" << endl;
hs->SetName(tit);
hs->SetTitle(tit);
cout << "name/title set" << endl;
Int_t nH=1;
TH2F * h=NULL;
if (hs->GetHists()) {
for (int i=0; i<3; i++) {
for (int i=0; i<nH; i++) {
if (hs->GetHists()->At(i)) {
cout << i << " " ;
(TH2F*)(hs->GetHists()->At(i))->Write();
// h->SetName(Form("h%d",i+1));
// h->SetTitle(Form("h%d",i+1));
// cerr<<h->GetEntries()<<" entries"<<endl;
// can->cd(i+1);
// h->Draw("colz");
// fout->cd();
// h->Write();
cout << i << " " ;
}
}
cout << " histos written " << endl;
} else
cout << "no hists in stack " << endl;
if(fout->IsOpen())fout->Close();
return;
}
/**********************************************************/
void readPixelCharge(std::string path,char* tit, std::string phase, std::string wtime, float OD, int sign,int runmin, int runmax, int x0, int y0){
char fname[1000];
sprintf(fname,"%s/%s_%s_%2.0fumSi_%s_f0_%%d.raw",path.c_str(),tit,phase.c_str(),OD,wtime.c_str());
char fnameout[1000];
sprintf(fnameout,"%s/%s_%s_%2.0fumSi_%s.root",path.c_str(),tit,phase.c_str(),OD,wtime.c_str());
// sprintf(fnameout,"%s/%s_phase%s_wtime%s_period0.2_OD%1.1f.root",path.c_str(),tit,phase.c_str(),wtime.c_str(),OD);
TFile * fout = new TFile(fnameout,"recreate");
for(int i=0; i<12; i++){
TH1F * hpix = SinglePixelHisto(fname,15000,-0.5,29999.5,runmin,runmax,x0+i,y0);
hpix->SetName(Form("h_%d_%d",x0+i,y0));
fout->cd();
hpix->Write();
}
return;
}
/*************************************************/
void LoopOnPixelCharge(std::string path,char* tit, std::string phase, std::string wtime, int sign,int runmin, int runmax, int x0, int y0){
float OD=700;
for(int i=0; i<6; i++){
readPixelCharge(path,tit,phase,wtime,OD,sign,runmin,runmax,x0,y0);
OD += 200;
}
return;
}
/**********************************************************/
void readPixelCorrelation(std::string path,char* tit, std::string frame, std::string phase, std::string wtime,int sign,int runmin, int runmax,int supercolumn){
int npx(400);
int npy(400);
char fname[1000];
int sc_width(25);
int sc_height(200);
int sc_number=supercolumn;
int xmin(0);
int xmax(0);
int ymin(0);
int ymax(0);
ADCMap * map = new ADCMap(npx,npy,sc_width,sc_height);
int ret = map->Init();
xmin=map->getXmin(sc_number);
xmax=map->getXmax(sc_number);
ymin=map->getYmin(sc_number);
ymax=map->getYmax(sc_number);
cerr<<"/**********************************/"<<endl;
cerr<<"Checking super column "<<sc_number<<endl;
cerr<<"x: "<<map->getXmin(sc_number)<<" - "<<map->getXmax(sc_number)<<endl;
cerr<<"y: "<<map->getYmin(sc_number)<<" - "<<map->getYmax(sc_number)<<endl;
cerr<<"/**********************************/"<<endl;
char fnameout[1000];
sprintf(fnameout,"%s/%s_wtime%s_%d_%d_Corr_SC%d.root",path.c_str(),tit,wtime.c_str(),runmin,runmax-1,sc_number);
TFile * fout = new TFile(fnameout,"recreate");
Int_t nbins(3000);
Float_t xlow(6999.5);
Float_t xup(9999.5);
TH2F * hcorr=NULL;
for(int i=0; i<25; i++){
phase.clear();
phase=(std::string)Form("%d",i*5);
cerr<<"check ADC phase "<<phase.c_str()<<endl;
sprintf(fname,"%s/%s_phase%s_wtime%s_period0.075_f0_%%d.raw",path.c_str(),tit,phase.c_str(),wtime.c_str());
// TH1F * h1 = new TH1F("h1",Form("ch%d",nchan1),nbins,xlow,xup);
// TH1F * h2 = new TH1F("h2",Form("ch%d",nchan2),nbins,xlow,xup);
hcorr= new TH2F(Form("hcorr_ph%s",phase.c_str()),Form("hcorr_ph%s",phase.c_str()),nbins,xlow,xup,nbins,xlow,xup);
DrawAllPixelCorrelation(fname,frame,runmin,runmax,xmin,xmax,ymin,ymax,sc_width,hcorr);
// h1->GetXaxis()->SetTitle("ADC");
// h2->GetXaxis()->SetTitle("ADC");
fout->cd();
// h1->Write();
// h2->Write();
hcorr->Write();
}
return;
fout->Close();
}
/****************************************************/
void readNoise(std::string path, char *tit, std::string flag, std::string phase, std::string wtime,float OD,int sign,int runmin, int runmax){
TFile * fout = new TFile(Form("%s/%s_%s_%s_noiseMap%s.root",path.c_str(),tit,phase.c_str(),wtime.c_str(),flag.c_str()),"recreate");
char fname[1000];
int nfiles=runmax-runmin;
// int runmin(382);
// int runmax(442);
// std::string target("Cu");
// THStack * hsnoise = NULL;
TH2F * hped=NULL;
TH2F * hnoise=NULL;
// for(int i=0; i<25; i++){
// phase.clear();
// phase=(std::string)Form("%d",i*5);
// sprintf(fname,"%s/%s_phase%s_wtime%s_period0.2_OD%1.1f_f0_%%d.raw",path.c_str(),tit,phase.c_str(),wtime.c_str(),OD);
// sprintf(fname,"%s/%s_phase%s_wtime%s_period0.2_f0_%%d.raw",path.c_str(),tit,phase.c_str(),wtime.c_str());
sprintf(fname,"%s/%s_%s_%s_f0_%%d.raw",path.c_str(),tit,phase.c_str(),wtime.c_str());
THStack * hsnoise=calcNoise(fname,flag,runmin,runmax,nfiles);
hsnoise->SetName(Form("%s_noiseMap_ph%s",tit,phase.c_str()));
fout->cd();
if (hsnoise->GetHists()) {
hped=(TH2F*)(hsnoise->GetHists()->At(0));
hped->SetName(Form("hped_ph%s",phase.c_str()));
hped->Write();
hnoise=(TH2F*)(hsnoise->GetHists()->At(1));
hnoise->SetName(Form("hnoise_ph%s",phase.c_str()));
hnoise->Write();
cout << " histos written for ADC Phase " << phase.c_str()<<endl;
} else
cout << "no hists in stack " << endl;
delete hsnoise;
// }
fout->Close();
return;
}
/************************************************************************/
void readFrames(std::string path, char *tit, std::string phase, std::string wtime,int sign,int runnum, int framemin, int framemax){
char fformat[1000];
char fname[1000];
sprintf(fformat,"%s/%s_phase%s_wtime%s_period0.075_f0_%d.raw",path.c_str(),tit,phase.c_str(),wtime.c_str(),runnum);
THStack * hs = DrawFrames(fformat,framemin,framemax);
int nframes=framemax-framemin;
TFile * fout = new TFile(Form("%s/Frames_phase%s_wtime%s_period0.075_f0_%d_fr%d-%d.root",path.c_str(),phase.c_str(),wtime.c_str(),runnum,framemin,framemax),"recreate");
TCanvas * c1= new TCanvas("c1","",800,600);
c1->Print(Form("%s/Frames_phase%s_wtime%s_period0.075_f0_%d_fr%d-%d.pdf[",path.c_str(),phase.c_str(),wtime.c_str(),runnum,framemin,framemax));
TH2F * h=NULL;
for(int hnum=0; hnum<nframes; hnum++){
h=(TH2F*)hs->GetHists()->At(hnum);
h->SetName(Form("h_%d",hnum+framemin));
h->SetTitle(Form("h_%d",hnum+framemin));
h->GetZaxis()->SetRangeUser(5000.,10000.);
h->Draw("colz");
c1->Print(Form("%s/Frames_phase%s_wtime%s_period0.075_f0_%d_fr%d-%d.pdf",path.c_str(),phase.c_str(),wtime.c_str(),runnum,framemin,framemax));
fout->cd();
h->Write();
}
c1->Print(Form("%s/Frames_phase%s_wtime%s_period0.075_f0_%d_fr%d-%d.pdf]",path.c_str(),phase.c_str(),wtime.c_str(),runnum,framemin,framemax));
fout->Close();
return;
}
/************************************************************************************/
void PlotSinglePixelHisto(std::string path, char *tit, std::string flag, std::string phase, std::string wtime,float OD,int sign,int runmin, int runmax, int x0, int y0){
char fname[1000];
sprintf(fname,"%s/%s_%s_%s_f0_%%d.raw",path.c_str(),tit,phase.c_str(),wtime.c_str());
int nbins(8000);
float xmin(-0.5);
float xmax(15999.5);
TH1F * h1 = SinglePixelHisto(fname,nbins,xmin,xmax,runmin,runmax,x0,y0);
h1->Draw("hist");
return;
}
// void raedNoiseDataN(char *tit, int sign=1){

66
slsDetectorCalibration/slsDetectorData.h
View File

@ -40,8 +40,8 @@ class slsDetectorData {
slsDetectorData(int npx, int npy, int dsize, int **dMap=NULL, dataType **dMask=NULL, int **dROI=NULL): nx(npx), ny(npy), dataSize(dsize) {
xmap=new int[nx*ny];
ymap=new int[nx*ny];
xmap=new int[dsize/sizeof(dataType)];
ymap=new int[dsize/sizeof(dataType)];
dataMask=new dataType*[ny];
for(int i = 0; i < ny; i++) {
@ -228,42 +228,44 @@ class slsDetectorData {
virtual dataType getChannel(char *data, int ix, int iy, int dr) {
dataType m=0;
uint64_t t;
int numBytes,divFactor,newix,pixelval;
//cout <<"ix:"<<ix<<" nx:"<<nx<<" iy:"<<ny<<" ny:"<<ny<<" datamap[iy][ix]:"<< dataMap[iy][ix] <<"datasize:"<< dataSize <<endl;
if (ix>=0 && ix<nx && iy>=0 && iy<ny && dataMap[iy][ix]>=0 && dataMap[iy][ix]<dataSize) {
m=dataMask[iy][ix];
/////All this stuff should go to eiger detector!!!!!!!!!!!!!!!
/* int numBytes,divFactor,newix,pixelval; */
numBytes = (nx * iy + ix);
divFactor=2;
if(dr == 4) divFactor = 16;
else if (dr == 8) divFactor = 8;
else if (dr == 16) divFactor = 4;
/* //cout <<"ix:"<<ix<<" nx:"<<nx<<" iy:"<<ny<<" ny:"<<ny<<" datamap[iy][ix]:"<< dataMap[iy][ix] <<"datasize:"<< dataSize <<endl; */
/* if (ix>=0 && ix<nx && iy>=0 && iy<ny && dataMap[iy][ix]>=0 && dataMap[iy][ix]<dataSize) { */
/* m=dataMask[iy][ix]; */
pixelval = numBytes % divFactor;
newix = ix - pixelval;
/* numBytes = (nx * iy + ix); */
/* divFactor=2; */
/* if(dr == 4) divFactor = 16; */
/* else if (dr == 8) divFactor = 8; */
/* else if (dr == 16) divFactor = 4; */
//cout <<"pixelval:"<<pixelval<<" newix:"<<newix<<endl;
//cout <<"64:"<< hex<<((uint64_t)(*((uint64_t*)(((char*)data)+(dataMap[iy][newix])))))<<endl;
t = (be64toh((uint64_t)(*((uint64_t*)(((char*)data)+(dataMap[iy][newix]))))));
//cout<<"t:"<<t<<endl;
/* pixelval = numBytes % divFactor; */
/* newix = ix - pixelval; */
}else
cprintf(RED,"outside limits\n");
/* //cout <<"pixelval:"<<pixelval<<" newix:"<<newix<<endl; */
/* //cout <<"64:"<< hex<<((uint64_t)(*((uint64_t*)(((char*)data)+(dataMap[iy][newix])))))<<endl; */
/* t = (be64toh((uint64_t)(*((uint64_t*)(((char*)data)+(dataMap[iy][newix])))))); */
/* //cout<<"t:"<<t<<endl; */
if(dr == 4)
//uint8_t value = t >> (pixelval*4); cout <<"value:"<< value << endl;
return ((t >> (pixelval*4)) & 0xf)^m;
else if(dr == 8)
//uint8_t value = t >> (pixelval*8); cout <<"value:"<< value << endl;
return ((t >> (pixelval*8)) & 0xff)^m;
else if(dr == 16){
//uint16_t value = t >> (pixelval*16); cout <<"value:"<< value << endl;
return ((t >> (pixelval*16)) & 0xffff)^m;
}else{
//uint32_t value = t >> (pixelval*32); cout <<"value:"<< value << endl;
return ((t >> (pixelval*32)) & 0xffffffff)^m;
}
/* }/\* else *\/ */
/* /\* cprintf(RED,"outside limits\n"); *\/ */
/* if(dr == 4) */
/* //uint8_t value = t >> (pixelval*4); cout <<"value:"<< value << endl; */
/* return ((t >> (pixelval*4)) & 0xf)^m; */
/* else if(dr == 8) */
/* //uint8_t value = t >> (pixelval*8); cout <<"value:"<< value << endl; */
/* return ((t >> (pixelval*8)) & 0xff)^m; */
/* else if(dr == 16){ */
/* //uint16_t value = t >> (pixelval*16); cout <<"value:"<< value << endl; */
/* return ((t >> (pixelval*16)) & 0xffff)^m; */
/* }else{ */
/* //uint32_t value = t >> (pixelval*32); cout <<"value:"<< value << endl; */
/* return ((t >> (pixelval*32)) & 0xffffffff)^m; */
/* } */
};
/**

2
slsDetectorCalibration/slsReceiverData.h
View File

@ -76,7 +76,7 @@ public:
while (dd<=(dsize-packetSize)) {
pnum=getPacketNumber(p);
fn=getFrameNumber(p);
//cout <<"pnum:"<<pnum<<" fn:"<<fn<<"\t";
cout <<"pnum:"<<pnum<<" fn:"<<fn<<"\t "<< np << endl;;
if (pnum<1 || pnum>nPackets) {
cout << "Bad packet number " << pnum << " frame "<< fn << endl;

178
slsDetectorCalibration/viewMoenchDRXRays.C Normal file
View File

@ -0,0 +1,178 @@
#define MYROOT1
#include <TH1D.h>
#include <TF1.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 <fstream>
#include "moench03CtbData.h"
//#include "moench03CommonMode.h"
//#include "singlePhotonDetector.h"
using namespace std;
THStack *viewMoenchDRXRays(int ix=70, int iy=88){
TF1 *poiss = new TF1("poiss", "[0]*TMath::Power(([1]/[2]),(x/[2]))*(TMath::Exp( ([1]/[2])))/TMath::Gamma((x/[2])+1)", 0, 5);
int thick[]={1700,1500,1300,1100,900,700};
int nt=6;
int nf=5;
char fname[1000];
char *data;
char tit[100];
TH1F *hh[nt], *hh3[nt], *hh5[nt];
TH2F *h2[nt], *hpix[nt];
Double_t val, val3, val5, val1;
int it=0;
Double_t ped[25];
TH1D *p;
THStack *hs=new THStack();
cout << nt << endl;
ifstream filebin;
moench03CtbData *decoder=new moench03CtbData();
sprintf(tit,"hpix_%dumSi_g1",thick[it]);
cout << tit << endl;
hpix[it]=new TH2F(tit,tit,2500,6000,16000,25,0,25);
hs->Add(hpix[it]);
sprintf(tit,"%dumSi_g1",thick[it]);
cout << tit << endl;
for (int iff=0; iff<nf; iff++) {
// cout << tit << " " << iff << endl;
sprintf(fname,"/mnt/moench/Moench03_MS_20150606/direct_beam_12.4keV_filter_scan/direct_beam_12.4keV_%s_400clk_f0_%d.raw",tit,iff);
filebin.open((const char *)(fname), ios::in | ios::binary);
if (filebin.is_open()) cout << "ok "<< fname << endl;
else cout << "could not open "<< fname << endl;
while ((data=decoder->readNextFrame(filebin))) {
for (int iiy=-2; iiy<3; iiy++)
for (int iix=-2; iix<3; iix++)
hpix[it]->Fill(decoder->getChannel(data,ix+iix,iy+iiy), (iiy+2)*5+iix+2);
delete [] data;
}
filebin.close();
cout << endl;
}
for (int iix=-2; iix<3; iix++) {
for (int iiy=-2; iiy<3; iiy++) {
cout << iix << " " << iiy << " " ;// <<endl;
p=hpix[0]->ProjectionX("p",(iiy+2)*5+iix+2+1,(2+iiy)*5+iix+2+1);
ped[(iiy+2)*5+iix+2]=p->GetBinCenter(p->GetMaximumBin());
cout << ped[(iiy+2)*5+iix+2] <<endl;
}
}
for (it=0; it<nt; it++) {
sprintf(tit,"hh_%dumSi_g1",thick[it]);
cout << tit << endl;
hh[it]=new TH1F(tit,tit,5000,0,10000);
hs->Add(hh[it]);
sprintf(tit,"hh3_%dumSi_g1",thick[it]);
hh3[it]=new TH1F(tit,tit,5000,0,30000);
hs->Add(hh3[it]);
sprintf(tit,"hh5_%dumSi_g1",thick[it]);
hh5[it]=new TH1F(tit,tit,5000,0,50000);
hs->Add(hh5[it]);
sprintf(tit,"%dumSi_g1",thick[it]);
cout << tit << endl;
hs->Add(hh[it]);
for (int iff=0; iff<nf; iff++) {
// cout << tit << " " << iff << endl;
sprintf(fname,"/mnt/moench/Moench03_MS_20150606/direct_beam_12.4keV_filter_scan/direct_beam_12.4keV_%s_400clk_f0_%d.raw",tit,iff);
filebin.open((const char *)(fname), ios::in | ios::binary);
if (filebin.is_open()) cout << "ok "<< fname << endl;
else cout << "could not open "<< fname << endl;
while ((data=decoder->readNextFrame(filebin))) {
cout << "-" ;
// for (int iy=0; iy<40; iy++)
// for (int ix=0; ix<350; ix++){
val1=0;
val3=0;
val5=0;
for (int iix=-2; iix<3; iix++) {
for (int iiy=-2; iiy<3; iiy++) {
// cout << iix << " " << iiy << " " ;// <<endl;
// p=hpix[]->ProjectionX("p",(iiy+2)*5+iix+2+1,(2+iiy)*5+iix+2+1);
// ped[it]=p->GetBinCenter(p->GetMaximumBin());
// cout << ped[it] <<endl;
val=decoder->getChannel(data,ix+iix,iy+iiy)-ped[(iiy+2)*5+iix+2];
if ((iix<-1 || iix>1) || (iiy<-1 || iiy>1))
val5+=val;
else if (iix!=0 || iiy!=0) {
val5+=val;
val3+=val;
} else {
val5+=val;
val3+=val;
val1+=val;
}
// if (iiy==1 && iix==0)
// h2[it]->Fill(val1,val);
}
}
hh5[it]->Fill(val5);
hh3[it]->Fill(val3);
hh[it]->Fill(val1);
delete [] data;
}
filebin.close();
cout << endl;
}
}
return hs;
}