20
20
21
21
#include " SpecUtils_config.h"
22
22
23
+ #include < string>
24
+ #include < vector>
23
25
#include < cstring>
24
26
#include < fstream>
25
- #include < iostream>
26
27
#include < numeric>
28
+ #include < sstream>
29
+ #include < iostream>
27
30
#include < stdexcept>
28
- #include < string>
29
- #include < vector>
30
31
31
32
#include " rapidxml/rapidxml.hpp"
32
33
33
34
#include " SpecUtils/DateTime.h"
34
- #include " SpecUtils/EnergyCalibration.h"
35
- #include " SpecUtils/RapidXmlUtils.hpp"
36
35
#include " SpecUtils/SpecFile.h"
36
+ #include " SpecUtils/ParseUtils.h"
37
37
#include " SpecUtils/StringAlgo.h"
38
+ #include " SpecUtils/RapidXmlUtils.hpp"
39
+ #include " SpecUtils/EnergyCalibration.h"
40
+
38
41
39
42
using namespace std ;
40
43
@@ -51,18 +54,25 @@ bool SpecFile::load_radiacode_file(const std::string& filename) {
51
54
if (!input.is_open ())
52
55
return false ;
53
56
54
- const bool success = load_from_radiacode (input);
57
+ bool success = load_from_radiacode (input);
55
58
59
+ if ( !success )
60
+ {
61
+ input.seekg (0 );
62
+ success = load_from_radiacode_spectrogram ( input );
63
+ }
64
+
56
65
if (success)
57
66
filename_ = filename;
58
67
59
68
return success;
60
- } // bool load_radiacode_file( const std::string &filename );
69
+ }// bool load_radiacode_file( const std::string &filename );
70
+
61
71
62
72
bool SpecFile::load_from_radiacode (std::istream& input) {
63
73
reset ();
64
74
65
- if ( !input.good ())
75
+ if ( !input.good () )
66
76
return false ;
67
77
68
78
std::unique_lock<std::recursive_mutex> scoped_lock (mutex_);
@@ -72,10 +82,10 @@ bool SpecFile::load_from_radiacode(std::istream& input) {
72
82
73
83
// Determine stream size
74
84
input.seekg (0 , ios::end);
75
- const size_t file_size = static_cast <size_t >(input.tellg () - start_pos);
85
+ size_t file_size = static_cast <size_t >(input.tellg () - start_pos);
76
86
input.seekg (start_pos);
77
87
78
- // The smallest valid 256 channel RadiaCode file I've been able to construct
88
+ // The smallest valid 256 channel RadiaCode XML file I've been able to construct
79
89
// is about 7KB. Typical 1024-channel foreground RC files are about 27KB going
80
90
// up to 31KB for files with many counts per channel. My largest real file
81
91
// with both foreground and background spectra is 59KB. My largest synthetic
@@ -102,21 +112,20 @@ bool SpecFile::load_from_radiacode(std::istream& input) {
102
112
103
113
// Look for some distinctive strings early in the file
104
114
// If they exist, this is probably a RadiaCode file.
105
- int signature_max_offset = 512 ;
106
- const auto fmtver_pos = filedata.find (" <FormatVersion>" );
107
- if (fmtver_pos == string::npos || fmtver_pos > signature_max_offset)
115
+ const size_t signature_max_offset = 512 ;
116
+ const string::size_type fmtver_pos = filedata.find (" <FormatVersion>" );
117
+ if ( (fmtver_pos == string::npos) || ( fmtver_pos > signature_max_offset) )
108
118
return false ;
109
119
110
- const auto dcr_pos = filedata.find (" <DeviceConfigReference>" );
111
- if (dcr_pos == string::npos || dcr_pos > signature_max_offset)
120
+ const string::size_type dcr_pos = filedata.find (" <DeviceConfigReference>" );
121
+ if ( (dcr_pos == string::npos) || ( dcr_pos > signature_max_offset) )
112
122
return false ;
113
123
114
- const auto device_model_pos = filedata.find (" RadiaCode-" );
115
- if (device_model_pos == string::npos ||
116
- device_model_pos > signature_max_offset)
124
+ const string::size_type device_model_pos = filedata.find (" RadiaCode-" );
125
+ if ( (device_model_pos == string::npos) || (device_model_pos > signature_max_offset) )
117
126
return false ;
118
127
119
- if ( device_model_pos < dcr_pos)
128
+ if ( device_model_pos < dcr_pos )
120
129
return false ;
121
130
122
131
try
@@ -372,4 +381,245 @@ bool SpecFile::load_from_radiacode(std::istream& input) {
372
381
return true ;
373
382
}// bool load_from_radiacode( std::istream &input )
374
383
384
+
385
+ bool SpecFile::load_from_radiacode_spectrogram ( std::istream& input )
386
+ {
387
+ reset ();
388
+
389
+ if ( !input.good () )
390
+ return false ;
391
+
392
+ std::unique_lock<std::recursive_mutex> scoped_lock (mutex_);
393
+
394
+ const istream::pos_type start_pos = input.tellg ();
395
+ input.unsetf (ios::skipws);
396
+
397
+ // Determine stream size
398
+ input.seekg (0 , ios::end);
399
+ size_t file_size = static_cast <size_t >(input.tellg () - start_pos);
400
+ input.seekg (start_pos);
401
+
402
+ // If under a kb, definitely not a valid file
403
+ if ( file_size < (1 *1024 ) )
404
+ return false ;
405
+
406
+ const size_t max_header_len = 512 ;
407
+ assert ( max_header_len < file_size );
408
+
409
+ // First we'll check the beginning of the file to make sure it looks like
410
+ // its probably a spectrogram file
411
+ string headerdata;
412
+ headerdata.resize (max_header_len + 1 );
413
+
414
+ input.read (&(headerdata[0 ]), static_cast <streamsize>(max_header_len));
415
+ headerdata[max_header_len] = 0 ; // jic.
416
+ input.clear ();
417
+ input.seekg (start_pos, ios::beg);
418
+
419
+ // Look for some distinctive strings early in the file - for the moment we'll be
420
+ // pretty restrictive in what we require to be there.
421
+ if ( (headerdata.find (" Spectrogram:" ) == string::npos)
422
+ || (headerdata.find (" Accumulation time:" ) == string::npos)
423
+ || (headerdata.find (" Timestamp:" ) == string::npos)
424
+ || (headerdata.find (" Time:" ) == string::npos)
425
+ || (headerdata.find (" Channels:" ) == string::npos) )
426
+ {
427
+ return false ;
428
+ }
429
+
430
+ try
431
+ {
432
+ // We'll start back over at the beginning
433
+ // The fields of the header are tab-seperated - we will rely on this
434
+ string header;
435
+ while ( safe_get_line (input, header, 10 *1024 ) && header.empty () )
436
+ {
437
+ }
438
+
439
+ auto get_header_field_str = [&headerdata]( const string &field, const bool required ) -> string {
440
+ const size_t pos = headerdata.find ( field + " :" );
441
+ if ( (pos == string::npos) && required )
442
+ throw logic_error ( " radiacode expected header field, '" + field + " ', not found" );
443
+ if ( pos == string::npos )
444
+ return " " ;
445
+
446
+ string value = headerdata.substr (pos + field.size () + 1 );
447
+ const size_t tab_pos = value.find (' \t ' );
448
+ if ( tab_pos != string::npos )
449
+ value = value.substr (0 ,tab_pos);
450
+ trim (value);
451
+ return value;
452
+ };// get_header_field_str lambda
453
+
454
+ const string name = get_header_field_str (" Spectrogram" , true );
455
+ const string time_str = get_header_field_str (" Time" , true );
456
+ const string timestamp_str = get_header_field_str (" Timestamp" , true );
457
+ // const string acc_time_str = get_header_field_str("Accumulation time", true);
458
+ const string channels_str = get_header_field_str (" Channels" , true );
459
+ const string serial_num = get_header_field_str (" Device serial" , false );
460
+ const string flags = get_header_field_str (" Flags" , false );
461
+ const string comment = get_header_field_str (" Comment" , false );
462
+
463
+ const time_point_t start_time = time_from_string ( time_str );
464
+
465
+ uint64_t timestamp = 0 ;
466
+ if ( !(stringstream (timestamp_str) >> timestamp) )
467
+ throw runtime_error ( " Unexpected timestamp format" );
468
+
469
+ size_t num_channels = 0 ;
470
+ if ( !(stringstream (channels_str) >> num_channels) || (num_channels < 16 ) || (num_channels > 4096 ) )
471
+ throw runtime_error ( " Invalid 'Channels' field" );
472
+
473
+ vector<shared_ptr<Measurement>> meass;
474
+ // We dont actually have energy calibration info...
475
+ auto energy_cal = make_shared<EnergyCalibration>();
476
+
477
+ int sample_num = 0 ;
478
+ uint64_t last_timestamp = timestamp;
479
+ size_t skipped_lines = 0 , total_lines = 0 ;
480
+ string line;
481
+ while ( safe_get_line (input, line, 64 *1024 ) )
482
+ {
483
+ total_lines += 1 ;
484
+
485
+ // We'll be pretty generous about allowing invalid lines
486
+ if ( (skipped_lines > 5 ) && (total_lines > 10 ) && (skipped_lines > (total_lines/10 )) )
487
+ throw runtime_error ( " To many invalid lines" );
488
+
489
+ trim ( line );
490
+ if ( line.empty () )
491
+ {
492
+ skipped_lines += 1 ;
493
+ continue ;
494
+ }
495
+
496
+ // The second line in the file starts with "Spectrum:", and looks like its probably an icon
497
+ // of the spectrum, and an example one I looked at was 12 kb, we'll skip it
498
+
499
+ if ( !isdigit ( (int )line.front () ) )
500
+ {
501
+ skipped_lines += 1 ;
502
+ continue ;
503
+ }
504
+
505
+ const size_t end_timestamp_pos = line.find (' \t ' );
506
+ if ( end_timestamp_pos == string::npos )
507
+ {
508
+ skipped_lines += 1 ;
509
+ continue ;
510
+ }
511
+
512
+ const string this_timestamp_str = line.substr (0 ,end_timestamp_pos);
513
+ uint64_t this_timestamp = 0 ;
514
+
515
+ if ( !(stringstream (this_timestamp_str) >> this_timestamp) )
516
+ {
517
+ skipped_lines += 1 ;
518
+ continue ;
519
+ }
520
+
521
+ const size_t end_nsecond = line.find (' \t ' , end_timestamp_pos + 1 );
522
+ if ( (end_nsecond == string::npos) || ((end_nsecond + 4 ) > line.size ()) )
523
+ {
524
+ skipped_lines += 1 ;
525
+ continue ;
526
+ }
527
+
528
+ const string num_seconds_str = line.substr (end_timestamp_pos + 1 , end_nsecond - end_timestamp_pos);
529
+ float num_seconds = 0 .0f ;
530
+ if ( !(stringstream (num_seconds_str) >> num_seconds) )
531
+ {
532
+ skipped_lines += 1 ;
533
+ continue ;
534
+ }
535
+
536
+ const char * const counts_start = line.c_str () + end_nsecond + 1 ;
537
+ const size_t counts_str_len = (line.c_str () + line.size ()) - counts_start;
538
+
539
+ auto channel_counts = make_shared<vector<float >>();
540
+ if ( !split_to_floats ( counts_start, counts_str_len, *channel_counts ) )
541
+ {
542
+ cerr << " Failed to split_to_floats" << endl;
543
+ }
544
+
545
+ if ( channel_counts->size () < 2 )
546
+ {
547
+ skipped_lines += 1 ;
548
+ continue ;
549
+ }
550
+
551
+ if ( channel_counts->size () > num_channels )
552
+ throw runtime_error ( " More channel counts than expected" );
553
+
554
+
555
+ vector<string> warnings;
556
+ channel_counts->resize ( num_channels, 0 .0f );
557
+
558
+ float real_time = 0 .0f ;
559
+ if ( last_timestamp <= this_timestamp )
560
+ {
561
+ real_time = num_seconds;
562
+ }else
563
+ {
564
+ real_time = 1.0E-8 * (this_timestamp - last_timestamp);
565
+ if ( fabs (real_time - num_seconds) > 1.5 )
566
+ {
567
+ warnings.push_back ( " Indeterminant real-time: timestamp implied " + to_string (real_time) + " seconds" );
568
+ real_time = num_seconds;
569
+ }
570
+ }
571
+ if ( (real_time < 0 .0f ) || IsInf (real_time) || IsNan (real_time) )
572
+ {
573
+ warnings.push_back ( " Real-time was negative, setting to zero." );
574
+ real_time = 0 .0f ;
575
+ }
576
+
577
+ last_timestamp = this_timestamp;
578
+ const double gamma_sum = std::accumulate ( begin (*channel_counts), end (*channel_counts), 0.0 );
579
+
580
+ auto meas = make_shared<Measurement>();
581
+ meas->real_time_ = real_time;
582
+ meas->live_time_ = real_time;
583
+ meas->gamma_counts_ = channel_counts;
584
+ meas->gamma_count_sum_ = gamma_sum;
585
+ meas->parse_warnings_ = warnings;
586
+ meas->energy_calibration_ = energy_cal;
587
+ meas->sample_number_ = sample_num;
588
+ meas->detector_name_ = " gamma" ;
589
+
590
+ if ( !is_special (start_time) && (this_timestamp > timestamp) )
591
+ {
592
+ const uint64_t nticks = (this_timestamp - timestamp);
593
+ meas->start_time_ = start_time + chrono::milliseconds (nticks/10000 );
594
+ }// if( !is_special(start_time) )
595
+
596
+
597
+ meass.push_back ( meas );
598
+
599
+ sample_num += 1 ;
600
+ }// while( safe_get_line(input, line, 64*1024) )
601
+
602
+ if ( meass.empty () )
603
+ throw runtime_error ( " No measurements" );
604
+
605
+ measurements_ = meass;
606
+
607
+ instrument_type_ = " Spectroscopic Personal Radiation Detector" ;
608
+ manufacturer_ = " Scan-Electronics" ;
609
+ detector_type_ = SpecUtils::DetectorType::RadiaCode;
610
+
611
+ cleanup_after_load ();
612
+ }catch ( std::exception & )
613
+ {
614
+ reset ();
615
+ input.clear ();
616
+ input.seekg (start_pos, ios::beg);
617
+
618
+ return false ;
619
+ }// try / catch
620
+
621
+ return true ;
622
+ }// bool parse_radiacode_spectrogram( std::istream& input )
623
+
624
+
375
625
}// namespace SpecUtils
0 commit comments