Improved support for multi-tile detector

cli/dumpTree.py

@@ -0,0 +1,188 @@
+#! /usr/bin/env python
+#
+# Read almost every field in the event tree.
+#
+
+from math import sqrt
+
+import numpy as np
+import fire
+import h5py
+
+from ROOT import TG4Event, TFile
+
+# Print the fields in a TG4PrimaryParticle object
+def printPrimaryParticle(depth, primaryParticle):
+    print(depth,"Class: ", primaryParticle.ClassName())
+    print(depth,"Track Id:", primaryParticle.GetTrackId())
+    print(depth,"Name:", primaryParticle.GetName())
+    print(depth,"PDG Code:",primaryParticle.GetPDGCode())
+    print(depth,"Momentum:",primaryParticle.GetMomentum().X(),
+          primaryParticle.GetMomentum().Y(),
+          primaryParticle.GetMomentum().Z(),
+          primaryParticle.GetMomentum().E(),
+          primaryParticle.GetMomentum().P(),
+          primaryParticle.GetMomentum().M())
+
+# Print the fields in an TG4PrimaryVertex object
+def printPrimaryVertex(depth, primaryVertex):
+    print(depth,"Class: ", primaryVertex.ClassName())
+    print(depth,"Position:", primaryVertex.GetPosition().X(),
+          primaryVertex.GetPosition().Y(),
+          primaryVertex.GetPosition().Z(),
+          primaryVertex.GetPosition().T())
+    print(depth,"Generator:",primaryVertex.GetGeneratorName())
+    print(depth,"Reaction:",primaryVertex.GetReaction())
+    print(depth,"Filename:",primaryVertex.GetFilename())
+    print(depth,"InteractionNumber:",primaryVertex.GetInteractionNumber())
+    depth = depth + ".."
+    for infoVertex in primaryVertex.Informational:
+        printPrimaryVertex(depth,infoVertex)
+    for primaryParticle in primaryVertex.Particles:
+        printPrimaryParticle(depth,primaryParticle)
+
+# Print the fields in a TG4TrajectoryPoint object
+def printTrajectoryPoint(depth, trajectoryPoint):
+    print(depth,"Class: ", trajectoryPoint.ClassName())
+    print(depth,"Position:", trajectoryPoint.GetPosition().X(),
+          trajectoryPoint.GetPosition().Y(),
+          trajectoryPoint.GetPosition().Z(),
+          trajectoryPoint.GetPosition().T())
+    print(depth,"Momentum:", trajectoryPoint.GetMomentum().X(),
+          trajectoryPoint.GetMomentum().Y(),
+          trajectoryPoint.GetMomentum().Z(),
+          trajectoryPoint.GetMomentum().Mag())
+    print(depth,"Process",trajectoryPoint.GetProcess())
+    print(depth,"Subprocess",trajectoryPoint.GetSubprocess())
+
+# Print the fields in a TG4Trajectory object
+def printTrajectory(depth, trajectory):
+    print(depth,"Class: ", trajectory.ClassName())
+    depth = depth + ".."
+    print(depth,"Track Id/Parent Id:",
+          trajectory.GetTrackId(),
+          trajectory.GetParentId())
+    print(depth,"Name:",trajectory.GetName())
+    print(depth,"PDG Code",trajectory.GetPDGCode())
+    print(depth,"Initial Momentum:",trajectory.GetInitialMomentum().X(),
+          trajectory.GetInitialMomentum().Y(),
+          trajectory.GetInitialMomentum().Z(),
+          trajectory.GetInitialMomentum().E(),
+          trajectory.GetInitialMomentum().P(),
+          trajectory.GetInitialMomentum().M())
+    for trajectoryPoint in trajectory.Points:
+        printTrajectoryPoint(depth,trajectoryPoint)
+
+# Print the fields in a TG4HitSegment object
+def printHitSegment(depth, hitSegment):
+    print(depth,"Class: ", hitSegment.ClassName())
+    print(depth,"Primary Id:", hitSegment.GetPrimaryId())
+    print(depth,"Energy Deposit:",hitSegment.GetEnergyDeposit())
+    print(depth,"Secondary Deposit:", hitSegment.GetSecondaryDeposit())
+    print(depth,"Track Length:",hitSegment.GetTrackLength())
+    print(depth,"Start:", hitSegment.GetStart().X(),
+          hitSegment.GetStart().Y(),
+          hitSegment.GetStart().Z(),
+          hitSegment.GetStart().T())
+    print(depth,"Stop:", hitSegment.GetStop().X(),
+          hitSegment.GetStop().Y(),
+          hitSegment.GetStop().Z(),
+          hitSegment.GetStop().T())
+    print(depth,"Contributor:", [contributor for contributor in hitSegment.Contrib])
+
+# Print the fields in a single element of the SegmentDetectors map.
+# The container name is the key, and the hitSegments is the value (a
+# vector of TG4HitSegment objects).
+def printSegmentContainer(depth, containerName, hitSegments):
+    print(depth,"Detector: ", containerName, hitSegments.size())
+    depth = depth + ".."
+    for hitSegment in hitSegments: printHitSegment(depth, hitSegment)
+
+# Read a file and dump it.
+def dump(input_file, output_file):
+
+    # The input file is generated in a previous test (100TestTree.sh).
+    inputFile = TFile(input_file)
+
+    # Get the input tree out of the file.
+    inputTree = inputFile.Get("EDepSimEvents")
+    print("Class:", inputTree.ClassName())
+
+    # Attach a brach to the events.
+    event = TG4Event()
+    inputTree.SetBranchAddress("Event",event)
+
+    # Read all of the events.
+    entries = inputTree.GetEntriesFast()
+
+    dtype = np.dtype([('eventID','u4'),('z_end','f4'),('trackID','u4'),
+                      ('tran_diff','f4'),('z_start','f4'),('x_end','f4'),
+                      ('y_end','f4'),('n_electrons','u4'),('pdgId','i4'),
+                      ('x_start','f4'),('y_start','f4'),('t_start','f4'),
+                      ('dx','f4'),('long_diff','f4'),('pixel_plane','u4'),
+                      ('t_end','f4'),('dEdx','f4'),('dE','f4'),('t','f4'),
+                      ('y','f4'),('x','f4'),('z','f4')])
+
+    segments = np.array([], dtype=dtype)
+
+    for jentry in range(entries):
+        print(jentry)
+        nb = inputTree.GetEntry(jentry)
+        if nb <= 0:
+            continue
+
+        print("Class: ", event.ClassName())
+        print("Event number:", event.EventId)
+
+        # Dump the primary vertices
+        # for primaryVertex in event.Primaries:
+        #    printPrimaryVertex("PP", primaryVertex)
+
+        # Dump the trajectories
+        trajectories = {'pdgId': [], 'trackId': []}
+        print("Number of trajectories ", len(event.Trajectories))
+
+        for trajectory in event.Trajectories:
+            trajectories['pdgId'].append(trajectory.GetPDGCode())
+            trajectories['trackId'].append(trajectory.GetTrackId())
+
+        # Dump the segment containers
+        print("Number of segment containers:", event.SegmentDetectors.size())
+
+        for containerName, hitSegments in event.SegmentDetectors:
+
+            for hitSegment in hitSegments:
+                segment = np.empty(1, dtype=dtype)
+                segment['eventID'] = jentry
+                segment['trackID'] = trajectories['trackId'][hitSegment.Contrib[0]]
+                segment['x_start'] = hitSegment.GetStart().X()/10
+                segment['y_start'] = hitSegment.GetStart().Y()/10
+                segment['z_start'] = hitSegment.GetStart().Z()/10
+                segment['x_end'] = hitSegment.GetStop().X()/10
+                segment['y_end'] = hitSegment.GetStop().Y()/10
+                segment['z_end'] = hitSegment.GetStop().Z()/10
+                segment['dE'] = hitSegment.GetEnergyDeposit()
+                segment['t'] = 0
+                segment['t_start'] = 0
+                segment['t_end'] = 0
+                xd = segment['x_end'] - segment['x_start']
+                yd = segment['y_end'] - segment['y_start']
+                zd = segment['z_end'] - segment['z_start']
+                dx = sqrt(xd**2 + yd**2 + zd**2)
+                segment['dx'] = dx
+                segment['x'] = (segment['x_start'] + segment['x_end'])/2.
+                segment['y'] = (segment['y_start'] + segment['y_end'])/2.
+                segment['z'] = (segment['z_start'] + segment['z_end'])/2.
+                segment['dEdx'] = hitSegment.GetEnergyDeposit()/dx if dx > 0 else 0
+                segment['pdgId'] = trajectories['pdgId'][hitSegment.Contrib[0]]
+                segment['n_electrons'] = 0
+                segment['long_diff'] = 0
+                segment['tran_diff'] = 0
+                segment['pixel_plane'] = 0
+                segments = np.hstack((segments, segment))
+
+    with h5py.File(output_file, 'w') as f:
+        f.create_dataset("segments", data=segments)
+
+if __name__ == "__main__":
+    fire.Fire(dump)

cli/simulate_pixels.py

@@ -90,17 +90,17 @@ def run_simulation(input_filename,
     RangePush("slicing_and_swapping")
     tracks = tracks[:n_tracks]
 
-    y_start = np.copy(tracks['y_start'] )
-    y_end = np.copy(tracks['y_end'])
-    y = np.copy(tracks['y'])
+    x_start = np.copy(tracks['x_start'] )
+    x_end = np.copy(tracks['x_end'])
+    x = np.copy(tracks['x'])
 
-    tracks['y_start'] = np.copy(tracks['z_start'])
-    tracks['y_end'] = np.copy(tracks['z_end'])
-    tracks['y'] = np.copy(tracks['z'])
+    tracks['x_start'] = np.copy(tracks['z_start'])
+    tracks['x_end'] = np.copy(tracks['z_end'])
+    tracks['x'] = np.copy(tracks['z'])
 
-    tracks['z_start'] = y_start
-    tracks['z_end'] = y_end
-    tracks['z'] = y
+    tracks['z_start'] = x_start
+    tracks['z_end'] = x_end
+    tracks['z'] = x
     RangePop()
 
     TPB = 256
@@ -124,39 +124,51 @@ def run_simulation(input_filename,
     RangePop()
     end_drifting = time()
     print(f" {end_drifting-start_drifting:.2f} s")
-    step = 200
+    step = 1
     adc_tot_list = cp.empty((0,fee.MAX_ADC_VALUES))
     adc_tot_ticks_list = cp.empty((0,fee.MAX_ADC_VALUES))
     MAX_TRACKS_PER_PIXEL = 5
-    backtracked_id_tot = cp.empty((0,fee.MAX_ADC_VALUES,MAX_TRACKS_PER_PIXEL))
+    backtracked_id_tot = cp.empty((0,fee.MAX_ADC_VALUES,MAX_TRACKS_PER_PIXEL,2))
     unique_pix_tot = cp.empty((0,2))
     tot_events = 0
-
+
+    tot_evids = np.unique(tracks['eventID'])
     # We divide the sample in portions that can be processed by the GPU
     tracks_batch_runtimes = []
-    for itrk in tqdm(range(0, tracks.shape[0], step), desc='Simulating pixels...'):
+    for itrk in tqdm(range(0, 5, step), desc='Simulating pixels...'):
         start_tracks_batch = time()
-        selected_tracks = tracks[itrk:itrk+step]
+        first_event = tot_evids[itrk]
+        last_event = tot_evids[min(itrk+step, tot_evids.shape[0]-1)]
+
+        if first_event == last_event:
+            last_event += 1
+
+        evt_tracks = tracks[(tracks['eventID']>=first_event) & (tracks['eventID']<last_event)]
+        selected_tracks = evt_tracks[:700]
 
         RangePush("event_id_map")
         # Here we build a map between tracks and event IDs
-        unique_eventIDs = cp.unique(selected_tracks['eventID'])
-        event_id_map = cp.searchsorted(unique_eventIDs,cp.asarray(selected_tracks['eventID']))
+        event_ids = selected_tracks['eventID']
+        unique_eventIDs = np.unique(event_ids)
+        event_id_map = np.searchsorted(unique_eventIDs,event_ids)
         RangePop()
 
         # We find the pixels intersected by the projection of the tracks on
         # the anode plane using the Bresenham's algorithm. We also take into
         # account the neighboring pixels, due to the transverse diffusion of the charges.
         RangePush("pixels_from_track")
-        longest_pix = ceil(max(selected_tracks["dx"])/consts.pixel_size[0])
-        max_radius = ceil(max(selected_tracks["tran_diff"])*5/consts.pixel_size[0])
+        longest_pix = ceil(max(selected_tracks["dx"])/consts.pixel_pitch)
+        max_radius = ceil(max(selected_tracks["tran_diff"])*5/consts.pixel_pitch)
         MAX_PIXELS = int((longest_pix*4+6)*max_radius*1.5)
         MAX_ACTIVE_PIXELS = int(longest_pix*1.5)
         active_pixels = cp.full((selected_tracks.shape[0], MAX_ACTIVE_PIXELS, 2), -1, dtype=np.int32)
         neighboring_pixels = cp.full((selected_tracks.shape[0], MAX_PIXELS, 2), -1, dtype=np.int32)
         n_pixels_list = cp.zeros(shape=(selected_tracks.shape[0]))
         threadsperblock = 128
         blockspergrid = ceil(selected_tracks.shape[0] / threadsperblock)
+
+        if not active_pixels.shape[1]:
+            continue
         pixels_from_track.get_pixels[blockspergrid,threadsperblock](selected_tracks,
                                                                     active_pixels,
                                                                     neighboring_pixels,
@@ -211,7 +223,7 @@ def run_simulation(input_filename,
         blockspergrid_y = ceil(signals.shape[1] / threadsperblock[1])
         blockspergrid_z = ceil(signals.shape[2] / threadsperblock[2])
         blockspergrid = (blockspergrid_x, blockspergrid_y, blockspergrid_z)
-        pixels_signals = cp.zeros((len(unique_pix), len(consts.time_ticks)*len(unique_eventIDs)*2))
+        pixels_signals = cp.zeros((len(unique_pix), len(consts.time_ticks)*len(unique_eventIDs)*3))
         detsim.sum_pixel_signals[blockspergrid,threadsperblock](pixels_signals,
                                                                 signals,
                                                                 track_starts,
@@ -220,7 +232,7 @@ def run_simulation(input_filename,
 
         RangePush("get_adc_values")
         # Here we simulate the electronics response (the self-triggering cycle) and the signal digitization
-        time_ticks = cp.linspace(0, len(unique_eventIDs)*consts.time_interval[1]*2, pixels_signals.shape[1]+1)
+        time_ticks = cp.linspace(0, len(unique_eventIDs)*consts.time_interval[1]*3, pixels_signals.shape[1]+1)
         integral_list = cp.zeros((pixels_signals.shape[0], fee.MAX_ADC_VALUES))
         adc_ticks_list = cp.zeros((pixels_signals.shape[0], fee.MAX_ADC_VALUES))
         TPB = 128
@@ -231,7 +243,7 @@ def run_simulation(input_filename,
                                     time_ticks,
                                     integral_list,
                                     adc_ticks_list,
-                                    consts.time_interval[1]*2*tot_events,
+                                    consts.time_interval[1]*3*tot_events,
                                     rng_states)
         adc_list = fee.digitize(integral_list)
         RangePop()
@@ -248,7 +260,7 @@ def run_simulation(input_filename,
         # Here we backtrack the ADC counts to the Geant4 tracks
         TPB = 128
         BPG = ceil(adc_list.shape[0] / TPB)
-        backtracked_id = cp.full((adc_list.shape[0], adc_list.shape[1], MAX_TRACKS_PER_PIXEL), -1)
+        backtracked_id = cp.full((adc_list.shape[0], adc_list.shape[1], MAX_TRACKS_PER_PIXEL, 2), -1)
         detsim.backtrack_adcs[BPG,TPB](selected_tracks,
                                        adc_list,
                                        adc_ticks_list,