diff --git a/src/include/stir/recon_buildblock/ProjMatrixByBinSPECTUB.h b/src/include/stir/recon_buildblock/ProjMatrixByBinSPECTUB.h
index eae1ecd495..51022e6ada 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,
@@ -97,7 +97,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 +146,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 +171,9 @@ class ProjMatrixByBinSPECTUB :
 
   bool already_setup;
 
+  static wm_da_type wm;
+  static 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 +208,10 @@ class ProjMatrixByBinSPECTUB :
   int maxszb;
 
 	
-  void compute_one_subset(const int kOS) const;
+  void compute_one_subset(const int kOS,
+                          wm_da_type& wm,
+                          wmh_type& wmh,
+                          float * Rrad) const;
   void delete_UB_SPECT_arrays();
   mutable std::vector<bool> subset_already_processed;
 };
diff --git a/src/include/stir/recon_buildblock/SPECTUB_Weight3d.h b/src/include/stir/recon_buildblock/SPECTUB_Weight3d.h
index 5b5f28d88d..bdc0848a49 100644
--- a/src/include/stir/recon_buildblock/SPECTUB_Weight3d.h
+++ b/src/include/stir/recon_buildblock/SPECTUB_Weight3d.h
@@ -51,20 +51,25 @@ void wm_calculation( const int kOS,
 					const bool *msk_2d,
 					const int maxszb,
 					const discrf_type *const gaussdens,
-					const int *const  NITEMS
+                    const int *const  NITEMS,
+                    wm_da_type& wm,
+                    wmh_type& wmh,
+                    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 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,
+                         wmh_type &wmh,
+                         float *Rrad);
 
 
 //... geometric component ............................................
@@ -96,7 +101,7 @@ void calc_psf_bin ( float center_psf, float binszcm, discrf_type const * const v
 
 void calc_att_path ( const bin_type& bin, const voxel_type& vox, const volume_type& vol, attpth_type *attpth );
 
-float calc_att ( const attpth_type *const attpth, const float *const attmap, int islc );
+float calc_att (const attpth_type *const attpth, const float *const attmap, int islc , 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..0164629f17 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,6 +18,7 @@
   \author Berta Marti Fuster
   \author Carles Falcon
   \author Kris Thielemans
+  \author Daniel Deidda
 */
 
 //#include "stir/ProjDataInterfile.h"
@@ -60,14 +62,14 @@
 #include "stir/recon_buildblock/SPECTUB_Weight3d.h"
 
 /* UB-SPECT global variables */
-namespace SPECTUB {
-  wm_da_type wm;
-  wmh_type wmh; 
-  float * Rrad;
-}
+//namespace SPECTUB {
+//  wm_da_type wm;
+//  wmh_type wmh;
+//  float * Rrad;
+//}
 
 START_NAMESPACE_STIR
-
+using namespace SPECTUB;
 
 const char * const 
 ProjMatrixByBinSPECTUB::registered_name =
@@ -658,7 +660,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;
 
@@ -734,7 +737,10 @@ delete_UB_SPECT_arrays()
 }
 void
 ProjMatrixByBinSPECTUB::
-compute_one_subset(const int kOS) const
+compute_one_subset(const int kOS,
+                   wm_da_type& wm,
+                   wmh_type& wmh,
+                   float * Rrad) const
 {
   using namespace SPECTUB;
 
@@ -807,7 +813,8 @@ compute_one_subset(const int kOS) const
 
   //... 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],
+                   wm, wmh, Rrad );
   info(boost::format("Weight matrix calculation done. time %1% (s)") % timer.value(),
        2);
 
@@ -868,8 +875,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,wm,wmh,Rrad);
       subset_already_processed[kOS]=true;
     }
   lor.erase();
diff --git a/src/recon_buildblock/SPECTUB_Weight3d.cxx b/src/recon_buildblock/SPECTUB_Weight3d.cxx
index 1348cf743d..aeb6195d5d 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, float *Rrad)
 {
 	
 	float weight;
@@ -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,
+                         float * Rrad)
 {
 	int   jp;
 	float eff;
@@ -926,7 +928,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;