Subversion Repositories f9daq

Compare Revisions

Ignore whitespace Rev 34 → Rev 35

/sipmScan/d2r.ini
1,72 → 1,72
# ch pad_center_x pad_center_y tdc_offset # sipm2, tdc modified 20/3/14 pos added 20/2/14
0 75600 5040 142.5
1 55440 5040 142.1
2 75600 15120 142.1
3 55440 15120 141.6
4 45360 25200 142.2
5 65520 25200 141.9
6 45360 35280 143
7 65520 35280 141.3
# ch pad_center_x pad_center_y tdc_offset # sipm2, tdc modified 27/3/14 pos added 20/2/14 threshold scan
0 75600 5040 143.5
1 55440 5040 143.5
2 75600 15120 143.5
3 55440 15120 143.5
4 45360 25200 143.5
5 65520 25200 143.5
6 45360 35280 143.5
7 65520 35280 143.5
 
8 75600 45360 141
9 55440 45360 140
10 75600 55440 140
11 55440 55440 140
12 45360 65520 140
13 65520 65520 140
14 45360 75600 140
15 65520 75600 140
8 75600 45360 143.5
9 55440 45360 142.5
10 75600 55440 142.5
11 55440 55440 142.5
12 45360 65520 142.5
13 65520 65520 142.5
14 45360 75600 142.5
15 65520 75600 142.5
 
16 65520 5040 142.7
17 45360 5040 142
18 65520 15120 142
19 45360 15120 142.1
20 55440 25200 142.1
21 75600 25200 141
22 55440 35280 142
23 75600 35280 142
16 65520 5040 144.5
17 45360 5040 144.5
18 65520 15120 144.5
19 45360 15120 144.5
20 55440 25200 144.5
21 75600 25200 144.5
22 55440 35280 144.5
23 75600 35280 143.5
 
24 65520 45360 142
25 45360 45360 141.5
26 65520 55440 141.8
27 45360 55440 141.9
28 55440 65520 141.3
29 75600 65520 141
30 55440 75600 141.7
31 75600 75600 145
24 65520 45360 143.5
25 45360 45360 143.5
26 65520 55440 143.5
27 45360 55440 143.5
28 55440 65520 143.5
29 75600 65520 143.5
30 55440 75600 143.5
31 75600 75600 146.5
 
32 5040 5040 141
33 25200 5040 141
34 5040 15120 141
35 25200 15120 141
36 35280 25200 141
37 15120 25200 141
38 35280 35280 141
39 15120 35280 140.1
32 5040 5040 143.5
33 25200 5040 143.5
34 5040 15120 143.5
35 25200 15120 143.5
36 35280 25200 143.5
37 15120 25200 143.5
38 35280 35280 144.5
39 15120 35280 144.5
 
40 5040 45360 142.8
41 25200 45360 142.7
42 5040 55440 143.7
43 25200 55440 142.7
44 35280 65520 142.7
45 15120 65520 142.4
46 35280 75600 142.1
47 15120 75600 141.9
40 5040 45360 144.5
41 25200 45360 144.5
42 5040 55440 144.5
43 25200 55440 144.5
44 35280 65520 144.5
45 15120 65520 144.5
46 35280 75600 144.5
47 15120 75600 144.5
 
48 15120 5040 143
49 35280 5040 140.9
50 15120 15120 142
51 35280 15120 141.5
52 25200 25200 140
53 5040 25200 141.3
54 25200 35280 142
55 5040 35280 141
48 15120 5040 144.5
49 35280 5040 144.5
50 15120 15120 144.5
51 35280 15120 144.5
52 25200 25200 143.5
53 5040 25200 142.5
54 25200 35280 142.5
55 5040 35280 142.5
 
56 15120 45360 140
57 35280 45360 140
58 15120 55440 140
59 35280 55440 140
60 25200 65520 140
61 5040 65520 140
62 25200 75600 140
63 5040 75600 140
56 15120 45360 142.5
57 35280 45360 142.5
58 15120 55440 142.5
59 35280 55440 142.5
60 25200 65520 142.5
61 5040 65520 142.5
62 25200 75600 142.5
63 5040 75600 142.5
/sipmScan/src/analysisThreshold.cpp
1,523 → 1,472
//////////////////////////////////////////
// Data to root conversion root script
//
// Contributors: Rok Pestotnik, Rok Dolenec, Dino Tahirovic
//
// 3/1/2014 TDC cut relative to tdc offset
//
 
#include "stdio.h"
#include "TROOT.h"
#include "TFile.h"
#include "TNtuple.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TF1.h"
#include "TMath.h"
#include "TStyle.h"
#include "TCanvas.h"
#include "TLine.h"
#include "zlib.h"
 
// ------------------------------------------------------------------------------
 
#define POSMARG 1000
 
#define READBUFFERLENGTH 10000
 
// data format
#define MAXDATA 16
#define NCH 64
#define TDC_BIN 1.0416 // 1 TDC bin in ns
#define MIKRO_BIN 0.49609/1000. //1 mikro step in mm; stage MM3MF
#define OFFSETX 5100 // Right edge of SiPM+Lightguide
#define OFFSETY 5200 // Lower edge of SiPM+Lightguide
 
#define RUNREC_ID 1
#define ENDREC_ID 2
#define POSREC_ID 3
#define EVTREC_ID 4
#define THRREC_ID 5
 
typedef struct {
unsigned int id,len;
unsigned int fver,time;
unsigned int nev,nch,ped,xy;
unsigned int thLow, thUp, thStep;
int nx,x0,dx,ny,y0,dy;
} RUNREC;
RUNREC *runrec;
RUNREC run;
 
typedef struct {
unsigned int id,len;
unsigned int time;
} ENDREC;
ENDREC *endrec;
 
typedef struct {
unsigned int id,len;
unsigned int time;
int ix;
int x;
int xset;
int iy;
int y;
int yset;
} POSREC;
POSREC *posrec;
POSREC pos;
 
typedef struct {
unsigned int id;
unsigned int len;
unsigned int nev;
} EVTREC;
EVTREC *evtrec;
 
typedef struct {
unsigned int id;
unsigned int len;
unsigned int threshold;
} THRREC;
THRREC *thrrec;
THRREC thr;
 
double padCenter[NCH][2];
 
int position(int, int, int);
 
// ------------------------------------------------------------------------------
 
int d2r(char* dfile0="test", int dbg=0, double tdcCut=5.0)
{
 
const double c_tdcOffset = +2.5; // ns
printf(" Data to root conversion program\nUsage:\nd2r(input file name <without .dat>, debug on/off, TDC cut +-[ns])\n\n");
char fullname[256];
char sbuff[256];
FILE *fp;
//Chanel information
double tdcOffset[NCH];
sprintf(fullname, "d2r.ini");
if( (fp=fopen(fullname, "rt")) == NULL )
printf("Cannot open pad centers file %s !!!\n", fullname);
else {
printf("Opened pad centers file %s\n", fullname);
char* result = fgets(sbuff,256, fp);
if (dbg) printf("Read buffer %s\n", result);
printf("%s", sbuff);
for(int i=0; i<NCH; i++) {
int channel;
int message = fscanf(fp, "%d %lf %lf %lf\n", &channel, &padCenter[i][0], &padCenter[i][1], &tdcOffset[i]);
if (dbg) printf("Read d2r.ini returned %d\n", message);
}
fclose(fp);
}
for(int i=0; i<NCH; i++) {
tdcOffset[i] += c_tdcOffset;
printf("%.2lf %.2lf %.2lf\n", padCenter[i][0], padCenter[i][1], tdcOffset[i]);
}
//TDC correction parameters
/*
double corpar[NCH][3]={ {-0.8, 25, 0},
{-0.8, 25, 0},
{-0.8, 25, 0},
{-0.8, 25, 0},
{-0.8, 25, 0},
{-0.8, 25, 0},
{-0.8, 25, 0},
{-0.8, 25, 0} };
sprintf(fullname, "data/%s_cor.txt", dfile0);
if( (fp=fopen(fullname, "rt")) == NULL )
printf("Cannot open parameter file %s !!!\n", fullname);
else {
printf("Opened parameter file %s\n", fullname);
for(int i=0; i<NCH; i++) {
fscanf(fp, "%lf %lf %lf\n", &corpar[i][0], &corpar[i][1], &corpar[i][2]);*/
/*// check if parameters make sense
if( (corpar[i][0] < (tdcmi-0.2*tdcmi)) || ((tdcma+0.2*tdcma) < corpar[i][0]) ||
(corpar[i][1] < 0) || (1e4*TDC_BIN < corpar[i][1]) ||
(1e4 < TMath::Abs(corpar[i][2])) ) {
printf("Warning: parameters for ch%d out of limits -> using default!\n", i);
corpar[i][0]=2200*TDC_BIN; corpar[i][1]=1000*TDC_BIN; corpar[i][2]=-100*TDC_BIN;
}*/
/* }
fclose(fp);
}*/
//for(int i=0; i<NCH; i++) printf("%.2lf %.2lf %.2lf\n", corpar[i][0], corpar[i][1], corpar[i][2]);
//histograms
char hname[256];
//double tdc;
TH2F *htdc;
TH2F* h_correctedTDC;
TH1F *hnhitsx[NCH], *hnhitsy[NCH];
TH2F *h2d[NCH];
TH2F *h_threshold;
TH2F *h_ch33;
TNtuple *nt;
//data buffer
unsigned int readbuf[READBUFFERLENGTH];
unsigned int buf[READBUFFERLENGTH];
//data file
gzFile dfp;
char dfile[256];
int ftype=0;
int fcount=1;
do {
switch(ftype++) {
case 0:
sprintf(dfile, "./data/%s_file%02d.dat", dfile0, fcount);
break;
case 1:
sprintf(dfile, "./data/%s_file%02d.dat.gz", dfile0, fcount);
break;
case 2:
sprintf(dfile, "./data/%s_file%02d.gz", dfile0, fcount);
break;
default:
printf(" Cannot find data file for %s !!!\n", dfile0);
return -1;
}
dfp=gzopen(dfile,"rb");
} while(!dfp);
printf("Opened data file %s\n", dfile);
//opens ROOT file
//TFile *rootfile;
char fnameroot[256];
sprintf(fnameroot, "root/%s.root", dfile0);
//rootfile = (TFile *) gROOT->FindObjectAny(dfile0);
//if (rootfile!=NULL) {printf("!!!\n");rootfile->Close();}
//rootfile = new TFile(fnameroot);
//if(rootfile) rootfile->Close();
TFile* rootfile = new TFile(fnameroot,"RECREATE",dfile0);
// -----------------------------------------------
// loop trough records
unsigned int rec_id, rec_len;
unsigned int ulsize = sizeof(unsigned int);
//unsigned int ulsize = 4;
if (dbg) printf("Size of unsigned int: %lu\n", sizeof(unsigned int));
int ceve=0;
int end_of_file = 0;
int ii;
int nint;
int nb;
int nSteps;
int status;
while(1) {
if(gzeof(dfp)) end_of_file = 1;
 
gzread(dfp, (voidp)&readbuf, 2*ulsize);
rec_id=readbuf[0];
rec_len=readbuf[1];
if(dbg) printf("-----------------------------------------------\n");
if(dbg) printf("[%d] rec_id = %d | rec_len = %u\n", ceve, rec_id, rec_len);
 
switch(rec_id)
{
case RUNREC_ID:
gzread(dfp, (voidp)&readbuf[2], (rec_len-2*ulsize));
runrec = (RUNREC*) readbuf;
run = *runrec;
if(dbg) {
printf("RUNREC_ID\n");
printf("id = %d, len = %d, time = %d\n", run.id, run.len, run.time);
printf("nev = %d, nch = %d\n", run.nev, run.nch);
printf("nx = %d, x0 = %d, dx = %d\n", run.nx, run.x0, run.dx);
printf("ny = %d, y0 = %d, dy = %d\n", run.ny, run.y0, run.dy);
printf("thLow = %d, thUp = %d, thStep = %d\n", run.thLow, run.thUp, run.thStep);
}
//create histograms
nt = new TNtuple("nt", "nt", "ch:x:y:tdc");
sprintf(hname, "htdc");
htdc = (TH2F*) gROOT->FindObject(hname);
if (htdc) delete htdc;
htdc = new TH2F("htdc","Raw TDC;TDC channel;SiPM channel",512,0,512,NCH,0,NCH);
h_correctedTDC = (TH2F*) gROOT->FindObject("h_correctedTDC");
if (h_correctedTDC) delete h_correctedTDC;
h_correctedTDC = new TH2F("h_correctedTDC","Corrected TDC;t [ns];SiPM channel",33, -16.5*TDC_BIN,16.5*TDC_BIN,NCH,0,NCH);
//TH1F* gsumV673A[NCH/16] = new TH1F(hn,hname,256,-0.5,255.5);
for(int i=0; i<NCH; i++) {
/*
sprintf(hname, "htdcpos%d", i);
htdcpos[i] = (TH1F*)gROOT->FindObject(hname);
if(htdcpos[i]) delete htdcpos[i];
htdcpos[i] = new TH1F(hname, hname, 512, 0, 512*TDC_BIN);
sprintf(hname, "htdc%d", i);
htdc[i] = (TH1F*)gROOT->FindObject(hname);
if(htdc[i]) delete htdc[i];
htdc[i] = new TH1F(hname, hname, 512, 0*TDC_BIN, 512*TDC_BIN);
*/
sprintf(hname, "hnhitsx%d", i);
hnhitsx[i] = (TH1F*)gROOT->FindObject(hname);
if(hnhitsx[i]) delete hnhitsx[i];
hnhitsx[i] = new TH1F(hname, hname, run.nx,
(run.x0 - OFFSETX - 0.5*run.dx)*MIKRO_BIN,
(run.x0 - OFFSETX + (run.nx-0.5)*run.dx)*MIKRO_BIN);
sprintf(hname, "hnhitsy%d", i);
hnhitsy[i] = (TH1F*)gROOT->FindObject(hname);
if(hnhitsy[i]) delete hnhitsy[i];
//hnhitsy[i] = new TH1F(hname, hname, run.ny, (run.y0-0.5*run.dy)*MIKRO_BIN, (run.y0+(run.ny-0.5)*run.dy)*MIKRO_BIN);
hnhitsy[i] = new TH1F(hname, hname, run.ny,
(run.y0 - 0.5*run.dy - OFFSETY)*MIKRO_BIN,
(run.y0 + (run.ny-0.5)*run.dy - OFFSETY)*MIKRO_BIN);
//hnhitsy[i] = new TH1F(hname, hname, 100, 0,100);
sprintf(hname, "h2d%d", i);
h2d[i] = (TH2F*)gROOT->FindObject(hname);
if(h2d[i]) delete h2d[i];
h2d[i] = new TH2F(hname, hname, run.nx,
(run.x0 - OFFSETX - 0.5*run.dx)*MIKRO_BIN,
(run.x0 - OFFSETX + (run.nx-0.5)*run.dx)*MIKRO_BIN,
run.ny,
(run.y0 - OFFSETY - 0.5*run.dy)*MIKRO_BIN,
(run.y0 - OFFSETY + (run.ny-0.5)*run.dy)*MIKRO_BIN);
}
if (h_threshold) delete h_threshold;
nSteps = (run.thUp - run.thLow)/double(run.thStep) + 1;
if (dbg) printf("nSteps %d\n", nSteps);
h_threshold = new TH2F("h_threshold","Threshold scan;SiPM ch;Threshold[mV]",64,-0.5,63.5,
nSteps,
run.thLow - 0.5*run.thStep,
run.thUp + 0.5*run.thStep);
//h_threshold = new TH2F("h_threshold","Threshold scan;SiPM ch;Threshold[mV]",64,-0.5,63.5,
// 101,995,2005);
if (h_ch33) delete h_ch33;
h_ch33 = new TH2F("h_ch33","ch. 33;x;y",100,20000,30000,100,0,10000);
break;
case POSREC_ID:
gzread(dfp, (voidp)&readbuf[2], (rec_len-2*ulsize));
posrec = (POSREC *) readbuf;
pos=*posrec;
if(dbg) {
printf("POSREC_ID\n");
printf("id = %d, len = %d, time = %d\n", posrec->id, posrec->len, posrec->time);
printf("ix = %d, x = %d, xset = %d\n", posrec->ix, posrec->x, posrec->xset);
printf("iy = %d, y = %d, yset = %d\n", posrec->iy, posrec->y, posrec->yset);
} else printf(" [%d,%d] %d, %d\n", pos.ix, pos.iy, pos.xset, pos.yset);
break;
case EVTREC_ID:
gzread(dfp, (voidp)&readbuf[2], ulsize); // last field of event record
evtrec = (EVTREC *) readbuf;
//evtrec->nev = buf[0];
//if (rec_len < 0 || rec_len > 10000) {
if (rec_len > READBUFFERLENGTH) {
printf("Len %u\n", rec_len);
return(0);
}
nb = rec_len - 3*ulsize; // no. of bytes to read
gzread(dfp, (voidp)&buf, nb);
if(dbg) {
printf("EVTREC_ID\n");
printf("id = %d, len = %d, nev = %d\n", evtrec->id, evtrec->len, evtrec->nev);
//for(int datai = 0; datai < NCH; datai++) printf("%u ", evtrec->data[datai]);
//printf("\n");
//for(int datai = NCH; datai < NCH+NCH; datai++) printf("%u ", evtrec->data[datai]);
//printf("\n");
break;
}
nint = nb / ulsize; // no. of subrecords
if (dbg) printf("No. of subrecords %d \n", nint);
ii=0;
while (ii<nint){
int recid = buf[ii++];
int len = buf[ii++];
if (dbg) printf("Buffer pointer %d\n", ii);
unsigned int *dbuf = (unsigned int *)&buf[ii];
//if (n%1000==0)
if (dbg) printf("%d 0x%03x Len=%d\n",evtrec->nev,recid,len);
//unsigned short edge;
//int nhits;
if (recid==0x140 || recid==0x141) {
for (int i=0; i<len; i++) {
int data = dbuf[i] & 0xFFFF ;
int edge_type = (dbuf[i]>>16)&0x1 ;
int overflow = (dbuf[i]>>17)&0x1 ;
int tdc_num = (dbuf[i]>>25)&0x1 ;
int channel = ((dbuf[i]>>18)&0x1F) | tdc_num<<5 ;
int ev_dat = (dbuf[i]>>23)&0x1 ;
int last_dat = (dbuf[i]>>30)&0x1 ;
int nval_dat = (dbuf[i]>>31)&0x1 ;
if (dbg){
if (ev_dat) printf("Event %d\n",data);
else printf("ch=%d edge=%d ev=%d data=%d last=%d nval=%d\n",
channel, edge_type,ev_dat,data,last_dat,nval_dat);
}
if (!ev_dat){
if (!edge_type && !overflow) {
htdc->Fill(data, channel);
if(dbg) printf("ch: %d tdc: %d\n", channel, data);
if (dbg) nt->Fill(channel, pos.ix, pos.iy, data);
double tdcmin=tdcOffset[channel] - tdcCut;
double tdcmax=tdcOffset[channel] + tdcCut;
double time = data*TDC_BIN;
if(time >= tdcmin && time <= tdcmax) {
h_correctedTDC->Fill((time - tdcOffset[channel]), channel);
hnhitsx[channel]->Fill((pos.xset - OFFSETX) * MIKRO_BIN);
hnhitsy[channel]->Fill((pos.yset - OFFSETY) * MIKRO_BIN);
h2d[channel]->Fill((pos.xset - OFFSETX) * MIKRO_BIN, (pos.yset - OFFSETY) * MIKRO_BIN);
//h_threshold->Fill(channel, thr.threshold);
if (position(pos.xset-OFFSETX, pos.yset-OFFSETY, channel)) {
h_ch33->Fill(pos.xset-OFFSETX, pos.yset-OFFSETY);
h_threshold->Fill(channel, thr.threshold);
}
}
//gV673A->Fill(data,channel);
//gsumV673A[channel/16]->Fill(data);
}
}
if (last_dat) break;
}
} // if (recid== 0x140 || recid== 0x141)
ii += len;
} //while
 
// events ------------------------------------------------------------------------------------------
// fill histograms
/*
for(int i=0; i<NCH; i++) {
//tdc=((double)evtrec->data[i])*TDC_BIN - tdcoffset[i];
//adc=(double)evtrec->data[i+NCH]-adcoffset[i];
if (gNtWrite) nt->Fill(i,pos.ix,pos.iy,tdc);
//if( (qdcmi < adc) && (adc < qdcma) ) {
htdc[i]->Fill(tdc);
//}
//hadc[i]->Fill(adc);
//hcor[i]->Fill(adc,tdc);
//if(adc > corpar[i][2])
//hctdc[i]->Fill( tdc - (corpar[i][0] + corpar[i][1]/TMath::Sqrt(adc - corpar[i][2])) );
 
if( (abs(padcenter[i][0] - pos.xset) < POSMARG) && (abs(padcenter[i][1] - pos.yset) < POSMARG) ) {
htdcpos[i]->Fill(tdc);
//hadcpos[i]->Fill(adc);
//hcorpos[i]->Fill(adc,tdc);
//if(adc > corpar[i][2])
//hctdcpos[i]->Fill( tdc - (corpar[i][0] + corpar[i][1]/TMath::Sqrt(adc - corpar[i][2])) );
}
if((tdcmi < tdc) && (tdc < tdcma)) {
//hadc_cut[i]->Fill(adc);
hnhitsx[i]->Fill((pos.xset - OFFSETX) * MIKRO_BIN);
hnhitsy[i]->Fill((pos.yset - OFFSETY) * MIKRO_BIN);
h2d[i]->Fill((pos.xset - OFFSETX) * MIKRO_BIN, (pos.yset - OFFSETY) * MIKRO_BIN);
//if(i==4) adc1=adc;
//if(i==5) adc2=adc;
}
} */
break;
case THRREC_ID:
status = gzread(dfp, (voidp)&readbuf[2], (rec_len-2*ulsize));
thrrec = (THRREC*) readbuf;
thr = *thrrec;
if (dbg) printf("THRREC id = %d len = %d threshold %d\n",
thrrec->id, thrrec->len, thrrec->threshold);
break;
case ENDREC_ID:
gzread(dfp, (voidp)&readbuf[2], (rec_len-2*ulsize));
endrec = (ENDREC *) readbuf;
if(dbg) {
printf("ENDREC_ID\n");
printf("id = %d, len = %d, time = %d\n", endrec->id, endrec->len, endrec->time);
} else printf(" ENDREC\n");
fcount++;
switch(ftype-1) {
case 0:
sprintf(dfile, "./data/%s_file%02d.dat", dfile0, fcount);
break;
case 1:
sprintf(dfile, "./data/%s_file%02d.dat.gz", dfile0, fcount);
break;
case 2:
sprintf(dfile, "./data/%s_file%02d.gz", dfile0, fcount);
break;
}
if(dfp) gzclose(dfp);
dfp=gzopen(dfile,"rb");
if(!dfp) {
printf(" Cannot open data file: %s ---> Exiting\n", dfile);
end_of_file = 1;
} else {
printf(" Opened data file: %s\n", dfile);
end_of_file = 0;
}
break;
default:
printf("switch(rec_id): default !!!\n");
end_of_file = 1;
break;
}
ceve++;
if( (ceve%50000) == 0) printf(" Current event = %d\n", ceve);
//if(dbg) if( ceve>dbg ) break;
if(end_of_file) break;
}
 
if(dfp) {
gzclose(dfp);
delete dfp;
}
if(dbg) return 1;
if(rootfile) {
nt->Write();
rootfile->Write();
printf("Saved to %s\n", fnameroot);
rootfile->Close();
delete rootfile;
}
return 1;
}
 
int position(int x, int y, int channel)
{
int flag = 0;
if ( (x > (padCenter[channel][0] - 5000)) && (x < (padCenter[channel][0] + 5000)) &&
(y > (padCenter[channel][1] - 5000)) && (y < (padCenter[channel][1] + 5000)) ) flag = 1;
return flag;
}
//////////////////////////////////////////
// Data to root conversion root script
//
// Contributors: Rok Pestotnik, Rok Dolenec, Dino Tahirovic
//
// 3/1/2014 TDC cut relative to tdc offset
// 5/3/2014 TH3F h_correctedTDC; commented 'delete' commands
// 12/3/2014 declarations of all histos before while{}
 
#include "stdio.h"
#include "TROOT.h"
#include "TFile.h"
#include "TNtuple.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TH3F.h"
#include "TF1.h"
#include "TMath.h"
#include "TStyle.h"
#include "TCanvas.h"
#include "TLine.h"
#include "zlib.h"
 
// ------------------------------------------------------------------------------
 
#define POSMARG 1000
 
#define READBUFFERLENGTH 10000
 
// data format
#define MAXDATA 16
#define NCH 64
#define TDC_BIN 1.0416 // 1 TDC bin in ns
#define MIKRO_BIN 0.49609/1000. //1 mikro step in mm; stage MM3MF
#define OFFSETX 5200 // Right edge of SiPM+Lightguide
#define OFFSETY 5400 // Lower edge of SiPM+Lightguide
 
#define RUNREC_ID 1
#define ENDREC_ID 2
#define POSREC_ID 3
#define EVTREC_ID 4
#define THRREC_ID 5
 
typedef struct {
unsigned int id,len;
unsigned int fver,time;
unsigned int nev,nch,ped,xy;
unsigned int thLow, thUp, thStep;
int nx,x0,dx,ny,y0,dy;
} RUNREC;
RUNREC *runrec;
RUNREC run;
 
typedef struct {
unsigned int id,len;
unsigned int time;
} ENDREC;
ENDREC *endrec;
 
typedef struct {
unsigned int id,len;
unsigned int time;
int ix;
int x;
int xset;
int iy;
int y;
int yset;
} POSREC;
POSREC *posrec;
POSREC pos;
 
typedef struct {
unsigned int id;
unsigned int len;
unsigned int nev;
} EVTREC;
EVTREC *evtrec;
 
typedef struct {
unsigned int id;
unsigned int len;
unsigned int threshold;
} THRREC;
THRREC *thrrec;
THRREC thr;
 
double padCenter[NCH][2];
 
int position(int, int, int);
 
// ------------------------------------------------------------------------------
 
int d2r(char* dfile0="test", int dbg=0, double tdcCut=5.0)
{
const double c_tdcOffset = 1; // ns
printf(" Data to root conversion program\nUsage:\nd2r(input file name <without.dat>, debug on/off, TDC cut +-[ns])\n\n");
char fullname[256];
char sbuff[256];
FILE *fp;
//Chanel information
double tdcOffset[NCH];
sprintf(fullname, "d2r.ini");
if( (fp=fopen(fullname, "rt")) == NULL )
printf("Cannot open pad centers file %s !!!\n", fullname);
else {
printf("Opened pad centers file %s\n", fullname);
char* result = fgets(sbuff,256, fp);
if (dbg) printf("Read buffer %s\n", result);
printf("%s", sbuff);
for(int i=0; i<NCH; i++) {
int channel;
int message = fscanf(fp, "%d %lf %lf %lf\n", &channel, &padCenter[i][0], &padCenter[i][1], &tdcOffset[i]);
if (dbg) printf("Read d2r.ini returned %d\n", message);
}
fclose(fp);
}
for(int i=0; i<NCH; i++) {
tdcOffset[i] += c_tdcOffset;
printf("%.2lf %.2lf %.2lf\n", padCenter[i][0], padCenter[i][1], tdcOffset[i]);
}
//data buffer
unsigned int readbuf[READBUFFERLENGTH];
unsigned int buf[READBUFFERLENGTH];
//data file
gzFile dfp;
char dfile[256];
int ftype=0;
int fcount=1;
do {
switch(ftype++) {
case 0:
sprintf(dfile, "./data/%s_file%02d.dat", dfile0, fcount);
break;
case 1:
sprintf(dfile, "./data/%s_file%02d.dat.gz", dfile0, fcount);
break;
case 2:
sprintf(dfile, "./data/%s_file%02d.gz", dfile0, fcount);
break;
default:
printf(" Cannot find data file for %s !!!\n", dfile0);
return -1;
}
dfp=gzopen(dfile,"rb");
} while(!dfp);
printf("Opened data file %s\n", dfile);
//opens ROOT file
char fnameroot[256];
sprintf(fnameroot, "root/%s.root", dfile0);
TFile rootfile(fnameroot,"RECREATE",dfile0);
if (rootfile.IsZombie()) {
std::cout << "Error opening file" << std::endl;
exit(-1);
}
if (rootfile.IsOpen()) std::cout << "ROOT file opened for writing." << std::endl;
// -----------------------------------------------
// loop trough records
unsigned int rec_id, rec_len;
unsigned int ulsize = sizeof(unsigned int);
//unsigned int ulsize = 4;
if (dbg) printf("Size of unsigned int: %lu\n", sizeof(unsigned int));
int ceve=0;
int end_of_file = 0;
int ii;
int nint;
int nb;
int status;
char hname[256];
TH2F* htdc = new TH2F("htdc",";TDC channel;SiPM channel",512,-0.5,511.5,NCH,-0.5,NCH-0.5);
TH3F* h_correctedTDC = new TH3F("h_correctedTDC",";SiPM channel; ASD threshold [V]; t [ns]",
NCH, -0.5, NCH-0.5,
101, 1.0, 2.0,
33, -16.5*TDC_BIN, 16.5*TDC_BIN);
TH1F* hnhitsx[NCH]; // move to 2d with (channel, position)
TH1F* hnhitsy[NCH]; //-||-
TH2F* h2d[NCH]; //-||-
TH2F* h_threshold = new TH2F("h_threshold",";SiPM ch;Threshold[V]",
64,-0.5,63.5,
101, 1.0, 2.0);
TH2F* h_ch33 = new TH2F("h_ch33","ch. 33;x;y",100,20000,30000,100,0,10000);
TNtuple* nt = new TNtuple("nt", "nt", "ch:x:y:tdc");
while(1) {
if(gzeof(dfp)) end_of_file = 1;
 
gzread(dfp, (voidp)&readbuf, 2*ulsize);
rec_id=readbuf[0];
rec_len=readbuf[1];
if(dbg) printf("-----------------------------------------------\n");
if(dbg) printf("[%d] rec_id = %d | rec_len = %u\n", ceve, rec_id, rec_len);
double nSteps = 0;
switch(rec_id)
{
case RUNREC_ID:
if (dbg) printf("RUNREC\n");
gzread(dfp, (voidp)&readbuf[2], (rec_len-2*ulsize));
runrec = (RUNREC *) readbuf;
run=*runrec;
if(dbg) {
printf("RUNREC_ID\n");
printf("id = %d, len = %d, time = %d\n", run.id, run.len, run.time);
printf("nev = %d, nch = %d\n", run.nev, run.nch);
printf("nx = %d, x0 = %d, dx = %d\n", run.nx, run.x0, run.dx);
printf("ny = %d, y0 = %d, dy = %d\n", run.ny, run.y0, run.dy);
printf("thLow = %d, thUp = %d, thStep = %d\n", run.thLow, run.thUp, run.thStep);
}
//create histograms
//nt = (TNtuple*) gROOT->FindObject("nt");
//if (nt) delete nt;
//sprintf(hname, "htdc");
//htdc = (TH2F*) gROOT->FindObject(hname);
//if (htdc) delete htdc;
//htdc = new TH2F("htdc",";TDC channel;SiPM channel",512,0,512,NCH,0,NCH);
h_correctedTDC = (TH3F*) gROOT->FindObject("h_correctedTDC");
if (h_correctedTDC) delete h_correctedTDC;
nSteps = (run.thUp - run.thLow)/double(run.thStep) + 1;
if (dbg) printf("nSteps %d\n", nSteps);
h_correctedTDC = new TH3F("h_correctedTDC",";SiPM channel; ASD threshold [V]; t [ns]",
NCH,
-0.5,
NCH-0.5,
nSteps,
(run.thLow - 0.5*run.thStep)/1000.0,
(run.thUp + 0.5*run.thStep)/1000.0,
2*tdcCut*TDC_BIN,
-tdcCut*TDC_BIN,
tdcCut*TDC_BIN);
//TH1F* gsumV673A[NCH/16] = new TH1F(hn,hname,256,-0.5,255.5);
for(int i=0; i<NCH; i++) {
/*
sprintf(hname, "htdcpos%d", i);
htdcpos[i] = (TH1F*)gROOT->FindObject(hname);
if(htdcpos[i]) delete htdcpos[i];
htdcpos[i] = new TH1F(hname, hname, 512, 0, 512*TDC_BIN);
sprintf(hname, "htdc%d", i);
htdc[i] = (TH1F*)gROOT->FindObject(hname);
if(htdc[i]) delete htdc[i];
htdc[i] = new TH1F(hname, hname, 512, 0*TDC_BIN, 512*TDC_BIN);
*/
sprintf(hname, "hnhitsx%d", i);
hnhitsx[i] = (TH1F*)gROOT->FindObject(hname);
if(hnhitsx[i]) delete hnhitsx[i];
hnhitsx[i] = new TH1F(hname, hname, run.nx,
(run.x0 - OFFSETX - 0.5*run.dx)*MIKRO_BIN,
(run.x0 - OFFSETX + (run.nx-0.5)*run.dx)*MIKRO_BIN);
sprintf(hname, "hnhitsy%d", i);
hnhitsy[i] = (TH1F*)gROOT->FindObject(hname);
if(hnhitsy[i]) delete hnhitsy[i];
//hnhitsy[i] = new TH1F(hname, hname, run.ny, (run.y0-0.5*run.dy)*MIKRO_BIN, (run.y0+(run.ny-0.5)*run.dy)*MIKRO_BIN);
hnhitsy[i] = new TH1F(hname, hname, run.ny,
(run.y0 - 0.5*run.dy - OFFSETY)*MIKRO_BIN,
(run.y0 + (run.ny-0.5)*run.dy - OFFSETY)*MIKRO_BIN);
//hnhitsy[i] = new TH1F(hname, hname, 100, 0,100);
sprintf(hname, "h2d%d", i);
h2d[i] = (TH2F*)gROOT->FindObject(hname);
if(h2d[i]) delete h2d[i];
h2d[i] = new TH2F(hname, hname, run.nx,
(run.x0 - OFFSETX - 0.5*run.dx)*MIKRO_BIN,
(run.x0 - OFFSETX + (run.nx-0.5)*run.dx)*MIKRO_BIN,
run.ny,
(run.y0 - OFFSETY - 0.5*run.dy)*MIKRO_BIN,
(run.y0 - OFFSETY + (run.ny-0.5)*run.dy)*MIKRO_BIN);
}
if (h_threshold) delete h_threshold;
h_threshold = new TH2F("h_threshold",";SiPM ch;Threshold[V]",64,-0.5,63.5,
nSteps,
(run.thLow - 0.5*run.thStep)/1000.0,
(run.thUp + 0.5*run.thStep)/1000.0);
//h_threshold = new TH2F("h_threshold","Threshold scan;SiPM ch;Threshold[mV]",64,-0.5,63.5,
// 101,995,2005);
if (h_ch33) delete h_ch33;
h_ch33 = new TH2F("h_ch33","ch. 33;x;y",100,20000,30000,100,0,10000);
 
if (dbg) printf("RUNREC: all histos created.\n");
break;
case POSREC_ID:
if (dbg) printf("POSREC\n");
gzread(dfp, (voidp)&readbuf[2], (rec_len-2*ulsize));
posrec = (POSREC *) readbuf;
pos=*posrec;
if(dbg) {
printf("POSREC_ID\n");
printf("id = %d, len = %d, time = %d\n", posrec->id, posrec->len, posrec->time);
printf("ix = %d, x = %d, xset = %d\n", posrec->ix, posrec->x, posrec->xset);
printf("iy = %d, y = %d, yset = %d\n", posrec->iy, posrec->y, posrec->yset);
} else printf(" [%d,%d] %d, %d\n", pos.ix, pos.iy, pos.xset, pos.yset);
break;
case EVTREC_ID:
gzread(dfp, (voidp)&readbuf[2], ulsize); // last field of event record
evtrec = (EVTREC *) readbuf;
//evtrec->nev = buf[0];
//if (rec_len < 0 || rec_len > 10000) {
if (rec_len > READBUFFERLENGTH) {
printf("Len %u\n", rec_len);
return(1);
}
nb = rec_len - 3*ulsize; // no. of bytes to read
gzread(dfp, (voidp)&buf, nb);
if(dbg) {
printf("EVTREC_ID\n");
printf("id = %d, len = %d, nev = %d\n", evtrec->id, evtrec->len, evtrec->nev);
//for(int datai = 0; datai < NCH; datai++) printf("%u ", evtrec->data[datai]);
//printf("\n");
//for(int datai = NCH; datai < NCH+NCH; datai++) printf("%u ", evtrec->data[datai]);
//printf("\n");
break;
}
nint = nb / ulsize; // no. of subrecords
if (dbg) printf("No. of subrecords %d \n", nint);
ii=0;
while (ii<nint){
int recid = buf[ii++];
int len = buf[ii++];
if (dbg) printf("Buffer pointer %d\n", ii);
unsigned int *dbuf = (unsigned int *)&buf[ii];
//if (n%1000==0)
if (dbg) printf("%d 0x%03x Len=%d\n",evtrec->nev,recid,len);
//unsigned short edge;
//int nhits;
if (recid==0x140 || recid==0x141) {
for (int i=0; i<len; i++) {
int data = dbuf[i] & 0xFFFF ;
int edge_type = (dbuf[i]>>16)&0x1 ;
int overflow = (dbuf[i]>>17)&0x1 ;
int tdc_num = (dbuf[i]>>25)&0x1 ;
int channel = ((dbuf[i]>>18)&0x1F) | tdc_num<<5 ;
int ev_dat = (dbuf[i]>>23)&0x1 ;
int last_dat = (dbuf[i]>>30)&0x1 ;
int nval_dat = (dbuf[i]>>31)&0x1 ;
if (dbg){
if (ev_dat) printf("Event %d\n",data);
else printf("ch=%d edge=%d ev=%d data=%d last=%d nval=%d\n",
channel, edge_type,ev_dat,data,last_dat,nval_dat);
}
if (!ev_dat){
if (!edge_type && !overflow) {
htdc->Fill(data, channel);
if (dbg) printf("ch: %d tdc: %d\n", channel, data);
if (dbg) nt->Fill(channel, pos.ix, pos.iy, data);
double tdcmin = tdcOffset[channel] - tdcCut;
double tdcmax = tdcOffset[channel] + tdcCut;
double time = data*TDC_BIN - tdcOffset[channel];
if(time >= -tdcCut and time <= tdcCut) {
h_correctedTDC->Fill(channel, thr.threshold/1000.0, time);
hnhitsx[channel]->Fill((pos.xset - OFFSETX) * MIKRO_BIN);
hnhitsy[channel]->Fill((pos.yset - OFFSETY) * MIKRO_BIN);
h2d[channel]->Fill((pos.xset - OFFSETX) * MIKRO_BIN, (pos.yset - OFFSETY) * MIKRO_BIN);
 
if (position(pos.xset-OFFSETX, pos.yset-OFFSETY, channel)) {
h_ch33->Fill(pos.xset-OFFSETX, pos.yset-OFFSETY);
h_threshold->Fill(channel, thr.threshold/1000.0);
}
}
//gV673A->Fill(data,channel);
//gsumV673A[channel/16]->Fill(data);
}
}
if (last_dat) break;
}
} // if (recid== 0x140 || recid== 0x141)
ii += len;
} //while
 
break;
case THRREC_ID:
status = gzread(dfp, (voidp)&readbuf[2], (rec_len-2*ulsize));
thrrec = (THRREC*) readbuf;
thr = *thrrec;
if (dbg) printf("THRREC id = %d len = %d threshold %d\n",
thrrec->id, thrrec->len, thrrec->threshold);
break;
case ENDREC_ID:
gzread(dfp, (voidp)&readbuf[2], (rec_len-2*ulsize));
endrec = (ENDREC *) readbuf;
if(dbg) {
printf("ENDREC_ID\n");
printf("id = %d, len = %d, time = %d\n", endrec->id, endrec->len, endrec->time);
} else printf(" ENDREC\n");
fcount++;
switch(ftype-1) {
case 0:
sprintf(dfile, "./data/%s_file%02d.dat", dfile0, fcount);
break;
case 1:
sprintf(dfile, "./data/%s_file%02d.dat.gz", dfile0, fcount);
break;
case 2:
sprintf(dfile, "./data/%s_file%02d.gz", dfile0, fcount);
break;
}
if(dfp) gzclose(dfp);
dfp=gzopen(dfile,"rb");
if(!dfp) {
printf(" Cannot open data file: %s ---> Exiting\n", dfile);
end_of_file = 1;
} else {
printf(" Opened data file: %s\n", dfile);
end_of_file = 0;
}
break;
default:
printf("switch(rec_id): default !!!\n");
end_of_file = 1;
break;
}
ceve++;
if ( (ceve%50000) == 0) printf(" Current event = %d\n", ceve);
//if(dbg) if( ceve>dbg ) break;
if (end_of_file) break;
}
 
if(dfp) {
gzclose(dfp);
//delete dfp;
}
//if(dbg) return 0;
if(rootfile.IsOpen()) {
nt->Write();
rootfile.Write();
printf("Saved to %s\n", fnameroot);
rootfile.Close();
//delete rootfile;
}
return 0;
}
 
int position(int x, int y, int channel)
{
int flag = 0;
if ( (x > (padCenter[channel][0] - 5000)) and (x < (padCenter[channel][0] + 5000)) and
(y > (padCenter[channel][1] - 5000)) and (y < (padCenter[channel][1] + 5000)) ) flag = 1;
return flag;
}
/sipmScan/src/analysisScan.cpp
31,7 → 31,7
#define TDC_BIN 1.0416 // 1 TDC bin in ns
#define MIKRO_BIN 0.49609/1000. //1 mikro step in mm; stage MM3MF
#define OFFSETX 5200 // Right edge of SiPM+Lightguide
#define OFFSETY 4800 // Lower edge of SiPM+Lightguide
#define OFFSETY 5300 // Lower edge of SiPM+Lightguide
 
#define RUNREC_ID 1
#define ENDREC_ID 2
371,11 → 371,11
if(time >= tdcmin && time <= tdcmax) {
h_correctedTDC->Fill((time - tdcOffset[channel]), channel);
hnhitsx[channel]->Fill((pos.xset - OFFSETX) * MIKRO_BIN);
hnhitsy[channel]->Fill((pos.yset - OFFSETY) * MIKRO_BIN);
h2d[channel]->Fill((pos.xset - OFFSETX) * MIKRO_BIN, (pos.yset - OFFSETY) * MIKRO_BIN);
hnhitsy[channel]->Fill((pos.yset - OFFSETY) * MIKRO_BIN);
//h_threshold->Fill(channel, thr.threshold);
//if (position(pos.xset-OFFSETX, pos.yset-OFFSETY, channel)) {
//}
if (position(pos.xset-OFFSETX, pos.yset-OFFSETY, channel)) {
h2d[channel]->Fill((pos.xset - OFFSETX) * MIKRO_BIN, (pos.yset - OFFSETY) * MIKRO_BIN);
}
}
//gV673A->Fill(data,channel);
//gsumV673A[channel/16]->Fill(data);
/sipmScan/examples/thresholdScan.C
0,0 → 1,513
#include "TROOT.h"
#include "TFile.h"
#include "TBenchmark.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TH3F.h"
#include "TCanvas.h"
#include "TStyle.h"
#include "TPad.h"
#include "TF1.h"
#include "TGraph.h"
#include "TSpectrum.h"
#include "stdio.h"
 
#include "include/RTUtil.h"
 
double getNoise(TH2F*, int, int);
 
int thresholdScan(char filename[256] = "test", char plopt[256]="th", int chXstart=0, int chXend=7, int chYstart=0, int chYend=7, bool debug = false)
{
const int c_nChannels = 64;
const double c_xOffset = 1; // mm
const double c_yOffset = 0.7;
int map[8][8]={{32,34,53,55,40,42,61,63},
{48,50,37,39,56,58,45,47},
{33,35,52,54,41,43,60,62},
{49,51,36,38,57,59,44,46},
{17,19,4,6,25,27,12,14},
{1,3,20,22,9,11,28,30},
{16,18,5,7,24,26,13,15},
{0,2,21,23,8,10,29,31}
};
RTSetStyle(gStyle);
// open the file with histograms
char fnameroot[256];
TFile* rootfile;
sprintf(fnameroot, "root/%s.root", filename);
rootfile = (TFile *) gROOT->FindObject(filename);
if(rootfile==NULL) rootfile = new TFile(fnameroot);
if(rootfile==NULL) {
printf("Cannot open root file %s!!!\n",fnameroot);
return(0);
}
 
if (strstr(plopt, "th") != NULL) {
TCanvas* canvas1 = new TCanvas("canvas1","",800,800);
canvas1->cd();
//gPad->SetLogz();
//canvas1->SetLogz();
TH2F* h_threshold = (TH2F*) rootfile->Get("h_threshold");
h_threshold->SetTitleOffset(1.5,"y");
//h_threshold->GetZaxis()->SetRangeUser(0,1000);
h_threshold->SetContour(20);
h_threshold->Draw("colz");
TCanvas* canvas2 = new TCanvas("canvas2","",1600,800);
canvas2->Divide(2);
TH1D* h_projection1 = h_threshold->ProjectionY("Ch 37",38,38);
canvas2->cd(1);
//h_projection1->GetYaxis()->SetRangeUser(0,1000);
//gPad->SetLogy();
h_projection1->Draw();
TH1D* h_projection2 = h_threshold->ProjectionY("Ch 38",39,39);
canvas2->cd(2);
//gPad->SetLogy();
h_projection2->Draw();
}
if(strstr(plopt, "tdc") != NULL) {
TCanvas *canvas21 = new TCanvas("canvas21","canvas21",1600,800);
TH3F* h0 = (TH3F*) rootfile->Get("h_correctedTDC");
Int_t binsX = h0->GetXaxis()->GetNbins();
Int_t minX = h0->GetXaxis()->GetFirst();
Int_t maxX = h0->GetXaxis()->GetLast()+1;
Int_t binsY = h0->GetYaxis()->GetNbins();
Int_t minY = h0->GetYaxis()->GetFirst();
Int_t maxY = h0->GetYaxis()->GetLast()+1;
Int_t binsZ = h0->GetZaxis()->GetNbins();
Int_t minZ = h0->GetZaxis()->GetFirst();
Int_t maxZ = h0->GetZaxis()->GetLast()+1;
Double_t xLowUser = h0->GetXaxis()->GetBinLowEdge(minX);
Double_t xUpUser = h0->GetXaxis()->GetBinUpEdge(maxX);
Double_t yLowUser = h0->GetYaxis()->GetBinLowEdge(minY);
Double_t yUpUser = h0->GetYaxis()->GetBinUpEdge(maxY);
canvas21->Divide(2);
//h0->Draw();
canvas21->cd(1);
//h0->Project3D("xy")->Draw("colz");
//canvas21->cd(2);
//h0->Project3D("yz")->Draw("colz");
//canvas21->cd(3);
TH2D* h_correctedTDC = (TH2D*) h0->Project3D("xz");
//h0->GetZaxis()->SetRangeUser(-5,5);
//h0->Project3D("xzo")->Draw("colz");
h_correctedTDC->SetTitle("; t [ns]; Channel");
h_correctedTDC->Draw("colz");
/*
TH1F* tdc1 = new TH1F("tdc1",";TDC [ns];Events",binsZ, minZ, maxZ);
//for(int j=minY; j<maxY; j++) {
for(int k=minZ; k<maxZ; k++) {
double signal = h0->GetBinContent(1,1,k);
tdc1->Fill(k,signal);
//}
}
tdc1->Draw();
*/
TH2D* h_timeWalk = (TH2D*) h0->Project3D("zy");
canvas21->cd(2);
h_timeWalk->SetTitle(";Threshold [V]; t [ns]");
h_timeWalk->Draw("colz");
}
if(strstr(plopt,"16") != NULL)
{
TH2F* h_threshold = (TH2F*) rootfile->Get("h_threshold");
TCanvas* canvas = new TCanvas("canvas","",1600,1600);
canvas->cd();
TVirtualPad *main = new TPad("main","main",0,0,1,1,10,1);
main->Draw();
main->Divide(4,4);
for(int i=chXstart; i<chXend; i++) {
for(int j=chYstart; j<chYend; j++) {
int channel = map[i][j];
TH1D* h_projection = h_threshold->ProjectionY("",channel+1,channel+1);
int canvasPosition = i-chXend+4*(chYend-j)+1;
printf(" %d ", canvasPosition);
main->cd(canvasPosition);
char name[128];
sprintf(name,"Channel %d",channel);
h_projection->SetTitle(name);
h_projection->DrawCopy();
}
}
printf("\n");
}
if( strstr(plopt, "all") != NULL ) {
TCanvas *canvas2 = new TCanvas("canvas2","Hits x;;",2000,2000);
TCanvas *canvas3 = new TCanvas("canvas3","Hits y;;",2000,2000);
canvas2->Divide(8,8);
canvas3->Divide(8,8);
TH1F* h_hitsx;
TH1F* h_hitsy;
for(int i=0; i<c_nChannels; i++) {
canvas2->cd(i+1);
char hname[128];
sprintf(hname, "hnhitsx%d", i);
h_hitsx = (TH1F*)rootfile->Get(hname);
h_hitsx->Draw();
canvas3->cd(i+1);
sprintf(hname, "hnhitsy%d", i);
h_hitsy = (TH1F*)rootfile->Get(hname);
h_hitsy->Draw();
}
}
if( strstr(plopt, "x") != NULL ) {
TCanvas *canvas10 = new TCanvas("canvas10","Ch x;;",500,500);
TH1F* h_hitsx;
canvas10->cd();
char hname[128];
sprintf(hname, "hnhitsx%d", chXstart);
h_hitsx = (TH1F*)rootfile->Get(hname);
h_hitsx->Draw("colz");
}
if( strstr(plopt, "y") != NULL ) {
TCanvas *canvas11 = new TCanvas("canvas11","Ch x;;",500,500);
TH1F* h_hitsy;
canvas11->cd();
char hname[128];
sprintf(hname, "hnhitsy%d", chXstart);
h_hitsy = (TH1F*)rootfile->Get(hname);
h_hitsy->Draw();
}
if( strstr(plopt, "share") != NULL ) {
/*TCanvas *canvas4 = new TCanvas("canvas1","canvas1",1000,1000);
int nChannels = chYend-chYstart+1;
int ncols = nChannels/2;
printf("nch %d nch\\2 %d\n", nChannels, ncols);
canvas4->Divide(2,ncols);
TH1F* h_hitsy;
for(int i=chYstart; i<=chYend; i++){
canvas4->cd(i-chYstart+1);
char hname[128];
int chPosition = map[0][i];
sprintf(hname, "hnhitsy%d", chPosition);
h_hitsy = (TH1F*)rootfile->Get(hname);
h_hitsy->Draw();
}*/
TCanvas *canvas4 = new TCanvas("canvas4","canvas4",500,500);
canvas4->cd();
for(int i=chXstart; i<=chXend; i++) {
TH1F* h_hitsx;
char hname[128];
int chPosition = map[i][chYstart];
sprintf(hname, "hnhitsx%d", chPosition);
h_hitsx = (TH1F*)rootfile->Get(hname);
h_hitsx->SetTitle("Scan X;x [mm]; Entries");
h_hitsx->GetYaxis()->SetTitleOffset(1.3);
h_hitsx->SetStats(0);
if (i == chXstart)
h_hitsx->Draw();
else {
h_hitsx->SetLineColor(i+1);
h_hitsx->Draw("same");
}
}
//sprintf(fullname, "ps/%s_Yshare.eps", filename);
//canvas4->SaveAs(fullname);
TCanvas *canvas5 = new TCanvas("canvas5","canvas5",500,500);
canvas5->cd();
for(int i=chYstart; i<=chYend; i++) {
TH1F* h_hitsy;
char hname[128];
int chPosition = map[chXstart][i];
sprintf(hname, "hnhitsy%d", chPosition);
h_hitsy = (TH1F*)rootfile->Get(hname);
h_hitsy->SetTitle("Scan Y;y [mm]; Entries");
h_hitsy->GetYaxis()->SetTitleOffset(1.3);
h_hitsy->SetStats(0);
if (i == chYstart)
h_hitsy->Draw();
else {
h_hitsy->SetLineColor(i+1);
h_hitsy->Draw("same");
}
}
//sprintf(fullname, "ps/%s_Yshare.eps", filename);
//canvas5->SaveAs(fullname);
}
/** Draws the signal from 8 channels in x-row
* for one specific y bin, so the background and cross-talk
* can be estimated.
* Draws also a 2d scan of these channels.
*/
if (strstr(plopt, "line") != NULL) {
TCanvas* canvas6 = new TCanvas("canvas6","canvas6",500,500);
canvas6->cd(0);
gStyle->SetOptStat(0);
TH2F* h0 = (TH2F*) rootfile->Get("h2d0");
Int_t binsX = h0->GetXaxis()->GetNbins();
Int_t minX = h0->GetXaxis()->GetFirst();
Int_t maxX = h0->GetXaxis()->GetLast()+1;
Int_t binsY = h0->GetYaxis()->GetNbins();
Int_t minY = h0->GetYaxis()->GetFirst();
Int_t maxY = h0->GetYaxis()->GetLast()+1;
Double_t xLowUser = h0->GetXaxis()->GetBinLowEdge(minX);
Double_t xUpUser = h0->GetXaxis()->GetBinUpEdge(maxX);
Double_t yLowUser = h0->GetYaxis()->GetBinLowEdge(minY);
Double_t yUpUser = h0->GetYaxis()->GetBinUpEdge(maxY);
if (debug) printf("xLow %f xUp %f\n",xLowUser,xUpUser);
//! 1-dimension position in x vs. hits
TH2F* h[8];
TH1F* h_line[8];
for(int j=0; j<8; j++) {
h_line[j] = new TH1F("h_line", "h_line", binsX, xLowUser, xUpUser);
}
for(int j=chXstart; j<=chXend; j++) {
int chPosition = map[j][chYstart];
char hname[128];
sprintf(hname, "h2d%d", chPosition);
int histogram = j;
h[histogram] = (TH2F *) rootfile->Get(hname);
int noise = getNoise(h[histogram], 1, 160);
for(int k=minX; k<=maxX; k++) {
int l=chYstart*20+12;
//for(int l=12; l<=16; l++) {
double signal = h[histogram]->GetBinContent(k,l);
//signal -= noise;
//signal /= 5*10000.0;
double eta = -log(1 - signal);
double x = k*(xUpUser-xLowUser)/double(binsX);
//double y = l*(yUpUser-yLowUser)/double(binsY);
h_line[j]->Fill(x-c_xOffset, signal);
//}
}
if (j == chXstart) {
h_line[j]->SetTitle("SiPM#2 w/o noise subtraction;x[mm];Hits");
//h_line[j]->GetYaxis()->SetRangeUser(-0.05, 0.3);
//h_line[j]->GetYaxis()->SetRangeUser(-50, 2500);
h_line[j]->Draw("");
}
else {
h_line[j]->SetLineColor(j+1);
h_line[j]->Draw("same");
}
}
//! 2d scan
TCanvas* canvas61 = new TCanvas("canvas61","canvas61",8*200,300);
canvas61->cd();
TH2F* h_corrected = new TH2F("h_corrected","h_corrected",binsX,xLowUser,xUpUser, binsY,yLowUser,yUpUser);
for(int i=chXstart; i<=chXend; i++) {
//int canvasPosition = nX*(i-chYstart)+(j-chXstart)+1;
//int canvasPosition = nX*(chYend-i)+chXstart+1;
//if (debug) printf("canvas %d\n",canvasPosition);
int chPosition = map[i][chYstart];
char hname[128];
sprintf(hname, "h2d%d", chPosition);
int histogram = i;
h[histogram] = (TH2F *) rootfile->Get(hname);
int noise = getNoise(h[histogram], 1, 100);
for(int k=minX; k<=maxX; k++) {
for(int l=minY; l<=maxY; l++) {
int signal = h[histogram]->GetBinContent(k,l); // detected
//p /= 10000.;
//double p0 = 1.0 - p; // events with zero photons
//double eta = (-log(p0) * p0 ) / (1-p0-0.00001);
//double eta = -log(p0);
//printf("p=%f p0=%f log(p0)=%f eta=%f\n",p,p0,log(p0),eta);
//double signal = ((p - noise) > 0.1) ? (p-noise) : 0.1;
double p = signal - noise;
p /= 10000.0;
double eta = -log(1 - p);
double x = k*(xUpUser-xLowUser)/double(binsX);
double y = l*(yUpUser-yLowUser)/double(binsY);
h_corrected->Fill(x-c_xOffset, y-c_yOffset, eta);
}
}
}
h_corrected->SetTitle("SiPM#2 n_pe = - ln(P0);x[mm];y[mm]");
h_corrected->GetZaxis()->SetRangeUser(-0.05,0.30);
h_corrected->Draw("colz");
}
/** Draws the sum of the channels
* Each channel is a 2d plot
* Intended for the study of 1 channel
*/
if (strstr(plopt, "2d") != NULL) {
int nX = chXend - chXstart + 1;
int nY = chYend - chYstart + 1;
TCanvas* canvas7 = new TCanvas("canvas7","canvas7", nX*400,nY*400);
printf("nx %d ny %d\n",nX,nY);
canvas7->Divide(nX,nY);
for(int i=chYstart; i<=chYend; i++) {
for(int j=chXstart; j<=chXend; j++) {
//int canvasPosition = nX*(i-chYstart)+(j-chXstart)+1;
int canvasPosition = nX*(chYend-i)+(j-chXstart)+1;
if (debug) printf("canvas %d\n",canvasPosition);
canvas7->cd(canvasPosition);
char hname[128];
int chPosition = map[j][i];
sprintf(hname, "h2d%d", chPosition);
TH2F* h_2d = (TH2F*)rootfile->Get(hname);
h_2d->Draw("colz");
} //x
}
// Number of photoelectrons - Poissonian correction
TCanvas* canvas8 = new TCanvas("canvas8","canvas8", 1000,1000);
canvas8->cd();
gStyle->SetOptStat(0);
char hname[128];
int chPosition = map[1][2];
sprintf(hname, "h2d%d", chPosition);
TH2F* h_2d = (TH2F*)rootfile->Get(hname);
 
Int_t binsX = h_2d->GetXaxis()->GetNbins();
Int_t minX = h_2d->GetXaxis()->GetFirst();
Int_t maxX = h_2d->GetXaxis()->GetLast()+1;
Int_t binsY = h_2d->GetYaxis()->GetNbins();
Int_t minY = h_2d->GetYaxis()->GetFirst();
Int_t maxY = h_2d->GetYaxis()->GetLast()+1;
Double_t xLowUser = h_2d->GetXaxis()->GetBinLowEdge(minX);
Double_t xUpUser = h_2d->GetXaxis()->GetBinUpEdge(maxX);
Double_t yLowUser = h_2d->GetYaxis()->GetBinLowEdge(minY);
Double_t yUpUser = h_2d->GetYaxis()->GetBinUpEdge(maxY);
if (debug) printf("xLow %f xUp %f\n",xLowUser,xUpUser);
TH2F* h_corrected = new TH2F("h_corrected","h_corrected",binsX, xLowUser, xUpUser, binsY, yLowUser, yUpUser);
double noise = getNoise(h_2d, 1, 89);
if(debug) printf("Noise = %f\n", noise);
for(int k=minX; k<=maxX; k++) {
for(int j=minY; j<=maxY; j++) {
double signal = h_2d->GetBinContent(k,j); // detected
//double p = ((signal - noise) > 1) ? (signal-noise) : 1;
double p = signal - noise;
p /= 10000.;
double p0 = 1.0 - p; // events with zero photons
//double eta = (-log(p0) * p0 ) / (1-p0-0.00001);
double eta = -log(p0); // constant of the poissonian statistics
if (debug) printf("p=%f p0=%f log(p0)=%f eta=%f\n",p,p0,log(p0),eta);
double x = k*(xUpUser-xLowUser)/double(binsX);
double y = j*(yUpUser-yLowUser)/double(binsY);
h_corrected->Fill(x+3,y+8, eta);
}
}
h_corrected->SetTitle("n_pe = - ln P(0);x[mm];y[mm]");
h_corrected->GetZaxis()->SetRangeUser(-0.05,0.3);
h_corrected->Draw("colz");
// collection efficiency
int nPoints =0;
double efficiency=0;
for (int i=18; i<=58; i++) {
for (int j=19; j<=59; j++) {
double signal = h_corrected->GetBinContent(i,j);
if(debug) printf("signal %f\n",signal);
efficiency += signal;
nPoints++;
}
}
printf("Signal sum = %f\n # of points = %d\n",efficiency,nPoints);
}
/** Draws the sum of channel signals
* Each channel is a 2d ('h2d') histogram
* Suitable for 8x8 chs scan
*/
if( strstr(plopt, "sum") != NULL ) {
int nX = chXend - chXstart + 1;
int nY = chYend - chYstart + 1;
TCanvas* canvas12 = new TCanvas("canvas12","c2",8*200, 8*200);
canvas12->cd();
gStyle->SetOptStat(0);
 
// final histogram parameters
TH2F* h0 = (TH2F*) rootfile->Get("h2d0");
Int_t binsX = h0->GetXaxis()->GetNbins();
Int_t minX = h0->GetXaxis()->GetFirst();
Int_t maxX = h0->GetXaxis()->GetLast()+1;
Int_t binsY = h0->GetYaxis()->GetNbins();
Int_t minY = h0->GetYaxis()->GetFirst();
Int_t maxY = h0->GetYaxis()->GetLast()+1;
Double_t xLowUser = h0->GetXaxis()->GetBinLowEdge(minX);
Double_t xUpUser = h0->GetXaxis()->GetBinUpEdge(maxX);
Double_t yLowUser = h0->GetYaxis()->GetBinLowEdge(minY);
Double_t yUpUser = h0->GetYaxis()->GetBinUpEdge(maxY);
if (debug) printf("xLow %f xUp %f\n",xLowUser,xUpUser);
TH2F* h_corrected = new TH2F("h_corrected","h_corrected",binsX,xLowUser,xUpUser, binsY,yLowUser,yUpUser);
TH2F* h[9];
// 2d histogram noise subtraction and poisson scaling
for(int i=chYstart; i<=chYend; i++) {
for(int j=chXstart; j<=chXend; j++) {
int chPosition = map[j][i];
char hname[128];
sprintf(hname, "h2d%d", chPosition);
int histogram = nX*(i-chYstart)+(j-chXstart);
h[histogram] = (TH2F *) rootfile->Get(hname);
int noise = getNoise(h[histogram], 1, 170);
if (debug) printf("noise: %d\n",noise);
for(int k=minX; k<=maxX; k++) {
for(int l=minY; l<=maxY; l++) {
int signal = h[histogram]->GetBinContent(k,l); // detected
//double p = ((signal - noise) > 0.1) ? (signal-noise) : 0.1;
double p = signal - noise;
p /= 10000.;
double p0 = 1.0 - p; // events with zero photons
//double eta = (-log(p0) * p0 ) / (1-p0-0.00001);
double eta = -log(p0);
//printf("p=%f p0=%f log(p0)=%f eta=%f\n",p,p0,log(p0),eta);
double x = k*(xUpUser-xLowUser)/double(binsX);
double y = l*(yUpUser-yLowUser)/double(binsY);
h_corrected->Fill(x-c_xOffset, y-c_yOffset, signal);
}
}
}
}
h_corrected->SetTitle("SiPM#2 n_p.e.;x[mm];y[mm]");
//h_corrected->GetZaxis()->SetRangeUser(-0.05,.30);
h_corrected->Draw("colz");
TCanvas* canvas13 = new TCanvas("canvas13","canvas13",600,300);
canvas13->Divide(2);
canvas13->cd(1);
//h[16]->Draw("colz");
canvas13->cd(2);
//h[8]->Draw("colz");
}
return(0);
}
 
/** Function calculates the noise from one channel
* it runs through the bins along x and returns the average value
*/
double getNoise(TH2F* histogram, int yStart, int yEnd)
{
double noise=0;
int count=0;
for(int j=yStart; j<yEnd; j++) {
double value = histogram->GetBinContent(j,2);
//if (noise < value) noise = value;
noise += value;
count++;
}
return (noise/double(count));
}
/sipmScan/examples/threshold64.C
0,0 → 1,43
 
{
int map[8][8]={{32,34,53,55,40,42,61,63},
{48,50,37,39,56,58,45,47},
{33,35,52,54,41,43,60,62},
{49,51,36,38,57,59,44,46},
{17,19,4,6,25,27,12,14},
{1,3,20,22,9,11,28,30},
{16,18,5,7,24,26,13,15},
{0,2,21,23,8,10,29,31}
};
TCanvas* canvas = new TCanvas("canvas","",1200,1200);
canvas->cd();
//canvas->Divide(4,4);
TVirtualPad *main = new TPad("main","main",0,0,1,1,10,1);
main->Draw();
main->Divide(4,4);
TH1F* h_thrProj[16];
for(int i=0; i<4; i++) {
for(int j=0; j<4; j++) {
int channel = map[i][j];
printf("%d ", channel);
TH1D* histogram = h_threshold->ProjectionY("",channel+1,channel+1);
int nBins = histogram->GetXaxis()->GetNbins();
char name[128];
int current = i*4+j;
sprintf(name,"histo%d",current);
//h_thrProj[i*4+j] = (TH1F*)gROOT->FindObject(name);
if (h_thrProj[i*4+j]) delete h_thrProj[i*4+j];
h_thrProj[current] = new TH1F(name, name, nBins,0,nBins);
for(int k=0; k< nBins; k++){
int signal = histogram->GetBinContent(k);
h_thrProj[current]->Fill(k,signal);
}
int canvasPosition = i+4*(3-j)+1;
main->cd(canvasPosition);
//printf("%d \n", canvasPosition);
h_thrProj[current]->Draw();
}
}
 
}
/sipmScan/examples/sipm.C
0,0 → 1,538
#include "TROOT.h"
#include "TFile.h"
#include "TBenchmark.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TCanvas.h"
#include "TStyle.h"
#include "TPad.h"
#include "TF1.h"
#include "TGraph.h"
#include "TSpectrum.h"
#include "stdio.h"
 
#include "include/RTUtil.h"
 
double getNoise(TH2F*, int, int);
 
int sipm(char filename[256] = "test", char plopt[256]="all", int chXstart=0, int chXend=7, int chYstart=0, int chYend=7, bool debug = false)
{
const int c_nChannels = 64;
const double c_xOffset = 0; // mm
const double c_yOffset = 0;
int map[8][8]={{32,34,53,55,40,42,61,63},
{48,50,37,39,56,58,45,47},
{33,35,52,54,41,43,60,62},
{49,51,36,38,57,59,44,46},
{17,19,4,6,25,27,12,14},
{1,3,20,22,9,11,28,30},
{16,18,5,7,24,26,13,15},
{0,2,21,23,8,10,29,31}
};
char fnameroot[256];
TFile* rootfile;
sprintf(fnameroot, "root/%s.root", filename);
rootfile = (TFile *) gROOT->FindObject(filename);
if(rootfile==NULL) rootfile = new TFile(fnameroot);
if(rootfile==NULL) {
printf("Cannot open root file %s!!!\n",fnameroot);
return(0);
}
 
// set draw style
gStyle->SetOptStat("ne");
gStyle->SetPalette(1, 0);
gStyle->SetPaperSize(TStyle::kA4);
gStyle->SetStatBorderSize(1);
gStyle->SetFrameBorderMode(0);
gStyle->SetFrameFillColor(0);
gStyle->SetCanvasBorderMode(0);
gStyle->SetPadBorderMode(0);
gStyle->SetPadColor(0);
gStyle->SetCanvasColor(0);
gStyle->SetStatColor(0);
gStyle->SetOptFit(11);
gStyle->SetOptStat();
gStyle->SetPadRightMargin(0.15);
gStyle->SetPadLeftMargin(0.12);
//gStyle->SetTitleYOffset(1.4);
 
if( strstr(plopt, "all") != NULL ) {
TCanvas *canvas2 = new TCanvas("canvas2","Hits x;;",2000,2000);
TCanvas *canvas3 = new TCanvas("canvas3","Hits y;;",2000,2000);
canvas2->Divide(8,8);
canvas3->Divide(8,8);
TH1F* h_hitsx;
TH1F* h_hitsy;
for(int i=0; i<c_nChannels; i++) {
canvas2->cd(i+1);
char hname[128];
sprintf(hname, "hnhitsx%d", i);
h_hitsx = (TH1F*)rootfile->Get(hname);
h_hitsx->Draw();
canvas3->cd(i+1);
sprintf(hname, "hnhitsy%d", i);
h_hitsy = (TH1F*)rootfile->Get(hname);
h_hitsy->Draw();
}
}
if( strstr(plopt, "x") != NULL ) {
TCanvas *canvas10 = new TCanvas("canvas10","Ch x;;",500,500);
TH1F* h_hitsx;
canvas10->cd();
char hname[128];
sprintf(hname, "hnhitsx%d", chXstart);
h_hitsx = (TH1F*)rootfile->Get(hname);
h_hitsx->Draw();
}
if( strstr(plopt, "y") != NULL ) {
TCanvas *canvas11 = new TCanvas("canvas11","Ch x;;",500,500);
TH1F* h_hitsy;
canvas11->cd();
char hname[128];
sprintf(hname, "hnhitsy%d", chXstart);
h_hitsy = (TH1F*)rootfile->Get(hname);
h_hitsy->Draw();
}
if( strstr(plopt, "share") != NULL ) {
/*TCanvas *canvas4 = new TCanvas("canvas1","canvas1",1000,1000);
int nChannels = chYend-chYstart+1;
int ncols = nChannels/2;
printf("nch %d nch\\2 %d\n", nChannels, ncols);
canvas4->Divide(2,ncols);
TH1F* h_hitsy;
for(int i=chYstart; i<=chYend; i++){
canvas4->cd(i-chYstart+1);
char hname[128];
int chPosition = map[0][i];
sprintf(hname, "hnhitsy%d", chPosition);
h_hitsy = (TH1F*)rootfile->Get(hname);
h_hitsy->Draw();
}*/
TCanvas *canvas4 = new TCanvas("canvas4","canvas4",500,500);
canvas4->cd();
for(int i=chXstart; i<=chXend; i++) {
TH1F* h_hitsx;
char hname[128];
int chPosition = map[i][chYstart];
sprintf(hname, "hnhitsx%d", chPosition);
h_hitsx = (TH1F*)rootfile->Get(hname);
h_hitsx->SetTitle("Scan X;x [mm]; Entries");
h_hitsx->GetYaxis()->SetTitleOffset(1.3);
h_hitsx->SetStats(0);
if (i == chXstart)
h_hitsx->Draw();
else {
h_hitsx->SetLineColor(i+1);
h_hitsx->Draw("same");
}
}
//sprintf(fullname, "ps/%s_Yshare.eps", filename);
//canvas4->SaveAs(fullname);
TCanvas *canvas5 = new TCanvas("canvas5","canvas5",500,500);
canvas5->cd();
for(int i=chYstart; i<=chYend; i++) {
TH1F* h_hitsy;
char hname[128];
int chPosition = map[chXstart][i];
sprintf(hname, "hnhitsy%d", chPosition);
h_hitsy = (TH1F*)rootfile->Get(hname);
h_hitsy->SetTitle("Scan Y;y [mm]; Entries");
h_hitsy->GetYaxis()->SetTitleOffset(1.3);
h_hitsy->SetStats(0);
if (i == chYstart)
h_hitsy->Draw();
else {
h_hitsy->SetLineColor(i+1);
h_hitsy->Draw("same");
}
}
//sprintf(fullname, "ps/%s_Yshare.eps", filename);
//canvas5->SaveAs(fullname);
}
/** Draw one channel in x and y
* Without and with the noise subtraction
*/
if (strstr(plopt, "noise") != NULL) {
TCanvas *canvas51 = new TCanvas("canvas51","canvas51",500,500);
canvas51->cd();
// Get the filled histogram from scan
TH1F* h_x;
char name[32];
sprintf(name, "hnhitsx%d", chXstart);
h_x = (TH1F*)rootfile->Get(name);
h_x->DrawCopy();
TH1F* h_y;
sprintf(name, "hnhitsy%d", chXstart);
h_y = (TH1F*)rootfile->Get(name);
// Create and fill corrected histogram
Int_t binsX = h_x->GetXaxis()->GetNbins();
if(debug) printf("nBins: %d\n", binsX);
//Double_t xLowUser = h_x->GetXaxis()->GetBinLowEdge(1);
//Double_t xUpUser = h_x->GetXaxis()->GetBinUpEdge(binsX);
Double_t xLowUser = h_x->GetXaxis()->GetXmin();
Double_t xUpUser = h_x->GetXaxis()->GetXmax();
TH1F* h_xCorrected = new TH1F("h_xCorrected",";;",binsX, xLowUser, xUpUser);
if (debug) printf("xLow %f xUp %f\n",xLowUser,xUpUser);
double noise = h_x->GetBinContent(2);
for (int jk = 0; jk < binsX; jk++) {
double signal = h_x->GetBinContent(jk+1);
signal -= noise;
double x = h_x->GetXaxis()->GetBinCenter(jk+1);
h_xCorrected->Fill(x, signal);
}
h_xCorrected->SetLineColor(2);
h_xCorrected->DrawCopy("same");
//Create and fill corrected histogram
Int_t binsY = h_y->GetXaxis()->GetNbins();
Double_t yLowUser = h_y->GetXaxis()->GetXmin();
Double_t yUpUser = h_y->GetXaxis()->GetXmax();
TH1F* h_yCorrected = new TH1F("h_yCorrected",";;",binsY, yLowUser, yUpUser);
noise = h_y->GetBinContent(2);
for (int jk = 0; jk < binsY; jk++) {
double signal = h_y->GetBinContent(jk+1);
signal -= noise;
double x = h_y->GetXaxis()->GetBinCenter(jk+1);
h_yCorrected->Fill(x, signal);
}
TCanvas *canvas52 = new TCanvas("canvas52","canvas52",500,500);
canvas52->cd();
h_y->DrawCopy();
h_yCorrected->SetLineColor(2);
h_yCorrected->DrawCopy("same");
}
/** Draws the signal from 8 channels in x-row
* for one specific y bin, so the background and cross-talk
* can be estimated.
* Draws also a 2d scan of these channels.
*/
if (strstr(plopt, "line") != NULL) {
TCanvas* canvas6 = new TCanvas("canvas6","canvas6",500,500);
canvas6->cd(0);
gStyle->SetOptStat(0);
TH2F* h0 = (TH2F*) rootfile->Get("h2d0");
Int_t binsX = h0->GetXaxis()->GetNbins();
Int_t minX = h0->GetXaxis()->GetFirst();
Int_t maxX = h0->GetXaxis()->GetLast()+1;
Int_t binsY = h0->GetYaxis()->GetNbins();
Int_t minY = h0->GetYaxis()->GetFirst();
Int_t maxY = h0->GetYaxis()->GetLast()+1;
Double_t xLowUser = h0->GetXaxis()->GetBinLowEdge(minX);
Double_t xUpUser = h0->GetXaxis()->GetBinUpEdge(maxX);
Double_t yLowUser = h0->GetYaxis()->GetBinLowEdge(minY);
Double_t yUpUser = h0->GetYaxis()->GetBinUpEdge(maxY);
if (debug) printf("xLow %f xUp %f\n",xLowUser,xUpUser);
//! 1-dimension position in x vs. hits
TH2F* h[8];
TH1F* h_line[8];
for(int j=0; j<8; j++) {
h_line[j] = new TH1F("h_line", "h_line", binsX, xLowUser, xUpUser);
}
for(int j=chXstart; j<=chXend; j++) {
int chPosition = map[j][chYstart];
char hname[128];
sprintf(hname, "h2d%d", chPosition);
int histogram = j;
h[histogram] = (TH2F *) rootfile->Get(hname);
int noise = getNoise(h[histogram], 1, 160);
for(int k=minX; k<=maxX; k++) {
int l=chYstart*20+12;
//for(int l=12; l<=16; l++) {
double signal = h[histogram]->GetBinContent(k,l);
//signal -= noise;
//signal /= 5*10000.0;
double eta = -log(1 - signal);
double x = xLowUser + k*(xUpUser-xLowUser)/double(binsX);
//double y = l*(yUpUser-yLowUser)/double(binsY);
h_line[j]->Fill(x-c_xOffset, signal);
//}
}
if (j == chXstart) {
h_line[j]->SetTitle("SiPM#2 w/o noise subtraction;x[mm];Hits");
//h_line[j]->GetYaxis()->SetRangeUser(-0.05, 0.3);
//h_line[j]->GetYaxis()->SetRangeUser(-50, 2500);
h_line[j]->Draw("");
}
else {
h_line[j]->SetLineColor(j+1);
h_line[j]->Draw("same");
}
}
//! 2d scan
TCanvas* canvas61 = new TCanvas("canvas61","canvas61",8*200,300);
canvas61->cd();
TH2F* h_corrected = new TH2F("h_corrected","h_corrected",binsX,xLowUser,xUpUser, binsY,yLowUser,yUpUser);
for(int i=chXstart; i<=chXend; i++) {
//int canvasPosition = nX*(i-chYstart)+(j-chXstart)+1;
//int canvasPosition = nX*(chYend-i)+chXstart+1;
//if (debug) printf("canvas %d\n",canvasPosition);
int chPosition = map[i][chYstart];
char hname[128];
sprintf(hname, "h2d%d", chPosition);
int histogram = i;
h[histogram] = (TH2F *) rootfile->Get(hname);
int noise = getNoise(h[histogram], 1, 100);
for(int k=minX; k<=maxX; k++) {
for(int l=minY; l<=maxY; l++) {
int signal = h[histogram]->GetBinContent(k,l); // detected
//p /= 10000.;
//double p0 = 1.0 - p; // events with zero photons
//double eta = (-log(p0) * p0 ) / (1-p0-0.00001);
//double eta = -log(p0);
//printf("p=%f p0=%f log(p0)=%f eta=%f\n",p,p0,log(p0),eta);
//double signal = ((p - noise) > 0.1) ? (p-noise) : 0.1;
double p = signal - noise;
p /= 10000.0;
double eta = -log(1 - p);
//double x = k*(xUpUser-xLowUser)/double(binsX);
//double y = l*(yUpUser-yLowUser)/double(binsY);
double x = h[histogram]->GetXaxis()->GetBinCenter(k);
double y = h[histogram]->GetYaxis()->GetBinCenter(l);
h_corrected->Fill(x-c_xOffset, y-c_yOffset, eta);
}
}
}
h_corrected->SetTitle("SiPM#2 n_pe = - ln(P0);x[mm];y[mm]");
h_corrected->GetZaxis()->SetRangeUser(-0.05,0.30);
h_corrected->Draw("colz");
}
/** Draws the sum of the channels
* Each channel is a 2d plot
* Intended for the study of 1 channel
*
* Input: scanned channel coordinates (x,y)
*/
if (strstr(plopt, "2d") != NULL) {
int nX = 3;
int nY = 3;
TCanvas* canvas7 = new TCanvas("canvas7","canvas7", nX*400,nY*400);
printf("nx %d ny %d\n",nX,nY);
canvas7->Divide(nX,nY);
for(int i=chXend-1; i<=chXend+1; i++) {
for(int j=chXstart-1; j<=chXstart+1; j++) {
//int canvasPosition = nX*(i-chYstart)+(j-chXstart)+1;
int canvasPosition = nX*(chXend+1-i) + (j-chXstart+1) +1;
if (debug) printf("canvas %d\n",canvasPosition);
canvas7->cd(canvasPosition);
char hname[128];
int chPosition = map[j][i];
sprintf(hname, "h2d%d", chPosition);
TH2F* h_2d = (TH2F*)rootfile->Get(hname);
h_2d->Draw("colz");
} //x
}
// Number of photoelectrons - Poissonian correction
TCanvas* canvas8 = new TCanvas("canvas8","canvas8", 1000,1000);
canvas8->cd();
gStyle->SetOptStat(0);
char hname[128];
int chPosition = map[chXstart][chXend];
sprintf(hname, "h2d%d", chPosition);
TH2F* h_2d = (TH2F*)rootfile->Get(hname);
 
Int_t binsX = h_2d->GetXaxis()->GetNbins();
Int_t minX = h_2d->GetXaxis()->GetFirst();
Int_t maxX = h_2d->GetXaxis()->GetLast()+1;
Int_t binsY = h_2d->GetYaxis()->GetNbins();
Int_t minY = h_2d->GetYaxis()->GetFirst();
Int_t maxY = h_2d->GetYaxis()->GetLast()+1;
Double_t xLowUser = h_2d->GetXaxis()->GetXmin();
Double_t xUpUser = h_2d->GetXaxis()->GetXmax();
Double_t yLowUser = h_2d->GetYaxis()->GetXmin();
Double_t yUpUser = h_2d->GetYaxis()->GetXmax();
if (debug) printf("xLow %f xUp %f\n",xLowUser,xUpUser);
TH2F* h_corrected = new TH2F("h_corrected","h_corrected",binsX, xLowUser, xUpUser, binsY, yLowUser, yUpUser);
//TH2F* h_corrected = new TH2F("h_corrected","h_corrected",binsX, minX, maxX, binsY, minY, maxY);
double noise = getNoise(h_2d, 1, 70);
if(debug) printf("Noise = %f\n", noise);
for(int k=minX; k<=maxX; k++) {
for(int j=minY; j<=maxY; j++) {
double signal = h_2d->GetBinContent(k,j); // detected
//double p = ((signal - noise) > 1) ? (signal-noise) : 1;
double p = signal - noise;
p /= 10000.;
double p0 = 1.0 - p; // events with zero photons
//double eta = (-log(p0) * p0 ) / (1-p0-0.00001);
double eta = -log(p0); // constant of the poissonian statistics
if (debug) printf("p=%f p0=%f log(p0)=%f eta=%f\n",p,p0,log(p0),eta);
//double x = xLowUser + k*(xUpUser - xLowUser) / double(binsX);
double x = h_2d->GetXaxis()->GetBinCenter(k);
//double y = yLowUser + j*(yUpUser-yLowUser)/double(binsY);
double y = h_2d->GetYaxis()->GetBinCenter(j);
if (debug) printf("x=%f y=%f\n", x, y);
h_corrected->Fill(x + c_xOffset, y + c_yOffset, eta);
}
}
//h_corrected->SetTitle("n_pe = - ln P(0);x[mm];y[mm]");
h_corrected->SetTitle(";x[mm];y[mm]");
h_corrected->GetZaxis()->SetRangeUser(-0.01,0.5);
h_corrected->SetContour(50);
h_corrected->Draw("colz");
// collection efficiency
int nPoints =0;
double efficiency=0;
for (int i=18; i<=58; i++) {
for (int j=19; j<=59; j++) {
double signal = h_corrected->GetBinContent(i,j);
if(debug) printf("signal %f\n",signal);
efficiency += signal;
nPoints++;
}
}
printf("Signal sum = %f\n # of points = %d\n",efficiency,nPoints);
}
/** Draws the sum of channel signals
* Each channel is a 2d ('h2d') histogram
* Suitable for 8x8 chs scan
*/
if( strstr(plopt, "sum") != NULL ) {
int nX = chXend - chXstart + 1;
int nY = chYend - chYstart + 1;
TCanvas* canvas12 = new TCanvas("canvas12","c2",8*200, 8*200);
canvas12->cd();
gStyle->SetOptStat(0);
 
// final histogram parameters
TH2F* h0 = (TH2F*) rootfile->Get("h2d0");
Int_t binsX = h0->GetXaxis()->GetNbins();
Int_t minX = h0->GetXaxis()->GetFirst();
Int_t maxX = h0->GetXaxis()->GetLast()+1;
Int_t binsY = h0->GetYaxis()->GetNbins();
Int_t minY = h0->GetYaxis()->GetFirst();
Int_t maxY = h0->GetYaxis()->GetLast()+1;
Double_t xLowUser = h0->GetXaxis()->GetBinLowEdge(minX);
Double_t xUpUser = h0->GetXaxis()->GetBinUpEdge(maxX);
Double_t yLowUser = h0->GetYaxis()->GetBinLowEdge(minY);
Double_t yUpUser = h0->GetYaxis()->GetBinUpEdge(maxY);
if (debug) printf("xLow %f xUp %f\n",xLowUser,xUpUser);
TH2F* h_corrected = new TH2F("h_corrected","h_corrected",binsX,xLowUser,xUpUser, binsY,yLowUser,yUpUser);
TH2F* h[nX*nY];
// 2d histogram noise subtraction and poisson scaling
for(int i=chYstart; i<=chYend; i++) {
for(int j=chXstart; j<=chXend; j++) {
int chPosition = map[j][i];
char hname[128];
sprintf(hname, "h2d%d", chPosition);
int histogram = nX*(i-chYstart)+(j-chXstart);
h[histogram] = (TH2F *) rootfile->Get(hname);
int noise = getNoise(h[histogram], 1, 170);
if (debug) printf("noise: %d\n",noise);
for(int k=minX; k<=maxX; k++) {
for(int l=minY; l<=maxY; l++) {
int signal = h[histogram]->GetBinContent(k,l); // detected
//double p = ((signal - noise) > 0.1) ? (signal-noise) : 0.1;
double p = signal - noise;
p /= 10000.;
double p0 = 1.0 - p; // events with zero photons
//double eta = (-log(p0) * p0 ) / (1-p0-0.00001);
double eta = -log(p0);
//printf("p=%f p0=%f log(p0)=%f eta=%f\n",p,p0,log(p0),eta);
double x = k*(xUpUser-xLowUser)/double(binsX);
double y = l*(yUpUser-yLowUser)/double(binsY);
h_corrected->Fill(x-c_xOffset, y-c_yOffset, eta);
}
}
}
}
h_corrected->SetTitle("SiPM#2 n_p.e.;x[mm];y[mm]");
h_corrected->GetZaxis()->SetRangeUser(-0.05,.30);
h_corrected->Draw("colz");
TCanvas* canvas13 = new TCanvas("canvas13","canvas13",600,300);
canvas13->Divide(2);
canvas13->cd(1);
h[16]->Draw("colz");
canvas13->cd(2);
h[8]->Draw("colz");
}
/** Draws the beam profile and fits it with error function
* on some background function
*/
if (strstr(plopt, "beam") != NULL) {
TCanvas* canvas9 = new TCanvas("canvas9","canvas9", 500,500);
canvas9->cd();
char hname[128];
sprintf(hname, "hnhitsx%d", 36);
TH1F* h_laser = (TH1F*)rootfile->Get(hname);
h_laser->Draw();
h_laser->SetStats(1);
TF1* err = new TF1("err","[0]+[1]*TMath::Erf((x-[2])/[3])",17.02,17.30);
err->SetParameter(0,2500);
err->SetParameter(1, h_laser->GetMaximum());
err->SetParameter(2, h_laser->GetBinCenter(h_laser->GetMaximumBin()));
err->SetParameter(3, 0.001);
h_laser->Fit(err,"qr");
h_laser->Fit(err,"lr");
double sigma = err->GetParameter(3);
printf("sigma = %2.0f um, FWHM = %2.0f um\n", sigma*1000, 2.35*sigma*1000);
}
if (strstr(plopt, "map") != NULL) {
TCanvas* canvas1 = new TCanvas("canvas1","canvas1",1000,1000);
canvas1->cd();
TH2I* h_map = new TH2I("h_map","h_map",8,-0.5,7.5,8,-0.5,7.5);
for (int i=7; i>=0; i--) {
//for (int j=7; j>=0; j--) printf("(%d, %d) ", j,i);
for (int j=7; j>=0; j--) {
printf("%2d (%d %d)\n", map[j][i], j*10080+5040, i*10080+5040);
h_map->Fill(j,i, map[j][i]);
}
printf("\n");
}
gStyle->SetOptStat(0);
h_map->Draw("text");
}
return(0);
}
 
/** Function calculates the noise from one channel
* it runs through the bins along x and returns the average value
*/
double getNoise(TH2F* histogram, int yStart, int yEnd)
{
double noise=0;
int count=0;
for(int j=yStart; j<yEnd; j++) {
double value = histogram->GetBinContent(j,2);
//if (noise < value) noise = value;
noise += value;
count++;
}
return (noise/double(count));
}
/sipmScan/examples/tdc.C
0,0 → 1,109
#include "TROOT.h"
#include "TFile.h"
#include "TBenchmark.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TCanvas.h"
#include "TStyle.h"
#include "TPad.h"
#include "TF1.h"
#include "TGraph.h"
#include "TSpectrum.h"
#include "stdio.h"
 
#include "include/RTUtil.h"
 
int tdc(char filename[256] = "test", int chX=0, int chY=0, double rangeLeft=-16, double rangeRight=16, bool debug = false)
{
//const int c_nChannels = 64;
//const double c_xOffset = 2.2; // mm
//const double c_yOffset = 2.3;
int map[8][8]={{32,34,53,55,40,42,61,63},
{48,50,37,39,56,58,45,47},
{33,35,52,54,41,43,60,62},
{49,51,36,38,57,59,44,46},
{17,19,4,6,25,27,12,14},
{1,3,20,22,9,11,28,30},
{16,18,5,7,24,26,13,15},
{0,2,21,23,8,10,29,31}
};
char fnameroot[256];
TFile* rootfile;
sprintf(fnameroot, "root/%s.root", filename);
rootfile = (TFile *) gROOT->FindObject(filename);
if(rootfile==NULL) rootfile = new TFile(fnameroot);
if(rootfile==NULL) {
printf("Cannot open root file %s!!!\n",fnameroot);
return(0);
}
 
// set draw style
gStyle->SetPalette(1, 0);
gStyle->SetPaperSize(TStyle::kA4);
gStyle->SetStatBorderSize(1);
gStyle->SetFrameBorderMode(0);
gStyle->SetFrameFillColor(0);
gStyle->SetCanvasBorderMode(0);
gStyle->SetPadBorderMode(0);
gStyle->SetPadColor(0);
gStyle->SetCanvasColor(0);
gStyle->SetStatColor(0);
gStyle->SetOptFit(11);
gStyle->SetOptStat(0);
gStyle->SetPadRightMargin(0.15);
gStyle->SetPadLeftMargin(0.12);
//gStyle->SetTitleYOffset(1.4);
TCanvas* canvas1 = new TCanvas("canvas1","canvas1",1000,1000);
TH2F* htdc = (TH2F*) rootfile->Get("htdc");
canvas1->cd();
htdc->Draw("colz");
TH2F* h_correctedTDC = (TH2F*) rootfile->Get("h_correctedTDC");
TCanvas* canvas2 = new TCanvas("canvas2","canvas2",800,800);
canvas2->cd();
h_correctedTDC->Draw("colz");
/*
canvas2->cd(2);
int binY = map[chX][chY];
TH1D* channelY = h_correctedTDC->ProjectionX("",binY+1,binY+1);
//channelY->SetStats(0);
char title[256];
sprintf(title,"Channel %d;t [ns];Events", binY);
channelY->SetTitle(title);
channelY->GetYaxis()->SetTitleOffset(1.7);
//TAxis* xAxis = h_correctedTDC->GetXaxis();
//int range = xAxis->GetBinUpEdge(xAxis->GetLast()+1);
//channelY->GetXaxis()->SetRangeUser(-range, range);
channelY->Draw();
 
TF1* f_gaus1 = new TF1("f_gaus1","[0] + gaus(1)", rangeLeft,rangeRight);
TF1* f_gaus2 = new TF1("f_gaus2","[0] + gaus(1) + gaus(4)",-8,8);
f_gaus1->SetParNames("Linear","Norm","#mu","#sigma");
f_gaus2->SetParNames("Linear","Norm1","Mean1","Sigma1","Norm2","Mean2","Sigma2");
Int_t n = channelY->GetMaximum();
Float_t mean = channelY->GetBinCenter(channelY->GetMaximumBin());
f_gaus1->SetParameters(channelY->GetMinimum(), n, mean, 2.0);
channelY->Fit(f_gaus1,"0");
channelY->Fit(f_gaus1,"r");
*/
//Fit in range +- 3 sigma
//channelY->Fit(f_gaus1,"r","", f_gaus1->GetParameter(2)-3*f_gaus1->GetParameter(3),
// f_gaus1->GetParameter(2)+3*f_gaus1->GetParameter(3));
 
//f_gaus2->SetParameters(-1, n, mean-0.1, 0.01,
// n, mean+0.1, 0.1);
//channelY->Fit(f_gaus2,"0");
//channelY->Fit(f_gaus2,"rl");
// set range to +- 3 sigma
//channelY->Fit(f_gaus,"lr","",f_gaus->GetParameter(2)-3*f_gaus->GetParameter(3),
// f_gaus->GetParameter(2)+3*f_gaus->GetParameter(3));
return (0);
}