WebSVN - f9daq - Rev 264 - /cvi/apps/petdemo/analyse.c

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

Subversion Repositories f9daq

(root)/cvi/apps/petdemo/analyse.c @ 360 - Rev 264