detectors/slsDetectorPackage
25
0
Fork 0
mirror of https://github.com/slsdetectorgroup/slsDetectorPackage.git synced 2025-04-26 16:20:03 +02:00
Code Issues Actions Packages Releases Activity

directory restuctured

Browse Source
This commit is contained in:
bergamaschi 2017-12-06 09:10:13 +01:00
parent 7f256c868b
commit 1bf7b4a6a1
78 changed files with 4212 additions and 8344 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

113
slsDetectorCalibration/IntMap.h
View File

@ -1,113 +0,0 @@
#ifndef INTMAP_H
#define INTMAP_H
#define N 60
#include <stdio.h>
#include <math.h>
class IntMap{
public:
IntMap(){};
~IntMap(){};
void Init(){
//lookup table of the output intensity for IR laser
//measurements performed by Dominic april 2014
//intensity[59]=intensity with 0 Optical Density
//intensity[0]=intensity with 5.9 Optical Density
intensity[0]=29;//5.9
intensity[1]=21;//5.8
intensity[2]=31;//5.7
intensity[3]=43;//5.6
intensity[4]=60;//5.5
intensity[5]=91;//5.4
intensity[6]=69;//5.3
intensity[7]=102;//5.2
intensity[8]=136;//5.1
intensity[9]=196;//5.0
intensity[10]=425;//4.9
intensity[11]=311;//4.8
intensity[12]=462;//4.7
intensity[13]=653;//4.6
intensity[14]=926;//4.5
intensity[15]=1423;//4.4
intensity[16]=1072;//4.3
intensity[17]=1592;//4.2
intensity[18]=2142;//4.1
intensity[19]=3085;//4.0
intensity[20]=729;//3.9
intensity[21]=533;//3.8
intensity[22]=793;//3.7
intensity[23]=1121;//3.6
intensity[24]=1588;//3.5
intensity[25]=2439;//3.4
intensity[26]=1842;//3.3
intensity[27]=2730;//3.2
intensity[28]=3663;//3.1
intensity[29]=5271;//3.0
intensity[30]=8102;//2.9
intensity[31]=5933;//2.8
intensity[32]=8789;//2.7
intensity[33]=12350;//2.6
intensity[34]=17358;//2.5
intensity[35]=26300;//2.4
intensity[36]=20029;//2.3
intensity[37]=29414;//2.2
intensity[38]=39202;//2.1
intensity[39]=55724;//2.0
intensity[40]=15697;//1.9
intensity[41]=11541;//1.8
intensity[42]=16976;//1.7
intensity[43]=23866;//1.6
intensity[44]=33478;//1.5
intensity[45]=50567;//1.4
intensity[46]=38552;//1.3
intensity[47]=56394;//1.2
intensity[48]=74897;//1.1
intensity[49]=106023;//1.0
intensity[50]=157384;//0.9
intensity[51]=117677;//0.8
intensity[52]=171101;//0.7
intensity[53]=236386;//0.6
intensity[54]=327248;//0.5
intensity[55]=492781;//0.4
intensity[56]=379641;//0.3
intensity[57]=546927;//0.2
intensity[58]=717203;//0.1
intensity[59]=1000000;//0.
return;
};
//_od is the total Optical Density
int getIntensity(float _od){
int _int(-1);
//these lines are to take into account rounding errors with floats
float hun_od = 100.*_od;
int Ihun_od = (int)round(hun_od);
float R_od =(float) Ihun_od/10.;
int I_od = (int)R_od;
if(I_od >-1 && I_od <60){
int index=59-I_od;
cerr<<index<<endl;
_int=intensity[index];
}else{
cerr<<"Optical density out of range!"<<endl;
}
return _int;
};
private:
int intensity[N];
};
#endif

194
slsDetectorCalibration/PlotGifs.C
View File

@ -1,194 +0,0 @@
#include "moench03ReadData.C"
/************************************************************************/
TH2F *readExactImage(char *fname, int iframe=0, int frperfile,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);
int framen(0);
moench03CtbData *decoder=new moench03CtbData();
char *buff=decoder->readNextFrame(filebin);
framen=decoder->getFrameNumber(buff);
int counter(0);
while(framen<iframe && counter<frperfile){
buff=decoder->readNextFrame(filebin);
framen=decoder->getFrameNumber(buff);
cerr<<"...";
if(framen%1000==0) cerr<<framen;
counter++;
}
if(counter<frperfile){
h2->SetName(Form("frame_%d",framen));
h2->SetTitle(Form("frame_%d",framen));
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);
}else{
cerr<<"frame number not found"<<endl;
}
return h2;
}
/************************************************************/
//fnamein filename: ..._f0_%d.raw
//runmin, runmax to calculate pedestals
//framen0, framen1 first and last frame you want to add to the gif
//frperfile number of frames per file
void PlotRawFrameGif(char * fnamein, int runmin, int runmax, int framen0,int framen1, int frperfile){
cerr<<"/***********************************/"<<endl;
cerr<<"calculating pedestals"<<endl;
TH2F * hp = calcPedestal(fnamein,runmin,runmax);
int filen = (int)(framen0/frperfile);
char fname[1000];
sprintf(fname,fnamein,filen);
cerr<<"/***********************************/"<<endl;
cerr<<"retrieving frame from"<<fname<<endl;
int fileframe0 = framen0%frperfile;
int fileframe1 = framen1%frperfile;
TImage * img = NULL;
TH2F * hf = NULL;
TCanvas * c1 =NULL;
for(int fileframe=fileframe0; fileframe<fileframe1; fileframe++){
hf=readExactImage(fname, fileframe,frperfile,hp);
delete img;
delete c1;
c1 = new TCanvas("c1","",800,600);
c1->cd();
hf->SetTitle(Form("Frame_%d",fileframe+framen0));
hf->SetName(Form("Frame_%d",fileframe+framen0));
hf->GetXaxis()->SetRangeUser(0,50);
hf->GetXaxis()->SetTitle("Column");
hf->GetYaxis()->SetRangeUser(240,290);
hf->GetYaxis()->SetTitle("Row");
hf->GetZaxis()->SetRangeUser(-50.,1300.);
hf->SetStats(kFALSE);
// c1->SetLogz();
hf->Draw("colz");
c1->Print(Form("/afs/psi/project/mythen/Marco/Pics/Fe_Raw3_%d.png",fileframe));
img = TImage::Open(Form("/afs/psi/project/mythen/Marco/Pics/Fe_Raw3_%d.png",fileframe));
if(fileframe<fileframe1-1){
img->WriteImage(Form("/afs/psi/project/mythen/Marco/Pics/Fe_Raw3_%d.gif+200",fileframe1-fileframe0));
}else{
img->WriteImage(Form("/afs/psi/project/mythen/Marco/Pics/Fe_Raw3_%d.gif++200++",fileframe1-fileframe0));
}
}
return;
}
/*********************************************************************/
TProfile2D * GetHitMap(TTree * treein, int framen,double zlow, double zup){
TProfile2D* map = new TProfile2D("map","",400,-0.5,399.5,400,-0.5,399.5,zlow,zup);
int x, y, iFrame;
double data[9];
TBranch * b_x = (TBranch*)treein->GetBranch("x");
TBranch * b_y = (TBranch*)treein->GetBranch("y");
TBranch * b_data = (TBranch*)treein->GetBranch("data");
TBranch * b_iFrame = (TBranch*)treein->GetBranch("iFrame");
b_x->SetAddress(&x);
b_y->SetAddress(&y);
b_data->SetAddress(data);
b_iFrame->SetAddress(&iFrame);
Int_t nEnt=treein->GetEntries();
for(Int_t i=0; i<nEnt; i++){
b_iFrame->GetEntry(i);
if(iFrame==framen){
b_x->GetEntry(i);
b_y->GetEntry(i);
b_data->GetEntry(i);
map->SetBinEntries(map->FindFixBin(x-1,y-1),1);
map->SetBinEntries(map->FindFixBin(x,y-1),1);
map->SetBinEntries(map->FindFixBin(x+1,y-1),1);
map->SetBinEntries(map->FindFixBin(x-1,y),1);
map->SetBinEntries(map->FindFixBin(x,y),1);
map->SetBinEntries(map->FindFixBin(x+1,y),1);
map->SetBinEntries(map->FindFixBin(x-1,y+1),1);
map->SetBinEntries(map->FindFixBin(x,y+1),1);
map->SetBinEntries(map->FindFixBin(x+1,y+1),1);
map->SetBinContent(map->FindFixBin(x-1,y-1),data[0]);
map->SetBinContent(map->FindFixBin(x,y-1),data[1]);
map->SetBinContent(map->FindFixBin(x+1,y-1),data[2]);
map->SetBinContent(map->FindFixBin(x-1,y),data[3]);
map->SetBinContent(map->FindFixBin(x,y),data[4]);
map->SetBinContent(map->FindFixBin(x+1,y),data[5]);
map->SetBinContent(map->FindFixBin(x-1,y+1),data[6]);
map->SetBinContent(map->FindFixBin(x,y+1),data[7]);
map->SetBinContent(map->FindFixBin(x+1,y+1),data[8]);
}
}
return map;
}
/*********************************************************/
/** creates an infinitely looping gif from clustered data
**/
void PlotClusterHitMapGif(std::string filename, std::string treename,int framen0,int framen1){
TFile * fin = new TFile(filename.c_str(),"read");
TTree * treein = (TTree*)fin->Get(treename.c_str());
TCanvas * c1 = NULL;
TProfile2D* hmap = NULL;
TImage * img = NULL;
for(int framen=framen0; framen<framen1; framen++){
delete c1;
c1 = new TCanvas("c1","",800,600);
c1->cd();
// c1->SetLogz();
hmap=GetHitMap(treein,framen,-50.,1300.);
hmap->SetName(Form("Frame_%d",framen));
hmap->SetTitle(Form("Frame_%d",framen));
hmap->GetXaxis()->SetRangeUser(0,50);
hmap->GetXaxis()->SetTitle("Column");
hmap->GetYaxis()->SetRangeUser(240,290);
hmap->GetYaxis()->SetTitle("Row");
hmap->GetZaxis()->SetRangeUser(-50.,1300.);
hmap->SetStats(kFALSE);
hmap->Draw("colz");
c1->Print(Form("/afs/psi/project/mythen/Marco/Pics/Fe_Cluster3_%d.png",framen));
img = TImage::Open(Form("/afs/psi/project/mythen/Marco/Pics/Fe_Cluster3_%d.png",framen));
if(framen<framen1-1){
img->WriteImage(Form("/afs/psi/project/mythen/Marco/Pics/Fe_Cluster3_%d.gif+200",framen1-framen0));
}else{
img->WriteImage(Form("/afs/psi/project/mythen/Marco/Pics/Fe_Cluster3_%d.gif++200++",framen1-framen0));
}
}
return;
}

365
slsDetectorCalibration/analogDetector.h
View File

@ -1,16 +1,21 @@
#ifndef ANALOGDETECTOR_H
#define ANALOGDETECTOR_H
//#include <mutex>
#include <pthread.h>
#include "slsDetectorData.h"
#include "pedestalSubtraction.h"
#include "commonModeSubtraction.h"
#include "tiffIO.h"
using namespace std;
#ifndef FRAMEMODE_DEF
#define FRAMEMODE_DEF
/**
enum to define the flags of the data set, which are needed to seect the type of processing it should undergo
enum to define the flags of the data set, which are needed to seect the type of processing it should undergo: frame, pedestal, flat
*/
enum frameMode { eFrame, ePedestal, eFlat };
#endif
@ -53,10 +58,17 @@ template <class dataType> class analogDetector {
}
}
image=new int[nx*ny];
xmin=0;
xmax=nx;
ymin=0;
ymax=ny;
fMode=ePedestal;
thr=0;
myFile=NULL;
fm=new pthread_mutex_t ;
};
/**
destructor. Deletes the cluster structure and the pdestalSubtraction array
destructor. Deletes the pdestalSubtraction array and the image
*/
virtual ~analogDetector() {for (int i=0; i<ny; i++) delete [] stat[i]; delete [] stat; delete [] image;};
@ -74,6 +86,16 @@ template <class dataType> class analogDetector {
gmap=orig->gmap;
cmSub=orig->cmSub;
id=orig->id;
xmin=orig->xmin;
xmax=orig->xmax;
ymin=orig->ymin;
ymax=orig->ymax;
thr=orig->thr;
// nSigma=orig->nSigma;
fMode=orig->fMode;
myFile=orig->myFile;
fm=orig->fm;
stat=new pedestalSubtraction*[ny];
for (int i=0; i<ny; i++) {
@ -92,8 +114,10 @@ template <class dataType> class analogDetector {
}
/**
clone. Must be virtual!
\returns a clone of the original analog detector
*/
virtual analogDetector *Clone() {
return new analogDetector(this);
@ -165,12 +189,13 @@ template <class dataType> class analogDetector {
return NULL;
}
/**
wrties a 32 bit tiff file of the size of the detector and contaning the gain map value, if any. If file doesn'e exist or gainmap is undefined, does not do anything.
writes a 32 bit tiff file of the size of the detector and contaning the gain map value, if any. If file doesn'e exist or gainmap is undefined, does not do anything.
\param imgname complete name of the file to be written
\returns NULL
\returns NULL if file writing didn't succeed, else a pointer
*/
void *writeGainMap(const char * imgname) {
float *gm=NULL;
void *ret;
if (gmap) {
gm=new float[nx*ny];
for (int ix=0; ix<nx; ix++) {
@ -178,19 +203,27 @@ template <class dataType> class analogDetector {
gm[iy*nx+ix]=gmap[iy*nx+ix];
}
}
WriteToTiff(gm, imgname, ny, nx);
ret=WriteToTiff(gm, imgname, ny, nx);
delete [] gm;
}
return NULL;
return ret;
}
/** resets the pedestalSubtraction array, the commonModeSubtraction and the image data*/
/** resets the pedestalSubtraction array and the commonModeSubtraction */
void newDataSet(){iframe=-1; for (int iy=0; iy<ny; iy++) for (int ix=0; ix<nx; ix++) stat[iy][ix].Clear(); if (cmSub) cmSub->Clear(); };
virtual void newDataSet(){
iframe=-1;
for (int iy=0; iy<ny; iy++)
for (int ix=0; ix<nx; ix++) {
stat[iy][ix].Clear();
image[iy*nx+ix]=0;
}
if (cmSub) cmSub->Clear();
};
/** resets the eventMask to undefined and the commonModeSubtraction */
void newFrame(){iframe++; if (cmSub) cmSub->newFrame();};
/** resets the commonModeSubtraction and increases the frame index */
virtual void newFrame(){iframe++; if (cmSub) cmSub->newFrame();};
/** sets the commonModeSubtraction algorithm to be used
@ -198,6 +231,11 @@ template <class dataType> class analogDetector {
\returns pointer to the actual common mode subtraction algorithm
*/
commonModeSubtraction *setCommonModeSubtraction(commonModeSubtraction *cm) {cmSub=cm; return cmSub;};
/**
gets the commonModeSubtraction algorithm to be used
\returns pointer to the actual common mode subtraction algorithm
*/
commonModeSubtraction *getCommonModeSubtraction() {return cmSub;};
/**
@ -213,6 +251,7 @@ template <class dataType> class analogDetector {
\param val value to be added
\param ix pixel x coordinate
\param iy pixel y coordinate
\param cm 1 adds the value to common mod, 0 skips it. Defaults to 0. - not properly implemented
*/
virtual void addToPedestal(double val, int ix, int iy=0){
if (ix>=0 && ix<nx && iy>=0 && iy<ny) {
@ -310,8 +349,8 @@ template <class dataType> class analogDetector {
*/
virtual void setPedestal(double *ped, double *rms=NULL, int m=-1){
double rr=0;
for (int ix=0; ix<nx; ix++) {
for (int iy=0; iy<ny; iy++) {
for (int ix=xmin; ix<xmax; ix++) {
for (int iy=ymin; iy<ymax; iy++) {
if (rms) rr=rms[iy*nx+ix];
stat[iy][ix].setPedestal(ped[iy*nx+ix],rr, m);
};
@ -331,9 +370,16 @@ template <class dataType> class analogDetector {
*/
virtual void setPedestalRMS(int ix, int iy, double rms=0){if (ix>=0 && ix<nx && iy>=0 && iy<ny) stat[iy][ix].setPedestalRMS(rms);};
/**
sets pedestal rms for all pixels
\param rms pointer to array of pedestal rms
*/
virtual void setPedestalRMS(double *rms){
for (int ix=0; ix<nx; ix++) {
for (int iy=0; iy<ny; iy++) {
for (int ix=xmin; ix<xmax; ix++) {
for (int iy=ymin; iy<ymax; iy++) {
stat[iy][ix].setPedestalRMS(rms[iy*nx+ix]);
};
};
@ -343,27 +389,33 @@ template <class dataType> class analogDetector {
/**
sets pedestal rms
\param ix pixel x coordinate
\param iy pixel y coordinate
\param rms value to set
write 32bit tiff file with detector image data
\param imgname file name to be written
\returns NULL if file writing didn't succed, otherwise a pointer
*/
virtual void *writeImage(const char * imgname) {
float *gm=NULL;
void *ret;
gm=new float[nx*ny];
for (int ix=0; ix<nx; ix++) {
for (int iy=0; iy<ny; iy++) {
gm[iy*nx+ix]=image[iy*nx+ix];
}
}
WriteToTiff(gm, imgname, ny, nx);
ret=WriteToTiff(gm, imgname, ny, nx);
delete [] gm;
return NULL;
return ret;
}
/**
write 32bit tiff file containing the pedestals
\param imgname file name to be written
\returns NULL if file writing didn't succed, otherwise a pointer
*/
virtual void *writePedestals(const char * imgname) {
float *gm=NULL;
void *ret;
gm=new float[nx*ny];
for (int ix=0; ix<nx; ix++) {
for (int iy=0; iy<ny; iy++) {
@ -373,11 +425,16 @@ template <class dataType> class analogDetector {
gm[iy*nx+ix]=stat[iy][ix].getPedestal();
}
}
WriteToTiff(gm, imgname, ny, nx);
ret=WriteToTiff(gm, imgname, ny, nx);
delete [] gm;
return NULL;
return ret;
}
/**
read 32bit tiff file containing the pedestals
\param imgname file name to be read
\returns 0 if file reading didn't succed, otherwise 1
*/
int readPedestals(const char * imgname) {
uint32 nnx, nny;
float *gm=ReadFromTiff( imgname, nny, nnx);
@ -398,6 +455,11 @@ template <class dataType> class analogDetector {
return NULL;
}
/**
read 32bit tiff file containing the image data
\param imgname file name to be read
\returns 0 if file reading didn't succed, otherwise 1
*/
int readImage(const char * imgname) {
uint32 nnx, nny;
float *gm=ReadFromTiff( imgname, nny, nnx);
@ -419,21 +481,36 @@ template <class dataType> class analogDetector {
}
/**
returns pointer to image data
\returns pointer to image data
*/
virtual int *getImage(){return image;};
/**
write 32bit tiff file containing the pedestals RMS
\param imgname file name to be written
\returns NULL if file writing didn't succed, otherwise a pointer
*/
void *writePedestalRMS(const char * imgname) {
float *gm=NULL;
void *ret;
gm=new float[nx*ny];
for (int ix=0; ix<nx; ix++) {
for (int iy=0; iy<ny; iy++) {
gm[iy*nx+ix]=stat[iy][ix].getPedestalRMS();
}
}
WriteToTiff(gm, imgname, ny, nx);
ret=WriteToTiff(gm, imgname, ny, nx);
delete [] gm;
return NULL;
return ret;
}
/**
read 32bit tiff file containing the pedestals RMS
\param imgname file name to be read
\returns 0 if file reading didn't succed, otherwise 1
*/
int readPedestalRMS(const char * imgname) {
uint32 nnx, nny;
float *gm=ReadFromTiff( imgname, nny, nnx);
@ -456,14 +533,19 @@ template <class dataType> class analogDetector {
/**
Adds all the data for each pixels in the selected region of interest to the pedestal
\param data pointer to the data
*/
virtual void addToPedestal(char *data) {
newFrame();
for (int ix=0; ix<nx; ix++) {
for (int iy=0; iy<ny; iy++) {
for (int ix=xmin; ix<xmax; ix++) {
for (int iy=ymin; iy<ymax; iy++) {
addToPedestal(data,ix,iy);
@ -474,6 +556,58 @@ template <class dataType> class analogDetector {
return ;
};
/**
Sets region of interest in which data should be processed
\param xmi minimum x. if -1 or out of range remains unchanged
\param xma maximum x. if -1 or out of range remains unchanged
\param ymi minimum y. if -1 or out of range remains unchanged
\param yma maximum y. if -1 or out of range remains unchanged
*/
void setROI(int xmi=-1, int xma=-1, int ymi=-1, int yma=-1) {
if (xmi>=0 && xmi<=nx) xmin=xmi;
if (xma>=0 && xma<=nx) xmax=xma;
if (xmax<xmin) {
xmi=xmin;
xmin=xmax;
xmax=xmi;
}
if (ymi>=0 && ymi<=ny) ymin=ymi;
if (yma>=0 && yma<=ny) ymax=yma;
if (ymax<ymin) {
ymi=ymin;
ymin=ymax;
ymax=ymi;
}
};
/**
Gets region of interest in which data are processed
\param xmi reference to minimum x.
\param xma reference to maximum x.
\param ymi reference to minimum y.
\param yma reference to maximum y.
*/
void getROI(int &xmi, int &xma, int &ymi, int &yma) {xmi=xmin; xma=xmax; ymi=ymin; yma=ymax;};
/**
Adds all the data for the selected pixel to the pedestal
\param data pointer to the data
\param ix pixel x coordinate
\param iy pixel y coordinate
*/
virtual void addToPedestal(char *data, int ix, int iy=0) {
@ -492,6 +626,13 @@ template <class dataType> class analogDetector {
};
/**
Subtracts pedestal from the data array in the region of interest
\param data pointer to the data
\param val pointer where the pedestal subtracted data should be added. If NULL, the internal image is used
\returns pointer to the pedestal subtracted data
*/
virtual double *subtractPedestal(char *data, double *val=NULL) {
@ -500,8 +641,8 @@ template <class dataType> class analogDetector {
if (val==NULL)
val=new double[nx*ny];
for (int ix=0; ix<nx; ix++) {
for (int iy=0; iy<ny; iy++) {
for (int ix=xmin; ix<xmax; ix++) {
for (int iy=ymin; iy<ymax; iy++) {
val[iy*nx+ix]+=subtractPedestal(data, ix, iy);
}
}
@ -511,6 +652,15 @@ template <class dataType> class analogDetector {
/**
Subtracts pedestal from the data for a selected pixel
\param data pointer to the data
\param ix pixel x coordinate
\param iy pixel y coordinate
\returns pedestal subtracted value
*/
virtual double subtractPedestal(char *data, int ix, int iy=0) {
double g=1.;
@ -528,52 +678,73 @@ template <class dataType> class analogDetector {
};
/**
sets threshold value for conversion into number of photons
\param t threshold to be set
\returns threshold value
*/
double setThreshold(double t){thr=t; return thr;};
virtual int getNPhotons(char *data, int ix, int iy=0, int thr=-1) {
double g=1.;
int v;
/**
gets threshold value for conversion into number of photons
\returns threshold value
*/
double getThreshold(){return thr;};
/**
converts the data into number of photons for the selected pixel
\param data pointer to the data
\param ix pixel x coordinate
\param iy pixel y coordinate
\returns converted number of photons. If no threshold is set, returns gain converted pedestal subtracted data.
*/
virtual int getNPhotons(char *data, int ix, int iy=0) {
int nph=0;
double v;
if (ix>=0 && ix<nx && iy>=0 && iy<ny) {
if (gmap) {
g=gmap[iy*nx+ix];
if (g==0) g=-1.;
}
if (thr<=0) thr=-1*thr*getPedestalRMS(ix,iy)/g;
/* if (det) */
/* return (dataSign*det->getValue(data, ix, iy)-getPedestal(ix,iy))/g/thr; */
/* else { */
//if (ix==30 && iy==50 )
v=subtractPedestal(data,ix,iy)/g/thr;
if (v>0)
v=subtractPedestal(data,ix,iy);
if (thr>0) {
v+=0.5*thr;
nph=v/thr;
return nph;
} else
return v;
//cout << v << endl;
return 0;
// }
}
return 0;
};
int *getNPhotons(char *data, double thr=-1, int *nph=NULL) {
/**
converts the data into number of photons for all pixels
\param data pointer to the data
\param nph pointer where the photons should added. If NULL,the internal image is used
\returns pointer to array containing the number of photons
*/
int *getNPhotons(char *data, int *nph=NULL) {
double val;
if (nph==NULL)
nph=image;//image=new int[nx*ny];
double tthr=thr;
nph=image;
newFrame();
for (int ix=0; ix<nx; ix++) {
for (int iy=0; iy<ny; iy++) {
nph[iy*nx+ix]+=getNPhotons(data, ix, iy, thr);
for (int ix=xmin; ix<xmax; ix++) {
for (int iy=ymin; iy<ymax; iy++) {
nph[iy*nx+ix]+=getNPhotons(data, ix, iy);
}
}
return nph;
}
virtual void clearImage(){ for (int ix=0; ix<nx; ix++) {
for (int iy=0; iy<ny; iy++) { image[iy*nx+ix]=0;}}};
/**
clears the image array
*/
virtual void clearImage(){
for (int ix=0; ix<nx; ix++) {
for (int iy=0; iy<ny; iy++) {
image[iy*nx+ix]=0;
}
}
};
/** sets/gets number of samples for moving average pedestal calculation
\param i number of samples to be set (0 or negative gets)
@ -588,6 +759,9 @@ template <class dataType> class analogDetector {
return stat[0][0].SetNPedestals();
};
/** gets number of samples for moving average pedestal calculation
\returns actual number of samples
*/
int GetNPedestals(int ix, int iy) {
if (ix>=0 && ix<nx && iy>=0 && iy<ny)
return stat[iy][ix].GetNPedestals();
@ -596,29 +770,74 @@ template <class dataType> class analogDetector {
};
double getROI(char *data, int xmin, int xmax, int ymin=0, int ymax=1) {
double val=0;
for (int ix=xmin; ix<xmax; ix++)
for (int iy=ymin; iy<ymax; iy++)
/** calculates the sum of photons in the specified region of interest.
\param data pointer to the data
\param xmi minimum x for the calculation. If -1 the minimum x of the predefined region of interest is used
\param xma maximum x for the calculation. If -1 the maximum x of the predefined region of interest is used
\param ymi minimum y for the calculation. If -1 the minimum y of the predefined region of interest is used
\param yma maximum y for the calculation. If -1 the maximum y of the predefined region of interest is used
\returns total number of photons in
*/
virtual int getTotalNumberOfPhotons(char *data, int xmi=-1, int xma=-1, int ymi=-1, int yma=-1) {
int val=0;
if (xmi<0) xmi=xmin;
if (xma<0) xma=xmax;
if (ymi<0) ymi=ymin;
if (yma<0) yma=ymax;
for (int ix=xmi; ix<xma; ix++)
for (int iy=ymi; iy<yma; iy++)
if (ix>=0 && ix<nx && iy>=0 && iy<ny)
val+=subtractPedestal(data, ix, iy);
val+=getNPhotons(data, ix, iy);
return val;
};
virtual void processData(char *data, frameMode i=eFrame, int *val=NULL) {
switch(i) {
/**
calculates the image converted into number of photons. If the frame mode is pedestal, it also it to the pdedestal subtraction.
\param data pointer to the data to be processed
\param val pointer of the data to be added to. If NULL, the internal image will be used
\param pointer to the processed data
\returns
*/
virtual void processData(char *data,int *val=NULL) {
switch(fMode) {
case ePedestal:
addToPedestal(data);
break;
default:
getNPhotons(data,-1,val);
getNPhotons(data,val);
}
};
virtual char *getInterpolation(){return NULL;};
/** sets the current frame mode for the detector
\param f frame mode to be set
\returns current frame mode
*/
frameMode setFrameMode(frameMode f) {fMode=f; return fMode;};
/** gets the current frame mode for the detector
\returns current frame mode
*/
frameMode getFrameMode() {return fMode;};
/** sets file pointer where to write the clusters to
\param f file pointer
\returns current file pointer
*/
FILE *setFilePointer(FILE *f){myFile=f; return myFile;};
/** gets file pointer where to write the clusters to
\returns current file pointer
*/
FILE *getFilePointer(){return myFile;};
void setMutex(pthread_mutex_t *m){fm=m;};
protected:
slsDetectorData<dataType> *det; /**< slsDetectorData to be used */
@ -631,7 +850,17 @@ template <class dataType> class analogDetector {
double *gmap;
int *image;
int id;
//int xmin, xmax, ymin, ymax;
//int xmin, xmax, ymin, ymax; int xmin; /**< minimum x of the region of interest */
int xmin; /**< minimum x of the region of interest */
int xmax; /**< maximum x of the region of interest */
int ymin;/**< minimum y of the region of interest */
int ymax;/**< maximum y of the region of interest */
double thr; /**< threshold to be used for conversion into number of photons */
// int nSigma; /**< number of sigma to be used for conversion into number of photons if threshold is undefined */
frameMode fMode; /**< current detector frame mode */
FILE *myFile; /**< file pointer to write to */
pthread_mutex_t *fm;
};
#endif

0
slsDetectorCalibration/Mythen3_01_jctbData.h → slsDetectorCalibration/dataStructures/Mythen3_01_jctbData.h
View File

0
slsDetectorCalibration/adcSar2_jctbData.h → slsDetectorCalibration/dataStructures/adcSar2_jctbData.h
View File

0
slsDetectorCalibration/chiptestBoardData.h → slsDetectorCalibration/dataStructures/chiptestBoardData.h
View File

0
slsDetectorCalibration/eigerHalfModuleData.h → slsDetectorCalibration/dataStructures/eigerHalfModuleData.h
View File

0
slsDetectorCalibration/gotthardModuleData.h → slsDetectorCalibration/dataStructures/gotthardModuleData.h
View File

0
slsDetectorCalibration/gotthardShortModuleData.h → slsDetectorCalibration/dataStructures/gotthardShortModuleData.h
View File

0
slsDetectorCalibration/jungfrau02Data.h → slsDetectorCalibration/dataStructures/jungfrau02Data.h
View File

0
slsDetectorCalibration/jungfrau10ModuleData.h → slsDetectorCalibration/dataStructures/jungfrau10ModuleData.h
View File

0
slsDetectorCalibration/moench02Ctb10GbData.h → slsDetectorCalibration/dataStructures/moench02Ctb10GbData.h
View File

0
slsDetectorCalibration/moench02CtbData.h → slsDetectorCalibration/dataStructures/moench02CtbData.h
View File

0
slsDetectorCalibration/moench02ModuleData.h → slsDetectorCalibration/dataStructures/moench02ModuleData.h
View File

0
slsDetectorCalibration/moench03Ctb10GbData.h → slsDetectorCalibration/dataStructures/moench03Ctb10GbData.h
View File

0
slsDetectorCalibration/moench03Ctb10GbT1Data.h → slsDetectorCalibration/dataStructures/moench03Ctb10GbT1Data.h
View File

0
slsDetectorCalibration/moench03CtbData.h → slsDetectorCalibration/dataStructures/moench03CtbData.h
View File

0
slsDetectorCalibration/moench03T1CtbData.h → slsDetectorCalibration/dataStructures/moench03T1CtbData.h
View File

285
slsDetectorCalibration/dataStructures/moench03T1ReceiverData.h Normal file
View File

@ -0,0 +1,285 @@
#ifndef MOENCH03T1ZMQDATA_H
#define MOENCH03T1ZMQDATA_H
#include "slsDetectorData.h"
/**
@short structure for a Detector Packet or Image Header
@li frameNumber is the frame number
@li expLength is the subframe number (32 bit eiger) or real time exposure time in 100ns (others)
@li packetNumber is the packet number
@li bunchId is the bunch id from beamline
@li timestamp is the time stamp with 10 MHz clock
@li modId is the unique module id (unique even for left, right, top, bottom)
@li xCoord is the x coordinate in the complete detector system
@li yCoord is the y coordinate in the complete detector system
@li zCoord is the z coordinate in the complete detector system
@li debug is for debugging purposes
@li roundRNumber is the round robin set number
@li detType is the detector type see :: detectorType
@li version is the version number of this structure format
*/
typedef struct {
uint64_t frameNumber; /**< is the frame number */
uint32_t expLength; /**< is the subframe number (32 bit eiger) or real time exposure time in 100ns (others) */
uint32_t packetNumber; /**< is the packet number */
uint64_t bunchId; /**< is the bunch id from beamline */
uint64_t timestamp; /**< is the time stamp with 10 MHz clock */
uint16_t modId; /**< is the unique module id (unique even for left, right, top, bottom) */
uint16_t xCoord; /**< is the x coordinate in the complete detector system */
uint16_t yCoord; /**< is the y coordinate in the complete detector system */
uint16_t zCoord; /**< is the z coordinate in the complete detector system */
uint32_t debug; /**< is for debugging purposes */
uint16_t roundRNumber; /**< is the round robin set number */
uint8_t detType; /**< is the detector type see :: detectorType */
uint8_t version; /**< is the version number of this structure format */
} sls_detector_header;
class moench03T1ReceiverData : public slsDetectorData<uint16_t> {
private:
int iframe;
int nadc;
int sc_width;
int sc_height;
const int nPackets; /**<number of UDP packets constituting one frame */
const int packetSize; /**< size of a udp packet */
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
*/
moench03T1ReceiverData(int npackets=40, int ps=8192): slsDetectorData<uint16_t>(400, 400, ps*npackets+sizeof(sls_detector_header)), packetSize(ps), nPackets(npackets) {
int nadc=32;
int sc_width=25;
int sc_height=200;
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,\
0,25,50,75,0,25,50,75};
int row, col;
int isample;
int iadc;
int ix, iy;
// int npackets=40;
int i;
int adc4(0);
for (int ip=0; ip<npackets; ip++) {
for (int is=0; is<128; is++) {
for (iadc=0; iadc<nadc; iadc++) {
i=128*ip+is;
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) {
row=199-i/sc_width;
} else {
row=200+i/sc_width;
}
dataMap[row][col]=sizeof(sls_detector_header)+(nadc*i+iadc)*2;//+16*(ip+1);
if (dataMap[row][col]<0 || dataMap[row][col]>=8192*40)
cout << "Error: pointer " << dataMap[row][col] << " out of range "<< endl;
}
}
}
}
int ipacket;
int ibyte;
int ii=0;
for (ibyte=0; ibyte<sizeof(sls_detector_header)/2; ibyte++){
xmap[ibyte]=-1;
ymap[ibyte]=-1;
}
int off=sizeof(sls_detector_header)/2;
for (ipacket=0; ipacket<npackets; ipacket++) {
for (ibyte=0; ibyte< 8192/2; ibyte++) {
i=ipacket*8208/2+ibyte;
isample=ii/nadc;
iadc=ii%nadc;
adc4 = (int)iadc/4;
ix=isample%sc_width;
iy=isample/sc_width;
if (adc4%2==0) {
xmap[i+off]=adc_nr[iadc]+ix;
ymap[i+off]=ny/2-1-iy;
} else {
xmap[i+off]=adc_nr[iadc]+ix;
ymap[i+off]=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
*/
/* class jfrau_packet_header_t { */
/* public: */
/* unsigned char reserved[4]; */
/* unsigned char packetNumber[1]; */
/* unsigned char frameNumber[3]; */
/* unsigned char bunchid[8]; */
/* }; */
int getFrameNumber(char *buff){return ((sls_detector_header*)buff)->frameNumber;};//*((int*)(buff+5))&0xffffff;};
/**
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){((sls_detector_header*)buff)->packetNumber;}//((*(((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 ff=-1, np=-1;
return readNextFrame(filebin, ff, np);
};
virtual char *readNextFrame(ifstream &filebin, int &ff) {
int np=-1;
return readNextFrame(filebin, ff, np);
};
virtual char *readNextFrame(ifstream &filebin, int& ff, int &np) {
char *data=new char[dataSize];
char *d=readNextFrame(filebin, ff, np, data);
if (d==NULL) {delete [] data; data=NULL;}
return data;
}
virtual char *readNextFrame(ifstream &filebin, int& ff, int &np, char *data) {
char *retval=0;
int nd;
int fnum = -1;
np=0;
int pn;
// cout << dataSize << endl;
if (ff>=0)
fnum=ff;
if (filebin.is_open()) {
if (filebin.read(data, dataSize) ){
ff=getFrameNumber(data);
np=getPacketNumber(data);
return data;
}
}
return NULL;
};
/**
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){
if (dsize<dataSize) ndata=dsize;
else ndata=dataSize;
return data;
}
int getPacketNumber(int x, int y) {return dataMap[y][x]/packetSize;};
};
#endif

0
slsDetectorCalibration/moench03T1ZmqData.h → slsDetectorCalibration/dataStructures/moench03T1ZmqData.h
View File

268
slsDetectorCalibration/dataStructures/moench03T1ZmqDataNew.h Normal file
View File

@ -0,0 +1,268 @@
#ifndef MOENCH03T1ZMQDATANEW_H
#define MOENCH03T1ZMQDATANEW_H
#include "slsDetectorData.h"
class moench03T1ZmqDataNew : public slsDetectorData<uint16_t> {
private:
int iframe;
int nadc;
int sc_width;
int sc_height;
const int nSamples;
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
*/
moench03T1ZmqDataNew(int ns=5000): slsDetectorData<uint16_t>(400, 400, ns*32*2), nSamples(ns) {
int nadc=32;
int sc_width=25;
int sc_height=200;
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,\
0,25,50,75,0,25,50,75};
int row, col;
int isample;
int iadc;
int ix, iy;
int npackets=40;
int i;
int adc4(0);
for (int ip=0; ip<npackets; ip++) {
for (int is=0; is<128; is++) {
for (iadc=0; iadc<nadc; iadc++) {
i=128*ip+is;
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) {
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]>=nSamples*2*32)
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< 8192/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;
if (isample<nSamples) {
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++;
// }
}
}
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
*/
/* class jfrau_packet_header_t { */
/* public: */
/* unsigned char reserved[4]; */
/* unsigned char packetNumber[1]; */
/* unsigned char frameNumber[3]; */
/* unsigned char bunchid[8]; */
/* }; */
int getFrameNumber(char *buff){return iframe;};//*((int*)(buff+5))&0xffffff;};
/**
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 0;}//((*(((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 ff=-1, np=-1;
return readNextFrame(filebin, ff, np);
};
virtual char *readNextFrame(ifstream &filebin, int &ff) {
int np=-1;
return readNextFrame(filebin, ff, np);
};
virtual char *readNextFrame(ifstream &filebin, int& ff, int &np) {
char *data=new char[32*2*nSamples];
char *d=readNextFrame(filebin, ff, np, data);
if (d==NULL) {delete [] data; data=NULL;}
return data;
}
virtual char *readNextFrame(ifstream &filebin, int& ff, int &np, char *data) {
char *retval=0;
int nd;
int fnum = -1;
np=0;
int pn;
if (ff>=0)
fnum=ff;
if (filebin.is_open()) {
if (filebin.read(data, 32*2*nSamples) ){
iframe++;
ff=iframe;
return data;
}
}
return NULL;
};
/**
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){
if (dsize<32*2*nSamples) ndata=dsize;
else ndata=32*2*nSamples;
return data;
}
int getPacketNumber(int x, int y) {return 0;};
};
#endif

0
slsDetectorCalibration/moench03TCtb10GbData.h → slsDetectorCalibration/dataStructures/moench03TCtb10GbData.h
View File

0
slsDetectorCalibration/moench03TCtbData.h → slsDetectorCalibration/dataStructures/moench03TCtbData.h
View File

31
slsDetectorCalibration/demoCreateTree.C
View File

@ -1,31 +0,0 @@
{
//.L moenchReadData.C
TFile *fout;
THStack *hs2N;
fout=new TFile("/scratch/outfile.root","RECREATE");
hs2N=moenchReadData("/data/moench_xbox_20140113/MoTarget_45kV_0_8mA_12us_120V_cds_g4_f00000%04d000_0.raw","dum",0,20,1500,-500,2500,1,0.,1,159,1,159, 0,1);
hs2N->SetName("cds_g4");
hs2N->SetTitle("cds_g4");
(TH2F*)(hs2N->GetHists()->At(0))->Write();
(TH2F*)(hs2N->GetHists()->At(1))->Write();
(TH2F*)(hs2N->GetHists()->At(2))->Write();
(TH2F*)(hs2N->GetHists()->At(3))->Write();
(TH2F*)(hs2N->GetHists()->At(4))->Write();
fout->Close();
}

457
slsDetectorCalibration/etaVEL/imageMacro.C
View File

@ -1,457 +0,0 @@
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=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");
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++) {
for (int iby=20*25; iby<hff->GetNbinsY()-20*25; iby++) {
hpixel->Fill((ibx-12)%25,(iby-12)%25,hff->GetBinContent(ibx+1,iby+1));
}
}
hpixel->Scale(1./hpixel->GetBinContent(13,13));
// new TCanvas();
// hpixel->Draw("colz");
// TH2D *hraw=(TH2D*)hff->Clone("hraw");
TH2D *hpix=(TH2D*)hff->Clone("hpix");
for (int ibx=0; ibx<hff->GetNbinsX(); ibx++) {
for (int iby=0; iby<hff->GetNbinsY(); iby++) {
hpix->SetBinContent(ibx+1,iby+1,hpixel->GetBinContent(hpixel->GetXaxis()->FindBin((ibx-12)%25),hpixel->GetXaxis()->FindBin((iby-12)%25)));
}
}
// return hpix;
//}
//void imageMacro(char *name,TH2D *hpix=NULL){
// hff->Divide(hpix);
// new TCanvas();
// hff->Draw("colz");
// if (hpix==NULL)
// hpix=makeNorm();
TH2D *hg;
char nn[1000];
if (strcmp(name,"blank")==NULL) {
hg=hff;
} else {
sprintf(nn,"/mnt/moench_data/tomcat_20160528_img/img_%s_eta_nb25.root", name);
// if (strcmp(name,"blank"))
TFile *fg=new TFile(nn);
// else
// TFile fg=gDirectory;
//TFile fg("/local_zfs_raid/tomcat_20160528/trees/img_grating_1d_eta_gmap.root");
//TH2D *hg=(TH2D*)fg.Get("imgHR");
// TFile fg("/local_zfs_raid/tomcat_20160528/trees/img_grating_1d_eta_nb25.root");
// 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);
// 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.;
int ip=0;
int max=0;
Double_t avg=0, rms=0;
for (int iby=0; iby<hg->GetNbinsY(); iby++) {
for (int ibx=0; ibx<hg->GetNbinsX(); ibx++) {
imageData[ip]=hg->GetBinContent(ibx+1,iby+1);
if (imageData[ip]>max) max=imageData[ip];
// if (imageData[ip]>3000) imageData[ip]=3000.;
avg+=((Double_t)imageData[ip])/(hg->GetNbinsY()*hg->GetNbinsX());
rms+=((Double_t)imageData[ip])*((Double_t)imageData[ip])/(hg->GetNbinsY()*hg->GetNbinsX());
ip++;
}
}
rms=TMath::Sqrt(rms-avg*avg);
ip=0;
for (int iby=0; iby<hg->GetNbinsY(); iby++) {
for (int ibx=0; ibx<hg->GetNbinsX(); ibx++) {
if (imageData[ip]>avg+nsigma*rms) imageData[ip]=avg+nsigma*rms;
// if (imageData[ip]>3000) imageData[ip]=3000.;
ip++;
}
}
cout << "MAXIMUM IS "<< max << endl;
cout << "AVERAGE IS "<< avg << endl;
cout << "RMS IS "<< rms << endl;
int nbx=hg->GetNbinsX();
int nby=hg->GetNbinsY();
Short_t *buffer=new Short_t[nbx];
// // cout << "Size of short int is "<<sizeof(char)<< endl;
// // cout << "width is "<<nbx<< endl;
// // cout << "height is "<<nby<< endl;
sprintf(nn,"%s_HR.bin", name);
ofstream myFile (nn, ios::out | ios::binary);
ip=0;
for (int iy=0; iy<nby; iy++) {
for (int ix=0; ix<nbx; ix++) {
buffer[ix]=imageData[ip]*65535/nsigma/avg;
ip++;
}
myFile.write((char*)buffer, nbx*sizeof(Short_t));
}
myFile.close();
TASImage *img=new TASImage ( "img",imageData, 400*25, 200*25);// TImagePalette * palette = 0
new TCanvas();
img->SetImageCompression(0);
img->SetImageQuality(TAttImage::kImgBest);
// img->Gray(kTRUE);
img->Draw();
sprintf(nn,"%s_HR.tiff", name);
img->WriteImage(nn, TImage::kTiff);
// new TCanvas();
// hg->SetMaximum(3000);
// hg->DrawCopy("colz");
hg->Rebin2D(25,25);
Double_t imageDataLR[200*400];
ip=0; max=0;avg=0; rms=0;
for (int iby=0; iby<hg->GetNbinsY(); iby++) {
for (int ibx=0; ibx<hg->GetNbinsX(); ibx++) {
imageDataLR[ip]=hg->GetBinContent(ibx+1,iby+1);
if (imageDataLR[ip]>max) max=imageDataLR[ip];
avg+=((Double_t)imageDataLR[ip])/(hg->GetNbinsY()*hg->GetNbinsX());
rms+=((Double_t)imageDataLR[ip])*((Double_t)imageDataLR[ip])/(hg->GetNbinsY()*hg->GetNbinsX());
ip++;
}
}
rms=TMath::Sqrt(rms-avg*avg);
ip=0;
for (int iby=0; iby<hg->GetNbinsY(); iby++) {
for (int ibx=0; ibx<hg->GetNbinsX(); ibx++) {
if (imageDataLR[ip]>avg+nsigma*rms) imageDataLR[ip]=avg+nsigma*rms;
// if (imageData[ip]>3000) imageData[ip]=3000.;
ip++;
}
}
cout << "MAXIMUM IS "<< max << endl;
cout << "AVERAGE IS "<< avg << endl;
cout << "RMS IS "<< rms << endl;
TASImage *imgLR=new TASImage ( "imgLR",imageDataLR, 400, 200);// TImagePalette * palette = 0
new TCanvas();
imgLR->SetImageCompression(0);
imgLR->SetImageQuality(TAttImage::kImgBest);
// imgLR->Gray(kTRUE);
imgLR->Draw();
sprintf(nn,"%s_LR.tiff", name);
imgLR->WriteImage(nn, TImage::kTiff);
int nbx1=hg->GetNbinsX();
int nby1=hg->GetNbinsY();
Short_t *buffer1=new Short_t[nbx1];
// // cout << "Size of short int is "<<sizeof(char)<< endl;
// // cout << "width is "<<nbx<< endl;
// // cout << "height is "<<nby<< endl;
sprintf(nn,"%s_LR.bin", name);
ofstream myFile1 (nn, ios::out | ios::binary);
ip=0;
for (int iy=0; iy<nby1; iy++) {
for (int ix=0; ix<nbx1; ix++) {
buffer1[ix]=imageDataLR[ip]*256/nsigma/avg;
ip++;
}
myFile1.write((char*)buffer1, nbx1*sizeof(Short_t));
}
myFile1.close();
cout << sizeof(Short_t) << endl;
//for grating2D the gratings are in the pixels (150-260)x(80-190)
// // hg->Draw("colz");
// int off=13;
// // hg->Divide(hpix);
// int ix,iy, sbx, sby, ibx2, iby2, sbx2, sby2;
// TH2F *hh13b=new TH2F("hh13b","hh13b",400*24,0,400,200*24,100,300);
// TH2F *hh13=new TH2F("hh13","hh13",400*23,0,400,200*23,100,300);
// TH2F *hh=new TH2F("hh","hh",400*21,0,400,200*21,100,300);
// TH2F *h1=new TH2F("h1","h1",400*23,0,400,200*23,100,300);
// for (int ibx=0; ibx<hg->GetNbinsX(); ibx++) {
// for (int iby=0; iby<hg->GetNbinsY(); iby++) {
// ix=(ibx-off)/25;
// iy=(iby-off)/25;
// sbx=(ibx-off)%25;
// sby=(iby-off)%25;
// sbx2=sbx-1;
// sby2=sby-1;
// if (sbx2<0) sbx2=0;
// if (sby2<0) sby2=0;
// if (sbx2>22) sbx2=22;
// if (sby2>22) sby2=22;
// ibx2=ix*23+(sbx2+off);
// iby2=iy*23+(sby2+off);
// hh13->Fill(hh13->GetXaxis()->GetBinCenter(ibx2+1),hh13->GetYaxis()->GetBinCenter(iby2+1),hg->GetBinContent(ibx+1,iby+1));
// // h1->Fill(h1->GetXaxis()->GetBinCenter(ibx2+1),h1->GetYaxis()->GetBinCenter(iby2+1),hpix->GetBinContent(ibx+1,iby+1));
// off=13;
// ix=(ibx-off)/25;
// iy=(iby-off)/25;
// sbx=(ibx-off)%25;
// sby=(iby-off)%25;
// sbx2=sbx-1;
// sby2=sby-1;
// // if (sbx2<0) sbx2=0;
// // if (sby2<0) sby2=0;
// ibx2=ix*24+(sbx2+off);
// iby2=iy*24+(sby2+off);
// if (sbx2<0 && sby2>=0) {
// ibx2=ix*24+(sbx2+off);
// hh13b->Fill(hh13b->GetXaxis()->GetBinCenter(ibx2+1),hh13b->GetYaxis()->GetBinCenter(iby2+1),0.5*hg->GetBinContent(ibx+1,iby+1));
// ibx2=ix*24+(sbx2+off)+1;
// hh13b->Fill(hh13b->GetXaxis()->GetBinCenter(ibx2+1),hh13b->GetYaxis()->GetBinCenter(iby2+1),0.5*hg->GetBinContent(ibx+1,iby+1));
// } else if (sby2<0 && sbx2>=0){
// iby2=iy*24+(sby2+off);
// hh13b->Fill(hh13b->GetXaxis()->GetBinCenter(ibx2+1),hh13b->GetYaxis()->GetBinCenter(iby2+1),0.5*hg->GetBinContent(ibx+1,iby+1));
// iby2=iy*24+(sby2+off)+1;
// hh13b->Fill(hh13b->GetXaxis()->GetBinCenter(ibx2+1),hh13b->GetYaxis()->GetBinCenter(iby2+1),0.5*hg->GetBinContent(ibx+1,iby+1));
// } else if (sby2<0 && sbx2<0){
// iby2=iy*24+(sby2+off);
// ibx2=ix*24+(sbx2+off);
// hh13b->Fill(hh13b->GetXaxis()->GetBinCenter(ibx2+1),hh13b->GetYaxis()->GetBinCenter(iby2+1),0.25*hg->GetBinContent(ibx+1,iby+1));
// iby2=iy*24+(sby2+off)+1;
// ibx2=ix*24+(sbx2+off);
// hh13b->Fill(hh13b->GetXaxis()->GetBinCenter(ibx2+1),hh13b->GetYaxis()->GetBinCenter(iby2+1),0.25*hg->GetBinContent(ibx+1,iby+1));
// iby2=iy*24+(sby2+off);
// ibx2=ix*24+(sbx2+off)+1;
// hh13b->Fill(hh13b->GetXaxis()->GetBinCenter(ibx2+1),hh13b->GetYaxis()->GetBinCenter(iby2+1),0.25*hg->GetBinContent(ibx+1,iby+1));
// iby2=iy*24+(sby2+off)+1;
// ibx2=ix*24+(sbx2+off)+1;
// hh13b->Fill(hh13b->GetXaxis()->GetBinCenter(ibx2+1),hh13b->GetYaxis()->GetBinCenter(iby2+1),0.25*hg->GetBinContent(ibx+1,iby+1));
// }else {
// hh13b->Fill(hh13b->GetXaxis()->GetBinCenter(ibx2+1),hh13b->GetYaxis()->GetBinCenter(iby2+1),hg->GetBinContent(ibx+1,iby+1));
// }
// ix=(ibx-off)/25;
// iy=(iby-off)/25;
// sbx=(ibx-off)%25;
// sby=(iby-off)%25;
// sbx2=sbx-2;
// sby2=sby-2;
// if (sbx2<0) sbx2=-1;
// if (sby2<0) sby2=-1;
// if (sbx2>20) sbx2=-1;
// if (sby2>20) sby2=-1;
// ibx2=ix*21+(sbx2+off);
// iby2=iy*21+(sby2+off);
// if (sbx2>=0 && sby2>=0) hh->Fill(hh->GetXaxis()->GetBinCenter(ibx2+1),hh->GetYaxis()->GetBinCenter(iby2+1),hg->GetBinContent(ibx+1,iby+1));
// }
// }
// new TCanvas();
// hg->GetXaxis()->SetRangeUser(84,87);
// hg->GetYaxis()->SetRangeUser(120,124);
// hg->Draw("colz");
// new TCanvas();
// // hh->Divide(h1);
// hh->GetXaxis()->SetRangeUser(84,87);
// hh->GetYaxis()->SetRangeUser(120,124);
// hh->Draw("colz");
// new TCanvas();
// // hh->Divide(h1);
// hh13b->GetXaxis()->SetRangeUser(84,87);
// hh13b->GetYaxis()->SetRangeUser(120,124);
// hh13b->Draw("colz");
// new TCanvas();
// // hh->Divide(h1);
// hh13->GetXaxis()->SetRangeUser(84,87);
// hh13->GetYaxis()->SetRangeUser(120,124);
// hh13->Draw("colz");
// // //TFile fsg("/local_zfs_raid/tomcat_20160528/trees/img_sample_grating_1d_eta_gmap.root");
// // //TH2D *hsg=(TH2D*)fsg.Get("sample_grating_1dHR");
// // TFile fsg("/local_zfs_raid/tomcat_20160528/trees/img_sample_grating_1d_eta.root");
// // TH2D *hsg=(TH2D*)fsg.Get("imgHR");
// // hsg->SetName("imgGratingSample");
// // hsg->Divide(hpix);
// // Double_t nf=hff->Integral(100,hff->GetNbinsX()-100,hff->GetYaxis()->FindBin(130), hff->GetYaxis()->FindBin(140))/((hff->GetNbinsX()-200)*250);
// // hff->Scale(1./nf);
// // // hg->Divide(hff);
// // //hsg->Divide(hff);
// // Double_t ng=hg->Integral(100,hg->GetNbinsX()-100,hg->GetYaxis()->FindBin(130), hg->GetYaxis()->FindBin(140))/((hg->GetNbinsX()-200)*250);
// // Double_t nsg=hsg->Integral(100,hsg->GetNbinsX()-100,hsg->GetYaxis()->FindBin(130), hsg->GetYaxis()->FindBin(140))/((hsg->GetNbinsX()-200)*250);
// // hg->Scale(1/ng);
// // hsg->Scale(1/nsg);
// // // new TCanvas();
// // // hg->SetMaximum(1.);
// // // hg->Draw("colz");
// // // new TCanvas();
// // // hsg->SetMaximum(1.);
// // // hsg->Draw("colz");
// // // TH1D *pg=hg->ProjectionX("pg",hg->GetYaxis()->FindBin(174.5)+1,hg->GetYaxis()->FindBin(175.5));
// // // TH1D *psg=hsg->ProjectionX("psg",hsg->GetYaxis()->FindBin(174.5)+1,hsg->GetYaxis()->FindBin(175.5));
// // // psg->SetLineColor(2);
// // // new TCanvas();
// // // pg->Draw("l");
// // // psg->Draw("l same");
// // int nbx=hg->GetNbinsX();
// // int nby=hg->GetNbinsY();
// // char buffer[nbx];
// // cout << "Size of short int is "<<sizeof(char)<< endl;
// // cout << "width is "<<nbx<< endl;
// // cout << "height is "<<nby<< endl;
// // ofstream myFile ("grating_1d.bin", ios::out | ios::binary);
// // for (int iy=0; iy<nby; iy++) {
// // for (int ix=0; ix<nbx; ix++) {
// // buffer[ix]=hg->GetBinContent(ix+1,iy+1);
// // }
// // myFile.write((char*)buffer, nbx*sizeof(char));
// // }
// // myFile.close();
// // nbx=hsg->GetNbinsX();
// // nby=hsg->GetNbinsY();
// // cout << "Size of short int is "<<sizeof(char)<< endl;
// // cout << "width is "<<nbx<< endl;
// // cout << "height is "<<nby<< endl;
// // ofstream myFile1 ("sample_grating_1d.bin", ios::out | ios::binary);
// // for (int iy=0; iy<nby; iy++) {
// // for (int ix=0; ix<nbx; ix++) {
// // buffer[ix]=hsg->GetBinContent(ix+1,iy+1);
// // }
// // myFile1.write((char*)buffer, nbx*sizeof(char));
// // }
// // myFile1.close();
// // NAME
// // raw2tiff - create a TIFF file from a raw data
// // SYNOPSIS
// // raw2tiff [ options ] input.raw output.tif
// // DESCRIPTION
// // raw2tiff converts a raw byte sequence into TIFF. By default, the TIFF image is created with data samples packed (PlanarConfiguration=1), compressed with the PackBits algorithm (Compression=32773),
// // and with each strip no more than 8 kilobytes. These characteristics can overridden, or explicitly specified with the options described below.
// // OPTIONS
// // -H number
// // size of input image file header in bytes (0 by default). This amount of data just will be skipped from the start of file while reading.
// // -w number
// // width of input image in pixels (can be guessed, see GUESSING THE IMAGE GEOMETRY below).
// // -l number
// // length of input image in lines (can be guessed, see GUESSING THE IMAGE GEOMETRY below).
// // -b number
// // number of bands in input image (1 by default).
// // -d data_type
// // type of samples in input image, where data_type may be:
// // byte 8-bit unsigned integer (default),
// // short 16-bit unsigned integer,
// // long 32-bit unsigned integer,
// // sbyte 8-bit signed integer,
// // sshort 16-bit signed integer,
// // slong 32-bit signed integer,
// // float 32-bit IEEE floating point,
// // double 64-bit IEEE floating point.
// // -i config
// // type of samples interleaving in input image, where config may be:
// // pixel pixel interleaved data (default),
// // band band interleaved data.
// // -p photo
// // photometric interpretation (color space) of the input image, where photo may be:
// // miniswhite white color represented with 0 value,
// // minisblack black color represented with 0 value (default),
// // rgb image has RGB color model,
// // cmyk image has CMYK (separated) color model,
// // ycbcr image has YCbCr color model,
// // cielab image has CIE L*a*b color model,
// // icclab image has ICC L*a*b color model,
// // itulab image has ITU L*a*b color model.
// // -s swap bytes fetched from the input file.
// // -L input data has LSB2MSB bit order (default).
// // -M input data has MSB2LSB bit order.
// // -c Specify a compression scheme to use when writing image data: -c none for no compression, -c packbits for the PackBits compression algorithm (the default), -c jpeg for the baseline JPEG compres-
// // sion algorithm, -c zip for the Deflate compression algorithm, and -c lzw for Lempel-Ziv & Welch.
// // -r number
// // Write data with a specified number of rows per strip; by default the number of rows/strip is selected so that each strip is approximately 8 kilobytes.
// //width is 10000
// //height is 5000
}

185
slsDetectorCalibration/etaVEL/imageMacroG2D.C
View File

@ -1,185 +0,0 @@
{
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("/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);
for (int ibx=10*25; ibx<hff->GetNbinsX()-10*25; ibx++) {
for (int iby=20*25; iby<hff->GetNbinsY()-20*25; iby++) {
hpixel->Fill((ibx-12)%25,(iby-12)%25,hff->GetBinContent(ibx+1,iby+1));
}
}
hpixel->Scale(1./hpixel->GetBinContent(13,13));
// new TCanvas();
// hpixel->Draw("colz");
TH2D *hraw=(TH2D*)hff->Clone("hraw");
TH2D *hpix=(TH2D*)hff->Clone("hpix");
for (int ibx=0; ibx<hff->GetNbinsX(); ibx++) {
for (int iby=0; iby<hff->GetNbinsY(); iby++) {
hpix->SetBinContent(ibx+1,iby+1,hpixel->GetBinContent(hpixel->GetXaxis()->FindBin((ibx-12)%25),hpixel->GetXaxis()->FindBin((iby-12)%25)));
}
}
hff->Divide(hpix);
// new TCanvas();
// hff->Draw("colz");
// TFile fg("/local_zfs_raid/tomcat_20160528/trees/img_grating_2d_eta.root");
// TH2D *hg=(TH2D*)fg.Get("imgHR");
TFile fg("/mnt/moench_data/tomcat_20160528_img/img_grating_2d_eta_nb25.root");
TH2D *hg=(TH2D*)fg.Get("grating_2dHR");
hg->SetName("imgGrating");
hg->Divide(hpix);
Double_t nf=hff->Integral(100,hff->GetNbinsX()-100,hff->GetYaxis()->FindBin(130), hff->GetYaxis()->FindBin(140))/((hff->GetNbinsX()-200)*250);
hff->Scale(1./nf);
// hg->Divide(hff);
//hsg->Divide(hff);
Double_t ng=hg->Integral(100,hg->GetNbinsX()-100,hg->GetYaxis()->FindBin(130), hg->GetYaxis()->FindBin(140))/((hg->GetNbinsX()-200)*250);
hg->Scale(1/ng);
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");
// TH1D *pg=hg->ProjectionX("pg",hg->GetYaxis()->FindBin(174.5)+1,hg->GetYaxis()->FindBin(175.5));
// TH1D *psg=hsg->ProjectionX("psg",hsg->GetYaxis()->FindBin(174.5)+1,hsg->GetYaxis()->FindBin(175.5));
// psg->SetLineColor(2);
// new TCanvas();
// pg->Draw("l");
// psg->Draw("l same");
int nbx=hg->GetNbinsX();
int nby=hg->GetNbinsY();
char buffer[nbx];
cout << "Size of short int is "<<sizeof(char)<< endl;
cout << "width is "<<nbx<< endl;
cout << "height is "<<nby<< endl;
ofstream myFile ("grating_2d.bin", ios::out | ios::binary);
for (int iy=0; iy<nby; iy++) {
for (int ix=0; ix<nbx; ix++) {
buffer[ix]=hg->GetBinContent(ix+1,iy+1);
}
myFile.write((char*)buffer, nbx*sizeof(char));
}
myFile.close();
// NAME
// raw2tiff - create a TIFF file from a raw data
// SYNOPSIS
// raw2tiff [ options ] input.raw output.tif
// DESCRIPTION
// raw2tiff converts a raw byte sequence into TIFF. By default, the TIFF image is created with data samples packed (PlanarConfiguration=1), compressed with the PackBits algorithm (Compression=32773),
// and with each strip no more than 8 kilobytes. These characteristics can overridden, or explicitly specified with the options described below.
// OPTIONS
// -H number
// size of input image file header in bytes (0 by default). This amount of data just will be skipped from the start of file while reading.
// -w number
// width of input image in pixels (can be guessed, see GUESSING THE IMAGE GEOMETRY below).
// -l number
// length of input image in lines (can be guessed, see GUESSING THE IMAGE GEOMETRY below).
// -b number
// number of bands in input image (1 by default).
// -d data_type
// type of samples in input image, where data_type may be:
// byte 8-bit unsigned integer (default),
// short 16-bit unsigned integer,
// long 32-bit unsigned integer,
// sbyte 8-bit signed integer,
// sshort 16-bit signed integer,
// slong 32-bit signed integer,
// float 32-bit IEEE floating point,
// double 64-bit IEEE floating point.
// -i config
// type of samples interleaving in input image, where config may be:
// pixel pixel interleaved data (default),
// band band interleaved data.
// -p photo
// photometric interpretation (color space) of the input image, where photo may be:
// miniswhite white color represented with 0 value,
// minisblack black color represented with 0 value (default),
// rgb image has RGB color model,
// cmyk image has CMYK (separated) color model,
// ycbcr image has YCbCr color model,
// cielab image has CIE L*a*b color model,
// icclab image has ICC L*a*b color model,
// itulab image has ITU L*a*b color model.
// -s swap bytes fetched from the input file.
// -L input data has LSB2MSB bit order (default).
// -M input data has MSB2LSB bit order.
// -c Specify a compression scheme to use when writing image data: -c none for no compression, -c packbits for the PackBits compression algorithm (the default), -c jpeg for the baseline JPEG compres-
// sion algorithm, -c zip for the Deflate compression algorithm, and -c lzw for Lempel-Ziv & Welch.
// -r number
// Write data with a specified number of rows per strip; by default the number of rows/strip is selected so that each strip is approximately 8 kilobytes.
//width is 10000
//height is 5000
}

173
slsDetectorCalibration/etaVEL/imageMacroSample.C
View File

@ -1,173 +0,0 @@
{
// 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");
// 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++) {
// for (int iby=20*25; iby<hff->GetNbinsY()-20*25; iby++) {
// hpixel->Fill((ibx-12)%25,(iby-12)%25,hff->GetBinContent(ibx+1,iby+1));
// }
// }
// hpixel->Scale(1./hpixel->GetBinContent(13,13));
// // new TCanvas();
// // hpixel->Draw("colz");
// TH2D *hraw=(TH2D*)hff->Clone("hraw");
// TH2D *hpix=(TH2D*)hff->Clone("hpix");
// for (int ibx=0; ibx<hff->GetNbinsX(); ibx++) {
// for (int iby=0; iby<hff->GetNbinsY(); iby++) {
// hpix->SetBinContent(ibx+1,iby+1,hpixel->GetBinContent(hpixel->GetXaxis()->FindBin((ibx-12)%25),hpixel->GetXaxis()->FindBin((iby-12)%25)));
// }
// }
// hff->Divide(hpix);
// Double_t nf=hff->Integral(100,hff->GetNbinsX()-100,hff->GetYaxis()->FindBin(130), hff->GetYaxis()->FindBin(140))/((hff->GetNbinsX()-200)*250);
// hff->Scale(1./nf);
// new TCanvas();
// hff->Draw("colz");
// TFile fg("/local_zfs_raid/tomcat_20160528/trees/img_grating_2d_eta.root");
// TH2D *hg=(TH2D*)fg.Get("imgHR");
TFile fg("/mnt/moench_data/tomcat_20160528_img/img_sample_eta_gmap_v2.root");
TH2D *hg=(TH2D*)fg.Get("sampleHR");
hg->SetName("imgSample");
// hg->Divide(hpix);
// hg->Divide(hff);
//hsg->Divide(hff);
Double_t ng=hg->Integral(100,hg->GetNbinsX()-100,hg->GetYaxis()->FindBin(130), hg->GetYaxis()->FindBin(140))/((hg->GetNbinsX()-200)*250);
hg->Scale(1/ng);
// new TCanvas();
// hsg->SetMaximum(1.);
// hsg->Draw("colz");
// TH1D *pg=hg->ProjectionX("pg",hg->GetYaxis()->FindBin(174.5)+1,hg->GetYaxis()->FindBin(175.5));
// TH1D *psg=hsg->ProjectionX("psg",hsg->GetYaxis()->FindBin(174.5)+1,hsg->GetYaxis()->FindBin(175.5));
// psg->SetLineColor(2);
// new TCanvas();
// pg->Draw("l");
// psg->Draw("l same");
TH2D *hpixel1=new TH2D("hpixel1","hpixel1",25,0,25,25,0,25);
for (int ibx=10*25; ibx<35*25; ibx++) {
for (int iby=25*25; iby<50*25; iby++) {
hpixel1->Fill((ibx-12)%25,(iby-12)%25,hg->GetBinContent(ibx+1,iby+1));
}
}
hpixel1->Scale(1./hpixel1->GetBinContent(13,13));
// new TCanvas();
// hpixel->Draw("colz");
TH2D *hpix1=(TH2D*)hg->Clone("hpix1");
for (int ibx=0; ibx<hg->GetNbinsX(); ibx++) {
for (int iby=0; iby<hg->GetNbinsY(); iby++) {
hpix1->SetBinContent(ibx+1,iby+1,hpixel1->GetBinContent(hpixel1->GetXaxis()->FindBin((ibx-12)%25),hpixel1->GetXaxis()->FindBin((iby-12)%25)));
}
}
hg->Divide(hpix1);
new TCanvas();
hg->SetMaximum(4);
hg->SetMinimum(0.5);
hg->Draw("colz");
int nbx=hg->GetNbinsX();
int nby=hg->GetNbinsY();
char buffer[nbx];
cout << "Size of short int is "<<sizeof(char)<< endl;
cout << "width is "<<nbx<< endl;
cout << "height is "<<nby<< endl;
ofstream myFile ("sample.bin", ios::out | ios::binary);
for (int iy=0; iy<nby; iy++) {
for (int ix=0; ix<nbx; ix++) {
buffer[ix]=hg->GetBinContent(ix+1,iy+1);
}
myFile.write((char*)buffer, nbx*sizeof(char));
}
myFile.close();
// NAME
// raw2tiff - create a TIFF file from a raw data
// SYNOPSIS
// raw2tiff [ options ] input.raw output.tif
// DESCRIPTION
// raw2tiff converts a raw byte sequence into TIFF. By default, the TIFF image is created with data samples packed (PlanarConfiguration=1), compressed with the PackBits algorithm (Compression=32773),
// and with each strip no more than 8 kilobytes. These characteristics can overridden, or explicitly specified with the options described below.
// OPTIONS
// -H number
// size of input image file header in bytes (0 by default). This amount of data just will be skipped from the start of file while reading.
// -w number
// width of input image in pixels (can be guessed, see GUESSING THE IMAGE GEOMETRY below).
// -l number
// length of input image in lines (can be guessed, see GUESSING THE IMAGE GEOMETRY below).
// -b number
// number of bands in input image (1 by default).
// -d data_type
// type of samples in input image, where data_type may be:
// byte 8-bit unsigned integer (default),
// short 16-bit unsigned integer,
// long 32-bit unsigned integer,
// sbyte 8-bit signed integer,
// sshort 16-bit signed integer,
// slong 32-bit signed integer,
// float 32-bit IEEE floating point,
// double 64-bit IEEE floating point.
// -i config
// type of samples interleaving in input image, where config may be:
// pixel pixel interleaved data (default),
// band band interleaved data.
// -p photo
// photometric interpretation (color space) of the input image, where photo may be:
// miniswhite white color represented with 0 value,
// minisblack black color represented with 0 value (default),
// rgb image has RGB color model,
// cmyk image has CMYK (separated) color model,
// ycbcr image has YCbCr color model,
// cielab image has CIE L*a*b color model,
// icclab image has ICC L*a*b color model,
// itulab image has ITU L*a*b color model.
// -s swap bytes fetched from the input file.
// -L input data has LSB2MSB bit order (default).
// -M input data has MSB2LSB bit order.
// -c Specify a compression scheme to use when writing image data: -c none for no compression, -c packbits for the PackBits compression algorithm (the default), -c jpeg for the baseline JPEG compres-
// sion algorithm, -c zip for the Deflate compression algorithm, and -c lzw for Lempel-Ziv & Welch.
// -r number
// Write data with a specified number of rows per strip; by default the number of rows/strip is selected so that each strip is approximately 8 kilobytes.
//width is 10000
//height is 5000
}

3109
slsDetectorCalibration/etaVEL/interpolationMacro.C
View File

File diff suppressed because it is too large Load Diff

11
slsDetectorCalibration/gMapDemo.C
View File

@ -1,11 +0,0 @@
{
//.L energyCalibration.cpp+
//.L gainMap.C+
TFile fin("/data/moench_xbox_20140113/MoTarget_45kV_0_8mA_120V_cds_g4.root");
TH2F *h2=fin.Get("h2");
TH2F *gMap=gainMap(h2,4);
gMap->Draw("colz");
}

228
slsDetectorCalibration/gainMap.C
View File

@ -1,228 +0,0 @@
#define MYROOT
#include <TH2F.h>
#include <TSpectrum.h>
#include <TH1D.h>
#include <TFile.h>
#include <TList.h>
#include <TPolyMarker.h>
#include <THStack.h>
#include <TF1.h>
#include <TGraphErrors.h>
#include <iostream>
#include <energyCalibration.h>
using namespace std;
#define FOPT "0"
TH2F *gainMap(TH2F *h2, float g) {
// int npx, npy;
int npx=160, npy=160;
// det->getDetectorSize(npx, npy);
TH2F *gMap=new TH2F("gmap",h2->GetTitle(), npx, -0.5, npx-0.5, npy, -0.5, npy-0.5);
Double_t ens[3]={0,8,17.5}, eens[3]={0.,0.1,0.1};
Double_t peaks[3], epeaks[3];
int ibin;
TH1D *px;
energyCalibration *enCal=new energyCalibration();
enCal->setPlotFlag(0);
// enCal->setChargeSharing(0);
enCal->setScanSign(1);
Double_t gain, off, egain, eoff;
TList *functions;
TPolyMarker *pm ;
Double_t *peakX, *peakY;
TSpectrum *s=new TSpectrum();
Double_t mypar[10], emypar[10];
Double_t prms, np;
int iit=0;
TGraph *glin;
Double_t peakdum, hpeakdum;
for (int ix=1; ix<npx-1; ix++) {
for (int iy=1; iy<npy-1; iy++) {
// cout << ix << " " << iy << " " << ibin << endl;
ibin=ix*npy+iy;
px=h2->ProjectionX("px",ibin+1,ibin+1);
prms=10*g;
iit=0;
np=s->Search(px,prms,"",0.2);
while (np !=2) {
if (np>2)
prms+=0.5*prms;
else
prms-=0.5*prms;
iit++;
if (iit>=10)
break;
np=s->Search(px,prms,"",0.2);
}
if (np!=2)
cout << "peak search could not converge " << ibin << endl;
if (np==2) {
pm=NULL;
functions=px->GetListOfFunctions();
if (functions)
pm = (TPolyMarker*)functions->FindObject("TPolyMarker");
if (pm) {
peakX=pm->GetX();
peakY=pm->GetY();
if (peakX[0]>peakX[1]) {
peakdum=peakX[0];
hpeakdum=peakY[0];
peakX[0]= peakX[1];
peakY[0]= peakY[1];
peakX[1]= peakdum;
peakY[1]= hpeakdum;
}
cout << "("<< ix << "," << iy << ") " << endl;
for (int ip=0; ip<np; ip++) {
// cout << ip << " " << peakX[ip] << " " << peakY[ip] << endl;
enCal->setFitRange(peakX[ip]-10*g,peakX[ip]+10*g);
mypar[0]=0;
mypar[1]=0;
mypar[2]=peakX[ip];
mypar[3]=g*10;
mypar[4]=peakY[ip];
mypar[5]=0;
enCal->setStartParameters(mypar);
enCal->fitSpectrum(px,mypar,emypar);
peaks[ip+1]=mypar[2];
epeaks[ip+1]=emypar[2];
}
peaks[0]=0;
epeaks[0]=1;
// for (int i=0; i<3; i++) cout << i << " " << ens[i] << " " << eens[i]<< " "<< peaks[i]<< " " << epeaks[i] << endl;
glin= enCal->linearCalibration(3,ens,eens,peaks,epeaks,gain,off,egain,eoff);
// cout << "Gain " << gain << " off " << off << endl;
if (off>-10 && off<10) {
gMap->SetBinContent(ix+1,iy+1,gain);
gMap->SetBinError(ix+1,iy+1,egain);
}
if (glin)
delete glin;
}
}
}
}
return gMap;
}
TH2F *noiseMap(TH2F *h2) {
// int npx, npy;
int npx=160, npy=160;
// det->getDetectorSize(npx, npy);
TH2F *nMap=new TH2F("nmap",h2->GetTitle(), npx, -0.5, npx-0.5, npy, -0.5, npy-0.5);
int ibin;
TH1D *px;
for (int ix=0; ix<npx; ix++) {
for (int iy=0; iy<npy; iy++) {
cout << ix << " " << iy << " " << ibin << endl;
ibin=ix*npy+iy;
px=h2->ProjectionX("px",ibin+1,ibin+1);
px->Fit("gaus", FOPT);
if (px->GetFunction("gaus")) {
nMap->SetBinContent(ix+1,iy+1,px->GetFunction("gaus")->GetParameter(2));
}
// delete px;
}
}
return nMap;
}
THStack *noiseHistos(char *tit) {
char fname[10000];
sprintf(fname,"/data/moench_xbox_20140116/noise_map_%s.root",tit);
TFile *fn=new TFile(fname);
TH2F *nmap=(TH2F*)fn->Get("nmap");
if (nmap==NULL) {
cout << "No noise map in file " << fname << endl;
return NULL;
}
sprintf(fname,"/data/moench_xbox_20140113/gain_map_%s.root",tit);
TFile *fg=new TFile(fname);
TH2F *gmap=(TH2F*)fg->Get("gmap");
if (gmap==NULL) {
cout << "No gain map in file " << fname << endl;
return NULL;
}
nmap->Divide(gmap);
nmap->Scale(1000./3.6);
THStack *hs=new THStack(tit,tit);
hs->SetTitle(tit);
TH1F *h;
char hname[100];
cout << tit << endl;
for (int is=0; is<4; is++) {
sprintf(hname,"h%ds",is+1);
h=new TH1F(hname,tit,500,0,500);
hs->Add(h);
// cout << hs->GetHists()->GetEntries() << endl;
for (int ix=40*is+2; ix<40*(is+1)-2; ix++) {
for (int iy=2; iy<158; iy++) {
if (ix<100 || ix>120)
h->Fill(nmap->GetBinContent(ix+1,iy+1));
}
}
cout << is+1 << "SC: " << "" << h->GetMean() << "+-" << h->GetRMS();
h->Fit("gaus","0Q");
h->SetLineColor(is+1);
if (h->GetFunction("gaus")) {
h->GetFunction("gaus")->SetLineColor(is+1);
cout << " or " << h->GetFunction("gaus")->GetParameter(1) << "+-" << h->GetFunction("gaus")->GetParError(1);
}
cout << endl;
}
// cout << hs->GetHists()->GetEntries() << endl;
return hs;
}

BIN
slsDetectorCalibration/interpolations/.__afs3A48 Executable file
View File

Binary file not shown.

93
slsDetectorCalibration/etaVEL/etaInterpolationBase.h → slsDetectorCalibration/interpolations/etaInterpolationBase.h
View File

@ -29,8 +29,8 @@ class etaInterpolationBase : public slsInterpolation {
#endif
#ifndef MYROOT1
heta=new int[nbeta*nbeta];
hhx=new int[nbeta*nbeta];
hhy=new int[nbeta*nbeta];
hhx=new float[nbeta*nbeta];
hhy=new float[nbeta*nbeta];
#endif
@ -51,10 +51,10 @@ class etaInterpolationBase : public slsInterpolation {
#ifndef MYROOT1
heta=new int[nbeta*nbeta];
memcpy(heta,orig->heta,nbeta*nbeta*sizeof(int));
hhx=new int[nbeta*nbeta];
memcpy(hhx,orig->hhx,nbeta*nbeta*sizeof(int));
hhy=new int[nbeta*nbeta];
memcpy(hhy,orig->hhy,nbeta*nbeta*sizeof(int));
hhx=new float[nbeta*nbeta];
memcpy(hhx,orig->hhx,nbeta*nbeta*sizeof(float));
hhy=new float[nbeta*nbeta];
memcpy(hhy,orig->hhy,nbeta*nbeta*sizeof(float));
#endif
@ -148,7 +148,7 @@ class etaInterpolationBase : public slsInterpolation {
if (nnx!=nny) {
cout << "different number of bins in x " << nnx << " and y " << nny<< " !"<< endl;
cout << "Aborting read"<< endl;
return NULL;
return 0;
}
nbeta=nnx;
if (gm) {
@ -159,8 +159,8 @@ class etaInterpolationBase : public slsInterpolation {
}
heta=new int[nbeta*nbeta];
hhx=new int[nbeta*nbeta];
hhy=new int[nbeta*nbeta];
hhx=new float[nbeta*nbeta];
hhy=new float[nbeta*nbeta];
for (int ix=0; ix<nbeta; ix++) {
for (int iy=0; iy<nbeta; iy++) {
@ -170,7 +170,7 @@ class etaInterpolationBase : public slsInterpolation {
delete [] gm;
return 1;
}
return NULL;
return 0;
};
@ -199,13 +199,13 @@ TH2D *gethhx()
#endif
#ifndef MYROOT1
int *gethhx()
float *gethhx()
{
// hhx->Scale((double)nSubPixels);
return hhx;
};
int *gethhy()
float *gethhy()
{
// hhy->Scale((double)nSubPixels);
return hhy;
@ -285,20 +285,25 @@ int *gethhx()
switch (corner)
{
case TOP_LEFT:
dX=-.99;
dY=+.99;
dX=-1.;//.99;
dY=0;//+1.;//.99;
//etay=1-etay;
break;
case TOP_RIGHT:
dX=+.99;
dY=+.99;
;
dX=0;//+1.;//.99;
dY=0;//+1.;//.99;
break;
case BOTTOM_LEFT:
dX=-.99;
dY=-.99;
//etax=1-etax;
//etay=1-etay;
dX=-1.;//99;
dY=-1.;//.99;
break;
case BOTTOM_RIGHT:
dX=+.99;
dY=-.99;
//etay=1-etay;
dX=0;//1.;//+.99;
dY=-1.;//-.99;
break;
default:
cout << "bad quadrant" << endl;
@ -316,27 +321,47 @@ int *gethhx()
#ifndef MYROOT1
ex=(etax-etamin)/etastep;
ey=(etay-etamin)/etastep;
if (ex<0) ex=0;
if (ex>=nbeta) ex=nbeta-1;
if (ey<0) ey=0;
if (ey>=nbeta) ey=nbeta-1;
if (ex<0) {
ex=0;
cout << "x*"<< ex << endl;
}
if (ex>=nbeta) {
ex=nbeta-1;
cout << "x?"<< ex << endl;
}
if (ey<0) {
ey=0;
cout << "y*"<< ey << endl;
}
if (ey>=nbeta) {
ey=nbeta-1;
cout << "y?"<< ey << endl;
}
xpos_eta=(((double)hhx[(ey*nbeta+ex)]))/((double)nSubPixels);
ypos_eta=(((double)hhy[(ey*nbeta+ex)]))/((double)nSubPixels);
xpos_eta=(((double)hhx[(ey*nbeta+ex)]))+dX ;///((double)nSubPixels);
ypos_eta=(((double)hhy[(ey*nbeta+ex)]))+dY ;///((double)nSubPixels);
//else
//return 0;
#endif
} else {
xpos_eta=-dX*0.25;
ypos_eta=-0.25*dY;
xpos_eta=0.5*dX+0.25;
ypos_eta=0.5*dY+0.25;
}
int_x=((double)x) + 0.5*dX + xpos_eta;
int_y=((double)y) + 0.5*dY + ypos_eta;
// cout << "***"<< x <<" " << y << " " << int_x << " " << int_y << endl;
// cout << etax << " " << ex << " " << etay << " " << ey << " " << xpos_eta << " " << int_x << " " << ypos_eta << " " << int_y << endl;
// int_x=((double)x) + 0.5*dX + xpos_eta;
// int_y=((double)y) + 0.5*dY + ypos_eta;
int_x=((double)x) + xpos_eta;
int_y=((double)y) + ypos_eta;
/* if (int_x<x-0.5 || int_y<y-0.5 ) { */
/* // cout << "***"<< x <<" " << y << " " << int_x << " " << int_y << endl; */
/* cout << corner << " X "<< x << " " << etax << " " << xpos_eta << " " << int_x <<" Y "<< y << " " << etay << " " << ypos_eta << " " <<int_y << endl; */
/* } */
//return 1;
}
@ -480,8 +505,8 @@ int *gethhx()
#endif
#ifndef MYROOT1
int *heta;
int *hhx;
int *hhy;
float *hhx;
float *hhy;
#endif
int nbeta;
double etamin, etamax, etastep;

0
slsDetectorCalibration/etaVEL/etaInterpolationGlobal.h → slsDetectorCalibration/interpolations/etaInterpolationGlobal.h
View File

30
slsDetectorCalibration/etaVEL/etaInterpolationPosXY.h → slsDetectorCalibration/interpolations/etaInterpolationPosXY.h
View File

@ -71,26 +71,48 @@ class etaInterpolationPosXY : public etaInterpolationBase{
ii=0;
for (int ibx=0; ibx<nbeta; ibx++) {
if (hix[ibx]>(ii+1)*tot_eta_x*bsize) ii++;
if (tot_eta_x==0) {
hhx[ibx+ib*nbeta]=((float)ibx)/((float)nbeta);
ii=(ibx)/nbeta;
} else //if (hix[ibx]>(ii+1)*tot_eta_x*bsize)
{
//ii++;
// cout << ib << " x " << ibx << " " << tot_eta_x << " " << (ii)*tot_eta_x*bsize << " " << ii << endl;
// }
#ifdef MYROOT1
hhx->SetBinContent(ibx+1,ib+1,ii);
#endif
#ifndef MYROOT1
hhx[ibx+ib*nbeta]=ii;
hhx[ibx+ib*nbeta]=hix[ibx]/((float)tot_eta_x);//ii;
#endif
}
}
/* if (ii!=(nSubPixels-1)) */
/* cout << ib << " x " << tot_eta_x << " " << (ii+1)*tot_eta_x*bsize << " " << ii << " " << hix[nbeta-1]<< endl; */
ii=0;
for (int ibx=0; ibx<nbeta; ibx++) {
if (hiy[ibx]>(ii+1)*tot_eta_y*bsize) ii++;
if (tot_eta_y==0) {
hhx[ibx+ib*nbeta]=((float)ibx)/((float)nbeta);
ii=(ibx*nSubPixels)/nbeta;
} else //if (hiy[ibx]>(ii+1)*tot_eta_y*bsize)
{
//ii++;
//cout << ib << " y " << ibx << " " << tot_eta_y << " "<< (ii)*tot_eta_y*bsize << " " << ii << endl;
//}
#ifdef MYROOT1
hhy->SetBinContent(ib+1,ibx+1,ii);
#endif
#ifndef MYROOT1
hhy[ib+ibx*nbeta]=ii;
hhy[ib+ibx*nbeta]=hiy[ibx]/((float)tot_eta_y);//ii;
#endif
}
}
/* if (ii!=(nSubPixels-1)) */
/* cout << ib << " y " << tot_eta_y << " " << (ii+1)*tot_eta_y*bsize << " " << ii << " " << hiy[nbeta-1]<< endl; */
// cout << "y " << nbeta << " " << (ii+1)*tot_eta_x*bsize << " " << ii << endl;

0
slsDetectorCalibration/etaVEL/EVELAlg.C → slsDetectorCalibration/interpolations/etaVEL/EVELAlg.C
View File

0
slsDetectorCalibration/etaVEL/EtaVEL.cpp → slsDetectorCalibration/interpolations/etaVEL/EtaVEL.cpp
View File

0
slsDetectorCalibration/etaVEL/EtaVEL.h → slsDetectorCalibration/interpolations/etaVEL/EtaVEL.h
View File

0
slsDetectorCalibration/etaVEL/EtaVELTr.py → slsDetectorCalibration/interpolations/etaVEL/EtaVELTr.py
View File

0
slsDetectorCalibration/etaVEL/interpolation_etaVEL.cpp → slsDetectorCalibration/interpolations/etaVEL/interpolation_etaVEL.cpp
View File

0
slsDetectorCalibration/etaVEL/interpolation_EtaVEL.h → slsDetectorCalibration/interpolations/interpolation_EtaVEL.h
View File

12
slsDetectorCalibration/etaVEL/linearInterpolation.h → slsDetectorCalibration/interpolations/linearInterpolation.h
View File

@ -98,18 +98,18 @@ class linearInterpolation : public slsInterpolation{
{
case TOP_LEFT:
dX=-1.;
dY=+1.;
dY=0;
break;
case TOP_RIGHT:
dX=+1.;
dY=+1.;
dX=0;
dY=0;
break;
case BOTTOM_LEFT:
dX=-1.;
dY=-1.;
break;
case BOTTOM_RIGHT:
dX=+1.;
dX=0;
dY=-1.;
break;
default:
@ -125,8 +125,8 @@ class linearInterpolation : public slsInterpolation{
xpos_eta=0;
xpos_eta=0;
}
int_x=((double)x) + 0.5*dX + xpos_eta;
int_y=((double)y) + 0.5*dY + ypos_eta;
int_x=((double)x) + dX + xpos_eta;
int_y=((double)y) + dY + ypos_eta;
return;
};

8
slsDetectorCalibration/etaVEL/noInterpolation.h → slsDetectorCalibration/interpolations/noInterpolation.h
View File

@ -40,11 +40,11 @@ class noInterpolation : public slsInterpolation{
};
virtual void getInterpolatedPosition(int x, int y, double etax, double etay, int corner, double &int_x, double &int_y)
{
return getInterpolatedPosition(x, y, NULL, int_x, int_y);
getInterpolatedPosition(x, y, NULL, int_x, int_y);
};
virtual void getInterpolatedPosition(int x, int y, double totquad,int quad,double *cl,double &etax, double &etay){
getInterpolatedPosition(x, y, NULL, etax, etay);
};
virtual int getInterpolatedPosition(int x, int y, double totquad,int quad,double *cl,double &etax, double &etay){
return getInterpolatedPosition(x, y, NULL, int_x, int_y);
}
//////////////////////////////////////////////////////////////////////////////////////
virtual void getPositionETA3(int x, int y, double *data, double &int_x, double &int_y)

3
slsDetectorCalibration/etaVEL/slsInterpolation.h → slsDetectorCalibration/interpolations/slsInterpolation.h
View File

@ -7,6 +7,7 @@
#include <TH2F.h>
#endif
#include <cstdlib>
#include "tiffIO.h"
#ifndef DEF_QUAD
#define DEF_QUAD
@ -18,7 +19,9 @@
UNDEFINED_QUADRANT=-1
};
#endif
#include <memory.h>
using namespace std;
//#ifdef MYROOT1
//: public TObject
//#endif

264
slsDetectorCalibration/jungfrauReadData.C
View File

@ -1,264 +0,0 @@
#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 "jungfrau02Data.h"
#include "jungfrau02CommonMode.h"
#include "singlePhotonDetector.h"
//#include "MovingStat.h"
using namespace std;
#define NC 48
#define NR 48
#define MY_DEBUG 1
#ifdef MY_DEBUG
#include <TCanvas.h>
#endif
/**
Loops over data file to find single photons, fills the tree (and writes it to file, although the root file should be opened before) and creates 1x1, 2x2, 3x3 cluster histograms with ADCu on the x axis, channel number (48*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
\param hitfinder if 0: performs pedestal subtraction, not hit finding; if 1: performs both pedestal subtraction and hit finding
\returns pointer to histo stack with cluster spectra
*/
THStack *jungfrauReadData(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){
/*
// read in calibration file
ifstream calfile("/home/l_msdetect/TriesteBeam2014/dummy data for scripts/CalibrationParametersTest_fake.txt");
if (calfile.is_open()==0){cout << "Unable to open calibration file!" << endl;}
int pix;
double of0,sl0,of1,sl1,of2,sl2;
double of_0[NC*NR], of_1[NC*NR], of_2[NC*NR], sl_0[NC*NR], sl_1[NC*NR], sl_2[NC*NR];
while (calfile >> pix >> of0 >> sl0 >> of1 >> sl1 >> of2 >> sl2){
of_0[pix]=of0;
sl_0[pix]=sl0;
of_1[pix]=of1;
sl_1[pix]=sl1;
of_2[pix]=of2;
sl_2[pix]=sl2; //if(pix==200) cout << "sl_2[200] " << sl_2[200] << endl;
}
calfile.close();
*/
double adc_value, num_photon;
jungfrau02Data *decoder=new jungfrau02Data(1,0,0);//(1,0,0); // (adc,offset,crosstalk) //(1,0,0) //(3,0,0) for readout of GB
jungfrau02CommonMode *cmSub=NULL;
if (cmsub)
cmSub=new jungfrau02CommonMode();
int nph=0;
int iev=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; // not needed for JF?
TH2F *hetaY; // not needed for JF?
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); // not needed for JF?
hetaY=new TH2F("hetaY",tit,nbins,-1,2,NC*NR,-0.5,NC*NR-0.5); // not needed for JF?
hs->Add(h2);
hs->Add(h3);
hs->Add(hetaX); // not needed for JF?
hs->Add(hetaY); // not needed for JF?
}
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);
int gainBits=decoder->getGainBits(buff,ix,iy); //XXX
#ifdef MY_DEBUG
if (hitfinder) {
if (iev%1000==0)
//he->SetBinContent(ix+1-xmin, iy+1-ymin, (int)thisEvent); // show single photon hits // original
//he->SetBinContent(ix+1-xmin, iy+1-ymin, cmSub->getCommonMode(ix,iy)); //show common mode!
he->SetBinContent(iy+1-ymin, ix+1-xmin, (int)gainBits); //show gain bits
//he->SetBinContent(ix+1-xmin, iy+1-ymin, (int)gainBits); // rows and columns reversed!!!
}
#endif
if (nf>1000) { // only start filing data after 1000 frames
// h1->Fill(decoder->getValue(buff,ix,iy), iy+NR*ix);
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 << "="; // one "=" for every frame that's been processed
delete [] buff;
}
cout << nph << endl; // number of photons found in file
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;
}

10
slsDetectorCalibration/makefile
View File

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

1117
slsDetectorCalibration/moench03ReadData.C
View File

File diff suppressed because it is too large Load Diff

496
slsDetectorCalibration/moench03ReadData10Gb.C
View File

@ -1,496 +0,0 @@
#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 "moench03Ctb10GbT1Data.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) {
moench03Ctb10GbT1Data *decoder=new moench03Ctb10GbT1Data();
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];
moench03Ctb10GbT1Data *decoder=new moench03Ctb10GbT1Data();
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];
moench03Ctb10GbT1Data *decoder=new moench03Ctb10GbT1Data();
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];
moench03Ctb10GbT1Data *decoder=new moench03Ctb10GbT1Data();
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;
moench03Ctb10GbT1Data *decoder=new moench03Ctb10GbT1Data();
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=NULL;//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;
// 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=-1;
Int_t np=-1;
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;
if (filebin.is_open())
;
else
cout << "Could not open file " << fname<< endl;
while ((buff=decoder->readNextFrame(filebin, iFrame))) {
//cout << (void*)buff << endl;
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 << buff << " " << iFrame << " " << 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));
// }
iFrame=-1;
np=-1;
// }
}
//////////////////////////////////////////////////////////
#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];
moench03Ctb10GbT1Data *decoder=new moench03Ctb10GbT1Data();
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;
}

33
slsDetectorCalibration/moenchExecutables/Makefile.cluster_finder Normal file
View File

@ -0,0 +1,33 @@
#CBFLIBDIR= /afs/psi.ch/project/sls_det_software/CBFlib-0.9.5/
#ZMQLIB=../slsReceiverSoftware/include
#LIBRARYCBF=$(CBFLIBDIR)/lib/*.o
INCDIR=-IslsDetectorCalibration -I. -IetaVEL -I../slsReceiverSoftware/include
#-I$(CBFLIBDIR)/include/
#LIBHDF5=
LDFLAG= -L/usr/lib64/ -lpthread -lm -lstdc++ -L. -pthread -lrt -L$(CBFLIBDIR)/lib/ -ltiff
#-L$(ZMQLIB) -lzmq
#-L../../bin
MAIN=moench03ClusterFinder.cpp
#-lhdf5
#DESTDIR?=../bin
all: moenchClusterFinder moenchMakeEta moenchInterpolation moenchNoInterpolation
moenchClusterFinder: $(MAIN) $(INCS) clean
g++ -o moenchClusterFinder $(MAIN) $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF) tiffIO.cpp -DSAVE_ALL
moenchMakeEta: moench03MakeEta.cpp $(INCS) clean
g++ -o moenchMakeEta moench03MakeEta.cpp $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF) tiffIO.cpp -DSAVE_ALL
moenchInterpolation: moench03Interpolation.cpp $(INCS) clean
g++ -o moenchInterpolation moench03Interpolation.cpp $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF) tiffIO.cpp -DSAVE_ALL
moenchNoInterpolation: moench03NoInterpolation.cpp $(INCS) clean
g++ -o moenchNoInterpolation moench03NoInterpolation.cpp $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF) tiffIO.cpp -DSAVE_ALL
clean:
rm -f moenchClusterFinder moenchMakeEta

0
slsDetectorCalibration/Makefile.moench → slsDetectorCalibration/moenchExecutables/Makefile.moench
View File

13
slsDetectorCalibration/Makefile.moench_eta → slsDetectorCalibration/moenchExecutables/Makefile.moench_eta
View File

@ -2,21 +2,22 @@
CBFLIBDIR= /home/l_msdetect/CBFlib-0.9.5
ZMQLIB=../slsReceiverSoftware/include
LIBRARYCBF=$(CBFLIBDIR)/lib/*.o
INCDIR=-IslsDetectorCalibration -I../slsReceiverSoftware/include -I$(CBFLIBDIR)/include/ -I. -IetaVEL
INCDIR=-IslsDetectorCalibration -I../slsReceiverSoftware/include -I$(CBFLIBDIR)/include/ -I.
#-IetaVEL
LIBHDF5=
LDFLAG= -L/usr/lib64/ -lpthread -lm -lstdc++ -L. -lzmq -pthread -lrt -L$(CBFLIBDIR)/lib/ -lhdf5 -ltiff -L$(ZMQLIB)
#-L../../bin
MAIN=moench03ZmqOnTheFlyEta.cpp
MAIN=moench03ClusterFinder.cpp
#DESTDIR?=../bin
all: moench03ZmqOnTheFlyEta
all: moench03ClusterFinder
moench03ZmqOnTheFlyEta: $(MAIN) $(INCS) clean
g++ -o moench03ZmqOnTheFlyEta $(MAIN) $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF) tiffIO.cpp -DSAVE_ALL
moench03ClusterFinder: $(MAIN) $(INCS) clean
g++ -o moenchClusterFinder $(MAIN) $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF) tiffIO.cpp -DSAVE_ALL
clean:
rm -f moench03ZmqOnTheFlyEta
rm -f moench03ClusterFinder

0
slsDetectorCalibration/Makefile.moench_zmq → slsDetectorCalibration/moenchExecutables/Makefile.moench_zmq
View File

30
slsDetectorCalibration/moenchExecutables/Makefile.phoenix Normal file
View File

@ -0,0 +1,30 @@
CBFLIBDIR= /afs/psi.ch/project/sls_det_software/CBFlib-0.9.5/
#ZMQLIB=../slsReceiverSoftware/include
LIBRARYCBF=$(CBFLIBDIR)/lib/*.o
INCDIR=-IslsDetectorCalibration -I. -IetaVEL -I../slsReceiverSoftware/include -I$(CBFLIBDIR)/include/
#LIBHDF5=
LDFLAG= -L/usr/lib64/ -lpthread -lm -lstdc++ -L. -pthread -lrt -L$(CBFLIBDIR)/lib/ -ltiff -lhdf5
#-L$(ZMQLIB) -lzmq
#-L../../bin
MAIN=moench03ClusterFinderPhoenix.cpp
#
#DESTDIR?=../bin
all: moenchClusterFinderPhoenix
#moenchMakeEtaPhoenix moenchInterpolationPhoenix
moenchClusterFinderPhoenix: $(MAIN) $(INCS) clean
g++ -o moenchClusterFinderPhoenix $(MAIN) $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF) tiffIO.cpp -DSAVE_ALL
moenchMakeEtaPhoenix: moench03MakeEtaPhoenix.cpp $(INCS) clean
g++ -o moenchMakeEta moench03MakeEta.cpp $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF) tiffIO.cpp -DSAVE_ALL
moenchInterpolationPhoenix: moench03InterpolationPhoenix.cpp $(INCS) clean
g++ -o moenchInterpolation moench03Interpolation.cpp $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF) tiffIO.cpp -DSAVE_ALL
clean:
rm -f moenchClusterFinderPhoenix moenchMakeEtaPhoenix moenchInterpolationPhoenix

25
slsDetectorCalibration/moenchExecutables/Makefile.tiff_to_th2f Normal file
View File

@ -0,0 +1,25 @@
#INCSROOT= receiverGui.h
#INCS= $(INCSROOT) moench03_receiver.h
#LINKDEF=receiverGuiLinkDef.h
CBFLIBDIR= /afs/psi.ch/project/sls_det_software/CBFlib-0.9.5/
ZMQLIB=../slsReceiverSoftware/include
LIBRARYCBF=$(CBFLIBDIR)/lib/*.o
INCDIR=-IslsDetectorCalibration -I. -IetaVEL -I../slsReceiverSoftware/include -I$(CBFLIBDIR)/include/ -I$(ROOTSYS)/include
#LIBHDF5=
LDFLAG= -L/usr/lib64/ -lpthread -lm -lstdc++ -L. -pthread -lrt -L$(CBFLIBDIR)/lib/ -lhdf5 -ltiff -L$(ZMQLIB) -lzmq
MAIN=tiff_to_th2f.cpp
all: tiff_to_th2f
tiff_to_th2f: $(MAIN) $(INCS)
g++ -o tiff_to_th2f $(MAIN) `root-config --cflags --glibs` -lMinuit -lm -ltiff -lstdc++ $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF) tiffIO.cpp
clean:
rm -f tiff_to_th2f

255
slsDetectorCalibration/moenchExecutables/moench03ClusterFinder.cpp Normal file
View File

@ -0,0 +1,255 @@
//#include "ansi.h"
#include <iostream>
//#include "moench03T1ZmqData.h"
#include "moench03T1ReceiverData.h"
// #include "interpolatingDetector.h"
//#include "etaInterpolationPosXY.h"
// #include "linearInterpolation.h"
// #include "noInterpolation.h"
#include "multiThreadedAnalogDetector.h"
#include "singlePhotonDetector.h"
//#include "interpolatingDetector.h"
#include <stdio.h>
#include <map>
#include <fstream>
#include <sys/stat.h>
#include <ctime>
using namespace std;
int main(int argc, char *argv[]) {
if (argc<7) {
cout << "Usage is " << argv[0] << "indir outdir fname runmin runmax pedfname" << endl;
}
int p=10000;
int fifosize=1000;
int nthreads=20;
int nsubpix=25;
int etabins=nsubpix*10;
double etamin=-1, etamax=2;
int csize=3;
int nx=400, ny=400;
int save=1;
int nsigma=5;
int nped=1000;
int ndark=100;
int ok;
int iprog=0;
moench03T1ReceiverData *decoder=new moench03T1ReceiverData();
//moench03T1ZmqData *decoder=new moench03T1ZmqData();
singlePhotonDetector *filter=new singlePhotonDetector(decoder,csize, nsigma, 1, 0, nped, 200);
// char tit[10000];
// filter->readPedestals("/scratch/ped_100.tiff");
// interp->readFlatField("/scratch/eta_100.tiff",etamin,etamax);
// cout << "filter "<< endl;
int size = 327680;////atoi(argv[3]);
int* image;
//int* image =new int[327680/sizeof(int)];
filter->newDataSet();
int ff, np;
int dsize=decoder->getDataSize();
cout << " data size is " << dsize;
char data[dsize];
ifstream filebin;
char *indir=argv[1];
char *outdir=argv[2];
char *fformat=argv[3];
int runmin=atoi(argv[4]);
int runmax=atoi(argv[5]);
char fname[10000];
char outfname[10000];
char imgfname[10000];
char pedfname[10000];
strcpy(pedfname,argv[6]);
char fn[10000];
std::time_t end_time;
FILE *of=NULL;
cout << "input directory is " << indir << endl;
cout << "output directory is " << outdir << endl;
cout << "fileformat is " << fformat << endl;
cout << "pedestal file is " << fformat << endl;
std::time(&end_time);
cout << std::ctime(&end_time) << endl;
filter->setFrameMode(eFrame);
// mt->setFrameMode(ePedestal);
cout << pedfname<< endl;
filebin.open((const char *)(pedfname), ios::in | ios::binary);
//filebin.open((const char *)(fname), ios::in | ios::binary);
// //open file
if (filebin.is_open()){
// //while read frame
cout << "pedestal file " << endl;
while (decoder->readNextFrame(filebin, ff, np,data)) {
// cout << ff << " " << np << endl;
// //push
// mt->pushData(buff);
// // //pop
//mt->nextThread();
// // // cout << " " << (void*)buff;
//mt->popFree(buff);
filter->processData(data);
}
filebin.close();
// //close file
// //join threads
// while (mt->isBusy()) {;}//wait until all data are processed from the queues
// cout << outfname << endl;
// filter->writePedestals(outfname);
// sprintf(outfname,"%s/%s_pedimg.tiff",outdir,fn);
// cout << outfname << endl;
// filter->writeImage(outfname);
// //mt->clearImage();
} else
cout << "Could not open "<< pedfname << " for reading " << endl;
// for (int ix=0; ix<400; ix++)
// for (int iy=0; iy<400; iy++)
// cout << ix << " " << iy << " " << filter->getPedestal(ix,iy) << " " << filter->getPedestalRMS(ix,iy) << endl;
char* buff;
multiThreadedAnalogDetector *mt=new multiThreadedAnalogDetector(filter,nthreads,fifosize);
// mt->setFrameMode(eFrame); //need to find a way to switch between flat and frames!
// mt->prepareInterpolation(ok);
mt->setFrameMode(eFrame);
mt->StartThreads();
mt->popFree(buff);
int ifr=0;
// //loop on files
// mt->setFrameMode(eFrame);
//mt->setFrameMode(eFlat);
for (int irun=runmin; irun<runmin+5; irun++) {
sprintf(fn,fformat,irun);
sprintf(fname,"%s/%s.raw",indir,fn);
filebin.open((const char *)(fname), ios::in | ios::binary);
// //open file
if (filebin.is_open()){
// //while read frame
while (decoder->readNextFrame(filebin, ff, np,buff)) {
// cout << "*"<<ifr++<<"*"<<ff<< endl;
// cout << ff << " " << np << endl;
// //push
mt->pushData(buff);
// // //pop
mt->nextThread();
// // // cout << " " << (void*)buff;
mt->popFree(buff);
}
// cout << "--" << endl;
filebin.close();
}
}
while (mt->isBusy()) {;}//wait until all data are processed from the queues
mt->clearImage();
for (int irun=runmin; irun<runmax; irun++) {
sprintf(fn,fformat,irun);
sprintf(fname,"%s/%s.raw",indir,fn);
sprintf(outfname,"%s/%s.clust",outdir,fn);
sprintf(imgfname,"%s/%s.tiff",outdir,fn);
std::time(&end_time);
cout << std::ctime(&end_time) << endl;
cout << fname << " " << outfname << " " << imgfname << endl;
filebin.open((const char *)(fname), ios::in | ios::binary);
// //open file
if (filebin.is_open()){
of=fopen(outfname,"w");
if (of) {
mt->setFilePointer(of);
// cout << "file pointer set " << endl;
} else {
cout << "Could not open "<< outfname << " for writing " << endl;
mt->setFilePointer(NULL);
return 1;
}
// //while read frame
while (decoder->readNextFrame(filebin, ff, np,buff)) {
// cout << "*"<<ifr++<<"*"<<ff<< endl;
// cout << ff << " " << np << endl;
// //push
mt->pushData(buff);
// // //pop
mt->nextThread();
// // // cout << " " << (void*)buff;
mt->popFree(buff);
}
// cout << "--" << endl;
filebin.close();
// //close file
// //join threads
while (mt->isBusy()) {;}//wait until all data are processed from the queues
if (of)
fclose(of);
mt->writeImage(imgfname);
mt->clearImage();
std::time(&end_time);
cout << std::ctime(&end_time) << endl;
} else
cout << "Could not open "<< fname << " for reading " << endl;
}
return 0;
}

192
slsDetectorCalibration/moenchExecutables/moench03ClusterFinderPhoenix.cpp Normal file
View File

@ -0,0 +1,192 @@
//#include "ansi.h"
#include <iostream>
//#include "moench03T1ZmqData.h"
//#include "moench03T1ReceiverData.h"
#include "moench03Ctb10GbT1Data.h"
// #include "interpolatingDetector.h"
//#include "etaInterpolationPosXY.h"
// #include "linearInterpolation.h"
// #include "noInterpolation.h"
#include "multiThreadedAnalogDetector.h"
#include "singlePhotonDetector.h"
//#include "interpolatingDetector.h"
#include <stdio.h>
#include <map>
#include <fstream>
#include <sys/stat.h>
#include <ctime>
using namespace std;
int main(int argc, char *argv[]) {
if (argc<6) {
cout << "Usage is " << argv[0] << "indir outdir fname runmin runmax " << endl;
return 0;
}
int p=10000;
int fifosize=1000;
int nthreads=20;
int nsubpix=25;
int etabins=nsubpix*10;
double etamin=-1, etamax=2;
int csize=3;
int nx=400, ny=400;
int save=1;
int nsigma=5;
int nped=1000;
int ndark=100;
int ok;
int iprog=0;
moench03Ctb10GbT1Data *decoder=new moench03Ctb10GbT1Data();
cout << "decoder" << endl;
//moench03T1ReceiverData *decoder=new moench03T1ReceiverData();
//moench03T1ZmqData *decoder=new moench03T1ZmqData();
singlePhotonDetector *filter=new singlePhotonDetector(decoder,csize, nsigma, 1, 0, nped, 200);
// char tit[10000];
cout << "filter" << endl;
// filter->readPedestals("/scratch/ped_100.tiff");
// interp->readFlatField("/scratch/eta_100.tiff",etamin,etamax);
// cout << "filter "<< endl;
int size = 327680;////atoi(argv[3]);
int* image;
//int* image =new int[327680/sizeof(int)];
filter->newDataSet();
cout << "dataset" << endl;
int ff, np;
int dsize=decoder->getDataSize();
cout << " data size is " << dsize << endl;
char data[dsize];
ifstream filebin;
char *indir=argv[1];
cout << "input directory is " << indir << endl;
char *outdir=argv[2];
cout << "output directory is " << outdir << endl;
char *fformat=argv[3];
cout << "fileformat is " << fformat << endl;
int runmin=atoi(argv[4]);
cout << "runmin : " << runmin << endl;
int runmax=atoi(argv[5]);
cout << "runmax : " << runmax << endl;
char fname[10000];
char outfname[10000];
char imgfname[10000];
char pedfname[10000];
char fn[10000];
std::time_t end_time;
FILE *of=NULL;
cout << "time " << endl;
std::time(&end_time);
cout << std::ctime(&end_time) << endl;
filter->setFrameMode(eFrame);
// mt->setFrameMode(ePedestal);
// for (int ix=0; ix<400; ix++)
// for (int iy=0; iy<400; iy++)
// cout << ix << " " << iy << " " << filter->getPedestal(ix,iy) << " " << filter->getPedestalRMS(ix,iy) << endl;
char* buff;
cout << "aa " << endl;
multiThreadedAnalogDetector *mt=new multiThreadedAnalogDetector(filter,nthreads,fifosize);
cout << "mt " << endl;
// mt->setFrameMode(eFrame); //need to find a way to switch between flat and frames!
// mt->prepareInterpolation(ok);
mt->setFrameMode(eFrame);
mt->StartThreads();
mt->popFree(buff);
int ifr=0;
// //loop on files
// mt->setFrameMode(eFrame);
//mt->setFrameMode(eFlat);
for (int irun=runmin; irun<runmax; irun++) {
sprintf(fn,fformat,irun);
sprintf(fname,"%s/%s.raw",indir,fn);
sprintf(outfname,"%s/%s.clust",outdir,fn);
sprintf(imgfname,"%s/%s.tiff",outdir,fn);
std::time(&end_time);
cout << std::ctime(&end_time) << endl;
cout << fname << " " << outfname << " " << imgfname << endl;
filebin.open((const char *)(fname), ios::in | ios::binary);
// //open file
if (filebin.is_open()){
of=fopen(outfname,"w");
if (of) {
mt->setFilePointer(of);
// cout << "file pointer set " << endl;
} else {
cout << "Could not open "<< outfname << " for writing " << endl;
mt->setFilePointer(NULL);
return 1;
}
// //while read frame
ff=-1;
while (decoder->readNextFrame(filebin, ff, np,buff)) {
// cout << "*"<<ifr++<<"*"<<ff<< endl;
// cout << ff << " " << np << endl;
// //push
mt->pushData(buff);
// // //pop
mt->nextThread();
// // // cout << " " << (void*)buff;
mt->popFree(buff);
ff=-1;
}
// cout << "--" << endl;
filebin.close();
// //close file
// //join threads
while (mt->isBusy()) {;}//wait until all data are processed from the queues
if (of)
fclose(of);
mt->writeImage(imgfname);
mt->clearImage();
std::time(&end_time);
cout << std::ctime(&end_time) << endl;
} else
cout << "Could not open "<< fname << " for reading " << endl;
}
return 0;
}

171
slsDetectorCalibration/moenchExecutables/moench03Interpolation.cpp Normal file
View File

@ -0,0 +1,171 @@
#include "ansi.h"
#include <iostream>
#include "moench03T1ZmqData.h"
#include "single_photon_hit.h"
#include "etaInterpolationPosXY.h"
using namespace std;
#define NC 400
#define NR 400
int main(int argc, char *argv[]) {
/**
* trial.o [socket ip] [starting port number] [outfname]
*
*/
if (argc<7) {
cout << "Wrong usage! Should be: "<< argv[0] << " infile " << " etafile outfile runmin runmax ns" << endl;
return 1;
}
char infname[10000];
char outfname[10000];
int runmin=atoi(argv[4]);
int runmax=atoi(argv[5]);
int nsubpix=atoi(argv[6]);
int etabins=1000;//nsubpix*2*100;
double etamin=-1, etamax=2;
int quad;
double sum, totquad;
double sDum[2][2];
double etax, etay, int_x, int_y;
int ok;
int ix, iy, isx, isy;
FILE *f=NULL;
single_photon_hit cl(3,3);
etaInterpolationPosXY *interp=new etaInterpolationPosXY(NC, NR, nsubpix, etabins, etamin, etamax);
interp->readFlatField(argv[2]);
interp->prepareInterpolation(ok);
int *img;
float *totimg=new float[NC*NR*nsubpix*nsubpix];
for (ix=0; ix<NC; ix++) {
for (iy=0; iy<NR; iy++) {
for (isx=0; isx<nsubpix; isx++) {
for (isy=0; isy<nsubpix; isy++) {
totimg[ix*nsubpix+isx+(iy*nsubpix+isy)*(NC*nsubpix)]=0;
}
}
}
}
#ifdef FF
float ff[nsubpix*nsubpix];
float *ffimg=new float[NC*NR*nsubpix*nsubpix];
float totff=0;
#endif
for (int irun=runmin; irun<runmax; irun++) {
sprintf(infname,argv[1],irun);
sprintf(outfname,argv[3],irun);
f=fopen(infname,"r");
if (f) {
while (cl.read(f)) {
quad=interp->calcQuad(cl.get_cluster(), sum, totquad, sDum);
if (sum>200 && sum<580) {
interp->getInterpolatedPosition(cl.x,cl.y, totquad,quad,cl.get_cluster(),int_x, int_y);
interp->addToImage(int_x, int_y);
}
}
fclose(f);
#ifdef FF
img=interp->getInterpolatedImage();
for (isx=0; isx<nsubpix; isx++) {
for (isy=0; isy<nsubpix; isy++) {
ff[isy*nsubpix+isx]=0;
}
}
totff=0;
for (ix=0; ix<NC; ix++) {
for (iy=0; iy<NR-100; iy++) {
for (isx=0; isx<nsubpix; isx++) {
for (isy=0; isy<nsubpix; isy++) {
ff[isy*nsubpix+isx]+=img[ix*nsubpix+isx+(iy*nsubpix+isy)*(NC*nsubpix)];
}
}
}
}
for (isx=0; isx<nsubpix; isx++) {
for (isy=0; isy<nsubpix; isy++) {
totff+=ff[isy*nsubpix+isx];
}
}
totff/=nsubpix*nsubpix;
if (totff) {
cout << "ff: " << totff << endl;
for (isx=0; isx<nsubpix; isx++) {
for (isy=0; isy<nsubpix; isy++) {
ff[isy*nsubpix+isx]/=totff;
cout << ff[isy*nsubpix+isx] << "\t";
}
cout << endl;
}
for (ix=0; ix<NC; ix++) {
for (iy=0; iy<NR; iy++) {
for (isx=0; isx<nsubpix; isx++) {
for (isy=0; isy<nsubpix; isy++) {
ffimg[ix*nsubpix+isx+(iy*nsubpix+isy)*(NC*nsubpix)]=img[ix*nsubpix+isx+(iy*nsubpix+isy)*(NC*nsubpix)]/ff[isy*nsubpix+isx];
totimg[ix*nsubpix+isx+(iy*nsubpix+isy)*(NC*nsubpix)]+=ffimg[ix*nsubpix+isx+(iy*nsubpix+isy)*(NC*nsubpix)];
}
}
}
}
WriteToTiff(ffimg, outfname,NC*nsubpix,NR*nsubpix);
}
#endif
#ifndef FF
interp->writeInterpolatedImage(outfname);
img=interp->getInterpolatedImage();
for (ix=0; ix<NC; ix++) {
for (iy=0; iy<NR; iy++) {
for (isx=0; isx<nsubpix; isx++) {
for (isy=0; isy<nsubpix; isy++) {
totimg[ix*nsubpix+isx+(iy*nsubpix+isy)*(NC*nsubpix)]+=img[ix*nsubpix+isx+(iy*nsubpix+isy)*(NC*nsubpix)];
}
}
}
}
#endif
interp->clearInterpolatedImage();
}
}
sprintf(outfname,argv[3],11111);
WriteToTiff(totimg, outfname,NC*nsubpix,NR*nsubpix);
return 0;
}

62
slsDetectorCalibration/moenchExecutables/moench03MakeEta.cpp Normal file
View File

@ -0,0 +1,62 @@
#include "ansi.h"
#include <iostream>
#include "moench03T1ZmqData.h"
#include "single_photon_hit.h"
#include "etaInterpolationPosXY.h"
using namespace std;
#define NC 400
#define NR 400
int main(int argc, char *argv[]) {
/**
* trial.o [socket ip] [starting port number] [outfname]
*
*/
int nsubpix=10;
int etabins=nsubpix*100;
double etamin=-1, etamax=2;
int quad;
double sum, totquad;
double sDum[2][2];
char fname[10000];
double etax, etay;
int runmin, runmax;
single_photon_hit cl(3,3);
if (argc<5) {
cout << "Wrong usage! Should be: "<< argv[0] << " infile " << " outfile runmin runmax" << endl;
return 1;
}
etaInterpolationPosXY *interp=new etaInterpolationPosXY(NR, NC, nsubpix, etabins, etamin, etamax);
runmin=atoi(argv[3]);
runmax=atoi(argv[4]);
FILE *f;
for (int i=runmin; i<runmax; i++) {
sprintf(fname,argv[1],i);
f=fopen(fname,"r");
if (f) {
cout << "*" << endl;
while (cl.read(f)) {
interp->calcQuad(cl.get_cluster(), sum, totquad, sDum);
if (sum>200 && sum<580 && cl.y<350)
interp->addToFlatField(cl.get_cluster(),etax, etay);
}
fclose(f);
interp->writeFlatField(argv[2]);
}
else cout << "could not open file " << fname << endl;
}
interp->writeFlatField(argv[2]);
return 0;
}

171
slsDetectorCalibration/moenchExecutables/moench03NoInterpolation.cpp Normal file
View File

@ -0,0 +1,171 @@
#include "ansi.h"
#include <iostream>
#include "moench03T1ZmqData.h"
#include "single_photon_hit.h"
#include "noInterpolation.h"
using namespace std;
#define NC 400
#define NR 400
int main(int argc, char *argv[]) {
/**
* trial.o [socket ip] [starting port number] [outfname]
*
*/
if (argc<7) {
cout << "Wrong usage! Should be: "<< argv[0] << " infile " << " etafile outfile runmin runmax ns" << endl;
return 1;
}
char infname[10000];
char outfname[10000];
int runmin=atoi(argv[4]);
int runmax=atoi(argv[5]);
int nsubpix=atoi(argv[6]);
int etabins=1000;//nsubpix*2*100;
double etamin=-1, etamax=2;
int quad;
double sum, totquad;
double sDum[2][2];
double etax, etay, int_x, int_y;
int ok;
int ix, iy, isx, isy;
FILE *f=NULL;
single_photon_hit cl(3,3);
// etaInterpolationPosXY *interp=new etaInterpolationPosXY(NC, NR, nsubpix, etabins, etamin, etamax);
noInterpolation *interp=new noInterpolation(NC, NR, nsubpix);
// interp->readFlatField(argv[2]);
// interp->prepareInterpolation(ok);
int *img;
float *totimg=new float[NC*NR*nsubpix*nsubpix];
for (ix=0; ix<NC; ix++) {
for (iy=0; iy<NR; iy++) {
for (isx=0; isx<nsubpix; isx++) {
for (isy=0; isy<nsubpix; isy++) {
totimg[ix*nsubpix+isx+(iy*nsubpix+isy)*(NC*nsubpix)]=0;
}
}
}
}
#ifdef FF
float ff[nsubpix*nsubpix];
float *ffimg=new float[NC*NR*nsubpix*nsubpix];
float totff=0;
#endif
for (int irun=runmin; irun<runmax; irun++) {
sprintf(infname,argv[1],irun);
sprintf(outfname,argv[3],irun);
f=fopen(infname,"r");
if (f) {
while (cl.read(f)) {
quad=interp->calcQuad(cl.get_cluster(), sum, totquad, sDum);
if (sum>200 && sum<580) {
interp->getInterpolatedPosition(cl.x,cl.y, totquad,quad,cl.get_cluster(),int_x, int_y);
interp->addToImage(int_x, int_y);
}
}
fclose(f);
#ifdef FF
img=interp->getInterpolatedImage();
for (isx=0; isx<nsubpix; isx++) {
for (isy=0; isy<nsubpix; isy++) {
ff[isy*nsubpix+isx]=0;
}
}
totff=0;
for (ix=0; ix<NC; ix++) {
for (iy=0; iy<NR-100; iy++) {
for (isx=0; isx<nsubpix; isx++) {
for (isy=0; isy<nsubpix; isy++) {
ff[isy*nsubpix+isx]+=img[ix*nsubpix+isx+(iy*nsubpix+isy)*(NC*nsubpix)];
}
}
}
}
for (isx=0; isx<nsubpix; isx++) {
for (isy=0; isy<nsubpix; isy++) {
totff+=ff[isy*nsubpix+isx];
}
}
totff/=nsubpix*nsubpix;
if (totff) {
cout << "ff: " << totff << endl;
for (isx=0; isx<nsubpix; isx++) {
for (isy=0; isy<nsubpix; isy++) {
ff[isy*nsubpix+isx]/=totff;
cout << ff[isy*nsubpix+isx] << "\t";
}
cout << endl;
}
for (ix=0; ix<NC; ix++) {
for (iy=0; iy<NR; iy++) {
for (isx=0; isx<nsubpix; isx++) {
for (isy=0; isy<nsubpix; isy++) {
ffimg[ix*nsubpix+isx+(iy*nsubpix+isy)*(NC*nsubpix)]=img[ix*nsubpix+isx+(iy*nsubpix+isy)*(NC*nsubpix)]/ff[isy*nsubpix+isx];
totimg[ix*nsubpix+isx+(iy*nsubpix+isy)*(NC*nsubpix)]+=ffimg[ix*nsubpix+isx+(iy*nsubpix+isy)*(NC*nsubpix)];
}
}
}
}
WriteToTiff(ffimg, outfname,NC*nsubpix,NR*nsubpix);
}
#endif
#ifndef FF
interp->writeInterpolatedImage(outfname);
img=interp->getInterpolatedImage();
for (ix=0; ix<NC; ix++) {
for (iy=0; iy<NR; iy++) {
for (isx=0; isx<nsubpix; isx++) {
for (isy=0; isy<nsubpix; isy++) {
totimg[ix*nsubpix+isx+(iy*nsubpix+isy)*(NC*nsubpix)]+=img[ix*nsubpix+isx+(iy*nsubpix+isy)*(NC*nsubpix)];
}
}
}
}
#endif
interp->clearInterpolatedImage();
}
}
sprintf(outfname,argv[3],11111);
WriteToTiff(totimg, outfname,NC*nsubpix,NR*nsubpix);
return 0;
}

0
slsDetectorCalibration/moench03OnTheFlyAnalysis.C → slsDetectorCalibration/moenchExecutables/moench03OnTheFlyAnalysis.C
View File

272
slsDetectorCalibration/moenchExecutables/moench03ZmqOnTheFly.cpp Normal file
View File

@ -0,0 +1,272 @@
#include "ZmqSocket.h"
#include "ansi.h"
#include <iostream>
#include "moench03T1ZmqData.h"
// #include "interpolatingDetector.h"
//#include "etaInterpolationPosXY.h"
// #include "linearInterpolation.h"
// #include "noInterpolation.h"
#include "multiThreadedAnalogDetector.h"
#include "singlePhotonDetector.h"
#include "interpolatingDetector.h"
using namespace std;
#define NC 400
#define NR 400
#define SLS_DETECTOR_JSON_HEADER_VERSION 0x2
int main(int argc, char *argv[]) {
/**
* trial.o [socket ip] [starting port number] [outfname]
*
*/
int p=10000;
int fifosize=1000;
int nthreads=20;
int nsubpix=25;
int etabins=nsubpix*10;
double etamin=-1, etamax=2;
int csize=3;
int nx=400, ny=400;
int save=1;
int nsigma=5;
int nped=1000;
int ndark=100;
int ok;
int iprog=0;
char outfname[10000];
moench03T1ZmqData *decoder=new moench03T1ZmqData();
// cout << "decoder "<< endl;
//etaInterpolationPosXY *interp=new etaInterpolationPosXY(nx, ny, nsubpix, etabins, etamin, etamax);
// linearInterpolation *interp=new linearInterpolation(NC, NR, nsubpix);
//noInterpolation *interp=new noInterpolation(NC, NR, nsubpix);
//interpolatingDetector *filter=new interpolatingDetector(decoder,interp, nsigma, 1, 0, nped, 100);
singlePhotonDetector *filter=new singlePhotonDetector(decoder,csize, nsigma, 1, 0, nped, 100);
char tit[10000];
// filter->readPedestals("/scratch/ped_100.tiff");
// interp->readFlatField("/scratch/eta_100.tiff",etamin,etamax);
// cout << "filter "<< endl;
int size = 327680;////atoi(argv[3]);
char* buff;
int* image;
//int* image =new int[327680/sizeof(int)];
filter->newDataSet();
multiThreadedAnalogDetector *mt=new multiThreadedAnalogDetector(filter,nthreads,fifosize);
//mt->setFrameMode(eFrame);
mt->setFrameMode(eFlat);
// mt->setFrameMode(eFrame); //need to find a way to switch between flat and frames!
// mt->prepareInterpolation(ok);
mt->StartThreads();
mt->popFree(buff);
// help
if (argc < 2 || (argc > 2)) {
cprintf(RED, "Help: %s [receive socket ip]\n", argv[0]);
return EXIT_FAILURE;
}
// receive parameters
bool send = false;
char* socketip = argv[1];
uint32_t portnum = 30001;//atoi(argv[2]);
// char fn[10000];
// strcpy(fn, argv[2]);
// send parameters if any
// char* socketip2 = 0;
// uint32_t portnum2 = 0;
// if (argc > 4) {
// send = true;
// socketip2 = argv[4];
// portnum2 = atoi(argv[5]);
// }
cout << "\nrx socket ip : " << socketip ; // <<
// "\nrx port num : " << portnum <<
// "\nsize : " << size;
// "\nfname : " << fn ;
// if (send) {
// cout << "\nsd socket ip : " << socketip2 <<
// "\nsd port num : " << portnum2;
// }
cout << endl;
// receive socket
ZmqSocket* zmqsocket = new ZmqSocket(socketip,portnum);
if (zmqsocket->IsError()) {
cprintf(RED, "Error: Could not create Zmq socket on port %d with ip %s\n", portnum, socketip);
delete zmqsocket;
return EXIT_FAILURE;
}
printf("Zmq Client at %s\n", zmqsocket->GetZmqServerAddress());
// send socket
// ZmqSocket* zmqsocket2 = 0;
// if (send) {
// zmqsocket2 = new ZmqSocket(portnum2, socketip2);
// if (zmqsocket2->IsError()) {
// bprintf(RED, "Error: Could not create Zmq socket server on port %d and ip %s\n", portnum2, socketip2);
// delete zmqsocket2;
// delete zmqsocket;
// return EXIT_FAILURE;
// }
// printf("Zmq Server started at %s\n", zmqsocket2->GetZmqServerAddress());
// }
// cout << "size " << 327680/sizeof(int) << endl;
// header variables
uint64_t acqIndex = -1;
uint64_t frameIndex = -1;
uint32_t subframeIndex = -1;
string filename = "";
int nnx, nny, nns;
int imsize=filter->getImageSize(nnx,nny,nns);
//int imsize=nx*ny;
int i_image=0;
cout << "Image size is "<< nnx << " " << nny << " " << nns << " " << imsize << endl;
int iframe=0;
int ff=0;
// infinite loop
while(1) {
// cprintf(GREEN, "Got ?\n");
// get header, (if dummy, fail is on parse error or end of acquisition)
if (!zmqsocket->ReceiveHeader(0, acqIndex, frameIndex, subframeIndex, filename)) {
cprintf(RED, "Acquisition finished\n");
cout << "Received " << ff << " frames for a total of "<< iframe << endl;
while (mt->isBusy()) {;}//wait until all data are processed from the queues
// stream dummy to socket2 to signal end of acquisition
// if (send) {
// zmqsocket2->SendData((char*)(mt->getImage()),imsize*sizeof(int));
// cprintf(BLUE, "Sent Interpolated image\n");
// // zmqsocket2->SendHeaderData(0, true, SLS_DETECTOR_JSON_HEADER_VERSION);
// // cprintf(RED, "Sent Dummy\n");
// }
// if (save) {
sprintf(tit,"%s_%05d.tiff",filename.c_str(), ff);
// cout << tit << endl;
mt->writeImage(tit);
mt->clearImage();
// }
// sprintf(tit,"%s_%05d_eta_i.tiff",filename.c_str(), ff);
// interp->writeFlatField(tit);
// sprintf(tit,"%s_%05d_eta.tiff",filename.c_str(), ff);
// cout << tit << endl;
// mt->writeFlatField(tit);
i_image++;
ff=0;
// dont get data
continue; //continue to not get out
}
if (ff==0)
cprintf(GREEN, "Start acquisition \n");
// cout << filename << endl;
// cprintf(GREEN, "Got Header \n");
// get data
// cprintf(GREEN, "Got Header\n");
//image=(int*)buff;
//cout << buff << endl;
int length = zmqsocket->ReceiveData(0, (int*)buff, size);
// int length = zmqsocket->ReceiveData(0, (int*)image, size);
// cprintf(GREEN, "Got Data\n");
//processing with image
//...
mt->pushData(buff);
mt->nextThread();
// cout << " " << (void*)buff;
mt->popFree(buff);
iframe++;
ff++;
// if (iframe%p==0) {
// while (mt->isBusy()) {;}//wait until all data are processed from the queues
// sprintf(tit,"tmp.tiff",filename.c_str(), i_image);
// mt->writeImage(tit);
// // mt->clearImage();
// // }
// //interp->writeFlatField(tit);
// // i_image++;
// // sprintf(tit,"%s_tmp.tiff",fn);
// // mt->writeImage(tit);
// // //mt->clearImage();
// // cout <<"*"<< iprog++ << endl;
// }
// //stream data from socket 2
// if (send) {
// zmqsocket2->SendHeaderData(0, false, SLS_DETECTOR_JSON_HEADER_VERSION,
// 0,0,0,acqIndex,frameIndex,(char*)"run", acqIndex, 0,0,0,0,0,0,0,0,0,0,0,1);
// cprintf(GREEN, "Sent Header\n");
// zmqsocket2->SendData((char*)image,length);
// cprintf(GREEN, "Sent Data\n");
// }
}// exiting infinite loop
delete zmqsocket;
// if (send)
// delete zmqsocket2;
// cout<<"Goodbye"<< endl;
return 0;
}

306
slsDetectorCalibration/moenchExecutables/moench03ZmqOnTheFlyEta.cpp Normal file
View File

@ -0,0 +1,306 @@
#include "ZmqSocket.h"
#include "ansi.h"
#include <iostream>
#include "moench03T1ZmqData.h"
// #include "interpolatingDetector.h"
#include "etaInterpolationPosXY.h"
// #include "linearInterpolation.h"
// #include "noInterpolation.h"
#include "multiThreadedAnalogDetector.h"
#include "singlePhotonDetector.h"
#include "interpolatingDetector.h"
using namespace std;
#define NC 400
#define NR 400
#define SLS_DETECTOR_JSON_HEADER_VERSION 0x2
int main(int argc, char *argv[]) {
/**
* trial.o [socket ip] [starting port number] [outfname]
*
*/
int p=1000; //1.5s with 15ms period
int fifosize=1000;
int nthreads=8;
int nsubpix=2;
int etabins=nsubpix*10;
double etamin=-1, etamax=2;
int csize=3;
int nx=400, ny=400;
int save=1;
int nsigma=5;
int nped=1000;
int ndark=100;
int ok;
int iprog=0;
char outfname[10000];
if (nsubpix>2) p=5000;
moench03T1ZmqData *decoder=new moench03T1ZmqData();
// cout << "decoder "<< endl;
etaInterpolationPosXY *interp=new etaInterpolationPosXY(nx, ny, nsubpix, etabins, etamin, etamax);
// linearInterpolation *interp=new linearInterpolation(NC, NR, nsubpix);
//noInterpolation *interp=new noInterpolation(NC, NR, nsubpix);
interpolatingDetector *filter=new interpolatingDetector(decoder,interp, nsigma, 1, 0, nped, 100);
// singlePhotonDetector *filter=new singlePhotonDetector(decoder,csize, nsigma, 1, 0, nped, 100);
char tit[10000];
// filter->readPedestals("/scratch/ped_100.tiff");
// interp->readFlatField("/scratch/eta_100.tiff",etamin,etamax);
// cout << "filter "<< endl;
int size = 327680;////atoi(argv[3]);
char* buff;
int* image;
//int* image =new int[327680/sizeof(int)];
filter->newDataSet();
multiThreadedAnalogDetector *mt=new multiThreadedAnalogDetector(filter,nthreads,fifosize);
//
if (nsubpix>2)
mt->setFrameMode(eFlat);
else
mt->setFrameMode(eFrame);
// mt->setFrameMode(eFrame); //need to find a way to switch between flat and frames!
// mt->prepareInterpolation(ok);
mt->StartThreads();
mt->popFree(buff);
// help
if (argc < 2 || (argc > 2)) {
cprintf(RED, "Help: %s [receive socket ip]\n", argv[0]);
return EXIT_FAILURE;
}
// receive parameters
bool send = false;
char* socketip = argv[1];
uint32_t portnum = 30001;//atoi(argv[2]);
// char fn[10000];
// strcpy(fn, argv[2]);
// send parameters if any
// char* socketip2 = 0;
// uint32_t portnum2 = 0;
// if (argc > 4) {
// send = true;
// socketip2 = argv[4];
// portnum2 = atoi(argv[5]);
// }
cout << "\nrx socket ip : " << socketip ; // <<
// "\nrx port num : " << portnum <<
// "\nsize : " << size;
// "\nfname : " << fn ;
// if (send) {
// cout << "\nsd socket ip : " << socketip2 <<
// "\nsd port num : " << portnum2;
// }
cout << endl;
// receive socket
ZmqSocket* zmqsocket = new ZmqSocket(socketip,portnum);
if (zmqsocket->IsError()) {
cprintf(RED, "Error: Could not create Zmq socket on port %d with ip %s\n", portnum, socketip);
delete zmqsocket;
return EXIT_FAILURE;
}
printf("Zmq Client at %s\n", zmqsocket->GetZmqServerAddress());
// send socket
// ZmqSocket* zmqsocket2 = 0;
// if (send) {
// zmqsocket2 = new ZmqSocket(portnum2, socketip2);
// if (zmqsocket2->IsError()) {
// bprintf(RED, "Error: Could not create Zmq socket server on port %d and ip %s\n", portnum2, socketip2);
// delete zmqsocket2;
// delete zmqsocket;
// return EXIT_FAILURE;
// }
// printf("Zmq Server started at %s\n", zmqsocket2->GetZmqServerAddress());
// }
// cout << "size " << 327680/sizeof(int) << endl;
// header variables
uint64_t acqIndex = -1;
uint64_t frameIndex = -1;
uint32_t subframeIndex = -1;
string filename = "";
int nnx, nny, nns;
int imsize=filter->getImageSize(nnx,nny,nns);
//int imsize=nx*ny;
int i_image=0;
cout << "Image size is "<< nnx << " " << nny << " " << nns << " " << imsize << endl;
int iframe=0;
int ff=0;
// infinite loop
while(1) {
// cprintf(GREEN, "Got ?\n");
// get header, (if dummy, fail is on parse error or end of acquisition)
if (!zmqsocket->ReceiveHeader(0, acqIndex, frameIndex, subframeIndex, filename)) {
cprintf(RED, "Acquisition finished\n");
cout << "Received " << ff << " frames for a total of "<< iframe << endl;
while (mt->isBusy()) {;}//wait until all data are processed from the queues
// stream dummy to socket2 to signal end of acquisition
// if (send) {
// zmqsocket2->SendData((char*)(mt->getImage()),imsize*sizeof(int));
// cprintf(BLUE, "Sent Interpolated image\n");
// // zmqsocket2->SendHeaderData(0, true, SLS_DETECTOR_JSON_HEADER_VERSION);
// // cprintf(RED, "Sent Dummy\n");
// }
// if (save) {
// if (iframe>p) {
sprintf(tit,"%s_%05d_%05d.tiff",filename.c_str(), ff, i_image);
cout << tit << endl;
mt->writeImage(tit);
sprintf(tit,"/home/l_msdetect/slsbl/x07mb/x07mbop/Data1/2017/11/06/tmp.tiff");
cout <<"*"<< iframe<< endl;
mt->writeImage(tit);
mt->clearImage();
//}
// }
//sprintf(tit,"%s_%05d_eta_i.tiff",filename.c_str(), ff);
//interp->writeFlatField(tit);
// sprintf(tit,"%s_%05d_eta_%d.tiff",filename.c_str(), ff, acqIndex);
// // cout << tit << endl;
// mt->writeFlatField(tit);
i_image++;
ff=0;
//filter->newDataSet();
// dont get data
continue; //continue to not get out
}
if (ff==0)
cprintf(GREEN, "Start acquisition %d\n",acqIndex);
// cout << filename << endl;
// cprintf(GREEN, "Got Header \n");
// get data
// cprintf(GREEN, "Got Header\n");
//image=(int*)buff;
//cout << buff << endl;
int length = zmqsocket->ReceiveData(0, (int*)buff, size);
// int length = zmqsocket->ReceiveData(0, (int*)image, size);
// cprintf(GREEN, "Got Data\n");
//processing with image
//...
mt->pushData(buff);
mt->nextThread();
// cout << " " << (void*)buff;
mt->popFree(buff);
iframe++;
if (nsubpix>2) {
if (iframe==p) {
mt->setFrameMode(eFrame);
mt->prepareInterpolation(ok);
sprintf(tit,"%s_%05d_eta.tiff",filename.c_str(), ff);
cout << tit << endl;
mt->writeFlatField(tit);
}
} //else {
// cout << iframe << " " << p << endl;
if (iframe%p==0) {
while (mt->isBusy()) {;}//wait until all data are processed from the queues
sprintf(tit,"/home/l_msdetect/slsbl/x07mb/x07mbop/Data1/2017/11/06/tmp.tiff");
cout <<"*"<< iframe<< endl;
mt->writeImage(tit);
// mt->clearImage();
}
//}
ff++;
// if (iframe%p==0) {
// while (mt->isBusy()) {;}//wait until all data are processed from the queues
// sprintf(tit,"tmp.tiff",filename.c_str(), i_image);
// mt->writeImage(tit);
// // mt->clearImage();
// // }
// //interp->writeFlatField(tit);
// // i_image++;
// // sprintf(tit,"%s_tmp.tiff",fn);
// // mt->writeImage(tit);
// // //mt->clearImage();
// // cout <<"*"<< iprog++ << endl;
// }
// //stream data from socket 2
// if (send) {
// zmqsocket2->SendHeaderData(0, false, SLS_DETECTOR_JSON_HEADER_VERSION,
// 0,0,0,acqIndex,frameIndex,(char*)"run", acqIndex, 0,0,0,0,0,0,0,0,0,0,0,1);
// cprintf(GREEN, "Sent Header\n");
// zmqsocket2->SendData((char*)image,length);
// cprintf(GREEN, "Sent Data\n");
// }
}// exiting infinite loop
delete zmqsocket;
// if (send)
// delete zmqsocket2;
// cout<<"Goodbye"<< endl;
return 0;
}

251
slsDetectorCalibration/moenchExecutables/moench03ZmqOnTheFlyInterpolating.cpp Normal file
View File

@ -0,0 +1,251 @@
#include "ZmqSocket.h"
#include "ansi.h"
#include <iostream>
#include "moench03T1ZmqData.h"
// #include "interpolatingDetector.h"
#include "etaInterpolationPosXY.h"
// #include "linearInterpolation.h"
// #include "noInterpolation.h"
#include "multiThreadedAnalogDetector.h"
//#include "singlePhotonDetector.h"
#include "interpolatingDetector.h"
using namespace std;
#define NC 400
#define NR 400
#define SLS_DETECTOR_JSON_HEADER_VERSION 0x2
int main(int argc, char *argv[]) {
/**
* trial.o [socket ip] [starting port number] [outfname]
*
*/
int p=1000;
int fifosize=1000;
int nthreads=6;
int nsubpix=25;
int etabins=nsubpix*10;
double etamin=-1, etamax=2;
int csize=3;
int nx=400, ny=400;
int save=1;
int nsigma=5;
int nped=1000;
int ndark=100;
int ok;
int iprog=0;
char outfname[10000];
moench03T1ZmqData *decoder=new moench03T1ZmqData();
// cout << "decoder "<< endl;
etaInterpolationPosXY *interp=new etaInterpolationPosXY(nx, ny, nsubpix, etabins, etamin, etamax);
// linearInterpolation *interp=new linearInterpolation(NC, NR, nsubpix);
//noInterpolation *interp=new noInterpolation(NC, NR, nsubpix);
// interpolatingDetector *filter=new interpolatingDetector(decoder,interp, nsigma, 1, 0, nped, 100);
singlePhotonDetector *filter=new singlePhotonDetector(decoder,csize, nsigma, 1, 0, nped, 100);
char tit[10000];
// filter->readPedestals("/scratch/ped_100.tiff");
// interp->readFlatField("/scratch/eta_100.tiff",etamin,etamax);
// cout << "filter "<< endl;
int size = 327680;////atoi(argv[3]);
char* buff;
int* image;
//int* image =new int[327680/sizeof(int)];
filter->newDataSet();
multiThreadedAnalogDetector *mt=new multiThreadedAnalogDetector(filter,nthreads,fifosize);
mt->setFrameMode(eFrame);
// mt->setFrameMode(eFrame); //need to find a way to switch between flat and frames!
// mt->prepareInterpolation(ok);
mt->StartThreads();
mt->popFree(buff);
// help
if (argc < 3 || (argc > 3)) {
cprintf(RED, "Help: %s [receive socket ip] [fname]\n", argv[0]);
return EXIT_FAILURE;
}
// receive parameters
bool send = false;
char* socketip = argv[1];
uint32_t portnum = 30001;//atoi(argv[2]);
char fn[10000];
strcpy(fn, argv[2]);
// send parameters if any
// char* socketip2 = 0;
// uint32_t portnum2 = 0;
// if (argc > 4) {
// send = true;
// socketip2 = argv[4];
// portnum2 = atoi(argv[5]);
// }
cout << "\nrx socket ip : " << socketip <<
// "\nrx port num : " << portnum <<
// "\nsize : " << size;
"\nfname : " << fn ;
// if (send) {
// cout << "\nsd socket ip : " << socketip2 <<
// "\nsd port num : " << portnum2;
// }
cout << endl;
// receive socket
ZmqSocket* zmqsocket = new ZmqSocket(socketip,portnum);
if (zmqsocket->IsError()) {
cprintf(RED, "Error: Could not create Zmq socket on port %d with ip %s\n", portnum, socketip);
delete zmqsocket;
return EXIT_FAILURE;
}
printf("Zmq Client at %s\n", zmqsocket->GetZmqServerAddress());
// send socket
// ZmqSocket* zmqsocket2 = 0;
// if (send) {
// zmqsocket2 = new ZmqSocket(portnum2, socketip2);
// if (zmqsocket2->IsError()) {
// bprintf(RED, "Error: Could not create Zmq socket server on port %d and ip %s\n", portnum2, socketip2);
// delete zmqsocket2;
// delete zmqsocket;
// return EXIT_FAILURE;
// }
// printf("Zmq Server started at %s\n", zmqsocket2->GetZmqServerAddress());
// }
cout << "size " << 327680/sizeof(int) << endl;
// header variables
uint64_t acqIndex = -1;
uint64_t frameIndex = -1;
uint32_t subframeIndex = -1;
string filename = "";
int nnx, nny, nns;
int imsize=filter->getImageSize(nnx,nny,nns);
//int imsize=nx*ny;
int i_image=0;
cout << "Image size is "<< nnx << " " << nny << " " << nns << " " << imsize << endl;
int iframe=0;
// infinite loop
while(1) {
// cprintf(GREEN, "Got ?\n");
// get header, (if dummy, fail is on parse error or end of acquisition)
if (!zmqsocket->ReceiveHeader(0, acqIndex, frameIndex, subframeIndex, filename)) {
cprintf(RED, "Acquisition finished\n");
cout << "Recieved " << iframe << " frames " << endl;
while (mt->isBusy()) {;}//wait until all data are processed from the queues
// stream dummy to socket2 to signal end of acquisition
// if (send) {
// zmqsocket2->SendData((char*)(mt->getImage()),imsize*sizeof(int));
// cprintf(BLUE, "Sent Interpolated image\n");
// // zmqsocket2->SendHeaderData(0, true, SLS_DETECTOR_JSON_HEADER_VERSION);
// // cprintf(RED, "Sent Dummy\n");
// }
// if (save) {
sprintf(tit,"%s_%d.tiff",fn, i_image);
mt->writeImage(tit);
mt->clearImage();
// }
//interp->writeFlatField(tit);
i_image++;
// dont get data
continue; //continue to not get out
}
// cprintf(GREEN, "Got Header \n");
// get data
// cprintf(GREEN, "Got Header\n");
//image=(int*)buff;
//cout << buff << endl;
int length = zmqsocket->ReceiveData(0, (int*)buff, size);
// int length = zmqsocket->ReceiveData(0, (int*)image, size);
// cprintf(GREEN, "Got Data\n");
//processing with image
//...
mt->pushData(buff);
mt->nextThread();
// cout << " " << (void*)buff;
mt->popFree(buff);
iframe++;
// if (iframe%p==0) {
// while (mt->isBusy()) {;}//wait until all data are processed from the queues
// sprintf(tit,"%s_tmp.tiff",fn);
// mt->writeImage(tit);
// //mt->clearImage();
// cout <<"*"<< iprog++ << endl;
// }
// //stream data from socket 2
// if (send) {
// zmqsocket2->SendHeaderData(0, false, SLS_DETECTOR_JSON_HEADER_VERSION,
// 0,0,0,acqIndex,frameIndex,(char*)"run", acqIndex, 0,0,0,0,0,0,0,0,0,0,0,1);
// cprintf(GREEN, "Sent Header\n");
// zmqsocket2->SendData((char*)image,length);
// cprintf(GREEN, "Sent Data\n");
// }
}// exiting infinite loop
delete zmqsocket;
// if (send)
// delete zmqsocket2;
// cout<<"Goodbye"<< endl;
return 0;
}

219
slsDetectorCalibration/moenchExecutables/moench03ZmqOnTheFlyInterpolation.cpp Normal file
View File

@ -0,0 +1,219 @@
#include "ZmqSocket.h"
#include "ansi.h"
#include <iostream>
#include "moench03T1ZmqData.h"
#include "interpolatingDetector.h"
#include "etaInterpolationPosXY.h"
#include "linearInterpolation.h"
#include "noInterpolation.h"
#include "multiThreadedDetector.h"
using namespace std;
#define SLS_DETECTOR_JSON_HEADER_VERSION 0x2
int main(int argc, char *argv[]) {
/**
* trial.o [socket ip] [starting port number] [image size] [send_socket ip] [send port number]
*
*/
int nthreads=3;
int nsubpix=5;
int etabins=550;
double etamin=-1, etamax=2;
int nx=400, ny=400;
int save=1;
int nsigma=5;
int nped=1000;
int ok;
//int* image = new int[(size/sizeof(int))]();
char* buff;
moench03T1ZmqData *decoder=new moench03T1ZmqData();
// cout << "decoder "<< endl;
etaInterpolationPosXY *interp=new etaInterpolationPosXY(nx, ny, nsubpix, etabins, etamin, etamax);
//linearInterpolation *interp=new linearInterpolation(NC, NR, nsubpix);
//noInterpolation *interp=new noInterpolation(NC, NR, nsubpix);
interpolatingDetector *filter=new interpolatingDetector(decoder,interp, nsigma, 1, 0, nped, 10);
filter->readPedestals("/scratch/ped_100.tiff");
interp->readFlatField("/scratch/eta_100.tiff",etamin,etamax);
cout << "filter "<< endl;
filter->newDataSet();
multiThreadedDetector *mt=new multiThreadedDetector(filter,nthreads,100);
mt->setFrameMode(eFrame); //need to find a way to switch between flat and frames!
mt->prepareInterpolation(ok);
mt->StartThreads();
mt->popFree(buff);
// help
if (argc < 4 || (argc > 4 && argc != 6)) {
cprintf(RED, "Help: ./trial [receive socket ip] [receive starting port number] [image size] [send_socket ip] [send starting port number]\n");
return EXIT_FAILURE;
}
// receive parameters
bool send = false;
char* socketip = argv[1];
uint32_t portnum = atoi(argv[2]);
int size = atoi(argv[3]);
// send parameters if any
char* socketip2 = 0;
uint32_t portnum2 = 0;
if (argc > 4) {
send = true;
socketip2 = argv[4];
portnum2 = atoi(argv[5]);
}
cout << "\nrx socket ip : " << socketip <<
"\nrx port num : " << portnum <<
"\nsize : " << size;
if (send) {
cout << "\nsd socket ip : " << socketip2 <<
"\nsd port num : " << portnum2;
}
cout << endl;
// receive socket
ZmqSocket* zmqsocket = new ZmqSocket(socketip,portnum);
if (zmqsocket->IsError()) {
cprintf(RED, "Error: Could not create Zmq socket on port %d with ip %s\n", portnum, socketip);
delete zmqsocket;
return EXIT_FAILURE;
}
printf("Zmq Client at %s\n", zmqsocket->GetZmqServerAddress());
// send socket
ZmqSocket* zmqsocket2 = 0;
if (send) {
zmqsocket2 = new ZmqSocket(portnum2, socketip2);
if (zmqsocket2->IsError()) {
bprintf(RED, "Error: Could not create Zmq socket server on port %d and ip %s\n", portnum2, socketip2);
delete zmqsocket2;
delete zmqsocket;
return EXIT_FAILURE;
}
printf("Zmq Server started at %s\n", zmqsocket2->GetZmqServerAddress());
}
// header variables
uint64_t acqIndex = -1;
uint64_t frameIndex = -1;
uint32_t subframeIndex = -1;
string filename = "";
int imsize=nx*ny*nsubpix*nsubpix;
int i_image=0;
// infinite loop
while(1) {
// get header, (if dummy, fail is on parse error or end of acquisition)
if (!zmqsocket->ReceiveHeader(0, acqIndex, frameIndex, subframeIndex, filename)) {
cprintf(RED, "Acquisition finished\n");
while (mt->isBusy()) {;}//wait until all data are processed from the queues
// stream dummy to socket2 to signal end of acquisition
if (send) {
zmqsocket2->SendData((char*)(mt->getInterpolatedImage()),imsize*sizeof(int));
cprintf(BLUE, "Sent Interpolated image\n");
// zmqsocket2->SendHeaderData(0, true, SLS_DETECTOR_JSON_HEADER_VERSION);
// cprintf(RED, "Sent Dummy\n");
}
if (save) {
char tit[10000];
sprintf(tit,"/scratch/int_image_%d.tiff",i_image);
mt->writeInterpolatedImage(tit);
mt->clearInterpolatedImage(tit);
}
i_image++;
// dont get data
break; //continue to not get out
}
cprintf(GREEN, "Got Header \n");
// get data
int length = zmqsocket->ReceiveData(0, (int*)buff, size);
cprintf(GREEN, "Got Data\n");
//processing with image
//...
mt->pushData(buff);
mt->nextThread();
// cout << " " << (void*)buff;
mt->popFree(buff);
// //stream data from socket 2
// if (send) {
// zmqsocket2->SendHeaderData(0, false, SLS_DETECTOR_JSON_HEADER_VERSION,
// 0,0,0,acqIndex,frameIndex,(char*)"run", acqIndex, 0,0,0,0,0,0,0,0,0,0,0,1);
// cprintf(GREEN, "Sent Header\n");
// zmqsocket2->SendData((char*)image,length);
// cprintf(GREEN, "Sent Data\n");
// }
}// exiting infinite loop
delete zmqsocket;
if (send)
delete zmqsocket2;
cout<<"Goodbye"<< endl;
return 0;
}

248
slsDetectorCalibration/moenchExecutables/moench03ZmqOnTheFlySinglePhoton.cpp Normal file
View File

@ -0,0 +1,248 @@
#include "ZmqSocket.h"
#include "ansi.h"
#include <iostream>
#include "moench03T1ZmqData.h"
// #include "interpolatingDetector.h"
// #include "etaInterpolationPosXY.h"
#include "linearInterpolation.h"
// #include "noInterpolation.h"
#include "multiThreadedAnalogDetector.h"
//#include "singlePhotonDetector.h"
#include "interpolatingDetector.h"
using namespace std;
#define NC 400
#define NR 400
#define SLS_DETECTOR_JSON_HEADER_VERSION 0x2
int main(int argc, char *argv[]) {
/**
* trial.o [socket ip] [starting port number] [outfname]
*
*/
int p=1000;
int fifosize=1000;
int nthreads=6;
int nsubpix=2;
// int etabins=550;
// double etamin=-1, etamax=2;
int csize=3;
int nx=400, ny=400;
int save=1;
int nsigma=5;
int nped=1000;
int ndark=100;
int ok;
//int* image = new int[(size/sizeof(int))]();
char* buff;
int iprog=0;
char outfname[10000];
moench03T1ZmqData *decoder=new moench03T1ZmqData();
// cout << "decoder "<< endl;
//etaInterpolationPosXY *interp=new etaInterpolationPosXY(nx, ny, nsubpix, etabins, etamin, etamax);
// linearInterpolation *interp=new linearInterpolation(NC, NR, nsubpix);
//noInterpolation *interp=new noInterpolation(NC, NR, nsubpix);
// interpolatingDetector *filter=new interpolatingDetector(decoder,interp, nsigma, 1, 0, nped, 10);
singlePhotonDetector *filter=new singlePhotonDetector(decoder,csize, nsigma, 1, 0, nped, 100);
char tit[10000];
// filter->readPedestals("/scratch/ped_100.tiff");
// interp->readFlatField("/scratch/eta_100.tiff",etamin,etamax);
// cout << "filter "<< endl;
filter->newDataSet();
multiThreadedAnalogDetector *mt=new multiThreadedAnalogDetector(filter,nthreads,fifosize);
mt->setFrameMode(eFrame); //need to find a way to switch between flat and frames!
// mt->prepareInterpolation(ok);
mt->StartThreads();
mt->popFree(buff);
// help
if (argc < 4 || (argc > 4 && argc != 6)) {
cprintf(RED, "Help: %s [receive socket ip] [receive starting port number] [fname]\n", argv[0]);
return EXIT_FAILURE;
}
// receive parameters
bool send = false;
char* socketip = argv[1];
uint32_t portnum = atoi(argv[2]);
int size = 327680;////atoi(argv[3]);
char fn[10000];
strcpy(fn, argv[3]);
// send parameters if any
// char* socketip2 = 0;
// uint32_t portnum2 = 0;
// if (argc > 4) {
// send = true;
// socketip2 = argv[4];
// portnum2 = atoi(argv[5]);
// }
cout << "\nrx socket ip : " << socketip <<
"\nrx port num : " << portnum <<
// "\nsize : " << size;
"\nfname : " << fn ;
// if (send) {
// cout << "\nsd socket ip : " << socketip2 <<
// "\nsd port num : " << portnum2;
// }
cout << endl;
// receive socket
ZmqSocket* zmqsocket = new ZmqSocket(socketip,portnum);
if (zmqsocket->IsError()) {
cprintf(RED, "Error: Could not create Zmq socket on port %d with ip %s\n", portnum, socketip);
delete zmqsocket;
return EXIT_FAILURE;
}
printf("Zmq Client at %s\n", zmqsocket->GetZmqServerAddress());
// send socket
// ZmqSocket* zmqsocket2 = 0;
// if (send) {
// zmqsocket2 = new ZmqSocket(portnum2, socketip2);
// if (zmqsocket2->IsError()) {
// bprintf(RED, "Error: Could not create Zmq socket server on port %d and ip %s\n", portnum2, socketip2);
// delete zmqsocket2;
// delete zmqsocket;
// return EXIT_FAILURE;
// }
// printf("Zmq Server started at %s\n", zmqsocket2->GetZmqServerAddress());
// }
// header variables
uint64_t acqIndex = -1;
uint64_t frameIndex = -1;
uint32_t subframeIndex = -1;
string filename = "";
int nnx, nny, nns;
int imsize=filter->getImageSize(nnx,nny,nns);
//int imsize=nx*ny;
int i_image=0;
int iframe=0;
// infinite loop
while(1) {
// get header, (if dummy, fail is on parse error or end of acquisition)
if (!zmqsocket->ReceiveHeader(0, acqIndex, frameIndex, subframeIndex, filename)) {
cprintf(RED, "Acquisition finished\n");
cout << "Recieved " << iframe << " frames " << endl;
while (mt->isBusy()) {;}//wait until all data are processed from the queues
cout << "threads done " << endl;
// stream dummy to socket2 to signal end of acquisition
// if (send) {
// zmqsocket2->SendData((char*)(mt->getImage()),imsize*sizeof(int));
// cprintf(BLUE, "Sent Interpolated image\n");
// // zmqsocket2->SendHeaderData(0, true, SLS_DETECTOR_JSON_HEADER_VERSION);
// // cprintf(RED, "Sent Dummy\n");
// }
// if (save) {
sprintf(tit,"%s_%d.tiff",fn, i_image);
cout << tit << endl;
mt->writeImage(tit);
cout << "wrote" << endl;
mt->clearImage();
// }
i_image++;
// dont get data
continue; //continue to not get out
}
// cprintf(GREEN, "Got Header \n");
// get data
int length = zmqsocket->ReceiveData(0, (int*)buff, size);
// cprintf(GREEN, "Got Data\n");
//processing with image
//...
mt->pushData(buff);
mt->nextThread();
// cout << " " << (void*)buff;
mt->popFree(buff);
iframe++;
if (iframe%p==0) {
while (mt->isBusy()) {;}//wait until all data are processed from the queues
sprintf(tit,"%s_tmp.tiff",fn);
mt->writeImage(tit);
//mt->clearImage();
cout <<"*"<< iprog++ << endl;
}
// //stream data from socket 2
// if (send) {
// zmqsocket2->SendHeaderData(0, false, SLS_DETECTOR_JSON_HEADER_VERSION,
// 0,0,0,acqIndex,frameIndex,(char*)"run", acqIndex, 0,0,0,0,0,0,0,0,0,0,0,1);
// cprintf(GREEN, "Sent Header\n");
// zmqsocket2->SendData((char*)image,length);
// cprintf(GREEN, "Sent Data\n");
// }
}// exiting infinite loop
delete zmqsocket;
// if (send)
// delete zmqsocket2;
cout<<"Goodbye"<< endl;
return 0;
}

209
slsDetectorCalibration/moenchExecutables/readClusters.C Normal file
View File

@ -0,0 +1,209 @@
#include "single_photon_hit.h"
//#include "etaVEL/etaInterpolationPosXY.h"
TH2F *readClusters(char *fname, int nx, int ny, TH2F *h2=NULL) {
FILE *f=fopen(fname,"r");
int iph=0;
int ns=4;
double px, py;
double left, right, top, bottom;
if (f) {
int x1,y1;
if (h2==NULL)
h2=new TH2F("h2",fname,nx, -0.5, nx-0.5, ny, -0.5, ny-0.5);
h2mult=new TH2F("h2mult",fname,nx, -0.5, nx-0.5, ny, -0.5, ny-0.5);
TH2F *hint=new TH2F("hint",fname,nx*ns, -0.5, nx-0.5, ny*ns, -0.5, ny-0.5);
TH2F *hff=new TH2F("hff","hff",ns, -0.5, 0.5, ns, -0.5, +0.5);
TH1F *hsp=new TH1F("hsp",fname,500,0,2000);
TCanvas *c=new TCanvas();
c->SetLogz(kTRUE);
h2->Draw("colz");
TCanvas *c1=new TCanvas();
hsp->Draw();
TCanvas *c2=new TCanvas();
hint->Draw("colz");
c2->SetLogz(kTRUE);
single_photon_hit cl(3,3);
double tot;
int w;
double phw=340, phs=62;
while (cl.read(f)) {
//cl.get_pixel(x1, y1);
//cout << cl.iframe << " " << cl.x << " " << cl.y << endl;
//if (cl.x>80 && cl.x<280 && cl.y>80 && cl.y<300) {
tot=0;
left=0;
right=0;
top=0;
bottom=0;
for (int ix=-1; ix<2; ix++)
for (int iy=-1; iy<2; iy++){
tot+=cl.get_data(ix,iy);
if (ix<0) left+=cl.get_data(ix,iy);
if (ix>0) right+=cl.get_data(ix,iy);
if (iy<0) bottom+=cl.get_data(ix,iy);
if (iy>0) top+=cl.get_data(ix,iy);
}
px=(-left+right)/tot;
py=(-bottom+top)/tot;
//max at 340
if (tot>200) {
w=(tot+3.5*phs)/phw;
} else
w=0;
if (w) {
hsp->Fill(tot);
if (w==1) {
h2->Fill(cl.x,cl.y,w);
hint->Fill(px+cl.x,py+cl.y,w);
if (cl.y<350)
hff->Fill(px,py,w);
}
h2mult->Fill(cl.x,cl.y,w);
//}
iph+=w;
if (iph%100000==0) {
// c->Modified();
// c->Update();
c1->Modified();;
c1->Update();
// c2->Modified();;
// c2->Update();
}
}
// if (iph>0.5E7) break;
}
fclose(f);
hff->Scale(hff->GetNbinsX()*hff->GetNbinsY()/hff->Integral());
TH2F *hint2=(TH2F*)hint->Clone("hint2");
double ff;
for (int ibx=0; ibx<hint->GetNbinsX(); ibx++) {
for (int iby=0; iby<hint->GetNbinsY(); iby++) {
ff=hff->GetBinContent(ibx%ns+1, iby%ns+1);
// cout << ibx << " " << iby << " " << ibx%ns << " " << iby%ns << " " << ff << endl;
if (ff>0)
hint2->SetBinContent(ibx+1, iby+1,hint->GetBinContent(ibx+1,iby+1)/ff);
}
}
return h2;
} else
cout << "could not open file " << fname << endl;
return NULL;
}
// TH2F *getEta(char *fname, int nx, int ny, TH2F *h2=NULL) {
// int ns=10;
// slsInterpolation *inte=new etaInterpolationPosXY(nx,ny,ns,
// FILE *f=fopen(fname,"r");
// int iph=0;
// int ns=4;
// double px, py;
// double left, right, top, bottom;
// if (f) {
// int x1,y1;
// if (h2==NULL)
// h2=new TH2F("h2",fname,nx, -0.5, nx-0.5, ny, -0.5, ny-0.5);
// h2mult=new TH2F("h2mult",fname,nx, -0.5, nx-0.5, ny, -0.5, ny-0.5);
// TH2F *hint=new TH2F("hint",fname,nx*ns, -0.5, nx-0.5, ny*ns, -0.5, ny-0.5);
// // TH2F *hff=new TH2F("hff","hff",ns, -0.5, 0.5, ns, -0.5, +0.5);
// TH1F *hsp=new TH1F("hsp",fname,500,0,2000);
// TCanvas *c=new TCanvas();
// c->SetLogz(kTRUE);
// h2->Draw("colz");
// TCanvas *c1=new TCanvas();
// hsp->Draw();
// TCanvas *c2=new TCanvas();
// hint->Draw("colz");
// c2->SetLogz(kTRUE);
// single_photon_hit cl(3,3);
// double tot;
// int w;
// double phw=340, phs=62;
// while (cl.read(f)) {
// //cl.get_pixel(x1, y1);
// //cout << cl.iframe << " " << cl.x << " " << cl.y << endl;
// //if (cl.x>80 && cl.x<280 && cl.y>80 && cl.y<300) {
// tot=0;
// left=0;
// right=0;
// top=0;
// bottom=0;
// for (int ix=-1; ix<2; ix++)
// for (int iy=-1; iy<2; iy++){
// tot+=cl.get_data(ix,iy);
// if (ix<0) left+=cl.get_data(ix,iy);
// if (ix>0) right+=cl.get_data(ix,iy);
// if (iy<0) bottom+=cl.get_data(ix,iy);
// if (iy>0) top+=cl.get_data(ix,iy);
// }
// px=(-left+right)/tot;
// py=(-bottom+top)/tot;
// //max at 340
// if (tot>200) {
// w=(tot+3.5*phs)/phw;
// } else
// w=0;
// if (w) {
// hsp->Fill(tot);
// if (w==1) {
// h2->Fill(cl.x,cl.y,w);
// hint->Fill(px+cl.x,py+cl.y,w);
// }
// h2mult->Fill(cl.x,cl.y,w);
// // if (cl.y<350)
// // hff->Fill(px,py,w);
// //}
// iph+=w;
// if (iph%100000==0) {
// // c->Modified();
// // c->Update();
// c1->Modified();;
// c1->Update();
// // c2->Modified();;
// // c2->Update();
// }
// }
// // if (iph>0.5E7) break;
// }
// fclose(f);
// // hff->Scale(hff->GetNbinsX()*hff->GetNbinsY()/hff->Integral());
// // TH2F *hint2=(TH2F*)hint->Clone("hint2");
// // double ff;
// // for (int ibx=0; ibx<hint->GetNbinsX(); ibx++) {
// // for (int iby=0; iby<hint->GetNbinsY(); iby++) {
// // ff=hff->GetBinContent(ibx%ns+1, iby%ns+1);
// // // cout << ibx << " " << iby << " " << ibx%ns << " " << iby%ns << " " << ff << endl;
// // if (ff>0)
// // hint2->SetBinContent(ibx+1, iby+1,hint->GetBinContent(ibx+1,iby+1)/ff);
// // }
// // }
// return h2;
// } else
// cout << "could not open file " << fname << endl;
// return NULL;
// }

57
slsDetectorCalibration/moenchExecutables/tiff_to_th2f.cpp Normal file
View File

@ -0,0 +1,57 @@
#include <TPaveText.h>
#include <TLegend.h>
#include <TF1.h>
#include <TGraphErrors.h>
#include <TH2F.h>
#include <TASImage.h>
#include <TImage.h>
#include <TFile.h>
#include <vector>
#include <string>
#include <sstream>
#include <iomanip>
#include <fstream>
#include "tiffIO.h"
#include<iostream>
using namespace std;
int main(int argc, char *argv[]) {
/**
* trial.o [socket ip] [starting port number] [outfname]
*
*/
if (argc<3) {
cout << "Wrong usage! Should be: "<< argv[0] << " infile " << " outfile " << endl;
return 1;
}
uint32 nx, ny;
float *data=ReadFromTiff(argv[1],nx, ny);
TH2F *h2=NULL;
if (data) {
TFile *fout=new TFile(argv[2],"RECREATE");
if (fout) {
h2=new TH2F("h2",argv[1],nx,0,nx,ny,0, ny);
for (int ix=0; ix<nx ; ix++) {
for (int iy=0; iy<ny ; iy++) {
h2->SetBinContent(ix+1, iy+1, data[ix+iy*nx]);
}
}
h2->Write();
fout->Close();
}
delete [] data;
}
return 0;
}

244
slsDetectorCalibration/moenchMakeTree.C
View File

@ -1,244 +0,0 @@
#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 "RunningStat.h"
#include "MovingStat.h"
#include "moench02ModuleData.h"
#include <TThread.h>
using namespace std;
//tree variables
int xC,yC,iFrameC;
double meC,sigC;
Double_t dataC[3][3];
TTree* tall;
typedef struct task_s{
char *fformat;
char *tname;
int runmin;
int runmax;
int sign;
} Task;
void setUpTree(char *tname){
tall=new TTree(tname,tname);
tall->Branch("iFrame",&iFrameC,"iframe/I");
tall->Branch("x",&xC,"x/I");
tall->Branch("y",&yC,"y/I");
tall->Branch("data",dataC,"data[3][3]/D");
tall->Branch("pedestal",&meC,"pedestal/D");
tall->Branch("rms",&sigC,"rms/D");
}
inline void storeEvent(int iF,int x,int y, Double_t data[][3], double me, double sig){
TThread::Lock();
xC = x; yC = y; iFrameC = iF;
memcpy(dataC,data,sizeof(Double_t)*3*3);
//cout << "X: " << x << " Y: " << y << endl;
/* for(int i = 0; i < 3; i++){
for(int j = 0; j < 3; j++){
cout << "i: " << i << " j: " << j << " dataC " << dataC[i][j] << " data " << data[i][j] << endl;
}
}*/
meC = me; sigC = sig;
tall->Fill();
TThread::UnLock();
}
void moenchMakeTree(char *fformat, char *tname, int runmin, int runmax, int sign=1) {
double nThSigma = 3.;
moench02ModuleData *decoder=new moench02ModuleData();
char *buff;
char fname[10000];
int nf=0;
int dum, nPhotons;
double me, sig, tot, maxNei, val, valNei;
MovingStat stat[160][160];
MovingStat nPhotonsStat;
ifstream filebin;
int nbg=20;
int ix, iy, ir, ic, mx,my;
Double_t data[3][3];
nPhotonsStat.Clear();
nPhotonsStat.SetN(1000);
for (ir=0; ir<160; ir++) {
for (ic=0; ic<160; ic++) {
stat[ir][ic].Clear();
stat[ir][ic].SetN(nbg);
}
}
for (int irun=runmin; irun<runmax; irun++) {
sprintf(fname,fformat,irun);
//cout << "process file " << fname; // cout.flush();
filebin.open((const char *)(fname), ios::in | ios::binary);
while ((buff=decoder->readNextFrame(filebin))) {
nPhotons=0;
//for (ix=0; ix<160; ix++){
for (ix=0; ix<40; ix++){
for (iy=0; iy<160; iy++) {
dum=0; //no hit
// tot=0;
me=stat[iy][ix].Mean();
sig=stat[iy][ix].StandardDeviation();
val=sign*(decoder->getChannelShort(buff,ix,iy)-me);
if (nf>nbg) {
me=stat[iy][ix].Mean();
sig=stat[iy][ix].StandardDeviation();
val=sign*decoder->getChannel(buff,ix,iy)-me;
// dum=0; //no hit
tot=0;
maxNei = 0;
if (val>nThSigma*sig)//{ //hit check if neighbors are higher
dum=1;
for (ir=-1; ir<2; ir++){
for (ic=-1; ic<2; ic++){
if ((ix+ic)>=0 && (ix+ic)<160 && (iy+ir)>=0 && (iy+ir)<160) {
valNei = sign*decoder->getChannel(buff,ix+ic,iy+ir)-stat[iy+ir][ix+ic].Mean();
if (sign*decoder->getChannel(buff,ix+ic,iy+ir)>(stat[iy+ir][ix+ic].Mean()+3.*stat[iy+ir][ix+ic].StandardDeviation())) dum=1; //is a hit or neighbour is a hit!
tot+=valNei;
data[ir+1][ic+1] = valNei;
if(valNei/stat[iy+ir][ix+ic].StandardDeviation() > maxNei){
maxNei = valNei/stat[iy+ir][ix+ic].StandardDeviation();
mx = ir;
my = ic;
}
}
}
}
// }
if (val<(-nThSigma*sig)) dum=2; //discard negative events!
if(dum == 1 && mx == 0 && my == 0){ // this is an event and we are in the center
storeEvent(nf,ix,iy,data,me,sig);
nPhotons++;
//cout << "BF X: " << ix << " Y: " << iy << " val: " << val << " tot: " << tot << " me: " << me << " sig: " << sig << endl;
}
//if (tot>9.*sig && dum == 1){ dum=3;}
//cout << dum;
}
//if (ix==20 && iy==40)
//cout << decoder->getChannelShort(buff,ix,iy)<< " " << val << " " << me << " " << sig << " " << dum << endl;
if ( dum==0) {
stat[iy][ix].Calc(decoder->getChannelShort(buff,ix,iy));
}
if (nf<nbg || dum==0)
stat[iy][ix].Calc(sign*decoder->getChannel(buff,ix,iy));
}
}
delete [] buff;
//cout << "="; cout.flush();
if(nf>nbg) nPhotonsStat.Calc((double)nPhotons);
nf++;
}
//cout << endl;
cout << "processed File " << fname << " done. Avg. Photons/Frame: " << nPhotonsStat.Mean() << " sig: " << nPhotonsStat.StandardDeviation() << " " << runmax-irun << " files need processing" << endl;
if (filebin.is_open())
filebin.close();
else
cout << "could not open file " << fname << endl;
}
delete decoder;
cout << "Read " << nf << " frames" << endl;
}
void *moenchMakeTreeTask(void *p){
Task *t = (Task *)p;
moenchMakeTree(t->fformat,t->tname,t->runmin,t->runmax,t->sign);
return 0;
}
//to compile: g++ -DMYROOT -g `root-config --cflags --glibs` -o moenchMakeTree moenchMakeTree.C
int main(int argc, char **argv){
if(argc != 6){ cout << "Usage: inFile outdir tname fileNrStart fileNrEnd" << endl; exit(-1); }
int nThreads = 10;
TThread *threads[nThreads];
char *inFile = argv[1];
char *outDir = argv[2];
char *tName = argv[3];
int start = atoi(argv[4]);
int end = atoi(argv[5]);
TFile *f;
char outfname[1000];
char threadName[1000];
sprintf(outfname,"%s/%s.root",outDir,tName);
f=new TFile(outfname,"RECREATE");
cout << "outputfile: " << outfname << endl;
setUpTree(tName);
for(int i = 0; i < nThreads; i++){
sprintf(threadName,"t%i",i);
Task *t = (Task *)malloc(sizeof(Task));
t->fformat = inFile;
t->tname = tName;
t->sign = 1.;
t->runmin = start + (end-start)/(nThreads)*i;
t->runmax = start + (end-start)/(nThreads)*(i+1);
if(i == nThreads - 1) t->runmax = end;
cout << "start thread " << i << " start: " << t->runmin << " end " << t->runmax << endl;
threads[i] = new TThread(threadName, moenchMakeTreeTask, t);
threads[i]->Run();
//moenchMakeTreeTask(t);
}
//TThread::Ps();
for(int i = 0; i < nThreads; i++){
threads[i]->Join();
}
tall->Write();
tall->Print();
f->Close();
}

238
slsDetectorCalibration/moenchReadData.C
View File

@ -1,238 +0,0 @@
#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 "moench02ModuleData.h"
#include "moenchCommonMode.h"
#include "singlePhotonDetector.h"
//#include "MovingStat.h"
using namespace std;
#define NC 160
#define NR 160
#define MY_DEBUG 1
#ifdef MY_DEBUG
#include <TCanvas.h>
#endif
/**
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 *moenchReadData(char *fformat, char *tit, int runmin, int runmax, int nbins=1500, int hmin=-500, int hmax=1000, int sign=1, double hc=0, int xmin=1, int xmax=NC-1, int ymin=1, int ymax=NR-1, int cmsub=0, int hitfinder=1) {
moench02ModuleData *decoder=new moench02ModuleData(hc);
moenchCommonMode *cmSub=NULL;
if (cmsub)
cmSub=new moenchCommonMode();
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;
}

52
slsDetectorCalibration/moenchReadDataMT.C
View File

@ -1,52 +0,0 @@
#include <TFile.h>
#include <TThread.h>
#include "moenchReadData.C"
typedef struct task_s{
char *fformat;
char *tname;
char *tdir;
int runmin;
int runmax;
int treeIndex;
} Task;
void *moenchMakeTreeTask(void *p){
TThread::Lock();
char fname[1000];
Task *t = (Task *)p;
sprintf(fname,"%s%s_%i.root",t->tdir,t->tname,t->treeIndex);
TFile *f = new TFile(fname,"RECREATE");
cout << "Call moenchReadData(" << t->fformat << "," << t->tname << "," << t->runmin<< "," << t->runmax <<")" << endl;
TThread::UnLock();
moenchReadData(t->fformat,t->tname,t->runmin,t->runmax);
f->Close();
return 0;
}
void moenchReadDataMT(char *fformat, char *tit, char *tdir, int runmin, int runoffset, int nThreads, int treeIndexStart=0){
char threadName[1000];
TThread *threads[nThreads];
for(int i = 0; i < nThreads; i++){
sprintf(threadName,"t%i",i);
Task *t = (Task *)malloc(sizeof(Task));
t->fformat = fformat;
t->tname = tit;
t->tdir = tdir;
t->runmin = runmin + i*runoffset;
t->runmax = runmin + (i+1)*runoffset - 1;
t->treeIndex = treeIndexStart + i;
cout << "start thread " << i << " start: " << t->runmin << " end " << t->runmax << endl;
threads[i] = new TThread(threadName, moenchMakeTreeTask, t);
threads[i]->Run();
}
for(int i = 0; i < nThreads; i++){
threads[i]->Join();
}
}

572
slsDetectorCalibration/multiThreadedAnalogDetector.h Normal file
View File

@ -0,0 +1,572 @@
#ifndef MULTITHREADED_ANALOG_DETECTOR_H
#define MULTITHREADED_ANALOG_DETECTOR_H
#include <vector>
#include <string>
#include <sstream>
#include <iomanip>
#include <fstream>
#include <stdio.h>
//#include <deque>
//#include <list>
//#include <queue>
#include <fstream>
#include <cstdlib>
#include <pthread.h>
//#include "analogDetector.h"
#include "singlePhotonDetector.h"
//#include "interpolatingDetector.h"
#include "circularFifo.h"
#include "slsInterpolation.h"
#include <unistd.h>
#include <stdio.h>
#include <stdlib.h>
//#include <mutex>
using namespace std;
class threadedAnalogDetector
{
public:
threadedAnalogDetector(analogDetector<uint16_t> *d, int fs=10000) {
char *mem, *mm;
det=d;
fifoFree=new CircularFifo<char>(fs);
fifoData=new CircularFifo<char>(fs);
mem=(char*)malloc(fs*det->getDataSize());
// cout << "data size is " << det->getDataSize()*fs << endl;
for (int i=0; i<fs; i++) {
mm=mem+i*det->getDataSize();
fifoFree->push(mm);
}
busy=0;
stop=1;
fMode=eFrame;
ff=NULL;
}
virtual int setFrameMode(int fm) {fMode=fm; if (fMode>=0) det->setFrameMode((frameMode)fMode);};
//fMode=fm; return fMode;}
/* void prepareInterpolation(int &ok) { */
/* cout << "-" << endl; */
/* det->prepareInterpolation(ok); */
/* }; */
virtual int *getImage() {
return det->getImage();
}
virtual int getImageSize(int &nnx, int &nny, int &ns) {return det->getImageSize(nnx, nny, ns);};
virtual int getDetectorSize(int &nnx, int &nny) {return det->getDetectorSize(nnx, nny);};
~threadedAnalogDetector() {StopThread(); free(mem); delete fifoFree; delete fifoData;}
/** Returns true if the thread was successfully started, false if there was an error starting the thread */
virtual bool StartThread()
{ stop=0;
return (pthread_create(&_thread, NULL, processData, this) == 0);
}
virtual void StopThread()
{ stop=1;
(void) pthread_join(_thread, NULL);
}
virtual bool pushData(char* &ptr) {
fifoData->push(ptr);
}
virtual bool popFree(char* &ptr) {
fifoFree->pop(ptr);
}
virtual int isBusy() {return busy;}
//protected:
/** Implement this method in your subclass with the code you want your thread to run. */
//virtual void InternalThreadEntry() = 0;
virtual void *writeImage(const char * imgname) {cout << "a" <<endl; return det->writeImage(imgname);};
virtual void clearImage(){det->clearImage();};
virtual void prepareInterpolation(int &ok){
slsInterpolation *interp=(slsInterpolation*)det->getInterpolation();
if (interp)
interp->prepareInterpolation(ok);
}
virtual int *getFlatField(){
slsInterpolation *interp=(slsInterpolation*)det->getInterpolation();
if (interp)
return interp->getFlatField();
else
return NULL;
}
virtual int *setFlatField(int *ff, int nb, double emin, double emax){
slsInterpolation *interp=(slsInterpolation*)det->getInterpolation();
if (interp)
return interp->setFlatField(ff, nb, emin, emax);
else
return NULL;
}
virtual int *getFlatField(int &nb, double emi, double ema){
slsInterpolation *interp=(slsInterpolation*)det->getInterpolation();
int *ff;
if (interp) {
ff=interp->getFlatField(nb,emi,ema);
// cout << "tdgff* ff has " << nb << " bins " << endl;
return ff;
} else
return NULL;
}
virtual char *getInterpolation() {
return det->getInterpolation();
}
virtual void setPedestal(double *ped, double *rms=NULL, int m=-1){det->setPedestal(ped,rms,m);};
virtual void setPedestalRMS(double *rms){ det->setPedestalRMS(rms);};
virtual double *getPedestal(double *ped=NULL){return det->getPedestal(ped);};
virtual double *getPedestalRMS(double *rms=NULL){ return det->getPedestalRMS(rms);};
/** sets file pointer where to write the clusters to
\param f file pointer
\returns current file pointer
*/
FILE *setFilePointer(FILE *f){return det->setFilePointer(f); };
/** gets file pointer where to write the clusters to
\returns current file pointer
*/
FILE *getFilePointer(){return det->getFilePointer();};
void setMutex(pthread_mutex_t *fmutex){det->setMutex(fmutex);};
private:
analogDetector<uint16_t> *det;
int fMode;
int *dataSize;
pthread_t _thread;
char *mem;
CircularFifo<char> *fifoFree;
CircularFifo<char> *fifoData;
int stop;
int busy;
char *data;
int *ff;
static void * processData(void * ptr) {
threadedAnalogDetector *This=((threadedAnalogDetector *)ptr);
return This->processData();
}
void * processData() {
busy=1;
while (!stop) {
if (fifoData->isEmpty()) {
busy=0;
usleep(100);
} else {
busy=1;
fifoData->pop(data); //blocking!
det->processData(data);
fifoFree->push(data);
}
}
return NULL;
}
};
class multiThreadedAnalogDetector
{
public:
multiThreadedAnalogDetector(analogDetector<uint16_t> *d, int n, int fs=1000) : nThreads(n), ithread(0) {
dd[0]=d;
if (nThreads==1)
dd[0]->setId(100);
else
dd[0]->setId(0);
for (int i=1; i<nThreads; i++) {
dd[i]=d->Clone();
dd[i]->setId(i);
}
for (int i=0; i<nThreads; i++) {
dets[i]=new threadedAnalogDetector(dd[i], fs);
}
int nnx, nny, ns;
int nn=dets[0]->getImageSize(nnx, nny,ns);
image=new int[nn];
ff=NULL;
ped=NULL;
}
~multiThreadedAnalogDetector() {
StopThreads();
for (int i=0; i<nThreads; i++)
delete dets[i];
for (int i=1; i<nThreads; i++)
delete dd[i];
delete [] image;
}
int setFrameMode(int fm) { int ret; for (int i=0; i<nThreads; i++) ret=dets[i]->setFrameMode(fm); return ret;};
int *getImage(int &nnx, int &nny, int &ns) {
int *img;
// int nnx, nny, ns;
int nn=dets[0]->getImageSize(nnx, nny, ns);
//for (i=0; i<nn; i++) image[i]=0;
for (int ii=0; ii<nThreads; ii++) {
// cout << ii << " " << dets[ii]->getImageSize(nnx, nny, ns) << " " << nnx << " " << nny << " " << ns << endl;
img=dets[ii]->getImage();
for (int i=0; i<nn; i++) {
if (ii==0)
image[i]=img[i];
else
image[i]+=img[i];
}
}
return image;
}
void clearImage() {
for (int ii=0; ii<nThreads; ii++) {
dets[ii]->clearImage();
}
}
virtual void *writeImage(const char * imgname) {
/* #ifdef SAVE_ALL */
/* for (int ii=0; ii<nThreads; ii++) { */
/* char tit[10000];cout << "m" <<endl; */
/* sprintf(tit,"/scratch/int_%d.tiff",ii); */
/* dets[ii]->writeImage(tit); */
/* } */
/* #endif */
int nnx, nny, ns;
getImage(nnx, nny, ns);
//int nnx, nny, ns;
int nn=dets[0]->getImageSize(nnx, nny, ns);
float *gm=new float[nn];
if (gm) {
for (int ix=0; ix<nn; ix++) {
// if (image[ix]>0) cout << ix << " " << image[ix]<< endl;
gm[ix]=image[ix];
}
//cout << "image " << nnx << " " << nny << endl;
WriteToTiff(gm,imgname ,nnx, nny);
delete [] gm;
} else cout << "Could not allocate float image " << endl;
return NULL;
}
virtual void StartThreads() {
for (int i=0; i<nThreads; i++)
dets[i]->StartThread();
}
virtual void StopThreads() {
for (int i=0; i<nThreads; i++)
dets[i]->StopThread();
}
int isBusy() {
int ret=0, ret1;
for (int i=0; i<nThreads; i++) {
ret1=dets[i]->isBusy();
ret|=ret1;
// if (ret1) cout << "thread " << i <<" still busy " << endl;
}
return ret;
}
bool pushData(char* &ptr) {
dets[ithread]->pushData(ptr);
}
bool popFree(char* &ptr) {
dets[ithread]->popFree(ptr);
}
int nextThread() {
ithread++;
if (ithread==nThreads) ithread=0;
return ithread;
}
virtual void prepareInterpolation(int &ok){
getFlatField(); //sum up all etas
setFlatField(); //set etas to all detectors
for (int i=0; i<nThreads; i++) {
dets[i]->prepareInterpolation(ok);
}
}
virtual int *getFlatField(){
int nb=0;
double emi, ema;
int *f0;
slsInterpolation* inte=(slsInterpolation*)dets[0]->getInterpolation();
if (inte) {
if (inte->getFlatField(nb,emi,ema)) {
if (ff) delete [] ff;
ff=new int[nb*nb];
for (int i=0; i<nThreads; i++) {
// cout << "mtgff* ff has " << nb << " bins " << endl;
inte=(slsInterpolation*)dets[i]->getInterpolation();
f0=inte->getFlatField();
if (f0) {
// cout << "ff " << i << endl;
for (int ib=0; ib<nb*nb; ib++) {
if (i==0)
ff[ib]=f0[ib];
else
ff[ib]+=f0[ib];
/* if (i==0 && f0[ib]>0) */
/* cout << i << " " << ib << " " << f0[ib] << " " << ff[ib] << endl; */
}
}
}
return ff;
}
}
return NULL;
}
virtual int *setFlatField(int *h=NULL, int nb=-1, double emin=1, double emax=0){
//int nb=0;
double emi, ema;
slsInterpolation* inte=(slsInterpolation*)dets[0]->getFlatField(nb,emi,ema);
if (inte) {
if (h==NULL) h=ff;
for (int i=0; i<nThreads; i++) {
dets[i]->setFlatField(h, nb, emin, emax);
}
}
return NULL;
};
void *writeFlatField(const char * imgname){
int nb=0;
double emi, ema;
slsInterpolation* inte=(slsInterpolation*)dets[0]->getFlatField(nb,emi,ema);
if (inte) {
if (getFlatField()) {
// cout << "mtwff* ff has " << nb << " bins " << endl;
float *gm=new float[nb*nb];
if (gm) {
for (int ix=0; ix<nb*nb; ix++) {
gm[ix]=ff[ix];
}
WriteToTiff(gm,imgname ,nb, nb);
delete [] gm;
} else cout << "Could not allocate float image " << endl;
}
} else
cout << "Not interpolating detector! " << endl;
return NULL;
};
void *readFlatField(const char * imgname, int nb=-1, double emin=1, double emax=0){
int* inte=(int*)dets[0]->getFlatField(nb,emin,emax);
if (inte) {
uint32 nnx;
uint32 nny;
float *gm=ReadFromTiff(imgname, nnx, nny);
if (ff) delete [] ff;
if (nnx>nb) nb=nnx;
if (nny>nb) nb=nny;
ff=new int[nb*nb];
for (int ix=0; ix<nb*nb; ix++) {
ff[ix]=gm[ix];
}
delete [] gm;
return setFlatField(ff,nb,emin,emax);
} else
cout << "Not interpolating detector! " << endl;
return NULL;
};
virtual double *getPedestal(){
int nx, ny;
dets[0]->getDetectorSize(nx,ny);
if (ped) delete [] ped;
ped=new double[nx*ny];
double *p0=new double[nx*ny];
for (int i=0; i<nThreads; i++) {
//inte=(slsInterpolation*)dets[i]->getInterpolation(nb,emi,ema);
p0=dets[i]->getPedestal(p0);
if (p0) {
for (int ib=0; ib<nx*ny; ib++) {
ped[ib]+=p0[ib];
}
}
delete [] p0;
}
return ped;
};
virtual double *setPedestal(double *h=NULL){
//int nb=0;
int nx, ny;
dets[0]->getDetectorSize(nx,ny);
if (h==NULL) h=ped;
for (int i=0; i<nThreads; i++) {
dets[i]->setPedestal(h);
}
return NULL;
};
void *writePedestal(const char * imgname){
int nx, ny;
dets[0]->getDetectorSize(nx,ny);
getPedestal();
float *gm=new float[nx*ny];
if (gm) {
for (int ix=0; ix<nx*ny; ix++) {
gm[ix]=ped[ix];
}
WriteToTiff(gm,imgname ,nx, ny);
delete [] gm;
} else cout << "Could not allocate float image " << endl;
return NULL;
};
void *readPedestal(const char * imgname, int nb=-1, double emin=1, double emax=0){
int nx, ny;
dets[0]->getDetectorSize(nx,ny);
uint32 nnx;
uint32 nny;
float *gm=ReadFromTiff(imgname, nnx, nny);
if (ped) delete [] ped;
if (nnx>nx) nx=nnx;
if (nny>ny) ny=nny;
ped=new double[nx*ny];
for (int ix=0; ix<nx*ny; ix++) {
ped[ix]=gm[ix];
}
delete [] gm;
return setPedestal();
};
/** sets file pointer where to write the clusters to
\param f file pointer
\returns current file pointer
*/
FILE *setFilePointer(FILE *f){
for (int i=0; i<nThreads; i++) {
dets[i]->setFilePointer(f);
//dets[i]->setMutex(&fmutex);
}
return dets[0]->getFilePointer();
};
/** gets file pointer where to write the clusters to
\returns current file pointer
*/
FILE *getFilePointer(){return dets[0]->getFilePointer();};
private:
bool stop;
const int nThreads;
threadedAnalogDetector *dets[20];
analogDetector<uint16_t> *dd[20];
int ithread;
int *image;
int *ff;
double *ped;
pthread_mutex_t fmutex;
};
#endif

459
slsDetectorCalibration/multiThreadedDetector.h
View File

@ -1,459 +0,0 @@
#include <vector>
#include <string>
#include <sstream>
#include <iomanip>
#include <fstream>
#include <stdio.h>
//#include <deque>
//#include <list>
//#include <queue>
#include <fstream>
#include <pthread.h>
#include "analogDetector.h"
#include "circularFifo.h"
include "etaVEL/slsInterpolation.h"
class threadedDetector
{
public:
threadedDetector(analogDetector *d, int fs=10000) {
char *mem, *mm;
det=d;
fifoFree=new CircularFifo<char>(fs);
fifoData=new CircularFifo<char>(fs);
mem=(char*)malloc(fs*det->getDataSize());
for (int i=0; i<fs; i++) {
mm=mem+i*det->getDataSize();
fifoFree->push(mm);
}
busy=0;
stop=1;
fMode=eFrame;
}
virtual int setFrameMode(int fm) {if (fMode>=0) fMode=fm; return fMode;}
/* void prepareInterpolation(int &ok) { */
/* cout << "-" << endl; */
/* det->prepareInterpolation(ok); */
/* }; */
virtual int *getImage() {
return det->getInterpolatedImage();
}
virtual int getImageSize(int &nnx, int &nny, int &ns) {return det->getImageSize(nnx, nny, ns);};
~threadedDetector() {StopThread(); free(mem); delete fifoFree; delete fifoData;}
/** Returns true if the thread was successfully started, false if there was an error starting the thread */
virtual bool StartThread()
{ stop=0;
return (pthread_create(&_thread, NULL, processData, this) == 0);
}
virtual void StopThread()
{ stop=1;
(void) pthread_join(_thread, NULL);
}
virtual bool pushData(char* &ptr) {
fifoData->push(ptr);
}
virtual bool popFree(char* &ptr) {
fifoFree->pop(ptr);
}
virtual int isBusy() {return busy;}
//protected:
/** Implement this method in your subclass with the code you want your thread to run. */
//virtual void InternalThreadEntry() = 0;
virtual void *writeImage(const char * imgname) {cout << "a" <<endl; return det->writeImage(imgname);};
virtual void clearImage(){det->clearImage();};
private:
interpolatingDetector *det;
int fMode;
int *dataSize;
pthread_t _thread;
char *mem;
CircularFifo<char> *fifoFree;
CircularFifo<char> *fifoData;
int stop;
int busy;
char *data;
static void * processData(void * ptr) {
threadedDetector *This=((threadedDetector *)ptr);
return This->processData();
}
void * processData() {
busy=1;
while (!stop) {
if (fifoData->isEmpty()) {
busy=0;
usleep(100);
} else {
busy=1;
fifoData->pop(data); //blocking!
det->processData(data,(frameMode)fMode);
fifoFree->push(data);
}
}
return NULL;
}
};
class multiThreadedDetector
{
public:
multiThreadedDetector(analogDetector *d, int n, int fs=1000) : nThreads(n), ithread(0) {
dd[0]=d;
if (nThreads==1)
dd[0]->setId(100);
else
dd[0]->setId(0);
for (int i=1; i<nThreads; i++) {
dd[i]=d->Clone();
dd[i]->setId(i);
}
for (int i=0; i<nThreads; i++) {
dets[i]=new threadedDetector(dd[i], fs);
}
int nnx, nny, ns;
int nn=dets[0]->getImageSize(nnx, nny,ns);
image=new int[nn];
}
~multiThreadedDetector() {
StopThreads();
for (int i=0; i<nThreads; i++)
delete dets[i];
for (int i=1; i<nThreads; i++)
delete dd[i];
delete [] image;
}
int setFrameMode(int fm) { int ret; for (int i=0; i<nThreads; i++) ret=dets[i]->setFrameMode(fm); return ret;};
int *getImage() {
int *img;
int nnx, nny, ns;
int nn=dets[0]->getImageSize(nnx, nny, ns);
//for (i=0; i<nn; i++) image[i]=0;
for (int ii=0; ii<nThreads; ii++) {
img=dets[ii]->getImage();
for (int i=0; i<nn; i++) {
if (ii==0)
image[i]=img[i];
else
image[i]+=img[i];
}
}
return image;
}
void clearImage() {
for (int ii=0; ii<nThreads; ii++) {
dets[ii]->clearImage();
}
}
void *writeImage(const char * imgname) {
#ifdef SAVE_ALL
for (int ii=0; ii<nThreads; ii++) {
char tit[10000];cout << "m" <<endl;
sprintf(tit,"/scratch/int_%d.tiff",ii);
dets[ii]->writeImage(tit);
}
#endif
getImage();
int nnx, nny, ns;
int nn=dets[0]->getImageSize(nnx, nny, ns);
float *gm=new float[ nn];
if (gm) {
for (int ix=0; ix<nn; ix++) {
gm[ix]=image[ix];
}
WriteToTiff(gm,imgname ,nnx, nny);
delete [] gm;
} else cout << "Could not allocate float image " << endl;
return NULL;
}
void StartThreads() {
for (int i=0; i<nThreads; i++)
dets[i]->StartThread();
}
void StopThreads() {
for (int i=0; i<nThreads; i++)
dets[i]->StopThread();
}
int isBusy() {
int ret=0, ret1;
for (int i=0; i<nThreads; i++) {
ret1=dets[ithread]->isBusy();
ret|=ret1;
if (ret1) cout << "thread " << i <<" still busy " << endl;
}
return ret;
}
bool pushData(char* &ptr) {
dets[ithread]->pushData(ptr);
}
bool popFree(char* &ptr) {
dets[ithread]->popFree(ptr);
}
int nextThread() {
ithread++;
if (ithread==nThreads) ithread=0;
return ithread;
}
virtual void prepareInterpolation(int &ok){
slsInterpolation *
};
private:
bool stop;
const int nThreads;
threadedDetector *dets[20];
interpolatingDetector *dd[20];
int ithread;
int *image;
};
/* void *moenchProcessFrame() { */
/* char fname[10000]; */
/* strcpy(fname,"/mnt/moench_data/m03-15_mufocustube/plant_40kV_10uA/m03-15_100V_g4hg_300us_dtf_0.raw"); */
/* int nph, nph1; */
/* single_photon_hit clusters[NR*NC]; */
/* // cout << "hits "<< endl; */
/* int etabins=550; */
/* double etamin=-1, etamax=2; */
/* int nsubpix=2; */
/* float *etah=new float[etabins*etabins]; */
/* // cout << "etah "<< endl; */
/* cout << "image size "<< nsubpix*nsubpix*NC*NR << endl; */
/* float *image=new float[nsubpix*nsubpix*NC*NR]; */
/* int *heta, *himage; */
/* moench03Ctb10GbT1Data *decoder=new moench03Ctb10GbT1Data(); */
/* // cout << "decoder "<< endl; */
/* etaInterpolationPosXY *interp=new etaInterpolationPosXY(NC, NR, nsubpix, etabins, etamin, etamax); */
/* // cout << "interp "<< endl; */
/* // linearInterpolation *interp=new linearInterpolation(NC, NR, nsubpix); */
/* interpolatingDetector *filter=new interpolatingDetector(decoder,interp, 5, 1, 0, 1000, 100); */
/* cout << "filter "<< endl; */
/* char *buff; */
/* int nf=0; */
/* int ok=0; */
/* ifstream filebin; */
/* std::time_t end_time; */
/* int iFrame=-1; */
/* int np=-1; */
/* filter->newDataSet(); */
/* nph=0; */
/* nph1=0; */
/* int iph; */
/* cout << "file name " << fname << endl; */
/* filebin.open((const char *)(fname), ios::in | ios::binary); */
/* if (filebin.is_open()) */
/* cout << "Opened file " << fname<< endl; */
/* else */
/* cout << "Could not open file " << fname<< endl; */
/* while ((buff=decoder->readNextFrame(filebin, iFrame)) && nf<1.5E4) { */
/* if (nf<9E3) */
/* ; */
/* else if (nf<1.1E4) { */
/* filter->processData(buff,eFlat); */
/* } else { */
/* if (ok==0) { */
/* cout << "**************************************************************************"<< endl; */
/* heta=interp->getFlatField(); */
/* // for (int ii=0; ii<etabins*etabins; ii++) { */
/* // etah[ii]=(float)heta[ii]; */
/* // } */
/* std::time(&end_time); */
/* cout << std::ctime(&end_time) << " " << nf << endl; */
/* // WriteToTiff(etah, "/scratch/eta.tiff", etabins, etabins); */
/* interp->prepareInterpolation(ok); */
/* cout << "**************************************************************************"<< endl; */
/* std::time(&end_time); */
/* cout << std::ctime(&end_time) << " " << nf << endl; */
/* } */
/* filter->processData(buff,eFrame); */
/* } */
/* nph+=nph1; */
/* if (nf%1000==0) { */
/* std::time(&end_time); */
/* cout << std::ctime(&end_time) << " " << nf << endl; */
/* } */
/* nf++; */
/* delete [] buff; */
/* iFrame=-1; */
/* } */
/* if (filebin.is_open()) */
/* filebin.close(); */
/* else */
/* cout << "could not open file " << fname << endl; */
/* himage=interp->getInterpolatedImage(); */
/* for (int ii=0; ii<nsubpix*nsubpix*NC*NR; ii++) { */
/* image[ii]=(float)himage[ii]; */
/* } */
/* WriteToTiff(image, "/scratch/int_image.tiff", nsubpix*NR, nsubpix*NC); */
/* delete [] etah; */
/* delete [] image; */
/* delete interp; */
/* delete decoder; */
/* cout << "Read " << nf << " frames" << endl; */
/* return NULL; */
/* } */
/* int main(int argc, char *argv[]){ */
/* moenchProcessFrame(); */
/* } */

136
slsDetectorCalibration/raedNoiseData.C
View File

@ -1,136 +0,0 @@
#include "moenchReadData.C"
void raedNoiseData(char *tit, int sign=1){
char fname[1000];
char f[1000];
TFile *fout;
THStack *hs2N;
sprintf(fname,"/data/moench_xbox_20140113/MoTarget_45kV_0_8mA_120V_%s_0.root",tit);
fout=new TFile(fname,"RECREATE");
sprintf(fname,"/data/moench_xbox_20140113/MoTarget_45kV_0_8mA_120V_%s_f00000%%04d000_0.raw",tit);
hs2N=moenchReadData(fname,0,3000,1500,-500,2500,1,0.,1,159,1,159, 0,1);
hs2N->SetName(tit);
hs2N->SetTitle(tit);
(TH2F*)(hs2N->GetHists()->At(0))->Write();
(TH2F*)(hs2N->GetHists()->At(1))->Write();
(TH2F*)(hs2N->GetHists()->At(2))->Write();
(TH2F*)(hs2N->GetHists()->At(3))->Write();
(TH2F*)(hs2N->GetHists()->At(4))->Write();
fout->Close();
}
void raedNoiseDataN(char *tit, int sign=1){
char fname[1000];
char f[1000];
TFile *fout;
THStack *hs2N;
sprintf(fname,"/data/moench_xbox_20140116/noise_%s.root",tit);
fout=new TFile(fname,"RECREATE");
sprintf(fname,"/data/moench_xbox_20140116/noise_%s_f00000%%04d000_0.raw",tit);
hs2N=moenchReadData(fname,tit,0,3000,1500,-500,2500,sign,0.,1,159,1,159, 0,0);
hs2N->SetName(tit);
hs2N->SetTitle(tit);
(TH2F*)(hs2N->GetHists()->At(0))->Write();
// (TH2F*)(hs2N->GetHists()->At(1))->Write();
// (TH2F*)(hs2N->GetHists()->At(2))->Write();
// (TH2F*)(hs2N->GetHists()->At(3))->Write();
// (TH2F*)(hs2N->GetHists()->At(4))->Write();
fout->Close();
}
void g4() {
raedNoiseData("cds_g4_low_gain");
raedNoiseData("cds_g4_sto1_only");
raedNoiseData("cds_g4_no sto");
}
void no_cds() {
raedNoiseData("cds_disable_low_gain",-1);
raedNoiseData("cds_disable_sto1_only",-1);
raedNoiseData("cds_disable_no sto",-1);
}
void all_gains() {
raedNoiseData("cds_g2");
raedNoiseData("cds_g2HC");
raedNoiseData("cds_g1_2");
raedNoiseData("cds_g2_3");
}
void all_low_gains() {
raedNoiseData("cds_g2_low_gain");
raedNoiseData("cds_g2HC_low_gain");
raedNoiseData("cds_g1_2_low_gain");
raedNoiseData("cds_g2_3_low_gain");
}
/*
clkdivider data
/data/moench_xbox_20140114/MoTarget_45kV_0_8mA_12us_120V_cds_g1_clkdiv17_f000000010000_0.raw
-rw-rw-r-- 1 l_msdetect l_msdetect 51440000 Jan 14 12:40 /data/moench_xbox_20140114/MoTarget_45kV_0_8mA_12us_120V_cds_g1_clkdiv25_f000000010000_0.raw
-rw-rw-r-- 1 l_msdetect l_msdetect 51440000 Jan 14 13:26 /data/moench_xbox_20140114/MoTarget_45kV_0_8mA_12us_120V_cds_g1_clkdiv35_f000000010000_0.raw
-rw-rw-r-- 1 l_msdetect l_msdetect 51440000 Jan 14 14:09 /data/moench_xbox_20140114/MoTarget_45kV_0_8mA_12us_120V_cds_g1_clkdiv50_f000000010000_0.raw
-rw-rw-r-- 1 l_msdetect l_msdetect 51440000 Jan 14 14:54 /data/moench_xbox_20140114/MoTarget_45kV_0_8mA_12us_120V_cds_g1_clkdiv70_f000000010000_0.raw
-rw-rw-r-- 1 l_msdetect l_msdetect 51440000 Jan 14 16:42 /data/moench_xbox_20140114/MoTarget_45kV_0_8mA_12us_120V_cds_g1_clkdiv110_f000000010000_0.raw
-rw-rw-r-- 1 l_msdetect l_msdetect 51440000 Jan 14 17:27 /data/moench_xbox_20140114/MoTarget_45kV_0_8mA_12us_120V_cds_g1_clkdiv170_f000000010000_0.raw
*/
/* oversampled data
-rw-rw-r-- 1 l_msdetect l_msdetect 51440000 Jan 14 18:12 /data/moench_xbox_20140114/MoTarget_45kV_0_8mA_12us_120V_os10_16rows_f000000010000_0.raw
-rw-rw-r-- 1 l_msdetect l_msdetect 51440000 Jan 14 18:47 /data/moench_xbox_20140114/MoTarget_45kV_0_8mA_12us_120V_os10_16rows_f000000010000_1.raw
-rw-rw-r-- 1 l_msdetect l_msdetect 51440000 Jan 14 19:22 /data/moench_xbox_20140114/MoTarget_45kV_0_8mA_12us_120V_os10_16rows_f000000010000_2.raw
-rw-rw-r-- 1 l_msdetect l_msdetect 51440000 Jan 14 20:02 /data/moench_xbox_20140114/MoTarget_45kV_0_8mA_12us_120V_os10_16rows_f000000010000_3.raw
-rw-rw-r-- 1 l_msdetect l_msdetect 51440000 Jan 14 20:41 /data/moench_xbox_20140114/MoTarget_45kV_0_8mA_12us_120V_os10_16rows_f000000010000_4.raw
-rw-rw-r-- 1 l_msdetect l_msdetect 51440000 Jan 14 21:16 /data/moench_xbox_20140114/MoTarget_45kV_0_8mA_12us_120V_os10_16rows_f000000010000_5.raw
-rw-rw-r-- 1 l_msdetect l_msdetect 51440000 Jan 14 21:56 /data/moench_xbox_20140114/MoTarget_45kV_0_8mA_12us_120V_os10_16rows_f000000010000_6.raw
-rw-rw-r-- 1 l_msdetect l_msdetect 51440000 Jan 14 22:35 /data/moench_xbox_20140114/MoTarget_45kV_0_8mA_12us_120V_os10_16rows_f000000010000_7.raw
-rw-rw-r-- 1 l_msdetect l_msdetect 51440000 Jan 14 23:11 /data/moench_xbox_20140114/MoTarget_45kV_0_8mA_12us_120V_os10_16rows_f000000010000_8.raw
-rw-rw-r-- 1 l_msdetect l_msdetect 51440000 Jan 14 23:50 /data/moench_xbox_20140114/MoTarget_45kV_0_8mA_12us_120V_os10_16rows_f000000010000_9.raw
*/

31
slsDetectorCalibration/readJungfrauData.C
View File

@ -1,31 +0,0 @@
{ // script to run Jungfrau data analysis
gROOT->ProcessLine(".L jungfrauReadData.C++");//");
TFile *fout; // output file
THStack *hs2N; // collection of objects
//fout=new TFile("/mnt/slitnas/datadir_jungfrau02/analysis_tests/tests_random.root","RECREATE"); //
fout=new TFile("/mnt/slitnas/datadir_jungfrau02/20140404_lowEXrays/test.root","RECREATE");
hs2N=jungfrauReadData("/mnt/slitnas/datadir_jungfrau02/20140404_lowEXrays/Ti_1000clk_13kV40mA_%d.bin","tit",0,86,3000,-499.5,2500.5,1,0,1,47,1,47,1,1); //1,1);
//hs2N=jungfrauReadData("/mnt/slitnas/datadir_jungfrau02/20140228_GainBits/GSdata_laser_%d.bin","tit",0,5,1500,-500,2500,1,0.,1,47,1,47,0,1);//,"/mnt/slitnas/datadir_jungfrau02/analysis_tests/CalibrationParametersTest_fake.txt"); //1,1); // data set test GB -> set (3,0,0) in junfrauReadData.C
//hs2N=jungfrauReadData("/mnt/slitnas/datadir_jungfrau02/digital_bit_adc_test/vb_1825_%d.bin","tit",0,7,1500,-500,2500,1,0.,1,47,1,47,0,1);//,"/mnt/slitnas/datadir_jungfrau02/analysis_tests/CalibrationParametersTest_fake.txt"); //1,1); // data set test GB -> set (3,0,0) in junfrauReadData.C
//hs2N->SetName("cds_g4");
//hs2N->SetTitle("cds_g4");
(TH2F*)(hs2N->GetHists()->At(0))->Write(); // write hists
(TH2F*)(hs2N->GetHists()->At(1))->Write();
(TH2F*)(hs2N->GetHists()->At(2))->Write();
(TH2F*)(hs2N->GetHists()->At(3))->Write();
(TH2F*)(hs2N->GetHists()->At(4))->Write();
//(TH2F*)(hs2N->GetHists()->At(5))->Write();
//(TH2F*)(hs2N->GetHists()->At(6))->Write();
fout->Close(); // uncomment
}

366
slsDetectorCalibration/readMoench03Data.C
View File

@ -1,366 +0,0 @@
#include "moench03ReadData.C"
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,"%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);
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<nH; i++) {
if (hs->GetHists()->At(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){
// char fname[1000];
// char f[1000];
// TFile *fout;
// THStack *hs2N;
// sprintf(fname,"/data/moench_xbox_20140116/noise_%s.root",tit);
// fout=new TFile(fname,"RECREATE");
// sprintf(fname,"/data/moench_xbox_20140116/noise_%s_f00000%%04d000_0.raw",tit);
// hs2N=moenchReadData(fname,tit,0,3000,1500,-500,2500,sign,0.,1,159,1,159, 0,0);
// hs2N->SetName(tit);
// hs2N->SetTitle(tit);
// (TH2F*)(hs2N->GetHists()->At(0))->Write();
// // (TH2F*)(hs2N->GetHists()->At(1))->Write();
// // (TH2F*)(hs2N->GetHists()->At(2))->Write();
// // (TH2F*)(hs2N->GetHists()->At(3))->Write();
// // (TH2F*)(hs2N->GetHists()->At(4))->Write();
// fout->Close();
// }
// void g4() {
// raedNoiseData("cds_g4_low_gain");
// raedNoiseData("cds_g4_sto1_only");
// raedNoiseData("cds_g4_no sto");
// }
// void no_cds() {
// raedNoiseData("cds_disable_low_gain",-1);
// raedNoiseData("cds_disable_sto1_only",-1);
// raedNoiseData("cds_disable_no sto",-1);
// }
// void all_gains() {
// raedNoiseData("cds_g2");
// raedNoiseData("cds_g2HC");
// raedNoiseData("cds_g1_2");
// raedNoiseData("cds_g2_3");
// }
// void all_low_gains() {
// raedNoiseData("cds_g2_low_gain");
// raedNoiseData("cds_g2HC_low_gain");
// raedNoiseData("cds_g1_2_low_gain");
// raedNoiseData("cds_g2_3_low_gain");
// }
// /*
// clkdivider data
// /data/moench_xbox_20140114/MoTarget_45kV_0_8mA_12us_120V_cds_g1_clkdiv17_f000000010000_0.raw
// -rw-rw-r-- 1 l_msdetect l_msdetect 51440000 Jan 14 12:40 /data/moench_xbox_20140114/MoTarget_45kV_0_8mA_12us_120V_cds_g1_clkdiv25_f000000010000_0.raw
// -rw-rw-r-- 1 l_msdetect l_msdetect 51440000 Jan 14 13:26 /data/moench_xbox_20140114/MoTarget_45kV_0_8mA_12us_120V_cds_g1_clkdiv35_f000000010000_0.raw
// -rw-rw-r-- 1 l_msdetect l_msdetect 51440000 Jan 14 14:09 /data/moench_xbox_20140114/MoTarget_45kV_0_8mA_12us_120V_cds_g1_clkdiv50_f000000010000_0.raw
// -rw-rw-r-- 1 l_msdetect l_msdetect 51440000 Jan 14 14:54 /data/moench_xbox_20140114/MoTarget_45kV_0_8mA_12us_120V_cds_g1_clkdiv70_f000000010000_0.raw
// -rw-rw-r-- 1 l_msdetect l_msdetect 51440000 Jan 14 16:42 /data/moench_xbox_20140114/MoTarget_45kV_0_8mA_12us_120V_cds_g1_clkdiv110_f000000010000_0.raw
// -rw-rw-r-- 1 l_msdetect l_msdetect 51440000 Jan 14 17:27 /data/moench_xbox_20140114/MoTarget_45kV_0_8mA_12us_120V_cds_g1_clkdiv170_f000000010000_0.raw
// */
// /* oversampled data
// -rw-rw-r-- 1 l_msdetect l_msdetect 51440000 Jan 14 18:12 /data/moench_xbox_20140114/MoTarget_45kV_0_8mA_12us_120V_os10_16rows_f000000010000_0.raw
// -rw-rw-r-- 1 l_msdetect l_msdetect 51440000 Jan 14 18:47 /data/moench_xbox_20140114/MoTarget_45kV_0_8mA_12us_120V_os10_16rows_f000000010000_1.raw
// -rw-rw-r-- 1 l_msdetect l_msdetect 51440000 Jan 14 19:22 /data/moench_xbox_20140114/MoTarget_45kV_0_8mA_12us_120V_os10_16rows_f000000010000_2.raw
// -rw-rw-r-- 1 l_msdetect l_msdetect 51440000 Jan 14 20:02 /data/moench_xbox_20140114/MoTarget_45kV_0_8mA_12us_120V_os10_16rows_f000000010000_3.raw
// -rw-rw-r-- 1 l_msdetect l_msdetect 51440000 Jan 14 20:41 /data/moench_xbox_20140114/MoTarget_45kV_0_8mA_12us_120V_os10_16rows_f000000010000_4.raw
// -rw-rw-r-- 1 l_msdetect l_msdetect 51440000 Jan 14 21:16 /data/moench_xbox_20140114/MoTarget_45kV_0_8mA_12us_120V_os10_16rows_f000000010000_5.raw
// -rw-rw-r-- 1 l_msdetect l_msdetect 51440000 Jan 14 21:56 /data/moench_xbox_20140114/MoTarget_45kV_0_8mA_12us_120V_os10_16rows_f000000010000_6.raw
// -rw-rw-r-- 1 l_msdetect l_msdetect 51440000 Jan 14 22:35 /data/moench_xbox_20140114/MoTarget_45kV_0_8mA_12us_120V_os10_16rows_f000000010000_7.raw
// -rw-rw-r-- 1 l_msdetect l_msdetect 51440000 Jan 14 23:11 /data/moench_xbox_20140114/MoTarget_45kV_0_8mA_12us_120V_os10_16rows_f000000010000_8.raw
// -rw-rw-r-- 1 l_msdetect l_msdetect 51440000 Jan 14 23:50 /data/moench_xbox_20140114/MoTarget_45kV_0_8mA_12us_120V_os10_16rows_f000000010000_9.raw
// */

203
slsDetectorCalibration/singlePhotonDetector.h
View File

@ -16,13 +16,14 @@
#ifndef EVTYPE_DEF
#define EVTYPE_DEF
/** enum to define the even types */
enum eventType {
PEDESTAL=0,
NEIGHBOUR=1,
PHOTON=2,
PHOTON_MAX=3,
NEGATIVE_PEDESTAL=4,
UNDEFINED_EVENT=-1
PEDESTAL=0, /** pedestal */
NEIGHBOUR=1, /** neighbour i.e. below threshold, but in the cluster of a photon */
PHOTON=2, /** photon i.e. above threshold */
PHOTON_MAX=3, /** maximum of a cluster satisfying the photon conditions */
NEGATIVE_PEDESTAL=4, /** negative value, will not be accounted for as pedestal in order to avoid drift of the pedestal towards negative values */
UNDEFINED_EVENT=-1 /** undefined */
};
#endif
@ -70,26 +71,31 @@ public analogDetector<uint16_t> {
if (ny==1)
clusterSizeY=1;
cluster=new single_photon_hit(clusterSize,clusterSizeY);
setClusterSize(csize);
// cluster=new single_photon_hit(clusterSize,clusterSizeY);
clusters=new single_photon_hit[nx*ny];
cluster=clusters;
setClusterSize(csize);
nphTot=0;
nphFrame=0;
};
/**
destructor. Deletes the cluster structure and the pdestalSubtraction array
destructor. Deletes the cluster structure, the pdestalSubtraction and the image array
*/
virtual ~singlePhotonDetector() {delete cluster; for (int i=0; i<ny; i++) delete [] eventMask[i]; delete [] eventMask; };
/* /\** resets the eventMask to undefined and the commonModeSubtraction *\/ */
/* void newFrame(){analogDetector::newFrame(); for (int iy=0; iy<ny; iy++) for (int ix=0; ix<nx; ix++) eventMask[iy][ix]=UNDEFINED_EVENT; }; */
/**
copy constructor
\param orig detector to be copied
*/
singlePhotonDetector(singlePhotonDetector *orig) : analogDetector<uint16_t>(orig) {
nDark=orig->nDark;
myFile=orig->myFile;
eventMask=new eventType*[ny];
for (int i=0; i<ny; i++) {
@ -100,16 +106,29 @@ public analogDetector<uint16_t> {
nSigma=orig->nSigma;
clusterSize=orig->clusterSize;
clusterSizeY=orig->clusterSizeY;
cluster=new single_photon_hit(clusterSize,clusterSizeY);
// cluster=new single_photon_hit(clusterSize,clusterSizeY);
clusters=new single_photon_hit[nx*ny];
cluster=clusters;
setClusterSize(clusterSize);
quad=UNDEFINED_QUADRANT;
tot=0;
quadTot=0;
gmap=orig->gmap;
nphTot=0;
nphFrame=0;
}
/**
duplicates the detector structure
\returns new single photon detector with same parameters
*/
virtual singlePhotonDetector *Clone() {
return new singlePhotonDetector(this);
}
@ -130,17 +149,29 @@ public analogDetector<uint16_t> {
clusterSize=n;
if (cluster)
delete cluster;
if (ny>1)
if (ny>clusterSize)
clusterSizeY=clusterSize;
cluster=new single_photon_hit(clusterSize,clusterSizeY);
else
clusterSizeY=1;
for (int ip=0; ip<nx*ny; ip++)
(clusters+ip)->set_cluster_size(clusterSize,clusterSizeY);
//cluster=new single_photon_hit(clusterSize,clusterSizeY);
}
return clusterSize;
}
virtual int *getNPhotons(char *data, double thr=-1, int *nph=NULL) {
/**
converts the image into number of photons
\param data pointer to data
\param nph pointer where to add the calculated photons. If NULL, the internal image will be used
\returns array with data converted into number of photons.
*/
virtual int *getNPhotons(char *data, int *nph=NULL) {
nphFrame=0;
double val;
if (nph==NULL)
nph=image;
@ -174,50 +205,25 @@ public analogDetector<uint16_t> {
return nph;
} else {
if (thr>0) {
cout << "threshold" << endl;
for (int ix=0; ix<nx; ix++) {
for (int iy=0; iy<ny; iy++) {
// cout << id << " " << ix << " " << iy << endl;
if (gmap) {
g=gmap[iy*nx+ix];
if (g==0) g=-1.;
}
//return UNDEFINED_EVENT;
for (int ix=xmin; ix<xmax; ix++) {
for (int iy=ymin; iy<ymax; iy++) {
val=subtractPedestal(data,ix,iy);
/* if (thr<=0) { */
/* tthr=nSigma*getPedestalRMS(ix,iy)/g; */
/* } */
/* */
nn=analogDetector<uint16_t>::getNPhotons(data,ix,iy,tthr);
nn=analogDetector<uint16_t>::getNPhotons(data,ix,iy);
nph[ix+nx*iy]+=nn;
rest[iy][ix]=(val-nn*tthr);
}
}
}
for (int ix=0; ix<nx; ix++) {
for (int iy=0; iy<ny; iy++) {
// }
for (int ix=xmin; ix<xmax; ix++) {
for (int iy=ymin; iy<ymax; iy++) {
// 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)/g;
if (ix==0 || iy==0)
rest[iy][ix]=subtractPedestal(data,ix,iy);
}
max=0;
tl=0;
tr=0;
@ -231,8 +237,6 @@ public analogDetector<uint16_t> {
if ((iy+ir)>=0 && (iy+ir)<ny && (ix+ic)>=0 && (ix+ic)<nx) {
//cluster->set_data(rest[iy+ir][ix+ic], ic, ir);
if (thr<=0 && ir>=0 && ic>=0 )
rest[iy+ir][ix+ic]=subtractPedestal(data,ix+ic,iy+ir);
v=rest[iy+ir][ix+ic];//cluster->get_data(ic,ir);
tot+=v;
@ -257,10 +261,7 @@ public analogDetector<uint16_t> {
}
}
}
if (rest[iy][ix]<=-tthr) {
eventMask[iy][ix]=NEGATIVE_PEDESTAL;
//if (cluster->get_data(0,0)>=max) {
} else if (rest[iy][ix]>=max) {
if (rest[iy][ix]>=max) {
if (bl>=br && bl>=tl && bl>=tr) {
quad=BOTTOM_LEFT;
quadTot=bl;
@ -275,35 +276,21 @@ public analogDetector<uint16_t> {
quadTot=tr;
}
if (max>tthr || tot>sqrt(ccy*ccs)*tthr || quadTot>sqrt(cy*cs)*tthr) {
/* cout << ix << " " << iy << " " << rest[iy][ix] <<" " << tot << " " << quadTot << endl; */
/* for (int ir=-(clusterSizeY/2); ir<(clusterSizeY/2)+1; ir++) { */
/* for (int ic=-(clusterSize/2); ic<(clusterSize/2)+1; ic++) */
/* cout << rest[iy+ir][ix+ic] << " " ; */
/* cout << endl; */
/* } */
eventMask[iy][ix]=PHOTON;
nph[ix+nx*iy]++;
// rest[iy][ix]-=tthr;
} //else if (thr<=0 ) {
//addToPedestal(data,ix,iy);
// }
}
if (thr<0 && eventMask[iy][ix]==PEDESTAL) {
addToPedestal(data,ix,iy);
}
}
}
}
// }
// }
nphFrame+=nph[ix+nx*iy];
nphTot+=nph[ix+nx*iy];
}
}
}
}
return nph;
} else return getClusters(data);
}
return NULL;
};
/** finds event type for pixel and fills cluster structure. The algorithm loops only if the evenMask for this pixel is still undefined.
@ -415,7 +402,16 @@ public analogDetector<uint16_t> {
};
int getClusters(char *data, single_photon_hit *clusters) {
/**
Loops in the region of interest to find the clusters
\param data pointer to the data structure
\returns number of clusters found
*/
int *getClusters(char *data) {
int nph=0;
@ -431,8 +427,8 @@ int getClusters(char *data, single_photon_hit *clusters) {
return 0;
}
newFrame();
for (int ix=0; ix<nx; ix++) {
for (int iy=0; iy<ny; iy++) {
for (int ix=xmin; ix<xmax; ix++) {
for (int iy=ymin; iy<ymax; iy++) {
max=0;
tl=0;
@ -502,6 +498,7 @@ int getClusters(char *data, single_photon_hit *clusters) {
(clusters+nph)->tot=tot;
(clusters+nph)->x=ix;
(clusters+nph)->y=iy;
(clusters+nph)->iframe=det->getFrameNumber(data);
(clusters+nph)->ped=getPedestal(ix,iy,0);
for (int ir=-(clusterSizeY/2); ir<(clusterSizeY/2)+1; ir++) {
for (int ic=-(clusterSize/2); ic<(clusterSize/2)+1; ic++) {
@ -509,7 +506,7 @@ int getClusters(char *data, single_photon_hit *clusters) {
}
}
nph++;
image[iy*nx+ix]++;
} else {
eventMask[iy][ix]=PHOTON;
@ -521,8 +518,12 @@ int getClusters(char *data, single_photon_hit *clusters) {
}
}
return nph;
nphFrame=nph;
nphTot+=nph;
//cout << nphFrame << endl;
// cout <<"**********************************"<< endl;
writeClusters();
return image;
};
@ -593,37 +594,59 @@ int getClusters(char *data, single_photon_hit *clusters) {
return tall;
};
#else
/** write cluster to filer*/
void writeCluster(FILE* myFile){cluster->write(myFile);};
/** write cluster to filer
\param f file pointer
*/
void writeCluster(FILE* f){cluster->write(f);};
/**
write clusters to file
\param f file pointer
\param clusters array of clusters structures
\param nph number of clusters to be written to file
*/
static void writeClusters(FILE *f, single_photon_hit *clusters, int nph){for (int i=0; i<nph; i++) (clusters+i)->write(f);};
void writeClusters(FILE *f){for (int i=0; i<nphFrame; i++) (clusters+i)->write(f);};
void writeClusters(){if (myFile) {
//cout << "++" << endl;
pthread_mutex_lock(fm);
for (int i=0; i<nphFrame; i++)
(clusters+i)->write(myFile);
pthread_mutex_unlock(fm);
//cout << "--" << endl;
}
};
#endif
virtual void processData(char *data, frameMode i=eFrame, int *val=NULL) {
virtual void processData(char *data, int *val=NULL) {
// cout << "sp" << endl;
switch(i) {
switch(fMode) {
case ePedestal:
addToPedestal(data);
break;
default:
getNPhotons(data,-1,val);
getNPhotons(data,val);
}
iframe++;
// cout << "done" << endl;
};
protected:
int nDark; /**< number of frames to be used at the beginning of the dataset to calculate pedestal without applying photon discrimination */
eventType **eventMask; /**< matrix of event type or each pixel */
double nSigma; /**< number of sigma parameter for photon discrimination */
int clusterSize; /**< cluster size in the x direction */
int clusterSizeY; /**< cluster size in the y direction i.e. 1 for strips, clusterSize for pixels */
single_photon_hit *cluster; /**< single photon hit data structure */
single_photon_hit *clusters; /**< single photon hit data structure */
quadrant quad; /**< quadrant where the photon is located */
double tot; /**< sum of the 3x3 cluster */
double quadTot; /**< sum of the maximum 2x2cluster */
int nphTot;
int nphFrame;
};

27
slsDetectorCalibration/single_photon_hit.h
View File

@ -33,13 +33,24 @@ class single_photon_hit {
/** binary write to file of all elements of the structure, except size of the cluster
\param myFile file descriptor
*/
void write(FILE *myFile) {fwrite((void*)this, 1, 3*sizeof(int)+4*sizeof(double)+sizeof(quad), myFile); fwrite((void*)data, 1, dx*dy*sizeof(double), myFile);};
size_t write(FILE *myFile) {
//fwrite((void*)this, 1, 3*sizeof(int)+4*sizeof(double)+sizeof(quad), myFile);
if (fwrite((void*)this, 1, 3*sizeof(int), myFile))
return fwrite((void*)data, 1, dx*dy*sizeof(double), myFile);
return 0;
};
/**
binary read from file of all elements of the structure, except size of the cluster. The structure is then filled with those args
\param myFile file descriptor
*/
void read(FILE *myFile) {fread((void*)this, 1, 3*sizeof(int)+4*sizeof(double)+sizeof(quad), myFile); fread((void*)data, 1, dx*dy*sizeof(double), myFile);};
size_t read(FILE *myFile) {
//fread((void*)this, 1, 3*sizeof(int)+4*sizeof(double)+sizeof(quad), myFile);
if (fread((void*)this, 1, 3*sizeof(int), myFile))
return fread((void*)data, 1, dx*dy*sizeof(double), myFile);
return 0;
};
/**
assign the value to the element of the cluster matrix, with relative coordinates where the center of the cluster is (0,0)
@ -49,6 +60,9 @@ class single_photon_hit {
*/
void set_data(double v, int ix, int iy=0){data[(iy+dy/2)*dx+ix+dx/2]=v;};
void set_cluster_size(int nx, int ny) {if (nx>0) dx=nx; if (ny>0) dy=ny; delete [] data; data=new double[dx*dy];};
void get_cluster_size(int &nx, int &ny) {nx=dx; ny=dy;};
void get_pixel(int &x1, int &y1) {x1=x; y1=y;};
/**
gets the value to the element of the cluster matrix, with relative coordinates where the center of the cluster is (0,0)
@ -58,16 +72,17 @@ class single_photon_hit {
*/
double get_data(int ix, int iy=0){return data[(iy+dy/2)*dx+ix+dx/2];};
double *get_cluster() {return data;};
int iframe; /**< frame number */
int x; /**< x-coordinate of the center of hit */
int y; /**< x-coordinate of the center of hit */
double rms; /**< noise of central pixel l -- at some point it can be removed*/
double ped; /**< pedestal of the central pixel -- at some point it can be removed*/
int iframe; /**< frame number */
const int dx; /**< size of data cluster in x */
const int dy; /**< size of data cluster in y */
quadrant quad; /**< quadrant where the photon is located */
double tot; /**< sum of the 3x3 cluster */
quadrant quad; /**< quadrant where the photon is located */
double quadTot; /**< sum of the maximum 2x2cluster */
int dx; /**< size of data cluster in x */
int dy; /**< size of data cluster in y */
double *data; /**< pointer to data */
};

18
slsDetectorCalibration/slsDetectorCalibration.doxy
View File

@ -80,8 +80,24 @@ LATEX_HIDE_INDICES = YES
PREDEFINED = __cplusplus
INPUT = analogDetector.h pedestalSubtraction.h MovingStat.h singlePhotonDetector.h interpolatingDetector.h tiffIO.h single_photon_hit.h slsDetectorData.h moench03Ctb10GbT1Data.h moench03TCtbData.h moench03T1CtbData.h moench03CtbData.h moench03Ctb10GbData.h moench03TCtb10GbData.h Mythen3_01_jctbData.h adcSar2_jctbData.h eigerHalfModuleData.h energyCalibration.h slsReceiverData.h gotthardModuleData.h gotthardShortModuleData.h jungfrau02Data.h jungfrau10ModuleData.h moench02Ctb10GbData.h moench02CtbData.h moench02ModuleData.h moench03CommonMode.h moenchCommonMode.h chiptestBoardData.h commonModeSubtraction.h RunningStat.h etaVEL/etaInterpolationBase.h etaVEL/slsInterpolation.h etaVEL/etaInterpolationPosXY.h etaVEL/linearInterpolation.h etaVEL/noInterpolation.h etaVEL/etaInterpolationGlobal.h etaVEL/interpolation_EtaVEL.h etaVEL/EtaVEL.h etaVEL/iterativeEtaInterpolation.h multiThreadedDetector.h moench03T1ZmqData.h
INPUT = analogDetector.h \
pedestalSubtraction.h \
MovingStat.h \
singlePhotonDetector.h \
interpolatingDetector.h tiffIO.h \
single_photon_hit.h \
dataStructures/slsDetectorData.h \
commonModeSubtraction.h \
RunningStat.h \
interpolations/etaInterpolationBase.h \
interpolations/slsInterpolation.h \
interpolations/etaInterpolationPosXY.h \
interpolations/linearInterpolation.h \
interpolations/noInterpolation.h \
interpolations/etaInterpolationGlobal.h \
multiThreadedAnalogDetector.h
OUTPUT_DIRECTORY = slsDetectorCalibrationDocs
#moench03Ctb10GbT1Data.h moench03TCtbData.h moench03T1CtbData.h moench03CtbData.h moench03Ctb10GbData.h moench03TCtb10GbData.h Mythen3_01_jctbData.h adcSar2_jctbData.h eigerHalfModuleData.h energyCalibration.h slsReceiverData.h gotthardModuleData.h gotthardShortModuleData.h jungfrau02Data.h jungfrau10ModuleData.h moench02Ctb10GbData.h moench02CtbData.h moench02ModuleData.h moench03CommonMode.h moenchCommonMode.h chiptestBoardData.h interpolation/etaVEL/interpolation_EtaVEL.h interpolation/etaVEL/EtaVEL.h etaVEL/iterativeEtaInterpolation.h dataStructures/moench03T1ZmqData.h

178
slsDetectorCalibration/viewMoenchDRXRays.C
View File

@ -1,178 +0,0 @@
#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;
}