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

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