Rev 1 | Blame | Compare with Previous | Last modification | View Log | RSS feed
/*vhodne datoteke:Mappingi:m16.map channel mappingmodules.map mapping modulovKalibracija:pedestals.dat ... adc sumaph511.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 3typedef 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 filetypedef 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 positionclass 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 generatorTH1I *m_crystalid[4];// kalibracijski histogramichar * m_modulesmapName; // ime datoteke s podatki o pozicijah modulovchar * m_channelsmapName; // ime datoteke s podatki o pozicijah elektronskih kanalov//char * m_sumpedestalsName;char * m_pedestalsName; // ime datoteke s podatki o pedestalihchar * m_photopeakName; // ime datoteke s podatki o vrhovihchar * m_calibrationrootName; // ime datoteke s kalibracijskimi histogramiTFile *m_calibration; // datoteke s kalibracijskimi histogramidouble m_dx; // pitch x kristalovdouble m_dy; // pitch y kristalovint m_nx; // stevilo kristalov v x smeriint m_ny; // stevilo kristalov v y smerifloat gPedestals[MAXCH];float gPeak[MAXCH];//float gSumPedestals[MAXCH];double m_peakscaling; // scaling za adcjedouble m_threshold; // threshold za upostevanje zadetkov na kanalihint 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 kristaluint 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 kristaladouble 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 kanalovfloat 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 processpublic:PETProjDataMgr *m_projector; // projektor za sinogramedouble m_rotation; // trenutna rotacija setupadouble m_threshold; //double gData[MAXCH]; // korigirani ADC jidouble gAdc[MAXCH]; // raw ADC jiint m_debug; // debug izpisiint m_write; // ce je ta flag nastavljen, se v root file izpisuje ntupleint 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) * scalingsum0 .. sum3 vsota energij po vseh kanalih fotopomnozevalkenh0 .. nh4 stevilo zadetkov na fotopomnozevalki (adc nad thresholdom)c0 .. c4 vsota energij za kanale nad thresholdompx0 .. px4 x koordinata COG na fotopomnozevalkipy0 .. py4 y koordinata COG na fotopomnozevalkimax0 .. max3 maximalna energija na kanalu na fotopomnozevalkia0 .. 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 thresholdomm_GlobalPosition = new TH2F("globalxy","Rekonstruirana koordinata zadetka",200,-80,80,200,-80,80); // reconstruirana koordinata v globalnem sistemuchar 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 adcjisprintf(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 PMTjusprintf(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 PMTjusprintf(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 PMTjusprintf(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 PMTjusprintf(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 PMTjusprintf(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 PMTjusprintf(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 thresholdomsprintf(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 pmtjusprintf(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 mappingTH1I *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 pmtjevint j2= (ipmt/2)*2+1-ipmt%2; // sosednja fotopomnozevalkafloat 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 fotopomnozevalkiint 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 fotopomnozevalkiposx[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 pmtjufloat 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 MAINint 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