@@ -103,13 +103,15 @@ def __init__(self, **params):
103103 self .packets_index_name = params .get ('packets_index_name' , self .packets_dset_name + '_index' )
104104 self .t0_dset_name = params .get ('t0_dset_name' )
105105 self .pedestal_file = params .get ('pedestal_file' , '' )
106+ self .gain_file = params .get ('gain_file' , '' )
106107 self .configuration_file = params .get ('configuration_file' , '' )
107108 self .pedestal_mv = params .get ('pedestal_mv' , self .default_pedestal_mv )
108109 self .vref_mv = params .get ('vref_mv' , self .default_vref_mv )
109110 self .vcm_mv = params .get ('vcm_mv' , self .default_vcm_mv )
110111 self .adc_counts = params .get ('adc_counts' , self .default_adc_counts )
111112 self .gain = params .get ('gain' , self .default_gain )
112113 self .adc_droop_calibration = params .get ('adc_droop_calibration' , False )
114+ self .elifetime_calibration = params .get ('elifetime_calibration' ,False )
113115 self .hit_ref = params .get ('hit_ref' , True )
114116
115117 #: ASIC ADC configuration lookup table
@@ -123,9 +125,14 @@ def __init__(self, **params):
123125 pedestal_mv = self .pedestal_mv
124126 ))
125127
128+ self .gains = defaultdict (lambda : dict (
129+ gain = self .gain
130+ ))
131+
126132 def init (self , source_name ):
127133 super (CalibHitBuilder , self ).init (source_name )
128134 self .load_pedestals ()
135+ self .load_gains ()
129136 self .load_configurations ()
130137
131138 def run (self , source_name , source_slice , cache ):
@@ -172,6 +179,7 @@ def run(self, source_name, source_slice, cache):
172179 packets_dset = self .packets_dset_name ,
173180 t0_dset = self .t0_dset_name ,
174181 pedestal_file = self .pedestal_file ,
182+ gain_file = self .gain_file ,
175183 configuration_file = self .configuration_file ,
176184 adc_droop_calibration = self .adc_droop_calibration
177185 )
@@ -236,9 +244,13 @@ def run(self, source_name, source_slice, cache):
236244 # This time we add the rollover period instead of subtracting.
237245 drift_t [after_sync_mask ] += resources ['RunData' ].rollover_ticks
238246
239- drift_d = drift_t * (resources ['LArData' ].v_drift * resources ['RunData' ].crs_ticks ) / units .cm # convert mm -> cm
247+ v_drift_arr = resources ['LArData' ].v_drift
248+ if len (v_drift_arr ) == 1 :
249+ v_drift = v_drift_arr [0 ] #Default vdrift
250+ else :
251+ v_drift = v_drift_arr [(packets_arr ['io_group' ]- 1 )// 2 ]
252+ drift_d = drift_t * (v_drift * resources ['RunData' ].crs_ticks ) / units .cm # convert mm -> cm
240253 x = resources ['Geometry' ].get_drift_coordinate (packets_arr ['io_group' ],packets_arr ['io_channel' ],drift_d )
241-
242254 ## true drift position pair
243255 #if has_mc_truth:
244256 # drift_t_true = packet_seg_bt_arr['t'] #us
@@ -265,6 +277,11 @@ def run(self, source_name, source_slice, cache):
265277 for unique_id in hit_uniqueid_str ])
266278 else :
267279 ped = np .full (len (hit_uniqueid_str ), self .pedestal_mv )
280+ if self .gain_file != '' :
281+ gain = np .array ([self .gains [unique_id ]['gain' ] for unique_id in hit_uniqueid_str ])
282+ else :
283+ gain = np .full (len (hit_uniqueid_str ), self .gain )
284+
268285 calib_hits_arr ['id' ] = calib_hits_slice .start + np .arange (n , dtype = int )
269286 calib_hits_arr ['x' ] = x
270287 #if has_mc_truth:
@@ -277,16 +294,20 @@ def run(self, source_name, source_slice, cache):
277294 calib_hits_arr ['io_channel' ] = packets_arr ['io_channel' ]
278295 calib_hits_arr ['chip_id' ] = packets_arr ['chip_id' ]
279296 calib_hits_arr ['channel_id' ] = packets_arr ['channel_id' ]
280- hits_charge = self .charge_from_dataword (packets_arr ['dataword' ], vref , vcm , ped , self .adc_counts , self . gain ) # ke-
297+ hits_charge = self .charge_from_dataword (packets_arr ['dataword' ], vref , vcm , ped , self .adc_counts , gain ) # ke-
281298 calib_hits_arr ['Q_raw' ] = hits_charge # ke-
282299 if self .adc_droop_calibration :
283- hits_charge_calibrated = self .charge_from_dataword_corrected (packets_arr ['dataword' ], packets_arr ['timestamp' ], hit_uniqueid , vref , vcm , ped , self .adc_counts , self . gain ) # ke-
284- calib_hits_arr ['Q' ] = hits_charge_calibrated # ke-
300+ hits_charge_calibrated = self .charge_from_dataword_corrected (packets_arr ['dataword' ], packets_arr ['timestamp' ], hit_uniqueid , vref , vcm , ped , self .adc_counts , gain ) # ke-
301+ calib_hits_arr ['Q' ] = hits_charge_calibrated # ke-
285302 else :
286- calib_hits_arr ['Q' ] = hits_charge
287-
303+ calib_hits_arr ['Q' ] = hits_charge # ke-
304+
305+
306+
288307 #FIXME supply more realistic dEdx in the recombination; also apply measured electron lifetime
289- calib_hits_arr ['E' ] = calib_hits_arr ['Q' ] * (1000 * units .e ) / resources ['LArData' ].ionization_recombination (mode = 2 ,dEdx = 2 ) * (resources ['LArData' ].ionization_w / units .MeV ) # MeV
308+ calib_hits_arr ['E' ] = calib_hits_arr ['Q' ] * (1000 * units .e ) / resources ['LArData' ].ionization_recombination (mode = 2 ,dEdx = 2 ) * (resources ['LArData' ].ionization_w / units .MeV ) # MeV
309+ if self .elifetime_calibration :
310+ calib_hits_arr ['E' ] /= resources ['LArData' ].charge_reduction_lifetime (t_drift = (drift_t * resources ['RunData' ].crs_ticks )) # ke- we change the drift_t to µs
290311 #if has_mc_truth:
291312 # true_recomb = resources['LArData'].ionization_recombination(mode=2,dEdx=packet_seg_bt_arr['dEdx'])
292313 # calib_hits_arr['E_true_recomb_elife'] = np.divide(hits_charge.reshape((hits_charge.shape[0],1)) * (1000 * units.e), true_recomb, out=np.zeros_like(true_recomb), where=true_recomb!=0) / resources['LArData'].charge_reduction_lifetime(t_drift=drift_t_true) * (resources['LArData'].ionization_w / units.MeV) # MeV
@@ -384,6 +405,12 @@ def load_pedestals(self):
384405 for key , value in json .load (infile ).items ():
385406 self .pedestal [key ] = value
386407
408+ def load_gains (self ):
409+ if self .gain_file != '' :
410+ with open (self .gain_file , 'r' ) as infile :
411+ for key , value in json .load (infile ).items ():
412+ self .gains [key ] = value
413+
387414 def load_configurations (self ):
388415 if self .configuration_file != '' and not resources ['RunData' ].is_mc :
389416 with open (self .configuration_file , 'r' ) as infile :
0 commit comments