diff --git a/documentation/release_6.4.htm b/documentation/release_6.4.htm
index efdddc2c09..6b918f9d43 100644
--- a/documentation/release_6.4.htm
+++ b/documentation/release_6.4.htm
@@ -40,6 +40,9 @@ <h2> Summary for end users (also to be read by developers)</h2>
   <h3>New functionality</h3>
   <h4>General</h4>
   <ul>
+    <li>
+      Two SAFIR scanners (Safir-I and Safir-II) with a Generic geometry were added to the list of scanners 
+      (see <a href=https://github.com/UCL/STIR/pull/1563>PR</a> for details).
   </ul>
   <h4>Python</h4>
   <ul>
@@ -47,6 +50,9 @@ <h4>Python</h4>
 
   <h4>Utilities</h4>
   <ul>
+    <li>
+      Two utilities were added to convert and trim sinograms in a generic geometry
+      (see <a href=https://github.com/UCL/STIR/pull/1563>PR</a> for details).
   </ul>
 
   <h3>Changed functionality</h3>
diff --git a/src/buildblock/Scanner.cxx b/src/buildblock/Scanner.cxx
index 6fbfe6ad9d..0da145d7ca 100644
--- a/src/buildblock/Scanner.cxx
+++ b/src/buildblock/Scanner.cxx
@@ -1146,6 +1146,70 @@ Scanner::Scanner(Type scanner_type)
       );
       break;
 
+      case SafirI:
+      set_params(SafirI, string_list("SafirI"),
+               24,    // num_rings_v
+               150,   // max_num_non_arccorrected_bins_v,
+               150,   // default_num_arccorrected_bins_v,
+               180,   // num_detectors_per_ring_v
+               64.05, // inner_ring_radius_v
+               5,     // average_depth_of_interaction_v
+               2.2,   // ring_spacing_v
+               1.1,   // bin_size_v
+               0,     // intrinsic_tilt_v
+               3,     // num_axial_blocks_per_bucket_v
+               1,     // num_transaxial_blocks_per_bucket_v
+               8,     // num_axial_crystals_per_block_v
+               15,    // num_transaxial_crystals_per_block_v
+               1,     // num_axial_crystals_per_singles_unit_v
+               1,     // num_transaxial_crystals_per_singles_unit_v
+               1,     // num_detector_layers_v
+               0.12,  // energy_resolution_v
+               511,   // reference_energy_v
+               (short int)1, //max_num_of_timing_poss_v,
+               0.F,   // size_timing_pos_v,
+               0.F,   // timing_resolution_v,
+               "",    // scanner_geometry_v
+               2.2,   // axial_crystal_spacing_v
+               2.2,   // transaxial_crystal_spacing_v
+               18.1,  // axial_block_spacing_v
+               33.6,  // transaxial_block_spacing_v
+               ""     // crystal_map_file_name_v
+              );
+      break;
+
+    case SafirII:
+      set_params(SafirII, string_list("SafirII"),
+               64,    // num_rings_v
+               150,   // max_num_non_arccorrected_bins_v,
+               150,   // default_num_arccorrected_bins_v,
+               180,   // num_detectors_per_ring_v
+               64.05, // inner_ring_radius_v
+               5,     // average_depth_of_interaction_v
+               2.2,   // ring_spacing_v
+               1.1,   // bin_size_v
+               0,     // intrinsic_tilt_v
+               8,     // num_axial_blocks_per_bucket_v
+               1,     // num_transaxial_blocks_per_bucket_v
+               8,     // num_axial_crystals_per_block_v
+               15,    // num_transaxial_crystals_per_block_v
+               1,     // num_axial_crystals_per_singles_unit_v
+               1,     // num_transaxial_crystals_per_singles_unit_v
+               1,     // num_detector_layers_v
+               0.12,  // energy_resolution_v
+               511,   // reference_energy_v
+               (short int)1, // max_num_of_timing_poss_v,
+               0.F,   // size_timing_pos_v,
+               0.F,   // timing_resolution_v,
+               "",    // scanner_geometry_v
+               2.2,   // axial_crystal_spacing_v
+               2.2,   // transaxial_crystal_spacing_v
+               18.1,  // axial_block_spacing_v
+               33.6,  // transaxial_block_spacing_v
+               ""     // crystal_map_file_name_v
+              );
+      break;
+
     case UPENN_5rings:
       set_params(UPENN_5rings,
                  string_list("UPENN_5rings"),
diff --git a/src/include/stir/Scanner.h b/src/include/stir/Scanner.h
index 5b2210c694..94759c66ca 100644
--- a/src/include/stir/Scanner.h
+++ b/src/include/stir/Scanner.h
@@ -175,6 +175,8 @@ class Scanner
     Allegro,
     GeminiTF,
     SAFIRDualRingPrototype,
+    SafirI,
+    SafirII,
     UPENN_5rings,
     UPENN_5rings_no_gaps,
     UPENN_6rings,
diff --git a/src/utilities/CMakeLists.txt b/src/utilities/CMakeLists.txt
index 05faf85500..960697f1e6 100644
--- a/src/utilities/CMakeLists.txt
+++ b/src/utilities/CMakeLists.txt
@@ -70,11 +70,14 @@ if (NOT MINI_STIR)
     find_ML_normfactors.cxx
     find_ML_singles_from_delayed.cxx
     find_normfactors_from_cylinder_data.cxx
+    find_normfactors_from_cylinder_data_SafirII.cxx
     find_recovery_coefficients_in_image_quality_phantom_nema_nu4.cxx
     write_sinogram_to_txt.cxx
     find_sum_projection_of_viewgram_and_sinogram.cxx
     separate_true_from_random_scatter_for_necr.cxx
     stir_timings.cxx
+    convert_projdata_types.cxx
+		trim_projdata.cxx
     )
 
   if (HAVE_ITK)
diff --git a/src/utilities/convert_projdata_types.cxx b/src/utilities/convert_projdata_types.cxx
new file mode 100644
index 0000000000..d254da003c
--- /dev/null
+++ b/src/utilities/convert_projdata_types.cxx
@@ -0,0 +1,95 @@
+//
+//
+/*!
+  \file
+  \ingroup utilities
+
+  \brief This program takes a projection data from one type and converts it to another type.
+  \author Parisa Khateri
+
+*/
+/*
+*/
+#include "stir/ProjData.h"
+#include "stir/IO/interfile.h"
+#include "stir/utilities.h"
+#include "stir/Bin.h"
+
+#include <fstream>
+#include <iostream>
+#include "stir/ProjDataFromStream.h"
+#include "stir/Viewgram.h"
+#include "stir/IO/read_from_file.h"
+#include "stir/SegmentByView.h"
+#include "stir/ProjDataInterfile.h"
+#include "stir/ProjDataInfo.h"
+#include "stir/LORCoordinates.h"
+
+#include "stir/GeometryBlocksOnCylindrical.h"
+#include "stir/DetectionPosition.h"
+#include "stir/CartesianCoordinate3D.h"
+#include "stir/DetectorCoordinateMap.h"
+#include <boost/make_shared.hpp>
+#include "stir/CPUTimer.h"
+#include "stir/shared_ptr.h"
+
+
+#ifndef STIR_NO_NAMESPACES
+using std::cerr;
+#endif
+
+USING_NAMESPACE_STIR
+
+
+
+int main(int argc, char *argv[])
+{
+  CPUTimer timer0;
+
+  if(argc<4)
+  {
+    cerr<<"Usage: " << argv[0] << " output_filename input_filename template_blk_projdata\n";
+    exit(EXIT_FAILURE);
+  }
+  std::string output_filename=argv[1];
+  shared_ptr<ProjData> in_pd_ptr = ProjData::read_from_file(argv[2]);
+  shared_ptr<ProjData> template_pd_ptr = ProjData::read_from_file(argv[3]);
+
+  shared_ptr<ProjDataInfo> in_pdi_ptr(in_pd_ptr->get_proj_data_info_sptr()->clone());
+  shared_ptr<ProjDataInfo> out_pdi_ptr(template_pd_ptr->get_proj_data_info_sptr()->clone());
+  ProjDataInterfile out_proj_data(template_pd_ptr->get_exam_info_sptr(), out_pdi_ptr, output_filename+".hs");
+  write_basic_interfile_PDFS_header(output_filename+".hs", out_proj_data);
+
+  timer0.start();
+
+  assert(in_pdi_ptr->get_min_segment_num()==-1*in_pdi_ptr->get_max_segment_num());
+  for (int seg=in_pdi_ptr->get_min_segment_num(); seg<=in_pdi_ptr->get_max_segment_num();++seg)
+  {
+    std::cout<<"seg_num = "<<seg<<"\n";
+    // keep sinograms out of the loop to avoid reallocations
+    // initialise to something because there's no default constructor
+    Viewgram<float> viewgram_blk = out_proj_data.get_empty_viewgram(out_proj_data.get_min_view_num(),seg);
+    Viewgram<float> viewgram_cyl = in_pd_ptr->get_empty_viewgram(in_pd_ptr->get_min_view_num(),seg);
+
+    for(int view=in_pdi_ptr->get_min_view_num(); view<=in_pdi_ptr->get_max_view_num();++view)
+    {
+      viewgram_blk = out_proj_data.get_empty_viewgram(view,seg);
+      viewgram_cyl = in_pd_ptr->get_viewgram(view,seg);
+
+      for(int ax=in_pdi_ptr->get_min_axial_pos_num(seg); ax<=in_pdi_ptr->get_max_axial_pos_num(seg);++ax)
+      {
+        for(int tang=in_pdi_ptr->get_min_tangential_pos_num(); tang<=in_pdi_ptr->get_max_tangential_pos_num(); ++tang)
+        {
+          viewgram_blk[ax][tang] = viewgram_cyl[ax][tang];
+        }
+      }
+      if (!(out_proj_data.set_viewgram(viewgram_blk)== Succeeded::yes))
+        warning("Error set_segment for projdata_symm %d\n", seg);
+    }
+  }
+
+  timer0.stop();
+  std::cerr << "\nConverting from cylindrical projdata to block took " << timer0.value() << "s CPU time.\n\n";
+
+  return EXIT_SUCCESS;
+}
diff --git a/src/utilities/find_normfactors_from_cylinder_data_SafirII.cxx b/src/utilities/find_normfactors_from_cylinder_data_SafirII.cxx
new file mode 100755
index 0000000000..43ef676cde
--- /dev/null
+++ b/src/utilities/find_normfactors_from_cylinder_data_SafirII.cxx
@@ -0,0 +1,363 @@
+/*
+   Copyright
+   */
+/*!
+
+  \file
+  \ingroup utilities
+
+  \brief Find normalisation factors given projection data of a cylinder (component-based method)
+
+  \author Parisa Khateri
+
+*/
+//#include "stir/Scanner.h"
+//#include "stir/stream.h"
+//#include "stir/CPUTimer.h"
+#include "stir/utilities.h"
+#include <iostream>
+#include <fstream>
+#include <string>
+#include <algorithm>
+#include "stir/shared_ptr.h"
+#include "stir/ProjData.h"
+#include "stir/CartesianCoordinate3D.h"
+#include "stir/LORCoordinates.h"
+#include "stir/Bin.h"
+#include "stir/ProjDataInterfile.h"
+#include "stir/IO/interfile.h"
+#include <limits>
+#include <tuple>
+
+USING_NAMESPACE_STIR
+
+
+std::tuple<int, int> get_crystalTangPos_from_LOR(int view_nr, int tang_pos) {
+	return {(view_nr + 180 + ((tang_pos - (abs(tang_pos) % 2)) / 2))%180, (view_nr + 270 - ((tang_pos + (abs(tang_pos) % 2)) / 2))%180};
+}
+
+std::tuple<int, int> get_crystalAxPos_from_LOR(int segment_nr, int ax_pos) {
+	if (segment_nr < 0)
+	{
+    	return  {ax_pos - segment_nr, ax_pos};
+	}
+	else
+	{
+    	return  {ax_pos, ax_pos + segment_nr};
+	}
+}
+
+
+int main(int argc, char **argv)
+{
+	if (argc!=4)
+	{//TODO ask the user to give the activity
+		std::cerr << "Usage: "<< argv[0]
+			<<" output_file_name_prefix cylider_measured_data cylinder_radius(mm)\n"
+			<<"only cylinder data are supported. The radius should be the radius of the measured cylinder data.\n"
+			<<"warning: mind the input order\n";
+		return EXIT_FAILURE;
+	}
+
+	shared_ptr<ProjData> cylinder_projdata_ptr = ProjData::read_from_file(argv[2]);
+	const std::string output_file_name = argv[1];
+	const float R = atof(argv[3]); // cylinder radius
+	if (R==0)
+	{
+		std::cerr << " Radius must be a float value\n"
+			<<"Usage: "<< argv[0]
+			<<" output_file_name_prefix cylider_measured_data cylinder_radius\n"
+			<<"warning: mind the input order\n";
+		return EXIT_FAILURE;
+	}
+
+	//output file
+	shared_ptr<ProjDataInfo> cylinder_pdi_ptr(cylinder_projdata_ptr->get_proj_data_info_sptr()->clone());
+
+	ProjDataInterfile output_projdata(cylinder_projdata_ptr->get_exam_info_sptr(), cylinder_pdi_ptr, output_file_name);
+	write_basic_interfile_PDFS_header(output_file_name, output_projdata);
+
+	CartesianCoordinate3D<float> c1, c2;
+	LORInAxialAndNoArcCorrSinogramCoordinates<float> lor;
+
+	// first find the average number of counts per LOR
+	std::cout << "Finding the Average NCounts per LOR \n";
+	float total_count = 0;
+	float min_count = std::numeric_limits<double>::max(); // minimum number of counts per LOR
+	float average_count = 0; //average number of counts per LOR in the active region
+	int num_active_LORs = 0; //number of LORs which pass through the cylinder
+
+	const int nTang_pos_max = 180;	
+	const int nAx_pos_max = 64;	
+
+	for (int seg = cylinder_projdata_ptr->get_min_segment_num(); seg <=cylinder_projdata_ptr->get_max_segment_num(); ++seg)
+	{
+		for (int view =cylinder_projdata_ptr->get_min_view_num(); view <=cylinder_projdata_ptr->get_max_view_num(); ++view)
+		{
+			Viewgram<float> cylinder_viewgram = cylinder_projdata_ptr->get_viewgram(view, seg);
+			
+			for (int ax = cylinder_projdata_ptr->get_min_axial_pos_num(seg); ax <=cylinder_projdata_ptr->get_max_axial_pos_num(seg); ++ax)
+			{
+				for (int tang = cylinder_projdata_ptr->get_min_tangential_pos_num(); tang <=cylinder_projdata_ptr->get_max_tangential_pos_num(); ++tang)
+				{
+					Bin bin(seg, view, ax, tang);
+					cylinder_projdata_ptr->get_proj_data_info_sptr()->get_LOR(lor, bin);
+					LORAs2Points<float> lor_as2points(lor);
+					LORAs2Points<float> intersection_coords;
+					if (find_LOR_intersections_with_cylinder(intersection_coords, lor_as2points, R) ==Succeeded::yes)
+					{ //this only succeeds if LOR is intersecting with the cylinder
+						float N_lor = cylinder_viewgram[ax][tang]; //counts seen by this lor
+						c1 = intersection_coords.p1();
+						c2 = intersection_coords.p2();
+						float c12 = sqrt( pow(c1.z()-c2.z(), 2) + pow(c1.y()-c2.y(), 2) + pow(c1.x()-c2.x(), 2) );  // length of intersection of lor with the cylinder
+						if (c12>0.5) // if LOR intersection is lager than 0.5 mm, check the count per LOR
+						{
+							float N_lor_corrected=N_lor/c12; // corrected for the length
+							total_count+=N_lor_corrected;
+							num_active_LORs+=1;
+							if (N_lor_corrected<min_count && N_lor_corrected!=0) min_count=N_lor_corrected;
+						}
+					}
+				}
+			}
+		}
+	}
+
+
+	average_count = total_count / num_active_LORs;
+	std::cout << "num_lor = "						<< num_active_LORs	<< "\n";
+	std::cout << "tot_count_per_length_unit = "		<< total_count		<< "\n";
+	std::cout << "average_count_per_length_unit = "	<< average_count	<< "\n";
+	std::cout << "non_zero_min_per_length_unit = "	<< min_count		<< "\n";
+
+	// Now find the Detector Efficiencies
+	std::array<std::array<float, nTang_pos_max>, nAx_pos_max> det_eff {{}};
+	std::array<std::array<float, nTang_pos_max>, nAx_pos_max> n_alive_channels {{}};
+
+	std::cout << "Finding the Detector Efficiencies \n";
+	for (int view =cylinder_projdata_ptr->get_min_view_num(); view <=cylinder_projdata_ptr->get_max_view_num(); ++view)
+	{
+		int seg = 0;
+		Viewgram<float> cylinder_viewgram = cylinder_projdata_ptr->get_viewgram(view, seg);
+		for (int ax = cylinder_projdata_ptr->get_min_axial_pos_num(seg); ax <=cylinder_projdata_ptr->get_max_axial_pos_num(seg); ++ax)
+		{
+			for (int tang = cylinder_projdata_ptr->get_min_tangential_pos_num(); tang <=cylinder_projdata_ptr->get_max_tangential_pos_num(); ++tang)
+			{
+				Bin bin(seg, view, ax, tang);
+				cylinder_projdata_ptr->get_proj_data_info_sptr()->get_LOR(lor, bin);
+				float N_lor = cylinder_viewgram[ax][tang]; //counts seen by this lor
+				if (N_lor > 1)
+				{
+					std::tuple<int, int> crystal_axial = get_crystalAxPos_from_LOR(seg, ax);
+					std::tuple<int, int> crystal_tang = get_crystalTangPos_from_LOR(view, tang);
+					det_eff[std::get<0>(crystal_axial)][std::get<0>(crystal_tang)] += N_lor;
+					det_eff[std::get<1>(crystal_axial)][std::get<1>(crystal_tang)] += N_lor;
+					n_alive_channels[std::get<0>(crystal_axial)][std::get<0>(crystal_tang)]++;
+					n_alive_channels[std::get<1>(crystal_axial)][std::get<1>(crystal_tang)]++;
+				}
+			}
+		}
+	}
+	
+	const int nSectors = 12;
+	float total_alive_channels = 0;
+	for (int i = 0; i < (nTang_pos_max/nSectors); i++) {
+		int dead_channels = 0;
+		float sum = 0;
+		for (int j = 0; j < nAx_pos_max; j++) {
+			for (int k = i; k < nTang_pos_max; k+=(nTang_pos_max/nSectors)) {
+				sum += det_eff[j][k];
+				total_alive_channels += n_alive_channels[j][k];
+			}
+		}
+		sum = sum / float(nSectors * nAx_pos_max);
+		total_alive_channels = total_alive_channels / float(nSectors * nAx_pos_max);
+		for (int j = 0; j < nAx_pos_max; j++) {
+			for (int k = i; k < nTang_pos_max; k+=(nTang_pos_max/nSectors)) {
+				float new_val = (det_eff[j][k] * total_alive_channels) / (sum * n_alive_channels[j][k]);
+				if (det_eff[j][k] < 10.0) {
+					std::cout << det_eff[j][k] << std::endl;
+					new_val = 0.00;
+				}
+				det_eff[j][k] = new_val;
+			}
+		}
+	}
+		
+	float mean_deteff = 0.0;
+	for (int j = 0; j < nAx_pos_max; j++) {
+		for (int k = 0; k < nTang_pos_max; k++) {
+			mean_deteff+=det_eff[j][j];
+		}
+	}
+	std::cout << "Mean Detector Efficiency" << mean_deteff/float(nAx_pos_max*nTang_pos_max) << std::endl;
+	
+	std::cout << "Apply detector efficiencies\n";
+	
+	total_count = 0;
+	min_count = std::numeric_limits<double>::max(); // minimum number of counts per LOR
+	average_count = 0; //average number of counts per LOR in the active region
+	num_active_LORs = 0; //number of LORs which pass through the cylinder
+
+	for (int seg = cylinder_projdata_ptr->get_min_segment_num(); seg <=cylinder_projdata_ptr->get_max_segment_num(); ++seg)
+	{
+		for (int view =cylinder_projdata_ptr->get_min_view_num(); view <=cylinder_projdata_ptr->get_max_view_num(); ++view)
+		{
+			Viewgram<float> cylinder_viewgram = cylinder_projdata_ptr->get_viewgram(view, seg);
+			
+			for (int ax = cylinder_projdata_ptr->get_min_axial_pos_num(seg); ax <=cylinder_projdata_ptr->get_max_axial_pos_num(seg); ++ax)
+			{
+				for (int tang = cylinder_projdata_ptr->get_min_tangential_pos_num(); tang <=cylinder_projdata_ptr->get_max_tangential_pos_num(); ++tang)
+				{
+					Bin bin(seg, view, ax, tang);
+					cylinder_projdata_ptr->get_proj_data_info_sptr()->get_LOR(lor, bin);
+					LORAs2Points<float> lor_as2points(lor);
+					LORAs2Points<float> intersection_coords;
+					if (find_LOR_intersections_with_cylinder(intersection_coords, lor_as2points, R) ==Succeeded::yes)
+					{ //this only succeeds if LOR is intersecting with the cylinder
+						float N_lor = cylinder_viewgram[ax][tang]; //counts seen by this lor
+						c1 = intersection_coords.p1();
+						c2 = intersection_coords.p2();
+						float c12 = sqrt( pow(c1.z()-c2.z(), 2) + pow(c1.y()-c2.y(), 2) + pow(c1.x()-c2.x(), 2) );  // length of intersection of lor with the cylinder
+						if (c12>0.5) // if LOR intersection is lager than 0.5 mm, check the count per LOR
+						{
+							std::tuple<int, int> crystal_axial = get_crystalAxPos_from_LOR(seg, ax);
+							std::tuple<int, int> crystal_tang = get_crystalTangPos_from_LOR(view, tang);
+							float N_lor_corrected=0.0; // corrected for the length and detector efficiencies
+							if (det_eff[std::get<0>(crystal_axial)][std::get<0>(crystal_tang)] * det_eff[std::get<1>(crystal_axial)][std::get<1>(crystal_tang)] != 0)
+							{
+								N_lor_corrected=N_lor/(c12 * det_eff[std::get<0>(crystal_axial)][std::get<0>(crystal_tang)] * det_eff[std::get<1>(crystal_axial)][std::get<1>(crystal_tang)]); // corrected for the length and detector efficiencies
+								num_active_LORs+=1;
+							}
+							total_count+=N_lor_corrected;
+							if (N_lor_corrected<min_count && N_lor_corrected!=0)
+							{
+								min_count=N_lor_corrected;
+							}
+						}
+					}
+				}
+			}
+		}
+	}
+	
+	average_count = total_count / num_active_LORs;
+	std::cout << "num_lor = "						<< num_active_LORs	<< "\n";
+	std::cout << "tot_count_per_length_unit = "		<< total_count		<< "\n";
+	std::cout << "average_count_per_length_unit = "	<< average_count	<< "\n";
+	std::cout << "non_zero_min_per_length_unit = "	<< min_count		<< "\n";
+
+	
+	
+	std::cout << "Get geometric coefficients\n";
+	const int nSeg_proj = 64;	
+	const int nView_proj = 90;	
+	const int nTang_pos_proj = 100;	
+	
+	if ((nTang_pos_proj-1)!=(cylinder_projdata_ptr->get_max_tangential_pos_num() - cylinder_projdata_ptr->get_min_tangential_pos_num()))
+	{//TODO ask the user to give the activity
+		std::cerr << "Error: Geometry Hardcoded for SAFIR-II. Check number of Tangential positions in code to projdata \n";
+		return EXIT_FAILURE;
+	}
+	if ((nSeg_proj-1)!=cylinder_projdata_ptr->get_max_segment_num())
+	{//TODO ask the user to give the activity
+		std::cerr << "Error: Geometry Hardcoded for SAFIR-II. Check number of segments in code to projdata \n";
+		return EXIT_FAILURE;
+	}
+
+	std::array<std::array<std::array<float, nTang_pos_proj>, nView_proj>, nSeg_proj> geom_coeff {{{}}};
+
+	for (int seg = 0; seg < nSeg_proj; ++seg)
+	{
+		for (int tang = -(nTang_pos_proj/2); tang < (nTang_pos_proj/2); ++tang)
+		{
+			for (int base_view = 0; base_view < (nView_proj / (nSectors / 2)); ++base_view)
+			{
+				num_active_LORs = 0;
+				float sum = 0;
+				for (int true_view = base_view; true_view < nView_proj; true_view+=(nView_proj / (nSectors / 2)))
+				{
+					int valid_segments[2] = {seg, -seg};
+					for (int true_seg : valid_segments)
+					{
+						Viewgram<float> cylinder_viewgram = cylinder_projdata_ptr->get_viewgram(true_view, true_seg);
+						for (int ax = cylinder_projdata_ptr->get_min_axial_pos_num(true_seg); ax <= cylinder_projdata_ptr->get_max_axial_pos_num(true_seg); ++ax)
+						{
+							Bin bin(true_seg, true_view, ax, tang);
+							cylinder_projdata_ptr->get_proj_data_info_sptr()->get_LOR(lor, bin);
+							LORAs2Points<float> lor_as2points(lor);
+							LORAs2Points<float> intersection_coords;
+
+							if (find_LOR_intersections_with_cylinder(intersection_coords, lor_as2points, R) ==Succeeded::yes)
+							{ //this only succeeds if LOR is intersecting with the cylinder
+
+								float N_lor = cylinder_viewgram[ax][tang]; //counts seen by this lor
+								c1 = intersection_coords.p1();
+								c2 = intersection_coords.p2();
+								float c12 = sqrt( pow(c1.z()-c2.z(), 2) + pow(c1.y()-c2.y(), 2) + pow(c1.x()-c2.x(), 2) );  // length of intersection of lor with the cylinder
+								
+								if (c12>0.5) // if LOR intersection is lager than 0.5 mm, check the count per LOR
+								{
+									if (N_lor>1) //if value is large enough
+									{
+										std::tuple<int, int> crystal_axial = get_crystalAxPos_from_LOR(true_seg, ax);
+										std::tuple<int, int> crystal_tang = get_crystalTangPos_from_LOR(true_view, tang);
+										float N_lor_corrected=0.0; // corrected for the length and detector efficiencies
+										if (det_eff[std::get<0>(crystal_axial)][std::get<0>(crystal_tang)] * det_eff[std::get<1>(crystal_axial)][std::get<1>(crystal_tang)] != 0)
+										{
+											float N_lor_corrected=N_lor/(c12 * det_eff[std::get<0>(crystal_axial)][std::get<0>(crystal_tang)] * det_eff[std::get<1>(crystal_axial)][std::get<1>(crystal_tang)]); // corrected for the length and detector efficiencies
+											num_active_LORs++;
+											sum += N_lor_corrected;
+										}
+									}
+								}
+							}
+						}
+					}
+				}
+
+				for (int true_view = base_view; true_view < nView_proj; true_view+=(nView_proj / (nSectors / 2)))
+				{
+					if (sum == 0)
+					{
+						geom_coeff[seg][true_view][tang+(nTang_pos_proj/2)] = 0.0000;
+					}
+					else
+					{
+						geom_coeff[seg][true_view][tang+(nTang_pos_proj/2)] = average_count / (sum / num_active_LORs);
+					}
+				}
+			}
+		}
+	}
+
+	std::cout << "Find Norm Factors per LOR\n";
+
+	for (int seg = cylinder_projdata_ptr->get_min_segment_num(); seg <= cylinder_projdata_ptr->get_max_segment_num(); ++seg)
+	{
+		for (int view =cylinder_projdata_ptr->get_min_view_num(); view <= cylinder_projdata_ptr->get_max_view_num(); ++view)
+		{
+			Viewgram<float> out_viewgram = cylinder_projdata_ptr->get_empty_viewgram(view, seg);
+			for (int ax = cylinder_projdata_ptr->get_min_axial_pos_num(seg); ax <= cylinder_projdata_ptr->get_max_axial_pos_num(seg); ++ax)
+			{
+				for (int tang = cylinder_projdata_ptr->get_min_tangential_pos_num(); tang <= cylinder_projdata_ptr->get_max_tangential_pos_num(); ++tang)
+				{
+					Bin bin(seg, view, ax, tang);
+					cylinder_projdata_ptr->get_proj_data_info_sptr()->get_LOR(lor, bin);
+
+					std::tuple<int, int> crystal_axial = get_crystalAxPos_from_LOR(seg, ax);
+					std::tuple<int, int> crystal_tang = get_crystalTangPos_from_LOR(view, tang);
+					if (geom_coeff[abs(seg)][view][tang+(nTang_pos_proj/2)] == 0 || det_eff[std::get<0>(crystal_axial)][std::get<0>(crystal_tang)] == 0 || det_eff[std::get<1>(crystal_axial)][std::get<1>(crystal_tang)] == 0)
+					{
+						out_viewgram[ax][tang] = 0.0;
+					}
+					else
+					{
+						out_viewgram[ax][tang] = geom_coeff[abs(seg)][view][tang+(nTang_pos_proj/2)] / (det_eff[std::get<0>(crystal_axial)][std::get<0>(crystal_tang)] * det_eff[std::get<1>(crystal_axial)][std::get<1>(crystal_tang)]);
+					}
+				}
+			}
+			output_projdata.set_viewgram(out_viewgram);
+		}
+	}
+}
diff --git a/src/utilities/trim_projdata.cxx b/src/utilities/trim_projdata.cxx
new file mode 100644
index 0000000000..4bac93a7f1
--- /dev/null
+++ b/src/utilities/trim_projdata.cxx
@@ -0,0 +1,253 @@
+//
+//
+/*
+    
+*/
+/*!
+
+  \file
+  \ingroup utilities
+  \brief Main program for trim_projdata
+
+  \author Parisa Khateri
+
+
+   \par Usage:
+   \code
+   trim_projdata [-t num_tang_poss_to_trim] \
+        output_filename input_projdata_name 
+   \endcode
+   \param num_tang_poss_to_trim has to be smaller than the available number
+      of tangential positions.
+
+  \par Example:
+  \code
+  trim_projdata -t 10 out in.hs
+  \endcode
+*/
+#include "stir/ProjData.h"
+#include "stir/shared_ptr.h"
+#include <string>
+#include <string.h>
+#include "stir/ProjDataFromStream.h"
+#include "stir/ProjDataInterfile.h"
+#include "stir/ProjDataInfoCylindrical.h"
+#include "stir/ProjDataInfoBlocksOnCylindrical.h"
+#include "stir/ProjDataInfoGeneric.h"
+#include "stir/Sinogram.h"
+#include "stir/Bin.h"
+#include "stir/round.h"
+#include <fstream>
+#include <algorithm>
+
+
+#ifndef STIR_NO_NAMESPACES
+using std::string;
+using std::cerr;
+#endif
+
+USING_NAMESPACE_STIR
+
+int main(int argc, char **argv)
+{
+  int num_tang_poss_to_trim = 0;
+  if (argc>1 && strcmp(argv[1], "-t")==0)
+    {
+      num_tang_poss_to_trim =  atoi(argv[2]);
+      argc -= 2; argv += 2;
+    }
+  if (argc > 5 || argc < 3 )
+    {
+      cerr << "Usage:\n"
+	   << argv[0] << " [-t num_tang_poss_to_trim] \\\n"
+	   << "\toutput_filename input_projdata_name \\\n"
+	   << "num_tang_poss_to_trim has to be smaller than the available number\n";
+      exit(EXIT_FAILURE);
+    }
+  const string  output_filename = argv[1];
+  shared_ptr<ProjData> in_projdata_ptr = ProjData::read_from_file(argv[2]);
+  
+  
+  if (in_projdata_ptr->get_num_tangential_poss() <=
+      num_tang_poss_to_trim)
+    error("trim_projdata: too large number of tangential positions to trim (%d)\n",
+    num_tang_poss_to_trim);
+  
+  if (in_projdata_ptr->get_proj_data_info_sptr()->get_scanner_ptr()->get_scanner_geometry() ==
+      "Cylindrical")
+  {
+    // const ProjDataInfoCylindrical * const in_projdata_info_cyl_ptr =
+      // dynamic_cast<ProjDataInfoCylindrical const * >(in_projdata_ptr->get_proj_data_info_sptr());
+    auto proj_data_info_ptr = in_projdata_ptr->get_proj_data_info_sptr();
+    auto in_projdata_info_cyl_ptr = dynamic_cast<ProjDataInfoCylindrical const * >(proj_data_info_ptr.get());
+    
+    if (in_projdata_info_cyl_ptr== NULL)
+      {
+        error("error converting to cylindrical projection data\n");
+      }
+    ProjDataInfoCylindrical * out_projdata_info_cyl_ptr = 
+      dynamic_cast<ProjDataInfoCylindrical * >
+      (in_projdata_info_cyl_ptr->clone());
+
+    out_projdata_info_cyl_ptr->
+      set_num_tangential_poss(in_projdata_info_cyl_ptr->get_num_tangential_poss() -
+            num_tang_poss_to_trim);
+    
+    shared_ptr<ProjDataInfo> out_projdata_info_ptr(out_projdata_info_cyl_ptr);
+    ProjDataInterfile out_projdata(in_projdata_ptr->get_exam_info_sptr(), 
+            out_projdata_info_ptr, output_filename, std::ios::out); 
+    
+    for (int seg = out_projdata.get_min_segment_num(); 
+             seg <= out_projdata.get_max_segment_num();
+             ++seg)
+    {
+      // keep sinograms out of the loop to avoid reallocations
+      // initialise to something because there's no default constructor
+      Sinogram<float> out_sino = 
+        out_projdata.get_empty_sinogram(out_projdata.get_min_axial_pos_num(seg),seg);
+      Sinogram<float> in_sino = 
+        in_projdata_ptr->get_empty_sinogram(in_projdata_ptr->get_min_axial_pos_num(seg),seg);
+
+      for (int ax = out_projdata.get_min_axial_pos_num(seg); 
+               ax <= out_projdata.get_max_axial_pos_num(seg);
+               ++ax )
+      {
+        out_sino= out_projdata.get_empty_sinogram(ax, seg);
+        in_sino = in_projdata_ptr->get_sinogram(ax, seg);
+        
+        {
+        for (int view=out_projdata.get_min_view_num();
+                 view <= out_projdata.get_max_view_num();
+                 ++view)
+          for (int tang=out_projdata.get_min_tangential_pos_num();
+                   tang <= out_projdata.get_max_tangential_pos_num();
+                   ++tang)
+            out_sino[view][tang] = in_sino[view][tang];
+        }
+        out_projdata.set_sinogram(out_sino);
+      }
+      
+    }
+      
+  }
+  else if (in_projdata_ptr->get_proj_data_info_sptr()->get_scanner_ptr()->get_scanner_geometry() ==
+            "BlocksOnCylindrical")
+  {
+    // const ProjDataInfoBlocksOnCylindrical * const in_projdata_info_blk_ptr =
+    //   dynamic_cast<ProjDataInfoBlocksOnCylindrical const * >(in_projdata_ptr->get_proj_data_info_sptr());
+    auto proj_data_info_ptr = in_projdata_ptr->get_proj_data_info_sptr();
+    auto in_projdata_info_blk_ptr = dynamic_cast<ProjDataInfoBlocksOnCylindrical const * >(proj_data_info_ptr.get());
+    if (in_projdata_info_blk_ptr== NULL)
+      {
+        error("error converting to BlocksOnCylindrical projection data\n");
+      }
+    ProjDataInfoBlocksOnCylindrical * out_projdata_info_blk_ptr = 
+      dynamic_cast<ProjDataInfoBlocksOnCylindrical * >
+      (in_projdata_info_blk_ptr->clone());
+
+    out_projdata_info_blk_ptr->
+      set_num_tangential_poss(in_projdata_info_blk_ptr->get_num_tangential_poss() -
+            num_tang_poss_to_trim);
+    
+    shared_ptr<ProjDataInfo> out_projdata_info_ptr(out_projdata_info_blk_ptr);
+    ProjDataInterfile out_projdata(in_projdata_ptr->get_exam_info_sptr(), 
+            out_projdata_info_ptr, output_filename, std::ios::out); 
+            
+  
+    for (int seg = out_projdata.get_min_segment_num(); 
+             seg <= out_projdata.get_max_segment_num();
+             ++seg)
+    {
+      // keep sinograms out of the loop to avoid reallocations
+      // initialise to something because there's no default constructor
+      Sinogram<float> out_sino = 
+        out_projdata.get_empty_sinogram(out_projdata.get_min_axial_pos_num(seg),seg);
+      Sinogram<float> in_sino = 
+        in_projdata_ptr->get_empty_sinogram(in_projdata_ptr->get_min_axial_pos_num(seg),seg);
+
+      for (int ax = out_projdata.get_min_axial_pos_num(seg); 
+               ax <= out_projdata.get_max_axial_pos_num(seg);
+               ++ax )
+      {
+        out_sino= out_projdata.get_empty_sinogram(ax, seg);
+        in_sino = in_projdata_ptr->get_sinogram(ax, seg);
+        
+        {
+        for (int view=out_projdata.get_min_view_num();
+                 view <= out_projdata.get_max_view_num();
+                 ++view)
+          for (int tang=out_projdata.get_min_tangential_pos_num();
+                   tang <= out_projdata.get_max_tangential_pos_num();
+                   ++tang)
+            out_sino[view][tang] = in_sino[view][tang];
+        }
+        out_projdata.set_sinogram(out_sino);
+      }
+      
+    }
+      
+
+  } 
+  else if (in_projdata_ptr->get_proj_data_info_sptr()->get_scanner_ptr()->get_scanner_geometry() ==
+            "Generic")
+  {
+    // const ProjDataInfoGeneric * const in_projdata_info_gen_ptr =
+    //   dynamic_cast<ProjDataInfoGeneric const * >(in_projdata_ptr->get_proj_data_info_sptr());
+    auto proj_data_info_ptr = in_projdata_ptr->get_proj_data_info_sptr();
+    auto in_projdata_info_gen_ptr = dynamic_cast<ProjDataInfoGeneric const * >(proj_data_info_ptr.get());
+    if (in_projdata_info_gen_ptr== NULL)
+      {
+        error("error converting to Generic projection data\n");
+      }
+    ProjDataInfoGeneric * out_projdata_info_gen_ptr = 
+      dynamic_cast<ProjDataInfoGeneric * >
+      (in_projdata_info_gen_ptr->clone());
+
+    out_projdata_info_gen_ptr->
+      set_num_tangential_poss(in_projdata_info_gen_ptr->get_num_tangential_poss() -
+            num_tang_poss_to_trim);
+    
+    shared_ptr<ProjDataInfo> out_projdata_info_ptr(out_projdata_info_gen_ptr);
+    ProjDataInterfile out_projdata(in_projdata_ptr->get_exam_info_sptr(), 
+            out_projdata_info_ptr, output_filename, std::ios::out); 
+            
+    for (int seg = out_projdata.get_min_segment_num(); 
+             seg <= out_projdata.get_max_segment_num();
+             ++seg)
+    {
+      // keep sinograms out of the loop to avoid reallocations
+      // initialise to something because there's no default constructor
+      Sinogram<float> out_sino = 
+        out_projdata.get_empty_sinogram(out_projdata.get_min_axial_pos_num(seg),seg);
+      Sinogram<float> in_sino = 
+        in_projdata_ptr->get_empty_sinogram(in_projdata_ptr->get_min_axial_pos_num(seg),seg);
+
+      for (int ax = out_projdata.get_min_axial_pos_num(seg); 
+               ax <= out_projdata.get_max_axial_pos_num(seg);
+               ++ax )
+      {
+        out_sino= out_projdata.get_empty_sinogram(ax, seg);
+        in_sino = in_projdata_ptr->get_sinogram(ax, seg);
+        
+        {
+        for (int view=out_projdata.get_min_view_num();
+                 view <= out_projdata.get_max_view_num();
+                 ++view)
+          for (int tang=out_projdata.get_min_tangential_pos_num();
+                   tang <= out_projdata.get_max_tangential_pos_num();
+                   ++tang)
+            out_sino[view][tang] = in_sino[view][tang];
+        }
+        out_projdata.set_sinogram(out_sino);
+      }
+      
+    }
+  
+  }
+  else
+  {
+    error("error the scanner geometry of projection data is not known\n");
+  }
+
+  return EXIT_SUCCESS;
+}