source: branches/Mars_McMismatchStudy/mcore/factofits.h@ 20051

Last change on this file since 20051 was 17304, checked in by lyard, 11 years ago
Adapted compilation on older Gcc
File size: 14.0 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 c.SetInt( "ZSHRINK", 1, "Catalog shrink factor");
227 c.End();
228
229 c.WriteHeader(*this);
230
231 const off_t here_I_am = tellp();
232
233 //go after the catalog to compress and write the table data
234 seekp(here_I_am + catalog_size);
235
236 //calculate RAWSUM
237 Checksum rawsum;
238 rawsum.add((char*)(fOffsetCalibration.data()), 1024*1440*sizeof(int16_t));
239#if GCC_VERSION < 40603
240 c.SetStr("RAWSUM", std::to_string((long long unsigned int)(rawsum.val())));
241#else
242 c.SetStr("RAWSUM", std::to_string(rawsum.val()));
243#endif
244
245 //compress data and calculate final, compressed size
246 const uint32_t compressed_header_size = sizeof(FITS::TileHeader) + sizeof(FITS::BlockHeader) + 1*sizeof(uint16_t);
247 std::vector<char> compressed_calib(1024*1440*2 + compressed_header_size + 8); //+8 for checksum;
248 char* data_start = compressed_calib.data() + compressed_header_size;
249 uint32_t compressed_size = compressHUFFMAN16(data_start, (char*)(fOffsetCalibration.data()), 1024*1440, 2, 1);;
250 compressed_size += compressed_header_size;
251
252 //Write tile header
253 FITS::TileHeader th(0, 0);
254 std::vector<uint16_t> seq(1, FITS::kFactHuffman16);
255 FITS::Compression bh(seq, FITS::kOrderByRow);
256 th.numRows = 1;
257 th.size = compressed_size;
258 bh.SetBlockSize(compressed_size-sizeof(FITS::TileHeader));
259 memcpy(compressed_calib.data(), &(th), sizeof(FITS::TileHeader));
260 bh.Memcpy(compressed_calib.data()+sizeof(FITS::TileHeader));
261
262 //calculate resulting compressed datasum
263 Checksum datasum;
264 memset(compressed_calib.data()+compressed_size, 0, 8-compressed_size%8);
265 datasum.add(compressed_calib.data(), compressed_size + 8-compressed_size%8);
266
267 //write the catalog !
268 seekp(here_I_am);
269
270 std::vector<uint64_t> catalog(2,0);
271 catalog[0] = compressed_size-sizeof(FITS::TileHeader);
272 catalog[1] = sizeof(FITS::TileHeader);
273
274 std::vector<char> swappedCatalog(catalog_size);
275 revcpy<sizeof(int64_t)>(swappedCatalog.data(), (char*)(catalog.data()), 2);//catalog_size);
276 datasum.add(swappedCatalog.data(), catalog_size);
277
278 write(swappedCatalog.data(), catalog_size);
279
280 //update relevant keywords
281 c.SetFloat("ZRATIO", (float)(1024*1440*2)/(float)(compressed_size));
282 c.SetInt("PCOUNT", compressed_size + catalog_size);
283
284#if GCC_VERSION < 40603
285 c.SetStr("DATASUM", std::to_string((long long unsigned int)(datasum.val())));
286#else
287 c.SetStr("DATASUM", std::to_string(datasum.val()));
288#endif
289
290 datasum += c.WriteHeader(*this);
291
292 c.SetStr("CHECKSUM", datasum.str());
293
294 c.WriteHeader(*this);
295
296 //write the compressed data
297 seekp(here_I_am + catalog_size);
298 write(compressed_calib.data(), compressed_size);
299
300 AlignTo2880Bytes();
301
302 return good();
303 }
304 ///Apply the drs offset calibration (overload of super-method)
305 virtual void DrsOffsetCalibrate(char* target_location)
306 {
307 if (!IsOffsetCalibration())
308 return;
309
310 const int16_t* startCell = reinterpret_cast<int16_t*>(target_location + fStartCellsOffset);
311 int16_t* data = reinterpret_cast<int16_t*>(target_location + fDataOffset);
312
313 for (uint32_t ch=0; ch<1440; ch++)
314 {
315 if (startCell[ch] < 0)
316 {
317 data += fNumSlices;
318 continue;
319 }
320
321 const int16_t modStart = startCell[ch]%1024;
322 const int16_t *off = fOffsetCalibration.data() + ch*1024;
323
324 const int16_t* cal = off+modStart;
325 const int16_t* end_stride = data+fNumSlices;
326
327 if (modStart+fNumSlices > 1024)
328 {
329 while (cal < off+1024)
330 *data++ -= *cal++;
331 cal = off;
332 }
333
334 while (data<end_stride)
335 *data++ -= *cal++;
336 }
337 }
338
339private:
340 /// Checks if the size of the input calibration is ok
341 bool VerifyCalibrationSize(uint32_t size)
342 {
343 if (size == 1440*1024)
344 return true;
345
346 std::ostringstream str;
347 str << "Cannot load calibration with anything else than 1440 pixels and 1024 samples per pixel. Got a total size of " << size;
348#ifdef __EXCEPTIONS
349 throw std::runtime_error(str.str());
350#else
351 gLog << ___err___ << "ERROR - " << str.str() << std::endl;
352 return false;
353#endif
354 }
355
356 //Offsets calibration stuff.
357 std::vector<int16_t> fOffsetCalibration; ///< The calibration itself
358 int32_t fStartCellsOffset; ///< Offset in bytes for the startcell data
359 int32_t fDataOffset; ///< Offset in bytes for the data
360 int32_t fNumSlices; ///< Number of samples per pixel per event
361
362}; //class factofits
363
364#endif /* FACTOFITS_H_ */
Note: See TracBrowser for help on using the repository browser.