Subversion Repositories f9daq

Rev

Rev 205 | Rev 208 | Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
101 f9daq 1
#include "TH1.h"
2
#include "TH2.h"
3
#include "TGraph.h"
4
#include "TCanvas.h"
5
#include "TStyle.h"
6
#include "TLegend.h"
204 f9daq 7
#include "TFile.h"
101 f9daq 8
//#include "TROOT.h"
9
 
10
#include <stdio.h>
180 f9daq 11
#include <time.h>
101 f9daq 12
 
13
#define DATA_PATH  "./data/"
14
#define CFRAC 0.4
15
#define DWIDTH 1024
180 f9daq 16
#define NCH 4
17
#include "drs.h"
18
 
19
ClassImp (drs)
20
 
21
drs::drs( const char *FileName, int nch, DrsChannel* d, int trgch, int nev, int updfrq)
101 f9daq 22
{
180 f9daq 23
        int debug=0;
24
        int i,ich,iev;
25
        int stat;
26
        char histname[128];
101 f9daq 27
//  char FileName[300]="./data/n2.dat";
180 f9daq 28
        FILE *fin;
101 f9daq 29
 
180 f9daq 30
        char recid[5]="SAMO";
31
        unsigned short *srecid = (unsigned short *) &recid[0];
175 f9daq 32
 
180 f9daq 33
        struct evrec {
34
                unsigned long evn;
35
                unsigned short tYear,tMonth,tDay,tHour,tMin,tSec,tmSec,rc;
36
        } evrec;
101 f9daq 37
 
180 f9daq 38
        float dtdata[NCH][DWIDTH];
39
        float tcdata[NCH][DWIDTH];
40
        float tccalib[NCH];
41
        float scaler[NCH];
42
        float *t;
43
        float *y;
44
        float *ay;
101 f9daq 45
 
180 f9daq 46
        short lwfm[NCH];
47
        int nwfm[NCH];
48
        int board;
49
        int version;
175 f9daq 50
 
180 f9daq 51
 
52
        unsigned short ichdat[NCH][DWIDTH];
53
        float chdat[NCH][DWIDTH];
54
        float achdat[NCH][DWIDTH];
55
        float chmax[NCH][2];
56
        float chmin[NCH][2];
57
        float letime[NCH],cftime[NCH],qdc[NCH];
58
 
59
 
101 f9daq 60
//  gROOT->Reset();
61
 
180 f9daq 62
        gStyle->SetOptStat(1111111);
101 f9daq 63
//  gStyle->SetPalette(52);
64
 
180 f9daq 65
        TH1F *h1_chmax[NCH];
66
        TH1F *h1_chmaxt[NCH];
67
        TH1F *h1_chmin[NCH];
68
        TH1F *h1_chmint[NCH];
69
        TH1F *h1_letime[NCH];
70
        TH1F *h1_cftime[NCH];
71
        TH1F *h1_cftdif[NCH];
72
        TH1F *h1_qdc[NCH];
73
        TH2F *h2_dpo[NCH];
74
        TH2F *h2_tdcadc[NCH];
75
        TH2F *h2_ctdcadc[NCH];
76
        TH2F *h2_ctalk;
77
        TGraph *gr_wfm[NCH];
204 f9daq 78
 
79
		TFile *froot = new  TFile(TString(FileName) + ".root","RECREATE");
80
 
180 f9daq 81
        TCanvas *c_wfm= new TCanvas("c_wfm","WFM",50,50,700,700);
82
        TLegend *leg = new TLegend(0.8, 0.8, 1, 1);
83
        c_wfm->Divide(1,1);
84
        c_wfm->cd(1);
101 f9daq 85
 
180 f9daq 86
        for (i=0; i<nch; i++) {
87
                sprintf(histname,"chmax_%02d",i);
204 f9daq 88
                h1_chmax[i]=new TH1F(histname,histname,1100,-0.88,0.88);
180 f9daq 89
                sprintf(histname,"chmaxt_%02d",i);
90
                h1_chmaxt[i]=new TH1F(histname,histname,DWIDTH,-0.1,204.7);
91
                sprintf(histname,"chmin_%02d",i);
204 f9daq 92
                h1_chmin[i]=new TH1F(histname,histname,1100,-0.55,0.55);
180 f9daq 93
                sprintf(histname,"chmint_%02d",i);
204 f9daq 94
                h1_chmint[i]=new TH1F(histname,histname,DWIDTH,-0.1,504.7);
180 f9daq 95
                sprintf(histname,"letime_%02d",i);
96
                h1_letime[i]=new TH1F(histname,histname,500,d[i].twin[0],d[i].twin[1]);
97
                sprintf(histname,"cftime_%02d",i);
98
                h1_cftime[i]=new TH1F(histname,histname,500,d[i].twin[0],d[i].twin[1]);
99
                sprintf(histname,"cftdif_%02d",i);
182 f9daq 100
                h1_cftdif[i]=new TH1F(histname,histname,2000,d[i].twin[0]-d[trgch].twin[0],d[i].twin[1]-d[trgch].twin[0]);
180 f9daq 101
                sprintf(histname,"qdc_%02d",i);
182 f9daq 102
                h1_qdc[i]=new TH1F(histname,histname,500,-0.5,4.5);
180 f9daq 103
                sprintf(histname,"dpo_%02d",i);
204 f9daq 104
                h2_dpo[i]=new TH2F(histname,histname,DWIDTH,-0.1,200.7,3000,-0.9,0.9);
180 f9daq 105
                sprintf(histname,"tdcadc_%02d",i);
182 f9daq 106
                h2_tdcadc[i]=new TH2F(histname,histname,250,-0.1,0.5,400,d[i].twin[0]-d[trgch].twin[0],d[i].twin[1]-d[trgch].twin[0]);
180 f9daq 107
                sprintf(histname,"ctdcadc_%02d",i);
182 f9daq 108
                h2_ctdcadc[i]=new TH2F(histname,histname,250,-0.1,0.5,400,d[i].twin[0]-d[trgch].twin[0],d[i].twin[1]-d[trgch].twin[0]);
180 f9daq 109
                gr_wfm[i]=new TGraph(DWIDTH);
204 f9daq 110
                gr_wfm[i]->SetMinimum(-0.9);
111
                gr_wfm[i]->SetMaximum(0.9);
180 f9daq 112
                gr_wfm[i]->SetLineColor(i+1);
113
                gr_wfm[i]->SetMarkerColor(i+1);
114
                sprintf(histname,"ch %d",i);
115
                leg->AddEntry(gr_wfm[i], histname, "l");
116
                if (i) gr_wfm[i]->Draw("LP"); else gr_wfm[i]->Draw("ALP");
117
                gr_wfm[i]->GetHistogram()->SetTitle("WaveForms");
118
        }
101 f9daq 119
 
180 f9daq 120
        leg->Draw();
121
        h2_ctalk=new TH2F("ctalk","ctalk",250,0.,0.5,250,0.,0.5);
204 f9daq 122
 
180 f9daq 123
        fin=fopen(FileName,"rb");
124
        if (fin==NULL) {
125
                printf("Error opening file %s\n",FileName);
126
                return;
127
        }
128
        printf("%s, %lu\n",recid,sizeof(evrec));
101 f9daq 129
 
130
 
180 f9daq 131
        for (ich=0; ich<nch; ich++) {
132
                nwfm[ich]=0;
133
                for (i=0; i<DWIDTH; i++) {
134
                        achdat[ich][i]=0;
135
                }
136
        }
137
        unsigned int trgcell = 0;
138
        int event = -1;
139
 
140
        iev=0;
141
        time_t t0,told;
142
        time(&t0);
143
        told=t0-1;
206 f9daq 144
		int posrec[9];
180 f9daq 145
        while(!feof(fin)) {
146
                if (iev>=nev) break;
147
                stat=fread(recid,1,4,fin);
148
                if (debug) {
149
                        if (recid[1]=='#') printf("%c%c%u\n",recid[0],recid[1],srecid[1]);
150
                        else printf("%s\n",recid );
151
                }
152
                switch (recid[0]) {
202 f9daq 153
 
205 f9daq 154
				        case 'P':{ // position record
155
									      int len;
206 f9daq 156
					       				  stat=fread(&len,1,4,fin);
157
										  if (len>0 && len<37) {
158
                                                     stat=fread(posrec,1,len,fin);
159
													 if (recid[3]=='R')
160
												        printf("Position  len=%d x=%d y=%d ix=%d iy=%d\n", len, posrec[0], posrec[1], posrec[3], posrec[4]);
161
												     else
162
                                                        printf("Position record buffer length %d\n", len);
163
										  } else {
164
													printf("Wrong buffer length %d\n", len);
165
										  }
202 f9daq 166
                        break;
205 f9daq 167
								}
180 f9daq 168
                case 'D': // DRS
169
                        version = atoi(&recid[3]);
170
                        break;
171
                case 'T': // TIME
172
                        switch (recid[1]) {
173
                        case 'I':
174
                                event=0;
175
                                if (version ==0) version++;
176
                                break;
177
                        case '#': trgcell = srecid[1]; break;
178
                        }
179
                        break;
180
                case 'B':
181
                        board = srecid[1];
182
                        break;
183
                case 'C':
184
                        sscanf(recid,"C%d",&ich);
185
                        ich=ich-1;
186
                        if (event) {
187
                                if (version>1) stat=fread(&scaler[ich],1,sizeof(float),fin);
188
                                stat=fread(ichdat[ich],1,sizeof(ichdat[ich]),fin);
189
                        } else {
190
 
191
                                stat=fread(dtdata[ich],1,sizeof(float)*DWIDTH,fin);
192
 
193
                        }
194
                        break;
195
                case 'E':
196
                        stat=fread(&evrec.evn,1,sizeof(evrec),fin);
197
                        iev++;
198
                        if (!version)  {
199
                                stat=fread(dtdata[0],1,sizeof(float)*DWIDTH,fin);
200
                        }
201
                        event=1;
202
                        break;
203
                default:
204
                        printf("Unknown header! %s\n", recid);
205
                        break;
206
                }
207
 
208
 
209
                if (recid[0] == 'E') {
210
                        for ( i=0; i<nch; i++) {
211
                                if (!lwfm[i]) continue;
182 f9daq 212
								double pulseheight =  (d[i].edge)? chmax[i][0]: chmin[i][0];
213
 
214
								if ((pulseheight >d[i].vcut[0])&&
215
								    (pulseheight <d[i].vcut[1]) ) {
216
 
217
                                        h1_cftdif[i]->Fill(cftime[i]-cftime[trgch]);
180 f9daq 218
                                }
182 f9daq 219
                                if (!d[i].edge) pulseheight = - pulseheight;
220
                                h2_tdcadc[i] ->Fill(pulseheight,letime[i]-cftime[trgch]);
221
                                h2_ctdcadc[i]->Fill(pulseheight,cftime[i]-cftime[trgch]);
222
 
223
								//if (i==1) printf("ph=%f %f\n", pulseheight, letime[i]-cftime[trgch]);
180 f9daq 224
                        }
182 f9daq 225
						double ph0 =  (d[0].edge)? chmax[0][0]: -chmin[0][0];
226
						double ph1 =  (d[1].edge)? chmax[1][0]: -chmin[1][0];
227
                        h2_ctalk->Fill(ph0,ph1);
180 f9daq 228
                }
229
                if (event<1) continue;
230
                if (recid[0]!='C') continue;
231
                if (ich>=nch) continue;
101 f9daq 232
 
180 f9daq 233
                lwfm[ich]=0;
234
                qdc[ich]=0;
235
                chmax[ich][0]=-1; chmax[ich][1]=-1;
236
                chmin[ich][0]=1; chmin[ich][1]=-1;
237
                letime[ich]=-1.; letime[ich]=-1.;
238
                cftime[ich]=-1.; cftime[ich]=-1.;
101 f9daq 239
 
240
 
180 f9daq 241
                lwfm[ich]=1;
242
                nwfm[ich]+=1;
101 f9daq 243
 
244
 
180 f9daq 245
 
246
                t  = &tcdata[ich][0];
247
                y  = &chdat[ich][0];
248
                ay = &achdat[ich][0];
249
 
250
                float sum=0;
251
                for (i=0; i<DWIDTH; i++) {
252
                        sum+= dtdata[ich][(i+trgcell)%DWIDTH];
253
                        if (version) t[i]=sum;
254
                        else t[i]=dtdata[0][i]; // old format
255
                }
256
                tccalib[ich]=t[(DWIDTH-trgcell)%DWIDTH];
257
 
258
                for (i=0; i<DWIDTH; i++) {
259
                        if (version) t[i]-=(tccalib[ich]-tccalib[0]); // assume ich 0 is connected and always first in the data
260
                        y[i]=-0.5+(float)ichdat[ich][i]/0xFFFF;
192 f9daq 261
                        if (version) y[i]+=evrec.rc/1000.;
180 f9daq 262
                        h2_dpo[ich]->Fill(t[i],y[i]);
263
                        gr_wfm[ich]->SetPoint(i,t[i],y[i]+0.05*(ich-1));
264
                }
265
 
266
 
267
                for (i=0; i<DWIDTH; i++) {
268
                        // accumulated waveforms
269
                        ay[i]+=y[i];
270
 
271
                        // signal minumum
272
                        if (y[i]>chmax[ich][0]) {
273
                                chmax[ich][0]=y[i];
274
                                chmax[ich][1]=t[i];
275
                        }
276
 
277
                        // signal maximum
278
                        if (y[i]<chmin[ich][0]) {
279
                                chmin[ich][0]=y[i];
280
                                chmin[ich][1]=t[i];
281
                        }
282
                        float t0  = t[i];
283
                        float t0p = t[i+1];
284
                        float t0n = t[i-1];
285
                        // charge in the range
286
                        if ((t0>d[ich].adcgate[0])&&(t0<d[ich].adcgate[1])) {
182 f9daq 287
                                float sign =  (d[ich].edge)? 1: -1;
288
								qdc[ich]+=sign*y[i]*(t0p-t0n)/2.;
180 f9daq 289
                        }
290
                }
291
 
292
                for (i=1; i<DWIDTH-1; i++) {
293
                        float t0  = t[i];
294
                        float t0p = t[i+1];
295
                        float t0n = t[i-1];
296
                        // leading edge in the range
297
                        if ((t0>d[ich].twin[0])&&(t0<d[ich].twin[1])) {
298
                                if ( ((!d[ich].edge)&&(y[i]<d[ich].threshold) ||
299
                                      ( d[ich].edge)&&(y[i]>d[ich].threshold ) )
300
                                     &&(letime[ich]<0.)) {
301
                                        letime[ich]=t0n+(t0-t0n)*(d[ich].threshold-y[i-1])/(y[i]-y[i-1]);
302
                                }
303
                                // constant fraction time in the range
204 f9daq 304
								if (cftime[ich]<0.){
305
								  // negative signals
306
                                  if ((!d[ich].edge) && (chmin[ich][0]<d[ich].threshold)&&(y[i]<chmin[ich][0]*d[ich].cfrac )){
307
									  cftime[ich]=t0n+(t0-t0n)*(chmin[ich][0]*d[ich].cfrac-y[i-1])/(y[i]-y[i-1]);
308
								  }
309
								  // positive signals
310
                                  if (( d[ich].edge) && (chmax[ich][0]>d[ich].threshold)&&(y[i]>chmax[ich][0]*d[ich].cfrac)){
311
									  cftime[ich]=t0n+(t0-t0n)*(chmax[ich][0]*d[ich].cfrac-y[i-1])/(y[i]-y[i-1]);
312
								  }
313
								}
180 f9daq 314
                        }
315
                }
316
 
317
                if ((chmax[ich][1]>d[ich].twin[0])&&(chmax[ich][1]<d[ich].twin[1])) {
318
                        h1_chmax[ich]->Fill(chmax[ich][0]);
319
                        if (chmax[ich][0]>d[ich].threshold) h1_chmaxt[ich]->Fill(chmax[ich][1]);
320
                }
321
                if ((chmin[ich][1]>d[ich].twin[0])&&(chmin[ich][1]<d[ich].twin[1])) {
322
                        h1_chmin[ich]->Fill(chmin[ich][0]);
323
                        if (chmin[ich][0]<d[ich].threshold) h1_chmint[ich]->Fill(chmin[ich][1]);
324
                }
325
                h1_qdc[ich]->Fill(qdc[ich]);
326
                h1_letime[ich]->Fill(letime[ich]);
327
                h1_cftime[ich]->Fill(cftime[ich]);
328
                time(&t0);
329
 
330
                if ((t0!=told) || ((updfrq>0)&&(!(iev%updfrq) ) ) ) {
331
                        c_wfm->Update();
332
                        for (i=0; i<DWIDTH; i++) {
333
                                ay[i]/=nwfm[ich];
334
                                gr_wfm[ich]->SetPoint(i,t[i],ay[i]+0.00*(ich-1));
335
                        }
336
                        printf("processing event %d\n",iev);
337
                        told=t0;
338
                }
175 f9daq 339
        }
101 f9daq 340
 
180 f9daq 341
        printf("number of events: %d\n",iev);
342
        fclose(fin);
101 f9daq 343
 
180 f9daq 344
        TCanvas* c_dpo= new TCanvas("c_dpo","DPO",750,50,700,700);
345
        c_dpo->Divide(2,2);
204 f9daq 346
        for (i=0; i<nch; i++) {c_dpo->cd(i+1); h2_dpo[i]->DrawCopy("colz");}
101 f9daq 347
 
180 f9daq 348
        TCanvas* c_cftdif= new TCanvas("c_cftdif","DPO",750,50,700,700);
349
        c_cftdif->Divide(2,2);
204 f9daq 350
        for (i=0; i<nch; i++) {c_cftdif->cd(i+1); h1_cftdif[i]->DrawCopy();}
351
		froot->Write();
352
		froot->Close();
180 f9daq 353
        return;
354
}
101 f9daq 355
 
180 f9daq 356
 
357
void drsana(int nev=1000, int updfrq=20, char *FileName="./data/x1.dat", int trgch=0){
358
 
359
//  float threshold[nch][2]={{0.,-0.02},{0.,0.},{0.,0.},{0.,-0.25}};
360
//  float twin[nch][2]={{120.,130.},{100.,150.},{100.,150.},{70.,90.}};
361
//  float adcgate[nch][2]={{12nch.,136.},{120.,140.},{120.,140.},{70.,90.}};
362
//  float threshold[4][2]={{0.,-0.05},{0.,-0.05},{0.,-0.25},{0.,-0.25}};
363
//  float twin[4][2]={{80.,100.},{80.,100.},{90.,110.},{50.,70.}};
364
//  float adcgate[4][2]={{80.,110.},{80.,110.},{90.,110.},{50.,70.}};
182 f9daq 365
        double threshold[4][2]={{0.,-0.05},{0.,-0.05},{0.,-0.25},{0.,-0.25}};
180 f9daq 366
        float twin[4][2]={{100.,120.},{100.,120.},{110.,130.},{60.,80.}};
367
        float adcgate[4][2]={{105.,125.},{105.,125.},{110.,130.},{60.,80.}};
182 f9daq 368
        double vcut[4][2]={{-0.49,-0.07},{-1,1},{-1,1},{-1,1}};
180 f9daq 369
        int edge[4]={1,1,1,1};// 0 -- rising, 1 -- falling
370
 
371
        DrsChannel c[2];
372
 
373
        c[0].cfrac = CFRAC;
374
        c[0].threshold = threshold[0][0];
375
        c[0].twin[0] = twin[0][0];
376
        c[0].twin[1] = twin[0][1];
377
 
378
        c[0].adcgate[0] = adcgate[0][0];
379
        c[0].adcgate[1] = adcgate[0][1];
380
 
381
        c[0].vcut[0] = vcut[0][0];
382
        c[0].vcut[1] = vcut[0][1];
383
 
384
        c[0].edge = edge[0];
385
 
386
 
387
        c[1].cfrac = CFRAC;
388
        c[1].threshold = threshold[1][0];
389
        c[1].twin[0] = twin[1][0];
390
        c[1].twin[1] = twin[1][1];
391
 
392
        c[1].adcgate[0] = adcgate[1][0];
393
        c[1].adcgate[1] = adcgate[1][1];
394
 
395
        c[1].vcut[0] = vcut[1][0];
396
        c[1].vcut[1] = vcut[1][1];
397
 
398
        c[1].edge = edge[1];
399
        drs * d = new drs(FileName, 2, c, trgch, nev, updfrq );
101 f9daq 400
}