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

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