Enable_pai flag in combination w/ gdml geometry and custom E-Field

sroux · January 21, 2023, 6:01pm

Dear Allpix Developers,

I am currently working with Allpix^2 for the first time (and Corryvreckan, for that matter:)) and am having some issues with the current configuration. I’d be very grateful if you’d share your insights with me.

The problem is the following:

setup of three ALPIDE detector planes that have a custom E-Field in my telescope. the config files are attached (omitting the custom field, I am not allowed to distribute it).

This is the setup roughly:

(the beam is replaced by a square Sr-90 source w/ opening angle 180deg)

everything seems to work well and relatively fast with my current configuration, as long as the enable_pai flag in [DepositionGeant4] in the master config is set to false
collimator in the setup was implemented via the gdml file
as soon as enable_pai flag (which is crucial for my analysis) to true, the simulation becomes impossibly slow

Can it be that the PAI model is somehow incompatible with a gdml- type material present in the setup?

If they indeed are not compatible, can you propose any alternatives regading the implementation of the collimator?

Thank you!

gdml:

<gdml xsi:noNamespaceSchemaLocation="http://service-spi.web.cern.ch/service-spi/app/releases/GDML/schema/gdml.xsd">
<define>
<rotation name="alignSurfXY" z="0.001"/>
<rotation name="alignSurfX" z="0.01"/>
<rotation name="alignSurfY" z="0.01"/>
</define>
<materials>
<material Z="13" name="G4_Al" state="solid">
<T unit="K" value="293.15"/>
<MEE unit="eV" value="166"/>
<D unit="g/cm3" value="2.699"/>
<atom unit="g/mole" value="26.9815"/>
</material>
<isotope N="204" Z="82" name="Pb204">
<atom unit="g/mole" value="203.973"/>
</isotope>
<isotope N="206" Z="82" name="Pb206">
<atom unit="g/mole" value="205.974"/>
</isotope>
<isotope N="207" Z="82" name="Pb207">
<atom unit="g/mole" value="206.976"/>
</isotope>
<isotope N="208" Z="82" name="Pb208">
<atom unit="g/mole" value="207.977"/>
</isotope>
<element name="Pb">
<fraction n="0.014" ref="Pb204"/>
<fraction n="0.241" ref="Pb206"/>
<fraction n="0.221" ref="Pb207"/>
<fraction n="0.524" ref="Pb208"/>
</element>
<material name="G4_Pb" state="solid">
<T unit="K" value="293.15"/>
<MEE unit="eV" value="823"/>
<D unit="g/cm3" value="11.35"/>
<fraction n="1" ref="Pb"/>
</material>
<isotope N="14" Z="7" name="N14">
<atom unit="g/mole" value="14.0031"/>
</isotope>
<isotope N="15" Z="7" name="N15">
<atom unit="g/mole" value="15.0001"/>
</isotope>
<element name="Nitrogen">
<fraction n="0.99632" ref="N14"/>
<fraction n="0.00368" ref="N15"/>
</element>
<isotope N="16" Z="8" name="O16">
<atom unit="g/mole" value="15.9949"/>
</isotope>
<isotope N="17" Z="8" name="O17">
<atom unit="g/mole" value="16.9991"/>
</isotope>
<isotope N="18" Z="8" name="O18">
<atom unit="g/mole" value="17.9992"/>
</isotope>
<element name="Oxygen">
<fraction n="0.99757" ref="O16"/>
<fraction n="0.00038" ref="O17"/>
<fraction n="0.00205" ref="O18"/>
</element>
<material name="Air20" state="gas">
<T unit="K" value="293"/>
<MEE unit="eV" value="85.7030667332999"/>
<D unit="g/cm3" value="0.001205"/>
<fraction n="0.7" ref="Nitrogen"/>
<fraction n="0.3" ref="Oxygen"/>
</material>
</materials>
<solids>
<box lunit="mm" name="solidAlCollimator0" x="300" y="100" z="10"/>
<tube aunit="deg" deltaphi="360" lunit="mm" name="solidAlCollHole" rmax="1" rmin="0" startphi="0" z="20"/>
<subtraction name="solidAlCollimator">
<first ref="solidAlCollimator0"/>
<second ref="solidAlCollHole"/>
</subtraction>
<tube aunit="deg" deltaphi="360" lunit="mm" name="solidLeadCollimator" rmax="30" rmin="2.5" startphi="0" z="10"/>
<tube aunit="deg" deltaphi="360" lunit="mm" name="solidSrContainer" rmax="30" rmin="20" startphi="0" z="100"/>
<box lunit="mm" name="World" x="400" y="400" z="300"/>
</solids>
<structure>
<volume name="logicAlCollimator">
<materialref ref="G4_Al"/>
<solidref ref="solidAlCollimator"/>
</volume>
<volume name="logicLeadCollimator">
<materialref ref="G4_Pb"/>
<solidref ref="solidLeadCollimator"/>
</volume>
<volume name="logicSrContainer">
<materialref ref="G4_Pb"/>
<solidref ref="solidSrContainer"/>
</volume>
<volume name="World">
<materialref ref="Air20"/>
<solidref ref="World"/>
<physvol copynumber="1" name="AlCollimator">
<volumeref ref="logicAlCollimator"/>
<position name="logicAlCollimator_pos" unit="mm" x="0" y="0" z="0"/>
<rotationref ref="alignSurfXY"/>
</physvol>
<physvol copynumber="1" name="LeadCollimator">
<volumeref ref="logicLeadCollimator"/>
<position name="LeadCollimator_pos" unit="mm" x="0" y="0" z="-10"/>
<rotationref ref="alignSurfX"/>
</physvol>
<physvol copynumber="1" name="SrContainer">
<volumeref ref="logicSrContainer"/>
<position name="SrContainer_pos" unit="mm" x="0" y="0" z="-65"/>
<rotationref ref="alignSurfY"/>
</physvol>
</volume>
</structure>
<setup name="Default" version="1.0">
<world ref="World"/>
</setup>
</gdml>

master config:

Allpix]
number_of_events = 10000
model_paths = "."
log_level = "WARNING"
detectors_file = "detector.conf"
multithreading = true
workers = 16
root_file = "modules_run95_customelectricfield_nominal.root"



[GeometryBuilderGeant4]


[DepositionGeant4]
physics_list = FTFP_BERT_LIV
enable_pai = true
particle_type = "Sr90"
source_energy = 0MeV
source_position = 0um 0um 0cm
source_type = "square"
square_side = 1mm
square_angle = 180deg
max_step_length = 10.0um
output_plots = true
#max_step_length = 20.0um


#physics_list = FTFP_BERT_EMY
#particle_type = "Sr90"
#source_energy = 0eV
#source_position = 0um 0um 0cm
#source_type = "point"
#max_step_length = 20.0um
#output_plots = true


################################################
#### custom electric field ###########################
###############################################

#Not for distribution

[ElectricFieldReader]
model = "custom"
log_level = "WARNING"

#############################################################################
###### 100 um thick alpide, but field is only up to 25 um              #######
#############################################################################

field_parameters = 1V/cm, 1V/cm, 1V/cm

################################################
#### other electric field reader parameters ####
################################################

### depletion_depth should match with the epitax_depth_z of alpide.py
depletion_depth = 25um
output_plots = 1
output_plots_project = y
output_plots_single_pixel = true

[GenericPropagation]
temperature = 293K
charge_per_step = 5
timestep_min = 0.5ps
timestep_max = 0.5ns
#timestep_start = 1ps
#spatial_precision = 0.1nm
integration_time = 20ns

#output_linegraphs = true
#output_plots_step = 100ps
#output_plots_align_pixels = true
#output_plots_use_pixel_units = true


[SimpleTransfer]
output_plots = 0

[DefaultDigitizer]
log_level = "WARNING"
electronics_noise  = 10e
threshold          = 120e
threshold_smearing = 5e
output_plots = 0

[DetectorHistogrammer]
log_level = "WARNING"
max_cluster_charge=4ke
output_plots = 1

[ROOTObjectWriter]
file_name = "data_run95_customelectricfield_nominal.root"
exclude = DepositedCharge, PropagatedCharge
output_plots = 0

geometry file:

[mydetector0]
type = "alpide"
position = 0um 0um 129mm
orientation_mode = "xyz"
orientation = 0deg 0deg 0deg

[mydetector1]
type = "alpide"
position = 0um 0um 84.64mm
orientation_mode = "xyz"
orientation = 0deg 0deg 0deg

[mydetector2]
type = "alpide"
position = 0um 0um 41mm
orientation_mode = "xyz"
orientation = 0deg 0deg 0deg

[collimator]
role = "passive"
type = "gdml"
#file_name = "pb_al_colimator.gdml"
#file_name = "assembly.gdml"
file_name = "al_collimator_pb_container.gdml"
position = 0cm 0cm 16mm
orientation = 0 0 0

alpide detector specification:

#Reference https://doi.org/10.1016/j.nima.2016.05.016
# Several dimensions taken from reference, some other are estimates for a realistic assembly
geometry = "pixel"
type = "monolithic"

number_of_pixels = 1024 512
pixel_size = 29.24um 26.88um

sensor_material = SILICON
### sensor_thickness should match with the sensor_depth_z from alpide.py
sensor_thickness = 39um
sensor_excess_top = 1200um
sensor_excess_bottom = 200um
sensor_excess_left = 30um
sensor_excess_right = 30um

chip_thickness = 25um

# PCB chip board with cut-out for the sensor
[support]
thickness = 1.6mm
size = 70mm 100mm
offset = 0 30.0mm
hole_size = 32mm 15mm
hole_offset = 0 -30mm
material = "g10"

# Kapton foils for shielding the sensor, 25um on both sides, see e.g. EUDET-type beam telescopes
[support]
thickness = 25um
size = 40mm 25mm
location = "absolute"
offset = 0 0 7.5mm
material = "kapton"
[support]
thickness = 25um
size = 40mm 25mm
location = "absolute"
offset = 0 0mm -7.5mm
material = "kapton"

simonspa · January 24, 2023, 6:06am

Hi @sroux

just a quick heads-up: we have looked into this and it is neither a problem with GDML nor with your custom e-field, but a bug in Geant4 that we have encountered previously (an issue in tracking when particles are going backwards in time slowly and essentially block the entire simulation).

With your help (the PAI switch) we might be able to give the Geant4 developers a better idea where to look for issues, and in the meantime we’re trying to figure out a short-term solution for you.

Stay tuned.

Simon

hwennlof · January 24, 2023, 2:49pm

Hi,

As Simon said, it does indeed seem to be a Geant4 problem, and specifically PAI. We’re trying to figure it out, and are in contact with Geant4 developers concerning this!
Looking at your setup with the visualiser, it seems that most decay particles end up in the collimator, and that’s where the simulation gets stuck. The opening angle of the source sadly has no effect when the particle used is a radioactively decaying one, as those decays are always isotropic, so a lot more than expected end up in the collimator.

If you don’t need the strontium-90 spectrum specifically, I would suggest running without collimator and with a beam of electrons with a given energy as a source for now, with a size given by the collimator hole.
If you do need the strontium-90, you can construct a spectrum using the General Particle Source (i.e. macro as source type in Allpix Squared). I attach an example of such a macro here (with no electrons with an energy below 100 keV generated; also note that the default unit for the General Particle Source is allegedly MeV when giving histograms to pull from, but the results I get with this indicate that the given energies are a bit low). By tweaking the maxtheta and mintheta, this can act as a collimator would. When aiming straight for the sensors, the passive material given by the PCB can also be removed without impacting the simulation results in any significant way. This can further reduce the risk of things getting stuck for a long time.

This is just a temporary workaround however, hopefully the PAI issue will be fixed soon!

Kind regards,
Håkan

sourceMacroExample.mac (28.9 KB)

sroux · January 31, 2023, 12:02pm

Dear Simon and Hakan,

Thanks a lot for your answers- it is unfortunate that we encounter this bug, for the meantime I will try to work w/ the source macro, as recommended.

We do indeed need the Sr-90 spectrum specifically.

And yes, you are right - in our setup it is inevitable that many electrons hit the passive materials, which made things very slow.

Best,
Sophie