Subversion Repositories f9daq

#include "TH1.h"

#include "TH2.h"

#include "TGraph.h"

#include "TCanvas.h"

#include "TStyle.h"

#include "TLegend.h"

//#include "TROOT.h"

#include <stdio.h>

#define DATA_PATH  "./data/"

#define CFRAC 0.4

#define DWIDTH 1024

void drs(int nev=1000, int updfrq=20, const char *FileName="drs000.dat", int trgch=0, int debug=0)

{

  int i,ich,iev;

  int stat;

  char histname[128];

//  char FileName[300]="./data/n2.dat";

  FILE *fin;

  char recid[5]="SAMO";

  unsigned short *srecid = (unsigned short *) &recid[0];

  struct evrec {

    unsigned long evn;

    unsigned short tYear,tMonth,tDay,tHour,tMin,tSec,tmSec,rc;

  } evrec;

  float dtdata[4][DWIDTH];

  float tcdata[4][DWIDTH];

  float tccalib[4];

  float scaler[4];

  float *t;

  float *y;

  float *ay;

  short lwfm[4];

  int nwfm[4];

  int board;

  int version;

  unsigned short ichdat[4][DWIDTH];

  float chdat[4][DWIDTH];

  float achdat[4][DWIDTH];

  float chmax[4][2];

  float chmin[4][2];

  float letime[4][2],cftime[4][2],qdc[4];

//  float threshold[4][2]={{0.,-0.02},{0.,0.},{0.,0.},{0.,-0.25}};

//  float twin[4][2]={{120.,130.},{100.,150.},{100.,150.},{70.,90.}};

//  float adcgate[4][2]={{124.,136.},{120.,140.},{120.,140.},{70.,90.}};

//  float threshold[4][2]={{0.,-0.05},{0.,-0.05},{0.,-0.25},{0.,-0.25}};

//  float twin[4][2]={{80.,100.},{80.,100.},{90.,110.},{50.,70.}};

//  float adcgate[4][2]={{80.,110.},{80.,110.},{90.,110.},{50.,70.}};

  float threshold[4][2]={{0.,-0.1  },{0.,-0.05 },{0.,-0.25 },{0.,-0.25}};

  float twin[4][2]     ={{100.,120.},{100.,120.},{40.,130. },{60.,80. }};

  float adcgate[4][2]  ={{105.,125.},{105.,125.},{110.,130.},{60.,80. }};

  float vcut[4][2]     ={{-1,1     },{-1,1     },{-1,1     },{-1,1    }};

//  gROOT->Reset();

  gStyle->SetOptStat(1111111);

//  gStyle->SetPalette(52);

  TH1F *h1_chmax[4];

  TH1F *h1_chmaxt[4];

  TH1F *h1_chmin[4];

  TH1F *h1_chmint[4];

  TH1F *h1_letime[4][2];

  TH1F *h1_cftime[4][2];

  TH1F *h1_cftdif[4];

  TH1F *h1_qdc[4];

  TH2F *h2_dpo[4];

  TH2F *h2_tdcadc[4];

  TH2F *h2_ctdcadc[4];

  TH2F *h2_ctalk;

  TGraph *gr_wfm[4];

  TCanvas *c_wfm= new TCanvas("c_wfm","WFM",50,50,700,700);

  TLegend *leg = new TLegend(0.8, 0.8, 1, 1);

  c_wfm->Divide(1,1);

  c_wfm->cd(1);

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

    sprintf(histname,"chmax_%02d",i);

    h1_chmax[i]=new TH1F(histname,histname,1100,-0.55,0.55);

    sprintf(histname,"chmaxt_%02d",i);

    h1_chmaxt[i]=new TH1F(histname,histname,DWIDTH,-0.1,204.7);

    sprintf(histname,"chmin_%02d",i);

    h1_chmin[i]=new TH1F(histname,histname,1100,-0.55,0.55);

    sprintf(histname,"chmint_%02d",i);

    h1_chmint[i]=new TH1F(histname,histname,DWIDTH,-0.1,204.7);

    sprintf(histname,"letimep_%02d",i);

    h1_letime[i][0]=new TH1F(histname,histname,500,twin[i][0],twin[i][1]);

    sprintf(histname,"letimen_%02d",i);

    h1_letime[i][1]=new TH1F(histname,histname,500,twin[i][0],twin[i][1]);

    sprintf(histname,"cftimep_%02d",i);

    h1_cftime[i][0]=new TH1F(histname,histname,500,twin[i][0],twin[i][1]);

    sprintf(histname,"cftimen_%02d",i);

    h1_cftime[i][1]=new TH1F(histname,histname,500,twin[i][0],twin[i][1]);

    sprintf(histname,"cftdif_%02d",i);

    h1_cftdif[i]=new TH1F(histname,histname,2000,twin[i][0],twin[i][1]);

    sprintf(histname,"qdc_%02d",i);

    h1_qdc[i]=new TH1F(histname,histname,500,-0.1,4.9);

    sprintf(histname,"dpo_%02d",i);

    h2_dpo[i]=new TH2F(histname,histname,DWIDTH,-0.1,204.7,1100,-0.55,0.55);

    sprintf(histname,"tdcadc_%02d",i);

    h2_tdcadc[i]=new TH2F(histname,histname,250,0.,0.5,400,twin[i][0],twin[i][1]);

    sprintf(histname,"ctdcadc_%02d",i);

    h2_ctdcadc[i]=new TH2F(histname,histname,250,0.,0.5,400,twin[i][0],twin[i][1]);

    gr_wfm[i]=new TGraph(DWIDTH);

    gr_wfm[i]->SetMinimum(-0.55);

    gr_wfm[i]->SetMaximum(0.55);

    gr_wfm[i]->SetLineColor(i+1);

    gr_wfm[i]->SetMarkerColor(i+1);

    sprintf(histname,"ch %d",i);

    leg->AddEntry(gr_wfm[i], histname, "l");

    if (i) gr_wfm[i]->Draw("LP"); else gr_wfm[i]->Draw("ALP");

  }

  leg->Draw();

  h2_ctalk=new TH2F("ctalk","ctalk",250,0.,0.5,250,0.,0.5);

  fin=fopen(FileName,"rb");

  if (fin==NULL) {

	  printf("Error opening file %s\n",FileName);

	  return;

  }

  printf("%s, %lu\n",recid,sizeof(evrec));

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

    nwfm[ich]=0;

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

      achdat[ich][i]=0;

    }

  }

  unsigned int trgcell = 0;

  int event = -1;

  iev=0;

  while(!feof(fin)){

    if (iev>=nev) break;

    stat=fread(recid,1,4,fin);

    if (debug){

      if (recid[1]=='#') printf("%c%c%u\n",recid[0],recid[1],srecid[1]);

      else  printf("%s\n",recid );

    }

    switch (recid[0]){

	   case 'D': // DRS

           version = atoi(&recid[3]);

           break;

       case 'T': // TIME

           switch (recid[1]){

             case 'I':

			   event=0;

               if (version ==0) version++;

   			   break;

             case '#': trgcell = srecid[1]; break;

           }

           break;

       case 'B':

           board = srecid[1];

           break;

       case 'C':

           sscanf(recid,"C%d",&ich);

           ich=ich-1;

           if (event) {

			 if (version>1) stat=fread(&scaler[ich],1,sizeof(float),fin);

             stat=fread(ichdat[ich],1,sizeof(ichdat[ich]),fin);

           } else {

             stat=fread(dtdata[ich],1,sizeof(float)*DWIDTH,fin);

           }

           break;

       case 'E':

           stat=fread(&evrec.evn,1,sizeof(evrec),fin);

           iev++;

		   if (!version)  {

			 stat=fread(dtdata[0],1,sizeof(float)*DWIDTH,fin);

		   }

           event=1;

           break;

       default:

           printf("Unknown header! %s\n", recid);

           break;

    }

    if (recid[0] == 'E'){

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

           if (!lwfm[i]) continue;

           if ((chmin[i][0]>vcut[i][0])&&(chmin[i][0]<vcut[i][1])) {

             h1_cftdif[i]->Fill(cftime[i][1]-cftime[trgch][1]+110.);

           }

//      h2_tdcadc[i]->Fill(-chmin[i][0],letime[i][1]-cftime[trgch][1]+72);

//      h2_ctdcadc[i]->Fill(-chmin[i][0],cftime[i][1]-cftime[trgch][1]+72);

           h2_tdcadc[i]->Fill(-chmin[0][0],letime[i][1]-cftime[trgch][1]+72);

           h2_ctdcadc[i]->Fill(-chmin[0][0],cftime[i][1]-cftime[trgch][1]+72);

         }

         h2_ctalk->Fill(-chmin[0][0],-chmin[1][0]);

    }

    if (event<1) continue;

    if (recid[0]!='C') continue;

    lwfm[ich]=0;

    qdc[ich]=0;

    chmax[ich][0]=-1; chmax[ich][1]=-1;

    chmin[ich][0]=1; chmin[ich][1]=-1;

    letime[ich][0]=-1.; letime[ich][1]=-1.;

    cftime[ich][0]=-1.; cftime[ich][1]=-1.;

    lwfm[ich]=1;

    nwfm[ich]+=1;

	t  = &tcdata[ich][0];

	y  = &chdat[ich][0];

	ay = &achdat[ich][0];

	float sum=0;

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

      sum+= dtdata[ich][(i+trgcell)%DWIDTH];

      if (version) t[i]=sum;

      else         t[i]=dtdata[0][i]; // old format

    }

    tccalib[ich]=t[(DWIDTH-trgcell)%DWIDTH];

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

      if (version)  t[i]-=(tccalib[ich]-tccalib[0]); // assume ich 0 is connected and always first in the data

      y[i]=-0.5+(float)ichdat[ich][i]/0xFFFF;

	  if (version) y[i]+=evrec.rc;

      h2_dpo[ich]->Fill(t[i],y[i]);

      gr_wfm[ich]->SetPoint(i,t[i],y[i]+0.05*(ich-1));

	}

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

	// accumulated waveforms

      ay[i]+=y[i];

	// signal minumum

      if (y[i]>chmax[ich][0]) {

        chmax[ich][0]=y[i];

        chmax[ich][1]=t[i];

      }

	// signal maximum

      if (y[i]<chmin[ich][0]) {

        chmin[ich][0]=y[i];

        chmin[ich][1]=t[i];

      }

      float t0  = t[i];

      float t0p = t[i+1];

      float t0n = t[i-1];

	// charge in the range

      if ((t0>adcgate[ich][0])&&(t0<adcgate[ich][1])) {

        qdc[ich]+=-y[i]*(t0p-t0n)/2;

      }

	}

	for (i=1;i<DWIDTH-1;i++) {

      float t0  = t[i];

      float t0p = t[i+1];

      float t0n = t[i-1];

    // leading edge in the range

      if ((t0>twin[ich][0])&&(t0<twin[ich][1])) {

        if ((y[i]<threshold[ich][1])&&(letime[ich][1]<0.)) {

           letime[ich][1]=t0n+(t0-t0n)*(threshold[ich][1]-y[i-1])/(y[i]-y[i-1]);

        }

	// constant fraction time in the range

        if ((chmin[ich][0]<threshold[ich][1])&&(y[i]<chmin[ich][0]*CFRAC)&&(cftime[ich][1]<0.)) {

            cftime[ich][1]=t0n+(t0-t0n)*(chmin[ich][0]*CFRAC-y[i-1])/(y[i]-y[i-1]);

        }

      }

    }

    if ((chmax[ich][1]>twin[ich][0])&&(chmax[ich][1]<twin[ich][1])) {

      h1_chmax[ich]->Fill(chmax[ich][0]);

      if (chmax[ich][0]>threshold[ich][0]) h1_chmaxt[ich]->Fill(chmax[ich][1]);

    }

    if ((chmin[ich][1]>twin[ich][0])&&(chmin[ich][1]<twin[ich][1])) {

      h1_chmin[ich]->Fill(chmin[ich][0]);

      if (chmin[ich][0]<threshold[ich][1]) h1_chmint[ich]->Fill(chmin[ich][1]);

    }

    h1_qdc[ich]->Fill(qdc[ich]);

    h1_letime[ich][1]->Fill(letime[ich][1]);

    h1_cftime[ich][1]->Fill(cftime[ich][1]);

    if (!(iev%updfrq)) {

      c_wfm->Update();

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

        ay[i]/=nwfm[ich];

        gr_wfm[ich]->SetPoint(i,t[i],ay[i]+0.05*(ich-1));

      }

      printf("processing event %d\n",iev);

    }

  }

  printf("number of events: %d\n",iev);

  fclose(fin);

  TCanvas* c_dpo= new TCanvas("c_dpo","DPO",750,50,700,700);

  c_dpo->Divide(2,2);

  for (i=0;i<4;i++) {c_dpo->cd(i+1);h2_dpo[i]->Draw("colz");}

  TCanvas* c_cftdif= new TCanvas("c_cftdif","DPO",750,50,700,700);

  c_cftdif->Divide(2,2);

  for (i=0;i<4;i++) {c_cftdif->cd(i+1);h1_cftdif[i]->Draw();}

  return;

}