Include sls_detector_calibration files
This commit is contained in:
@ -1,11 +1,11 @@
|
||||
// to analyse the backplane pulsing data per module
|
||||
// changes by VH 210906: to eliminate hardcoded absolute paths, uses location of the analysis root files as additional input argument (accordingly changed in filename_creator.sh)
|
||||
|
||||
#include "../sls_detector_calibration/jungfrauCommonHeader.h"
|
||||
#include "../sls_detector_calibration/jungfrauCommonFunctions.h"
|
||||
#include "sls_detector_calibration/jungfrauCommonHeader.h"
|
||||
#include "sls_detector_calibration/jungfrauCommonFunctions.h"
|
||||
|
||||
#include "../sls_detector_calibration/jungfrauFile.C"
|
||||
#include "../sls_detector_calibration/jungfrauPedestal.C"
|
||||
#include "sls_detector_calibration/jungfrauFile.C"
|
||||
#include "sls_detector_calibration/jungfrauPedestal.C"
|
||||
|
||||
#include "TGraphErrors.h"
|
||||
#include "TF1.h"
|
||||
|
||||
@ -2,12 +2,12 @@
|
||||
// changes by VH 210906: to eliminate hardcoded absolute paths, uses location of the analysis root files as additional input argument (accordingly changed in filename_creator.sh)
|
||||
|
||||
|
||||
#include "../sls_detector_calibration/jungfrauCommonHeader.h"
|
||||
#include "../sls_detector_calibration/jungfrauCommonFunctions.h"
|
||||
#include "sls_detector_calibration/jungfrauCommonHeader.h"
|
||||
#include "sls_detector_calibration/jungfrauCommonFunctions.h"
|
||||
|
||||
#include "../sls_detector_calibration/jungfrauFile.C"
|
||||
#include "../sls_detector_calibration/jungfrauPedestal.C"
|
||||
#include "../sls_detector_calibration/jungfrauPixelMask.C"
|
||||
#include "sls_detector_calibration/jungfrauFile.C"
|
||||
#include "sls_detector_calibration/jungfrauPedestal.C"
|
||||
#include "sls_detector_calibration/jungfrauPixelMask.C"
|
||||
|
||||
#include "TGraph.h"
|
||||
#include "TGraphErrors.h"
|
||||
|
||||
@ -4,15 +4,15 @@
|
||||
// and save peak position and uncertainty
|
||||
// changes by VH 210906: to eliminate hardcoded absolute paths, uses location of the analysis root files as additional input argument (accordingly changed in filename_creator.sh)
|
||||
|
||||
#include "../sls_detector_calibration/jungfrauCommonHeader.h"
|
||||
#include "../sls_detector_calibration/jungfrauCommonFunctions.h"
|
||||
#include "sls_detector_calibration/jungfrauCommonHeader.h"
|
||||
#include "sls_detector_calibration/jungfrauCommonFunctions.h"
|
||||
|
||||
#include "../sls_detector_calibration/jungfrauFile.C"
|
||||
#include "../sls_detector_calibration/jungfrauPixelMask.C"
|
||||
#include "../sls_detector_calibration/jungfrauPedestal.C"
|
||||
#include "sls_detector_calibration/jungfrauFile.C"
|
||||
#include "sls_detector_calibration/jungfrauPixelMask.C"
|
||||
#include "sls_detector_calibration/jungfrauPedestal.C"
|
||||
|
||||
#include "../sls_detector_calibration/energyCalibration.h"
|
||||
#include "../sls_detector_calibration/energyCalibration.cpp"
|
||||
#include "sls_detector_calibration/energyCalibration.h"
|
||||
#include "sls_detector_calibration/energyCalibration.cpp"
|
||||
|
||||
#include "TGraph.h"
|
||||
#include "TGraphErrors.h"
|
||||
|
||||
145
sls_detector_calibration/MovingStat.h
Executable file
145
sls_detector_calibration/MovingStat.h
Executable file
@ -0,0 +1,145 @@
|
||||
#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), u(0) {}
|
||||
|
||||
/**
|
||||
clears the moving average number of samples parameter, mean and standard deviation
|
||||
*/
|
||||
void Clear()
|
||||
{
|
||||
m_n = 0;
|
||||
m_newM=0;
|
||||
m_newM2=0;
|
||||
u = 0;
|
||||
}
|
||||
|
||||
void ResetUpdates()
|
||||
{
|
||||
u = 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) {
|
||||
u++;
|
||||
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 );
|
||||
}
|
||||
|
||||
/** returns the number of times Calc was called
|
||||
\returns the number of times Calc was called
|
||||
*/
|
||||
inline double GetNUpdates() const
|
||||
{
|
||||
return u;
|
||||
}
|
||||
|
||||
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 */
|
||||
int u; /**< number of times Calc was called (number of updates) */
|
||||
};
|
||||
#endif
|
||||
798
sls_detector_calibration/energyCalibration.cpp
Normal file
798
sls_detector_calibration/energyCalibration.cpp
Normal file
@ -0,0 +1,798 @@
|
||||
#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;
|
||||
// Gaussian exponent
|
||||
if (par[3]!=0) {
|
||||
arg=sign*(x[0]-par[2])/par[3];
|
||||
}
|
||||
// the Gaussian
|
||||
f=TMath::Exp(-1*arg*arg/2.);
|
||||
// Gaussian + error function
|
||||
f=f+par[5]/2.*(TMath::Erfc(arg/(TMath::Sqrt(2.))));
|
||||
// Gaussian + error function + pedestal
|
||||
return par[4]*f+pedestal(x,par);
|
||||
}
|
||||
|
||||
Double_t energyCalibrationFunctions::gaussChargeSharingKb(Double_t *x, Double_t *par) {
|
||||
Double_t f, arg=0,argb=0;
|
||||
// Gaussian exponent
|
||||
if (par[3]!=0) {
|
||||
arg=sign*(x[0]-par[2])/par[3];
|
||||
argb=sign*(x[0]-(par[6]*par[2]))/par[3]; // using absolute kb mean might seem better but like this the ratio can be fixed
|
||||
}
|
||||
// the Gaussian
|
||||
f=TMath::Exp(-1*arg*arg/2.);
|
||||
f=f+par[7]*(TMath::Exp(-1*argb*argb/2.));
|
||||
// Gaussian + error function
|
||||
f=f+par[5]/2.*(TMath::Erfc(arg/(TMath::Sqrt(2.))));
|
||||
f=f+par[7]*par[5]/2.*(TMath::Erfc(argb/(TMath::Sqrt(2.))));
|
||||
// Gaussian + error function + pedestal
|
||||
return par[4]*f+pedestal(x,par);
|
||||
}
|
||||
|
||||
Double_t energyCalibrationFunctions::gaussChargeSharingKaDoublet(Double_t *x, Double_t *par) {
|
||||
Double_t f, f2, arg=0, arg2=0;
|
||||
// Gaussian exponent
|
||||
if (par[3]!=0) {
|
||||
arg=sign*(x[0]-par[2])/par[3];
|
||||
arg2=sign*(x[0]-par[6])/par[3];
|
||||
}
|
||||
// the Gaussian
|
||||
f=TMath::Exp(-1*arg*arg/2.);
|
||||
f2=TMath::Exp(-1*arg2*arg2/2.);
|
||||
// Gaussian + error function
|
||||
f=f+par[5]/2.*(TMath::Erfc(arg/(TMath::Sqrt(2.))));
|
||||
f2=f2+par[5]/2.*(TMath::Erfc(arg/(TMath::Sqrt(2.)))); // shouldn't this be arg2?
|
||||
// Gaussian + error function + pedestal
|
||||
return par[4]*f+par[7]*f2+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::spectrumkb(Double_t *x, Double_t *par) {
|
||||
return gaussChargeSharingKb(x,par);
|
||||
}
|
||||
|
||||
Double_t energyCalibrationFunctions::spectrumkadoublet(Double_t *x, Double_t *par) {
|
||||
return gaussChargeSharingKaDoublet(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),
|
||||
kb_mean(-1),
|
||||
kb_frac(-1),
|
||||
mean2(-1),
|
||||
ampl2(-1),
|
||||
fscurve(NULL),
|
||||
fspectrum(NULL),
|
||||
fspectrumkb(NULL),
|
||||
fspectrumkadoublet(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");
|
||||
fspectrumkb=new TF1("fspectrumkb",funcs,&energyCalibrationFunctions::spectrumkb,0,1000,8,"energyCalibrationFunctions","spectrumkb");
|
||||
fspectrumkb->SetParNames("Background Pedestal","Background slope", "Peak position","Noise RMS", "Number of Photons","Charge Sharing Pedestal","kb mean","kb frac");
|
||||
|
||||
fspectrumkadoublet=new TF1("fspectrumkadoublet",funcs,&energyCalibrationFunctions::spectrumkadoublet,0,1000,8,"energyCalibrationFunctions","spectrumkadoublet");
|
||||
fspectrumkadoublet->SetParNames("Background Pedestal","Background slope", "Peak position","Noise RMS", "Number of Photons","Charge Sharing Pedestal","ka2 mean","n2");
|
||||
|
||||
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);
|
||||
fspectrumkb->FixParameter(ip, val);
|
||||
fspectrumkadoublet->FixParameter(ip, val);
|
||||
}
|
||||
|
||||
|
||||
void energyCalibration::releaseParameter(int ip){
|
||||
|
||||
fscurve->ReleaseParameter(ip);
|
||||
fspectrum->ReleaseParameter(ip);
|
||||
fspectrumkb->ReleaseParameter(ip);
|
||||
fspectrumkadoublet->ReleaseParameter(ip);
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
energyCalibration::~energyCalibration(){
|
||||
#ifdef MYROOT
|
||||
delete fscurve;
|
||||
delete fspectrum;
|
||||
delete fspectrumkb;
|
||||
delete fspectrumkadoublet;
|
||||
#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::setStartParametersKb(Double_t *par){
|
||||
bg_offset=par[0];
|
||||
bg_slope=par[1];
|
||||
flex=par[2];
|
||||
noise=par[3];
|
||||
ampl=par[4];
|
||||
cs_slope=par[5];
|
||||
kb_mean=par[6];
|
||||
kb_frac=par[7];
|
||||
//fit_min = 400; // used for soleil flat field
|
||||
//fit_max = 800;
|
||||
}
|
||||
|
||||
void energyCalibration::setStartParametersKaDoublet(Double_t *par){
|
||||
bg_offset=par[0];
|
||||
bg_slope=par[1];
|
||||
flex=par[2];
|
||||
noise=par[3];
|
||||
ampl=par[4];
|
||||
cs_slope=par[5];
|
||||
mean2=par[6];
|
||||
ampl2=par[7];
|
||||
//fit_min = 400; // used for soleil flat field
|
||||
//fit_max = 800;
|
||||
}
|
||||
|
||||
|
||||
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);
|
||||
fspectrumkb->ReleaseParameter(1);
|
||||
fspectrumkadoublet->ReleaseParameter(1);
|
||||
} else {
|
||||
fscurve->FixParameter(5,0);
|
||||
fspectrum->FixParameter(1,0);
|
||||
fspectrumkb->FixParameter(1,0);
|
||||
fspectrumkadoublet->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);
|
||||
|
||||
}
|
||||
|
||||
void energyCalibration::initFitFunctionKb(TF1 *fun, TH1 *h1) {
|
||||
|
||||
Double_t min=fit_min, max=fit_max;
|
||||
|
||||
Double_t mypar[8];
|
||||
|
||||
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;
|
||||
|
||||
if (kb_mean==-1)
|
||||
mypar[6]=0;
|
||||
else
|
||||
mypar[6]=kb_mean;
|
||||
|
||||
if (kb_frac==-1)
|
||||
mypar[7]=0;
|
||||
else
|
||||
mypar[7]=kb_frac;
|
||||
|
||||
fun->SetParameters(mypar);
|
||||
|
||||
fun->SetRange(min,max);
|
||||
|
||||
}
|
||||
|
||||
void energyCalibration::initFitFunctionKaDoublet(TF1 *fun, TH1 *h1) {
|
||||
|
||||
Double_t min=fit_min, max=fit_max;
|
||||
|
||||
Double_t mypar[8];
|
||||
|
||||
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;
|
||||
|
||||
if (mean2==-1)
|
||||
mypar[6]=0;
|
||||
else
|
||||
mypar[6]=mean2;
|
||||
|
||||
if (ampl2==-1)
|
||||
mypar[7]=0;
|
||||
else
|
||||
mypar[7]=ampl2;
|
||||
|
||||
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::fitFunctionKb(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<8; ip++) {
|
||||
emypar[ip]=fitfun->GetParError(ip);
|
||||
}
|
||||
return fitfun;
|
||||
}
|
||||
|
||||
TF1* energyCalibration::fitFunctionKaDoublet(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<8; 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);
|
||||
}
|
||||
|
||||
TF1* energyCalibration::fitSpectrumKb(TH1 *h1, Double_t *mypar, Double_t *emypar) {
|
||||
initFitFunctionKb(fspectrumkb,h1);
|
||||
return fitFunctionKb(fspectrumkb, h1, mypar, emypar);
|
||||
}
|
||||
|
||||
TF1* energyCalibration::fitSpectrumKaDoublet(TH1 *h1, Double_t *mypar, Double_t *emypar) {
|
||||
initFitFunctionKaDoublet(fspectrumkadoublet,h1);
|
||||
return fitFunctionKaDoublet(fspectrumkadoublet, 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
|
||||
|
||||
|
||||
|
||||
|
||||
469
sls_detector_calibration/energyCalibration.h
Normal file
469
sls_detector_calibration/energyCalibration.h
Normal file
@ -0,0 +1,469 @@
|
||||
#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);
|
||||
Double_t gaussChargeSharingKb(Double_t *x, Double_t *par);
|
||||
Double_t gaussChargeSharingKaDoublet(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);
|
||||
Double_t spectrumkb(Double_t *x, Double_t *par);
|
||||
Double_t spectrumkadoublet(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);
|
||||
void setStartParametersKb(Double_t *par);
|
||||
void setStartParametersKaDoublet(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);
|
||||
TF1 *fitSpectrumKb(TH1 *h1, Double_t *mypar, Double_t *emypar);
|
||||
TF1 *fitSpectrumKaDoublet(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);
|
||||
void initFitFunctionKb(TF1 *fun, TH1 *h1);
|
||||
void initFitFunctionKaDoublet(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);
|
||||
TF1 *fitFunctionKb(TF1 *fun, TH1 *h1, Double_t *mypar, Double_t *emypar);
|
||||
TF1 *fitFunctionKaDoublet(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 */
|
||||
Double_t kb_mean;
|
||||
Double_t kb_frac;
|
||||
Double_t mean2;
|
||||
Double_t ampl2;
|
||||
|
||||
TF1 *fscurve; /**< function with which the s-curve will be fitted */
|
||||
|
||||
TF1 *fspectrum; /**< function with which the spectrum will be fitted */
|
||||
TF1 *fspectrumkb; /**< function with which the spectrum will be fitted */
|
||||
TF1 *fspectrumkadoublet; /**< 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
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
344
sls_detector_calibration/jungfrauCommonFunctions.h
Normal file
344
sls_detector_calibration/jungfrauCommonFunctions.h
Normal file
@ -0,0 +1,344 @@
|
||||
#ifndef JUNGFRAUCOMMONFUNCTIONS_H
|
||||
#define JUNGFRAUCOMMONFUNCTIONS_H
|
||||
|
||||
#include <algorithm>
|
||||
|
||||
int chipOfPixel(int i) {
|
||||
int chip = 0;
|
||||
if (i/1024 >= 256 && i%1024 >= 0 && i%1024 < 256) {
|
||||
chip = 1;
|
||||
} else if (i/1024 >= 256 && i%1024 >= 256 && i%1024 < 512) {
|
||||
chip = 2;
|
||||
} else if (i/1024 >= 256 && i%1024 >= 512 && i%1024 < 768) {
|
||||
chip = 3;
|
||||
} else if (i/1024 >= 256 && i%1024 >= 768 && i%1024 < 1024) {
|
||||
chip = 4;
|
||||
} else if (i/1024 < 256 && i%1024 >= 0 && i%1024 < 256) {
|
||||
chip = 5;
|
||||
} else if (i/1024 < 256 && i%1024 >= 256 && i%1024 < 512) {
|
||||
chip = 6;
|
||||
} else if (i/1024 < 256 && i%1024 >= 512 && i%1024 < 768) {
|
||||
chip = 7;
|
||||
} else if (i/1024 < 256 && i%1024 >= 768 && i%1024 < 1024) {
|
||||
chip = 8;
|
||||
} else {
|
||||
cout << "problem" << endl;
|
||||
}
|
||||
return chip;
|
||||
}
|
||||
|
||||
|
||||
int supercolumnOfPixel(int i) {
|
||||
int sc = ((i-((i/256)*256))/64) + 1;
|
||||
return sc;
|
||||
}
|
||||
|
||||
|
||||
int ADCOfPixel(int i) {
|
||||
int adc = -1;
|
||||
if (i%1024 < 256) {
|
||||
adc = 0;
|
||||
} else if (i%1024 < 512) {
|
||||
adc = 1;
|
||||
} else if (i%1024 < 768) {
|
||||
adc = 2;
|
||||
} else {
|
||||
adc = 3;
|
||||
}
|
||||
return adc;
|
||||
}
|
||||
|
||||
|
||||
void LoadPaletteFalse() {
|
||||
|
||||
// jet
|
||||
Double_t stops[5] = { 0.00, 0.34, 0.61, 0.84, 1.00 };
|
||||
Double_t red[5] = { 0.00, 0.00, 0.87, 1.00, 0.51 };
|
||||
Double_t green[5] = { 0.00, 0.81, 1.00, 0.20, 0.00 };
|
||||
Double_t blue[5] = { 0.51, 1.00, 0.12, 0.00, 0.00 };
|
||||
TColor::CreateGradientColorTable(5, stops, red, green, blue, 255, 1.0);
|
||||
gStyle->SetNumberContours(255);
|
||||
|
||||
}
|
||||
|
||||
|
||||
void UsePaletteViridis() {
|
||||
|
||||
// viridis
|
||||
Double_t stops[9] = {0.0000, 0.1250, 0.2500, 0.3750, 0.5000, 0.6250, 0.7500, 0.8750, 1.0000};
|
||||
Double_t red[9] = {26./255., 51./255., 43./255., 33./255., 28./255., 35./255., 74./255., 144./255., 246./255.};
|
||||
Double_t green[9] = {9./255., 24./255., 55./255., 87./255., 118./255., 150./255., 180./255., 200./255., 222./255.};
|
||||
Double_t blue[9] = {30./255., 96./255., 112./255., 114./255., 112./255., 101./255., 72./255., 35./255., 0./255.};
|
||||
TColor::CreateGradientColorTable(9, stops, red, green, blue, 255, 1.0);
|
||||
gStyle->SetNumberContours(255);
|
||||
|
||||
}
|
||||
|
||||
|
||||
void UsePaletteDBR() {
|
||||
|
||||
// dark body radiator
|
||||
Double_t stops[9] = {0.0000, 0.1250, 0.2500, 0.3750, 0.5000, 0.6250, 0.7500, 0.8750, 1.0000};
|
||||
Double_t red[9] = {0./255., 45./255., 99./255., 156./255., 212./255., 230./255., 237./255., 234./255., 242./255.};
|
||||
Double_t green[9] = {0./255., 0./255., 0./255., 45./255., 101./255., 168./255., 238./255., 238./255., 243./255.};
|
||||
Double_t blue[9] = {0./255., 1./255., 1./255., 3./255., 9./255., 8./255., 11./255., 95./255., 230./255.};
|
||||
TColor::CreateGradientColorTable(9, stops, red, green, blue, 255, 1.0);
|
||||
gStyle->SetNumberContours(255);
|
||||
|
||||
}
|
||||
|
||||
|
||||
void jungfrauStyle()
|
||||
{
|
||||
gROOT->SetStyle("Plain"); /*Default white background for all plots*/
|
||||
/*set bkg color of all to white*/
|
||||
gStyle->SetCanvasColor(kWhite);
|
||||
gStyle->SetFrameFillColor(kWhite);
|
||||
gStyle->SetStatColor(kWhite);
|
||||
gStyle->SetPadColor(kWhite);
|
||||
gStyle->SetFillColor(10);
|
||||
gStyle->SetTitleFillColor(kWhite);
|
||||
|
||||
/* SetPaperSize wants width & height in cm: A4 is 20,26 & US is 20,24*/
|
||||
gStyle->SetPaperSize(20, 26);
|
||||
/* No yellow border around histogram*/
|
||||
gStyle->SetDrawBorder(0);
|
||||
/* remove border of canvas*/
|
||||
gStyle->SetCanvasBorderMode(0);
|
||||
/* remove border of pads*/
|
||||
gStyle->SetPadBorderMode(0);
|
||||
gStyle->SetFrameBorderMode(0);
|
||||
gStyle->SetLegendBorderSize(0);
|
||||
|
||||
/* default text size*/
|
||||
gStyle->SetTextSize(0.05);
|
||||
gStyle->SetTitleSize(0.07,"xyz");
|
||||
gStyle->SetLabelSize(0.06,"xyz");
|
||||
/* title offset: distance between given text and axis, here x,y,z*/
|
||||
gStyle->SetLabelOffset(0.015,"xyz");
|
||||
gStyle->SetTitleOffset(1.1,"yz");
|
||||
gStyle->SetTitleOffset(1.0,"x");
|
||||
|
||||
/* Use visible font for all text*/
|
||||
int font = 42;
|
||||
gStyle->SetTitleFont(font);
|
||||
gStyle->SetTitleFontSize(0.06);
|
||||
gStyle->SetStatFont(font);
|
||||
gStyle->SetStatFontSize(0.07);
|
||||
gStyle->SetTextFont(font);
|
||||
gStyle->SetLabelFont(font,"xyz");
|
||||
gStyle->SetTitleFont(font,"xyz");
|
||||
gStyle->SetTitleBorderSize(0);
|
||||
gStyle->SetStatBorderSize(1);
|
||||
gStyle->SetLegendFont(font);
|
||||
|
||||
/* big marker points*/
|
||||
gStyle->SetMarkerStyle(1);
|
||||
gStyle->SetLineWidth(2);
|
||||
gStyle->SetMarkerSize(1.2);
|
||||
/*set palette in 2d histogram to nice and colorful one*/
|
||||
gStyle->SetPalette(1,0);
|
||||
|
||||
/*No title for histograms*/
|
||||
gStyle->SetOptTitle(0);
|
||||
/* show the errors on the stat box */
|
||||
gStyle->SetOptStat(0);
|
||||
/* show errors on fitted parameters*/
|
||||
gStyle->SetOptFit(0);
|
||||
/* number of decimals used for errors*/
|
||||
gStyle->SetEndErrorSize(5);
|
||||
|
||||
/* set line width to 2 by default so that histograms are visible when printed small
|
||||
idea: emphasize the data, not the frame around*/
|
||||
gStyle->SetHistLineWidth(2);
|
||||
gStyle->SetFrameLineWidth(2);
|
||||
gStyle->SetFuncWidth(2);
|
||||
gStyle->SetHistLineColor(kBlack);
|
||||
gStyle->SetFuncColor(kRed);
|
||||
gStyle->SetLabelColor(kBlack,"xyz");
|
||||
|
||||
//set the margins
|
||||
gStyle->SetPadBottomMargin(0.2);
|
||||
gStyle->SetPadTopMargin(0.05);
|
||||
gStyle->SetPadRightMargin(0.05);
|
||||
gStyle->SetPadLeftMargin(0.2);
|
||||
|
||||
//set the number of divisions to show
|
||||
gStyle->SetNdivisions(506, "xy");
|
||||
|
||||
//turn off xy grids
|
||||
gStyle->SetPadGridX(0);
|
||||
gStyle->SetPadGridY(0);
|
||||
|
||||
//set the tick mark style
|
||||
gStyle->SetPadTickX(1);
|
||||
gStyle->SetPadTickY(1);
|
||||
|
||||
gStyle->SetCanvasDefW(800);
|
||||
gStyle->SetCanvasDefH(700);
|
||||
|
||||
LoadPaletteFalse();
|
||||
|
||||
gROOT->ForceStyle();
|
||||
}
|
||||
|
||||
// types of pixels
|
||||
bool isQuad(int i) {
|
||||
if ( (i/NC == 255 || i/NC == 256) && (i%NC == 255 || i%NC == 256 || i%NC == 511 || i%NC == 512 || i%NC == 767 || i%NC == 768) ) {
|
||||
return true;
|
||||
} else {
|
||||
return false;
|
||||
}
|
||||
}
|
||||
|
||||
bool isDouble(int i) {
|
||||
if(i/NC == 255 || i/NC == 256 || i%NC == 255 || i%NC == 256 || i%NC == 511 || i%NC == 512 || i%NC == 767 || i%NC == 768) {
|
||||
if (not isQuad(i)) {
|
||||
return true;
|
||||
} else {
|
||||
return false;
|
||||
}
|
||||
} else {
|
||||
return false;
|
||||
}
|
||||
}
|
||||
|
||||
bool isNextToDouble(int i) {
|
||||
if(i/NC == 254 || i/NC == 257 || i%NC == 254 || i%NC == 257 || i%NC == 510 || i%NC == 513 || i%NC == 766 || i%NC == 769) {
|
||||
if (not isDouble(i) && not isQuad(i)) {
|
||||
return true;
|
||||
} else {
|
||||
return false;
|
||||
}
|
||||
} else {
|
||||
return false;
|
||||
}
|
||||
}
|
||||
|
||||
bool isEdge(int i) {
|
||||
if (i%NC == 0 || i%NC == 1023 || i/NC == 0 || i/NC == 511) {
|
||||
if (not isDouble(i) && not isNextToDouble(i)) {
|
||||
return true;
|
||||
} else {
|
||||
return false;
|
||||
}
|
||||
} else {
|
||||
return false;
|
||||
}
|
||||
}
|
||||
|
||||
bool isInnerEdge(int i) {
|
||||
if (i%NC == 1 || i%NC == 1022 || i/NC == 1 || i/NC == 510) {
|
||||
if (not isEdge(i) && not isDouble(i) && not isNextToDouble(i)) {
|
||||
return true;
|
||||
} else {
|
||||
return false;
|
||||
}
|
||||
} else {
|
||||
return false;
|
||||
}
|
||||
}
|
||||
|
||||
bool isBulk(int i) {
|
||||
if (not isEdge(i) && not isInnerEdge(i) && not isDouble(i) && not isQuad(i) && not isNextToDouble(i)) {
|
||||
return true;
|
||||
} else {
|
||||
return false;
|
||||
}
|
||||
}
|
||||
|
||||
double highestPointBeforeSwitching(const vector<double> &lower_filter, const vector<double> &higher_filter) {
|
||||
|
||||
double highest_point_before_switching;
|
||||
|
||||
if (higher_filter.size() > 0) {
|
||||
// find the highest value in lower_filter that is lower than all entries in higher_filter
|
||||
highest_point_before_switching = *min_element(lower_filter.begin(),lower_filter.end())-1;
|
||||
|
||||
double lowest_entry_in_higher_filter = *min_element(higher_filter.begin(),higher_filter.end());
|
||||
for(vector<double>::const_iterator it = lower_filter.begin(); it != lower_filter.end(); ++it) {
|
||||
if (*it < lowest_entry_in_higher_filter) {
|
||||
if (*it > highest_point_before_switching) {
|
||||
highest_point_before_switching = *it;
|
||||
}
|
||||
}
|
||||
}
|
||||
} else {
|
||||
highest_point_before_switching = *max_element(lower_filter.begin(),lower_filter.end());
|
||||
}
|
||||
|
||||
return highest_point_before_switching;
|
||||
|
||||
}
|
||||
|
||||
double lowestPointAfterSwitching(const vector<double> &higher_filter, const vector<double> &lower_filter) {
|
||||
|
||||
double lowest_point_after_switching;
|
||||
|
||||
if (lower_filter.size() > 0) {
|
||||
// find the lowest value in higher_filter that is higher than all entries in lower_filter
|
||||
lowest_point_after_switching = *max_element(higher_filter.begin(),higher_filter.end())+1;
|
||||
|
||||
double highest_entry_in_lower_filter = *max_element(lower_filter.begin(),lower_filter.end());
|
||||
for(vector<double>::const_iterator it = higher_filter.begin(); it != higher_filter.end(); ++it) {
|
||||
if (*it > highest_entry_in_lower_filter) {
|
||||
if (*it < lowest_point_after_switching) {
|
||||
lowest_point_after_switching = *it;
|
||||
}
|
||||
}
|
||||
}
|
||||
} else {
|
||||
lowest_point_after_switching = *min_element(higher_filter.begin(),higher_filter.end());
|
||||
}
|
||||
|
||||
return lowest_point_after_switching;
|
||||
|
||||
}
|
||||
|
||||
double highestPointBeforeSaturation(const vector<double> &filter, const vector<double> &adc_g2) {
|
||||
|
||||
// first look if saturation ever occurs
|
||||
bool saturation = false;
|
||||
for(vector<double>::const_iterator it_adc = adc_g2.begin(); it_adc != adc_g2.end(); ++it_adc) {
|
||||
if (*it_adc == 0) {
|
||||
saturation = true;
|
||||
}
|
||||
}
|
||||
|
||||
if (saturation == false) {
|
||||
|
||||
double highest_point_before_saturation = *max_element(filter.begin(),filter.end());
|
||||
return highest_point_before_saturation;
|
||||
|
||||
} else {
|
||||
|
||||
// find the lowest filter value which saturates
|
||||
double lowest_point_which_saturates = *max_element(filter.begin(),filter.end());
|
||||
|
||||
vector<double>::const_iterator it_adc = adc_g2.begin();
|
||||
|
||||
for(vector<double>::const_iterator it_filter = filter.begin(); it_filter != filter.end(); ++it_filter) {
|
||||
if (*it_adc == 0) {
|
||||
if (*it_filter < lowest_point_which_saturates) {
|
||||
lowest_point_which_saturates = *it_filter;
|
||||
}
|
||||
}
|
||||
++it_adc;
|
||||
}
|
||||
|
||||
// find the highest filter value which is below saturation
|
||||
double highest_point_before_saturation = *min_element(filter.begin(),filter.end());
|
||||
for(vector<double>::const_iterator it_filter = filter.begin(); it_filter != filter.end(); ++it_filter) {
|
||||
if (*it_filter < lowest_point_which_saturates) {
|
||||
if (*it_filter > highest_point_before_saturation) {
|
||||
highest_point_before_saturation = *it_filter;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return highest_point_before_saturation;
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
#endif /* JUNGFRAUCOMMONFUNCTIONS_H */
|
||||
14
sls_detector_calibration/jungfrauCommonHeader.h
Normal file
14
sls_detector_calibration/jungfrauCommonHeader.h
Normal file
@ -0,0 +1,14 @@
|
||||
#include<iostream> // std::cout
|
||||
|
||||
// LoadPalette
|
||||
#include "TColor.h"
|
||||
#include "TStyle.h"
|
||||
// jungfrauStyle
|
||||
#include "TROOT.h"
|
||||
|
||||
#define NC 1024
|
||||
#define NR 512
|
||||
#define NCH NC*NR
|
||||
#define NSC 16
|
||||
|
||||
using namespace std;
|
||||
150
sls_detector_calibration/jungfrauFile.C
Normal file
150
sls_detector_calibration/jungfrauFile.C
Normal file
@ -0,0 +1,150 @@
|
||||
//#include "jungfrauCommonHeader.h"
|
||||
|
||||
#include <fstream> // std::ifstream
|
||||
#include <stdint.h> // uint16_t
|
||||
#include <ctime> // clock
|
||||
#include <iomanip> // std::setprecision
|
||||
|
||||
#include <time.h>
|
||||
#include <sys/time.h>
|
||||
|
||||
class oneframe {
|
||||
|
||||
public:
|
||||
uint64_t framenumber;
|
||||
uint64_t bunchid;
|
||||
uint16_t imagedata[NCH];
|
||||
|
||||
};
|
||||
|
||||
class jungfrauFile {
|
||||
|
||||
private:
|
||||
|
||||
public:
|
||||
ifstream filebin; // file descriptor
|
||||
oneframe thisframe;
|
||||
long long int filebinSize; // need 64 bit
|
||||
clock_t filebinStart;
|
||||
clock_t filebinStop;
|
||||
struct timeval tss,tsss; //for timing
|
||||
|
||||
jungfrauFile() { // constructor
|
||||
|
||||
cout << "jungfrauFile constructed" << endl;
|
||||
//cout << "sizeof(thisframe) " << sizeof(thisframe) << endl;
|
||||
//cout << "address of data array " << &(thisframe.imagedata)<< endl;
|
||||
|
||||
}
|
||||
|
||||
uint16_t * getFrameDataHandle() {
|
||||
|
||||
return (uint16_t*)&(thisframe.imagedata);
|
||||
|
||||
}
|
||||
|
||||
bool readNextFrame() {
|
||||
|
||||
if (filebin.is_open()) {
|
||||
if (!filebin.eof()) {
|
||||
filebin.read(((char*)&thisframe),sizeof(thisframe));
|
||||
|
||||
std::cout << std::fixed << std::setprecision(2)
|
||||
<< "\r [" << std::string(int(filebin.tellg()*10 / float(filebinSize)), '#')
|
||||
<< std::string(10 + 1 - int(filebin.tellg()*10 / float(filebinSize)), ' ')
|
||||
<< "] " << 100 * filebin.tellg() / float(filebinSize) << "%";
|
||||
|
||||
if (filebin.good()) {return true;} // necessary to stop the frame 'after' the last frame returning true
|
||||
else {
|
||||
filebinStop = clock();
|
||||
double elapsed_secs = double(filebinStop - filebinStart) / CLOCKS_PER_SEC;
|
||||
cout << " end of file reached in " << elapsed_secs << "s" << endl;
|
||||
return false;
|
||||
}
|
||||
|
||||
} else {
|
||||
filebinStop = clock();
|
||||
double elapsed_secs = double(filebinStop - filebinStart) / CLOCKS_PER_SEC;
|
||||
cout << " end of file reached in " << elapsed_secs << "s" << endl;
|
||||
|
||||
return false;
|
||||
}
|
||||
} else {
|
||||
cout << "file not open" << endl;
|
||||
return false;
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
int currentFrameNumber() {
|
||||
|
||||
return int(thisframe.framenumber);
|
||||
|
||||
}
|
||||
|
||||
uint64_t currentBunchID() {
|
||||
|
||||
return uint64_t(thisframe.bunchid);
|
||||
|
||||
}
|
||||
|
||||
/*****************/
|
||||
/* VH 2022-05-05 */
|
||||
/*****************/
|
||||
uint64_t currentSCnumber() {
|
||||
|
||||
return ( ( thisframe.bunchid >> 8 )&0xf );
|
||||
|
||||
}
|
||||
/*****************/
|
||||
/* VH 2022-05-05 */
|
||||
/*****************/
|
||||
|
||||
void rewind() {
|
||||
if (filebin.is_open()) {
|
||||
cout << "rewinding file" << endl;
|
||||
filebin.clear();
|
||||
filebin.seekg(0);
|
||||
} else {
|
||||
cout << "file was not open" << endl;
|
||||
}
|
||||
}
|
||||
|
||||
void open(char *fformat, int fileindex) {
|
||||
|
||||
char fname[1000];
|
||||
sprintf(fname,fformat,fileindex);
|
||||
cout << "file name " << fname << endl;
|
||||
if (filebin.is_open()) {
|
||||
cout << "WARNING: a file is already open" << endl;
|
||||
cout << "WARNING: closing this file" << endl;
|
||||
filebin.close();
|
||||
}
|
||||
cout << "opening file" << endl;
|
||||
filebin.open((const char *)(fname), ios::in | ios::binary);
|
||||
filebin.seekg(0,filebin.end);
|
||||
filebinSize = filebin.tellg();
|
||||
filebin.seekg(0,filebin.beg);
|
||||
filebinStart = clock();
|
||||
}
|
||||
|
||||
void close() {
|
||||
|
||||
if (filebin.is_open()) {
|
||||
filebin.close();
|
||||
cout << "closed file" << endl;
|
||||
} else {
|
||||
cout << "file was not open" << endl;
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
bool isOpen() {
|
||||
if (filebin.is_open()) {
|
||||
return true;
|
||||
} else {
|
||||
return false;
|
||||
}
|
||||
}
|
||||
|
||||
};
|
||||
217
sls_detector_calibration/jungfrauPedestal.C
Normal file
217
sls_detector_calibration/jungfrauPedestal.C
Normal file
@ -0,0 +1,217 @@
|
||||
#include "MovingStat.h"
|
||||
|
||||
#include "TCanvas.h"
|
||||
#include "TH1F.h"
|
||||
#include "TH2F.h"
|
||||
|
||||
class jungfrauPedestal {
|
||||
|
||||
private:
|
||||
MovingStat stat [NCH];
|
||||
|
||||
public:
|
||||
TH1F* pedestals1d;
|
||||
TH2F* pedestals2d;
|
||||
TH1F* pedestalvars1d;
|
||||
TH2F* pedestalvars2d;
|
||||
TH2F* pedestalsvars;
|
||||
|
||||
jungfrauPedestal() {
|
||||
cout << "jungfrauPedestal constructed" << endl;
|
||||
}
|
||||
|
||||
bool addFrameToPedestalCalculation(uint16_t *imagedata) {
|
||||
for (int i = 0; i<NCH; i++) {
|
||||
stat[i].Calc(double(imagedata[i]&0x3fff));
|
||||
}
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
bool addG0FrameToPedestalCalculation(uint16_t *imagedata) {
|
||||
for (int i = 0; i<NCH; i++) {
|
||||
if (((imagedata[i]&0xc000) >> 14) == 0) {
|
||||
stat[i].Calc(double(imagedata[i]&0x3fff));
|
||||
}
|
||||
}
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
bool addG0PixelToPedestalCalculation(int i, double adc_double) {
|
||||
stat[i].Calc(adc_double);
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
bool addG1FrameToPedestalCalculation(uint16_t *imagedata) {
|
||||
for (int i = 0; i<NCH; i++) {
|
||||
//if ((((imagedata[i]&0xc000) >> 14) == 1)) {
|
||||
if ((((imagedata[i]&0xc000) >> 14) == 1)||(((imagedata[i]&0xc000) >> 14) == 2)) { //to fix FW error during JF11 FW develop. needed for modules 243-379
|
||||
stat[i].Calc(double(imagedata[i]&0x3fff));
|
||||
}
|
||||
}
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
bool addG2FrameToPedestalCalculation(uint16_t *imagedata) {
|
||||
for (int i = 0; i<NCH; i++) {
|
||||
if (((imagedata[i]&0xc000) >> 14) == 3) {
|
||||
stat[i].Calc(double(imagedata[i]&0x3fff));
|
||||
}
|
||||
}
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
// single threshold
|
||||
bool addG0FrameToPedestalCalculationWThreshold(uint16_t *imagedata, jungfrauPedestal *pedestalObj, uint16_t threshold) {
|
||||
for (int i = 0; i<NCH; i++) {
|
||||
if (((imagedata[i]&0xc000) >> 14) == 0) {
|
||||
if (fabs((imagedata[i]&0x3fff) - pedestalObj->pedestalOfChannel(i)) < threshold) {
|
||||
stat[i].Calc(double(imagedata[i]&0x3fff));
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
// array of thresholds
|
||||
bool addG0FrameToPedestalCalculationWThreshold(uint16_t *imagedata, jungfrauPedestal *pedestalObj, double *thresholds) {
|
||||
for (int i = 0; i<NCH; i++) {
|
||||
if (((imagedata[i]&0xc000) >> 14) == 0) {
|
||||
if (fabs((imagedata[i]&0x3fff) - pedestalObj->pedestalOfChannel(i)) < thresholds[i]) {
|
||||
stat[i].Calc(double(imagedata[i]&0x3fff));
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
double pedestalOfChannel(int ichannel) {
|
||||
return stat[ichannel].Mean();
|
||||
}
|
||||
|
||||
double semOfChannel(int ichannel) {
|
||||
return stat[ichannel].StandardDeviation() / sqrt(stat[ichannel].NumDataValues());
|
||||
}
|
||||
|
||||
double rmsOfChannel(int ichannel) {
|
||||
return stat[ichannel].StandardDeviation();
|
||||
}
|
||||
|
||||
double varianceOfChannel(int ichannel) {
|
||||
return stat[ichannel].Variance();
|
||||
}
|
||||
|
||||
bool pedestalData(double* pedestals) {
|
||||
for (int i = 0; i < NCH; i++) {
|
||||
pedestals[i] = stat[i].Mean();
|
||||
}
|
||||
return true;
|
||||
}
|
||||
|
||||
bool pedestalData(uint16_t* pedestals) {
|
||||
for (int i = 0; i < NCH; i++) {
|
||||
pedestals[i] = uint16_t(round(stat[i].Mean()));
|
||||
}
|
||||
return true;
|
||||
}
|
||||
|
||||
bool pedestalRMSData(double* pedestalRMSs) {
|
||||
for (int i = 0; i < NCH; i++) {
|
||||
pedestalRMSs[i] = stat[i].StandardDeviation();
|
||||
}
|
||||
return true;
|
||||
}
|
||||
|
||||
int pedestalNFrames() {
|
||||
return stat[0].GetN();
|
||||
}
|
||||
|
||||
int pedestalUpdates(int ichannel) {
|
||||
return stat[ichannel].GetNUpdates();
|
||||
}
|
||||
|
||||
void pedestalResetUpdates() {
|
||||
for (int i=0; i < NCH; i++){
|
||||
stat[i].ResetUpdates();
|
||||
}
|
||||
}
|
||||
|
||||
void pedestalClear() {
|
||||
for (int i=0; i < NCH; i++){
|
||||
stat[i].Clear();
|
||||
}
|
||||
}
|
||||
|
||||
bool pedestalSetNFrames(int nframes) {
|
||||
for (int i=0; i < NCH; i++){
|
||||
stat[i].SetN(nframes);
|
||||
}
|
||||
return true;
|
||||
}
|
||||
|
||||
void savePedestals(uint16_t* pedestals, string fileName) {
|
||||
ofstream pedefile (fileName.c_str());
|
||||
if (pedefile.is_open()) {
|
||||
cout << "saving pedestals" << endl;
|
||||
for (int i=0; i<NCH; i++){
|
||||
pedefile << pedestals[i] << " " ;
|
||||
}
|
||||
pedefile.close();
|
||||
} else {
|
||||
cout << "Unable to open file" << endl;
|
||||
}
|
||||
}
|
||||
|
||||
void savePedestals(double* pedestals, string fileName) {
|
||||
ofstream pedefile (fileName.c_str());
|
||||
if (pedefile.is_open()) {
|
||||
cout << "saving pedestals" << endl;
|
||||
for (int i=0; i<NCH; i++){
|
||||
pedefile << pedestals[i] << " " ;
|
||||
}
|
||||
pedefile.close();
|
||||
} else {
|
||||
cout << "Unable to open file" << endl;
|
||||
}
|
||||
}
|
||||
|
||||
bool readPedestals(uint16_t* pedestals, string fileName) {
|
||||
ifstream pedefile(fileName.c_str());
|
||||
if (pedefile.is_open()) {
|
||||
cout << "reading pedestal file" << endl;
|
||||
for (int i = 0; i < NCH; i++) {
|
||||
pedefile >> pedestals[i];
|
||||
}
|
||||
return true;
|
||||
} else {
|
||||
cout << "Unable to open file" << endl;
|
||||
return false;
|
||||
}
|
||||
}
|
||||
|
||||
bool readPedestals(double* pedestals, string fileName) {
|
||||
ifstream pedefile(fileName.c_str());
|
||||
if (pedefile.is_open()) {
|
||||
cout << "reading pedestal file" << endl;
|
||||
for (int i = 0; i < NCH; i++) {
|
||||
pedefile >> pedestals[i];
|
||||
}
|
||||
return true;
|
||||
} else {
|
||||
cout << "Unable to open file" << endl;
|
||||
return false;
|
||||
}
|
||||
}
|
||||
|
||||
bool setPedestalOfChannel(int ichannel, uint16_t set_pede) {
|
||||
stat[ichannel].Calc(double(set_pede));
|
||||
return true;
|
||||
}
|
||||
|
||||
};
|
||||
245
sls_detector_calibration/jungfrauPixelMask.C
Executable file
245
sls_detector_calibration/jungfrauPixelMask.C
Executable file
@ -0,0 +1,245 @@
|
||||
|
||||
#include <algorithm> // std::count
|
||||
#include "TCanvas.h"
|
||||
#include "TH2F.h"
|
||||
|
||||
class jungfrauPixelMask {
|
||||
|
||||
private:
|
||||
|
||||
public:
|
||||
TH2F* pixel_mask2d;
|
||||
TCanvas* maskcanvas;
|
||||
|
||||
jungfrauPixelMask() {
|
||||
cout << "jungfrauPixelMask constructed" << endl;
|
||||
}
|
||||
|
||||
void initialisePixelMask(bool* mask) {
|
||||
std::fill_n(mask, NCH, true); // initialise with true (not masked)
|
||||
cout << "jungfrauPixelMask initialised" << endl;
|
||||
}
|
||||
|
||||
int getNMasked(bool* mask) {
|
||||
return count(mask, mask + NCH, false);
|
||||
}
|
||||
|
||||
void maskPixel(int pixel, bool *mask) {
|
||||
mask[pixel]= false;
|
||||
}
|
||||
|
||||
void maskChip(int chip, bool* mask) {
|
||||
if (chip < 1 || chip > 8){
|
||||
cout << "expected chip number 1-8" << endl;
|
||||
}
|
||||
for (int i=0; i<NCH; i++) {
|
||||
if (chipOfPixel(i) == chip) {
|
||||
mask[i] = false;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void maskSupercolumn(int chip, int sc, bool* mask) {
|
||||
if (chip < 1 || chip > 8 || sc < 1 || sc > 4){
|
||||
cout << "expected chip number 1-8 and supercolumn number 1-4" << endl;
|
||||
}
|
||||
for (int i=0; i<NCH; i++) {
|
||||
if (chipOfPixel(i) == chip){
|
||||
if (supercolumnOfPixel(i) == sc) {
|
||||
mask[i] = false;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void maskColumn(int chip, int c, bool* mask) {
|
||||
if (chip < 1 || chip > 8 || c < 0 || c >= NC){
|
||||
cout << "expected chip number 1-8 and column number 0--(NC-1)" << endl;
|
||||
}
|
||||
for (int i=0; i<NCH; i++) {
|
||||
if (chipOfPixel(i) == chip){
|
||||
if (i%NC == c){
|
||||
mask[i] = false;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void maskChipEdges(bool* mask) {
|
||||
for (int i=0; i<NCH; i++) {
|
||||
if (i%1024 == 0 ||
|
||||
i%1024 == 255 ||
|
||||
i%1024 == 256 ||
|
||||
i%1024 == 511 ||
|
||||
i%1024 == 512 ||
|
||||
i%1024 == 767 ||
|
||||
i%1024 == 768 ||
|
||||
i%1024 == 1023 ||
|
||||
i/1024 == 0 ||
|
||||
i/1024 == 255 ||
|
||||
i/1024 == 256 ||
|
||||
i/1024 == 511) {
|
||||
mask[i] = false;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void maskNeighbours(bool* mask) {
|
||||
bool tmp_mask [NCH];
|
||||
std::fill_n(tmp_mask, NCH, true);
|
||||
|
||||
for (int i=0; i<NCH; i++) {
|
||||
if (mask[i] == false) {
|
||||
tmp_mask[i] = false;
|
||||
|
||||
if (i<1024) { // bottom
|
||||
if (i==0) { // bottom left corner
|
||||
tmp_mask[i+1] = false;
|
||||
tmp_mask[i+1024] = false;
|
||||
tmp_mask[i+1024+1] = false;
|
||||
} else if (i==1023) { // bottom right corner
|
||||
tmp_mask[i-1] = false;
|
||||
tmp_mask[i+1024] = false;
|
||||
tmp_mask[i+1024-1] = false;
|
||||
} else { // bottom centre
|
||||
tmp_mask[i-1] = false;
|
||||
tmp_mask[i+1] = false;
|
||||
tmp_mask[i+1024-1] = false;
|
||||
tmp_mask[i+1024] = false;
|
||||
tmp_mask[i+1024+1] = false;
|
||||
}
|
||||
} else if (i>=523264){ // top
|
||||
if (i==523264){ // top left corner
|
||||
tmp_mask[i+1] = false;
|
||||
tmp_mask[i-1024] = false;
|
||||
tmp_mask[i-1024+1] = false;
|
||||
} else if (i==524287){ // top right corner
|
||||
tmp_mask[i-1] = false;
|
||||
tmp_mask[i-1024] = false;
|
||||
tmp_mask[i-1024-1] = false;
|
||||
} else { // top center
|
||||
tmp_mask[i+1] = false;
|
||||
tmp_mask[i-1] = false;
|
||||
tmp_mask[i-1024+1] = false;
|
||||
tmp_mask[i-1024] = false;
|
||||
tmp_mask[i-1024-1] = false;
|
||||
}
|
||||
} else if ((i%1024)==0) { // left hand side center
|
||||
tmp_mask[i+1] = false;
|
||||
tmp_mask[i+1024] = false;
|
||||
tmp_mask[i+1024+1] = false;
|
||||
tmp_mask[i-1024] = false;
|
||||
tmp_mask[i-1024+1] = false;
|
||||
} else if ((i%1024)==1023) { // right hand side center
|
||||
tmp_mask[i-1] = false;
|
||||
tmp_mask[i+1024] = false;
|
||||
tmp_mask[i+1024-1] = false;
|
||||
tmp_mask[i-1024] = false;
|
||||
tmp_mask[i-1024-1] = false;
|
||||
} else { // center
|
||||
tmp_mask[i-1] = false;
|
||||
tmp_mask[i+1] = false;
|
||||
tmp_mask[i+1024-1] = false;
|
||||
tmp_mask[i+1024] = false;
|
||||
tmp_mask[i+1024+1] = false;
|
||||
tmp_mask[i-1024-1] = false;
|
||||
tmp_mask[i-1024] = false;
|
||||
tmp_mask[i-1024+1] = false;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
for (int i = 0; i < NCH; i++) {
|
||||
mask[i] = tmp_mask[i];
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
void maskIfGainNot(int expected, uint16_t *imagedata, bool* mask) {
|
||||
// if gain of a pixel is not as expected, mask pixel
|
||||
for (int i = 0; i<NCH; i++) {
|
||||
if (((imagedata[i]&0xc000) >> 14) != expected) {
|
||||
mask[i] = false;
|
||||
//cout << "gain is " << ((imagedata[i]&0xc000) >> 14) << " not " << expected << " as expected so masking" << endl;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void maskFromPedestal(bool* mask, double* pedestals, double* pedestalRMSs) {
|
||||
// if pedestal RMS is outside limits, mask pixel
|
||||
for (int i = 0; i<NCH; i++) {
|
||||
if (pedestalRMSs[i] > 100 || pedestals[i] < 1000 || pedestals[i] > 4000) {
|
||||
mask[i] = false;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void maskIfPedeRMSG0GreaterThan(bool* mask, double* pedestalG0RMSs, int val) {
|
||||
for (int i = 0; i < NCH; i++) {
|
||||
if (pedestalG0RMSs[i] > val) {
|
||||
mask[i] = false;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void maskIfPedestalStep(uint16_t *imagedata, uint16_t* last_frame, bool* mask, int step) {
|
||||
for (int i = 0; i < NCH; i++) {
|
||||
if ( fabs(last_frame[i] - (imagedata[i]&0x3fff)) > step) {
|
||||
mask[i] = false;
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
void plotPixelMask(bool* mask, string fileName) {
|
||||
pixel_mask2d = new TH2F("pixel_mask2d","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);
|
||||
for (int i = 0; i < NCH; i++) {
|
||||
if (mask[i] == false) {
|
||||
pixel_mask2d->Fill(i%NC,i/NC,1);
|
||||
}
|
||||
}
|
||||
|
||||
maskcanvas = new TCanvas("maskcanvas","",150,10,800,400);
|
||||
maskcanvas->SetLeftMargin(0.1);
|
||||
maskcanvas->SetRightMargin(0.13);
|
||||
maskcanvas->SetTopMargin(0.08);
|
||||
maskcanvas->SetBottomMargin(0.15);
|
||||
|
||||
pixel_mask2d->GetXaxis()->SetTitle("Column");
|
||||
pixel_mask2d->GetYaxis()->SetTitle("Row");
|
||||
pixel_mask2d->GetYaxis()->SetTitleOffset(0.7);
|
||||
pixel_mask2d->Draw("colz");
|
||||
maskcanvas->SaveAs(fileName.c_str());
|
||||
maskcanvas->Close();
|
||||
pixel_mask2d->Delete();
|
||||
}
|
||||
|
||||
void saveMask(bool* mask, string fileName) {
|
||||
ofstream maskfile (fileName.c_str());
|
||||
if (maskfile.is_open()) {
|
||||
cout << "saving pixel mask" << endl;
|
||||
for (int i=0; i<NCH; i++){
|
||||
maskfile << mask[i] << " " ;
|
||||
}
|
||||
maskfile.close();
|
||||
} else {
|
||||
cout << "Unable to open file" << endl;
|
||||
}
|
||||
}
|
||||
|
||||
bool readMask(bool* mask, string fileName) {
|
||||
ifstream maskfile(fileName.c_str());
|
||||
if (maskfile.is_open()) {
|
||||
cout << "reading pixel mask" << endl;
|
||||
for (int i = 0; i < NCH; ++i) {
|
||||
maskfile >> mask[i];
|
||||
}
|
||||
return true;
|
||||
} else {
|
||||
cout << "Unable to open pixel mask" << endl;
|
||||
return false;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
};
|
||||
Reference in New Issue
Block a user