00001
00012 #include "FitsController.h"
00013
00014 #include "FitsFile.h"
00015 #include "FitsNTuple.h"
00016
00017 #include "controllers/CutController.h"
00018 #include "datasrcs/DataSourceController.h"
00019 #include "datasrcs/TupleCut.h"
00020 #include "pattern/string_convert.h"
00021 #include "plotters/PlotterBase.h"
00022
00023 #include <algorithm>
00024 #include <fstream>
00025
00026
00027 #include <cassert>
00028
00029 using std::ifstream;
00030 using std::runtime_error;
00031 using std::string;
00032 using std::vector;
00033
00034 using namespace hippodraw;
00035
00038 typedef std::map < std::string, FitsFile * >::iterator FileMapIterator;
00039
00040 FitsController * FitsController::s_instance = 0;
00041
00042 FitsController *
00043 FitsController::
00044 instance ()
00045 {
00046 if ( s_instance == 0 ) {
00047 s_instance = new FitsController ();
00048 }
00049 return s_instance;
00050 }
00051
00052 FitsController::
00053 ~FitsController()
00054 {
00055 }
00056
00057 const std::string &
00058 FitsController::
00059 version () const
00060 {
00061 float version = 0.0;
00062 fits_get_version ( & version );
00063
00064 m_version = String::convert ( version );
00065
00066 return m_version;
00067 }
00068
00069 FitsFile *
00070 FitsController::
00071 openFile ( const std::string & name )
00072 {
00073 FitsFile * file = 0;
00074 FileMapIterator first = m_file_map.find ( name );
00075
00076 if ( first == m_file_map.end() ) {
00077 ifstream test ( name.c_str (), std::ios::in );
00078 if ( test.is_open () == false ) {
00079 string what ( "FitsController: File `" );
00080 what += name;
00081 what += "' not found";
00082 throw std::runtime_error ( what );
00083 }
00084
00085 file = new FitsFile ( name );
00086 int status = file -> status();
00087
00088 if ( status != 0 ) {
00089 string what ( "cfitsio: Error opening file `" );
00090 what += name;
00091 what += "'";
00092 delete file;
00093 throw runtime_error ( what );
00094 }
00095 m_file_map [ name ] = file;
00096 }
00097 else {
00098 file = first -> second;
00099 }
00100
00101 return file;
00102 }
00103
00106 void
00107 FitsController::
00108 closeFile ( const std::string & )
00109 {
00110
00111
00112
00113
00114
00115
00116 }
00117
00118 const vector < string > &
00119 FitsController::
00120 getNTupleNames ( const std::string & file_name )
00121 {
00122 m_ntuple_names.clear ();
00123 FitsFile * file = openFile ( file_name );
00124 if ( file != 0 ) {
00125 file -> fillHDUNames ( m_ntuple_names );
00126 }
00127
00128 return m_ntuple_names;
00129 }
00130
00131
00132 DataSource *
00133 FitsController::
00134 createNTuple ( const std::string & filename,
00135 const std::string & name )
00136 {
00137 int index = -1;
00138
00139 const vector < string > & names
00140 = getNTupleNames ( filename );
00141 unsigned int size = names.size ();
00142
00143 for ( unsigned int i = 0; i < size; i++ ) {
00144 if ( name == names[i] ) {
00145 index = i;
00146 break;
00147 }
00148 }
00149
00150 if ( index == -1 ) {
00151 string what ( "FitsController: File `" );
00152 what += filename;
00153 what += "' ";
00154 what += "does not have extensison with name `";
00155 what += name;
00156 what += "'";
00157 throw runtime_error ( what );
00158 }
00159
00160 return createNTuple ( filename, name, index );
00161 }
00162
00163 DataSource *
00164 FitsController::
00165 createNTuple ( const std::string & filename,
00166 const std::string & name,
00167 int index )
00168 {
00169 DataSource * ntuple = 0;
00170
00171 string ds_name = filename;
00172 ds_name += ": ";
00173 ds_name += name;
00174
00175 FitsFile * file = openFile ( filename );
00176 if ( file == NULL ) return NULL;
00177
00178 int retval = file -> moveToHDU ( index );
00179 if ( retval == 0 ) {
00180 try {
00181 ntuple = new FitsNTuple ( file );
00182 }
00183 catch ( const runtime_error & e ) {
00184 throw e;
00185 }
00186
00187 ntuple -> setTitle ( name );
00188 ntuple -> setName ( ds_name );
00189 DataSourceController * controller = DataSourceController::instance();
00190 controller -> registerNTuple ( ds_name, ntuple );
00191 controller -> registerDataSourceFile ( ntuple );
00192 }
00193
00194 return ntuple;
00195 }
00196
00197 void
00198 FitsController::
00199 checkForImage ( PlotterBase * plotter, const DataSource & source )
00200 {
00201 const FitsFile * file = 0;
00202 try {
00203 const FitsNTuple & ntuple
00204 = dynamic_cast < const FitsNTuple & > ( source );
00205 file = ntuple.getFile ();
00206 }
00207 catch ( ... ) {
00208
00209 }
00210 if ( file == 0 ) return;
00211
00212 FitsFileBase::HduType type = file -> getHduType ();
00213 if ( type != FitsFileBase::Image ) return;
00214
00215 vector < long > sizes;
00216 file -> fillAxisSizes ( sizes );
00217 if ( sizes.size () > 3 ) {
00218 string what ( "FitsController: greater then 3D images not supported" );
00219 throw std::runtime_error ( what );
00220 }
00221
00222 plotter -> setNumberOfBins ( Axes::X, sizes[0] );
00223 plotter -> setNumberOfBins ( Axes::Y, sizes[1] );
00224
00225 vector < double > deltas;
00226 file -> fillImageDeltas ( deltas );
00227 plotter -> setBinWidth ( Axes::X, deltas[0] );
00228 plotter -> setBinWidth ( Axes::Y, deltas[1] );
00229
00230 vector < double > ref_values;
00231 file -> fillRefPixelValues ( ref_values );
00232 vector < int > ref_indices;
00233 file -> fillRefPixelIndices ( ref_indices );
00234
00235
00236 double pixToRef;
00237 if ((file ->pixCenter ())==true) pixToRef=0.5;
00238 else pixToRef=0.0;
00239
00240 double x_orig = - deltas[0] * ( ref_indices[0]-1+pixToRef ) + ref_values[0];
00241 double y_orig = - deltas[1] * ( ref_indices[1]-1+pixToRef ) + ref_values[1];
00242 plotter ->setOffset ( Axes::X, x_orig );
00243 plotter ->setOffset ( Axes::Y, y_orig );
00244
00245
00246
00247 Range range;
00248 if ( deltas[0] < 0. ) {
00249 range.setRange ( x_orig + sizes[0] * deltas[0], x_orig, -deltas[0] );
00250 }
00251 else {
00252 range.setRange ( x_orig, x_orig + (sizes[0] ) * deltas[0], 1. );
00253 }
00254 plotter -> setRange ( Axes::X, range, false, true );
00255 range.setLow ( y_orig );
00256 range.setLength ( ( sizes[1] ) * deltas[1] );
00257 plotter -> setRange ( Axes::Y, range, false, true );
00258
00259 plotter -> matrixTranspose ( true );
00260
00261 bool yes = file -> isHammerAitoff ();
00262 if ( yes ) {
00263 plotter ->setAspectRatio ( 2.0 );
00264 plotter ->setFitsTransform ("HammerAito");
00265 }
00266 }
00267
00270 void
00271 FitsController::
00272 writeNTupleToFile ( const DataSource * ntuple,
00273 const std::string & filename )
00274 {
00275 FitsFile * file = new FitsFile ( filename, true );
00276
00277 long row = static_cast<long> ( ntuple -> rows() );
00278 int column = static_cast<int> (ntuple -> columns());
00279 const vector <string> labels = ntuple -> getLabels();
00280 vector < vector < int > > shapes;
00281
00282 for ( int i = 0; i < column; i++ ){
00283 vector < int > shape;
00284 ntuple -> fillShape ( shape, i );
00285 shapes.push_back ( shape );
00286 }
00287 string extname;
00288 const string fullname = ntuple -> getName();
00289 string::size_type pos1 = fullname.find_last_of ( ":" );
00290 string::size_type pos2 = fullname.find_last_of ( "/" );
00291
00292 if ( pos1 != string::npos && pos1!=1 ) {
00293 extname = fullname.substr ( pos1+2 );
00294 }
00295 else if (pos2 != string::npos) {
00296 extname = "<no name>";
00297 }
00298 else {
00299 extname = fullname;
00300 }
00301
00302 file -> writeHDU ( row, column, labels, shapes, extname );
00303 for ( int i = 0; i < column; i++ ) {
00304 const vector < double > & col = ntuple -> getColumn (i );
00305 file -> writeColumn ( i, col );
00306 }
00307
00308 file -> writeCloseFile ();
00309
00310 delete file;
00311
00312 }
00313
00314 void
00315 FitsController::
00316 writeNTupleToFile ( const std::string & name, const std::string & filename )
00317 {
00318 DataSourceController * controller = DataSourceController::instance ();
00319 DataSource * ntuple
00320 = controller -> findDataSource ( name );
00321 if ( ntuple == 0 ) return;
00322
00323 writeNTupleToFile ( ntuple, filename );
00324 }
00325
00326 void
00327 FitsController::
00328 writeImageToFile ( unsigned int x, unsigned int y,
00329 const std::vector <double> & data,
00330 const std::string & filename )
00331 {
00332
00333 FitsFile * file = new FitsFile ( filename, true );
00334
00335 long xx = static_cast<long> ( x );
00336 long yy = static_cast<long> ( y );
00337
00338 double xlow = 0.0 - static_cast<double> (x) / 2 ;
00339 double ylow = 0.0 - static_cast<double> (y) / 2;
00340
00341 file -> writeImageHDU ( xx,yy );
00342 file -> writeRefPixelValues ( xlow, ylow );
00343 file -> writePix ( xx, yy, data );
00344
00345 file -> writeCloseFile ();
00346
00347 delete file;
00348
00349 }
00350
00351 size_t
00352 FitsController::
00353 calcColumnWidth ( const DataSource * source, unsigned int column ) const
00354 {
00355 vector < int > shape;
00356 source -> fillShape ( shape, column );
00357 size_t rank = shape.size ();
00358 size_t count = 1;
00359
00360 for ( size_t i = 1; i < rank; i++ ) {
00361 count *= shape[i];
00362 }
00363 return count;
00364 }
00365
00366 int
00367 FitsController::
00368 writeNTupleToFile ( const DataSource * ds,
00369 const std::string & filename,
00370 const std::string & dsname,
00371 const std::vector < std::string > & column_list,
00372 const std::vector < const TupleCut * > & cut_list)
00373 {
00374 if ( column_list.empty() ) return 1;
00375
00376 FitsFile * file = new FitsFile ( filename, true );
00377
00378 unsigned int columnNumber = column_list.size();
00379 unsigned int size = ds->rows();
00380
00381 vector < bool > acceptArray;
00382 CutController::fillAcceptedRows ( acceptArray, ds, cut_list );
00383
00384 long rows = std::count ( acceptArray.begin(), acceptArray.end(), true );
00385 int columns = static_cast<int> (columnNumber);
00386
00387
00388
00389 vector < int > col_indices ( columnNumber );
00390
00391 for ( unsigned int i = 0; i < columnNumber; i++ ) {
00392 const string & label = column_list [ i ];
00393 int index = ds -> indexOf ( label );
00394 if ( index < 0 ) {
00395 ds -> throwIfInvalidLabel ( label );
00396 }
00397 col_indices [i] = index;
00398 }
00399
00400 vector < size_t > sizes;
00401 vector < vector < int > > shapes;
00402 for ( int i = 0; i < columns; i++ ) {
00403 unsigned int size = calcColumnWidth ( ds, col_indices[i] );
00404 sizes.push_back ( size );
00405 vector < int > shape;
00406 ds -> fillShape ( shape, col_indices[i] );
00407 shapes.push_back ( shape );
00408 }
00409
00410 file -> writeHDU ( rows, columns, column_list,
00411 shapes, dsname );
00412
00413
00414 for ( unsigned int k = 0; k < columnNumber; k++ ) {
00415 std::vector < double > tempColumn;
00416 if ( sizes [ k ] == 1 ) {
00417 tempColumn.reserve (size );
00418 }
00419 else {
00420 tempColumn.reserve ( size * sizes[k] );
00421 }
00422 for ( unsigned int i = 0; i<size; i++ ) {
00423 if ( acceptArray[i]==true ) {
00424 int index = col_indices[k];
00425 if ( sizes [ k ] == 1 ) {
00426 tempColumn.push_back ( ds -> valueAtNoCache (i, index) );
00427 }
00428 else {
00429 double * array = ds -> doubleArrayAt ( i, index );
00430 for ( std::size_t l = 0; l < sizes [ k ]; l++ ) {
00431 tempColumn.push_back ( array[l] );
00432 }
00433 }
00434 }
00435 }
00436 file -> writeColumn ( k, tempColumn );
00437
00438 tempColumn.clear();
00439 }
00440
00441 file -> writeCloseFile ();
00442
00443 delete file;
00444
00445 return 0;
00446 }