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

slsdetectorcalibration is a repository on its own

Browse Source
This commit is contained in:
Maliakal Dhanya 2014-07-16 15:34:58 +02:00
parent bbb0a05fca
commit 4814fd56c1
26 changed files with 0 additions and 4185 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

130
slsReceiverSoftware/slsDetectorCalibration/MovingStat.h
View File

@ -1,130 +0,0 @@
#ifndef MOVINGSTAT_H
#define MOVINGSTAT_H
#include <math.h>
class MovingStat
{
/** @short approximated moving average structure */
public:
/** constructor
\param nn number of samples parameter to be used
*/
MovingStat(int nn=1000) : n(nn), m_n(0) {}
/**
clears the moving average number of samples parameter, mean and standard deviation
*/
void Clear()
{
m_n = 0;
m_newM=0;
m_newM2=0;
}
/** sets number of samples parameter
\param i number of samples parameter to be set
*/
void SetN(int i) {if (i>=1) n=i;};
/**
gets number of samples parameter
\returns actual number of samples parameter
*/
int GetN() {return n;};
/** calculates the moving average i.e. adds if number of elements is lower than number of samples parameter, pushes otherwise
\param x value to calculate the moving average
*/
inline void Calc(double x) {
if (m_n<n) Add(x);
else Push(x);
}
/** adds the element to the accumulated average and standard deviation
\param x value to add
*/
inline void Add(double x) {
m_n++;
if (m_n == 1)
{
m_newM = x;
m_newM2 = x*x;
} else {
m_newM = m_newM + x;
m_newM2 = m_newM2 + x*x;
}
}
inline void Push(double x)
{
/** adds the element to the accumulated average and squared mean, while subtracting the current value of the average and squared average
\param x value to push
*/
if (m_n == 0)
{
m_newM = x;
m_newM2 = x*x;
m_n++;
} else {
m_newM = m_newM + x - m_newM/m_n;
m_newM2 = m_newM2 + x*x - m_newM2/m_n;
}
}
/** returns the current number of elements of the moving average
\returns returns the current number of elements of the moving average
*/
int NumDataValues() const
{
return m_n;
}
/** returns the mean, 0 if no elements are inside
\returns returns the mean
*/
inline double Mean() const
{
return (m_n > 0) ? m_newM/m_n : 0.0;
}
/** returns the squared mean, 0 if no elements are inside
\returns returns the squared average
*/
double M2() const
{
return ( (m_n > 1) ? m_newM2/m_n : 0.0 );
}
/** returns the variance, 0 if no elements are inside
\returns returns the variance
*/
inline double Variance() const
{
return ( (m_n > 1) ? (M2()-Mean()*Mean()) : 0.0 );
}
/** returns the standard deviation, 0 if no elements are inside
\returns returns the standard deviation
*/
inline double StandardDeviation() const
{
return ( (Variance() > 0) ? sqrt( Variance() ) : -1 );
}
private:
int n; /**< number of samples parameter */
int m_n; /**< current number of elements */
double m_newM; /**< accumulated average */
double m_newM2; /**< accumulated squared average */
};
#endif

55
slsReceiverSoftware/slsDetectorCalibration/RunningStat.h
View File

@ -1,55 +0,0 @@
class RunningStat
{
public:
RunningStat() : m_n(0) {}
void Clear()
{
m_n = 0;
}
void Push(double x)
{
m_n++;
// See Knuth TAOCP vol 2, 3rd edition, page 232
if (m_n == 1)
{
m_oldM = m_newM = x;
m_oldS = 0.0;
}
else
{
m_newM = m_oldM + (x - m_oldM)/m_n;
m_newS = m_oldS + (x - m_oldM)*(x - m_newM);
// set up for next iteration
m_oldM = m_newM;
m_oldS = m_newS;
}
}
int NumDataValues() const
{
return m_n;
}
double Mean() const
{
return (m_n > 0) ? m_newM : 0.0;
}
double Variance() const
{
return ( (m_n > 1) ? m_newS/(m_n - 1) : 0.0 );
}
double StandardDeviation() const
{
return sqrt( Variance() );
}
private:
int m_n;
double m_oldM, m_newM, m_oldS, m_newS;
};

89
slsReceiverSoftware/slsDetectorCalibration/chiptestBoardData.h
View File

@ -1,89 +0,0 @@
#ifndef CHIPTESTDATA_H
#define CHIPTESTDATA_H
#include "slsDetectorData.h"
class chiptestBoardData : public slsDetectorData<uint16_t> {
public:
/**
chiptestBoard data structure. Works for data acquired using the chiptestBoard.
Inherits and implements slsDetectorData.
Constructor (no error checking if datasize and offsets are compatible!)
\param npx number of pixels in the x direction
\param npy number of pixels in the y direction (1 for strips)
\param nadc number of adcs
\param offset offset at the beginning of the pattern
\param dMap array of size nx*ny storing the pointers to the data in the dataset (as offset)
\param dMask Array of size nx*ny storing the polarity of the data in the dataset (should be 0 if no inversion is required, 0xffffffff is inversion is required)
\param dROI Array of size nx*ny. The elements are 1s if the channel is good or in the ROI, 0 is bad or out of the ROI. NULL (default) means all 1s.
*/
chiptestBoardData(int npx, int npy, int nadc, int offset, int **dMap=NULL, uint16_t **dMask=NULL, int **dROI=NULL): slsDetectorData<uint16_t>(npx, npy, nadc*(npx*npy)+offset, dMap, dMask, dROI), nAdc(nadc), offSize(offset), iframe(0) {}; // should be? nadc*(npx*npy+offset)
/**
Returns the frame number for the given dataset. Virtual func: works for slsDetectorReceiver data (also for each packet), but can be overloaded.
\param buff pointer to the dataset
\returns frame number
*/
virtual int getFrameNumber(char *buff){(void)buff; return iframe;};
/**
Loops over a memory slot until a complete frame is found (i.e. all packets 0 to nPackets, same frame number). Can be overloaded for different kind of detectors!
\param data pointer to the memory to be analyzed
\param ndata size of frame returned
\param dsize size of the memory slot to be analyzed
\returns always return the pointer to data (no frame loss!)
*/
virtual char *findNextFrame(char *data, int &ndata, int dsize) {ndata=dsize;setDataSize(dsize); return data;};
/**
Loops over a file stream until a complete frame is found (i.e. all packets 0 to nPackets, same frame number). Can be overloaded for different kind of detectors!
\param filebin input file stream (binary)
\returns pointer to the first packet of the last good frame, NULL if no frame is found or last frame is incomplete
*/
virtual char *readNextFrame(ifstream &filebin) {
int afifo_length=0;
uint16_t *afifo_cont;
if (filebin.is_open()) {
if (filebin.read((char*)&afifo_length,sizeof(uint32_t))) {
setDataSize(afifo_length*nAdc*sizeof(uint16_t));
afifo_cont=new uint16_t[afifo_length*nAdc];
if (filebin.read((char*)afifo_cont,afifo_length*sizeof(uint16_t)*nAdc)) {
iframe++;
return (char*)afifo_cont;
} else {
delete [] afifo_cont;
return NULL;
}
} else {
return NULL;
}
}
return NULL;
};
private:
const int nAdc; /**<number of ADC read out */
const int offSize; /**< offset at the beginning of the frame (depends on the pattern) */
int iframe; /**< frame number (calculated in software! not in the data)*/
};
#endif

82
slsReceiverSoftware/slsDetectorCalibration/commonModeSubtraction.h
View File

@ -1,82 +0,0 @@
#ifndef COMMONMODESUBTRACTION_H
#define COMMONMODESUBTRACTION_H
#include "MovingStat.h"
class 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
*/
commonModeSubtraction(int nn=1000, int iroi=1) : cmStat(NULL), cmPed(NULL), nCm(NULL), nROI(iroi) {cmStat=new MovingStat[nROI]; for (int i=0; i<nROI; i++) cmStat[i].SetN(nn); cmPed=new double[nROI]; nCm=new double[nROI];};
/** destructor - deletes the moving average(s) and the sum of pedestals calculator(s) */
virtual ~commonModeSubtraction() {delete [] cmStat; delete [] cmPed; delete [] nCm;};
/** clears the moving average and the sum of pedestals calculation - virtual func*/
virtual void Clear(){
for (int i=0; i<nROI; i++) {
cmStat[i].Clear();
nCm[i]=0;
cmPed[i]=0;
}};
/** adds the average of pedestals to the moving average and reinitializes the calculation of the sum of pedestals for all ROIs. - virtual func*/
virtual void newFrame(){
for (int i=0; i<nROI; i++) {
if (nCm[i]>0) cmStat[i].Calc(cmPed[i]/nCm[i]);
nCm[i]=0;
cmPed[i]=0;
}};
/** adds the pixel to the sum of pedestals -- virtual func must be overloaded to define the regions of interest
\param val value to add
\param ix pixel x coordinate
\param iy pixel y coordinate
*/
virtual void addToCommonMode(double val, int ix=0, int iy=0) {
(void) ix; (void) iy;
//if (isc>=0 && isc<nROI) {
cmPed[0]+=val;
nCm[0]++;//}
};
/** gets the common mode i.e. the difference between the current average sum of pedestals mode and the average pedestal -- virtual func must be overloaded to define the regions of interest
\param ix pixel x coordinate
\param iy pixel y coordinate
\return the difference between the current average sum of pedestals and the average pedestal
*/
virtual double getCommonMode(int ix=0, int iy=0) {
(void) ix; (void) iy;
if (nCm[0]>0) return cmPed[0]/nCm[0]-cmStat[0].Mean();
return 0;};
protected:
MovingStat *cmStat; /**<array of moving average of the pedestal average per region of interest */
double *cmPed; /**< array storing the sum of pedestals per region of interest */
double *nCm; /**< array storing the number of pixels currently contributing to the pedestals */
const int nROI; /**< constant parameter for number of regions on which the common mode should be calculated separately e.g. supercolumns */
};
#endif

31
slsReceiverSoftware/slsDetectorCalibration/demoCreateTree.C
View File

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

85
slsReceiverSoftware/slsDetectorCalibration/doxy.config
View File

@ -1,85 +0,0 @@
# If the EXTRACT_ALL tag is set to YES doxygen will assume all entities in
# documentation are documented, even if no documentation was available.
# Private class members and static file members will be hidden unless
# the EXTRACT_PRIVATE and EXTRACT_STATIC tags are set to YES
EXTRACT_ALL = YES
# If the EXTRACT_PRIVATE tag is set to YES all private members of a class
# will be included in the documentation.
EXTRACT_PRIVATE = NO
# If the EXTRACT_STATIC tag is set to YES all static members of a file
# will be included in the documentation.
EXTRACT_STATIC = YES
# If the EXTRACT_LOCAL_CLASSES tag is set to YES classes (and structs)
# defined locally in source files will be included in the documentation.
# If set to NO only classes defined in header files are included.
EXTRACT_LOCAL_CLASSES = YES
# This flag is only useful for Objective-C code. When set to YES local
# methods, which are defined in the implementation section but not in
# the interface are included in the documentation.
# If set to NO (the default) only methods in the interface are included.
EXTRACT_LOCAL_METHODS = YES
# If this flag is set to YES, the members of anonymous namespaces will be
# extracted and appear in the documentation as a namespace called
# 'anonymous_namespace{file}', where file will be replaced with the base
# name of the file that contains the anonymous namespace. By default
# anonymous namespace are hidden.
EXTRACT_ANON_NSPACES = NO
# If the HIDE_UNDOC_MEMBERS tag is set to YES, Doxygen will hide all
# undocumented members of documented classes, files or namespaces.
# If set to NO (the default) these members will be included in the
# various overviews, but no documentation section is generated.
# This option has no effect if EXTRACT_ALL is enabled.
HIDE_UNDOC_MEMBERS = NO
# If the HIDE_UNDOC_CLASSES tag is set to YES, Doxygen will hide all
# undocumented classes that are normally visible in the class hierarchy.
# If set to NO (the default) these classes will be included in the various
# overviews. This option has no effect if EXTRACT_ALL is enabled.
HIDE_UNDOC_CLASSES = NO
# If the HIDE_FRIEND_COMPOUNDS tag is set to YES, Doxygen will hide all
# friend (class|struct|union) declarations.
# If set to NO (the default) these declarations will be included in the
# documentation.
HIDE_FRIEND_COMPOUNDS = NO
INTERNAL_DOCS = NO
SHOW_INCLUDE_FILES = NO
SHOW_FILES = NO
SHOW_NAMESPACES = NO
COMPACT_LATEX = YES
PAPER_TYPE = a4
PDF_HYPERLINKS = YES
USE_PDFLATEX = YES
LATEX_HIDE_INDICES = YES
PREDEFINED = __cplusplus
INPUT = MovingStat.h slsDetectorData.h slsReceiverData.h moench02ModuleData.h pedestalSubtraction.h commonModeSubtraction.h moenchCommonMode.h singlePhotonDetector.h energyCalibration.h moenchReadData.C single_photon_hit.h chiptestBoardData.h jungfrau02Data.h jungfrauReadData.C jungfrau02CommonMode.h
OUTPUT_DIRECTORY = docs

527
slsReceiverSoftware/slsDetectorCalibration/energyCalibration.cpp
View File

@ -1,527 +0,0 @@
#include "energyCalibration.h"
#ifdef __CINT
#define MYROOT
#endif
#ifdef MYROOT
#include <TMath.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TGraphErrors.h>
#endif
#include <iostream>
#define max(a,b) ((a) > (b) ? (a) : (b))
#define min(a,b) ((a) < (b) ? (a) : (b))
#define ELEM_SWAP(a,b) { register int t=(a);(a)=(b);(b)=t; }
using namespace std;
#ifdef MYROOT
Double_t energyCalibrationFunctions::pedestal(Double_t *x, Double_t *par) {
return par[0]-par[1]*sign*x[0];
}
Double_t energyCalibrationFunctions::gaussChargeSharing(Double_t *x, Double_t *par) {
Double_t f, arg=0;
if (par[3]!=0) arg=sign*(x[0]-par[2])/par[3];
f=TMath::Exp(-1*arg*arg/2.);
f=f+par[5]/2.*(TMath::Erfc(arg/(TMath::Sqrt(2.))));
return par[4]*f+pedestal(x,par);
}
Double_t energyCalibrationFunctions::gaussChargeSharingPixel(Double_t *x, Double_t *par) {
Double_t f;
if (par[3]<=0 || par[2]*(*x)<=0 || par[5]<0 || par[4]<=0) return 0;
Double_t pp[3];
pp[0]=0;
pp[1]=par[2];
pp[2]=par[3];
f=(par[5]-par[6]*(TMath::Log(*x/par[2])))*erfBox(x,pp);
f+=par[4]*TMath::Gaus(*x, par[2], par[3], kTRUE);
return f+pedestal(x,par);
}
Double_t energyCalibrationFunctions::erfBox(Double_t *z, Double_t *par) {
Double_t m=par[0];
Double_t M=par[1];
if (par[0]>par[1]) {
m=par[1];
M=par[0];
}
if (m==M)
return 0;
if (par[2]<=0) {
if (*z>=m && *z<=M)
return 1./(M-m);
else
return 0;
}
return (TMath::Erfc((z[0]-M)/par[2])-TMath::Erfc((z[0]-m)/par[2]))*0.5/(M-m);
}
// basic erf function
Double_t energyCalibrationFunctions::erfFunction(Double_t *x, Double_t *par) {
double arg=0;
if (par[1]!=0) arg=(par[0]-x[0])/par[1];
return ((par[2]/2.*(1+TMath::Erf(sign*arg/(TMath::Sqrt(2))))));
};
Double_t energyCalibrationFunctions::erfFunctionChargeSharing(Double_t *x, Double_t *par) {
Double_t f;
f=erfFunction(x, par+2)*(1+par[5]*(par[2]-x[0]))+par[0]-par[1]*x[0]*sign;
return f;
};
Double_t energyCalibrationFunctions::erfFuncFluo(Double_t *x, Double_t *par) {
Double_t f;
f=erfFunctionChargeSharing(x, par)+erfFunction(x, par+6)*(1+par[9]*(par[6]-x[0]));
return f;
};
#endif
double energyCalibrationFunctions::median(double *x, int n){
// sorts x into xmed array and returns median
// n is number of values already in the xmed array
double xmed[n];
int k,i,j;
for (i=0; i<n; i++) {
k=0;
for (j=0; j<n; j++) {
if(*(x+i)>*(x+j))
k++;
if (*(x+i)==*(x+j)) {
if (i>j)
k++;
}
}
xmed[k]=*(x+i);
}
k=n/2;
return xmed[k];
}
int energyCalibrationFunctions::quick_select(int arr[], int n){
int low, high ;
int median;
int middle, ll, hh;
low = 0 ; high = n-1 ; median = (low + high) / 2;
for (;;) {
if (high <= low) /* One element only */
return arr[median] ;
if (high == low + 1) { /* Two elements only */
if (arr[low] > arr[high])
ELEM_SWAP(arr[low], arr[high]) ;
return arr[median] ;
}
/* Find median of low, middle and high items; swap into position low */
middle = (low + high) / 2;
if (arr[middle] > arr[high]) ELEM_SWAP(arr[middle], arr[high]) ;
if (arr[low] > arr[high]) ELEM_SWAP(arr[low], arr[high]) ;
if (arr[middle] > arr[low]) ELEM_SWAP(arr[middle], arr[low]) ;
/* Swap low item (now in position middle) into position (low+1) */
ELEM_SWAP(arr[middle], arr[low+1]) ;
/* Nibble from each end towards middle, swapping items when stuck */
ll = low + 1;
hh = high;
for (;;) {
do ll++; while (arr[low] > arr[ll]) ;
do hh--; while (arr[hh] > arr[low]) ;
if (hh < ll)
break;
ELEM_SWAP(arr[ll], arr[hh]) ;
}
/* Swap middle item (in position low) back into correct position */
ELEM_SWAP(arr[low], arr[hh]) ;
/* Re-set active partition */
if (hh <= median)
low = ll;
if (hh >= median)
high = hh - 1;
}
}
int energyCalibrationFunctions::kth_smallest(int *a, int n, int k){
register int i,j,l,m ;
register double x ;
l=0 ; m=n-1 ;
while (l<m) {
x=a[k] ;
i=l ;
j=m ;
do {
while (a[i]<x) i++ ;
while (x<a[j]) j-- ;
if (i<=j) {
ELEM_SWAP(a[i],a[j]) ;
i++ ; j-- ;
}
} while (i<=j) ;
if (j<k) l=i ;
if (k<i) m=j ;
}
return a[k] ;
}
#ifdef MYROOT
Double_t energyCalibrationFunctions::spectrum(Double_t *x, Double_t *par) {
return gaussChargeSharing(x,par);
}
Double_t energyCalibrationFunctions::spectrumPixel(Double_t *x, Double_t *par) {
return gaussChargeSharingPixel(x,par);
}
Double_t energyCalibrationFunctions::scurve(Double_t *x, Double_t *par) {
return erfFunctionChargeSharing(x,par);
}
Double_t energyCalibrationFunctions::scurveFluo(Double_t *x, Double_t *par) {
return erfFuncFluo(x,par);
}
#endif
energyCalibration::energyCalibration() :
#ifdef MYROOT
fit_min(-1),
fit_max(-1),
bg_offset(-1),
bg_slope(-1),
flex(-1),
noise(-1),
ampl(-1),
cs_slope(-1),
fscurve(NULL),
fspectrum(NULL),
#endif
funcs(NULL),
plot_flag(1), // fit parameters output to screen
cs_flag(1)
{
#ifdef MYROOT
funcs=new energyCalibrationFunctions();
fscurve=new TF1("fscurve",funcs,&energyCalibrationFunctions::scurve,0,1000,6,"energyCalibrationFunctions","scurve");
fscurve->SetParNames("Background Offset","Background Slope","Inflection Point","Noise RMS", "Number of Photons","Charge Sharing Slope");
fspectrum=new TF1("fspectrum",funcs,&energyCalibrationFunctions::spectrum,0,1000,6,"energyCalibrationFunctions","spectrum");
fspectrum->SetParNames("Background Pedestal","Background slope", "Peak position","Noise RMS", "Number of Photons","Charge Sharing Pedestal");
fspixel=new TF1("fspixel",funcs,&energyCalibrationFunctions::spectrumPixel,0,1000,7,"energyCalibrationFunctions","spectrumPixel");
fspixel->SetParNames("Background Pedestal","Background slope", "Peak position","Noise RMS", "Number of Photons","Charge Sharing Pedestal","Corner");
#endif
}
void energyCalibration::fixParameter(int ip, Double_t val){
fscurve->FixParameter(ip, val);
fspectrum->FixParameter(ip, val);
}
void energyCalibration::releaseParameter(int ip){
fscurve->ReleaseParameter(ip);
fspectrum->ReleaseParameter(ip);
}
energyCalibration::~energyCalibration(){
#ifdef MYROOT
delete fscurve;
delete fspectrum;
#endif
}
#ifdef MYROOT
TH1F* energyCalibration::createMedianHistogram(TH2F* h2, int ch0, int nch, int direction) {
if (h2==NULL || nch==0)
return NULL;
double *x=new double[nch];
TH1F *h1=NULL;
double val=-1;
if (direction==0) {
h1=new TH1F("median","Median",h2->GetYaxis()->GetNbins(),h2->GetYaxis()->GetXmin(),h2->GetYaxis()->GetXmax());
for (int ib=0; ib<h1->GetXaxis()->GetNbins(); ib++) {
for (int ich=0; ich<nch; ich++) {
x[ich]=h2->GetBinContent(ch0+ich+1,ib+1);
}
val=energyCalibrationFunctions::median(x, nch);
h1->SetBinContent(ib+1,val);
}
} else if (direction==1) {
h1=new TH1F("median","Median",h2->GetXaxis()->GetNbins(),h2->GetXaxis()->GetXmin(),h2->GetXaxis()->GetXmax());
for (int ib=0; ib<h1->GetYaxis()->GetNbins(); ib++) {
for (int ich=0; ich<nch; ich++) {
x[ich]=h2->GetBinContent(ib+1,ch0+ich+1);
}
val=energyCalibrationFunctions::median(x, nch);
h1->SetBinContent(ib+1,val);
}
}
delete [] x;
return h1;
}
void energyCalibration::setStartParameters(Double_t *par){
bg_offset=par[0];
bg_slope=par[1];
flex=par[2];
noise=par[3];
ampl=par[4];
cs_slope=par[5];
}
void energyCalibration::getStartParameters(Double_t *par){
par[0]=bg_offset;
par[1]=bg_slope;
par[2]=flex;
par[3]=noise;
par[4]=ampl;
par[5]=cs_slope;
}
#endif
int energyCalibration::setChargeSharing(int p) {
if (p>=0) {
cs_flag=p;
#ifdef MYROOT
if (p) {
fscurve->ReleaseParameter(5);
fspectrum->ReleaseParameter(1);
} else {
fscurve->FixParameter(5,0);
fspectrum->FixParameter(1,0);
}
#endif
}
return cs_flag;
}
#ifdef MYROOT
void energyCalibration::initFitFunction(TF1 *fun, TH1 *h1) {
Double_t min=fit_min, max=fit_max;
Double_t mypar[6];
if (max==-1)
max=h1->GetXaxis()->GetXmax();
if (min==-1)
min=h1->GetXaxis()->GetXmin();
if (bg_offset==-1)
mypar[0]=0;
else
mypar[0]=bg_offset;
if (bg_slope==-1)
mypar[1]=0;
else
mypar[1]=bg_slope;
if (flex==-1)
mypar[2]=(min+max)/2.;
else
mypar[2]=flex;
if (noise==-1)
mypar[3]=0.1;
else
mypar[3]=noise;
if (ampl==-1)
mypar[4]=h1->GetBinContent(h1->GetXaxis()->FindBin(0.5*(max+min)));
else
mypar[4]=ampl;
if (cs_slope==-1)
mypar[5]=0;
else
mypar[5]=cs_slope;
fun->SetParameters(mypar);
fun->SetRange(min,max);
}
TF1* energyCalibration::fitFunction(TF1 *fun, TH1 *h1, Double_t *mypar, Double_t *emypar) {
TF1* fitfun;
char fname[100];
strcpy(fname, fun->GetName());
if (plot_flag) {
h1->Fit(fname,"R0Q");
} else
h1->Fit(fname,"R0Q");
fitfun= h1->GetFunction(fname);
fitfun->GetParameters(mypar);
for (int ip=0; ip<6; ip++) {
emypar[ip]=fitfun->GetParError(ip);
}
return fitfun;
}
TF1* energyCalibration::fitSCurve(TH1 *h1, Double_t *mypar, Double_t *emypar) {
initFitFunction(fscurve,h1);
return fitFunction(fscurve, h1, mypar, emypar);
}
TF1* energyCalibration::fitSpectrum(TH1 *h1, Double_t *mypar, Double_t *emypar) {
initFitFunction(fspectrum,h1);
return fitFunction(fspectrum, h1, mypar, emypar);
}
TGraphErrors* energyCalibration::linearCalibration(int nscan, Double_t *en, Double_t *een, Double_t *fl, Double_t *efl, Double_t &gain, Double_t &off, Double_t &egain, Double_t &eoff) {
TGraphErrors *gr;
Double_t mypar[2];
gr = new TGraphErrors(nscan,en,fl,een,efl);
if (plot_flag) {
gr->Fit("pol1");
gr->SetMarkerStyle(20);
} else
gr->Fit("pol1","0Q");
TF1 *fitfun= gr->GetFunction("pol1");
fitfun->GetParameters(mypar);
egain=fitfun->GetParError(1);
eoff=fitfun->GetParError(0);
gain=funcs->setScanSign()*mypar[1];
off=mypar[0];
return gr;
}
TGraphErrors* energyCalibration::calibrate(int nscan, Double_t *en, Double_t *een, TH1F **h1, Double_t &gain, Double_t &off, Double_t &egain, Double_t &eoff, int integral) {
TH1F *h;
Double_t mypar[6], emypar[6];
Double_t fl[nscan], efl[nscan];
for (int ien=0; ien<nscan; ien++) {
h=h1[ien];
if (integral)
fitSCurve(h,mypar,emypar);
else
fitSpectrum(h,mypar,emypar);
fl[ien]=mypar[2];
efl[ien]=emypar[2];
}
return linearCalibration(nscan,en,een,fl,efl,gain,off, egain, eoff);
}
#endif
Fit data is empty

452
slsReceiverSoftware/slsDetectorCalibration/energyCalibration.h
View File

@ -1,452 +0,0 @@
#ifndef ENERGYCALIBRATION_H
#define ENERGYCALIBRATION_H
#ifdef __CINT__
#define MYROOT
#endif
#ifdef G__ROOT
#define MYROOT
#endif
#ifdef __MAKECINT__
#define MYROOT
#endif
#ifdef ROOT_VERSION
#define MYROOT
#endif
#define MYROOT
#ifdef MYROOT
#include <TROOT.h>
#include <TF1.h>
class TH1F;
class TH2F;
class TGraphErrors;
#endif
using namespace std;
const double conven=1000./3.6; /**< electrons/keV */
const double el=1.67E-4; /**< electron charge in fC */
/**
\mainpage Common Root library for SLS detectors data analysis
*
* \section intro_sec Introduction
We know very well s-curves etc. but at the end everybody uses different functions ;-).
* \subsection mot_sec Motivation
It would be greate to use everybody the same functions...
*/
/**
*
*
@libdoc The energiCalibration class contains all the necessary functions for s-curve fitting and linear calibration of the threshold.
*
* @short Energy calibration functions
* @author Anna Bergamaschi
* @version 0.1alpha
*/
/**
class containing all the possible energy calibration functions (scurves with and without charge sharing, gaussian spectrum with and without charge sharing, possibility of chosing the sign of the X-axis)
*/
class energyCalibrationFunctions {
public:
energyCalibrationFunctions(int s=-1) {setScanSign(s);};
/** sets scan sign
\param s can be 1 (energy and x-axis have the same direction) or -1 (energy and x-axis have opposite directions) otherwise gets
\returns current scan sign can be 1 (energy and x-axis have the same direction) or -1 (energy and x-axis have opposite directions)
*/
int setScanSign(int s=0) {if (s==1 || s==-1) sign=s; return sign;};;
#ifdef MYROOT
/**
Gaussian Function with charge sharing pedestal
par[0] is the absolute height of the background pedestal
par[1] is the slope of the background pedestal
*/
Double_t pedestal(Double_t *x, Double_t *par);
/**
Gaussian Function with charge sharing pedestal
par[0] is the absolute height of the background pedestal
par[1] is the slope of the background pedestal
par[2] is the gaussian peak position
par[3] is the RMS of the gaussian (and of the pedestal)
par[4] is the height of the function
par[5] is the fractional height of the charge sharing pedestal (scales with par[3])
*/
Double_t gaussChargeSharing(Double_t *x, Double_t *par);
/**
Gaussian Function with charge sharing pedestal
par[0] is the absolute height of the background pedestal
par[1] is the slope of the background pedestal
par[2] is the gaussian peak position
par[3] is the RMS of the gaussian (and of the pedestal)
par[4] is the height of the function
par[5] is the fractional height of the charge sharing pedestal (scales with par[3])
*/
Double_t gaussChargeSharingPixel(Double_t *x, Double_t *par);
/**
Basic erf function
par[0] is the inflection point
par[1] is the RMS
par[2] is the amplitude
*/
Double_t erfFunction(Double_t *x, Double_t *par) ;
Double_t erfBox(Double_t *z, Double_t *par);
/** Erf function with charge sharing slope
par[0] is the pedestal
par[1] is the slope of the pedestal
par[2] is the inflection point
par[3] is the RMS
par[4] is the amplitude
par[5] is the angual coefficient of the charge sharing slope (scales with par[3])
*/
Double_t erfFunctionChargeSharing(Double_t *x, Double_t *par);
/** Double Erf function with charge sharing slope
par[0] is the pedestal
par[1] is the slope of the pedestal
par[2] is the inflection point of the first energy
par[3] is the RMS of the first energy
par[4] is the amplitude of the first energy
par[5] is the angual coefficient of the charge sharing slope of the first energy (scales with par[3])
par[6] is the inflection point of the second energy
par[7] is the RMS of the second energy
par[8] is the amplitude of the second energy
par[9] is the angual coefficient of the charge sharing slope of the second energy (scales with par[8])
*/
Double_t erfFuncFluo(Double_t *x, Double_t *par);
/**
static function Gaussian with charge sharing pedestal with the correct scan sign
par[0] is the absolute height of the background pedestal
par[1] is the slope of the pedestal
par[2] is the gaussian peak position
par[3] is the RMS of the gaussian (and of the pedestal)
par[4] is the height of the function
par[5] is the fractional height of the charge sharing pedestal (scales with par[4]
*/
Double_t spectrum(Double_t *x, Double_t *par);
/**
static function Gaussian with charge sharing pedestal with the correct scan sign
par[0] is the absolute height of the background pedestal
par[1] is the slope of the pedestal
par[2] is the gaussian peak position
par[3] is the RMS of the gaussian (and of the pedestal)
par[4] is the height of the function
par[5] is the fractional height of the charge sharing pedestal (scales with par[4]
*/
Double_t spectrumPixel(Double_t *x, Double_t *par);
/** Erf function with charge sharing slope with the correct scan sign
par[0] is the pedestal
par[1] is the slope of the pedestal
par[2] is the inflection point
par[3] is the RMS
par[4] is the amplitude
par[5] is the angual coefficient of the charge sharing slope (scales with par[3])
*/
Double_t scurve(Double_t *x, Double_t *par);
/** Double Erf function with charge sharing slope
par[0] is the pedestal
par[1] is the slope of the pedestal
par[2] is the inflection point of the first energy
par[3] is the RMS of the first energy
par[4] is the amplitude of the first energy
par[5] is the angual coefficient of the charge sharing slope of the first energy (scales with par[3])
par[6] is the inflection point of the second energy
par[7] is the RMS of the second energy
par[8] is the amplitude of the second energy
par[9] is the angual coefficient of the charge sharing slope of the second energy (scales with par[8])
*/
Double_t scurveFluo(Double_t *x, Double_t *par);
#endif
/** Calculates the median of an array of n elements */
static double median(double *x, int n);
/** Calculates the median of an array of n elements (swaps the arrays!)*/
static int quick_select(int arr[], int n);
/** Calculates the median of an array of n elements (swaps the arrays!)*/
static int kth_smallest(int *a, int n, int k);
private:
int sign;
};
/**
class alowing the energy calibration of photon counting and anlogue detectors
*/
class energyCalibration {
public:
/**
default constructor - creates the function with which the s-curves will be fitted
*/
energyCalibration();
/**
default destructor - deletes the function with which the s-curves will be fitted
*/
~energyCalibration();
/** sets plot flag
\param p plot flag (-1 gets, 0 unsets, >0 plot)
\returns current plot flag
*/
int setPlotFlag(int p=-1) {if (p>=0) plot_flag=p; return plot_flag;};
/** sets scan sign
\param s can be 1 (energy and x-axis have the same direction) or -1 (energy and x-axis have opposite directions) otherwise gets
\returns current scan sign can be 1 (energy and x-axis have the same direction) or -1 (energy and x-axis have opposite directions)
*/
int setScanSign(int s=0) {return funcs->setScanSign(s);};
/** sets plot flag
\param p plot flag (-1 gets, 0 unsets, >0 plot)
\returns current plot flag
*/
int setChargeSharing(int p=-1);
void fixParameter(int ip, Double_t val);
void releaseParameter(int ip);
#ifdef MYROOT
/**
Creates an histogram with the median of nchannels starting from a specified one. the direction on which it is mediated can be selected (defaults to x=0)
\param h2 2D histogram on which the median will be calculated
\param ch0 starting channel
\param nch number of channels to be mediated
\param direction can be either 0 (x, default) or 1 (y)
\returns a TH1F histogram with the X-axis as a clone of the h2 Y (if direction=0) or X (if direction=0) axis, and on the Y axis the median of the counts of the mediated channels f h2
*/
static TH1F* createMedianHistogram(TH2F* h2, int ch0, int nch, int direction=0);
/** sets the s-curve fit range
\param mi minimum of the fit range (-1 is histogram x-min)
\param ma maximum of the fit range (-1 is histogram x-max)
*/
void setFitRange(Double_t mi, Double_t ma){fit_min=mi; fit_max=ma;};
/** gets the s-curve fit range
\param mi reference for minimum of the fit range (-1 is histogram x-min)
\param ma reference for maximum of the fit range (-1 is histogram x-max)
*/
void getFitRange(Double_t &mi, Double_t &ma){mi=fit_min; ma=fit_max;};
/** set start parameters for the s-curve function
\param par parameters, -1 sets to auto-calculation
par[0] is the pedestal
par[1] is the slope of the pedestal
par[2] is the inflection point
par[3] is the RMS
par[4] is the amplitude
par[5] is the angual coefficient of the charge sharing slope (scales with par[3]) -- always positive
*/
void setStartParameters(Double_t *par);
/** get start parameters for the s-curve function
\param par parameters, -1 means auto-calculated
par[0] is the pedestal
par[1] is the slope of the pedestal
par[2] is the inflection point
par[3] is the RMS
par[4] is the amplitude
par[5] is the angual coefficient of the charge sharing slope (scales with par[3]) -- always positive
*/
void getStartParameters(Double_t *par);
/**
fits histogram with the s-curve function
\param h1 1d-histogram to be fitted
\param mypar pointer to fit parameters array
\param emypar pointer to fit parameter errors
\returns the fitted function - can be used e.g. to get the Chi2 or similar
*/
TF1 *fitSCurve(TH1 *h1, Double_t *mypar, Double_t *emypar);
/**
fits histogram with the spectrum
\param h1 1d-histogram to be fitted
\param mypar pointer to fit parameters array
\param emypar pointer to fit parameter errors
\returns the fitted function - can be used e.g. to get the Chi2 or similar
*/
TF1 *fitSpectrum(TH1 *h1, Double_t *mypar, Double_t *emypar);
/**
calculates gain and offset for the set of inflection points
\param nscan number of energy scans
\param en array of energies (nscan long)
\param een array of errors on energies (nscan long) - can be NULL!
\param fl array of inflection points (nscan long)
\param efl array of errors on the inflection points (nscan long)
\param gain reference to gain resulting from the fit
\param off reference to offset resulting from the fit
\param egain reference to error on the gain resulting from the fit
\param eoff reference to the error on the offset resulting from the fit
\returns graph energy vs inflection point
*/
TGraphErrors* linearCalibration(int nscan, Double_t *en, Double_t *een, Double_t *fl, Double_t *efl, Double_t &gain, Double_t &off, Double_t &egain, Double_t &eoff);
/**
calculates gain and offset for the set of energy scans
\param nscan number of energy scans
\param en array of energies (nscan long)
\param een array of errors on energies (nscan long) - can be NULL!
\param h1 array of TH1
\param gain reference to gain resulting from the fit
\param off reference to offset resulting from the fit
\param egain reference to error on the gain resulting from the fit
\param eoff reference to the error on the offset resulting from the fit
\returns graph energy vs inflection point
*/
TGraphErrors* calibrateScurves(int nscan, Double_t *en, Double_t *een, TH1F **h1, Double_t &gain, Double_t &off, Double_t &egain, Double_t &eoff){return calibrate(nscan, en, een, h1, gain, off, egain, eoff, 1);};
/**
calculates gain and offset for the set of energy spectra
\param nscan number of energy scans
\param en array of energies (nscan long)
\param een array of errors on energies (nscan long) - can be NULL!
\param h1 array of TH1
\param gain reference to gain resulting from the fit
\param off reference to offset resulting from the fit
\param egain reference to error on the gain resulting from the fit
\param eoff reference to the error on the offset resulting from the fit
\returns graph energy vs peak
*/
TGraphErrors* calibrateSpectra(int nscan, Double_t *en, Double_t *een, TH1F **h1, Double_t &gain, Double_t &off, Double_t &egain, Double_t &eoff){return calibrate(nscan, en, een, h1, gain, off, egain, eoff, 0);};
#endif
private:
#ifdef MYROOT
/**
calculates gain and offset for the set of energies
\param nscan number of energy scans
\param en array of energies (nscan long)
\param een array of errors on energies (nscan long) - can be NULL!
\param h1 array of TH1
\param gain reference to gain resulting from the fit
\param off reference to offset resulting from the fit
\param egain reference to error on the gain resulting from the fit
\param eoff reference to the error on the offset resulting from the fit
\param integral 1 is an s-curve set (default), 0 spectra
\returns graph energy vs peak/inflection point
*/
TGraphErrors* calibrate(int nscan, Double_t *en, Double_t *een, TH1F **h1, Double_t &gain, Double_t &off, Double_t &egain, Double_t &eoff, int integral=1);
/**
Initializes the start parameters and the range of the fit depending on the histogram characteristics and/or on the start parameters specified by the user
\param fun pointer to function to be initialized
\param h1 histogram from which to extract the range and start parameters, if not already specified by the user
*/
void initFitFunction(TF1 *fun, TH1 *h1);
/**
Performs the fit according to the flags specified and returns the fitted function
\param fun function to fit
\param h1 histogram to fit
\param mypar pointer to fit parameters array
\param emypar pointer to fit parameter errors
\returns the fitted function - can be used e.g. to get the Chi2 or similar
*/
TF1 *fitFunction(TF1 *fun, TH1 *h1, Double_t *mypar, Double_t *emypar);
#endif
#ifdef MYROOT
Double_t fit_min; /**< minimum of the s-curve fitting range, -1 is histogram x-min */
Double_t fit_max; /**< maximum of the s-curve fitting range, -1 is histogram x-max */
Double_t bg_offset; /**< start value for the background pedestal */
Double_t bg_slope; /**< start value for the background slope */
Double_t flex; /**< start value for the inflection point */
Double_t noise; /**< start value for the noise */
Double_t ampl; /**< start value for the number of photons */
Double_t cs_slope; /**< start value for the charge sharing slope */
TF1 *fscurve; /**< function with which the s-curve will be fitted */
TF1 *fspectrum; /**< function with which the spectrum will be fitted */
TF1 *fspixel; /**< function with which the spectrum will be fitted */
#endif
energyCalibrationFunctions *funcs;
int plot_flag; /**< 0 does not plot, >0 plots (flags?) */
int cs_flag; /**< 0 functions without charge sharing contribution, >0 with charge sharing contribution */
};
#endif

11
slsReceiverSoftware/slsDetectorCalibration/gMapDemo.C
View File

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

228
slsReceiverSoftware/slsDetectorCalibration/gainMap.C
View File

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

148
slsReceiverSoftware/slsDetectorCalibration/gotthardModuleData.h
View File

@ -1,148 +0,0 @@
#ifndef GOTTHARDMODULEDATA_H
#define GOTTHARDMODULEDATA_H
#include "slsReceiverData.h"
#define FRAMEMASK 0xFFFFFFFE
#define PACKETMASK 1
#define FRAMEOFFSET 0x1
class gotthardModuleData : public slsReceiverData<uint16_t> {
public:
/**
Implements the slsReceiverData structure for the gotthard read out by a module i.e. using the slsReceiver
(1x1280 pixels, 2 packets 1286 large etc.)
\param c crosstalk parameter for the output buffer
*/
gotthardModuleData(double c=0): slsReceiverData<uint16_t>(xpixels, ypixels, npackets, buffersize), xtalk(c) {
uint16_t **dMask;
int **dMap;
int ix, iy;
int initial_offset = 2;
int offset = initial_offset;
dMask=new uint16_t*[ypixels];
dMap=new int*[ypixels];
for (int i = 0; i < ypixels; i++) {
dMap[i] = new int[xpixels];
dMask[i] = new uint16_t[xpixels];
}
for(ix=0; ix<ypixels; ++ix)
for(iy=0; iy<xpixels; ++iy)
dMask[ix][iy] = 0x0;
for(ix=0; ix<ypixels; ++ix){
if(ix == (ypixels/2)){
offset += initial_offset;//2
offset++;
}
for(iy=0; iy<xpixels; ++iy){
dMap[ix][iy] = offset;
offset++;
}
}
setDataMap(dMap);
setDataMask(dMask);
};
/**
Returns the frame number for the given dataset.
\param buff pointer to the dataset
\returns frame number
*/
int getFrameNumber(char *buff){
int np=(*(int*)buff);
//gotthards frame header must be incremented
++np;
//packet index should be 1 or 2
return ((np&FRAMEMASK)>>FRAMEOFFSET);
};
/**
gets the packets number (last packet is labelled with 0 and is replaced with 40)
\param buff pointer to the memory
\returns packet number
*/
int getPacketNumber(char *buff){
int np=(*(int*)buff);
//gotthards frame header must be incremented
++np;
//packet index should be 1 or 2
return ((np&PACKETMASK)+1);
};
/**
returns the pixel value as double correcting for the output buffer crosstalk
\param data pointer to the memory
\param ix coordinate in the x direction
\param iy coordinate in the y direction
\returns channel value as double
*/
double getValue(char *data, int ix, int iy=0) {
//check how it is for gotthard
if (xtalk==0)
return slsDetectorData<uint16_t>::getValue(data, ix, iy);
else
return slsDetectorData<uint16_t>::getValue(data, ix, iy)-xtalk*slsDetectorData<uint16_t>::getValue(data, ix-1, iy);
};
/** sets the output buffer crosstalk correction parameter
\param c output buffer crosstalk correction parameter to be set
\returns current value for the output buffer crosstalk correction parameter
*/
double setXTalk(double c) {xtalk=c; return xtalk;}
/** gets the output buffer crosstalk parameter
\returns current value for the output buffer crosstalk correction parameter
*/
double getXTalk() {return xtalk;}
private:
double xtalk; /**<output buffer crosstalk correction parameter */
const static int xpixels = 1280;
const static int ypixels = 1;
const static int npackets = 2;
const static int buffersize = 1286;
};
#endif

127
slsReceiverSoftware/slsDetectorCalibration/gotthardShortModuleData.h
View File

@ -1,127 +0,0 @@
#ifndef GOTTHARDSHORTMODULEDATA_H
#define GOTTHARDSHORTMODULEDATA_H
#include "slsReceiverData.h"
class gotthardShortModuleData : public slsReceiverData<uint16_t> {
public:
/**
Implements the slsReceiverData structure for the gotthard short read out by a module i.e. using the slsReceiver
(1x256 pixels, 1 packet 256 large etc.)
\param c crosstalk parameter for the output buffer
*/
gotthardShortModuleData(double c=0): slsReceiverData<uint16_t>(xpixels, ypixels, npackets, buffersize), xtalk(c){
uint16_t **dMask;
int **dMap;
int ix, iy;
int offset = 2;
dMask=new uint16_t*[ypixels];
dMap=new int*[ypixels];
for (int i = 0; i < ypixels; i++) {
dMap[i] = new int[xpixels];
dMask[i] = new uint16_t[xpixels];
}
for(ix=0; ix<ypixels; ++ix)
for(iy=0; iy<xpixels; ++iy)
dMask[ix][iy] = 0x0;
for(ix=0; ix<ypixels; ++ix)
for(iy=0; iy<xpixels; ++iy){
dMap[ix][iy] = offset;
offset++;
}
setDataMap(dMap);
setDataMask(dMask);
};
/**
Returns the frame number for the given dataset.
\param buff pointer to the dataset
\returns frame number
*/
int getFrameNumber(char *buff){
return (*(int*)buff);
};
/**
gets the packets number (last packet is labelled with 0 and is replaced with 40)
\param buff pointer to the memory
\returns packet number
*/
int getPacketNumber(char *buff){
return 1;
};
/**
returns the pixel value as double correcting for the output buffer crosstalk
\param data pointer to the memory
\param ix coordinate in the x direction
\param iy coordinate in the y direction
\returns channel value as double
*/
double getValue(char *data, int ix, int iy=0) {
//check how it is for gotthard
if (xtalk==0)
return slsDetectorData<uint16_t>::getValue(data, ix, iy);
else
return slsDetectorData<uint16_t>::getValue(data, ix, iy)-xtalk*slsDetectorData<uint16_t>::getValue(data, ix-1, iy);
};
/** sets the output buffer crosstalk correction parameter
\param c output buffer crosstalk correction parameter to be set
\returns current value for the output buffer crosstalk correction parameter
*/
double setXTalk(double c) {xtalk=c; return xtalk;}
/** gets the output buffer crosstalk parameter
\returns current value for the output buffer crosstalk correction parameter
*/
double getXTalk() {return xtalk;}
private:
double xtalk; /**<output buffer crosstalk correction parameter */
const static int xpixels = 256;
const static int ypixels = 1;
const static int npackets = 1;
const static int buffersize = 518;
};
#endif

156
slsReceiverSoftware/slsDetectorCalibration/jungfrau02Data.h
View File

@ -1,156 +0,0 @@
#ifndef JUNGFRAU02DATA_H
#define JUNGFRAU02DATA_H
#include "chiptestBoardData.h"
class jungfrau02Data : public chiptestBoardData {
public:
/**
Implements the chiptestBoardData structure for the jungfrau02 prototype
(48x48 pixels, ADC2 for analog output, 2 more ADCs used for readout of digital bits, pixel offset configurable.)
\param c crosstalk parameter for the output buffer
*/
jungfrau02Data(int nadc, int offset, double c=0): chiptestBoardData(48, 48, nadc, offset),
xtalk(c) {
uint16_t **dMask;
int **dMap;
dMask=new uint16_t*[48];
dMap=new int*[48];
for (int i = 0; i < 48; i++) {
dMap[i] = new int[48];
dMask[i] = new uint16_t[48];
}
for (int iy=0; iy<48; iy++) {
for (int ix=0; ix<48; ix++) {
dMap[ix][iy]=offset+(iy*48+ix)*nAdc;//dMap[iy][ix]=offset+(iy*48+ix)*nAdc;
// cout << ix << " " << iy << " " << dMap[ix][iy] << endl;
}
}
cout << (0,0) << " " << dMap[0][0] << endl;
cout << (47,47) << " " << dMap[47][47] << endl;
for (int ix=0; ix<48; ix++) {
for (int iy=0; iy<48; iy++)
dMask[iy][ix]=0x0;
setDataMap(dMap);
setDataMask(dMask);
}
};
/**
Returns the value of the selected channel for the given dataset. Since the ADC is only 14bit wide, only bit 0-13 are occupied. If the gain bits are read out, they are returned in bit 14-15.
\param data pointer to the dataset (including headers etc)
\param ix pixel number in the x direction
\param iy pixel number in the y direction
\returns data for the selected channel, with inversion if required
*/
virtual uint16_t getChannel(char *data, int ix, int iy=0) {
uint16_t m=0, d=0;
if (ix>=0 && ix<nx && iy>=0 && iy<ny && dataMap[iy][ix]>=0 && dataMap[iy][ix]<dataSize) {
m=dataMask[iy][ix];
d=*(((uint16_t*)(data))+dataMap[iy][ix]);
// cout << ix << " " << iy << " " << (uint16_t*)data << " " << ((uint16_t*)(data))+dataMap[iy][ix] << " " << d << endl;
if (nAdc==3) {
//cout << "Gain bit!" << endl;
if (*((uint16_t*)(data)+dataMap[iy][ix]+2)>8000) //exchange if gainBits==2 is returned!
d|=(1<<14); // gain bit 1
if (*((uint16_t*)(data)+dataMap[iy][ix]+1)>8000) //exchange if gainBits==2 is returned!
d|=(1<<15); // gain bit 0
}
}
return d^m;
};
/**
returns the pixel value as double correcting for the output buffer crosstalk
\param data pointer to the memory
\param ix coordinate in the x direction
\param iy coordinate in the y direction
\returns channel value as double
*/
double getValue(char *data, int ix, int iy=0) {
// cout << "##" << (void*)data << " " << ix << " " <<iy << endl;
uint16_t d=getChannel(data, ix, iy);
// cout << d << " " << (d&0x3fff) << endl;
if (xtalk==0 || ix==0)
return (double)(getChannel(data, ix, iy)&0x3fff);
else
return (double)(getChannel(data, ix, iy)&0x3fff)-xtalk* (double)(getChannel(data, ix, iy)&0x3fff);
};
/**
returns the gain bit value, i.e. returns 0 if(GB0==0 && GB1==0), returns 1 if(GB0==1 && GB1==0), returns 3 if(GB0==1 && GB1==1), if it returns 2 -> the gain bits are read out the wrong way around (i.e. GB0 and GB1 have to be reversed!)
\param data pointer to the memory
\param ix coordinate in the x direction
\param iy coordinate in the y direction
\returns gain bits as int
*/
int getGainBits(char *data, int ix, int iy=0) {
uint16_t d=getChannel(data, ix, iy);
return ((d&0xc000)>>14);
};
//*/
/** sets the output buffer crosstalk correction parameter
\param c output buffer crosstalk correction parameter to be set
\returns current value for the output buffer crosstalk correction parameter
*/
double setXTalk(double c) {xtalk=c; return xtalk;}
/** gets the output buffer crosstalk parameter
\returns current value for the output buffer crosstalk correction parameter
*/
double getXTalk() {return xtalk;}
private:
double xtalk; /**<output buffer crosstalk correction parameter */
};
#endif

264
slsReceiverSoftware/slsDetectorCalibration/jungfrauReadData.C
View File

@ -1,264 +0,0 @@
#include <TH1D.h>
#include <TH2D.h>
#include <TPad.h>
#include <TDirectory.h>
#include <TEntryList.h>
#include <TFile.h>
#include <TMath.h>
#include <TTree.h>
#include <TChain.h>
#include <THStack.h>
#include <TCanvas.h>
#include <stdio.h>
//#include <deque>
//#include <list>
//#include <queue>
#include <fstream>
#include "jungfrau02Data.h"
#include "jungfrau02CommonMode.h"
#include "singlePhotonDetector.h"
//#include "MovingStat.h"
using namespace std;
#define NC 48
#define NR 48
#define MY_DEBUG 1
#ifdef MY_DEBUG
#include <TCanvas.h>
#endif
/**
Loops over data file to find single photons, fills the tree (and writes it to file, although the root file should be opened before) and creates 1x1, 2x2, 3x3 cluster histograms with ADCu on the x axis, channel number (48*x+y) on the y axis.
\param fformat file name format
\param tit title of the tree etc.
\param runmin minimum run number
\param runmax max run number
\param nbins number of bins for spectrum hists
\param hmin histo minimum for spectrum hists
\param hmax histo maximum for spectrum hists
\param sign sign of the spectrum to find hits
\param hc readout correlation coefficient with previous pixel
\param xmin minimum x coordinate
\param xmax maximum x coordinate
\param ymin minimum y coordinate
\param ymax maximum y coordinate
\param cmsub enable commonmode subtraction
\param hitfinder if 0: performs pedestal subtraction, not hit finding; if 1: performs both pedestal subtraction and hit finding
\returns pointer to histo stack with cluster spectra
*/
THStack *jungfrauReadData(char *fformat, char *tit, int runmin, int runmax, int nbins=1500, int hmin=-500, int hmax=1000, int sign=1, double hc=0, int xmin=1, int xmax=NC-1, int ymin=1, int ymax=NR-1, int cmsub=0, int hitfinder=1){
/*
// read in calibration file
ifstream calfile("/home/l_msdetect/TriesteBeam2014/dummy data for scripts/CalibrationParametersTest_fake.txt");
if (calfile.is_open()==0){cout << "Unable to open calibration file!" << endl;}
int pix;
double of0,sl0,of1,sl1,of2,sl2;
double of_0[NC*NR], of_1[NC*NR], of_2[NC*NR], sl_0[NC*NR], sl_1[NC*NR], sl_2[NC*NR];
while (calfile >> pix >> of0 >> sl0 >> of1 >> sl1 >> of2 >> sl2){
of_0[pix]=of0;
sl_0[pix]=sl0;
of_1[pix]=of1;
sl_1[pix]=sl1;
of_2[pix]=of2;
sl_2[pix]=sl2; //if(pix==200) cout << "sl_2[200] " << sl_2[200] << endl;
}
calfile.close();
*/
double adc_value, num_photon;
jungfrau02Data *decoder=new jungfrau02Data(1,0,0);//(1,0,0); // (adc,offset,crosstalk) //(1,0,0) //(3,0,0) for readout of GB
jungfrau02CommonMode *cmSub=NULL;
if (cmsub)
cmSub=new jungfrau02CommonMode();
int nph=0;
int iev=0;
singlePhotonDetector<uint16_t> *filter=new singlePhotonDetector<uint16_t>(decoder, 3, 5, sign, cmSub);
char *buff;
char fname[10000];
int nf=0;
eventType thisEvent=PEDESTAL;
// int iframe;
// double *data, ped, sigma;
// data=decoder->getCluster();
TH2F *h2;
TH2F *h3;
TH2F *hetaX; // not needed for JF?
TH2F *hetaY; // not needed for JF?
THStack *hs=new THStack("hs",fformat);
TH2F *h1=new TH2F("h1",tit,nbins,hmin-0.5,hmax-0.5,NC*NR,-0.5,NC*NR-0.5);
hs->Add(h1);
if (hitfinder) {
h2=new TH2F("h2",tit,nbins,hmin-0.5,hmax-0.5,NC*NR,-0.5,NC*NR-0.5);
h3=new TH2F("h3",tit,nbins,hmin-0.5,hmax-0.5,NC*NR,-0.5,NC*NR-0.5);
hetaX=new TH2F("hetaX",tit,nbins,-1,2,NC*NR,-0.5,NC*NR-0.5); // not needed for JF?
hetaY=new TH2F("hetaY",tit,nbins,-1,2,NC*NR,-0.5,NC*NR-0.5); // not needed for JF?
hs->Add(h2);
hs->Add(h3);
hs->Add(hetaX); // not needed for JF?
hs->Add(hetaY); // not needed for JF?
}
ifstream filebin;
int ix=20, iy=20, ir, ic;
Int_t iFrame;
TTree *tall;
if (hitfinder)
tall=filter->initEventTree(tit, &iFrame);
#ifdef MY_DEBUG
TCanvas *myC;
TH2F *he;
TCanvas *cH1;
TCanvas *cH2;
TCanvas *cH3;
if (hitfinder) {
myC=new TCanvas();
he=new TH2F("he","Event Mask",xmax-xmin, xmin, xmax, ymax-ymin, ymin, ymax);
he->SetStats(kFALSE);
he->Draw("colz");
cH1=new TCanvas();
cH1->SetLogz();
h1->Draw("colz");
cH2=new TCanvas();
cH2->SetLogz();
h2->Draw("colz");
cH3=new TCanvas();
cH3->SetLogz();
h3->Draw("colz");
}
#endif
filter->newDataSet();
for (int irun=runmin; irun<=runmax; irun++){
sprintf(fname,fformat,irun);
cout << "file name " << fname << endl;
filebin.open((const char *)(fname), ios::in | ios::binary);
nph=0;
while ((buff=decoder->readNextFrame(filebin))) {
if (hitfinder) {
filter->newFrame();
//calculate pedestals and common modes
if (cmsub) {
// cout << "cm" << endl;
for (ix=xmin-1; ix<xmax+1; ix++)
for (iy=ymin-1; iy<ymax+1; iy++) {
thisEvent=filter->getEventType(buff, ix, iy,0);
}
}
}
// cout << "new frame " << endl;
for (ix=xmin-1; ix<xmax+1; ix++)
for (iy=ymin-1; iy<ymax+1; iy++) {
// cout << ix << " " << iy << endl;
thisEvent=filter->getEventType(buff, ix, iy, cmsub);
int gainBits=decoder->getGainBits(buff,ix,iy); //XXX
#ifdef MY_DEBUG
if (hitfinder) {
if (iev%1000==0)
//he->SetBinContent(ix+1-xmin, iy+1-ymin, (int)thisEvent); // show single photon hits // original
//he->SetBinContent(ix+1-xmin, iy+1-ymin, cmSub->getCommonMode(ix,iy)); //show common mode!
he->SetBinContent(iy+1-ymin, ix+1-xmin, (int)gainBits); //show gain bits
//he->SetBinContent(ix+1-xmin, iy+1-ymin, (int)gainBits); // rows and columns reversed!!!
}
#endif
if (nf>1000) { // only start filing data after 1000 frames
// h1->Fill(decoder->getValue(buff,ix,iy), iy+NR*ix);
h1->Fill(filter->getClusterTotal(1), iy+NR*ix);
if (hitfinder) {
if (thisEvent==PHOTON_MAX ) {
nph++;
h2->Fill(filter->getClusterTotal(2), iy+NR*ix);
h3->Fill(filter->getClusterTotal(3), iy+NR*ix);
iFrame=decoder->getFrameNumber(buff);
tall->Fill();
}
}
}
}
//////////////////////////////////////////////////////////
#ifdef MY_DEBUG
if (iev%1000==0) {
myC->Modified();
myC->Update();
cH1->Modified();
cH1->Update();
cH2->Modified();
cH2->Update();
cH3->Modified();
cH3->Update();
}
iev++;
#endif
nf++;
cout << "="; // one "=" for every frame that's been processed
delete [] buff;
}
cout << nph << endl; // number of photons found in file
if (filebin.is_open())
filebin.close();
else
cout << "Could not open file :( ... " << fname << endl;
}
if (hitfinder)
tall->Write(tall->GetName(),TObject::kOverwrite);
delete decoder;
cout << "Read " << nf << " frames." << endl;
return hs;
}

137
slsReceiverSoftware/slsDetectorCalibration/moench02ModuleData.h
View File

@ -1,137 +0,0 @@
#ifndef MOENCH02MODULEDATA_H
#define MOENCH02MODULEDATA_H
#include "slsReceiverData.h"
class moench02ModuleData : public slsReceiverData<uint16_t> {
public:
/**
Implements the slsReceiverData structure for the moench02 prototype read out by a module i.e. using the slsReceiver
(160x160 pixels, 40 packets 1286 large etc.)
\param c crosstalk parameter for the output buffer
*/
moench02ModuleData(double c=0): slsReceiverData<uint16_t>(160, 160, 40, 1286),
xtalk(c) {
uint16_t **dMask;
int **dMap;
int ix, iy;
dMask=new uint16_t*[160];
dMap=new int*[160];
for (int i = 0; i < 160; i++) {
dMap[i] = new int[160];
dMask[i] = new uint16_t[160];
}
for (int isc=0; isc<4; isc++) {
for (int ip=0; ip<10; ip++) {
for (int ir=0; ir<16; ir++) {
for (int ic=0; ic<40; ic++) {
ix=isc*40+ic;
iy=ip*16+ir;
dMap[iy][ix]=1286*(isc*10+ip)+2*ir*40+2*ic+4;
// cout << ix << " " << iy << " " << dMap[ix][iy] << endl;
}
}
}
}
for (ix=0; ix<120; ix++) {
for (iy=0; iy<160; iy++)
dMask[iy][ix]=0x3fff;
}
for (ix=120; ix<160; ix++) {
for (iy=0; iy<160; iy++)
dMask[iy][ix]=0x0;
}
setDataMap(dMap);
setDataMask(dMask);
};
/**
gets the packets number (last packet is labelled with 0 and is replaced with 40)
\param buff pointer to the memory
\returns packet number
*/
int getPacketNumber(char *buff){
int np=(*(int*)buff)&0xff;
if (np==0)
np=40;
return np;
};
/**
returns the pixel value as double correcting for the output buffer crosstalk
\param data pointer to the memory
\param ix coordinate in the x direction
\param iy coordinate in the y direction
\returns channel value as double
*/
double getValue(char *data, int ix, int iy=0) {
// cout << "##" << (void*)data << " " << ix << " " <<iy << endl;
if (xtalk==0 || ix%40==0)
return slsDetectorData<uint16_t>::getValue(data, ix, iy);
else
return slsDetectorData<uint16_t>::getValue(data, ix, iy)-xtalk*slsDetectorData<uint16_t>::getValue(data, ix-1, iy);
};
/** sets the output buffer crosstalk correction parameter
\param c output buffer crosstalk correction parameter to be set
\returns current value for the output buffer crosstalk correction parameter
*/
double setXTalk(double c) {xtalk=c; return xtalk;}
/** gets the output buffer crosstalk parameter
\returns current value for the output buffer crosstalk correction parameter
*/
double getXTalk() {return xtalk;}
private:
double xtalk; /**<output buffer crosstalk correction parameter */
};
#endif

45
slsReceiverSoftware/slsDetectorCalibration/moenchCommonMode.h
View File

@ -1,45 +0,0 @@
#ifndef MOENCHCOMMONMODE_H
#define MOENCHCOMMONMODE_H
#include "commonModeSubtraction.h"
class moenchCommonMode : public commonModeSubtraction {
/** @short class to calculate the common mode noise for moench02 i.e. on 4 supercolumns separately */
public:
/** constructor - initalizes a commonModeSubtraction with 4 different regions of interest
\param nn number of samples for the moving average
*/
moenchCommonMode(int nn=1000) : commonModeSubtraction(nn,4){} ;
/** add value to common mode as a function of the pixel value, subdividing the region of interest in the 4 supercolumns of 40 columns each;
\param val value to add to the common mode
\param ix pixel coordinate in the x direction
\param iy pixel coordinate in the y direction
*/
virtual void addToCommonMode(double val, int ix=0, int iy=0) {
(void) iy;
int isc=ix/40;
if (isc>=0 && isc<nROI) {
cmPed[isc]+=val;
nCm[isc]++;
}
};
/**returns common mode value as a function of the pixel value, subdividing the region of interest in the 4 supercolumns of 40 columns each;
\param ix pixel coordinate in the x direction
\param iy pixel coordinate in the y direction
\returns common mode value
*/
virtual double getCommonMode(int ix=0, int iy=0) {
(void) iy;
int isc=ix/40;
if (isc>=0 && isc<nROI) {
if (nCm[isc]>0) return cmPed[isc]/nCm[isc]-cmStat[isc].Mean();
}
return 0;
};
};
#endif

244
slsReceiverSoftware/slsDetectorCalibration/moenchMakeTree.C
View File

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

236
slsReceiverSoftware/slsDetectorCalibration/moenchReadData.C
View File

@ -1,236 +0,0 @@
#include <TH1D.h>
#include <TH2D.h>
#include <TPad.h>
#include <TDirectory.h>
#include <TEntryList.h>
#include <TFile.h>
#include <TMath.h>
#include <TTree.h>
#include <TChain.h>
#include <THStack.h>
#include <TCanvas.h>
#include <stdio.h>
//#include <deque>
//#include <list>
//#include <queue>
#include <fstream>
#include "moench02ModuleData.h"
#include "moenchCommonMode.h"
#include "singlePhotonDetector.h"
//#include "MovingStat.h"
using namespace std;
#define NC 160
#define NR 160
#define MY_DEBUG 1
#ifdef MY_DEBUG
#include <TCanvas.h>
#endif
/**
Loops over data file to find single photons, fills the tree (and writes it to file, althoug the root file should be opened before) and creates 1x1, 2x2, 3x3 cluster histograms with ADCu on the x axis, channel number (160*x+y) on the y axis.
\param fformat file name format
\param tit title of the tree etc.
\param runmin minimum run number
\param runmax max run number
\param nbins number of bins for spectrum hists
\param hmin histo minimum for spectrum hists
\param hmax histo maximum for spectrum hists
\param sign sign of the spectrum to find hits
\param hc readout correlation coefficient with previous pixel
\param xmin minimum x coordinate
\param xmax maximum x coordinate
\param ymin minimum y coordinate
\param ymax maximum y coordinate
\param cmsub enable commonmode subtraction
\returns pointer to histo stack with cluster spectra
*/
THStack *moenchReadData(char *fformat, char *tit, int runmin, int runmax, int nbins=1500, int hmin=-500, int hmax=1000, int sign=1, double hc=0, int xmin=1, int xmax=NC-1, int ymin=1, int ymax=NR-1, int cmsub=0, int hitfinder=1) {
moench02ModuleData *decoder=new moench02ModuleData(hc);
moenchCommonMode *cmSub=NULL;
if (cmsub)
cmSub=new moenchCommonMode();
int nph=0;
singlePhotonDetector<uint16_t> *filter=new singlePhotonDetector<uint16_t>(decoder, 3, 5, sign, cmSub);
char *buff;
char fname[10000];
int nf=0;
eventType thisEvent=PEDESTAL;
// int iframe;
// double *data, ped, sigma;
// data=decoder->getCluster();
TH2F *h2;
TH2F *h3;
TH2F *hetaX;
TH2F *hetaY;
THStack *hs=new THStack("hs",fformat);
TH2F *h1=new TH2F("h1",tit,nbins,hmin-0.5,hmax-0.5,NC*NR,-0.5,NC*NR-0.5);
hs->Add(h1);
if (hitfinder) {
h2=new TH2F("h2",tit,nbins,hmin-0.5,hmax-0.5,NC*NR,-0.5,NC*NR-0.5);
h3=new TH2F("h3",tit,nbins,hmin-0.5,hmax-0.5,NC*NR,-0.5,NC*NR-0.5);
hetaX=new TH2F("hetaX",tit,nbins,-1,2,NC*NR,-0.5,NC*NR-0.5);
hetaY=new TH2F("hetaY",tit,nbins,-1,2,NC*NR,-0.5,NC*NR-0.5);
hs->Add(h2);
hs->Add(h3);
hs->Add(hetaX);
hs->Add(hetaY);
}
ifstream filebin;
int ix=20, iy=20, ir, ic;
Int_t iFrame;
TTree *tall;
if (hitfinder)
tall=filter->initEventTree(tit, &iFrame);
#ifdef MY_DEBUG
TCanvas *myC;
TH2F *he;
TCanvas *cH1;
TCanvas *cH2;
TCanvas *cH3;
if (hitfinder) {
myC=new TCanvas();
he=new TH2F("he","Event Mask",xmax-xmin, xmin, xmax, ymax-ymin, ymin, ymax);
he->SetStats(kFALSE);
he->Draw("colz");
cH1=new TCanvas();
cH1->SetLogz();
h1->Draw("colz");
cH2=new TCanvas();
cH2->SetLogz();
h2->Draw("colz");
cH3=new TCanvas();
cH3->SetLogz();
h3->Draw("colz");
}
#endif
filter->newDataSet();
for (int irun=runmin; irun<runmax; irun++) {
sprintf(fname,fformat,irun);
cout << "file name " << fname << endl;
filebin.open((const char *)(fname), ios::in | ios::binary);
nph=0;
while ((buff=decoder->readNextFrame(filebin))) {
if (hitfinder) {
filter->newFrame();
//calculate pedestals and common modes
if (cmsub) {
// cout << "cm" << endl;
for (ix=xmin-1; ix<xmax+1; ix++)
for (iy=ymin-1; iy<ymax+1; iy++) {
thisEvent=filter->getEventType(buff, ix, iy,0);
}
}
}
// cout << "new frame " << endl;
for (ix=xmin-1; ix<xmax+1; ix++)
for (iy=ymin-1; iy<ymax+1; iy++) {
// cout << ix << " " << iy << endl;
thisEvent=filter->getEventType(buff, ix, iy, cmsub);
#ifdef MY_DEBUG
if (hitfinder) {
if (iev%1000==0)
he->SetBinContent(ix+1-xmin, iy+1-ymin, (int)thisEvent);
}
#endif
if (nf>1000) {
h1->Fill(filter->getClusterTotal(1), iy+NR*ix);
if (hitfinder) {
if (thisEvent==PHOTON_MAX ) {
nph++;
h2->Fill(filter->getClusterTotal(2), iy+NR*ix);
h3->Fill(filter->getClusterTotal(3), iy+NR*ix);
iFrame=decoder->getFrameNumber(buff);
tall->Fill();
}
}
}
}
//////////////////////////////////////////////////////////
#ifdef MY_DEBUG
if (iev%1000==0) {
myC->Modified();
myC->Update();
cH1->Modified();
cH1->Update();
cH2->Modified();
cH2->Update();
cH3->Modified();
cH3->Update();
}
iev++;
#endif
nf++;
cout << "=" ;
delete [] buff;
}
cout << nph << endl;
if (filebin.is_open())
filebin.close();
else
cout << "could not open file " << fname << endl;
}
if (hitfinder)
tall->Write(tall->GetName(),TObject::kOverwrite);
delete decoder;
cout << "Read " << nf << " frames" << endl;
return hs;
}

52
slsReceiverSoftware/slsDetectorCalibration/moenchReadDataMT.C
View File

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

48
slsReceiverSoftware/slsDetectorCalibration/pedestalSubtraction.h
View File

@ -1,48 +0,0 @@
#ifndef PEDESTALSUBTRACTION_H
#define PEDESTALSUBTRACTION_H
#include "MovingStat.h"
class pedestalSubtraction {
/** @short class defining the pedestal subtraction based on an approximated moving average */
public:
/** constructor
\param nn number of samples to calculate the moving average (defaults to 1000)
*/
pedestalSubtraction (int nn=1000) : stat(nn) {};
/** virtual destructorr
*/
virtual ~pedestalSubtraction() {};
/** clears the moving average */
virtual void Clear() {stat.Clear();}
/** adds the element to the moving average
\param val value to be added
*/
virtual void addToPedestal(double val){stat.Calc(val);};
/** returns the average value of the pedestal
\returns mean of the moving average
*/
virtual double getPedestal(){return stat.Mean();};
/** returns the standard deviation of the moving average
\returns standard deviation of the moving average
*/
virtual double getPedestalRMS(){return stat.StandardDeviation();};
/**sets/gets the number of samples for the moving average
\param i number of elements for the moving average. If -1 (default) or negative, gets.
\returns actual number of samples for the moving average
*/
virtual int SetNPedestals(int i=-1) {if (i>0) stat.SetN(i); return stat.GetN();};
private:
MovingStat stat; /**< approximated moving average struct */
};
#endif

136
slsReceiverSoftware/slsDetectorCalibration/raedNoiseData.C
View File

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

31
slsReceiverSoftware/slsDetectorCalibration/readJungfrauData.C
View File

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

387
slsReceiverSoftware/slsDetectorCalibration/singlePhotonDetector.h
View File

@ -1,387 +0,0 @@
#ifndef SINGLEPHOTONDETECTOR_H
#define SINGLEPHOTONDETECTOR_H
#include "slsDetectorData.h"
#include "single_photon_hit.h"
#include "pedestalSubtraction.h"
#include "commonModeSubtraction.h"
//#define MYROOT1
#ifdef MYROOT1
#include <TTree.h>
#endif
#include <iostream>
using namespace std;
enum eventType {
PEDESTAL=0,
NEIGHBOUR=1,
PHOTON=2,
PHOTON_MAX=3,
NEGATIVE_PEDESTAL=4,
UNDEFINED_EVENT=-1
};
enum quadrant {
TOP_LEFT=0,
TOP_RIGHT=1,
BOTTOM_LEFT=2,
BOTTOM_RIGHT=3,
UNDEFINED_QUADRANT=-1
};
template <class dataType>
class singlePhotonDetector {
/** @short class to perform pedestal subtraction etc. and find single photon clusters for an analog detector */
public:
/**
Constructor (no error checking if datasize and offsets are compatible!)
\param d detector data structure to be used
\param csize cluster size (should be an odd number). Defaults to 3
\param nsigma number of rms to discriminate from the noise. Defaults to 5
\param sign 1 if photons are positive, -1 if negative
\param cm common mode subtraction algorithm, if any. Defaults to NULL i.e. none
\param nped number of samples for pedestal averaging
\param nd number of dark frames to average as pedestals without photon discrimination at the beginning of the measurement
*/
singlePhotonDetector(slsDetectorData<dataType> *d,
int csize=3,
double nsigma=5,
int sign=1,
commonModeSubtraction *cm=NULL,
int nped=1000,
int nd=100) : det(d), nx(0), ny(0), stat(NULL), cmSub(cm), nDark(nd), eventMask(NULL),nSigma (nsigma), clusterSize(csize), clusterSizeY(csize), cluster(NULL), iframe(-1), dataSign(sign), quad(UNDEFINED_QUADRANT), tot(0), quadTot(0) {
det->getDetectorSize(nx,ny);
stat=new pedestalSubtraction*[ny];
eventMask=new eventType*[ny];
for (int i=0; i<ny; i++) {
stat[i]=new pedestalSubtraction[nx];
stat[i]->SetNPedestals(nped);
eventMask[i]=new eventType[nx];
}
if (ny==1)
clusterSizeY=1;
cluster=new single_photon_hit(clusterSize,clusterSizeY);
setClusterSize(csize);
};
/**
destructor. Deletes the cluster structure and the pdestalSubtraction array
*/
virtual ~singlePhotonDetector() {delete cluster; for (int i=0; i<ny; i++) delete [] stat[i]; delete [] stat;};
/** resets the pedestalSubtraction array and the commonModeSubtraction */
void newDataSet(){iframe=-1; for (int iy=0; iy<ny; iy++) for (int ix=0; ix<nx; ix++) stat[iy][ix].Clear(); if (cmSub) cmSub->Clear(); };
/** resets the eventMask to undefined and the commonModeSubtraction */
void newFrame(){iframe++; for (int iy=0; iy<ny; iy++) for (int ix=0; ix<nx; ix++) eventMask[iy][ix]=UNDEFINED_EVENT; if (cmSub) cmSub->newFrame();};
/** 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;};
/**
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
*/
virtual void addToPedestal(double val, int ix, int iy){
// cout << "*"<< ix << " " << iy << " " << val << endl;
if (ix>=0 && ix<nx && iy>=0 && iy<ny) {
// cout << ix << " " << iy << " " << val << endl;
stat[iy][ix].addToPedestal(val);
if (cmSub && det->isGood(ix, iy) )
cmSub->addToCommonMode(val, ix, iy);
};
};
/**
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)
*/
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;};
/**
gets pedestal rms (i.e. noise)
\param ix pixel x coordinate
\param iy pixel y coordinate
*/
double getPedestalRMS(int ix, int iy){if (ix>=0 && ix<nx && iy>=0 && iy<ny) return stat[iy][ix].getPedestalRMS();else return -1;};
/** sets/gets number of rms threshold to detect photons
\param n number of sigma to be set (0 or negative gets)
\returns actual number of sigma parameter
*/
double setNSigma(double n=-1){if (n>0) nSigma=n; return nSigma;}
/** sets/gets cluster size
\param n cluster size to be set, (0 or negative gets). If even is incremented by 1.
\returns actual cluster size
*/
int setClusterSize(int n=-1){
if (n>0 && n!=clusterSize) {
if (n%2==0)
n+=1;
clusterSize=n;
if (cluster)
delete cluster;
if (ny>1)
clusterSizeY=clusterSize;
cluster=new single_photon_hit(clusterSize,clusterSizeY);
}
return clusterSize;
}
/** finds event type for pixel and fills cluster structure. The algorithm loops only if the evenMask for this pixel is still undefined.
if pixel or cluster around it are above threshold (nsigma*pedestalRMS) cluster is filled and pixel mask is PHOTON_MAX (if maximum in cluster) or NEIGHBOUR; If PHOTON_MAX, the elements of the cluster are also set as NEIGHBOURs in order to speed up the looping
if below threshold the pixel is either marked as PEDESTAL (and added to the pedestal calculator) or NEGATIVE_PEDESTAL is case it's lower than -threshold, otherwise the pedestal average would drift to negative values while it should be 0.
/param data pointer to the data
/param ix pixel x coordinate
/param iy pixel y coordinate
/param cm enable(1)/disable(0) common mode subtraction (if defined).
/returns event type for the given pixel
*/
eventType getEventType(char *data, int ix, int iy, int cm=0) {
// eventType ret=PEDESTAL;
double max=0, tl=0, tr=0, bl=0,br=0, v;
// cout << iframe << endl;
tot=0;
quadTot=0;
quad=UNDEFINED_QUADRANT;
if (iframe<nDark) {
if (cm==0) {
// cout << "=" << endl;
addToPedestal(det->getValue(data, ix, iy),ix,iy);
cluster->set_data(dataSign*(det->getValue(data, ix, iy)-getPedestal(ix,iy,cm)), 0,0 );
// cout << "=" << endl;
}
return UNDEFINED_EVENT;
}
// if (eventMask[iy][ix]==UNDEFINED) {
eventMask[iy][ix]=PEDESTAL;
cluster->x=ix;
cluster->y=iy;
cluster->rms=getPedestalRMS(ix,iy);
cluster->ped=getPedestal(ix,iy, cm);
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(dataSign*(det->getValue(data, ix+ic, iy+ir)-getPedestal(ix+ic,iy+ir,cm)), ic, ir );
v=cluster->get_data(ic,ir);
tot+=v;
if (ir<=0 && ic<=0)
bl+=v;
if (ir<=0 && ic>=0)
br+=v;
if (ir>=0 && ic<=0)
tl+=v;
if (ir>=0 && ic>=0)
tr+=v;
if (cluster->get_data(ic,ir)>max) {
max=v;
}
if (ir==0 && ic==0) {
if (v>nSigma*cluster->rms) {
eventMask[iy][ix]=PHOTON;
} else if (cluster->get_data(ic,ir)<-nSigma*cluster->rms)
eventMask[iy][ix]=NEGATIVE_PEDESTAL;
}
}
}
}
if (eventMask[iy][ix]!=PHOTON && tot>sqrt(clusterSizeY*clusterSize)*nSigma*cluster->rms) {
eventMask[iy][ix]=NEIGHBOUR;
} else if (eventMask[iy][ix]==PHOTON) {
if (cluster->get_data(0,0)>=max) {
eventMask[iy][ix]=PHOTON_MAX;
}
} else if (eventMask[iy][ix]==PEDESTAL) {
if (cm==0)
addToPedestal(det->getValue(data, ix, iy),ix,iy);
}
if (bl>=br && bl>=tl && bl>=tr) {
quad=BOTTOM_LEFT;
quadTot=bl;
} else if (br>=bl && br>=tl && br>=tr) {
quad=BOTTOM_RIGHT;
quadTot=br;
} else if (tl>=br && tl>=bl && tl>=tr) {
quad=TOP_LEFT;
quadTot=tl;
} else if (tr>=bl && tr>=tl && tr>=br) {
quad=TOP_RIGHT;
quadTot=tr;
}
return eventMask[iy][ix];
};
/**<
retrurns the total signal in a cluster
\param size cluser size should be 1,2 or 3
\returns cluster center if size=1, sum of the maximum quadrant if size=2, total of the cluster if size=3 or anything else
*/
double getClusterTotal(int size) {
switch (size) {
case 1:
return getClusterElement(0,0);
case 2:
return quadTot;
default:
return tot;
};
};
/**<
retrurns the quadrant with maximum signal
\returns quadrant where the cluster is located */
quadrant getQuadrant() {return quad;};
/** 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) {int ix=0, iy=0; if (i>0) for (ix=0; ix<nx; ix++) for (iy=0; iy<ny; iy++) stat[iy][ix].SetNPedestals(i); return stat[0][0].SetNPedestals();};
/** returns value for cluster element in relative coordinates
\param ic x coordinate (center is (0,0))
\param ir y coordinate (center is (0,0))
\returns cluster element
*/
double getClusterElement(int ic, int ir=0){return cluster->get_data(ic,ir);};
/** returns event mask for the given pixel
\param ic x coordinate (center is (0,0))
\param ir y coordinate (center is (0,0))
\returns event mask enum for the given pixel
*/
eventType getEventMask(int ic, int ir=0){return eventMask[ir][ic];};
#ifdef MYROOT1
/** generates a tree and maps the branches
\param tname name for the tree
\param iFrame pointer to the frame number
\returns returns pointer to the TTree
*/
TTree *initEventTree(char *tname, int *iFrame=NULL) {
TTree* tall=new TTree(tname,tname);
if (iFrame)
tall->Branch("iFrame",iFrame,"iframe/I");
else
tall->Branch("iFrame",&(cluster->iframe),"iframe/I");
tall->Branch("x",&(cluster->x),"x/I");
tall->Branch("y",&(cluster->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");
return tall;
};
#else
/** write cluster to filer*/
void writeCluster(FILE* myFile){cluster->write(myFile);};
#endif
private:
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 */
commonModeSubtraction *cmSub;/**< commonModeSubtraction class */
int nDark; /**< number of frames to be used at the beginning of the dataset to calculate pedestal without applying photon discrimination */
eventType **eventMask; /**< matrix of event type or each pixel */
double nSigma; /**< number of sigma parameter for photon discrimination */
int clusterSize; /**< cluster size in the x direction */
int clusterSizeY; /**< cluster size in the y direction i.e. 1 for strips, clusterSize for pixels */
single_photon_hit *cluster; /**< single photon hit data structure */
int iframe; /**< frame number (not from file but incremented within the dataset every time newFrame is called */
int dataSign; /**< sign of the data i.e. 1 if photon is positive, -1 if negative */
quadrant quad; /**< quadrant where the photon is located */
double tot; /**< sum of the 3x3 cluster */
double quadTot; /**< sum of the maximum 2x2cluster */
};
#endif

62
slsReceiverSoftware/slsDetectorCalibration/single_photon_hit.h
View File

@ -1,62 +0,0 @@
#ifndef SINGLE_PHOTON_HIT_H
#define SINGLE_PHOTON_HIT_h
typedef double double32_t;
typedef float float32_t;
typedef int int32_t;
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, int ny=1): dx(nx), dy(ny) {data=new double[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
*/
void write(FILE *myFile) {fwrite((void*)this, 1, 3*sizeof(int)+2*sizeof(double), myFile); fwrite((void*)data, 1, dx*dy*sizeof(double), myFile);};
/**
binary read from file of all elements of the structure, except size of the cluster. The structure is then filled with those args
\param myFile file descriptor
*/
void read(FILE *myFile) {fread((void*)this, 1, 3*sizeof(int)+2*sizeof(double), myFile); fread((void*)data, 1, dx*dy*sizeof(double), myFile);};
/**
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;};
/**
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 x; /**< x-coordinate of the center of hit */
int y; /**< x-coordinate of the center of hit */
double rms; /**< noise of central pixel l -- at some point it can be removed*/
double ped; /**< pedestal of the central pixel -- at some point it can be removed*/
int iframe; /**< frame number */
double *data; /**< pointer to data */
const int dx; /**< size of data cluster in x */
const int dy; /**< size of data cluster in y */
};
#endif

251
slsReceiverSoftware/slsDetectorCalibration/slsDetectorData.h
View File

@ -1,251 +0,0 @@
#ifndef SLSDETECTORDATA_H
#define SLSDETECTORDATA_H
#include <inttypes.h>
#include <iostream>
#include <fstream>
using namespace std;
template <class dataType>
class slsDetectorData {
public:
/**
General slsDetectors data structure. Works for data acquired using the slsDetectorReceiver. Can be generalized to other detectors (many virtual funcs).
Constructor (no error checking if datasize and offsets are compatible!)
\param npx number of pixels in the x direction
\param npy number of pixels in the y direction (1 for strips)
\param dsize size of the data
\param dMap array of size nx*ny storing the pointers to the data in the dataset (as offset)
\param dMask Array of size nx*ny storing the polarity of the data in the dataset (should be 0 if no inversion is required, 0xffffffff is inversion is required)
\param dROI Array of size nx*ny. The elements are 1s if the channel is good or in the ROI, 0 is bad or out of the ROI. NULL (default) means all 1s.
*/
slsDetectorData(int npx, int npy, int dsize, int **dMap=NULL, dataType **dMask=NULL, int **dROI=NULL): nx(npx), ny(npy), dataSize(dsize) {
dataMask=new dataType*[ny];
for(int i = 0; i < ny; i++) {
dataMask[i] = new dataType[nx];
}
dataMap=new int*[ny];
for(int i = 0; i < ny; i++) {
dataMap[i] = new int[nx];
}
dataROIMask=new int*[ny];
for(int i = 0; i < ny; i++) {
dataROIMask[i] = new int[nx];
for (int j=0; j<nx; j++)
dataROIMask[i][j]=1;
}
setDataMap(dMap);
setDataMask(dMask);
setDataROIMask(dROI);
};
virtual ~slsDetectorData() {
for(int i = 0; i < ny; i++) {
delete [] dataMap[i];
delete [] dataMask[i];
delete [] dataROIMask[i];
}
delete [] dataMap;
delete [] dataMask;
delete [] dataROIMask;
}
/**
defines the data map (as offset) - no error checking if datasize and offsets are compatible!
\param dMap array of size nx*ny storing the pointers to the data in the dataset (as offset). If NULL (default),the data are arranged as if read out row by row (dataMap[iy][ix]=(iy*nx+ix)*sizeof(dataType);)
*/
void setDataMap(int **dMap=NULL) {
if (dMap==NULL) {
for (int iy=0; iy<ny; iy++)
for (int ix=0; ix<nx; ix++)
dataMap[iy][ix]=(iy*nx+ix)*sizeof(dataType);
} else {
for (int iy=0; iy<ny; iy++)
for (int ix=0; ix<nx; ix++)
dataMap[iy][ix]=dMap[iy][ix];
}
};
/**
defines the data mask i.e. the polarity of the data
\param dMask Array of size nx*ny storing the polarity of the data in the dataset (should be 0 if no inversion is required, 0xffffffff is inversion is required)
*/
void setDataMask(dataType **dMask=NULL){
if (dMask!=NULL) {
for (int iy=0; iy<ny; iy++)
for (int ix=0; ix<nx; ix++)
dataMask[iy][ix]=dMask[iy][ix];
} else {
for (int iy=0; iy<ny; iy++)
for (int ix=0; ix<nx; ix++)
dataMask[iy][ix]=0;
}
};
/**
defines the region of interest and/or the bad channels mask
\param dROI Array of size nx*ny. The lements are 1s if the channel is good or in the ROI, 0 is bad or out of the ROI. NULL (default) means all 1s.
*/
void setDataROIMask(int **dROI=NULL){
if (dROI!=NULL) {
for (int iy=0; iy<ny; iy++)
for (int ix=0; ix<nx; ix++)
dataROIMask[iy][ix]=dROI[iy][ix];
} else {
for (int iy=0; iy<ny; iy++)
for (int ix=0; ix<nx; ix++)
dataROIMask[iy][ix]=1;
}
};
/**
Define bad channel or roi mask for a single channel
\param ix channel x coordinate
\param iy channel y coordinate (1 for strips)
\param i 1 if pixel is good (or in the roi), 0 if bad
\returns 1 if pixel is good, 0 if it's bad, -1 if pixel is out of range
*/
int setGood(int ix, int iy, int i=1) { if (ix>=0 && ix<nx && iy>=0 && iy<ny) dataROIMask[iy][ix]=i; return isGood(ix,iy);};
/**
Define bad channel or roi mask for a single channel
\param ix channel x coordinate
\param iy channel y coordinate (1 for strips)
\returns 1 if pixel is good, 0 if it's bad, -1 if pixel is out of range
*/
int isGood(int ix, int iy) { if (ix>=0 && ix<nx && iy>=0 && iy<ny) return dataROIMask[iy][ix]; else return -1;};
/**
Returns detector size in x,y
\param npx reference to number of channels in x
\param npy reference to number of channels in y (will be 1 for strips)
\returns total number of channels
*/
int getDetectorSize(int &npx, int &npy){npx=nx; npy=ny; return nx*ny;};
/** Returns the size of the data frame */
int getDataSize() {return dataSize;};
/** changes the size of the data frame */
int setDataSize(int d) {dataSize=d; return dataSize;};
/**
Returns the value of the selected channel for the given dataset. Virtual function, can be overloaded.
\param data pointer to the dataset (including headers etc)
\param ix pixel number in the x direction
\param iy pixel number in the y direction
\returns data for the selected channel, with inversion if required
*/
virtual dataType getChannel(char *data, int ix, int iy=0) {
dataType m=0, d=0;
if (ix>=0 && ix<nx && iy>=0 && iy<ny && dataMap[iy][ix]>=0 && dataMap[iy][ix]<dataSize) {
m=dataMask[iy][ix];
d=*((dataType*)(data+dataMap[iy][ix]));
}
return d^m;
};
/**
Returns the value of the selected channel for the given dataset as double.
\param data pointer to the dataset (including headers etc)
\param ix pixel number in the x direction
\param iy pixel number in the y direction
\returns data for the selected channel, with inversion if required as double
*/
virtual double getValue(char *data, int ix, int iy=0) {return (double)getChannel(data, ix, iy);};
/**
Returns the frame number for the given dataset. Purely virtual func.
\param buff pointer to the dataset
\returns frame number
*/
virtual int getFrameNumber(char *buff)=0;
/**
Returns the packet number for the given dataset. purely virtual func
\param buff pointer to the dataset
\returns packet number number
virtual int getPacketNumber(char *buff)=0;
*/
/**
Loops over a memory slot until a complete frame is found (i.e. all packets 0 to nPackets, same frame number). purely virtual func
\param data pointer to the memory to be analyzed
\param ndata reference to the amount of data found for the frame, in case the frame is incomplete at the end of the memory slot
\param dsize size of the memory slot to be analyzed
\returns pointer to the beginning of the last good frame (might be incomplete if ndata smaller than dataSize), or NULL if no frame is found
*/
virtual char *findNextFrame(char *data, int &ndata, int dsize)=0;
/**
Loops over a file stream until a complete frame is found (i.e. all packets 0 to nPackets, same frame number). Can be overloaded for different kind of detectors!
\param filebin input file stream (binary)
\returns pointer to the begin of the last good frame, NULL if no frame is found or last frame is incomplete
*/
virtual char *readNextFrame(ifstream &filebin)=0;
protected:
const int nx; /**< Number of pixels in the x direction */
const int ny; /**< Number of pixels in the y direction */
int dataSize; /**<size of the data constituting one frame */
int **dataMap; /**< Array of size nx*ny storing the pointers to the data in the dataset (as offset)*/
dataType **dataMask; /**< Array of size nx*ny storing the polarity of the data in the dataset (should be 0 if no inversion is required, 0xffffffff is inversion is required) */
int **dataROIMask; /**< Array of size nx*ny 1 if channel is good (or in the ROI), 0 if bad channel (or out of ROI) */
};
#endif

171
slsReceiverSoftware/slsDetectorCalibration/slsReceiverData.h
View File

@ -1,171 +0,0 @@
#ifndef SLSRECEIVERDATA_H
#define SLSRECEIVERDATA_H
#include "slsDetectorData.h"
template <class dataType>
class slsReceiverData : public slsDetectorData<dataType> {
public:
/**
slsReceiver data structure. Works for data acquired using the slsDetectorReceiver subdivided in different packets with headers and footers.
Inherits and implements slsDetectorData.
Constructor (no error checking if datasize and offsets are compatible!)
\param npx number of pixels in the x direction
\param npy number of pixels in the y direction (1 for strips)
\param np number of packets
\param psize packets size
\param dMap array of size nx*ny storing the pointers to the data in the dataset (as offset)
\param dMask Array of size nx*ny storing the polarity of the data in the dataset (should be 0 if no inversion is required, 0xffffffff is inversion is required)
\param dROI Array of size nx*ny. The elements are 1s if the channel is good or in the ROI, 0 is bad or out of the ROI. NULL (default) means all 1s.
*/
slsReceiverData(int npx, int npy, int np, int psize, int **dMap=NULL, dataType **dMask=NULL, int **dROI=NULL): slsDetectorData<dataType>(npx, npy, np*psize, dMap, dMask, dROI), nPackets(np), packetSize(psize) {};
/**
Returns the frame number for the given dataset. Virtual func: works for slsDetectorReceiver data (also for each packet), but can be overloaded.
\param buff pointer to the dataset
\returns frame number
*/
virtual int getFrameNumber(char *buff){return ((*(int*)buff)&(0xffffff00))>>8;};
/**
Returns the packet number for the given dataset. Virtual func: works for slsDetectorReceiver packets, but can be overloaded.
\param buff pointer to the dataset
\returns packet number number
*/
virtual int getPacketNumber(char *buff){return (*(int*)buff)&0xff;};
/**
Loops over a memory slot until a complete frame is found (i.e. all packets 0 to nPackets, same frame number). Can be overloaded for different kind of detectors!
\param data pointer to the memory to be analyzed
\param ndata size of frame returned
\param dsize size of the memory slot to be analyzed
\returns pointer to the first packet of the last good frame (might be incomplete if npackets lower than the number of packets), or NULL if no frame is found
*/
virtual char *findNextFrame(char *data, int &ndata, int dsize) {
char *retval=NULL, *p=data;
int dd=0;
int fn, fnum=-1, np=0, pnum=-1;
while (dd<=(dsize-packetSize)) {
pnum=getPacketNumber(p);
fn=getFrameNumber(p);
if (pnum<1 || pnum>nPackets) {
cout << "Bad packet number " << pnum << " frame "<< fn << endl;
retval=NULL;
np=0;
} else if (pnum==1) {
retval=p;
if (np>0)
/*cout << "*Incomplete frame number " << fnum << endl;*/
np=0;
fnum=fn;
} else if (fn!=fnum) {
if (fnum!=-1) {
/* cout << " **Incomplete frame number " << fnum << " pnum " << pnum << " " << getFrameNumber(p) << endl;*/
retval=NULL;
}
np=0;
}
p+=packetSize;
dd+=packetSize;
np++;
// cout << pnum << " " << fn << " " << np << " " << dd << " " << dsize << endl;
if (np==nPackets)
if (pnum==nPackets) {
// cout << "Frame found!" << endl;
break;
} else {
cout << "Too many packets for this frame! "<< fnum << " " << pnum << endl;
retval=NULL;
}
}
if (np<nPackets) {
if (np>0)
cout << "Too few packets for this frame! "<< fnum << " " << pnum << endl;
}
ndata=np*packetSize;
// cout << "return " << ndata << endl;
return retval;
};
/**
Loops over a file stream until a complete frame is found (i.e. all packets 0 to nPackets, same frame number). Can be overloaded for different kind of detectors!
\param filebin input file stream (binary)
\returns pointer to the first packet of the last good frame, NULL if no frame is found or last frame is incomplete
*/
virtual char *readNextFrame(ifstream &filebin) {
char *data=new char[packetSize*nPackets];
char *retval=0;
int np=0, nd;
if (filebin.is_open()) {
while (filebin.read(data+np*packetSize,packetSize)) {
if (np==(nPackets-1)) {
retval=findNextFrame(data,nd,packetSize*nPackets);
np=nd/packetSize;
// cout << np << endl;
if (retval==data && np==nPackets) {
// cout << "-" << endl;
return data;
} else if (np>nPackets) {
cout << "too many packets!!!!!!!!!!" << endl;
delete [] data;
return NULL;
} else if (retval!=NULL) {
// cout << "+" << endl;;
for (int ip=0; ip<np; ip++)
memcpy(data+ip*packetSize,retval+ip*packetSize,packetSize);
}
} else if (np>nPackets) {
cout << "*******too many packets!!!!!!!!!!" << endl;
delete [] data;
return NULL;
} else {
// cout << "." << endl;;
np++;
}
}
}
delete [] data;
return NULL;
};
private:
const int nPackets; /**<number of UDP packets constituting one frame */
const int packetSize; /**< size of a udp packet */
};
#endif