Multiple Energy Window Code

recon_test_pack/root_header.hroot

@@ -34,6 +34,14 @@ number of crystals_Z := 4
 
 Singles readout depth := 1
 exclude scattered events := ${EXCLUDE_SCATTERED}
+
+;
+; If set only a certain order of scatter will be kept
+exclude scattered events := ${MAXIMUM_ORDER_OF_SCATTER}
+
+;
+; If set the random events will be skipped
+
 exclude random events := ${EXCLUDE_RANDOM}
 offset (num of detectors) := 0 
 

src/IO/InputStreamFromROOTFile.cxx

@@ -40,12 +40,13 @@ InputStreamFromROOTFile()
 InputStreamFromROOTFile::
 InputStreamFromROOTFile(std::string filename,
                         std::string chain_name,
-                        bool exclude_scattered, bool exclude_randoms,
-                        float low_energy_window, float up_energy_window,
+                        int maximum_order_of_scatter, bool exclude_randoms,
+                        float low_energy_window_1, float up_energy_window_1,
+                        float low_energy_window_2, float up_energy_window_2,
                         int offset_dets)
     : filename(filename), chain_name(chain_name),
-      exclude_scattered(exclude_scattered), exclude_randoms(exclude_randoms),
-      low_energy_window(low_energy_window), up_energy_window(up_energy_window), offset_dets(offset_dets)
+      maximum_order_of_scatter(maximum_order_of_scatter), exclude_randoms(exclude_randoms),
+     low_energy_window_1(low_energy_window_1), up_energy_window_1(up_energy_window_1),low_energy_window_2(low_energy_window_2), up_energy_window_2(up_energy_window_2), offset_dets(offset_dets)
 {
     set_defaults();
     reset();
@@ -56,10 +57,12 @@ InputStreamFromROOTFile::set_defaults()
 {
     starting_stream_position = 0;
     singles_readout_depth = -1;
-    exclude_scattered = false;
+    maximum_order_of_scatter = 0;
     exclude_randoms = false;
-    low_energy_window = 0.f;
-    up_energy_window = 1000.f;
+    low_energy_window_1 = 0.f;
+    up_energy_window_1 = 1000.f;
+    low_energy_window_2 = 0.f;
+    up_energy_window_2 = 1000.f;
     read_optional_root_fields=false;
 }
 
@@ -69,11 +72,13 @@ InputStreamFromROOTFile::initialise_keymap()
     this->parser.add_key("name of data file", &this->filename);
     this->parser.add_key("Singles readout depth", &this->singles_readout_depth);
     this->parser.add_key("name of input TChain", &this->chain_name);
-    this->parser.add_key("exclude scattered events", &this->exclude_scattered);
+    this->parser.add_key("maximum order of scatter", &this->maximum_order_of_scatter);
     this->parser.add_key("exclude random events", &this->exclude_randoms);
     this->parser.add_key("offset (num of detectors)", &this->offset_dets);
-    this->parser.add_key("low energy window (keV)", &this->low_energy_window);
-    this->parser.add_key("upper energy window (keV)", &this->up_energy_window);
+    this->parser.add_key("low energy window 1 (MeV)", &this->low_energy_window_1);
+    this->parser.add_key("upper energy window 1 (MeV)", &this->up_energy_window_1);
+    this->parser.add_key("low energy window 2 (MeV)", &this->low_energy_window_2);
+    this->parser.add_key("upper energy window 2 (MeV)", &this->up_energy_window_2);
     this->parser.add_key("read optional ROOT fields", &this->read_optional_root_fields);
 }
 

src/IO/InputStreamFromROOTFileForCylindricalPET.cxx

@@ -38,8 +38,9 @@ InputStreamFromROOTFileForCylindricalPET(std::string _filename,
                                          int submodule_repeater_x, int submodule_repeater_y, int submodule_repeater_z,
                                          int module_repeater_x, int module_repeater_y, int module_repeater_z,
                                          int rsector_repeater,
-                                         bool _exclude_scattered, bool _exclude_randoms,
-                                         float _low_energy_window, float _up_energy_window,
+                                         int _maximum_order_of_scatter, bool _exclude_randoms,
+                                         float _low_energy_window_1, float _up_energy_window_1,
+                                         float _low_energy_window_2, float _up_energy_window_2,
                                          int _offset_dets):
     base_type(),
     crystal_repeater_x(crystal_repeater_x), crystal_repeater_y(crystal_repeater_y), crystal_repeater_z(crystal_repeater_z),
@@ -51,10 +52,12 @@ InputStreamFromROOTFileForCylindricalPET(std::string _filename,
 
     filename = _filename;
     chain_name = _chain_name;
-    exclude_scattered = _exclude_scattered;
+    maximum_order_of_scatter = _maximum_order_of_scatter;
     exclude_randoms = _exclude_randoms;
-    low_energy_window = _low_energy_window;
-    up_energy_window = _up_energy_window;
+    low_energy_window_1 = _low_energy_window_1;
+    up_energy_window_1 = _up_energy_window_1;
+    low_energy_window_2 = _low_energy_window_2;
+    up_energy_window_2 = _up_energy_window_2;
     offset_dets = _offset_dets;
 
     half_block = module_repeater_y * submodule_repeater_y * crystal_repeater_y / 2  - 1;
@@ -78,14 +81,14 @@ get_next_record(CListRecordROOT& record)
 
         current_position ++ ;
 
-        if ( (this->comptonphantom1 > 0 || this->comptonphantom2 > 0) && this->exclude_scattered )
+        if (comptonphantom1 + comptonphantom2 > maximum_order_of_scatter)
             continue;
         if ( (this->eventID1 != this->eventID2) && this->exclude_randoms)
             continue;
-        if (this->energy1 < this->low_energy_window ||
-                 this->energy1 > this->up_energy_window ||
-                 this->energy2 < this->low_energy_window ||
-                 this->energy2 > this->up_energy_window)
+        if (this->energy1 < this->low_energy_window_1 ||
+             this->energy1 > this->up_energy_window_1 ||
+             this->energy2 < this->low_energy_window_2||
+             this->energy2 > this->up_energy_window_2)
             continue;
 
         break;
@@ -122,7 +125,8 @@ get_next_record(CListRecordROOT& record)
             record.init_from_data(ring1, ring2,
                                   crystal1, crystal2,
                                   time1, time2,
-                                  eventID1, eventID2);
+                                  eventID1, eventID2,
+                                  energy1,energy2);
 }
 
 std::string

src/IO/InputStreamFromROOTFileForECATPET.cxx

@@ -37,8 +37,9 @@ InputStreamFromROOTFileForECATPET(std::string _filename,
                                   std::string _chain_name,
                                   int crystal_repeater_x, int crystal_repeater_y, int crystal_repeater_z,
                                   int block_repeater_y, int block_repeater_z,
-                                  bool _exclude_scattered, bool _exclude_randoms,
-                                  float _low_energy_window, float _up_energy_window,
+                                  int _maximum_order_of_scater, bool _exclude_randoms,
+                                  float _low_energy_window_1, float _up_energy_window_1,
+                                  float _low_energy_window_2, float _up_energy_window_2,
                                   int _offset_dets):
     base_type(),
     crystal_repeater_x(crystal_repeater_x), crystal_repeater_y(crystal_repeater_y), crystal_repeater_z(crystal_repeater_z),
@@ -48,10 +49,12 @@ InputStreamFromROOTFileForECATPET(std::string _filename,
 
     filename = _filename;
     chain_name = _chain_name;
-    exclude_scattered = _exclude_scattered;
+    maximum_order_of_scatter = _maximum_order_of_scater;
     exclude_randoms = _exclude_randoms;
-    low_energy_window = _low_energy_window;
-    up_energy_window = _up_energy_window;
+    low_energy_window_1 = _low_energy_window_1;
+    up_energy_window_1 = _up_energy_window_1;
+    low_energy_window_2 = _low_energy_window_2;
+    up_energy_window_2 = _up_energy_window_2;
     offset_dets = _offset_dets;
 
     half_block = crystal_repeater_y / 2  - 1;
@@ -74,14 +77,14 @@ get_next_record(CListRecordROOT& record)
 
         current_position ++ ;
 
-        if ( (comptonphantom1 > 0 || comptonphantom2 > 0) && exclude_scattered )
+        if (comptonphantom1 + comptonphantom2 > maximum_order_of_scatter)
             continue;
         if ( eventID1 != eventID2 && exclude_randoms )
             continue;
-        if (energy1 < low_energy_window ||
-                 energy1 > up_energy_window ||
-                 energy2 < low_energy_window ||
-                 energy2 > up_energy_window)
+        if (this->energy1 < this->low_energy_window_1 ||
+             this->energy1 > this->up_energy_window_1 ||
+             this->energy2 < this->low_energy_window_2||
+             this->energy2 > this->up_energy_window_2)
             continue;
 
         break;
@@ -112,7 +115,8 @@ get_next_record(CListRecordROOT& record)
             record.init_from_data(ring1, ring2,
                                   crystal1, crystal2,
                                   time1, time2,
-                                  eventID1, eventID2);
+                                  eventID1, eventID2,
+                                  energy1,energy2);
 }
 
 std::string

src/IO/InterfileHeader.cxx

@@ -102,6 +102,7 @@ void MinimalInterfileHeader::set_version_specific_keys()
 InterfileHeader::InterfileHeader()
   : MinimalInterfileHeader()
 {
+
   number_format_values.push_back("bit");
   number_format_values.push_back("ascii");
   number_format_values.push_back("signed integer");
@@ -171,6 +172,7 @@ InterfileHeader::InterfileHeader()
 
 
 
+
   add_key("name of data file", &data_file_name);
   add_key("originating system", &exam_info_sptr->originating_system);
   ignore_key("GENERAL DATA");
@@ -221,15 +223,16 @@ InterfileHeader::InterfileHeader()
   // support for Louvain la Neuve's extension of 3.3
   add_key("quantification units", &lln_quantification_units);
 
-  add_key("number of energy windows",
-    KeyArgument::INT,	(KeywordProcessor)&InterfileHeader::read_num_energy_windows,&num_energy_windows);
+  add_key("energy window pair", KeyArgument::LIST_OF_INTS, (KeywordProcessor)&InterfileHeader::en_window_pair_set, &energy_window_pair);
+  add_key("number of energy windows", KeyArgument::INT,	(KeywordProcessor)&InterfileHeader::read_num_energy_windows,&num_energy_windows);
   add_vectorised_key("energy window lower level", &lower_en_window_thresholds);
   add_vectorised_key("energy window upper level", &upper_en_window_thresholds);
 
   bed_position_horizontal = 0.F;
   add_key("start horizontal bed position (mm)", &bed_position_horizontal);
   bed_position_vertical = 0.F;
   add_key("start vertical bed position (mm)", &bed_position_vertical);
+
 }
 
 void InterfileHeader::set_version_specific_keys()
@@ -239,16 +242,17 @@ void InterfileHeader::set_version_specific_keys()
     {
       info("Setting energy window keys as in STIR3.0");
       // only a single energy window, and non-vectorised
-      remove_key("energy window lower level");
+      /*remove_key("energy window lower level");
       remove_key("energy window upper level");
       add_key("energy window lower level", &lower_en_window_thresholds[0]);
-      add_key("energy window upper level", &upper_en_window_thresholds[0]);
+      add_key("energy window upper level", &upper_en_window_thresholds[0]);*/
     }
 }
 
 // MJ 17/05/2000 made bool
 bool InterfileHeader::post_processing()
 {
+
   if(type_of_data_index<0)
     {
       warning("Interfile Warning: 'type_of_data' keyword required");
@@ -365,19 +369,45 @@ bool InterfileHeader::post_processing()
                lln_quantification_units);
     }      
   } // lln_quantification_units
-  if (num_energy_windows>0)
-    {
-      if (num_energy_windows>1)
-        warning("Currently only reading the first energy window.");
-      if (upper_en_window_thresholds[0] > 0 && lower_en_window_thresholds[0] > 0 )
-        {
-          exam_info_sptr->set_high_energy_thres(upper_en_window_thresholds[0]);
-          exam_info_sptr->set_low_energy_thres(lower_en_window_thresholds[0]);
-        }
-    }
+
+
+ //set the lower and the higher energy thresholds for all the energy windows available. Default: 1.
+
+  //set the number of energy windows
+
+  exam_info_sptr->set_num_energy_windows(num_energy_windows);
+
+   //set the number of energy window pair
+  if (energy_window_pair.size() > 0) {
+
+      if (energy_window_pair.size() != 2)
+          error("should have two.");
+      if (energy_window_pair[0] < 0)
+          error("first window should be >= 0.");
+      if (energy_window_pair[1] < 0)
+          error("second window should be >= 0.");
+      if (energy_window_pair[0] > num_energy_windows)
+          error("The selected window %d exceeds the  number of energy windows %d.\n",energy_window_pair[0],num_energy_windows);
+      if (energy_window_pair[1] > num_energy_windows)
+          error("The selected window %d exceeds the  number of energy windows %d.\n",energy_window_pair[1],num_energy_windows);
+
+      exam_info_sptr->set_energy_window_pair(energy_window_pair);
+  }
+
+    //set the high and low energy window threshold
+
+      for (int i = 0; i < num_energy_windows; ++i)
+      {
+
+          if (upper_en_window_thresholds[i] >0 && lower_en_window_thresholds[i] >0)
+            {
+              exam_info_sptr->set_high_energy_thres(upper_en_window_thresholds[i],i);
+              exam_info_sptr->set_low_energy_thres(lower_en_window_thresholds[i],i);
+             }
+       }
 
   exam_info_sptr->time_frame_definitions = 
-    TimeFrameDefinitions(image_relative_start_times, image_durations);
+  TimeFrameDefinitions(image_relative_start_times, image_durations);
 
   return false;
 
@@ -393,6 +423,14 @@ void InterfileHeader::read_matrix_info()
 
 }
 
+void InterfileHeader::en_window_pair_set()
+{
+
+    energy_window_pair.resize(2);
+    set_variable();
+
+
+}
 void InterfileHeader::read_num_energy_windows()
 {
   set_variable();
@@ -420,7 +458,7 @@ void InterfileHeader::set_type_of_data()
               &PET_data_type_values);
       ignore_key("process status");
       ignore_key("IMAGE DATA DESCRIPTION");
-      // TODO rename keyword 
+      // TODO rename keyword
       add_vectorised_key("data offset in bytes", &data_offset_each_dataset);
 
     }
@@ -499,6 +537,7 @@ bool InterfileImageHeader::post_processing()
   if (InterfileHeader::post_processing() == true)
     return true;
 
+
   if (PET_data_type_values[PET_data_type_index] != "Image")
     { warning("Interfile error: expecting an image\n");  return true; }
 
@@ -633,6 +672,9 @@ void InterfilePDFSHeader::resize_segments_and_set()
 
 }
 
+
+
+
 int InterfilePDFSHeader::find_storage_order()
 {
 
@@ -1312,12 +1354,12 @@ bool InterfilePDFSHeader::post_processing()
 	      scanner_ptr_from_file->parameter_info().c_str());
     }
 
-  
+
   // float azimuthal_angle_sampling =_PI/num_views;
-  
-  
-   
-  
+
+
+
+
   if (is_arccorrected)
     {
       if (effective_central_bin_size_in_cm <= 0)