source: trunk/Mars/mcore/factofits.h@ 17262

Last change on this file since 17262 was 17253, checked in by lyard, 12 years ago
zofits with catalog shrinking and compressed calibration table
File size: 13.9 KB
Line 
1/*
2 * factofits.h
3 *
4 * Created on: Oct 16, 2013
5 * Author: lyard
6 */
7
8#ifndef FACTOFITS_H_
9#define FACTOFITS_H_
10
11#include "zofits.h"
12#include "DrsCalib.h"
13
14#ifndef __MARS__
15namespace std
16{
17#else
18using namespace std;
19#endif
20
21class factofits : public zofits
22{
23
24 public:
25
26 /// constructors
27 factofits(uint32_t numTiles=1000,
28 uint32_t rowPerTile=100,
29 uint64_t maxUsableMem=0) : zofits(numTiles, rowPerTile, maxUsableMem)
30 {
31 fStartCellsOffset = -1;
32 fDataOffset = -1;
33 }
34
35 factofits(const char* fname,
36 uint32_t numTiles=1000,
37 uint32_t rowPerTile=100,
38 uint64_t maxUsableMem=0) : zofits(fname, numTiles, rowPerTile, maxUsableMem)
39 {
40 fStartCellsOffset = -1;
41 fDataOffset = -1;
42 }
43
44 virtual ~factofits()
45 {
46 }
47
48 /// whether or not a calibration was given to the file writer
49 virtual bool IsOffsetCalibration()
50 {
51 return (fOffsetCalibration.size() != 0);
52 }
53
54 ///assign a given drs offset calibration
55 void SetDrsCalibration(const vector<float> &calib)
56 {
57 VerifyCalibrationSize(calib.size());
58
59 if (!IsOffsetCalibration())
60 fOffsetCalibration.resize(1440*1024);
61
62 for (uint32_t i=0; i<1440*1024; i++)
63 fOffsetCalibration[i] = (int16_t)(calib[i]*4096/2000);
64 }
65
66 ///assign a given drs offset calibration
67 void SetDrsCalibration(const vector<int16_t>& vec)
68 {
69 VerifyCalibrationSize(vec.size());
70
71 if (!IsOffsetCalibration())
72 fOffsetCalibration.resize(1440*1024);
73
74 for (uint32_t i=0; i<1440*1024; i++)
75 fOffsetCalibration[i] = vec[i];
76 }
77
78 ///assign a given drs offset calibration
79 void SetDrsCalibration(const DrsCalibration& drs)
80 {
81 if (drs.fNumOffset==0)
82 return;
83
84 VerifyCalibrationSize(drs.fOffset.size());
85
86 if (!IsOffsetCalibration())
87 fOffsetCalibration.resize(1440*1024);
88
89 for (uint32_t i=0; i<1024*1440; i++)
90 fOffsetCalibration[i] = drs.fOffset[i]/drs.fNumOffset;
91 }
92
93 ///Overload of the super function
94 bool WriteTableHeader(const char* name="DATA")
95 {
96 if (!zofits::WriteTableHeader(name))
97 return false;
98
99 if (IsOffsetCalibration())
100 {//retrieve the column storing the start cell offsets, if required.
101
102 for (auto it=fRealColumns.begin(); it!=fRealColumns.end(); it++)//Table.cols.begin(); it!= fTable.cols.end(); it++)
103 {
104 if (it->col.name == "StartCellData")
105 fStartCellsOffset = it->col.offset;
106
107 if (it->col.name == "Data")
108 {
109 fNumSlices = it->col.num;
110 fDataOffset = it->col.offset;
111 if (fNumSlices % 1440 != 0)
112 {
113#ifdef __EXCEPTIONS
114 throw runtime_error("Number of data samples not a multiple of 1440.");
115#else
116 gLog << ___err___ << "ERROR - Number of data samples not a multiple of 1440. Doing it uncalibrated." << endl;
117#endif
118 fOffsetCalibration.resize(0);
119 }
120 fNumSlices /= 1440;
121 }
122 }
123
124 if (fStartCellsOffset < 0)
125 {
126#ifdef __EXCEPTIONS
127 throw runtime_error("FACT Calibration requested, but \"StartCellData\" column not found.");
128#else
129 gLog << ___err___ << "ERROR - FACT Calibration requested, but \"StartCellData\" column not found. Doing it uncalibrated." << endl;
130#endif
131 //throw away the calibration data
132 fOffsetCalibration.resize(0);
133 }
134
135 if (fDataOffset < 0)
136 {
137#ifdef __EXCEPTIONS
138 throw runtime_error("FACT Calibration requested, but \"Data\" column not found.");
139#else
140 gLog << ___err___ << "ERROR - FACT Calibration requested, but \"Data\" column not found. Doing it uncalibrated." << endl;
141#endif
142 //throw away the calibration data
143 fOffsetCalibration.resize(0);
144 }
145 }
146
147 return true;
148 }
149
150 ///Uncompressed version of the DrsCalibration table
151 /* virtual bool WriteDrsOffsetsTable()
152 {
153 if (!IsOffsetCalibration())
154 return false;
155
156 ofits c;
157 c.SetStr("XTENSION", "BINTABLE" , "binary table extension");
158 c.SetInt("BITPIX" , 8 , "8-bit bytes");
159 c.SetInt("NAXIS" , 2 , "2-dimensional binary table");
160 c.SetInt("NAXIS1" , 1024*1440*2 , "width of table in bytes");
161 c.SetInt("NAXIS2" , 1 , "number of rows in table");
162 c.SetInt("PCOUNT" , 0 , "size of special data area");
163 c.SetInt("GCOUNT" , 1 , "one data group (required keyword)");
164 c.SetInt("TFIELDS" , 1 , "number of fields in each row");
165 c.SetStr("CHECKSUM", "0000000000000000" , "Checksum for the whole HDU");
166 c.SetStr("DATASUM" , " 0" , "Checksum for the data block");
167 c.SetStr("EXTNAME" , "ZDrsCellOffsets" , "name of this binary table extension");
168 c.SetStr("TTYPE1" , "OffsetCalibration" , "label for field 1");
169 c.SetStr("TFORM1" , "1474560I" , "data format of field: 2-byte INTEGER");
170 c.End();
171
172 vector<char> swappedOffsets;
173 swappedOffsets.resize(1024*1440*sizeof(int16_t));
174 revcpy<sizeof(int16_t)>(swappedOffsets.data(), (char*)(fOffsetCalibration.data()), 1024*1440);
175
176 Checksum datasum;
177 datasum.add(swappedOffsets.data(), sizeof(int16_t)*1024*1440);
178
179 ostringstream dataSumStr;
180 dataSumStr << datasum.val();
181 c.SetStr("DATASUM", dataSumStr.str());
182
183 datasum += c.WriteHeader(*this);
184
185 const off_t here_I_am = tellp();
186
187 c.SetStr("CHECKSUM", datasum.str());
188 c.WriteHeader(*this);
189
190 seekp(here_I_am);
191
192 write(swappedOffsets.data(), swappedOffsets.size());
193
194 AlignTo2880Bytes();
195
196 return good();
197 }*/
198 ///Actually write the drs calibration table
199 virtual bool WriteDrsOffsetsTable()
200 {
201 if (!IsOffsetCalibration())
202 return false;
203
204 const uint32_t catalog_size = sizeof(int64_t)*2;
205
206 ofits c;
207 c.SetStr("XTENSION", "BINTABLE" , "binary table extension");
208 c.SetInt("BITPIX" , 8 , "8-bit bytes");
209 c.SetInt("NAXIS" , 2 , "2-dimensional binary table");
210 c.SetInt("NAXIS1" , catalog_size , "width of table in bytes");
211 c.SetInt("NAXIS2" , 1 , "number of rows in table");
212 c.SetInt("PCOUNT" , 0 , "size of special data area");
213 c.SetInt("GCOUNT" , 1 , "one data group (required keyword)");
214 c.SetInt("TFIELDS" , 1 , "number of fields in each row");
215 c.SetStr("CHECKSUM", "0000000000000000" , "Checksum for the whole HDU");
216 c.SetStr("DATASUM" , " 0" , "Checksum for the data block");
217 c.SetStr("EXTNAME" , "ZDrsCellOffsets" , "name of this binary table extension");
218 c.SetStr("TTYPE1" , "OffsetCalibration" , "label for field 1");
219 c.SetStr("ZFORM1" , "1474560I" , "data format of field: 2-byte INTEGER");
220 c.SetStr("TFORM1" , "1QB" , "data format of variable length bytes");
221 c.SetStr("ZCTYP1" , "FACT" , "Compression type FACT");
222
223 c.SetBool( "ZTABLE", true, "Table is compressed");
224 c.SetInt( "ZNAXIS1", 1024*1440*2, "Width of uncompressed rows");
225 c.SetInt( "ZNAXIS2", 1, "Number of uncompressed rows");
226 c.SetInt( "ZPCOUNT", 0, "");
227 c.SetInt( "ZHEAPPTR", catalog_size, "");
228 c.SetInt( "ZTILELEN", 1, "Number of rows per tile");
229 c.SetInt( "THEAP", catalog_size, "");
230 c.SetStr( "RAWSUM", " 0", "Checksum of raw little endian data");
231 c.SetFloat("ZRATIO", 0, "Compression ratio");
232 c.SetInt( "ZSHRINK", 1, "Catalog shrink factor");
233 c.End();
234
235 c.WriteHeader(*this);
236
237 const off_t here_I_am = tellp();
238
239 //go after the catalog to compress and write the table data
240 seekp(here_I_am + catalog_size);
241
242 //calculate RAWSUM
243 Checksum rawsum;
244 rawsum.add((char*)(fOffsetCalibration.data()), 1024*1440*sizeof(int16_t));
245 ostringstream sum_str;
246 sum_str << rawsum.val();
247 c.SetStr("RAWSUM", sum_str.str());
248
249 //compress data and calculate final, compressed size
250 const uint32_t compressed_header_size = sizeof(TileHeader) + sizeof(BlockHeader) + 1*sizeof(uint16_t);
251 vector<char> compressed_calib(1024*1440*2 + compressed_header_size + 8); //+8 for checksum;
252 char* data_start = compressed_calib.data() + compressed_header_size;
253 uint32_t compressed_size = compressHUFFMAN16(data_start, (char*)(fOffsetCalibration.data()), 1024*1440, 2, 1);;
254 compressed_size += compressed_header_size;
255
256 //Write tile header
257 TileHeader th;
258 vector<uint16_t> seq(1, kFactHuffman16);
259 Compression bh(seq, kOrderByRow);
260 th.numRows = 1;
261 th.size = compressed_size;
262 bh.SetBlockSize(compressed_size-sizeof(TileHeader));
263 memcpy(compressed_calib.data(), &(th), sizeof(TileHeader));
264 bh.Memcpy(compressed_calib.data()+sizeof(TileHeader));
265
266 //calculate resulting compressed datasum
267 Checksum datasum;
268 memset(compressed_calib.data()+compressed_size, 0, 8-compressed_size%8);
269 datasum.add(compressed_calib.data(), compressed_size + 8-compressed_size%8);
270
271 //write the catalog !
272 seekp(here_I_am);
273
274 vector<uint64_t> catalog(2,0);
275 catalog[0] = compressed_size-sizeof(TileHeader);
276 catalog[1] = sizeof(TileHeader);
277
278 vector<char> swappedCatalog(catalog_size);
279 revcpy<sizeof(int64_t)>(swappedCatalog.data(), (char*)(catalog.data()), 2);//catalog_size);
280 datasum.add(swappedCatalog.data(), catalog_size);
281
282 write(swappedCatalog.data(), catalog_size);
283
284 //update relevant keywords
285 c.SetFloat("ZRATIO", (float)(1024*1440*2)/(float)(compressed_size));
286 c.SetInt("PCOUNT", compressed_size + catalog_size);
287
288 sum_str.str("");
289 sum_str << datasum.val();
290 c.SetStr("DATASUM", sum_str.str());
291
292 datasum += c.WriteHeader(*this);
293
294 c.SetStr("CHECKSUM", datasum.str());
295
296 c.WriteHeader(*this);
297
298 //write the compressed data
299 seekp(here_I_am + catalog_size);
300 write(compressed_calib.data(), compressed_size);
301
302 AlignTo2880Bytes();
303
304 return good();
305 }
306 ///Apply the drs offset calibration (overload of super-method)
307 virtual void DrsOffsetCalibrate(char* target_location)
308 {
309 if (IsOffsetCalibration())
310 {
311 const int16_t* startCell = reinterpret_cast<int16_t*>(target_location + fStartCellsOffset);
312 int16_t* data = reinterpret_cast<int16_t*>(target_location + fDataOffset);
313
314 for (uint32_t ch=0; ch<1440; ch++)
315 {
316 if (startCell[ch] < 0)
317 {
318 data += fNumSlices;
319 continue;
320 }
321
322 const int16_t modStart = startCell[ch]%1024;
323 const int16_t *off = fOffsetCalibration.data() + ch*1024;
324
325 const int16_t* cal = off+modStart;
326 const int16_t* end_stride = data+fNumSlices;
327
328 if (modStart+fNumSlices > 1024)
329 {
330 while (cal < off+1024)
331 *data++ -= *cal++;
332 cal = off;
333 }
334
335 while (data<end_stride)
336 *data++ -= *cal++;
337 }
338 }
339 }
340
341private:
342 /// Checks if the size of the input calibration is ok
343 void VerifyCalibrationSize(uint32_t size)
344 {
345 if (size != 1440*1024)
346 {
347 ostringstream str;
348 str << "Cannot load calibration with anything else than 1440 pixels and 1024 samples per pixel. Got a total size of " << size;
349#ifdef __EXCEPTIONS
350 throw runtime_error(str.str());
351#else
352 gLog << ___err___ << "ERROR - " << str.str() << endl;
353#endif
354 }
355 }
356
357 //Offsets calibration stuff.
358 vector<int16_t> fOffsetCalibration; ///< The calibration itself
359 int32_t fStartCellsOffset; ///< Offset in bytes for the startcell data
360 int32_t fDataOffset; ///< Offset in bytes for the data
361 int32_t fNumSlices; ///< Number of samples per pixel per event
362
363}; //class factofits
364
365#ifndef __MARS
366}; //namespace std
367#endif
368
369#endif /* FACTOFITS_H_ */
Note: See TracBrowser for help on using the repository browser.