Subversion Repositories f9daq

/*

vhodne datoteke:

Mappingi:

m16.map channel mapping

modules.map mapping modulov

Kalibracija:

pedestals.dat  ... adc suma

ph511.dat      ... adc fotovrhov vrhov

// sumpedestals.dat vsota adc po modulih

./readdata -c config.xml -i input.dat -o output -n nevents -w writentuple -d debuglevel

*/

#include <stdlib.h>

#include <stdio.h>

#include "TNtuple.h"

#include "TH1F.h"

#include "TH2F.h"

#include "TF1.h"

#include "TCanvas.h"

#include "TStyle.h"

#include "TFile.h"

#include "TMath.h"

#include "TRandom3.h"

#include <zlib.h>

#include <iostream>

#include <fstream>

#include <stdlib.h>

#include <stdio.h>

#include <string.h>

#include <assert.h>

#include <vector>

#include <libxml/tree.h>

#include <libxml/parser.h>

#include <libxml/xpath.h>

#include <libxml/xpathInternals.h>

#include "PETProjDataMgr.h"

// en krog je 4790 korakov

#define NSTEPS 4790

// nt->Draw("px1/max1:py1/max1>>hmy(200,-1,4,200,-1,4)","max1>500","colz")

//#define MAXCH 64

#define MAXCH 72

#define RUNREC_ID                       1

#define EVTREC_ID                       2

#define POSREC_ID                       3

typedef struct {

  unsigned int num_events,num_channels,pedestal,xy;

  int nx,x0,dx,ny,y0,dy;

} RUNREC;

RUNREC *runrec;  // start header: appears once at the start of the file

typedef struct {

  int mikro_pos_x,set_pos_x;

  int num_iter_y,mikro_pos_y,set_pos_y;

} POSREC;

POSREC *posrec;  // position header: appears at every change of position

class Channel {

public:

  Channel(){};

  ~Channel(){};

  int ix;

  int iy;

  int idx;

  void SetIxIy(int a,int b){ix=a;iy=b;idx=ix+4*iy;};

};

class Module {

public:

  Module(){};

  ~Module(){};

  double r;

  double phi;

  void SetRPhi(double rr, double pphi){r=rr;phi=pphi*TMath::Pi()/180.;};

};

class Geometry {

public:

  ~Geometry(){ m_calibration->Close(); };

  Channel ch[16];

  Module  module[4];

  TRandom3 *m_rnd;   // random generator

  TH1I     *m_crystalid[4];// kalibracijski histogrami

  char * m_modulesmapName;  // ime datoteke s podatki o pozicijah modulov

  char * m_channelsmapName; // ime datoteke s podatki o pozicijah elektronskih kanalov

    //char * m_sumpedestalsName;

  char * m_pedestalsName;   // ime datoteke s podatki o pedestalih

  char * m_photopeakName;   // ime datoteke s podatki o vrhovih

  char * m_calibrationrootName;   // ime datoteke s kalibracijskimi histogrami

  TFile *m_calibration; //  datoteke s kalibracijskimi histogrami

  double m_dx;        // pitch x kristalov

  double m_dy;        // pitch y kristalov

  int    m_nx;        // stevilo kristalov v x smeri

  int    m_ny;        // stevilo kristalov v y smeri

  float gPedestals[MAXCH];

  float gPeak[MAXCH];

  //float gSumPedestals[MAXCH];

  double m_peakscaling;   // scaling za adcje

  double m_threshold;    // threshold za upostevanje zadetkov na kanalih

  int    m_debug;

  //---------------------------------------------------

  int SetThreshold(int threshold){

    m_threshold=threshold;

    return 0;

  };

  double GetThreshold(){

    return m_threshold;

  }

  int readfile(const char *fname, float *x, int nmax, int defaultvalue=0){

    int id;

    float ix;

    std::ifstream f(fname);

    if (f.is_open()){

      std::cout << "Reading ... " <<  fname <<  std::endl;

      while (f.good()){

	f >> id >> ix;

	if (id<nmax)  x[id]=ix;

	if (m_debug) std::cout << fname << " " << id << " " << ix <<" "<< std::endl;

      }

      f.close();

    } else {

      std::cout << "Cannot read " <<  fname <<  std::endl;

      std::cout << "Default value will be used :" << defaultvalue  <<  std::endl;

      for (int i=0;i<nmax;i++){

        x[i]=defaultvalue;

      }

    }

    return 0;

  };

  void ReadChannelMap(const char *fname){

    int id,ix,iy;

    char line[400];

    std::ifstream f(fname);

    if (f.is_open()){

      std::cout << "Reading ... " <<  fname <<  std::endl;

      f.getline(line,400);

      while (f.good()){

	f >> id >> ix >> iy;

	if (id<16) ch[id].SetIxIy(ix,iy);

	if (m_debug) std::cout << fname << " " << id << " " << ix <<" "<< iy << std::endl;

      }

      f.close();

    } else {

       std::cout << "Cannot read " <<  fname <<  std::endl;

    }

  };

  void ReadModuleMap(const char *fname){

    int id;

    double r,phi;

    char line[400];

    std::ifstream f(fname);

    if (f.is_open()){

      std::cout << "Reading ... " <<  fname <<  std::endl;

      f.getline(line,400);

      while (f.good()){

	f >> id >> r >> phi;

	if (id<4) module[id].SetRPhi(r,phi);

	if (m_debug) std::cout << fname << " " << id << " " << r <<" "<< phi << std::endl;

      }

      f.close();

    } else {

       std::cout << "Cannot read " <<  fname <<  std::endl;

    }

  };

  int getvalue(xmlXPathContextPtr xpathCtx, const char* xpathExpr, std::vector<xmlChar *> &retval) {

    xmlXPathObjectPtr xpathObj;

    xmlNodeSetPtr nodes;

    xmlChar * ret=NULL;

    int size;

    int i;

    assert(xpathExpr);

    /* Evaluate xpath expression */

    xpathObj = xmlXPathEvalExpression(BAD_CAST xpathExpr, xpathCtx);

    if(xpathObj == NULL) {

      fprintf(stderr,"Error: unable to evaluate xpath expression \"%s\"\n", xpathExpr);

      return retval.size();

    }

    nodes  = xpathObj->nodesetval;

    size   = (nodes) ? nodes->nodeNr : 0;

    for(i = 0; i< size; i++) {

      ret = xmlNodeGetContent(nodes->nodeTab[i]);

      if(ret) {

        retval.push_back(ret);

        printf("[%d] %s         Return value: %s\n", i, xpathExpr, ret);

      }

    }

    /* Cleanup of XPath data */

    xmlXPathFreeObject(xpathObj);

    return retval.size();

  }

  char *  getvalue(xmlXPathContextPtr xpathCtx, const char* xpathExpr, int which) {

    xmlXPathObjectPtr xpathObj;

    xmlNodeSetPtr nodes;

    xmlChar * ret=NULL;

    int size;

    int i;

    assert(xpathExpr);

    xpathObj = xmlXPathEvalExpression(BAD_CAST xpathExpr, xpathCtx);

    if(xpathObj == NULL) {

      fprintf(stderr,"Error: unable to evaluate xpath expression \"%s\"\n", xpathExpr);

      return NULL;

    }

    nodes  = xpathObj->nodesetval;

    size   = (nodes) ? nodes->nodeNr : 0;

    for(i = 0; i< size; i++) {

      ret = xmlNodeGetContent(nodes->nodeTab[i]);

      if(ret) {

        printf("%s=%s\n",  xpathExpr, ret);

        if (which == i) break;

      }

    }

    xmlXPathFreeObject(xpathObj);

    return (char *) ret;

  }

  int readxmlconfig(const char *fname){

    xmlInitParser();

    LIBXML_TEST_VERSION

    /* Load XML document */

    xmlDocPtr doc = xmlParseFile(fname);

    if (doc == NULL) {

      fprintf(stderr, "Error: unable to parse file \"%s\"\n", fname);

      return(-1);

    } else {

      std::cout << "Reading ... " <<  fname <<  std::endl;

    }

    /* Create xpath evaluation context */

    xmlXPathContextPtr xpathCtx = xmlXPathNewContext(doc);

    if(xpathCtx == NULL) {

      fprintf(stderr,"Error: unable to create new XPath context\n");

      xmlFreeDoc(doc);

      return(-1);

    }

    m_pedestalsName       =  getvalue(xpathCtx, "//pedestals"   , 0);

    //m_sumpedestalsName    =  getvalue(xpathCtx, "//sumpedestals", 0);

    m_photopeakName       =  getvalue(xpathCtx, "//photopeak"   , 0);

    m_modulesmapName      =  getvalue(xpathCtx, "//modules"     , 0);

    m_channelsmapName     =  getvalue(xpathCtx, "//channels"    , 0);

    m_calibrationrootName =  getvalue(xpathCtx, "//channelcalibration", 0);

    m_threshold           =  atoi(getvalue(xpathCtx, "//adcthreshold", 0));

    m_nx                   = atoi(getvalue(xpathCtx, "//nofcrystalsx", 0));

    m_ny                   = atoi(getvalue(xpathCtx, "//nofcrystalsx", 0));

    m_dx                   = atof(getvalue(xpathCtx, "//crystalpitchx", 0));

    m_dy                   = atof(getvalue(xpathCtx, "//crystalpitchy", 0));

    xmlXPathFreeContext(xpathCtx);

    xmlFreeDoc(doc);

    xmlCleanupParser();

    xmlMemoryDump();

    return 0;

  };

  double GetEnergy(int ch, int adc){

   return (adc-gPedestals[ch])/(gPeak[ch]-gPedestals[ch])*m_peakscaling;

  }

  int GetRealPosition(int ipmt, float px,float py,float &cx,float &cy){

      // iz CoG izracuna pozicijo na kristalu

    int binx = m_crystalid[ipmt]->GetXaxis()->FindBin(px);

    int biny = m_crystalid[ipmt]->GetYaxis()->FindBin(py);

    int crystalid  = m_crystalid[ipmt]->GetBinContent(binx,biny);

    cx= ( crystalid % m_nx - m_nx/2. + 0.5 ) * m_dx;

    cy= ( crystalid / m_nx - m_ny/2. + 0.5 ) * m_dy;

    // razmazi pozicijo po povrsini kristala

    double rndx = (m_rnd->Rndm()-0.5) * m_dx;

    double rndy = (m_rnd->Rndm()-0.5) * m_dy;

    cx+= rndx;

    cy+= rndy;

    return crystalid;

  }

  TVector3 GetGlobalPosition(int ipmt, double angle, float cx, float cy){

      double phi = angle + module[ipmt].phi;

      double &r   = module[ipmt].r;

      double sinphi = sin(phi);

      double cosphi = cos(phi);

      TVector3 rv( cosphi* r + sinphi * cx , sinphi* r - cosphi * cx, cy);

      //if (m_debug) printf(  "[%d] phi=%f %f global ( %f , %f , %f )\n" ,ipmt,  module[ipmt].phi,angle, rv.x(), rv.y(), rv.z());

      return  rv;

  }

  Geometry(const char *fnameconfig, int debug){

    m_debug = debug;

    readxmlconfig(fnameconfig);

    ReadModuleMap(m_modulesmapName);

    ReadChannelMap(m_channelsmapName);

    std::cout << "Reading ... " <<  m_calibrationrootName <<  std::endl;

    m_calibration = new TFile(m_calibrationrootName);

    for (int i=0; i<4;i++) {

      char hn[256];

      sprintf(hn,"pmt1%d_calib",i);

      m_crystalid[i]= (TH1I *) m_calibration->Get(hn);

      m_crystalid[i]->ls();

    }

    //readfile(m_sumpedestalsName ,gSumPedestals,4);

    m_peakscaling=3000;

    readfile(m_pedestalsName,gPedestals,4*16, 0);

    readfile(m_photopeakName,gPeak,4*16, m_peakscaling);

    m_rnd= new TRandom3();

  };

};

class readdata {

private:

  TNtuple* gNtuple;

  TH1F* m_Adc[MAXCH];

  TH1F* m_AdcCut[MAXCH];

  TH2F* m_Adc_vs_Nhits[MAXCH];

  TH2F* m_Adc_vs_Sum[MAXCH];

  TH2F* m_Adc_vs_Sum_Uncorected[MAXCH];

  TH2F* m_Nhits;

  TH1F* m_AdcSum[6];

  TH2F* m_CenterOfGravity[6];

  TH2F* m_ReconstructedPosition[6];

  TH2F* m_CenterOfGravityforChAboveThreshold[6];

  TH2F* m_GlobalPosition;

  TH1F* m_AdcSumCluster[6];

  TH2F* m_MaxAdc[4];

  TH2F* m_SumAdc[4];

  Geometry *m_geo;    // pozicije modulov in kanalov

  float gSum[6];

  float gRawSum[6];

  float gMax[6];

  float gSumCluster[6];

  float gNtdata[MAXCH*3];

  int   gNabove[5];

  int m_neve; // number of events to process

public:

  PETProjDataMgr *m_projector;  // projektor za sinograme

  double m_rotation;   // trenutna rotacija setupa

  double m_threshold;  //

  double gData[MAXCH]; // korigirani ADC ji

  double gAdc[MAXCH];  // raw ADC ji

  int m_debug;         // debug izpisi

  int m_write;         // ce je ta flag nastavljen, se v root file izpisuje ntuple

  int m_coincidences;  // stevilo koincidenc

  //---------------------------------------------------

  int Initialize(){

    if (m_write) {

      char varlist[1024], tmplist[1024];

      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");

      for (int i=0;i<MAXCH;i++) {

        sprintf(tmplist,"%s",varlist);

        sprintf(varlist,"%s:a%d",tmplist,i);

      }

    /*

      energija .. (adc -pedestal) * scaling

      sum0 .. sum3 vsota energij po vseh kanalih fotopomnozevalke

      nh0  .. nh4   stevilo zadetkov na fotopomnozevalki (adc nad thresholdom)

      c0   .. c4 vsota energij za kanale nad thresholdom

      px0  .. px4 x koordinata COG na fotopomnozevalki

      py0  .. py4 y koordinata COG na fotopomnozevalki

      max0 .. max3 maximalna energija na kanalu na fotopomnozevalki

      a0   .. a63 energija na i-tem kanalu

    */

      gNtuple = new TNtuple("nt","Pet RAW data",varlist);

      printf ("#%s#\n",varlist);

    }

    m_Nhits = new TH2F("hnabove","Stevilo zadetkov nad nivojem diskriminacije",4,-0.5,3.5, 16,+0.5,16.5 );  // stevilo hitov nad thresholdom

    m_GlobalPosition = new TH2F("globalxy","Rekonstruirana koordinata zadetka",200,-80,80,200,-80,80);  // reconstruirana koordinata v globalnem sistemu

    char hname[0xFF];

    char hn[0xFF];

    for (int i=0;i<MAXCH;i++){

      sprintf(hname,"Raw ADC Ch. %d ;ADC;N",i);

      sprintf(hn,"ach%d",i);

      m_Adc[i]   = new TH1F(hn,hname,4000,-0.5,3999.5);       // osnovni adcji

      sprintf(hname,"cADC za dogodke z manj kot x hitov na pmtju Ch. %d ;ADC;N",i);

      sprintf(hn,"cutch%d",i);

      m_AdcCut[i]   = new TH1F(hn,hname,4000,-0.5,3999.5);    // adcji za zadetke z manj kot 7 hiti na PMTju

      sprintf(hname,"cADC proti stevilu zadetkov na pmtju Ch. %d ;ADC;N",i);

      sprintf(hn,"singlech%d",i);

      m_Adc_vs_Nhits[i]   = new TH2F(hn,hname,200,0,4000, 17,-0.5,16.5); // adc proti stevilu zadetkov na PMTju

      sprintf(hname,"cADC proti vsoti kanalov na pmtju,  Ch. %d ;ADC;ADCsum",i);

      sprintf(hn,"corrch%d",i);

      m_Adc_vs_Sum[i]   = new TH2F(hn,hname,200,0,4000, 200,0,12000 );  // raw adc proti vsoti kanalov na PMTju

      sprintf(hname,"Raw ADC proti vsoti kanalov na pmtju,  Ch. %d ;ADC;ADCsum",i);

      sprintf(hn,"adcvssumch%d",i);

      m_Adc_vs_Sum_Uncorected[i]   = new TH2F(hn,hname,100,0,3500, 100,5000,12000 );  // adc proti vsoti kanalov na PMTju

    }

    for (int i=0;i<4;i++) {

      sprintf(hname,"Vsota cADC na PMTju %d ;ADC;N",i);

      sprintf(hn,"pmt%d",i);

      m_AdcSum[i]   = new TH1F(hn,hname,200,0,20000); // vsota adcjev na PMTju

      sprintf(hname,"Vsota cADC na PMTju za dogodke z manj kot x zadetki %d ;ADC;N",i);

      sprintf(hn,"clusterpmt%d",i);

      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

      sprintf(hname,"Center naboja CoG PMT %d ;x;y",i);

      sprintf(hn,"pmt1%d",i);

      m_CenterOfGravity[i]   =  new TH2F(hn,hname,200,-0.25,3.25,200,-0.25,3.25); // center naboja na PMTju

      sprintf(hname,"Center naboja CoG za zadetke nbad thresholdom PMT %d ;x;y",i);

      sprintf(hn,"pmt2%d",i);

      m_CenterOfGravityforChAboveThreshold[i]   = new TH2F(hn,hname,200,0,3,200,0,3);  // center naboja za zadetke nad thresholdom

      sprintf(hname,"Rekonstruirana koordinata PMT %d ;x(mm);y(mm)",i);

      sprintf(hn,"pmt3%d",i);

      m_ReconstructedPosition[i]   = new TH2F(hn,hname,200,-12,12,200,-12,12); // korigirana koordinata na pmtju

      sprintf(hname,"SumAdc PMT %d ;SumAdc(a.u.);crystal ID",i);

      sprintf(hn,"sumadc%d",i);

      m_SumAdc[i]   = new TH2F(hn,hname,200,0,12000,81,-0.5,80.5); // SumADC spekter za razlicne kristale

    }

    // store mapping

    TH1I *hx= new TH1I("mapping", "ch mapping", 16, -0.5,15.5);

    for (int j=0;j<16;j++) hx->Fill(j,m_geo->ch[j].idx);

    return 0;

  };

  int FillHistograms(){

	for (int i=0;i<MAXCH;i++) {  // zanka cez vse elektronske kanale;

	  float x=gData[i];

	  int ipmt= i/16;

	  float s=gSum[ipmt];

	  // int idx= 2*(ipmt)+64;

	  // const  float fSlope1=2;

          // const  float fSlope2=3;

	  //  if ( ((s-gSumPedestals[ipmt] ) > fSlope1*(x-gPedestals[i])) && ((s-gSumPedestals[ipmt] ) < fSlope2*(x-gPedestals[ipmt])) ) {

	  if (gNabove[ipmt]<7){

            m_Adc_vs_Sum[i]->Fill(x,s);

            m_Adc_vs_Sum_Uncorected[i]->Fill(gAdc[i],gRawSum[ipmt]);

            m_AdcCut[i]->Fill(x);

	  }

	}

	TVector3 r_coincidence[2];

	int  f_coincidence[2]={0,0};

	for (int ipmt=0;ipmt<4;ipmt++){ // zanka preko pmtjev

	  int j2= (ipmt/2)*2+1-ipmt%2;  // sosednja fotopomnozevalka

          float posx[2]={0,0};

          float posy[2]={0,0};

          float sum[2]={0,0};

          for (int ich=0;ich<16;ich++){  // zanka preko elektronskih kanalov na fotopomnozevalki

            int ch= ich+ipmt*16;

            if (gMax[ipmt]>m_threshold) {

	      posx[0]+= gData[ch]*m_geo->ch[ich].ix;

	      posy[0]+= gData[ch]*m_geo->ch[ich].iy;

	      sum[0] += gData[ch];

	      if (gData[ch]> 0.2*gMax[ipmt]){  // pri racunanju pozicije upostevaj le kanale, ki imajo vrednost vecjo od ratio*maksimalna na tisti fotopomnozevalki

		posx[1]+= gData[ch]*m_geo->ch[ich].ix;

		posy[1]+= gData[ch]*m_geo->ch[ich].iy;

		sum[1] += gData[ch];

	      }

	    }

	  }

          if (m_write) {

            gNtdata[12+ipmt]=posx[0]/sum[0];

            gNtdata[16+ipmt]=posy[0]/sum[0];

            gNtdata[20+ipmt]=(sum[1]>0)?posx[1]/sum[1]:-10;

            gNtdata[24+ipmt]=(sum[1]>0)?posy[1]/sum[1]:-10;

            gNtdata[28+ipmt]=sum[0];

          }

	  if (gMax[ipmt]>gMax[j2]){ // izberi fotopomnozevalko, ki ima vec zadetkov kot soseda TRG= (P1|P2) & (P3|P4)

	    if ( sum[0] > 0  ) {

	      float px=posx[0]/sum[0];

	      float py=posy[0]/sum[0];

	      m_CenterOfGravity[ipmt]->Fill(px,py);

	      float cx=0; //koordinata na pmtju

	      float cy=0;

	      int crystalID = m_geo->GetRealPosition(ipmt,px,py,cx,cy);

              m_ReconstructedPosition[ipmt]->Fill(cx,cy);

              m_SumAdc[ipmt]->Fill(gSum[ipmt], crystalID );

              if (m_debug > 2) printf("ipmt=%d cx=%f cy=%f\n",ipmt,cx,cy);

              r_coincidence[ipmt/2] =  m_geo->GetGlobalPosition(ipmt, m_rotation, cx, cy) ;

              f_coincidence[ipmt/2] =1;

              m_GlobalPosition->Fill( r_coincidence[ipmt/2].x(),r_coincidence[ipmt/2].y());

	    }

	    if ( sum[1] > 0  ) {

	      float px=posx[1]/sum[1];

	      float py=posy[1]/sum[1];

              m_CenterOfGravityforChAboveThreshold[ipmt]->Fill(px,py);

	    }

	  }

	}

	for (int i=0;i<MAXCH;i++) if (gData[i]>m_threshold) { m_Adc_vs_Nhits[i]->Fill(gData[i],gNabove[i/16]); }

	for (int i=0;i<MAXCH;i++) if (gData[i]>m_threshold && gNabove[i/16]<5 ) { gSumCluster[i/16]+=gData[i]; }

	for (int ipmt=0;ipmt<4;ipmt++) {

	  m_AdcSum[ipmt]->Fill(gSum[ipmt]);

	  m_AdcSumCluster[ipmt]->Fill(gSumCluster[ipmt]);

	}

	for (int i=0;i<4;i++) m_Nhits->Fill(i,gNabove[i]);

	if (m_write) {

	  for (int i=0;i<4;i++) gNtdata[i]=gSum[i];

	  for (int i=0;i<4;i++) gNtdata[4+i]=gNabove[i];

	  for (int i=0;i<4;i++) gNtdata[8+i]=gSumCluster[i];

	  for (int i=0;i<MAXCH;i++) gNtdata[32+i]=gData[i];

	  gNtuple->Fill(gNtdata);

	}

	//

	// coincidences

	//

	if (f_coincidence[0]&&f_coincidence[1]){

          m_coincidences++;

          if (m_debug > 1) {

            printf(  "Coincidence ( %f , %f , %f ) ( %f , %f , %f )\n" ,

            r_coincidence[0].x(), r_coincidence[0].y(), r_coincidence[0].z(),

            r_coincidence[1].x(), r_coincidence[1].y(), r_coincidence[1].z() );

          }

          m_projector->AddEvent( r_coincidence[0] , r_coincidence[1]);

        }

    return 0;

  };

int Process( const char *fnamein, int nevetoprocess){

    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);

    }

    const int bsize=20000;

    unsigned int *buf = new unsigned int[bsize];

    int nr=0;

    int hdr[4];

    int nread=0;

    while (!gzeof(fp) ){

      int nitems  = gzread(fp,hdr,   sizeof(int)*4 );

      if (nitems !=4*sizeof(int)) break;

      int recid   = hdr[0];

      int nb      = hdr[1] - 4*sizeof(int);

      if ( nb > bsize) {

	nb=bsize;

	printf ("Error bsize %d nb=%d\n",bsize,nb);

      }

      // int nint = nb/sizeof(int);

      int n=hdr[3];

      nitems=gzread(fp, buf,  nb);

      nr++;

      if (nitems!=nb){  printf("Read Error nb=%d nitems=%d\n",nb,nitems); continue; }

      switch (recid){

      case RUNREC_ID:

	runrec=(RUNREC *) (&buf[0]);

	printf("RUNREC RECID        = %u\n",hdr[0]);

	printf("RUNREC Length       = %u\n",hdr[1]);

	printf("RUNREC File version = %u\n",hdr[2]);

	printf("RUNREC Time = %u\n",hdr[3]);

	printf("Number of events per step   = %u\n",runrec->num_events);

	printf("Number of channels measured = %u\n",runrec->num_channels);

	printf("Pedestal              = %u\n",runrec->pedestal);

	printf("Scan type = %u :: 0 -> single data scan :: 1 -> XY position scan\n",runrec->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);

	break;

      case POSREC_ID:

	posrec=(POSREC *) (&buf[0]);

	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);

	m_rotation =  posrec->set_pos_x*2*TMath::Pi()/NSTEPS;

	break;

      case EVTREC_ID:

	break;

      }

      if (recid!=EVTREC_ID)  continue;

      if (nread++==nevetoprocess) break;

      int idx=1;

      int neve=buf[0]/2;

      //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

      for (int nev=0;nev<neve;nev++){

	int len=buf[idx];

	int sb =buf[idx+1];

	unsigned int *pbuf=&buf[idx+2];

	if (sb!=0xffab) {

	  printf("0x%04x!0xffab len=%d\n",sb,len);

	  break;

	}

	// postavi na nic

#define InitArrayWithAValue(arr,n,x) {for (int i=0;i<n;i++) arr[i]=x;}

	InitArrayWithAValue( gSum       , 4    , 0);

	InitArrayWithAValue( gRawSum    , 4    , 0);

	InitArrayWithAValue( gMax       , 4    , 0);

	InitArrayWithAValue( gSumCluster, 4    , 0);

	InitArrayWithAValue( gNabove    , 4    , 0);

	InitArrayWithAValue( gData      , MAXCH, 0);

	InitArrayWithAValue( gAdc       , MAXCH, 0);

        //------------------------------------------------------------

	for (int i0=0;i0<len-2;i0+=2) {

	  int data0  = pbuf[i0];

	  int data1  = pbuf[i0+1];

	  int geo    = (data1 >> 11) & 0x1f;

	  int ch     = (data1&0x1f)  | (geo<<5);

	  int dtype  = (data1>>9)&0x3;

	  int adc    =  data0&0xfff;

	  switch (dtype) {

          case 0x0:

            if (ch<MAXCH) {

              m_Adc[ch]->Fill(adc);

              int ipmt = ch/16;

              gAdc[ch]=adc;

              gRawSum[ipmt]+=adc;

              if (ch<64) gData[ch]=m_geo->GetEnergy(ch,adc); else gData[ch]=adc;

              gSum[ipmt]+=gData[ch];

              if (gData[ch]   >gMax[ipmt] )     gMax[ipmt]= gData[ch];

              if (gData[ch]   >m_threshold )    gNabove[ipmt]++;

            }

            break;

          case 0x10:

          case 0x11:

          case 0x01:

            break;

	  }

	};// for (int i0=0;i0<len-2;i0+=2)

        //------------------------------------------------------------

	idx+=len+1;

	FillHistograms();

      } // for (int nev=0;nev<neve;nev++)

      //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

    } //  while (!gzeof(fp) )

    gzclose(fp);

    delete buf;

    return nr;

    }

  readdata(const char *fnameconfig, std::vector<TString> &fnamein, const char *fnameout, int write2root=0, int nevetoread=-1, int debug=0 ){

    m_debug = debug;

    m_write = write2root;

    m_neve  = nevetoread;

    m_geo = new Geometry(fnameconfig,m_debug);

    m_threshold=m_geo->GetThreshold();

    char rootname[0xFF]; sprintf(rootname,"%s.root",fnameout);

    TFile *rootFile = new TFile(rootname,"RECREATE");

    Initialize();

    printf("Konec inicializacije ...\n");

    m_projector=PETProjDataMgr::GetInstance();

    m_projector->SetDebug(m_debug);

    m_projector->SetRingDiameter(2*m_geo->module[0].r);

    m_projector->SetFilename(fnameout);

    m_coincidences=0;

    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);

      else if (neve>0)   nr += Process(fnamein[i].Data(), neve);

    }

    printf("End...\nreaddata::Number of events read=%d\nreaddata::Number of coincidences=%d\n",nr,m_coincidences);

    m_projector->WriteInterfile();

    rootFile->Write();

    rootFile->Close();

  }

  ~readdata(){

   if (m_geo) delete m_geo;

  };

};

#ifdef MAIN

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

  std::cout << "Usage:" << argv[0] <<  " -c config.xml -i input.dat -o output -n nevents -w writentuple -d debuglevel "<< std::endl;

  //-----------------------------------------------------------------

  char configfile[256]="ini/config.xml";

  char outputfile[256]="/tmp/petfmf";

  int nevetoread      =-1;

  int writentuple     =0;

  int verbose         =0;

  char c;

  std::vector<TString> inputfilelist;

  while ((c=getopt (argc, argv, "c:i:o:n:w:d:")) != -1){

    switch(c){

    case 'c': sprintf(configfile,"%s",optarg); break;

    case 'i': 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;

    default: abort();

    }

  }

  std::cout << "Used: " << argv[0] <<  " -c " << configfile << " -o " << outputfile << " -n " << nevetoread << " -w " << writentuple << " -d  " << verbose ;

  for (unsigned int i=0;i<  inputfilelist.size() ;i++) std::cout << " -i " << inputfilelist[i];

  std::cout << std::endl;

  //-----------------------------------------------------------------

  new readdata(configfile, inputfilelist,outputfile, writentuple, nevetoread, verbose);

  std::cout << "Output data files " << outputfile <<  std::endl;

  return 0;

}

#endif