Subversion Repositories f9daq

#include <stdlib.h>

#include <stdio.h>

#include <stdint.h>

#include "daq.h"

#include "H1D.h"

#include "H2D.h"

#include "PETProjDataMgr.h"

double gAdc[MAXCH];  // raw ADC ji

float  gNtdata[MAXCH*3];

int m_write=0;

double gData[MAXCH]; // korigirani ADC ji

float gSum[MAXPMT];

float gRawSum[MAXPMT];

float gMax[MAXPMT];

float gSumCluster[MAXPMT];

int   gNabove[MAXPMT];

double GetEnergy(int ch, int adc) {

        return (adc-conf.apedestals[ch])/(conf.apeak[ch]-conf.apedestals[ch])*conf.peakscaling;

}

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

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

        double r   = conf.module[ipmt].r;

        double sinphi = sin(phi);

        double cosphi = cos(phi);

        gposition[0]= cosphi* r + sinphi * cx ;

        gposition[1]= sinphi* r - cosphi * cx ;

        gposition[2]= 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  0;

}

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

        // iz CoG izracuna pozicijo na kristalu

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

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

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

        *cx= ( crystalid % conf.nofcrystalsx - conf.nofcrystalsy/2. + 0.5 ) * conf.crystalpitchx;

        *cy= ( crystalid / conf.nofcrystalsx - conf.nofcrystalsx/2. + 0.5 ) * conf.crystalpitchy;

        // razmazi pozicijo po povrsini kristala

        double rndx = Random(-0.5,0.5)* conf.crystalpitchx;

        double rndy = Random(-0.5,0.5)* conf.crystalpitchy;

        *cx+= rndx;

        *cy+= rndy;

        return crystalid;

}

int HistogramsInit() {

        if (conf.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 = H2D_Init(301,"hnabove","Stevilo zadetkov nad nivojem diskriminacije",4,-0.5,3.5, 16,+0.5,16.5 );  // stevilo hitov nad thresholdom

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

        char hname[0xFF];

        char hn[0xFF];

        int i;

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

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

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

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

                m_Adc[i]   = H1D_Init(i,hn,hname,400,500-0.5,899.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]   = H1D_Init(i+MAXCH, 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]   = H2D_Init(i,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]   = H2D_Init(i+MAXCH,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]   =H2D_Init(i+MAXCH*2,hn,hname,100,0,3500, 100,5000,12000 );  // adc proti vsoti kanalov na PMTju

        }

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

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

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

                m_AdcSum[i]   =H1D_Init(200+i,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]   = H1D_Init(216+i,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]   =H2D_Init(200+i,hn,hname,200,-1.25,4.25,200,-1.25,4.25); // center naboja na PMTju

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

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

                m_CenterOfGravityforChAboveThreshold[i] = H2D_Init(216+i,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]   = H2D_Init(232+i,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]   =H2D_Init(248+i,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 HistogramsWrite(const char *fname){

        FILE *fp = fopen(fname, "wb");

        if (fp){

          int i;

                for (i=0;i<H1DMAX;i++){ if (H1D_Exist(i)) H1D_Write2File(i, fp); }

                for (i=0;i<H2DMAX;i++){ if (H2D_Exist(i)) H2D_Write2File(i, fp); }

                for (i=0;i<H3DMAX;i++){ if (H3D_Exist(i)) H3D_Write2File(i, fp); }

                fclose(fp);

        }

        return 0;

}

int HistogramsFill() {

        int i=0, ipmt;

        for (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) {

                        H2D_Fill(m_Adc_vs_Sum[i],x,s,1);

                        H2D_Fill(m_Adc_vs_Sum_Uncorected[i], gAdc[i],gRawSum[ipmt],1);

                        H1D_Fill(m_AdcCut[i],x,1);

                }

        }

        HVector3 r_coincidence[2];

        for (int k=0; k<2; k++) for (int i1=0; i1<3; i1++)  r_coincidence[k].x[i1]=0;

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

        for (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]>conf.adcthreshold) {

                                posx[0]+= gData[ch]*conf.channel[ich].ix;

                                posy[0]+= gData[ch]*conf.channel[ich].iy;

                                sum[0] += gData[ch];

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

                                        posx[1]+= gData[ch]*conf.channel[ich].ix;

                                        posy[1]+= gData[ch]*conf.channel[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];

                                H2D_Fill(m_CenterOfGravity[ipmt],px,py,1);

                                float cx=0; //koordinata na pmtju

                                float cy=0;

                                int crystalID = GetRealPosition(ipmt,px,py,&cx,&cy);

                                H2D_Fill(m_ReconstructedPosition[ipmt],cx,cy,1);

                                H2D_Fill(m_SumAdc[ipmt],gSum[ipmt], crystalID,1 );

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

                                GetGlobalPosition(ipmt, conf.rotation, cx, cy,r_coincidence[ipmt/2].x) ;

                                f_coincidence[ipmt/2] =1;

                                H2D_Fill(m_GlobalPosition,r_coincidence[ipmt/2].x[0],r_coincidence[ipmt/2].x[1],1);

                        }

                        if ( sum[1] > 0  ) {

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

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

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

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

                                // if (px0<0.25 || px0>2.75 ||  py0< 0.25 || py0> 2.75)

                                // if  (((py0>0.75 && py0<1.25) || (py0>1.75 && py0<2.25))  &&  (px0< 0.25 || px0>2.75 ) )

                                H2D_Fill(m_CenterOfGravityforChAboveThreshold[ipmt],px,py,1);

                        }

                }

        }

        for (i=0; i<MAXCH; i++) if (gData[i]>conf.adcthreshold) {

                        H2D_Fill(m_Adc_vs_Nhits[i],gData[i],gNabove[i/16],1);

                }

        for (i=0; i<MAXCH; i++) if (gData[i]>conf.adcthreshold && gNabove[i/16]<5 ) {

                        gSumCluster[i/16]+=gData[i];

                }

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

                H1D_Fill(m_AdcSum[ipmt],gSum[ipmt],1);

                H1D_Fill(m_AdcSumCluster[ipmt],gSumCluster[ipmt],1);

        }

        for (i=0; i<4; i++) H1D_Fill(m_Nhits,i,gNabove[i]);

        if (m_write) {

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

                        gNtdata[i]=gSum[i];

                        gNtdata[4+i]=gNabove[i];

                        gNtdata[8+i]=gSumCluster[i];

                }

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

                //gNtuple->Fill(gNtdata);

        }

        //

        // coincidences

        //

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

                conf.coincidences++;

                if (debug > 1) {

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

                                                         r_coincidence[0].x[0], r_coincidence[0].x[1], r_coincidence[0].x[2],

                                                         r_coincidence[1].x[0], r_coincidence[1].x[1], r_coincidence[1].x[2] );

                }

                AddEvent( r_coincidence[0] , r_coincidence[1]);

        }

        return 0;

};

int analyse(int nb, uint32_t *buf) {

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

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

                // postavi na nic

                InitArrayWithAValue( gSum       , MAXPMT    , 0);

                InitArrayWithAValue( gRawSum    , MAXPMT    , 0);

                InitArrayWithAValue( gMax       , MAXPMT    , 0);

                InitArrayWithAValue( gSumCluster, MAXPMT    , 0);

                InitArrayWithAValue( gNabove    , MAXPMT    , 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;

                        if (dtype == 0 & ch<MAXCH) {

                                if (debug)      printf("%d ADC %d %d\n", geo, ch, adc);

                                if (ch<MAXCH) H1D_Fill(ch, adc,1);

                                int ipmt = ch/16;

                                gAdc[ch]=adc;

                                gRawSum[ipmt]+=adc;

                                gData[ch]=GetEnergy(ch,adc);

                                gSum[ipmt]+=gData[ch];

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

                                if (gData[ch]   >conf.adcthreshold )    gNabove[ipmt]++;

                                //#### m_Adc[ch]->Fill(adc);

                        }

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

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

                idx+=len+1;

                HistogramsFill();

                if (m_write) {

                        gNtdata[0]=0;

                        for (int i=0; i<MAXCH; i++) gNtdata[i+1]=gAdc[i];

                        //####gNtuple->Fill(gNtdata);

                }

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

        return 0;

}