mirror of
https://github.com/slsdetectorgroup/slsDetectorPackage.git
synced 2025-05-02 19:00:05 +02:00
quite moenchTree creator for first nanoscope beamtime
git-svn-id: file:///afs/psi.ch/project/sls_det_software/svn/slsDetectorCalibration@47 113b152e-814d-439b-b186-022a431db7b5
This commit is contained in:
l_msdetect
parent
fb40e2236e
commit
0f881677d2
@ -447,10 +447,12 @@ TF1* energyCalibration::fitFunction(TF1 *fun, TH1 *h1, Double_t *mypar, Double_t
|
||||
|
||||
|
||||
fitfun= h1->GetFunction(fname);
|
||||
if (fitfun) {
|
||||
fitfun->GetParameters(mypar);
|
||||
for (int ip=0; ip<6; ip++) {
|
||||
emypar[ip]=fitfun->GetParError(ip);
|
||||
}
|
||||
}
|
||||
return fitfun;
|
||||
}
|
||||
|
||||
|
||||
@ -16,22 +16,25 @@ using namespace std;
|
||||
|
||||
#define FOPT "0"
|
||||
|
||||
TH2F *gainMap(TH2F *h2, float g) {
|
||||
THStack *gainMap(TH2F *h2, float g) {
|
||||
// int npx, npy;
|
||||
int npx=160, npy=160;
|
||||
// det->getDetectorSize(npx, npy);
|
||||
|
||||
THStack *hs=new THStack();
|
||||
|
||||
TH2F *gMap=new TH2F("gmap",h2->GetTitle(), npx, -0.5, npx-0.5, npy, -0.5, npy-0.5);
|
||||
TH2F *nMap=new TH2F("nmap",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];
|
||||
|
||||
hs->Add(gMap);
|
||||
hs->Add(nMap);
|
||||
|
||||
Double_t ens[1]={20.}, eens[1]={20.};
|
||||
Double_t peaks[1], epeaks[1];
|
||||
|
||||
int ibin;
|
||||
TH1D *px;
|
||||
|
||||
|
||||
energyCalibration *enCal=new energyCalibration();
|
||||
enCal->setPlotFlag(0);
|
||||
// enCal->setChargeSharing(0);
|
||||
@ -50,94 +53,48 @@ TH2F *gainMap(TH2F *h2, float g) {
|
||||
TGraph *glin;
|
||||
Double_t peakdum, hpeakdum;
|
||||
|
||||
for (int ix=1; ix<npx-1; ix++) {
|
||||
for (int iy=1; iy<npy-1; iy++) {
|
||||
int ix=20;
|
||||
int iy=40;
|
||||
|
||||
|
||||
for ( ix=1; ix<npx-1; ix++) {
|
||||
for ( 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);
|
||||
enCal->setFitRange(50,3000);
|
||||
//enCal->setChargeSharing(0);
|
||||
if (px) {
|
||||
enCal->fixParameter(0,0);
|
||||
enCal->fixParameter(1,0);
|
||||
enCal->fixParameter(5,0);
|
||||
mypar[0]=0;
|
||||
mypar[1]=0;
|
||||
mypar[2]=peakX[ip];
|
||||
mypar[3]=g*10;
|
||||
mypar[4]=peakY[ip];
|
||||
mypar[2]=px->GetBinCenter(px->GetMaximumBin());
|
||||
mypar[3]=10;
|
||||
mypar[4]=px->GetMaximum();
|
||||
mypar[5]=0;
|
||||
|
||||
|
||||
|
||||
enCal->setStartParameters(mypar);
|
||||
enCal->fitSpectrum(px,mypar,emypar);
|
||||
|
||||
|
||||
peaks[ip+1]=mypar[2];
|
||||
epeaks[ip+1]=emypar[2];
|
||||
cout << ix << " " << iy << " " << mypar[2] << endl;
|
||||
}
|
||||
|
||||
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 (mypar[2]>0) {
|
||||
gMap->SetBinContent(ix+1,iy+1,mypar[2]/ens[0]);
|
||||
gMap->SetBinError(ix+1,iy+1,emypar[2]/ens[0]);
|
||||
nMap->SetBinContent(ix+1,iy+1,mypar[3]);
|
||||
nMap->SetBinError(ix+1,iy+1,emypar[3]);
|
||||
}
|
||||
if (glin)
|
||||
delete glin;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
return gMap;
|
||||
return hs;
|
||||
}
|
||||
|
||||
|
||||
TH2F *noiseMap(TH2F *h2, TH2F *gMap=NULL) {
|
||||
TH2F *noiseMap(TH2F *h2) {
|
||||
// int npx, npy;
|
||||
int npx=160, npy=160;
|
||||
// det->getDetectorSize(npx, npy);
|
||||
@ -152,90 +109,37 @@ TH2F *noiseMap(TH2F *h2, TH2F *gMap=NULL) {
|
||||
cout << ix << " " << iy << " " << ibin << endl;
|
||||
ibin=h2->GetYaxis()->FindBin(ix*npy+iy);
|
||||
px=h2->ProjectionX("px",ibin,ibin);
|
||||
px->Fit("gaus", FOPT);
|
||||
px->Fit("gaus", FOPT,"",-100,100);
|
||||
if (px->GetFunction("gaus")) {
|
||||
if (px->GetFunction("gaus")->GetParameter(1)>-5 && px->GetFunction("gaus")->GetParameter(1)<5)
|
||||
nMap->SetBinContent(ix+1,iy+1,px->GetFunction("gaus")->GetParameter(2));
|
||||
}
|
||||
// delete px;
|
||||
}
|
||||
}
|
||||
if (gMap) {
|
||||
nMap->Divide(gMap);
|
||||
nMap->Scale(1000./3.6);
|
||||
}
|
||||
|
||||
|
||||
return nMap;
|
||||
}
|
||||
|
||||
|
||||
THStack *noiseHistos(char *tit) {
|
||||
char fname[10000];
|
||||
|
||||
sprintf(fname,"/data/moench_xbox_201401_trees/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_201401_trees/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;
|
||||
|
||||
|
||||
}
|
||||
|
||||
THStack *noiseHistos(TH2F *nmap) {
|
||||
THStack *noiseHistos(TH2F *nmap, TH2F *gmap=NULL) {
|
||||
|
||||
|
||||
char tit[1000];
|
||||
|
||||
|
||||
if (nmap==NULL) {
|
||||
cout << "No noise map" << endl;
|
||||
|
||||
return NULL;
|
||||
}
|
||||
|
||||
if (gmap) {
|
||||
nmap->Divide(gmap);
|
||||
nmap->Scale(1000./3.6);
|
||||
}
|
||||
|
||||
|
||||
strcpy(tit,nmap->GetTitle());
|
||||
THStack *hs=new THStack(tit,tit);
|
||||
hs->SetTitle(tit);
|
||||
@ -260,7 +164,7 @@ THStack *noiseHistos(TH2F *nmap) {
|
||||
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 << " or " << h->GetFunction("gaus")->GetParameter(1) << "+-" << h->GetFunction("gaus")->GetParameter(2);
|
||||
}
|
||||
cout << endl;
|
||||
}
|
||||
|
||||
@ -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
|
||||
@ -79,6 +79,7 @@ THStack *moenchReadData(char *fformat, char *tit, int runmin, int runmax, int nb
|
||||
TH2F *h3;
|
||||
TH2F *hetaX;
|
||||
TH2F *hetaY;
|
||||
TH2D *clustHist;
|
||||
|
||||
THStack *hs=new THStack("hs",fformat);
|
||||
|
||||
@ -117,6 +118,9 @@ THStack *moenchReadData(char *fformat, char *tit, int runmin, int runmax, int nb
|
||||
|
||||
|
||||
#ifdef MY_DEBUG
|
||||
quadrant quad;
|
||||
tall->Branch("q",&quad,"q/I");
|
||||
|
||||
|
||||
TCanvas *myC;
|
||||
TH2F *he;
|
||||
@ -124,6 +128,9 @@ THStack *moenchReadData(char *fformat, char *tit, int runmin, int runmax, int nb
|
||||
TCanvas *cH2;
|
||||
TCanvas *cH3;
|
||||
|
||||
int quadrants[5];
|
||||
for(int i = 0; i < 5; i++){ quadrants[i] = 0; }
|
||||
|
||||
if (hitfinder) {
|
||||
myC=new TCanvas();
|
||||
he=new TH2F("he","Event Mask",xmax-xmin, xmin, xmax, ymax-ymin, ymin, ymax);
|
||||
@ -138,6 +145,8 @@ THStack *moenchReadData(char *fformat, char *tit, int runmin, int runmax, int nb
|
||||
cH3=new TCanvas();
|
||||
cH3->SetLogz();
|
||||
h3->Draw("colz");
|
||||
|
||||
clustHist= new TH2D("clustHist","clustHist",3,-1.5,1.5,3,-1.5,1.5);
|
||||
}
|
||||
#endif
|
||||
filter->newDataSet();
|
||||
@ -145,7 +154,7 @@ THStack *moenchReadData(char *fformat, char *tit, int runmin, int runmax, int nb
|
||||
|
||||
for (int irun=runmin; irun<runmax; irun++) {
|
||||
sprintf(fname,fformat,irun);
|
||||
cout << "file name " << fname << endl;
|
||||
cout << "file name " << fname << " ( " << (((double)(irun-runmin))*100./((double)(runmax-runmin))) << "% )" << endl;
|
||||
filebin.open((const char *)(fname), ios::in | ios::binary);
|
||||
nph=0;
|
||||
while ((buff=decoder->readNextFrame(filebin))) {
|
||||
@ -189,8 +198,16 @@ THStack *moenchReadData(char *fformat, char *tit, int runmin, int runmax, int nb
|
||||
h3->Fill(filter->getClusterTotal(3), iy+NR*ix);
|
||||
iFrame=decoder->getFrameNumber(buff);
|
||||
|
||||
tall->Fill();
|
||||
|
||||
#ifdef MY_DEBUG
|
||||
for(int cx=-1; cx <2; cx++)
|
||||
for(int cy=-1; cy < 2; cy++)
|
||||
clustHist->Fill(cx,cy,filter->getClusterElement(cx,cy));
|
||||
|
||||
quad = filter->getQuadrant();
|
||||
quadrants[quad]++;
|
||||
#endif
|
||||
tall->Fill();
|
||||
|
||||
}
|
||||
|
||||
@ -217,10 +234,10 @@ THStack *moenchReadData(char *fformat, char *tit, int runmin, int runmax, int nb
|
||||
#endif
|
||||
nf++;
|
||||
|
||||
cout << "=" ;
|
||||
//cout << "=" ;
|
||||
delete [] buff;
|
||||
}
|
||||
cout << nph << endl;
|
||||
//cout << nph << endl;
|
||||
if (filebin.is_open())
|
||||
filebin.close();
|
||||
else
|
||||
@ -233,6 +250,14 @@ THStack *moenchReadData(char *fformat, char *tit, int runmin, int runmax, int nb
|
||||
|
||||
delete decoder;
|
||||
cout << "Read " << nf << " frames" << endl;
|
||||
|
||||
#ifdef MY_DEBUG
|
||||
cout << "quadrants: " ;
|
||||
for(int i = 0; i<5; i++) cout << i << ": " << quadrants[i] << " || " ;
|
||||
cout << endl;
|
||||
cout << "Read Events " << nph << endl;
|
||||
#endif
|
||||
|
||||
return hs;
|
||||
}
|
||||
|
||||
|
||||
@ -9,23 +9,30 @@ typedef struct task_s{
|
||||
int runmin;
|
||||
int runmax;
|
||||
int treeIndex;
|
||||
int tNumber;
|
||||
double xTalkFactor;
|
||||
} Task;
|
||||
|
||||
double hc = 0.1; // read - out crosstalk
|
||||
|
||||
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;
|
||||
double xTalkC = t->xTalkFactor; //((double)t->tNumber) * 0.01; Calibrated with data in /data/moench_trieste_20140313_trees/xtalkScan
|
||||
cout << "Call moenchReadData(" << t->fformat << "," << t->tname << "," << t->runmin<< "," << t->runmax <<") with xTalkC: " << xTalkC << endl;
|
||||
TThread::UnLock();
|
||||
moenchReadData(t->fformat,t->tname,t->runmin,t->runmax);
|
||||
moenchReadData(t->fformat,t->tname,t->runmin,t->runmax, 1500, -500, 1000, 1, xTalkC);
|
||||
if(f && f->IsOpen()){
|
||||
f->Close();
|
||||
}
|
||||
return 0;
|
||||
}
|
||||
|
||||
|
||||
void moenchReadDataMT(char *fformat, char *tit, char *tdir, int runmin, int runoffset, int nThreads, int treeIndexStart=0){
|
||||
void moenchReadDataMT(char *fformat, char *tit, char *tdir, int runmin, int runoffset, int nThreads, int treeIndexStart=0, double xTalkFactor=0.044){
|
||||
char threadName[1000];
|
||||
TThread *threads[nThreads];
|
||||
for(int i = 0; i < nThreads; i++){
|
||||
@ -37,6 +44,8 @@ void moenchReadDataMT(char *fformat, char *tit, char *tdir, int runmin, int runo
|
||||
t->runmin = runmin + i*runoffset;
|
||||
t->runmax = runmin + (i+1)*runoffset - 1;
|
||||
t->treeIndex = treeIndexStart + i;
|
||||
t->tNumber = i;
|
||||
t->xTalkFactor = xTalkFactor;
|
||||
cout << "start thread " << i << " start: " << t->runmin << " end " << t->runmax << endl;
|
||||
threads[i] = new TThread(threadName, moenchMakeTreeTask, t);
|
||||
threads[i]->Run();
|
||||
@ -46,14 +55,16 @@ void moenchReadDataMT(char *fformat, char *tit, char *tdir, int runmin, int runo
|
||||
threads[i]->Join();
|
||||
}
|
||||
|
||||
cout << " ( DONE ) " << endl;
|
||||
|
||||
}
|
||||
|
||||
|
||||
|
||||
//to compile: g++ -DMYROOT -DMYROOT1 -g `root-config --cflags --glibs` -o moenchReadDataMT moenchReadDataMT.C
|
||||
int main(int argc, char **argv){
|
||||
if(argc != 8){
|
||||
cout << "Usage: " << argv[0] << " fformat tit tdir runmin runoffset nThreads treeIndexStart" << endl;
|
||||
if(argc < 8){
|
||||
cout << "Usage: " << argv[0] << " fformat tit tdir runmin runoffset nThreads treeIndexStart [xTalkFactor]" << endl;
|
||||
exit(-1);
|
||||
}
|
||||
|
||||
@ -64,8 +75,12 @@ int main(int argc, char **argv){
|
||||
int runoffset = atoi(argv[5]);
|
||||
int nThreads = atoi(argv[6]);
|
||||
int treeIndexStart = atoi(argv[7]);
|
||||
double xTalkFactor = 0.044;
|
||||
|
||||
moenchReadDataMT(fformat, tit, tdir,runmin,runoffset,nThreads,treeIndexStart);
|
||||
if(argc == 9)
|
||||
xTalkFactor = atof(argv[8]);
|
||||
|
||||
moenchReadDataMT(fformat, tit, tdir,runmin,runoffset,nThreads,treeIndexStart, xTalkFactor);
|
||||
|
||||
}
|
||||
|
||||
|
||||
@ -14,12 +14,12 @@ void raedNoiseData(char *tit, int sign=1){
|
||||
TFile *fout;
|
||||
THStack *hs2N;
|
||||
|
||||
sprintf(fname,"/data/moench_xbox_20140113/MoTarget_45kV_0_8mA_120V_%s.root",tit);
|
||||
sprintf(fname,"/data/moench_trieste_calibration_trees/flat_20keV_%s.root",tit);
|
||||
fout=new TFile(fname,"RECREATE");
|
||||
|
||||
sprintf(fname,"/data/moench_xbox_20140113/MoTarget_45kV_0_8mA_120V_%s_f00000%%04d000_0.raw",tit);
|
||||
sprintf(fname,"/data/moench_trieste_calibration/flat_20keV_%s_f00000%%04d000_0.raw",tit);
|
||||
|
||||
hs2N=moenchReadData(fname,tit,0,3000,1500,-500,2500,sign,0.,1,159,1,159, 0,1);
|
||||
hs2N=moenchReadData(fname,tit,0,2000,1500,-500,2500,sign,0.,1,159,1,159, 0,1);
|
||||
hs2N->SetName(tit);
|
||||
hs2N->SetTitle(tit);
|
||||
(TH2F*)(hs2N->GetHists()->At(0))->Write();
|
||||
@ -100,10 +100,10 @@ void raedNoiseDataN(char *tit, int sign=1){
|
||||
TFile *fout;
|
||||
THStack *hs2N;
|
||||
|
||||
sprintf(fname,"/data/moench_xbox_20140116/noise_%s.root",tit);
|
||||
sprintf(fname,"/data/moench_noise_20140327_trees/noise_%s.root",tit);
|
||||
fout=new TFile(fname,"RECREATE");
|
||||
|
||||
sprintf(fname,"/data/moench_xbox_20140116/noise_%s_f00000%%04d000_0.raw",tit);
|
||||
sprintf(fname,"/data/moench_noise_20140327/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);
|
||||
|
||||
@ -9,7 +9,7 @@
|
||||
#include "commonModeSubtraction.h"
|
||||
|
||||
|
||||
//#define MYROOT1
|
||||
#define MYROOT1
|
||||
|
||||
#ifdef MYROOT1
|
||||
#include <TTree.h>
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user