Subversion Repositories f9daq

Rev

Blame | Last modification | View Log | RSS feed

// g++  -I`root-config  --incdir` sa02read.C `root-config  --libs` -lz -o sa02read
#ifdef WIN32
#pragma warning( disable : 4996  )
#endif
#include <stdlib.h>
#include <stdio.h>
#include <stdint.h>
#include <iostream>
#include <fstream>
#include <zlib.h>
#include <iostream>
#include <vector>
#include <string>
#include "TString.h"
#include "TFile.h"
#include "TNtuple.h"
#include "TH3F.h"
#include "TH2F.h"
#include "TH1F.h"
#include "TH3D.h"
#include "TH1D.h"
#include "TGraph.h"
#include "TROOT.h"
#include "TObject.h"

#include "SA02_DEF.h"
#include "dataio.h"
#include "H3D.h"
#include "H2D.h"
#include "H1D.h"
#ifdef WIN32
#  include "Tchar.h"
#  include "XGetopt.h"
# else /* WIN32 */
#  include <getopt.h>
#endif  /* WIN32 */

int gposx=-1;
int gposy=-1;
char gmapfile[256]="";
int gMapped=0;
int gInitialized=0;
unsigned int gEid2HapdXy[144];
int gBoardNumber=0;
int gIsNsteps=1;
int gScanType=0;

MONREC monrec;
RUNREC runrec;
POSREC posrec;
ITERATORREC iteratorrec;
DATREC datrec;
RUNINFO runinfo;
EVTREC evtrec;
H3D h3d;
H2D h2d;
H1D h1d;
int V729_first=1;
//--------------------------------------------------------------------
int ReadMapFile(char *fname) {

  int chip, ch,x,y;
  std::ifstream din;
  din.open (fname);
  for (int i=0; i<144; i++) {
    din >> chip >> ch >> x>> y;
    x = 11 -x;
    y = 11 -y;
    gEid2HapdXy[chip*36+ch]=x*12+y;
  }
  din.close();
  return 1;
}

//--------------------------------------------------------------------
//--------------------------------------------------------------------
class readdata {

private:
  TGraph *m_GrTime;
  TGraph *m_GrData;
  TGraph *m_GrEvent;
  TNtuple* m_Ntuple;
  TNtuple* m_thresh;
  TH1F* m_NHit;
  TH2F* m_Scan;
  TH2F* m_Histo;
  TH2F* m_ScanSet;
  TH3F* m_Scan3d;
  TH2F* m_xy[144];
  TH2F* m_xysum;
  TTree * m_ttmon;
  TNtuple *m_ntpos;
  TTree *m_ntrun;
public:
  int m_neve; // number of events to process
  int m_ny;
  int m_dy;
  int m_y0;
  int m_val;
  int m_time;
  int m_cureve;
  int m_y;
  int m_towrite;
//public:
  int m_mask;          // maska za status register
  int m_shft;          // shift za status register
  int m_debug;         // debug izpisi
  int m_write;         // ce je ta flag nastavljen, se v root file izpisuje ntuple
  unsigned int m_bitmask ;

//--------------------------------------------------------------------
#define MAXCH 144

  int Initialize(int nx,int ny,double ny0,double ny1) {

    printf("Initialize %d %d %2.2f %2.2f\n",nx,ny,ny0,ny1);
    if (m_write) {
      m_Ntuple = new TNtuple("nt","SA02 RAW data","time:eve:x:y");
      m_GrTime = new TGraph(ny);
      m_GrTime->SetTitle("grtime");
      m_GrEvent= new TGraph(ny);
      m_GrEvent->SetTitle("grevent");
      m_GrData = new TGraph(ny);
      m_GrData->SetTitle("grdata");
    }
    m_NHit = new TH1F("NHit","N hit per event",nx+1,-0.5,nx+0.5);
    m_Scan  = new TH2F("Scan","Scan;SA02 ch;Threshold value",nx,-0.5,nx-0.5, ny,ny0,ny1 );
//     m_Scan  = new TH2F("Scan","Scan;SA02 ch;Threshold value",nx,-0.5,nx-0.5, (ny+1.2452)/0.0024,ny0,ny1 );
    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 );
//     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 );
    m_Histo = new TH2F("Histo","Histo;SA02 ch;bit",nx,-0.5,nx-0.5, 8,-0.5,7.5 );
    m_ScanSet  = new TH2F("ScanSet","Scan;SA02 ch;set value",nx,-0.5,nx-0.5, ny,ny0,ny1 );
    m_thresh = new TNtuple("thresh","bo","ch:y:xl:yl");
    return 1;
  }

//--------------------------------------------------------------------
  int FillHistograms(unsigned int *data, double yset, double y) {

    static int trgn = 0;
    int nhit=0;

    for (int ch=0; ch<144; ch++) {
      int id=(17-ch/8)+gBoardNumber*18;
      int shft=(ch%8)*4;
      for (int i=0; i<4; i++) {
        if (data[id] & (1<<(i+shft)) )  {
          if ( m_Scan3d ) m_Scan3d->Fill(ch,i,y);
          if ( m_Histo  ) m_Histo->Fill(ch,i);
        }
      }
      if (data[id] & (m_bitmask<<shft) ) {
        nhit++;
        if (m_Scan) m_Scan->Fill(ch,y);
        if (m_ScanSet) m_ScanSet->Fill(ch,yset);
        unsigned int ii= gEid2HapdXy[ch];
        if (m_xy[ii]) {
          double xs;
          double ys;
          if (runrec.fver) {
            xs = posrec.ix*runrec.dx+runrec.x0;
            ys = posrec.iy*runrec.dy+runrec.y0;
          } else {
            xs = posrec.xset;
            ys = posrec.yset;
          }
          if (m_xysum) m_xysum->Fill( xs, ys);
//!          m_xy[ii]->Fill(xs, ys);
          m_xy[ii]->Fill(posrec.xset, posrec.yset);
//          printf("%3d c=%d ch=%2d x=%6d y=%6d\n",ch, ii/36,ii%36,posrec.set_pos_x, posrec.set_pos_y);
        } else {
          printf("error m_xy[%d] does not exist\n",ii);
        }
        m_thresh->Fill(ch,y,trgn,1);
      }
    }
    if (m_NHit) m_NHit->Fill(nhit);
    trgn++;
    return 1;
  }

//--------------------------------------------------------------------
  int ProcessData(unsigned int *rdata,int count,int neve) {

    int nc=0;
    for (int j=0; j<count;) {
      unsigned int recid   = rdata[j++];
      unsigned int len     = rdata[j++];
      if (m_debug) printf(" recid=0x%0x len=%d pos=%d(maxpos %d) val=%d\n",recid, len, j, count, rdata[j]);
      if (m_debug>2 && len>2) printf("\n");
      if(j+(int)len>=count) return nc;
      switch (recid) {
        case 3:
          if (len !=18) printf(" recid=0x%0x len=%d pos=%d(maxpos %d) val=%d\n",recid, len, j, count, rdata[j]);
          nc+=FillHistograms(&rdata[j], m_y, m_val);
          if (nc==neve) return nc;
          break;
        default:
          m_val=(rdata[j] >> m_shft) & m_mask ;
          if (m_write) {
            m_Ntuple->Fill(m_time,m_cureve,m_y,m_val);
            m_GrTime->SetPoint(m_cureve,m_time,m_val);
            m_GrEvent->SetPoint(m_cureve,m_cureve,m_val);
            m_GrData->SetPoint(m_cureve,m_y,m_val);
          }
          break;
      }
      if (m_debug>2) {
        for (unsigned int i=0; i<len; i++) {
          if (j<count) printf("0x%0x ", rdata[j]);
          if (i%4==3) printf("\n");
          j++;
        }
        printf("\n");
      } else j+=len;
//      if(j+4>=count) return nc;
    }

    return nc;
  }

//--------------------------------------------------------------------
  double H3DGetBinContentFrom( H3D *h ,unsigned int atx, unsigned int aty, unsigned int atz) {

    unsigned int idx;
    if (!h) return 0;
    if (h->nx <= atx)  return 0;
    if (h->ny <= aty)  return 0;
    if (h->nz <= atz)  return 0;
    idx = atx+(aty+atz*h->ny)*h->nx;
    if (idx*sizeof(double)  < h->size ) {
//    if (h->data[idx]>0)  printf("xy %d %d %g\n",atx,aty,h->data[idx]);
      return h->data[idx];
    }
    return 0;
  }

//--------------------------------------------------------------------
  double H2DGetBinContentFrom( H2D *h ,unsigned int atx, unsigned int aty) {

    unsigned int idx;
    if (!h) return 0;
    if (h->nx <= atx)  return 0;
    if (h->ny <= aty)  return 0;
    idx = atx+aty*h->nx;
    if (idx*sizeof(double)  < h->size ) {
//    if (h->data[idx]>0)  printf("xy %d %d %g\n",atx,aty,h->data[idx]);
      return h->data[idx];
    }
    return 0;
  }

//--------------------------------------------------------------------
  double H1DGetBinContentFrom( H1D *h ,unsigned int atx) {

    if (!h) return 0;
    if (h->nx <= atx)  return 0;
    if (atx*sizeof(double)  < h->size ) {
      return h->data[atx];
    }
    return 0;
  }

//--------------------------------------------------------------------
  TH2D *hwf[4];
  int V729_init() {
    printf(" V729_init\n");
    const double maxy=512;
//    const double maxy=512;
    const int ny=100;
//    const int h0=10;
    int len=100;

    TObject *hobj = gROOT->FindObject("gwf0");
    if (hobj !=NULL) hobj->Write();
    hobj = gROOT->FindObject("gwf1");
    if (hobj !=NULL) hobj->Write();
    hobj = gROOT->FindObject("gwf2");
    if (hobj !=NULL) hobj->Write();
    hobj = gROOT->FindObject("gwf3");
    if (hobj !=NULL) hobj->Write();
    hwf[0]= new TH2D("gwf0","CAENV_729 ch0",len,-0.5,len-0.5,ny,-maxy,maxy);
    hwf[1]= new TH2D("gwf1","CAENV_729 ch1",len,-0.5,len-0.5,ny,-maxy,maxy);
    hwf[2]= new TH2D("gwf2","CAENV_729 ch2",len,-0.5,len-0.5,ny,-maxy,maxy);
    hwf[3]= new TH2D("gwf3","CAENV_729 ch3",len,-0.5,len-0.5,ny,-maxy,maxy);
//    h1d=new TH1F("h1d","Charge;signal charge(a.u.);N",500,-maxy,5*maxy);
    return 0;
  }

//--------------------------------------------------------------------
  int V729_plot(uint16_t *buf ) {
    
        if (V729_first){
          V729_init();
      V729_first=0;       
        }
    int start=310;
    int end=335;
//    const double maxy=2048;
//    const double maxy=512;
//    const double maxy=512;
//    const int ny=100;
//    const int h0=10;
    int i;
    uint32_t *hdr = (uint32_t *)buf;
    int nb= hdr[1]-2*sizeof(uint32_t);
    int nitems = nb /sizeof(uint16_t);
    unsigned int recid=hdr[0];
    uint16_t *data = &buf[4];
    double x,y;
//    double baseline=0;
    double charge=0;
    while (nitems>0) {
      int id = data[0]-1;
      int len= data[1]/sizeof(uint16_t)-2;
      if (m_debug)  printf("\t->> recid=%d nb=%d nitems=%d id=%d len=%d\n", recid, hdr[1],nitems,id,len);
      if (id>3) return 0;
      data+=2;
//      for (i=0; i<200;i++) baseline += data[i];
//      baseline /=200.;
      for (i=0; i<len; i++) {
        if (m_debug >2) printf("ADC %d %d\n", i, data[i]);
        x=i;
        y=data[i];
//        y-=baseline;
        if (i>start && i<end) charge += y;

        hwf[id]->Fill(x,y-2000);
      }
      data+=len;
      nitems=nitems - len-2;
//      if (id==2) h1d->Fill(charge);
    }
    return 0;
  }

//--------------------------------------------------------------------
  int Process( const char *fnamein, int nevetoprocess, int histoonly) {

    int nxh,nyh;
    gzFile fp=gzopen(fnamein,"r");
    printf("Process....\n");
    if (!fp) {
      printf("File %s cannot be opened!\n", fnamein);
      return 0;
    } else {
      printf("Processing file %s\n", fnamein);
    }
    unsigned int bsize=1000000; // initial buffer size
    unsigned int *buf = new unsigned int[bsize];
    unsigned int *data= buf+4; // data shift
    int nr=0;
    HEADERREC *hdr = (HEADERREC *) buf;
//    int nread=0;
    while (!gzeof(fp) ) {
      int nitems  = gzread(fp,buf,sizeof(HEADERREC));
      if (nitems !=sizeof(HEADERREC)) break;
      int recid   = hdr->id;
      int nb      = hdr->len - sizeof(HEADERREC);
      if (hdr->len > bsize*sizeof(int)) {
        printf ("[INFO] Requested data size=%u bigger than buffer size %u\n",hdr->len/sizeof(uint32_t), bsize);
        bsize=2*hdr->len/sizeof(int);
        printf("[INFO] Increasing data buffer to %d elements\n", bsize);
        delete buf;
        buf = new unsigned int[bsize];
        buf[0] = recid;
        buf[1] = nb + sizeof(HEADERREC);
        hdr = (HEADERREC *) buf;
        data = buf +4;
      }
      nitems=gzread(fp,((unsigned char *)buf)+sizeof(HEADERREC),nb);
      nr++;
      if (m_debug || nr<30) printf("[RECID=%u] Length=%u items=%u\n",hdr->id,hdr->len,hdr->len/sizeof(uint32_t));
      if (nitems!=nb) {
        printf("Read Error nb=%d nitems=%d. Exiting ...\n",nb,nitems);
        break;
//        exit(0);
      }
      switch (recid) {
        case POSREC_ID: {
          memcpy(&posrec.id, buf, sizeof(posrec));
          printf("POSRECID = %u pos= %d %d\n",posrec.id, posrec.ix,posrec.iy);
          //if (!histoonly)  V729_init();
          if (m_debug || nr<20) {
            printf("Length = %u\n",posrec.len);
            printf("Time = %u\n",posrec.time);
            printf("Iteration X = %d\n",posrec.ix);
            printf("MIKRO Position X = %d\n",posrec.x);
            printf("Set position X = %d\n",posrec.xset);
            printf("Iteration Y = %d\n",posrec.iy);
            printf("MIKRO Position Y = %d\n",posrec.y);
            printf("Set position Y = %d\n",posrec.yset);
          }
                m_ntpos->Fill(posrec.time-runrec.time, posrec.ix,posrec.x,posrec.xset, posrec.iy,posrec.y,posrec.yset);
        }
        break;
                
                
                case MONREC_ID : {
                        memcpy(&monrec.id, buf, sizeof(monrec));

                        //m_ttmon->Branch("monrec",&monrec,"id/i:len:time:imon[6]/I:vmon[6]:status[6]");
                        m_ttmon->Branch("id",&monrec.id,"id/i");
                        m_ttmon->Branch("len",&monrec.len,"len/i");
                        m_ttmon->Branch("time",&monrec.time,"time/i");
                        m_ttmon->Branch("imon",&monrec.imon[0],"imon[24]/I");
                        m_ttmon->Branch("iset",&monrec.iset[0],"iset[24]/I");
                        m_ttmon->Branch("vmon",&monrec.vmon[0],"vmon[24]/I");
                        m_ttmon->Branch("vset",&monrec.vset[0],"vset[24]/I");
                        m_ttmon->Branch("status",&monrec.status[0],"status[24]/I");
                        m_ttmon->Fill();  

                        break;
                }
                
        case RUNREC_ID: {
                        
                  memcpy(&runrec.id, buf, sizeof(runrec));
                        m_ntrun->Branch("time",&runrec.time,"time/i");
                        m_ntrun->Branch("fver",&runrec.fver,"fver/i");
                        m_ntrun->Branch("nev",&runrec.nev,"nev/i");
                        m_ntrun->Branch("ped",&runrec.ped,"ped/i");
                        m_ntrun->Branch("nx",&runrec.nx,"nx/I");
                        m_ntrun->Branch("x0",&runrec.x0,"x0/I");
                        m_ntrun->Branch("dx",&runrec.dx,"dx/I");
                        m_ntrun->Branch("ny",&runrec.ny,"ny/I");
                        m_ntrun->Branch("y0",&runrec.y0,"y0/I");
                        m_ntrun->Branch("dy",&runrec.dy,"dy/I");
                        m_ntrun->Fill();
                                
          printf("***************RECID = %u\n",runrec.id);
          printf("Length = %u\n",runrec.len);
          printf("File version = %u\n",runrec.fver);
          printf("Time = %u\n",runrec.time);
          printf("Number of events per step = %u\n",runrec.nev);
          printf("Pedestal = %u\n",runrec.ped);
//          printf("Scan type = %u :: 0 -> single data scan :: 1 -> XY position scan\n",scanrec.xy);
          printf("Number of steps in X = %d\n",runrec.nx);
          printf("Start position X = %d\n",runrec.x0);
          printf("Step size direction X = %d\n",runrec.dx);
          printf("Number of steps in Y = %d\n",runrec.ny);
          printf("Start position Y = %d\n",runrec.y0);
          printf("Step size direction Y = %d\n",runrec.dy);
          if (gIsNsteps) {
            nxh=runrec.nx;
            nyh=runrec.ny;
          } else {
            nxh=(runrec.nx-1)/runrec.dx+1;
            nyh=(runrec.ny-1)/runrec.dy+1;
          }  
          for (int i=0;i<144;i++) {
            char hname[256];
            if (gMapped) sprintf(hname,"ch_%02d_%02d",i%12,i/12);
            else sprintf(hname,"xy_%02d_%02d",i/36,i%36);
            m_xy[i]=new TH2F(hname,hname,nxh,runrec.x0-runrec.dx*0.5,runrec.x0+runrec.dx*(nxh-0.5),
                                         nyh,runrec.y0-runrec.dy*0.5,runrec.y0+runrec.dy*(nyh-0.5));
          }
          m_xysum= (TH2F *) (m_xy[0]->Clone("xy"));
          m_xysum->SetTitle("xy all channels");
        }
        break;
        case RUNINFO_ID: {
          memcpy(&runinfo.id, buf, sizeof(runinfo));
          if (!gInitialized) {
            m_y0  = runinfo.x0;
            m_dy  = runinfo.dx;
            m_ny  = runinfo.nx;
            m_towrite = runinfo.writemode;
            if (m_towrite ==3) m_y0=runinfo.chip*36+runinfo.ch;
            printf("RUNINFO x0=%d nx=%d dx=%d\n", runinfo.x0,runinfo.dx,runinfo.nx);
            gInitialized=Initialize(144,m_ny,m_y0-0.5,m_y0-0.5+m_dy*m_ny);
          }
        }
        break;
        case EVTREC_ID:
          memcpy(&evtrec.id, buf, sizeof(evtrec));
          m_time   = evtrec.time;
          m_cureve = evtrec.nev;
          m_y = m_y0 +m_cureve * m_dy;
//          nread++;
          break;
//        case ITERATORREC_ID:
//          memcpy(&iteratorrec.id, buf, sizeof(iteratorrec));
//          printf("ITERATORREC level=%d n=%d value=%f step=%f\n", iteratorrec.level,iteratorrec.n,iteratorrec.value, iteratorrec.step);
//          break;
        case DATREC_ID:
          memcpy(&datrec.id, buf, sizeof(datrec));
          printf("DATREC chip=%d channel=%d cmd=%d data=%d retval=%d\n", datrec.chip,datrec.channel,datrec.cmd, datrec.data, datrec.retval);
          break;
        case CAENV729REC_ID: {
          if (!histoonly) V729_plot( (uint16_t*)buf );
        }
        break;
        case H3D_ID: {
          int hlen= sizeof(h3d)-sizeof(double *);
          memcpy(&h3d.id, buf, hlen);
          h3d.data = new double[h3d.size/sizeof(double)];
//          int align=0;
//          if (sizeof(double *)==4) align=0;
//          memcpy(h3d.data, ((char *) buf)+hlen-align, h3d.size);
          memcpy(h2d.data, ((char *) buf)+hlen, h2d.size);
          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 );
//          printf("H3D sizeof(H2D)=%lu len-datasize=%d len=%lu datasize=%d\t",hlen,h3d.len-h3d.size,h3d.len,h3d.size);
//          printf("H3D sz=%d\n",sizeof(double *));
          H3D  *h = &h3d;
          char hname[100];
          sprintf(hname,"%s_%d",h->name,datrec.data);
          TObject *hobj = gROOT->FindObject(hname);
          if (hobj !=NULL) {
            printf("TH3D duplicate name %s!\n",hname);
            hobj->Write();
            delete hobj;
//            break;
          }
          if (h->stepx==0) h->stepx=1;
          if (h->stepy==0) h->stepy=1;
          if (h->stepz==0) h->stepz=1;
          TH3D *h3 = new TH3D(hname, h->title,
                              h->nx, h->minx-0.5*h->stepx,h->minx+(h->nx-0.5)*h->stepx,
                              h->ny, h->miny-0.5*h->stepy,h->miny+(h->ny-0.5)*h->stepy,
                              h->nz, h->minz-0.5*h->stepz,h->minz+(h->nz-0.5)*h->stepz);
          h3->SetTitle(h->title);
          h3->GetXaxis()->SetTitle(h->titlex);
          h3->GetYaxis()->SetTitle(h->titley);
          h3->GetZaxis()->SetTitle(h->titlez);
          if (m_debug) {
            printf("TH3D name=%s %s %s %s %s",hname, h->title,h->titlex, h->titley, h->titlez);
            printf("\tnx=%d %f %f \tny=%d %f %f\tnz=%d %f %f\n",
                   h->nx, h->minx-0.5*h->stepx,h->minx+(h->nx-0.5)*h->stepx,
                   h->ny, h->miny-0.5*h->stepy,h->miny+(h->ny-0.5)*h->stepy,
                   h->nz, h->minz-0.5*h->stepz,h->minz+(h->nz-0.5)*h->stepz);
          }
          for (unsigned int ix=0; ix<h->nx; ix++) {
            for (unsigned int iy=0; iy<h->ny; iy++) {
              for (unsigned int iz=0; iz<h->nz; iz++) {
                double val = H3DGetBinContentFrom(h,ix,iy,iz);
                if (m_debug>2) printf("%d %d %d %g\n", ix,iy,iz,val);
                h3->SetBinContent(ix+1,iy+1,iz+1,val);
              }
            }
          }
          if (m_debug) {
            h3->ls();
            h3->Print();
            printf("TH3D Sum %f\n", h3->GetMaximum());
          }
        }
        break;
        case H2D_ID: {
          if(gScanType==1) break;
          int hlen= sizeof(h2d)-sizeof(double *);
          memcpy(&h2d.id, buf, hlen);
          h2d.data = new double[h2d.size/sizeof(double)];
//          int align=0;
//          if (sizeof(double *)==4) align=0;
//          memcpy(h2d.data, ((char *) buf)+hlen-align, h2d.size);
          memcpy(h2d.data, ((char *) buf)+hlen, h2d.size);
          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 );
//          printf("H2D sizeof(H2D)=%lu len-datasize=%d len=%lu datasize=%d\t",hlen,h2d.len-h2d.size,h2d.len,h2d.size);
//          printf("H2D sz=%d\n",sizeof(double *));
          H2D  *h = &h2d;
          char hname[100];
          sprintf(hname,"%s_%d",h->name,datrec.data);
          TObject *hobj = gROOT->FindObject(hname);
          if (hobj !=NULL) {
            printf("TH2D duplicate name %s!\n",hname);
            hobj->Write();
            delete hobj;
//            break;
          }
          if (h->stepx==0) h->stepx=1;
          if (h->stepy==0) h->stepy=1;
          TH2D *h2 = new TH2D(hname, h->title,
                              h->nx, h->minx-0.5*h->stepx,h->minx+(h->nx-0.5)*h->stepx,
                              h->ny, h->miny-0.5*h->stepy,h->miny+(h->ny-0.5)*h->stepy);
          h2->SetTitle(h->title);
          h2->GetXaxis()->SetTitle(h->titlex);
          h2->GetYaxis()->SetTitle(h->titley);
          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,
                                h->nx, h->minx-0.5*h->stepx,h->minx+(h->nx-0.5)*h->stepx,
                                h->ny, h->miny-0.5*h->stepy,h->miny+(h->ny-0.5)*h->stepy);
          for (unsigned int ix=0; ix<h->nx; ix++) {
            for (unsigned int iy=0; iy<h->ny; iy++) {
              double val = H2DGetBinContentFrom(h,ix,iy);
              if (m_debug>2) printf("%d %d %g\n", ix,iy,val);
              h2->SetBinContent(ix+1,iy+1,val);
            }
          }
          if (m_debug) {
            h2->ls();
            h2->Print();
            printf("TH2D Sum %f\n", h2->GetMaximum());
          }
        }
        break;
        case H1D_ID: {
          int hlen= sizeof(h1d)-sizeof(double *);
          memcpy(&h1d.id, buf, hlen);
          h1d.data = new double[h1d.size/sizeof(double)];
//          int align=0;
//          if (sizeof(double *)==4) align=0;
//          memcpy(h1d.data, ((char *) buf)+hlen-align, h1d.size);
          memcpy(h1d.data, ((char *) buf)+hlen, h1d.size);
          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 );
          H1D  *h = &h1d;
          char hname[100];
          sprintf(hname,"%s_%d",h->name,datrec.data);
          TObject *hobj = gROOT->FindObject(hname);
          if (hobj !=NULL) {
            printf("TH1D duplicate name %s!\n",hname);
            hobj->Write();
            delete hobj;
//            break;
          }
          if (h->stepx==0) h->stepx=1;
          TH1D *h1 = new TH1D(hname, h->title, h->nx, h->minx-0.5*h->stepx,h->minx+(h->nx-0.5)*h->stepx);
          h1->SetTitle(h->title);
          h1->GetXaxis()->SetTitle(h->titlex);
          h1->GetYaxis()->SetTitle(h->titley);
//          h1->GetXaxis()->SetTitle("timebin");
//          h1->GetYaxis()->SetTitle(h->titlex);
          if (m_debug) printf("TH1D name=%s %s %s %s\tnx=%d %g %g\n",hname, h->title,h->titlex, h->titley,
                                h->nx, h->minx-0.5*h->stepx,h->minx+(h->nx-0.5)*h->stepx);
          for (unsigned int ix=0; ix<h->nx; ix++) {
            double val = H1DGetBinContentFrom(h,ix);
            if (m_debug>2) printf("%d %g\n", ix,val);
            h1->SetBinContent(ix+1,val);
          }
          if (m_debug) {
            h1->ls();
            h1->Print();
            printf("TH1D Sum %f\n", h1->GetMaximum());
          }
        }
        break;
        default:
          printf("Unknown record 0x%x Exiting .... \n", recid);
                  gzclose(fp);
          delete buf;
          return nr;
          
      }
      if (recid!=EVTREC_ID)  continue;
      if (histoonly) continue;
      if (m_debug) printf("EVTREC %d histo=%d\n",recid,histoonly);
      if (gposx>=0 && gposx!=posrec.ix ) continue;
      if (gposy>=0 && gposy!=posrec.iy ) continue;
      int nc=ProcessData(data,nb/sizeof(unsigned int)-2,nevetoprocess);
      if (m_debug) printf("EVTREC Events %d nb=%u\n",nc,nb/sizeof(unsigned int)-2);
//      if (nread++==nevetoprocess) break;
    } //  while (!gzeof(fp))
    gzclose(fp);
    delete buf;
    return nr;
  }

//--------------------------------------------------------------------
  readdata( std::vector<TString> &fnamein,const char *fnameout,int write2root=0,
            int nevetoread=-1,unsigned int bmask=0xff,int mask=0xFFFFFFFF,int shift=0,
            int neve=-1,int debug=0,int histoonly=0) {

    m_debug = debug;
    m_write = write2root;
    m_neve  = nevetoread;
    m_mask  = mask;
    m_shft =  shift;
    m_towrite = 0;
    m_bitmask = bmask;
    char rootname[0xFF];
    sprintf(rootname,"%s",fnameout);
    TFile *rootFile = new TFile(rootname,"RECREATE");
        m_ntpos = new TNtuple("pos","Position record", "time:ix:x:xset:iy:y:yset");
        m_ntrun = new TTree("run","Run record");
        m_ttmon = new TTree("Monitor","Voltage, Current and Status Monitor");

//    Initialize();
    printf("Konec inicializacije ...\n");
    if (m_debug) printf("readdata histo=%d\n",histoonly);
    int nr = 0;
    for (unsigned int i=0; i<  fnamein.size() ; i++) {
      int neve = m_neve -nr;
      if (m_neve == -1)  nr += Process(fnamein[i].Data(), m_neve, histoonly);
      else if (neve>0)   nr += Process(fnamein[i].Data(), neve, histoonly);
    }
    printf("End...\nreaddata::Number of events read=%d required=%d\n",nr,m_neve);
    if (m_write) {
      m_GrTime->Write();
      m_GrData->Write();
      m_GrEvent->Write();
    }
    rootFile->Write();
    if (m_debug) rootFile->ls();
    rootFile->Close();
  }

//--------------------------------------------------------------------
  ~readdata() { };
};

//--------------------------------------------------------------------
//-----------------------------------------------------------------
#ifdef MAIN

int main(int argc, char **argv) {

  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;
  char outputfile[256]="/tmp/sa02asic.root";
  int nevetoread      =-1;
  int writentuple     =0;
  int verbose         =0;
  int shift=0;
  int histoonly=0;
  Int_t mask =0xFFFFFFFF;
  unsigned int bitmask=0xf;
  char c;
  for (int i=0; i<144; i++) gEid2HapdXy[i]=i;
  std::vector<TString> inputfilelist;
  while ((c=getopt (argc, argv, "l:i:I:o:n:w:d:m:s:x:y:f:b:h:t:")) != -1) {
    switch (c) {
      case 'I': // runrec.nx = xmax-xmin - for magnet tests in folder 2014
        gIsNsteps=0;
      case 'i': // runrec.nx == x steps - standard
        inputfilelist.push_back(TString(optarg));
        break;
      case 'o':
        sprintf(outputfile,"%s",optarg);
        break;
      case 'n':
        nevetoread = atoi(optarg)    ;
        break;
      case 'w':
        writentuple = atoi(optarg);
        break;
      case 'd':
        verbose = atoi(optarg);
        break;
      case 'l':
        gBoardNumber = atoi(optarg);
        break;
      case 's':
        shift = strtoul(optarg,NULL,0);
        break;
      case 'm':
        mask  = strtoul(optarg,NULL,0);
        break;
      case 'b':
        bitmask  = strtoul(optarg,NULL,0);
        break;
      case 'x':
        gposx = atoi(optarg)    ;
        break;
      case 'y':
        gposy = atoi(optarg)    ;
        break;
      case 'h':
        histoonly = atoi(optarg)    ;
        break;
      case 'f':
        sprintf(gmapfile,"%s", optarg)    ;
        gMapped=ReadMapFile(gmapfile);
        break;
      case 't':
        gScanType = atoi(optarg);
        break;
      default:
        abort();
    }
  }
  std::cout << "Used: " << argv[0]  << " -o " << outputfile
            << " -n " << nevetoread << " -w " << writentuple
            << " -d  " << verbose  << " -m  0x" << std::hex << mask
            << " -s  0x" << std::hex << shift  << " -l  " << gBoardNumber << " -h  " << histoonly ;
  for (unsigned int i=0; i<  inputfilelist.size() ; i++) std::cout << " -i " << inputfilelist[i];
  std::cout << std::endl;
  if (!inputfilelist.size()) {
    std::cout << "No input file specified" << std::endl;
  } else {
    new readdata(inputfilelist,outputfile,writentuple,nevetoread,bitmask,mask,shift,-1,verbose,histoonly);
    std::cout << "Output data files " << outputfile <<  std::endl;
  }
  return 0;
}
#endif /* MAIN */