Subversion Repositories f9daq

Rev

Blame | Last modification | View Log | RSS feed

  1.  
  2. #include <stdlib.h>
  3. #include <stdio.h>
  4. #include <stdint.h>
  5. #include "daq.h"
  6. #include "H1D.h"
  7. #include "H2D.h"
  8. #include "PETProjDataMgr.h"
  9.  
  10. double gAdc[MAXCH];  // raw ADC ji
  11. float  gNtdata[MAXCH*3];
  12. int m_write=0;
  13.  
  14. double gData[MAXCH]; // korigirani ADC ji
  15.  
  16. float gSum[MAXPMT];
  17. float gRawSum[MAXPMT];
  18. float gMax[MAXPMT];
  19. float gSumCluster[MAXPMT];
  20. int   gNabove[MAXPMT];
  21.  
  22.  
  23.  
  24.  
  25.  
  26. double GetEnergy(int ch, int adc) {
  27.         return (adc-conf.apedestals[ch])/(conf.apeak[ch]-conf.apedestals[ch])*conf.peakscaling;
  28. }
  29.  
  30.  
  31.  
  32. int  GetGlobalPosition(int ipmt, double angle, float cx, float cy, float *gposition) {
  33.  
  34.         double phi = angle + conf.module[ipmt].phi;
  35.         double r   = conf.module[ipmt].r;
  36.         double sinphi = sin(phi);
  37.         double cosphi = cos(phi);
  38.  
  39.         gposition[0]= cosphi* r + sinphi * cx ;
  40.         gposition[1]= sinphi* r - cosphi * cx ;
  41.         gposition[2]= cy;
  42.  
  43.         //if (m_debug) printf(  "[%d] phi=%f %f global ( %f , %f , %f )\n" ,ipmt,  module[ipmt].phi,angle, rv.x(), rv.y(), rv.z());
  44.         return  0;
  45. }
  46.  
  47.  
  48. int GetRealPosition(int ipmt, float px,float py,float *cx,float *cy) {
  49.         // iz CoG izracuna pozicijo na kristalu
  50.  
  51.         int binx = 0; //####m_crystalid[ipmt]->GetXaxis()->FindBin(px);
  52.         int biny = 0; //####m_crystalid[ipmt]->GetYaxis()->FindBin(py);
  53.         int crystalid  = 0;//####m_crystalid[ipmt]->GetBinContent(binx,biny);
  54.         *cx= ( crystalid % conf.nofcrystalsx - conf.nofcrystalsy/2. + 0.5 ) * conf.crystalpitchx;
  55.         *cy= ( crystalid / conf.nofcrystalsx - conf.nofcrystalsx/2. + 0.5 ) * conf.crystalpitchy;
  56.  
  57.  
  58.         // razmazi pozicijo po povrsini kristala
  59.         double rndx = Random(-0.5,0.5)* conf.crystalpitchx;
  60.         double rndy = Random(-0.5,0.5)* conf.crystalpitchy;
  61.         *cx+= rndx;
  62.         *cy+= rndy;
  63.         return crystalid;
  64. }
  65.  
  66.  
  67. int HistogramsInit() {
  68.  
  69.         if (conf.write) {
  70.                 char varlist[1024], tmplist[1024];
  71.                 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");
  72.                 for (int i=0; i<MAXCH; i++) {
  73.                         sprintf(tmplist,"%s",varlist);
  74.                         sprintf(varlist,"%s:a%d",tmplist,i);
  75.                 }
  76.                 /*
  77.                         energija .. (adc -pedestal) * scaling
  78.  
  79.                         sum0 .. sum3 vsota energij po vseh kanalih fotopomnozevalke
  80.                         nh0  .. nh4   stevilo zadetkov na fotopomnozevalki (adc nad thresholdom)
  81.                         c0   .. c4 vsota energij za kanale nad thresholdom
  82.                         px0  .. px4 x koordinata COG na fotopomnozevalki
  83.                         py0  .. py4 y koordinata COG na fotopomnozevalki
  84.                         max0 .. max3 maximalna energija na kanalu na fotopomnozevalki
  85.                         a0   .. a63 energija na i-tem kanalu
  86.                 */
  87.                 //gNtuple = new TNtuple("nt","Pet RAW data",varlist);
  88.                 printf ("#%s#\n",varlist);
  89.         }
  90.         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
  91.         m_GlobalPosition =H2D_Init(302,"globalxy","Rekonstruirana koordinata zadetka",200,-80,80,200,-80,80);  // reconstruirana koordinata v globalnem sistemu
  92.  
  93.  
  94.         char hname[0xFF];
  95.         char hn[0xFF];
  96.         int i;
  97.         for (i=0; i<MAXCH; i++) {
  98.                 sprintf(hname,"Raw ADC Ch. %d ;ADC;N",i);
  99.                 sprintf(hn,"ach%d",i);
  100.                 //m_Adc[i]   = H1D_Init(i,hn,hname,4000,-0.5,3999.5);       // osnovni adcji
  101.                 m_Adc[i]   = H1D_Init(i,hn,hname,400,500-0.5,899.5);       // osnovni adcji
  102.                 sprintf(hname,"cADC za dogodke z manj kot x hitov na pmtju Ch. %d ;ADC;N",i);
  103.                 sprintf(hn,"cutch%d",i);
  104.                 m_AdcCut[i]   = H1D_Init(i+MAXCH, hn,hname,4000,-0.5,3999.5);    // adcji za zadetke z manj kot 7 hiti na PMTju
  105.                 sprintf(hname,"cADC proti stevilu zadetkov na pmtju Ch. %d ;ADC;N",i);
  106.                 sprintf(hn,"singlech%d",i);
  107.                 m_Adc_vs_Nhits[i]   = H2D_Init(i,hn,hname,200,0,4000, 17,-0.5,16.5); // adc proti stevilu zadetkov na PMTju
  108.                 sprintf(hname,"cADC proti vsoti kanalov na pmtju,  Ch. %d ;ADC;ADCsum",i);
  109.                 sprintf(hn,"corrch%d",i);
  110.                 m_Adc_vs_Sum[i]   = H2D_Init(i+MAXCH,hn,hname,200,0,4000, 200,0,12000 );  // raw adc proti vsoti kanalov na PMTju
  111.                 sprintf(hname,"Raw ADC proti vsoti kanalov na pmtju,  Ch. %d ;ADC;ADCsum",i);
  112.                 sprintf(hn,"adcvssumch%d",i);
  113.                 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
  114.         }
  115.  
  116.         for (i=0; i<MAXPMT; i++) {
  117.                 sprintf(hname,"Vsota cADC na PMTju %d ;ADC;N",i);
  118.                 sprintf(hn,"pmt%d",i);
  119.                 m_AdcSum[i]   =H1D_Init(200+i,hn,hname,200,0,20000); // vsota adcjev na PMTju
  120.                 sprintf(hname,"Vsota cADC na PMTju za dogodke z manj kot x zadetki %d ;ADC;N",i);
  121.                 sprintf(hn,"clusterpmt%d",i);
  122.                 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
  123.  
  124.                 sprintf(hname,"Center naboja CoG PMT %d ;x;y",i);
  125.                 sprintf(hn,"pmt1%d",i);
  126.                 m_CenterOfGravity[i]   =H2D_Init(200+i,hn,hname,200,-1.25,4.25,200,-1.25,4.25); // center naboja na PMTju
  127.  
  128.                 sprintf(hname,"Center naboja CoG za zadetke nad thresholdom PMT %d ;x;y",i);
  129.                 sprintf(hn,"pmt2%d",i);
  130.                 m_CenterOfGravityforChAboveThreshold[i] = H2D_Init(216+i,hn,hname,200,0,3,200,0,3);  // center naboja za zadetke nad thresholdom
  131.  
  132.                 sprintf(hname,"Rekonstruirana koordinata PMT %d ;x(mm);y(mm)",i);
  133.                 sprintf(hn,"pmt3%d",i);
  134.                 m_ReconstructedPosition[i]   = H2D_Init(232+i,hn,hname,200,-12,12,200,-12,12); // korigirana koordinata na pmtju
  135.  
  136.                 sprintf(hname,"SumAdc PMT %d ;SumAdc(a.u.);crystal ID",i);
  137.                 sprintf(hn,"sumadc%d",i);
  138.                 m_SumAdc[i]   =H2D_Init(248+i,hn,hname,200,0,12000,81,-0.5,80.5); // SumADC spekter za razlicne kristale
  139.         }
  140.  
  141.         // store mapping
  142.         //TH1I *hx= new TH1I("mapping", "ch mapping", 16, -0.5,15.5);
  143.         //for (int j=0;j<16;j++) hx->Fill(j,m_geo->ch[j].idx);
  144.  
  145.         return 0;
  146. };
  147.  
  148. int HistogramsWrite(const char *fname){
  149.         FILE *fp = fopen(fname, "wb");
  150.        
  151.         if (fp){
  152.           int i;
  153.                 for (i=0;i<H1DMAX;i++){ if (H1D_Exist(i)) H1D_Write2File(i, fp); }
  154.                 for (i=0;i<H2DMAX;i++){ if (H2D_Exist(i)) H2D_Write2File(i, fp); }
  155.                 for (i=0;i<H3DMAX;i++){ if (H3D_Exist(i)) H3D_Write2File(i, fp); }
  156.                 fclose(fp);
  157.        
  158.         }
  159.  
  160.         return 0;
  161. }
  162.  
  163.  
  164. int HistogramsFill() {
  165.         int i=0, ipmt;
  166.         for (i=0; i<MAXCH; i++) { // zanka cez vse elektronske kanale;
  167.                 float x=gData[i];
  168.                 int ipmt= i/16;
  169.                 float s=gSum[ipmt];
  170.                 // int idx= 2*(ipmt)+64;
  171.                 // const  float fSlope1=2;
  172.                 // const  float fSlope2=3;
  173.                 //  if ( ((s-gSumPedestals[ipmt] ) > fSlope1*(x-gPedestals[i])) && ((s-gSumPedestals[ipmt] ) < fSlope2*(x-gPedestals[ipmt])) ) {
  174.                 if (gNabove[ipmt]<7) {
  175.                         H2D_Fill(m_Adc_vs_Sum[i],x,s,1);
  176.                         H2D_Fill(m_Adc_vs_Sum_Uncorected[i], gAdc[i],gRawSum[ipmt],1);
  177.                         H1D_Fill(m_AdcCut[i],x,1);
  178.                 }
  179.         }
  180.  
  181.         HVector3 r_coincidence[2];
  182.         for (int k=0; k<2; k++) for (int i1=0; i1<3; i1++)  r_coincidence[k].x[i1]=0;
  183.         int  f_coincidence[2]= {0,0};
  184.  
  185.         for (ipmt=0; ipmt<4; ipmt++) { // zanka preko pmtjev
  186.                 int j2= (ipmt/2)*2+1-ipmt%2;  // sosednja fotopomnozevalka
  187.  
  188.                 float posx[2]= {0,0};
  189.                 float posy[2]= {0,0};
  190.                 float sum[2]= {0,0};
  191.                 for (int ich=0; ich<16; ich++) { // zanka preko elektronskih kanalov na fotopomnozevalki
  192.                         int ch= ich+ipmt*16;
  193.                         if (gMax[ipmt]>conf.adcthreshold) {
  194.                                 posx[0]+= gData[ch]*conf.channel[ich].ix;
  195.                                 posy[0]+= gData[ch]*conf.channel[ich].iy;
  196.                                 sum[0] += gData[ch];
  197.  
  198.                                 if (gData[ch]> 0.1*gMax[ipmt]) { // pri racunanju pozicije upostevaj le kanale, ki imajo vrednost vecjo od ratio*maksimalna na tisti fotopomnozevalki ********
  199.                                         posx[1]+= gData[ch]*conf.channel[ich].ix;
  200.                                         posy[1]+= gData[ch]*conf.channel[ich].iy;
  201.                                         sum[1] += gData[ch];
  202.                                 }
  203.                         }
  204.                 }
  205.  
  206.                 if (m_write) {
  207.                         gNtdata[12+ipmt]=posx[0]/sum[0];
  208.                         gNtdata[16+ipmt]=posy[0]/sum[0];
  209.                         gNtdata[20+ipmt]=(sum[1]>0)?posx[1]/sum[1]:-10;
  210.                         gNtdata[24+ipmt]=(sum[1]>0)?posy[1]/sum[1]:-10;
  211.                         gNtdata[28+ipmt]=sum[0];
  212.                 }
  213.                 if (gMax[ipmt]>gMax[j2]) { // izberi fotopomnozevalko, ki ima vec zadetkov kot soseda TRG= (P1|P2) & (P3|P4)
  214.  
  215.                         if ( sum[0] > 0  ) {
  216.                                 float px=posx[0]/sum[0];
  217.                                 float py=posy[0]/sum[0];
  218.  
  219.                                 H2D_Fill(m_CenterOfGravity[ipmt],px,py,1);
  220.  
  221.                                 float cx=0; //koordinata na pmtju
  222.                                 float cy=0;
  223.                                 int crystalID = GetRealPosition(ipmt,px,py,&cx,&cy);
  224.                                 H2D_Fill(m_ReconstructedPosition[ipmt],cx,cy,1);
  225.                                 H2D_Fill(m_SumAdc[ipmt],gSum[ipmt], crystalID,1 );
  226.                                 if (debug > 2) printf("ipmt=%d cx=%f cy=%f\n",ipmt,cx,cy);
  227.                                 GetGlobalPosition(ipmt, conf.rotation, cx, cy,r_coincidence[ipmt/2].x) ;
  228.                                 f_coincidence[ipmt/2] =1;
  229.                                 H2D_Fill(m_GlobalPosition,r_coincidence[ipmt/2].x[0],r_coincidence[ipmt/2].x[1],1);
  230.  
  231.                         }
  232.                         if ( sum[1] > 0  ) {
  233.  
  234.                                 float px=posx[1]/sum[1];
  235.                                 float py=posy[1]/sum[1];
  236.                                 // float px0=posx[0]/sum[0];
  237.                                 // float py0=posy[0]/sum[0];
  238.                                 // if (px0<0.25 || px0>2.75 ||  py0< 0.25 || py0> 2.75)
  239.                                 // if  (((py0>0.75 && py0<1.25) || (py0>1.75 && py0<2.25))  &&  (px0< 0.25 || px0>2.75 ) )
  240.                                 H2D_Fill(m_CenterOfGravityforChAboveThreshold[ipmt],px,py,1);
  241.                         }
  242.  
  243.                 }
  244.         }
  245.  
  246.         for (i=0; i<MAXCH; i++) if (gData[i]>conf.adcthreshold) {
  247.                         H2D_Fill(m_Adc_vs_Nhits[i],gData[i],gNabove[i/16],1);
  248.                 }
  249.         for (i=0; i<MAXCH; i++) if (gData[i]>conf.adcthreshold && gNabove[i/16]<5 ) {
  250.                         gSumCluster[i/16]+=gData[i];
  251.                 }
  252.         for (ipmt=0; ipmt<4; ipmt++) {
  253.                 H1D_Fill(m_AdcSum[ipmt],gSum[ipmt],1);
  254.                 H1D_Fill(m_AdcSumCluster[ipmt],gSumCluster[ipmt],1);
  255.         }
  256.  
  257.         for (i=0; i<4; i++) H1D_Fill(m_Nhits,i,gNabove[i]);
  258.  
  259.         if (m_write) {
  260.                 for (i=0; i<4; i++) {
  261.                         gNtdata[i]=gSum[i];
  262.                         gNtdata[4+i]=gNabove[i];
  263.                         gNtdata[8+i]=gSumCluster[i];
  264.                 }
  265.                 for (i=0; i<MAXCH; i++) gNtdata[32+i]=gData[i];
  266.                 //gNtuple->Fill(gNtdata);
  267.         }
  268.         //
  269.         // coincidences
  270.         //
  271.         if (f_coincidence[0]&&f_coincidence[1]) {
  272.                 conf.coincidences++;
  273.                 if (debug > 1) {
  274.                         printf(  "Coincidence ( %f , %f , %f ) ( %f , %f , %f )\n" ,
  275.                                                          r_coincidence[0].x[0], r_coincidence[0].x[1], r_coincidence[0].x[2],
  276.                                                          r_coincidence[1].x[0], r_coincidence[1].x[1], r_coincidence[1].x[2] );
  277.                 }
  278.                 AddEvent( r_coincidence[0] , r_coincidence[1]);
  279.         }
  280.         return 0;
  281. };
  282.  
  283.  
  284.  
  285. int analyse(int nb, uint32_t *buf) {
  286.  
  287.         int idx=1;
  288.         int neve=buf[0]/2;
  289.         //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  290.         for (int nev=0; nev<neve; nev++) {
  291.                 int len=buf[idx];
  292.                 int sb =buf[idx+1];
  293.                 unsigned int *pbuf=&buf[idx+2];
  294.                 if (sb!=0xffab) {
  295.                         printf("0x%04x!0xffab len=%d\n",sb,len);
  296.                         break;
  297.                 }
  298.                 // postavi na nic
  299. #define InitArrayWithAValue(arr,n,x) {for (int i=0;i<n;i++) arr[i]=x;}
  300.  
  301.                 //------------------------------------------------------------
  302.                 // postavi na nic
  303.  
  304.                 InitArrayWithAValue( gSum       , MAXPMT    , 0);
  305.                 InitArrayWithAValue( gRawSum    , MAXPMT    , 0);
  306.                 InitArrayWithAValue( gMax       , MAXPMT    , 0);
  307.                 InitArrayWithAValue( gSumCluster, MAXPMT    , 0);
  308.                 InitArrayWithAValue( gNabove    , MAXPMT    , 0);
  309.                 InitArrayWithAValue( gData      , MAXCH, 0);
  310.                 InitArrayWithAValue( gAdc       , MAXCH, 0);
  311.                 for (int i0=0; i0<len-2; i0+=2) {
  312.  
  313.                         int data0  = pbuf[i0];
  314.                         int data1  = pbuf[i0+1];
  315.                         int geo    = (data1 >> 11) & 0x1f;
  316.                         int ch     = (data1&0x1f)  | (geo<<5);
  317.                         int dtype  = (data1>>9)&0x3;
  318.                         int adc    =  data0&0xfff;
  319.                         if (dtype == 0 & ch<MAXCH) {
  320.                                 if (debug)      printf("%d ADC %d %d\n", geo, ch, adc);
  321.                                 if (ch<MAXCH) H1D_Fill(ch, adc,1);
  322.                                 int ipmt = ch/16;
  323.  
  324.                                 gAdc[ch]=adc;
  325.                                 gRawSum[ipmt]+=adc;
  326.  
  327.                                 gData[ch]=GetEnergy(ch,adc);
  328.  
  329.                                 gSum[ipmt]+=gData[ch];
  330.                                 if (gData[ch]   >gMax[ipmt] )     gMax[ipmt]= gData[ch];
  331.                                 if (gData[ch]   >conf.adcthreshold )    gNabove[ipmt]++;
  332.                                 //#### m_Adc[ch]->Fill(adc);
  333.                         }
  334.                 };// for (int i0=0;i0<len-2;i0+=2)
  335.                 //------------------------------------------------------------
  336.  
  337.                 idx+=len+1;
  338.                 HistogramsFill();
  339.                 if (m_write) {
  340.                         gNtdata[0]=0;
  341.                         for (int i=0; i<MAXCH; i++) gNtdata[i+1]=gAdc[i];
  342.                         //####gNtuple->Fill(gNtdata);
  343.                 }
  344.  
  345.         } // for (int nev=0;nev<neve;nev++)
  346.  
  347.  
  348.  
  349.  
  350.         return 0;
  351. }
  352.