added resamplig capabilities to convert_voltage2adc

grand/sim/detector/adc.py

@@ -3,6 +3,7 @@
 """
 import numpy as np
 import logging
+import scipy.fft as sf
 
 logger = logging.getLogger(__name__)
 
@@ -21,6 +22,39 @@ def __init__(self):
         self.max_bit_value = 8192 # 14 bit ADC;  2 x 2^13 bits for negative and positive ADC values
         self.max_voltage   = 9e5  # [µV]; saturation voltage of ADC (absolute value)
 
+
+    def downsample(self,
+                  voltage_trace,input_sampling_rate_mhz):
+        '''
+        downsamples the voltage trace to the target sampling rate
+        
+        Arguments
+        ---------
+        `voltage_trace`
+        type        : np.ndarray[double]
+        units       : uV 
+        description : Array of voltage traces, with shape (N_du,3,N_samples)
+                                    
+        Returns
+        -------
+        `downsampled_voltage_trace`
+        type        : np.ndarray[double]
+        units       : uV
+        description : Array of downsamplef voltage traces, with shape (N_du,3,N_samples)
+        '''
+        if self.sampling_rate != input_sampling_rate_mhz : 
+          #compute the fft
+          voltage_trace_f=sf.rfft(voltage_trace)
+          #compute new number of points
+          ratio=(self.sampling_rate/input_sampling_rate_mhz)        
+          m=int(np.shape(voltage_trace)[2]*ratio)
+          logger.info(f"resampling the voltage from {input_sampling_rate_mhz} to an ADC of {self.sampling_rate} MHz")        
+          downsampled_voltage_trace=sf.irfft(voltage_trace,m)*ratio
+        else:
+          downsampled_voltage_trace=voltage_trace
+
+        return downsampled_voltage_trace
+
 
     def _digitize(self,
                   voltage_trace):
@@ -51,7 +85,6 @@ def _digitize(self,
         adc_trace = np.trunc(adc_trace).astype(int)
 
         return adc_trace
-
 
     def _saturate(self,
                   adc_trace):
@@ -79,7 +112,6 @@ def _saturate(self,
 
         return saturated_adc_trace
 
-
     def process(self,
                 voltage_trace,
                 noise_trace=None):
@@ -107,7 +139,8 @@ def process(self,
         description : Array of ADC traces with shape (N_du,3,N_samples)
         '''
 
-        assert isinstance(voltage_trace,np.ndarray)
+        assert isinstance(voltage_trace,np.ndarray)       
+
 
         adc_trace = self._digitize(voltage_trace)
 
@@ -122,4 +155,4 @@ def process(self,
         # Make sure the saturation occurs AFTER adding noise
         adc_trace = self._saturate(adc_trace)
 
-        return adc_trace
+        return adc_trace

scripts/convert_efield2efield.py

@@ -511,7 +511,7 @@ def impz(b, a):
        #to resample we use fourier interpolation, becouse it seems to be better than scipy.decimate (points are closer to the original trace)
        vout = sf.irfft(vout_f, m)*ratio #renormalize the amplitudes
 
-       if(event_idx==0):
+       if(event_idx==0 and PLOT):
          plt.scatter(np.arange(0,len(vout[6][0]))/ratio,vout[6][0],label="sampled",c="red")
 
 

scripts/convert_voltage2adc.py

@@ -147,7 +147,12 @@ def manage_args():
                         type=str,
                         default=None,
                         help='Path to directory containing files with measured noise in GrandRoot format (TADC). Default adds no noise.')
-
+    parser.add_argument(
+                        "--target_sampling_rate_mhz",
+                        type=float,
+                        default=0,
+                        help="Target sampling rate of the data in Mhz",
+    )       
     parser.add_argument('-s',
                         '--seed',
                         type=int,
@@ -170,25 +175,28 @@ def manage_args():
 
     #-#-#- Get parser arguments -#-#-#
     args      = manage_args()
-    f_input   = args.in_file
+    f_input_dir   = args.in_file
     f_output  = args.out_file
     noise_dir = args.noise_dir
 
+    f_input_file=glob.glob(f_input_dir+"/voltage_*_L0_*.root")[0]
+
     if f_output == None:
-        f_output = f_input.replace('voltage','adc')
+        f_output = f_input_file.replace('voltage','adc')
     if noise_dir == None:
         noise_trace = None
 
     manage_log.create_output_for_logger(args.verbose,log_stdout=True)
     logger.info( manage_log.string_begin_script() )
     #-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#
 
-    logger.info(f'Converting voltage traces from {f_input} to ADC traces')
+    logger.info(f'Converting voltage traces from {f_input_file} to ADC traces')
 
     #-#-#- Load TVoltage -#-#-#
-    df       = rt.DataFile(f_input)
+    df       = rt.DataDirectory(f_input_dir)
     tvoltage = df.tvoltage
     entries  = tvoltage.get_number_of_entries()
+    trun = df.trun
 
     #-#-#- Prepare TADC -#-#-#
     if os.path.exists(f_output):
@@ -199,6 +207,7 @@ def manage_args():
 
     #-#-#- Initiate ADC object and RNG -#-#-#
     adc = ADC()
+
     rng = np.random.default_rng(args.seed)
     if noise_dir is not None:
         logger.info(f'Set RNG seed to {args.seed}')
@@ -211,6 +220,19 @@ def manage_args():
         tvoltage.get_entry(entry)
         voltage_trace = np.array(tvoltage.trace)
 
+        event_number = tvoltage.event_number
+        run_number = tvoltage.run_number
+        trun.get_run(run_number)
+        event_dus_indices = tvoltage.get_dus_indices_in_run(trun)
+        dt_ns = np.asarray(trun.t_bin_size)[event_dus_indices] # sampling time in ns, sampling freq = 1e9/dt_ns. 
+        f_samp_mhz = 1e3/dt_ns                                 # MHz  
+        input_sampling_rate_mhz = f_samp_mhz[0]                # and here we asume all sampling rates are the same!. In any case, we are asuming all the ADCs are the same...        
+
+        #-#-#- Downsample if needed -#-#-# (this could be added to the "process" method to hide it from the public, and add input_sampling_rate as input to process.
+        if( input_sampling_rate_mhz != adc.sampling_rate):         
+           voltage_trace=adc.downsample(voltage_trace,input_sampling_rate_mhz)
+
+
         #-#-#- Get noise trace if requested -#-#-#
         if noise_dir is not None:
             noise_trace = get_noise_trace(noise_dir,