PeaksWorkspace.cpp 31.4 KB
Newer Older
Ruth Mikkelson's avatar
Ruth Mikkelson committed
1
2
3
4
5
#include "MantidKernel/Logger.h"
#include "MantidKernel/Unit.h"
#include "MantidKernel/DateAndTime.h"
#include "MantidAPI/ColumnFactory.h"
#include "MantidGeometry/Quat.h"
6
#include "MantidDataObjects/TableWorkspace.h"
Ruth Mikkelson's avatar
Ruth Mikkelson committed
7
8
9
10
11
12
13
14
15
16
#include "MantidAPI/WorkspaceProperty.h"
#include "MantidAPI/WorkspaceFactory.h"
#include "MantidDataObjects/PeaksWorkspace.h"
#include "MantidGeometry/V3D.h"
#include "MantidAPI/AlgorithmFactory.h"
#include "MantidAPI/MatrixWorkspace.h"
#include "boost/shared_ptr.hpp"
#include <algorithm>
#include <string>
#include <ostream>
17
#include <exception>
Ruth Mikkelson's avatar
Ruth Mikkelson committed
18
19
20
21
22
23
24
25
#include <iostream>
#include <fstream>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

namespace Mantid
{
26
27
28
29
namespace DataObjects
{
  /// Register the workspace as a type
  DECLARE_WORKSPACE(PeaksWorkspace );
Ruth Mikkelson's avatar
Ruth Mikkelson committed
30

31
  double pi = 3.1415926535;
32
/*//moved to header file
33
  std::vector< std::string > DetNames;
Ruth Mikkelson's avatar
Ruth Mikkelson committed
34

35
  std::vector< double* > DetInfo;
Ruth Mikkelson's avatar
Ruth Mikkelson committed
36

37
  Kernel::DateAndTime  C_experimentDate;
38
  */
39
  Kernel::Logger& PeaksWorkspace::g_log = Kernel::Logger::get("PeaksWorkspace");
Ruth Mikkelson's avatar
Ruth Mikkelson committed
40
41


42
43
44
45
46
47
  /** Constructor. Create a table with all the required columns.
   *
   * @return PeaksWorkspace object
   */
  PeaksWorkspace::PeaksWorkspace():TableWorkspace( 0 )
  {
Ruth Mikkelson's avatar
Ruth Mikkelson committed
48

49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
    TableWorkspace::addColumn( std::string( "V3D" )  ,  std::string( "hkl" ) );
    TableWorkspace::addColumn( std::string( "V3D" ) ,  std::string( "position" ) );
    TableWorkspace::addColumn( std::string( "V3D" ) ,  std::string( "sample_orientation" ) );
    TableWorkspace::addColumn( std::string( "int" ) ,  std::string( "reflag" ) );
    TableWorkspace::addColumn( std::string( "int" ) ,  std::string( "bank" ) );
    TableWorkspace::addColumn( std::string( "int" ) ,  std::string( "run" ) );
    TableWorkspace::addColumn( std::string( "double" ) ,  std::string( "ipeak" ) );
    TableWorkspace::addColumn( std::string( "double" ) ,  std::string( "inti" ) );
    TableWorkspace::addColumn( std::string( "double" ) ,  std::string( "sigi" ) );
    TableWorkspace::addColumn( std::string( "double" ) ,  std::string( "moncount" ) );
    TableWorkspace::addColumn( std::string( "double" ) ,  std::string( "row" ) );
    TableWorkspace::addColumn( std::string( "double" ) ,  std::string( "col" ) );
    TableWorkspace::addColumn( std::string( "double" ) ,  std::string( "chan" ) );
    TableWorkspace::addColumn( std::string( "double" ) ,  std::string( "L1" ) );
    TableWorkspace::addColumn( std::string( "double" ) ,  std::string( "L2" ) );
    TableWorkspace::addColumn( std::string( "double" ) ,  std::string( "time" ) );
    TableWorkspace::addColumn( std::string( "double" ) ,  std::string( "t_offset" ) );

  }


  /** Destructor */
  PeaksWorkspace::~PeaksWorkspace()
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
  {
    ClearDeleteCalibrationData();
  }

  void PeaksWorkspace:: ClearDeleteCalibrationData()
  {
    DetNames.clear();// elements std::string so will delete self??

    for( size_t n=DetInfo.size(); n >0;)
    {
      double* val = DetInfo[n-1];
      DetInfo.pop_back();
      delete val;

      n=DetInfo.size();
    }
  }
89
90

  void PeaksWorkspace::initialize( double      L1 ,               //m
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
                                      double      time_offset ,      //microseconds
                                      std::string Facility ,
                                      std::string Instrument ,
                                      Kernel::DateAndTime  experimentDate ,
                                      std::string version ,

                                      std::vector<std::string> &PanelNames,
                                      std::vector< double*> &PanelInfo
                                      )
     {
       C_L1 = L1;
       C_time_offset = time_offset;
       C_Facility = Facility;
       C_Instrument = Instrument;
       C_version = version;
       C_experimentDate = experimentDate;
107
       ClearDeleteCalibrationData();
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
       for( int i=0; i<PanelNames.size() ; i++)
         DetNames.push_back( std::string(PanelNames[i]));

       for( int i=0; i<PanelInfo.size(); i++)
       {
         double* data = PanelInfo[i];
         double* data1 = new double[15];
         for( int j=0; j<15;j++)
         {
           data1[j] = data[j];
           if( (j>=2) &&( j <=8))
             data1[j] *=100;
         }
         DetInfo.push_back(data1);
       }
     }

     void PeaksWorkspace::initialize(  std::string DetCalFileName)
     {
127
       ClearDeleteCalibrationData();
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143

       try
       {
         std::ifstream in( DetCalFileName.c_str() );

         std::string s = readHeader( in );
         in.close();
       }catch( char * str)
       {
         std::cout <<"Exception reading detector info "<<str <<std::endl;
       }

     }

  void PeaksWorkspace::addPeak( const Geometry::V3D position ,
      const double time ,
144
145
146
147
148
149
150
151
152
153
      const Geometry::V3D hkl ,
      const Geometry::V3D sample_orientation ,  //radians ,  phi,chi,omega
      const int  reflag ,
      const int  runNum ,
      const double monCount ,
      const int bankName ,
      const double PeakCellCount ,
      const double row ,
      const double col ,
      const double chan ,
154
      const double L2 ,
155
156
157
158
159
160
161
      const double PeakIntegrateCount ,
      const double PeakIntegrateError
  )
  {
    int i = rowCount();
    insertRow( i );
    double T = time;
Ruth Mikkelson's avatar
Ruth Mikkelson committed
162

163
    try
Ruth Mikkelson's avatar
Ruth Mikkelson committed
164
165
    {

166
167
168
169
      getRef< double >( std::string( "ipeak" ) ,  i ) = PeakCellCount;
      getRef< double >( std::string( "inti" ) ,  i ) = PeakIntegrateCount;
      getRef< double >( std::string( "sigi" ) ,  i ) = PeakIntegrateError;
      getRef< double >( std::string( "moncount" ) ,  i ) = monCount;
Ruth Mikkelson's avatar
Ruth Mikkelson committed
170

171
172
      getRef< Geometry::V3D >( std::string( "position" ) ,  i ) =
          Geometry::V3D( position );
Ruth Mikkelson's avatar
Ruth Mikkelson committed
173

174
175
      getRef< Geometry::V3D >( std::string( "sample_orientation" ) ,  i ) =
          Geometry::V3D( sample_orientation );
Ruth Mikkelson's avatar
Ruth Mikkelson committed
176

177
178
      getRef< Geometry::V3D >( std::string( "hkl" ) ,  i ) =
          Geometry::V3D( hkl );
Ruth Mikkelson's avatar
Ruth Mikkelson committed
179

180
181
182
183
      getRef< double >( std::string( "L1" ) ,  i ) = C_L1;
      getRef< double >( std::string( "time" ) ,  i ) = time;
      getRef< double >( std::string( "t_offset" ) ,  i ) = C_time_offset;
      getRef< int >( std::string( "run" ) ,  i ) = runNum;
Ruth Mikkelson's avatar
Ruth Mikkelson committed
184

185
      getRef< double >( std::string( "chan" ) ,  i ) = chan;
186
      getRef< double >( std::string( "L2" ) ,  i ) = L2;
187
188
      getRef< int >( std::string( "bank" ) ,  i ) =
          ( bankName );
Ruth Mikkelson's avatar
Ruth Mikkelson committed
189

190
191
192
      getRef< double >( std::string( "row" ) ,  i ) = row;
      getRef< double >( std::string( "col" ) ,  i ) = col;
      getRef< int >( std::string( "reflag" ) ,  i ) = reflag;
Ruth Mikkelson's avatar
Ruth Mikkelson committed
193

194
195
196
197
    }catch( char * str )
    {
      std::cout <<  "exception occurred " <<  str <<  "\n";
    }
Ruth Mikkelson's avatar
Ruth Mikkelson committed
198

199
  }
Ruth Mikkelson's avatar
Ruth Mikkelson committed
200

201

202
  /**
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
   * @return the number of peaks
   */
  int PeaksWorkspace::getNumberPeaks() const
  {
    return rowCount();
  }

  /** Remove a peak
   *
   * @param peakNum :: peak index to remove
   */
  void PeaksWorkspace::removePeak( int peakNum )
  {
    removeRow( peakNum );
  }
Ruth Mikkelson's avatar
Ruth Mikkelson committed
218
219


220
221
222
223
  //returns true if strictly less than
  bool compareAsVecs( std::vector< int > arg1, std::vector< int > arg2 )
  {
    if( arg1.size() < 3 )
Ruth Mikkelson's avatar
Ruth Mikkelson committed
224
    {
225
226
227
228
      if( arg2.size() < 3 )
        return true;
      else
        return false;
Ruth Mikkelson's avatar
Ruth Mikkelson committed
229
230
    }

231
232
    if( arg2.size() < 3 )
      return false;   //hopefully bad values go to back
Ruth Mikkelson's avatar
Ruth Mikkelson committed
233

234
235
    int r1 = arg1[0];
    int r2 = arg2[0];
Ruth Mikkelson's avatar
Ruth Mikkelson committed
236

237
238
    if( r1 < r2 )
      return true;
Ruth Mikkelson's avatar
Ruth Mikkelson committed
239

240
241
    if( r1 > r2 )
      return false;
Ruth Mikkelson's avatar
Ruth Mikkelson committed
242

243
244
245
246
247
248
249
250
    r1 = arg1[1];
    r2 = arg2[1];

    if( r1 < r2 )
      return true;

    return false;
  }
Ruth Mikkelson's avatar
Ruth Mikkelson committed
251
252
253



254
255
256
257
258
259
260
261
262
  /** Write an ISAW-style .peaks file
   *
   * @param filename :: filename to write to
   */
  void PeaksWorkspace::write( std::string filename )
  {
    try
    {
      std::ofstream out( filename.c_str() );
263
264
      std::string date =C_experimentDate.to_ISO8601_string();

265
      out <<  "Version: " <<  C_version <<  " Facility: " <<  C_Facility  ;
266
267
      out <<  " Instrument: " <<  C_Instrument <<  " Date: " ;
      out<< std::setw(date.length())<<date<< std::endl;
Ruth Mikkelson's avatar
Ruth Mikkelson committed
268

269
      out <<  "6        L1    T0_SHIFT" <<  std::endl;
270
      out << "7 "<< std::setw( 11 )  ;
271
272
273
      out <<   std::setprecision( 4 ) <<  std::fixed <<  ( C_L1*100 ) ;
      out << std::setw( 12 ) <<  std::setprecision( 3 ) <<  std::fixed  ;
      out << ( C_time_offset ) <<  std::endl;
Ruth Mikkelson's avatar
Ruth Mikkelson committed
274

275
276
277
      out <<  "4 DETNUM  NROWS  NCOLS   WIDTH   HEIGHT     DEPTH   DETD  "
          <<  " CenterX   CenterY   CenterZ   BaseX    BaseY    BaseZ    "
          <<  "  UpX      UpY      UpZ" <<  std::endl;
Ruth Mikkelson's avatar
Ruth Mikkelson committed
278

279
280
281
282
      for(size_t i = 0 ; i < DetNames.size() ; i++ )
      {
        double* info = DetInfo[ i ];
        out <<  "5" <<  std::setw( 7 ) <<  DetNames[ i ];
Ruth Mikkelson's avatar
Ruth Mikkelson committed
283

284
285
286
        out <<  std::setw( 7 ) <<  (int)info[ 0 ] <<  std::setw( 7 );
        out  <<  (int)info[ 1 ]<<  std::setw( 9 ) <<  std::setprecision( 4 );
        out  <<  std::fixed <<  info[ 2 ];
Ruth Mikkelson's avatar
Ruth Mikkelson committed
287

288
289
290
291
        for( int j = 3 ; j < 15 ; j++ )
        {
          out <<  std::setw( 9 ) <<  std::setprecision( 4 ) <<  std::fixed;
          out <<  info[ j ];
Ruth Mikkelson's avatar
Ruth Mikkelson committed
292

293
294
        }
        out  <<  std::endl;
Ruth Mikkelson's avatar
Ruth Mikkelson committed
295

296
      }
Ruth Mikkelson's avatar
Ruth Mikkelson committed
297
298


299
      std::vector< std::vector< int > > SortKeys;
Ruth Mikkelson's avatar
Ruth Mikkelson committed
300

301
302
303
      for( int i = 0 ; i < rowCount() ; i++ )
      {
        std::vector<int > dat ;
Ruth Mikkelson's avatar
Ruth Mikkelson committed
304

305
306
        dat.push_back(getRef< int >( "run" ,  i ) );
        dat.push_back( getRef< int >( "bank" ,  i ) );
Ruth Mikkelson's avatar
Ruth Mikkelson committed
307

308
309
310
        dat.push_back( i );
        SortKeys.push_back( dat );
      }
Ruth Mikkelson's avatar
Ruth Mikkelson committed
311

312
      stable_sort( SortKeys.begin() ,  SortKeys.end() ,  compareAsVecs );
Ruth Mikkelson's avatar
Ruth Mikkelson committed
313

314
315
      int LastRun = -1;;
      int LastDet;
Ruth Mikkelson's avatar
Ruth Mikkelson committed
316

317
318
319
320
321
      int thisRun ,  thisReflag;
      double this_h ,  this_k ,  this_l ,   this_row ,  this_col ,
      this_chan ,  this_L2 ,  this_az ,  this_polar ,  this_time ,
      this_inti ,  this_ipk ,   this_sigi ,  this_chi ,  this_phi ,
      this_omega ,   this_monCt ,  this_tOffset ,  this_L1;
Ruth Mikkelson's avatar
Ruth Mikkelson committed
322

323
324
325
326
      int thisDet;
      Geometry::V3D thisSampleOrientation;
      int seqNum = 0;
      std::vector< std::vector< int > >::iterator it;
Ruth Mikkelson's avatar
Ruth Mikkelson committed
327

328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
      for ( it = SortKeys.begin(); it != SortKeys.end(); ++it )
      {
        seqNum++;
        std::vector<int > elt = *it;
        int peakNum =elt[2];

        this_row =
            getRef< double  >( std::string( "row" ) ,  peakNum );
        Geometry::V3D hkl = getRef< Geometry::V3D  > ( std::string( "hkl" ) ,
            peakNum );
        this_h = hkl.X();
        this_k = hkl.Y();
        this_l = hkl.Z();
        this_col = getRef < double >( std::string( "col" ) ,  peakNum );
        this_chan = getRef< double >( std::string( "chan" ) ,  peakNum );
        this_ipk = getRef< double >( "ipeak" ,  peakNum );
        this_inti = getRef< double >( "inti" ,  peakNum );
        this_sigi = getRef< double >( "sigi" ,  peakNum );

        this_monCt = getRef< double >( std::string( "moncount" ) ,  peakNum );
        this_L1 = getRef< double >( "L1" ,  peakNum );
        this_tOffset = getRef< double >( "t_offset" ,  peakNum );
        this_time = getRef< double >( "time" ,  peakNum );
        thisRun = getRef< int >( "run" ,  peakNum );
        thisReflag = Int( peakNum ,  3 );
        thisDet = getRef< int >( "bank" ,  peakNum );

        Geometry::V3D position = getRef< Geometry::V3D >( "position" ,
            peakNum );

        double x = position.X();
        double y = position.Y();
        double z = position.Z();
        this_L2 = sqrt( x*x + y*y + z*z );
        this_polar = acos(  z/this_L2 );
        this_az = atan2( y ,  x );


        thisSampleOrientation = getRef<Geometry::V3D >( "sample_orientation" ,
            peakNum );
        double mult = 180/3.1415926535897932384626433832795;
        this_chi = thisSampleOrientation.Y()*mult;
        this_phi = thisSampleOrientation.X()*mult;
        this_omega = thisSampleOrientation.Z()*mult;

        if( thisRun != LastRun  || LastDet != thisDet )
        {
          out <<  "0 NRUN DETNUM    CHI    PHI  OMEGA MONCNT" <<  std::endl;
          out <<  "1" <<  std::setw( 5 ) <<  thisRun <<  std::setw( 7 ) <<
              std::right <<  thisDet;

          out  <<  std::setw( 7 ) <<  std::fixed <<  std::setprecision( 2 )
          <<  this_chi;

          out  <<  std::setw( 7 ) <<  std::fixed <<  std::setprecision( 2 )
          <<  this_phi;

          out  <<  std::setw( 7 ) <<  std::fixed <<  std::setprecision( 2 )
          <<  this_omega;
          out  <<  std::setw( 7 ) <<  (int)( .5 + this_monCt ) <<  std::endl;
          LastRun = thisRun;
          LastDet = thisDet;
          out <<  "2   SEQN    H    K    L     COL     ROW    CHAN       L2  "
              <<  "2_THETA       AZ        WL        D   IPK      INTI   "
              << "SIGI RFLG" <<  std::endl;
        }
394

395
396
397
        int h = (int)floor( this_h+.5) , //assume when set via UB those
            k = (int)floor( this_k +.5 ) , //not close enuf get 0's
            l = (int)floor( this_l +.5 );
Ruth Mikkelson's avatar
Ruth Mikkelson committed
398

399

400
401
402
403
        out <<  "3" <<  std::setw( 7 ) <<  seqNum <<  std::setw( 5 ) <<  h
            <<  std::setw( 5 ) <<  k <<  std::setw( 5 ) <<  l;
        out <<  std::setw( 8 ) <<  std::fixed << std::setprecision( 2 )
        << this_col;
Ruth Mikkelson's avatar
Ruth Mikkelson committed
404

405
406
        out << std::setw( 8 ) << std::fixed << std::setprecision( 2 )
        << this_row;
Ruth Mikkelson's avatar
Ruth Mikkelson committed
407

408
409
        out << std::setw( 8 ) << std::fixed << std::setprecision( 2 )
        << this_chan;
Ruth Mikkelson's avatar
Ruth Mikkelson committed
410
411


412
413
        out << std::setw( 9 ) << std::fixed << std::setprecision( 3 )
        << ( this_L2*100 );
Ruth Mikkelson's avatar
Ruth Mikkelson committed
414

415
416
        out << std::setw( 9 ) << std::fixed << std::setprecision( 5 )
        << this_polar;
Ruth Mikkelson's avatar
Ruth Mikkelson committed
417

418
419
        out << std::setw( 9 ) << std::fixed << std::setprecision( 5 )
        << this_az;
Ruth Mikkelson's avatar
Ruth Mikkelson committed
420
421


422
423
424
425
426
        std::vector< double >xx ,  yy;
        xx.push_back( this_time );
        Kernel::Units::Wavelength wl;
        wl.fromTOF( xx ,  yy ,  this_L1 ,  this_L2 ,  this_polar ,
            0 ,  1.2 ,  1.2 );
Ruth Mikkelson's avatar
Ruth Mikkelson committed
427

428
429
        out << std::setw( 10 ) << std::fixed << std::setprecision( 6 )
        << xx[ 0 ];
Ruth Mikkelson's avatar
Ruth Mikkelson committed
430

431
432
        xx[ 0 ] = this_time;
        Kernel::Units::dSpacing dsp;
433

434
435
        dsp.fromTOF( xx ,  yy ,  this_L1 ,  this_L2 ,  this_polar ,
            0 ,  1.2 ,  1.2 );
Ruth Mikkelson's avatar
Ruth Mikkelson committed
436

437
438
        out << std::setw( 9 ) << std::fixed << std::setprecision( 4 )
        <<  xx[ 0 ];
Ruth Mikkelson's avatar
Ruth Mikkelson committed
439

440
441
        out << std::setw( 6 ) << (int)this_ipk << std::setw( 10 )
        << std::fixed << std::setprecision( 2 ) << this_inti;
Ruth Mikkelson's avatar
Ruth Mikkelson committed
442

443
444
        out << std::setw( 7 ) << std::fixed << std::setprecision( 2 )
        << this_sigi << std::setw( 5 ) << thisReflag << std::endl;
Ruth Mikkelson's avatar
Ruth Mikkelson committed
445

446
      }
Ruth Mikkelson's avatar
Ruth Mikkelson committed
447

448
449
      out.flush();
      out.close();
Ruth Mikkelson's avatar
Ruth Mikkelson committed
450

451
    }catch( char *s )
Ruth Mikkelson's avatar
Ruth Mikkelson committed
452
    {
453
454
      g_log.error(std::string(s));
      throw( std::string(s));
Ruth Mikkelson's avatar
Ruth Mikkelson committed
455
    }
456
457
   catch( std::exception &e)
   {
458
459
     g_log.error(e.what());
          throw( e.what());
460
461
   }catch(...)
   {
462
463
     g_log.error("Unknown Error");
     throw( "Unknown Error");
464
   }
465
  }
Ruth Mikkelson's avatar
Ruth Mikkelson committed
466
467


468
469
470
471
472
473
  /** Remove all the peaks from the workspace */
  void PeaksWorkspace::removeAllPeaks()
  {
    for( int i = rowCount() - 1 ; i >= 0 ; i-- )
      removeRow( i );
  }
Ruth Mikkelson's avatar
Ruth Mikkelson committed
474
475


476
477
478
479
480
481
482
483
484
485
  //stops at spaces and \n. strips
  std::string getWord( std::ifstream &in ,  bool consumeEOL )
  {
    std::string s;
    char c = 0;
    if( in.good() )
      for(  c = in.get() ; c == ' ' && in.good() ; c = in.get() )
      {}
    else
      return std::string();
Ruth Mikkelson's avatar
Ruth Mikkelson committed
486

487
488
489
    if( c == '\n' )
    {
      if( !consumeEOL )
Ruth Mikkelson's avatar
Ruth Mikkelson committed
490
491
        in.putback( c );

492
      return std::string();
Ruth Mikkelson's avatar
Ruth Mikkelson committed
493
494
    }

495
    s.push_back( c );
Ruth Mikkelson's avatar
Ruth Mikkelson committed
496

497
    if( in.good() )
498
      for(  c = in.get() ; in.good()&&c != ' ' && c != '\n' && c != '\r' ; c = in.get() )
499
        s.push_back( c );
Ruth Mikkelson's avatar
Ruth Mikkelson committed
500

501
    if( ((c == '\n') || (c == '\r')) && !consumeEOL )
502
      in.putback( c );
Ruth Mikkelson's avatar
Ruth Mikkelson committed
503

504
505
    return s;
  }
Ruth Mikkelson's avatar
Ruth Mikkelson committed
506
507


508
509
510
511
  void readToEndOfLine( std::ifstream& in ,  bool ConsumeEOL )
  {
    while( in.good() && getWord( in ,  false ).length() > 0  )
      getWord( in ,  false );
Ruth Mikkelson's avatar
Ruth Mikkelson committed
512

513
514
    if( !ConsumeEOL )
      return ;
Ruth Mikkelson's avatar
Ruth Mikkelson committed
515

516
517
    getWord( in ,  true );
  }
Ruth Mikkelson's avatar
Ruth Mikkelson committed
518
519


520
521
522
523
524
525
  /** Reads the header of a .peaks file
   *
   * @param in :: stream of the input file
   * @return TODO: I don't know what here
   */
  std::string PeaksWorkspace::readHeader( std::ifstream& in )
526
  {
Ruth Mikkelson's avatar
Ruth Mikkelson committed
527

528
    std::string r = getWord( in ,  false );
529

530
531
    if( r.length() < 1 )
      throw std::logic_error( std::string( "No first line of Peaks file" ) );
Ruth Mikkelson's avatar
Ruth Mikkelson committed
532

533
534
535
    if( r.compare( std::string( "Version:" ) ) != 0 )
      throw std::logic_error(
          std::string( "No Version: on first line of Peaks file" ) );
Ruth Mikkelson's avatar
Ruth Mikkelson committed
536

537
    C_version = getWord( in ,  false );
Ruth Mikkelson's avatar
Ruth Mikkelson committed
538

539
540
    if( C_version.length() < 1 )
      throw  std::logic_error( std::string( "No Version for Peaks file" ) );
Ruth Mikkelson's avatar
Ruth Mikkelson committed
541

542
    C_Facility = getWord( in ,  false );
Ruth Mikkelson's avatar
Ruth Mikkelson committed
543

544
545
546
    if( C_Facility.length() < 1  )
      throw  std::logic_error(
          std::string( "No Facility tag for Peaks file" ) );
Ruth Mikkelson's avatar
Ruth Mikkelson committed
547

548
549
550
    if( C_Facility.compare( std::string( "Facility:" ) ) != 0 )
      throw  std::logic_error(
          std::string( "No Facility tag for Peaks file" ) );
Ruth Mikkelson's avatar
Ruth Mikkelson committed
551

552
    C_Facility = getWord( in ,  false );
Ruth Mikkelson's avatar
Ruth Mikkelson committed
553

554
555
    if( C_Facility.length() < 1 )
      throw  std::logic_error( std::string( "No Facility in Peaks file" ) );
Ruth Mikkelson's avatar
Ruth Mikkelson committed
556

557
    C_Instrument = getWord( in ,  false );
Ruth Mikkelson's avatar
Ruth Mikkelson committed
558

559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
    if( C_Instrument.length() < 1 )
      throw  std::logic_error(
          std::string( "No Instrument tag for Peaks file" ) );

    if( C_Instrument.compare( std::string( "Instrument:" ) ) != 0 )
      throw  std::logic_error(
          std::string( "No Instrument tag for Peaks file" ) );

    C_Instrument = getWord( in ,  false );

    if( C_Instrument.length() < 1 )
      throw std::logic_error(
          std::string( "No Instrument for Peaks file" ) );

    std::string date = getWord( in ,  false );
Ruth Mikkelson's avatar
Ruth Mikkelson committed
574

575
576
577
578
579
580
581
582
    if(date .length()< 1)
    {
      Kernel::DateAndTime x = Kernel::DateAndTime::get_current_time();
      C_experimentDate.set_from_time_t(x.to_time_t());
    }

    else
      if( date.compare( std::string( "Date:" ) ) == 0 )
Ruth Mikkelson's avatar
Ruth Mikkelson committed
583
584
585
586
587
588
      {
        date = getWord( in ,  false );
        C_experimentDate = Kernel::DateAndTime( date );
      }


589
590
    while( getWord( in ,  false ).length() > 0 )
      getWord( in ,  false );
Ruth Mikkelson's avatar
Ruth Mikkelson committed
591

592
    // Now get info on detector banks
Ruth Mikkelson's avatar
Ruth Mikkelson committed
593

594
595
596
597
598
    getWord( in ,  true );
    std::string s = getWord( in ,  false );
    if( s.compare( "6" ) != 0 )
      throw std::logic_error(
          std::string( "Peaks File has no L0 and T0 header info" ) );
Ruth Mikkelson's avatar
Ruth Mikkelson committed
599

600
601
    readToEndOfLine( in ,  true );
    s = getWord( in ,  false );
Ruth Mikkelson's avatar
Ruth Mikkelson committed
602

603
604
605
    if( s.compare( "7" ) != 0 )
      throw std::logic_error(
          std::string( "Peaks File has no L0 and T0 info" ) );
Ruth Mikkelson's avatar
Ruth Mikkelson committed
606

607
608
    C_L1 = strtod( getWord( in ,  false ).c_str() ,  0 )/100;
    C_time_offset = strtod( getWord( in ,  false ).c_str() ,  0 );
Ruth Mikkelson's avatar
Ruth Mikkelson committed
609
610


611
612
613
614
615
    readToEndOfLine( in ,  true );
    for( s = getWord( in ,  false ) ; s.length() < 1 ; getWord( in ,  true ) )
    {
      s = getWord( in ,  false );
    }
Ruth Mikkelson's avatar
Ruth Mikkelson committed
616

617
618
    if( s.compare( "4" ) != 0 )
      throw std::logic_error( std::string( "Peaks File has no bank descriptor line" ) );
Ruth Mikkelson's avatar
Ruth Mikkelson committed
619

620
621
622
623
624
    readToEndOfLine( in ,  true );
    s = getWord( in ,  false );
    while( s.length() < 1 && in.good() )
    {
      s = getWord( in ,  true );
Ruth Mikkelson's avatar
Ruth Mikkelson committed
625
      s = getWord( in ,  false );
626
    }
Ruth Mikkelson's avatar
Ruth Mikkelson committed
627

628
629
    if(  !in.good() )
      throw std::runtime_error( std::string( "Peaks File Quit too fast" ) );
Ruth Mikkelson's avatar
Ruth Mikkelson committed
630

631
632
633
634
635
636
637
638
639
    if( s.compare( std::string( "5" ) ) != 0 )
      throw std::logic_error( "No information on banks given" );//log this only


    bool err = false;
    while( s.compare( std::string( "5" ) ) == 0  )
    {
      double* data = new double[ 15 ];//Use boost pointers here??? does gc
      s = getWord( in ,  false );
Ruth Mikkelson's avatar
Ruth Mikkelson committed
640

641
642
643
644
      if( s.length() < 1 )
        err = true;
      else
        DetNames.push_back( s );
Ruth Mikkelson's avatar
Ruth Mikkelson committed
645

646
      for( int i = 0 ; i < 15 && !err ; i++ )
Ruth Mikkelson's avatar
Ruth Mikkelson committed
647
648
649
650
651
      {
        s = getWord( in ,  false );
        if( s.length() < 1 )
          err = true;
        else
652
          data[ i ] = strtod( s.c_str() ,  0 );
Ruth Mikkelson's avatar
Ruth Mikkelson committed
653

654
      }
Ruth Mikkelson's avatar
Ruth Mikkelson committed
655

656
657
658
659
660
      if( err )
        s = std::string( "2" );
      else
      {
        DetInfo.push_back( data );
Ruth Mikkelson's avatar
Ruth Mikkelson committed
661

662
663
        for( s = getWord( in ,  false ) ; s.length() > 0 ; s = getWord( in ,  false ) )
        {}
Ruth Mikkelson's avatar
Ruth Mikkelson committed
664

665
666
667
668
        s = getWord( in ,  true );
        s = getWord( in ,  false );//no blank lines will be allowed so far
      }
    }
Ruth Mikkelson's avatar
Ruth Mikkelson committed
669
670


671
672
    if( err )
      return std::string();
Ruth Mikkelson's avatar
Ruth Mikkelson committed
673

674
675
    return s;
  }
Ruth Mikkelson's avatar
Ruth Mikkelson committed
676
677
678



679
680
681
682
683
684
685
686
687
  std::string readPeak( std::string lastStr ,  std::ifstream& in ,
      double& h , double& k ,  double& l ,  double& col ,
      double& row , double& chan ,  double& L2 ,
      double&  ScatAng , double& Az ,  double& wl ,
      double& D ,  double& IPK , double& Inti ,
      double& SigI ,  int &iReflag )
  {

    std::string s = lastStr;
Ruth Mikkelson's avatar
Ruth Mikkelson committed
688

689
    if( s.length() < 1 && in.good() )//blank line
Ruth Mikkelson's avatar
Ruth Mikkelson committed
690
    {
691
      readToEndOfLine( in ,  true );
Ruth Mikkelson's avatar
Ruth Mikkelson committed
692

693
694
      s = getWord( in ,  false );;
    }
Ruth Mikkelson's avatar
Ruth Mikkelson committed
695

696
697
    if( s.length() < 1 )
      return 0;
Ruth Mikkelson's avatar
Ruth Mikkelson committed
698

699
700
701
702
703
704
705
706
707

    if( s.compare( "2" ) == 0 )
    {
      readToEndOfLine( in ,  true );

      for( s = getWord( in ,  false ) ; s.length() < 1 && in.good() ;
          s = getWord( in ,  true ) )
      {
        s = getWord( in ,  false );
Ruth Mikkelson's avatar
Ruth Mikkelson committed
708
      }
709
    }
Ruth Mikkelson's avatar
Ruth Mikkelson committed
710

711
712
    if( s.length() < 1 )
      return 0;
Ruth Mikkelson's avatar
Ruth Mikkelson committed
713

714
715
    if( s.compare( "3" ) != 0 )
      return s;
Ruth Mikkelson's avatar
Ruth Mikkelson committed
716

717
    int seqNum = atoi( getWord( in ,  false ).c_str() );
Ruth Mikkelson's avatar
Ruth Mikkelson committed
718

719
720
721
    h = strtod( getWord( in ,  false ).c_str() ,  0 ) ;
    k = strtod( getWord( in ,  false ).c_str() ,  0 ) ;
    l = strtod( getWord( in ,  false ).c_str() ,  0 ) ;
Ruth Mikkelson's avatar
Ruth Mikkelson committed
722

723
724
725
726
727
    col = strtod( getWord( in ,  false ).c_str() ,  0 ) ;
    row = strtod( getWord( in ,  false ).c_str() ,  0 ) ;
    chan = strtod( getWord( in ,  false ).c_str() ,  0 ) ;
    L2 = strtod( getWord( in ,  false ).c_str() ,  0 )/100 ;
    ScatAng = strtod( getWord( in ,  false ).c_str() ,  0 ) ;
Ruth Mikkelson's avatar
Ruth Mikkelson committed
728

729
730
731
732
    Az = strtod( getWord( in , false ).c_str() , 0 ) ;
    wl = strtod( getWord( in , false ).c_str() , 0 ) ;
    D = strtod( getWord( in , false ).c_str() , 0 ) ;
    IPK = strtod( getWord( in , false ).c_str() , 0 ) ;
Ruth Mikkelson's avatar
Ruth Mikkelson committed
733

734
735
736
    Inti = strtod( getWord( in , false ).c_str() , 0 ) ;
    SigI = strtod( getWord( in , false ).c_str() , 0 ) ;
    iReflag = atoi( getWord( in , false ).c_str() ) ;
Ruth Mikkelson's avatar
Ruth Mikkelson committed
737
738


739
740
741
    readToEndOfLine( in ,  true );
    return getWord( in , false );
  }
Ruth Mikkelson's avatar
Ruth Mikkelson committed
742
743


744
745
746
747
748
749
  std::string readPeakBlockHeader( std::string lastStr ,  std::ifstream &in ,
      int&run ,int& detName ,
      double&chi , double&phi ,
      double&omega , double&monCount )
  {
    std::string s = lastStr;
Ruth Mikkelson's avatar
Ruth Mikkelson committed
750

751
752
753
754
755
    if( s.length() < 1 && in.good() )//blank line
    {
      readToEndOfLine( in ,  true );
      s = getWord( in , false );
    }
Ruth Mikkelson's avatar
Ruth Mikkelson committed
756

757
758
759
760
761
    if( s.length() < 1 )
      return std::string();

    if( s.compare( "0" ) == 0 )
    {
Ruth Mikkelson's avatar
Ruth Mikkelson committed
762
      readToEndOfLine( in ,  true );
763
764
765
766
767
768
      s = getWord( in , false );
      while( s.length() < 1 )
      {
        readToEndOfLine( in ,  true );
        s = getWord( in , false );
      }
Ruth Mikkelson's avatar
Ruth Mikkelson committed
769
770
    }

771
772
    if( s.compare( std::string( "1" ) ) != 0 )
      return s;
Ruth Mikkelson's avatar
Ruth Mikkelson committed
773

774
775
776
777
    run = atoi( getWord( in , false ).c_str() );
    detName =atoi( getWord( in , false ).c_str());
    chi = strtod( getWord( in , false ).c_str() , 0 );
    phi = strtod( getWord( in , false ).c_str() , 0 );
Ruth Mikkelson's avatar
Ruth Mikkelson committed
778

779
780
781
    omega = strtod( getWord( in , false ).c_str() , 0 );
    monCount = strtod( getWord( in , false ).c_str() , 0 );
    readToEndOfLine( in ,  true );
Ruth Mikkelson's avatar
Ruth Mikkelson committed
782

783
    return getWord( in , false );
Ruth Mikkelson's avatar
Ruth Mikkelson committed
784

785
786
  }
  //special formatted file Use clear peaks to no append
Ruth Mikkelson's avatar
Ruth Mikkelson committed
787

788
789
790
791
792
  /** Append the peaks from a .peaks file into the workspace
   *
   * @param filename :: path to the .peaks file
   */
  void PeaksWorkspace::append( std::string filename )
793
  {
794
795
    try
    {
Ruth Mikkelson's avatar
Ruth Mikkelson committed
796

797
      std::ifstream in( filename.c_str() );
Ruth Mikkelson's avatar
Ruth Mikkelson committed
798

799
      std::string s = readHeader( in );
Ruth Mikkelson's avatar
Ruth Mikkelson committed
800

801
802
803
      if( !in.good() || s.length() < 1 )
        throw std::runtime_error(
            std::string( "End of Peaks file before peaks" ) );
Ruth Mikkelson's avatar
Ruth Mikkelson committed
804

805
806
807
      if( s.compare( std::string( "0" ) ) != 0 )
        throw std::logic_error(
            std::string( "No header for Peak segments" ) );
Ruth Mikkelson's avatar
Ruth Mikkelson committed
808

809
810
      readToEndOfLine( in ,  true );
      s = getWord( in , false );
Ruth Mikkelson's avatar
Ruth Mikkelson committed
811

812
813
      int run , Reflag, detName;
      double chi , phi , omega , monCount;
Ruth Mikkelson's avatar
Ruth Mikkelson committed
814

815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
      double h , k , l , col , row , chan , L2 , ScatAng , Az , wl , D ,
      IPK , Inti , SigI;
      double x , y , z , r , time;

      while( in.good() )
      {
        s = readPeakBlockHeader( s ,  in  , run , detName , chi , phi ,
            omega , monCount );

        s = readPeak( s , in , h , k , l , col , row , chan , L2 ,
            ScatAng , Az , wl , D , IPK , Inti , SigI , Reflag );

        chi *= pi/180;
        phi *= pi/180;
        omega *= pi/180;

        z = L2*cos( ScatAng );
        r = sqrt( L2*L2 - z*z );
        x = r*cos( Az );
        y = r*sin( Az );
        //Would use the V3D spherical stuff, but it seemed weird
        std::vector< double >xx , yy;
        xx.push_back( wl );

        Kernel::Units::Wavelength WL;
        WL.toTOF( xx , yy , C_L1 , L2 , ScatAng , 12 , 12.1 , 12.1 );
        time = xx[ 0 ];


        this->addPeak( Geometry::V3D( x , y , z ) , time ,
            Geometry::V3D( h , k , l ) ,
            Geometry::V3D( phi , chi , omega ) ,
            Reflag , run , monCount , detName , IPK ,
            row , col , chan , L2, Inti , SigI );
      }
  }catch( std::exception & xe)
  {
852
853
    g_log.error( xe.what());
    throw (std::logic_error(xe.what()));
854
855
856
  }catch( char* str)
  {

857
858
859

    g_log.error( std::string(str));
    throw (std::logic_error(std::string(str)));
860
861
862
  }catch(...)
  {

863
864
    g_log.error( "Unknown Error");
        throw ("Unknown Error");
865
  }
Ruth Mikkelson's avatar
Ruth Mikkelson committed
866

867
  }
Ruth Mikkelson's avatar
Ruth Mikkelson committed
868

869
870
871
872
873
874

  /** Get the d spacing of the given peak
   *
   * @param peakNum :: index of the peak
   * @return double, d_spacing
   */
875
876
877
878
879
880
  double  PeaksWorkspace::get_dspacing( int peakNum )
  {
    double time = cell< double >(peakNum, ItimeCol);
    double L1   =cell< double >(peakNum, IL1Col );
    Geometry::V3D position =
        cell< Geometry::V3D >( peakNum , IpositionCol);
Ruth Mikkelson's avatar
Ruth Mikkelson committed
881

882
    double rho, polar,az;
883
884
885
    position.getSpherical( rho, polar, az);
    polar *= pi/180;
    double L2=rho;
886
887
    std::vector< double >xx , yy;
    xx.push_back( time );
Ruth Mikkelson's avatar
Ruth Mikkelson committed
888

889
890
891
    Kernel::Units::dSpacing dsp;
    dsp.fromTOF( xx , yy , L1 , L2 , polar , 12 , 12.1 , 12.1 );
    return xx[ 0 ];
Ruth Mikkelson's avatar
Ruth Mikkelson committed
892

893
  }
Ruth Mikkelson's avatar
Ruth Mikkelson committed
894

895
896
897
898
899
900
  double  PeaksWorkspace::get_wavelength( int peakNum )
  {
    double time = cell< double >(peakNum, ItimeCol);
    double L1   =cell< double >(peakNum, IL1Col );
    Geometry::V3D position = cell< Geometry::V3D >( peakNum , IpositionCol);
    double rho, polar,az;
901
902
903
    position.getSpherical( rho, polar, az);
    polar *= pi/180;
    double L2=rho;
904
905
    std::vector< double >xx , yy;
    xx.push_back( time );
Ruth Mikkelson's avatar
Ruth Mikkelson committed
906

907
908
909
910
    Kernel::Units::Wavelength wl;
    wl.fromTOF( xx , yy , L1 , L2 , polar , 12 , 12.1 , 12.1 );
    return xx[ 0 ];
  }
Ruth Mikkelson's avatar
Ruth Mikkelson committed
911

912

913
914
915
  //I believe their |Q| = 2pi/d.  This is what is returned.
  double  PeaksWorkspace::get_Qmagnitude( int peakNum )
  {
Ruth Mikkelson's avatar
Ruth Mikkelson committed
916
917


918
919
920
921
    double time = cell< double >(peakNum, ItimeCol);
    double L1   =cell< double >(peakNum, IL1Col );
    Geometry::V3D position = cell< Geometry::V3D >( peakNum , IpositionCol);
    double rho, polar,az;
922
923
924
    position.getSpherical( rho, polar, az);
    polar *= pi/180;
    double L2=rho;
925
926
    std::vector< double >xx , yy;
    xx.push_back( time );
Ruth Mikkelson's avatar
Ruth Mikkelson committed
927

928
929
    Kernel::Units::MomentumTransfer Q;
    Q.fromTOF( xx , yy , L1 , L2 , polar , 12 , 12.1 , 12.1 );
930
    return xx[ 0 ]/2/pi;
931
  }
Ruth Mikkelson's avatar
Ruth Mikkelson committed
932

933
934
935
  Geometry::V3D    PeaksWorkspace::get_Qlab( int peakNum )
  {
    double MagQ = get_Qmagnitude( peakNum);
Ruth Mikkelson's avatar
Ruth Mikkelson committed
936
937


938
939
940
941
942
    Geometry::V3D position =
        cell< Geometry::V3D >( peakNum , IpositionCol);
    position =Geometry::V3D(position);
    position /=position.norm();
    position.setZ( position.Z()-1);
943
944
    double nrm= position.norm();
    position *=MagQ/nrm;
945
    return position;
Ruth Mikkelson's avatar
Ruth Mikkelson committed
946

947

948
  }
Ruth Mikkelson's avatar
Ruth Mikkelson committed
949

950
951
  Geometry::V3D     PeaksWorkspace::get_QXtal( int peakNum )
  {
Ruth Mikkelson's avatar
Ruth Mikkelson committed
952

953
    Geometry::V3D Qvec= Geometry::V3D( get_Qlab( peakNum) );
Ruth Mikkelson's avatar
Ruth Mikkelson committed
954

955
956
    Geometry::V3D sampOrient = getSampleOrientation( peakNum );

957
958
    //phi, chi, omega
    Geometry::Quat rot;
959
    sampOrient *=180/pi;
960
    rot.setAngleAxis( -sampOrient[2],Geometry::V3D( 0,1,0));
Ruth Mikkelson's avatar
Ruth Mikkelson committed
961

962
    rot.rotate( Qvec );
Ruth Mikkelson's avatar
Ruth Mikkelson committed
963

964
965
    rot.setAngleAxis( -sampOrient[1] ,Geometry::V3D( 0,0,1));
    rot.rotate( Qvec );
Ruth Mikkelson's avatar
Ruth Mikkelson committed
966

967
    rot.setAngleAxis( -sampOrient[0],Geometry::V3D( 0,1,0));
Ruth Mikkelson's avatar
Ruth Mikkelson committed
968

969
    rot.rotate( Qvec );
Ruth Mikkelson's avatar
Ruth Mikkelson committed
970
971


972
973
    return Qvec;
  }
Ruth Mikkelson's avatar
Ruth Mikkelson committed
974

975
976
  Geometry::V3D   PeaksWorkspace::get_hkl( int peakNum )
  {
Ruth Mikkelson's avatar
Ruth Mikkelson committed
977

978
979
    return Geometry::V3D( cell< Geometry::V3D >(peakNum , IhklCol ));
  }
Ruth Mikkelson's avatar
Ruth Mikkelson committed
980

981
982
983
984
  double     PeaksWorkspace::get_row(int peakNum )
  {
    return cell< double >( peakNum , IPeakRowCol );
  }
Ruth Mikkelson's avatar
Ruth Mikkelson committed
985

986
987
988
989
990
991
  double     PeaksWorkspace::get_ipk( int peakNum )
  {

    return cell< double >( peakNum , IPeakIntensityCol );
  }

992
993
  double     PeaksWorkspace::get_column( int peakNum )
  {
Ruth Mikkelson's avatar
Ruth Mikkelson committed
994

995
996
    return cell< double >( peakNum , IPeakColCol );
  }
Ruth Mikkelson's avatar
Ruth Mikkelson committed
997

998
999
  double     PeaksWorkspace::get_time_channel( int peakNum )
  {
Ruth Mikkelson's avatar
Ruth Mikkelson committed
1000