Subversion Repositories f9daq

Rev

Rev 1 | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
1 f9daq 1
/*
2
 
3
vhodne datoteke:
4
 
5
Mappingi:
6
m16.map channel mapping
7
modules.map mapping modulov
8
 
9
Kalibracija:
10
pedestals.dat  ... adc suma
11
ph511.dat      ... adc fotovrhov vrhov
12
 
13
 
14
// sumpedestals.dat vsota adc po modulih
15
 
16
./readdata -c config.xml -i input.dat -o output -n nevents -w writentuple -d debuglevel
17
 
18
*/
19
#include <stdlib.h>
20
#include <stdio.h>
21
#include "TNtuple.h"
22
#include "TH1F.h"
23
#include "TH2F.h"
24
#include "TF1.h"
25
#include "TCanvas.h"
26
#include "TStyle.h"
27
#include "TFile.h"
28
#include "TMath.h"
29
#include "TRandom3.h"
30
#include <zlib.h>
31
#include <iostream>
32
#include <fstream>
33
 
34
 
35
#include <stdlib.h>
36
#include <stdio.h>
37
#include <string.h>
38
#include <assert.h>
39
#include <vector>
40
 
41
#include <libxml/tree.h>
42
#include <libxml/parser.h>
43
#include <libxml/xpath.h>
44
#include <libxml/xpathInternals.h>
45
 
46
#include "PETProjDataMgr.h"
47
 
48
 
49
// en krog je 4790 korakov
50
#define NSTEPS 4790
51
// nt->Draw("px1/max1:py1/max1>>hmy(200,-1,4,200,-1,4)","max1>500","colz")
52
//#define MAXCH 64
53
#define MAXCH 72
54
 
55
 
56
#define RUNREC_ID                       1
57
#define EVTREC_ID                       2
58
#define POSREC_ID                       3
59
 
60
typedef struct {
61
  unsigned int num_events,num_channels,pedestal,xy;
62
  int nx,x0,dx,ny,y0,dy;
63
} RUNREC;
64
RUNREC *runrec;  // start header: appears once at the start of the file
65
 
66
typedef struct {
67
  int mikro_pos_x,set_pos_x;
68
  int num_iter_y,mikro_pos_y,set_pos_y;
69
} POSREC;
70
POSREC *posrec;  // position header: appears at every change of position
71
 
72
class Channel {
73
public:
74
  Channel(){};
75
  ~Channel(){};
76
  int ix;
77
  int iy;
78
  int idx;
79
  void SetIxIy(int a,int b){ix=a;iy=b;idx=ix+4*iy;};
80
};
81
 
82
class Module {
83
public:
84
  Module(){};
85
  ~Module(){};
86
  double r;
87
  double phi;
88
  void SetRPhi(double rr, double pphi){r=rr;phi=pphi*TMath::Pi()/180.;};
89
};
90
 
91
class Geometry {
92
 
93
public:
94
 
95
 
96
  ~Geometry(){ m_calibration->Close(); };
97
  Channel ch[16];
98
  Module  module[4];
99
 
100
 
101
  TRandom3 *m_rnd;   // random generator
102
  TH1I     *m_crystalid[4];// kalibracijski histogrami
103
  char * m_modulesmapName;  // ime datoteke s podatki o pozicijah modulov
104
  char * m_channelsmapName; // ime datoteke s podatki o pozicijah elektronskih kanalov
105
    //char * m_sumpedestalsName;
106
  char * m_pedestalsName;   // ime datoteke s podatki o pedestalih
107
  char * m_photopeakName;   // ime datoteke s podatki o vrhovih
108
  char * m_calibrationrootName;   // ime datoteke s kalibracijskimi histogrami
109
  TFile *m_calibration; //  datoteke s kalibracijskimi histogrami
110
 
111
  double m_dx;        // pitch x kristalov
112
  double m_dy;        // pitch y kristalov
113
  int    m_nx;        // stevilo kristalov v x smeri
114
  int    m_ny;        // stevilo kristalov v y smeri
115
 
116
  float gPedestals[MAXCH];
117
  float gPeak[MAXCH];
118
 
119
  //float gSumPedestals[MAXCH];
120
  double m_peakscaling;   // scaling za adcje
121
  double m_threshold;    // threshold za upostevanje zadetkov na kanalih
122
  int    m_debug;
123
 
124
 
125
  //---------------------------------------------------
126
  int SetThreshold(int threshold){
127
    m_threshold=threshold;
128
    return 0;
129
  };
130
  double GetThreshold(){
131
    return m_threshold;
132
  }
133
 
134
 
135
  int readfile(const char *fname, float *x, int nmax, int defaultvalue=0){
136
    int id;
137
    float ix;
138
    std::ifstream f(fname);
139
    if (f.is_open()){
140
      std::cout << "Reading ... " <<  fname <<  std::endl;
141
      while (f.good()){
142
	f >> id >> ix;
143
	if (id<nmax)  x[id]=ix;
144
	if (m_debug) std::cout << fname << " " << id << " " << ix <<" "<< std::endl;
145
      }
146
      f.close();
147
    } else {
148
      std::cout << "Cannot read " <<  fname <<  std::endl;
149
      std::cout << "Default value will be used :" << defaultvalue  <<  std::endl;
150
      for (int i=0;i<nmax;i++){
151
        x[i]=defaultvalue;
152
      }
153
    }
154
    return 0;
155
  };
156
 
157
  void ReadChannelMap(const char *fname){
158
    int id,ix,iy;
159
    char line[400];
160
    std::ifstream f(fname);
161
    if (f.is_open()){
162
      std::cout << "Reading ... " <<  fname <<  std::endl;
163
      f.getline(line,400);
164
      while (f.good()){
165
	f >> id >> ix >> iy;
166
	if (id<16) ch[id].SetIxIy(ix,iy);
167
	if (m_debug) std::cout << fname << " " << id << " " << ix <<" "<< iy << std::endl;
168
      }
169
      f.close();
170
    } else {
171
       std::cout << "Cannot read " <<  fname <<  std::endl;
172
    }
173
 
174
  };
175
 
176
  void ReadModuleMap(const char *fname){
177
    int id;
178
    double r,phi;
179
    char line[400];
180
    std::ifstream f(fname);
181
    if (f.is_open()){
182
      std::cout << "Reading ... " <<  fname <<  std::endl;
183
      f.getline(line,400);
184
      while (f.good()){
185
	f >> id >> r >> phi;
186
	if (id<4) module[id].SetRPhi(r,phi);
187
	if (m_debug) std::cout << fname << " " << id << " " << r <<" "<< phi << std::endl;
188
      }
189
      f.close();
190
    } else {
191
       std::cout << "Cannot read " <<  fname <<  std::endl;
192
    }
193
 
194
  };
195
 
196
  int getvalue(xmlXPathContextPtr xpathCtx, const char* xpathExpr, std::vector<xmlChar *> &retval) {
197
 
198
    xmlXPathObjectPtr xpathObj;
199
    xmlNodeSetPtr nodes;
200
    xmlChar * ret=NULL;
201
    int size;
202
    int i;
203
    assert(xpathExpr);
204
    /* Evaluate xpath expression */
205
    xpathObj = xmlXPathEvalExpression(BAD_CAST xpathExpr, xpathCtx);
206
    if(xpathObj == NULL) {
207
      fprintf(stderr,"Error: unable to evaluate xpath expression \"%s\"\n", xpathExpr);
208
      return retval.size();
209
    }
210
 
211
    nodes  = xpathObj->nodesetval;
212
    size   = (nodes) ? nodes->nodeNr : 0;
213
    for(i = 0; i< size; i++) {
214
      ret = xmlNodeGetContent(nodes->nodeTab[i]);
215
      if(ret) {
216
        retval.push_back(ret);
217
        printf("[%d] %s         Return value: %s\n", i, xpathExpr, ret);
218
      }
219
    }
220
 
221
    /* Cleanup of XPath data */
222
    xmlXPathFreeObject(xpathObj);
223
 
224
    return retval.size();
225
  }
226
 
227
  char *  getvalue(xmlXPathContextPtr xpathCtx, const char* xpathExpr, int which) {
228
 
229
    xmlXPathObjectPtr xpathObj;
230
    xmlNodeSetPtr nodes;
231
    xmlChar * ret=NULL;
232
    int size;
233
    int i;
234
    assert(xpathExpr);
235
 
236
    xpathObj = xmlXPathEvalExpression(BAD_CAST xpathExpr, xpathCtx);
237
    if(xpathObj == NULL) {
238
      fprintf(stderr,"Error: unable to evaluate xpath expression \"%s\"\n", xpathExpr);
239
      return NULL;
240
    }
241
 
242
    nodes  = xpathObj->nodesetval;
243
    size   = (nodes) ? nodes->nodeNr : 0;
244
    for(i = 0; i< size; i++) {
245
      ret = xmlNodeGetContent(nodes->nodeTab[i]);
246
      if(ret) {
247
        printf("%s=%s\n",  xpathExpr, ret);
248
        if (which == i) break;
249
      }
250
    }
251
 
252
    xmlXPathFreeObject(xpathObj);
253
    return (char *) ret;
254
 
255
  }
256
 
257
  int readxmlconfig(const char *fname){
258
    xmlInitParser();
259
    LIBXML_TEST_VERSION
260
 
261
    /* Load XML document */
262
    xmlDocPtr doc = xmlParseFile(fname);
263
    if (doc == NULL) {
264
      fprintf(stderr, "Error: unable to parse file \"%s\"\n", fname);
265
      return(-1);
266
    } else {
267
      std::cout << "Reading ... " <<  fname <<  std::endl;
268
    }
269
 
270
    /* Create xpath evaluation context */
271
    xmlXPathContextPtr xpathCtx = xmlXPathNewContext(doc);
272
    if(xpathCtx == NULL) {
273
      fprintf(stderr,"Error: unable to create new XPath context\n");
274
      xmlFreeDoc(doc);
275
      return(-1);
276
    }
277
 
278
    m_pedestalsName       =  getvalue(xpathCtx, "//pedestals"   , 0);
279
    //m_sumpedestalsName    =  getvalue(xpathCtx, "//sumpedestals", 0);
280
    m_photopeakName       =  getvalue(xpathCtx, "//photopeak"   , 0);
281
    m_modulesmapName      =  getvalue(xpathCtx, "//modules"     , 0);
282
    m_channelsmapName     =  getvalue(xpathCtx, "//channels"    , 0);
283
    m_calibrationrootName =  getvalue(xpathCtx, "//channelcalibration", 0);
284
    m_threshold           =  atoi(getvalue(xpathCtx, "//adcthreshold", 0));
285
    m_nx                   = atoi(getvalue(xpathCtx, "//nofcrystalsx", 0));
286
    m_ny                   = atoi(getvalue(xpathCtx, "//nofcrystalsx", 0));
287
    m_dx                   = atof(getvalue(xpathCtx, "//crystalpitchx", 0));
288
    m_dy                   = atof(getvalue(xpathCtx, "//crystalpitchy", 0));
289
 
290
    xmlXPathFreeContext(xpathCtx);
291
    xmlFreeDoc(doc);
292
    xmlCleanupParser();
293
    xmlMemoryDump();
294
    return 0;
295
  };
296
 
297
 
298
  double GetEnergy(int ch, int adc){
299
   return (adc-gPedestals[ch])/(gPeak[ch]-gPedestals[ch])*m_peakscaling;
300
  }
301
  int GetRealPosition(int ipmt, float px,float py,float &cx,float &cy){
302
      // iz CoG izracuna pozicijo na kristalu
303
 
304
    int binx = m_crystalid[ipmt]->GetXaxis()->FindBin(px);
305
    int biny = m_crystalid[ipmt]->GetYaxis()->FindBin(py);
306
    int crystalid  = m_crystalid[ipmt]->GetBinContent(binx,biny);
307
    cx= ( crystalid % m_nx - m_nx/2. + 0.5 ) * m_dx;
308
    cy= ( crystalid / m_nx - m_ny/2. + 0.5 ) * m_dy;
309
 
310
 
311
    // razmazi pozicijo po povrsini kristala
312
    double rndx = (m_rnd->Rndm()-0.5) * m_dx;
313
    double rndy = (m_rnd->Rndm()-0.5) * m_dy;
314
    cx+= rndx;
315
    cy+= rndy;
316
    return crystalid;
317
  }
318
 
319
  TVector3 GetGlobalPosition(int ipmt, double angle, float cx, float cy){
320
 
321
      double phi = angle + module[ipmt].phi;
322
      double &r   = module[ipmt].r;
323
      double sinphi = sin(phi);
324
      double cosphi = cos(phi);
325
 
326
      TVector3 rv( cosphi* r + sinphi * cx , sinphi* r - cosphi * cx, cy);
327
      //if (m_debug) printf(  "[%d] phi=%f %f global ( %f , %f , %f )\n" ,ipmt,  module[ipmt].phi,angle, rv.x(), rv.y(), rv.z());
328
      return  rv;
329
  }
330
 
331
  Geometry(const char *fnameconfig, int debug){
332
    m_debug = debug;
333
    readxmlconfig(fnameconfig);
334
 
335
    ReadModuleMap(m_modulesmapName);
336
    ReadChannelMap(m_channelsmapName);
337
 
338
    std::cout << "Reading ... " <<  m_calibrationrootName <<  std::endl;
339
    m_calibration = new TFile(m_calibrationrootName);
340
    for (int i=0; i<4;i++) {
341
      char hn[256];
342
      sprintf(hn,"pmt1%d_calib",i);
343
      m_crystalid[i]= (TH1I *) m_calibration->Get(hn);
344
      m_crystalid[i]->ls();
345
    }
346
 
347
    //readfile(m_sumpedestalsName ,gSumPedestals,4);
348
    m_peakscaling=3000;
349
 
350
    readfile(m_pedestalsName,gPedestals,4*16, 0);
351
    readfile(m_photopeakName,gPeak,4*16, m_peakscaling);
352
 
353
    m_rnd= new TRandom3();
354
  };
355
 
356
};
357
 
358
 
359
class readdata {
360
private:
361
 
362
  TNtuple* gNtuple;
363
 
364
 
365
  TH1F* m_Adc[MAXCH];
366
  TH1F* m_AdcCut[MAXCH];
367
  TH2F* m_Adc_vs_Nhits[MAXCH];
368
  TH2F* m_Adc_vs_Sum[MAXCH];
369
  TH2F* m_Adc_vs_Sum_Uncorected[MAXCH];
370
  TH2F* m_Nhits;
371
  TH1F* m_AdcSum[6];
372
  TH2F* m_CenterOfGravity[6];
373
  TH2F* m_ReconstructedPosition[6];
374
  TH2F* m_CenterOfGravityforChAboveThreshold[6];
375
  TH2F* m_GlobalPosition;
376
  TH1F* m_AdcSumCluster[6];
377
  TH2F* m_MaxAdc[4];
378
  TH2F* m_SumAdc[4];
379
 
380
  Geometry *m_geo;    // pozicije modulov in kanalov
381
  float gSum[6];
382
  float gRawSum[6];
383
  float gMax[6];
384
  float gSumCluster[6];
385
  float gNtdata[MAXCH*3];
386
  int   gNabove[5];
387
  int m_neve; // number of events to process
388
 
389
public:
390
  PETProjDataMgr *m_projector;  // projektor za sinograme
391
 
392
  double m_rotation;   // trenutna rotacija setupa
393
  double m_threshold;  //
394
 
395
  double gData[MAXCH]; // korigirani ADC ji
396
  double gAdc[MAXCH];  // raw ADC ji
397
  int m_debug;         // debug izpisi
398
  int m_write;         // ce je ta flag nastavljen, se v root file izpisuje ntuple
399
  int m_coincidences;  // stevilo koincidenc
400
 
401
  //---------------------------------------------------
402
  int Initialize(){
403
 
404
    if (m_write) {
405
      char varlist[1024], tmplist[1024];
406
      sprintf( varlist,"sum0:sum1:sum2:sum3:nh0:nh1:nh2:nh3:c0:c1:c2:c3:px0:px1:px2:px3:py0:py1:py2:py3:rx0:rx1:rx2:rx3:ry0:ry1:ry2:ry3:max0:max1:max2:max3");
407
      for (int i=0;i<MAXCH;i++) {
408
        sprintf(tmplist,"%s",varlist);
409
        sprintf(varlist,"%s:a%d",tmplist,i);
410
      }
411
    /*
412
      energija .. (adc -pedestal) * scaling
413
 
414
      sum0 .. sum3 vsota energij po vseh kanalih fotopomnozevalke
415
      nh0  .. nh4   stevilo zadetkov na fotopomnozevalki (adc nad thresholdom)
416
      c0   .. c4 vsota energij za kanale nad thresholdom
417
      px0  .. px4 x koordinata COG na fotopomnozevalki
418
      py0  .. py4 y koordinata COG na fotopomnozevalki
419
      max0 .. max3 maximalna energija na kanalu na fotopomnozevalki
420
      a0   .. a63 energija na i-tem kanalu
421
    */
422
      gNtuple = new TNtuple("nt","Pet RAW data",varlist);
423
      printf ("#%s#\n",varlist);
424
    }
425
    m_Nhits = new TH2F("hnabove","Stevilo zadetkov nad nivojem diskriminacije",4,-0.5,3.5, 16,+0.5,16.5 );  // stevilo hitov nad thresholdom
426
    m_GlobalPosition = new TH2F("globalxy","Rekonstruirana koordinata zadetka",200,-80,80,200,-80,80);  // reconstruirana koordinata v globalnem sistemu
427
 
428
 
429
    char hname[0xFF];
430
    char hn[0xFF];
431
    for (int i=0;i<MAXCH;i++){
432
      sprintf(hname,"Raw ADC Ch. %d ;ADC;N",i);
433
      sprintf(hn,"ach%d",i);
434
      m_Adc[i]   = new TH1F(hn,hname,4000,-0.5,3999.5);       // osnovni adcji
435
      sprintf(hname,"cADC za dogodke z manj kot x hitov na pmtju Ch. %d ;ADC;N",i);
436
      sprintf(hn,"cutch%d",i);
437
      m_AdcCut[i]   = new TH1F(hn,hname,4000,-0.5,3999.5);    // adcji za zadetke z manj kot 7 hiti na PMTju
438
      sprintf(hname,"cADC proti stevilu zadetkov na pmtju Ch. %d ;ADC;N",i);
439
      sprintf(hn,"singlech%d",i);
440
      m_Adc_vs_Nhits[i]   = new TH2F(hn,hname,200,0,4000, 17,-0.5,16.5); // adc proti stevilu zadetkov na PMTju
441
      sprintf(hname,"cADC proti vsoti kanalov na pmtju,  Ch. %d ;ADC;ADCsum",i);
442
      sprintf(hn,"corrch%d",i);
443
      m_Adc_vs_Sum[i]   = new TH2F(hn,hname,200,0,4000, 200,0,12000 );  // raw adc proti vsoti kanalov na PMTju
444
      sprintf(hname,"Raw ADC proti vsoti kanalov na pmtju,  Ch. %d ;ADC;ADCsum",i);
445
      sprintf(hn,"adcvssumch%d",i);
446
      m_Adc_vs_Sum_Uncorected[i]   = new TH2F(hn,hname,100,0,3500, 100,5000,12000 );  // adc proti vsoti kanalov na PMTju
447
    }
448
 
449
    for (int i=0;i<4;i++) {
450
      sprintf(hname,"Vsota cADC na PMTju %d ;ADC;N",i);
451
      sprintf(hn,"pmt%d",i);
452
      m_AdcSum[i]   = new TH1F(hn,hname,200,0,20000); // vsota adcjev na PMTju
453
      sprintf(hname,"Vsota cADC na PMTju za dogodke z manj kot x zadetki %d ;ADC;N",i);
454
      sprintf(hn,"clusterpmt%d",i);
455
      m_AdcSumCluster[i]   = new TH1F(hn,hname,200,0,20000);  // vsota adcjev na pmtju, ki so nad threshholdom,  za dogodke z manj kot 5 zadetki na PMTju
456
 
457
      sprintf(hname,"Center naboja CoG PMT %d ;x;y",i);
458
      sprintf(hn,"pmt1%d",i);
459
      m_CenterOfGravity[i]   =  new TH2F(hn,hname,200,-0.25,3.25,200,-0.25,3.25); // center naboja na PMTju
460
 
461
      sprintf(hname,"Center naboja CoG za zadetke nbad thresholdom PMT %d ;x;y",i);
462
      sprintf(hn,"pmt2%d",i);
463
      m_CenterOfGravityforChAboveThreshold[i]   = new TH2F(hn,hname,200,0,3,200,0,3);  // center naboja za zadetke nad thresholdom
464
 
465
      sprintf(hname,"Rekonstruirana koordinata PMT %d ;x(mm);y(mm)",i);
466
      sprintf(hn,"pmt3%d",i);
467
      m_ReconstructedPosition[i]   = new TH2F(hn,hname,200,-12,12,200,-12,12); // korigirana koordinata na pmtju
468
 
469
      sprintf(hname,"SumAdc PMT %d ;SumAdc(a.u.);crystal ID",i);
470
      sprintf(hn,"sumadc%d",i);
471
      m_SumAdc[i]   = new TH2F(hn,hname,200,0,12000,81,-0.5,80.5); // SumADC spekter za razlicne kristale
472
    }
473
 
474
    // store mapping
475
    TH1I *hx= new TH1I("mapping", "ch mapping", 16, -0.5,15.5);
476
    for (int j=0;j<16;j++) hx->Fill(j,m_geo->ch[j].idx);
477
 
478
    return 0;
479
  };
480
 
481
  int FillHistograms(){
482
	for (int i=0;i<MAXCH;i++) {  // zanka cez vse elektronske kanale;
483
	  float x=gData[i];
484
	  int ipmt= i/16;
485
	  float s=gSum[ipmt];
486
	  // int idx= 2*(ipmt)+64;
487
	  // const  float fSlope1=2;
488
          // const  float fSlope2=3;
489
	  //  if ( ((s-gSumPedestals[ipmt] ) > fSlope1*(x-gPedestals[i])) && ((s-gSumPedestals[ipmt] ) < fSlope2*(x-gPedestals[ipmt])) ) {
490
	  if (gNabove[ipmt]<7){
491
            m_Adc_vs_Sum[i]->Fill(x,s);
492
            m_Adc_vs_Sum_Uncorected[i]->Fill(gAdc[i],gRawSum[ipmt]);
493
            m_AdcCut[i]->Fill(x);
494
	  }
495
	}
496
 
497
	TVector3 r_coincidence[2];
498
	int  f_coincidence[2]={0,0};
499
 
500
	for (int ipmt=0;ipmt<4;ipmt++){ // zanka preko pmtjev
501
	  int j2= (ipmt/2)*2+1-ipmt%2;  // sosednja fotopomnozevalka
502
 
503
          float posx[2]={0,0};
504
          float posy[2]={0,0};
505
          float sum[2]={0,0};
506
          for (int ich=0;ich<16;ich++){  // zanka preko elektronskih kanalov na fotopomnozevalki
507
            int ch= ich+ipmt*16;
508
            if (gMax[ipmt]>m_threshold) {
509
	      posx[0]+= gData[ch]*m_geo->ch[ich].ix;
510
	      posy[0]+= gData[ch]*m_geo->ch[ich].iy;
511
	      sum[0] += gData[ch];
512
 
513
	      if (gData[ch]> 0.2*gMax[ipmt]){  // pri racunanju pozicije upostevaj le kanale, ki imajo vrednost vecjo od ratio*maksimalna na tisti fotopomnozevalki
514
		posx[1]+= gData[ch]*m_geo->ch[ich].ix;
515
		posy[1]+= gData[ch]*m_geo->ch[ich].iy;
516
		sum[1] += gData[ch];
517
	      }
518
	    }
519
	  }
520
 
521
          if (m_write) {
522
            gNtdata[12+ipmt]=posx[0]/sum[0];
523
            gNtdata[16+ipmt]=posy[0]/sum[0];
524
            gNtdata[20+ipmt]=(sum[1]>0)?posx[1]/sum[1]:-10;
525
            gNtdata[24+ipmt]=(sum[1]>0)?posy[1]/sum[1]:-10;
526
            gNtdata[28+ipmt]=sum[0];
527
          }
528
	  if (gMax[ipmt]>gMax[j2]){ // izberi fotopomnozevalko, ki ima vec zadetkov kot soseda TRG= (P1|P2) & (P3|P4)
529
 
530
	    if ( sum[0] > 0  ) {
531
	      float px=posx[0]/sum[0];
532
	      float py=posy[0]/sum[0];
533
 
534
	      m_CenterOfGravity[ipmt]->Fill(px,py);
535
 
536
	      float cx=0; //koordinata na pmtju
537
	      float cy=0;
538
	      int crystalID = m_geo->GetRealPosition(ipmt,px,py,cx,cy);
539
              m_ReconstructedPosition[ipmt]->Fill(cx,cy);
540
              m_SumAdc[ipmt]->Fill(gSum[ipmt], crystalID );
541
              if (m_debug > 2) printf("ipmt=%d cx=%f cy=%f\n",ipmt,cx,cy);
542
              r_coincidence[ipmt/2] =  m_geo->GetGlobalPosition(ipmt, m_rotation, cx, cy) ;
543
              f_coincidence[ipmt/2] =1;
544
              m_GlobalPosition->Fill( r_coincidence[ipmt/2].x(),r_coincidence[ipmt/2].y());
545
 
546
	    }
547
	    if ( sum[1] > 0  ) {
548
 
549
	      float px=posx[1]/sum[1];
550
	      float py=posy[1]/sum[1];
551
              m_CenterOfGravityforChAboveThreshold[ipmt]->Fill(px,py);
552
	    }
553
 
554
	  }
555
	}
556
 
557
	for (int i=0;i<MAXCH;i++) if (gData[i]>m_threshold) { m_Adc_vs_Nhits[i]->Fill(gData[i],gNabove[i/16]); }
558
	for (int i=0;i<MAXCH;i++) if (gData[i]>m_threshold && gNabove[i/16]<5 ) { gSumCluster[i/16]+=gData[i]; }
559
	for (int ipmt=0;ipmt<4;ipmt++) {
560
	  m_AdcSum[ipmt]->Fill(gSum[ipmt]);
561
	  m_AdcSumCluster[ipmt]->Fill(gSumCluster[ipmt]);
562
	}
563
 
564
	for (int i=0;i<4;i++) m_Nhits->Fill(i,gNabove[i]);
565
 
566
	if (m_write) {
567
	  for (int i=0;i<4;i++) gNtdata[i]=gSum[i];
568
	  for (int i=0;i<4;i++) gNtdata[4+i]=gNabove[i];
569
	  for (int i=0;i<4;i++) gNtdata[8+i]=gSumCluster[i];
570
	  for (int i=0;i<MAXCH;i++) gNtdata[32+i]=gData[i];
571
	  gNtuple->Fill(gNtdata);
572
	}
573
	//
574
	// coincidences
575
	//
576
	if (f_coincidence[0]&&f_coincidence[1]){
577
          m_coincidences++;
578
          if (m_debug > 1) {
579
            printf(  "Coincidence ( %f , %f , %f ) ( %f , %f , %f )\n" ,
580
            r_coincidence[0].x(), r_coincidence[0].y(), r_coincidence[0].z(),
581
            r_coincidence[1].x(), r_coincidence[1].y(), r_coincidence[1].z() );
582
          }
583
          m_projector->AddEvent( r_coincidence[0] , r_coincidence[1]);
584
        }
585
    return 0;
586
  };
587
 
588
 
589
int Process( const char *fnamein, int nevetoprocess){
590
    gzFile fp=gzopen(fnamein,"r");
591
    printf("Process....\n");
592
    if (!fp) {
593
      printf("File %s cannot be opened!\n", fnamein);
594
      return 0;
595
    } else {
596
      printf("Processing file %s\n", fnamein);
597
    }
598
    const int bsize=20000;
599
 
600
    unsigned int *buf = new unsigned int[bsize];
601
    int nr=0;
602
    int hdr[4];
603
 
604
    int nread=0;
605
    while (!gzeof(fp) ){
606
      int nitems  = gzread(fp,hdr,   sizeof(int)*4 );
607
      if (nitems !=4*sizeof(int)) break;
608
      int recid   = hdr[0];
609
      int nb      = hdr[1] - 4*sizeof(int);
610
      if ( nb > bsize) {
611
	nb=bsize;
612
	printf ("Error bsize %d nb=%d\n",bsize,nb);
613
      }
614
      // int nint = nb/sizeof(int);
615
      int n=hdr[3];
616
      nitems=gzread(fp, buf,  nb);
617
      nr++;
618
      if (nitems!=nb){  printf("Read Error nb=%d nitems=%d\n",nb,nitems); continue; }
619
 
620
      switch (recid){
621
      case RUNREC_ID:
622
	runrec=(RUNREC *) (&buf[0]);
623
	printf("RUNREC RECID        = %u\n",hdr[0]);
624
	printf("RUNREC Length       = %u\n",hdr[1]);
625
	printf("RUNREC File version = %u\n",hdr[2]);
626
 
627
	printf("RUNREC Time = %u\n",hdr[3]);
628
	printf("Number of events per step   = %u\n",runrec->num_events);
629
	printf("Number of channels measured = %u\n",runrec->num_channels);
630
	printf("Pedestal              = %u\n",runrec->pedestal);
631
	printf("Scan type = %u :: 0 -> single data scan :: 1 -> XY position scan\n",runrec->xy);
632
	printf("Number of steps in  X = %d\n",runrec->nx);
633
	printf("Start position      X = %d\n",runrec->x0);
634
	printf("Step size direction X = %d\n",runrec->dx);
635
	printf("Number of steps in  Y = %d\n",runrec->ny);
636
	printf("Start position      Y = %d\n",runrec->y0);
637
	printf("Step size direction Y = %d\n",runrec->dy);
638
 
639
	break;
640
      case POSREC_ID:
641
	posrec=(POSREC *) (&buf[0]);
642
	printf("POSREC nx=%d , posx=%d setx=%d posy=%d sety=%d angle(deg)= %3.2f \n",n,posrec->mikro_pos_x,posrec->set_pos_x,posrec->mikro_pos_y,posrec->set_pos_y, posrec->set_pos_x*360./NSTEPS);
643
	m_rotation =  posrec->set_pos_x*2*TMath::Pi()/NSTEPS;
644
	break;
645
      case EVTREC_ID:
646
	break;
647
 
648
      }
649
 
650
      if (recid!=EVTREC_ID)  continue;
651
      if (nread++==nevetoprocess) break;
652
      int idx=1;
653
      int neve=buf[0]/2;
654
      //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
655
      for (int nev=0;nev<neve;nev++){
656
	int len=buf[idx];
657
	int sb =buf[idx+1];
658
	unsigned int *pbuf=&buf[idx+2];
659
	if (sb!=0xffab) {
660
	  printf("0x%04x!0xffab len=%d\n",sb,len);
661
	  break;
662
	}
663
	// postavi na nic
664
#define InitArrayWithAValue(arr,n,x) {for (int i=0;i<n;i++) arr[i]=x;}
665
	InitArrayWithAValue( gSum       , 4    , 0);
666
	InitArrayWithAValue( gRawSum    , 4    , 0);
667
	InitArrayWithAValue( gMax       , 4    , 0);
668
	InitArrayWithAValue( gSumCluster, 4    , 0);
669
	InitArrayWithAValue( gNabove    , 4    , 0);
670
	InitArrayWithAValue( gData      , MAXCH, 0);
671
	InitArrayWithAValue( gAdc       , MAXCH, 0);
672
        //------------------------------------------------------------
673
	for (int i0=0;i0<len-2;i0+=2) {
674
 
675
	  int data0  = pbuf[i0];
676
	  int data1  = pbuf[i0+1];
677
	  int geo    = (data1 >> 11) & 0x1f;
678
	  int ch     = (data1&0x1f)  | (geo<<5);
679
	  int dtype  = (data1>>9)&0x3;
680
	  int adc    =  data0&0xfff;
681
 
682
	  switch (dtype) {
683
          case 0x0:
684
            if (ch<MAXCH) {
685
              m_Adc[ch]->Fill(adc);
686
 
687
              int ipmt = ch/16;
688
 
689
              gAdc[ch]=adc;
690
              gRawSum[ipmt]+=adc;
691
 
692
              if (ch<64) gData[ch]=m_geo->GetEnergy(ch,adc); else gData[ch]=adc;
693
 
694
              gSum[ipmt]+=gData[ch];
695
              if (gData[ch]   >gMax[ipmt] )     gMax[ipmt]= gData[ch];
696
              if (gData[ch]   >m_threshold )    gNabove[ipmt]++;
697
            }
698
            break;
699
          case 0x10:
700
          case 0x11:
701
          case 0x01:
702
            break;
703
	  }
704
 
705
	};// for (int i0=0;i0<len-2;i0+=2)
706
        //------------------------------------------------------------
707
 
708
	idx+=len+1;
709
	FillHistograms();
710
      } // for (int nev=0;nev<neve;nev++)
711
      //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
712
    } //  while (!gzeof(fp) )
713
    gzclose(fp);
714
    delete buf;
715
    return nr;
716
    }
717
 
718
  readdata(const char *fnameconfig, std::vector<TString> &fnamein, const char *fnameout, int write2root=0, int nevetoread=-1, int debug=0 ){
719
    m_debug = debug;
720
    m_write = write2root;
721
    m_neve  = nevetoread;
722
 
723
    m_geo = new Geometry(fnameconfig,m_debug);
724
    m_threshold=m_geo->GetThreshold();
725
 
726
    char rootname[0xFF]; sprintf(rootname,"%s.root",fnameout);
727
    TFile *rootFile = new TFile(rootname,"RECREATE");
728
    Initialize();
729
    printf("Konec inicializacije ...\n");
730
 
731
    m_projector=PETProjDataMgr::GetInstance();
732
    m_projector->SetDebug(m_debug);
733
    m_projector->SetRingDiameter(2*m_geo->module[0].r);
734
    m_projector->SetFilename(fnameout);
735
 
736
    m_coincidences=0;
737
 
738
    int nr = 0;
739
    for (unsigned int i=0;i<  fnamein.size() ;i++) {
740
      int neve = m_neve -nr;
741
      if (m_neve == -1)  nr += Process(fnamein[i].Data(), m_neve);
742
      else if (neve>0)   nr += Process(fnamein[i].Data(), neve);
743
    }
744
    printf("End...\nreaddata::Number of events read=%d\nreaddata::Number of coincidences=%d\n",nr,m_coincidences);
745
 
746
    m_projector->WriteInterfile();
747
    rootFile->Write();
748
    rootFile->Close();
749
  }
750
 
751
  ~readdata(){
752
   if (m_geo) delete m_geo;
753
  };
754
 
755
};
756
 
757
#ifdef MAIN
758
 
759
 
760
int main(int argc, char **argv){
761
 
762
  std::cout << "Usage:" << argv[0] <<  " -c config.xml -i input.dat -o output -n nevents -w writentuple -d debuglevel "<< std::endl;
763
  //-----------------------------------------------------------------
764
  char configfile[256]="ini/config.xml";
765
  char outputfile[256]="/tmp/petfmf";
766
  int nevetoread      =-1;
767
  int writentuple     =0;
768
  int verbose         =0;
769
  char c;
770
 
771
  std::vector<TString> inputfilelist;
772
  while ((c=getopt (argc, argv, "c:i:o:n:w:d:")) != -1){
773
    switch(c){
774
    case 'c': sprintf(configfile,"%s",optarg); break;
775
    case 'i': inputfilelist.push_back(TString(optarg));  break;
776
    case 'o': sprintf(outputfile,"%s",optarg); break;
777
    case 'n': nevetoread = atoi(optarg)    ;   break;
778
    case 'w': writentuple = atoi(optarg);      break;
779
    case 'd': verbose = atoi(optarg);          break;
780
    default: abort();
781
    }
782
  }
783
  std::cout << "Used: " << argv[0] <<  " -c " << configfile << " -o " << outputfile << " -n " << nevetoread << " -w " << writentuple << " -d  " << verbose ;
784
  for (unsigned int i=0;i<  inputfilelist.size() ;i++) std::cout << " -i " << inputfilelist[i];
785
  std::cout << std::endl;
786
  //-----------------------------------------------------------------
787
  new readdata(configfile, inputfilelist,outputfile, writentuple, nevetoread, verbose);
788
 
789
  std::cout << "Output data files " << outputfile <<  std::endl;
790
  return 0;
791
}
792
 
793
#endif