mirror of
https://github.com/slsdetectorgroup/slsDetectorPackage.git
synced 2025-04-22 03:40:04 +02:00
some modifications and demos
git-svn-id: file:///afs/psi.ch/project/sls_det_software/svn/slsDetectorCalibration@31 113b152e-814d-439b-b186-022a431db7b5
This commit is contained in:
l_msdetect
parent
f3e58f95ea
commit
dd26dccde4
31
slsDetectorCalibration/demoCreateTree.C
Normal file
31
slsDetectorCalibration/demoCreateTree.C
Normal file
@ -0,0 +1,31 @@
|
||||
{
|
||||
|
||||
//.L moenchReadData.C
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
TFile *fout;
|
||||
THStack *hs2N;
|
||||
|
||||
fout=new TFile("/scratch/outfile.root","RECREATE");
|
||||
|
||||
hs2N=moenchReadData("/data/moench_xbox_20140113/MoTarget_45kV_0_8mA_120V_cds_g4_f00000%%04d000_0.raw",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();
|
||||
|
||||
|
||||
|
||||
}
|
||||
|
||||
@ -289,7 +289,9 @@ energyCalibration::~energyCalibration(){
|
||||
|
||||
|
||||
|
||||
TH1F* energyCalibration::createMedianHistogram(TH2F* h2, int ch0, int nch) {
|
||||
|
||||
|
||||
TH1F* energyCalibration::createMedianHistogram(TH2F* h2, int ch0, int nch, int direction) {
|
||||
|
||||
if (h2==NULL || nch==0)
|
||||
return NULL;
|
||||
@ -299,9 +301,8 @@ TH1F* energyCalibration::createMedianHistogram(TH2F* h2, int ch0, int nch) {
|
||||
|
||||
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);
|
||||
@ -309,10 +310,20 @@ TH1F* energyCalibration::createMedianHistogram(TH2F* h2, int ch0, int nch) {
|
||||
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;
|
||||
|
||||
|
||||
|
||||
}
|
||||
|
||||
|
||||
@ -480,6 +491,7 @@ TGraphErrors* energyCalibration::linearCalibration(int nscan, Double_t *en, Doub
|
||||
eoff=fitfun->GetParError(0);
|
||||
|
||||
gain=funcs->setScanSign()*mypar[1];
|
||||
|
||||
off=mypar[0];
|
||||
|
||||
return gr;
|
||||
|
||||
@ -251,7 +251,15 @@ class energyCalibration {
|
||||
|
||||
#ifdef MYROOT
|
||||
|
||||
static TH1F* createMedianHistogram(TH2F* h2, int ch0, int nch);
|
||||
/**
|
||||
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
|
||||
|
||||
11
slsDetectorCalibration/gMapDemo.C
Normal file
11
slsDetectorCalibration/gMapDemo.C
Normal file
@ -0,0 +1,11 @@
|
||||
{
|
||||
//.L energyCalibration.cpp+
|
||||
//.L gainMap.C+
|
||||
TFile fin("/data/moench_xbox_20140113/MoTarget_45kV_0_8mA_120V_cds_g4.root");
|
||||
TH2F *h2=fin.Get("h2");
|
||||
TH2F *gMap=gainMap(h2,4);
|
||||
gMap->Draw("colz");
|
||||
|
||||
|
||||
|
||||
}
|
||||
228
slsDetectorCalibration/gainMap.C
Normal file
228
slsDetectorCalibration/gainMap.C
Normal file
@ -0,0 +1,228 @@
|
||||
#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;
|
||||
|
||||
|
||||
}
|
||||
@ -26,7 +26,7 @@ using namespace std;
|
||||
#define NR 160
|
||||
|
||||
|
||||
//#define MY_DEBUG 1
|
||||
#define MY_DEBUG 1
|
||||
#ifdef MY_DEBUG
|
||||
#include <TCanvas.h>
|
||||
#endif
|
||||
|
||||
@ -1,21 +1,24 @@
|
||||
#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.root",tit);
|
||||
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,tit,0,3000,1500,-500,2500,sign,0.,1,159,1,159, 0);
|
||||
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();
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user