Add pedestal variations for FSD

[diaza]

This PR adds pedestal variations for FSD.

To test the pedestal implementation, I take the difference of the ADC values for each hit with and without pedestal variation. Given this equation, that value should be consistent with (pedestal_from_file - 580.)*256/(v_ref-v_cm). Where pedestal_from_file is the pedestal extracted from the inputted pedestal file for that given pixel, and 580 is the default pedestal voltage. Below is the code used, and the outputs

with open('/global/common/software/dune/inputs/FSD/pedestals/FSD_pedestals_20241112.json', 'r') as data_file:
    data_swapped = json.load(data_file)

default = f2['packets']['io_group','io_channel', 'chip_id', 'channel_id', 'dataword'][mask2].astype(np.dtype([('io_group', 'int64'), ('io_channel', 'int64'), ('chip_id', 'int64'), ('channel_id', 'int64'), ('dataword', 'int64')]))
default_tile_id = ((default['io_channel']) - 1) //4 +1 + (default['io_group']-1)*10
default_unique_id = ((default['io_group']*10000+default_tile_id)*1000 + default['chip_id'])*100 + default['channel_id']

varied = f3['packets']['io_group','io_channel', 'chip_id', 'channel_id', 'dataword'][mask3].astype(np.dtype([('io_group', 'int64'), ('io_channel', 'int64'), ('chip_id', 'int64'), ('channel_id', 'int64'), ('dataword', 'int64')]))
varied_tile_id = ((varied['io_channel']) - 1) //4 +1 + (varied['io_group']-1)*10
varied_unique_id = ((varied['io_group']*10000+varied_tile_id)*1000 + varied['chip_id'])*100 + varied['channel_id']

n=len(varied)
mv_diff = np.zeros(n)
adc_diff = np.zeros(n)
adc_counts = 256
v_ref = 1568.
v_cm=478.
n_notfound=0
n_maxed = 0
for i in range(n):

    if not varied['dataword'][i]<255:
        n_maxed+=1
        continue
    try:
        ped = data_swapped[str(varied_unique_id[i])]
        mv_diff[i] = (ped['pedestal_mv'] - 580.)*256/(v_ref-v_cm)
        
    except KeyError:
        n_notfound+=1
        mv_diff[i] = 0
        
    adc_diff[i] = varied['dataword'][i] - default['dataword'][i]

print("Number of total hits: ", n)
print("Number of hits without pedestal in file: ", n_notfound)
print("Number of satured hits: ", n_maxed)

print("Are all extracted adc differences within expectation?") 
print(np.all( np.abs(mv_diff - adc_diff) < 1))

Number of total hits:  318948
Number of hits without pedestal in file:  2156
Number of satured hits:  103
Are all extracted adc differences within expectation?
True

[diaza]

I also did a check with the flow outputs. I flowed a file with the new pedestals. With the output, I looped over each hit and manually calculate what Q should be given the adc value and the pedestal from the pedestal file. I check that this calculated Q matches with the output from flow.

with open('/global/common/software/dune/inputs/FSD/pedestals/FSD_pedestals_20241112.json', 'r') as f:
    data = json.load(f)

def charge_from_dataword(dw, vref, vcm, ped, adc_counts, gain):
        return (dw / adc_counts * (vref - vcm) + vcm - ped) / gain

mask = ~dm['charge/calib_prompt_hits/data']['is_disabled']
packets = dm['charge/calib_prompt_hits','/charge/packets'][['dataword', 'io_group', 'io_channel', 'chip_id', 'channel_id']][:,0][mask]
hits = dm['charge/calib_prompt_hits/data'][mask]

diff = list()
notfound = list()
n_notfound = 0
unique_ids = list()
for (hit, packet) in zip(hits, packets):

    (adc, io_group, io_channel, chip_id, channel_id) = packet[['dataword', 'io_group', 'io_channel', 'chip_id', 'channel_id']]
    tile_id = 10*(io_group-1) + (io_channel-1)//4 +1

    unique_id = (int(io_group)*1000_000_000
                      + int(tile_id)*100_000
                      + int(chip_id)*100
                      + int(channel_id))
    unique_ids.append(unique_id)

    try:
        ped = data[str(unique_id)]['pedestal_mv']
    except KeyError:
        ped= 580
        n_notfound+=1
        notfound.append(unique_id)
        
    charge = charge_from_dataword(dw = adc, vref = 1568., vcm=478.1, ped = ped, adc_counts=256, gain=4.522)
    
    if charge - hit['Q'] != 0:
        print((io_group, io_channel, chip_id, channel_id))
        print(tile_id)
        print(i, unique_id, charge, hit['Q'])
        raise ValueError("fuck")
    
    diff.append(charge -  hit['Q'])

print("Number of hits: {}".format(len(packets)))
print("Number of unique ids: {}".format(np.unique(np.array(unique_ids)).shape))
print("Number of hits without corresponding id in pedestal file: {}".format(n_notfound))
print("Number of unique ids without corresponding id in pedestal file: {}".format(np.unique(np.array(notfound)).shape))
      
print("Are valid hits all consistent?", np.all(np.array(diff) == 0.))

Number of hits: 315570
Number of unique ids: (211301,)
Number of hits without corresponding id in pedestal file: 217
Number of unique ids without corresponding id in pedestal file: (140,)
Are valid hits all consistent? True

[diaza]

(The checks above are copied from the checks I did with the 2x2 pedestal implementations here and here.)

[diaza]

Note @cuddandr, I think the gains used in larndsim are currently different than the gains in flow.

[cuddandr]

Note @cuddandr, I think the gains used in larndsim are currently different than the gains in flow.

Where are those even defined in flow? @jaafar-chakrani updated the FSD gain in larnd-sim so I can see them being different from flow right now.

[diaza]

@cuddandr
In flow, they should be defined in CalibHitBuilderMC.yaml. You can see them already defined in CalibHitBuilderData.yaml. We might also need to update the voltages in larnd-sim, I notice that they're different in fsd.yaml compared to the ones in CalibHitBuilderData.yaml. You know which ones are the correct ones?

[cuddandr]

Of course not. But the voltage values in flow are more recent than the ones in larnd-sim, which are the original values when that file was first made.

So I'd be inclined to copy the values from CalibHitBuilderData to CalibHitBuilderMC in flow and to fsd.yaml in larnd-sim.

[cuddandr]

Everything looks good to me. I also updated the voltages in fsd.yaml to match what is used in flow after conferring with Jaafar. With these changes larnd-sim still runs fine for me.