source: trunk/FACT++/src/fitsdump.cc@ 20115

Last change on this file since 20115 was 20092, checked in by tbretz, 4 years ago
Made consitsent return codes.
File size: 39.6 KB
Line 
1//****************************************************************
2//****************************************************************
3/** @class FitsDumper
4
5 @brief Dumps contents of fits tables to stdout or a file
6
7 */
8 //****************************************************************
9#include "Configuration.h"
10
11#include <float.h>
12
13#include <map>
14#include <fstream>
15
16#include <boost/regex.hpp>
17
18#include "tools.h"
19#include "Time.h"
20#include "factfits.h"
21
22#ifdef HAVE_ROOT
23#include "TFile.h"
24#include "TTree.h"
25#include "TFormula.h"
26#endif
27
28using namespace std;
29
30struct MyColumn
31{
32 string name;
33
34 fits::Table::Column col;
35
36 uint32_t first;
37 uint32_t last;
38
39 void *ptr;
40};
41
42struct minMaxStruct
43{
44 double min;
45 double max;
46 long double average;
47 long double squared;
48 long numValues;
49 minMaxStruct() : min(FLT_MAX), max(-FLT_MAX), average(0), squared(0), numValues(0) { }
50
51 void add(long double val)
52 {
53 average += val;
54 squared += val*val;
55
56 if (val<min)
57 min = val;
58
59 if (val>max)
60 max = val;
61
62 numValues++;
63 }
64};
65
66
67class FitsDumper : public factfits
68{
69private:
70 string fFilename;
71 bool fStream;
72
73 // Convert CCfits::ValueType into a human readable string
74 string ValueTypeToStr(char type) const;
75
76 /// Lists all columns of an open file
77 void List();
78 void ListFileContent() const;
79 void ListHeader(const string& filename);
80 void ListKeywords(ostream &);
81
82 vector<MyColumn> InitColumns(vector<string> list);
83 vector<MyColumn> InitColumnsRoot(vector<string> &list);
84
85 double GetDouble(const MyColumn &, size_t) const;
86 int64_t GetInteger(const MyColumn &, size_t) const;
87 string Format(const string &fmt, const double &val) const;
88 string Format(const string &fmt, const MyColumn &, size_t) const;
89
90 ///Display the selected columns values VS time
91 int Dump(ostream &, const vector<string> &, const vector<MyColumn> &, const string &, size_t, size_t, const string &, Configuration &);
92 void DumpRoot(ostream &, const vector<string> &, const string &, size_t, size_t, const string &);
93 void DumpMinMax(ostream &, const vector<MyColumn> &, size_t, size_t, bool);
94 void DumpStats(ostream &, const vector<MyColumn> &, const string &, size_t, size_t);
95
96public:
97 FitsDumper(const string &fname, const string &tablename);
98
99 ///Configures the fitsLoader from the config file and/or command arguments.
100 int Exec(Configuration& conf);
101};
102
103// --------------------------------------------------------------------------
104//
105//! Constructor
106//! @param out
107//! the ostream where to redirect the outputs
108//
109FitsDumper::FitsDumper(const string &fname, const string &tablename) : factfits(fname, tablename), fFilename(fname), fStream(false)
110{
111}
112
113string FitsDumper::ValueTypeToStr(char type) const
114{
115 switch (type)
116 {
117 case 'L': return "bool(8)";
118 case 'A': return "char(8)";
119 case 'B': return "byte(8)";
120 case 'I': return "short(16)";
121 case 'J': return "int(32)";
122 case 'K': return "int(64)";
123 case 'E': return "float(32)";
124 case 'D': return "double(64)";
125 default:
126 return "unknown";
127 }
128}
129
130void FitsDumper::List()
131{
132 const fits::Table::Columns &fColMap = GetColumns();
133
134 cout << endl;
135
136 cout << "\nFile: " << fFilename << "\n";
137
138 cout << " " << Get<std::string>("EXTNAME", "");
139 cout << " [" << GetNumRows() << "]\n";
140
141 for (auto it=fColMap.cbegin(); it!=fColMap.cend(); it++)
142 {
143 cout << " " << setw(3) << it->second.id;
144 cout << "|" << setw(10) << ValueTypeToStr(it->second.type);
145 cout << " " << it->first << "[" << it->second.num << "]";
146 cout << " / " << it->second.comment;
147 cout << " (" << it->second.unit << ")" << endl;
148 }
149
150 cout << endl;
151}
152
153void FitsDumper::ListKeywords(ostream &fout)
154{
155 const fits::Table::Keys &fKeyMap = GetKeys();
156
157 for (auto it=fKeyMap.begin(); it != fKeyMap.end(); it++)
158 {
159 fout << "## " << ::left << setw(8) << it->first << "= ";
160
161 if (it->second.type=='T')
162 fout << ::left << setw(20) << ("'"+it->second.value+"'");
163 else
164 fout << ::right << setw(20) << it->second.value;
165
166 if (!it->second.comment.empty())
167 fout << " / " << it->second.comment;
168 fout << '\n';
169 }
170
171 fout << flush;
172}
173
174void FitsDumper::ListFileContent() const
175{
176 const std::vector<std::string> &tables = GetTables();
177
178 cout << "File " << fFilename << " has " << tables.size() << " table(s): " << endl;
179 for (auto it=tables.begin(); it!=tables.end(); it++)
180 cout << " * " << *it << endl;
181}
182
183void FitsDumper::ListHeader(const string& filename)
184{
185 ostream fout(cout.rdbuf());
186
187 ofstream sout;
188 if (filename!="-")
189 {
190 sout.open(filename);
191 if (!sout)
192 {
193 cerr << "Cannot open output stream " << filename << ": " << strerror(errno) << endl;
194 return;
195 }
196 fout.rdbuf(sout.rdbuf());
197 }
198
199 const fits::Table::Keys &fKeyMap = GetKeys();
200
201 fout << "\nTable: " << fKeyMap.find("EXTNAME")->second.value << " (rows=" << GetNumRows() << ")\n";
202 if (fKeyMap.find("COMMENT") != fKeyMap.end())
203 fout << "Comment: \t" << fKeyMap.find("COMMENT")->second.value << "\n";
204
205 ListKeywords(fout);
206 fout << endl;
207}
208
209vector<MyColumn> FitsDumper::InitColumns(vector<string> names)
210{
211 static const boost::regex expr("([[:word:].]+)(\\[([[:digit:]]+)?(:)?([[:digit:]]+)?\\])?");
212
213 const fits::Table::Columns &fColMap = GetColumns();
214
215 if (names.empty())
216 for (auto it=fColMap.begin(); it!=fColMap.end(); it++)
217 if (it->second.num>0)
218 names.push_back(it->first);
219
220 vector<MyColumn> vec;
221
222 for (auto it=names.begin(); it!=names.end(); it++)
223 {
224 boost::smatch what;
225 if (!boost::regex_match(*it, what, expr, boost::match_extra))
226 {
227 cerr << "Couldn't parse expression '" << *it << "' " << endl;
228 return vector<MyColumn>();
229 }
230
231 const string name = what[1];
232
233 const auto iter = fColMap.find(name);
234 if (iter==fColMap.end())
235 {
236 cerr << "ERROR - Column '" << name << "' not found in table." << endl;
237 return vector<MyColumn>();
238 }
239
240 const fits::Table::Column &col = iter->second;
241
242 const string val0 = what[3];
243 const string delim = what[4];
244 const string val1 = what[5];
245
246 const uint32_t first = atol(val0.c_str());
247 const uint32_t last = (val0.empty() && delim.empty()) ? col.num-1 : (val1.empty() ? first : atoi(val1.c_str()));
248
249 if (first>=col.num)
250 {
251 cerr << "ERROR - First index " << first << " for column " << name << " exceeds number of elements " << col.num << endl;
252 return vector<MyColumn>();
253 }
254
255 if (last>=col.num)
256 {
257 cerr << "ERROR - Last index " << last << " for column " << name << " exceeds number of elements " << col.num << endl;
258 return vector<MyColumn>();
259 }
260
261 if (first>last)
262 {
263 cerr << "ERROR - Last index " << last << " for column " << name << " exceeds first index " << first << endl;
264 return vector<MyColumn>();
265 }
266
267 MyColumn mycol;
268
269 mycol.name = name;
270 mycol.col = col;
271 mycol.first = first;
272 mycol.last = last;
273 mycol.ptr = SetPtrAddress(name);
274
275 vec.push_back(mycol);
276 }
277
278 return vec;
279}
280
281double FitsDumper::GetDouble(const MyColumn &it, size_t i) const
282{
283 switch (it.col.type)
284 {
285 case 'A':
286 return reinterpret_cast<const char*>(it.ptr)[i];
287
288 case 'L':
289 return reinterpret_cast<const bool*>(it.ptr)[i];
290
291 case 'B':
292 return (unsigned int)reinterpret_cast<const uint8_t*>(it.ptr)[i];
293
294 case 'I':
295 return reinterpret_cast<const int16_t*>(it.ptr)[i];
296
297 case 'J':
298 return reinterpret_cast<const int32_t*>(it.ptr)[i];
299
300 case 'K':
301 return reinterpret_cast<const int64_t*>(it.ptr)[i];
302
303 case 'E':
304 return reinterpret_cast<const float*>(it.ptr)[i];
305
306 case 'D':
307 return reinterpret_cast<const double*>(it.ptr)[i];
308 }
309
310 return 0;
311}
312
313int64_t FitsDumper::GetInteger(const MyColumn &it, size_t i) const
314{
315 switch (it.col.type)
316 {
317 case 'A':
318 return reinterpret_cast<const char*>(it.ptr)[i];
319
320 case 'L':
321 return reinterpret_cast<const bool*>(it.ptr)[i];
322
323 case 'B':
324 return (unsigned int)reinterpret_cast<const uint8_t*>(it.ptr)[i];
325
326 case 'I':
327 return reinterpret_cast<const int16_t*>(it.ptr)[i];
328
329 case 'J':
330 return reinterpret_cast<const int32_t*>(it.ptr)[i];
331
332 case 'K':
333 return reinterpret_cast<const int64_t*>(it.ptr)[i];
334
335 case 'E':
336 return reinterpret_cast<const float*>(it.ptr)[i];
337
338 case 'D':
339 return reinterpret_cast<const double*>(it.ptr)[i];
340 }
341
342 return 0;
343}
344
345#pragma GCC diagnostic push
346#pragma GCC diagnostic ignored "-Wformat-nonliteral"
347string FitsDumper::Format(const string &format, const MyColumn &col, size_t i) const
348{
349 switch (*format.rbegin())
350 {
351 case 'd':
352 case 'i':
353 case 'o':
354 case 'u':
355 case 'x':
356 case 'X':
357 return Tools::Form(format.c_str(), GetDouble(col, i));
358
359 case 'e':
360 case 'E':
361 case 'f':
362 case 'F':
363 case 'g':
364 case 'G':
365 case 'a':
366 case 'A':
367 return Tools::Form(format.c_str(), GetInteger(col, i));
368
369 case 'h':
370 {
371 string rc = Tools::Scientific(GetDouble(col, i));
372 *remove_if(rc.begin(), rc.end(), ::isspace)=0;
373 return rc;
374 }
375 }
376
377 return "";
378}
379
380string FitsDumper::Format(const string &format, const double &val) const
381{
382 switch (*format.rbegin())
383 {
384 case 'd':
385 case 'i':
386 case 'o':
387 case 'u':
388 case 'x':
389 case 'X':
390 return Tools::Form(format.c_str(), int64_t(val));
391
392 case 'e':
393 case 'E':
394 case 'f':
395 case 'F':
396 case 'g':
397 case 'G':
398 case 'a':
399 case 'A':
400 return Tools::Form(format.c_str(), val);
401
402 case 'h':
403 {
404 string rc = Tools::Scientific(val);
405 *remove_if(rc.begin(), rc.end(), ::isspace)=0;
406 return rc;
407 }
408 }
409
410 return "";
411}
412#pragma GCC diagnostic pop
413
414
415#ifdef HAVE_ROOT
416template<typename T>
417void SetBranch(TTree *ttree, const string &bname, const MyColumn &col)
418{
419 TBranch *b = ttree->GetBranch(bname.c_str());
420
421 T *address = reinterpret_cast<T*>(col.ptr)+col.first;
422 if (b)
423 b->SetAddress(address);
424 else
425 ttree->Branch(bname.c_str(), address);
426}
427#endif
428
429// --------------------------------------------------------------------------
430//
431//! Perform the actual dump, based on the current parameters
432//
433#ifdef HAVE_ROOT
434int FitsDumper::Dump(ostream &fout, const vector<string> &format, const vector<MyColumn> &cols, const string &filter, size_t first, size_t limit, const string &filename, Configuration &conf)
435#else
436int FitsDumper::Dump(ostream &fout, const vector<string> &format, const vector<MyColumn> &cols, const string &, size_t first, size_t limit, const string &filename, Configuration &)
437#endif
438{
439 const fits::Table::Keys &fKeyMap = GetKeys();
440
441#ifdef HAVE_ROOT
442 TFormula select;
443 if (!filter.empty() && select.Compile(filter.c_str()))
444 throw runtime_error("Syntax Error: TFormula::Compile failed for '"+filter+"'");
445#endif
446
447 fout << "## --------------------------------------------------------------------------\n";
448 fout << "## Fits file: \t" << fFilename << '\n';
449 if (filename!="-")
450 fout << "## File: \t" << filename << '\n';
451 fout << "## Table: \t" << fKeyMap.find("EXTNAME")->second.value << '\n';
452 fout << "## NumRows: \t" << GetNumRows() << '\n';
453 fout << "## Comment: \t" << ((fKeyMap.find("COMMENT") != fKeyMap.end()) ? fKeyMap.find("COMMENT")->second.value : "") << '\n';
454#ifdef HAVE_ROOT
455 if (!filter.empty())
456 fout << "## Selection: \t" << select.GetExpFormula() << '\n';
457#endif
458 fout << "## --------------------------------------------------------------------------\n";
459 ListKeywords(fout);
460 fout << "## --------------------------------------------------------------------------\n";
461 fout << "#\n";
462
463#ifdef HAVE_ROOT
464 const bool rootfile = filename.rfind(".root")==filename.size()-5;
465
466 TFile *tfile = 0;
467 if (rootfile)
468 {
469 tfile = new TFile(filename.c_str(), conf.Get<bool>("update")?"UPDATE":(conf.Get<bool>("overwrite")?"RECREATE":"CREATE"), fFilename.c_str(), conf.Get<bool>("force"));
470 if (tfile->IsZombie())
471 {
472 cerr << "Could not open root file '" << filename << "'" << endl;
473 return 8;
474 }
475 }
476
477 // FIXME: Add header values
478 TTree *ttree = 0;
479 if (tfile)
480 {
481 const string tname = fKeyMap.find("EXTNAME")->second.value;
482 tfile->GetObject(tname.c_str(), ttree);
483 if (!ttree)
484 ttree = new TTree(tname.c_str(), Tools::Form("From N=%d [cnt=%d] {filter=%s}", first, limit, filter.c_str()).c_str());
485 }
486#endif
487
488 size_t num = 0;
489 for (auto it=cols.begin(); it!=cols.end(); it++)
490 {
491 fout << "# " << it->name;
492
493 if (it->first==it->last)
494 {
495 if (it->first!=0)
496 fout << "[" << it->first << "]";
497 }
498 else
499 fout << "[" << it->first << ":" << it->last << "]";
500
501 if (!it->col.unit.empty())
502 fout << ": " << it->col.unit;
503 fout << '\n';
504
505 num += it->last-it->first+1;
506 }
507 fout << "#" << endl;
508
509 // -----------------------------------------------------------------
510
511#ifdef HAVE_ROOT
512 if (ttree)
513 {
514 fout << "## --------------------------------------------------------------------------\n";
515
516 fout << "#\n";
517 fout << "# Tree: " << fKeyMap.find("EXTNAME")->second.value << "\n";
518 fout << "#\n";
519
520 for (auto it=cols.begin(); it!=cols.end(); it++)
521 {
522 string branch = it->name;
523 if (it->last!=it->first)
524 branch += "["+to_string(it->last-it->first+1)+"]";
525
526 fout << "# Branch: " << branch << "\n";
527
528 switch (it->col.type)
529 {
530 case 'B':
531 SetBranch<unsigned char>(ttree, branch, *it);
532 break;
533 case 'L':
534 SetBranch<bool>(ttree, branch, *it);
535 break;
536 case 'I':
537 SetBranch<int16_t>(ttree, branch, *it);
538 break;
539 case 'J':
540 SetBranch<int32_t>(ttree, branch, *it);
541 break;
542 case 'K':
543 SetBranch<int64_t>(ttree, branch, *it);
544 break;
545 case 'E':
546 SetBranch<float>(ttree, branch, *it);
547 break;
548 case 'D':
549 SetBranch<double>(ttree, branch, *it);
550 break;
551 default:
552 ;
553 }
554 }
555
556 fout << "#\n";
557 fout << "## --------------------------------------------------------------------------\n";
558}
559#endif
560
561
562 // -----------------------------------------------------------------
563
564
565
566#ifdef HAVE_ROOT
567 vector<Double_t> data(num+1);
568#endif
569
570 const size_t last = limit ? first + limit : size_t(-1);
571
572 while (GetRow(first++, !fStream))
573 {
574 const size_t row = GetRow();
575 if ((!fStream && row==GetNumRows()) || row==last)
576 break;
577
578 size_t p = 0;
579
580#ifdef HAVE_ROOT
581 data[p++] = first-1;
582#endif
583
584 ostringstream sout;
585 sout.precision(fout.precision());
586 sout.flags(fout.flags());
587
588 uint32_t col = 0;
589 for (auto it=cols.begin(); it!=cols.end(); it++, col++)
590 {
591 string msg;
592 for (uint32_t i=it->first; i<=it->last; i++, p++)
593 {
594 if (col<format.size())
595 sout << Format("%"+format[col], *it, i) << " ";
596 else
597 {
598 switch (it->col.type)
599 {
600 case 'A':
601 msg += reinterpret_cast<const char*>(it->ptr)[i];
602 break;
603 case 'B':
604 sout << (unsigned int)reinterpret_cast<const unsigned char*>(it->ptr)[i] << " ";
605 break;
606 case 'L':
607 sout << reinterpret_cast<const bool*>(it->ptr)[i] << " ";
608 break;
609 case 'I':
610 sout << reinterpret_cast<const int16_t*>(it->ptr)[i] << " ";
611 break;
612 case 'J':
613 sout << reinterpret_cast<const int32_t*>(it->ptr)[i] << " ";
614 break;
615 case 'K':
616 sout << reinterpret_cast<const int64_t*>(it->ptr)[i] << " ";
617 break;
618 case 'E':
619 sout << reinterpret_cast<const float*>(it->ptr)[i] << " ";
620 break;
621 case 'D':
622 sout << reinterpret_cast<const double*>(it->ptr)[i] << " ";
623 break;
624 default:
625 ;
626 }
627 }
628#ifdef HAVE_ROOT
629 if (!filter.empty())
630 data[p] = GetDouble(*it, i);
631#endif
632 }
633
634 if (it->col.type=='A')
635 sout << "'" << msg.c_str() << "' ";
636 }
637#ifdef HAVE_ROOT
638 if (!filter.empty() && select.EvalPar(0, data.data())<0.5)
639 continue;
640
641 if (ttree)
642 {
643 ttree->Fill();
644 continue;
645 }
646#endif
647
648 fout << sout.str() << endl;
649 }
650
651#ifdef HAVE_ROOT
652 if (tfile)
653 {
654 fout << "#\n";
655 fout << "# " << ttree->GetEntries() << " rows filled into tree.\n";
656 fout << "#\n";
657 fout << "## --------------------------------------------------------------------------\n";
658 fout << endl;
659
660 ttree->Write("", TObject::kWriteDelete);
661 delete tfile;
662 }
663#endif
664
665 return 0;
666}
667
668vector<MyColumn> FitsDumper::InitColumnsRoot(vector<string> &names)
669{
670 static const boost::regex expr("[^\\[]([[:word:].]+)(\\[([[:digit:]]+)\\])?");
671
672 const fits::Table::Columns &cols = GetColumns();
673
674 vector<MyColumn> vec;
675
676 for (auto it=names.begin(); it!=names.end(); it++)
677 {
678 if (it->empty())
679 continue;
680
681 *it = ' '+*it;
682
683 string::const_iterator ibeg = it->begin();
684 string::const_iterator iend = it->end();
685
686 boost::smatch what;
687 while (boost::regex_search(ibeg, iend, what, expr, boost::match_extra))
688 {
689 const string all = what[0];
690 const string name = what[1];
691 const size_t idx = atol(string(what[3]).c_str());
692
693 // Check if found colum is valid
694 const auto ic = cols.find(name);
695 if (ic==cols.end())
696 {
697 ibeg++;
698 //cout << "Column '" << name << "' does not exist." << endl;
699 //return vector<MyColumn>();
700 continue;
701 }
702 if (idx>=ic->second.num)
703 {
704 cout << "Column '" << name << "' has no index " << idx << "." << endl;
705 return vector<MyColumn>();
706 }
707
708 // find index if column already exists
709 size_t p = 0;
710 for (; p<vec.size(); p++)
711 if (vec[p].name==name)
712 break;
713
714 const string id = '['+to_string(p)+']';
715
716 // Replace might reallocate the memory. Therefore, we cannot use what[0].first
717 // directly but have to store the offset
718 const size_t offset = what[0].first - it->begin();
719
720 it->replace(ibeg-it->begin()+what.position(1), what.length()-1, id);
721
722 ibeg = it->begin() + offset + id.size();
723 iend = it->end();
724
725 if (p<vec.size())
726 continue;
727
728 // Column not found, add new column
729 MyColumn mycol;
730
731 mycol.name = name;
732 mycol.col = ic->second;
733 mycol.first = idx;
734 mycol.last = idx;
735 mycol.ptr = SetPtrAddress(name);
736
737 vec.push_back(mycol);
738 }
739 }
740
741 ostringstream id;
742 id << '[' << vec.size() << ']';
743
744 for (auto it=names.begin(); it!=names.end(); it++)
745 {
746 while (1)
747 {
748 auto p = it->find_first_of('#');
749 if (p==string::npos)
750 break;
751
752 it->replace(p, 1, id.str());
753 }
754 }
755
756 //cout << endl;
757 //for (size_t i=0; i<vec.size(); i++)
758 // cout << "val[" << i << "] = " << vec[i].name << '[' << vec[i].first << ']' << endl;
759 //cout << endl;
760
761 return vec;
762}
763
764#ifdef HAVE_ROOT
765void FitsDumper::DumpRoot(ostream &fout, const vector<string> &cols, const string &filter, size_t first, size_t limit, const string &filename)
766#else
767void FitsDumper::DumpRoot(ostream &, const vector<string> &, const string &, size_t, size_t, const string &)
768#endif
769{
770#ifdef HAVE_ROOT
771 vector<string> names(cols);
772 names.insert(names.begin(), filter);
773
774 const vector<MyColumn> vec = InitColumnsRoot(names);
775 if (vec.empty())
776 return;
777
778 vector<TFormula> form(names.size());
779
780 auto ifo = form.begin();
781 for (auto it=names.begin(); it!=names.end(); it++, ifo++)
782 {
783 if (!it->empty() && ifo->Compile(it->c_str()))
784 throw runtime_error("Syntax Error: TFormula::Compile failed for '"+*it+"'");
785 }
786
787 const fits::Table::Keys &fKeyMap = GetKeys();
788
789 fout << "## --------------------------------------------------------------------------\n";
790 fout << "## Fits file: \t" << fFilename << '\n';
791 if (filename!="-")
792 fout << "## File: \t" << filename << '\n';
793 fout << "## Table: \t" << fKeyMap.find("EXTNAME")->second.value << '\n';
794 fout << "## NumRows: \t" << GetNumRows() << '\n';
795 fout << "## Comment: \t" << ((fKeyMap.find("COMMENT") != fKeyMap.end()) ? fKeyMap.find("COMMENT")->second.value : "") << '\n';
796 fout << "## --------------------------------------------------------------------------\n";
797 ListKeywords(fout);
798 fout << "## --------------------------------------------------------------------------\n";
799 fout << "##\n";
800 if (!filter.empty())
801 fout << "## Selection: " << form[0].GetExpFormula() << "\n##\n";
802
803 size_t num = 0;
804 for (auto it=vec.begin(); it!=vec.end(); it++, num++)
805 {
806 fout << "## [" << num << "] = " << it->name;
807
808 if (it->first==it->last)
809 {
810 if (it->first!=0)
811 fout << "[" << it->first << "]";
812 }
813 else
814 fout << "[" << it->first << ":" << it->last << "]";
815
816 if (!it->col.unit.empty())
817 fout << ": " << it->col.unit;
818 fout << '\n';
819 }
820 fout << "##\n";
821 fout << "## --------------------------------------------------------------------------\n";
822 fout << "#\n";
823
824 fout << "# ";
825 for (auto it=form.begin()+1; it!=form.end(); it++)
826 fout << " \"" << it->GetExpFormula() << "\"";
827 fout << "\n#" << endl;
828
829 // -----------------------------------------------------------------
830
831 vector<Double_t> data(vec.size()+1);
832
833 const size_t last = limit ? first + limit : size_t(-1);
834
835 while (GetRow(first++, !fStream))
836 {
837 const size_t row = GetRow();
838 if ((!fStream && row==GetNumRows()) || row==last)
839 break;
840
841 size_t p = 0;
842 for (auto it=vec.begin(); it!=vec.end(); it++, p++)
843 data[p] = GetDouble(*it, it->first);
844
845 data[p] = first;
846
847 if (!filter.empty() && form[0].EvalPar(0, data.data())<0.5)
848 continue;
849
850 for (auto iform=form.begin()+1; iform!=form.end(); iform++)
851 fout << iform->EvalPar(0, data.data()) << " ";
852
853 fout << endl;
854 }
855#endif
856}
857
858
859void FitsDumper::DumpMinMax(ostream &fout, const vector<MyColumn> &cols, size_t first, size_t limit, bool fNoZeroPlease)
860{
861 vector<minMaxStruct> statData(cols.size());
862
863 // Loop over all columns in our list of requested columns
864 const size_t last = limit ? first + limit : size_t(-1);
865
866 while (GetRow(first++, !fStream))
867 {
868 const size_t row = GetRow();
869 if ((!fStream && row==GetNumRows()) || row==last)
870 break;
871
872 auto statsIt = statData.begin();
873
874 for (auto it=cols.begin(); it!=cols.end(); it++, statsIt++)
875 {
876 if ((it->name=="UnixTimeUTC" || it->name=="PCTime") && it->first==0 && it->last==1)
877 {
878 const uint32_t *val = reinterpret_cast<const uint32_t*>(it->ptr);
879 if (fNoZeroPlease && val[0]==0 && val[1]==0)
880 continue;
881
882 statsIt->add(Time(val[0], val[1]).Mjd());
883 continue;
884 }
885
886 for (uint32_t i=it->first; i<=it->last; i++)
887 {
888 const double cValue = GetDouble(*it, i);
889
890 if (fNoZeroPlease && cValue == 0)
891 continue;
892
893 statsIt->add(cValue);
894 }
895 }
896 }
897
898 // okay. So now I've got ALL the data, loaded.
899 // let's do the summing and averaging in a safe way (i.e. avoid overflow
900 // of variables as much as possible)
901 auto statsIt = statData.begin();
902 for (auto it=cols.begin(); it!=cols.end(); it++, statsIt++)
903 {
904 fout << "\n[" << it->name << ':' << it->first;
905 if (it->first!=it->last)
906 fout << ':' << it->last;
907 fout << "]\n";
908
909 if (statsIt->numValues == 0)
910 {
911 fout << "Min: -\nMax: -\nAvg: -\nRms: -" << endl;
912 continue;
913 }
914
915 const long &num = statsIt->numValues;
916
917 long double &avg = statsIt->average;
918 long double &rms = statsIt->squared;
919
920 avg /= num;
921 rms /= num;
922 rms += avg*avg;
923 rms = rms<0 ? 0 : sqrt(rms);
924
925 fout << "Min: " << statsIt->min << '\n';
926 fout << "Max: " << statsIt->max << '\n';
927 fout << "Avg: " << avg << '\n';
928 fout << "Rms: " << rms << endl;
929 }
930}
931
932template<typename T>
933void displayStats(vector<char> &array, ostream& out)
934{
935 const size_t numElems = array.size()/sizeof(T);
936 if (numElems == 0)
937 {
938 out << "Min: -\nMax: -\nMed: -\nAvg: -\nRms: -" << endl;
939 return;
940 }
941
942 T *val = reinterpret_cast<T*>(array.data());
943
944 sort(val, val+numElems);
945
946 out << "Min: " << double(val[0]) << '\n';
947 out << "Max: " << double(val[numElems-1]) << '\n';
948
949 if (numElems%2 == 0)
950 out << "Med: " << (double(val[numElems/2-1]) + double(val[numElems/2]))/2 << '\n';
951 else
952 out << "Med: " << double(val[numElems/2]) << '\n';
953
954 long double avg = 0;
955 long double rms = 0;
956 for (uint32_t i=0;i<numElems;i++)
957 {
958 const long double v = val[i];
959 avg += v;
960 rms += v*v;
961 }
962
963 avg /= numElems;
964 rms /= numElems;
965 rms -= avg*avg;
966 rms = rms<0 ? 0 : sqrt(rms);
967
968
969 out << "Avg: " << avg << '\n';
970 out << "Rms: " << rms << endl;
971}
972
973#ifdef HAVE_ROOT
974void FitsDumper::DumpStats(ostream &fout, const vector<MyColumn> &cols, const string &filter, size_t first, size_t limit)
975#else
976void FitsDumper::DumpStats(ostream &fout, const vector<MyColumn> &cols, const string &, size_t first, size_t limit)
977#endif
978{
979#ifdef HAVE_ROOT
980 TFormula select;
981 if (!filter.empty() && select.Compile(filter.c_str()))
982 throw runtime_error("Syntax Error: TFormula::Compile failed for '"+filter+"'");
983#endif
984
985 // Loop over all columns in our list of requested columns
986 vector<vector<char>> statData;
987
988 const size_t rows = limit==0 || GetNumRows()<limit ? GetNumRows() : limit;
989
990 for (auto it=cols.begin(); it!=cols.end(); it++)
991 statData.emplace_back(vector<char>(it->col.size*rows*(it->last-it->first+1)));
992
993#ifdef HAVE_ROOT
994 size_t num = 0;
995 for (auto it=cols.begin(); it!=cols.end(); it++)
996 num += it->last-it->first+1;
997
998 vector<Double_t> data(num+1);
999#endif
1000
1001 // Loop over all columns in our list of requested columns
1002 const size_t last = limit ? first + limit : size_t(-1);
1003
1004 uint64_t counter = 0;
1005
1006 while (GetRow(first++, !fStream))
1007 {
1008 const size_t row = GetRow();
1009 if ((!fStream && row==GetNumRows()) || row==last)
1010 break;
1011
1012#ifdef HAVE_ROOT
1013 if (!filter.empty())
1014 {
1015 size_t p = 0;
1016
1017 data[p++] = first-1;
1018
1019 for (auto it=cols.begin(); it!=cols.end(); it++)
1020 for (uint32_t i=it->first; i<=it->last; i++, p++)
1021 data[p] = GetDouble(*it, i);
1022
1023 if (select.EvalPar(0, data.data())<0.5)
1024 continue;
1025 }
1026#endif
1027
1028 auto statsIt = statData.begin();
1029 for (auto it=cols.begin(); it!=cols.end(); it++, statsIt++)
1030 {
1031 const char *src = reinterpret_cast<const char*>(it->ptr);
1032 const size_t sz = (it->last-it->first+1)*it->col.size;
1033 memcpy(statsIt->data()+counter*sz, src+it->first*it->col.size, sz);
1034 }
1035
1036 counter++;
1037 }
1038
1039 auto statsIt = statData.begin();
1040 for (auto it=cols.begin(); it!=cols.end(); it++, statsIt++)
1041 {
1042 fout << "\n[" << it->name << ':' << it->first;
1043 if (it->last!=it->first)
1044 fout << ':' << it->last;
1045 fout << "]\n";
1046
1047 const size_t sz = (it->last-it->first+1)*it->col.size;
1048 statsIt->resize(counter*sz);
1049
1050 switch (it->col.type)
1051 {
1052 case 'L':
1053 displayStats<bool>(*statsIt, fout);
1054 break;
1055 case 'B':
1056 displayStats<char>(*statsIt, fout);
1057 break;
1058 case 'I':
1059 displayStats<int16_t>(*statsIt, fout);
1060 break;
1061 case 'J':
1062 displayStats<int32_t>(*statsIt, fout);
1063 break;
1064 case 'K':
1065 displayStats<int64_t>(*statsIt, fout);
1066 break;
1067 case 'E':
1068 displayStats<float>(*statsIt, fout);
1069 break;
1070 case 'D':
1071 displayStats<double>(*statsIt, fout);
1072 break;
1073 default:
1074 ;
1075 }
1076 }
1077}
1078
1079// --------------------------------------------------------------------------
1080//
1081//! Retrieves the configuration parameters
1082//! @param conf
1083//! the configuration object
1084//
1085int FitsDumper::Exec(Configuration& conf)
1086{
1087 fStream = conf.Get<bool>("stream");
1088
1089 if (conf.Get<bool>("autocal"))
1090 resetCalibration();
1091
1092 if (conf.Get<bool>("list"))
1093 List();
1094
1095 if (conf.Get<bool>("filecontent"))
1096 ListFileContent();
1097
1098 if (conf.Get<bool>("header"))
1099 ListHeader(conf.Get<string>("outfile"));
1100
1101
1102 if (conf.Get<bool>("header") || conf.Get<bool>("list") || conf.Get<bool>("filecontent"))
1103 return 0;
1104
1105 // ------------------------------------------------------------
1106
1107 if (conf.Get<bool>("minmax") && conf.Get<bool>("stat"))
1108 {
1109 cerr << "Invalid combination of options: cannot do stats and minmax." << endl;
1110 return 1;
1111 }
1112 if (conf.Get<bool>("stat") && conf.Get<bool>("nozero"))
1113 {
1114 cerr << "Invalid combination of options: nozero only works with minmax." << endl;
1115 return 2;
1116 }
1117
1118 if (conf.Get<bool>("scientific") && conf.Get<bool>("fixed"))
1119 {
1120 cerr << "Switched --scientific and --fixed are mutually exclusive." << endl;
1121 return 3;
1122 }
1123
1124 if (conf.Has("%") && conf.Has("%%"))
1125 {
1126 cerr << "Switched --% and --%% are mutually exclusive." << endl;
1127 return 4;
1128 }
1129
1130 // ------------------------------------------------------------
1131
1132 string filename = conf.Get<string>("outfile");
1133
1134 if (!filename.empty() && filename[0]=='.')
1135 {
1136 if (fFilename.rfind(".fits")==fFilename.size()-5)
1137 filename.insert(0, fFilename.substr(0, fFilename.size()-5));
1138 if (fFilename.rfind(".fits.fz")==fFilename.size()-8)
1139 filename.insert(0, fFilename.substr(0, fFilename.size()-8));
1140 }
1141
1142 ostream fout(cout.rdbuf());
1143
1144 ofstream sout;
1145 if (filename!="-" && filename.rfind(".root")!=filename.length()-5)
1146 {
1147 sout.open(filename);
1148 if (!sout)
1149 {
1150 cerr << "Cannot open output stream " << filename << ": " << strerror(errno) << endl;
1151 return 5;
1152 }
1153 fout.rdbuf(sout.rdbuf());
1154 }
1155
1156 fout.precision(conf.Get<int>("precision"));
1157 if (conf.Get<bool>("fixed"))
1158 fout << fixed;
1159 if (conf.Get<bool>("scientific"))
1160 fout << scientific;
1161
1162 const string filter = conf.Has("filter") ? conf.Get<string>("filter") : "";
1163 const size_t first = conf.Get<size_t>("first");
1164 const size_t limit = conf.Get<size_t>("limit");
1165
1166#ifdef HAVE_ROOT
1167 if (conf.Get<bool>("root"))
1168 {
1169 DumpRoot(fout, conf.Vec<string>("col"), filter, first, limit, filename);
1170 return 0;
1171 }
1172#endif
1173
1174 const vector<string> format = conf.Vec<string>("%");
1175 for (auto it=format.begin(); it<format.end(); it++)
1176 {
1177 static const boost::regex expr("-?[0-9]*[.]?[0-9]*[diouxXeEfFgGaAh]");
1178
1179 boost::smatch what;
1180 if (!boost::regex_match(*it, what, expr, boost::match_extra))
1181 {
1182 cerr << "Format '" << *it << "' not supported." << endl;
1183 return 6;
1184 }
1185 }
1186
1187 const vector<MyColumn> cols = InitColumns(conf.Vec<string>("col"));
1188 if (cols.empty())
1189 return 7;
1190
1191 if (conf.Get<bool>("minmax"))
1192 {
1193 DumpMinMax(fout, cols, first, limit, conf.Get<bool>("nozero"));
1194 return 0;
1195 }
1196
1197 if (conf.Get<bool>("stat"))
1198 {
1199 DumpStats(fout, cols, filter, first, limit);
1200 return 0;
1201 }
1202
1203 return Dump(fout, format, cols, filter, first, limit, filename, conf);
1204}
1205
1206void PrintUsage()
1207{
1208 cout <<
1209 "fitsdump is a tool to dump data from a FITS table as ascii.\n"
1210 "\n"
1211 "Usage: fitsdump [OPTIONS] fitsfile col col ... \n"
1212 " or: fitsdump [OPTIONS]\n"
1213 "\n"
1214 "Addressing a column:\n"
1215 " ColumnName: Will address all fields of a column\n"
1216 " ColumnName[n]: Will address the n-th field of a column (starts with 0)\n"
1217 " ColumnName[n1:n2]: Will address all fields between n1 and including n2\n"
1218#ifdef HAVE_ROOT
1219 "\n"
1220 "Selecting a column:\n"
1221 " Commandline option: --filter\n"
1222 " Explanation: Such a selection is evaluated using TFormula, hence, every "
1223 "mathematical operation allowed in TFormula is allowed there, too. "
1224 "The reference is the column index as printed in the output stream, "
1225 "starting with 1. The index 0 is reserved for the row number.\n"
1226#endif
1227 ;
1228 cout << endl;
1229}
1230
1231void PrintHelp()
1232{
1233#ifdef HAVE_ROOT
1234 cout <<
1235 "\n\n"
1236 "Examples:\n"
1237 "In --root mode, fitsdump support TFormula's syntax for all columns and the filter "
1238 "You can then refer to a column or a (single) index of the column just by its name "
1239 "If the index is omitted, 0 is assumed. Note that the [x:y] syntax in this mode is "
1240 "not supported\n"
1241 "\n"
1242 " fitsdump Zd --filter=\"[0]>20 && cos([1])*TMath::RadToDeg()<45\"\n"
1243 "\n"
1244 "The columns can also be addressed with their names\n"
1245 "\n"
1246 " fitsdump -r \"(Zd+Err)*TMath::DegToRad()\" --filter=\"[0]<25 && [1]<0.05\"\n"
1247 "\n"
1248 "is identical to\n"
1249 "\n"
1250 " fitsdump -r \"(Zd[0]+Err[0])*TMath::DegToRad()\" --filter=\"[0]<25 && [1]<0.05\"\n"
1251 "\n"
1252 "A special placeholder exists for the row number\n"
1253 "\n"
1254 " fitsdump -r \"#\" --filter=\"#>10 && #<100\"\n"
1255 "\n"
1256 "To format a single column you can do\n"
1257 "\n"
1258 " fitsdump col1 -%.1f col2 -%d\n"
1259 "\n"
1260 "A special format is provided converting to 'human readable format'\n"
1261 "\n"
1262 " fitsdump col1 -%h\n"
1263 "\n"
1264 "To convert a fits table to a root file, use --outfile, for example\n\n"
1265 " fitsdump --outfile outfile.root [--compression N] [--overwrite|--update]\n\n"
1266 "Note that this option can not be combined with the root mode (--root). "
1267 "According to the FITS standard, all columns (except bytes) are written as "
1268 "signed values.\n"
1269 "\n";
1270 cout << endl;
1271#endif
1272}
1273
1274
1275void SetupConfiguration(Configuration& conf)
1276{
1277 po::options_description configs("Fitsdump options");
1278 configs.add_options()
1279 ("filecontent", po_switch(), "List the number of tables in the file, along with their name")
1280 ("header,h", po_switch(), "Dump header of given table")
1281 ("list,l", po_switch(), "List all tables and columns in file")
1282 ("fitsfile", var<string>()
1283#if BOOST_VERSION >= 104200
1284 ->required()
1285#endif
1286 , "Name of FITS file")
1287 ("col,c", vars<string>(), "List of columns to dump\narg is a list of columns, separated by a space.\nAdditionnally, a list of sub-columns can be added\ne.g. Data[3] will dump sub-column 3 of column Data\nData[3:4] will dump sub-columns 3 and 4\nOmitting this argument dump the entire column\nnote: all indices start at zero")
1288 ("outfile,o", var<string>("-"), "Name of output file (-:/dev/stdout), if outfile starts with a dot, the extension .fits and .fits.fz are replaced automatically.")
1289 ("precision,p", var<int>(20), "Precision of ofstream")
1290 ("stat,s", po_switch(), "Perform statistics instead of dump")
1291 ("minmax,m", po_switch(), "Calculates min and max of data")
1292 ("nozero,z", po_switch(), "skip 0 values for stats")
1293 ("fixed", po_switch(), "Switch output stream to floating point values in fixed-point notation")
1294 ("scientific", po_switch(), "Switch output stream to floating point values in scientific notation")
1295 ("%,%", vars<string>(), "Format for the output (currently not available in root-mode)")
1296 ("force", po_switch(), "Force reading the fits file even if END key is missing")
1297 ("stream", po_switch(), "Stream reading the fits file (ignore the row count in the header)")
1298 ("first", var<size_t>(size_t(0)), "First number of row to read")
1299 ("limit", var<size_t>(size_t(0)), "Limit for the maximum number of rows to read (0=unlimited)")
1300 ("tablename,t", var<string>(""), "Name of the table to open. If not specified, first binary table is opened")
1301 ("autocal", po_switch(), "This option can be used to skip the application of the noise calibration when reading back a compressed fits file. This is identical in using zfits instead of factfits.")
1302#ifdef HAVE_ROOT
1303 ("root,r", po_switch(), "Enable root mode")
1304 ("filter,f", var<string>(""), "Filter to restrict the selection of events (e.g. '[0]>10 && [0]<20'; does not work with stat and minmax yet)")
1305 ("overwrite", po_switch(), "Force overwriting an existing root file ('RECREATE')")
1306 ("update", po_switch(), "Update an existing root file ('UPDATE') or an existing tree with more data.")
1307 ("compression", var<uint16_t>(1), "zlib compression level for the root file")
1308#endif
1309 ;
1310
1311 po::positional_options_description p;
1312 p.add("fitsfile", 1); // The first positional options
1313 p.add("col", -1); // All others
1314
1315 conf.AddOptions(configs);
1316 conf.SetArgumentPositions(p);
1317}
1318
1319int main(int argc, const char** argv)
1320{
1321 Configuration conf(argv[0]);
1322 conf.SetPrintUsage(PrintUsage);
1323 SetupConfiguration(conf);
1324
1325 if (!conf.DoParse(argc, argv, PrintHelp))
1326 return -1;
1327
1328 if (!conf.Has("fitsfile"))
1329 {
1330 cerr << "Filename required." << endl;
1331 return -1;
1332 }
1333
1334 FitsDumper loader(conf.Get<string>("fitsfile"), conf.Get<string>("tablename"));
1335 if (!loader)
1336 {
1337 cerr << "ERROR - Opening " << conf.Get<string>("fitsfile");
1338 cerr << " failed: " << strerror(errno) << endl;
1339 return -1;
1340 }
1341
1342 return loader.Exec(conf);
1343}
Note: See TracBrowser for help on using the repository browser.