diff --git a/cosipy/response/FullDetectorResponse.py b/cosipy/response/FullDetectorResponse.py
index 2c4f1d36..e4a10050 100644
--- a/cosipy/response/FullDetectorResponse.py
+++ b/cosipy/response/FullDetectorResponse.py
@@ -198,7 +198,7 @@ def _open_rsp(cls, filename, Spectrumfile=None,norm="Linear" ,single_pixel = Fal
                     try :
                         norm = str(line[1])
                     except :
-                        print(f"norm not found in the file ! We assume {norm}")
+                        logger.info(f"norm not found in the file ! We assume {norm}")
 
                     if norm =="Linear" :
                         emin = int(line[2])
@@ -255,12 +255,12 @@ def _open_rsp(cls, filename, Spectrumfile=None,norm="Linear" ,single_pixel = Fal
         assert nevents_sim != 0,"number of simulated events is 0 !" 
         
         
-        print("normalisation is {0}".format(norm))
+        logger.info("normalisation is {0}".format(norm))
         if sparse == None :
-            print("Sparse paramater not found in the file : We assume this is a non sparse matrice !")
+            logger.info("Sparse paramater not found in the file : We assume this is a non sparse matrice !")
             sparse = False
         else :
-            print("Sparse matrice ? {0}".format(sparse))
+            logger.info("Sparse matrice ? {0}".format(sparse))
         edges = ()
 
         for axis_edges, axis_type in zip(axes_edges, axes_types):
@@ -298,8 +298,8 @@ def _open_rsp(cls, filename, Spectrumfile=None,norm="Linear" ,single_pixel = Fal
 
             # If fast method fails, use long method, which should work in all cases.
             except:
-                print("Initial attempt failed.")
-                print("Using long method...")
+                logger.info("Initial attempt failed.")
+                logger.info("Using long method...")
                 nlines = sum(1 for _ in gzip.open(filename,"rt"))
                 
             # Preallocate arrays
@@ -308,7 +308,7 @@ def _open_rsp(cls, filename, Spectrumfile=None,norm="Linear" ,single_pixel = Fal
 
             # Calculate the memory usage in Gigabytes
             memory_size = ((nlines * data.itemsize)+(axes.ndim*nlines*coords.itemsize))/(1024*1024*1024)
-            print(f"Estimated RAM you need to read the file : {memory_size} GB")
+            logger.info(f"Estimated RAM you need to read the file : {memory_size} GB")
 
     
                 
@@ -320,7 +320,7 @@ def _open_rsp(cls, filename, Spectrumfile=None,norm="Linear" ,single_pixel = Fal
 
             # Calculate the memory usage in Gigabytes
             memory_size = (nlines * data.itemsize)/(1024*1024*1024)
-            print(f"Estimated RAM you need to read the file : {memory_size} GB")
+            logger.info(f"Estimated RAM you need to read the file : {memory_size} GB")
 
 
         # Loop
@@ -381,7 +381,7 @@ def _open_rsp(cls, filename, Spectrumfile=None,norm="Linear" ,single_pixel = Fal
                     if binLine :
                         #check we have same number of bin than values read
                         if len(line)!=nbins :
-                            print("nb of bin content read ({0}) != nb of bins {1}".format(len(line),nbins))
+                            logger.info("nb of bin content read ({0}) != nb of bins {1}".format(len(line),nbins))
                             sys.exit()
                         
                         for i in tqdm(range(nbins), desc="Processing", unit="bin"):
@@ -393,7 +393,7 @@ def _open_rsp(cls, filename, Spectrumfile=None,norm="Linear" ,single_pixel = Fal
 
                         break
         
-        print("response file read ! Now we create the histogram and weight in order to "+ 
+        logger.info("response file read ! Now we create the histogram and weight in order to "+ 
                 "get the effective area")
         # create histpy histogram
 
@@ -424,11 +424,11 @@ def _open_rsp(cls, filename, Spectrumfile=None,norm="Linear" ,single_pixel = Fal
             gauss_int = 0.5 * (1 + erf( (edges[0]-Gauss_mean)/(4*np.sqrt(2)) ) ) + 0.5 * (1 + erf( (edges[1]-Gauss_mean)/(4*np.sqrt(2)) ) )
             
             assert gauss_int == 1, "The gaussian spectrum is not fully contained in this single bin !"
-            print("Only one bin so we will use the Mono normalisation")
+            logger.info("Only one bin so we will use the Mono normalisation")
             norm ="Mono"
 
         if Spectrumfile is not None and norm=="file":
-            print("normalisation : spectrum file")
+            logger.info("normalisation : spectrum file")
             # From spectrum file
             spec = pd.read_csv(Spectrumfile, sep=" ")
             spec = spec.iloc[:-1]
@@ -443,7 +443,7 @@ def _open_rsp(cls, filename, Spectrumfile=None,norm="Linear" ,single_pixel = Fal
             nperchannel_norm = hspec[:]
 
         elif norm=="powerlaw":
-            print("normalisation : powerlaw with index {0} with energy range [{1}-{2}]keV".format(alpha,emin,emax))
+            logger.info("normalisation : powerlaw with index {0} with energy range [{1}-{2}]keV".format(alpha,emin,emax))
             # From powerlaw
 
             e_lo = dr.axes['Ei'].lower_bounds
@@ -466,11 +466,11 @@ def _open_rsp(cls, filename, Spectrumfile=None,norm="Linear" ,single_pixel = Fal
                 nperchannel_norm = (e_hi**a - e_lo**a) / (emax**a - emin**a)            
 
         elif norm =="Linear" :
-            print("normalisation : linear with energy range [{0}-{1}]".format(emin,emax))
+            logger.info("normalisation : linear with energy range [{0}-{1}]".format(emin,emax))
             nperchannel_norm = ewidth / (emax-emin)
             
         elif norm=="Mono" :
-            print("normalisation : mono")
+            logger.info("normalisation : mono")
 
             nperchannel_norm = np.array([1.])
         
@@ -1163,7 +1163,7 @@ def command_dump():
                 apar.error(f"Argument '{option}' not valid for 'dump' command")
 
             if args.output is None:
-                print(result)
+                logger.info(result)
             else:
                 logger.info(f"Saving result to {Path(args.output).resolve()}")
                 f = open(args.output, 'a')
diff --git a/cosipy/spacecraftfile/SpacecraftFile.py b/cosipy/spacecraftfile/SpacecraftFile.py
index 914ad94d..ce2e3d69 100644
--- a/cosipy/spacecraftfile/SpacecraftFile.py
+++ b/cosipy/spacecraftfile/SpacecraftFile.py
@@ -956,7 +956,7 @@ def plot_rmf(self, file_name = None, save_name = None, dpi = 300):
 
         # Read the MATRIX extension
         matrix_ext = self.rmf['MATRIX']
-        #print(repr(matrix_hdu.header[:60]))
+        #logger.info(repr(matrix_hdu.header[:60]))
         energy_low = matrix_ext.data["ENERG_LO"] # energy bin lower edges for measured energies
         energy_high = matrix_ext.data["ENERG_HI"] # energy bin higher edges for measured energies
         data = matrix_ext.data
@@ -988,13 +988,13 @@ def plot_rmf(self, file_name = None, save_name = None, dpi = 300):
         energy_all_edges = np.append(energy_low,energy_high[-1])
         #bin_edges = np.array([incident_energy_bins,incident_energy_bins]) # doesn't work
         bin_edges = np.vstack((energy_all_edges, energy_all_edges))
-        #print(bin_edges)
+        #logger.info(bin_edges)
 
         self.probability = []
         for i in np.arange(10):
             for j in np.arange(10):
                 self.probability.append(rmf_matrix[i][j])
-        #print(type(probability))
+        #logger.info(type(probability))
 
         plt.hist2d(x=x_center_coords,y=y_center_coords,weights=self.probability,bins=bin_edges, norm=LogNorm())
         plt.xscale('log')
diff --git a/cosipy/test_data/20280301_2s.ori b/cosipy/test_data/20280301_2s.ori
new file mode 100644
index 00000000..1d95981f
--- /dev/null
+++ b/cosipy/test_data/20280301_2s.ori
@@ -0,0 +1,5 @@
+Type OrientationsGalactic 
+OG 1835487300.0 68.44034002307066 44.61117227945379 -21.559659976929343 44.61117227945379
+OG 1835487301.0 68.38920658776064 44.642124027021296 -21.610793412239364 44.642124027021296
+OG 1835487302.0 68.3380787943012 44.67309722321497 -21.661921205698793 44.67309722321497
+EN
\ No newline at end of file
diff --git a/cosipy/test_data/ei_cds_array.npy b/cosipy/test_data/ei_cds_array.npy
new file mode 100644
index 00000000..731157ba
Binary files /dev/null and b/cosipy/test_data/ei_cds_array.npy differ
diff --git a/cosipy/test_data/ts_map_bkg.h5 b/cosipy/test_data/ts_map_bkg.h5
new file mode 100644
index 00000000..b1283704
Binary files /dev/null and b/cosipy/test_data/ts_map_bkg.h5 differ
diff --git a/cosipy/test_data/ts_map_signal.h5 b/cosipy/test_data/ts_map_signal.h5
new file mode 100644
index 00000000..a19b9009
Binary files /dev/null and b/cosipy/test_data/ts_map_signal.h5 differ
diff --git a/cosipy/test_data/ts_map_src_bkg.h5 b/cosipy/test_data/ts_map_src_bkg.h5
new file mode 100644
index 00000000..ae8e5c81
Binary files /dev/null and b/cosipy/test_data/ts_map_src_bkg.h5 differ
diff --git a/cosipy/ts_map/TSMap.py b/cosipy/ts_map/TSMap.py
index 8048b04d..c5017aa5 100644
--- a/cosipy/ts_map/TSMap.py
+++ b/cosipy/ts_map/TSMap.py
@@ -12,6 +12,8 @@
 
 import astropy.io.fits as fits
 
+import logging
+logger = logging.getLogger(__name__)
 
 class TSMap:
 
@@ -175,8 +177,8 @@ def ts_grid_data(self):
             for j in range(self.log_like.axes['dec'].nbins):
         
                 # progress
-                print(f"\rra = {i:2d}/{self.log_like.axes['ra'].nbins}   ", end = "")
-                print(f"dec = {j:2d}/{self.log_like.axes['dec'].nbins}   ", end = "")
+                logger.info(f"\rra = {i:2d}/{self.log_like.axes['ra'].nbins}   ", end = "")
+                logger.info(f"dec = {j:2d}/{self.log_like.axes['dec'].nbins}   ", end = "")
         
                 # changing the position parameters
                 # converting rad to deg due to ra and dec in 3ML PointSource
@@ -237,10 +239,10 @@ def print_best_fit(self):
         else:
             self.best_ra = (self.ts.axes['ra'].centers[self.argmax[0]]) * (180/np.pi) # deg
         self.best_dec = self.ts.axes['dec'].centers[self.argmax[1]] * (180/np.pi) # deg
-        print(f"Best fit position: RA = {self.best_ra} deg, Dec = {self.best_dec} deg")
+        logger.info(f"Best fit position: RA = {self.best_ra} deg, Dec = {self.best_dec} deg")
         
         # convert to significance based on Wilk's theorem
-        print(f"Expected significance: {stats.norm.isf(stats.chi2.sf(self.ts_max, df = 2)):.1f} sigma")
+        logger.info(f"Expected significance: {stats.norm.isf(stats.chi2.sf(self.ts_max, df = 2)):.1f} sigma")
         
     def save_ts(self, output_file_name):
 
diff --git a/cosipy/ts_map/__init__.py b/cosipy/ts_map/__init__.py
index 129fa2c2..03f650a1 100644
--- a/cosipy/ts_map/__init__.py
+++ b/cosipy/ts_map/__init__.py
@@ -1,2 +1,3 @@
 from .TSMap import TSMap
-from .fast_ts_fit import FastTSMap
\ No newline at end of file
+from .fast_ts_fit import FastTSMap
+from .fast_norm_fit import FastNormFit
\ No newline at end of file
diff --git a/cosipy/ts_map/fast_ts_fit.py b/cosipy/ts_map/fast_ts_fit.py
index 9ed71b6f..323f7003 100644
--- a/cosipy/ts_map/fast_ts_fit.py
+++ b/cosipy/ts_map/fast_ts_fit.py
@@ -19,6 +19,9 @@
 import gc
 import matplotlib.pyplot as plt
 
+import logging
+logger = logging.getLogger(__name__)
+
 class FastTSMap():
     
     def __init__(self, data, bkg_model, response_path, orientation = None, cds_frame = "local", scheme = "RING"):
@@ -236,14 +239,14 @@ def get_ei_cds_array(hypothesis_coord, energy_channel, response_path, spectrum,
                                                                         quiet = True)
             #time_coord_convert_end = time.time()
             #time_coord_convert_used = time_coord_convert_end - time_coord_convert_start
-            #print(f"The time used for coordinate conversion is {time_coord_convert_used}s.")
+            #logger.info(f"The time used for coordinate conversion is {time_coord_convert_used}s.")
 
             #time_dwell_start = time.time()
             # get the dwell time map: the map of the time spent on each pixel in the local frame
             dwell_time_map = orientation.get_dwell_map(response = response_path)
             #time_dwell_end = time.time()
             #time_dwell_used = time_dwell_end - time_dwell_start
-            #print(f"The time used for dwell time map is {time_dwell_used}s.")
+            #logger.info(f"The time used for dwell time map is {time_dwell_used}s.")
             
             #time_psr_start = time.time()
             # convolve the response with the dwell_time_map to get the point source response
@@ -251,7 +254,7 @@ def get_ei_cds_array(hypothesis_coord, energy_channel, response_path, spectrum,
                 psr = response.get_point_source_response(dwell_time_map)
             #time_psr_end = time.time()
             #time_psr_used = time_psr_end - time_psr_start
-            #print(f"The time used for psr is {time_psr_used}s.")
+            #logger.info(f"The time used for psr is {time_psr_used}s.")
 
         elif cds_frame == "galactic":
             
@@ -272,7 +275,7 @@ def get_ei_cds_array(hypothesis_coord, energy_channel, response_path, spectrum,
 
         #time_cds_end = time.time()
         #time_cds_used = time_cds_end - time_cds_start
-        #print(f"The time used for cds is {time_cds_used}s.")
+        #logger.info(f"The time used for cds is {time_cds_used}s.")
         
         return ei_cds_array
     
@@ -304,7 +307,7 @@ def fast_ts_fit(hypothesis_coord,
         cds_frame : str
             "local" or "galactic", it's the Compton data space (CDS) frame of the data, bkg_model and the response. In other words, they should have the same cds frame .
         ts_nside : int
-            The nside of the ts map.
+            The nside of the ts map.   
         ts_scheme : str
             The scheme of the Ts map.
 
@@ -315,11 +318,9 @@ def fast_ts_fit(hypothesis_coord,
         """
         
         start_fast_ts_fit = time.time()
-        
-        # get the pix number of the ts map
-        data_array = np.zeros(hp.nside2npix(ts_nside))
-        ts_temp = HealpixMap(data = data_array, scheme = ts_scheme, coordsys = "galactic")
-        pix = ts_temp.ang2pix(hypothesis_coord)
+
+        # get the indices of the pixels to fit
+        pix = hp.ang2pix(nside = ts_nside, theta = hypothesis_coord.l.deg, phi = hypothesis_coord.b.deg, lonlat = True)
         
         # get the expected counts in the flattened cds array
         start_ei_cds_array = time.time()
@@ -342,7 +343,7 @@ def fast_ts_fit(hypothesis_coord,
         return [pix, result[0], result[1], result[2], result[3], result[4], time_ei_cds_array, time_fit, time_fast_ts_fit]
 
         
-    def parallel_ts_fit(self, hypothesis_coords, energy_channel, spectrum, ts_scheme = "RING", start_method = "fork", cpu_cores = None):
+    def parallel_ts_fit(self, hypothesis_coords, energy_channel, spectrum, ts_scheme = "RING", start_method = "fork", cpu_cores = None, ts_nside = None):
         
         """
         Perform parallel computation on all the hypothesis coordinates.
@@ -350,7 +351,7 @@ def parallel_ts_fit(self, hypothesis_coords, energy_channel, spectrum, ts_scheme
         Parameters
         ----------
         hypothesis_coords : list
-            A list of the hypothesis coordinates
+            A list of the hypothesis coordinates to fit
         energy_channel : list
             the energy channel you want to use: [lower_channel, upper_channel]. lower_channel is inclusive while upper_channel is exclusive.
         spectrum : astromodels.functions
@@ -361,6 +362,8 @@ def parallel_ts_fit(self, hypothesis_coords, energy_channel, spectrum, ts_scheme
             The starting method of the parallel computation (the default is "fork", which implies using the fork method to start parallel computation).
         cpu_cores : int, optional
             The number of cpu cores you wish to use for the parallel computation (the default is `None`, which implies using all the available number of cores -1 to perform the parallel computation).
+        ts_nside : int, optional
+            The nside of the ts map. This must be given if the number of hypothesis_coords isn't equal to the number of pixels of the total ts map, which means that you fit only a portion of the total ts map. (the default is `None`, which means that you fit the full ts map). 
         
         Returns
         -------
@@ -368,8 +371,9 @@ def parallel_ts_fit(self, hypothesis_coords, energy_channel, spectrum, ts_scheme
             The result of the ts fit over all the hypothesis coordinates.
         """
         
-        # decide the ts_nside from the list of hypothesis coordinates
-        ts_nside = hp.npix2nside(len(hypothesis_coords))
+        # decide the ts_nside from the list of hypothesis coordinates if not given
+        if ts_nside == None:
+            ts_nside = hp.npix2nside(len(hypothesis_coords))
         
         # get the flattened data_cds_array
         data_cds_array = FastTSMap.get_cds_array(self._data, energy_channel).flatten()
@@ -390,10 +394,10 @@ def parallel_ts_fit(self, hypothesis_coords, energy_channel, spectrum, ts_scheme
             # if you don't specify the number of cpu cores to use or the specified number of cpu cores is the same as the total number of cores you have
             # it will use the [total_cores - 1] number of cores to run the parallel computation.
             cores = total_cores - 1
-            print(f"You have total {total_cores} CPU cores, using {cores} CPU cores for parallel computation.")
+            logger.info(f"You have total {total_cores} CPU cores, using {cores} CPU cores for parallel computation.")
         else:
             cores = cpu_cores
-            print(f"You have total {total_cores} CPU cores, using {cores} CPU cores for parallel computation.")
+            logger.info(f"You have total {total_cores} CPU cores, using {cores} CPU cores for parallel computation.")
 
         start = time.time()
         multiprocessing.set_start_method(start_method, force = True)
@@ -409,7 +413,7 @@ def parallel_ts_fit(self, hypothesis_coords, energy_channel, spectrum, ts_scheme
         
         elapsed_seconds = end - start
         elapsed_minutes = elapsed_seconds/60
-        print(f"The time used for the parallel TS map computation is {elapsed_minutes} minutes")
+        logger.info(f"The time used for the parallel TS map computation is {elapsed_minutes} minutes")
         
         results = np.array(results)  # turn to a numpy array
         results = results[results[:, 0].argsort()]  # arrange the order by the pixel numbering
@@ -447,6 +451,7 @@ def _plot_ts(ts_array, skycoord = None, containment = None, save_plot = False, s
             lon = skycoord.l.deg
             lat = skycoord.b.deg
 
+
         # get the ts value
         m_ts = ts_array
         
@@ -462,14 +467,14 @@ def _plot_ts(ts_array, skycoord = None, containment = None, save_plot = False, s
             hp.mollview(m_ts[:], max = max_ts, min = max_ts-critical, title = f"Containment {percentage}%", coord = "G", hold = True) 
         elif containment == None:
             hp.mollview(m_ts[:], coord = "G", hold = True) 
-            
-        
+
         if skycoord != None:
             hp.projscatter(lon, lat, marker = "x", linewidths = 0.5, lonlat=True, coord = "G", label = f"True location at l={lon}, b={lat}", color = "fuchsia")
         hp.projscatter(0, 0, marker = "o", linewidths = 0.5, lonlat=True, coord = "G", color = "red")
         hp.projtext(350, 0, "(l=0, b=0)", lonlat=True, coord = "G", color = "red")
 
         if save_plot == True:
+
             fig.savefig(Path(save_dir)/save_name, dpi = dpi)
 
         return
@@ -498,6 +503,7 @@ def plot_ts(self, ts_array = None, skycoord = None, containment = None, save_plo
         if ts_array is not None:
 
             FastTSMap._plot_ts(ts_array = ts_array, skycoord = skycoord, containment = containment, 
+
                                save_plot = save_plot, save_dir = save_dir, save_name = save_name, dpi = dpi)
 
         else:
@@ -533,6 +539,6 @@ def show_memory_info(hint):
     
         info = p.memory_full_info()
         memory = info.uss / 1024. / 1024
-        print('{} memory used: {} MB'.format(hint, memory))
+        logger.info('{} memory used: {} MB'.format(hint, memory))
         
     
\ No newline at end of file
diff --git a/tests/spacecraftfile/__init__.py b/tests/spacecraftfile/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/tests/spacecraftfile/test_spacecraftfile.py b/tests/spacecraftfile/test_spacecraftfile.py
new file mode 100644
index 00000000..f4e78c60
--- /dev/null
+++ b/tests/spacecraftfile/test_spacecraftfile.py
@@ -0,0 +1,425 @@
+from cosipy import test_data
+from pytest import approx
+from cosipy import SpacecraftFile
+import numpy as np
+import astropy.units as u
+from astropy.coordinates import SkyCoord
+from astropy.io import fits
+import os
+from pathlib import Path
+from astropy.time import Time 
+
+def test_get_time():
+
+    ori_path = test_data.path / "20280301_first_10sec.ori"
+    
+    ori = SpacecraftFile.parse_from_file(ori_path)
+    
+    assert np.allclose(ori.get_time().value,
+                       [1835478000.0, 1835478001.0, 1835478002.0,
+                        1835478003.0, 1835478004.0, 1835478005.0,
+                        1835478006.0, 1835478007.0, 1835478008.0,
+                        1835478009.0, 1835478010.0])
+
+
+def test_get_time_delta():
+
+    ori_path = test_data.path / "20280301_first_10sec.ori"
+    ori = SpacecraftFile.parse_from_file(ori_path)
+    time_delta = ori.get_time_delta()
+    time_delta.format = "sec"
+
+    assert np.allclose(time_delta.value, np.array([1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 
+                                                   1.000000, 1.000000, 1.000000, 1.000000, 1.000000]))
+
+
+
+def test_get_attitude():
+
+    ori_path = test_data.path / "20280301_first_10sec.ori"
+    ori = SpacecraftFile.parse_from_file(ori_path)
+    
+    attitude = ori.get_attitude()
+
+    matrix = np.array([[[0.215904, -0.667290, -0.712818],
+                        [0.193436, 0.744798, -0.638638],
+                        [0.957062, 0.000000, 0.289883]],
+                       
+                       [[0.216493, -0.667602, -0.712347],
+                        [0.194127, 0.744518, -0.638754],
+                        [0.956789, 0.000000, 0.290783]],
+                       
+                       [[0.217081, -0.667914, -0.711875],
+                        [0.194819, 0.744238, -0.638870],
+                        [0.956515, -0.000000, 0.291683]],
+                       
+                       [[0.217669, -0.668227, -0.711402],
+                        [0.195511, 0.743958, -0.638985],
+                        [0.956240, 0.000000, 0.292582]],
+                       
+                       [[0.218255, -0.668539, -0.710929],
+                        [0.196204, 0.743677, -0.639100],
+                        [0.955965, 0.000000, 0.293481]],
+                       
+                       [[0.218841, -0.668852, -0.710455],
+                        [0.196897, 0.743396, -0.639214],
+                        [0.955688, -0.000000, 0.294380]],
+                       
+                       [[0.219426, -0.669165, -0.709980],
+                        [0.197590, 0.743114, -0.639327],
+                        [0.955411, 0.000000, 0.295279]],
+                       
+                       [[0.220010, -0.669477, -0.709504],
+                        [0.198284, 0.742833, -0.639440],
+                        [0.955133, -0.000000, 0.296177]],
+                       
+                       [[0.220594, -0.669790, -0.709027],
+                        [0.198978, 0.742551, -0.639552],
+                        [0.954854, 0.000000, 0.297075]],
+                       
+                       [[0.221176, -0.670103, -0.708550],
+                        [0.199673, 0.742268, -0.639663],
+                        [0.954574, -0.000000, 0.297973]],
+                       
+                       [[0.221758, -0.670416, -0.708072],
+                        [0.200368, 0.741986, -0.639773],
+                        [0.954294, -0.000000, 0.298871]]])
+    
+    assert np.allclose(attitude.as_matrix(), matrix)
+
+
+def test_get_target_in_sc_frame():
+
+    ori_path = test_data.path / "20280301_first_10sec.ori"
+    ori = SpacecraftFile.parse_from_file(ori_path)
+
+    target_name = "Crab"
+    target_coord = SkyCoord(l=184.5551, b = -05.7877, unit = (u.deg, u.deg), frame = "galactic")
+
+    path_in_sc = ori.get_target_in_sc_frame(target_name, target_coord)
+
+    assert np.allclose(path_in_sc.lon.deg, 
+                       np.array([118.393522, 118.425255, 118.456868, 118.488362, 118.519735, 
+                                 118.550989, 118.582124, 118.613139, 118.644035, 118.674813, 118.705471]))
+
+    assert np.allclose(path_in_sc.lat.deg, 
+                       np.array([46.733430, 46.687559, 46.641664, 46.595745, 46.549801, 46.503833,
+                                 46.457841, 46.411825, 46.365785, 46.319722, 46.273634]))
+
+
+def test_get_dwell_map():
+
+    response_path =test_data.path / "test_full_detector_response.h5"
+    ori_path = test_data.path / "20280301_first_10sec.ori"
+    ori = SpacecraftFile.parse_from_file(ori_path)
+    
+    target_name = "Crab"
+    target_coord = SkyCoord(l=184.5551, b = -05.7877, unit = (u.deg, u.deg), frame = "galactic")
+    path_in_sc = ori.get_target_in_sc_frame(target_name, target_coord)
+    
+    dwell_map = ori.get_dwell_map(response = response_path)
+    
+    assert np.allclose(dwell_map[:].value, 
+                       np.array([1.895057, 7.615584, 0.244679, 0.244679, 0.000000, 0.000000, 
+                                0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000]))
+
+
+def test_get_psr_rsp():
+
+    response_path = test_data.path / "test_full_detector_response.h5"
+    ori_path = test_data.path / "20280301_first_10sec.ori"
+    ori = SpacecraftFile.parse_from_file(ori_path)
+    
+    target_name = "Crab"
+    target_coord = SkyCoord(l=184.5551, b = -05.7877, unit = (u.deg, u.deg), frame = "galactic")
+    path_in_sc = ori.get_target_in_sc_frame(target_name, target_coord)
+    
+    dwell_map = ori.get_dwell_map(response = response_path)
+    
+    Ei_edges, Ei_lo, Ei_hi, Em_edges, Em_lo, Em_hi, areas, matrix = ori.get_psr_rsp()
+    
+    assert np.allclose(Ei_edges, 
+                       np.array([150., 220., 325., 480., 520., 765., 1120., 1650., 2350., 3450., 5000.]))
+    
+    assert np.allclose(Ei_lo, 
+                       np.array([150., 220., 325., 480., 520., 765., 1120., 1650., 2350., 3450.]))
+    
+    assert np.allclose(Ei_hi, 
+                       np.array([220.,  325.,  480.,  520.,  765., 1120., 1650., 2350., 3450., 5000.]))
+    
+    assert np.allclose(Em_edges, 
+                       np.array([150.,  220.,  325.,  480.,  520.,  765., 1120., 1650., 2350., 3450., 5000.]))
+    
+    assert np.allclose(Em_lo, 
+                       np.array([150.,  220.,  325.,  480.,  520.,  765., 1120., 1650., 2350., 3450.]))
+    
+    assert np.allclose(Em_hi, 
+                       np.array([220.,  325.,  480.,  520.,  765., 1120., 1650., 2350., 3450., 5000.]))
+    
+    assert np.allclose(areas, 
+                       np.array([0.06089862, 0.4563752 , 1.1601573 , 1.6237522 , 2.0216975 , 
+                                 2.2039971 , 2.0773466 , 1.7005537 , 1.1626455 , 0.80194914]))
+    
+    assert np.allclose(matrix, 
+                       np.array([[9.80996430e-01, 4.68325317e-02, 1.82471890e-02, 9.86817386e-03,
+                                  5.82037494e-03, 3.47572053e-03, 2.80415593e-03, 3.13903880e-03,
+                                  4.89909900e-03, 6.68705115e-03],
+                                  [1.90035217e-02, 9.44634676e-01, 1.28470331e-01, 9.38407257e-02,
+                                  4.32382338e-02, 2.23877952e-02, 1.63043533e-02, 1.73287615e-02,
+                                  2.80312393e-02, 3.78256924e-02],
+                                   [0.00000000e+00, 8.53277557e-03, 8.48568857e-01, 2.18858123e-01,
+                                  1.85861006e-01, 7.39495233e-02, 4.45922092e-02, 4.06639054e-02,
+                                  6.96888119e-02, 9.27841067e-02],
+                                   [0.00000000e+00, 0.00000000e+00, 4.71363496e-03, 6.62667990e-01,
+                                  6.19757064e-02, 2.71992888e-02, 1.51670892e-02, 1.46367634e-02,
+                                  3.69769707e-02, 7.03022778e-02],
+                                   [0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 1.47649962e-02,
+                                  7.00923026e-01, 2.60504693e-01, 9.65307504e-02, 7.03864172e-02,
+                                  1.15635686e-01, 1.53913230e-01],
+                                   [0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
+                                  2.18164618e-03, 6.11085474e-01, 2.28024259e-01, 9.29291621e-02,
+                                  1.14003479e-01, 1.54005408e-01],
+                                   [0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
+                                  0.00000000e+00, 1.39757351e-03, 5.95472097e-01, 2.54652113e-01,
+                                  1.32362068e-01, 1.71157718e-01],
+                                   [0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
+                                  0.00000000e+00, 0.00000000e+00, 1.10507896e-03, 5.05610526e-01,
+                                  2.00507417e-01, 1.41500503e-01],
+                                   [0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
+                                  0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 6.53312833e-04,
+                                  2.97714621e-01, 1.26633704e-01],
+                                   [0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
+                                  0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
+                                  1.80651987e-04, 4.51902114e-02]]))
+    
+def test_get_arf():
+    
+    response_path = test_data.path / "test_full_detector_response.h5"
+    ori_path = test_data.path / "20280301_first_10sec.ori"
+    ori = SpacecraftFile.parse_from_file(ori_path)
+    
+    target_name = "Crab"
+    target_coord = SkyCoord(l=184.5551, b = -05.7877, unit = (u.deg, u.deg), frame = "galactic")
+    
+    path_in_sc = ori.get_target_in_sc_frame(target_name, target_coord)
+    
+    dwell_map = ori.get_dwell_map(response = response_path)
+    
+    _ = ori.get_psr_rsp()
+    
+    ori.get_arf(out_name = "test")
+    
+    fits_file = fits.open("test.arf")
+    
+    assert np.allclose(fits_file[1].data.field("ENERG_LO"), 
+                       np.array([150.,  220.,  325.,  480.,  520.,  765., 1120., 1650., 2350., 3450.]))
+    
+    assert np.allclose(fits_file[1].data.field("ENERG_HI"), 
+                       np.array([220.,  325.,  480.,  520.,  765., 1120., 1650., 2350., 3450., 5000.]))   
+    
+    assert np.allclose(fits_file[1].data.field("SPECRESP"), 
+                       np.array([0.06089862, 0.4563752 , 1.1601573 , 1.6237522 , 2.0216975 , 
+                                 2.2039971 , 2.0773466 , 1.7005537 , 1.1626455 , 0.80194914]))   
+    
+    os.remove("test.arf")
+
+def test_get_rmf():
+    
+    response_path = test_data.path / "test_full_detector_response.h5"
+    ori_path = test_data.path / "20280301_first_10sec.ori"
+    ori = SpacecraftFile.parse_from_file(ori_path)
+    
+    target_name = "Crab"
+    target_coord = SkyCoord(l=184.5551, b = -05.7877, unit = (u.deg, u.deg), frame = "galactic")
+    
+    path_in_sc = ori.get_target_in_sc_frame(target_name, target_coord)
+    
+    dwell_map = ori.get_dwell_map(response = response_path)
+    
+    _ = ori.get_psr_rsp()
+    
+    ori.get_rmf(out_name = "test")
+    
+    fits_file = fits.open("test.rmf")
+    
+    assert np.allclose(fits_file[1].data.field("ENERG_LO"), 
+                       np.array([150.,  220.,  325.,  480.,  520.,  765., 1120., 1650., 2350., 3450.]))
+    
+    assert np.allclose(fits_file[1].data.field("ENERG_HI"), 
+                       np.array([220.,  325.,  480.,  520.,  765., 1120., 1650., 2350., 3450., 5000.]))   
+    
+    assert np.allclose(fits_file[1].data.field("N_GRP"), 
+                       np.array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1]))  
+    
+    matrix_flattened = []
+    for i in fits_file[1].data.field("MATRIX"):
+        matrix_flattened += i.tolist()
+    
+    assert np.allclose(matrix_flattened,
+                        [0.9809964299201965,
+                         0.019003521651029587,
+                         0.046832531690597534,
+                         0.9446346759796143,
+                         0.008532775565981865,
+                         0.01824718900024891,
+                         0.12847033143043518,
+                         0.848568856716156,
+                         0.0047136349603533745,
+                         0.009868173860013485,
+                         0.09384072571992874,
+                         0.21885812282562256,
+                         0.662667989730835,
+                         0.014764996245503426,
+                         0.005820374935865402,
+                         0.043238233774900436,
+                         0.1858610063791275,
+                         0.06197570636868477,
+                         0.7009230256080627,
+                         0.00218164618127048,
+                         0.003475720528513193,
+                         0.02238779515028,
+                         0.07394952327013016,
+                         0.027199288830161095,
+                         0.26050469279289246,
+                         0.6110854744911194,
+                         0.0013975735055282712,
+                         0.0028041559271514416,
+                         0.01630435325205326,
+                         0.04459220916032791,
+                         0.01516708917915821,
+                         0.09653075039386749,
+                         0.22802425920963287,
+                         0.5954720973968506,
+                         0.001105078961700201,
+                         0.0031390388030558825,
+                         0.017328761518001556,
+                         0.04066390544176102,
+                         0.014636763371527195,
+                         0.07038641721010208,
+                         0.0929291620850563,
+                         0.25465211272239685,
+                         0.5056105256080627,
+                         0.000653312832582742,
+                         0.004899099003523588,
+                         0.0280312392860651,
+                         0.0696888118982315,
+                         0.03697697073221207,
+                         0.11563568562269211,
+                         0.11400347948074341,
+                         0.13236206769943237,
+                         0.20050741732120514,
+                         0.29771462082862854,
+                         0.0001806519867386669,
+                         0.006687051150947809,
+                         0.03782569244503975,
+                         0.0927841067314148,
+                         0.07030227780342102,
+                         0.1539132297039032,
+                         0.15400540828704834,
+                         0.17115771770477295,
+                         0.14150050282478333,
+                         0.12663370370864868,
+                         0.04519021138548851])
+    
+    os.remove("test.rmf")
+
+
+def test_get_pha():
+
+    response_path = test_data.path / "test_full_detector_response.h5"
+    ori_path = test_data.path / "20280301_first_10sec.ori"
+    ori = SpacecraftFile.parse_from_file(ori_path)
+    
+    target_name = "Crab"
+    target_coord = SkyCoord(l=184.5551, b = -05.7877, unit = (u.deg, u.deg), frame = "galactic")
+    
+    path_in_sc = ori.get_target_in_sc_frame(target_name, target_coord)
+    dwell_map = ori.get_dwell_map(response = response_path)
+    _ = ori.get_psr_rsp()
+    ori.get_arf(out_name = "test")
+    ori.get_rmf(out_name = "test")
+    
+    counts = np.array([0.01094232, 0.04728866, 0.06744612, 0.01393708, 0.05420688,
+                       0.03141498, 0.01818584, 0.00717219, 0.00189568, 0.00010503])*1000
+    
+    errors = np.sqrt(counts)
+    
+    ori.get_pha(src_counts=counts, errors=errors, exposure_time=10)
+    
+    os.remove("test.arf")
+    os.remove("test.rmf")
+    
+    fits_file = fits.open("test.pha")
+    os.remove("test.pha")
+    
+    assert np.allclose(fits_file[1].data.field("CHANNEL"), 
+                           np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9]))
+    
+    assert np.allclose(fits_file[1].data.field("COUNTS"), 
+                           np.array([10, 47, 67, 13, 54, 31, 18,  7,  1,  0]))
+    
+    assert np.allclose(fits_file[1].data.field("STAT_ERR"), 
+                           np.array([3, 6, 8, 3, 7, 5, 4, 2, 1, 0]))
+
+def test_plot_arf():
+
+    response_path = test_data.path / "test_full_detector_response.h5"
+    ori_path = test_data.path / "20280301_first_10sec.ori"
+    ori = SpacecraftFile.parse_from_file(ori_path)
+    
+    target_name = "Crab"
+    target_coord = SkyCoord(l=184.5551, b = -05.7877, unit = (u.deg, u.deg), frame = "galactic")
+    
+    path_in_sc = ori.get_target_in_sc_frame(target_name, target_coord)
+    dwell_map = ori.get_dwell_map(response = response_path)
+    _ = ori.get_psr_rsp()
+    ori.get_arf(out_name = "test")
+    
+    ori.plot_arf()
+    
+    assert Path("Effective_area_for_test.png").exists()
+    
+    os.remove("test.arf")
+    os.remove("Effective_area_for_test.png")
+
+def test_plot_rmf():
+
+    response_path = test_data.path / "test_full_detector_response.h5"
+    ori_path = test_data.path / "20280301_first_10sec.ori"
+    ori = SpacecraftFile.parse_from_file(ori_path)
+    
+    target_name = "Crab"
+    target_coord = SkyCoord(l=184.5551, b = -05.7877, unit = (u.deg, u.deg), frame = "galactic")
+    
+    path_in_sc = ori.get_target_in_sc_frame(target_name, target_coord)
+    dwell_map = ori.get_dwell_map(response = response_path)
+    _ = ori.get_psr_rsp()
+    ori.get_rmf(out_name = "test")
+    
+    ori.plot_rmf()
+    
+    assert Path("Redistribution_matrix_for_test.png").exists()
+    
+    os.remove("test.rmf")
+    os.remove("Redistribution_matrix_for_test.png")
+
+def test_source_interval():
+
+    response_path = test_data.path / "test_full_detector_response.h5"
+    ori_path = test_data.path / "20280301_first_10sec.ori"
+    ori = SpacecraftFile.parse_from_file(ori_path)
+
+    new_ori = ori.source_interval(Time(ori._load_time[0]+0.1, format = "unix"), 
+                                  Time(ori._load_time[0]+2.1, format = "unix"))
+
+    assert np.allclose(new_ori._load_time, 
+                       np.array([1.835478e+09, 1.835478e+09, 1.835478e+09, 1.835478e+09]))
+
+    assert np.allclose(new_ori._x_direction.flatten(), 
+                       np.array([41.86062093, 73.14368765, 41.88225011, 73.09517927, 
+                                 41.90629597, 73.0412838 , 41.9087019 , 73.03589454]))
+
+    assert np.allclose(new_ori._z_direction.flatten(), 
+                       np.array([221.86062093,  16.85631235, 221.88225011,  16.90482073,
+                                221.90629597,  16.9587162 , 221.9087019 ,  16.96410546]))
+    
diff --git a/tests/ts_map/__init__.py b/tests/ts_map/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/tests/ts_map/test_fast_norm_fit.py b/tests/ts_map/test_fast_norm_fit.py
new file mode 100644
index 00000000..2c2c9746
--- /dev/null
+++ b/tests/ts_map/test_fast_norm_fit.py
@@ -0,0 +1,35 @@
+from cosipy import test_data
+from cosipy import FastTSMap
+from histpy import Histogram
+import numpy as np
+from cosipy.ts_map import FastNormFit as fnf
+
+
+def test_solve():
+
+    # read the signal+background
+    src_bkg_path = test_data.path / "ts_map_src_bkg.h5"
+    src_bkg = Histogram.open(src_bkg_path)
+
+    # read the background
+    bkg_path = test_data.path / "ts_map_bkg.h5"
+    bkg = Histogram.open(bkg_path)
+
+    # get the cds arrays of src_bkg and bkg
+    src_bkg_cds_array = FastTSMap.get_cds_array(src_bkg, [0,10])
+    bkg_cds_array = FastTSMap.get_cds_array(bkg, [0,10])
+
+    # read the cds array of expectation
+    ei_path = test_data.path / "ei_cds_array.npy"
+    ei_cds_array = np.load(ei_path)
+
+
+    # calculate the ts value
+    fit = fnf(max_iter=1000)
+    result = fit.solve(src_bkg_cds_array, bkg_cds_array, ei_cds_array)
+
+    assert result[0] == 187.3360310655543
+
+    assert result[1] == 0.02119470713546078
+
+    assert result[2] == 0.0055665881497504646
\ No newline at end of file
diff --git a/tests/ts_map/test_fast_ts_map.py b/tests/ts_map/test_fast_ts_map.py
new file mode 100644
index 00000000..1b721cc6
--- /dev/null
+++ b/tests/ts_map/test_fast_ts_map.py
@@ -0,0 +1,44 @@
+from cosipy import test_data
+from pytest import approx
+from threeML import Powerlaw
+from cosipy import FastTSMap, SpacecraftFile
+from histpy import Histogram
+import numpy as np
+from astropy.coordinates import SkyCoord
+import astropy.units as u
+
+
+def test_parallel_ts_fit():
+
+    src_bkg_path = test_data.path / "ts_map_src_bkg.h5"
+    bkg_path = test_data.path / "ts_map_bkg.h5"
+    response_path = test_data.path / "test_full_detector_response.h5"
+
+    orientation_path = test_data.path / "20280301_2s.ori"
+    ori = SpacecraftFile.parse_from_file(orientation_path)
+
+    src_bkg = Histogram.open(src_bkg_path).project(['Em', 'PsiChi', 'Phi'])
+    bkg = Histogram.open(bkg_path).project(['Em', 'PsiChi', 'Phi'])
+
+    ts = FastTSMap(data = src_bkg, bkg_model = bkg, orientation = ori, 
+                   response_path = response_path, cds_frame = "local", scheme = "RING")
+
+    hypothesis_coords = FastTSMap.get_hypothesis_coords(nside = 1)
+    
+    index = -2.2
+    K = 10 / u.cm / u.cm / u.s / u.keV
+    piv = 100 * u.keV
+    spectrum = Powerlaw()
+    spectrum.index.value = index
+    spectrum.K.value = K.value
+    spectrum.piv.value = piv.value 
+    spectrum.K.unit = K.unit
+    spectrum.piv.unit = piv.unit
+
+    ts_results = ts.parallel_ts_fit(hypothesis_coords = hypothesis_coords, energy_channel = [2,3], spectrum = spectrum, ts_scheme = "RING", cpu_cores = 2)
+
+    assert np.allclose(ts_results[:,1],
+                       [51.30709447, 51.16302889, 51.11429069, 
+                        51.19306142, 51.27823575, 51.30579709, 
+                        51.09094512, 51.10914182, 51.30271261, 
+                        51.27412572, 51.15872202, 51.29249638])
\ No newline at end of file