diff --git a/documentation/release_5.2.htm b/documentation/release_5.2.htm
index 82ac1762cc..ba6d283e49 100644
--- a/documentation/release_5.2.htm
+++ b/documentation/release_5.2.htm
@@ -31,6 +31,9 @@ <h3>Bug fixes</h3>
 
 <h3>New functionality</h3>
 <ul>
+<li>Global variables in SPECTUB have been substituted by class members, such that multiple SPECTUB projectors can be used.
+    <br /><a href="https://github.com/UCL/STIR/pull/1169/">PR #1169</a>.
+  </li>
   <li>Scatter estimation is now smoothed in axial direction for BlocksOnCylindrical scanners.
     <br /><a href="https://github.com/UCL/STIR/pull/1172/">PR #1172</a>.
   </li>
diff --git a/src/include/stir/recon_buildblock/ProjMatrixByBinSPECTUB.h b/src/include/stir/recon_buildblock/ProjMatrixByBinSPECTUB.h
index eae1ecd495..3ca43a9683 100644
--- a/src/include/stir/recon_buildblock/ProjMatrixByBinSPECTUB.h
+++ b/src/include/stir/recon_buildblock/ProjMatrixByBinSPECTUB.h
@@ -75,7 +75,7 @@ class Bin;
 End Projection Matrix By Bin SPECT UB Parameters:=
 \endverbatim
 */
-
+//using namespace SPECTUB;
 class ProjMatrixByBinSPECTUB : 
   public RegisteredParsingObject<
 	      ProjMatrixByBinSPECTUB,
@@ -90,6 +90,9 @@ class ProjMatrixByBinSPECTUB :
   //! Default constructor (calls set_defaults())
   ProjMatrixByBinSPECTUB();
 
+  // disable copy-constructor as currently unsafe to copy due to bare pointers
+  ProjMatrixByBinSPECTUB(const ProjMatrixByBinSPECTUB&) = delete;
+
   //! Destructor (deallocates UB SPECT memory)
   ~ProjMatrixByBinSPECTUB();
 
@@ -97,7 +100,7 @@ class ProjMatrixByBinSPECTUB :
   virtual void set_up(		 
 		      const shared_ptr<const ProjDataInfo>& proj_data_info_ptr,
                       const shared_ptr<const DiscretisedDensity<3,float> >& density_info_ptr // TODO should be Info only
-                      );
+);
 
   bool get_keep_all_views_in_cache() const;
   //! Enable keeping the matrix in memory
@@ -146,7 +149,6 @@ class ProjMatrixByBinSPECTUB :
 
   void
     set_resolution_model(const float collimator_sigma_0_in_mm, const float collimator_slope_in_mm, const bool full_3D = true);
-
  private:
 
   // parameters that will be parsed
@@ -172,6 +174,9 @@ class ProjMatrixByBinSPECTUB :
 
   bool already_setup;
 
+  mutable SPECTUB::wm_da_type wm;
+  mutable SPECTUB::wmh_type wmh; // this could be an arry of wmh_type for each index
+  float * Rrad;
 
   virtual void 
     calculate_proj_matrix_elems_for_one_bin(
@@ -206,7 +211,8 @@ class ProjMatrixByBinSPECTUB :
   int maxszb;
 
 	
-  void compute_one_subset(const int kOS) const;
+  void compute_one_subset(const int kOS,
+                          const float *Rrad) const;
   void delete_UB_SPECT_arrays();
   mutable std::vector<bool> subset_already_processed;
 };
diff --git a/src/include/stir/recon_buildblock/SPECTUB_Tools.h b/src/include/stir/recon_buildblock/SPECTUB_Tools.h
index 2f7d738ab0..467aae3e36 100644
--- a/src/include/stir/recon_buildblock/SPECTUB_Tools.h
+++ b/src/include/stir/recon_buildblock/SPECTUB_Tools.h
@@ -319,16 +319,16 @@ typedef struct
 //... functions from wmtools_SPECT.cpp .........................................
 
 
-void write_wm_FC ();                               // to write double array weight matrix 
+void write_wm_FC (wm_da_type &wm);                               // to write double array weight matrix
 
-void write_wm_hdr ();                              // to write header of a matrix
+void write_wm_hdr (wm_da_type &wm, wmh_type &wmh);                              // to write header of a matrix
 
-void write_wm_STIR ();                             // to write matrix in STIR format
+void write_wm_STIR (SPECTUB::wm_da_type& wm);                             // to write matrix in STIR format
 
 
-void index_calc ( int *indexs );                   // to calculate projection index order in subsets 
+void index_calc (int *indexs , wmh_type &wmh);                   // to calculate projection index order in subsets
 
-void read_Rrad ();								   // to read variable rotation radius from a text file (1 radius per line)
+void read_Rrad (float *Rrad, wmh_type &wmh);								   // to read variable rotation radius from a text file (1 radius per line)
 
 //
 //void col_params ( collim_type *COL );              // to fill collimator structure
@@ -336,21 +336,21 @@ void read_Rrad ();								   // to read variable rotation radius from a text fil
 //void read_col_params ( collim_type *COL);          // to read collimator parameters from a file
 
 
-void fill_ang ( angle_type *ang );				   // to fill angle structure
+void fill_ang (angle_type *ang , wmh_type &wmh, const float *Rrad);				   // to fill angle structure
 
-void generate_msk ( bool *msk_3d, bool *msk_2d, float *att, volume_type *vol); // to create a boolean mask for wm (no weights outside the msk)
+void generate_msk (bool *msk_3d, bool *msk_2d, float *att, volume_type *vol, wmh_type &wmh); // to create a boolean mask for wm (no weights outside the msk)
 
-void read_msk_file ( bool * msk );                 // to read mask from a file
+void read_msk_file (bool * msk , wmh_type &wmh);                 // to read mask from a file
 
 
-void read_att_map ( float *attmap );               // to read attenuation map from a file
+void read_att_map (float *attmap , wmh_type &wmh);               // to read attenuation map from a file
 
 
-int  max_psf_szb ( angle_type *ang );
+int  max_psf_szb (angle_type *ang , wmh_type &wmh);
 
 float calc_sigma_h ( voxel_type vox, collim_type COL);
 
-float calc_sigma_v ( voxel_type vox, collim_type COL);
+float calc_sigma_v (voxel_type vox, collim_type COL, wmh_type &wmh);
 
 
 char *itoa ( int n, char *s);                      // to conver integer to ascii
diff --git a/src/include/stir/recon_buildblock/SPECTUB_Weight3d.h b/src/include/stir/recon_buildblock/SPECTUB_Weight3d.h
index 5b5f28d88d..8f15d93467 100644
--- a/src/include/stir/recon_buildblock/SPECTUB_Weight3d.h
+++ b/src/include/stir/recon_buildblock/SPECTUB_Weight3d.h
@@ -23,68 +23,62 @@
   Berta Marti Fuster, Carles Falcon, Charalampos Tsoumpas, Lefteris Livieratos, Pablo Aguiar, Albert Cot, Domenec Ros and Kris Thielemans,
   <i>Integration of advanced 3D SPECT modeling into the open-source STIR framework</i>,
   Med. Phys. 40, 092502 (2013); http://dx.doi.org/10.1118/1.4816676
-
-  \todo Variables wm, wmh and Rrad are currently global variables. This means that this code would be very dangerous
-  in a parallel setting.
 */
 
 namespace SPECTUB {
 
-/* Global variables: the matrix etc TODO 
-   Need to be defined elsewhere
-
-   A lot of the functions below modify these variables.
-*/
-  extern wm_da_type wm; //! weight (or probability) matrix
-  extern wmh_type wmh;  //! information to construct wm 
-  extern float * Rrad;  //! radii per view
-
-
-void wm_calculation( const int kOS,
-					const angle_type *const ang, 
-					voxel_type vox, 
-					bin_type bin, 
-					const volume_type& vol, 
-					const proj_type& prj, 
-					const float *attmap,
-					const bool *msk_3d,
-					const bool *msk_2d,
-					const int maxszb,
-					const discrf_type *const gaussdens,
-					const int *const  NITEMS
-					);
+void wm_calculation(const int kOS,
+                    const SPECTUB::angle_type *const ang,
+                    SPECTUB::voxel_type vox,
+                    bin_type bin,
+                    const SPECTUB::volume_type& vol,
+                    const proj_type& prj,
+                    const float *attmap,
+                    const bool *msk_3d,
+                    const bool *msk_2d,
+                    const int maxszb,
+                    const SPECTUB::discrf_type *const gaussdens,
+                    const int *const  NITEMS,
+                    SPECTUB::wm_da_type& wm,
+                    SPECTUB::wmh_type& wmh,
+                    const float *Rrad
+                    );
 
 void wm_size_estimation (int kOS,
-						 const angle_type * const ang, 
-						 voxel_type vox, 
-						 bin_type bin, 
-						 const volume_type& vol, 
-						 const proj_type& prj, 
-						 const bool * const msk_3d,
-						 const bool *const msk_2d,
-						 const int maxszb,
-						 const discrf_type * const gaussdens,
-						 int *NITEMS);
+                         const SPECTUB::angle_type * const ang,
+                         SPECTUB::voxel_type vox,
+                         bin_type bin,
+                         const SPECTUB::volume_type& vol,
+                         const SPECTUB::proj_type& prj,
+                         const bool * const msk_3d,
+                         const bool *const msk_2d,
+                         const int maxszb,
+                         const SPECTUB::discrf_type * const gaussdens,
+                         int *NITEMS,
+                         SPECTUB::wmh_type &wmh,
+                         const float *Rrad);
 
 
 //... geometric component ............................................
 
-void calc_gauss ( discrf_type *gaussdens );             
+void calc_gauss ( SPECTUB::discrf_type *gaussdens );
 
-void calc_vxprj ( angle_type *ang );
+void calc_vxprj ( SPECTUB::angle_type *ang );
 
-void voxel_projection ( voxel_type *vox, float * eff, float lngcmd2 );
+void voxel_projection (SPECTUB::voxel_type *vox, float * eff, float lngcmd2 , SPECTUB::wmh_type &wmh);
 
-void fill_psf_no( psf2da_type *psf, psf1d_type *psf1d_h, const voxel_type& vox, const angle_type * const ang, float szdx );
+void fill_psf_no(SPECTUB::psf2da_type *psf, SPECTUB::psf1d_type *psf1d_h, const SPECTUB::voxel_type& vox, const angle_type * const ang, float szdx , SPECTUB::wmh_type &wmh);
 
-void fill_psf_2d( psf2da_type *psf, psf1d_type *psf1d_h, const voxel_type& vox, discrf_type const*const gaussdens, float szdx );
+void fill_psf_2d( SPECTUB::psf2da_type *psf, SPECTUB::psf1d_type *psf1d_h, const SPECTUB::voxel_type& vox, SPECTUB::discrf_type const*const gaussdens, float szdx, SPECTUB::wmh_type &wmh);
 
-void fill_psf_3d(psf2da_type *psf,
-                 psf1d_type *psf1d_h,
-                 psf1d_type *psf1d_v,
-                 const voxel_type& vox, discrf_type const * const gaussdens, float szdx, float thdx, float thcmd2 );
+void fill_psf_3d(SPECTUB::psf2da_type *psf,
+                 SPECTUB::psf1d_type *psf1d_h,
+                 SPECTUB::psf1d_type *psf1d_v,
+                 const SPECTUB::voxel_type& vox, SPECTUB::discrf_type const * const gaussdens,
+                 float szdx, float thdx, float thcmd2,
+                 wmh_type &wmh);
 
-void calc_psf_bin ( float center_psf, float binszcm, discrf_type const * const vxprj, psf1d_type *psf );
+void calc_psf_bin (float center_psf, float binszcm, SPECTUB::discrf_type const * const vxprj, SPECTUB::psf1d_type *psf , SPECTUB::wmh_type &wmh);
 
 
 //... attenuation...................................................
@@ -94,9 +88,9 @@ void calc_psf_bin ( float center_psf, float binszcm, discrf_type const * const v
 
 //void size_attpth_full( psf2da_type *psf, voxel_type vox, volume_type vol, float *att, angle_type *ang );
 
-void calc_att_path ( const bin_type& bin, const voxel_type& vox, const volume_type& vol, attpth_type *attpth );
+void calc_att_path ( const bin_type& bin, const SPECTUB::voxel_type& vox, const SPECTUB::volume_type& vol, SPECTUB::attpth_type *attpth );
 
-float calc_att ( const attpth_type *const attpth, const float *const attmap, int islc );
+float calc_att (const SPECTUB::attpth_type *const attpth, const float *const attmap, int islc , SPECTUB::wmh_type &wmh);
 
 int comp_dist ( float dx, float dy, float dz, float dlast );
 
diff --git a/src/recon_buildblock/ProjMatrixByBinSPECTUB.cxx b/src/recon_buildblock/ProjMatrixByBinSPECTUB.cxx
index c6ae155aa2..8516d6de89 100644
--- a/src/recon_buildblock/ProjMatrixByBinSPECTUB.cxx
+++ b/src/recon_buildblock/ProjMatrixByBinSPECTUB.cxx
@@ -2,6 +2,7 @@
     Copyright (C) 2013, Institute for Bioengineering of Catalonia
     Copyright (C) Biomedical Image Group (GIB), Universitat de Barcelona, Barcelona, Spain.
     Copyright (C) 2013-2014, 2019, 2020 University College London
+    Copyright (C) 2023 National Physical Laboratory
     This file is part of STIR.
 
     SPDX-License-Identifier: Apache-2.0
@@ -17,16 +18,14 @@
   \author Berta Marti Fuster
   \author Carles Falcon
   \author Kris Thielemans
+  \author Daniel Deidda
 */
 
-//#include "stir/ProjDataInterfile.h"
 #include "stir/recon_buildblock/ProjMatrixByBinSPECTUB.h"
 #include "stir/recon_buildblock/TrivialDataSymmetriesForBins.h"
 #include "stir/ProjDataInfoCylindricalArcCorr.h"
-//#include "stir/KeyParser.h"
 #include "stir/IO/read_from_file.h"
 #include "stir/ProjDataInfo.h"
-//#include "stir/utilities.h"
 #include "stir/VoxelsOnCartesianGrid.h"
 #include "stir/Succeeded.h"
 #include "stir/is_null_ptr.h"
@@ -53,22 +52,14 @@
 #include <stdlib.h>
 #include <math.h>
 #include <ctype.h>
-//#include <time.h>
 
 //... user defined libraries .............................................................
 
 #include "stir/recon_buildblock/SPECTUB_Weight3d.h"
 
-/* UB-SPECT global variables */
-namespace SPECTUB {
-  wm_da_type wm;
-  wmh_type wmh; 
-  float * Rrad;
-}
-
+using namespace SPECTUB;
 START_NAMESPACE_STIR
 
-
 const char * const 
 ProjMatrixByBinSPECTUB::registered_name =
 "SPECT UB";
@@ -248,8 +239,6 @@ set_up(
     }
 #endif
 
-  using namespace SPECTUB;
-
    const VoxelsOnCartesianGrid<float> * image_info_ptr =
       dynamic_cast<const VoxelsOnCartesianGrid<float>*> (density_info_ptr.get());
 
@@ -502,12 +491,12 @@ set_up(
 	//... to sort angles into subsets ......................................
 
 	prj.order = new int [ prj.Nang ];
-	index_calc( prj.order );
+    index_calc( prj.order, wmh);
 
 	//... to fill ang structure ............................................
 
 	ang = new angle_type [ prj.Nang ];		
-	fill_ang( ang );			   
+    fill_ang( ang,wmh, Rrad );
 
 	//... to fill high resolution discrete distribution functions ..............
 
@@ -574,12 +563,12 @@ set_up(
                 // we do this to avoid using its own read_msk_file
                 wmh.do_msk_file = false;
 	        wmh.do_msk_att = true;
-                generate_msk( msk_3d, msk_2d, mask_from_file, &vol);
+                generate_msk( msk_3d, msk_2d, mask_from_file, &vol, wmh);
                 delete[] mask_from_file;
               }
             else
               {
-		generate_msk( msk_3d, msk_2d, attmap, &vol);
+        generate_msk( msk_3d, msk_2d, attmap, &vol, wmh);
               }
           }
 	else msk_2d = msk_3d = NULL;
@@ -591,7 +580,7 @@ set_up(
 
 	//... setting PSF maximum size (in bins) and memory allocation for PSF values .......
 
-	this->maxszb = max_psf_szb( ang );  // maximum PSF size (horizontal component of PSF)
+    this->maxszb = max_psf_szb( ang, wmh );  // maximum PSF size (horizontal component of PSF)
 	NITEMS = new int * [prj.NOS];
 	for (int kOS=0; kOS<prj.NOS; ++kOS) {
 	  NITEMS[kOS] = new int [ wm.NbOS ];
@@ -658,7 +647,8 @@ set_up(
 
 		//... size estimations ........................................................
 
-		wm_size_estimation ( kOS,  ang, vox, bin, vol, prj, msk_3d, msk_2d, maxszb, &gaussdens, NITEMS[kOS] );
+        wm_size_estimation ( kOS,  ang, vox, bin, vol, prj, msk_3d, msk_2d, maxszb, &gaussdens, NITEMS[kOS],
+                             wmh,Rrad );
 
 		//cout << "\nwm_SPECT. Size estimation done. time (s): " << double( clock()-ini )/CLOCKS_PER_SEC <<std::endl;
 
@@ -686,7 +676,6 @@ delete_UB_SPECT_arrays()
   if (!this->already_setup)
     return;
   //... freeing matrix memory....................................
-  using namespace SPECTUB;
   delete [] Rrad;
 
   if ( !wmh.do_psf ){
@@ -734,9 +723,9 @@ delete_UB_SPECT_arrays()
 }
 void
 ProjMatrixByBinSPECTUB::
-compute_one_subset(const int kOS) const
+compute_one_subset(const int kOS,
+                   const float *Rrad) const
 {
-  using namespace SPECTUB;
 
   CPUTimer timer;
   timer.start();
@@ -744,12 +733,12 @@ compute_one_subset(const int kOS) const
 
   //... to fill wmh fields related to the subset ..................................
 
-  wmh.subset_ind = kOS;
+  this->wmh.subset_ind = kOS;
 
   for ( int i = 0 ; i < prj.NangOS ; i ++ ){
 
-    wmh.index[ i ] = prj.order[ i + kOS * prj.NangOS ];
-    wmh.Rrad [ i ] = Rrad[ wmh.index[ i ] ];
+    this->wmh.index[ i ] = prj.order[ i + kOS * prj.NangOS ];
+    this->wmh.Rrad [ i ] = Rrad[ this->wmh.index[ i ] ];
   }
 
   //... NITEMS initialization  ......................
@@ -765,26 +754,26 @@ compute_one_subset(const int kOS) const
 
   int ne = 0;
 
-  for ( int i = 0 ; i < wmh.prj.NbOS ; i++ ) ne += NITEMS[kOS][ i ];
+  for ( int i = 0 ; i < this->wmh.prj.NbOS ; i++ ) ne += NITEMS[kOS][ i ];
 
   //... size information ....................................................................
 
   info(boost::format("total number of non-zero weights in this view: %1%, estimated size: %2% MB") 
        % ne
-       % ( wm.do_save_STIR ?  (ne + 10* prj.NbOS)/104857.6 : ne/131072),
+       % ( this->wm.do_save_STIR ?  (ne + 10* prj.NbOS)/104857.6 : ne/131072),
        2);
 
   //... memory allocation for wm float arrays ...................................
 
-  for( int i = 0 ; i < wmh.prj.NbOS ; i++ ){
+  for( int i = 0 ; i < this->wmh.prj.NbOS ; i++ ){
 
-    if ( ( wm.val[ i ] = new (std::nothrow) float [ NITEMS[kOS][ i ] ]) == NULL) 
+    if ( ( this->wm.val[ i ] = new (std::nothrow) float [ NITEMS[kOS][ i ] ]) == NULL)
       {
         //error_wm_SPECT( 200, "wm.val[][]" );
         error("Error allocating space to store values for SPECTUB matrix");
       }
 
-    if ( ( wm.col[ i ] = new (std::nothrow) int   [ NITEMS[kOS][ i ] ]) == NULL) 
+    if ( ( this->wm.col[ i ] = new (std::nothrow) int   [ NITEMS[kOS][ i ] ]) == NULL)
       {
         //error_wm_SPECT( 200, "wm.col[]" );
         error("Error allocating space to store column indices for SPECTUB matrix");
@@ -793,46 +782,47 @@ compute_one_subset(const int kOS) const
 
   //... to initialize wm to zero ......................
 
-  for ( int i = 0 ; i < wm.NbOS ; i++ ){
+  for ( int i = 0 ; i < this->wm.NbOS ; i++ ){
 
-    wm.ne[ i ] = 0;
+    this->wm.ne[ i ] = 0;
 
     for( int j = 0 ; j < NITEMS[kOS][ i ] ; j++ ){
 
-      wm.val[ i ][ j ] = (float)0.;
-      wm.col[ i ][ j ] = 0;
+      this->wm.val[ i ][ j ] = (float)0.;
+      this->wm.col[ i ][ j ] = 0;
     }
   }
-  wm.ne[ wm.NbOS ] = 0;
+  this->wm.ne[ this->wm.NbOS ] = 0;
 
   //... wm calculation for this subset ...........................
 
-  wm_calculation ( kOS, ang, vox, bin, vol, prj, attmap, msk_3d, msk_2d, maxszb, &gaussdens, NITEMS[kOS] );
+  wm_calculation ( kOS, ang, vox, bin, vol, prj, attmap, msk_3d, msk_2d, maxszb, &gaussdens, NITEMS[kOS],
+                   this->wm, this->wmh, Rrad );
   info(boost::format("Weight matrix calculation done. time %1% (s)") % timer.value(),
        2);
 
   //... fill lor .........................
 
-  for( int j = 0 ; j < wm.NbOS ; j++ ){
+  for( int j = 0 ; j < this->wm.NbOS ; j++ ){
     ProjMatrixElemsForOneBin lor;
     Bin bin;
     bin.segment_num()=0;	
-    bin.view_num()=wm.na [ j ];	
-    bin.axial_pos_num()=wm.ns [ j ];	
-    bin.tangential_pos_num()=wm.nb [ j ];	
+    bin.view_num()=this->wm.na [ j ];
+    bin.axial_pos_num()=this->wm.ns [ j ];
+    bin.tangential_pos_num()=this->wm.nb [ j ];
     bin.set_bin_value(0);
     lor.set_bin(bin);
 
-    lor.reserve(wm.ne[ j ]);
-    for ( int i = 0 ; i < wm.ne[ j ] ; i++ ){
+    lor.reserve(this->wm.ne[ j ]);
+    for ( int i = 0 ; i < this->wm.ne[ j ] ; i++ ){
 
       const ProjMatrixElemsForOneBin::value_type 
-        elem(Coordinate3D<int>(wm.nz[ wm.col[ j ][ i ] ],wm.ny[ wm.col[ j ][ i ] ],wm.nx[ wm.col[ j ][ i ] ]), wm.val[ j ][ i ]);      
+        elem(Coordinate3D<int>(this->wm.nz[ this->wm.col[ j ][ i ] ],this->wm.ny[ this->wm.col[ j ][ i ] ],this->wm.nx[ this->wm.col[ j ][ i ] ]), this->wm.val[ j ][ i ]);
       lor.push_back( elem);	
     }
 
-    delete [] wm.val[ j ];
-    delete [] wm.col[ j ];
+    delete [] this->wm.val[ j ];
+    delete [] this->wm.col[ j ];
 
     this->cache_proj_matrix_elems_for_one_bin(lor);
   }
@@ -868,8 +858,8 @@ calculate_proj_matrix_elems_for_one_bin(ProjMatrixElemsForOneBin& lor
 	  subset_already_processed.assign(prj.NOS,false);
 	}
       info(boost::format("Computing matrix elements for view %1%") % view_num,
-        2);
-      compute_one_subset(kOS);
+        2);// potentially pass a wm, wmh[threadh] not sure if works then in setup we need an array of wmh
+      compute_one_subset(kOS,Rrad);
       subset_already_processed[kOS]=true;
     }
   lor.erase();
diff --git a/src/recon_buildblock/SPECTUB_Tools.cxx b/src/recon_buildblock/SPECTUB_Tools.cxx
index 989f893fdc..e71a9f0e70 100644
--- a/src/recon_buildblock/SPECTUB_Tools.cxx
+++ b/src/recon_buildblock/SPECTUB_Tools.cxx
@@ -22,7 +22,6 @@
 #include <boost/math/constants/constants.hpp>
 
 using namespace std;
-//using std::string;
 
 //... user defined libraries .......................................
 
@@ -45,18 +44,11 @@ namespace SPECTUB {
 #define DELIMITER1 '#' //delimiter character in input parameter text file
 #define DELIMITER2 '%' //delimiter character in input parameter text file
 
-//... global variables ..............................................
-
-extern wm_da_type wm;
-extern wmh_type wmh; 
-extern float * Rrad;
-
-
 //=============================================================================
 //=== write_wm_FC =============================================================
 //=============================================================================
 
-void write_wm_FC()
+void write_wm_FC(SPECTUB::wm_da_type& wm)
 {
 	FILE *fid;
 	
@@ -110,7 +102,7 @@ void write_wm_FC()
 //=== write_wm_hdr ============================================================
 //=============================================================================
 
-void write_wm_hdr()
+void write_wm_hdr(SPECTUB::wm_da_type& wm, SPECTUB::wmh_type& wmh)
 {
 	ofstream stream1( wm.fn_hdr.c_str() );
 	if( !stream1 ) error_wmtools_SPECT( 31, wm.fn_hdr );  
@@ -192,7 +184,7 @@ void write_wm_hdr()
 //=== write_wm_STIR ===========================================================
 //=============================================================================
 
-void write_wm_STIR()
+void write_wm_STIR(wm_da_type &wm)
 {
 	int seg_num = 0;             // segment number for STIR matrix (always zero)
 	FILE *fid;
@@ -228,7 +220,7 @@ void write_wm_STIR()
 //=== index_calc ==============================================================
 //=============================================================================
 
-void index_calc ( int *indexs )
+void index_calc ( int *indexs , SPECTUB::wmh_type& wmh)
 {
 	if ( wmh.prj.NOS == 1 ){
 		for ( int i = 0 ; i < wmh.prj.Nang ; i++ ){
@@ -332,7 +324,7 @@ void index_calc ( int *indexs )
 //=== read rotation radius ==================================================
 //=============================================================================
 
-void read_Rrad()
+void read_Rrad(float * Rrad, SPECTUB::wmh_type& wmh)
 {
 	string line;
 	ifstream stream1( wmh.Rrad_fn.c_str() );
@@ -529,7 +521,7 @@ void read_Rrad()
 //=== calc_sigma_v =========================================================
 //==========================================================================
 
-float calc_sigma_v( voxel_type	vox, collim_type COL)
+float calc_sigma_v( voxel_type	vox, collim_type COL, SPECTUB::wmh_type& wmh)
 {
 	float sigma;
 	if ( COL.do_fb ){
@@ -547,7 +539,7 @@ float calc_sigma_v( voxel_type	vox, collim_type COL)
 //=== fill_ang ================================================================
 //=============================================================================
 
-void fill_ang ( angle_type *ang )
+void fill_ang ( angle_type *ang, SPECTUB::wmh_type& wmh,const float * Rrad)
 {
 	float DX    = (float) 0.5 / wmh.psfres ;
 	float dg2rd = boost::math::constants::pi<float>() / (float)180. ;
@@ -614,7 +606,7 @@ void fill_ang ( angle_type *ang )
 //=== generate msk ============================================================
 //=============================================================================
 
-void generate_msk ( bool *msk_3d, bool *msk_2d, float *attmap, volume_type * vol )
+void generate_msk ( bool *msk_3d, bool *msk_2d, float *attmap, volume_type * vol,SPECTUB::wmh_type& wmh )
 {
 	//... initialzation of msk to true .........................
 	
@@ -672,7 +664,7 @@ void generate_msk ( bool *msk_3d, bool *msk_2d, float *attmap, volume_type * vol
 		else {
 			//... to read a mask from a (int) file ....................
 			
-			if ( wmh.do_msk_file ) read_msk_file( msk_3d );             
+            if ( wmh.do_msk_file ) read_msk_file( msk_3d, wmh);
 		}
 	}
 
@@ -717,7 +709,7 @@ void generate_msk ( bool *msk_3d, bool *msk_2d, float *attmap, volume_type * vol
 //=== read_mask file ==========================================================
 //=============================================================================
 
-void read_msk_file( bool *msk )
+void read_msk_file( bool *msk, SPECTUB::wmh_type& wmh )
 {
 	FILE *fid;
 	int *aux;
@@ -739,7 +731,7 @@ void read_msk_file( bool *msk )
 //=== read_att_map ============================================================
 //=============================================================================
 
-void read_att_map( float *attmap )
+void read_att_map( float *attmap,SPECTUB::wmh_type& wmh )
 {
 	FILE *fid;
 	if ( ( fid = fopen( wmh.att_fn.c_str() , "rb") ) == NULL ) error_wmtools_SPECT ( 124, wmh.att_fn );
@@ -763,7 +755,7 @@ void read_att_map( float *attmap )
 //=== max_psf_szb ==========================================================
 //==========================================================================
 
-int max_psf_szb( angle_type *ang )
+int max_psf_szb( angle_type *ang, SPECTUB::wmh_type& wmh)
 { 
 	int maxszb;
 	float Rrad_max = ang[0].Rrad;
@@ -803,7 +795,7 @@ int max_psf_szb( angle_type *ang )
 		maxszb = (int)floor( wmh.maxsigm * (float)2. * sig_h_max_cm / wmh.prj.szcm ) + 3;
 	
 		if ( wmh.do_psf_3d ){
-			float sig_v_max_cm = calc_sigma_v( vox, wmh.COL );
+            float sig_v_max_cm = calc_sigma_v( vox, wmh.COL, wmh );
 			int maxszb_v = (int)floor( wmh.maxsigm * (float)2. * sig_v_max_cm / wmh.prj.thcm ) + 3;
 			maxszb = maxim( maxszb , maxszb_v ); 	
 		}
diff --git a/src/recon_buildblock/SPECTUB_Weight3d.cxx b/src/recon_buildblock/SPECTUB_Weight3d.cxx
index 1348cf743d..c254004a5b 100644
--- a/src/recon_buildblock/SPECTUB_Weight3d.cxx
+++ b/src/recon_buildblock/SPECTUB_Weight3d.cxx
@@ -45,18 +45,18 @@ using namespace std;
 //=== wm_calculation =======================================================
 //==========================================================================
 
-void wm_calculation( const int kOS,
-					const angle_type *const ang, 
-					voxel_type vox, 
-				        bin_type bin, 
-					const volume_type& vol, 
-					const proj_type& prj, 
-					const float *attmap,
-					const bool *msk_3d,
-					const bool *msk_2d,
-					const int maxszb,
-					const discrf_type *const gaussdens,
-		     const int *const  NITEMS)
+void wm_calculation(const int kOS,
+                    const angle_type * const ang,
+                    voxel_type vox,
+                        bin_type bin,
+                    const volume_type &vol,
+                    const proj_type& prj,
+                    const float *attmap,
+                    const bool *msk_3d,
+                    const bool *msk_2d,
+                    const int maxszb,
+                    const discrf_type * const gaussdens,
+             const int *const  NITEMS, wm_da_type &wm, wmh_type &wmh, const float *Rrad)
 {
 	
 	float weight;
@@ -174,7 +174,7 @@ void wm_calculation( const int kOS,
 				
 				//... to project voxels onto the detection plane and to calculate other distances .....
 				
-                voxel_projection( &vox , &eff , prj.lngcmd2 );
+                voxel_projection( &vox , &eff , prj.lngcmd2, wmh );
 				
 				//... setting PSF to zero	.........................................	
 				
@@ -186,13 +186,13 @@ void wm_calculation( const int kOS,
 				
 				//... correction for PSF ..............................
 				
-				if ( !wmh.do_psf  )	fill_psf_no ( &psf, &psf1d_h, vox, &ang[ ka ], bin.szdx );
+                if ( !wmh.do_psf  )	fill_psf_no ( &psf, &psf1d_h, vox, &ang[ ka ], bin.szdx, wmh);
 				
 				else{
 					
-					if ( wmh.do_psf_3d ) fill_psf_3d ( &psf, &psf1d_h, &psf1d_v, vox, gaussdens, bin.szdx, bin.thdx, bin.thcmd2 );
+                    if ( wmh.do_psf_3d ) fill_psf_3d ( &psf, &psf1d_h, &psf1d_v, vox, gaussdens, bin.szdx, bin.thdx, bin.thcmd2, wmh );
 					
-					else fill_psf_2d ( &psf, &psf1d_h, vox, gaussdens, bin.szdx );
+                    else fill_psf_2d ( &psf, &psf1d_h, vox, gaussdens, bin.szdx, wmh);
 				}
 				
 				//... correction for attenuation .................................................
@@ -232,7 +232,7 @@ void wm_calculation( const int kOS,
 						if ( !msk_3d[ vox.iv ] ) continue;
 					}
 					
-					if ( wmh.do_att && !wmh.do_full_att ) coeff_att = calc_att( &attpth[ 0 ], attmap , vox.islc );
+                    if ( wmh.do_att && !wmh.do_full_att ) coeff_att = calc_att( &attpth[ 0 ], attmap , vox.islc, wmh);
 					
 					//... weight matrix values calculation .......................................
 					
@@ -248,7 +248,7 @@ void wm_calculation( const int kOS,
 
 						jp = k * prj.Nbp + ks * prj.Nbin + psf.ib[ ie ];
 						
-						if ( wmh.do_full_att ) coeff_att = calc_att( &attpth[ ie ], attmap, vox.islc );
+                        if ( wmh.do_full_att ) coeff_att = calc_att( &attpth[ ie ], attmap, vox.islc, wmh );
 						
 						weight = psf.val[ ie ] * eff * coeff_att ;
                         
@@ -312,7 +312,9 @@ void wm_size_estimation (int kOS,
 						 const bool *const msk_2d,
 						 const int maxszb,
 						 const discrf_type * const gaussdens,
-						 int *NITEMS)
+                         int *NITEMS,
+                         wmh_type& wmh,
+                         const float * Rrad)
 {
 	int   jp;
 	float eff;
@@ -380,7 +382,7 @@ void wm_size_estimation (int kOS,
 				
 				//... to project voxels onto the detection plane and to calculate other distances .....
 				
-                voxel_projection( &vox , &eff , prj.lngcmd2 );
+                voxel_projection( &vox , &eff , prj.lngcmd2, wmh);
 				
 				//... setting PSF to zero	.........................................	
 				
@@ -391,13 +393,13 @@ void wm_size_estimation (int kOS,
 				
 				//... correction for PSF ..............................	
 				
-				if ( !wmh.do_psf  )	fill_psf_no ( &psf, &psf1d_h, vox, &ang[ ka ], bin.szdx );
+                if ( !wmh.do_psf  )	fill_psf_no ( &psf, &psf1d_h, vox, &ang[ ka ], bin.szdx, wmh);
 				
 				else{
 					
-					if ( wmh.do_psf_3d ) fill_psf_3d ( &psf, &psf1d_h, &psf1d_v, vox, gaussdens, bin.szdx, bin.thdx, bin.thcmd2 );
+                    if ( wmh.do_psf_3d ) fill_psf_3d ( &psf, &psf1d_h, &psf1d_v, vox, gaussdens, bin.szdx, bin.thdx, bin.thcmd2, wmh );
 					
-					else fill_psf_2d ( &psf, &psf1d_h, vox, gaussdens, bin.szdx );
+                    else fill_psf_2d ( &psf, &psf1d_h, vox, gaussdens, bin.szdx, wmh);
 				}
 
 				
@@ -526,7 +528,7 @@ void calc_vxprj( angle_type *ang )
 //=== voxel_projection =====================================================
 //==========================================================================
 
-void voxel_projection ( voxel_type *vox, float * eff, float lngcmd2)
+void voxel_projection ( voxel_type *vox, float * eff, float lngcmd2, wmh_type &wmh)
 {
 	
 	if ( wmh.COL.do_fb ){				// fan_beam
@@ -557,7 +559,12 @@ void voxel_projection ( voxel_type *vox, float * eff, float lngcmd2)
 //=== fill_psf_no ==========================================================
 //==========================================================================
 
-void fill_psf_no( psf2da_type *psf, psf1d_type * psf1d_h, const voxel_type& vox, angle_type const *const ang , float szdx )
+void fill_psf_no( psf2da_type *psf,
+                  psf1d_type * psf1d_h,
+                  const voxel_type& vox,
+                  angle_type const *const ang ,
+                  float szdx,
+                  wmh_type &wmh)
 {
 	psf1d_h->sgmcm   = vox.szcm;
 
@@ -571,7 +578,7 @@ void fill_psf_no( psf2da_type *psf, psf1d_type * psf1d_h, const voxel_type& vox,
 	psf1d_h->lngcmd2 = psf1d_h->lngcm / (float)2.;
 	psf1d_h->efres   = ang->vxprj.res * psf1d_h->sgmcm;  // to resize discretization resolution once applied sgmcm
 	
-	calc_psf_bin( vox.xd0, wmh.prj.szcm, &ang->vxprj, psf1d_h );
+    calc_psf_bin( vox.xd0, wmh.prj.szcm, &ang->vxprj, psf1d_h, wmh);
     
     for ( int ie = 0 ; ie < psf1d_h->Nib ; ie++ ){
         
@@ -586,7 +593,7 @@ void fill_psf_no( psf2da_type *psf, psf1d_type * psf1d_h, const voxel_type& vox,
 //=== fill_psf_2d ==========================================================
 //==========================================================================
 
-void fill_psf_2d( psf2da_type *psf, psf1d_type * psf1d_h, const voxel_type& vox, discrf_type const* const gaussdens, float szdx )
+void fill_psf_2d( psf2da_type *psf, psf1d_type * psf1d_h, const voxel_type& vox, discrf_type const* const gaussdens, float szdx,wmh_type &wmh )
 {
  
     psf1d_h->sgmcm   = calc_sigma_h( vox, wmh.COL );
@@ -597,7 +604,7 @@ void fill_psf_2d( psf2da_type *psf, psf1d_type * psf1d_h, const voxel_type& vox,
 	
 	psf1d_h->efres   = gaussdens->res * psf1d_h->sgmcm ;
 	
-	calc_psf_bin( vox.xd0, wmh.prj.szcm, gaussdens, psf1d_h );
+    calc_psf_bin( vox.xd0, wmh.prj.szcm, gaussdens, psf1d_h, wmh);
     
     for ( int ie = 0 ; ie < psf1d_h->Nib ; ie++ ){
         
@@ -616,7 +623,10 @@ void fill_psf_2d( psf2da_type *psf, psf1d_type * psf1d_h, const voxel_type& vox,
 void fill_psf_3d (psf2da_type *psf,
                   psf1d_type *psf1d_h,
                   psf1d_type *psf1d_v,
-                  const voxel_type& vox, discrf_type const * const gaussdens, float szdx, float thdx, float thcmd2 )
+                  const voxel_type& vox,
+                  discrf_type const * const gaussdens,
+                  float szdx, float thdx, float thcmd2,
+                  SPECTUB::wmh_type &wmh)
 {
 	
 	//... horizontal component ...........................
@@ -636,11 +646,11 @@ void fill_psf_3d (psf2da_type *psf,
 	
 	//... calculation of the horizontal component of psf ...................
 	
-	calc_psf_bin( vox.xd0, wmh.prj.szcm, gaussdens, psf1d_h );
+    calc_psf_bin( vox.xd0, wmh.prj.szcm, gaussdens, psf1d_h, wmh);
 
 	//... vertical component ..............................
 	
-	psf1d_v->sgmcm   = calc_sigma_v( vox, wmh.COL);
+    psf1d_v->sgmcm   = calc_sigma_v( vox, wmh.COL, wmh);
 	psf1d_v->lngcmd2 = psf1d_v->sgmcm * wmh.maxsigm;
 	psf1d_v->lngcm   = psf1d_v->lngcmd2 * (float)2.;
 	psf1d_v->di      = min( (int) floor( thdx / psf1d_v->sgmcm ), gaussdens->lng - 1 ) ;
@@ -655,7 +665,7 @@ void fill_psf_3d (psf2da_type *psf,
 	
 	//... calculation of the vertical component of psf ....................
 	
-	calc_psf_bin( thcmd2, wmh.prj.thcm, gaussdens, psf1d_v );
+    calc_psf_bin( thcmd2, wmh.prj.thcm, gaussdens, psf1d_v, wmh);
 	
     //... mixing and setting PSF area to 1 (to correct for tail truncation of Gaussian function) .....
 	
@@ -699,7 +709,8 @@ void fill_psf_3d (psf2da_type *psf,
 void calc_psf_bin (float center_psf,
 				   float binszcm,
 				   discrf_type const * const vxprj,
-				   psf1d_type *psf)
+                   psf1d_type *psf,
+                   SPECTUB::wmh_type &wmh)
 {
 	float weight, preval;
 
@@ -926,7 +937,8 @@ int comp_dist( float dx,
 //=== cal_att =================================================================
 //=============================================================================
 
-float calc_att( const attpth_type *const attpth, const float *const attmap , int nsli ){
+float calc_att( const attpth_type *const attpth, const float *const attmap , int nsli,
+                wmh_type& wmh){
 	
 	float att_coef = (float)0.;
 	int iv;