Skip to content

N
NeuRad_tests
Commit bf84b293 authored by Vratislav Chudoba's avatar Vratislav Chudoba
Browse files
Options

Branch with class RawEvent added in read_binary_DRS application.

parent 7b3acf5f
Hide whitespace changes
Inline Side-by-side
Showing with 303 additions and 280 deletions
convertDRS4/read_binary.cpp
View file @ bf84b293
...@@ -4,11 +4,11 @@ ...@@ -4,11 +4,11 @@
Date: July 30th, 2014 Date: July 30th, 2014
Purpose: Example file to read binary data saved by DRSOsc. Purpose: Example file to read binary data saved by DRSOsc.
Compile and run it with: Compile and run it with:
gcc -o read_binary read_binary.cpp gcc -o read_binary read_binary.cpp
./read_binary <filename> ./read_binary <filename>
This program assumes that a pulse from a signal generator is split This program assumes that a pulse from a signal generator is split
...@@ -17,7 +17,7 @@ ...@@ -17,7 +17,7 @@
for time measurements. for time measurements.
$Id: read_binary.cpp 22290 2016-04-27 14:51:37Z ritt $ $Id: read_binary.cpp 22290 2016-04-27 14:51:37Z ritt $
*/ */
#include <stdio.h> #include <stdio.h>
#include <fcntl.h> #include <fcntl.h>
...@@ -31,299 +31,311 @@ ...@@ -31,299 +31,311 @@
#include "TTree.h" #include "TTree.h"
#include "TH1F.h" #include "TH1F.h"
//our code
#include "../dataClasses/RawData.h"
typedef struct { typedef struct {
char tag[3]; char tag[3];
char version; char version;
} FHEADER; } FHEADER;
typedef struct { typedef struct {
char time_header[4]; char time_header[4];
} THEADER; } THEADER;
typedef struct { typedef struct {
char bn[2]; char bn[2];
unsigned short board_serial_number; unsigned short board_serial_number;
} BHEADER; } BHEADER;
typedef struct { typedef struct {
char event_header[4]; char event_header[4];
unsigned int event_serial_number; unsigned int event_serial_number;
unsigned short year; unsigned short year;
unsigned short month; unsigned short month;
unsigned short day; unsigned short day;
unsigned short hour; unsigned short hour;
unsigned short minute; unsigned short minute;
unsigned short second; unsigned short second;
unsigned short millisecond; unsigned short millisecond;
unsigned short range; // range center in mV unsigned short range; // range center in mV
} EHEADER; } EHEADER;
typedef struct { typedef struct {
char tc[2]; char tc[2];
unsigned short trigger_cell; unsigned short trigger_cell;
} TCHEADER; } TCHEADER;
typedef struct { typedef struct {
char c[1]; char c[1];
char cn[3]; char cn[3];
} CHEADER; } CHEADER;
/*-----------------------------------------------------------------------------*/ /*-----------------------------------------------------------------------------*/
int main(int argc, const char * argv[]) int main(int argc, const char * argv[])
{ {
FHEADER fh; FHEADER fh;
THEADER th; THEADER th;
BHEADER bh; BHEADER bh;
EHEADER eh; EHEADER eh;
TCHEADER tch; TCHEADER tch;
CHEADER ch; CHEADER ch;
unsigned int scaler; unsigned int scaler;
unsigned short voltage[1024]; unsigned short voltage[1024];
double waveform[16][4][1024], time[16][4][1024]; double waveform[16][4][1024], time[16][4][1024];
float bin_width[16][4][1024]; float bin_width[16][4][1024];
int i, j, b, chn, n, chn_index, n_boards; int i, j, b, chn, n, chn_index, n_boards;
double t1, t2, t3, t4, dt, dt34; double t1, t2, t3, t4, dt, dt34;
char filename[256]; //for input binary file char filename[256]; //for input binary file
char outroot[256]; //for output root file char outroot[256]; //for output root file
int ndt; int ndt;
double threshold, sumdt, sumdt2; double threshold, sumdt, sumdt2;
if (argc == 3) { if (argc == 3) {
strcpy(filename, argv[1]); strcpy(filename, argv[1]);
strcpy(outroot, argv[2]); strcpy(outroot, argv[2]);
} }
else if (argc == 2) { else if (argc == 2) {
printf("Error: both input binary file and output root file should be specified!\n"); printf("Error: both input binary file and output root file should be specified!\n");
return -1; return -1;
} }
else { else {
printf("Error: input binary file and output root file should be specified!\n"); printf("Error: input binary file and output root file should be specified!\n");
return -1; return -1;
} }
// ---------------for ROOT // ---------------for ROOT
TFile* rfile = new TFile(outroot, "RECREATE"); TFile* rfile = new TFile(outroot, "RECREATE");
TTree* rtree = new TTree("rtree", "tree for drs4 analysis"); 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("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 // rtree->Branch("t2", &t2, "t2/D"); //br for time of threshold crossing signal in 2 ch
int ncell; int ncell;
const int ncellMax = 1030; const int ncellMax = 1030;
double amp_ch1[ncellMax], time_ch1[ncellMax]; //variable size array double amp_ch1[ncellMax], time_ch1[ncellMax]; //variable size array
//------for other channels //------for other channels
// double amp_ch2[ncellMax], time_ch2[ncellMax]; // double amp_ch2[ncellMax], time_ch2[ncellMax];
rtree->Branch("ncell", &ncell, "ncell/I"); rtree->Branch("ncell", &ncell, "ncell/I");
rtree->Branch("amp_ch1", amp_ch1, "amp_ch1[ncell]/D"); rtree->Branch("amp_ch1", amp_ch1, "amp_ch1[ncell]/D");
rtree->Branch("time_ch1", time_ch1, "time_ch1[ncell]/D"); rtree->Branch("time_ch1", time_ch1, "time_ch1[ncell]/D");
//------for other channels RawData *event = new RawData();
// rtree->Branch("amp_ch2", amp_ch2, "amp_ch2[ncell]/D"); rtree->Bronch("rawEvent", "RawData", &event);
// rtree->Branch("time_ch2", time_ch2, "time_ch2[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");
// open the binary waveform file
FILE *f = fopen(filename, "r"); //----------------
if (f == NULL) {
printf("Cannot find file \'%s\'\n", filename); // open the binary waveform file
return 0; FILE *f = fopen(filename, "r");
} if (f == NULL) {
printf("Cannot find file \'%s\'\n", filename);
// read file header return 0;
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); // read file header
return 0; 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);
if (fh.version != '2') { return 0;
printf("Found invalid file version \'%c\' in file \'%s\', should be \'2\', aborting.\n", fh.version, filename); }
return 0;
} if (fh.version != '2') {
printf("Found invalid file version \'%c\' in file \'%s\', should be \'2\', aborting.\n", fh.version, filename);
// read time header return 0;
fread(&th, sizeof(th), 1, f); }
if (memcmp(th.time_header, "TIME", 4) != 0) {
printf("Invalid time header in file \'%s\', aborting.\n", filename); // read time header
return 0; fread(&th, sizeof(th), 1, f);
} if (memcmp(th.time_header, "TIME", 4) != 0) {
printf("Invalid time header in file \'%s\', aborting.\n", filename);
for (b = 0 ; ; b++) { return 0;
// read board header }
fread(&bh, sizeof(bh), 1, f);
if (memcmp(bh.bn, "B#", 2) != 0) { for (b = 0 ; ; b++) {
// probably event header found // read board header
fseek(f, -4, SEEK_CUR); fread(&bh, sizeof(bh), 1, f);
break; if (memcmp(bh.bn, "B#", 2) != 0) {
} // probably event header found
fseek(f, -4, SEEK_CUR);
printf("Found data for board #%d\n", bh.board_serial_number); break;
}
// read time bin widths
memset(bin_width[b], sizeof(bin_width[0]), 0); printf("Found data for board #%d\n", bh.board_serial_number);
for (chn=0 ; chn<5 ; chn++) {
fread(&ch, sizeof(ch), 1, f); // read time bin widths
if (ch.c[0] != 'C') { memset(bin_width[b], sizeof(bin_width[0]), 0);
// event header found for (chn=0 ; chn<5 ; chn++) {
fseek(f, -4, SEEK_CUR); fread(&ch, sizeof(ch), 1, f);
break; if (ch.c[0] != 'C') {
} // event header found
i = ch.cn[2] - '0' - 1; fseek(f, -4, SEEK_CUR);
printf("Found timing calibration for channel #%d\n", i+1); break;
fread(&bin_width[b][i][0], sizeof(float), 1024, f); }
/*my printf 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]); */ printf("bin width %d \n", bin_width[b][i][10]); */
// fix for 2048 bin mode: double channel // fix for 2048 bin mode: double channel
if (bin_width[b][i][1023] > 10 || bin_width[b][i][1023] < 0.01) { if (bin_width[b][i][1023] > 10 || bin_width[b][i][1023] < 0.01) {
for (j=0 ; j<512 ; j++) for (j=0 ; j<512 ; j++)
bin_width[b][i][j+512] = bin_width[b][i][j]; bin_width[b][i][j+512] = bin_width[b][i][j];
/*my printf /*my printf
printf("bin width %d \n", bin_width[b][i][j+512]); */ printf("bin width %d \n", bin_width[b][i][j+512]); */
} }
} }
} }
n_boards = b; n_boards = b;
// initialize statistics // initialize statistics
ndt = 0; ndt = 0;
sumdt = sumdt2 = 0; sumdt = sumdt2 = 0;
// loop over all events in the data file // loop over all events in the data file
for (n=0 ; ; n++) { for (n=0 ; ; n++) {
// read event header // read event header
i = (int)fread(&eh, sizeof(eh), 1, f); i = (int)fread(&eh, sizeof(eh), 1, f);
if (i < 1) if (i < 1)
break; break;
if ( !(eh.event_serial_number%100) ) { if ( !(eh.event_serial_number%100) ) {
printf("Found event #%d\n", eh.event_serial_number); printf("Found event #%d\n", eh.event_serial_number);
} }
// loop over all boards in data file // loop over all boards in data file
for (b=0 ; b<n_boards ; b++) { for (b=0 ; b<n_boards ; b++) {
// read board header // read board header
fread(&bh, sizeof(bh), 1, f); fread(&bh, sizeof(bh), 1, f);
if (memcmp(bh.bn, "B#", 2) != 0) { if (memcmp(bh.bn, "B#", 2) != 0) {
printf("Invalid board header in file \'%s\', aborting.\n", filename); printf("Invalid board header in file \'%s\', aborting.\n", filename);
return 0; return 0;
} }
// read trigger cell // read trigger cell
fread(&tch, sizeof(tch), 1, f); fread(&tch, sizeof(tch), 1, f);
if (memcmp(tch.tc, "T#", 2) != 0) { if (memcmp(tch.tc, "T#", 2) != 0) {
printf("Invalid trigger cell header in file \'%s\', aborting.\n", filename); printf("Invalid trigger cell header in file \'%s\', aborting.\n", filename);
return 0; return 0;
} }
if (n_boards > 1) if (n_boards > 1)
printf("Found data for board #%d\n", bh.board_serial_number); printf("Found data for board #%d\n", bh.board_serial_number);
// reach channel data // reach channel data
for (chn=0 ; chn<4 ; chn++) { for (chn=0 ; chn<4 ; chn++) {
// read channel header // read channel header
fread(&ch, sizeof(ch), 1, f); fread(&ch, sizeof(ch), 1, f);
if (ch.c[0] != 'C') { if (ch.c[0] != 'C') {
// event header found // event header found
fseek(f, -4, SEEK_CUR); fseek(f, -4, SEEK_CUR);
break; break;
} }
chn_index = ch.cn[2] - '0' - 1; chn_index = ch.cn[2] - '0' - 1;
// printf("print channel %d \n",chn_index); // printf("print channel %d \n",chn_index);
fread(&scaler, sizeof(int), 1, f); fread(&scaler, sizeof(int), 1, f);
fread(voltage, sizeof(short), 1024, f); fread(voltage, sizeof(short), 1024, f);
for (i=0 ; i<1024 ; i++) { for (i=0 ; i<1024 ; i++) {
// convert data to volts // 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 waveform[b][chn_index][i] = (voltage[i] / 65536. + eh.range/1000.0 - 0.5); //calculation of amplitudes values for each cell
//for ROOT //for ROOT
ncell = i; ncell = i;
if(chn_index == 0) {amp_ch1[i] = waveform[b][chn_index][i];} if(chn_index == 0) {
// if(chn_index == 1) {amp_ch2[i] = waveform[b][chn_index][i];} amp_ch1[i] = waveform[b][chn_index][i];
event->SetAmp(waveform[b][chn_index][i], i);
// calculate time for this cell }
for (j=0,time[b][chn_index][i]=0 ; j<i ; j++){ // if(chn_index == 1) {amp_ch2[i] = waveform[b][chn_index][i];}
time[b][chn_index][i] += bin_width[b][chn_index][(j+tch.trigger_cell) % 1024];
} // calculate time for this cell
} for (j=0,time[b][chn_index][i]=0 ; j<i ; j++){
} // end of the channel loop (chn) time[b][chn_index][i] += bin_width[b][chn_index][(j+tch.trigger_cell) % 1024];
}
// align cell #0 of all channels }
t1 = time[b][0][(1024-tch.trigger_cell) % 1024]; } // end of the channel loop (chn)
//my print;
// printf("t1 %1.6lf \n",time[b][0][(1024-tch.trigger_cell) % 1024]); // align cell #0 of all channels
for (chn=1 ; chn<4 ; chn++) { t1 = time[b][0][(1024-tch.trigger_cell) % 1024];
t2 = time[b][chn][(1024-tch.trigger_cell) % 1024]; //my print;
//adding channels 3 and 4 // printf("t1 %1.6lf \n",time[b][0][(1024-tch.trigger_cell) % 1024]);
t3 = time[b][chn][(1024-tch.trigger_cell) % 1024]; for (chn=1 ; chn<4 ; chn++) {
t4 = time[b][chn][(1024-tch.trigger_cell) % 1024]; t2 = time[b][chn][(1024-tch.trigger_cell) % 1024];
//my prinf //adding channels 3 and 4
//printf("t2 %f for %d %d %d \n",time[b][chn][(1024-tch.trigger_cell) % 1024], b, chn, i); t3 = time[b][chn][(1024-tch.trigger_cell) % 1024];
dt = t1 - t2; t4 = time[b][chn][(1024-tch.trigger_cell) % 1024];
dt34 = t3 - t4; //my prinf
for (i=0 ; i<1024 ; i++) { //printf("t2 %f for %d %d %d \n",time[b][chn][(1024-tch.trigger_cell) % 1024], b, chn, i);
time[b][chn][i] += dt; //each element of time gets dt correction dt = t1 - t2;
//my print; dt34 = t3 - t4;
// printf("time %1.6lf for %d %d %d \n",time[b][chn][i], b, chn, i); 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 }
t1 = t2 = t3 = t4 = 0;
//for ROOT threshold = -0.045; //my threshold, used to be 0.3
for(i=0 ; i<1024 ; i++) {
time_ch1[i] = time[b][0][i];
// time_ch2[i] = time[b][1][i]; //for ROOT
} for(i=0 ; i<1024 ; i++) {
time_ch1[i] = time[b][0][i];
// find peak in channel 1 above threshold // time_ch2[i] = time[b][1][i];
for (i=0 ; i<1022 ; i++) { }
if (waveform[b][0][i] < threshold && waveform[b][0][i+1] >= threshold) { // find peak in channel 1 above 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]; for (i=0 ; i<1022 ; i++) {
//my prinf
//printf("t1 recalc %1.6lf %d \n",t1, i); if (waveform[b][0][i] < threshold && waveform[b][0][i+1] >= threshold) {
break; 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]; // find peak in channel 2 above threshold
//my prinf for (i=0 ; i<1022 ; i++) {
//printf("t2 recalc %1.6lf \n",t2); if (waveform[b][1][i] < threshold && waveform[b][1][i+1] >= threshold) {
break; 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; // calculate distance of peaks with statistics
sumdt += dt; if (t1 > 0 && t2 > 0) {
sumdt2 += dt*dt; ndt++;
} dt = t2 - t1;
} //end of the boards loop sumdt += dt;
sumdt2 += dt*dt;
rtree->Fill(); }
} // end of the events loop } //end of the boards loop
// print statistics rtree->Fill();
printf("dT = %1.3lfns +- %1.1lfps\n", sumdt/ndt, 1000*sqrt(1.0/(ndt-1)*(sumdt2-1.0/ndt*sumdt*sumdt))); } // end of the events loop
rfile->Write(); // print statistics
rfile->Close(); printf("dT = %1.3lfns +- %1.1lfps\n", sumdt/ndt, 1000*sqrt(1.0/(ndt-1)*(sumdt2-1.0/ndt*sumdt*sumdt)));
return 1;
rfile->Write();
rfile->Close();
delete event;
return 1;
} }
dataClasses/RawData.cpp
View file @ bf84b293
...@@ -20,7 +20,7 @@ RawData::~RawData() { ...@@ -20,7 +20,7 @@ RawData::~RawData() {
void RawData::Reset() { void RawData::Reset() {
for (Int_t i = 0; i < 1024; i++) { for (Int_t i = 0; i < NCELLS; i++) {
Amp[i] = 0; Amp[i] = 0;
Time[i] = 0; Time[i] = 0;
} }
...@@ -29,7 +29,7 @@ void RawData::Reset() { ...@@ -29,7 +29,7 @@ void RawData::Reset() {
} }
void RawData::Print() { void RawData::Print() {
for (Int_t i = 0; i < 1024; i++) { for (Int_t i = 0; i < NCELLS; i++) {
cout << Amp[i] << endl; cout << Amp[i] << endl;
// Time[i] = 0; // Time[i] = 0;
} }
......
dataClasses/RawData.h
View file @ bf84b293
...@@ -15,11 +15,13 @@ ...@@ -15,11 +15,13 @@
using std::cout; using std::cout;
using std::endl; using std::endl;
#define NCELLS 1024
class RawData { class RawData {
private: private:
Double_t Amp[1024]; Double_t Amp[NCELLS];
Double_t Time[1024]; Double_t Time[NCELLS];
public: public:
RawData(); RawData();
...@@ -36,6 +38,15 @@ public: ...@@ -36,6 +38,15 @@ public:
return Time; return Time;
} }
void SetAmp(Double_t a, Int_t i) {
if (i >=NCELLS) {
cout << "chren'" << endl;
return;
}
Amp[i] = a;
return;
}
void Print(); void Print();
}; };
......
makefile
View file @ bf84b293
...@@ -59,11 +59,11 @@ libData.so: $(DATAOBJS) ...@@ -59,11 +59,11 @@ libData.so: $(DATAOBJS)
@echo 'Finished building target: $@' @echo 'Finished building target: $@'
@echo ' ' @echo ' '
read_binary_DRS4: $(CONVERTDRS4)/read_binary.cpp read_binary_DRS4: libData.so $(CONVERTDRS4)/read_binary.cpp
@echo 'Building target: $@' @echo 'Building target: $@'
@echo 'Invoking: GCC C++ Linker' @echo 'Invoking: GCC C++ Linker'
# $(CC) -L $(ROOTLIBS) -shared -o"libData.so" $(DATAOBJS) $(DATALIBS) # $(CC) -L $(ROOTLIBS) -shared -o"libData.so" $(DATAOBJS) $(DATALIBS)
$(CC) -o $(CONVERTDRS4)/read_binary_DRS4 $(CONVERTDRS4)/read_binary.cpp -lm `root-config --cflags --libs` $(CC) -o $(CONVERTDRS4)/read_binary_DRS4 $(CONVERTDRS4)/read_binary.cpp -lm `root-config --cflags --libs` -L $(PWD) -lData -Wl,-rpath,$(PWD)
@echo 'Finished building target: $@' @echo 'Finished building target: $@'
@echo ' ' @echo ' '
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment