diff --git a/convertDRS4/read_binary.cpp b/convertDRS4/read_binary.cpp index 1b66e0907e593fa67b53ad5abeaabad2ee5ccd48..ac0f190a22d1d1ce7b54e93c6230abcce21026a4 100644 --- a/convertDRS4/read_binary.cpp +++ b/convertDRS4/read_binary.cpp @@ -4,11 +4,11 @@ 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 This program assumes that a pulse from a signal generator is split @@ -17,7 +17,7 @@ for time measurements. $Id: read_binary.cpp 22290 2016-04-27 14:51:37Z ritt $ -*/ + */ #include #include @@ -31,299 +31,311 @@ #include "TTree.h" #include "TH1F.h" +//our code +#include "../dataClasses/RawData.h" + typedef struct { - char tag[3]; - char version; + char tag[3]; + char version; } FHEADER; typedef struct { - char time_header[4]; + char time_header[4]; } THEADER; typedef struct { - char bn[2]; - unsigned short board_serial_number; + 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 + 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; + char tc[2]; + unsigned short trigger_cell; } TCHEADER; typedef struct { - char c[1]; - char cn[3]; + 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]; //for input binary file - char outroot[256]; //for output root file - - int ndt; - double threshold, sumdt, sumdt2; - - if (argc == 3) { - strcpy(filename, argv[1]); - strcpy(outroot, argv[2]); - } - else if (argc == 2) { - printf("Error: both input binary file and output root file should be specified!\n"); - return -1; - } - else { - printf("Error: input binary file and output root file should be specified!\n"); - return -1; - } - -// ---------------for ROOT - - TFile* rfile = new TFile(outroot, "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"); - -//---------------- - - // 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 + 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]; //for input binary file + char outroot[256]; //for output root file + + int ndt; + double threshold, sumdt, sumdt2; + + if (argc == 3) { + strcpy(filename, argv[1]); + strcpy(outroot, argv[2]); + } + else if (argc == 2) { + printf("Error: both input binary file and output root file should be specified!\n"); + return -1; + } + else { + printf("Error: input binary file and output root file should be specified!\n"); + return -1; + } + + // ---------------for ROOT + + TFile* rfile = new TFile(outroot, "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"); + + RawData *event = new RawData(); + rtree->Bronch("rawEvent", "RawData", &event); + + //------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); + 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 + // 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 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= 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; + } + } + } + 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 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]; + event->SetAmp(waveform[b][chn_index][i], 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= 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(); + + delete event; + + return 1; } diff --git a/data/dataDSR4/analysis_07_1.root b/data/dataDSR4/analysis_07_1.root index a02252f0f4c0182996abcff0496f4a0e055c9b6a..ab55ba18c0c18f5334c1632cd6b0892cb1f6d941 100644 Binary files a/data/dataDSR4/analysis_07_1.root and b/data/dataDSR4/analysis_07_1.root differ diff --git a/dataClasses/RawData.cpp b/dataClasses/RawData.cpp index 0d9a50fc3e0f6702c60b0389beaee9876f46e15b..29994d88fd645956e1a3e955e2b51c4b58b27f3f 100644 --- a/dataClasses/RawData.cpp +++ b/dataClasses/RawData.cpp @@ -20,7 +20,7 @@ RawData::~RawData() { void RawData::Reset() { - for (Int_t i = 0; i < 1024; i++) { + for (Int_t i = 0; i < NCELLS; i++) { Amp[i] = 0; Time[i] = 0; } @@ -29,7 +29,7 @@ void RawData::Reset() { } void RawData::Print() { - for (Int_t i = 0; i < 1024; i++) { + for (Int_t i = 0; i < NCELLS; i++) { cout << Amp[i] << endl; // Time[i] = 0; } diff --git a/dataClasses/RawData.h b/dataClasses/RawData.h index 52200a277abc5fa248199f82ec75cbe096e0a64b..da26666e006a45ddf37177704b68402793a52240 100644 --- a/dataClasses/RawData.h +++ b/dataClasses/RawData.h @@ -15,11 +15,13 @@ using std::cout; using std::endl; +#define NCELLS 1024 + class RawData { private: - Double_t Amp[1024]; - Double_t Time[1024]; + Double_t Amp[NCELLS]; + Double_t Time[NCELLS]; public: RawData(); @@ -36,6 +38,15 @@ public: return Time; } + void SetAmp(Double_t a, Int_t i) { + if (i >=NCELLS) { + cout << "chren'" << endl; + return; + } + Amp[i] = a; + return; + } + void Print(); }; diff --git a/makefile b/makefile index 1c9e420904e5e0b9c8f526c01098a9c163592582..ff16c52f84eece69110f928a59a4afd6889bbc9f 100755 --- a/makefile +++ b/makefile @@ -59,11 +59,11 @@ libData.so: $(DATAOBJS) @echo 'Finished building target: $@' @echo ' ' -read_binary_DRS4: $(CONVERTDRS4)/read_binary.cpp +read_binary_DRS4: libData.so $(CONVERTDRS4)/read_binary.cpp @echo 'Building target: $@' @echo 'Invoking: GCC C++ Linker' # $(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 ' '