improved readNextFrame for moench 10G streaming; implemented many photons finder

This commit is contained in:
bergamaschi 2017-05-30 12:21:48 +02:00
parent 73ec3903bf
commit 5dbfbdb82b
9 changed files with 524 additions and 235 deletions

View File

@ -29,11 +29,14 @@ class MovingStat
/**
clears the moving average number of samples parameter, mean and standard deviation
*/
void Set(double val)
void Set(double val, double rms=0)
{
m_n = n;
m_newM=val*n;
m_newM2=val*val*n;
if (rms<=0)
m_newM2=val*val*n;
else
m_newM2=(n*rms*rms+m_newM*m_newM/n);
}

View File

@ -44,8 +44,9 @@ class mythen3_01_jctbData : public slsDetectorData<short unsigned int> {
static short unsigned int* mythen03_frame(char *ptr, int dr=24, int nch=64*3, int off=5) {
// off=0;
int iarg;
int64_t word;
int64_t word, *wp;
short unsigned int* val=new short unsigned int[nch];
int bit[64];
int nb=2;
@ -53,30 +54,53 @@ class mythen3_01_jctbData : public slsDetectorData<short unsigned int> {
int idr=0;
int ib=0;
int iw=0;
int ii=0;
bit[0]=19;
bit[1]=8;
idr=0;
for (ib=0; ib<nch; ib++) {
val[ib]=0;
}
wp=(int64_t*)ptr;
for (iw=0; iw<nch/nb; iw) {
word=*((int64_t*)(ptr+8));
if (ioff<off) ioff++;
else {
if (idr<16) {
for (ib=0; ib<nb; ib++) {
if (word&(1<<bit[ib])) val[iw+nch*ib/nb]|=(1<<idr);
}//end for()
}
idr++;
if (idr==dr) {
idr=0;
iw++;
}//end if()
}//end else()
word=*wp;;
if (ioff<off) {
ioff++;
cout <<"*";
} else {
if (idr<16) {
for (ib=0; ib<nb; ib++) {
if (word&(1<<bit[ib])) {
cout << "+" ;
val[iw+nch*(ib/nb)]|=(1<<idr);
} else {
cout << "-" ;
}
}//end for()
}
idr++;
if (idr==dr) {
idr=0;
// cout << dec << " " << iw << " " << val[iw] << " " << val[iw+nch/2] << endl;
cout <<dec << iw<<endl;
iw++;
}//end if()
}//end else()
wp+=1;
ii++;
}//end for
cout << "Decoded "<<ii << " samples"<< endl;
cout << "Should be "<< nch/nb*dr+off << " samples"<< endl;
return val;
}
}
virtual int setFrameNumber(int f=0) {if (f>=0) frameNumber=f; return frameNumber; };
virtual int setDynamicRange(int d=-1) {if (d>0 && d<=24) dynamicRange=d; return dynamicRange;};

View File

@ -1,11 +1,12 @@
void imageMacro(char *name){
TH2D *imageMacro(char *name) {
//TH2D *makeNorm(){
//TFile ff("/local_zfs_raid/tomcat_20160528/trees/img_blank_eta_gmap.root");
//TH2D *hff=(TH2D*)ff.Get("imgHR");
TFile ff("/local_zfs_raid/tomcat_20160528/trees/img_blank_eta_nb25.root");
TFile *ff=new TFile("/mnt/moench_data/tomcat_20160528_img/img_blank_eta_nb25.root");
// TFile ff("/local_zfs_raid/tomcat_20160528/trees/img_blank_eta_gcorr_nb25.root");
TH2D *hff=(TH2D*)ff.Get("blankHR");
TH2D *hff=(TH2D*)ff->Get("blankHR");
hff->SetName("imgBlank");
TH2D *hpixel=new TH2D("hpixel","hpixel",25,0,25,25,0,25);
for (int ibx=10*25; ibx<hff->GetNbinsX()-10*25; ibx++) {
@ -45,9 +46,9 @@ void imageMacro(char *name){
} else {
sprintf(nn,"/local_zfs_raid/tomcat_20160528/trees/img_%s_eta_nb25.root", name);
sprintf(nn,"/mnt/moench_data/tomcat_20160528_img/img_%s_eta_nb25.root", name);
// if (strcmp(name,"blank"))
TFile fg(nn);
TFile *fg=new TFile(nn);
// else
// TFile fg=gDirectory;
@ -57,13 +58,21 @@ void imageMacro(char *name){
// TFile fg("/local_zfs_raid/tomcat_20160528/trees/img_sample_eta_nb25.root");
// TFile fg("/local_zfs_raid/tomcat_20160528/trees/img_grating_1d_eta_gcorr_nb25.root");
sprintf(nn,"%sHR",name);
TH2D *hg=(TH2D*)fg.Get(nn);
TH2D *hg=(TH2D*)fg->Get(nn);
// hg->SetName("imgGrating");
//hg->Divide(hff);
}
if (hpix)
hg->Divide(hpix);
new TCanvas();
hg->Draw("colz");
return hg;
}
void imageMacro(TH2D *hg){
Double_t imageData[200*400*25*25];
const int nsigma=5.;

View File

@ -1,7 +1,31 @@
{
TColor::InitializeColors();
const Int_t NRGBs = 5;
const Int_t NCont = 255;//90;
Double_t stops[NRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 };
Double_t red[NRGBs] = { 0.00, 0.00, 0.87, 1.00, 0.51 };
Double_t green[NRGBs] = { 0.00, 0.81, 1.00, 0.20, 0.00 };
Double_t blue[NRGBs] = { 0.51, 1.00, 0.12, 0.00, 0.00 };
Double_t gray[NRGBs] = { 1., 0.34, 0.61, 0.84, 1.00};
Double_t zero[NRGBs] = { 0., 0.0,0.0, 0.0, 0.00};
//TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
TColor::CreateGradientColorTable(NRGBs, stops, gray, gray, gray, NCont);
gStyle->SetNumberContours(NCont);
gStyle->SetPadTopMargin(0);
gStyle->SetPadRightMargin(0);
gStyle->SetPadBottomMargin(0);
gStyle->SetPadLeftMargin(0);
gROOT->ForceStyle();
// TFile ff("/local_zfs_raid/tomcat_20160528/trees/img_blank_eta.root");
// TH2D *hff=(TH2D*)ff.Get("imgHR");
TFile ff("/local_zfs_raid/tomcat_20160528/trees/img_blank_eta_gmap_v2.root");
TFile ff("/mnt/moench_data/tomcat_20160528_img/img_blank_eta_nb25.root");
TH2D *hff=(TH2D*)ff.Get("blankHR");
hff->SetName("imgBlank");
TH2D *hpixel=new TH2D("hpixel","hpixel",25,0,25,25,0,25);
@ -32,7 +56,7 @@
// TFile fg("/local_zfs_raid/tomcat_20160528/trees/img_grating_2d_eta.root");
// TH2D *hg=(TH2D*)fg.Get("imgHR");
TFile fg("/local_zfs_raid/tomcat_20160528/trees/img_grating_2d_eta_gmap_v2.root");
TFile fg("/mnt/moench_data/tomcat_20160528_img/img_grating_2d_eta_nb25.root");
TH2D *hg=(TH2D*)fg.Get("grating_2dHR");
hg->SetName("imgGrating");
@ -50,9 +74,21 @@
hg->Scale(1/ng);
new TCanvas();
hg->SetMaximum(1.);
hg->Draw("colz");
new TCanvas("c1","c1",800,800);
hg->SetMaximum(2.);
Double_t xmin=hg->GetXaxis()->GetXmin();
Double_t xmax=hg->GetXaxis()->GetXmax();
Double_t ymin=hg->GetYaxis()->GetXmin();
Double_t ymax=hg->GetYaxis()->GetXmax();
hg->SetTitle(";mm;mm;Normalized intensity (a.u.)");
hg->GetXaxis()->Set(hg->GetXaxis()->GetNbins(),xmin*0.025-180.*0.025+0.1,xmax*0.025-180.*0.025+0.1);
hg->GetYaxis()->Set(hg->GetYaxis()->GetNbins(),ymin*0.025-180.*0.025-1.5,ymax*0.025-180*0.025-1.5);
hg->GetXaxis()->SetRangeUser(0,0.5);
hg->GetYaxis()->SetRangeUser(0,0.5);
hg->Draw("col");
hg->SetStats(kFALSE);
// new TCanvas();
// hsg->SetMaximum(1.);
// hsg->Draw("colz");

View File

@ -35,7 +35,7 @@
// TFile fg("/local_zfs_raid/tomcat_20160528/trees/img_grating_2d_eta.root");
// TH2D *hg=(TH2D*)fg.Get("imgHR");
TFile fg("/local_zfs_raid/tomcat_20160528/trees/img_sample_eta_gmap_v2.root");
TFile fg("/mnt/moench_data/tomcat_20160528_img/img_sample_eta_gmap_v2.root");
TH2D *hg=(TH2D*)fg.Get("sampleHR");
hg->SetName("imgSample");

View File

@ -6,12 +6,13 @@
class moench03Ctb10GbData : public slsReceiverData<uint16_t> {
private:
protected:
int iframe;
int nadc;
int sc_width;
int sc_height;
public:
@ -165,58 +166,116 @@ class moench03Ctb10GbData : public slsReceiverData<uint16_t> {
/* }; */
virtual char *readNextFrame(ifstream &filebin, int& fnum) {
char *data=new char[packetSize*nPackets];
char *retval=0;
int np=0, nd;
fnum = -1;
int pn;
char aa[8224];
char *packet=(char *)aa;
virtual char *readNextFrame(ifstream &filebin, int& fnum, char *data=NULL) {
int dd=0;
if (data==NULL) {
data=new char[packetSize*nPackets];
dd=1;
}
char *retval=0;
int np=0, nd;
fnum = -1;
int pn;
char aa[8224];
char *packet=(char *)aa;
int place;
if (filebin.is_open()) {
if (filebin.is_open()) {
place=filebin.tellg();
while(filebin.read((char*)packet, 8208) && np<nPackets){
pn=getPacketNumber(packet);
if (pn==1 && fnum<0)
fnum= getFrameNumber(packet);
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++;
// cout <<getFrameNumber(packet)<< " fn: " << fnum << "\t pn: " << pn << " np " << np << endl;
if (fnum>=0) {
if (getFrameNumber(packet) !=fnum) {
cout << "-"<<endl;
filebin.seekg(place);
if (np==0){
if (dd)
delete [] data;
return NULL;
} else {
return data;
}
}
}
memcpy(data+(pn-1)*packetSize, packet, packetSize);
np++;
if (np==nPackets)
return data;
if (np==0){
delete [] data;
return NULL;
}
}
place=filebin.tellg();
}
};
}
if (np==0){
cout << "-"<<endl;
filebin.seekg(place);
if (dd)
delete [] data;
return NULL;
}
//filebin.seekg(place);
return data;
};
/* virtual char *readNextFrame(ifstream &filebin, int& fnum) { */
/* char *data=new char[packetSize*nPackets]; */
/* char *retval=0; */
/* int np=0, nd; */
/* fnum = -1; */
/* int pn; */
/* char aa[8224]; */
/* char *packet=(char *)aa; */
/* if (filebin.is_open()) { */
/* while(filebin.read((char*)packet, 8208) && np<nPackets){ */
/* pn=getPacketNumber(packet); */
/* if (pn==1 && fnum<0) */
/* fnum= getFrameNumber(packet); */
/* // cout << "fn: " << fnum << "\t pn: " << pn << endl; */
/* if (fnum>=0) { */
/* if (getFrameNumber(packet) !=fnum) { */
/* if (np==0){ */
/* delete [] data; */
/* return NULL; */
/* } else */
/* return data; */
/* } */
/* memcpy(data+(pn-1)*packetSize, packet, packetSize); */
/* np++; */
/* } */
/* } */
/* } */
/* if (np==0){ */
/* delete [] data; */
/* return NULL; */
/* } */
/* }; */
int getPacketNumber(int x, int y) {return dataMap[y][x]/8208;};
virtual char *readNextFrame(ifstream &filebin) {
int fnum;
return readNextFrame(filebin, fnum);
};
};

View File

@ -1,17 +1,12 @@
#ifndef MOENCH03TCTB10GBDATA_H
#define MOENCH03TCTB10GBDATA_H
#include "slsReceiverData.h"
#include "moench03Ctb10GbData.h"
class moench03TCtb10GbData : public slsReceiverData<uint16_t> {
class moench03TCtb10GbData : public moench03Ctb10GbData { //slsReceiverData<uint16_t> {
private:
int iframe;
int nadc;
int sc_width;
int sc_height;
public:
@ -26,8 +21,9 @@ class moench03TCtb10GbData : public slsReceiverData<uint16_t> {
// moench03TCtb10GbData(int ns=5000): slsDetectorData<uint16_t>(400, 400, 8208*40, NULL, NULL) , nadc(32), sc_width(25), sc_height(200) {
moench03TCtb10GbData(int ns=5000): slsReceiverData<uint16_t>(400, 400, 40, 8208), nadc(32), sc_width(25), sc_height(200) {
moench03TCtb10GbData(int ns=5000): moench03Ctb10GbData(ns) {//slsReceiverData<uint16_t>(400, 400, 40, 8208), nadc(32), sc_width(25), sc_height(200) {
// cout <<"constructor"<< endl;
// nadc=32;
int adc_nr[32]={300,325,350,375,300,325,350,375, \
200,225,250,275,200,225,250,275,\
100,125,150,175,100,125,150,175,\
@ -47,11 +43,13 @@ class moench03TCtb10GbData : public slsReceiverData<uint16_t> {
for (iadc=0; iadc<nadc; iadc++) {
i=128*ip+is;
adc4=(int)iadc/4;
//cout << i << endl;
adc4=(int)iadc/4;
if (i<sc_width*sc_height) {
// for (int i=0; i<sc_width*sc_height; i++) {
col=adc_nr[iadc]+(i%sc_width);
if (adc4%2==0) {
// if (iadc<(nadc/2)) {
row=199-i/sc_width;
} else {
row=200+i/sc_width;
@ -63,181 +61,221 @@ class moench03TCtb10GbData : public slsReceiverData<uint16_t> {
}
}
}
// cout <<"map"<< endl;
int ipacket;
int ibyte;
int ii=0;
for (int ipacket=0; ipacket<npackets; ipacket++) {
for (int ibyte=0; ibyte< 8208/2; ibyte++) {
i=ipacket*8208/2+ibyte;
if (ibyte<8) {
//header!
xmap[i]=-1;
ymap[i]=-1;
} else {
// ii=ibyte+128*32*ipacket;
isample=ii/nadc;
iadc=ii%nadc;
adc4 = (int)iadc/4;
i=ipacket*8208/2+ibyte;
// cout << i << " ";
if (ibyte<8) {
//header!
xmap[i]=-1;
ymap[i]=-1;
} else {
// ii=ibyte+128*32*ipacket;
// cout << nadc << endl;
isample=ii/nadc;
iadc=ii%nadc;
ix=isample%sc_width;
iy=isample/sc_width;
if (adc4%2==0) {
xmap[i]=adc_nr[iadc]+ix;
ymap[i]=ny/2-1-iy;
adc4 = (int)iadc/4;
ix=isample%sc_width;
iy=isample/sc_width;
// cout << isample << " " << iadc << " " << ix << " " << iy ;//<< endl;
if (adc4%2==0) {//iadc<(nadc/2)) {
//if (iadc<(nadc/2)) {
xmap[i]=adc_nr[iadc]+ix;
ymap[i]=ny/2-1-iy;
} else {
xmap[i]=adc_nr[iadc]+ix;
ymap[i]=ny/2+iy;
}
} else {
xmap[i]=adc_nr[iadc]+ix;
ymap[i]=ny/2+iy;
}
ii++;
}
ii++;
}
}
}
/* int ipacket; */
/* int ibyte; */
/* int ii=0; */
/* for (int ipacket=0; ipacket<npackets; ipacket++) { */
/* for (int ibyte=0; ibyte< 8208/2; ibyte++) { */
/* i=ipacket*8208/2+ibyte; */
/* cout << i << endl; */
/* if (ibyte<8) { */
/* //header! */
/* xmap[i]=-1; */
/* ymap[i]=-1; */
/* } else { */
/* ii=ibyte+128*32*ipacket; */
/* isample=ii/nadc; */
/* iadc=ii%nadc; */
/* adc4 = (int)iadc/4; */
/* ix=isample%sc_width; */
/* iy=isample/sc_width; */
/* if (adc4%2==0) { */
/* xmap[i]=adc_nr[iadc]+ix; */
/* ymap[i]=ny/2-1-iy; */
/* } else { */
/* xmap[i]=adc_nr[iadc]+ix; */
/* ymap[i]=ny/2+iy; */
/* } */
/* ii++; */
/* } */
/* } */
/* } */
// cout <<"done"<< endl;
iframe=0;
// cout << "data struct created" << endl;
};
/**
/* /\** */
Returns the frame number for the given dataset. Purely virtual func.
\param buff pointer to the dataset
\returns frame number
*/
int getFrameNumber(char *buff){return *((int*)buff)&0xffffffff;};
/**
Returns the packet number for the given dataset. purely virtual func
\param buff pointer to the dataset
\returns packet number number
*/
int getPacketNumber(char *buff){return ((*(((int*)(buff+4))))&0xff)+1;};
/* /\** */
/* Loops over a memory slot until a complete frame is found (i.e. all packets 0 to nPackets, same frame number). purely virtual func */
/* \param data pointer to the memory to be analyzed */
/* \param ndata reference to the amount of data found for the frame, in case the frame is incomplete at the end of the memory slot */
/* \param dsize size of the memory slot to be analyzed */
/* \returns pointer to the beginning of the last good frame (might be incomplete if ndata smaller than dataSize), or NULL if no frame is found */
/* Returns the frame number for the given dataset. Purely virtual func. */
/* \param buff pointer to the dataset */
/* \returns frame number */
/* *\/ */
/* virtual char *findNextFrame(char *data, int &ndata, int dsize){ndata=dsize; setDataSize(dsize); return data;}; */
/* /\** */
/* int getFrameNumber(char *buff){return *((int*)buff)&0xffffffff;}; */
/* /\** */
/* Returns the packet number for the given dataset. purely virtual func */
/* \param buff pointer to the dataset */
/* \returns packet number number */
/* Loops over a file stream until a complete frame is found (i.e. all packets 0 to nPackets, same frame number). Can be overloaded for different kind of detectors! */
/* \param filebin input file stream (binary) */
/* \returns pointer to the begin of the last good frame, NULL if no frame is found or last frame is incomplete */
/* *\/ */
/* virtual char *readNextFrame(ifstream &filebin){ */
/* // int afifo_length=0; */
/* uint16_t *afifo_cont; */
/* int ib=0; */
/* if (filebin.is_open()) { */
/* afifo_cont=new uint16_t[dataSize/2]; */
/* while (filebin.read(((char*)afifo_cont)+ib,2)) { */
/* ib+=2; */
/* if (ib==dataSize) break; */
/* } */
/* if (ib>0) { */
/* iframe++; */
/* // cout << ib << "-" << endl; */
/* return (char*)afifo_cont; */
/* } else { */
/* delete [] afifo_cont; */
/* return NULL; */
/* } */
/* } */
/* return NULL; */
/* }; */
/* int getPacketNumber(char *buff){return ((*(((int*)(buff+4))))&0xff)+1;}; */
/* /\* /\\** *\/ */
/* /\* Loops over a memory slot until a complete frame is found (i.e. all packets 0 to nPackets, same frame number). purely virtual func *\/ */
/* /\* \param data pointer to the memory to be analyzed *\/ */
/* /\* \param ndata reference to the amount of data found for the frame, in case the frame is incomplete at the end of the memory slot *\/ */
/* /\* \param dsize size of the memory slot to be analyzed *\/ */
/* /\* \returns pointer to the beginning of the last good frame (might be incomplete if ndata smaller than dataSize), or NULL if no frame is found *\/ */
/* /\* *\\/ *\/ */
/* /\* virtual char *findNextFrame(char *data, int &ndata, int dsize){ndata=dsize; setDataSize(dsize); return data;}; *\/ */
virtual char *readNextFrame(ifstream &filebin, int& ff, int &np) {
char *data=new char[packetSize*nPackets];
char *retval=0;
int nd;
np=0;
int pn;
char aa[8224]={0};
char *packet=(char *)aa;
int fnum = -1;
/* /\* /\\** *\/ */
if (ff>=0)
fnum=ff;
/* /\* 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 *\/ */
if (filebin.is_open()) {
/* /\* *\\/ *\/ */
/* /\* 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; *\/ */
/* /\* }; *\/ */
cout << "+";
/* virtual char *readNextFrame(ifstream &filebin, int& ff, int &np) { */
/* char *data=new char[packetSize*nPackets]; */
/* char *retval=0; */
/* int nd; */
/* np=0; */
/* int pn; */
/* char aa[8224]={0}; */
/* char *packet=(char *)aa; */
/* int fnum = -1; */
while(filebin.read((char*)packet, 8208)){
/* if (ff>=0) */
/* fnum=ff; */
pn=getPacketNumber(packet);
/* if (filebin.is_open()) { */
if (fnum<0)
fnum= getFrameNumber(packet);
if (fnum>=0) {
if (getFrameNumber(packet) !=fnum) {
/* cout << "+"; */
if (np==0){
cout << "-";
delete [] data;
return NULL;
} else
filebin.seekg(-8208,ios_base::cur);
/* while(filebin.read((char*)packet, 8208)){ */
return data;
}
/* pn=getPacketNumber(packet); */
memcpy(data+(pn-1)*packetSize, packet, packetSize);
np++;
if (np==nPackets)
break;
if (pn==nPackets)
break;
}
}
/* if (fnum<0) */
/* fnum= getFrameNumber(packet); */
}
/* if (fnum>=0) { */
/* if (getFrameNumber(packet) !=fnum) { */
if (np==0){
cout << "?";
delete [] data;
return NULL;
}// else if (np<nPackets)
// cout << "Frame " << fnum << " lost " << nPackets-np << " packets " << endl;
/* if (np==0){ */
/* cout << "-"; */
/* delete [] data; */
/* return NULL; */
/* } else */
/* filebin.seekg(-8208,ios_base::cur); */
return data;
/* return data; */
/* } */
};
/* memcpy(data+(pn-1)*packetSize, packet, packetSize); */
/* np++; */
/* if (np==nPackets) */
/* break; */
/* if (pn==nPackets) */
/* break; */
/* } */
/* } */
int getPacketNumber(int x, int y) {return dataMap[y][x]/8208;};
/* } */
virtual char *readNextFrame(ifstream &filebin) {
int fnum=-1, np;
return readNextFrame(filebin, fnum, np);
};
/* if (np==0){ */
/* cout << "?"; */
/* delete [] data; */
/* return NULL; */
/* }// else if (np<nPackets) */
/* // cout << "Frame " << fnum << " lost " << nPackets-np << " packets " << endl; */
/* return data; */
/* }; */
/* int getPacketNumber(int x, int y) {return dataMap[y][x]/8208;}; */
/* virtual char *readNextFrame(ifstream &filebin) { */
/* int fnum=-1, np; */
/* return readNextFrame(filebin, fnum, np); */
/* }; */
};

View File

@ -40,7 +40,7 @@ class pedestalSubtraction {
virtual int SetNPedestals(int i=-1) {if (i>0) stat.SetN(i); return stat.GetN();};
/** sets the moving average */
virtual void setPedestal(double val) {stat.Set(val);}
virtual void setPedestal(double val, double rms=0) {stat.Set(val, rms);}

View File

@ -31,6 +31,8 @@ using namespace std;
UNDEFINED_EVENT=-1
};
#ifndef DEF_QUAD
#define DEF_QUAD
enum quadrant {
TOP_LEFT=0,
TOP_RIGHT=1,
@ -38,7 +40,7 @@ using namespace std;
BOTTOM_RIGHT=3,
UNDEFINED_QUADRANT=-1
};
#endif
template <class dataType>
@ -158,7 +160,7 @@ class singlePhotonDetector {
\param iy pixel y coordinate
\param val value to set
*/
virtual void setPedestal(int ix, int iy, double val){if (ix>=0 && ix<nx && iy>=0 && iy<ny) stat[iy][ix].setPedestal(val);};
virtual void setPedestal(int ix, int iy, double val, double rms=0){if (ix>=0 && ix<nx && iy>=0 && iy<ny) stat[iy][ix].setPedestal(val,rms);};
@ -190,6 +192,121 @@ class singlePhotonDetector {
int *getNPhotons(char *data, double thr=-1) {
double val;
int *nph=new int[nx*ny];
double rest[ny][nx];
int cy=(clusterSizeY+1)/2;
int cs=(clusterSize+1)/2;
int ccs=clusterSize;
int ccy=clusterSizeY;
double tthr=thr;
int nn;
double max=0, tl=0, tr=0, bl=0,br=0, v;
if (thr>=0) {
cy=1;
cs=1;
ccs=1;
ccy=1;
}
for (int ix=0; ix<nx; ix++) {
for (int iy=0; iy<ny; iy++) {
val=dataSign*(det->getValue(data, ix, iy)-getPedestal(ix,iy,0));
if (thr<=0) tthr=nSigma*getPedestalRMS(ix,iy);
nn=val/tthr;
if (nn>0) {
rest[iy][ix]=val-nn*tthr;
nph[ix+nx*iy]=nn;
} else {
rest[iy][ix]=val;
nph[ix+nx*iy]=0;
}
}
}
for (int ix=clusterSize/2; ix<clusterSize/2-1; ix++) {
for (int iy=clusterSizeY/2; iy<ny-clusterSizeY/2; iy++) {
eventMask[iy][ix]=PEDESTAL;
if (thr<=0) tthr=nSigma*getPedestalRMS(ix,iy);
max=0;
tl=0;
tr=0;
bl=0;
br=0;
for (int ir=-(clusterSizeY/2); ir<(clusterSizeY/2)+1; ir++) {
for (int ic=-(clusterSize/2); ic<(clusterSize/2)+1; ic++) {
if ((iy+ir)>=0 && (iy+ir)<ny && (ix+ic)>=0 && (ix+ic)<nx) {
//cluster->set_data(rest[iy+ir][ix+ic], ic, ir);
v=rest[iy+ir][ix+ic];//cluster->get_data(ic,ir);
tot+=v;
if (ir<=0 && ic<=0)
bl+=v;
if (ir<=0 && ic>=0)
br+=v;
if (ir>=0 && ic<=0)
tl+=v;
if (ir>=0 && ic>=0)
tr+=v;
if (v>max) {
max=v;
}
if (ir==0 && ic==0) {
if (v>tthr) {
eventMask[iy][ix]=PHOTON;
}
}
}
}
}
//if (cluster->get_data(0,0)>=max) {
if (rest[iy][ix]>=max) {
if (bl>=br && bl>=tl && bl>=tr) {
quad=BOTTOM_LEFT;
quadTot=bl;
} else if (br>=bl && br>=tl && br>=tr) {
quad=BOTTOM_RIGHT;
quadTot=br;
} else if (tl>=br && tl>=bl && tl>=tr) {
quad=TOP_LEFT;
quadTot=tl;
} else if (tr>=bl && tr>=tl && tr>=br) {
quad=TOP_RIGHT;
quadTot=tr;
}
if (rest[iy][ix]>tthr || tot>sqrt(ccy*ccs)*tthr || quadTot>sqrt(cy*cs)*tthr) {
nph[ix+nx*iy]++;
rest[iy][ix]-=tthr;
}
}
}
}
return nph;
}
/** finds event type for pixel and fills cluster structure. The algorithm loops only if the evenMask for this pixel is still undefined.
if pixel or cluster around it are above threshold (nsigma*pedestalRMS) cluster is filled and pixel mask is PHOTON_MAX (if maximum in cluster) or NEIGHBOUR; If PHOTON_MAX, the elements of the cluster are also set as NEIGHBOURs in order to speed up the looping
@ -207,6 +324,11 @@ class singlePhotonDetector {
double max=0, tl=0, tr=0, bl=0,br=0, v;
// cout << iframe << endl;
int cy=(clusterSizeY+1)/2;
int cs=(clusterSize+1)/2;
tot=0;
quadTot=0;
quad=UNDEFINED_QUADRANT;
@ -250,31 +372,17 @@ class singlePhotonDetector {
if (ir>=0 && ic>=0)
tr+=v;
if (cluster->get_data(ic,ir)>max) {
if (v>max) {
max=v;
}
if (ir==0 && ic==0) {
if (v>nSigma*cluster->rms) {
eventMask[iy][ix]=PHOTON;
} else if (cluster->get_data(ic,ir)<-nSigma*cluster->rms)
if (v<-nSigma*cluster->rms)
eventMask[iy][ix]=NEGATIVE_PEDESTAL;
}
}
}
}
if (eventMask[iy][ix]!=PHOTON && tot>sqrt(clusterSizeY*clusterSize)*nSigma*cluster->rms) {
eventMask[iy][ix]=NEIGHBOUR;
} else if (eventMask[iy][ix]==PHOTON) {
if (cluster->get_data(0,0)>=max) {
eventMask[iy][ix]=PHOTON_MAX;
}
} else if (eventMask[iy][ix]==PEDESTAL) {
if (cm==0)
addToPedestal(det->getValue(data, ix, iy),ix,iy);
}
if (bl>=br && bl>=tl && bl>=tr) {
quad=BOTTOM_LEFT;
quadTot=bl;
@ -289,6 +397,18 @@ class singlePhotonDetector {
quadTot=tr;
}
if (max>nSigma*cluster->rms || tot>sqrt(clusterSizeY*clusterSize)*nSigma*cluster->rms || quadTot>cy*cs*nSigma*cluster->rms) {
if (cluster->get_data(0,0)>=max) {
eventMask[iy][ix]=PHOTON_MAX;
} else {
eventMask[iy][ix]=PHOTON;
}
} else if (eventMask[iy][ix]==PEDESTAL) {
if (cm==0)
addToPedestal(det->getValue(data, ix, iy),ix,iy);
}
return eventMask[iy][ix];