Subversion Repositories f9daq

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