source: fact/tools/pyscripts/sandbox/dneise/fact_compress/c++/comp2.cpp@ 14269

Last change on this file since 14269 was 14268, checked in by neise, 12 years ago
evolving
File size: 12.3 KB
Line 
1#include <iostream>
2#include <fstream>
3#include <vector>
4#include <string.h>
5#include <math.h>
6
7// I will use the fits class by TB for reading the interesting data from the files
8
9#include "izstream.h"
10
11#define HAVE_ZLIB
12#include "fits.h"
13
14#include "helper.hxx"
15
16using namespace std;
17
18ifstream::pos_type size;
19char * memblock;
20
21void make_groups( vector<int>& group_starts , vector<int>& group_length, short * sizes, int len_sizes, int roi);
22long get_event_size( vector<int>& gs, vector<int>& gl, short * sizes);
23
24int main (int argc, char * argv[]) {
25
26 char * data_file_name = "20120712_004.fits";
27 char * calib_file_name ="20120712_003.drs.fits.gz";
28 char * out_file_name = NULL;
29 out_file_name = new char[strlen(data_file_name)+4];
30 strcpy(out_file_name, data_file_name);
31 strcpy(out_file_name+strlen(data_file_name),".fc");
32 out_file_name[strlen(data_file_name)+3] = 0;
33
34
35 //=========================================================================
36 // CALIBRATION CONSTANTS
37 //=========================================================================
38 fits * calib =new fits(calib_file_name);
39 calib->PrintKeys();
40 calib->PrintColumns();
41
42 int size_of_offset = 1440*1024;
43 float * offset_mV = new float[size_of_offset];
44 short * offset = new short[size_of_offset];
45 calib->SetPtrAddress("BaselineMean",offset_mV, size_of_offset);
46
47 calib->GetNextRow();
48 for (int i =0 ; i<size_of_offset; i++)
49 {
50 offset[i] = short(offset_mV[i] / 2000. * 4096 - 0.5);
51 }
52 cout << endl;
53 delete[] offset_mV;
54 delete calib;
55 calib = NULL;
56 //=========================================================================
57 // END OF CALIBRATION CONSTANTS
58 //=========================================================================
59
60
61 //=========================================================================
62 // DATA FILE - with FITS class
63 //=========================================================================
64
65 fits * datafits = new fits(data_file_name);
66
67 datafits->PrintKeys();
68 datafits->PrintColumns();
69
70 int nevents = (int)datafits->GetNumRows();
71 int roi = (int)datafits->GetInt("NROI");
72 int npixel = (int)datafits->GetInt("NPIX");
73 int event_size = npixel * roi;
74
75
76
77 short *event = new short[event_size];
78 short *sc = new short[npixel];
79 datafits->SetPtrAddress("Data",event, event_size );
80 datafits->SetPtrAddress("StartCellData",sc, npixel );
81
82 short first_word;
83 short *diffs = new short[event_size-1];
84 short *sizes = new short[event_size-1];
85
86 //=========================================================================
87 // DATA FILE - with ISTREAM
88 //=========================================================================
89 ifstream datafile (data_file_name, ios::in|ios::binary);
90 ofstream out (out_file_name, ios::out|ios::binary|ios::trunc);
91 helper h(&datafile, &out);
92
93
94 if (datafile.is_open() && out.is_open())
95 {
96 // create our own header
97// out.write( "FACT_COMPRESS ", 80);
98// out.write( "END. ", 80);
99
100// h.copy_ascii_header();
101
102 for ( int event_id = 0 ; event_id < nevents; event_id++)
103 {
104// h.copy_event_header();
105 printf("last event header @: 0x%10X \n", h.event_header_start);
106
107
108 // jump over npixel * roi * shorts
109 datafile.seekg( datafile.tellg() + long(npixel * roi * sizeof(short)) , ios::beg);
110 // since the data is stored big-endian in the FITS file, I can't memcopy it ..
111 // have to use the methods of the fits-class
112
113 datafits->GetNextRow();
114
115 // perform a little DRS calibration
116 for ( int pix=0; pix<npixel; pix++){
117 short *pevt = event+pix*roi;
118 short cell = sc[pix];
119 short *poff = offset+pix*1024;
120
121 for ( int sl=0 ; sl<roi; sl++)
122 {
123 pevt[sl] -= poff[(sl+cell)%1024];
124 }
125
126 }
127
128
129 first_word = event[0];
130 //event_size = 100;
131
132 short d;
133 for (int i=0; i<event_size-1; i++)
134 {
135 diffs[i] = event[i+1]-event[i];
136 d = diffs[i];
137 for (int b =4; b<20; b++)
138 {
139 if ( (d >= -1*(1<<b)) && (d <= (1<<b)-1))
140 {
141 sizes[i] = b+1;
142 break;
143 }
144 }
145 }
146
147 vector<int> gs;
148 vector<int> gl;
149 make_groups( gs , gl, sizes, event_size-1, event_size);
150
151 long cevts = get_event_size( gs,gl, sizes);
152
153 cout << cevts << endl;
154 cout << "number of groups:" << gs.size() << endl;
155
156 // write the length of the event
157 out.write( (char*) &cevts, sizeof(cevts) );
158 // write the event to disk
159 out.write( (char*)&first_word, sizeof(first_word));
160 unsigned short header;
161 for (int g=0; g<gs.size(); g++){
162
163 if (g<50)
164 // printf("group: %3d length: %4d -- size: %2d ", g, gl[g], sizes[gs[g]]);
165
166 header = (unsigned short)gl[g] & 0x0FFF;
167 header |= (sizes[gs[g]]-5) << 12;
168
169 out.write( (char*)&header, sizeof(header));
170 int group_size = ceil((gl[g]*sizes[gs[g]])/8.);
171 if (g<50)
172 // printf("group_size: %4d ", group_size);
173
174 int size = sizes[gs[g]];
175 unsigned char * buffer = new unsigned char[ group_size ];
176 for (int i=0; i<group_size; i++) buffer[i] = 0;
177 unsigned short *data = (unsigned short *)diffs;
178
179 int *bs = new int[gl[g]];
180 int *Bs = new int[gl[g]];
181
182 for (int i=0; i<gl[g] ;i++)
183 {
184 int bit_pos = ((i+1)*size)-1;
185 Bs[i] = bit_pos/8;
186 bs[i] = bit_pos%8;
187 }
188
189 for (int i=0; i<gl[g] ;i++)
190 {
191 int ii = i+gs[i];
192 unsigned char dat = (unsigned char)(diffs[ii]);
193 if (bs[i]+1 >= size) // diff entry fits in one byte
194 {
195 buffer[Bs[i]] |= dat <<(bs[i]+1-size);
196 }
197 else // entry is distributed over two bytes
198 {
199 buffer[Bs[i]] |= dat>>(size - bs[i]+1);
200 buffer[Bs[i-1]] |= dat<<(8-size - bs[i]+1);
201 }
202 }
203 delete[] bs;
204 delete[] Bs;
205 if (g<50)
206 {
207 int tellp = out.tellp();
208// printf("start_addres: 0x%8X\n", tellp);
209 }
210
211
212 out.write( (char*)buffer, group_size);
213 delete[] buffer;
214
215 }
216
217
218
219 return 0;
220/*
221
222
223 // measure how large the stuff, is
224 int cs =0;
225 // for all groups write header and group
226 short header = 0;
227 int diff_index = 0;
228 for ( int i = 0 ; i < (int)groups.size() ; i++)
229 {
230 // write header
231 if (groups[i] == 8)
232 {
233 header = short(group_sizes[i]);
234 }
235 else
236 {
237 header = 0x8000 | short(group_sizes[i] );
238 }
239 cs += sizeof(short);
240
241 // write group
242 if (groups[i] == 8)
243 {
244 for (int j = 0; j<group_sizes[i]; j++)
245 {
246 cs += 1;
247 }
248 }
249 else
250 {
251 for (int j = 0; j<group_sizes[i]; j++)
252 {
253 cs += 2;
254 }
255 }
256 }
257 //cout << "compressed size:" << cs << endl;
258 out.write( (char*)&cs, sizeof(cs));
259
260 cout << "sizeof(cs):" << sizeof(cs) << endl;
261
262 // write the first short
263 out.write(memblock, 1*sizeof(short) );
264 printf("first short: 0x%X \t %d \n", ( (short*)memblock )[0], ( (short*)memblock )[0] );
265 // for all groups write header and group
266 header = 0;
267 diff_index = 0;
268 for ( int i = 0 ; i < (int)groups.size() ; i++)
269 {
270 // write header
271 if (groups[i] == 8)
272 {
273 header = short(group_sizes[i]);
274 }
275 else
276 {
277 header = 0x8000 | short(group_sizes[i] );
278 }
279 out.write( (char*)&header, sizeof(short));
280
281 // write group
282 if (groups[i] == 8)
283 {
284 for (int j = 0; j<group_sizes[i]; j++)
285 {
286 out.write( (char*)&(diffs[diff_index++]), 1);
287 }
288 }
289 else
290 {
291 for (int j = 0; j<group_sizes[i]; j++)
292 {
293 out.write( (char*)&(diffs[diff_index++]), 2);
294 }
295 }
296 }
297
298 //out.write(memblock, data_size*sizeof(short) );
299
300 delete[] memblock;
301 delete[] diffs;
302 delete[] sizes;
303*/
304 //}
305
306 }
307 cout << "finished with 1000 events." << endl;
308 long after_address = datafile.tellg();
309 datafile.seekg (0, ios::end);
310 long end = datafile.tellg();
311 datafile.seekg (after_address, ios::beg);
312
313 cout << "between last event and end:" << end - after_address << endl;
314 // read the last piece...
315 const int rest_size = end-after_address;
316 char * memblock2 = new char [rest_size];
317 datafile.read(memblock2, rest_size);
318 cout << "first char in memblock: " << int(memblock2[0]) << endl;
319 char lastchar = memblock2[0];
320 for (int i =0 ; i<rest_size; i++)
321 {
322 if (memblock2[i] != lastchar)
323 {
324 cout << "new char at: " << i << " with value:" << int(memblock[0]) << endl;
325 lastchar = memblock2[i];
326 }
327 }
328
329 out.write(memblock2, rest_size);
330 delete[] memblock2;
331
332
333 datafile.close();
334 out.close();
335 } //end of if file open
336} //end of main
337
338void make_groups( vector<int>& group_starts , vector<int>& group_length, short * sizes, int len_sizes, int roi)
339{
340
341 short * diff_sizes = new short[len_sizes-1];
342
343 for (int i=0; i<len_sizes-1; i++)
344 {
345 diff_sizes[i] = sizes[i+1] - sizes[i];
346 }
347
348 vector<int> gs;
349 vector<int> gl;
350
351
352 gs.push_back(0);
353 for (int i=0; i<len_sizes-1; i++)
354 {
355 if (diff_sizes[i] != 0)
356 {
357 gs.push_back(i+1);
358 }
359 }
360 gs.push_back(roi-1);
361
362 // make group_length the same length as group_start
363 gl.reserve( gs.size()-1 );
364 for (int i=0; i<gs.size()-1; i++)
365 {
366 gl[i] = gs[i+1] - gs[i];
367 }
368
369 gs.erase(gs.begin()+gs.size()-1);
370
371 group_starts.swap(gs);
372 group_length.swap(gl);
373}
374
375long get_event_size( vector<int>& gs, vector<int>& gl, short * sizes)
376{
377 long bits;
378 // eventsize in bytes
379 long evts = 2; // the first word
380
381 // gs and gl have equal length
382 for (int i=0; i<gs.size(); i++)
383 {
384 evts += 2; // group header
385 bits = sizes[gs[i]] * gl[i];
386 evts += ceil( bits/8.);
387
388 }
389 return evts;
390}
Note: See TracBrowser for help on using the repository browser.