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

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