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

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