NcFITSIO.cxx

Go to the documentation of this file.

 
#include "NcFITSIO.h"
#include "Riostream.h"
 
ClassImp(NcFITSIO); // Class implementation to enable ROOT I/O

NcFITSIO::NcFITSIO(const char* name,const char* title) : TNamed(name,title)
{
 
 Initialize();
}

NcFITSIO::~NcFITSIO()
{
 
 Reset();
}

NcFITSIO::NcFITSIO(const NcFITSIO& q) : TNamed(q)
{
 
 Initialize();
 fFilename=q.fFilename;
 fFilenameFilter=q.fFilenameFilter;
 
 if (!LoadHeaderInfo())
 {
  cout << endl;
  cout << " *" << ClassName() << "::NcFITSIO* Failure in copy constructor." << endl;
 }
}

void NcFITSIO::Initialize()
{
 
 fFilename="";           // The (full path) name of the FITS file on the computer system
 fFilenameFilter="";     // The FITS filename with the HDU selection filter
 fInput=0;               // Pointer to the FITS input file
 fOutput=0;              // Pointer to the FITS output file
 fType=kImageHDU;        // The HDU type
 fExtensionName="";      // Extension name of the HDU
 fExtensionNumber=0;     // Extension number of the HDU (0=PRIMARY)
 fNkeys=0;               // The number of HDU keywords
 fKeyNames=0;            // The HDU key names
 fKeyValues=0;           // The HDU key values
 fComments=0;            // The HDU (key) comments
 fNrows=0;               // The number of table rows
 fNcolumns=0;            // The number of table columns
 fColumnNames=0;         // The names of the table columns
 fColumnTypes=0;         // The types of the table columns
 fColumnLayers=0;        // The number of layers of the table column
 fSizes=0;               // Image sizes in each dimension (when fType == kImageHDU)
}

void NcFITSIO::Reset()
{
 
 if (fKeyNames)
 {
  delete[] fKeyNames;
  fKeyNames=0;
 }
 if (fKeyValues)
 {
  delete[] fKeyValues;
  fKeyValues=0;
 }
 if (fComments)
 {
  delete[] fComments;
  fComments=0;
 }
 if (fColumnNames)
 {
  delete[] fColumnNames;
  fColumnNames=0;
 }
 if (fColumnTypes)
 {
  delete[] fColumnTypes;
  fColumnTypes=0;
 }
 if (fColumnLayers)
 {
  delete[] fColumnLayers;
  fColumnLayers=0;
 }
 if (fSizes)
 {
  delete fSizes;
  fSizes=0;
 }
 
 // Reset all parameters
 Initialize();
}

Bool_t NcFITSIO::OpenInputFile(TString specs)
{
 
 TString file=gSystem->ExpandPathName(specs.Data());
   
 Initialize();
   
 // Store both the filename with the optional HDU selection filter appended and the plain filename
 fFilenameFilter=file;
 fFilename=StripFilter(fFilenameFilter);
 
 Bool_t good=LoadHeaderInfo();
 
 if (good)
 {
  cout << endl;
  cout << " *" << ClassName() << "::OpenInputFile* FITS file specs: " << fFilenameFilter << endl;
 }
 else
 {
  cout << endl;
  cout << " *" << ClassName() << "::OpenInputFile* Could not open " << fFilenameFilter << endl;
  Reset();
 }
 
 return good;  
}

TString NcFITSIO::StripFilter(TString filename) const
{
 
 TString s=filename;
 Int_t idx=s.Index("[",1,0,TString::kExact);
 if (idx>=0) s.Resize(idx);
 
 return s;
}

Bool_t NcFITSIO::LoadHeaderInfo()
{
 
 Bool_t good=kTRUE;
 fInput=0;
 int status=0;
 
 if (fKeyNames)
 {
  delete[] fKeyNames;
  fKeyNames=0;
 }
 if (fKeyValues)
 {
  delete[] fKeyValues;
  fKeyValues=0;
 }
 if (fComments)
 {
  delete[] fComments;
  fComments=0;
 }
 
 // Open the FITS file as specified via OpenInputFile() or SelectHDU()
 fits_open_file(&fInput,fFilenameFilter.Data(),READONLY,&status);
 
 if (status)
 {
  cout << endl;
  cout << " *" << ClassName() << "::LoadHeaderInfo* Could not open " << fFilenameFilter << endl;
  Reset();
  good=kFALSE;
  return good;
 }
 
 // Read the number of this HDU
 int hdunum;
 fits_get_hdu_num(fInput,&hdunum);
 fExtensionNumber=Int_t(hdunum)-1;
 
 // Read the type of this HDU
 int hdutype;
 fits_get_hdu_type(fInput,&hdutype,&status);
 
 if (status)
 {
  cout << endl;
  cout << " *" << ClassName() << "::LoadHeaderInfo* Could not retrieve the HDU type." << endl;
  fits_close_file(fInput,&status);
  Reset();
  good=kFALSE;
  return good;
 }
 
 fType=kTableHDU;
 if (hdutype==IMAGE_HDU) fType=kImageHDU;
 
 // Read the HDU header records
 int nkeys,morekeys;
 fits_get_hdrspace(fInput,&nkeys,&morekeys,&status);
 
 if (status)
 {
  cout << endl;
  cout << " *" << ClassName() << "::LoadHeaderInfo* Could not retrieve the HDU header space." << endl;
  fits_close_file(fInput,&status);
  Reset();
  good=kFALSE;
  return good;
 }
 
 // Store the HDU header information
 fKeyNames=new TString[nkeys];
 fKeyValues=new TString[nkeys];
 fComments=new TString[nkeys];
 
 char keyname[FLEN_KEYWORD+1];
 char keyvalue[FLEN_VALUE+1];
 char comment[FLEN_COMMENT+1];
 fNkeys=0;
 for (Int_t i=1; i<=nkeys; i++)
 {
  fits_read_keyn(fInput,i,keyname,keyvalue,comment,&status);
 
  if (status)
  {
   cout << endl;
   cout << " *" << ClassName() << "::LoadHeaderInfo* Could not retrieve info of HDU key record "<< i << endl;
   fits_close_file(fInput,&status);
   Reset(); 
   good=kFALSE;
   return good;
  }
 
  fKeyNames[i-1]=keyname;
  fKeyValues[i-1]=keyvalue;
  fComments[i-1]=comment;
 
  fNkeys++;
 
  // Obtain the extension name
  fExtensionName="[";
   if (fKeyNames[i-1]=="EXTNAME")
  {
   fExtensionName+=fKeyValues[i-1];
   fExtensionName.ReplaceAll("'",""); // Remove the single quotes
  }
  else
  {
   fExtensionName+=fExtensionNumber;
  }
  fExtensionName+="]";
 }
 
 // Obtain row and column information in case of Table data
 if (fType==kTableHDU)
 {
  long nrows;
  fits_get_num_rows(fInput,&nrows,&status);
 
  if (status)
  {
   cout << endl;
   cout << " *" << ClassName() << "::LoadHeaderInfo* Could not retrieve number of table rows." << endl;
   fits_close_file(fInput,&status);
   Reset(); 
   good=kFALSE;
   return good;
  }
      
  fNrows=Int_t(nrows);
      
  int ncols;
  fits_get_num_cols(fInput,&ncols,&status);
 
  if (status)
  {
   cout << endl;
   cout << " *" << ClassName() << "::LoadHeaderInfo* Could not retrieve number of table columns." << endl;
   fits_close_file(fInput,&status);
   Reset(); 
   good=kFALSE;
   return good;
  }
 
  // Store the Table column information
  fNcolumns=Int_t(ncols);
  fColumnNames=new TString[ncols];
  fColumnTypes=new eColumnTypes[ncols];
  fColumnLayers=new Int_t[ncols];
            
  // Read the column names
  status=0;
  char colname[80];     
  int jcol;
  fits_get_colname(fInput,CASEINSEN,(char*)"*",colname,&jcol,&status);
 
  if (status==COL_NOT_FOUND)
  {
   cout << endl;
   cout << " *" << ClassName() << "::LoadHeaderInfo* Could not find any table column." << endl;
   fits_close_file(fInput,&status);
   Reset(); 
   good=kFALSE;
   return good;
  }
 
  if (jcol>0 && jcol<=ncols) fColumnNames[jcol-1]=colname;
 
  while (status==COL_NOT_UNIQUE)
  {
   fits_get_colname(fInput,CASEINSEN,(char*)"*",colname,&jcol,&status);
   if (jcol>0 && jcol<=ncols) fColumnNames[jcol-1]=colname;
  }
  if (status && status!=COL_NOT_FOUND)
  {
   cout << endl;
   cout << " *" << ClassName() << "::LoadHeaderInfo* Could not retrieve name of table column " << jcol << endl;
   fits_close_file(fInput,&status);
   Reset(); 
   good=kFALSE;
   return good;
  }
 
  // Read the column data types
  status=0;
  int typecode;
  long repeat,width;
  Int_t dim=1; // Dimension of the stored column elements
  for (jcol=1; jcol<=fNcolumns; jcol++)
  {
   fits_get_coltype(fInput,jcol,&typecode,&repeat,&width,&status);
 
   if (status)
   {
    cout << endl;
    cout << " *" << ClassName() << "::LoadHeaderInfo* Could not retrieve data type of table column " << jcol << endl;
    fits_close_file(fInput,&status);
    Reset(); 
    good=kFALSE;
    return good;
   }
 
   // Determine the dimension of the stored column elements
   dim=Int_t(repeat);
   if (typecode==TSTRING) dim=Int_t(repeat/width);
   if (dim<=0) dim=1;
 
   fColumnLayers[jcol-1]=dim;
 
   if (typecode==TSTRING)
   {
    fColumnTypes[jcol-1]=kString;
    if (dim>1) fColumnTypes[jcol-1]=kStringArray;
   }
   else if (typecode==TLOGICAL)
   {
    fColumnTypes[jcol-1]=kLogical;
    if (dim>1) fColumnTypes[jcol-1]=kLogicalArray;
   }
   else if (typecode==TCOMPLEX || typecode==TDBLCOMPLEX)
   {
    fColumnTypes[jcol-1]=kComplexNumber;
    if (dim>1) fColumnTypes[jcol-1]=kComplexArray;
   }
   else
   {
    fColumnTypes[jcol-1]=kRealNumber;
    if (dim>1) fColumnTypes[jcol-1]=kRealArray;
   }
  }
 } // End of Table info
 
 // Obtain dimension and size information in case of Image data
 if (fType==kImageHDU)
 {
  // The number of dimensions of the N-dimensional data 
  status=0;
  int ndims=0;
  fits_get_img_dim(fInput,&ndims,&status);
 
  if (status)
  {
   cout << endl;
   cout << " *" << ClassName() << "::LoadHeaderInfo* Could not retrieve the number of dimensions of the Image." << endl;
   fits_close_file(fInput,&status);
   Reset(); 
   good=kFALSE;
   return good;
  }
 
  // Empty Image
  if (ndims<=0)
  {
   fSizes=new TArrayI();
   return good;
  }
 
  // The size of each dimension
  long* dimsizes=new long[ndims];
  fits_get_img_size(fInput,ndims,dimsizes,&status);
 
  if (status)
  {
   cout << endl;
   cout << " *" << ClassName() << "::LoadHeaderInfo* Could not retrieve the sizes of the dimensions of the Image." << endl;
   fits_close_file(fInput,&status);
   Reset(); 
   good=kFALSE;
   delete[] dimsizes;
   return good;
  }
 
  // Store the size of each dimension
  fSizes=new TArrayI(ndims);
  for (Int_t i= 0; i<ndims; i++)
  {
   fSizes->SetAt(dimsizes[i],i);
  }
 
  delete [] dimsizes;
 } // End of Image info
 
 return good;
}

Bool_t NcFITSIO::SelectHDU(TString extname)
{
 
 fFilenameFilter=fFilename;
 fFilenameFilter+=extname;
 
 Bool_t good=LoadHeaderInfo();
 
 if (good)
 {
  cout << endl;
  cout << " *" << ClassName() << "::SelectHDU* Current selection: " << fFilenameFilter << endl;
 }
 else
 {
  cout << endl;
  cout << " *" << ClassName() << "::SelectHDU* Could not select " << fFilenameFilter << endl;
  Reset();
 }
 
 return good;
}

Bool_t NcFITSIO::SelectHDU(Int_t extnumber)
{
 
 TString extname="[";
 extname+=extnumber;
 extname+="]";
 
 Bool_t good=SelectHDU(extname);
 
 return good;
}

TString NcFITSIO::GetKeywordValue(TString keyname,Int_t mode)
{
 
 
 TString value="";
 
 for (Int_t i=0; i<fNkeys; i++)
 {
  if (fKeyNames[i]==keyname || (mode && fKeyNames[i].Contains(keyname.Data())))
  {
   value=fKeyValues[i];
   break;
  } 
 }
 
 value.ReplaceAll("'","");          // Remove the single quotes
 value=value.Strip(value.kBoth,' '); // Remove leading and trailing blanks
 
 return value;
}

Bool_t NcFITSIO::IsTable() const
{
 
 
 Bool_t flag=kFALSE;
 if (fType==kTableHDU) flag=kTRUE;
 
 return flag; 
}

Int_t NcFITSIO::GetHDUCount() const
{
 
 Int_t nhdus=0;
 
 // Open a local mode of the FITS file without any filters
 fitsfile* fp=0;
 int status=0;
 fits_open_file(&fp,fFilename.Data(),READONLY,&status);
 
 if (status)
 {
  cout << " *" << ClassName() << "::GetHDUCount* Could not open file : " << fFilename << endl;
  fits_close_file(fp,&status);
  return 0;
 }
 
 int n=0;
 fits_get_num_hdus(fp,&n,&status);
 
 if (status)
 {
  cout << " *" << ClassName() << "::GetHDUCount* Could not read the number of HDUs" << endl;
  fits_close_file(fp,&status);
  return 0;
 }
 
 nhdus=n;
 
 // Close this local "unfiltered" file mode
 fits_close_file(fp,&status);
 
 return nhdus; 
}

Int_t NcFITSIO::GetTableNrows() const
{
 
 return fNrows;
}

Int_t NcFITSIO::GetTableNcolumns() const
{
 
 return fNcolumns;
}

Int_t NcFITSIO::GetColumnNumber(TString colname,Int_t mode) const
{
   
 Int_t colnum=0;
 Int_t match=0;
 for (Int_t i=1; i<=fNcolumns; i++)
 {
  if (!mode && fColumnNames[i-1]==colname) match=1;
  if (mode && fColumnNames[i-1].Contains(colname)) match=1;
  if (match)
  {
   colnum=i;
   break;
  }
 }
 return colnum;
}

TString NcFITSIO::GetColumnName(Int_t colnum) const
{
 
 TString name="";
 
 if (fType!=kTableHDU || colnum<1 || colnum>fNcolumns) return name;
 
 name=fColumnNames[colnum-1];
 return name;
}

Int_t NcFITSIO::GetTableCell(Double_t &val,Int_t row,Int_t col,Int_t layer)
{
 
 TArrayD arr;
 Int_t ndim=GetTableCell(arr,row,col);
 
 if (!ndim || layer<1 || layer>ndim)
 {
  val=0;
  return 0;
 }
 
 val=arr[layer-1];
 return ndim;
}

Int_t NcFITSIO::GetTableCell(Double_t &val,Int_t row,TString colname,Int_t layer,Int_t mode)
{
 
 Int_t col=GetColumnNumber(colname,mode);
 
 Int_t ndim=GetTableCell(val,row,col);
 
 return ndim;
}

Int_t NcFITSIO::GetTableCell(TArrayD &arr,Int_t row,Int_t col)
{
 
 arr.Set(0);
 
 if (row<=0 || row>fNrows || col<=0 || col>fNcolumns) return 0;
 
 if (fType!=kTableHDU || fColumnTypes[col-1]==kString || fColumnTypes[col-1]==kStringArray) return 0;
 
 // Obtain the number of (array) elements stored in this (row,col) cell.
 Int_t dim=0;
 int colnum=col;
 long long rownum=row;
 long repeat=0;
 long offset=0;
 int status=0;
 fits_read_descript(fInput,colnum,rownum,&repeat,&offset,&status);
 
 if (!status) // The column was created to hold a variable number of elements
 {
  dim=repeat;
 }
 else // The column was created to hold a fixed number of elements
 {
  dim=fColumnLayers[col-1];
 }
 
 if (dim<=0) return 0;
 
 // Read the cell contents into an array
 double nulval=0;
 int anynul=0;
 double* array=new double[dim];
 status=0;
 if (fColumnTypes[col-1]==kLogical || fColumnTypes[col-1]==kLogicalArray)
 {
  bool* barray=new bool[dim];
  bool bnulval=0;
  fits_read_col(fInput,TLOGICAL,col,row,1,dim,&bnulval,barray,&anynul,&status);
  for (Int_t i=0; i<dim; i++)
  {
   array[i]=double(barray[i]);
  }
  delete [] barray;
 }
 else
 {
  fits_read_col(fInput,TDOUBLE,col,row,1,dim,&nulval,array,&anynul,&status);
 }
 
 if (status)
 {
  cout << endl;
  cout << " *" << ClassName() << "::GetTableCell* Could not retrieve data type of table cell [" 
               << row << "," << col << "]." << endl;
  fits_close_file(fInput,&status);
  Reset();
  delete [] array; 
  return 0;
 }
 
 // Copy the cell values into the output array
 arr.Set(dim);
 for (Int_t i=0; i<dim; i++)
 {
  arr.SetAt(array[i],i);
 }
 
 delete [] array;
 return dim;
}

Int_t NcFITSIO::GetTableCell(TArrayD &arr,Int_t row,TString colname,Int_t mode)
{
 
 Int_t col=GetColumnNumber(colname,mode);
 
 Int_t ndim=GetTableCell(arr,row,col);
 
 return ndim;
}

Int_t NcFITSIO::GetTableCell(TString &str,Int_t row,Int_t col,Int_t layer)
{
 
 TString* arr=0;
 Int_t ndim=GetTableCell(arr,row,col);
 if (!ndim || layer<1 || layer>ndim)
 {
  str="";
  ndim=0;
 }
 else
 {
  str=arr[layer-1];
 }
 
 delete[] arr;
 return ndim;
}

Int_t NcFITSIO::GetTableCell(TString &str,Int_t row,TString colname,Int_t layer,Int_t mode)
{
 
 Int_t col=GetColumnNumber(colname,mode);
 
 Int_t ndim=GetTableCell(str,row,col,layer);
 
 return ndim;
}

Int_t NcFITSIO::GetTableCell(TString* &arr,Int_t row,Int_t col)
{
 
 if (!arr) arr=new TString[1];
 arr[0]="";
 
 if (row<=0 || row>fNrows || col<=0 || col>fNcolumns) return 0;
 
 if (fType!=kTableHDU) return 0;
 
 // Obtain the number of (array) elements stored in this (row,col) cell.
 Int_t dim=0;
 int colnum=col;
 long long rownum=row;
 long repeat=0;
 long offset=0;
 int status=0;
 fits_read_descript(fInput,colnum,rownum,&repeat,&offset,&status);
 
 if (!status) // The column was created to hold a variable number of elements
 {
  dim=repeat;
 }
 else // The column was created to hold a fixed number of elements
 {
  dim=fColumnLayers[col-1];
 }
 
 if (dim<=0) return 0;
 
 // Retrieve the character width for this column
 int dispwidth=0;
 status=0;
 fits_get_col_display_width(fInput,col,&dispwidth,&status);
 
 if (status)
 {
  cout << endl;
  cout << " *" << ClassName() << "::GetTableCell* Could not retrieve string width of table cell [" 
               << row << "," << col << "]." << endl;
  fits_close_file(fInput,&status);
  Reset();
  return 0;
 }
 
 if (dispwidth<=0) dispwidth=1;
 
 // Read the cell contents into an array
 char* nulval=(char*)"";
 int anynul=0;
 char** array=new char*[dim];
 for (Int_t i=0; i<dim; i++)
 {
  array[i]=new char[dispwidth+1];
 }
 fits_read_col(fInput,TSTRING,col,row,1,dim,nulval,array,&anynul,&status);
 
 if (status)
 {
  cout << endl;
  cout << " *" << ClassName() << "::GetTableCell* Could not retrieve string of table cell [" 
               << row << "," << col << "]." << endl;
  fits_close_file(fInput,&status);
  Reset();
  delete [] array; 
  return 0;
 }
 
 // Copy the data into the output array
 delete[] arr;
 arr=new TString[dim];
 for (Int_t j=0; j<dim; j++)
 {
  arr[j]=array[j];
 }
 
 delete[] array;
 return dim;
}

Int_t NcFITSIO::GetTableCell(TString* &arr,Int_t row,TString colname,Int_t mode)
{
 
 Int_t col=GetColumnNumber(colname,mode);
 
 Int_t ndim=GetTableCell(arr,row,col);
 
 return ndim;
}

Int_t NcFITSIO::GetTableCell(TObjArray &arr,Int_t row,Int_t col)
{
 
 TString* temp=0;
 Int_t ndim=GetTableCell(temp,row,col);
 
 arr.Clear();
 arr.SetOwner();
 
 TObjString* sobj;
 for (Int_t i=0; i<ndim; i++)
 {
  sobj=new TObjString();
  sobj->SetString(temp[i]);
  arr.Add(sobj);
 }
 
 delete[] temp;
 
 return ndim;
}

Int_t NcFITSIO::GetTableCell(TObjArray &arr,Int_t row,TString colname,Int_t mode)
{
 
 Int_t col=GetColumnNumber(colname,mode);
 
 Int_t ndim=GetTableCell(arr,row,col);
 
 return ndim;
}

Int_t NcFITSIO::GetTableColumn(TArrayD &arr,Int_t col,Int_t rstart,Int_t rend,Int_t layer)
{
 
 if (!rend) rend=fNrows;
 
 arr.Set(0);
 
 if (rstart<=0 || rstart>fNrows || rend<=0 || rend>fNrows || col<=0 || col>fNcolumns || layer<1) return 0;
 
 if (fType!=kTableHDU || fColumnTypes[col-1]==kString || fColumnTypes[col-1]==kStringArray) return 0;
 
 Int_t ndim=(rend-rstart)+1;
 arr.Set(ndim);
 
 Double_t val=0;
 Int_t dim=0;
 Int_t n=0;
 for (Int_t irow=rstart; irow<=rend; irow++)
 {
  dim=GetTableCell(val,irow,col);
 
  if (!dim || layer>dim)
  {
   arr.Set(0);
   return 0;
  }
 
  arr.AddAt(val,n);
  n++;
 }
 
 return n;
}

Int_t NcFITSIO::GetTableColumn(TArrayD &arr,TString colname,Int_t rstart,Int_t rend,Int_t layer,Int_t mode)
{
 
 Int_t col=GetColumnNumber(colname,mode);
 
 Int_t n=GetTableColumn(arr,col,rstart,rend,layer);
 
 return n;
}

Int_t NcFITSIO::GetTableColumn(TString* &arr,Int_t col,Int_t rstart,Int_t rend,Int_t layer)
{
 
 if (!rend) rend=fNrows;
 
 if (!arr) arr=new TString[1];
 arr[0]="";
 
 if (rstart<=0 || rstart>fNrows || rend<=0 || rend>fNrows || col<=0 || col>fNcolumns || layer<1) return 0;
 
 if (fType!=kTableHDU) return 0;
 
 Int_t ndim=(rend-rstart)+1;
 
 delete[] arr;
 arr=new TString[ndim];
 
 TString str;
 Int_t dim=0;
 Int_t n=0;
 for (Int_t irow=rstart; irow<=rend; irow++)
 {
  dim=GetTableCell(str,irow,col,layer);
 
  if (!dim || layer>dim)
  {
   delete[] arr;
   arr=new TString[1];
   arr[0]="";
   return 0;
  }
 
  arr[n]=str;
  n++;
 }
 
 return n;
}

Int_t NcFITSIO::GetTableColumn(TString* &arr,TString colname,Int_t rstart,Int_t rend,Int_t layer,Int_t mode)
{
 
 Int_t col=GetColumnNumber(colname,mode);
 
 Int_t n=GetTableColumn(arr,col,rstart,rend,layer);
 
 return n;
}

Int_t NcFITSIO::GetTableColumn(TObjArray &arr,Int_t col,Int_t rstart,Int_t rend,Int_t layer)
{
 
 if (!rend) rend=fNrows;
 
 arr.Clear();
 arr.SetOwner();
 
 if (rstart<=0 || rstart>fNrows || rend<=0 || rend>fNrows || col<=0 || col>fNcolumns || layer<1) return 0;
 
 if (fType!=kTableHDU) return 0;
 
 TString str;
 Int_t dim=0;
 TObjString* sx;
 Int_t n=0;
 for (Int_t irow=rstart; irow<=rend; irow++)
 {
  dim=GetTableCell(str,irow,col,layer);
 
  if (!dim || layer>dim)
  {
   arr.Clear();
   return 0;
  }
 
  sx=new TObjString();
  sx->SetString(str.Data());
  arr.Add(sx);
  n++;
 }
 
 return n;
}

Int_t NcFITSIO::GetTableColumn(TObjArray &arr,TString colname,Int_t rstart,Int_t rend,Int_t layer,Int_t mode)
{
 
 Int_t col=GetColumnNumber(colname,mode);
 
 Int_t n=GetTableColumn(arr,col,rstart,rend,layer);
 
 return n;
}

Int_t NcFITSIO::GetImageDimension(Int_t i) const
{
 
 if (!fSizes) return 0;
 
 Int_t dim=fSizes->GetSize();
 
 if (!i) return dim;
 
 if (i<1 || i>dim) return 0;
 
 dim=fSizes->At(i-1);
 
 return dim;
}

Int_t NcFITSIO::GetImageLayer(TASImage &im,Int_t layer,Double_t* thres,Double_t max)
{
 
 // Set the image empty
 im.SetImage(0,0);
 
 TArrayD arr;
 if (!LoadLayer(arr,layer)) return 0;
 
 Int_t ndim1=fSizes->At(0);
 Int_t npix=arr.GetSize();
 
 if (npix)
 {
  if (thres) ApplyPixelThreshold(arr,*thres);
  if (max>0) RescalePixels(arr,max);
  im.SetImage(arr,ndim1);
 }
 
 return npix;
}

Int_t NcFITSIO::GetImageLayer(TMatrixD &m,Int_t layer,Double_t* thres,Double_t max)
{
 
 // Clear the matrix elements
 m.Clear();
 
 TArrayD arr;
 if (!LoadLayer(arr,layer)) return 0;
 
 Int_t ndim1=fSizes->At(0);
 Int_t ndim2=fSizes->At(1);
 Int_t npix=arr.GetSize();
 
 if (npix)
 {
  if (thres) ApplyPixelThreshold(arr,*thres);
  if (max>0) RescalePixels(arr,max);
  m.ResizeTo(ndim2,ndim1);
  Int_t i=0;
  for (Int_t jrow=0; jrow<ndim2; jrow++)
  {
   for (Int_t jcol=0; jcol<ndim1; jcol++)
   {
    m(jrow,jcol)=arr.At(i);
    i++;
   }
  }
 }
 
 return npix;
}

Int_t NcFITSIO::GetImageLayer(TH2D &his,Int_t layer,Double_t* thres,Double_t max)
{
 
 // Clear the histogram
 his.Reset();
 TString title="Histogram of layer ";
 title+=layer;
 his.SetTitle(title.Data());
 
 TMatrixD m;
 Int_t npix=GetImageLayer(m,layer,thres,max);
 
 if (!npix) return 0;
 
 Int_t nrows=m.GetNrows();
 Int_t ncols=m.GetNcols();
 
 // Set the binning and range of the histogram axes
 his.SetBins(ncols+1,0,ncols,nrows+1,0,nrows);
 
 // Fill the histogram cells with the corresponding matrix contents
 Double_t val=0;
 for (Int_t icol=0; icol<ncols; icol++)
 {
  for (Int_t irow=0; irow<nrows; irow++)
  {
   val=m(irow,icol);
   his.SetBinContent(icol+1,irow+1,val);
  }
 }
 
 return npix;
}

UInt_t NcFITSIO::GetImageArray(TArrayD &arr,TArrayI ifirst,TArrayI ilast,TArrayI incr)
{
 
 arr.Set(0);
 
 if (fType!=kImageHDU || !fSizes) return 0;
 
 Int_t ndims=fSizes->GetSize();
 
 if (ndims<1 || ifirst.GetSize()<ndims || ilast.GetSize()<ndims || incr.GetSize()<ndims) return 0;
 
 // Read the pixels of the specified (sub)set
 long* fpixel=new long[ndims];
 long* lpixel=new long[ndims];
 long* inc=new long[ndims];
 long long npixels=1;
 int status=0;
 
 Int_t istart=0;
 Int_t iend=0;
 Int_t istep=0;
 for (Int_t i=0; i<ndims; i++)
 {
  istart=ifirst.At(i);
  iend=ilast.At(i);
  istep=incr.At(i);
 
  if (istart<1 || iend<1 || istep<1 || iend<istart)
  {
   cout << endl;
   cout << " *" << ClassName() << "::GetImageArray* Inconsistent ifirst, ilast or incr input array(s)." << endl;
   fits_close_file(fInput,&status);
   Reset();
   delete[] fpixel;
   delete[] lpixel;
   delete[] inc;
   return 0;
  }
 
  fpixel[i]=istart;
  lpixel[i]=iend;
  inc[i]=istep;
  if (iend>istart) npixels*=1+(iend-istart)/istep;
 }
 
 if (!npixels) npixels=1;
 double* pixels=new double[npixels];
 
 double nulval=0;
 int anynul;
 status=0;
 fits_read_subset(fInput,TDOUBLE,fpixel,lpixel,inc,(void*)&nulval,(void*)pixels,&anynul,&status);
 
 if (status)
 {
  cout << endl;
  cout << " *" << ClassName() << "::GetImageArray* Could not read pixel data." << endl;
  fits_close_file(fInput,&status);
  Reset();
  delete[] fpixel;
  delete[] lpixel;
  delete[] inc;
  delete[] pixels;
  return 0;
 }
 
 // Store the pixel data in the output array
 arr.Set(npixels);
 Double_t val=0;
 for (long long i=0; i<npixels; i++)
 {
  val=pixels[i];
  arr.SetAt(val,i);
 }
 
 UInt_t npix=npixels;
 delete[] pixels;
 return npix;
}

UInt_t NcFITSIO::GetImageArray(TArrayD &arr,TArrayI ifirst,UInt_t npix)
{
 
 arr.Set(0);
 
 if (fType!=kImageHDU || !fSizes || npix<1) return 0;
 
 Int_t ndims=GetImageDimension();
 
 if (ndims<1 || ifirst.GetSize()<ndims) return 0;
 
 ULong_t nmax=1;
 for (Int_t i=1; i<=ndims; i++)
 {
  nmax*=GetImageDimension(i);
 } 
 
 if (npix>nmax) return 0;
 
 // Read the pixels of the specified (sub)set
 long* fpixel=new long[ndims];
 int status=0;
 
 Int_t istart=0;
 for (Int_t i=0; i<ndims; i++)
 {
  istart=ifirst.At(i);
 
  if (istart<1)
  {
   cout << endl;
   cout << " *" << ClassName() << "::GetImageArray* Inconsistent ifirst input array." << endl;
   fits_close_file(fInput,&status);
   Reset();
   delete[] fpixel;
   return 0;
  }
 
  fpixel[i]=istart;
 }
 
 long long npixels=npix;
 double* pixels=new double[npixels];
 double nulval=0;
 int anynul;
 status=0;
 fits_read_pix(fInput,TDOUBLE,fpixel,npixels,(void*)&nulval,(void*)pixels,&anynul,&status);
 
 if (status)
 {
  cout << endl;
  cout << " *" << ClassName() << "::GetImageArray* Could not read pixel data." << endl;
  fits_close_file(fInput,&status);
  Reset();
  delete[] fpixel;
  delete[] pixels;
  return 0;
 }
 
 // Store the pixel data in the output array
 arr.Set(npixels);
 Double_t val=0;
 for (long long i=0; i<npixels; i++)
 {
  val=pixels[i];
  arr.SetAt(val,i);
 }
 
 UInt_t nread=npixels;
 delete[] pixels;
 return nread;
}

void NcFITSIO::ListTable(Int_t width,Int_t rstart,Int_t rend,Int_t cstart,Int_t cend)
{
   
 if (fType != kTableHDU)
 {
  cout << " *" << ClassName() <<"::ListTable* This is not a table HDU." << endl;
  return;
 }
   
 if (!width)
 {
  cout << " *" << ClassName() <<"::ListTable* Invalid argument : width=" << width << endl;
  return;
 }
 
 if (!rend) rend=fNrows;
 
 if (!cend) cend=fNcolumns;
 
 if (rstart<=0 || rstart>fNrows || rend<=0 || rend>fNrows ||
     cstart<=0 || cstart>fNcolumns || cend<=0 || cend>fNcolumns)
 {
  cout << " *" << ClassName() <<"::ListTable* Invalid input rstart=" << rstart << " rend=" << rend
       << " cstart=" << cstart << " cend=" << cend << endl;
  return;
 }
 
 TString type;
 TString name;
 if (width<0) // List description of all table columns
 {  
  width=-width;
 
  cout << endl;
  cout << " *" << ClassName() << "::ListTable* Table column description for col=["
       << cstart << "," << cend << "] (name width is " << width << " characters)."<< endl;
  cout << endl;
 
  for (Int_t i=cstart; i<=cend; i++)
  {
   if (fColumnTypes[i-1]==kString) type="STRING";
   if (fColumnTypes[i-1]==kStringArray) type="STRING ARRAY";
   if (fColumnTypes[i-1]==kRealNumber) type="REAL NUMBER";
   if (fColumnTypes[i-1]==kRealArray) type="REAL ARRAY";
   if (fColumnTypes[i-1]==kComplexNumber) type="COMPLEX NUMBER";
   if (fColumnTypes[i-1]==kComplexArray) type="COMPLEX ARRAY";
   if (fColumnTypes[i-1]==kLogical) type="LOGICAL";
   if (fColumnTypes[i-1]==kLogicalArray) type="LOGICAL ARRAY";
 
   name=fColumnNames[i-1];
   name=name.Strip(name.kBoth,' ');
   name.Resize(width);
 
   cout << " " << name << " : " << type;
   if (type.Contains("ARRAY")) cout << "[" << fColumnLayers[i-1] << "]";
   cout << endl;
  }
  return;
 }
 
 if (width>0) // List the (selected) table contents
 {
 
  cout << endl;
  cout << " *" << ClassName() << "::ListTable* Table contents for row=[" << rstart << "," << rend << "] and"
       << " col=[" << cstart << "," << cend << "] (column width is " << width << " characters)."<< endl;
  cout << endl;
 
  // List the header with the (selected) column names
  if (width<7) width=7; // Minimal width needed for scientific notation
  Int_t prec=width-7;   // Adjust the precision
  TString str;
  Int_t nchars=0;
  for (Int_t col=cstart; col<=cend; col++)
  {
   str=fColumnNames[col-1];
   str=str.Strip(str.kBoth,' ');
   str.Resize(width);
   if (col==cstart) cout << " ";
   printf("%s| ",str.Data());
   nchars+=(width+2);
  }
  cout << endl;
  if (nchars>0) nchars--; // Don't count the space at the end
  cout << " ";
  while(nchars--)
  {
   cout << "-";
  }
  cout << endl;
   
  // List the (selected) row content
  Double_t val=0;
  Int_t ndim=0;
  for (Int_t row=rstart; row<=rend; row++)
  {
   for (Int_t col=cstart; col<=cend; col++)
   {
    if (col==cstart) cout << " ";
    if (fColumnTypes[col-1]==kString || fColumnTypes[col-1]==kStringArray || fColumnTypes[col-1]==kLogical || fColumnTypes[col-1]==kLogicalArray)
    {
     ndim=GetTableCell(str,row,col);
     if (!ndim) str="---";
     str=str.Strip(str.kBoth,' ');
     str.Resize(width);
     printf("%s",str.Data());
    }
    else // Numerical value
    {
     str="---";
     str.Resize(width);
     ndim=GetTableCell(val,row,col);
     if (ndim)
     {
      printf("%-*.*g",width,prec,val);
     }
     else
     {
      printf("%s",str.Data());
     }
    }
    cout << "| ";
   } // end of column loop
   cout << endl;
  } // end of row loop
 }
}

void NcFITSIO::ListHDUHeader() const
{
 
 cout << endl;
 cout << " *" << ClassName() << "::ListHDUHeader* The current HDU header record " << fExtensionName  << endl;
 cout << endl;
 
 for (Int_t i=0; i<fNkeys; i++)
 {
  printf(" %-8s = %-s",fKeyNames[i].Data(),fKeyValues[i].Data());
  if (fComments[i].Length()>0) printf(" / %-s",fComments[i].Data());
  printf("\n");
 }
}

void NcFITSIO::ListFileHeader(Int_t mode) const
{
 
 // Open a local mode of the FITS file without any filters and start at the primary HDU
 fitsfile* fp=0;
 int status=0;
 fits_open_file(&fp,fFilename.Data(),READONLY,&status);
 
 if (status)
 {
  cout << " *" << ClassName() << "::ListFileHeader* Could not open file : " << fFilename << endl;
  fits_close_file(fp,&status);
  return;
 }
 
 cout << endl;
 if (!mode)
 {
  cout << " *" << ClassName() << "::ListFileHeader* Short summary of the FITS file header information" << endl;
 }
 else
 {
  cout << " *" << ClassName() << "::ListFileHeader* Full FITS file header information" << endl;
 }
 cout << endl;
 
 // Obtain the number of HDUs present in the FITS file
 int nhdus=0;
 fits_get_num_hdus(fp,&nhdus,&status);
 
 if (status)
 {
  cout << " *" << ClassName() << "::ListFileHeader* Could not read the number of HDUs" << endl;
  fits_close_file(fp,&status);
  return;
 }
 
 cout << " Total number of HDUs : " << nhdus << endl;
 
 // Loop over all available HDUs
 int hdutype=0;
 int nkeys=0;
 int morekeys=0;
 char keyname[FLEN_KEYWORD+1];
 char keyvalue[FLEN_VALUE+1];
 char comment[FLEN_COMMENT+1];
 TString exttype;
 TString extname;
 TString* keynames=0;
 TString* keyvalues=0;
 TString* comments=0;
 for (Int_t jhdu=1; jhdu<=nhdus; jhdu++)
 {
  // Read the HDU type
  fits_get_hdu_type(fp,&hdutype,&status);
 
  if (status)
  {
   cout << " *" << ClassName() << "::ListFileHeader* Could not read the type of HDU ["<< (jhdu-1) << "]" << endl;
   fits_close_file(fp,&status);
   return;
  }
 
  exttype="unknown";
  if (hdutype==IMAGE_HDU) exttype="IMAGE";
  if (hdutype==ASCII_TBL) exttype="ASCII-TABLE";
  if (hdutype==BINARY_TBL) exttype="BINARY-TABLE";
 
  //Read the HDU header space information
  fits_get_hdrspace(fp,&nkeys,&morekeys,&status);
 
  if (status)
  {
   cout << " *" << ClassName() << "::ListFileHeader* Could not read the header space of HDU ["<< (jhdu-1) << "]" << endl;
   fits_close_file(fp,&status);
   return;
  }
 
  keynames=new TString[nkeys];
  keyvalues=new TString[nkeys];
  comments=new TString[nkeys];
 
  // Read the description of the HDU keys
  extname="";
  for (Int_t i=1; i<=nkeys; i++)
  {
   fits_read_keyn(fp,i,keyname,keyvalue,comment,&status);
 
   if (status)
   {
    cout << " *" << ClassName() << "::ListFileHeader* Could not read key number " << i << " of HDU ["<< (jhdu-1) << "]" << endl;
    fits_close_file(fp,&status);
    if (keynames) delete[] keynames;
    if (keyvalues) delete[] keyvalues;
    if (comments) delete[] comments;
    return;
   }
 
   keynames[i-1]=keyname;
   keyvalues[i-1]=keyvalue;
   comments[i-1]=comment;
 
   // Obtain the extension name if specified in the HDU
   if (keynames[i-1]=="EXTNAME")
   {
    extname=keyvalues[i-1];
    extname.ReplaceAll("'",""); // Remove all single quotes
    extname=extname.Strip(extname.kBoth,' '); // Remove leading and trailing blanks
    extname.Prepend("[");
    extname+="]";
   }
  }
 
  // List the HDU extension number, type and (optional) name
  cout << " [" << (jhdu-1) << "] " << exttype << " " << extname << endl;
 
  // List the contents of the HDU header record
  if (mode)
  {
   for (Int_t i=0; i<nkeys; i++)
   {
    printf(" %-8s = %-s",keynames[i].Data(),keyvalues[i].Data());
    if (comments[i].Length()>0) printf(" / %-s",comments[i].Data());
    printf("\n");
   }
   cout << endl;
  }
 
  if (keynames) delete[] keynames;
  if (keyvalues) delete[] keyvalues;
  if (comments) delete[] comments;
 
  // Move to the next HDU
  if (jhdu<nhdus)
  {
   fits_movrel_hdu(fp,1,&hdutype,&status);
 
   if (status)
   {
    cout << " *" << ClassName() << "::ListFileHeader* Could not move to HDU [" << jhdu << "]" << endl;
    fits_close_file(fp,&status);
    return;
   }
  }
 }
 
 // Close this loal "unfiltered" file mode
 fits_close_file(fp,&status);
 return;
}

Int_t NcFITSIO::LoadLayer(TArrayD &arr,Int_t layer)
{
 
 arr.Set(0);
 
 if (fType!=kImageHDU || !fSizes || layer<1) return 0;
 
 Int_t ndims=GetImageDimension();
 
 // Check whether the data dimensions are consistent with a (layered) image. 
 if (ndims<2) return 0;
 
 Int_t ndim1=GetImageDimension(1);
 Int_t ndim2=GetImageDimension(2);
 Int_t ndim3=GetImageDimension(3);
 
 // Check whether the layer number is within bounds
 if ((ndims==2 && layer>1) || (ndims>2 && layer>ndim3)) return 0;
 
 // Read the pixels of the specified layer
 TArrayI ifirst(ndims);
 TArrayI ilast(ndims);
 TArrayI incr(ndims);
 for (Int_t i=0; i<ndims; i++)
 {
  ifirst.SetAt(1,i);
  ilast.SetAt(1,i);
  incr.SetAt(1,i);
 }
 
 // Read a full layer content
 ilast.SetAt(ndim1,0);
 ilast.SetAt(ndim2,1);
 
 // The selected layer
 if (ndim3)
 {
  ifirst.SetAt(layer,2);
  ilast.SetAt(layer,2);
 }
 
 // Read the layer into the data storage
 // For a (nrow,ncol) matrix, the FITS parameters correspond to ncol=NAXIS1 (=width) and nrow=NAXIS2 (=height). 
 // The FITS image data start at the lower left corner and end at the right upper corner and follow a
 // pixel array storage with the following sequence convention for pixel coordinates (ix,iy) :
 // (1,1),(2,1),(3,1),...(NAXIS1,1),(1,2),(2,2),(3,2),...(NAXIS1,2),...(1,NAXIS2),(2,NAXIS2),...(NAXIS1,NAXIS2).
 // This implies that for a matrix interpretation the row numbering should be considered as inverted,
 // such that (row,col)=(1,1) indicates the lower left corner instead of the (usual) upper left corner.
 // With this convention, the upper left corner is (nrow,1), the lower left corner is (1,ncol)
 // and the upper right corner is (nrow,ncol).     
 GetImageArray(arr,ifirst,ilast,incr);
 
 Int_t npix=arr.GetSize();
 return npix;
}

void NcFITSIO::ApplyPixelThreshold(TArrayD &arr,Double_t thres)
{
 
 UInt_t npix=arr.GetSize();
 
 for (UInt_t i=0; i<npix; i++)
 {
  if (arr.At(i) < thres) arr.SetAt(0,i); 
 }
}

void NcFITSIO::RescalePixels(TArrayD &arr,Double_t max)
{
 
 UInt_t npix=arr.GetSize();
 
 if (!npix) return;
 
 Double_t minval=0;
 Double_t maxval=0;
 Double_t val=0;
 for (UInt_t i=0; i<npix; i++)
 {
  val=arr.At(i);
 
  if (!i)
  {
   minval=val;
   maxval=val;
   continue;
  }
 
  if (val<minval) minval=val; 
  if (val>maxval) maxval=val; 
 }
 
 Double_t range=maxval-minval;
 Double_t fact=-1;
 if (range>0) fact=max/range;
 
 for (UInt_t i=0; i<npix; i++)
 {
  if (fact<0)
  {
   arr.SetAt(max,i);
  }
  else
  {
   val=arr.At(i);
   val=fact*(val-minval);
   arr.SetAt(val,i);
  }
 }
}

TObject* NcFITSIO::Clone(const char* name) const
{
 
 NcFITSIO* q=new NcFITSIO(*this);
 if (name)
 {
  if (strlen(name)) q->SetName(name);
 }
 return q;
}