Subversion Repositories f9daq

Rev

Blame | Last modification | View Log | RSS feed

  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 */
  785.