detectors/slsDetectorPackage
25
0
Fork 0
mirror of https://github.com/slsdetectorgroup/slsDetectorPackage.git synced 2025-04-19 18:40:01 +02:00
Code Issues Actions Packages Releases Activity
slsDetectorPackage/slsDetectorCalibration/analogDetector.h
Erik Fröjdh 61f290441e
Local loop indicies in analogDetector (#709) 
* making ix and iy local variables in analogDetector
2023-05-03 09:18:09 +02:00

1390 lines
44 KiB
C++
Raw Blame History

// SPDX-License-Identifier: LGPL-3.0-or-other
// Copyright (C) 2021 Contributors to the SLS Detector Package
#ifndef ANALOGDETECTOR_H
#define ANALOGDETECTOR_H
//#include <mutex>
#include "commonModeSubtractionNew.h"
#include "ghostSummation.h"
#include "pedestalSubtraction.h"
#include "slsDetectorData.h"
#include "slsInterpolation.h"
#include "sls/tiffIO.h"
#include <pthread.h>
#ifdef ROOTSPECTRUM
#include <TASImage.h>
#include <TF1.h>
#include <TFile.h>
#include <TGraphErrors.h>
#include <TH2F.h>
#include <TImage.h>
#include <TLegend.h>
#include <TPaveText.h>
#endif
using namespace std;
#ifndef FRAMEMODE_DEF
#define FRAMEMODE_DEF
/**
enum to define the flags of the data set, which are needed to seect the type of
processing it should undergo: frame, pedestal, flat
*/
enum frameMode { eFrame, ePedestal, eFlat, eRaw };
/**
enum to define the detector mode
*/
enum detectorMode { eAnalog, ePhotonCounting, eInterpolating };
#endif
template <class dataType> class analogDetector {
/** @short class to perform pedestal subtraction etc. for an analog detector
*/
public:
/**
Constructor (no error checking if datasize and offsets are compatible!)
\param d detector data structure to be used - if null it is assumed that
the data are in ordered ip=iy*nx+ix \param sign is the sign of the data
\param nped number of samples for pedestal averaging
\param cm common mode subtraction algorithm, if any. Defaults to NULL
i.e. none \param nnx detector size in x - must be specified if no data
structure is defined, otherwise defaults to the size of the data
structure. \param nny detector size in y - must be specified if no data
structure is defined, otherwise defaults to the size of the data
structure. \param gm pointer to tha gain map matrix
*/
analogDetector(slsDetectorData<dataType> *d, int sign = 1,
commonModeSubtraction *cm = NULL, int nped = 1000,
int nnx = -1, int nny = -1, double *gm = NULL,
ghostSummation<dataType> *gs = NULL)
: det(d), nx(nnx), ny(nny), stat(NULL), cmSub(cm), dataSign(sign),
iframe(-1), gmap(gm), ghSum(gs), id(0) {
if (det)
det->getDetectorSize(nx, ny);
stat = new pedestalSubtraction *[ny];
/* pedMean=new double*[ny]; */
/* pedVariance=new double*[ny]; */
for (int i = 0; i < ny; i++) {
stat[i] = new pedestalSubtraction[nx];
/* pedMean[i]=new double[nx]; */
/* pedVariance[i]=new double[nx]; */
for (int ix = 0; ix < nx; ++ix) {
stat[i][ix].SetNPedestals(nped);
/* stat[i][ix].setPointers(&(pedMean[iy][ix]),&(pedVariance[iy][ix]));
*/
}
}
image = new int[nx * ny];
xmin = 0;
xmax = nx;
ymin = 0;
ymax = ny;
fMode = ePedestal;
thr = 0;
myFile = NULL;
#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 [] pedMean[i]; */
/* delete [] pedVariance[i]; */
};
/* delete [] pedMean; */
/* delete [] pedVariance; */
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
\param orig analog Detector structure to be cloned
*/
analogDetector(analogDetector *orig) {
/* copy construction from orig*/
det = orig->det;
nx = orig->nx;
ny = orig->ny;
dataSign = orig->dataSign;
iframe = orig->iframe;
gmap = orig->gmap;
// cmSub=orig->cmSub;
id = orig->id;
xmin = orig->xmin;
xmax = orig->xmax;
ymin = orig->ymin;
ymax = orig->ymax;
thr = orig->thr;
// nSigma=orig->nSigma;
fMode = orig->fMode;
myFile = orig->myFile;
stat = new pedestalSubtraction *[ny];
/* pedMean=new double*[ny]; */
/* pedVariance=new double*[ny]; */
for (int i = 0; i < ny; i++) {
stat[i] = new pedestalSubtraction[nx];
/* pedMean[i]=new double[nx]; */
/* pedVariance[i]=new double[nx]; */
}
int nped = orig->SetNPedestals();
// cout << nped << " " << orig->getPedestal(ix,iy) <<
// orig->getPedestalRMS(ix,iy) << endl;
for (int iy = 0; iy < ny; ++iy) {
for (int ix = 0; ix < nx; ++ix) {
// stat[iy][ix].setPointers(&(pedMean[iy][ix]),&(pedVariance[iy][ix]));
stat[iy][ix].SetNPedestals(nped);
setPedestal(ix, iy, orig->getPedestal(ix, iy),
orig->getPedestalRMS(ix, iy),
orig->GetNPedestals(ix, iy));
/* if (ix==50 && iy==50) */
/* cout << "clone ped " << " " << ix << " " << iy << " " <<
* getPedestal(ix,iy) << " " << getPedestalRMS(ix,iy)<< " " <<
* GetNPedestals(ix,iy) << endl; */
}
}
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
if (orig->cmSub) {
cmSub = (orig->cmSub)->Clone();
cout << "cloning cm" << endl;
} else
cmSub = NULL;
if (orig->ghSum) {
ghSum = (orig->ghSum)->Clone();
cout << "cloning gs" << endl;
} else
ghSum = NULL;
}
/**
clone. Must be virtual!
\returns a clone of the original analog detector
*/
virtual analogDetector *Clone() = 0;
/**
Gives an id to the structure. For debugging purposes in case of
multithreading. \param i is to be set \returns current id
*/
int setId(int i) {
id = i;
return id;
};
/**
Returns id of the structure. For debugging purposes in case of
multithreading. \returns current id
*/
int getId() { return id; };
/**
Returns data size of the detector data structure
\returns data size of the detector data structurein bytes
*/
int getDataSize() { return det->getDataSize(); };
/**
Returns data size of the detector image matrix
\param nnx reference to image size in x
\param nny reference to image size in y
\param nns reference to number of subpixels for interpolating detector,
will always be 1 in this case \returns number of pixels of the detector
image
*/
virtual int getImageSize(int &nnx, int &nny, int &nnsx, int &nnsy) {
nnx = nx;
nny = ny;
nnsx = 1;
nnsy = 1;
return nx * ny;
};
/**
Returns data size of the detector image matrix
\param nnx reference to pixel size in x
\param nny reference to pixel size in y
\returns number of pixels of the detector image
*/
virtual int getDetectorSize(int &nnx, int &nny) {
nnx = nx;
nny = ny;
return nx * ny;
};
/**
set gain map
\param gm pointer to gain map matrix to be set - NULL unsets
\returns pointer to current gain map
*/
double *setGainMap(double *gm) {
gmap = gm;
return gmap;
};
/**
return gain map
\returns pointer to current gain map
*/
double *getGainMap() { return gmap; };
/**
reads a 32 bit tiff file of the size of the detector and sets its values
as gain map for the detector. If file does not exist returns NULL, but
does not change gainmap compared to previous settings. \param imgname
complete name of the file containing the gain map data \returns pointer
to current gain map is file reading succeeded, NULL is file reading
didn't work.
*/
double *readGainMap(const char *imgname) {
uint32_t nnx, nny;
float *gm = ReadFromTiff(imgname, nny, nnx);
if (gm) {
if (gmap)
delete[] gmap;
gmap = new double[nnx * nny];
for (int iy = 0; iy < static_cast<int>(nny); ++iy) {
for (int ix = 0; ix < static_cast<int>(nnx); ++ix) {
gmap[iy * nnx + ix] = gm[iy * nnx + ix];
// cout << gmap[iy*nnx+ix] << " " ;
}
}
return gmap;
}
return NULL;
}
/**
writes a 32 bit tiff file of the size of the detector and contaning the
gain map value, if any. If file doesn'e exist or gainmap is undefined,
does not do anything. \param imgname complete name of the file to be
written \returns NULL if file writing didn't succeed, else a pointer
*/
void *writeGainMap(const char *imgname) {
float *gm = NULL;
void *ret;
if (gmap) {
gm = new float[nx * ny];
for (int iy = 0; iy < ny; ++iy) {
for (int ix = 0; ix < nx; ++ix) {
gm[iy * nx + ix] = gmap[iy * nx + ix];
}
}
ret = WriteToTiff(gm, imgname, ny, nx);
delete[] gm;
}
return ret;
}
/** resets the pedestalSubtraction array, the commonModeSubtraction and the
* image data*/
virtual void newDataSet() {
iframe = -1;
for (int iy = 0; iy < ny; ++iy)
for (int ix = 0; ix < nx; ++ix) {
stat[iy][ix].Clear();
image[iy * nx + ix] = 0;
}
if (cmSub)
cmSub->Clear();
#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 */
virtual void newFrame() {
iframe++;
if (cmSub)
cmSub->newFrame();
det->newFrame();
};
/** resets the commonModeSubtraction and increases the frame index */
virtual void newFrame(char *data) {
iframe++;
if (cmSub)
cmSub->newFrame();
det->newFrame();
// det->getData(data);
calcGhost(data);
// cout << getId() << " Calc ghost " << getGhost(15,15) << endl;
};
/** sets the commonModeSubtraction algorithm to be used
\param cm commonModeSubtraction algorithm to be used (NULL unsets)
\returns pointer to the actual common mode subtraction algorithm
*/
commonModeSubtraction *setCommonModeSubtraction(commonModeSubtraction *cm) {
cmSub = cm;
return cmSub;
};
/**
gets the commonModeSubtraction algorithm to be used
\returns pointer to the actual common mode subtraction algorithm
*/
commonModeSubtraction *getCommonModeSubtraction() { return cmSub; };
ghostSummation<dataType> *getGhostSummation() { return ghSum; };
ghostSummation<dataType> *setGhostSummation(ghostSummation<dataType> *gs) {
ghSum = gs;
return ghSum;
};
/**
sets the sign of the data
\param sign 1 means positive values for photons, -1 negative, 0 gets
\returns current sign for the data
*/
int setDataSign(int sign = 0) {
if (sign == 1 || sign == -1)
dataSign = sign;
return dataSign;
};
/**
adds value to pedestal (and common mode) for the given pixel
\param val value to be added
\param ix pixel x coordinate
\param iy pixel y coordinate
\param cm 1 adds the value to common mod, 0 skips it. Defaults to 0. -
not properly implemented
*/
virtual void addToPedestal(double val, int ix, int iy, int cm = 0) {
if (ix >= 0 && ix < nx && iy >= 0 && iy < ny) {
// cout << val << " " ;
if (cmSub && cm > 0) {
val -= getCommonMode(ix, iy);
}
// cout << val << " " ;
val += getGhost(ix, iy);
// cout << val ;
// cout << endl;
stat[iy][ix].addToPedestal(val);
/* 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);
}
return 0;
}
virtual void addToCommonMode(char *data) {
// cout << "+"<< getId() << endl;
if (cmSub) {
// cout << "*" << endl;
for (int iy = ymin; iy < ymax; ++iy) {
for (int ix = xmin; ix < xmax; ++ix) {
// if (getNumpedestals(ix,iy)>0)
// if (det->isGood(ix,iy)) {
addToCommonMode(data, ix, iy);
// cout << ":";
// }
}
}
// cout << "cm " << getCommonMode(0,0) << " " << getCommonMode(1,0)
// << endl;
}
}
virtual void addToCommonMode(char *data, int ix, int iy = 0) {
// cout <<".";
if (cmSub) {
// cout <<":";
if (det)
if (det->isGood(ix, iy) == 0)
return;
if (getNumpedestals(ix, iy) > 0) {
double val;
if (det) {
/* if (det->getChannel(data, ix, iy)>=0x3fff) */
/* cout << ix << " " << iy << " " << det->getChannel(data,
* ix, iy) <<endl; */
val = (dataSign * det->getValue(data, ix, iy) -
getPedestal(ix, iy, 0));
} else
val =
(((double *)data)[iy * nx + ix] - getPedestal(ix, iy));
val += getGhost(ix, iy);
cmSub->addToCommonMode(val, ix, iy);
// cout << ":";
}
}
}
/**
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() + getCommonMode(ix, iy);
// return pedMean[iy][ix]+getCommonMode(ix,iy);
}
// return pedMean[iy][ix];
return stat[iy][ix].getPedestal();
} else
return -1;
};
/**
gets pedestal rms (i.e. noise)
\param ix pixel x coordinate
\param iy pixel y coordinate
\returns pedestal rms
*/
virtual double getPedestalRMS(int ix, int iy) {
double g = 1;
if (ix >= 0 && ix < nx && iy >= 0 && iy < ny) {
if (gmap) {
g = gmap[iy * nx + ix];
if (g == 0)
g = -1.;
}
// return sqrt(pedVariance[iy][ix])/g;
return stat[iy][ix].getPedestalRMS() / g; // divide by gain?
}
return -1;
};
virtual int getNumpedestals(int ix, int iy) {
if (ix >= 0 && ix < nx && iy >= 0 && iy < ny)
return stat[iy][ix].getNumpedestals();
return -1;
};
/**
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(double *ped) {
if (ped == NULL) {
ped = new double[nx * ny];
}
for (int iy = 0; iy < ny; ++iy) {
for (int ix = 0; ix < nx; ++ix) {
ped[iy * nx + ix] = stat[iy][ix].getPedestal();
// cout << ped[iy*nx+ix] << " " ;
}
}
return ped;
};
/**
gets pedestal rms (i.e. noise)
\param ix pixel x coordinate
\param iy pixel y coordinate
\returns pedestal rms
*/
virtual double *getPedestalRMS(double *ped = NULL) {
if (ped == NULL) {
ped = new double[nx * ny];
}
for (int iy = 0; iy < ny; ++iy) {
for (int ix = 0; ix < nx; ++ix) {
ped[iy * nx + ix] = stat[iy][ix].getPedestalRMS();
}
}
return ped;
};
/**
sets pedestal
\param ix pixel x coordinate
\param iy pixel y coordinate
\param val value to set
\param rms rms to be set if any, defaults to 0
\param m number of pedestal samples to be set or the moving stat
structure is any, defaults to 0
*/
virtual void setPedestal(int ix, int iy, double val, double rms = 0,
int m = -1) {
if (ix >= 0 && ix < nx && iy >= 0 && iy < ny)
stat[iy][ix].setPedestal(val, rms, m);
};
/**
sets pedestal
\param ix pixel x coordinate
\param iy pixel y coordinate
\param val value to set
\param rms rms to be set if any, defaults to 0
\param m number of pedestal samples to be set or the moving stat
structure is any, defaults to 0
*/
virtual void setPedestal(double *ped, double *rms = NULL, int m = -1) {
double rr = 0;
for (int iy = ymin; iy < ymax; ++iy) {
for (int ix = xmin; ix < xmax; ++ix) {
if (rms)
rr = rms[iy * nx + ix];
stat[iy][ix].setPedestal(ped[iy * nx + ix], rr, m);
// cout << ix << " " << iy << " " << ped[iy*nx+ix] << " " <<
// stat[iy][ix].getPedestal() << endl;
};
};
}
/**
sets pedestal rms
\param ix pixel x coordinate
\param iy pixel y coordinate
\param rms value to set
*/
virtual void setPedestalRMS(int ix, int iy, double rms = 0) {
if (ix >= 0 && ix < nx && iy >= 0 && iy < ny)
stat[iy][ix].setPedestalRMS(rms);
};
/**
sets pedestal rms for all pixels
\param rms pointer to array of pedestal rms
*/
virtual void setPedestalRMS(double *rms) {
for (int iy = ymin; iy < ymax; ++iy) {
for (int ix = xmin; ix < xmax; ++ix) {
stat[iy][ix].setPedestalRMS(rms[iy * nx + ix]);
};
};
}
virtual void calcGhost(char *data, int ix, int iy = 1) {
if (ghSum) {
ghSum->calcGhost(data, ix, iy);
}
};
virtual void calcGhost(char *data) {
if (ghSum) {
ghSum->calcGhost(data);
// cout << getId() << " @ " << ghSum->getXTalk() << " " <<
// ghSum->getGhost(15,15) << endl;
}
};
virtual double getGhost(int ix, int iy) {
if (ghSum)
return ghSum->getGhost(ix, iy);
return 0;
};
/**
write 32bit tiff file with detector image data
\param imgname file name to be written
\returns NULL if file writing didn't succed, otherwise a pointer
*/
virtual void *writeImage(const char *imgname) {
float *gm = NULL;
void *ret;
#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 iy = 0; iy < ny; ++iy) {
for (int ix = 0; ix < nx; ++ix) {
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
\returns NULL if file writing didn't succed, otherwise a pointer
*/
virtual void *writePedestals(const char *imgname) {
float *gm = NULL;
void *ret;
gm = new float[nx * ny];
#ifdef ROOTSPECTRUM
TH2F *hmap =
new TH2F("hmap", "hmap", nx, -0.5, nx - 0.5, ny, -0.5, ny - 0.5);
#endif
for (int iy = 0; iy < ny; ++iy) {
for (int ix = 0; ix < nx; ++ix) {
/* 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;
#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;
}
/**
read 32bit tiff file containing the pedestals
\param imgname file name to be read
\returns 0 if file reading didn't succed, otherwise 1
*/
int readPedestals(const char *imgname) {
uint32_t nnx, nny;
float *gm = ReadFromTiff(imgname, nny, nnx);
if (nnx > nx)
nnx = nx;
if (nny > ny)
nny = ny;
if (gm) {
for (int iy = 0; iy < nny; ++iy) {
for (int ix = 0; ix < nnx; ++ix) {
stat[iy][ix].setPedestal(gm[iy * nx + ix], -1, -1);
}
}
delete[] gm;
return 1;
}
return 0;
}
/**
read 32bit tiff file containing the image data
\param imgname file name to be read
\returns 0 if file reading didn't succed, otherwise 1
*/
int readImage(const char *imgname) {
uint32_t nnx, nny;
float *gm = ReadFromTiff(imgname, nny, nnx);
if (nnx > nx)
nnx = nx;
if (nny > ny)
nny = ny;
if (gm) {
for (int iy = 0; iy < nny; ++iy) {
for (int ix = 0; ix < nnx; ++ix) {
image[iy * nx + ix] = gm[iy * nx + ix];
}
}
delete[] gm;
return 1;
}
return 0;
}
/**
returns pointer to image data
\returns pointer to image data
*/
virtual int *getImage() { return image; };
/**
write 32bit tiff file containing the pedestals RMS
\param imgname file name to be written
\returns NULL if file writing didn't succed, otherwise a pointer
*/
void *writePedestalRMS(const char *imgname) {
float *gm = NULL;
void *ret;
gm = new float[nx * ny];
for (int iy = 0; iy < ny; ++iy) {
for (int ix = 0; ix < nx; ++ix) {
gm[iy * nx + ix] = stat[iy][ix].getPedestalRMS();
}
}
ret = WriteToTiff(gm, imgname, ny, nx);
delete[] gm;
return ret;
}
/**
read 32bit tiff file containing the pedestals RMS
\param imgname file name to be read
\returns 0 if file reading didn't succed, otherwise 1
*/
int readPedestalRMS(const char *imgname) {
uint32_t nnx, nny;
float *gm = ReadFromTiff(imgname, nny, nnx);
if (nnx > nx)
nnx = nx;
if (nny > ny)
nny = ny;
if (gm) {
for (uint32_t iy = 0; iy < nny; ++iy) {
for (uint32_t ix = 0; ix < nnx; ++ix) {
stat[iy][ix].setPedestalRMS(gm[iy * nx + ix]);
}
}
delete[] gm;
return 1;
}
return 0;
}
/**
Adds all the data for each pixels in the selected region of interest to
the pedestal \param data pointer to the data
*/
virtual void addToPedestal(char *data, int cm = 0) {
// cout << "add to pedestal " << endl;
newFrame(data);
// calcGhost(data);
if (cmSub && cm) {
// cout <<",";
addToCommonMode(data);
}
// cout << xmin << " " << xmax << endl;
// cout << ymin << " " << ymax << endl;
for (int iy = ymin; iy < ymax; ++iy) {
for (int ix = xmin; ix < xmax; ++ix) {
if (det->isGood(ix, iy)) {
// addToPedestal(data,ix,iy,1);
addToPedestal(data, ix, iy, cm);
/* if (ix==50 && iy==50) */
/* cout<< "*ped* " << id << " " << ix << " " << iy << " "
* << det->getChannel(data,ix,iy) << " " <<
* stat[iy][ix].getPedestal() << " " <<
* stat[iy][ix].getPedestalRMS() <<" " << stat[iy][ix].
* GetNPedestals() << endl; */
// if (ix==10 && iy==10)
// cout <<ix << " " << iy << " " << getPedestal(ix,iy)<<
// endl;
#ifdef ROOTSPECTRUM
subtractPedestal(data, ix, iy, cm);
#endif
}
}
}
return;
};
/**
Sets region of interest in which data should be processed
\param xmi minimum x. if -1 or out of range remains unchanged
\param xma maximum x. if -1 or out of range remains unchanged
\param ymi minimum y. if -1 or out of range remains unchanged
\param yma maximum y. if -1 or out of range remains unchanged
*/
void setROI(int xmi = -1, int xma = -1, int ymi = -1, int yma = -1) {
if (xmi >= 0 && xmi <= nx)
xmin = xmi;
if (xma >= 0 && xma <= nx)
xmax = xma;
if (xmax < xmin) {
xmi = xmin;
xmin = xmax;
xmax = xmi;
}
if (ymi >= 0 && ymi <= ny)
ymin = ymi;
if (yma >= 0 && yma <= ny)
ymax = yma;
if (ymax < ymin) {
ymi = ymin;
ymin = ymax;
ymax = ymi;
}
#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
};
/**
Gets region of interest in which data are processed
\param xmi reference to minimum x.
\param xma reference to maximum x.
\param ymi reference to minimum y.
\param yma reference to maximum y.
*/
void getROI(int &xmi, int &xma, int &ymi, int &yma) {
xmi = xmin;
xma = xmax;
ymi = ymin;
yma = ymax;
};
/**
Adds all the data for the selected pixel to the pedestal
\param data pointer to the data
\param ix pixel x coordinate
\param iy pixel y coordinate
*/
virtual void addToPedestal(char *data, int ix, int iy, int cm = 0) {
double val;
if (ix >= 0 && ix < nx && iy >= 0 && iy < ny) {
if (det)
val = dataSign * det->getValue(data, ix, iy);
else
val = ((double *)data)[iy * nx + ix];
// cout << val << endl;
/* if (ix==10 && iy==10) */
/* cout << ix << " " << iy << " " << val ; */
/* if (ix==100 && iy==100) */
/* cout << ix << " " << iy << " " << val; */
addToPedestal(val, ix, iy);
// cout << val << endl;
/* if (ix==10 && iy==10) */
/* cout <<" " << getPedestal(ix,iy)<< endl; */
/* if (ix==100 && iy==100) */
/* cout << " " << getPedestal(ix,iy)<< endl; */
}
return;
};
/**
Subtracts pedestal from the data array in the region of interest
\param data pointer to the data
\param val pointer where the pedestal subtracted data should be added.
If NULL, the internal image is used \returns pointer to the pedestal
subtracted data
*/
// virtual int *subtractPedestal(char *data, int *val=NULL) {
virtual int *subtractPedestal(char *data, int *val = NULL, int cm = 0) {
newFrame(data);
if (val == NULL)
val = image; // new double[nx*ny];
// calcGhost(data);
for (int iy = ymin; iy < ymax; ++iy) {
for (int ix = xmin; ix < xmax; ++ix) {
if (det->isGood(ix, iy))
val[iy * nx + ix] += subtractPedestal(data, ix, iy, cm);
}
}
return val;
};
/**
Subtracts pedestal from the data for a selected pixel
\param data pointer to the data
\param ix pixel x coordinate
\param iy pixel y coordinate
\returns pedestal subtracted value
*/
virtual double subtractPedestal(char *data, int ix, int iy = 0,
int cm = 0) {
double g = 1.;
double val = 0;
if (ix >= 0 && ix < nx && iy >= 0 && iy < ny) {
if (gmap) {
g = gmap[iy * nx + ix];
if (g == 0)
g = -1.;
}
if (det) {
/* if (det->getChannel(data, ix, iy)>=0x3fff) */
/* cout << ix << " " << iy << " " << det->getChannel(data, ix,
* iy) <<endl; */
val = (dataSign * det->getValue(data, ix, iy) -
getPedestal(ix, iy, cm)) /
g;
} else
val =
(((double *)data)[iy * nx + ix] - getPedestal(ix, iy)) / g;
// if (val>=0.5*thr)
// cout << val << " " << (dataSign*det->getValue(data, ix,
// iy)-getPedestal(ix,iy,cm)) << " " << g << " ";
val += getGhost(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; i++ ix)
for (int iiy = -4; iiy < 5; i++ iy) {
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;
}
return val;
};
/**
sets threshold value for conversion into number of photons
\param t threshold to be set
\returns threshold value
*/
double setThreshold(double t) {
thr = t;
return thr;
};
virtual int setClusterSize(int n = -1) = 0;
/**
gets threshold value for conversion into number of photons
\returns threshold value
*/
double getThreshold() { return thr; };
/**
converts the data into number of photons for the selected pixel
\param data pointer to the data
\param ix pixel x coordinate
\param iy pixel y coordinate
\returns converted number of photons. If no threshold is set, returns gain
converted pedestal subtracted data.
*/
int convertToPhotons(char *data, int ix, int iy = 0) {
int nph = 0;
double v;
if (ix >= 0 && ix < nx && iy >= 0 && iy < ny) {
v = subtractPedestal(data, ix, iy, 1);
/* // 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;
/* if (ix ==10 && iy<20) */
/* cout << det->getValue(data,ix,iy) << " " <<
* stat[iy][ix].getPedestal() << " " << getCommonMode(ix,iy) <<
* " " << getPedestal(ix,iy,0)<< " " << getPedestal(ix,iy,1) <<
* endl; */
// cout << subtractPedestal(data,ix,iy,0) << " " <<
// subtractPedestal(data,ix,iy,1) << " " << nph << endl;
if (nph > 0) {
// cout << " " << nph << endl;
return nph;
}
return 0;
}
return v;
}
return 0;
};
/**
converts the data into number of photons for all pixels
\param data pointer to the data
\param nph pointer where the photons should added. If NULL,the internal
image is used \returns pointer to array containing the number of photons
*/
virtual int *getNPhotons(char *data, int *nph = NULL) {
if (nph == NULL)
nph = image;
newFrame(data);
/* cout << fMode << endl; */
/* switch(fMode) { */
/* case eRaw: */
/* cout << "raw" << endl; */
/* break; */
/* default: */
/* cout << "analog" << endl; */
/* } */
// calcGhost(data);
addToCommonMode(data);
for (int iy = ymin; iy < ymax; ++iy) {
for (int ix = xmin; ix < xmax; ++ix) {
switch (fMode) {
case eRaw:
// cout << "raw" << endl;
nph[iy * nx + ix] = det->getChannel(data, ix, iy);
break;
default:
if (det->isGood(ix, iy))
nph[iy * nx + ix] += convertToPhotons(data, ix, iy);
}
}
}
return nph;
}
/**
clears the image array
*/
virtual void clearImage() {
for (int iy = 0; iy < ny; ++iy) {
for (int ix = 0; ix < nx; ++ix) {
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
\param i number of samples to be set (0 or negative gets)
\returns actual number of samples
*/
int SetNPedestals(int i = -1) {
if (i > 0)
for (int iy = 0; iy < ny; ++iy)
for (int ix = 0; ix < nx; ++ix)
stat[iy][ix].SetNPedestals(i);
return stat[0][0].SetNPedestals();
};
/** gets number of samples for moving average pedestal calculation
\returns actual number of samples
*/
int GetNPedestals(int ix, int iy) {
if (ix >= 0 && ix < nx && iy >= 0 && iy < ny)
return stat[iy][ix].GetNPedestals();
else
return -1;
};
/** calculates the sum of photons in the specified region of interest.
\param data pointer to the data
\param xmi minimum x for the calculation. If -1 the minimum x of the
predefined region of interest is used \param xma maximum x for the
calculation. If -1 the maximum x of the predefined region of interest is
used \param ymi minimum y for the calculation. If -1 the minimum y of the
predefined region of interest is used \param yma maximum y for the
calculation. If -1 the maximum y of the predefined region of interest is
used
\returns total number of photons in
*/
virtual int getTotalNumberOfPhotons(char *data, int xmi = -1, int xma = -1,
int ymi = -1, int yma = -1) {
int val = 0;
if (xmi < 0)
xmi = xmin;
if (xma < 0)
xma = xmax;
if (ymi < 0)
ymi = ymin;
if (yma < 0)
yma = ymax;
newFrame(data);
// calcGhost(data);
addToCommonMode(data);
for (int iy = ymi; iy < yma; ++iy)
for (int ix = xmi; ix < xma; ++ix)
if (det->isGood(ix, iy)) {
if (ix >= 0 && ix < nx && iy >= 0 && iy < ny)
val += convertToPhotons(data, ix, iy);
}
return val;
};
/**
calculates the image converted into number of photons. If the frame mode
is pedestal, it also it to the pdedestal subtraction. \param data pointer
to the data to be processed \param val pointer of the data to be added
to. If NULL, the internal image will be used \param pointer to the
processed data \returns
*/
virtual void processData(char *data, int *val = NULL) {
switch (fMode) {
case ePedestal:
// cout << "analog ped " << endl;
addToPedestal(data, 1);
break;
default:
// cout << "analog frame" << endl;
// subtractPedestal(data);
analogDetector<dataType>::getNPhotons(data);
}
};
// virtual char *getInterpolation(){return NULL;};
/** sets the current frame mode for the detector
\param f frame mode to be set
\returns current frame mode
*/
frameMode setFrameMode(frameMode f) {
fMode = f;
return fMode;
};
/** gets the current frame mode for the detector
\returns current frame mode
*/
frameMode getFrameMode() { return fMode; };
// enum detectorMode { eAnalog, ePhotonCounting, eInterpolating };
/** sets the current detector mode for the detector
\param f detector mode to be set
\returns current detector mode
*/
detectorMode setDetectorMode(detectorMode f) {
dMode = f;
return dMode;
};
/** gets the current detector mode for the detector
\returns current detector mode
*/
detectorMode getDetectorMode() { return dMode; };
/** sets file pointer where to write the clusters to
\param f file pointer
\returns current file pointer
*/
FILE *setFilePointer(FILE *f) {
myFile = f;
return myFile;
};
/** gets file pointer where to write the clusters to
\returns current file pointer
*/
FILE *getFilePointer() { return myFile; };
virtual slsInterpolation *getInterpolation() { return NULL; };
virtual slsInterpolation *setInterpolation(slsInterpolation *ii) {
return NULL;
}
virtual double setNSigma(double n) { return 0; };
virtual void setEnergyRange(double emi, double ema) { ; };
protected:
slsDetectorData<dataType> *det; /**< slsDetectorData to be used */
int nx; /**< Size of the detector in x direction */
int ny; /**< Size of the detector in y direction */
pedestalSubtraction **stat; /**< pedestalSubtraction class */
/* double **pedMean; /\**< pedestalSubtraction class *\/ */
/* double **pedVariance; /\**< pedestalSubtraction class *\/ */
commonModeSubtraction *cmSub; /**< commonModeSubtraction class */
int dataSign; /**< sign of the data i.e. 1 if photon is positive, -1 if
negative */
int iframe; /**< frame number (not from file but incremented within the
dataset every time newFrame is called */
double *gmap;
ghostSummation<dataType> *ghSum; /**< ghostSummation class */
int *image;
int id;
// int xmin, xmax, ymin, ymax; int xmin; /**< minimum x of the region of
// interest */
int xmin; /**< minimum x of the region of interest */
int xmax; /**< maximum x of the region of interest */
int ymin; /**< minimum y of the region of interest */
int ymax; /**< maximum y of the region of interest */
double
thr; /**< threshold to be used for conversion into number of photons */
// int nSigma; /**< number of sigma to be used for conversion into number
// of photons if threshold is undefined */
frameMode fMode; /**< current detector frame mode */
detectorMode dMode; /**< 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
};
#endif
View Git Blame Copy Permalink