source: fact/tools/pyscripts/sandbox/dneise/fact_compress/c++/comp1.cpp@ 16858

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