adding drs4 raw and analysis files

convertDRS4/README.txt

 Application for convertation of binary files acquired by DRS4 to ROOT file.
+
+1. g++ -o read_binary read_binary.cpp -lm `root-config --cflags --libs`
+This compiles read_binary.cpp file and makes object file read_binary.
+2. ./read_binary ./data/rawDataDSR4/file_name.dat
+With the help of the object file it is possible to convert binary file .dat to file .root (which we call raw /home/dariak/NeuRad_tests/data/rawDataDSR4)
+3. read_root.C 
+This script reads raw file and make another file .root with simple analysis /home/dariak/NeuRad_tests/data/dataDSR4

convertDRS4/read_binary.cpp

+/*
+   Name:           read_binary.cpp
+   Created by:     Stefan Ritt <stefan.ritt@psi.ch>
+   Date:           July 30th, 2014
+
+   Purpose:        Example file to read binary data saved by DRSOsc.
+ 
+   Compile and run it with:
+ 
+      gcc -o read_binary read_binary.cpp
+ 
+      ./read_binary <filename>
+
+   This program assumes that a pulse from a signal generator is split
+   and fed into channels #1 and #2. It then calculates the time difference
+   between these two pulses to show the performance of the DRS board
+   for time measurements.
+
+   $Id: read_binary.cpp 22290 2016-04-27 14:51:37Z ritt $
+*/
+
+#include <stdio.h>
+#include <fcntl.h>
+#include <unistd.h>
+#include <string.h>
+#include <math.h>
+
+//for ROOT
+#include "TROOT.h"
+#include "TFile.h"
+#include "TTree.h"
+#include "TH1F.h"
+
+typedef struct {
+   char           tag[3];
+   char           version;
+} FHEADER;
+
+typedef struct {
+   char           time_header[4];
+} THEADER;
+
+typedef struct {
+   char           bn[2];
+   unsigned short board_serial_number;
+} BHEADER;
+
+typedef struct {
+   char           event_header[4];
+   unsigned int   event_serial_number;
+   unsigned short year;
+   unsigned short month;
+   unsigned short day;
+   unsigned short hour;
+   unsigned short minute;
+   unsigned short second;
+   unsigned short millisecond;
+   unsigned short range;	// range center in mV
+} EHEADER;
+
+typedef struct {
+   char           tc[2];
+   unsigned short trigger_cell;
+} TCHEADER;
+
+typedef struct {
+   char           c[1];
+   char           cn[3];
+} CHEADER;
+
+/*-----------------------------------------------------------------------------*/
+
+int main(int argc, const char * argv[])
+{
+   FHEADER  fh;
+   THEADER  th;
+   BHEADER  bh;
+   EHEADER  eh;
+   TCHEADER tch;
+   CHEADER  ch;
+   
+   unsigned int scaler;
+   unsigned short voltage[1024];
+   double waveform[16][4][1024], time[16][4][1024];
+   float bin_width[16][4][1024];
+   int i, j, b, chn, n, chn_index, n_boards;
+   double t1, t2, t3, t4, dt, dt34;
+   char filename[256];
+
+   int ndt;
+   double threshold, sumdt, sumdt2;
+
+
+//for ROOT
+
+   TFile* rfile = new TFile("/home/dariak/NeuRad_tests/data/rawDataDSR4/NeuRad_test_08_2.root", "RECREATE");
+   TTree* rtree = new TTree("rtree", "tree for drs4 analysis");
+   //rtree->Branch("t1", &t1, "t1/D"); //br for time of threshold crossing signal in 1 ch
+  // rtree->Branch("t2", &t2, "t2/D"); //br for time of threshold crossing signal in 2 ch
+   int ncell;
+   const int ncellMax = 1030;
+   double amp_ch1[ncellMax], time_ch1[ncellMax];	//variable size array 
+//------for other channels 
+//   double amp_ch2[ncellMax], time_ch2[ncellMax];
+
+   rtree->Branch("ncell", &ncell, "ncell/I"); 
+   rtree->Branch("amp_ch1", amp_ch1, "amp_ch1[ncell]/D"); 
+   rtree->Branch("time_ch1", time_ch1, "time_ch1[ncell]/D"); 
+
+//------for other channels 
+//   rtree->Branch("amp_ch2", amp_ch2, "amp_ch2[ncell]/D"); 
+//  rtree->Branch("time_ch2", time_ch2, "time_ch2[ncell]/D"); 
+
+   if (argc > 1)
+      strcpy(filename, argv[1]);
+   else {
+      printf("Usage: read_binary <filename>\n");
+      return 0;
+   }
+   
+   // open the binary waveform file
+   FILE *f = fopen(filename, "r");
+   if (f == NULL) {
+      printf("Cannot find file \'%s\'\n", filename);
+      return 0;
+   }
+
+   // read file header
+   fread(&fh, sizeof(fh), 1, f);
+   if (fh.tag[0] != 'D' || fh.tag[1] != 'R' || fh.tag[2] != 'S') {
+      printf("Found invalid file header in file \'%s\', aborting.\n", filename);
+      return 0;
+   }
+   
+   if (fh.version != '2') {
+      printf("Found invalid file version \'%c\' in file \'%s\', should be \'2\', aborting.\n", fh.version, filename);
+      return 0;
+   }
+
+   // read time header
+   fread(&th, sizeof(th), 1, f);
+   if (memcmp(th.time_header, "TIME", 4) != 0) {
+      printf("Invalid time header in file \'%s\', aborting.\n", filename);
+      return 0;
+   }
+
+   for (b = 0 ; ; b++) {
+      // read board header
+      fread(&bh, sizeof(bh), 1, f);
+      if (memcmp(bh.bn, "B#", 2) != 0) {
+         // probably event header found
+         fseek(f, -4, SEEK_CUR);
+         break;
+      }
+      
+      printf("Found data for board #%d\n", bh.board_serial_number);
+      
+      // read time bin widths
+      memset(bin_width[b], sizeof(bin_width[0]), 0);
+      for (chn=0 ; chn<5 ; chn++) {
+         fread(&ch, sizeof(ch), 1, f);
+         if (ch.c[0] != 'C') {
+            // event header found
+            fseek(f, -4, SEEK_CUR);
+            break;
+         }
+         i = ch.cn[2] - '0' - 1;
+         printf("Found timing calibration for channel #%d\n", i+1);
+         fread(&bin_width[b][i][0], sizeof(float), 1024, f);
+	/*my printf
+		printf("bin width %d \n", bin_width[b][i][10]); */
+         // fix for 2048 bin mode: double channel
+         if (bin_width[b][i][1023] > 10 || bin_width[b][i][1023] < 0.01) {
+            for (j=0 ; j<512 ; j++)
+               bin_width[b][i][j+512] = bin_width[b][i][j];
+	/*my printf
+		printf("bin width %d \n", bin_width[b][i][j+512]); */
+         }
+      }
+   }
+   n_boards = b;
+   
+   // initialize statistics
+   ndt = 0;
+   sumdt = sumdt2 = 0;
+   
+   // loop over all events in the data file
+   for (n=0 ; ; n++) {
+      // read event header
+      i = (int)fread(&eh, sizeof(eh), 1, f);
+      if (i < 1)
+         break;
+      
+      if ( !(eh.event_serial_number%100) ) {
+	 printf("Found event #%d\n", eh.event_serial_number);
+      }
+      
+      // loop over all boards in data file
+      for (b=0 ; b<n_boards ; b++) {
+         
+         // read board header
+         fread(&bh, sizeof(bh), 1, f);
+         if (memcmp(bh.bn, "B#", 2) != 0) {
+            printf("Invalid board header in file \'%s\', aborting.\n", filename);
+            return 0;
+         }
+         
+         // read trigger cell
+         fread(&tch, sizeof(tch), 1, f);
+         if (memcmp(tch.tc, "T#", 2) != 0) {
+            printf("Invalid trigger cell header in file \'%s\', aborting.\n", filename);
+            return 0;
+         }
+
+         if (n_boards > 1)
+            printf("Found data for board #%d\n", bh.board_serial_number);
+         
+         // reach channel data
+         for (chn=0 ; chn<4 ; chn++) {
+            
+            // read channel header
+            fread(&ch, sizeof(ch), 1, f);
+            if (ch.c[0] != 'C') {
+               // event header found
+               fseek(f, -4, SEEK_CUR);
+               break;
+            }
+            chn_index = ch.cn[2] - '0' - 1;
+	//	printf("print channel %d \n",chn_index);
+            fread(&scaler, sizeof(int), 1, f);
+            fread(voltage, sizeof(short), 1024, f);
+            
+            for (i=0 ; i<1024 ; i++) {
+               // convert data to volts
+               waveform[b][chn_index][i] = (voltage[i] / 65536. + eh.range/1000.0 - 0.5); //calculation of amplitudes values for each cell
+	
+		//for ROOT
+	       ncell = i;
+	       if(chn_index == 0) {amp_ch1[i] = waveform[b][chn_index][i];}
+	     //  if(chn_index == 1) {amp_ch2[i] = waveform[b][chn_index][i];}
+
+               // calculate time for this cell
+               for (j=0,time[b][chn_index][i]=0 ; j<i ; j++){
+                  time[b][chn_index][i] += bin_width[b][chn_index][(j+tch.trigger_cell) % 1024];
+	       }
+            }	    
+         } // end of the channel loop (chn)
+         
+         // align cell #0 of all channels
+         t1 = time[b][0][(1024-tch.trigger_cell) % 1024];
+	//my print;
+	// printf("t1 %1.6lf \n",time[b][0][(1024-tch.trigger_cell) % 1024]);
+         for (chn=1 ; chn<4 ; chn++) {
+            t2 = time[b][chn][(1024-tch.trigger_cell) % 1024];
+//adding channels 3 and 4
+            t3 = time[b][chn][(1024-tch.trigger_cell) % 1024];
+            t4 = time[b][chn][(1024-tch.trigger_cell) % 1024];
+	//my prinf
+	//printf("t2 %f for %d %d %d \n",time[b][chn][(1024-tch.trigger_cell) % 1024], b, chn, i);
+            dt = t1 - t2;
+            dt34 = t3 - t4;
+            for (i=0 ; i<1024 ; i++) {
+               time[b][chn][i] += dt;  //each element of time gets dt correction
+		//my print;
+	      // printf("time %1.6lf for %d %d %d \n",time[b][chn][i], b, chn, i);
+	    }
+		
+         }
+         
+         t1 = t2 = t3 = t4 = 0;
+         threshold = -0.045; //my threshold, used to be 0.3
+
+
+	    //for ROOT
+	 for(i=0 ; i<1024 ; i++) {
+	    time_ch1[i] = time[b][0][i];
+	   // time_ch2[i] = time[b][1][i];
+	 }
+
+             // find peak in channel 1 above threshold
+         for (i=0 ; i<1022 ; i++) {
+
+            if (waveform[b][0][i] < threshold && waveform[b][0][i+1] >= threshold) {
+               t1 = (threshold-waveform[b][0][i])/(waveform[b][0][i+1]-waveform[b][0][i])*(time[b][0][i+1]-time[b][0][i])+time[b][0][i];
+		//my prinf
+		//printf("t1 recalc %1.6lf %d \n",t1, i);
+               break;
+            }
+
+	 }
+         
+         // find peak in channel 2 above threshold
+         for (i=0 ; i<1022 ; i++) {
+            if (waveform[b][1][i] < threshold && waveform[b][1][i+1] >= threshold) {
+               t2 = (threshold-waveform[b][1][i])/(waveform[b][1][i+1]-waveform[b][1][i])*(time[b][1][i+1]-time[b][1][i])+time[b][1][i];
+		//my prinf
+		//printf("t2 recalc %1.6lf \n",t2);
+               break;
+            }
+	 }
+
+         // calculate distance of peaks with statistics
+         if (t1 > 0 && t2 > 0) {
+            ndt++;
+            dt = t2 - t1;
+            sumdt += dt;
+            sumdt2 += dt*dt;
+         }
+      } //end of the boards loop
+
+      rtree->Fill();
+   } // end of the events loop
+   
+   // print statistics
+   printf("dT = %1.3lfns +- %1.1lfps\n", sumdt/ndt, 1000*sqrt(1.0/(ndt-1)*(sumdt2-1.0/ndt*sumdt*sumdt)));
+
+   rfile->Write();
+   rfile->Close();
+   return 1;
+}
+

convertDRS4/read_root.C

+void read_root() 
+{
+	TFile *f = new TFile("/home/dariak/NeuRad_tests/data/rawDataDSR4/NeuRad_test_07_1.root");
+	TTree *tr = (TTree*)f->Get("rtree");
+	tr->SetMakeClass(1);
+
+//	TArrayD ampl1;
+//	TArrayD time1;
+//	ampl1.Set(ncellMax);
+//	time1.Set(ncellMax);
+
+	const int ncellMax = 1030;
+	Double_t ampl1[ncellMax];
+	Double_t time1[ncellMax];
+
+	TFile *fw = new TFile("/home/dariak/NeuRad_tests/data/dataDSR4/analysis_07_1.root", "RECREATE");
+	TTree *tw = new TTree("drs4analysis", "title of drs4 analysis tree");
+	Double_t minAmplitude1;
+	tw->Branch("minAmplitude1", &minAmplitude1, "minAmplitude1/D");
+/*	Double_t ampl1[1024];
+	Double_t time1[1024];*/
+
+	tr->SetBranchAddress("amp_ch1", ampl1);
+	tr->SetBranchAddress("time_ch1", time1);
+
+	const Long64_t nEntries = tr->GetEntries();
+	cout <<"Number of entries: "<<nEntries<<endl;
+
+//----for pictures 
+/*	TCanvas *c1 = new TCanvas();
+	c1->Divide(3,4);
+
+	for(Int_t i=12; i<24; i++) {
+		c1->cd(i-11);
+		tr->GetEntry(i);
+		TGraph *gr1 = new TGraph(1023, time1, ampl1);
+		gr1->SetTitle("Signal shape for one event");
+		gr1->GetXaxis();
+		gr1->Draw("Al*");
+	}
+
+	tr->GetEntry(428);
+	TGraph *gr1 = new TGraph(1023, time1, ampl1);
+	gr1->SetTitle("Signal shape for one event");
+	gr1->Draw("Al*");*/
+//-----
+	Double_t minAmpl1 = 0.; 
+	
+	for(Int_t i=0; i<tr->GetEntries(); i++) {
+		tr->GetEntry(i);
+		minAmpl1 = ampl1[0];
+
+		//find minimum by hand 
+		for(Int_t j=0; j<1023; j++) {	
+			//cout<<ampl1[j]<<" "<<j<<endl;
+			if(ampl1[j] < minAmpl1) {
+				minAmpl1 = ampl1[j];
+			}
+			/*cout<<"Time "<<i<<" "<<j<<" "<<time1[j]<<endl;
+			//cout<<"Amplitude "<<i<<" "<<j<<" "<<ampl1[j]<<endl;
+			cout<<endl;*/
+		}
+		//fitting 90 percent of rising edge
+		for(Int_t k=0; k<1023; k++) {
+			if(ampl1[k] < 0.1*minAmpl1 && ampl1[k] > 0.9*minAmpl1) { //we have negative signals, amplitude between 10 and 90 proc from minimum
+				/*cout<<"rjgnfre"<<endl;
+				TGraph *gr2 = new TGraph(1023, time1, ampl1);
+				gr2->SetTitle("stupido");
+				gr2->Fit("pol1");
+				gr2->Draw("Al*");*/
+				
+			}
+		}
+		//getting integral
+
+		TGraph *gr3 = new TGraph(1023, time1, ampl1);
+		//for(Int_t k=120; k<180; k++) {
+		gr3->Integral(120,180);
+		cout<<"INtegral "<<gr3->Integral(120,180)<<endl;
+		//}
+		minAmplitude1 = minAmpl1;
+		cout<<"Min amplitude "<<minAmpl1<<" for entry "<<i<<endl;
+		
+
+		tw->Fill();
+	}
+	fw->cd();
+	tw->Write();
+	fw->Close();
+}