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

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