detectors/slsDetectorPackage
25
0
Fork 0
mirror of https://github.com/slsdetectorgroup/slsDetectorPackage.git synced 2025-05-12 07:30:03 +02:00
Code Issues Actions Packages Releases Activity

common mode subtraction implemented

Browse Source 
git-svn-id: file:///afs/psi.ch/project/sls_det_software/svn/slsDetectorCalibration@9 113b152e-814d-439b-b186-022a431db7b5
This commit is contained in:
bergamaschi 2013-11-25 11:03:04 +00:00
parent 37196f14ca
commit 4f7c6f633f
2 changed files with 135 additions and 41 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

65
slsDetectorCalibration/moench02ModuleData.h
View File

@ -109,8 +109,55 @@ class moench02ModuleData : public slsDetectorData {
}; };
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 (getEventType(ix, iy, hc, tc, 1)==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]);
}
};
double getCommonMode(int ix, int iy) {
int isc=ix/40;
if (nCm[isc]>0)
return cmPed[isc]/nCm[isc]-cmStat[isc].Mean();
else
return 0;
};
char *readNextFrame(ifstream &filebin) { char *readNextFrame(ifstream &filebin) {
if (oldbuff) if (oldbuff)
delete [] oldbuff; delete [] oldbuff;
@ -223,10 +270,7 @@ class moench02ModuleData : public slsDetectorData {
} else if (v<-nSigma*getPedestalRMS(ix,iy)) } else if (v<-nSigma*getPedestalRMS(ix,iy))
ret=NEGATIVE_PEDESTAL; ret=NEGATIVE_PEDESTAL;
} }
} }
} }
} }
@ -238,14 +282,19 @@ class moench02ModuleData : public slsDetectorData {
} }
cx=ix;
cy=iy;
return ret; return ret;
}; };
double getClusterElement(int ic, int ir){return cluster[ic+1][ir+1];}; double getClusterElement(int ic, int ir, int cm=0){
if (cm) cluster[ic+1][ir+1]-getCommonMode(cx,cy);
else return cluster[ic+1][ir+1];
};
double *getCluster(){return &cluster[0][0];}; double *getCluster(){return &cluster[0][0];};
@ -271,8 +320,12 @@ class moench02ModuleData : public slsDetectorData {
char *buff; char *buff;
char *oldbuff; char *oldbuff;
int pedSub; int pedSub;
int cmSub;
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;
}; };

103
slsDetectorCalibration/moenchReadData.C
View File

@ -16,15 +16,35 @@
#include <fstream> #include <fstream>
#include "moench02ModuleData.h" #include "moench02ModuleData.h"
//#include "MovingStat.h"
using namespace std; using namespace std;
#define NC 160 #define NC 160
#define NR 160 #define NR 160
/**
char *fformat, file name format
char *tit, title of the tree etc.
int runmin, minimum run number
int runmax, max run number
int nbins=1500, number of bins for spectrum hists
int hmin=-500, minimum for spectrum hists
int hmax=1000, maximum for spectrum hists
int sign=1, sign of the spectrum to find hits
double hc=0, readout correlation coefficient with previous pixel
double tc=0, time correlation coefficient with previous frame (case of bad reset)
int xmin=0, minimum x coordinate
int xmax=NC, maximum x coordinate
int ymin=0, minimum y coordinate
int ymax=NR, maximum y coordinate
int pedsub=1, enable pedestal subtraction
int cmsub=0 enable commonmode subtraction
*/
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, double tc=0, int xmin=0, int xmax=NC, int ymin=0, int ymax=NR, int pedsub=1) {
moench02ModuleData *decoder=new moench02ModuleData(); moench02ModuleData *decoder=new moench02ModuleData();
char *buff; char *buff;
@ -33,7 +53,6 @@ THStack *moenchReadData(char *fformat, char *tit, int runmin, int runmax, int nb
double oldval; double oldval;
int nf=0; int nf=0;
moench02ModuleData::eventType thisEvent=moench02ModuleData::PEDESTAL; moench02ModuleData::eventType thisEvent=moench02ModuleData::PEDESTAL;
// int iframe; // int iframe;
@ -82,6 +101,8 @@ THStack *moenchReadData(char *fformat, char *tit, int runmin, int runmax, int nb
double tot, tl, tr, bl, br, v; double tot, tl, tr, bl, br, v;
int scmin=xmin/40, scmax=(xmax-1)/40+1;
// 6% x-talk from previous pixel // 6% x-talk from previous pixel
// 12% x-talk from previous frame // 12% x-talk from previous frame
@ -93,19 +114,39 @@ THStack *moenchReadData(char *fformat, char *tit, int runmin, int runmax, int nb
while ((buff=decoder->readNextFrame(filebin))) { while ((buff=decoder->readNextFrame(filebin))) {
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);
}
}
}
for (ix=xmin-1; ix<xmax+1; ix++) for (ix=xmin-1; ix<xmax+1; ix++)
for (iy=ymin-1; iy<ymax+1; iy++) { for (iy=ymin-1; iy<ymax+1; iy++) {
//////////////////////////////////////////////////////// ////////////////////////////////////////////////////////
if (pedsub) { if (pedsub) {
if (nf>100) {
thisEvent= decoder->getEventType(ix, iy, hc, tc, 1);
} if (nf>100) {
thisEvent= decoder->getEventType(ix, iy, hc, tc, 1);
}
if (thisEvent==moench02ModuleData::PEDESTAL) {
decoder->addToPedestal( decoder->getChannelShort(ix, iy, hc, tc), ix, iy); if (thisEvent==moench02ModuleData::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);
}
/*********************************************************** /***********************************************************
@ -118,14 +159,14 @@ Add here the function that you want to call: fill histos, make trees etc.
tr=0; tr=0;
bl=0; bl=0;
br=0; br=0;
h1->Fill(decoder->getClusterElement(0,0), iy+NR*ix); h1->Fill(decoder->getClusterElement(0,0, cmsub), iy+NR*ix);
if (thisEvent==moench02ModuleData::PHOTON_MAX ) { if (thisEvent==moench02ModuleData::PHOTON_MAX ) {
for (ir=-1; ir<2; ir++) { for (ir=-1; ir<2; ir++) {
for (ic=-1; ic<2; ic++) { for (ic=-1; ic<2; ic++) {
v=decoder->getClusterElement(ic,ir); v=decoder->getClusterElement(ic,ir,cmsub);
data[ic+1][ir+1]=v; data[ic+1][ir+1]=v;
tot+=v; tot+=v;
if (ir<1) { if (ir<1) {
@ -154,34 +195,34 @@ Add here the function that you want to call: fill histos, make trees etc.
if (bl>br && bl>tl && bl>tr) { if (bl>br && bl>tl && bl>tr) {
h2->Fill(bl, iy+NR*ix); h2->Fill(bl, iy+NR*ix);
if (bl>0) { if (bl>0) {
hetaX->Fill((decoder->getClusterElement(0,0)+decoder->getClusterElement(0,-1))/bl,iy+NR*ix); hetaX->Fill((decoder->getClusterElement(0,0,cmsub)+decoder->getClusterElement(0,-1,cmsub))/bl,iy+NR*ix);
hetaY->Fill((decoder->getClusterElement(0,0)+decoder->getClusterElement(-1,0))/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)+decoder->getClusterElement(0,-1))/bl,iy+NR*(ix-1)); 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)+decoder->getClusterElement(-1,0))/bl,(iy-1)+NR*ix); hetaY->Fill(1.-(decoder->getClusterElement(0,0,cmsub)+decoder->getClusterElement(-1,0,cmsub))/bl,(iy-1)+NR*ix);
} }
} else if (br>bl && br>tl && br>tr) { } else if (br>bl && br>tl && br>tr) {
h2->Fill(br, iy+NR*ix); h2->Fill(br, iy+NR*ix);
if (br>0) { if (br>0) {
hetaX->Fill((decoder->getClusterElement(0,0)+decoder->getClusterElement(0,-1))/br,iy+NR*ix); hetaX->Fill((decoder->getClusterElement(0,0,cmsub)+decoder->getClusterElement(0,-1,cmsub))/br,iy+NR*ix);
hetaY->Fill((decoder->getClusterElement(0,0)+decoder->getClusterElement(1,0))/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)+decoder->getClusterElement(0,-1))/br,iy+NR*(ix+1)); 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)+decoder->getClusterElement(1,0))/br,iy-1+NR*ix); hetaY->Fill(1.-(decoder->getClusterElement(0,0,cmsub)+decoder->getClusterElement(1,0,cmsub))/br,iy-1+NR*ix);
} }
} else if (tl>br && tl>bl && tl>tr) { } else if (tl>br && tl>bl && tl>tr) {
h2->Fill(tl, iy+NR*ix); h2->Fill(tl, iy+NR*ix);
if (tl>0) { if (tl>0) {
hetaX->Fill((decoder->getClusterElement(0,0)+decoder->getClusterElement(0,1))/tl,iy+NR*ix); hetaX->Fill((decoder->getClusterElement(0,0,cmsub)+decoder->getClusterElement(0,1,cmsub))/tl,iy+NR*ix);
hetaY->Fill((decoder->getClusterElement(0,0)+decoder->getClusterElement(-1,0))/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)+decoder->getClusterElement(0,1))/tl,iy+NR*(ix-1)); 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)+decoder->getClusterElement(-1,0))/tl,iy+1+NR*ix); hetaY->Fill(1.-(decoder->getClusterElement(0,0,cmsub)+decoder->getClusterElement(-1,0,cmsub))/tl,iy+1+NR*ix);
} }
} else if (tr>bl && tr>tl && tr>br) { } else if (tr>bl && tr>tl && tr>br) {
h2->Fill(tr, iy+NR*ix); h2->Fill(tr, iy+NR*ix);
if (tr>0) { if (tr>0) {
hetaX->Fill((decoder->getClusterElement(0,0)+decoder->getClusterElement(0,1))/tr,iy+NR*ix); hetaX->Fill((decoder->getClusterElement(0,0,cmsub)+decoder->getClusterElement(0,1,cmsub))/tr,iy+NR*ix);
hetaY->Fill((decoder->getClusterElement(0,0)+decoder->getClusterElement(1,0))/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)+decoder->getClusterElement(0,1))/tr,iy+NR*(ix+1)); 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)+decoder->getClusterElement(1,0))/tr,iy+1+NR*ix); hetaY->Fill(1.-(decoder->getClusterElement(0,0,cmsub)+decoder->getClusterElement(1,0,cmsub))/tr,iy+1+NR*ix);
} }
@ -207,17 +248,17 @@ Add here the function that you want to call: fill histos, make trees etc.
h1->Fill(decoder->getChannelShort(ix, iy), iy+NR*ix); h1->Fill(decoder->getChannelShort(ix, iy), iy+NR*ix);
} }
////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////
} }
cout << "=" ; cout << "=" ;
nf++; nf++;
} }
cout << endl; cout << endl;
if (filebin.is_open()) if (filebin.is_open())
filebin.close(); filebin.close();
else else
cout << "could not open file " << fname << endl; cout << "could not open file " << fname << endl;
} }
if (pedsub) if (pedsub)
tall->Write(tall->GetName(),TObject::kOverwrite); tall->Write(tall->GetName(),TObject::kOverwrite);
@ -226,5 +267,5 @@ Add here the function that you want to call: fill histos, make trees etc.
delete decoder; delete decoder;
cout << "Read " << nf << " frames" << endl; cout << "Read " << nf << " frames" << endl;
return hs; return hs;
} }