-
Notifications
You must be signed in to change notification settings - Fork 0
/
b2r.C
183 lines (162 loc) · 6.67 KB
/
b2r.C
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
// Save binary output from Struck ADCs to ROOT trees.
// The binary file info. is fetched from its index file generated by idx.C
void b2r(const char* index_file = "index.csv")
{
cout<<"Open "<<index_file<<endl;
ifstream index(index_file);
if (index.is_open()==false) { cout<<"Failed! Quit"<<endl; return; }
string line, hash, tmp, binary_file, unit, bit;
index>>hash>>tmp>>tmp>>tmp>>binary_file;
cout<<"Open "<<binary_file<<endl;
ifstream input(binary_file, ios::binary);
if (index.is_open()==false) { cout<<"Failed! Quit"<<endl; return; }
cout<<"Get sampling rate: ";
int nm=21, nc=16; float sampling_rate=250, dt;
index>>hash>>tmp>>tmp>>tmp>>tmp>>tmp>>hash>>tmp>>tmp>>sampling_rate>>unit;
cout<<sampling_rate<<" "<<unit<<endl;
if (unit=="MHz") dt=1000/sampling_rate; // dt unit: ns
else dt=1/sampling_rate; // unit: GHz
cout<<"Get number of modules used: ";
index>>hash>>tmp>>tmp>>tmp>>tmp>>nm;
cout<<nm<<endl;
cout<<"Get number of channels per module: ";
index>>hash>>tmp>>tmp>>tmp>>tmp>>tmp>>nc;
cout<<nc<<endl;
cout<<"Get channel status (0: empty, 1: used):"<<endl;
for (int i=0; i<3+nm; i++) getline(index,line); // skip lines about ch Id
bool used[21][16] = {0}; // whether a channel is used
for (int m=0; m<nm; m++) {
index>>hash;
for (int c=0; c<nc; c++) { index>>used[m][c]; cout<<used[m][c]<<" "; }
cout<<endl;
}
cout<<"Get channel sync requirement (1: sync, 0: no need to sync):"<<endl;
for (int i=0; i<2; i++) getline(index,line); // skip 2 lines
bool sync[21][16] = {0}; // whether events in a channel are synced
for (int m=0; m<nm; m++) {
index>>hash;
for (int c=0; c<nc; c++) { index>>sync[m][c]; cout<<sync[m][c]<<" "; }
cout<<endl;
}
cout<<"Get channel format bits:"<<endl;
for (int i=0; i<2; i++) getline(index,line); // skip 2 lines
bitset<4> format[21][16];
for (int m=0; m<nm; m++) {
index>>hash;
for (int c=0; c<nc; c++) {
index>>bit; istringstream bit_stream(bit);
bit_stream>>format[m][c];
cout<<format[m][c]<<" ";
}
cout<<endl;
}
cout<<"Get spill positions and sizes: ";
for (int i=0; i<6; i++) getline(index,line); // skip 6 lines
int nspill=0; // total number of spills in file
int pos[20000] = {0}; // position of each spill in file (max 20000 spills)
int min_size[20000] = {0}; // minimal size of a spill in all channels in a card
while (getline(index, line)) {
nspill=stoi(line.substr(0,7));
pos[nspill]=stoi(line.substr(9,19));
min_size[nspill]=stoi(line.substr(21,28));
}
nspill++; cout<<nspill<<" spills to be read"<<endl<<endl;
// generate root file name
TString root_file(index_file); root_file.ReplaceAll("csv","root");
cout<<"\nCreate "<<root_file<<endl;
TFile *output = new TFile(root_file, "recreate");
int sum[8]={0}; // accumulator sums
bool pu; int ih, max, vbt, vat; float es, em, h;
int n, ns; unsigned long long ts; float s[99999]={0}, t[99999];
for (int i=0; i<99999; i++) t[i]=i*dt;
cout<<"Create a tree per channel"<<endl<<endl;
TTree* tree[21][16]={0};
for (int m=0; m<nm; m++) {
for (int c=0; c<nc; c++) {
if (used[m][c]==false) continue; // don't create tree for empty channel
tree[m][c] = new TTree(Form("t%d",m*nc+c),
Form("SIS3316 channel %d",m*nc+c));
tree[m][c]->Branch("n", &n, "n/I"); // number of waveform samples
tree[m][c]->Branch("s", s, "s[n]/F"); // waveform samples
tree[m][c]->Branch("t", t, "t[n]/F"); // time [ns] of waveform samples
tree[m][c]->Branch("ts", &ts, "ts/l"); // 48-bit event timestamp
tree[m][c]->Branch("pu", &pu, "pu/O"); // pile-up flag
if (format[m][c][0]==1) {
tree[m][c]->Branch("ns", &ns, "ns/I"); // number of accumulator sums
tree[m][c]->Branch("sum", sum, "sum[ns]/I"); // accumulator sums
tree[m][c]->Branch("h", &h, "h/F"); // height of a waveform
tree[m][c]->Branch("ih", &ih, "ih/I"); // sample index of the highest point
}
if (format[m][c][2]==1) {
tree[m][c]->Branch("max", &max, "max/I"); // MAW maximum value
tree[m][c]->Branch("vbt", &vbt, "vbt/I"); // MAW value before trigger
tree[m][c]->Branch("vat", &vat, "vat/I"); // MAW value after trigger
}
if (format[m][c][3]==1) {
tree[m][c]->Branch("es", &es, "es/F"); // start energy value
tree[m][c]->Branch("em", &em, "em/F"); // max energy value
}
}
}
cout<<"Convert binary input to ROOT trees..."<<endl;
unsigned int word[100]; // read not more than 100 words at a time
char* byte = (char*) word; // index in unit of byte instead of word
for (int ispill=0; ispill<nspill; ispill++) { // loop over spills
input.seekg(pos[ispill], ios::beg); // jump to the start of a spill
for (int m=0; m<nm; m++) { // loop over modules
input.seekg(8, ios::cur); // skip module header (8 bytes)
for (int c=0; c<nc; c++) { // loop over channels
input.read(byte,32); // get channel header (8 words)
int nwords=word[7]; // size of a spill of data for a channel (words)
int dw = nwords-min_size[ispill];
if (sync[m][c]) nwords=min_size[ispill];
while (nwords>0) { // loop over events
input.read(byte,8); nwords-=2; // get event header
ts = *word&0xffff0000; // upper timestamp (started at 16th bit)
ts = ts<<16; // move it up by 16 bit again
// combine upper (16) and lower (32) bits
ts = ts | static_cast<unsigned long long>(word[1]);
if (format[m][c][0]==1) { // peak + accumulator sums (gate 1~6)
input.read(byte,28); nwords-=7;
h = float(*word&0xffff); ih = int((*word&0xffff0000)>>16);
ns=6; sum[0]=word[1]&0xffffff;
for (int i=1; i<6; i++) sum[i]=word[i+1];
}
if (format[m][c][1]==1) { // accumulator sums from gate 7 & 8
input.read(byte,8); nwords-=2;
ns=8; sum[6]=word[0]; sum[7]=word[1];
}
if (format[m][c][2]==1) { // moving average window (MAW) info.
input.read(byte,12); nwords-=3;
max=word[0]; vbt=word[1]; vat=word[2];
}
if (format[m][c][3]==1) { // energy values
input.read(byte,8); nwords-=2;
es=word[0]; em=word[1];
}
input.read(byte,4); nwords-=1;
n=2*((*word)&0x3ffffff); // number of waveform samples
if ((byte[3]&0x4)>0) pu=1; else pu=0; // pile-up flag
for (int i=0; i<n; i+=2) {
input.read(byte,4); nwords-=1; // read two samples
s[i]=*word&0xffff; s[i+1]=(*word&0xffff0000)>>16;
}
if (used[m][c]) tree[m][c]->Fill();
}
if (sync[m][c] && dw>0) input.seekg(dw*4,ios::cur);
}
}
cout<<ispill+1<<"/"<<nspill<<" spills processed"<<endl;
}
input.close();
cout<<"\nWrite trees to "<<output->GetName()<<endl;
for (int m=0; m<nm; m++) {
for (int c=0; c<nc; c++) {
if (used[m][c]==false) continue;
tree[m][c]->Write("", TObject::kOverwrite);
printf("# of events in channel %2d: %8llu\n",
m*nc+c, tree[m][c]->GetEntries());
}
}
output->Close();
}