Subversion Repositories f9daq

Rev

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

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