source: trunk/FACT++/src/corruptFileFixer.cc@ 19285

Last change on this file since 19285 was 18544, checked in by lyard, 9 years ago
Added -ugly- program to recover raw data files corrupted due to computer crash
File size: 12.0 KB
Line 
1/*
2 * corruptFileFixer.cc
3 *
4 * Takes a compressed fits file as an input, and copy the "good" part of it to a <filename>_recovered.fits.fz
5 * No tstart or tstop nor checksums are updated. For this the script fixHeaderKeys.sh should be used along with the ftools
6 *
7 * How to compile me:
8 *
9 * setenv PATH /gpfs/fact/swdev/root_v5.32.00/bin:/usr/lib64/qt-3.3/bin:/usr/local/bin:/bin:/usr/bin:/opt/dell/srvadmin/bin
10 *
11 * setenv LD_LIBRARY_PATH /lib64:/home/isdc/lyard/FACT++/.libs:/gpfs/fact/swdev/root_v5.32.00/lib:/swdev_nfs/FACT++/.libs:
12 *
13 * g++ -o corruptZFitsFixer --std=c++0x src/corruptFileFixer.cc -I./externals -DHAVE_ZLIB
14 *
15 *
16 * Created on: Aug 23, 2016
17 * Author: lyard
18 */
19
20#include <iostream>
21#include <fstream>
22#include <list>
23
24#include "factfits.h"
25
26using namespace std;
27
28
29 template<size_t N>
30 void revcpy(char *dest, const char *src, int num)
31 {
32 const char *pend = src + num*N;
33 for (const char *ptr = src; ptr<pend; ptr+=N, dest+=N)
34 std::reverse_copy(ptr, ptr+N, dest);
35 }
36
37 string format_integer(unsigned long value)
38 {
39 ostringstream str;
40 str << " ";
41 if (value < 10000000000) str << " ";
42 if (value < 1000000000) str << " ";
43 if (value < 100000000) str << " ";
44 if (value < 10000000) str << " ";
45 if (value < 1000000) str << " ";
46 if (value < 100000) str << " ";
47 if (value < 10000) str << " ";
48 if (value < 1000) str << " ";
49 if (value < 100) str << " ";
50 if (value < 10) str << " ";
51 str << value << " ";
52 return str.str();
53 }
54
55
56int main(int argc, char** argv)
57{
58 if (argc != 2)
59 {
60 cout << "Only one argument please: the input filename" << endl;
61 return -1;
62 }
63
64 ifstream input(argv[1], ios::in | ios::binary);
65
66 if (!input)
67 {
68 cout << "Impossible to open file named \"" << argv[1] << "\"" << endl;
69 return -2;
70 }
71
72 //First look for the calibration table and make sure that it is complete otherwise there is nothing that we can do
73 //read data by 80 chars chunks
74 char* fits_buffer = new char[81];
75 fits_buffer[80] = 0;
76 unsigned int n_lines_read = 0;
77 unsigned int num_start_tables = 0;
78 unsigned int num_end_tables = 0;
79 streampos start_of_table_header = 0;
80
81 //read fits "rows" until we arrive at the end of the data table header, i.e. 3 END and 2 XTENSION
82 while (num_start_tables != 2 || num_end_tables != 3)
83 {
84 input.read(fits_buffer, 80);
85 if (input.eof()) break;
86 n_lines_read++;
87 if (!strcmp(fits_buffer, "XTENSION= 'BINTABLE' / binary table extension "))
88 {
89 num_start_tables ++;
90 if (num_start_tables == 2)
91 start_of_table_header = input.tellg() - streampos(80);
92 }
93 if (!strcmp(fits_buffer, "END "))
94 num_end_tables++;
95 //look for relevant header keywords to be updated and remember their address in the file
96 if (num_start_tables == 2)
97 {//display the broken table header, for informatio
98 //TODO make sure that this is a data run, otherwise constants used below will not work
99 cout << fits_buffer << endl;
100 }
101 }
102
103 if (num_start_tables != 2 || num_end_tables != 3)
104 {
105 cout << "Could not reach the end of the data table header: nothing could be recovered... sorry." << endl;
106 return -2;
107 }
108
109 //we are now at the end of the data table header. Progress until the start of the data, i.e. move to 36 - (n_lines_read%36)
110 int i_end = (36-(n_lines_read%36));
111
112 for (int i=0;i<i_end; i++)
113 {
114 input.read(fits_buffer, 80);
115 if (input.eof()) break;
116 n_lines_read++;
117 }
118
119 streampos catalog_beginning = input.tellg();
120 size_t catalog_reserved_size = 4800000;
121 //we passed the header padding: skip the data and reserved catalog space
122 input.seekg(catalog_beginning + (streampos)(catalog_reserved_size));
123 n_lines_read += 60000;
124
125 //get a list of valid tiles
126 list<FITS::TileHeader> good_tiles;
127 //remember where the data starts
128 streampos data_beginning = input.tellg();
129 streampos previous_tile_begin = input.tellg();
130 //read the very first tile header
131 input.read(fits_buffer, 80);
132
133 if (memcmp(fits_buffer, "TILE", 4))
134 {
135 cout << "Compressed data does not start by string TILE: nothing could be recovered, sorry." << endl;
136 return -4;
137 }
138
139 FITS::TileHeader* thead = reinterpret_cast<FITS::TileHeader*>(fits_buffer);
140 FITS::TileHeader previous_tile = *thead;
141 //previous_tile has the previous tile header. Move forward, and remember the previous valid tile every time we find a new valid header
142 while (!input.eof())
143 {
144 //skip previous tile
145 input.seekg(previous_tile_begin + (streampos)(thead->size));
146 streampos this_tile_start = input.tellg();
147 input.read(fits_buffer, 80);
148 if (input.eof())
149 break;
150 if (fits_buffer[0] != 'T' || fits_buffer[1] != 'I' || fits_buffer[2] != 'L' || fits_buffer[3] != 'E')
151 break;
152
153 //we've found a new valid tile header: remember the previous one
154 good_tiles.push_back(previous_tile);
155 previous_tile = *thead;
156 previous_tile_begin = this_tile_start;
157 }
158
159 unsigned int num_tiles_recovered = good_tiles.size();
160 //done.
161 cout << "We have found " << num_tiles_recovered << " valid tiles. Recovering catalog now " << endl;
162
163 std::vector<std::vector<std::pair<int64_t, int64_t> > > catalog;
164
165 FITS::BlockHeader bhead;
166
167 streamoff offset_in_heap = 0;
168 streamoff valid_data = 0;
169 unsigned int num_cols = 9;
170 input.close();
171 input.open(argv[1], ios::in | ios::binary);
172 input.seekg(data_beginning);
173 for (unsigned int i=0;i<num_tiles_recovered;i++)
174 {
175 input.read(fits_buffer, sizeof(FITS::TileHeader));
176 if (!input.good())
177 break;
178 if (memcmp(thead->id, "TILE", 4))
179 break;
180
181 catalog.emplace_back();
182 offset_in_heap += sizeof(FITS::TileHeader);
183
184 //skip through the columns
185 for (unsigned int i=0;i<num_cols;i++)
186 {
187 input.read((char*)(&bhead), sizeof(FITS::BlockHeader));
188 if (!input.good())
189 {
190 break;
191 }
192 catalog.back().emplace_back((int64_t)(bhead.size), offset_in_heap);
193 offset_in_heap += bhead.size;
194 input.seekg(data_beginning + offset_in_heap);
195 }
196
197 //at the very last, 0 size time-something column
198 catalog.back().emplace_back(0,0);
199
200 if (!input.good())
201 {
202 catalog.pop_back();
203 break;
204 }
205 valid_data = offset_in_heap;
206 }
207
208 if (catalog.size() != num_tiles_recovered)
209 cout << "Notice: some apparently OK tiles are in fact corrupted: could only recover " << catalog.size() << " tiles." << endl;
210
211 string recovered_filename(argv[1]);
212 recovered_filename = recovered_filename + ".recovered";
213
214 cout << "Catalog recovered. Now writing " << recovered_filename << endl;
215
216 ifstream test_input(recovered_filename.c_str(), ios::in | ios::binary);
217
218 if (test_input)
219 {
220 cout << "Error: output file already exists. Aborting. " << endl;
221 return -10;
222 }
223
224 ofstream output(recovered_filename.c_str(), ios::out | ios::binary);
225 input.close();
226 input.open(argv[1], ios::in | ios::binary);
227
228 cout << "Writing calibration table...";
229 cout.flush();
230 //first skip the early section of the file: basic fits table, and calibration table
231 while (input.tellg() != start_of_table_header)
232 {
233 input.read(fits_buffer, 80);
234 output.write(fits_buffer, 80);
235 }
236
237 cout << "done." << endl << "Writing updated table header...";
238 cout.flush();
239
240 ostringstream updated_key;
241 //calculate updated header keys
242 unsigned long theap = 160*num_tiles_recovered;
243 unsigned long pcount = valid_data + catalog_reserved_size - theap;
244
245 while (input.tellg() != catalog_beginning)
246 {
247 input.read(fits_buffer, 80);
248 if (!memcmp(fits_buffer, "NAXIS2", 6))
249 {
250 updated_key.str("");
251 updated_key << "NAXIS2 =" << format_integer(num_tiles_recovered) << "/ number of rows in table ";
252 output.write(updated_key.str().c_str(), 80);
253 continue;
254 }
255 if (!memcmp(fits_buffer, "PCOUNT", 6))
256 {
257 updated_key.str("");
258 updated_key << "PCOUNT =" << format_integer(pcount) << "/ size of special data area ";
259 output.write(updated_key.str().c_str(), 80);
260 continue;
261 }
262 if (!memcmp(fits_buffer, "ZNAXIS2", 7))
263 {
264 updated_key.str("");
265 updated_key << "ZNAXIS2 =" << format_integer(num_tiles_recovered) << "/ Number of uncompressed rows ";
266 output.write(updated_key.str().c_str(), 80);
267 continue;
268 }
269 if (!memcmp(fits_buffer, "ZHEAPPTR", 8))
270 {
271 updated_key.str("");
272 updated_key << "ZHEAPPTR=" << format_integer(catalog_reserved_size) << " ";
273 output.write(updated_key.str().c_str(), 80);
274 continue;
275 }
276 if (!memcmp(fits_buffer, "THEAP", 5))
277 {
278 updated_key.str("");
279 updated_key << "THEAP =" << format_integer(theap) << " ";
280 output.write(updated_key.str().c_str(), 80);
281 continue;
282 }
283
284 output.write(fits_buffer, 80);
285 }
286cout << "num tiles recovered: " << num_tiles_recovered << endl;
287cout << "pcount: " << pcount << endl;
288cout << "catalog_reserved_size: " << catalog_reserved_size << endl;
289cout << "theap: " << theap << endl;
290cout << "offset in heap: " << offset_in_heap << endl;
291cout << "valid data: " << valid_data << endl;
292 cout << "done." << endl << "Writing updated catalog...";
293 cout.flush();
294
295 //write the catalog itself
296 vector<char> swapped_catalog(catalog_reserved_size);
297 unsigned int shift = 0;
298 for (auto it=catalog.cbegin(); it!=catalog.cend(); it++)
299 {
300 revcpy<sizeof(uint64_t)>(swapped_catalog.data() + shift, (char*)(it->data()), (num_cols+1)*2);
301 shift += 160;
302 }
303
304 if (catalog.size() < 30000)
305 memset(swapped_catalog.data()+shift, 0, catalog_reserved_size - shift);
306
307 output.write(swapped_catalog.data(), catalog_reserved_size);
308
309 cout << "done." << endl << "Writing recovered data...";
310 cout.flush();
311
312 //write the actual data
313 input.seekg(input.tellg() + (streampos)(catalog_reserved_size));
314cout << "starting writing data at " << output.tellp() << endl;
315 unsigned int actual_data_size = pcount + theap - catalog_reserved_size;
316 unsigned int i=0;
317 for (;i<actual_data_size-80;i+=80)
318 {
319 input.read(fits_buffer, 80);
320 output.write(fits_buffer, 80);
321 }
322
323 input.read(fits_buffer, actual_data_size%80);
324 output.write(fits_buffer, actual_data_size%80);
325cout << "ended writing data at " << output.tellp() << endl;
326 cout << "done." << endl << "Writing FITS padding...";
327 cout.flush();
328cout << "Current extra chars: " << output.tellp()%(80*36) << endl;
329cout << "Will write " << 2880 - output.tellp()%(2880) << " filling bytes while at byte " << output.tellp() << endl;
330 //eventually write the fits padding
331 if (output.tellp()%(80*36) > 0)
332 {
333 std::vector<char> filler(2880-output.tellp()%(2880), 0);
334 output.write(filler.data(), filler.size());
335 }
336 output.close();
337 cout << "All done !." << endl;
338 cout << "TSTOP is also most likely invalid, so:" << endl;
339 cout << " - run /swdev_nfs/FACT++/fitsdump <filename> -c UnixTimeUTC --minmax --nozero" << endl;
340 cout << " - update header key using e.g. fv" << endl;
341 cout << "Checksums are now invalid: please run fchecksum <filename> update+ datasum+ " << endl;
342 return 0;
343}
344
345
346
Note: See TracBrowser for help on using the repository browser.