Subversion Repositories f9daq

Rev

Rev 101 | Rev 180 | Go to most recent revision | Blame | Compare with Previous | Last modification | View Log | RSS feed

  1. #include "TH1.h"
  2. #include "TH2.h"
  3. #include "TGraph.h"
  4. #include "TCanvas.h"
  5. #include "TStyle.h"
  6. #include "TLegend.h"
  7. //#include "TROOT.h"
  8.  
  9. #include <stdio.h>
  10.  
  11. #define DATA_PATH  "./data/"
  12. #define CFRAC 0.4
  13. #define DWIDTH 1024
  14. void drs(int nev=1000, int updfrq=20, const char *FileName="drs000.dat", int trgch=0, int debug=0)
  15. {
  16.   int i,ich,iev;
  17.   int stat;
  18.   char histname[128];
  19. //  char FileName[300]="./data/n2.dat";
  20.   FILE *fin;
  21.  
  22.   char recid[5]="SAMO";
  23.   unsigned short *srecid = (unsigned short *) &recid[0];
  24.  
  25.   struct evrec {
  26.     unsigned long evn;
  27.     unsigned short tYear,tMonth,tDay,tHour,tMin,tSec,tmSec,rc;
  28.   } evrec;
  29.  
  30.   float dtdata[4][DWIDTH];
  31.   float tcdata[4][DWIDTH];
  32.   float tccalib[4];
  33.   float scaler[4];
  34.   float *t;
  35.   float *y;
  36.   float *ay;
  37.  
  38.   short lwfm[4];
  39.   int nwfm[4];
  40.   int board;
  41.   int version;
  42.  
  43.  
  44.   unsigned short ichdat[4][DWIDTH];
  45.   float chdat[4][DWIDTH];
  46.   float achdat[4][DWIDTH];
  47.   float chmax[4][2];
  48.   float chmin[4][2];
  49.   float letime[4][2],cftime[4][2],qdc[4];
  50. //  float threshold[4][2]={{0.,-0.02},{0.,0.},{0.,0.},{0.,-0.25}};
  51. //  float twin[4][2]={{120.,130.},{100.,150.},{100.,150.},{70.,90.}};
  52. //  float adcgate[4][2]={{124.,136.},{120.,140.},{120.,140.},{70.,90.}};
  53. //  float threshold[4][2]={{0.,-0.05},{0.,-0.05},{0.,-0.25},{0.,-0.25}};
  54. //  float twin[4][2]={{80.,100.},{80.,100.},{90.,110.},{50.,70.}};
  55. //  float adcgate[4][2]={{80.,110.},{80.,110.},{90.,110.},{50.,70.}};
  56.   float threshold[4][2]={{0.,-0.1  },{0.,-0.05 },{0.,-0.25 },{0.,-0.25}};
  57.   float twin[4][2]     ={{100.,120.},{100.,120.},{40.,130. },{60.,80. }};
  58.   float adcgate[4][2]  ={{105.,125.},{105.,125.},{110.,130.},{60.,80. }};
  59.   float vcut[4][2]     ={{-1,1     },{-1,1     },{-1,1     },{-1,1    }};
  60.  
  61. //  gROOT->Reset();
  62.  
  63.   gStyle->SetOptStat(1111111);
  64. //  gStyle->SetPalette(52);
  65.  
  66.   TH1F *h1_chmax[4];
  67.   TH1F *h1_chmaxt[4];
  68.   TH1F *h1_chmin[4];
  69.   TH1F *h1_chmint[4];
  70.   TH1F *h1_letime[4][2];
  71.   TH1F *h1_cftime[4][2];
  72.   TH1F *h1_cftdif[4];
  73.   TH1F *h1_qdc[4];
  74.   TH2F *h2_dpo[4];
  75.   TH2F *h2_tdcadc[4];
  76.   TH2F *h2_ctdcadc[4];
  77.   TH2F *h2_ctalk;
  78.   TGraph *gr_wfm[4];
  79.   TCanvas *c_wfm= new TCanvas("c_wfm","WFM",50,50,700,700);
  80.   TLegend *leg = new TLegend(0.8, 0.8, 1, 1);
  81.   c_wfm->Divide(1,1);
  82.   c_wfm->cd(1);
  83.  
  84.   for (i=0;i<4;i++) {
  85.     sprintf(histname,"chmax_%02d",i);
  86.     h1_chmax[i]=new TH1F(histname,histname,1100,-0.55,0.55);
  87.     sprintf(histname,"chmaxt_%02d",i);
  88.     h1_chmaxt[i]=new TH1F(histname,histname,DWIDTH,-0.1,204.7);
  89.     sprintf(histname,"chmin_%02d",i);
  90.     h1_chmin[i]=new TH1F(histname,histname,1100,-0.55,0.55);
  91.     sprintf(histname,"chmint_%02d",i);
  92.     h1_chmint[i]=new TH1F(histname,histname,DWIDTH,-0.1,204.7);
  93.     sprintf(histname,"letimep_%02d",i);
  94.     h1_letime[i][0]=new TH1F(histname,histname,500,twin[i][0],twin[i][1]);
  95.     sprintf(histname,"letimen_%02d",i);
  96.     h1_letime[i][1]=new TH1F(histname,histname,500,twin[i][0],twin[i][1]);
  97.     sprintf(histname,"cftimep_%02d",i);
  98.     h1_cftime[i][0]=new TH1F(histname,histname,500,twin[i][0],twin[i][1]);
  99.     sprintf(histname,"cftimen_%02d",i);
  100.     h1_cftime[i][1]=new TH1F(histname,histname,500,twin[i][0],twin[i][1]);
  101.     sprintf(histname,"cftdif_%02d",i);
  102.     h1_cftdif[i]=new TH1F(histname,histname,2000,twin[i][0],twin[i][1]);
  103.     sprintf(histname,"qdc_%02d",i);
  104.     h1_qdc[i]=new TH1F(histname,histname,500,-0.1,4.9);
  105.     sprintf(histname,"dpo_%02d",i);
  106.     h2_dpo[i]=new TH2F(histname,histname,DWIDTH,-0.1,204.7,1100,-0.55,0.55);
  107.     sprintf(histname,"tdcadc_%02d",i);
  108.     h2_tdcadc[i]=new TH2F(histname,histname,250,0.,0.5,400,twin[i][0],twin[i][1]);
  109.     sprintf(histname,"ctdcadc_%02d",i);
  110.     h2_ctdcadc[i]=new TH2F(histname,histname,250,0.,0.5,400,twin[i][0],twin[i][1]);
  111.     gr_wfm[i]=new TGraph(DWIDTH);
  112.     gr_wfm[i]->SetMinimum(-0.55);
  113.     gr_wfm[i]->SetMaximum(0.55);
  114.     gr_wfm[i]->SetLineColor(i+1);
  115.     gr_wfm[i]->SetMarkerColor(i+1);
  116.     sprintf(histname,"ch %d",i);
  117.     leg->AddEntry(gr_wfm[i], histname, "l");
  118.     if (i) gr_wfm[i]->Draw("LP"); else gr_wfm[i]->Draw("ALP");
  119.    
  120.   }
  121.   leg->Draw();
  122.   h2_ctalk=new TH2F("ctalk","ctalk",250,0.,0.5,250,0.,0.5);
  123.  
  124.   fin=fopen(FileName,"rb");
  125.   if (fin==NULL) {
  126.           printf("Error opening file %s\n",FileName);
  127.           return;
  128.   }      
  129.   printf("%s, %lu\n",recid,sizeof(evrec));
  130.  
  131.  
  132.   for (ich=0;ich<4;ich++){
  133.     nwfm[ich]=0;
  134.     for (i=0;i<DWIDTH;i++){
  135.       achdat[ich][i]=0;
  136.     }
  137.   }
  138.   unsigned int trgcell = 0;
  139.   int event = -1;
  140.  
  141.   iev=0;
  142.   while(!feof(fin)){
  143.     if (iev>=nev) break;
  144.     stat=fread(recid,1,4,fin);
  145.     if (debug){
  146.       if (recid[1]=='#') printf("%c%c%u\n",recid[0],recid[1],srecid[1]);
  147.       else  printf("%s\n",recid );
  148.     }
  149.     switch (recid[0]){
  150.            case 'D': // DRS
  151.            version = atoi(&recid[3]);
  152.            break;
  153.        case 'T': // TIME
  154.            switch (recid[1]){
  155.              case 'I':
  156.                            event=0;
  157.                if (version ==0) version++;                       
  158.                            break;
  159.              case '#': trgcell = srecid[1]; break;
  160.            }
  161.            break;
  162.        case 'B':
  163.            board = srecid[1];
  164.            break;
  165.        case 'C':
  166.            sscanf(recid,"C%d",&ich);
  167.            ich=ich-1;
  168.            if (event) {
  169.                          if (version>1) stat=fread(&scaler[ich],1,sizeof(float),fin);
  170.              stat=fread(ichdat[ich],1,sizeof(ichdat[ich]),fin);
  171.            } else {
  172.                          
  173.              stat=fread(dtdata[ich],1,sizeof(float)*DWIDTH,fin);
  174.              
  175.            }
  176.            break;
  177.        case 'E':
  178.            stat=fread(&evrec.evn,1,sizeof(evrec),fin);
  179.            iev++;
  180.                    if (!version)  {
  181.                          stat=fread(dtdata[0],1,sizeof(float)*DWIDTH,fin);
  182.                    }  
  183.            event=1;
  184.            break;
  185.        default:
  186.            printf("Unknown header! %s\n", recid);
  187.            break;
  188.     }
  189.  
  190.  
  191.     if (recid[0] == 'E'){
  192.          for ( i=0;i<4;i++){
  193.            if (!lwfm[i]) continue;
  194.            if ((chmin[i][0]>vcut[i][0])&&(chmin[i][0]<vcut[i][1])) {
  195.              h1_cftdif[i]->Fill(cftime[i][1]-cftime[trgch][1]+110.);
  196.            }
  197. //      h2_tdcadc[i]->Fill(-chmin[i][0],letime[i][1]-cftime[trgch][1]+72);
  198. //      h2_ctdcadc[i]->Fill(-chmin[i][0],cftime[i][1]-cftime[trgch][1]+72);
  199.            h2_tdcadc[i]->Fill(-chmin[0][0],letime[i][1]-cftime[trgch][1]+72);
  200.            h2_ctdcadc[i]->Fill(-chmin[0][0],cftime[i][1]-cftime[trgch][1]+72);
  201.          }
  202.          h2_ctalk->Fill(-chmin[0][0],-chmin[1][0]);
  203.     }
  204.     if (event<1) continue;
  205.     if (recid[0]!='C') continue;
  206.  
  207.    
  208.     lwfm[ich]=0;
  209.     qdc[ich]=0;
  210.     chmax[ich][0]=-1; chmax[ich][1]=-1;
  211.     chmin[ich][0]=1; chmin[ich][1]=-1;
  212.     letime[ich][0]=-1.; letime[ich][1]=-1.;
  213.     cftime[ich][0]=-1.; cftime[ich][1]=-1.;
  214.    
  215.  
  216.     lwfm[ich]=1;
  217.     nwfm[ich]+=1;
  218.  
  219.      
  220.          
  221.         t  = &tcdata[ich][0];
  222.         y  = &chdat[ich][0];
  223.         ay = &achdat[ich][0];
  224.  
  225.         float sum=0;
  226.     for (i=0;i<DWIDTH;i++) {
  227.       sum+= dtdata[ich][(i+trgcell)%DWIDTH];
  228.       if (version) t[i]=sum;
  229.       else         t[i]=dtdata[0][i]; // old format        
  230.     }
  231.     tccalib[ich]=t[(DWIDTH-trgcell)%DWIDTH];
  232.  
  233.     for (i=0;i<DWIDTH;i++) {
  234.       if (version)  t[i]-=(tccalib[ich]-tccalib[0]); // assume ich 0 is connected and always first in the data          
  235.       y[i]=-0.5+(float)ichdat[ich][i]/0xFFFF;
  236.           if (version) y[i]+=evrec.rc;
  237.       h2_dpo[ich]->Fill(t[i],y[i]);
  238.       gr_wfm[ich]->SetPoint(i,t[i],y[i]+0.05*(ich-1));
  239.         }
  240.          
  241.          
  242.         for (i=0;i<DWIDTH;i++) {       
  243.         // accumulated waveforms
  244.       ay[i]+=y[i];
  245.                
  246.         // signal minumum
  247.       if (y[i]>chmax[ich][0]) {
  248.         chmax[ich][0]=y[i];
  249.         chmax[ich][1]=t[i];
  250.       }
  251.                
  252.         // signal maximum
  253.       if (y[i]<chmin[ich][0]) {
  254.         chmin[ich][0]=y[i];
  255.         chmin[ich][1]=t[i];
  256.       }
  257.       float t0  = t[i];
  258.       float t0p = t[i+1];
  259.       float t0n = t[i-1];
  260.         // charge in the range
  261.       if ((t0>adcgate[ich][0])&&(t0<adcgate[ich][1])) {
  262.         qdc[ich]+=-y[i]*(t0p-t0n)/2;
  263.       }
  264.         }
  265.        
  266.         for (i=1;i<DWIDTH-1;i++) {
  267.       float t0  = t[i];
  268.       float t0p = t[i+1];
  269.       float t0n = t[i-1];              
  270.     // leading edge in the range
  271.       if ((t0>twin[ich][0])&&(t0<twin[ich][1])) {
  272.         if ((y[i]<threshold[ich][1])&&(letime[ich][1]<0.)) {
  273.            letime[ich][1]=t0n+(t0-t0n)*(threshold[ich][1]-y[i-1])/(y[i]-y[i-1]);                          
  274.         }
  275.         // constant fraction time in the range
  276.         if ((chmin[ich][0]<threshold[ich][1])&&(y[i]<chmin[ich][0]*CFRAC)&&(cftime[ich][1]<0.)) {
  277.             cftime[ich][1]=t0n+(t0-t0n)*(chmin[ich][0]*CFRAC-y[i-1])/(y[i]-y[i-1]);                          
  278.         }
  279.       }
  280.     }
  281.      
  282.     if ((chmax[ich][1]>twin[ich][0])&&(chmax[ich][1]<twin[ich][1])) {
  283.       h1_chmax[ich]->Fill(chmax[ich][0]);
  284.       if (chmax[ich][0]>threshold[ich][0]) h1_chmaxt[ich]->Fill(chmax[ich][1]);
  285.     }
  286.     if ((chmin[ich][1]>twin[ich][0])&&(chmin[ich][1]<twin[ich][1])) {
  287.       h1_chmin[ich]->Fill(chmin[ich][0]);
  288.       if (chmin[ich][0]<threshold[ich][1]) h1_chmint[ich]->Fill(chmin[ich][1]);
  289.     }
  290.     h1_qdc[ich]->Fill(qdc[ich]);
  291.     h1_letime[ich][1]->Fill(letime[ich][1]);
  292.     h1_cftime[ich][1]->Fill(cftime[ich][1]);
  293.  
  294.     if (!(iev%updfrq)) {
  295.       c_wfm->Update();
  296.       for (i=0;i<DWIDTH;i++) {
  297.         ay[i]/=nwfm[ich];
  298.         gr_wfm[ich]->SetPoint(i,t[i],ay[i]+0.05*(ich-1));
  299.       }
  300.       printf("processing event %d\n",iev);
  301.     }
  302.   }
  303.  
  304.   printf("number of events: %d\n",iev);
  305.   fclose(fin);
  306.  
  307.   TCanvas* c_dpo= new TCanvas("c_dpo","DPO",750,50,700,700);
  308.   c_dpo->Divide(2,2);
  309.   for (i=0;i<4;i++) {c_dpo->cd(i+1);h2_dpo[i]->Draw("colz");}
  310.  
  311.   TCanvas* c_cftdif= new TCanvas("c_cftdif","DPO",750,50,700,700);
  312.   c_cftdif->Divide(2,2);
  313.   for (i=0;i<4;i++) {c_cftdif->cd(i+1);h1_cftdif[i]->Draw();}
  314.   return;
  315. }
  316.