Logo ROOT  
Reference Guide
TFITS.cxx
Go to the documentation of this file.
1// @(#)root/graf2d:$Id$
2// Author: Claudi Martinez, July 19th 2010
3
4/*************************************************************************
5 * Copyright (C) 1995-2010, Rene Brun and Fons Rademakers. *
6 * All rights reserved. *
7 * *
8 * For the licensing terms see $ROOTSYS/LICENSE. *
9 * For the list of contributors see $ROOTSYS/README/CREDITS. *
10 *************************************************************************/
11
12/// \defgroup fitsio FITS file
13/// \brief Interface to FITS file.
14/// \ingroup Graphics2D
15///
16/// TFITS is an interface that lets you reading Flexible Image Transport System
17/// (FITS) files, which are generally used in astronomy. This file format
18/// was standardized 1981 and today is still widely used among professional
19/// and amateur astronomers. FITS is not only an image file, but also
20/// it can contain spectrums, data tables, histograms, and multidimensional
21/// data. Furthermore, FITS data can be described itself by containing
22/// human-readable information that let us to interpret the data within
23/// the FITS file. For example, a FITS could contain a 3D data cube,
24/// but an additional description would tell us that we must read it, for
25/// example, as a 3-layer image.
26///
27/// TFITS requires CFITSIO library to be installed on your system. It
28/// is currently maintained by NASA/GSFC and can be downloaded from
29/// [NASA/GSFC web site](http://fits.gsfc.nasa.gov), as well as documentation.
30///
31/// Using this interface is easy and straightforward. There is only 1 class
32/// called "TFITSHDU" which has several methods to extract data from a
33/// FITS file, more specifically, from an HDU within the file. An HDU, or
34/// Header Data Unit, is a chunk of data with a header containing several
35/// "keyword = value" tokens. The header describes the structure of data
36/// within the HDU. An HDU can be of two types: an "image HDU" or a "table
37/// HDU". The former can be any kind of multidimensional array of real numbers,
38/// by which the name "image" may be confusing: you can store an image, but
39/// you can also store a N-dimensional data cube. On the other hand, table
40/// HDUs are sets of several rows and columns (a.k.a fields) which contain
41/// generic data, as strings, real or complex numbers and even arrays.
42///
43/// Please have a look to the tutorials ($ROOTSYS/tutorials/fitsio/) to see
44/// some examples. IMPORTANT: to run tutorials it is required that
45/// you change the current working directory of ROOT (CINT) shell to the
46/// tutorials directory. Example:
47/// ~~~ {.cpp}
48/// root [1] gSystem->ChangeDirectory("tutorials/fitsio")
49/// root [1] .x FITS_tutorial1.C
50/// ~~~
51/// LIST OF TODO
52/// - Support for complex values within data tables
53/// - Support for reading arrays from table cells
54/// - Support for grouping
55///
56/// IMPLEMENTATION NOTES:
57///
58/// CFITSIO library uses standard C types ('int', 'long', ...). To avoid
59/// confusion, the same types are used internally by the access methods.
60/// However, class's fields are ROOT-defined types.
61
62/** \class TFITSHDU
63\ingroup fitsio
64
65FITS file interface class
66
67TFITSHDU is a class that allows extracting images and data from FITS files and contains
68several methods to manage them.
69*/
70
71#include "TFITS.h"
72#include "TROOT.h"
73#include "TImage.h"
74#include "TArrayI.h"
75#include "TArrayD.h"
76#include "TH1D.h"
77#include "TH2D.h"
78#include "TH3D.h"
79#include "TVectorD.h"
80#include "TMatrixD.h"
81#include "TObjArray.h"
82#include "TObjString.h"
83#include "TCanvas.h"
84#include "TMath.h"
85
86#include "fitsio.h"
87#include <stdlib.h>
88
90
91////////////////////////////////////////////////////////////////////////////////
92/// Clean path from possible filter and put the result in 'dst'.
93
94void TFITSHDU::CleanFilePath(const char *filepath_with_filter, TString &dst)
95{
96 dst = filepath_with_filter;
97
98 Ssiz_t ndx = dst.Index("[", 1, 0, TString::kExact);
99 if (ndx != kNPOS) {
100 dst.Resize(ndx);
101 }
102}
103
104
105////////////////////////////////////////////////////////////////////////////////
106/// TFITSHDU constructor from file path with HDU selection filter.
107/// Please refer to CFITSIO manual for more information about
108/// HDU selection filters.
109///
110/// Examples:
111/// - `TFITSHDU("/path/to/myfile.fits")`: just open the PRIMARY HDU
112/// - `TFITSHDU("/path/to/myfile.fits[1]")`: open HDU #1
113/// - `TFITSHDU("/path/to/myfile.fits[PICS]")`: open HDU called 'PICS'
114/// - `TFITSHDU("/path/to/myfile.fits[ACQ][EXPOSURE > 5]")`: open the (table) HDU called 'ACQ' and
115/// selects the rows that have column 'EXPOSURE'
116/// greater than 5.
117
118TFITSHDU::TFITSHDU(const char *filepath_with_filter)
119{
121
122 fFilePath = filepath_with_filter;
123 CleanFilePath(filepath_with_filter, fBaseFilePath);
124
125 if (kFALSE == LoadHDU(fFilePath)) {
127 throw -1;
128 }
129}
130
131////////////////////////////////////////////////////////////////////////////////
132/// TFITSHDU constructor from filepath and extension number.
133
134TFITSHDU::TFITSHDU(const char *filepath, Int_t extension_number)
135{
137 CleanFilePath(filepath, fBaseFilePath);
138
139 //Add "by extension number" filter
140 fFilePath.Form("%s[%d]", fBaseFilePath.Data(), extension_number);
141
142 if (kFALSE == LoadHDU(fFilePath)) {
144 throw -1;
145 }
146}
147
148////////////////////////////////////////////////////////////////////////////////
149/// TFITSHDU constructor from filepath and extension name.
150
151TFITSHDU::TFITSHDU(const char *filepath, const char *extension_name)
152{
154 CleanFilePath(filepath, fBaseFilePath);
155
156 //Add "by extension number" filter
157 fFilePath.Form("%s[%s]", fBaseFilePath.Data(), extension_name);
158
159
160 if (kFALSE == LoadHDU(fFilePath)) {
162 throw -1;
163 }
164}
165
166////////////////////////////////////////////////////////////////////////////////
167/// TFITSHDU destructor.
168
170{
172}
173
174////////////////////////////////////////////////////////////////////////////////
175/// Release internal resources.
176
178{
179 if (fRecords) delete [] fRecords;
180
181 if (fType == kImageHDU) {
182 if (fSizes) delete fSizes;
183 if (fPixels) delete fPixels;
184 } else {
185 if (fColumnsInfo) {
186 if (fCells) {
187 for (Int_t i = 0; i < fNColumns; i++) {
188 if (fColumnsInfo[i].fType == kString) {
189 //Deallocate character arrays allocated for kString columns
190 Int_t offset = i * fNRows;
191 for (Int_t row = 0; row < fNRows; row++) {
192 delete [] fCells[offset+row].fString;
193 }
194 } else if (fColumnsInfo[i].fType == kRealArray) {
195 //Deallocate character arrays allocated for kString columns
196 Int_t offset = i * fNRows;
197 for (Int_t row = 0; row < fNRows; row++) {
198 delete [] fCells[offset+row].fRealArray;
199 }
200 } else if (fColumnsInfo[i].fType == kRealVector) {
201 // Deallocate character arrays allocated for variable-length array columns
202 Int_t offset = i * fNRows;
203 for (Int_t row = 0; row < fNRows; row++) {
204 delete fCells[offset + row].fRealVector;
205 }
206 }
207 }
208
209 delete [] fCells;
210 }
211
212 delete [] fColumnsInfo;
213 }
214
215
216 }
217}
218
219////////////////////////////////////////////////////////////////////////////////
220/// Do some initializations.
221
223{
224 fRecords = 0;
225 fPixels = 0;
226 fSizes = 0;
227 fColumnsInfo = 0;
228 fNColumns = fNRows = 0;
229 fCells = 0;
230}
231
232////////////////////////////////////////////////////////////////////////////////
233/// Load HDU from fits file satisfying the specified filter.
234/// Returns kTRUE if success. Otherwise kFALSE.
235/// If filter == "" then the primary array is selected
236
238{
239 fitsfile *fp=0;
240 int status = 0;
241 char errdescr[FLEN_STATUS+1];
242
243 // Open file with filter
244 fits_open_file(&fp, filepath_filter.Data(), READONLY, &status);
245 if (status) goto ERR;
246
247 // Read HDU number
248 int hdunum;
249 fits_get_hdu_num(fp, &hdunum);
250 fNumber = Int_t(hdunum);
251
252 // Read HDU type
253 int hdutype;
254 fits_get_hdu_type(fp, &hdutype, &status);
255 if (status) goto ERR;
256 fType = (hdutype == IMAGE_HDU) ? kImageHDU : kTableHDU;
257
258 //Read HDU header records
259 int nkeys, morekeys;
260 char keyname[FLEN_KEYWORD+1];
261 char keyvalue[FLEN_VALUE+1];
262 char comment[FLEN_COMMENT+1];
263
264 fits_get_hdrspace(fp, &nkeys, &morekeys, &status);
265 if (status) goto ERR;
266
267 fRecords = new struct HDURecord[nkeys];
268
269 for (int i = 1; i <= nkeys; i++) {
270 fits_read_keyn(fp, i, keyname, keyvalue, comment, &status);
271 if (status) goto ERR;
272 fRecords[i-1].fKeyword = keyname;
273 fRecords[i-1].fValue = keyvalue;
275 }
276
277 fNRecords = Int_t(nkeys);
278
279 //Set extension name
280 fExtensionName = "PRIMARY"; //Default
281 for (int i = 0; i < nkeys; i++) {
282 if (fRecords[i].fKeyword == "EXTNAME") {
284 break;
285 }
286 }
287
288 //Read HDU's data
289 if (fType == kImageHDU) {
290 Info("LoadHDU", "The selected HDU contains an Image Extension");
291
292 int param_ndims=0;
293 long *param_dimsizes;
294
295 //Read image number of dimensions
296 fits_get_img_dim(fp, &param_ndims, &status);
297 if (status) goto ERR;
298 if (param_ndims > 0) {
299 //Read image sizes in each dimension
300 param_dimsizes = new long[param_ndims];
301 fits_get_img_size(fp, param_ndims, param_dimsizes, &status);
302 if (status) {
303 delete [] param_dimsizes;
304 goto ERR;
305 }
306
307 fSizes = new TArrayI(param_ndims);
308 fSizes = new TArrayI(param_ndims);
309 for (int i = 0; i < param_ndims; i++) { //Use for loop to copy values instead of passing array to constructor, since 'Int_t' size may differ from 'long' size
310 fSizes->SetAt(param_dimsizes[i], i);
311 }
312
313 delete [] param_dimsizes;
314
315 //Read pixels
316 int anynul;
317 long *firstpixel = new long[param_ndims];
318 double nulval = 0;
319 long npixels = 1;
320
321 for (int i = 0; i < param_ndims; i++) {
322 npixels *= (long) fSizes->GetAt(i); //Compute total number of pixels
323 firstpixel[i] = 1; //Set first pixel to read from.
324 }
325
326 double *pixels = new double[npixels];
327
328 fits_read_pix(fp, TDOUBLE, firstpixel, npixels,
329 (void *) &nulval, (void *) pixels, &anynul, &status);
330
331 if (status) {
332 delete [] firstpixel;
333 delete [] pixels;
334 goto ERR;
335 }
336
337 fPixels = new TArrayD(npixels, pixels);
338
339 delete [] firstpixel;
340 delete [] pixels;
341
342 } else {
343 //Null array
344 fSizes = new TArrayI();
345 fPixels = new TArrayD();
346 }
347 } else {
348 Info("LoadHDU", "The selected HDU contains a Table Extension");
349
350 // Get table's number of rows and columns
351 long table_rows;
352 int table_cols;
353
354 fits_get_num_rows(fp, &table_rows, &status);
355 if (status) goto ERR;
356
357 fNRows = Int_t(table_rows);
358
359 fits_get_num_cols(fp, &table_cols, &status);
360 if (status) goto ERR;
361
362 fNColumns = Int_t(table_cols);
363
364 // Allocate column info array
365 fColumnsInfo = new struct Column[table_cols];
366
367 // Read column names
368 char colname[80];
369 int colnum;
370
371 fits_get_colname(fp, CASEINSEN, (char*) "*", colname, &colnum, &status);
372 while (status == COL_NOT_UNIQUE)
373 {
374 fColumnsInfo[colnum-1].fName = colname;
375 fits_get_colname(fp, CASEINSEN, (char*) "*", colname, &colnum, &status);
376 }
377 if (status != COL_NOT_FOUND) goto ERR;
378 status = 0;
379
380 //Allocate cells
381 fCells = new union Cell [table_rows * table_cols];
382
383 // Read columns
384 int typecode;
385 long repeat, width;
386 Int_t cellindex;
387
388
389 for (colnum = 0, cellindex = 0; colnum < fNColumns; colnum++) {
390 fits_get_coltype(fp, colnum+1, &typecode, &repeat, &width, &status);
391
392 if (status) goto ERR;
393
394 if ((typecode == TDOUBLE) || (typecode == TSHORT) || (typecode == TLONG)
395 || (typecode == TFLOAT) || (typecode == TLOGICAL) || (typecode == TBIT)
396 || (typecode == TBYTE) || (typecode == TSTRING)) {
397
398 if (typecode == TSTRING) {
399 // this column contains strings
400 fColumnsInfo[colnum].fType = kString;
401
402 int dispwidth=0;
403 fits_get_col_display_width(fp, colnum+1, &dispwidth, &status);
404 if (status) goto ERR;
405
406
407 char *nulval = (char*) "";
408 int anynul=0;
409 char **array;
410
411 if (dispwidth <= 0) {
412 dispwidth = 1;
413 }
414
415 array = new char* [table_rows];
416 for (long row = 0; row < table_rows; row++) {
417 array[row] = new char[dispwidth+1]; //also room for end null!
418 }
419
420 if (repeat > 0) {
421 fits_read_col(fp, TSTRING, colnum+1, 1, 1, table_rows, nulval, array, &anynul, &status);
422 if (status) {
423 for (long row = 0; row < table_rows; row++) {
424 delete [] array[row];
425 }
426 delete [] array;
427 goto ERR;
428 }
429
430 } else {
431 //No elements: set dummy
432 for (long row = 0; row < table_rows; row++) {
433 strlcpy(array[row], "-",dispwidth+1);
434 }
435 }
436
437 //Save values
438 for (long row = 0; row < table_rows; row++) {
439 fCells[cellindex++].fString = array[row];
440 }
441
442 delete [] array; //Delete temporal array holding pointer to strings, but not delete strings themselves!
443
444
445 } else {
446 // this column contains either a number or a fixed-length array
447 double nulval = 0;
448 int anynul=0;
449
450 fColumnsInfo[colnum].fDim = (Int_t) repeat;
451
452 double *array = 0;
453 char *arrayl = 0;
454
455 if (repeat > 0) {
456
457 if (typecode == TLOGICAL) {
458 arrayl = new char[table_rows * repeat];
459 fits_read_col(fp, TLOGICAL, colnum + 1, 1, 1, table_rows * repeat, &nulval, arrayl, &anynul,
460 &status);
461 if (status) {
462 delete[] arrayl;
463 goto ERR;
464 }
465 } else {
466 array = new double[table_rows * repeat]; // Hope you got a big machine! Ask China otherwise :-)
467 fits_read_col(fp, TDOUBLE, colnum + 1, 1, 1, table_rows * repeat, &nulval, array, &anynul,
468 &status);
469 if (status) {
470 delete[] array;
471 goto ERR;
472 }
473 }
474
475 } else {
476 // No elements: set dummy
477 array = new double[table_rows];
478 for (long row = 0; row < table_rows; row++) {
479 array[row] = 0.0;
480 }
481 }
482
483 // Save values
484 if (repeat == 1) {
485 // this column contains scalars
487 if (typecode == TLOGICAL) {
488 for (long row = 0; row < table_rows; row++) {
489 int temp = (signed char)arrayl[row];
490 fCells[cellindex++].fRealNumber = (double)temp;
491 }
492 delete[] arrayl;
493 } else {
494 for (long row = 0; row < table_rows; row++) {
495 fCells[cellindex++].fRealNumber = array[row];
496 }
497 delete[] array;
498 }
499 } else if (repeat > 1) {
500 // this column contains fixed-length arrays
501 fColumnsInfo[colnum].fType = kRealArray;
502 if (typecode == TLOGICAL) {
503 for (long row = 0; row < table_rows; row++) {
504 double *vec = new double[repeat];
505 long offset = row * repeat;
506 for (long component = 0; component < repeat; component++) {
507 int temp = (signed char)arrayl[offset++];
508 vec[component] = (double)temp;
509 }
510 fCells[cellindex++].fRealArray = vec;
511 }
512 delete[] arrayl;
513 } else {
514 for (long row = 0; row < table_rows; row++) {
515 double *vec = new double[repeat];
516 long offset = row * repeat;
517 for (long component = 0; component < repeat; component++) {
518 vec[component] = array[offset++];
519 }
520 fCells[cellindex++].fRealArray = vec;
521 }
522 delete[] array;
523 }
524 }
525 }
526 } else if (typecode < 1) {
527 // this column contains variable-length arrays
529
530 // null variables needed by the fits_read_col
531 double nulval = 0;
532 int anynul=0;
533
534 // loop over rows to derive the total memory requirement
535 fColumnsInfo[colnum].fRowStart.assign(table_rows+1, 0);
536 fColumnsInfo[colnum].fVarLengths.assign(table_rows, 0);
537 fColumnsInfo[colnum].fRowStart[0] = 0;
538
539 for (long row = 0; row < table_rows; row++) {
540 long offset = 0;
541 long repeat_row = 0;
542
543 fits_read_descript(fp, colnum+1, row+1, &repeat_row, &offset, &status);
544 // store the starting of each row and the number of elements it contains
545 fColumnsInfo[colnum].fRowStart[row+1] = fColumnsInfo[colnum].fRowStart[row] + repeat_row;
546 fColumnsInfo[colnum].fVarLengths[row] = repeat_row;
547 }
548
549 for (long row = 0; row < table_rows; row++) {
550 // number of elements in the cell we want to read, i.e.
551 // number of elements in the variable-length array
552 int nelements = fColumnsInfo[colnum].fRowStart[row + 1] - fColumnsInfo[colnum].fRowStart[row];
553 // size of the variable-length array
554 const int size = fColumnsInfo[colnum].fVarLengths[row];
555 // vector to store the results
556 TArrayD *vec = new TArrayD(size);
557 // variable-length array have negative DATATYPE
558 int abstype = TMath::Abs(typecode);
559
560 // define the array to load the data with the CFITSIO functions
561 // a fixed arrays is needed as argument to the fits_read_col function
562 // so a vector is defined and then its `.data()` pointer is passed
563 // to fits_read_col
564 //
565 // TODO: add all type cases, for now only short, long, float and double are considered
566 //
567 if (abstype == 21) {
568 std::vector<short> data;
569 data.resize(size);
570 fits_read_col(fp, abstype, colnum + 1, row + 1, 1, nelements, &nulval, data.data(), &anynul, &status);
571 for (int i = 0; i < size; i++)
572 vec->SetAt(data[i], i);
573 } else if (abstype == 41) {
574 std::vector<int> data;
575 data.resize(size);
576 fits_read_col(fp, abstype, colnum + 1, row + 1, 1, nelements, &nulval, data.data(), &anynul, &status);
577 for (int i = 0; i < size; i++)
578 vec->SetAt(data[i], i);
579 } else if (abstype == 42) {
580 std::vector<float> data;
581 data.resize(size);
582 fits_read_col(fp, abstype, colnum + 1, row + 1, 1, nelements, &nulval, data.data(), &anynul, &status);
583 for (int i = 0; i < size; i++)
584 vec->SetAt(data[i], i);
585 } else if (abstype == 82) {
586 std::vector<double> data;
587 data.resize(size);
588 fits_read_col(fp, abstype, colnum + 1, row + 1, 1, nelements, &nulval, data.data(), &anynul, &status);
589 for (int i = 0; i < size; i++)
590 vec->SetAt(data[i], i);
591 } else {
592 Error("LoadHDU", "The variable-length array type in column %d is unknown", colnum + 1);
593 }
594 // place the vector storing the variable-length array in the corresponding cell
595 fCells[cellindex++].fRealVector = vec;
596 }
597 } else {
598 Warning("LoadHDU", "error opening FITS file. Column type %d is currently not supported", typecode);
599 }
600 }
601
602 if (hdutype == ASCII_TBL) {
603 // ASCII table
604
605 } else {
606 // Binary table
607 }
608 }
609
610 // Close file
611 fits_close_file(fp, &status);
612 return kTRUE;
613
614ERR:
615 fits_get_errstatus(status, errdescr);
616 Warning("LoadHDU", "error opening FITS file. Details: %s", errdescr);
617 status = 0;
618 if (fp) fits_close_file(fp, &status);
619 return kFALSE;
620}
621
622////////////////////////////////////////////////////////////////////////////////
623/// Get record by keyword.
624
625struct TFITSHDU::HDURecord* TFITSHDU::GetRecord(const char *keyword)
626{
627 for (int i = 0; i < fNRecords; i++) {
628 if (fRecords[i].fKeyword == keyword) {
629 return &fRecords[i];
630 }
631 }
632 return 0;
633}
634
635////////////////////////////////////////////////////////////////////////////////
636/// Get the value of a given keyword. Return "" if not found.
637
639{
640 HDURecord *rec = GetRecord(keyword);
641 if (rec) {
642 return rec->fValue;
643 } else {
644 return *(new TString(""));
645 }
646}
647
648////////////////////////////////////////////////////////////////////////////////
649/// Print records.
650
652{
653 for (int i = 0; i < fNRecords; i++) {
654 if (fRecords[i].fComment.Length() > 0) {
655 printf("%-10s = %s / %s\n", fRecords[i].fKeyword.Data(), fRecords[i].fValue.Data(), fRecords[i].fComment.Data());
656 } else {
657 printf("%-10s = %s\n", fRecords[i].fKeyword.Data(), fRecords[i].fValue.Data());
658 }
659 }
660}
661
662////////////////////////////////////////////////////////////////////////////////
663/// Print HDU's parent file's metadata.
664
666{
667 fitsfile *fp=0;
668 int status = 0;
669 char errdescr[FLEN_STATUS+1];
670 int hducount, extnum;
671 int hdutype = IMAGE_HDU;
672 const char *exttype;
673 char extname[FLEN_CARD]="PRIMARY"; //room enough
674 int verbose = (opt[0] ? 1 : 0);
675
676 // Open file with no filters: current HDU will be the primary one.
677 fits_open_file(&fp, fBaseFilePath.Data(), READONLY, &status);
678 if (status) goto ERR;
679
680 // Read HDU count
681 fits_get_num_hdus(fp, &hducount, &status);
682 if (status) goto ERR;
683 printf("Total: %d HDUs\n", hducount);
684
685 extnum = 0;
686 while(hducount) {
687 // Read HDU type
688 fits_get_hdu_type(fp, &hdutype, &status);
689 if (status) goto ERR;
690
691 if (hdutype == IMAGE_HDU) {
692 exttype="IMAGE";
693 } else if (hdutype == ASCII_TBL) {
694 exttype="ASCII TABLE";
695 } else {
696 exttype="BINARY TABLE";
697 }
698
699 //Read HDU header records
700 int nkeys, morekeys;
701 char keyname[FLEN_KEYWORD+1];
702 char keyvalue[FLEN_VALUE+1];
703 char comment[FLEN_COMMENT+1];
704
705 fits_get_hdrspace(fp, &nkeys, &morekeys, &status);
706 if (status) goto ERR;
707
708 struct HDURecord *records = new struct HDURecord[nkeys];
709
710 for (int i = 1; i <= nkeys; i++) {
711 fits_read_keyn(fp, i, keyname, keyvalue, comment, &status);
712 if (status) {
713 delete [] records;
714 goto ERR;
715 }
716
717 records[i-1].fKeyword = keyname;
718 records[i-1].fValue = keyvalue;
719 records[i-1].fComment = comment;
720
721 if (strcmp(keyname, "EXTNAME") == 0) {
722 //Extension name
723 strlcpy(extname, keyvalue,FLEN_CARD);
724 }
725 }
726
727 //HDU info
728 printf(" [%d] %s (%s)\n", extnum, exttype, extname);
729
730 //HDU records
731 if (verbose) {
732 for (int i = 0; i < nkeys; i++) {
733 if (comment[0]) {
734 printf(" %-10s = %s / %s\n", records[i].fKeyword.Data(), records[i].fValue.Data(), records[i].fComment.Data());
735 } else {
736 printf(" %-10s = %s\n", records[i].fKeyword.Data(), records[i].fValue.Data());
737 }
738 }
739 }
740 printf("\n");
741
742 delete [] records;
743
744 hducount--;
745 extnum++;
746 if (hducount){
747 fits_movrel_hdu(fp, 1, &hdutype, &status);
748 if (status) goto ERR;
749 }
750 }
751
752 // Close file
753 fits_close_file(fp, &status);
754 return;
755
756ERR:
757 fits_get_errstatus(status, errdescr);
758 Warning("PrintFileMetadata", "error opening FITS file. Details: %s", errdescr);
759 status = 0;
760 if (fp) fits_close_file(fp, &status);
761}
762
763////////////////////////////////////////////////////////////////////////////////
764/// Print column information
765
767{
768 if (fType != kTableHDU) {
769 Warning("PrintColumnInfo", "this is not a table HDU.");
770 return;
771 }
772
773 for (Int_t i = 0; i < fNColumns; i++) {
774 switch (fColumnsInfo[i].fType){
775 case kString:
776 printf("%-20s : %s\n", fColumnsInfo[i].fName.Data(), "STRING");
777 break;
778 case kRealNumber:
779 printf("%-20s : %s\n", fColumnsInfo[i].fName.Data(), "REAL NUMBER");
780 break;
781 case kRealArray:
782 printf("%-20s : %s\n", fColumnsInfo[i].fName.Data(), "FIXED-LENGTH ARRAY");
783 break;
784 case kRealVector:
785 printf("%-20s : %s\n", fColumnsInfo[i].fName.Data(), "VARIABLE-LENGTH ARRAY");
786 break;
787 }
788 }
789}
790
791////////////////////////////////////////////////////////////////////////////////
792/// Print full table contents
793
795{
796 int printed_chars;
797
798 if (fType != kTableHDU) {
799 Warning("PrintColumnInfo", "this is not a table HDU.");
800 return;
801 }
802
803 // check that the table does not contain fixed or variable length arrays
804 // in that case is not possible to print a flattened table
805 for (Int_t col = 0; col < fNColumns; col++){
806 if (fColumnsInfo[col].fType == kRealArray) {
807 Warning("PrintColumnInfo", "The table contains column with fixed-length arrays and cannot be flattened for printing.");
808 return;
809 }
810 else if (fColumnsInfo[col].fType == kRealVector) {
811 Warning("PrintColumnInfo", "The table contains column with variable-length arrays and cannot be flattened for printing.");
812 return;
813 }
814 }
815
816 // Dump header
817 putchar('\n');
818 printed_chars = 0;
819 for (Int_t col = 0; col < fNColumns; col++) {
820 printed_chars += printf("%-10s| ", fColumnsInfo[col].fName.Data());
821 }
822 putchar('\n');
823 while(printed_chars--) {
824 putchar('-');
825 }
826 putchar('\n');
827
828 // Dump rows
829 for (Int_t row = 0; row < fNRows; row++) {
830 for (Int_t col = 0; col < fNColumns; col++) {
831 if (fColumnsInfo[col].fType == kString) {
832 printf("%-10s", fCells[col * fNRows + row].fString);
833 } else if (fColumnsInfo[col].fType == kRealNumber) {
834 printed_chars = printf("%.2lg", fCells[col * fNRows + row].fRealNumber);
835 printed_chars -= 10;
836 while (printed_chars < 0) {
837 putchar(' ');
838 printed_chars++;
839 }
840 }
841
842 if (col <= fNColumns - 1) printf("| ");
843 }
844 printf("\n");
845 }
846}
847
848////////////////////////////////////////////////////////////////////////////////
849/// Print metadata.
850/// Currently supported options:
851///
852/// - "" : print HDU record data
853/// - "F" : print FITS file's extension names, numbers and types
854/// - "F+": print FITS file's extension names and types and their record data
855/// - "T" : print column information when HDU is a table
856/// - "T+" : print full table (columns header and rows)
857
858void TFITSHDU::Print(const Option_t *opt) const
859{
860 if ((opt[0] == 'F') || (opt[0] == 'f')) {
861 PrintFileMetadata((opt[1] == '+') ? "+" : "");
862 } else if ((opt[0] == 'T') || (opt[0] == 't')) {
863 if (opt[1] == '+') {
864 PrintFullTable("");
865 } else {
866 PrintColumnInfo("");
867 }
868
869 } else {
871 }
872}
873
874
875////////////////////////////////////////////////////////////////////////////////
876/// Read image HDU as a displayable image. Return 0 if conversion cannot be done.
877/// If the HDU seems to be a multilayer image, 'layer' parameter can be used
878/// to retrieve the specified layer (starting from 0)
879
881{
882 if (fType != kImageHDU) {
883 Warning("ReadAsImage", "this is not an image HDU.");
884 return 0;
885 }
886
887 if (((fSizes->GetSize() != 2) && (fSizes->GetSize() != 3) && (fSizes->GetSize() != 4)) || ((fSizes->GetSize() == 4) && (fSizes->GetAt(3) > 1))) {
888 Warning("ReadAsImage", "could not convert image HDU to image because it has %d dimensions.", fSizes->GetSize());
889 return 0;
890 }
891
892 Int_t width, height;
893 UInt_t pixels_per_layer;
894
895 width = Int_t(fSizes->GetAt(0));
896 height = Int_t(fSizes->GetAt(1));
897
898 pixels_per_layer = UInt_t(width) * UInt_t(height);
899
900 if ( ((fSizes->GetSize() == 2) && (layer > 0))
901 || (((fSizes->GetSize() == 3) || (fSizes->GetSize() == 4)) && (layer >= fSizes->GetAt(2)))) {
902
903 Warning("ReadAsImage", "layer out of bounds.");
904 return 0;
905 }
906
907 // Get the maximum and minimum pixel values in the layer to auto-stretch pixels
908 Double_t maxval = 0, minval = 0;
909 UInt_t i;
910 Double_t pixvalue;
911 Int_t offset = layer * pixels_per_layer;
912
913 for (i = 0; i < pixels_per_layer; i++) {
914 pixvalue = fPixels->GetAt(offset + i);
915
916 if (pixvalue > maxval) {
917 maxval = pixvalue;
918 }
919
920 if ((i == 0) || (pixvalue < minval)) {
921 minval = pixvalue;
922 }
923 }
924
925 //Build the image stretching pixels into a range from 0.0 to 255.0
926 //TImage *im = new TImage(width, height);
927 TImage *im = TImage::Create();
928 if (!im) return 0;
929 TArrayD *layer_pixels = new TArrayD(pixels_per_layer);
930
931
932 if (maxval == minval) {
933 //plain image
934 for (i = 0; i < pixels_per_layer; i++) {
935 layer_pixels->SetAt(255.0, i);
936 }
937 } else {
938 Double_t factor = 255.0 / (maxval-minval);
939 for (i = 0; i < pixels_per_layer; i++) {
940 pixvalue = fPixels->GetAt(offset + i);
941 layer_pixels->SetAt(factor * (pixvalue-minval), i) ;
942 }
943 }
944
945 if (pal == 0) {
946 // Default palette: grayscale palette
947 pal = new TImagePalette(256);
948 for (i = 0; i < 256; i++) {
949 pal->fPoints[i] = ((Double_t)i)/255.0;
950 pal->fColorAlpha[i] = 0xffff;
951 pal->fColorBlue[i] = pal->fColorGreen[i] = pal->fColorRed[i] = i << 8;
952 }
953 pal->fPoints[0] = 0;
954 pal->fPoints[255] = 1.0;
955
956 im->SetImage(*layer_pixels, UInt_t(width), pal);
957
958 delete pal;
959
960 } else {
961 im->SetImage(*layer_pixels, UInt_t(width), pal);
962 }
963
964 delete layer_pixels;
965
966 return im;
967}
968
969////////////////////////////////////////////////////////////////////////////////
970/// If the HDU is an image, draw the first layer of the primary array
971/// To set a title to the canvas, pass it in "opt"
972
974{
975 if (fType != kImageHDU) {
976 Warning("Draw", "cannot draw. This is not an image HDU.");
977 return;
978 }
979
980 TImage *im = ReadAsImage(0, 0);
981 if (im) {
983 Int_t height = Int_t(fSizes->GetAt(1));
984 TString cname, ctitle;
985 cname.Form("%sHDU", this->GetName());
986 ctitle.Form("%d x %d", width, height);
987 new TCanvas(cname, ctitle, width, height);
988 im->Draw();
989 }
990}
991
992
993////////////////////////////////////////////////////////////////////////////////
994/// Read image HDU as a matrix. Return 0 if conversion cannot be done
995/// If the HDU seems to be a multilayer image, 'layer' parameter can be used
996/// to retrieve the specified layer (starting from 0) in matrix form.
997/// Options (value of 'opt'):
998/// "S": stretch pixel values to a range from 0.0 to 1.0
999
1001{
1002 if (fType != kImageHDU) {
1003 Warning("ReadAsMatrix", "this is not an image HDU.");
1004 return 0;
1005 }
1006
1007 if (((fSizes->GetSize() != 2) && (fSizes->GetSize() != 3) && (fSizes->GetSize() != 4)) || ((fSizes->GetSize() == 4) && (fSizes->GetAt(3) > 1))) {
1008 Warning("ReadAsMatrix", "could not convert image HDU to image because it has %d dimensions.", fSizes->GetSize());
1009 return 0;
1010 }
1011
1012
1013 if ( ((fSizes->GetSize() == 2) && (layer > 0))
1014 || (((fSizes->GetSize() == 3) || (fSizes->GetSize() == 4)) && (layer >= fSizes->GetAt(2)))) {
1015 Warning("ReadAsMatrix", "layer out of bounds.");
1016 return 0;
1017 }
1018
1019 Int_t width, height;
1020 UInt_t pixels_per_layer;
1021 Int_t offset;
1022 UInt_t i;
1023 TMatrixD *mat=0;
1024
1025 width = Int_t(fSizes->GetAt(0));
1026 height = Int_t(fSizes->GetAt(1));
1027
1028 pixels_per_layer = UInt_t(width) * UInt_t(height);
1029 offset = layer * pixels_per_layer;
1030
1031 double *layer_pixels = new double[pixels_per_layer];
1032
1033 if ((opt[0] == 'S') || (opt[0] == 's')) {
1034 //Stretch
1035 // Get the maximum and minimum pixel values in the layer to auto-stretch pixels
1036 Double_t factor, maxval=0, minval=0;
1037 Double_t pixvalue;
1038 for (i = 0; i < pixels_per_layer; i++) {
1039 pixvalue = fPixels->GetAt(offset + i);
1040
1041 if (pixvalue > maxval) {
1042 maxval = pixvalue;
1043 }
1044
1045 if ((i == 0) || (pixvalue < minval)) {
1046 minval = pixvalue;
1047 }
1048 }
1049
1050 if (maxval == minval) {
1051 //plain image
1052 for (i = 0; i < pixels_per_layer; i++) {
1053 layer_pixels[i] = 1.0;
1054 }
1055 } else {
1056 factor = 1.0 / (maxval-minval);
1057 mat = new TMatrixD(height, width);
1058
1059 for (i = 0; i < pixels_per_layer; i++) {
1060 layer_pixels[i] = factor * (fPixels->GetAt(offset + i) - minval);
1061 }
1062 }
1063
1064 } else {
1065 //No stretching
1066 mat = new TMatrixD(height, width);
1067
1068 for (i = 0; i < pixels_per_layer; i++) {
1069 layer_pixels[i] = fPixels->GetAt(offset + i);
1070 }
1071 }
1072
1073 if (mat) {
1074 // mat->Use(height, width, layer_pixels);
1075 memcpy(mat->GetMatrixArray(), layer_pixels, pixels_per_layer*sizeof(double));
1076 }
1077
1078 delete [] layer_pixels;
1079 return mat;
1080}
1081
1082
1083////////////////////////////////////////////////////////////////////////////////
1084/// Read image HDU as a histogram. Return 0 if conversion cannot be done.
1085/// The returned object can be TH1D, TH2D or TH3D depending on data dimensionality.
1086/// Please, check condition (returnedValue->IsA() == TH*D::Class()) to
1087/// determine the object class.
1088///
1089/// NOTE: do not confuse with image histogram! This function interprets
1090/// the array as a histogram. It does not compute the histogram of pixel
1091/// values of an image! Here "pixels" are interpreted as number of entries.
1092
1094{
1095 if (fType != kImageHDU) {
1096 Warning("ReadAsHistogram", "this is not an image HDU.");
1097 return 0;
1098 }
1099
1100 TH1 *result=0;
1101
1102 if ((fSizes->GetSize() != 1) && (fSizes->GetSize() != 2) && (fSizes->GetSize() != 3)) {
1103 Warning("ReadAsHistogram", "could not convert image HDU to histogram because it has %d dimensions.", fSizes->GetSize());
1104 return 0;
1105 }
1106
1107 if (fSizes->GetSize() == 1) {
1108 //1D
1109 UInt_t Nx = UInt_t(fSizes->GetAt(0));
1110 UInt_t x;
1111
1112 TH1D *h = new TH1D("", "", Int_t(Nx), 0, Nx-1);
1113
1114 for (x = 0; x < Nx; x++) {
1116 if (nentries < 0) nentries = 0; //Crop negative values
1117 h->Fill(x, nentries);
1118 }
1119
1120 result = h;
1121
1122 } else if (fSizes->GetSize() == 2) {
1123 //2D
1124 UInt_t Nx = UInt_t(fSizes->GetAt(0));
1125 UInt_t Ny = UInt_t(fSizes->GetAt(1));
1126 UInt_t x,y;
1127
1128 TH2D *h = new TH2D("", "", Int_t(Nx), 0, Nx-1, Int_t(Ny), 0, Ny-1);
1129
1130 for (y = 0; y < Ny; y++) {
1131 UInt_t offset = y * Nx;
1132 for (x = 0; x < Nx; x++) {
1133 Int_t nentries = Int_t(fPixels->GetAt(offset + x));
1134 if (nentries < 0) nentries = 0; //Crop negative values
1135 h->Fill(x,y, nentries);
1136 }
1137 }
1138
1139 result = h;
1140
1141 } else if (fSizes->GetSize() == 3) {
1142 //3D
1143 UInt_t Nx = UInt_t(fSizes->GetAt(0));
1144 UInt_t Ny = UInt_t(fSizes->GetAt(1));
1145 UInt_t Nz = UInt_t(fSizes->GetAt(2));
1146 UInt_t x,y,z;
1147
1148 TH3D *h = new TH3D("", "", Int_t(Nx), 0, Nx-1, Int_t(Ny), 0, Ny-1, Int_t(Nz), 0, Nz-1);
1149
1150
1151 for (z = 0; z < Nz; z++) {
1152 UInt_t offset1 = z * Nx * Ny;
1153 for (y = 0; y < Ny; y++) {
1154 UInt_t offset2 = y * Nx;
1155 for (x = 0; x < Nx; x++) {
1156 Int_t nentries = Int_t(fPixels->GetAt(offset1 + offset2 + x));
1157 if (nentries < 0) nentries = 0; //Crop negative values
1158 h->Fill(x, y, z, nentries);
1159 }
1160 }
1161 }
1162
1163 result = h;
1164 }
1165
1166 return result;
1167}
1168
1169////////////////////////////////////////////////////////////////////////////////
1170/// Get a row from the image HDU when it's a 2D array.
1171
1173{
1174 if (fType != kImageHDU) {
1175 Warning("GetArrayRow", "this is not an image HDU.");
1176 return 0;
1177 }
1178
1179 if (fSizes->GetSize() != 2) {
1180 Warning("GetArrayRow", "could not get row from HDU because it has %d dimensions.", fSizes->GetSize());
1181 return 0;
1182 }
1183
1184 UInt_t i, offset, W,H;
1185
1186 W = UInt_t(fSizes->GetAt(0));
1187 H = UInt_t(fSizes->GetAt(1));
1188
1189
1190 if (row >= H) {
1191 Warning("GetArrayRow", "index out of bounds.");
1192 return 0;
1193 }
1194
1195 offset = W * row;
1196 double *v = new double[W];
1197
1198 for (i = 0; i < W; i++) {
1199 v[i] = fPixels->GetAt(offset+i);
1200 }
1201
1202 TVectorD *vec = new TVectorD(W, v);
1203
1204 delete [] v;
1205
1206 return vec;
1207}
1208
1209////////////////////////////////////////////////////////////////////////////////
1210/// Get a column from the image HDU when it's a 2D array.
1211
1213{
1214 if (fType != kImageHDU) {
1215 Warning("GetArrayColumn", "this is not an image HDU.");
1216 return 0;
1217 }
1218
1219 if (fSizes->GetSize() != 2) {
1220 Warning("GetArrayColumn", "could not get row from HDU because it has %d dimensions.", fSizes->GetSize());
1221 return 0;
1222 }
1223
1224 UInt_t i, W,H;
1225
1226 W = UInt_t(fSizes->GetAt(0));
1227 H = UInt_t(fSizes->GetAt(1));
1228
1229
1230 if (col >= W) {
1231 Warning("GetArrayColumn", "index out of bounds.");
1232 return 0;
1233 }
1234
1235 double *v = new double[H];
1236
1237 for (i = 0; i < H; i++) {
1238 v[i] = fPixels->GetAt(W*i+col);
1239 }
1240
1241 TVectorD *vec = new TVectorD(H, v);
1242
1243 delete [] v;
1244
1245 return vec;
1246}
1247
1248
1249////////////////////////////////////////////////////////////////////////////////
1250///Get column number given its name
1251
1253{
1254 Int_t colnum;
1255 for (colnum = 0; colnum < fNColumns; colnum++) {
1256 if (fColumnsInfo[colnum].fName == colname) {
1257 return colnum;
1258 }
1259 }
1260 return -1;
1261}
1262
1263////////////////////////////////////////////////////////////////////////////////
1264/// Get a string-typed column from a table HDU given its column index (>=0).
1265
1267{
1268 if (fType != kTableHDU) {
1269 Warning("GetTabStringColumn", "this is not a table HDU.");
1270 return 0;
1271 }
1272
1273 if ((colnum < 0) || (colnum >= fNColumns)) {
1274 Warning("GetTabStringColumn", "column index out of bounds.");
1275 return 0;
1276 }
1277
1278 if (fColumnsInfo[colnum].fType != kString) {
1279 Warning("GetTabStringColumn", "attempting to read a column that is not of type 'kString'.");
1280 return 0;
1281 }
1282
1283 Int_t offset = colnum * fNRows;
1284
1285 TObjArray *res = new TObjArray();
1286 for (Int_t row = 0; row < fNRows; row++) {
1287 res->Add(new TObjString(fCells[offset + row].fString));
1288 }
1289
1290 return res;
1291}
1292
1293////////////////////////////////////////////////////////////////////////////////
1294/// Get a string-typed column from a table HDU given its name
1295
1297{
1298 if (fType != kTableHDU) {
1299 Warning("GetTabStringColumn", "this is not a table HDU.");
1300 return 0;
1301 }
1302
1303
1304 Int_t colnum = GetColumnNumber(colname);
1305
1306 if (colnum == -1) {
1307 Warning("GetTabStringColumn", "column not found.");
1308 return 0;
1309 }
1310
1311 if (fColumnsInfo[colnum].fType != kString) {
1312 Warning("GetTabStringColumn", "attempting to read a column that is not of type 'kString'.");
1313 return 0;
1314 }
1315
1316 Int_t offset = colnum * fNRows;
1317
1318 TObjArray *res = new TObjArray();
1319 for (Int_t row = 0; row < fNRows; row++) {
1320 res->Add(new TObjString(fCells[offset + row].fString));
1321 }
1322
1323 return res;
1324}
1325
1326////////////////////////////////////////////////////////////////////////////////
1327/// Get a real number-typed column from a table HDU given its column index (>=0).
1328
1330{
1331 if (fType != kTableHDU) {
1332 Warning("GetTabRealVectorColumn", "this is not a table HDU.");
1333 return 0;
1334 }
1335
1336 if ((colnum < 0) || (colnum >= fNColumns)) {
1337 Warning("GetTabRealVectorColumn", "column index out of bounds.");
1338 return 0;
1339 }
1340
1341 if (fColumnsInfo[colnum].fType == kRealArray) {
1342 Warning("GetTabRealVectorColumn", "attempting to read a column whose cells have embedded fixed-length arrays");
1343 Info("GetTabRealVectorColumn", "Use GetTabRealVectorCells() or GetTabRealVectorCell() instead.");
1344 return 0;
1345 } else if (fColumnsInfo[colnum].fType == kRealVector) {
1346 Warning("GetTabRealVectorColumn", "attempting to read a column whose cells have embedded variable-length arrays");
1347 Info("GetTabRealVectorColumn", "Use GetTabVarLengthCell() instead.");
1348 return 0;
1349 }
1350
1351 Int_t offset = colnum * fNRows;
1352
1353 Double_t *arr = new Double_t[fNRows];
1354
1355 for (Int_t row = 0; row < fNRows; row++) {
1356 arr[row] = fCells[offset + row].fRealNumber;
1357 }
1358
1359 TVectorD *res = new TVectorD();
1360 res->Use(fNRows, arr);
1361
1362 return res;
1363}
1364
1365////////////////////////////////////////////////////////////////////////////////
1366/// Get a real number-typed column from a table HDU given its name
1367
1369{
1370 if (fType != kTableHDU) {
1371 Warning("GetTabRealVectorColumn", "this is not a table HDU.");
1372 return 0;
1373 }
1374
1375 Int_t colnum = GetColumnNumber(colname);
1376
1377 if (colnum == -1) {
1378 Warning("GetTabRealVectorColumn", "column not found.");
1379 return 0;
1380 }
1381
1382 if (fColumnsInfo[colnum].fType == kRealArray) {
1383 Warning("GetTabRealVectorColumn", "attempting to read a column whose cells have embedded fixed-length arrays");
1384 Info("GetTabRealVectorColumn", "Use GetTabRealVectorCells() or GetTabRealVectorCell() instead.");
1385 return 0;
1386 } else if (fColumnsInfo[colnum].fType == kRealVector) {
1387 Warning("GetTabRealVectorColumn", "attempting to read a column whose cells have embedded variable-length arrays");
1388 Info("GetTabRealVectorColumn", "Use GetTabVarLengthCell() instead.");
1389 return 0;
1390 }
1391
1392 Int_t offset = colnum * fNRows;
1393
1394 Double_t *arr = new Double_t[fNRows];
1395
1396 for (Int_t row = 0; row < fNRows; row++) {
1397 arr[row] = fCells[offset + row].fRealNumber;
1398 }
1399
1400 TVectorD *res = new TVectorD();
1401 res->Use(fNRows, arr);
1402
1403 return res;
1404}
1405
1406////////////////////////////////////////////////////////////////////////////////
1407/// Change to another HDU given by "filter".
1408/// The parameter "filter" will be appended to the
1409/// FITS file's base path. For example:
1410/// hduObject.Change("[EVENTS][TIME > 5]");
1411/// Please, see documentation of TFITSHDU(const char *filepath_with_filter) constructor
1412/// for further information.
1413
1414Bool_t TFITSHDU::Change(const char *filter)
1415{
1416 TString tmppath;
1417 tmppath.Form("%s%s", fBaseFilePath.Data(), filter);
1418
1421
1422 if (kFALSE == LoadHDU(tmppath)) {
1423 //Failed! Restore previous hdu
1424 Warning("Change", "error changing HDU. Restoring the previous one...");
1425
1428
1429 if (kFALSE == LoadHDU(fFilePath)) {
1430 Warning("Change", "could not restore previous HDU. This object is no longer reliable!!");
1431 }
1432 return kFALSE;
1433 }
1434
1435 //Set new full path
1436 fFilePath = tmppath;
1437 return kTRUE;
1438}
1439
1440////////////////////////////////////////////////////////////////////////////////
1441/// Change to another HDU given by extension_number
1442
1444{
1445 TString tmppath;
1446 tmppath.Form("[%d]", extension_number);
1447
1448 return Change(tmppath.Data());
1449}
1450
1451////////////////////////////////////////////////////////////////////////////////
1452/// Get a collection of real vectors embedded in cells along a given column from a table HDU. colnum >= 0.
1453
1455{
1456 if (fType != kTableHDU) {
1457 Warning("GetTabRealVectorCells", "this is not a table HDU.");
1458 return 0;
1459 }
1460
1461 if ((colnum < 0) || (colnum >= fNColumns)) {
1462 Warning("GetTabRealVectorCells", "column index out of bounds.");
1463 return 0;
1464 }
1465
1466 if (fColumnsInfo[colnum].fType == kRealVector) {
1467 Warning("GetTabRealVectorCells", "attempting to read a column whose cells have embedded variable-length arrays");
1468 Info("GetTabRealVectorCells", "Use GetTabVarLengthCell() instead.");
1469 return 0;
1470 }
1471
1472 Int_t offset = colnum * fNRows;
1473
1474 TObjArray *res = new TObjArray();
1475 Int_t dim = fColumnsInfo[colnum].fDim;
1476
1477 for (Int_t row = 0; row < fNRows; row++) {
1478 TVectorD *v = new TVectorD();
1479 v->Use(dim, fCells[offset + row].fRealArray);
1480 res->Add(v);
1481 }
1482
1483 //Make the collection to own the allocated TVectorD objects, so when
1484 //destroying the collection, the vectors will be destroyed too.
1485 res->SetOwner(kTRUE);
1486
1487 return res;
1488}
1489
1490////////////////////////////////////////////////////////////////////////////////
1491/// Get a collection of real vectors embedded in cells along a given column from a table HDU by name
1492
1494{
1495 if (fType != kTableHDU) {
1496 Warning("GetTabRealVectorCells", "this is not a table HDU.");
1497 return 0;
1498 }
1499
1500 Int_t colnum = GetColumnNumber(colname);
1501
1502 if (colnum == -1) {
1503 Warning("GetTabRealVectorCells", "column not found.");
1504 return 0;
1505 }
1506
1507 return GetTabRealVectorCells(colnum);
1508}
1509
1510////////////////////////////////////////////////////////////////////////////////
1511/// Get a real array (with fixed or variable-length) embedded in a cell given by (row>=0, column>=0)
1512
1514{
1515 if (fType != kTableHDU) {
1516 Warning("GetTabRealVectorCell", "this is not a table HDU.");
1517 return 0;
1518 }
1519
1520 if ((colnum < 0) || (colnum >= fNColumns)) {
1521 Warning("GetTabRealVectorCell", "column index out of bounds.");
1522 return 0;
1523 }
1524
1525 if ((rownum < 0) || (rownum >= fNRows)) {
1526 Warning("GetTabRealVectorCell", "row index out of bounds.");
1527 return 0;
1528 }
1529
1530 if (fColumnsInfo[colnum].fType == kRealVector) {
1531 Warning("GetTabRealVectorCells", "attempting to read a column whose cells have embedded variable-length arrays");
1532 Info("GetTabRealVectorCells", "Use GetTabVarLengthCell() instead.");
1533 return 0;
1534 }
1535
1536 TVectorD *v = new TVectorD();
1537 v->Use(fColumnsInfo[colnum].fDim, fCells[(colnum * fNRows) + rownum].fRealArray);
1538 return v;
1539}
1540
1541////////////////////////////////////////////////////////////////////////////////
1542/// Get a real vector embedded in a cell given by (row>=0, column name)
1543
1544TVectorD *TFITSHDU::GetTabRealVectorCell(Int_t rownum, const char *colname)
1545{
1546 if (fType != kTableHDU) {
1547 Warning("GetTabRealVectorCell", "this is not a table HDU.");
1548 return 0;
1549 }
1550
1551 Int_t colnum = GetColumnNumber(colname);
1552
1553 if (colnum == -1) {
1554 Warning("GetTabRealVectorCell", "column not found.");
1555 return 0;
1556 }
1557
1558 return GetTabRealVectorCell(rownum, colnum);
1559}
1560
1561////////////////////////////////////////////////////////////////////////////////
1562/// Get the name of a column given its index (column>=0).
1563/// In case of error the column name is "".
1564
1566{
1567 static TString noName;
1568 if (fType != kTableHDU) {
1569 Error("GetColumnName", "this is not a table HDU.");
1570 return noName;
1571 }
1572
1573 if ((colnum < 0) || (colnum >= fNColumns)) {
1574 Error("GetColumnName", "column index out of bounds.");
1575 return noName;
1576 }
1577 return fColumnsInfo[colnum].fName;
1578}
1579
1580////////////////////////////////////////////////////////////////////////////////
1581/// Get the variable-length array contained in a cell given by (row>=0, column name)
1582
1584
1585 if (fType != kTableHDU) {
1586 Warning("GetTabVarLengthVectorCell", "this is not a table HDU.");
1587 return 0;
1588 }
1589
1590 if ((colnum < 0) || (colnum >= fNColumns)) {
1591 Warning("GetTabVarLengthVectorCell", "column index out of bounds.");
1592 return 0;
1593 }
1594
1595 if ((rownum < 0) || (rownum >= fNRows)) {
1596 Warning("GetTabVarLengthVectorCell", "row index out of bounds.");
1597 return 0;
1598 }
1599
1600 return fCells[(colnum * fNRows) + rownum].fRealVector;
1601}
1602
1603////////////////////////////////////////////////////////////////////////////////
1604/// Get the variable-length array contained in a cell given by (row>=0, column name)
1605
1607{
1608 if (fType != kTableHDU) {
1609 Warning("GetTabVarLengthVectorCell", "this is not a table HDU.");
1610 return 0;
1611 }
1612
1613 Int_t colnum = GetColumnNumber(colname);
1614
1615 if (colnum == -1) {
1616 Warning("GetTabVarLengthVectorCell", "column not found.");
1617 return 0;
1618 }
1619
1620 return GetTabVarLengthVectorCell(rownum, colnum);
1621}
double
Definition: Converters.cxx:921
long
Definition: Converters.cxx:858
#define h(i)
Definition: RSha256.hxx:106
const Ssiz_t kNPOS
Definition: RtypesCore.h:113
int Int_t
Definition: RtypesCore.h:43
unsigned int UInt_t
Definition: RtypesCore.h:44
const Bool_t kFALSE
Definition: RtypesCore.h:90
double Double_t
Definition: RtypesCore.h:57
const Bool_t kTRUE
Definition: RtypesCore.h:89
const char Option_t
Definition: RtypesCore.h:64
#define ClassImp(name)
Definition: Rtypes.h:361
include TDocParser_001 C image html pict1_TDocParser_001 png width
Definition: TDocParser.cxx:121
int nentries
Definition: THbookFile.cxx:89
TMatrixT< Double_t > TMatrixD
Definition: TMatrixDfwd.h:22
TVectorT< Double_t > TVectorD
Definition: TVectorDfwd.h:22
Array of doubles (64 bits per element).
Definition: TArrayD.h:27
void SetAt(Double_t v, Int_t i)
Definition: TArrayD.h:51
Double_t GetAt(Int_t i) const
Definition: TArrayD.h:45
Array of integers (32 bits per element).
Definition: TArrayI.h:27
void SetAt(Double_t v, Int_t i)
Definition: TArrayI.h:51
Double_t GetAt(Int_t i) const
Definition: TArrayI.h:45
Int_t GetSize() const
Definition: TArray.h:47
The Canvas class.
Definition: TCanvas.h:27
virtual void SetOwner(Bool_t enable=kTRUE)
Set whether this collection is the owner (enable==true) of its content.
FITS file interface class.
Definition: TFITS.h:35
union Cell * fCells
Table cells (when fType == kTableHDU).
Definition: TFITS.h:90
void PrintColumnInfo(const Option_t *) const
Print column information.
Definition: TFITS.cxx:766
void _initialize_me()
Do some initializations.
Definition: TFITS.cxx:222
void Draw(Option_t *opt="")
If the HDU is an image, draw the first layer of the primary array To set a title to the canvas,...
Definition: TFITS.cxx:973
void PrintFullTable(const Option_t *) const
Print full table contents.
Definition: TFITS.cxx:794
TString & GetKeywordValue(const char *keyword)
Get the value of a given keyword. Return "" if not found.
Definition: TFITS.cxx:638
@ kImageHDU
Definition: TFITS.h:43
@ kTableHDU
Definition: TFITS.h:44
Int_t fNRows
Number of rows (when fType == kTableHDU)
Definition: TFITS.h:89
const TString & GetColumnName(Int_t colnum)
Get the name of a column given its index (column>=0).
Definition: TFITS.cxx:1565
Int_t fNRecords
Number of records.
Definition: TFITS.h:81
TString fFilePath
Path to HDU's file including filter.
Definition: TFITS.h:78
struct Column * fColumnsInfo
Information about columns (when fType == kTableHDU)
Definition: TFITS.h:87
TFITSHDU(const char *filepath_with_filter)
TFITSHDU constructor from file path with HDU selection filter.
Definition: TFITS.cxx:118
Int_t fNColumns
Number of columns (when fType == kTableHDU)
Definition: TFITS.h:88
TImage * ReadAsImage(Int_t layer=0, TImagePalette *pal=0)
Read image HDU as a displayable image.
Definition: TFITS.cxx:880
TVectorD * GetTabRealVectorColumn(Int_t colnum)
Get a real number-typed column from a table HDU given its column index (>=0).
Definition: TFITS.cxx:1329
@ kString
Definition: TFITS.h:48
@ kRealArray
Definition: TFITS.h:50
@ kRealVector
Definition: TFITS.h:51
@ kRealNumber
Definition: TFITS.h:49
TArrayI * fSizes
Image sizes in each dimension (when fType == kImageHDU)
Definition: TFITS.h:85
Bool_t Change(const char *filter)
Change to another HDU given by "filter".
Definition: TFITS.cxx:1414
TMatrixD * ReadAsMatrix(Int_t layer=0, Option_t *opt="")
Read image HDU as a matrix.
Definition: TFITS.cxx:1000
Int_t fNumber
HDU number (1=PRIMARY)
Definition: TFITS.h:84
void PrintHDUMetadata(const Option_t *opt="") const
Print records.
Definition: TFITS.cxx:651
~TFITSHDU()
TFITSHDU destructor.
Definition: TFITS.cxx:169
struct HDURecord * fRecords
HDU metadata records.
Definition: TFITS.h:80
TVectorD * GetArrayRow(UInt_t row)
Get a row from the image HDU when it's a 2D array.
Definition: TFITS.cxx:1172
void _release_resources()
Release internal resources.
Definition: TFITS.cxx:177
void PrintFileMetadata(const Option_t *opt="") const
Print HDU's parent file's metadata.
Definition: TFITS.cxx:665
static void CleanFilePath(const char *filepath_with_filter, TString &dst)
Clean path from possible filter and put the result in 'dst'.
Definition: TFITS.cxx:94
struct HDURecord * GetRecord(const char *keyword)
Get record by keyword.
Definition: TFITS.cxx:625
void Print(const Option_t *opt="") const
Print metadata.
Definition: TFITS.cxx:858
TObjArray * GetTabStringColumn(Int_t colnum)
Get a string-typed column from a table HDU given its column index (>=0).
Definition: TFITS.cxx:1266
TString fBaseFilePath
Path to HDU's file excluding filter.
Definition: TFITS.h:79
TH1 * ReadAsHistogram()
Read image HDU as a histogram.
Definition: TFITS.cxx:1093
TObjArray * GetTabRealVectorCells(Int_t colnum)
Get a collection of real vectors embedded in cells along a given column from a table HDU....
Definition: TFITS.cxx:1454
TArrayD * fPixels
Image pixels (when fType == kImageHDU)
Definition: TFITS.h:86
TVectorD * GetArrayColumn(UInt_t col)
Get a column from the image HDU when it's a 2D array.
Definition: TFITS.cxx:1212
enum EHDUTypes fType
HDU type.
Definition: TFITS.h:82
TVectorD * GetTabRealVectorCell(Int_t rownum, Int_t colnum)
Get a real array (with fixed or variable-length) embedded in a cell given by (row>=0,...
Definition: TFITS.cxx:1513
Bool_t LoadHDU(TString &filepath_filter)
Load HDU from fits file satisfying the specified filter.
Definition: TFITS.cxx:237
TString fExtensionName
Extension Name.
Definition: TFITS.h:83
TArrayD * GetTabVarLengthVectorCell(Int_t rownum, Int_t colnum)
Get the variable-length array contained in a cell given by (row>=0, column name)
Definition: TFITS.cxx:1583
Int_t GetColumnNumber(const char *colname)
Get column number given its name.
Definition: TFITS.cxx:1252
1-D histogram with a double per channel (see TH1 documentation)}
Definition: TH1.h:614
The TH1 histogram class.
Definition: TH1.h:56
2-D histogram with a double per channel (see TH1 documentation)}
Definition: TH2.h:292
3-D histogram with a double per channel (see TH1 documentation)}
Definition: TH3.h:304
A class to define a conversion from pixel values to pixel color.
Definition: TAttImage.h:33
UShort_t * fColorRed
[fNumPoints] red color at each anchor point
Definition: TAttImage.h:38
Double_t * fPoints
[fNumPoints] value of each anchor point [0..1]
Definition: TAttImage.h:37
UShort_t * fColorGreen
[fNumPoints] green color at each anchor point
Definition: TAttImage.h:39
UShort_t * fColorBlue
[fNumPoints] blue color at each anchor point
Definition: TAttImage.h:40
UShort_t * fColorAlpha
[fNumPoints] alpha at each anchor point
Definition: TAttImage.h:41
An abstract interface to image processing library.
Definition: TImage.h:29
virtual void SetImage(const Double_t *, UInt_t, UInt_t, TImagePalette *=0)
Definition: TImage.h:116
static TImage * Create()
Create an image.
Definition: TImage.cxx:36
virtual const Element * GetMatrixArray() const
Definition: TMatrixT.h:222
TString fName
Definition: TNamed.h:32
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
An array of TObjects.
Definition: TObjArray.h:37
void Add(TObject *obj)
Definition: TObjArray.h:74
Collectable string class.
Definition: TObjString.h:28
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
Definition: TObject.cxx:877
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:891
virtual void Draw(Option_t *option="")
Default Draw method for all objects.
Definition: TObject.cxx:195
virtual void Info(const char *method, const char *msgfmt,...) const
Issue info message.
Definition: TObject.cxx:865
Basic string class.
Definition: TString.h:131
Ssiz_t Length() const
Definition: TString.h:405
const char * Data() const
Definition: TString.h:364
void Resize(Ssiz_t n)
Resize the string. Truncate or add blanks as necessary.
Definition: TString.cxx:1095
@ kExact
Definition: TString.h:263
void Form(const char *fmt,...)
Formats a string using a printf style format descriptor.
Definition: TString.cxx:2289
Ssiz_t Index(const char *pat, Ssiz_t i=0, ECaseCompare cmp=kExact) const
Definition: TString.h:634
TVectorT< Element > & Use(Int_t lwb, Int_t upb, Element *data)
Use the array data to fill the vector lwb..upb].
Definition: TVectorT.cxx:350
Double_t y[n]
Definition: legend1.C:17
Double_t x[n]
Definition: legend1.C:17
#define H(x, y, z)
static const std::string comment("comment")
Short_t Abs(Short_t d)
Definition: TMathBase.h:120
std::vector< int > fRowStart
Definition: TFITS.h:66
enum EColumnTypes fType
Definition: TFITS.h:62
TString fName
Definition: TFITS.h:61
Int_t fDim
Definition: TFITS.h:63
std::vector< int > fVarLengths
Definition: TFITS.h:67
TString fKeyword
Definition: TFITS.h:55
TString fValue
Definition: TFITS.h:56
TString fComment
Definition: TFITS.h:57
TArrayD * fRealVector
Definition: TFITS.h:74
Double_t * fRealArray
Definition: TFITS.h:73
Double_t fRealNumber
Definition: TFITS.h:72
Char_t * fString
Definition: TFITS.h:71