Merge branch 'developer' into esrfchanges

This commit is contained in:
maliakal_d 2018-05-02 16:57:03 +02:00
commit 5b01b4cbd8
45 changed files with 2716 additions and 1467 deletions

View File

@ -10,6 +10,17 @@
#include "tiffIO.h"
#ifdef ROOTSPECTRUM
#include <TPaveText.h>
#include <TLegend.h>
#include <TF1.h>
#include <TGraphErrors.h>
#include <TH2F.h>
#include <TASImage.h>
#include <TImage.h>
#include <TFile.h>
#endif
using namespace std;
#ifndef FRAMEMODE_DEF
@ -66,11 +77,30 @@ template <class dataType> class analogDetector {
thr=0;
myFile=NULL;
fm=new pthread_mutex_t ;
#ifdef ROOTSPECTRUM
hs=new TH2F("hs","hs",2000,-100,10000,nx*ny,-0.5,nx*ny-0.5);
#ifdef ROOTCLUST
hs3=new TH2F("hs3","hs3",2000,-100,3*3*10000,nx*ny,-0.5,nx*ny-0.5);
hs5=new TH2F("hs5","hs5",2000,-100,5*5*10000,nx*ny,-0.5,nx*ny-0.5);
hs7=new TH2F("hs7","hs7",2000,-100,7*7*10000,nx*ny,-0.5,nx*ny-0.5);
hs9=new TH2F("hs9","hs9",2000,-100,9*9*10000,nx*ny,-0.5,nx*ny-0.5);
#endif
#endif
};
/**
destructor. Deletes the pdestalSubtraction array and the image
*/
virtual ~analogDetector() {for (int i=0; i<ny; i++) delete [] stat[i]; delete [] stat; delete [] image;};
virtual ~analogDetector() {for (int i=0; i<ny; i++) delete [] stat[i]; delete [] stat; delete [] image;
#ifdef ROOTSPECTRUM
delete hs;
#ifdef ROOTCLUST
delete hs3;
delete hs5;
delete hs7;
delete hs9;
#endif
#endif
};
/**
constructor cloning another analog detector
@ -111,7 +141,16 @@ template <class dataType> class analogDetector {
}
}
image=new int[nx*ny];
#ifdef ROOTSPECTRUM
hs=(TH2F*)(orig->hs)->Clone();//new TH2F("hs","hs",(orig->hs)-getNbins,-500,9500,nx*ny,-0.5,nx*ny-0.5);
#ifdef ROOTCLUST
hs3=(TH2F*)(orig->hs3)->Clone();//new TH2F("hs","hs",(orig->hs)-getNbins,-500,9500,nx*ny,-0.5,nx*ny-0.5);
hs5=(TH2F*)(orig->hs5)->Clone();//new TH2F("hs","hs",(orig->hs)-getNbins,-500,9500,nx*ny,-0.5,nx*ny-0.5);
hs7=(TH2F*)(orig->hs7)->Clone();//new TH2F("hs","hs",(orig->hs)-getNbins,-500,9500,nx*ny,-0.5,nx*ny-0.5);
hs9=(TH2F*)(orig->hs9)->Clone();//new TH2F("hs","hs",(orig->hs)-getNbins,-500,9500,nx*ny,-0.5,nx*ny-0.5);
#endif
#endif
}
@ -220,6 +259,15 @@ template <class dataType> class analogDetector {
image[iy*nx+ix]=0;
}
if (cmSub) cmSub->Clear();
#ifdef ROOTSPECTRUM
hs->Reset();
#ifdef ROOTCLUST
hs3->Reset();//new TH2F("hs","hs",(orig->hs)-getNbins,-500,9500,nx*ny,-0.5,nx*ny-0.5);
hs5->Reset();//new TH2F("hs","hs",(orig->hs)-getNbins,-500,9500,nx*ny,-0.5,nx*ny-0.5);
hs7->Reset();//new TH2F("hs","hs",(orig->hs)-getNbins,-500,9500,nx*ny,-0.5,nx*ny-0.5);
hs9->Reset();//new TH2F("hs","hs",(orig->hs)-getNbins,-500,9500,nx*ny,-0.5,nx*ny-0.5);
#endif
#endif
};
/** resets the commonModeSubtraction and increases the frame index */
@ -253,24 +301,60 @@ template <class dataType> class analogDetector {
\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) {
virtual void addToPedestal(double val, int ix, int iy=0, int cm=0){
if (ix>=0 && ix<nx && iy>=0 && iy<ny) {
if (cmSub && cm>0) {
val-= getCommonMode(ix, iy);
}
stat[iy][ix].addToPedestal(val);
if (cmSub) {
if (det) if (det->isGood(ix, iy)==0) return;
cmSub->addToCommonMode(val, ix, iy);
};
/* if (cmSub && cm>0) { */
/* if (det) if (det->isGood(ix, iy)==0) return; */
/* cmSub->addToCommonMode(val, ix, iy); */
/* }; */
};
}
/**
double getCommonMode(int ix, int iy) {
if (cmSub) return cmSub->getCommonMode(ix, iy);
else return 0;
}
virtual void addToCommonMode(char *data){
if (cmSub) {
for (int ix=xmin; ix<xmax; ix++) {
for (int iy=ymin; iy<ymax; iy++) {
// if (getNumpedestals(ix,iy)>0)
addToCommonMode(data, ix, iy);
}
}
//cout << "cm " << getCommonMode(0,0) << " " << getCommonMode(1,0) << endl;
}
}
virtual void addToCommonMode(char *data, int ix, int iy=0){
if (cmSub) {
if (det) if (det->isGood(ix, iy)==0) return;
if (getNumpedestals(ix,iy)>0){
cmSub->addToCommonMode(subtractPedestal(data,ix,iy,0), ix, iy);
// cout << ix << " " <<subtractPedestal(data,ix,iy,0) << endl;
}
}
}
/**
gets pedestal (and common mode)
\param ix pixel x coordinate
\param iy pixel y coordinate
\param cm 0 (default) without common mode subtraction, 1 with common mode subtraction (if defined)
\returns pedestal value
*/
virtual double getPedestal(int ix, int iy, int cm=0){if (ix>=0 && ix<nx && iy>=0 && iy<ny) if (cmSub && cm>0) return stat[iy][ix].getPedestal()-cmSub->getCommonMode(); else return stat[iy][ix].getPedestal(); else return -1;};
virtual double getPedestal(int ix, int iy, int cm=0){
if (ix>=0 && ix<nx && iy>=0 && iy<ny)
if (cmSub && cm>0)
return stat[iy][ix].getPedestal()+getCommonMode(ix,iy);
else return stat[iy][ix].getPedestal();
else return -1;
};
/**
gets pedestal rms (i.e. noise)
@ -278,9 +362,17 @@ template <class dataType> class analogDetector {
\param iy pixel y coordinate
\returns pedestal rms
*/
virtual double getPedestalRMS(int ix, int iy){if (ix>=0 && ix<nx && iy>=0 && iy<ny) return stat[iy][ix].getPedestalRMS();else return -1;};
virtual double getPedestalRMS(int ix, int iy){
if (ix>=0 && ix<nx && iy>=0 && iy<ny)
return stat[iy][ix].getPedestalRMS();
else return -1;
};
virtual int getNumpedestals(int ix, int iy){
if (ix>=0 && ix<nx && iy>=0 && iy<ny)
return stat[iy][ix].getNumpedestals();
else return -1;
};
/**
gets pedestal (and common mode)
\param ix pixel x coordinate
@ -327,7 +419,6 @@ template <class dataType> class analogDetector {
/**
sets pedestal
@ -396,17 +487,45 @@ template <class dataType> class analogDetector {
virtual void *writeImage(const char * imgname) {
float *gm=NULL;
void *ret;
#ifdef ROOTSPECTRUM
TH2F *hmap=new TH2F("hmap","hmap",nx, -0.5,nx-0.5, ny, -0.5, ny-0.5);
#endif
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];
#ifdef ROOTSPECTRUM
hmap->SetBinContent(ix+1, iy+1,image[iy*nx+ix]);
#endif
}
}
ret=WriteToTiff(gm, imgname, ny, nx);
delete [] gm;
#ifdef ROOTSPECTRUM
char rootfn[10000];
sprintf(rootfn,"%s.root",imgname);
TFile *f=new TFile(rootfn,"RECREATE");
hs->Write("hs");
#ifdef ROOTCLUST
hs3->Write("hs3");
hs5->Write("hs5");
hs7->Write("hs7");
hs9->Write("hs9");
#endif
hmap->Write("hmap");
f->Close();
delete f;
delete hmap;
#endif
return ret;
}
#ifdef ROOTSPECTRUM
TH2F *getSpectrum(){return hs;};
#endif
/**
write 32bit tiff file containing the pedestals
\param imgname file name to be written
@ -416,18 +535,43 @@ template <class dataType> class analogDetector {
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++) {
if (cmSub)
gm[iy*nx+ix]=stat[iy][ix].getPedestal()-cmSub->getCommonMode();
else
gm[iy*nx+ix]=stat[iy][ix].getPedestal();
}
gm=new float[nx*ny];
#ifdef ROOTSPECTRUM
TH2F *hmap=new TH2F("hmap","hmap",nx, -0.5,nx-0.5, ny, -0.5, ny-0.5);
#endif
for (int ix=0; ix<nx; ix++) {
for (int iy=0; iy<ny; iy++) {
/* if (cmSub) */
/* gm[iy*nx+ix]=stat[iy][ix].getPedestal()-cmSub->getCommonMode(); */
/* else */
gm[iy*nx+ix]=stat[iy][ix].getPedestal();
#ifdef ROOTSPECTRUM
hmap->SetBinContent(ix+1, iy+1,gm[iy*nx+ix]);
#endif
}
ret=WriteToTiff(gm, imgname, ny, nx);
delete [] gm;
return ret;
}
ret=WriteToTiff(gm, imgname, ny, nx);
delete [] gm;
#ifdef ROOTSPECTRUM
char rootfn[10000];
sprintf(rootfn,"%s.root",imgname);
TFile *f=new TFile(rootfn,"RECREATE");
hs->Write("hs");
#ifdef ROOTCLUST
hs3->Write("hs3");
hs5->Write("hs5");
hs7->Write("hs7");
hs9->Write("hs9");
#endif
hmap->Write("hmap");
f->Close();
delete f;
delete hmap;
#endif
return ret;
}
/**
@ -538,18 +682,24 @@ template <class dataType> class analogDetector {
\param data pointer to the data
*/
virtual void addToPedestal(char *data) {
virtual void addToPedestal(char *data, int cm=0) {
newFrame();
if (cmSub) {
addToCommonMode(data);
}
// cout << xmin << " " << xmax << endl;
// cout << ymin << " " << ymax << endl;
for (int ix=xmin; ix<xmax; ix++) {
for (int iy=ymin; iy<ymax; iy++) {
addToPedestal(data,ix,iy,1);
addToPedestal(data,ix,iy);
#ifdef ROOTSPECTRUM
subtractPedestal(data,ix,iy,cm);
#endif
}
}
@ -585,6 +735,20 @@ template <class dataType> class analogDetector {
}
#ifdef ROOTSPECTRUM
delete hs;
hs=new TH2F("hs","hs",2000,-100,10000,(xmax-xmin)*(ymax-ymin),-0.5,(xmax-xmin)*(ymax-ymin)-0.5);
#ifdef ROOTCLUST
delete hs3;
hs3=new TH2F("hs3","hs3",2000,-100,3*3*10000,(xmax-xmin)*(ymax-ymin),-0.5,(xmax-xmin)*(ymax-ymin)-0.5);
delete hs5;
hs5=new TH2F("hs5","hs5",2000,-100,5*5*10000,(xmax-xmin)*(ymax-ymin),-0.5,(xmax-xmin)*(ymax-ymin)-0.5);
delete hs7;
hs7=new TH2F("hs7","hs7",2000,-100,7*7*10000,(xmax-xmin)*(ymax-ymin),-0.5,(xmax-xmin)*(ymax-ymin)-0.5);
delete hs9;
hs9=new TH2F("hs9","hs9",2000,-100,9*9*10000,(xmax-xmin)*(ymax-ymin),-0.5,(xmax-xmin)*(ymax-ymin)-0.5);
#endif
#endif
};
/**
@ -608,7 +772,7 @@ template <class dataType> class analogDetector {
*/
virtual void addToPedestal(char *data, int ix, int iy=0) {
virtual void addToPedestal(char *data, int ix, int iy=0, int cm=0) {
double val;
@ -619,7 +783,15 @@ template <class dataType> class analogDetector {
else
val=((double*)data)[iy*nx+ix];
/* if (ix==10 && iy==10) */
/* cout << ix << " " << iy << " " << val ; */
/* if (ix==100 && iy==100) */
/* cout << ix << " " << iy << " " << val; */
addToPedestal(val,ix,iy);
/* if (ix==10 && iy==10) */
/* cout <<" " << getPedestal(ix,iy)<< endl; */
/* if (ix==100 && iy==100) */
/* cout << " " << getPedestal(ix,iy)<< endl; */
}
return ;
@ -631,19 +803,19 @@ template <class dataType> class analogDetector {
\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) {
*/
// virtual int *subtractPedestal(char *data, int *val=NULL) {
virtual int *subtractPedestal(char *data, int *val=NULL, int cm=0) {
newFrame();
if (val==NULL)
val=new double[nx*ny];
val=image;//new double[nx*ny];
for (int ix=xmin; ix<xmax; ix++) {
for (int iy=ymin; iy<ymax; iy++) {
val[iy*nx+ix]+=subtractPedestal(data, ix, iy);
val[iy*nx+ix]+=subtractPedestal(data, ix, iy,cm);
}
}
return val;
@ -662,8 +834,9 @@ template <class dataType> class analogDetector {
virtual double subtractPedestal(char *data, int ix, int iy=0) {
virtual double subtractPedestal(char *data, int ix, int iy=0, int cm=0) {
double g=1.;
double val;
if (ix>=0 && ix<nx && iy>=0 && iy<ny) {
if (gmap) {
g=gmap[iy*nx+ix];
@ -671,10 +844,41 @@ template <class dataType> class analogDetector {
}
if (det)
return (dataSign*det->getValue(data, ix, iy)-getPedestal(ix,iy))/g;
val= (dataSign*det->getValue(data, ix, iy)-getPedestal(ix,iy,cm))/g;
else
return (((double*)data)[iy*nx+ix]-getPedestal(ix,iy))/g;
}
val= (((double*)data)[iy*nx+ix]-getPedestal(ix,iy))/g;
#ifdef ROOTSPECTRUM
hs->Fill(val,(iy-ymin)*(xmax-xmin)+(ix-xmin));
#ifdef ROOTCLUST
double v3=0,v5=0,v7=0,v9=0;
for (int iix=-4; iix<5; iix++)
for (int iiy=-4; iiy<5; iiy++) {
if (det)
val= (dataSign*det->getValue(data, ix+iix, iy+iiy)-getPedestal(ix+iix,iy+iiy,cm))/g;
else
val= (((double*)data)[(iy+iiy)*nx+ix+iix]-getPedestal(ix+iix,iy+iiy,cm))/g;
if (iix>-4 && iiy>-4 && iix<4 && iiy<4) {
if (iix>-3 && iiy>-3 && iix<3 && iiy<3){
if (iix>-2 && iiy>-2 && iix<2 && iiy<2){
v3+=val;
}
v5+=val;
}
v7+=val;
}
v9+=val;
}
hs3->Fill(v3,(iy-ymin)*(xmax-xmin)+(ix-xmin));
hs5->Fill(v5,(iy-ymin)*(xmax-xmin)+(ix-xmin));
hs7->Fill(v7,(iy-ymin)*(xmax-xmin)+(ix-xmin));
hs9->Fill(v9,(iy-ymin)*(xmax-xmin)+(ix-xmin));
#endif
#endif
return val;
}
};
@ -703,12 +907,20 @@ template <class dataType> class analogDetector {
double v;
if (ix>=0 && ix<nx && iy>=0 && iy<ny) {
v=subtractPedestal(data,ix,iy);
/* // cout << v << " " ; */
/* #ifdef ROOTSPECTRUM */
/* // cout << (iy-ymin)*(xmax-xmin)+(ix-xmin) << endl; */
/* hs->Fill(v,(iy-ymin)*(xmax-xmin)+(ix-xmin)); */
/* #endif */
if (thr>0) {
v+=0.5*thr;
nph=v/thr;
return nph;
} else
return v;
if (nph>0)
return nph;
return 0;
}
return v;
}
return 0;
};
@ -722,15 +934,18 @@ template <class dataType> class analogDetector {
int *getNPhotons(char *data, int *nph=NULL) {
double val;
if (nph==NULL)
nph=image;
newFrame();
for (int ix=xmin; ix<xmax; ix++) {
for (int iy=ymin; iy<ymax; iy++) {
nph[iy*nx+ix]+=getNPhotons(data, ix, iy);
if (nph==NULL)
nph=image;
newFrame();
addToCommonMode(data);
for (int ix=xmin; ix<xmax; ix++) {
for (int iy=ymin; iy<ymax; iy++) {
nph[iy*nx+ix]+=getNPhotons(data, ix, iy);
}
}
}
return nph;
return nph;
}
@ -744,6 +959,23 @@ template <class dataType> class analogDetector {
image[iy*nx+ix]=0;
}
}
#ifdef ROOTSPECTRUM
//cout << "reset histogram " << endl;
if (hs)
hs->Reset();
#ifdef ROOTCLUST
if (hs3)
hs3->Reset();
if (hs5)
hs5->Reset();
if (hs7)
hs7->Reset();
if (hs9)
hs9->Reset();
#endif
//cout << "done " << endl;
#endif
};
/** sets/gets number of samples for moving average pedestal calculation
@ -806,10 +1038,12 @@ template <class dataType> class analogDetector {
virtual void processData(char *data,int *val=NULL) {
switch(fMode) {
case ePedestal:
// cout << "ped " << endl;
addToPedestal(data);
break;
default:
getNPhotons(data,val);
//subtractPedestal(data);
getNPhotons(data);
}
};
@ -831,13 +1065,13 @@ template <class dataType> class analogDetector {
\param f file pointer
\returns current file pointer
*/
FILE *setFilePointer(FILE *f){myFile=f; return myFile;};
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;};
FILE *getFilePointer(){return myFile;};
void setMutex(pthread_mutex_t *m){fm=m;};
protected:
slsDetectorData<dataType> *det; /**< slsDetectorData to be used */
@ -859,7 +1093,15 @@ template <class dataType> class analogDetector {
// 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 */
#ifdef ROOTSPECTRUM
TH2F *hs;
#ifdef ROOTCLUST
TH2F *hs3;
TH2F *hs5;
TH2F *hs7;
TH2F *hs9;
#endif
#endif
pthread_mutex_t *fm;
};

View File

@ -34,6 +34,9 @@ public:
#ifdef BCHIP074_BCHIP075
cout << "This is a bchip074-bchip075 system " << endl;
#endif
uint16_t **dMask;
@ -59,6 +62,11 @@ public:
dMap[0][ix] = 1280*2+2*offset+ipix*2;//dataSize-2-ix;//+2*offset;
// dMap[0][ix] = 2*ipix+offset*(imod+1)+1280*2*imod;
dMask[0][ix] = 0x0;
#ifdef BCHIP074_BCHIP075
int ibad=ix/2+1280*imod;
if ((ibad>=128*4 && ibad<128*5) || (ibad>=9*128 && ibad<10*128) || (ibad>=(1280+128*4) && ibad<ibad>=(1280+128*6)))
dataROIMask[0][ix]=0;
#endif
}
setDataMap(dMap);
@ -78,7 +86,7 @@ public:
int getFrameNumber(char *buff){if (offset>=sizeof(sls_detector_header)) return ((sls_detector_header*)buff)->frameNumber; return iframe;};//*((int*)(buff+5))&0xffffff;};
/**
gets the packets number (last packet is labelled with 0 and is replaced with 40)

View File

@ -0,0 +1,30 @@
#ifndef GOTTHARDDOUBLECOMMONMODESUBTRACTION_H
#define GOTTHARDDOUBLECOMMONMODESUBTRACTION_H
#include "commonModeSubtractionNew.h"
class gotthardDoubleModuleCommonModeSubtraction : public commonModeSubtraction {
/** @short class to calculate the common mode of the pedestals based on an approximated moving average*/
public:
/** constructor
\param nn number of samples for the moving average to calculate the average common mode
\param iroi number of regions on which one can calculate the common mode separately. Defaults to 1 i.e. whole detector
*/
gotthardDoubleModuleCommonModeSubtraction(int ns=3) : commonModeSubtraction(2, ns) {};
/**
gets the common mode ROI for pixel ix, iy
*/
virtual int getROI(int ix, int iy){return ix%2;};
};
#endif

View File

@ -1,215 +0,0 @@
#include <vector>
#include <string>
#include <sstream>
#include <iomanip>
#include <fstream>
#include <stdio.h>
//#include <deque>
//#include <list>
//#include <queue>
#include <fstream>
#include <cstring>
#include "gotthardModuleDataNew.h"
#include "gotthardDoubleModuleDataNew.h"
#include "singlePhotonDetector.h"
//#include "interpolatingDetector.h"
//#include "linearInterpolation.h"
#include "multiThreadedAnalogDetector.h"
#include <ctime>
#define NC 1280
#define NR 1
#include "tiffIO.h"
void *gotthardProcessFrame() {
if (argc < 3 ) {
cprintf(RED, "Help: ./trial [receive socket ip] [receive starting port number] [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 = 32*2*5000;//atoi(argv[3]);
// send parameters if any
char* socketip2 = 0;
uint32_t portnum2 = 0;
if (argc > 3) {
send = true;
socketip2 = argv[3];
portnum2 = atoi(argv[4]);
}
cout << "\nrx socket ip : " << socketip <<
"\nrx port num : " << portnum ;
if (send) {
cout << "\nsd socket ip : " << socketip2 <<
"\nsd port num : " << portnum2;
}
cout << endl;
char fname0[10000], fname1[10000];
char fformat[10000];
int fifosize=1000;
strcpy(fformat,"/external_pool/gotthard_data/datadir_gotthardI/bchip074075/20170731/Xray/xray_15kV_200uA_5us_d%d_f000000000000_0.raw");
sprintf(fname0,fformat,0,0);
sprintf(fname1,fformat,1,1);
int nthreads=3;
int nph, nph1;
// single_photon_hit clusters[NR*NC];
// cout << "hits "<< endl;
int etabins=550;
double etamin=-1, etamax=2;
int nsubpix=1;
float *etah=new float[etabins*etabins];
// cout << "etah "<< endl;
cout << "image size "<< nsubpix*nsubpix*NC*NR << endl;
int *heta, *himage;
gotthardModuleDataNew *decoder=new gotthardModuleDataNew();
gotthardDoubleModuleDataNew *det=new gotthardDoubleModuleDataNew();
// cout << "decoder "<< endl;
// etaInterpolationPosXY *interp=new etaInterpolationPosXY(NC, NR, nsubpix, etabins, etamin, etamax);
// cout << "interp "<< endl;
// linearInterpolation *interp=new linearInterpolation(NC, NR, nsubpix);
//noInterpolation *interp=new noInterpolation(NC, NR, nsubpix);
// interp->readFlatField("/scratch/eta_100.tiff",etamin,etamax);
// interpolatingDetector *filter0=new interpolatingDetector(decoder,interp, 5, 1, 0, 1000, 10);
// interpolatingDetector *filter1=new interpolatingDetector(decoder,interp, 5, 1, 0, 1000, 10);
//filter->readPedestals("/scratch/ped_100.tiff");
//cout << "filter "<< endl;
singlePhotonDetector *filter=new singlePhotonDetector(det,3, 5, 1, 0, 1000, 200);
filter->setFrameMode(eFrame);
char *buff;//[2*(48+1280*2)];
char *buff0;
char *buff1;
multiThreadedAnalogDetector *mt=new multiThreadedAnalogDetector(filter,nthreads,fifosize);
mt->setFrameMode(eFrame);
mt->StartThreads();
mt->popFree(buff);
buff0=buff;
buff1=buff+48+1280*2;
int photons[1280*2];
int nf=0;
int ok=0;
ifstream filebin0,filebin1;
std::time_t end_time;
int16_t dout[1280*2];
int iFrame=-1;
int np=-1;
nph=0;
nph1=0;
//int np;
int iph;
int data_ready=1;
int *image;
// filter->setROI(0,512,0,1);
filebin0.open((const char *)(fname0), ios::in | ios::binary);
filebin1.open((const char *)(fname1), ios::in | ios::binary);
if (filebin0.is_open() && filebin1.is_open()) {
cout << "Opened file " << fname0<< endl;
cout << "Opened file " << fname1<< endl;
// mt->setFrameMode(eFrame);
// mt->prepareInterpolation(ok);
// mt->StartThreads();
// mt->popFree(buff);
nf=0;
iFrame=-1;
while ((decoder->readNextFrame(filebin0, iFrame, np, buff0)) && (decoder->readNextFrame(filebin1, iFrame, np, buff1))) {
//filter->processData(buff, photons);
// cout << nf << " " << decoder->getFrameNumber(buff0) << " " << decoder->getFrameNumber(buff1) << " " << filter->getPhFrame() << " " << filter->getPhTot() << endl;
// for (int i=0; i<1280*2; i++) {
// filter->addToPedestal(buff,i,0);
// dout[i]=filter->subtractPedestal(buff,i,0);
// if (nf>10 && i<512)
// if (i%2)
// cout << nf << " " << i << " " << filter->getPedestal(i,0) << " " << det->getValue(buff,i,0) << " " << decoder->getValue(buff1,1280-1-i/2,0)<< " " << dout[i] << endl;
// else
// cout << nf << " " << i << " " << filter->getPedestal(i,0) << " " << det->getValue(buff,i,0) << " " << decoder->getValue(buff0,i/2,0)<< " " << dout[i] << endl;
// }
mt->pushData(buff);
mt->nextThread();
mt->popFree(buff);
buff0=buff;
buff1=buff+48+1280*2;
nf++;
// cout << id << " " << nf << endl;
if (nf%10000==0) {
while (mt->isBusy()) {;}
image=filter->getImage();
if (image) {
cout << nf << "*****************" << endl;
for (int i=0; i<512; i++) {
cout << image[i] << "\t";
}
cout << endl;
}
filter->clearImage();
std::time(&end_time);
cout << std::ctime(&end_time) << " " << nf << endl;
}
iFrame=-1;
}
filebin0.close();
filebin1.close();
}
else
cout << "Could not open file " << fname0<< " or " << fname1 << endl;
return NULL;
}
int main(int argc, char *argv[]){
gotthardProcessFrame();
}

View File

@ -14,8 +14,11 @@
#include <fstream>
#include <cstring>
//#define BCHIP074_BCHIP075
#include "gotthardModuleDataNew.h"
#include "gotthardDoubleModuleDataNew.h"
#include "gotthardDoubleModuleCommonModeSubtractionNew.h"
#include "singlePhotonDetector.h"
//#include "interpolatingDetector.h"
@ -27,8 +30,7 @@
#define NC 1280
#define NR 1
#include "tiffIO.h"
//#include "tiffIO.h"
@ -40,19 +42,6 @@ int main(int argc, char *argv[]){
int fifosize=1000;
int nthreads=1;
int nph, nph1;
@ -69,16 +58,22 @@ int main(int argc, char *argv[]){
offset=48;
#endif
//commonModeSubtraction *cm=NULL;
gotthardDoubleModuleCommonModeSubtraction *cm=new gotthardDoubleModuleCommonModeSubtraction();
gotthardModuleDataNew *decoder=new gotthardModuleDataNew();
gotthardDoubleModuleDataNew *det=new gotthardDoubleModuleDataNew(offset);
singlePhotonDetector *filter=new singlePhotonDetector(det,3, 5, 1, 0, 1000, 100);
singlePhotonDetector *filter=new singlePhotonDetector(det,3, 5, 1, cm, 1000, 100);
// analogDetector<uint16_t> *filter=new analogDetector<uint16_t>(det, 1, cm, 1000);
// analogDetector<uint16_t> *filter_nocm=new analogDetector<uint16_t>(det, 1, NULL, 1000);
filter->setROI(0,2560,0,1);
char *buff;//[2*(48+1280*2)];
char *buff0;
char *buff1;
multiThreadedAnalogDetector *mt=new multiThreadedAnalogDetector(filter,nthreads,fifosize);
mt->setFrameMode(eFrame);
// mt->setFrameMode(eFrame);
// mt->setFrameMode(ePedestal);
mt->StartThreads();
mt->popFree(buff);
buff0=buff;
@ -210,7 +205,8 @@ int main(int argc, char *argv[]){
string filename0 = "";
int eoa0=0;
int eoa1=0;
@ -230,7 +226,9 @@ int main(int argc, char *argv[]){
char ofname[10000];
char fn0[10000], fn1[10000];
FILE *fout=NULL;
FILE *fclust=NULL;
for (int i=0; i<nnx; i++)
dout[i]=0;
char fname0[10000], fname1[10000];
@ -279,89 +277,204 @@ int main(int argc, char *argv[]){
int end_of_acquisition;
while(1) {
end_of_acquisition=0;
eoa0=0;
eoa1=0;
// cout << "Receive header " << nf << endl;
if (!zmqsocket0->ReceiveHeader(0, acqIndex0, frameIndex0, subframeIndex0, filename0, fileindex0)) {
cout << "************************************************************************** packet0!*****************************"<< endl;
// cout << "************************************************************************** packet0!*****************************"<< endl;
eoa0=1;
end_of_acquisition++;
}
if (!zmqsocket1->ReceiveHeader(0, acqIndex1, frameIndex1, subframeIndex1, filename1, fileindex1)) {
cout << "************************************************************************** packet1!*****************************"<< endl;
//cout << "************************************************************************** packet1!*****************************"<< endl;
eoa1=1;
end_of_acquisition++;
}
// if ((!zmqsocket0->ReceiveHeader(0, acqIndex0, frameIndex0, subframeIndex0, filename0, fileindex0)) && (!zmqsocket1->ReceiveHeader(0, acqIndex1, frameIndex1, subframeIndex1, filename1, fileindex1))){
if (end_of_acquisition) {
cout << "************************************************************************** END OF FRAME" << end_of_acquisition << " !*****************************"<< endl;
// return 0;
// if ((!zmqsocket0->ReceiveHeader(0, acqIndex0, frameIndex0, subframeIndex0, filename0, fileindex0)) && (!zmqsocket1->ReceiveHeader(0, acqIndex1, frameIndex1, subframeIndex1, filename1, fileindex1))){
if (end_of_acquisition==0) {
if (acqIndex0!=acqIndex1)
cout << "different acquisition indexes " << acqIndex0 << " and " << acqIndex1 << endl;
if (frameIndex0!=frameIndex1)
cout << "different frame indexes " << frameIndex0 << " and " << frameIndex1 << endl;
while (frameIndex0<frameIndex1) {
cout << "aligning det 0 " << endl;
length = zmqsocket0->ReceiveData(0, buff0, size/2);
if (!zmqsocket0->ReceiveHeader(0, acqIndex0, frameIndex0, subframeIndex0, filename0, fileindex0)) {
end_of_acquisition++;
eoa0=1;
break;
}
}
while (frameIndex1<frameIndex0) {
cout << "aligning det 1 " << endl;
length = zmqsocket1->ReceiveData(0, buff1, size/2);
if (!zmqsocket1->ReceiveHeader(0, acqIndex1, frameIndex1, subframeIndex1, filename1, fileindex1)) {
end_of_acquisition++;
eoa1=1;
break;
}
}
}
if (eoa0!=eoa1) {
while (eoa0<1) {
length = zmqsocket0->ReceiveData(0, buff0, size/2);
if (!zmqsocket0->ReceiveHeader(0, acqIndex0, frameIndex0, subframeIndex0, filename0, fileindex0)) {
end_of_acquisition++;
eoa0=1;
}
}
while (eoa1<1) {
length = zmqsocket1->ReceiveData(0, buff1, size/2);
if (!zmqsocket1->ReceiveHeader(0, acqIndex1, frameIndex1, subframeIndex1, filename1, fileindex1)) {
end_of_acquisition++;
eoa1=1;
}
}
}
if (end_of_acquisition) {
// cout << "************************************************************************** END OF FRAME" << end_of_acquisition << " !*****************************"<< endl;
// return 0;
sprintf(ofname,"%s_%d.ph",fn0,irun);
while (mt->isBusy()) {;}
image=filter->getImage();
if (image) {
//fout=fopen(ofname,"w");
cout << nf << "*****************" << endl;
fout=fopen(ofname,"w");
cout << nf << "*****************" << endl;
for (int i=0; i<2560; i++) {
// // fprintf(fout,"%d %d\n",i,image[i]);
dout[i]=image[i];
fprintf(fout,"%d %d\n",i,image[i]);
dout[i]=image[i];
if (dout[i]<0)
dout[i]=0;
}
// // fclose(fout);;
fclose(fout);
}
if (send) {
strcpy(fname0,filename0.c_str());
strcpy(fname1,filename1.c_str());
// zmqsocket2->SendHeaderData(0, false, SLS_DETECTOR_JSON_HEADER_VERSION,16,fileindex,400,400,400*400, acqIndex,frameIndex,fname, acqIndex, 0,0,0,0,0,0,0,0,0,0,0,1);
zmqsocket2->SendHeaderData(0, false, SLS_DETECTOR_JSON_HEADER_VERSION,0,0,0,0,0, 0,0,fname0, 0, 0,0,0,0,0,0,0,0,0,0,0,1);
zmqsocket3->SendHeaderData(0, false, SLS_DETECTOR_JSON_HEADER_VERSION,0,0,0,0,0, 0,0,fname1, 0, 0,0,0,0,0,0,0,0,0,0,0,1);
zmqsocket2->SendData((char*)dout,size/2);
zmqsocket3->SendData(((char*)dout)+size/2,size/2);
// cprintf(GREEN, "Sent Data\n");
zmqsocket2->SendHeaderData(0, false, SLS_DETECTOR_JSON_HEADER_VERSION,0,0,0,0,0, 0,0,fn0, 0, 0,0,0,0,0,0,0,0,0,0,0,1);
zmqsocket3->SendHeaderData(0, false, SLS_DETECTOR_JSON_HEADER_VERSION,0,0,0,0,0, 0,0,fn1, 0, 0,0,0,0,0,0,0,0,0,0,0,1);
zmqsocket2->SendData((char*)dout,size/2);
zmqsocket3->SendData(((char*)dout)+size/2,size/2);
// // cprintf(GREEN, "Sent Data\n");
zmqsocket2->SendHeaderData(0, true, SLS_DETECTOR_JSON_HEADER_VERSION);
zmqsocket3->SendHeaderData(0, true, SLS_DETECTOR_JSON_HEADER_VERSION);
cprintf(RED, "Received %d frames\n", nf);
}
//mt->setFrameMode(eFrame);
filter->clearImage();
// std::time(&end_time);
// cout << std::ctime(&end_time) << " " << nf << endl;
fclose(fclust);
fclust=NULL;
continue;
}
if (acqIndex0!=acqIndex1)
cout << "different acquisition indexes " << acqIndex0 << " and " << acqIndex1 << endl;
if (frameIndex0!=frameIndex1)
cout << "different frame indexes " << frameIndex0 << " and " << frameIndex1 << endl;
if (fclust==NULL) {
strcpy(fn0,filename0.c_str());
strcpy(fn1,filename1.c_str());
sprintf(ofname,"%s_%d.clust",fn0,irun);
fclust=fopen(ofname,"w");
while (mt->isBusy()) {;}
mt->setFilePointer(fclust);
}
// strcpy(fn0,filename0.c_str());
// strcpy(fn1,filename1.c_str());
// cout << "Receive data " << nf << endl;
length = zmqsocket0->ReceiveData(0, buff0, size/2);
length += zmqsocket1->ReceiveData(0, buff1, size/2);
irun=fileindex0;
sprintf(ofname,"%s_%d.ph",filename0.c_str(),fileindex0);
// // if (nf>100)
// // mt->setFrameMode(eFrame);
// //filter->clearImage();
#endif
mt->pushData(buff);
mt->nextThread();
// cout << "==" << nf << endl;
// while (mt->isBusy()) {;}
// image=filter->getImage();
// if (image) {
// for (int i=0; i<2560; i++) {
// // if (i<512)
// // fprintf(fout,"%d %d\n",i,image[i]);
// dout[i]=filter->subtractPedestal(buff,i,0,1);//image[i];//filter->getPedestal(i,0);//
// if (dout[i]<0)
// dout[i]=0;
// // cout << i << " " << image[i] << " " << dout[i] << endl;
// }
// }
// if (send) {
// strcpy(fname0,filename0.c_str());
// strcpy(fname1,filename1.c_str());
// // zmqsocket2->SendHeaderData(0, false, SLS_DETECTOR_JSON_HEADER_VERSION,16,fileindex,400,400,400*400, acqIndex,frameIndex,fname, acqIndex, 0,0,0,0,0,0,0,0,0,0,0,1);
// zmqsocket2->SendHeaderData(0, false, SLS_DETECTOR_JSON_HEADER_VERSION,0,0,0,0,0, 0,0,fname0, 0, 0,0,0,0,0,0,0,0,0,0,0,1);
// zmqsocket3->SendHeaderData(0, false, SLS_DETECTOR_JSON_HEADER_VERSION,0,0,0,0,0, 0,0,fname1, 0, 0,0,0,0,0,0,0,0,0,0,0,1);
// zmqsocket2->SendData((char*)dout,size/2);
// zmqsocket3->SendData(((char*)dout)+size/2,size/2);
// // cprintf(GREEN, "Sent Data\n");
// // zmqsocket2->SendHeaderData(0, true, SLS_DETECTOR_JSON_HEADER_VERSION);
// // zmqsocket3->SendHeaderData(0, true, SLS_DETECTOR_JSON_HEADER_VERSION);
// // cprintf(RED, "Received %d frames\n", nf);
// }
mt->popFree(buff);
buff0=buff;
buff1=buff+offset*2+1280*2;
nf++;
@ -370,6 +483,8 @@ int main(int argc, char *argv[]){
#ifndef ZMQ
while (mt->isBusy()) {;}
image=filter->getImage();
if (image) {
@ -378,7 +493,7 @@ int main(int argc, char *argv[]){
for (int i=0; i<512; i++) {
fprintf(fout,"%d %d\n",i,image[i]);
}
fclose(fout);;
fclose(fout);
}
filter->clearImage();

View File

@ -0,0 +1,436 @@
#ifndef ETA2_INTERPOLATION_BASE_H
#define ETA2_INTERPOLATION_BASE_H
#ifdef MYROOT1
#include <TObject.h>
#include <TTree.h>
#include <TH2D.h>
#include <TH2F.h>
#endif
#include "etaInterpolationBase.h"
class eta2InterpolationBase : public virtual etaInterpolationBase {
public:
eta2InterpolationBase(int nx=400, int ny=400, int ns=25, int nb=-1, double emin=1, double emax=0) : etaInterpolationBase(nx,ny, ns, nb, emin, emax) {
// cout << "e2ib " << nb << " " << emin << " " << emax << endl;
if (nbeta<=0) {
nbeta=nSubPixels*10;
}
if (etamin>=etamax) {
etamin=-1;
etamax=2;
cout << ":" <<endl;
}
etastep=(etamax-etamin)/nbeta;
#ifdef MYROOT1
delete heta;
delete hhx;
delete hhy;
heta=new TH2D("heta","heta",nbeta,etamin,etamax,nbeta,etamin,etamax);
hhx=new TH2D("hhx","hhx",nbeta,etamin,etamax,nbeta,etamin,etamax);
hhy=new TH2D("hhy","hhy",nbeta,etamin,etamax,nbeta,etamin,etamax);
#endif
#ifndef MYROOT1
delete [] heta;
delete [] hhx;
delete [] hhy;
heta=new int[nbeta*nbeta];
hhx=new float[nbeta*nbeta];
hhy=new float[nbeta*nbeta];
#endif
// cout << nbeta << " " << etamin << " " << etamax << endl;
};
eta2InterpolationBase(eta2InterpolationBase *orig): etaInterpolationBase(orig){ };
virtual eta2InterpolationBase* Clone()=0;/* {
return new eta2InterpolationBase(this);
};
*/
//////////////////////////////////////////////////////////////////////////////
//////////// /*It return position hit for the event in input */ //////////////
virtual void getInterpolatedPosition(int x, int y, int *data, double &int_x, double &int_y)
{
double sDum[2][2];
double tot, totquad;
double etax,etay;
int corner;
corner=calcQuad(data, tot, totquad, sDum);
if (nSubPixels>2)
calcEta(totquad, sDum, etax, etay);
getInterpolatedPosition(x,y,etax,etay,corner,int_x,int_y);
return;
};
virtual void getInterpolatedPosition(int x, int y, double *data, double &int_x, double &int_y)
{
double sDum[2][2];
double tot, totquad;
double etax,etay;
int corner;
corner=calcQuad(data, tot, totquad, sDum);
if (nSubPixels>2)
calcEta(totquad, sDum, etax, etay);
getInterpolatedPosition(x,y,etax,etay,corner,int_x,int_y);
return;
};
virtual void getInterpolatedPosition(int x, int y, double totquad,int quad,double *cl,double &int_x, double &int_y) {
double cc[2][2];
double *cluster[3];
int xoff, yoff;
cluster[0]=cl;
cluster[1]=cl+3;
cluster[2]=cl+6;
switch (quad) {
case BOTTOM_LEFT:
xoff=0;
yoff=0;
break;
case BOTTOM_RIGHT:
xoff=1;
yoff=0;
break;
case TOP_LEFT:
xoff=0;
yoff=1;
break;
case TOP_RIGHT:
xoff=1;
yoff=1;
break;
default:
;
}
double etax, etay;
if (nSubPixels>2) {
cc[0][0]=cluster[yoff][xoff];
cc[1][0]=cluster[yoff+1][xoff];
cc[0][1]=cluster[yoff][xoff+1];
cc[1][1]=cluster[yoff+1][xoff+1];
calcEta(totquad,cc,etax,etay);
}
return getInterpolatedPosition(x,y,etax, etay,quad,int_x,int_y);
}
virtual void getInterpolatedPosition(int x, int y, double totquad,int quad,int *cl,double &int_x, double &int_y) {
double cc[2][2];
int *cluster[3];
int xoff, yoff;
cluster[0]=cl;
cluster[1]=cl+3;
cluster[2]=cl+6;
switch (quad) {
case BOTTOM_LEFT:
xoff=0;
yoff=0;
break;
case BOTTOM_RIGHT:
xoff=1;
yoff=0;
break;
case TOP_LEFT:
xoff=0;
yoff=1;
break;
case TOP_RIGHT:
xoff=1;
yoff=1;
break;
default:
;
}
double etax, etay;
if (nSubPixels>2) {
cc[0][0]=cluster[yoff][xoff];
cc[1][0]=cluster[yoff+1][xoff];
cc[0][1]=cluster[yoff][xoff+1];
cc[1][1]=cluster[yoff+1][xoff+1];
calcEta(totquad,cc,etax,etay);
}
return getInterpolatedPosition(x,y,etax, etay,quad,int_x,int_y);
}
virtual void getInterpolatedPosition(int x, int y, double etax, double etay, int corner, double &int_x, double &int_y)
{
double xpos_eta=0,ypos_eta=0;
double dX,dY;
int ex,ey;
switch (corner)
{
case TOP_LEFT:
dX=-1.;
dY=0;
break;
case TOP_RIGHT:
;
dX=0;
dY=0;
break;
case BOTTOM_LEFT:
dX=-1.;
dY=-1.;
break;
case BOTTOM_RIGHT:
dX=0;
dY=-1.;
break;
default:
cout << "bad quadrant" << endl;
dX=0.;
dY=0.;
}
if (nSubPixels>2) {
#ifdef MYROOT1
xpos_eta=(hhx->GetBinContent(hhx->GetXaxis()->FindBin(etax),hhy->GetYaxis()->FindBin(etay)))/((double)nSubPixels);
ypos_eta=(hhy->GetBinContent(hhx->GetXaxis()->FindBin(etax),hhy->GetYaxis()->FindBin(etay)))/((double)nSubPixels);
#endif
#ifndef MYROOT1
ex=(etax-etamin)/etastep;
ey=(etay-etamin)/etastep;
if (ex<0) {
cout << "x*"<< ex << endl;
ex=0;
}
if (ex>=nbeta) {
cout << "x?"<< ex << endl;
ex=nbeta-1;
}
if (ey<0) {
cout << "y*"<< ey << endl;
ey=0;
}
if (ey>=nbeta) {
cout << "y?"<< ey << endl;
ey=nbeta-1;
}
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=0.5*dX+0.25;
ypos_eta=0.5*dY+0.25;
}
int_x=((double)x) + xpos_eta+0.5;
int_y=((double)y) + ypos_eta+0.5;
}
virtual int addToFlatField(double totquad,int quad,int *cl,double &etax, double &etay) {
double cc[2][2];
int *cluster[3];
int xoff, yoff;
cluster[0]=cl;
cluster[1]=cl+3;
cluster[2]=cl+6;
switch (quad) {
case BOTTOM_LEFT:
xoff=0;
yoff=0;
break;
case BOTTOM_RIGHT:
xoff=1;
yoff=0;
break;
case TOP_LEFT:
xoff=0;
yoff=1;
break;
case TOP_RIGHT:
xoff=1;
yoff=1;
break;
default:
;
}
cc[0][0]=cluster[yoff][xoff];
cc[1][0]=cluster[yoff+1][xoff];
cc[0][1]=cluster[yoff][xoff+1];
cc[1][1]=cluster[yoff+1][xoff+1];
/* cout << cl[0] << " " << cl[1] << " " << cl[2] << endl; */
/* cout << cl[3] << " " << cl[4] << " " << cl[5] << endl; */
/* cout << cl[6] << " " << cl[7] << " " << cl[8] << endl; */
/* cout <<"******"<<totquad << " " << quad << endl; */
/* cout << cc[0][0]<< " " << cc[0][1] << endl; */
/* cout << cc[1][0]<< " " << cc[1][1] << endl; */
//calcMyEta(totquad,quad,cl,etax, etay);
calcEta(totquad, cc,etax, etay);
// cout <<"******"<< etax << " " << etay << endl;
return addToFlatField(etax,etay);
}
virtual int addToFlatField(double totquad,int quad,double *cl,double &etax, double &etay) {
double cc[2][2];
double *cluster[3];
int xoff, yoff;
cluster[0]=cl;
cluster[1]=cl+3;
cluster[2]=cl+6;
switch (quad) {
case BOTTOM_LEFT:
xoff=0;
yoff=0;
break;
case BOTTOM_RIGHT:
xoff=1;
yoff=0;
break;
case TOP_LEFT:
xoff=0;
yoff=1;
break;
case TOP_RIGHT:
xoff=1;
yoff=1;
break;
default:
;
}
cc[0][0]=cluster[yoff][xoff];
cc[1][0]=cluster[yoff+1][xoff];
cc[0][1]=cluster[yoff][xoff+1];
cc[1][1]=cluster[yoff+1][xoff+1];
/* cout << cl[0] << " " << cl[1] << " " << cl[2] << endl; */
/* cout << cl[3] << " " << cl[4] << " " << cl[5] << endl; */
/* cout << cl[6] << " " << cl[7] << " " << cl[8] << endl; */
/* cout <<"******"<<totquad << " " << quad << endl; */
/* cout << cc[0][0]<< " " << cc[0][1] << endl; */
/* cout << cc[1][0]<< " " << cc[1][1] << endl; */
//calcMyEta(totquad,quad,cl,etax, etay);
calcEta(totquad, cc,etax, etay);
// cout <<"******"<< etax << " " << etay << endl;
return addToFlatField(etax,etay);
}
//////////////////////////////////////////////////////////////////////////////////////
virtual int addToFlatField(double *cluster, double &etax, double &etay){
double sDum[2][2];
double tot, totquad;
int corner;
corner=calcQuad(cluster, tot, totquad, sDum);
double xpos_eta,ypos_eta;
double dX,dY;
calcEta(totquad, sDum, etax, etay);
return addToFlatField(etax,etay);
};
virtual int addToFlatField(int *cluster, double &etax, double &etay){
double sDum[2][2];
double tot, totquad;
int corner;
corner=calcQuad(cluster, tot, totquad, sDum);
double xpos_eta,ypos_eta;
double dX,dY;
calcEta(totquad, sDum, etax, etay);
return addToFlatField(etax,etay);
};
virtual int addToFlatField(double etax, double etay){
#ifdef MYROOT1
heta->Fill(etax,etay);
#endif
#ifndef MYROOT1
int ex,ey;
ex=(etax-etamin)/etastep;
ey=(etay-etamin)/etastep;
if (ey<nbeta && ex<nbeta && ex>=0 && ey>=0)
heta[ey*nbeta+ex]++;
#endif
return 0;
};
/* protected: */
/* #ifdef MYROOT1 */
/* TH2D *heta; */
/* TH2D *hhx; */
/* TH2D *hhy; */
/* #endif */
/* #ifndef MYROOT1 */
/* int *heta; */
/* float *hhx; */
/* float *hhy; */
/* #endif */
/* int nbeta; */
/* double etamin, etamax, etastep; */
};
#endif

View File

@ -0,0 +1,294 @@
#ifndef ETA3_INTERPOLATION_BASE_H
#define ETA3_INTERPOLATION_BASE_H
#ifdef MYROOT1
#include <TObject.h>
#include <TTree.h>
#include <TH2D.h>
#include <TH2F.h>
#endif
#include "etaInterpolationBase.h"
class eta3InterpolationBase : public virtual etaInterpolationBase {
public:
eta3InterpolationBase(int nx=400, int ny=400, int ns=25, int nb=-1, double emin=1, double emax=0) : etaInterpolationBase(nx, ny, ns, nb, emin, emax) {
// cout << "e3ib " << nb << " " << emin << " " << emax << endl;
if (nbeta<=0) {
nbeta=nSubPixels*10;
}
if (etamin>=etamax) {
etamin=-1;
etamax=1;
}
etastep=(etamax-etamin)/nbeta;
#ifdef MYROOT1
delete heta;
delete hhx;
delete hhy;
heta=new TH2D("heta","heta",nbeta,etamin,etamax,nbeta,etamin,etamax);
hhx=new TH2D("hhx","hhx",nbeta,etamin,etamax,nbeta,etamin,etamax);
hhy=new TH2D("hhy","hhy",nbeta,etamin,etamax,nbeta,etamin,etamax);
#endif
#ifndef MYROOT1
delete [] heta;
delete [] hhx;
delete [] hhy;
heta=new int[nbeta*nbeta];
hhx=new float[nbeta*nbeta];
hhy=new float[nbeta*nbeta];
#endif
// cout << nbeta << " " << etamin << " " << etamax << endl;
};
eta3InterpolationBase(eta3InterpolationBase *orig): etaInterpolationBase(orig){ };
virtual eta3InterpolationBase* Clone()=0;
// virtual void prepareInterpolation(int &ok){};
//////////////////////////////////////////////////////////////////////////////
//////////// /*It return position hit for the event in input */ //////////////
virtual void getInterpolatedPosition(int x, int y, int *data, double &int_x, double &int_y)
{
double tot, totquad;
double etax,etay;
int corner=calcEta3(data,etax,etay, totquad);
getInterpolatedPosition(x,y,etax,etay,corner,int_x,int_y);
return;
};
virtual void getInterpolatedPosition(int x, int y, double *data, double &int_x, double &int_y)
{
double sDum[2][2];
double tot, totquad;
double etax,etay;
int corner=calcEta3(data,etax,etay, totquad);
getInterpolatedPosition(x,y,etax,etay,corner,int_x,int_y);
return;
};
virtual void getInterpolatedPosition(int x, int y, double totquad,int quad,double *cl,double &int_x, double &int_y) {
double etax, etay;
if (nSubPixels>2) {
calcEta3(cl,etax,etay, totquad);
}
return getInterpolatedPosition(x,y,etax, etay,quad,int_x,int_y);
}
virtual void getInterpolatedPosition(int x, int y, double totquad,int quad,int *cl,double &int_x, double &int_y) {
double etax, etay;
if (nSubPixels>2) {
calcEta3(cl,etax,etay, totquad);
}
return getInterpolatedPosition(x,y,etax, etay,quad,int_x,int_y);
}
virtual void getInterpolatedPosition(int x, int y, double etax, double etay, int corner, double &int_x, double &int_y)
{
double xpos_eta=0,ypos_eta=0;
int ex,ey;
if (nSubPixels>2) {
#ifdef MYROOT1
xpos_eta=(hhx->GetBinContent(hhx->GetXaxis()->FindBin(etax),hhy->GetYaxis()->FindBin(etay)))/((double)nSubPixels);
ypos_eta=(hhy->GetBinContent(hhx->GetXaxis()->FindBin(etax),hhy->GetYaxis()->FindBin(etay)))/((double)nSubPixels);
#endif
#ifndef MYROOT1
ex=(etax-etamin)/etastep;
ey=(etay-etamin)/etastep;
if (ex<0) {
/* cout << etax << " " << etamin << " "; */
/* cout << "3x*"<< ex << endl; */
ex=0;
}
if (ex>=nbeta) {
/* cout << etax << " " << etamin << " "; */
/* cout << "3x?"<< ex << endl; */
ex=nbeta-1;
}
if (ey<0) {
/* cout << etay << " " << etamin << " "; */
/* cout << "3y*"<< ey << endl; */
ey=0;
}
if (ey>=nbeta) {
/* cout << etay << " " << etamin << " "; */
/* cout << "3y?"<< ey << endl; */
ey=nbeta-1;
}
xpos_eta=(((double)hhx[(ey*nbeta+ex)]));///((double)nSubPixels);
ypos_eta=(((double)hhy[(ey*nbeta+ex)]));///((double)nSubPixels);
#endif
} else {
switch (corner) {
case BOTTOM_LEFT:
xpos_eta=-0.25;
ypos_eta=-0.25;
break;
case BOTTOM_RIGHT:
xpos_eta=0.25;
ypos_eta=-0.25;
break;
case TOP_LEFT:
xpos_eta=-0.25;
ypos_eta=0.25;
break;
case TOP_RIGHT:
xpos_eta=0.25;
ypos_eta=0.25;
break;
default:
xpos_eta=0;
ypos_eta=0;
}
}
int_x=((double)x) + xpos_eta;
int_y=((double)y) + ypos_eta;
// int_x=5. + xpos_eta;
// int_y=5. + ypos_eta;
}
/* ///////////////////////////////////////////////////////////////////////////////////////////////// */
/* virtual void getPositionETA3(int x, int y, double *data, double &int_x, double &int_y) */
/* { */
/* double sDum[2][2]; */
/* double tot, totquad; */
/* double eta3x,eta3y; */
/* double ex,ey; */
/* calcQuad(data, tot, totquad, sDum); */
/* calcEta3(data,eta3x, eta3y,tot); */
/* double xpos_eta,ypos_eta; */
/* #ifdef MYROOT1 */
/* xpos_eta=((hhx->GetBinContent(hhx->GetXaxis()->FindBin(eta3x),hhy->GetYaxis()->FindBin(eta3y))))/((double)nSubPixels); */
/* ypos_eta=((hhy->GetBinContent(hhx->GetXaxis()->FindBin(eta3x),hhy->GetYaxis()->FindBin(eta3y))))/((double)nSubPixels); */
/* #endif */
/* #ifndef MYROOT1 */
/* ex=(eta3x-etamin)/etastep; */
/* ey=(eta3y-etamin)/etastep; */
/* if (ex<0) ex=0; */
/* if (ex>=nbeta) ex=nbeta-1; */
/* if (ey<0) ey=0; */
/* if (ey>=nbeta) ey=nbeta-1; */
/* xpos_eta=(((double)hhx[(int)(ey*nbeta+ex)]))/((double)nSubPixels); */
/* ypos_eta=(((double)hhy[(int)(ey*nbeta+ex)]))/((double)nSubPixels); */
/* #endif */
/* int_x=((double)x) + xpos_eta; */
/* int_y=((double)y) + ypos_eta; */
/* return; */
/* }; */
virtual int addToFlatField(double totquad,int quad,int *cl,double &etax, double &etay) {
calcEta3(cl, etax, etay, totquad);
return addToFlatField(etax,etay);
}
virtual int addToFlatField(double totquad,int quad,double *cl,double &etax, double &etay) {
calcEta3(cl, etax, etay, totquad);
return addToFlatField(etax,etay);
}
//////////////////////////////////////////////////////////////////////////////////////
virtual int addToFlatField(double *cluster, double &etax, double &etay){
double totquad;
calcEta3(cluster, etax, etay, totquad);
return addToFlatField(etax,etay);
};
virtual int addToFlatField(int *cluster, double &etax, double &etay){
double totquad;
calcEta3(cluster, etax, etay, totquad);
return addToFlatField(etax,etay);
};
virtual int addToFlatField(double etax, double etay){
#ifdef MYROOT1
heta->Fill(etax,etay);
#endif
#ifndef MYROOT1
int ex,ey;
ex=(etax-etamin)/etastep;
ey=(etay-etamin)/etastep;
if (ey<nbeta && ex<nbeta && ex>=0 && ey>=0)
heta[ey*nbeta+ex]++;
#endif
return 0;
};
/* protected: */
/* #ifdef MYROOT1 */
/* TH2D *heta; */
/* TH2D *hhx; */
/* TH2D *hhy; */
/* #endif */
/* #ifndef MYROOT1 */
/* int *heta; */
/* float *hhx; */
/* float *hhy; */
/* #endif */
/* int nbeta; */
/* double etamin, etamax, etastep; */
};
#endif

View File

@ -181,53 +181,73 @@ class etaInterpolationAdaptiveBins : public etaInterpolationPosXY {
double avg=tot_eta/((double)(nSubPixels*nSubPixels));
cout << "total eta entries is :"<< tot_eta << " avg: "<< avg << endl;
cout << "Start " << endl;
double old_diff=calcDiff(avg, hhx, hhy), new_diff=old_diff+1;
cout << " diff= " << new_diff << endl;
double old_diff=calcDiff(avg, hhx, hhy), new_diff=old_diff+1, best_diff=old_diff;
cout << " diff= " << old_diff << endl;
int iint=0;
float *newhhx=new float[nbeta*nbeta]; //profile x
float *newhhy=new float[nbeta*nbeta]; //profile y
float *besthhx=hhx; //profile x
float *besthhy=hhy; //profile y
while (iint<nint) {
cout << "Iteration " << iint << endl;
iterate(newhhx,newhhy);
new_diff=calcDiff(avg, newhhx, newhhy);
cout << " diff= " << new_diff << endl;
#ifdef SAVE_ALL
for (int ii=0; ii<etabins*etabins; ii++) {
etah[ii]=newhhx[ii];
if (etah[ii]>1 || etah[ii]<0 ) cout << "***"<< ii << etah[ii] << endl;
/* #ifdef SAVE_ALL */
/* for (int ii=0; ii<etabins*etabins; ii++) { */
/* etah[ii]=newhhx[ii]; */
/* if (etah[ii]>1 || etah[ii]<0 ) cout << "***"<< ii << etah[ii] << endl; */
}
sprintf(tit,"/scratch/neweta_hhx_%d.tiff",iint);
WriteToTiff(etah, tit, etabins, etabins);
/* } */
/* sprintf(tit,"/scratch/neweta_hhx_%d.tiff",iint); */
/* WriteToTiff(etah, tit, etabins, etabins); */
for (int ii=0; ii<etabins*etabins; ii++) {
etah[ii]=newhhy[ii];
if (etah[ii]>1 || etah[ii]<0 ) cout << "***"<< ii << etah[ii] << endl;
/* for (int ii=0; ii<etabins*etabins; ii++) { */
/* etah[ii]=newhhy[ii]; */
/* if (etah[ii]>1 || etah[ii]<0 ) cout << "***"<< ii << etah[ii] << endl; */
/* } */
/* sprintf(tit,"/scratch/neweta_hhy_%d.tiff",iint); */
/* WriteToTiff(etah, tit, etabins, etabins); */
/* #endif */
if (new_diff<best_diff) {
best_diff=new_diff;
besthhx=newhhx;
besthhy=newhhy;
}
sprintf(tit,"/scratch/neweta_hhy_%d.tiff",iint);
WriteToTiff(etah, tit, etabins, etabins);
#endif
// if (new_diff<old_diff) {
delete [] hhx;
delete [] hhy;
hhx=newhhx;
hhy=newhhy;
newhhx=new float[nbeta*nbeta]; //profile x
newhhy=new float[nbeta*nbeta]; //profile y
old_diff=new_diff;
// } /* else { */
if (hhx!=besthhx)
delete [] hhx;
if (hhy!=besthhy)
delete [] hhy;
hhx=newhhx;
hhy=newhhy;
newhhx=new float[nbeta*nbeta]; //profile x
newhhy=new float[nbeta*nbeta]; //profile y
old_diff=new_diff;
//} /* else { */
/* cout << "Difference not decreasing after "<< iint << " iterations (" << old_diff << " < " << new_diff << ")"<< endl; */
/* break; */
/* } */
iint++;
iint++;
}
delete [] newhhx;
delete [] newhhy;
delete [] newhhx;
delete [] newhhy;
if (hhx!=besthhx)
delete [] hhx;
if (hhy!=besthhy)
delete [] hhy;
hhx=besthhx;
hhy=besthhy;

View File

@ -14,10 +14,16 @@
class etaInterpolationBase : public slsInterpolation {
public:
etaInterpolationBase(int nx=400, int ny=400, int ns=25, int nb=-1, double emin=1, double emax=0) : slsInterpolation(nx,ny,ns), hhx(NULL), hhy(NULL), heta(NULL),nbeta(nb),etamin(emin), etamax(emax) {
if (nb<=0)
nbeta=nSubPixels*10;
etaInterpolationBase(int nx=400, int ny=400, int ns=25, int nb=-1, double emin=1, double emax=0) : slsInterpolation(nx,ny,ns), hhx(NULL), hhy(NULL), heta(NULL), nbeta(nb), etamin(emin), etamax(emax) {
// cout << "eb " << nb << " " << emin << " " << emax << endl;
// cout << nb << " " << etamin << " " << etamax << endl;
if (nbeta<=0) {
cout << "aaa:" <<endl;
nbeta=nSubPixels*10;
}
if (etamin>=etamax) {
cout << "aaa:" <<endl;
etamin=-1;
etamax=2;
}
@ -34,6 +40,7 @@ class etaInterpolationBase : public slsInterpolation {
#endif
//cout << nbeta << " " << etamin << " " << etamax << endl;
};
etaInterpolationBase(etaInterpolationBase *orig): slsInterpolation(orig){
@ -60,11 +67,11 @@ class etaInterpolationBase : public slsInterpolation {
};
virtual etaInterpolationBase* Clone() {
virtual etaInterpolationBase* Clone()=0;/*{
return new etaInterpolationBase(this);
};
*/
@ -179,8 +186,6 @@ class etaInterpolationBase : public slsInterpolation {
#endif
virtual void prepareInterpolation(int &ok){};
/* ////////////////////////////////////////////////////////////////////////////// */
@ -214,277 +219,35 @@ float *gethhx()
//////////////////////////////////////////////////////////////////////////////
//////////// /*It return position hit for the event in input */ //////////////
virtual void getInterpolatedPosition(int x, int y, double *data, double &int_x, double &int_y)
{
double sDum[2][2];
double tot, totquad;
double etax,etay;
int corner;
corner=calcQuad(data, tot, totquad, sDum);
if (nSubPixels>2)
calcEta(totquad, sDum, etax, etay);
getInterpolatedPosition(x,y,etax,etay,corner,int_x,int_y);
/* virtual void getInterpolatedPosition(int x, int y, int *data, double &int_x, double &int_y)=0; */
/* virtual void getInterpolatedPosition(int x, int y, double *data, double &int_x, double &int_y)=0; */
/* virtual void getInterpolatedPosition(int x, int y, double totquad,int quad,double *cl,double &int_x, double &int_y)=0; */
/* virtual void getInterpolatedPosition(int x, int y, double totquad,int quad,int *cl,double &int_x, double &int_y)=0; */
/* virtual void getInterpolatedPosition(int x, int y, double etax, double etay, int corner, double &int_x, double &int_y)=0; */
return;
};
/* virtual int addToFlatField(double totquad,int quad,int *cl,double &etax, double &etay)=0; */
/* virtual int addToFlatField(double totquad,int quad,double *cl,double &etax, double &etay)=0; */
/* virtual int addToFlatField(double *cluster, double &etax, double &etay)=0; */
/* virtual int addToFlatField(int *cluster, double &etax, double &etay)=0; */
virtual void getInterpolatedPosition(int x, int y, double totquad,int quad,double *cl,double &int_x, double &int_y) {
double cc[2][2];
double *cluster[3];
int xoff, yoff;
cluster[0]=cl;
cluster[1]=cl+3;
cluster[2]=cl+6;
switch (quad) {
case BOTTOM_LEFT:
xoff=0;
yoff=0;
break;
case BOTTOM_RIGHT:
xoff=1;
yoff=0;
break;
case TOP_LEFT:
xoff=0;
yoff=1;
break;
case TOP_RIGHT:
xoff=1;
yoff=1;
break;
default:
;
}
double etax, etay;
if (nSubPixels>2) {
cc[0][0]=cluster[yoff][xoff];
cc[1][0]=cluster[yoff+1][xoff];
cc[0][1]=cluster[yoff][xoff+1];
cc[1][1]=cluster[yoff+1][xoff+1];
calcEta(totquad,cc,etax,etay);
}
return getInterpolatedPosition(x,y,etax, etay,quad,int_x,int_y);
}
virtual void getInterpolatedPosition(int x, int y, double etax, double etay, int corner, double &int_x, double &int_y)
{
double xpos_eta=0,ypos_eta=0;
double dX,dY;
int ex,ey;
switch (corner)
{
case TOP_LEFT:
dX=-1.;
dY=0;
break;
case TOP_RIGHT:
;
dX=0;
dY=0;
break;
case BOTTOM_LEFT:
dX=-1.;
dY=-1.;
break;
case BOTTOM_RIGHT:
dX=0;
dY=-1.;
break;
default:
cout << "bad quadrant" << endl;
dX=0.;
dY=0.;
}
if (nSubPixels>2) {
#ifdef MYROOT1
xpos_eta=(hhx->GetBinContent(hhx->GetXaxis()->FindBin(etax),hhy->GetYaxis()->FindBin(etay)))/((double)nSubPixels);
ypos_eta=(hhy->GetBinContent(hhx->GetXaxis()->FindBin(etax),hhy->GetYaxis()->FindBin(etay)))/((double)nSubPixels);
#endif
#ifndef MYROOT1
ex=(etax-etamin)/etastep;
ey=(etay-etamin)/etastep;
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)]))+dX ;///((double)nSubPixels);
ypos_eta=(((double)hhy[(ey*nbeta+ex)]))+dY ;///((double)nSubPixels);
//else
//return 0;
#endif
} else {
xpos_eta=0.5*dX+0.25;
ypos_eta=0.5*dY+0.25;
}
int_x=((double)x) + xpos_eta+0.5;
int_y=((double)y) + ypos_eta+0.5;
}
/////////////////////////////////////////////////////////////////////////////////////////////////
virtual void getPositionETA3(int x, int y, double *data, double &int_x, double &int_y)
{
double sDum[2][2];
double tot, totquad;
double eta3x,eta3y;
double ex,ey;
calcQuad(data, tot, totquad, sDum);
calcEta3(data,eta3x, eta3y,tot);
double xpos_eta,ypos_eta;
#ifdef MYROOT1
xpos_eta=((hhx->GetBinContent(hhx->GetXaxis()->FindBin(eta3x),hhy->GetYaxis()->FindBin(eta3y))))/((double)nSubPixels);
ypos_eta=((hhy->GetBinContent(hhx->GetXaxis()->FindBin(eta3x),hhy->GetYaxis()->FindBin(eta3y))))/((double)nSubPixels);
#endif
#ifndef MYROOT1
ex=(eta3x-etamin)/etastep;
ey=(eta3y-etamin)/etastep;
if (ex<0) ex=0;
if (ex>=nbeta) ex=nbeta-1;
if (ey<0) ey=0;
if (ey>=nbeta) ey=nbeta-1;
xpos_eta=(((double)hhx[(int)(ey*nbeta+ex)]))/((double)nSubPixels);
ypos_eta=(((double)hhy[(int)(ey*nbeta+ex)]))/((double)nSubPixels);
#endif
int_x=((double)x) + xpos_eta;
int_y=((double)y) + ypos_eta;
return;
};
virtual int addToFlatField(double totquad,int quad,double *cl,double &etax, double &etay) {
double cc[2][2];
double *cluster[3];
int xoff, yoff;
cluster[0]=cl;
cluster[1]=cl+3;
cluster[2]=cl+6;
switch (quad) {
case BOTTOM_LEFT:
xoff=0;
yoff=0;
break;
case BOTTOM_RIGHT:
xoff=1;
yoff=0;
break;
case TOP_LEFT:
xoff=0;
yoff=1;
break;
case TOP_RIGHT:
xoff=1;
yoff=1;
break;
default:
;
}
cc[0][0]=cluster[yoff][xoff];
cc[1][0]=cluster[yoff+1][xoff];
cc[0][1]=cluster[yoff][xoff+1];
cc[1][1]=cluster[yoff+1][xoff+1];
/* cout << cl[0] << " " << cl[1] << " " << cl[2] << endl; */
/* cout << cl[3] << " " << cl[4] << " " << cl[5] << endl; */
/* cout << cl[6] << " " << cl[7] << " " << cl[8] << endl; */
/* cout <<"******"<<totquad << " " << quad << endl; */
/* cout << cc[0][0]<< " " << cc[0][1] << endl; */
/* cout << cc[1][0]<< " " << cc[1][1] << endl; */
//calcMyEta(totquad,quad,cl,etax, etay);
calcEta(totquad, cc,etax, etay);
// cout <<"******"<< etax << " " << etay << endl;
return addToFlatField(etax,etay);
}
//////////////////////////////////////////////////////////////////////////////////////
virtual int addToFlatField(double *cluster, double &etax, double &etay){
double sDum[2][2];
double tot, totquad;
int corner;
corner=calcQuad(cluster, tot, totquad, sDum);
double xpos_eta,ypos_eta;
double dX,dY;
calcEta(totquad, sDum, etax, etay);
return addToFlatField(etax,etay);
};
virtual int addToFlatField(double etax, double etay){
int ex,ey;
#ifdef MYROOT1
heta->Fill(etax,etay);
#endif
#ifndef MYROOT1
int ex,ey;
ex=(etax-etamin)/etastep;
ey=(etay-etamin)/etastep;
// cout << etax << " " << ex << " " << etay << " " << ey << " " << ey*nbeta+ex << endl;
if (ey<nbeta && ex<nbeta && ex>=0 && ey>=0)
heta[ey*nbeta+ex]++;
// cout << "*"<< etax << " " << etay << endl;
/* cout << etax << " " << etay << " " << ex << " " << ey << " " << ey*nbeta+ex << endl; */
/* cout <<"********"<< endl << endl ; */
#endif
return 0;
};
};
// virtual void prepareInterpolation(int &ok)=0;
protected:
#ifdef MYROOT1

View File

@ -2,24 +2,27 @@
#define ETA_INTERPOLATION_POSXY_H
#include "tiffIO.h"
//#include "tiffIO.h"
#include "etaInterpolationBase.h"
#include "eta2InterpolationBase.h"
#include "eta3InterpolationBase.h"
class etaInterpolationPosXY : public etaInterpolationBase{
class etaInterpolationPosXY : public virtual etaInterpolationBase{
public:
etaInterpolationPosXY(int nx=400, int ny=400, int ns=25, int nb=-1, double emin=1, double emax=0) : etaInterpolationBase(nx,ny,ns, nb, emin,emax){};
etaInterpolationPosXY(int nx=400, int ny=400, int ns=25, int nb=-1, double emin=1, double emax=0) : etaInterpolationBase(nx,ny,ns, nb, emin,emax){
// cout << "epxy " << nb << " " << emin << " " << emax << endl; cout << nbeta << " " << etamin << " " << etamax << endl;
};
etaInterpolationPosXY(etaInterpolationPosXY *orig): etaInterpolationBase(orig){};
etaInterpolationPosXY(etaInterpolationPosXY *orig): etaInterpolationBase(orig) {};
virtual etaInterpolationPosXY* Clone() {
virtual etaInterpolationPosXY* Clone()=0;/** {
return new etaInterpolationPosXY(this);
};
};*/
virtual void prepareInterpolation(int &ok)
{
cout <<"?"<< endl;
ok=1;
#ifdef MYROOT1
if (hhx) delete hhx;
@ -114,10 +117,6 @@ class etaInterpolationPosXY : public etaInterpolationBase{
// cout << "y " << nbeta << " " << (ii+1)*tot_eta_x*bsize << " " << ii << endl;
}
#ifdef SAVE_ALL
@ -150,4 +149,25 @@ class etaInterpolationPosXY : public etaInterpolationBase{
};
class eta2InterpolationPosXY : public virtual eta2InterpolationBase, public virtual etaInterpolationPosXY {
public:
eta2InterpolationPosXY(int nx=400, int ny=400, int ns=25, int nb=-1, double emin=1, double emax=0) : etaInterpolationBase(nx,ny,ns, nb, emin,emax),eta2InterpolationBase(nx,ny,ns, nb, emin,emax),etaInterpolationPosXY(nx,ny,ns, nb, emin,emax){
// cout << "e2pxy " << nb << " " << emin << " " << emax << endl;
};
eta2InterpolationPosXY(eta2InterpolationPosXY *orig): etaInterpolationBase(orig), etaInterpolationPosXY(orig) {};
virtual eta2InterpolationPosXY* Clone() { return new eta2InterpolationPosXY(this);};
};
class eta3InterpolationPosXY : public virtual eta3InterpolationBase, public virtual etaInterpolationPosXY {
public:
eta3InterpolationPosXY(int nx=400, int ny=400, int ns=25, int nb=-1, double emin=1, double emax=0) : etaInterpolationBase(nx,ny,ns, nb, emin,emax),eta3InterpolationBase(nx,ny,ns, nb, emin,emax), etaInterpolationPosXY(nx,ny,ns, nb, emin,emax){
cout << "e3pxy " << nbeta << " " << etamin << " " << etamax << " " << nSubPixels<< endl;
};
eta3InterpolationPosXY(eta3InterpolationPosXY *orig): etaInterpolationBase(orig), etaInterpolationPosXY(orig) {};
virtual eta3InterpolationPosXY* Clone() { return new eta3InterpolationPosXY(this);};
};
#endif

View File

@ -0,0 +1,417 @@
#ifndef ETA_INTERPOLATION_RANDOMBINS_H
#define ETA_INTERPOLATION_RANDOMBINS_H
#include "tiffIO.h"
//#include "etaInterpolationBase.h"
#include "etaInterpolationPosXY.h"
#include <cstdlib>
#include <algorithm>
//#include <math>
#include <cmath> // std::abs
#define PI 3.14159265
#define TWOPI 2.*PI
using namespace std;
class etaInterpolationRandomBins : public etaInterpolationPosXY {
private:
double calcDiff(double avg, float *hx, float *hy) {
double p_tot=0;
double diff=0;
double bsize=1./nSubPixels;
for (int ipx=0; ipx<nSubPixels; ipx++) {
for (int ipy=0; ipy<nSubPixels; ipy++) {
p_tot=0;
for (int ibx=0; ibx<nbeta; ibx++) {
for (int iby=0; iby<nbeta; iby++) {
if ( hx[ibx+iby*nbeta]>=((ipx)*bsize) && hx[ibx+iby*nbeta]<((ipx+1)*bsize) && hy[ibx+iby*nbeta]>=((ipy)*bsize) && hy[ibx+iby*nbeta]<((ipy+1)*bsize)) {
p_tot+=heta[ibx+iby*nbeta];
}
}
}
// cout << p_tot << " \t ";
diff+=(p_tot-avg)*(p_tot-avg);
}
// cout << "\n";
}
return diff;
}
double iterate(float *newhhx, float *newhhy, double avg) {
double bsize=1./nSubPixels;
double hy[nbeta]; //profile y
double hx[nbeta]; //profile x
double hix[nbeta]; //integral of projection x
double hiy[nbeta]; //integral of projection y
double tot_eta_x=0;
double tot_eta_y=0;
int p0;
int vx[(nSubPixels+1)*(nSubPixels+1)], vy[(nSubPixels+1)*(nSubPixels+1)];
int arrx[nSubPixels+1], arry[nSubPixels+1];
int bad=1;
int isby, isbx;
int ii=0;
// using default comparison (operator <):
// std::sort (myvector.begin(), myvector.begin()+4); //(12 32 45 71)26 80 53 33
for (isby=0; isby<(nSubPixels+1)/2+1; isby++) {
for (isbx=0; isbx<(nSubPixels+1)/2+1; isbx++) {
p0=isby*(nSubPixels+1)+isbx;
// for (int iv=0; iv<(nSubPixels+1)*(nSubPixels+1); iv++) {
if (isbx==0) {
vy[p0]=isby*nbeta/nSubPixels;
vx[p0]=0;
} else if ( isby==0 ) {
vy[p0]=0;
vx[p0]=isbx*nbeta/nSubPixels;
}
else {
vy[p0]=rand()%(nbeta/2);
vx[p0]=rand()%(nbeta/2);
if (nSubPixels%2==0 && isbx==nSubPixels/2)
vx[p0]=nbeta/2;
if (nSubPixels%2==0 && isby==nSubPixels/2 )
vy[p0]=nbeta/2;
}
// cout << "(" << vx[p0] << " , " << vy[p0] << " ) \t" ;
// }
}
//cout << endl;
}
// cout << "rand" << endl;
while (bad) {
for (isby=0; isby<(nSubPixels+1)/2+1; isby++) {
for (isbx=0; isbx<(nSubPixels+1)/2+1; isbx++) {
arrx[isbx]=vx[isby*(nSubPixels+1)+isbx];
arry[isbx]=vy[isbx*(nSubPixels+1)+isby];
//cout << isbx << " " << arrx[isbx] << " " << isby << " " << arry[isbx] << endl;
}
sort(arrx,arrx+(nSubPixels+1)/2+1);
sort(arry,arry+(nSubPixels+1)/2+1);
// cout << "*****"<< endl;
// cout << endl;
for (int isbx=0; isbx<(nSubPixels+1)/2+1; isbx++) {
vx[isby*(nSubPixels+1)+isbx]=arrx[isbx];
vy[isbx*(nSubPixels+1)+isby]=arry[isbx];
vx[(nSubPixels-isby)*(nSubPixels+1)+(nSubPixels-isbx)]=nbeta-arrx[isbx];
vy[(nSubPixels-isbx)*(nSubPixels+1)+(nSubPixels-isby)]=nbeta-arry[isbx];
vx[isby*(nSubPixels+1)+(nSubPixels-isbx)]=nbeta-arrx[isbx];
vy[isbx*(nSubPixels+1)+(nSubPixels-isby)]=arry[isbx];
vx[(nSubPixels-isby)*(nSubPixels+1)+(isbx)]=arrx[isbx];
vy[(nSubPixels-isbx)*(nSubPixels+1)+(isby)]=nbeta-arry[isbx];
}
}
/* for (isby=0; isby<nSubPixels+1; isby++) { */
/* for (isbx=0; isbx<nSubPixels+1; isbx++) { */
/* cout << "("<< vx[isby*(nSubPixels+1)+isbx] << " " << vy[isby*(nSubPixels+1)+isbx] << ")\t";//<< endl; */
/* } */
/* cout << endl; */
/* } */
bad=0;
for (isby=1; isby<(nSubPixels+1)/2+1; isby++) {
for (isbx=1; isbx<(nSubPixels+1)/2+1; isbx++) {
if (heta[vx[isby*(nSubPixels+1)+isbx]+vy[isby*(nSubPixels+1)+isbx]*nbeta]<avg*(nSubPixels*nSubPixels)/(nbeta*nbeta)) {
// cout << ii << " " << isbx << " " << isby << " " << vx[isby*(nSubPixels+1)+isbx] << " " << vy[isby*(nSubPixels+1)+isbx] << " " << heta[vx[isby*(nSubPixels+1)+isbx]+vy[isby*(nSubPixels+1)+isbx]*nbeta] << endl;
if (nSubPixels%2==0 && isbx==nSubPixels/2)
;
else
vx[isby*(nSubPixels+1)+isbx]=rand()%(nbeta/2);
if (nSubPixels%2==0 && isbx==nSubPixels/2)
;
else
vy[isby*(nSubPixels+1)+isbx]=rand()%(nbeta/2);
if (bad==0)
ii++;
bad=1;
// break;
}
}
//if (bad) break;
}
// cout << "sort" << endl;
}
cout << ii << " sub iteractions " << avg*(nSubPixels*nSubPixels)/(nbeta*nbeta) << endl;
double m,q;
int in_quad;
int p[4];
int p1x,p2x, p1y, p2y;
// cout << nbeta << endl;
double angle;
double dtheta,theta1,theta2;
for (int ibx=0; ibx<nbeta; ibx++) {
for (int iby=0; iby<nbeta; iby++) {
in_quad=0;
/* if (ibx==0) */
/* isbx=0; */
/* else */
/* isbx= (newhhx[ibx-1+iby*nbeta])/bsize-1; */
/* if (isbx<0) isbx=0; */
/* if (isbx>nSubPixels-1) isbx=nSubPixels-1; */
/* if (iby==0) */
/* isby=0; */
/* else */
/* isby= (newhhx[ibx+(iby-1)*nbeta])/bsize-1; */
/* if (isby<0) isbx=0; */
/* if (isby>nSubPixels-1) isby=nSubPixels-1; */
/* // cout << isbx << " " << isby << endl; */
for (isby=0; isby<nSubPixels; isby++) {
for (isbx=0; isbx<nSubPixels; isbx++) {
// cout << ibx << " " << iby << " " << isbx << " " << isby << endl;
p[0]=isby*(nSubPixels+1)+isbx;
p[1]=isby*(nSubPixels+1)+isbx+1;
p[2]=(isby+1)*(nSubPixels+1)+isbx+1;
p[3]=(isby+1)*(nSubPixels+1)+isbx;
angle=0;
for (int i=0;i<4;i++) {
p1x = vx[p[i]] - ibx;
p1y = vy[p[i]] - iby;
p2x = vx[p[(i+1)%4]] - ibx;
p2y = vy[p[(i+1)%4]] - iby;
theta1 = atan2(p1y,p1x);
theta2 = atan2(p2y,p2x);
dtheta = theta2 - theta1;
while (dtheta > PI)
dtheta -= TWOPI;
while (dtheta < -PI)
dtheta += TWOPI;
angle += dtheta;
}
if (abs((double)angle) < PI)
in_quad=0;
else
in_quad=1;
if (in_quad) {
newhhx[ibx+iby*nbeta]=bsize*((double)isbx);
newhhy[ibx+iby*nbeta]=bsize*((double)isby);
break;
}
}
if (in_quad) break;
}
}
}
// cout << "hist" << endl;
return calcDiff(avg, newhhx, newhhy);
}
public:
etaInterpolationRandomBins(int nx=400, int ny=400, int ns=25, int nb=-1, double emin=1, double emax=0) : etaInterpolationPosXY(nx,ny,ns, nb, emin,emax){};
etaInterpolationRandomBins(etaInterpolationRandomBins *orig): etaInterpolationPosXY(orig){};
virtual etaInterpolationRandomBins* Clone() {
return new etaInterpolationRandomBins(this);
};
virtual void prepareInterpolation(int &ok)
{
ok=1;
cout << "Adaptive bins" << endl;
///*Eta Distribution Rebinning*///
double bsize=1./nSubPixels; //precision
// cout<<"nPixelsX = "<<nPixelsX<<" nPixelsY = "<<nPixelsY<<" nSubPixels = "<<nSubPixels<<endl;
double tot_eta=0;
for (int ip=0; ip<nbeta*nbeta; ip++)
tot_eta+=heta[ip];
if (tot_eta<=0) {ok=0; return;};
int ii=0;
int nint=1000;
double thr=1./((double)nSubPixels);
double avg=tot_eta/((double)(nSubPixels*nSubPixels));
cout << "total eta entries is :"<< tot_eta << " avg: "<< avg << endl;
cout << "Start " << endl;
double old_diff=-1, new_diff=-1;
// cout << " diff= " << new_diff << endl;
etaInterpolationPosXY::prepareInterpolation(ok);
old_diff=calcDiff(avg, hhx, hhy);
cout << " diff= " << old_diff << endl;
int iint=0;
float *newhhx=new float[nbeta*nbeta]; //profile x
float *newhhy=new float[nbeta*nbeta]; //profile y
int igood=0, ibad=0;
#ifdef SAVE_ALL
int etabins=nbeta;
float *etah=new float[nbeta*nbeta];
char tit[1000];
#endif
while (iint<nint) {
cout << "Iteration " << iint << endl;
new_diff=iterate(newhhx,newhhy, avg);
//new_diff=calcDiff(avg, newhhx, newhhy);
cout << " diff= " << new_diff << " ( " << old_diff<< " ) " << endl;
/* #ifdef SAVE_ALL */
/* for (int ii=0; ii<etabins*etabins; ii++) { */
/* etah[ii]=newhhx[ii]; */
/* if (etah[ii]>1 || etah[ii]<0 ) cout << "***"<< ii << etah[ii] << endl; */
/* } */
/* sprintf(tit,"/scratch/randeta_hhx_%d.tiff",iint); */
/* WriteToTiff(etah, tit, etabins, etabins); */
/* for (int ii=0; ii<etabins*etabins; ii++) { */
/* etah[ii]=newhhy[ii]; */
/* if (etah[ii]>1 || etah[ii]<0 ) cout << "***"<< ii << etah[ii] << endl; */
/* } */
/* sprintf(tit,"/scratch/randeta_hhy_%d.tiff",iint); */
/* WriteToTiff(etah, tit, etabins, etabins); */
/* #endif */
if (new_diff<old_diff) {
cout << "******************** GOOD! ***********************"<< endl;
delete [] hhx;
delete [] hhy;
igood++;
hhx=newhhx;
hhy=newhhy;
newhhx=new float[nbeta*nbeta]; //profile x */
newhhy=new float[nbeta*nbeta]; //profile y */
old_diff=new_diff;
} else
ibad++;
iint++;
}
delete [] newhhx;
delete [] newhhy;
cout << "performed " << iint << " iterations of which " << igood << " positive " << endl;
/* #ifdef SAVE_ALL */
/* for (int ii=0; ii<etabins*etabins; ii++) { */
/* etah[ii]=hhx[ii]; */
/* } */
/* sprintf(tit,"/scratch/eta_hhx_%d.tiff",id); */
/* WriteToTiff(etah, tit, etabins, etabins); */
/* for (int ii=0; ii<etabins*etabins; ii++) { */
/* etah[ii]=hhy[ii]; */
/* } */
/* sprintf(tit,"/scratch/eta_hhy_%d.tiff",id); */
/* WriteToTiff(etah, tit, etabins, etabins); */
/* for (int ii=0; ii<etabins*etabins; ii++) { */
/* etah[ii]=heta[ii]; */
/* } */
/* sprintf(tit,"/scratch/eta_%d.tiff",id); */
/* WriteToTiff(etah, tit, etabins, etabins); */
/* delete [] etah; */
/* #endif */
return ;
}
};
#endif

View File

@ -38,13 +38,31 @@ class noInterpolation : public slsInterpolation{
return ;
};
virtual void getInterpolatedPosition(int x, int y, int *data, double &int_x, double &int_y)
{
return getInterpolatedPosition(x, y, (double*)NULL, int_x, int_y);
}
virtual void getInterpolatedPosition(int x, int y, double etax, double etay, int corner, double &int_x, double &int_y)
{
getInterpolatedPosition(x, y, NULL, int_x, int_y);
getInterpolatedPosition(x, y, (double*)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);
getInterpolatedPosition(x, y, (double*)NULL, etax, etay);
};
virtual void getInterpolatedPosition(int x, int y, double totquad,int quad,int *cl,double &etax, double &etay){
getInterpolatedPosition(x, y, (double*)NULL, etax, etay);
};
//////////////////////////////////////////////////////////////////////////////////////
virtual void getPositionETA3(int x, int y, double *data, double &int_x, double &int_y)
@ -55,13 +73,27 @@ class noInterpolation : public slsInterpolation{
return ;
};
virtual void getPositionETA3(int x, int y, int *data, double &int_x, double &int_y)
{
return getPositionETA3(x, y, (double*)NULL, int_x, int_y);
};
//////////////////////////////////////////////////////////////////////////////////////
virtual int addToFlatField(double *cluster, double &etax, double &etay){return 0;};
virtual int addToFlatField(int *cluster, double &etax, double &etay){return 0;};
virtual int addToFlatField(double etax, double etay){return 0;};
virtual int addToFlatField(double totquad,int quad,double *cl,double &etax, double &etay){return 0;};
virtual int addToFlatField(double totquad,int quad,int *cl,double &etax, double &etay){return 0;};
protected:
;
// TRandom *eventGenerator;

View File

@ -8,7 +8,10 @@
#endif
#include <cstdlib>
#ifndef MY_TIFF_IO_H
#include "tiffIO.h"
#endif
#ifndef DEF_QUAD
#define DEF_QUAD
enum quadrant {
@ -88,7 +91,7 @@ hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny);
void *writeInterpolatedImage(const char * imgname) {
cout << "!" <<endl;
//cout << "!" <<endl;
float *gm=NULL;
gm=new float[ nSubPixels* nSubPixels* nPixelsX*nPixelsY];
if (gm) {
@ -107,8 +110,10 @@ hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny);
//return position inside the pixel for the given photon
virtual void getInterpolatedPosition(int x, int y, double *data, double &int_x, double &int_y)=0;
virtual void getInterpolatedPosition(int x, int y, int *data, double &int_x, double &int_y)=0;
//return position inside the pixel for the given photon
virtual void getInterpolatedPosition(int x, int y, double etax, double etay, int quad, double &int_x, double &int_y)=0;
virtual void getInterpolatedPosition(int x, int y, double totquad,int quad,int *cluster,double &etax, double &etay)=0;
virtual void getInterpolatedPosition(int x, int y, double totquad,int quad,double *cluster,double &etax, double &etay)=0;
@ -138,7 +143,9 @@ hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny);
#endif
#ifndef MYROOT1
virtual int *addToImage(double int_x, double int_y){ int iy=nSubPixels*int_y; int ix=nSubPixels*int_x;
virtual int *addToImage(double int_x, double int_y){
int iy=((double)nSubPixels)*int_y;
int ix=((double)nSubPixels)*int_x;
// cout << int_x << " " << int_y << " " << " " << ix << " " << iy << " " << ix+iy*nPixelsX*nSubPixels << endl;
if (ix>=0 && ix<(nPixelsX*nSubPixels) && iy<(nSubPixels*nPixelsY) && iy>=0 )(*(hint+ix+iy*nPixelsX*nSubPixels))+=1;
return hint;
@ -147,9 +154,10 @@ hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny);
virtual int addToFlatField(double *cluster, double &etax, double &etay)=0;
virtual int addToFlatField(double etax, double etay)=0;
virtual int addToFlatField(int *cluster, double &etax, double &etay)=0;
virtual int addToFlatField(double totquad,int quad,int *cl,double &etax, double &etay)=0;
virtual int addToFlatField(double totquad,int quad,double *cluster,double &etax, double &etay)=0;
virtual int addToFlatField(double etax, double etay)=0;
#ifdef MYROOT1
virtual TH2D *getFlatField(){return NULL;};
@ -167,6 +175,14 @@ hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny);
//virtual void Streamer(TBuffer &b);
static int calcQuad(int *cl, double &sum, double &totquad, double sDum[2][2]){
double cli[3*3];//=new int[3*3];
for (int i=0; i<9; i++)
cli[i]=cl[i];
return calcQuad(cli, sum, totquad, sDum);
}
static int calcQuad(double *cl, double &sum, double &totquad, double sDum[2][2]){
@ -175,48 +191,64 @@ hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny);
cluster[0]=cl;
cluster[1]=cl+3;
cluster[2]=cl+6;
sum=0;
double sumBL=0;
double sumTL=0;
double sumBR=0;
double sumTR=0;
int xoff=0, yoff=0;
for (int ix=0; ix<3; ix++) {
for (int iy=0; iy<3; iy++) {
sum+=cluster[iy][ix];
if (ix<=1 && iy<=1) sumBL+=cluster[iy][ix];
if (ix<=1 && iy>=1) sumTL+=cluster[iy][ix];
if (ix>=1 && iy<=1) sumBR+=cluster[iy][ix];
if (ix>=1 && iy>=1) sumTR+=cluster[iy][ix];
}
}
sum = cluster[0][0] + cluster[1][0] + cluster[2][0] + cluster[0][1] + cluster[1][1] + cluster[2][1] + cluster[0][2] + cluster[1][2] + cluster[2][2];
double sumBL = cluster[0][0] + cluster[1][0] + cluster[0][1] + cluster[1][1]; //2 ->BL
double sumTL = cluster[1][0] + cluster[2][0] + cluster[2][1] + cluster[1][1]; //0 ->TL
double sumBR = cluster[0][1] + cluster[0][2] + cluster[1][2] + cluster[1][1]; //3 ->BR
double sumTR = cluster[1][2] + cluster[2][1] + cluster[2][2] + cluster[1][1]; //1 ->TR
double sumMax = 0;
double t, r;
// if(sumTL >= sumMax){
sDum[0][0] = cluster[0][0]; sDum[1][0] = cluster[1][0];
sDum[0][1] = cluster[0][1]; sDum[1][1] = cluster[1][1];
/* sDum[0][0] = cluster[0][0]; sDum[1][0] = cluster[1][0]; */
/* sDum[0][1] = cluster[0][1]; sDum[1][1] = cluster[1][1]; */
corner = BOTTOM_LEFT;
sumMax=sumBL;
// }
if(sumTL >= sumMax){
sDum[0][0] = cluster[1][0]; sDum[1][0] = cluster[2][0];
sDum[0][1] = cluster[1][1]; sDum[1][1] = cluster[2][1];
totquad=sumBL;
if(sumTL >= totquad){
/* sDum[0][0] = cluster[1][0]; sDum[1][0] = cluster[2][0]; */
/* sDum[0][1] = cluster[1][1]; sDum[1][1] = cluster[2][1]; */
corner = TOP_LEFT;
sumMax=sumTL;
totquad=sumTL;
xoff=0;
yoff=1;
}
if(sumBR >= sumMax){
sDum[0][0] = cluster[0][1]; sDum[1][0] = cluster[1][1];
sDum[0][1] = cluster[0][2]; sDum[1][1] = cluster[1][2];
if(sumBR >= totquad){
/* sDum[0][0] = cluster[0][1]; sDum[1][0] = cluster[1][1]; */
/* sDum[0][1] = cluster[0][2]; sDum[1][1] = cluster[1][2]; */
xoff=1;
yoff=0;
corner = BOTTOM_RIGHT;
sumMax=sumBR;
totquad=sumBR;
}
if(sumTR >= sumMax){
sDum[0][0] = cluster[1][1]; sDum[1][0] = cluster[2][1];
sDum[0][1] = cluster[1][2]; sDum[1][1] = cluster[2][2];
if(sumTR >= totquad){
xoff=1;
yoff=1;
/* sDum[0][0] = cluster[1][1]; sDum[1][0] = cluster[2][1]; */
/* sDum[0][1] = cluster[1][2]; sDum[1][1] = cluster[2][2]; */
corner = TOP_RIGHT;
sumMax=sumTR;
totquad=sumTR;
}
totquad=sumMax;
for (int ix=0; ix<2; ix++) {
for (int iy=0; iy<2; iy++) {
sDum[iy][ix] = cluster[iy+yoff][ix+xoff];
}
}
return corner;
@ -235,7 +267,6 @@ hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny);
}
static int calcEta(double *cl, double &etax, double &etay, double &sum, double &totquad, double sDum[2][2]) {
int corner = calcQuad(cl,sum,totquad,sDum);
calcEta(totquad, sDum, etax, etay);
@ -244,6 +275,14 @@ hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny);
}
static int calcEta(int *cl, double &etax, double &etay, double &sum, double &totquad, double sDum[2][2]) {
int corner = calcQuad(cl,sum,totquad,sDum);
calcEta(totquad, sDum, etax, etay);
return corner;
}
static int calcEtaL(double totquad, int corner, double sDum[2][2], double &etax, double &etay){
double t,r, toth, totv;
if (totquad>0) {
@ -294,6 +333,13 @@ hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny);
return corner;
}
static int calcEtaL(int *cl, double &etax, double &etay, double &sum, double &totquad, double sDum[2][2]) {
int corner = calcQuad(cl,sum,totquad,sDum);
calcEtaL(totquad, corner, sDum, etax, etay);
return corner;
}
static int calcEtaC3(double *cl, double &etax, double &etay, double &sum, double &totquad, double sDum[2][2]){
@ -306,23 +352,95 @@ hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny);
static int calcEtaC3(int *cl, double &etax, double &etay, double &sum, double &totquad, double sDum[2][2]){
int corner = calcQuad(cl,sum,totquad,sDum);
calcEta(sum, sDum, etax, etay);
return corner;
}
static int calcEta3(double *cl, double &etax, double &etay, double &sum) {
double l,r,t,b;
sum=cl[0]+cl[1]+cl[2]+cl[3]+cl[4]+cl[5]+cl[6]+cl[7]+cl[8];
double l=0,r=0,t=0,b=0, val;
sum=0;
// int quad;
for (int ix=0; ix<3; ix++) {
for (int iy=0; iy<3; iy++) {
val=cl[iy+3*ix];
sum+=val;
if (iy==0) l+=val;
if (iy==2) r+=val;
if (ix==0) b+=val;
if (ix==2) t+=val;
}
}
if (sum>0) {
l=cl[0]+cl[3]+cl[6];
r=cl[2]+cl[5]+cl[8];
b=cl[0]+cl[1]+cl[2];
t=cl[6]+cl[7]+cl[8];
etax=(-l+r)/sum;
etay=(-b+t)/sum;
}
/* if (etax<-1 || etax>1 || etay<-1 || etay>1) { */
/* cout << "**********" << etax << " " << etay << endl; */
/* for (int ix=0; ix<3; ix++) { */
/* for (int iy=0; iy<3; iy++) { */
/* cout << cl[iy+3*ix] << "\t" ; */
/* } */
/* cout << endl; */
/* } */
/* cout << sum << " " << l << " " << r << " " << t << " " << b << endl; */
/* } */
if (etax>=0 && etay>=0)
return TOP_RIGHT;
if (etax<0 && etay>=0)
return TOP_LEFT;
if (etax<0 && etay<0)
return BOTTOM_LEFT;
return BOTTOM_RIGHT;
}
static int calcEta3(int *cl, double &etax, double &etay, double &sum) {
double cli[9];
for (int ix=0; ix<9; ix++) cli[ix]=cl[ix];
return calcEta3(cli, etax, etay, sum);
}
static int calcMyEta(double totquad, int quad, double *cl, double &etax, double &etay) {
double l,r,t,b, sum;
int yoff;
switch (quad) {
case BOTTOM_LEFT:
case BOTTOM_RIGHT:
yoff=0;
break;
case TOP_LEFT:
case TOP_RIGHT:
yoff=1;
break;
default:
;
}
l=cl[0+yoff*3]+cl[0+yoff*3+3];
r=cl[2+yoff*3]+cl[2+yoff*3+3];
b=cl[0+yoff*3]+cl[1+yoff*3]*cl[2+yoff*3];
t=cl[0+yoff*3+3]+cl[1+yoff*3+3]*cl[0+yoff*3+3];
sum=t+b;
if (sum>0) {
etax=(-l+r)/sum;
etay=(+t)/sum;
}
return -1;
}
static int calcMyEta(double totquad, int quad, double *cl, double &etax, double &etay) {
static int calcMyEta(double totquad, int quad, int *cl, double &etax, double &etay) {
double l,r,t,b, sum;
int yoff;
switch (quad) {
@ -367,6 +485,21 @@ hint=new TH2F("hint","hint",ns*nx, 0, nx, ns*ny, 0, ny);
}
static int calcEta3X(int *cl, double &etax, double &etay, double &sum) {
double l,r,t,b;
sum=cl[0]+cl[1]+cl[2]+cl[3]+cl[4]+cl[5]+cl[6]+cl[7]+cl[8];
if (sum>0) {
l=cl[3];
r=cl[5];
b=cl[1];
t=cl[7];
etax=(-l+r)/sum;
etay=(-b+t)/sum;
}
return -1;
}

View File

@ -11,23 +11,29 @@ MAIN=moench03ClusterFinder.cpp
#-lhdf5
#DESTDIR?=../bin
all: moenchClusterFinder moenchMakeEta moenchInterpolation moenchNoInterpolation
all: moenchClusterFinder moenchMakeEta moenchInterpolation moenchNoInterpolation moenchPhotonCounter moenchAnalog
moenchClusterFinder: $(MAIN) $(INCS) clean
g++ -o moenchClusterFinder $(MAIN) $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF) -DSAVE_ALL -DOLDDATA
g++ -o moenchClusterFinder $(MAIN) $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF) -DSAVE_ALL -DOLDDATA #-DNEWRECEIVER
moenchMakeEta: moench03MakeEta.cpp $(INCS) clean
g++ -o moenchMakeEta moench03MakeEta.cpp $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF) -DSAVE_ALL
moenchInterpolation: moench03Interpolation.cpp $(INCS) clean
g++ -o moenchInterpolation moench03Interpolation.cpp $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF) -DSAVE_ALL
g++ -o moenchInterpolation moench03Interpolation.cpp $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF) -DSAVE_ALL
moenchNoInterpolation: moench03NoInterpolation.cpp $(INCS) clean
g++ -o moenchNoInterpolation moench03NoInterpolation.cpp $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF) -DSAVE_ALL
moenchPhotonCounter: moenchPhotonCounter.cpp $(INCS) clean
g++ -o moenchPhotonCounter moenchPhotonCounter.cpp $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF) -DSAVE_ALL -DNEWRECEIVER
moenchAnalog: moenchPhotonCounter.cpp $(INCS) clean
g++ -o moenchAnalog moenchPhotonCounter.cpp $(LDFLAG) $(INCDIR) $(LIBHDF5) $(LIBRARYCBF) -DSAVE_ALL -DNEWRECEIVER -DANALOG
clean:
rm -f moenchClusterFinder moenchMakeEta moenchInterpolation moenchNoInterpolation
rm -f moenchClusterFinder moenchMakeEta moenchInterpolation moenchNoInterpolation moenchPhotonCounter

View File

@ -119,53 +119,11 @@ int main(int argc, char *argv[]) {
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;
@ -176,9 +134,7 @@ int main(int argc, char *argv[]) {
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);
@ -187,46 +143,7 @@ int main(int argc, char *argv[]) {
cout << "mt " << endl;
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);
@ -259,6 +176,8 @@ int main(int argc, char *argv[]) {
mt->nextThread();
// // // cout << " " << (void*)buff;
mt->popFree(buff);
ifr++;
if (ifr%10000==0) cout << ifr << " " << ff << endl;
ff=-1;
}
cout << "--" << endl;

View File

@ -2,53 +2,101 @@
#include "ansi.h"
#include <iostream>
#include "moench03T1ZmqData.h"
#include "single_photon_hit.h"
//#include "moench03T1ZmqData.h"
//#define DOUBLE_SPH
//#define MANYFILES
// #include "etaInterpolationPosXY.h"
#include "etaInterpolationAdaptiveBins.h"
#ifdef DOUBLE_SPH
#include "single_photon_hit_double.h"
#endif
#ifndef DOUBLE_SPH
#include "single_photon_hit.h"
#endif
#include "etaInterpolationPosXY.h"
#include "noInterpolation.h"
//#include "etaInterpolationAdaptiveBins.h"
//#include "etaInterpolationRandomBins.h"
using namespace std;
#define NC 400
#define NR 400
#define XTALK
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;
if (argc<9) {
cout << "Wrong usage! Should be: "<< argv[0] << " infile " << " etafile outfile runmin runmax ns cmin cmax" << endl;
return 1;
}
char infname[10000];
char fname[10000];
char outfname[10000];
int runmin=atoi(argv[4]);
int runmax=atoi(argv[5]);
int nsubpix=atoi(argv[6]);
float cmin=atof(argv[7]);
float cmax=atof(argv[8]);
int etabins=1000;//nsubpix*2*100;
double etamin=-1, etamax=2;
double eta3min=-2, eta3max=2;
int quad;
double sum, totquad;
double sDum[2][2];
double etax, etay, int_x, int_y;
double eta3x, eta3y, int3_x, int3_y, noint_x, noint_y;
int ok;
int f0=-1;
int ix, iy, isx, isy;
int nframes=0, lastframe=-1;
double d_x, d_y, res=5, xx, yy;
#ifdef MANYFILES
int ff=1000;
#endif
int nph=0, badph=0, totph=0;
FILE *f=NULL;
#ifdef DOUBLE_SPH
single_photon_hit_double cl(3,3);
#endif
#ifndef DOUBLE_SPH
single_photon_hit cl(3,3);
// etaInterpolationPosXY *interp=new etaInterpolationPosXY(NC, NR, nsubpix, etabins, etamin, etamax);
etaInterpolationAdaptiveBins *interp=new etaInterpolationAdaptiveBins (NC, NR, nsubpix, etabins, etamin, etamax);
#endif
#ifdef XTALK
int old_val[3][3];
int new_val[3][3];
double xcorr=0.04;
// int ix=0;
#endif
eta2InterpolationPosXY *interp=new eta2InterpolationPosXY(NC, NR, nsubpix, etabins, etamin, etamax);
eta3InterpolationPosXY *interp3=new eta3InterpolationPosXY(NC, NR, nsubpix, etabins, eta3min, eta3max);
noInterpolation *dummy=new noInterpolation(NC, NR, nsubpix);
noInterpolation *nointerp=new noInterpolation(NC, NR, nsubpix);
noInterpolation *mult=new noInterpolation(NC, NR, nsubpix);
//etaInterpolationAdaptiveBins *interp=new etaInterpolationAdaptiveBins (NC, NR, nsubpix, etabins, etamin, etamax);
//etaInterpolationRandomBins *interp=new etaInterpolationRandomBins (NC, NR, nsubpix, etabins, etamin, etamax);
//#ifndef FF
cout << "read ff " << argv[2] << endl;
interp->readFlatField(argv[2]);
sprintf(fname,"%s_eta2.tiff",argv[2]);
interp->readFlatField(fname);
interp->prepareInterpolation(ok);
sprintf(fname,"%s_eta3.tiff",argv[2]);
interp3->readFlatField(fname);
interp3->prepareInterpolation(ok);
//#endif
int *img;
@ -70,20 +118,99 @@ int main(int argc, char *argv[]) {
#endif
for (int irun=runmin; irun<runmax; irun++) {
int irun;
for (irun=runmin; irun<runmax; irun++) {
sprintf(infname,argv[1],irun);
#ifndef MANYFILES
sprintf(outfname,argv[3],irun);
#endif
f=fopen(infname,"r");
if (f) {
cout << infname << endl;
nframes=0;
f0=-1;
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);
totph++;
if (lastframe!=cl.iframe) {
lastframe=cl.iframe;
// cout << cl.iframe << endl;
// f0=cl.iframe;
if (nframes==0) f0=lastframe;
nframes++;
#ifdef MANYFILES
if (nframes%ff==0) {
sprintf(outfname,argv[3],irun,nframes-ff);
cout << outfname << endl;
interp->writeInterpolatedImage(outfname);
interp->clearInterpolatedImage();
}
#endif
}
#ifdef XTALK
if ((cl.x+1)%25!=0) {
for (int ix=-1; ix<2; ix++) {
for (int iy=-1; iy<2; iy++) {
old_val[iy+1][ix+1]=cl.get_data(ix,iy);
if (ix>=0) {
new_val[iy+1][ix+1]=old_val[iy+1][ix+1]-old_val[iy+1][ix]*xcorr;
cl.set_data(new_val[iy+1][ix+1],ix,iy);
}
}
}
}
#endif
quad=interp->calcQuad(cl.get_cluster(), sum, totquad, sDum);
if (sum>cmin && totquad/sum>0.8 && totquad/sum<1.2 && totquad<cmax && quad<cmax) {
nph++;
// if (sum>200 && sum<580) {
// interp->getInterpolatedPosition(cl.x,cl.y, totquad,quad,cl.get_cluster(),int_x, int_y);
interp->getInterpolatedPosition(cl.x,cl.y, cl.get_cluster(),int_x, int_y);
interp->addToImage(int_x, int_y);
interp3->getInterpolatedPosition(cl.x,cl.y, cl.get_cluster(),int3_x, int3_y);
interp3->addToImage(int3_x, int3_y);
nointerp->getInterpolatedPosition(cl.x,cl.y, cl.get_cluster(),noint_x, noint_y);
nointerp->addToImage(noint_x, noint_y);
d_x= (int_x-int3_x)*25.;
d_y= (int_y-int3_y)*25.;
dummy->calcEta(totquad, sDum, etax, etay);
xx=int_x;
yy=int_y;
if (etax<0.1 || etax>0.9) xx=int3_x;
if (etay<0.1 || etay>0.9) yy=int3_y;
dummy->addToImage(xx,yy);
if (d_x>res || d_x<-res || d_y>res || d_y<-res) {
badph++;
// cout << "delta (um): "<< d_x << " " << d_y << " " << cl.x << " " << cl.y << endl;
// cout << sum << " " << totquad << " " << etax << " "<< etay << endl;
// //cout<< int_x << " " << int_y << " " << int3_x << " " << int3_y << endl;
}
mult->addToImage(noint_x, noint_y);
if (nph%1000000==0) cout << nph << endl;
if (nph%100000000==0) {
sprintf(outfname,"%s_inteta2.tiff", argv[3]);
interp->writeInterpolatedImage(outfname);
sprintf(outfname,"%s_inteta3.tiff", argv[3]);
interp3->writeInterpolatedImage(outfname);
sprintf(outfname,"%s_mix.tiff", argv[3]);
dummy->writeInterpolatedImage(outfname);
sprintf(outfname,"%s_noint.tiff", argv[3]);
nointerp->writeInterpolatedImage(outfname);
sprintf(outfname,"%s_mult.tiff", argv[3]);
mult->writeInterpolatedImage(outfname);
}
} else {
mult->getInterpolatedPosition(cl.x,cl.y, cl.get_cluster(),int_x, int_y);
for (int imult=0; imult<2.*sum/(cmax+cmin); imult++) mult->addToImage(int_x, int_y);
}
}
fclose(f);
#ifdef FF
@ -147,8 +274,22 @@ int main(int argc, char *argv[]) {
}
#endif
#ifndef FF
#ifdef MANYFILES
sprintf(outfname,argv[3],irun,nframes-ff);
cout << outfname << endl;
#endif
sprintf(outfname,"%s_inteta2.tiff", argv[3]);
interp->writeInterpolatedImage(outfname);
img=interp->getInterpolatedImage();
sprintf(outfname,"%s_inteta3.tiff", argv[3]);
interp3->writeInterpolatedImage(outfname);
sprintf(outfname,"%s_mix.tiff", argv[3]);
dummy->writeInterpolatedImage(outfname);
sprintf(outfname,"%s_mult.tiff", argv[3]);
mult->writeInterpolatedImage(outfname);
#ifndef MANYFILES
img=interp->getInterpolatedImage();
for (ix=0; ix<NC; ix++) {
for (iy=0; iy<NR; iy++) {
for (isx=0; isx<nsubpix; isx++) {
@ -158,18 +299,24 @@ int main(int argc, char *argv[]) {
}
}
}
#endif
#endif
cout << "Read " << nframes << " frames (first frame: " << f0 << " last frame: " << lastframe << " delta:" << lastframe-f0 << ")"<<endl;
interp->clearInterpolatedImage();
interp3->clearInterpolatedImage();
dummy->clearInterpolatedImage();
mult->clearInterpolatedImage();
}
} else
cout << "could not open file " << infname << endl;
}
cout << irun << " " << runmax << endl;
#ifndef MANYFILES
sprintf(outfname,argv[3],11111);
WriteToTiff(totimg, outfname,NC*nsubpix,NR*nsubpix);
#endif
cout << "Filled " << nph << " (/"<< totph <<") of which " << badph << " badly interpolated " << endl;
return 0;
}

View File

@ -11,6 +11,7 @@ using namespace std;
#define NC 400
#define NR 400
#define XTALK
int main(int argc, char *argv[]) {
/**
@ -20,6 +21,7 @@ int main(int argc, char *argv[]) {
int nsubpix=10;
int etabins=nsubpix*100;
double etamin=-1, etamax=2;
double eta3min=-2, eta3max=2;
int quad;
double sum, totquad;
double sDum[2][2];
@ -27,16 +29,29 @@ int main(int argc, char *argv[]) {
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;
int iph=0;
if (argc<7) {
cout << "Wrong usage! Should be: "<< argv[0] << " infile " << " outfile runmin runmax cmin cmax" << endl;
return 1;
}
etaInterpolationPosXY *interp=new etaInterpolationPosXY(NR, NC, nsubpix, etabins, etamin, etamax);
eta2InterpolationPosXY *interp2=new eta2InterpolationPosXY(NR, NC, nsubpix, etabins, etamin, etamax);
cout << "###########"<< endl;
eta3InterpolationPosXY *interp3=new eta3InterpolationPosXY(NR, NC, nsubpix, etabins, eta3min, eta3max);
// cout << eta3min << " " << eta3max << endl;
runmin=atoi(argv[3]);
runmax=atoi(argv[4]);
double cmin=atof(argv[5]); //200
double cmax=atof(argv[6]); //3000
#ifdef XTALK
int old_val[3][3];
int new_val[3][3];
double xcorr=0.04;
// int ix=0;
#endif
FILE *f;
for (int i=runmin; i<runmax; i++) {
@ -45,18 +60,55 @@ int main(int argc, char *argv[]) {
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);
#ifdef XTALK
if ((cl.x+1)%25!=0) {
for (int ix=-1; ix<2; ix++) {
for (int iy=-1; iy<2; iy++) {
old_val[iy+1][ix+1]=cl.get_data(ix,iy);
if (ix>=0) {
new_val[iy+1][ix+1]=old_val[iy+1][ix+1]-old_val[iy+1][ix]*xcorr;
cl.set_data(new_val[iy+1][ix+1],ix,iy);
}
}
}
}
#endif
quad=interp2->calcQuad(cl.get_cluster(), sum, totquad, sDum);
if (sum>cmin && totquad/sum>0.8 && totquad/sum<1.2 && totquad<cmax && quad<cmax) {
/* Cross talk corrections !!! */
// for (int ix=0; ix<2; ix++) {
interp2->addToFlatField(cl.get_cluster(),etax, etay);
// if (etax>0.49 && etax<0.51 && etay>0.49 && etay<0.51 ) {
// cout << cl.y << " " << cl.x << " " << quad << " "<< totquad << " " <<sum << endl;
// }
interp3->addToFlatField(cl.get_cluster(),etax, etay);
iph++;
if (iph%1000000==0) cout << iph << endl;
if (iph%100000000==0) {
sprintf(fname,"%s_eta2.tiff",argv[2]);
interp2->writeFlatField(fname);
sprintf(fname,"%s_eta3.tiff",argv[2]);
interp3->writeFlatField(fname);
}
// if (iph>1E8) break;
}
// }
}
fclose(f);
interp->writeFlatField(argv[2]);
}
else cout << "could not open file " << fname << endl;
}
interp->writeFlatField(argv[2]);
sprintf(fname,"%s_eta2.tiff",argv[2]);
interp2->writeFlatField(fname);
sprintf(fname,"%s_eta3.tiff",argv[2]);
interp3->writeFlatField(fname);
return 0;
}

View File

@ -1,160 +0,0 @@
#include <vector>
#include <string>
#include <sstream>
#include <iomanip>
#include <fstream>
#include <stdio.h>
//#include <deque>
//#include <list>
//#include <queue>
#include <fstream>
#include "moench03Ctb10GbT1Data.h"
#include "interpolatingDetector.h"
#include "etaInterpolationPosXY.h"
#include "linearInterpolation.h"
#include "noInterpolation.h"
#include "multiThreadedDetector.h"
#include <ctime>
#define NC 400
#define NR 400
#include "tiffIO.h"
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 nthreads=3;
int nph, nph1;
single_photon_hit clusters[NR*NC];
// cout << "hits "<< endl;
int etabins=550;
double etamin=-1, etamax=2;
int nsubpix=4;
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);
//noInterpolation *interp=new noInterpolation(NC, NR, nsubpix);
interp->readFlatField("/scratch/eta_100.tiff",etamin,etamax);
interpolatingDetector *filter=new interpolatingDetector(decoder,interp, 5, 1, 0, 1000, 10);
filter->readPedestals("/scratch/ped_100.tiff");
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();
multiThreadedDetector *mt=new multiThreadedDetector(filter,nthreads,100);
nph=0;
nph1=0;
//int np;
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;
mt->setFrameMode(eFrame);
mt->prepareInterpolation(ok);
mt->StartThreads();
mt->popFree(buff);
while ((decoder->readNextFrame(filebin, iFrame, np, buff)) && nf<1.5E4) {
if (nf<9E3)
;
else {
// if (nf>1.1E4 && ok==0) {
// mt->prepareInterpolation(ok);
// mt->setFrameMode(eFrame);
// //ok=1;
// }
mt->pushData(buff);
mt->nextThread();
// cout << " " << (void*)buff;
mt->popFree(buff);
// 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;
mt->StopThreads();
char tit[10000];
sprintf(tit,"/scratch/int_image_mt%d.tiff",nthreads);
mt->writeInterpolatedImage(tit);
// delete [] etah;
// delete interp;
//delete decoder;
//cout << "Read " << nf << " frames" << endl;
return NULL;
}
int main(int argc, char *argv[]){
moenchProcessFrame();
}

View File

@ -1,267 +0,0 @@
#include "../single_photon_hit.h"
//#include "etaVEL/etaInterpolationPosXY.h"
#include <TH1F.h>
#include <TH2F.h>
#include <TCanvas.h>
#include <iostream>
using namespace std;
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);
TH1F *hf=new TH1F("hf","hf",1000,0,10E6);
//TH2F *hff=new TH2F("hff","hff",ns, -0.5, 0.5, ns, -0.5, +0.5);
TH1F *hsp=new TH1F("hsp",fname,500,0,10000);
// TH1F *hsp1=new TH1F("hsp1",fname,500,0,10000);
// TH1F *hsp2=new TH1F("hsp2",fname,500,0,1000);
// TH1F *hsp3=new TH1F("hsp3",fname,500,0,1000);
// hsp1->SetLineColor(2);
// hsp2->SetLineColor(3);
// hsp3->SetLineColor(4);
TCanvas *c=new TCanvas();
h2->Draw("colz");
TCanvas *c1=new TCanvas();
hsp->Draw();
c1->SetLogy(kTRUE);
// hsp1->Draw("same");
// hsp2->Draw("same");
// hsp3->Draw("same");
TCanvas *c2=new TCanvas();
hf->Draw();
// hint->Draw("colz");
//c2->SetLogz(kTRUE);
single_photon_hit cl(3,3);
double tot;
int w;
double phw=340, phs=62;
int f0=-1;
double tl, bl, tr, br, qt;
int iimage=0;
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;*/
tl=0; tr=0; bl=0; br=0;
for (int ix=-1; ix<2; ix++)
for (int iy=-1; iy<2; iy++){
tot+=cl.get_data(ix,iy);
if (ix<=0 && iy<=0) bl+=cl.get_data(ix,iy);
if (ix<=0 && iy>=0) tl+=cl.get_data(ix,iy);
if (ix>=0 && iy<=0) br+=cl.get_data(ix,iy);
if (ix>=0 && iy>=0) tr+=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);
}
qt=bl;
if (br>qt) qt=br;
if (tl>qt) qt=tl;
if (tr>qt) qt=tr;
/*
px=(-left+right)/tot;
py=(-bottom+top)/tot;*/
//max at 340
//if (tot>200) {
w=1;
if (qt>1000) {
if (qt/tot>0.8 && qt/tot<1.2){
if (f0<0)
f0=cl.iframe;
hf->Fill(cl.iframe-f0);
// if (qt>540) w++;
// if (qt>820) w++;
//(tot+3.5*phs)/phw;
//} else
//w=0;
// if (w) {
// if (cl.y<350) {
// if (cl.y<100 || cl.y>300) {
// if (cl.x>150 && cl.x<250 && cl.y>200 && cl.y<250)
// hsp1->Fill(qt);
// else
hsp->Fill(qt);
// hsp2->Fill(cl.get_data(0,0));
// } else {
// hsp1->Fill(qt);
// hsp3->Fill(cl.get_data(0,0));
// }
//if (cl.x>160 && cl.x<260 && cl.y>30 && cl.y<80 && tot>0)
//if (w==1) {
// if (w==2)
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>1E7)
// break;
//}
// 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;
// }

View File

@ -51,8 +51,11 @@ class pedestalSubtraction {
/** sets the moving average */
virtual void setPedestalRMS(double rms) {stat.SetRMS(rms);}
/**sets/gets the number of samples for the moving average
\returns actual number of samples for the moving average
*/
virtual int getNumpedestals() {return stat.NumDataValues();};
private:
MovingStat stat; /**< approximated moving average struct */

View File

@ -58,7 +58,7 @@ public analogDetector<uint16_t> {
int sign=1,
commonModeSubtraction *cm=NULL,
int nped=1000,
int nd=100, int nnx=-1, int nny=-1, double *gm=NULL) : analogDetector<uint16_t>(d, sign, cm, nped, nnx, nny, gm), nDark(nd), eventMask(NULL),nSigma (nsigma), clusterSize(csize), clusterSizeY(csize), cluster(NULL), quad(UNDEFINED_QUADRANT), tot(0), quadTot(0) {
int nd=100, int nnx=-1, int nny=-1, double *gm=NULL) : analogDetector<uint16_t>(d, sign, cm, nped, nnx, nny, gm), nDark(nd), eventMask(NULL),nSigma (nsigma), clusterSize(csize), clusterSizeY(csize), clusters(NULL), quad(UNDEFINED_QUADRANT), tot(0), quadTot(0) {
@ -74,7 +74,7 @@ public analogDetector<uint16_t> {
// cluster=new single_photon_hit(clusterSize,clusterSizeY);
clusters=new single_photon_hit[nx*ny];
cluster=clusters;
// cluster=clusters;
setClusterSize(csize);
nphTot=0;
nphFrame=0;
@ -82,7 +82,7 @@ public analogDetector<uint16_t> {
/**
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; };
virtual ~singlePhotonDetector() {delete [] clusters; for (int i=0; i<ny; i++) delete [] eventMask[i]; delete [] eventMask; };
@ -109,7 +109,7 @@ public analogDetector<uint16_t> {
// cluster=new single_photon_hit(clusterSize,clusterSizeY);
clusters=new single_photon_hit[nx*ny];
cluster=clusters;
// cluster=clusters;
setClusterSize(clusterSize);
@ -147,8 +147,8 @@ public analogDetector<uint16_t> {
if (n%2==0)
n+=1;
clusterSize=n;
if (cluster)
delete cluster;
// if (clusters)
// delete [] clusters;
if (ny>clusterSize)
clusterSizeY=clusterSize;
else
@ -192,29 +192,37 @@ public analogDetector<uint16_t> {
double max=0, tl=0, tr=0, bl=0,br=0, v;
if (thr>=0) {
int cm=0;
if (cmSub) cm=1;
if (thr>0) {
cy=1;
cs=1;
ccs=1;
ccy=1;
}
if (iframe<nDark) {
//cout << "ped " << iframe << endl;
// cout << "ped " << iframe << endl;
//this already adds to common mode
addToPedestal(data);
return nph;
} else {
if (thr>0) {
newFrame();
if (cmSub) {
addToCommonMode(data);
}
for (int ix=xmin; ix<xmax; ix++) {
for (int iy=ymin; iy<ymax; iy++) {
val=subtractPedestal(data,ix,iy);
val=subtractPedestal(data,ix,iy, cm);
nn=analogDetector<uint16_t>::getNPhotons(data,ix,iy);
nn=val/tthr;//analogDetector<uint16_t>::getNPhotons(data,ix,iy);
nph[ix+nx*iy]+=nn;
rest[iy][ix]=(val-nn*tthr);
nphFrame+=nn;
nphTot+=nn;
}
}
// }
@ -235,10 +243,10 @@ public analogDetector<uint16_t> {
for (int ir=-(clusterSizeY/2); ir<(clusterSizeY/2)+1; ir++) {
for (int ic=-(clusterSize/2); ic<(clusterSize/2)+1; ic++) {
if ((iy+ir)>=0 && (iy+ir)<ny && (ix+ic)>=0 && (ix+ic)<nx) {
//cluster->set_data(rest[iy+ir][ix+ic], ic, ir);
//clusters->set_data(rest[iy+ir][ix+ic], ic, ir);
v=rest[iy+ir][ix+ic];//cluster->get_data(ic,ir);
v=rest[iy+ir][ix+ic];//clusters->get_data(ic,ir);
tot+=v;
if (ir<=0 && ic<=0)
@ -278,14 +286,15 @@ public analogDetector<uint16_t> {
if (max>tthr || tot>sqrt(ccy*ccs)*tthr || quadTot>sqrt(cy*cs)*tthr) {
eventMask[iy][ix]=PHOTON;
nph[ix+nx*iy]++;
nphFrame+=nph[ix+nx*iy];
nphTot+=nph[ix+nx*iy];
nphFrame++;
nphTot++;
}
}
}
}
cout << iframe << " " << nph << endl;
// cout << iframe << " " << nphFrame << " " << nphTot << endl;
//cout << iframe << " " << nph << endl;
} else return getClusters(data, nph);
}
return NULL;
@ -329,10 +338,10 @@ public analogDetector<uint16_t> {
eventMask[iy][ix]=PEDESTAL;
cluster->x=ix;
cluster->y=iy;
cluster->rms=getPedestalRMS(ix,iy);
cluster->ped=getPedestal(ix,iy, cm);
clusters->x=ix;
clusters->y=iy;
clusters->rms=getPedestalRMS(ix,iy);
clusters->ped=getPedestal(ix,iy, cm);
for (int ir=-(clusterSizeY/2); ir<(clusterSizeY/2)+1; ir++) {
@ -341,8 +350,8 @@ public analogDetector<uint16_t> {
v=subtractPedestal(data, ix+ic, iy+ir);
cluster->set_data(v, ic, ir);
// v=cluster->get_data(ic,ir);
clusters->set_data(v, ic, ir);
// v=clusters->get_data(ic,ir);
tot+=v;
if (ir<=0 && ic<=0)
bl+=v;
@ -357,7 +366,7 @@ public analogDetector<uint16_t> {
max=v;
}
if (ir==0 && ic==0) {
if (v<-nSigma*cluster->rms)
if (v<-nSigma*clusters->rms)
eventMask[iy][ix]=NEGATIVE_PEDESTAL;
}
}
@ -378,8 +387,8 @@ public analogDetector<uint16_t> {
quadTot=tr;
}
if (max>nSigma*cluster->rms || tot>sqrt(clusterSizeY*clusterSize)*nSigma*cluster->rms || quadTot>cy*cs*nSigma*cluster->rms) {
if (cluster->get_data(0,0)>=max) {
if (max>nSigma*clusters->rms || tot>sqrt(clusterSizeY*clusterSize)*nSigma*clusters->rms || quadTot>cy*cs*nSigma*clusters->rms) {
if (clusters->get_data(0,0)>=max) {
eventMask[iy][ix]=PHOTON_MAX;
} else {
eventMask[iy][ix]=PHOTON;
@ -420,7 +429,8 @@ int *getClusters(char *data, int *ph=NULL) {
int ir, ic;
double max=0, tl=0, tr=0, bl=0,br=0, *v, vv;
int cm=0;
if (cmSub) cm=1;
if (ph==NULL)
ph=image;
@ -429,6 +439,13 @@ int *getClusters(char *data, int *ph=NULL) {
return 0;
}
newFrame();
if (cm)
addToCommonMode(data);
for (int ix=xmin; ix<xmax; ix++) {
for (int iy=ymin; iy<ymax; iy++) {
@ -447,14 +464,14 @@ int *getClusters(char *data, int *ph=NULL) {
(clusters+nph)->rms=getPedestalRMS(ix,iy);
cluster=clusters+nph;
// cluster=clusters+nph;
for (int ir=-(clusterSizeY/2); ir<(clusterSizeY/2)+1; ir++) {
for (int ic=-(clusterSize/2); ic<(clusterSize/2)+1; ic++) {
if ((iy+ir)>=iy && (iy+ir)<ny && (ix+ic)>=ix && (ix+ic)<nx) {
val[iy+ir][ix+ic]=subtractPedestal(data,ix+ic,iy+ir);
val[iy+ir][ix+ic]=subtractPedestal(data,ix+ic,iy+ir, cm);
}
v=&(val[iy+ir][ix+ic]);
@ -473,7 +490,7 @@ int *getClusters(char *data, int *ph=NULL) {
if (ir==0 && ic==0) {
if (*v<-nSigma*cluster->rms)
if (*v<-nSigma*(clusters+nph)->rms)
eventMask[iy][ix]=NEGATIVE_PEDESTAL;
}
@ -494,7 +511,7 @@ int *getClusters(char *data, int *ph=NULL) {
(clusters+nph)->quadTot=tr;
}
if (max>nSigma*cluster->rms || tot>sqrt(clusterSizeY*clusterSize)*nSigma*cluster->rms || ((clusters+nph)->quadTot)>sqrt(cy*cs)*nSigma*cluster->rms) {
if (max>nSigma*(clusters+nph)->rms || tot>sqrt(clusterSizeY*clusterSize)*nSigma*(clusters+nph)->rms || ((clusters+nph)->quadTot)>sqrt(cy*cs)*nSigma*(clusters+nph)->rms) {
if (val[iy][ix]>=max) {
eventMask[iy][ix]=PHOTON_MAX;
(clusters+nph)->tot=tot;
@ -507,6 +524,7 @@ int *getClusters(char *data, int *ph=NULL) {
(clusters+nph)->set_data(val[iy+ir][ix+ic],ic,ir);
}
}
// cout << (clusters+nph)->iframe << " " << ix << " " << nph << " " << tot << " " << (clusters+nph)->quadTot << endl;
nph++;
image[iy*nx+ix]++;
@ -514,7 +532,7 @@ int *getClusters(char *data, int *ph=NULL) {
eventMask[iy][ix]=PHOTON;
}
} else if (eventMask[iy][ix]==PEDESTAL) {
addToPedestal(data,ix,iy);
addToPedestal(data,ix,iy,cm);
}
@ -559,7 +577,7 @@ int *getClusters(char *data, int *ph=NULL) {
\param ir y coordinate (center is (0,0))
\returns cluster element
*/
double getClusterElement(int ic, int ir=0){return cluster->get_data(ic,ir);};
double getClusterElement(int ic, int ir=0){return clusters->get_data(ic,ir);};
/** returns event mask for the given pixel
\param ic x coordinate (center is (0,0))
@ -581,25 +599,25 @@ int *getClusters(char *data, int *ph=NULL) {
if (iFrame)
tall->Branch("iFrame",iFrame,"iframe/I");
else
tall->Branch("iFrame",&(cluster->iframe),"iframe/I");
tall->Branch("iFrame",&(clusters->iframe),"iframe/I");
tall->Branch("x",&(cluster->x),"x/I");
tall->Branch("y",&(cluster->y),"y/I");
tall->Branch("x",&(clusters->x),"x/I");
tall->Branch("y",&(clusters->y),"y/I");
char tit[100];
sprintf(tit,"data[%d]/D",clusterSize*clusterSizeY);
tall->Branch("data",cluster->data,tit);
tall->Branch("pedestal",&(cluster->ped),"pedestal/D");
tall->Branch("rms",&(cluster->rms),"rms/D");
tall->Branch("tot",&(cluster->tot),"tot/D");
tall->Branch("quadTot",&(cluster->quadTot),"quadTot/D");
tall->Branch("quad",&(cluster->quad),"quad/I");
tall->Branch("data",clusters->data,tit);
tall->Branch("pedestal",&(clusters->ped),"pedestal/D");
tall->Branch("rms",&(clusters->rms),"rms/D");
tall->Branch("tot",&(clusters->tot),"tot/D");
tall->Branch("quadTot",&(clusters->quadTot),"quadTot/D");
tall->Branch("quad",&(clusters->quad),"quad/I");
return tall;
};
#else
/** write cluster to filer
\param f file pointer
*/
void writeCluster(FILE* f){cluster->write(f);};
void writeCluster(FILE* f){clusters->write(f);};
/**
write clusters to file
@ -643,7 +661,7 @@ void writeClusters(FILE *f){for (int i=0; i<nphFrame; i++) (clusters+i)->write(f
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 *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 */

View File

@ -1,5 +1,7 @@
#ifndef SINGLE_PHOTON_HIT_H
#define SINGLE_PHOTON_HIT_H
#include <stdio.h>
#include <stdint.h>
typedef double double32_t;
typedef float float32_t;
@ -26,7 +28,7 @@ class single_photon_hit {
\param nx cluster size in x direction
\param ny cluster size in y direction (defaults to 1 for 1D detectors)
*/
single_photon_hit(int nx=3, int ny=3): dx(nx), dy(ny) {data=new double[dx*dy];};
single_photon_hit(int nx=3, int ny=3): dx(nx), dy(ny) {data=new int[dx*dy];};
~single_photon_hit(){delete [] data;}; /**< destructor, deletes the data array */
@ -35,8 +37,9 @@ class single_photon_hit {
*/
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);
if (fwrite((void*)this, 1, sizeof(int)+2*sizeof(int16_t), myFile))
return fwrite((void*)data, 1, dx*dy*sizeof(int), myFile);
return 0;
};
@ -47,8 +50,8 @@ class single_photon_hit {
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);
if (fread((void*)this, 1, sizeof(int)+2*sizeof(int16_t), myFile))
return fread((void*)data, 1, dx*dy*sizeof(int), myFile);
return 0;
};
@ -60,7 +63,7 @@ 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 set_cluster_size(int nx, int ny) {if (nx>0) dx=nx; if (ny>0) dy=ny; delete [] data; data=new int[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;};
@ -71,11 +74,11 @@ class single_photon_hit {
\returns value of the cluster element
*/
double get_data(int ix, int iy=0){return data[(iy+dy/2)*dx+ix+dx/2];};
double *get_cluster() {return data;};
int *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 */
int16_t x; /**< x-coordinate of the center of hit */
int16_t 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*/
double tot; /**< sum of the 3x3 cluster */
@ -83,7 +86,7 @@ class single_photon_hit {
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 */
int *data; /**< pointer to data */
};

View File

@ -0,0 +1,98 @@
#ifndef SINGLE_PHOTON_HIT_H
#define SINGLE_PHOTON_HIT_H
#include <stdio.h>
#include <stdint.h>
typedef double double32_t;
typedef float float32_t;
typedef int int32_t;
#ifndef DEF_QUAD
#define DEF_QUAD
enum quadrant {
TOP_LEFT=0,
TOP_RIGHT=1,
BOTTOM_LEFT=2,
BOTTOM_RIGHT=3,
UNDEFINED_QUADRANT=-1
};
#endif
class single_photon_hit_double {
/** @short Structure for a single photon hit */
public:
/** constructor, instantiates the data array -- all class elements are public!
\param nx cluster size in x direction
\param ny cluster size in y direction (defaults to 1 for 1D detectors)
*/
single_photon_hit_double(int nx=3, int ny=3): dx(nx), dy(ny) {
data=new double[dx*dy];
};
~single_photon_hit_double(){delete [] data;}; /**< destructor, deletes the data array */
/** binary write to file of all elements of the structure, except size of the cluster
\param myFile file descriptor
*/
size_t write(FILE *myFile) {
//fwrite((void*)this, 1, 3*sizeof(int)+4*sizeof(double)+sizeof(quad), myFile);
if (fwrite((void*)this, 1, sizeof(int)+2*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
*/
size_t read(FILE *myFile) {
//fread((void*)this, 1, 3*sizeof(int)+4*sizeof(double)+sizeof(quad), myFile);
if (fread((void*)this, 1, sizeof(int)+2*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)
\param v value to be set
\param ix coordinate x within the cluster (center is (0,0))
\param iy coordinate y within the cluster (center is (0,0))
*/
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)
\param ix coordinate x within the cluster (center is (0,0))
\param iy coordinate y within the cluster (center is (0,0))
\returns value of the cluster element
*/
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*/
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 */
};
#endif

View File

@ -0,0 +1,93 @@
#ifndef SINGLE_PHOTON_HIT_H
#define SINGLE_PHOTON_HIT_H
#include <stdio.h>
#include <stdint.h>
typedef double double32_t;
typedef float float32_t;
typedef int int32_t;
#ifndef DEF_QUAD
#define DEF_QUAD
enum quadrant {
TOP_LEFT=0,
TOP_RIGHT=1,
BOTTOM_LEFT=2,
BOTTOM_RIGHT=3,
UNDEFINED_QUADRANT=-1
};
#endif
class single_photon_hit {
/** @short Structure for a single photon hit */
public:
/** constructor, instantiates the data array -- all class elements are public!
\param nx cluster size in x direction
\param ny cluster size in y direction (defaults to 1 for 1D detectors)
*/
single_photon_hit(int nx=3, int ny=3): dx(nx), dy(ny) {data=new int[dx*dy];};
~single_photon_hit(){delete [] data;}; /**< destructor, deletes the data array */
/** binary write to file of all elements of the structure, except size of the cluster
\param myFile file descriptor
*/
size_t write(FILE *myFile) {
//fwrite((void*)this, 1, 3*sizeof(int)+4*sizeof(double)+sizeof(quad), myFile);
if (fwrite((void*)this, 1, sizeof(int)+2*sizeof(int16_t), myFile))
return fwrite((void*)data, 1, dx*dy*sizeof(int), 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
*/
size_t read(FILE *myFile) {
//fread((void*)this, 1, 3*sizeof(int)+4*sizeof(double)+sizeof(quad), myFile);
if (fread((void*)this, 1, sizeof(int)+2*sizeof(int16_t), myFile))
return fread((void*)data, 1, dx*dy*sizeof(int), 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)
\param v value to be set
\param ix coordinate x within the cluster (center is (0,0))
\param iy coordinate y within the cluster (center is (0,0))
*/
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 int[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)
\param ix coordinate x within the cluster (center is (0,0))
\param iy coordinate y within the cluster (center is (0,0))
\returns value of the cluster element
*/
double get_data(int ix, int iy=0){return data[(iy+dy/2)*dx+ix+dx/2];};
int *get_cluster() {return data;};
int iframe; /**< frame number */
int16_t x; /**< x-coordinate of the center of hit */
int16_t 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*/
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 */
int *data; /**< pointer to data */
};
#endif

View File

@ -298,7 +298,7 @@ class slsDetectorData {
*/
virtual int getFrameNumber(char *buff)=0;
/**
Returns the packet number for the given dataset. purely virtual func

View File

@ -1,5 +1,6 @@
#ifndef MY_TIFF_IO_H
#include "tiffIO.h"
#endif
#include<iostream>
using namespace std;
// #undef cbf_failnez

View File

@ -1,9 +1,18 @@
Path: slsDetectorsPackage/slsDetectorGui
URL: origin git@github.com:slsdetectorgroup/slsDetectorPackage.git
Repository Root: origin git@github.com:slsdetectorgroup/slsDetectorPackage.git
<<<<<<< HEAD
Repsitory UUID: ab06c33107ecfeb4741d49407903ff80286cf75b
Revision: 492
Branch: developer
Last Changed Author: Anna_Bergamaschi
Last Changed Rev: 3731
Last Changed Date: 2018-03-15 12:27:06.000000002 +0100 ./src/qTabMeasurement.cpp
=======
Repsitory UUID: fe2ba8621b33bc00f51c8cd4d33e98b7d4ebfa41
Revision: 493
Branch: developer
Last Changed Author: Dhanya_Thattil
Last Changed Rev: 3747
Last Changed Date: 2018-03-27 17:30:53.000000002 +0200 ./include/qTabMeasurement.h
>>>>>>> 7cd35f24b87501374fbaf45693a2adf16dfae3e3

View File

@ -1,4 +1,5 @@
#define GITURL "git@github.com:slsdetectorgroup/slsDetectorPackage.git"
#define GITREPUUID "fe2ba8621b33bc00f51c8cd4d33e98b7d4ebfa41"
#define GITAUTH "Dhanya_Thattil"
#define GITREV 0x3747

View File

@ -1,9 +1,18 @@
Path: slsDetectorsPackage/slsDetectorSoftware
URL: origin git@github.com:slsdetectorgroup/slsDetectorPackage.git
Repository Root: origin git@github.com:slsdetectorgroup/slsDetectorPackage.git
<<<<<<< HEAD
Repsitory UUID: ab06c33107ecfeb4741d49407903ff80286cf75b
Revision: 1846
Branch: developer
Last Changed Author: Anna_Bergamaschi
Last Changed Rev: 3731
Last Changed Date: 2018-03-15 12:30:34.000000002 +0100 ./threadFiles/ThreadPool.o
=======
Repsitory UUID: b8bdbf4da61f95b88893b02ddabc2491b16fa10f
Revision: 1852
Branch: developer
Last Changed Author: Dhanya_Thattil
Last Changed Rev: 3746
Last Changed Date: 2018-03-27 10:47:02.000000002 +0200 ./slsDetector/slsDetector.o
>>>>>>> 7cd35f24b87501374fbaf45693a2adf16dfae3e3

View File

@ -3,8 +3,10 @@
#define CSP0 0x20200000
#define MEM_SIZE 0x100000
#define MEM_MAP_SHIFT 1
#include <sys/types.h>
#include <sys/types.h>
int mapCSP0(void);
u_int16_t bus_r16(u_int32_t offset);

View File

@ -1289,6 +1289,7 @@ int64_t setExposureTime(int64_t value){
/* time is in ns */
if (value!=-1)
value*=(1E-3*clkDivider[0]);//(1E-9*CLK_FREQ);
setPatternWaitTime(0,value);
return set64BitReg(value,SET_EXPTIME_LSB_REG, SET_EXPTIME_MSB_REG)/(1E-3*clkDivider[0]);//(1E-9*CLK_FREQ);
}

View File

@ -1,9 +1,9 @@
Path: slsDetectorsPackage/slsDetectorSoftware/jctbDetectorServer
URL: origin git@github.com:slsdetectorgroup/slsDetectorPackage.git
Repository Root: origin git@github.com:slsdetectorgroup/slsDetectorPackage.git
Repsitory UUID: 3da6a6df6556312f7467407a8b5691bdc478424e
Revision: 21
Repsitory UUID: c35203ebfb35ba605eaff601ce46e4d29c1f5690
Revision: 24
Branch: developer
Last Changed Author: Dhanya_Thattil
Last Changed Rev: 3597
Last Changed Date: 2018-02-07 10:44:13.000000002 +0100 ./server_funcs.h
Last Changed Author: Anna_Bergamaschi
Last Changed Rev: 3761
Last Changed Date: 2018-04-27 14:46:10.176661554 +0200 ./server_funcs.h

View File

@ -1,6 +1,6 @@
#define GITURL "git@github.com:slsdetectorgroup/slsDetectorPackage.git"
#define GITREPUUID "3da6a6df6556312f7467407a8b5691bdc478424e"
#define GITAUTH "Dhanya_Thattil"
#define GITREV 0x3597
#define GITDATE 0x20180207
#define GITREPUUID "c35203ebfb35ba605eaff601ce46e4d29c1f5690"
#define GITAUTH "Anna_Bergamaschi"
#define GITREV 0x3761
#define GITDATE 0x20180427
#define GITBRANCH "developer"

282
slsDetectorSoftware/jctbDetectorServer/registers_m.h Executable file → Normal file
View File

@ -20,7 +20,7 @@
#define FPGA_INIT_ADDR 0xb0000000
//#ifdef JUNGFRAU_DHANYA
#define POWER_ON_REG 0x5e<<11
#define POWER_ON_REG 0x5e << MEM_MAP_SHIFT
// Pwr_I2C_SDA <= PowerReg_s(1) when PowerReg_s(3)='1' else 'Z';
// Pwr_I2C_SCL <= PowerReg_s(0) when PowerReg_s(2)='1' else 'Z';
@ -29,71 +29,71 @@
#define PWR_I2C_SCL_EN_BIT 2
#define PWR_I2C_SDA_EN_BIT 3
#define POWER_STATUS_REG 41<<11
#define POWER_STATUS_REG 41 << MEM_MAP_SHIFT
#define ADCREG1 0x08
#define ADCREG2 0x14//20
#define ADCREG3 0x4
#define ADCREG4 0x5
#define ADCREG_VREFS 24
#define DBIT_PIPELINE_REG 89<<11 //0x59 same PATTERN_N_LOOP2_REG
#define MEM_MACHINE_FIFOS_REG 79<<11 //from gotthard
#define CONFGAIN_REG 93<<11 //from gotthard
#define ADC_PIPELINE_REG 66<<11 //0x42 same as ADC_OFFSET_REG
#define DBIT_PIPELINE_REG 89 << MEM_MAP_SHIFT //0x59 same PATTERN_N_LOOP2_REG
#define MEM_MACHINE_FIFOS_REG 79 << MEM_MAP_SHIFT //from gotthard
#define CONFGAIN_REG 93 << MEM_MAP_SHIFT //from gotthard
#define ADC_PIPELINE_REG 66 << MEM_MAP_SHIFT //0x42 same as ADC_OFFSET_REG
//#endif
//#define ADC_OFFSET_REG 93<<11 //same as DAQ_REG
#define ADC_INVERSION_REG 67<<11
//#define ADC_OFFSET_REG 93 << MEM_MAP_SHIFT //same as DAQ_REG
#define ADC_INVERSION_REG 67 << MEM_MAP_SHIFT
#define DAC_REG 64<<11//0x17<<11// control the dacs
#define DAC_REG 64 << MEM_MAP_SHIFT//0x17 << MEM_MAP_SHIFT// control the dacs
//ADC
#define ADC_WRITE_REG 65<<11//0x18<<11
//#define ADC_SYNC_REG 66<<11//0x19<<11
//#define HV_REG 67<<11//0x20<<11
#define ADC_WRITE_REG 65 << MEM_MAP_SHIFT//0x18 << MEM_MAP_SHIFT
//#define ADC_SYNC_REG 66 << MEM_MAP_SHIFT//0x19 << MEM_MAP_SHIFT
//#define HV_REG 67 << MEM_MAP_SHIFT//0x20 << MEM_MAP_SHIFT
//#define MUTIME_REG 0x1a<<11
//#define MUTIME_REG 0x1a << MEM_MAP_SHIFT
//temperature
#define TEMP_IN_REG 0x1b<<11
#define TEMP_OUT_REG 0x1c<<11
#define TEMP_IN_REG 0x1b << MEM_MAP_SHIFT
#define TEMP_OUT_REG 0x1c << MEM_MAP_SHIFT
//configure MAC
#define TSE_CONF_REG 0x1d<<11
#define ENET_CONF_REG 0x1e<<11
//#define WRTSE_SHAD_REG 0x1f<<11
#define TSE_CONF_REG 0x1d << MEM_MAP_SHIFT
#define ENET_CONF_REG 0x1e << MEM_MAP_SHIFT
//#define WRTSE_SHAD_REG 0x1f << MEM_MAP_SHIFT
//HV
#define DUMMY_REG 68<<11//0x21<<11
#define FPGA_VERSION_REG 0<<11 //0x22<<11
#define PCB_REV_REG 0<<11
#define FIX_PATT_REG 1<<11 //0x23<<11
#define CONTROL_REG 79<<11//0x24<<11
#define STATUS_REG 2<<11 //0x25<<11
#define CONFIG_REG 77<<11//0x26<<11
#define EXT_SIGNAL_REG 78<<11// 0x27<<11
//#define FPGA_SVN_REG 0x29<<11
#define DUMMY_REG 68 << MEM_MAP_SHIFT//0x21 << MEM_MAP_SHIFT
#define FPGA_VERSION_REG 0 << MEM_MAP_SHIFT //0x22 << MEM_MAP_SHIFT
#define PCB_REV_REG 0 << MEM_MAP_SHIFT
#define FIX_PATT_REG 1 << MEM_MAP_SHIFT //0x23 << MEM_MAP_SHIFT
#define CONTROL_REG 79 << MEM_MAP_SHIFT//0x24 << MEM_MAP_SHIFT
#define STATUS_REG 2 << MEM_MAP_SHIFT //0x25 << MEM_MAP_SHIFT
#define CONFIG_REG 77 << MEM_MAP_SHIFT//0x26 << MEM_MAP_SHIFT
#define EXT_SIGNAL_REG 78 << MEM_MAP_SHIFT// 0x27 << MEM_MAP_SHIFT
//#define FPGA_SVN_REG 0x29 << MEM_MAP_SHIFT
#define CHIP_OF_INTRST_REG 0x2A<<11
#define CHIP_OF_INTRST_REG 0x2A << MEM_MAP_SHIFT
//FIFO
#define LOOK_AT_ME_REG 3<<11 //0x28<<11
#define SYSTEM_STATUS_REG 4<<11
#define LOOK_AT_ME_REG 3 << MEM_MAP_SHIFT //0x28 << MEM_MAP_SHIFT
#define SYSTEM_STATUS_REG 4 << MEM_MAP_SHIFT
#define FIFO_DATA_REG 6<<11
#define FIFO_STATUS_REG 7<<11
#define FIFO_DATA_REG 6 << MEM_MAP_SHIFT
#define FIFO_STATUS_REG 7 << MEM_MAP_SHIFT
// constant FifoDigitalInReg_c : integer := 60;
#define FIFO_DIGITAL_DATA_LSB_REG 60<<11
#define FIFO_DIGITAL_DATA_MSB_REG 61<<11
#define FIFO_DIGITAL_DATA_LSB_REG 60 << MEM_MAP_SHIFT
#define FIFO_DIGITAL_DATA_MSB_REG 61 << MEM_MAP_SHIFT
#define FIFO_DATA_REG_OFF 0x50<<11 ///////
#define FIFO_DATA_REG_OFF 0x50 << MEM_MAP_SHIFT ///////
//to read back dac registers
//#define MOD_DACS1_REG 0x65<<11
//#define MOD_DACS2_REG 0x66<<11
//#define MOD_DACS3_REG 0x67<<11
//#define MOD_DACS1_REG 0x65 << MEM_MAP_SHIFT
//#define MOD_DACS2_REG 0x66 << MEM_MAP_SHIFT
//#define MOD_DACS3_REG 0x67 << MEM_MAP_SHIFT
//user entered
@ -102,135 +102,135 @@
#define GET_ACTUAL_TIME_LSB_REG 16<<11
#define GET_ACTUAL_TIME_MSB_REG 17<<11
#define GET_ACTUAL_TIME_LSB_REG 16 << MEM_MAP_SHIFT
#define GET_ACTUAL_TIME_MSB_REG 17 << MEM_MAP_SHIFT
#define GET_MEASUREMENT_TIME_LSB_REG 38<<11
#define GET_MEASUREMENT_TIME_MSB_REG 39<<11
#define GET_MEASUREMENT_TIME_LSB_REG 38 << MEM_MAP_SHIFT
#define GET_MEASUREMENT_TIME_MSB_REG 39 << MEM_MAP_SHIFT
#define SET_DELAY_LSB_REG 96<<11 //0x68<<11
#define SET_DELAY_MSB_REG 97<<11 //0x69<<11
#define GET_DELAY_LSB_REG 18<<11//0x6a<<11
#define GET_DELAY_MSB_REG 19<<11//0x6b<<11
#define SET_DELAY_LSB_REG 96 << MEM_MAP_SHIFT //0x68 << MEM_MAP_SHIFT
#define SET_DELAY_MSB_REG 97 << MEM_MAP_SHIFT //0x69 << MEM_MAP_SHIFT
#define GET_DELAY_LSB_REG 18 << MEM_MAP_SHIFT//0x6a << MEM_MAP_SHIFT
#define GET_DELAY_MSB_REG 19 << MEM_MAP_SHIFT//0x6b << MEM_MAP_SHIFT
#define SET_CYCLES_LSB_REG 98<<11//0x6c<<11
#define SET_CYCLES_MSB_REG 99<<11//0x6d<<11
#define GET_CYCLES_LSB_REG 20<<11//0x6e<<11
#define GET_CYCLES_MSB_REG 21<<11//0x6f<<11
#define SET_CYCLES_LSB_REG 98 << MEM_MAP_SHIFT//0x6c << MEM_MAP_SHIFT
#define SET_CYCLES_MSB_REG 99 << MEM_MAP_SHIFT//0x6d << MEM_MAP_SHIFT
#define GET_CYCLES_LSB_REG 20 << MEM_MAP_SHIFT//0x6e << MEM_MAP_SHIFT
#define GET_CYCLES_MSB_REG 21 << MEM_MAP_SHIFT//0x6f << MEM_MAP_SHIFT
#define SET_FRAMES_LSB_REG 100<<11//0x70<<11
#define SET_FRAMES_MSB_REG 101<<11//0x71<<11
#define GET_FRAMES_LSB_REG 22<<11//0x72<<11
#define GET_FRAMES_MSB_REG 23<<11//0x73<<11
#define SET_FRAMES_LSB_REG 100 << MEM_MAP_SHIFT//0x70 << MEM_MAP_SHIFT
#define SET_FRAMES_MSB_REG 101 << MEM_MAP_SHIFT//0x71 << MEM_MAP_SHIFT
#define GET_FRAMES_LSB_REG 22 << MEM_MAP_SHIFT//0x72 << MEM_MAP_SHIFT
#define GET_FRAMES_MSB_REG 23 << MEM_MAP_SHIFT//0x73 << MEM_MAP_SHIFT
#define SET_PERIOD_LSB_REG 102<<11//0x74<<11
#define SET_PERIOD_MSB_REG 103<<11//0x75<<11
#define GET_PERIOD_LSB_REG 24<<11//0x76<<11
#define GET_PERIOD_MSB_REG 25<<11//0x77<<11
#define SET_PERIOD_LSB_REG 102 << MEM_MAP_SHIFT//0x74 << MEM_MAP_SHIFT
#define SET_PERIOD_MSB_REG 103 << MEM_MAP_SHIFT//0x75 << MEM_MAP_SHIFT
#define GET_PERIOD_LSB_REG 24 << MEM_MAP_SHIFT//0x76 << MEM_MAP_SHIFT
#define GET_PERIOD_MSB_REG 25 << MEM_MAP_SHIFT//0x77 << MEM_MAP_SHIFT
//#define PATTERN_WAIT0_TIME_REG_LSB 114<<11
//#define PATTERN_WAIT0_TIME_REG_MSB 115<<11
#define SET_EXPTIME_LSB_REG 114<<11//0x78<<11
#define SET_EXPTIME_MSB_REG 115<<11//0x79<<11
#define GET_EXPTIME_LSB_REG 26<<11//0x7a<<11
#define GET_EXPTIME_MSB_REG 27<<11//0x7b<<11
//#define PATTERN_WAIT0_TIME_REG_LSB 114 << MEM_MAP_SHIFT
//#define PATTERN_WAIT0_TIME_REG_MSB 115 << MEM_MAP_SHIFT
#define SET_EXPTIME_LSB_REG 114 << MEM_MAP_SHIFT//0x78 << MEM_MAP_SHIFT
#define SET_EXPTIME_MSB_REG 115 << MEM_MAP_SHIFT//0x79 << MEM_MAP_SHIFT
#define GET_EXPTIME_LSB_REG 26 << MEM_MAP_SHIFT//0x7a << MEM_MAP_SHIFT
#define GET_EXPTIME_MSB_REG 27 << MEM_MAP_SHIFT//0x7b << MEM_MAP_SHIFT
#define SET_GATES_LSB_REG 106<<11//0x7c<<11
#define SET_GATES_MSB_REG 107<<11//0x7d<<11
#define GET_GATES_LSB_REG 28<<11//0x7e<<11
#define GET_GATES_MSB_REG 29<<11//0x7f<<11
#define SET_GATES_LSB_REG 106 << MEM_MAP_SHIFT//0x7c << MEM_MAP_SHIFT
#define SET_GATES_MSB_REG 107 << MEM_MAP_SHIFT//0x7d << MEM_MAP_SHIFT
#define GET_GATES_LSB_REG 28 << MEM_MAP_SHIFT//0x7e << MEM_MAP_SHIFT
#define GET_GATES_MSB_REG 29 << MEM_MAP_SHIFT//0x7f << MEM_MAP_SHIFT
#define DATA_IN_LSB_REG 30<<11
#define DATA_IN_MSB_REG 31<<11
#define DATA_IN_LSB_REG 30 << MEM_MAP_SHIFT
#define DATA_IN_MSB_REG 31 << MEM_MAP_SHIFT
#define PATTERN_OUT_LSB_REG 32<<11
#define PATTERN_OUT_MSB_REG 33<<11
#define PATTERN_OUT_LSB_REG 32 << MEM_MAP_SHIFT
#define PATTERN_OUT_MSB_REG 33 << MEM_MAP_SHIFT
#define FRAMES_FROM_START_LSB_REG 34<<11
#define FRAMES_FROM_START_MSB_REG 35<<11
#define FRAMES_FROM_START_LSB_REG 34 << MEM_MAP_SHIFT
#define FRAMES_FROM_START_MSB_REG 35 << MEM_MAP_SHIFT
#define FRAMES_FROM_START_PG_LSB_REG 36<<11
#define FRAMES_FROM_START_PG_MSB_REG 37<<11
#define FRAMES_FROM_START_PG_LSB_REG 36 << MEM_MAP_SHIFT
#define FRAMES_FROM_START_PG_MSB_REG 37 << MEM_MAP_SHIFT
#define SLOW_ADC_REG 43<<11
#define SLOW_ADC_REG 43 << MEM_MAP_SHIFT
#define PLL_PARAM_REG 80<<11//0x37<<11
#define PLL_PARAM_OUT_REG 5<<11 //0x38<<11
#define PLL_CNTRL_REG 81<<11//0x34<<11
#define PLL_PARAM_REG 80 << MEM_MAP_SHIFT//0x37 << MEM_MAP_SHIFT
#define PLL_PARAM_OUT_REG 5 << MEM_MAP_SHIFT //0x38 << MEM_MAP_SHIFT
#define PLL_CNTRL_REG 81 << MEM_MAP_SHIFT//0x34 << MEM_MAP_SHIFT
#ifdef NEW_GBE_INTERFACE
#define GBE_PARAM_OUT_REG 40<<11
#define GBE_PARAM_REG 69<<11
#define GBE_CNTRL_REG 70<<11
#define GBE_PARAM_OUT_REG 40 << MEM_MAP_SHIFT
#define GBE_PARAM_REG 69 << MEM_MAP_SHIFT
#define GBE_CNTRL_REG 70 << MEM_MAP_SHIFT
#else
#define RX_UDP_AREG 69<<11 //rx_udpip_AReg_c : integer:= 69; *\/
#define UDPPORTS_AREG 70<<11// udpports_AReg_c : integer:= 70; *\/
#define RX_UDPMACL_AREG 71<<11//rx_udpmacL_AReg_c : integer:= 71; *\/
#define RX_UDPMACH_AREG 72<<11//rx_udpmacH_AReg_c : integer:= 72; *\/
#define DETECTORMACL_AREG 73<<11//detectormacL_AReg_c : integer:= 73; *\/
#define DETECTORMACH_AREG 74<<11//detectormacH_AReg_c : integer:= 74; *\/
#define DETECTORIP_AREG 75<<11//detectorip_AReg_c : integer:= 75; *\/
#define IPCHKSUM_AREG 76<<11//ipchksum_AReg_c : integer:= 76; *\/ */
#define RX_UDP_AREG 69 << MEM_MAP_SHIFT //rx_udpip_AReg_c : integer:= 69; *\/
#define UDPPORTS_AREG 70 << MEM_MAP_SHIFT// udpports_AReg_c : integer:= 70; *\/
#define RX_UDPMACL_AREG 71 << MEM_MAP_SHIFT//rx_udpmacL_AReg_c : integer:= 71; *\/
#define RX_UDPMACH_AREG 72 << MEM_MAP_SHIFT//rx_udpmacH_AReg_c : integer:= 72; *\/
#define DETECTORMACL_AREG 73 << MEM_MAP_SHIFT//detectormacL_AReg_c : integer:= 73; *\/
#define DETECTORMACH_AREG 74 << MEM_MAP_SHIFT//detectormacH_AReg_c : integer:= 74; *\/
#define DETECTORIP_AREG 75 << MEM_MAP_SHIFT//detectorip_AReg_c : integer:= 75; *\/
#define IPCHKSUM_AREG 76 << MEM_MAP_SHIFT//ipchksum_AReg_c : integer:= 76; *\/ */
#endif
#define PATTERN_CNTRL_REG 82<<11
#define PATTERN_LIMITS_AREG 83<<11
#define PATTERN_CNTRL_REG 82 << MEM_MAP_SHIFT
#define PATTERN_LIMITS_AREG 83 << MEM_MAP_SHIFT
#define PATTERN_LOOP0_AREG 84<<11
#define PATTERN_N_LOOP0_REG 85<<11
#define PATTERN_LOOP0_AREG 84 << MEM_MAP_SHIFT
#define PATTERN_N_LOOP0_REG 85 << MEM_MAP_SHIFT
#define PATTERN_LOOP1_AREG 86<<11
#define PATTERN_N_LOOP1_REG 87<<11
#define PATTERN_LOOP1_AREG 86 << MEM_MAP_SHIFT
#define PATTERN_N_LOOP1_REG 87 << MEM_MAP_SHIFT
#define PATTERN_LOOP2_AREG 88<<11
#define PATTERN_N_LOOP2_REG 89<<11
#define PATTERN_LOOP2_AREG 88 << MEM_MAP_SHIFT
#define PATTERN_N_LOOP2_REG 89 << MEM_MAP_SHIFT
#define PATTERN_WAIT0_AREG 90<<11
#define PATTERN_WAIT1_AREG 91<<11
#define PATTERN_WAIT2_AREG 92<<11
#define PATTERN_WAIT0_AREG 90 << MEM_MAP_SHIFT
#define PATTERN_WAIT1_AREG 91 << MEM_MAP_SHIFT
#define PATTERN_WAIT2_AREG 92 << MEM_MAP_SHIFT
//#define DAQ_REG 93<<11 //unused
#define NSAMPLES_REG 93<<11
//#define DAQ_REG 93 << MEM_MAP_SHIFT //unused
#define NSAMPLES_REG 93 << MEM_MAP_SHIFT
#define HV_REG 95<<11
#define HV_REG 95 << MEM_MAP_SHIFT
#define PATTERN_IOCTRL_REG_LSB 108<<11
#define PATTERN_IOCTRL_REG_MSB 109<<11
#define PATTERN_IOCTRL_REG_LSB 108 << MEM_MAP_SHIFT
#define PATTERN_IOCTRL_REG_MSB 109 << MEM_MAP_SHIFT
#define PATTERN_IOCLKCTRL_REG_LSB 110<<11
#define PATTERN_IOCLKCTRL_REG_MSB 111<<11
#define PATTERN_IN_REG_LSB 112<<11
#define PATTERN_IN_REG_MSB 113<<11
#define PATTERN_WAIT0_TIME_REG_LSB 114<<11
#define PATTERN_WAIT0_TIME_REG_MSB 115<<11
#define PATTERN_WAIT1_TIME_REG_LSB 116<<11
#define PATTERN_WAIT1_TIME_REG_MSB 117<<11
#define PATTERN_WAIT2_TIME_REG_LSB 118<<11
#define PATTERN_WAIT2_TIME_REG_MSB 119<<11
#define PATTERN_IOCLKCTRL_REG_LSB 110 << MEM_MAP_SHIFT
#define PATTERN_IOCLKCTRL_REG_MSB 111 << MEM_MAP_SHIFT
#define PATTERN_IN_REG_LSB 112 << MEM_MAP_SHIFT
#define PATTERN_IN_REG_MSB 113 << MEM_MAP_SHIFT
#define PATTERN_WAIT0_TIME_REG_LSB 114 << MEM_MAP_SHIFT
#define PATTERN_WAIT0_TIME_REG_MSB 115 << MEM_MAP_SHIFT
#define PATTERN_WAIT1_TIME_REG_LSB 116 << MEM_MAP_SHIFT
#define PATTERN_WAIT1_TIME_REG_MSB 117 << MEM_MAP_SHIFT
#define PATTERN_WAIT2_TIME_REG_LSB 118 << MEM_MAP_SHIFT
#define PATTERN_WAIT2_TIME_REG_MSB 119 << MEM_MAP_SHIFT
//#define DAC_REG_OFF 120
//#define DAC_0_1_VAL_REG 120<<11
//#define DAC_2_3_VAL_REG 121<<11
//#define DAC_4_5_VAL_REG 122<<11
//#define DAC_6_7_VAL_REG 123<<11
//#define DAC_8_9_VAL_REG 124<<11
//#define DAC_10_11_VAL_REG 125<<11
//#define DAC_12_13_VAL_REG 126<<11
//#define DAC_14_15_VAL_REG 127<<11
#define DAC_VAL_REG 121<<11
#define DAC_NUM_REG 122<<11
#define DAC_VAL_OUT_REG 42<<11
#define ADC_LATCH_DISABLE_REG 120<<11
//#define DAC_0_1_VAL_REG 120 << MEM_MAP_SHIFT
//#define DAC_2_3_VAL_REG 121 << MEM_MAP_SHIFT
//#define DAC_4_5_VAL_REG 122 << MEM_MAP_SHIFT
//#define DAC_6_7_VAL_REG 123 << MEM_MAP_SHIFT
//#define DAC_8_9_VAL_REG 124 << MEM_MAP_SHIFT
//#define DAC_10_11_VAL_REG 125 << MEM_MAP_SHIFT
//#define DAC_12_13_VAL_REG 126 << MEM_MAP_SHIFT
//#define DAC_14_15_VAL_REG 127 << MEM_MAP_SHIFT
#define DAC_VAL_REG 121 << MEM_MAP_SHIFT
#define DAC_NUM_REG 122 << MEM_MAP_SHIFT
#define DAC_VAL_OUT_REG 42 << MEM_MAP_SHIFT
#define ADC_LATCH_DISABLE_REG 120 << MEM_MAP_SHIFT
@ -241,27 +241,27 @@
/* registers defined in FPGA */
#define GAIN_REG 0
//#define FLOW_CONTROL_REG 0x11<<11
//#define FLOW_STATUS_REG 0x12<<11
//#define FRAME_REG 0x13<<11
//#define FLOW_CONTROL_REG 0x11 << MEM_MAP_SHIFT
//#define FLOW_STATUS_REG 0x12 << MEM_MAP_SHIFT
//#define FRAME_REG 0x13 << MEM_MAP_SHIFT
#define MULTI_PURPOSE_REG 0
//#define TIME_FROM_START_REG 0x16<<11
//#define TIME_FROM_START_REG 0x16 << MEM_MAP_SHIFT
#define ROI_REG 0 // 0x35<<11
#define OVERSAMPLING_REG 0 // 0x36<<11
#define MOENCH_CNTR_REG 0 // 0x31<<11
#define MOENCH_CNTR_OUT_REG 0 // 0x33<<11
#define MOENCH_CNTR_CONF_REG 0 // 0x32<<11
#define ROI_REG 0 // 0x35 << MEM_MAP_SHIFT
#define OVERSAMPLING_REG 0 // 0x36 << MEM_MAP_SHIFT
#define MOENCH_CNTR_REG 0 // 0x31 << MEM_MAP_SHIFT
#define MOENCH_CNTR_OUT_REG 0 // 0x33 << MEM_MAP_SHIFT
#define MOENCH_CNTR_CONF_REG 0 // 0x32 << MEM_MAP_SHIFT
//image
#define DARK_IMAGE_REG 0 // 0x81<<11
#define GAIN_IMAGE_REG 0 // 0x82<<11
#define DARK_IMAGE_REG 0 // 0x81 << MEM_MAP_SHIFT
#define GAIN_IMAGE_REG 0 // 0x82 << MEM_MAP_SHIFT
//counter block memory
#define COUNTER_MEMORY_REG 0 // 0x85<<11
#define COUNTER_MEMORY_REG 0 // 0x85 << MEM_MAP_SHIFT
//not used

7
slsDetectorSoftware/jctbDetectorServer/server_funcs.c Executable file → Normal file
View File

@ -7,8 +7,9 @@
#include "slow_adc.h"
#include "registers_m.h"
#include "gitInfoMoench.h"
#include "blackfin.h"
#define FIFO_DATA_REG_OFF 0x50<<11
#define FIFO_DATA_REG_OFF 0x50 << MEM_MAP_SHIFT
// Global variables
@ -866,7 +867,7 @@ int write_register(int file_des) {
if(ret!=FAIL){
address=(addr<<11);
address=(addr << MEM_MAP_SHIFT);
if((address==FIFO_DATA_REG_OFF)||(address==CONTROL_REG))
ret = bus_w16(address,val);
else
@ -932,7 +933,7 @@ int read_register(int file_des) {
//#endif
if(ret!=FAIL){
address=(addr<<11);
address=(addr << MEM_MAP_SHIFT);
if((address==FIFO_DATA_REG_OFF)||(address==CONTROL_REG))
retval=bus_r16(address);
else

View File

@ -1,6 +1,7 @@
#include "firmware_funcs.h"
#include "registers_m.h"
#include "server_defs.h"
#include "blackfin.h"
int prepareSlowADCSeq() {

View File

@ -13,7 +13,7 @@ INSTMODE= 0777
SRCS= server.c server_funcs.c communication_funcs.c firmware_funcs.c mcb_funcs.c trimming_funcs.c sharedmemory.c
OBJS= $(SRCS:%.c=%.o)
VFLAGS=
VFLAGS=
#-DVERBOSE
#-DVERYVERBOSE
CFLAGS+= -Wall -DC_ONLY -DMCB_FUNCS -DDACS_INT $(VFLAGS)
@ -32,7 +32,7 @@ versioning:
$(PROGS): $(OBJS)
echo $(OBJS)
$(CC) $(LDFLAGS) $^ $(LDLIBS) $(CFLAGS) -o $@
$(CC) $(LDFLAGS) $^ $(LDLIBS) $(CFLAGS) -o $@ -DVERBOSE
install: $(PROGS)
$(INSTALL) -d $(INSTDIR)

View File

@ -17,6 +17,7 @@
#include <sys/stat.h>
#include <stdlib.h>
//#define VERBOSE
//for memory mapping
u_int64_t CSP0BASE;
@ -702,9 +703,9 @@ int setNMod(int n) {
shiftfifo=SHIFTFIFO;
#ifdef VERBOSE
//#ifdef VERBOSE
printf("SetNMod called arg %d -- dr %d shiftfifo %d\n",n,dynamicRange,shiftfifo);
#endif
//#endif
if (n>=0 && n<=ntot) {
nModX=n;
@ -713,12 +714,16 @@ int setNMod(int n) {
reg=bus_r(FIFO_COUNTR_REG_OFF+(ififo<<shiftfifo));
if (ififo<n*NCHIP) {
if (reg&FIFO_DISABLED_BIT) {
bus_w(FIFO_CNTRL_REG_OFF+(ififo<<shiftfifo), FIFO_DISABLE_TOGGLE_BIT);
#ifdef VERBOSE
if (bus_r(FIFO_COUNTR_REG_OFF+(ififo<<shiftfifo))&FIFO_DISABLED_BIT) {
// if (bus_r(FIFO_COUNTR_REG_OFF+(ififo<<shiftfifo))&FIFO_DISABLED_BIT) {
printf("Fifo %d is %x (nm %d nc %d addr %08x)",ififo,reg, (reg&FIFO_NM_MASK)>>FIFO_NM_OFF, (reg&FIFO_NC_MASK)>>FIFO_NC_OFF, FIFO_COUNTR_REG_OFF+(ififo<<shiftfifo));
printf(" enabling %08x\n",bus_r(FIFO_COUNTR_REG_OFF+(ififo<<shiftfifo)));
}
// }
#endif
}
//#ifdef VERBOSE
@ -729,10 +734,10 @@ int setNMod(int n) {
if ((reg&FIFO_ENABLED_BIT)) {
bus_w(FIFO_CNTRL_REG_OFF+(ififo<<shiftfifo), FIFO_DISABLE_TOGGLE_BIT);
#ifdef VERBOSE
if ((bus_r(FIFO_COUNTR_REG_OFF+(ififo<<shiftfifo))&FIFO_ENABLED_BIT)) {
// if ((bus_r(FIFO_COUNTR_REG_OFF+(ififo<<shiftfifo))&FIFO_ENABLED_BIT)) {
printf("Fifo %d is %x (nm %d nc %d addr %08x)",ififo,reg, (reg&FIFO_NM_MASK)>>FIFO_NM_OFF, (reg&FIFO_NC_MASK)>>FIFO_NC_OFF, FIFO_COUNTR_REG_OFF+(ififo<<shiftfifo));
printf(" disabling %08x\n",bus_r(FIFO_COUNTR_REG_OFF+(ififo<<shiftfifo)));
}
// }
#endif
}
//#ifdef VERBOSE

View File

@ -1,9 +1,9 @@
Path: slsDetectorsPackage/slsDetectorSoftware/mythenDetectorServer
URL: origin git@git.psi.ch:sls_detectors_software/sls_detector_software.git/mythenDetectorServer
Repository Root: origin git@git.psi.ch:sls_detectors_software/sls_detector_software.git
Repsitory UUID: 8aceb5d4b0ca6bd95a11b53e7a799b463b92d51b
Revision: 106
Branch: developer
Last Changed Author: Dhanya_Maliakal
Last Changed Rev: 334
Last Changed Date: 2016-08-12 11:08:03 +0200
URL:
Repository Root:
Repsitory UUID:
Revision: 0
Branch:
Last Changed Author: _
Last Changed Rev: 0
Last Changed Date: 2018-01-30 10:42:45.000000000 +0100 ./Makefile

View File

@ -1,6 +1,6 @@
#define GITURL "git@git.psi.ch:sls_detectors_software/sls_detector_software.git/mythenDetectorServer"
#define GITREPUUID "8aceb5d4b0ca6bd95a11b53e7a799b463b92d51b"
#define GITAUTH "Dhanya_Maliakal"
#define GITREV 0x334
#define GITDATE 0x20160812
#define GITBRANCH "blabla"
#define GITURL "URL:"
#define GITREPUUID "UUID:"
#define GITAUTH "_"
#define GITREV 0x0
#define GITDATE 0x20180130
#define GITBRANCH "Branch:"

View File

@ -1,6 +1,8 @@
#define GITURL "git@github.com:slsdetectorgroup/slsDetectorPackage.git"
#define GITREPUUID "b8bdbf4da61f95b88893b02ddabc2491b16fa10f"
#define GITAUTH "Dhanya_Thattil"
#define GITREV 0x3746
#define GITDATE 0x20180327
#define GITBRANCH "developer"

View File

@ -1,9 +1,18 @@
Path: slsDetectorsPackage/slsReceiverSoftware
URL: origin git@github.com:slsdetectorgroup/slsDetectorPackage.git
Repository Root: origin git@github.com:slsdetectorgroup/slsDetectorPackage.git
<<<<<<< HEAD
Repsitory UUID: ab06c33107ecfeb4741d49407903ff80286cf75b
Revision: 765
Branch: developer
Last Changed Author: Anna_Bergamaschi
Last Changed Rev: 3731
Last Changed Date: 2018-03-15 12:27:06.000000002 +0100 ./src/slsReceiverTCPIPInterface.cpp
=======
Repsitory UUID: b8bdbf4da61f95b88893b02ddabc2491b16fa10f
Revision: 767
Branch: developer
Last Changed Author: Dhanya_Thattil
Last Changed Rev: 3746
Last Changed Date: 2018-03-27 10:43:44.000000002 +0200 ./src/slsReceiverTCPIPInterface.cpp
>>>>>>> 7cd35f24b87501374fbaf45693a2adf16dfae3e3

View File

@ -1,4 +1,5 @@
#define GITURL "git@github.com:slsdetectorgroup/slsDetectorPackage.git"
#define GITREPUUID "b8bdbf4da61f95b88893b02ddabc2491b16fa10f"
#define GITAUTH "Dhanya_Thattil"
#define GITREV 0x3746