Subversion Repositories f9daq

Rev

Details | Last modification | View Log | RSS feed

Rev Author Line No. Line
326 f9daq 1
// g++  -I`root-config  --incdir` sa02read.C `root-config  --libs` -lz -o sa02read
2
#ifdef WIN32
3
#pragma warning( disable : 4996  )
4
#endif
5
#include <stdlib.h>
6
#include <stdio.h>
7
#include <stdint.h>
8
#include <iostream>
9
#include <fstream>
10
#include <zlib.h>
11
#include <iostream>
12
#include <vector>
13
#include <string>
14
#include "TString.h"
15
#include "TFile.h"
16
#include "TNtuple.h"
17
#include "TH3F.h"
18
#include "TH2F.h"
19
#include "TH1F.h"
20
#include "TH3D.h"
21
#include "TH1D.h"
22
#include "TGraph.h"
23
#include "TROOT.h"
24
#include "TObject.h"
25
 
26
#include "SA02_DEF.h"
27
#include "dataio.h"
28
#include "H3D.h"
29
#include "H2D.h"
30
#include "H1D.h"
31
#ifdef WIN32
32
#  include "Tchar.h"
33
#  include "XGetopt.h"
34
# else /* WIN32 */
35
#  include <getopt.h>
36
#endif  /* WIN32 */
37
 
38
int gposx=-1;
39
int gposy=-1;
40
char gmapfile[256]="";
41
int gMapped=0;
42
int gInitialized=0;
43
unsigned int gEid2HapdXy[144];
44
int gBoardNumber=0;
45
int gIsNsteps=1;
46
int gScanType=0;
47
 
48
MONREC monrec;
49
RUNREC runrec;
50
POSREC posrec;
51
ITERATORREC iteratorrec;
52
DATREC datrec;
53
RUNINFO runinfo;
54
EVTREC evtrec;
55
H3D h3d;
56
H2D h2d;
57
H1D h1d;
58
int V729_first=1;
59
//--------------------------------------------------------------------
60
int ReadMapFile(char *fname) {
61
 
62
  int chip, ch,x,y;
63
  std::ifstream din;
64
  din.open (fname);
65
  for (int i=0; i<144; i++) {
66
    din >> chip >> ch >> x>> y;
67
    x = 11 -x;
68
    y = 11 -y;
69
    gEid2HapdXy[chip*36+ch]=x*12+y;
70
  }
71
  din.close();
72
  return 1;
73
}
74
 
75
//--------------------------------------------------------------------
76
//--------------------------------------------------------------------
77
class readdata {
78
 
79
private:
80
  TGraph *m_GrTime;
81
  TGraph *m_GrData;
82
  TGraph *m_GrEvent;
83
  TNtuple* m_Ntuple;
84
  TNtuple* m_thresh;
85
  TH1F* m_NHit;
86
  TH2F* m_Scan;
87
  TH2F* m_Histo;
88
  TH2F* m_ScanSet;
89
  TH3F* m_Scan3d;
90
  TH2F* m_xy[144];
91
  TH2F* m_xysum;
92
  TTree * m_ttmon;
93
  TNtuple *m_ntpos;
94
  TTree *m_ntrun;
95
public:
96
  int m_neve; // number of events to process
97
  int m_ny;
98
  int m_dy;
99
  int m_y0;
100
  int m_val;
101
  int m_time;
102
  int m_cureve;
103
  int m_y;
104
  int m_towrite;
105
//public:
106
  int m_mask;          // maska za status register
107
  int m_shft;          // shift za status register
108
  int m_debug;         // debug izpisi
109
  int m_write;         // ce je ta flag nastavljen, se v root file izpisuje ntuple
110
  unsigned int m_bitmask ;
111
 
112
//--------------------------------------------------------------------
113
#define MAXCH 144
114
 
115
  int Initialize(int nx,int ny,double ny0,double ny1) {
116
 
117
    printf("Initialize %d %d %2.2f %2.2f\n",nx,ny,ny0,ny1);
118
    if (m_write) {
119
      m_Ntuple = new TNtuple("nt","SA02 RAW data","time:eve:x:y");
120
      m_GrTime = new TGraph(ny);
121
      m_GrTime->SetTitle("grtime");
122
      m_GrEvent= new TGraph(ny);
123
      m_GrEvent->SetTitle("grevent");
124
      m_GrData = new TGraph(ny);
125
      m_GrData->SetTitle("grdata");
126
    }
127
    m_NHit = new TH1F("NHit","N hit per event",nx+1,-0.5,nx+0.5);
128
    m_Scan  = new TH2F("Scan","Scan;SA02 ch;Threshold value",nx,-0.5,nx-0.5, ny,ny0,ny1 );
129
//     m_Scan  = new TH2F("Scan","Scan;SA02 ch;Threshold value",nx,-0.5,nx-0.5, (ny+1.2452)/0.0024,ny0,ny1 );
130
    m_Scan3d  = new TH3F("Scan3d","Scan3d;SA02 ch;bit;Threshold value",nx,-0.5,nx-0.5,8,-0.5,7.5,ny,ny0,ny1 );
131
//     m_Scan3d  = new TH3F("Scan3d","Scan3d;SA02 ch;bit;Threshold value",nx,-0.5,nx-0.5,  8,-0.5,7.5 ,(ny+1.2452)/0.0024,ny0,ny1 );
132
    m_Histo = new TH2F("Histo","Histo;SA02 ch;bit",nx,-0.5,nx-0.5, 8,-0.5,7.5 );
133
    m_ScanSet  = new TH2F("ScanSet","Scan;SA02 ch;set value",nx,-0.5,nx-0.5, ny,ny0,ny1 );
134
    m_thresh = new TNtuple("thresh","bo","ch:y:xl:yl");
135
    return 1;
136
  }
137
 
138
//--------------------------------------------------------------------
139
  int FillHistograms(unsigned int *data, double yset, double y) {
140
 
141
    static int trgn = 0;
142
    int nhit=0;
143
 
144
    for (int ch=0; ch<144; ch++) {
145
      int id=(17-ch/8)+gBoardNumber*18;
146
      int shft=(ch%8)*4;
147
      for (int i=0; i<4; i++) {
148
        if (data[id] & (1<<(i+shft)) )  {
149
          if ( m_Scan3d ) m_Scan3d->Fill(ch,i,y);
150
          if ( m_Histo  ) m_Histo->Fill(ch,i);
151
        }
152
      }
153
      if (data[id] & (m_bitmask<<shft) ) {
154
        nhit++;
155
        if (m_Scan) m_Scan->Fill(ch,y);
156
        if (m_ScanSet) m_ScanSet->Fill(ch,yset);
157
        unsigned int ii= gEid2HapdXy[ch];
158
        if (m_xy[ii]) {
159
          double xs;
160
          double ys;
161
          if (runrec.fver) {
162
            xs = posrec.ix*runrec.dx+runrec.x0;
163
            ys = posrec.iy*runrec.dy+runrec.y0;
164
          } else {
165
            xs = posrec.xset;
166
            ys = posrec.yset;
167
          }
168
          if (m_xysum) m_xysum->Fill( xs, ys);
169
//!          m_xy[ii]->Fill(xs, ys);
170
          m_xy[ii]->Fill(posrec.xset, posrec.yset);
171
//          printf("%3d c=%d ch=%2d x=%6d y=%6d\n",ch, ii/36,ii%36,posrec.set_pos_x, posrec.set_pos_y);
172
        } else {
173
          printf("error m_xy[%d] does not exist\n",ii);
174
        }
175
        m_thresh->Fill(ch,y,trgn,1);
176
      }
177
    }
178
    if (m_NHit) m_NHit->Fill(nhit);
179
    trgn++;
180
    return 1;
181
  }
182
 
183
//--------------------------------------------------------------------
184
  int ProcessData(unsigned int *rdata,int count,int neve) {
185
 
186
    int nc=0;
187
    for (int j=0; j<count;) {
188
      unsigned int recid   = rdata[j++];
189
      unsigned int len     = rdata[j++];
190
      if (m_debug) printf(" recid=0x%0x len=%d pos=%d(maxpos %d) val=%d\n",recid, len, j, count, rdata[j]);
191
      if (m_debug>2 && len>2) printf("\n");
192
      if(j+(int)len>=count) return nc;
193
      switch (recid) {
194
        case 3:
195
          if (len !=18) printf(" recid=0x%0x len=%d pos=%d(maxpos %d) val=%d\n",recid, len, j, count, rdata[j]);
196
          nc+=FillHistograms(&rdata[j], m_y, m_val);
197
          if (nc==neve) return nc;
198
          break;
199
        default:
200
          m_val=(rdata[j] >> m_shft) & m_mask ;
201
          if (m_write) {
202
            m_Ntuple->Fill(m_time,m_cureve,m_y,m_val);
203
            m_GrTime->SetPoint(m_cureve,m_time,m_val);
204
            m_GrEvent->SetPoint(m_cureve,m_cureve,m_val);
205
            m_GrData->SetPoint(m_cureve,m_y,m_val);
206
          }
207
          break;
208
      }
209
      if (m_debug>2) {
210
        for (unsigned int i=0; i<len; i++) {
211
          if (j<count) printf("0x%0x ", rdata[j]);
212
          if (i%4==3) printf("\n");
213
          j++;
214
        }
215
        printf("\n");
216
      } else j+=len;
217
//      if(j+4>=count) return nc;
218
    }
219
 
220
    return nc;
221
  }
222
 
223
//--------------------------------------------------------------------
224
  double H3DGetBinContentFrom( H3D *h ,unsigned int atx, unsigned int aty, unsigned int atz) {
225
 
226
    unsigned int idx;
227
    if (!h) return 0;
228
    if (h->nx <= atx)  return 0;
229
    if (h->ny <= aty)  return 0;
230
    if (h->nz <= atz)  return 0;
231
    idx = atx+(aty+atz*h->ny)*h->nx;
232
    if (idx*sizeof(double)  < h->size ) {
233
//    if (h->data[idx]>0)  printf("xy %d %d %g\n",atx,aty,h->data[idx]);
234
      return h->data[idx];
235
    }
236
    return 0;
237
  }
238
 
239
//--------------------------------------------------------------------
240
  double H2DGetBinContentFrom( H2D *h ,unsigned int atx, unsigned int aty) {
241
 
242
    unsigned int idx;
243
    if (!h) return 0;
244
    if (h->nx <= atx)  return 0;
245
    if (h->ny <= aty)  return 0;
246
    idx = atx+aty*h->nx;
247
    if (idx*sizeof(double)  < h->size ) {
248
//    if (h->data[idx]>0)  printf("xy %d %d %g\n",atx,aty,h->data[idx]);
249
      return h->data[idx];
250
    }
251
    return 0;
252
  }
253
 
254
//--------------------------------------------------------------------
255
  double H1DGetBinContentFrom( H1D *h ,unsigned int atx) {
256
 
257
    if (!h) return 0;
258
    if (h->nx <= atx)  return 0;
259
    if (atx*sizeof(double)  < h->size ) {
260
      return h->data[atx];
261
    }
262
    return 0;
263
  }
264
 
265
//--------------------------------------------------------------------
266
  TH2D *hwf[4];
267
  int V729_init() {
268
    printf(" V729_init\n");
269
    const double maxy=512;
270
//    const double maxy=512;
271
    const int ny=100;
272
//    const int h0=10;
273
    int len=100;
274
 
275
    TObject *hobj = gROOT->FindObject("gwf0");
276
    if (hobj !=NULL) hobj->Write();
277
    hobj = gROOT->FindObject("gwf1");
278
    if (hobj !=NULL) hobj->Write();
279
    hobj = gROOT->FindObject("gwf2");
280
    if (hobj !=NULL) hobj->Write();
281
    hobj = gROOT->FindObject("gwf3");
282
    if (hobj !=NULL) hobj->Write();
283
    hwf[0]= new TH2D("gwf0","CAENV_729 ch0",len,-0.5,len-0.5,ny,-maxy,maxy);
284
    hwf[1]= new TH2D("gwf1","CAENV_729 ch1",len,-0.5,len-0.5,ny,-maxy,maxy);
285
    hwf[2]= new TH2D("gwf2","CAENV_729 ch2",len,-0.5,len-0.5,ny,-maxy,maxy);
286
    hwf[3]= new TH2D("gwf3","CAENV_729 ch3",len,-0.5,len-0.5,ny,-maxy,maxy);
287
//    h1d=new TH1F("h1d","Charge;signal charge(a.u.);N",500,-maxy,5*maxy);
288
    return 0;
289
  }
290
 
291
//--------------------------------------------------------------------
292
  int V729_plot(uint16_t *buf ) {
293
 
294
        if (V729_first){
295
          V729_init();
296
      V729_first=0;      
297
        }
298
    int start=310;
299
    int end=335;
300
//    const double maxy=2048;
301
//    const double maxy=512;
302
//    const double maxy=512;
303
//    const int ny=100;
304
//    const int h0=10;
305
    int i;
306
    uint32_t *hdr = (uint32_t *)buf;
307
    int nb= hdr[1]-2*sizeof(uint32_t);
308
    int nitems = nb /sizeof(uint16_t);
309
    unsigned int recid=hdr[0];
310
    uint16_t *data = &buf[4];
311
    double x,y;
312
//    double baseline=0;
313
    double charge=0;
314
    while (nitems>0) {
315
      int id = data[0]-1;
316
      int len= data[1]/sizeof(uint16_t)-2;
317
      if (m_debug)  printf("\t->> recid=%d nb=%d nitems=%d id=%d len=%d\n", recid, hdr[1],nitems,id,len);
318
      if (id>3) return 0;
319
      data+=2;
320
//      for (i=0; i<200;i++) baseline += data[i];
321
//      baseline /=200.;
322
      for (i=0; i<len; i++) {
323
        if (m_debug >2) printf("ADC %d %d\n", i, data[i]);
324
        x=i;
325
        y=data[i];
326
//        y-=baseline;
327
        if (i>start && i<end) charge += y;
328
 
329
        hwf[id]->Fill(x,y-2000);
330
      }
331
      data+=len;
332
      nitems=nitems - len-2;
333
//      if (id==2) h1d->Fill(charge);
334
    }
335
    return 0;
336
  }
337
 
338
//--------------------------------------------------------------------
339
  int Process( const char *fnamein, int nevetoprocess, int histoonly) {
340
 
341
    int nxh,nyh;
342
    gzFile fp=gzopen(fnamein,"r");
343
    printf("Process....\n");
344
    if (!fp) {
345
      printf("File %s cannot be opened!\n", fnamein);
346
      return 0;
347
    } else {
348
      printf("Processing file %s\n", fnamein);
349
    }
350
    unsigned int bsize=1000000; // initial buffer size
351
    unsigned int *buf = new unsigned int[bsize];
352
    unsigned int *data= buf+4; // data shift
353
    int nr=0;
354
    HEADERREC *hdr = (HEADERREC *) buf;
355
//    int nread=0;
356
    while (!gzeof(fp) ) {
357
      int nitems  = gzread(fp,buf,sizeof(HEADERREC));
358
      if (nitems !=sizeof(HEADERREC)) break;
359
      int recid   = hdr->id;
360
      int nb      = hdr->len - sizeof(HEADERREC);
361
      if (hdr->len > bsize*sizeof(int)) {
362
        printf ("[INFO] Requested data size=%u bigger than buffer size %u\n",hdr->len/sizeof(uint32_t), bsize);
363
        bsize=2*hdr->len/sizeof(int);
364
        printf("[INFO] Increasing data buffer to %d elements\n", bsize);
365
        delete buf;
366
        buf = new unsigned int[bsize];
367
        buf[0] = recid;
368
        buf[1] = nb + sizeof(HEADERREC);
369
        hdr = (HEADERREC *) buf;
370
        data = buf +4;
371
      }
372
      nitems=gzread(fp,((unsigned char *)buf)+sizeof(HEADERREC),nb);
373
      nr++;
374
      if (m_debug || nr<30) printf("[RECID=%u] Length=%u items=%u\n",hdr->id,hdr->len,hdr->len/sizeof(uint32_t));
375
      if (nitems!=nb) {
376
        printf("Read Error nb=%d nitems=%d. Exiting ...\n",nb,nitems);
377
        break;
378
//        exit(0);
379
      }
380
      switch (recid) {
381
        case POSREC_ID: {
382
          memcpy(&posrec.id, buf, sizeof(posrec));
383
          printf("POSRECID = %u pos= %d %d\n",posrec.id, posrec.ix,posrec.iy);
384
          //if (!histoonly)  V729_init();
385
          if (m_debug || nr<20) {
386
            printf("Length = %u\n",posrec.len);
387
            printf("Time = %u\n",posrec.time);
388
            printf("Iteration X = %d\n",posrec.ix);
389
            printf("MIKRO Position X = %d\n",posrec.x);
390
            printf("Set position X = %d\n",posrec.xset);
391
            printf("Iteration Y = %d\n",posrec.iy);
392
            printf("MIKRO Position Y = %d\n",posrec.y);
393
            printf("Set position Y = %d\n",posrec.yset);
394
          }
395
                m_ntpos->Fill(posrec.time-runrec.time, posrec.ix,posrec.x,posrec.xset, posrec.iy,posrec.y,posrec.yset);
396
        }
397
        break;
398
 
399
 
400
                case MONREC_ID : {
401
                        memcpy(&monrec.id, buf, sizeof(monrec));
402
 
403
                        //m_ttmon->Branch("monrec",&monrec,"id/i:len:time:imon[6]/I:vmon[6]:status[6]");
404
                        m_ttmon->Branch("id",&monrec.id,"id/i");
405
                        m_ttmon->Branch("len",&monrec.len,"len/i");
406
                        m_ttmon->Branch("time",&monrec.time,"time/i");
407
                        m_ttmon->Branch("imon",&monrec.imon[0],"imon[24]/I");
408
                        m_ttmon->Branch("iset",&monrec.iset[0],"iset[24]/I");
409
                        m_ttmon->Branch("vmon",&monrec.vmon[0],"vmon[24]/I");
410
                        m_ttmon->Branch("vset",&monrec.vset[0],"vset[24]/I");
411
                        m_ttmon->Branch("status",&monrec.status[0],"status[24]/I");
412
                        m_ttmon->Fill();  
413
 
414
                        break;
415
                }
416
 
417
        case RUNREC_ID: {
418
 
419
                  memcpy(&runrec.id, buf, sizeof(runrec));
420
                        m_ntrun->Branch("time",&runrec.time,"time/i");
421
                        m_ntrun->Branch("fver",&runrec.fver,"fver/i");
422
                        m_ntrun->Branch("nev",&runrec.nev,"nev/i");
423
                        m_ntrun->Branch("ped",&runrec.ped,"ped/i");
424
                        m_ntrun->Branch("nx",&runrec.nx,"nx/I");
425
                        m_ntrun->Branch("x0",&runrec.x0,"x0/I");
426
                        m_ntrun->Branch("dx",&runrec.dx,"dx/I");
427
                        m_ntrun->Branch("ny",&runrec.ny,"ny/I");
428
                        m_ntrun->Branch("y0",&runrec.y0,"y0/I");
429
                        m_ntrun->Branch("dy",&runrec.dy,"dy/I");
430
                        m_ntrun->Fill();
431
 
432
          printf("***************RECID = %u\n",runrec.id);
433
          printf("Length = %u\n",runrec.len);
434
          printf("File version = %u\n",runrec.fver);
435
          printf("Time = %u\n",runrec.time);
436
          printf("Number of events per step = %u\n",runrec.nev);
437
          printf("Pedestal = %u\n",runrec.ped);
438
//          printf("Scan type = %u :: 0 -> single data scan :: 1 -> XY position scan\n",scanrec.xy);
439
          printf("Number of steps in X = %d\n",runrec.nx);
440
          printf("Start position X = %d\n",runrec.x0);
441
          printf("Step size direction X = %d\n",runrec.dx);
442
          printf("Number of steps in Y = %d\n",runrec.ny);
443
          printf("Start position Y = %d\n",runrec.y0);
444
          printf("Step size direction Y = %d\n",runrec.dy);
445
          if (gIsNsteps) {
446
            nxh=runrec.nx;
447
            nyh=runrec.ny;
448
          } else {
449
            nxh=(runrec.nx-1)/runrec.dx+1;
450
            nyh=(runrec.ny-1)/runrec.dy+1;
451
          }  
452
          for (int i=0;i<144;i++) {
453
            char hname[256];
454
            if (gMapped) sprintf(hname,"ch_%02d_%02d",i%12,i/12);
455
            else sprintf(hname,"xy_%02d_%02d",i/36,i%36);
456
            m_xy[i]=new TH2F(hname,hname,nxh,runrec.x0-runrec.dx*0.5,runrec.x0+runrec.dx*(nxh-0.5),
457
                                         nyh,runrec.y0-runrec.dy*0.5,runrec.y0+runrec.dy*(nyh-0.5));
458
          }
459
          m_xysum= (TH2F *) (m_xy[0]->Clone("xy"));
460
          m_xysum->SetTitle("xy all channels");
461
        }
462
        break;
463
        case RUNINFO_ID: {
464
          memcpy(&runinfo.id, buf, sizeof(runinfo));
465
          if (!gInitialized) {
466
            m_y0  = runinfo.x0;
467
            m_dy  = runinfo.dx;
468
            m_ny  = runinfo.nx;
469
            m_towrite = runinfo.writemode;
470
            if (m_towrite ==3) m_y0=runinfo.chip*36+runinfo.ch;
471
            printf("RUNINFO x0=%d nx=%d dx=%d\n", runinfo.x0,runinfo.dx,runinfo.nx);
472
            gInitialized=Initialize(144,m_ny,m_y0-0.5,m_y0-0.5+m_dy*m_ny);
473
          }
474
        }
475
        break;
476
        case EVTREC_ID:
477
          memcpy(&evtrec.id, buf, sizeof(evtrec));
478
          m_time   = evtrec.time;
479
          m_cureve = evtrec.nev;
480
          m_y = m_y0 +m_cureve * m_dy;
481
//          nread++;
482
          break;
483
//        case ITERATORREC_ID:
484
//          memcpy(&iteratorrec.id, buf, sizeof(iteratorrec));
485
//          printf("ITERATORREC level=%d n=%d value=%f step=%f\n", iteratorrec.level,iteratorrec.n,iteratorrec.value, iteratorrec.step);
486
//          break;
487
        case DATREC_ID:
488
          memcpy(&datrec.id, buf, sizeof(datrec));
489
          printf("DATREC chip=%d channel=%d cmd=%d data=%d retval=%d\n", datrec.chip,datrec.channel,datrec.cmd, datrec.data, datrec.retval);
490
          break;
491
        case CAENV729REC_ID: {
492
          if (!histoonly) V729_plot( (uint16_t*)buf );
493
        }
494
        break;
495
        case H3D_ID: {
496
          int hlen= sizeof(h3d)-sizeof(double *);
497
          memcpy(&h3d.id, buf, hlen);
498
          h3d.data = new double[h3d.size/sizeof(double)];
499
//          int align=0;
500
//          if (sizeof(double *)==4) align=0;
501
//          memcpy(h3d.data, ((char *) buf)+hlen-align, h3d.size);
502
          memcpy(h2d.data, ((char *) buf)+hlen, h2d.size);
503
          if (m_debug) printf("H3D %s nx=%d minx=%f stepx=%f ny=%d miny=%f stepy=%f  nz=%d minz=%f stepz=%f  min=%f max=%f size=%d title='%s' x='%s' y='%s' z='%s'  \n", h3d.name, h3d.nx, h3d.minx, h3d.stepx, h3d.ny, h3d.miny, h3d.stepy,h3d.nz, h3d.minz, h3d.stepz, h3d.min, h3d.max, h3d.size, h3d.title, h3d.titlex, h3d.titley, h3d.titlez );
504
//          printf("H3D sizeof(H2D)=%lu len-datasize=%d len=%lu datasize=%d\t",hlen,h3d.len-h3d.size,h3d.len,h3d.size);
505
//          printf("H3D sz=%d\n",sizeof(double *));
506
          H3D  *h = &h3d;
507
          char hname[100];
508
          sprintf(hname,"%s_%d",h->name,datrec.data);
509
          TObject *hobj = gROOT->FindObject(hname);
510
          if (hobj !=NULL) {
511
            printf("TH3D duplicate name %s!\n",hname);
512
            hobj->Write();
513
            delete hobj;
514
//            break;
515
          }
516
          if (h->stepx==0) h->stepx=1;
517
          if (h->stepy==0) h->stepy=1;
518
          if (h->stepz==0) h->stepz=1;
519
          TH3D *h3 = new TH3D(hname, h->title,
520
                              h->nx, h->minx-0.5*h->stepx,h->minx+(h->nx-0.5)*h->stepx,
521
                              h->ny, h->miny-0.5*h->stepy,h->miny+(h->ny-0.5)*h->stepy,
522
                              h->nz, h->minz-0.5*h->stepz,h->minz+(h->nz-0.5)*h->stepz);
523
          h3->SetTitle(h->title);
524
          h3->GetXaxis()->SetTitle(h->titlex);
525
          h3->GetYaxis()->SetTitle(h->titley);
526
          h3->GetZaxis()->SetTitle(h->titlez);
527
          if (m_debug) {
528
            printf("TH3D name=%s %s %s %s %s",hname, h->title,h->titlex, h->titley, h->titlez);
529
            printf("\tnx=%d %f %f \tny=%d %f %f\tnz=%d %f %f\n",
530
                   h->nx, h->minx-0.5*h->stepx,h->minx+(h->nx-0.5)*h->stepx,
531
                   h->ny, h->miny-0.5*h->stepy,h->miny+(h->ny-0.5)*h->stepy,
532
                   h->nz, h->minz-0.5*h->stepz,h->minz+(h->nz-0.5)*h->stepz);
533
          }
534
          for (unsigned int ix=0; ix<h->nx; ix++) {
535
            for (unsigned int iy=0; iy<h->ny; iy++) {
536
              for (unsigned int iz=0; iz<h->nz; iz++) {
537
                double val = H3DGetBinContentFrom(h,ix,iy,iz);
538
                if (m_debug>2) printf("%d %d %d %g\n", ix,iy,iz,val);
539
                h3->SetBinContent(ix+1,iy+1,iz+1,val);
540
              }
541
            }
542
          }
543
          if (m_debug) {
544
            h3->ls();
545
            h3->Print();
546
            printf("TH3D Sum %f\n", h3->GetMaximum());
547
          }
548
        }
549
        break;
550
        case H2D_ID: {
551
          if(gScanType==1) break;
552
          int hlen= sizeof(h2d)-sizeof(double *);
553
          memcpy(&h2d.id, buf, hlen);
554
          h2d.data = new double[h2d.size/sizeof(double)];
555
//          int align=0;
556
//          if (sizeof(double *)==4) align=0;
557
//          memcpy(h2d.data, ((char *) buf)+hlen-align, h2d.size);
558
          memcpy(h2d.data, ((char *) buf)+hlen, h2d.size);
559
          if (m_debug) printf("H2D nx=%d minx=%f stepx=%f ny=%d miny=%f stepy=%f min=%f max=%f size=%d title='%s' x='%s' y='%s' \n", h2d.nx, h2d.minx, h2d.stepx, h2d.ny, h2d.miny, h2d.stepy, h2d.min, h2d.max, h2d.size, h2d.title, h2d.titlex, h2d.titley );
560
//          printf("H2D sizeof(H2D)=%lu len-datasize=%d len=%lu datasize=%d\t",hlen,h2d.len-h2d.size,h2d.len,h2d.size);
561
//          printf("H2D sz=%d\n",sizeof(double *));
562
          H2D  *h = &h2d;
563
          char hname[100];
564
          sprintf(hname,"%s_%d",h->name,datrec.data);
565
          TObject *hobj = gROOT->FindObject(hname);
566
          if (hobj !=NULL) {
567
            printf("TH2D duplicate name %s!\n",hname);
568
            hobj->Write();
569
            delete hobj;
570
//            break;
571
          }
572
          if (h->stepx==0) h->stepx=1;
573
          if (h->stepy==0) h->stepy=1;
574
          TH2D *h2 = new TH2D(hname, h->title,
575
                              h->nx, h->minx-0.5*h->stepx,h->minx+(h->nx-0.5)*h->stepx,
576
                              h->ny, h->miny-0.5*h->stepy,h->miny+(h->ny-0.5)*h->stepy);
577
          h2->SetTitle(h->title);
578
          h2->GetXaxis()->SetTitle(h->titlex);
579
          h2->GetYaxis()->SetTitle(h->titley);
580
          if (m_debug) printf("TH2D name=%s %s %s %s\tnx=%d %g %g \tny=%d %g %g\n",hname, h->title,h->titlex, h->titley,
581
                                h->nx, h->minx-0.5*h->stepx,h->minx+(h->nx-0.5)*h->stepx,
582
                                h->ny, h->miny-0.5*h->stepy,h->miny+(h->ny-0.5)*h->stepy);
583
          for (unsigned int ix=0; ix<h->nx; ix++) {
584
            for (unsigned int iy=0; iy<h->ny; iy++) {
585
              double val = H2DGetBinContentFrom(h,ix,iy);
586
              if (m_debug>2) printf("%d %d %g\n", ix,iy,val);
587
              h2->SetBinContent(ix+1,iy+1,val);
588
            }
589
          }
590
          if (m_debug) {
591
            h2->ls();
592
            h2->Print();
593
            printf("TH2D Sum %f\n", h2->GetMaximum());
594
          }
595
        }
596
        break;
597
        case H1D_ID: {
598
          int hlen= sizeof(h1d)-sizeof(double *);
599
          memcpy(&h1d.id, buf, hlen);
600
          h1d.data = new double[h1d.size/sizeof(double)];
601
//          int align=0;
602
//          if (sizeof(double *)==4) align=0;
603
//          memcpy(h1d.data, ((char *) buf)+hlen-align, h1d.size);
604
          memcpy(h1d.data, ((char *) buf)+hlen, h1d.size);
605
          if (m_debug) printf("H1D nx=%d minx=%f stepx=%f min=%f max=%f size=%d title='%s' x='%s' y='%s' \n", h1d.nx, h1d.minx, h1d.stepx, h1d.min, h1d.max, h1d.size, h1d.title, h1d.titlex, h1d.titley );
606
          H1D  *h = &h1d;
607
          char hname[100];
608
          sprintf(hname,"%s_%d",h->name,datrec.data);
609
          TObject *hobj = gROOT->FindObject(hname);
610
          if (hobj !=NULL) {
611
            printf("TH1D duplicate name %s!\n",hname);
612
            hobj->Write();
613
            delete hobj;
614
//            break;
615
          }
616
          if (h->stepx==0) h->stepx=1;
617
          TH1D *h1 = new TH1D(hname, h->title, h->nx, h->minx-0.5*h->stepx,h->minx+(h->nx-0.5)*h->stepx);
618
          h1->SetTitle(h->title);
619
          h1->GetXaxis()->SetTitle(h->titlex);
620
          h1->GetYaxis()->SetTitle(h->titley);
621
//          h1->GetXaxis()->SetTitle("timebin");
622
//          h1->GetYaxis()->SetTitle(h->titlex);
623
          if (m_debug) printf("TH1D name=%s %s %s %s\tnx=%d %g %g\n",hname, h->title,h->titlex, h->titley,
624
                                h->nx, h->minx-0.5*h->stepx,h->minx+(h->nx-0.5)*h->stepx);
625
          for (unsigned int ix=0; ix<h->nx; ix++) {
626
            double val = H1DGetBinContentFrom(h,ix);
627
            if (m_debug>2) printf("%d %g\n", ix,val);
628
            h1->SetBinContent(ix+1,val);
629
          }
630
          if (m_debug) {
631
            h1->ls();
632
            h1->Print();
633
            printf("TH1D Sum %f\n", h1->GetMaximum());
634
          }
635
        }
636
        break;
637
        default:
638
          printf("Unknown record 0x%x Exiting .... \n", recid);
639
                  gzclose(fp);
640
          delete buf;
641
          return nr;
642
 
643
      }
644
      if (recid!=EVTREC_ID)  continue;
645
      if (histoonly) continue;
646
      if (m_debug) printf("EVTREC %d histo=%d\n",recid,histoonly);
647
      if (gposx>=0 && gposx!=posrec.ix ) continue;
648
      if (gposy>=0 && gposy!=posrec.iy ) continue;
649
      int nc=ProcessData(data,nb/sizeof(unsigned int)-2,nevetoprocess);
650
      if (m_debug) printf("EVTREC Events %d nb=%u\n",nc,nb/sizeof(unsigned int)-2);
651
//      if (nread++==nevetoprocess) break;
652
    } //  while (!gzeof(fp))
653
    gzclose(fp);
654
    delete buf;
655
    return nr;
656
  }
657
 
658
//--------------------------------------------------------------------
659
  readdata( std::vector<TString> &fnamein,const char *fnameout,int write2root=0,
660
            int nevetoread=-1,unsigned int bmask=0xff,int mask=0xFFFFFFFF,int shift=0,
661
            int neve=-1,int debug=0,int histoonly=0) {
662
 
663
    m_debug = debug;
664
    m_write = write2root;
665
    m_neve  = nevetoread;
666
    m_mask  = mask;
667
    m_shft =  shift;
668
    m_towrite = 0;
669
    m_bitmask = bmask;
670
    char rootname[0xFF];
671
    sprintf(rootname,"%s",fnameout);
672
    TFile *rootFile = new TFile(rootname,"RECREATE");
673
        m_ntpos = new TNtuple("pos","Position record", "time:ix:x:xset:iy:y:yset");
674
        m_ntrun = new TTree("run","Run record");
675
        m_ttmon = new TTree("Monitor","Voltage, Current and Status Monitor");
676
 
677
//    Initialize();
678
    printf("Konec inicializacije ...\n");
679
    if (m_debug) printf("readdata histo=%d\n",histoonly);
680
    int nr = 0;
681
    for (unsigned int i=0; i<  fnamein.size() ; i++) {
682
      int neve = m_neve -nr;
683
      if (m_neve == -1)  nr += Process(fnamein[i].Data(), m_neve, histoonly);
684
      else if (neve>0)   nr += Process(fnamein[i].Data(), neve, histoonly);
685
    }
686
    printf("End...\nreaddata::Number of events read=%d required=%d\n",nr,m_neve);
687
    if (m_write) {
688
      m_GrTime->Write();
689
      m_GrData->Write();
690
      m_GrEvent->Write();
691
    }
692
    rootFile->Write();
693
    if (m_debug) rootFile->ls();
694
    rootFile->Close();
695
  }
696
 
697
//--------------------------------------------------------------------
698
  ~readdata() { };
699
};
700
 
701
//--------------------------------------------------------------------
702
//-----------------------------------------------------------------
703
#ifdef MAIN
704
 
705
int main(int argc, char **argv) {
706
 
707
  std::cout << "Usage:" << argv[0] <<  "  -i input.dat -o output -n nevents -w writentuple -d debuglevel -s shift -m mask -l boardnumber -h histoonly -f electronicmap.map "<< std::endl;
708
  char outputfile[256]="/tmp/sa02asic.root";
709
  int nevetoread      =-1;
710
  int writentuple     =0;
711
  int verbose         =0;
712
  int shift=0;
713
  int histoonly=0;
714
  Int_t mask =0xFFFFFFFF;
715
  unsigned int bitmask=0xf;
716
  char c;
717
  for (int i=0; i<144; i++) gEid2HapdXy[i]=i;
718
  std::vector<TString> inputfilelist;
719
  while ((c=getopt (argc, argv, "l:i:I:o:n:w:d:m:s:x:y:f:b:h:t:")) != -1) {
720
    switch (c) {
721
      case 'I': // runrec.nx = xmax-xmin - for magnet tests in folder 2014
722
        gIsNsteps=0;
723
      case 'i': // runrec.nx == x steps - standard
724
        inputfilelist.push_back(TString(optarg));
725
        break;
726
      case 'o':
727
        sprintf(outputfile,"%s",optarg);
728
        break;
729
      case 'n':
730
        nevetoread = atoi(optarg)    ;
731
        break;
732
      case 'w':
733
        writentuple = atoi(optarg);
734
        break;
735
      case 'd':
736
        verbose = atoi(optarg);
737
        break;
738
      case 'l':
739
        gBoardNumber = atoi(optarg);
740
        break;
741
      case 's':
742
        shift = strtoul(optarg,NULL,0);
743
        break;
744
      case 'm':
745
        mask  = strtoul(optarg,NULL,0);
746
        break;
747
      case 'b':
748
        bitmask  = strtoul(optarg,NULL,0);
749
        break;
750
      case 'x':
751
        gposx = atoi(optarg)    ;
752
        break;
753
      case 'y':
754
        gposy = atoi(optarg)    ;
755
        break;
756
      case 'h':
757
        histoonly = atoi(optarg)    ;
758
        break;
759
      case 'f':
760
        sprintf(gmapfile,"%s", optarg)    ;
761
        gMapped=ReadMapFile(gmapfile);
762
        break;
763
      case 't':
764
        gScanType = atoi(optarg);
765
        break;
766
      default:
767
        abort();
768
    }
769
  }
770
  std::cout << "Used: " << argv[0]  << " -o " << outputfile
771
            << " -n " << nevetoread << " -w " << writentuple
772
            << " -d  " << verbose  << " -m  0x" << std::hex << mask
773
            << " -s  0x" << std::hex << shift  << " -l  " << gBoardNumber << " -h  " << histoonly ;
774
  for (unsigned int i=0; i<  inputfilelist.size() ; i++) std::cout << " -i " << inputfilelist[i];
775
  std::cout << std::endl;
776
  if (!inputfilelist.size()) {
777
    std::cout << "No input file specified" << std::endl;
778
  } else {
779
    new readdata(inputfilelist,outputfile,writentuple,nevetoread,bitmask,mask,shift,-1,verbose,histoonly);
780
    std::cout << "Output data files " << outputfile <<  std::endl;
781
  }
782
  return 0;
783
}
784
#endif /* MAIN */