From 8a9a448f374cfb6c0daad40bbf8b94c64f36cf39 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9?= Date: Tue, 10 Jun 2025 19:43:10 +0200 Subject: [PATCH 01/26] added vacuum as a possible material for generic collimators (i.e. for empty space in FLUKA) --- xcoll/scattering_routines/fluka/generic_prototype.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/xcoll/scattering_routines/fluka/generic_prototype.py b/xcoll/scattering_routines/fluka/generic_prototype.py index e8010240c..f75efe3a8 100644 --- a/xcoll/scattering_routines/fluka/generic_prototype.py +++ b/xcoll/scattering_routines/fluka/generic_prototype.py @@ -23,7 +23,8 @@ def xcoll_to_fluka_material(material): "c": " CARBON", "si": " SILICON", "cu": " COPPER", - "mogr": "AC150GPH" + "mogr": "AC150GPH", + "vacuum": "VACUUM" } if material not in material_dict.keys(): From 2faff4bd31ec9c4385cc98651f1772b893a99f43 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9?= Date: Tue, 10 Jun 2025 19:50:01 +0200 Subject: [PATCH 02/26] Preliminary version of beam-beam interaction in Xsyute --- xcoll/scattering_routines/fluka/includes.py | 42 ++++++++++++++++++--- 1 file changed, 37 insertions(+), 5 deletions(-) diff --git a/xcoll/scattering_routines/fluka/includes.py b/xcoll/scattering_routines/fluka/includes.py index 92daa1f68..c74d4ceab 100644 --- a/xcoll/scattering_routines/fluka/includes.py +++ b/xcoll/scattering_routines/fluka/includes.py @@ -43,10 +43,13 @@ def get_include_files(particle_ref, include_files=[], *, verbose=True, lower_mom return_leptons=None, return_neutrinos=None, return_protons=None, return_neutrons=None, return_ions=None, return_exotics=None, return_all=None, return_neutral=None, use_crystals=False, **kwargs): + + bb_int = kwargs.get('bb_int', None) + this_include_files = include_files.copy() # Required default include files if 'include_settings_beam.inp' not in [file.name for file in this_include_files]: - this_include_files.append(_beam_include_file(particle_ref)) + this_include_files.append(_beam_include_file(particle_ref, bb_int=bb_int)) if 'include_settings_physics.inp' in [file.name for file in this_include_files]: if photon_lower_momentum_cut is not None: raise ValueError("Physics include file already provided. Cannot change " @@ -217,10 +220,11 @@ def _assignmat_include_file(): return filename -def _beam_include_file(particle_ref): +def _beam_include_file(particle_ref, bb_int=0): filename = FsPath("include_settings_beam.inp").resolve() pdg_id = particle_ref.pdg_id[0] momentum_cut = particle_ref.p0c[0] / 1e9 * 1.05 + hi_prope = "*" if pdg_id in fluka_names: name = fluka_names[pdg_id] @@ -233,7 +237,6 @@ def _beam_include_file(particle_ref): raise ValueError(f"Reference particle {get_name_from_pdg_id(pdg_id)} not " + "supported by FLUKA.") beam = f"BEAM {format_fluka_float(momentum_cut)}{50*' '}{name}" - template = f"""\ ****************************************************************************** * BEAM SETTINGS * @@ -245,12 +248,39 @@ def _beam_include_file(particle_ref): * BEAMPOS * -* Only asking for loss map and touches map as in Sixtrack +""" + if bb_int: + bb_int_template = f"""\ +* ..+....1....+....2....+....3....+....4....+....5....+....6....+....7.. +* W(1):interaction type, W(2):Index IR, W(3):Length of IR[cm], W(4):SigmaZ +SOURCE 1. 1. 0. 0. +SOURCE 89. 90. 91. 0. & +* ..+....1....+....2....+....3....+....4....+....5....+....6....+....7.. +*USRICALL 0.511232 -0.026532 0.0148447 0.0020532 2.3e-6 BBCOL_H +*USRICALL 0.59894741-0.0002269-0.0001212-0.00035714.8293E-10 BBCOL_H +* USRICALL 0.511232 -0.026532 0.0148447 0.0020532 2.3e-6 BBCOL_H +USRICALL {format_fluka_float(bb_int['theta2'])}{format_fluka_float(bb_int['xs'])}{format_fluka_float(bb_int['sigma_p_x2'])}{format_fluka_float(bb_int['sigma_p_y2'])}{Z:8}.0{A:8}.0BBEAMCOL +USRICALL {format_fluka_float(bb_int['betx'])}{format_fluka_float(bb_int['alfx'])}{format_fluka_float(bb_int['dx'])}{format_fluka_float(bb_int['dpx'])}{format_fluka_float(bb_int['rms_emx'])} BBCOL_H +USRICALL {format_fluka_float(bb_int['bety'])}{format_fluka_float(bb_int['alfy'])}{format_fluka_float(bb_int['dy'])}{format_fluka_float(bb_int['dpy'])}{format_fluka_float(bb_int['rms_emy'])} BBCOL_V +* ..+....1....+....2....+....3....+....4....+....5....+....6....+....7.. +*USRICALL #OFFSET BBCOL_O +USRICALL 0.0 0.0 0.0 0.0 0.0 0.0BBCOL_O +*USRICALL #SIGDPP BBCOL_L +USRICALL 0.0 0.0 0.0 0.0 0.0 0.0BBCOL_L +* +""" + template += bb_int_template + else: + template_source += """\ +* Only asking for loss map and touches map as in Xsuite * ..+....1....+....2....+....3....+....4....+....5....+....6....+....7.. SOURCE 87. 88. 1. SOURCE 89. 90. 91. 0. -1. 10.& -*SOURCE 0.0 0.0 97.0 1.0 96.0 1.0&& +SOURCE 0.0 0.0 97.0 1.0 96.0 1.0&& """ + + template += template_source + with filename.open('w') as fp: fp.write(template) return filename @@ -311,6 +341,8 @@ def _physics_include_file(*, verbose, lower_momentum_cut, photon_lower_momentum_ PHYSICS 1.D+5 1.D+5 1.D+5 1.D+5 1.D+5 1.D+5PEATHRES PHYSICS 1. 0.005 0.15 2.0 2.0 3.0IONSPLIT PHYSICS 2. EM-DISSO +* beam-beam collisoins +PHYSICS 8000.0 LIMITS *THRESHOL 0.0 0.0 8000.0 8000.0 0.0 0.0 """ with filename.open('w') as fp: From 796470d83e310581e052a87b20a5946b1c842a7e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9?= Date: Wed, 11 Jun 2025 09:45:13 +0200 Subject: [PATCH 03/26] minor bug in includes for bbint and source cards --- xcoll/scattering_routines/fluka/includes.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/xcoll/scattering_routines/fluka/includes.py b/xcoll/scattering_routines/fluka/includes.py index c74d4ceab..0a7a2841e 100644 --- a/xcoll/scattering_routines/fluka/includes.py +++ b/xcoll/scattering_routines/fluka/includes.py @@ -250,7 +250,7 @@ def _beam_include_file(particle_ref, bb_int=0): * """ if bb_int: - bb_int_template = f"""\ + template += f"""\ * ..+....1....+....2....+....3....+....4....+....5....+....6....+....7.. * W(1):interaction type, W(2):Index IR, W(3):Length of IR[cm], W(4):SigmaZ SOURCE 1. 1. 0. 0. @@ -269,9 +269,8 @@ def _beam_include_file(particle_ref, bb_int=0): USRICALL 0.0 0.0 0.0 0.0 0.0 0.0BBCOL_L * """ - template += bb_int_template else: - template_source += """\ + template += f"""\ * Only asking for loss map and touches map as in Xsuite * ..+....1....+....2....+....3....+....4....+....5....+....6....+....7.. SOURCE 87. 88. 1. @@ -279,7 +278,6 @@ def _beam_include_file(particle_ref, bb_int=0): SOURCE 0.0 0.0 97.0 1.0 96.0 1.0&& """ - template += template_source with filename.open('w') as fp: fp.write(template) From b713b34651e19255dcbcf5993929111f568408f1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9?= Date: Fri, 27 Jun 2025 16:53:33 +0200 Subject: [PATCH 04/26] Reminder that the beam-beam interaction routine currently uses this 'empty collimators'. Maybe including dimensions for generic collimation would be desirable --- xcoll/scattering_routines/fluka/generic_prototype.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/xcoll/scattering_routines/fluka/generic_prototype.py b/xcoll/scattering_routines/fluka/generic_prototype.py index f75efe3a8..26b097ddd 100644 --- a/xcoll/scattering_routines/fluka/generic_prototype.py +++ b/xcoll/scattering_routines/fluka/generic_prototype.py @@ -213,6 +213,8 @@ def _body_file(fedb_tag, length, width, height, **kwargs): template_tank = f"""\ RPP {fedb_tag}_T -28 28 -28 28 -{length*100 + 5} {length*100 + 5} RPP {fedb_tag}_I -28 28 -28 28 -{length*100 + 5} {length*100 + 5} +*RPP {fedb_tag}_T -28 28 -28 28 -{length*100 + 1e-12} {length*100 + 1e-12} +*RPP {fedb_tag}_I -28 28 -28 28 -{length*100 + 1e-12} {length*100 + 1e-12} """ tank_file = _write_file("bodies", f"generic_{fedb_tag}_T.bodies", template_tank) From 5c070479ade243783c9c7183b2aab85aa50f37aa Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9?= Date: Fri, 27 Jun 2025 17:09:51 +0200 Subject: [PATCH 05/26] added fluka routine for Beam-beam interaction. TODO: Add example --- xcoll/scattering_routines/fluka/includes.py | 29 ++++++++++++--------- 1 file changed, 17 insertions(+), 12 deletions(-) diff --git a/xcoll/scattering_routines/fluka/includes.py b/xcoll/scattering_routines/fluka/includes.py index 33c877335..4248cb1b3 100644 --- a/xcoll/scattering_routines/fluka/includes.py +++ b/xcoll/scattering_routines/fluka/includes.py @@ -219,9 +219,9 @@ def _assignmat_include_file(): raise ValueError(f"Crystal {name} has multiple positions in the prototypes file. " + "(too many dependant aseemblies).") pos = crystal.dependant_assemblies[0].fluka_position - pos1 = format_fluka_float(pos[3]) + pos1 = format_fluka_float(pos[3]-50.0+1e-4) pos2 = format_fluka_float(pos[4]) - pos3 = format_fluka_float(pos[5]) + pos3 = format_fluka_float(pos[5]+1e-4) template += f"""\ * ..+....1....+....2....+....3....+....4....+....5....+....6....+....7.. CRYSTAL {name:>8}{format_fluka_float(bang)}{format_fluka_float(l)} 0.0 0.0 300.0 110 @@ -235,7 +235,7 @@ def _assignmat_include_file(): return filename -def _beam_include_file(particle_ref, bb_int=0): +def _beam_include_file(particle_ref, bb_int=False): filename = FsPath("include_settings_beam.inp").resolve() pdg_id = particle_ref.pdg_id[0] momentum_cut = particle_ref.p0c[0] / 1e9 * 1.05 @@ -267,20 +267,25 @@ def _beam_include_file(particle_ref, bb_int=0): if bb_int: template += f"""\ * ..+....1....+....2....+....3....+....4....+....5....+....6....+....7.. -* W(1):interaction type, W(2):Index IR, W(3):Length of IR[cm], W(4):SigmaZ -SOURCE 1. 1. 0. 0. +* W(1): interaction type, W(2): Index IR, W(3): Length of IR[cm], W(4): SigmaZ +* Interaction type: 1.0 (Inelastic), 10.0 (Elastic), 100.0 (EMD) +* SigmaZ: sigma_z (cm) for the Gaussian sampling of the, collision position around the center of the insertion +SOURCE {format_fluka_float(bb_int['int_type'])} 1. 0. 0. SOURCE 89. 90. 91. 0. & * ..+....1....+....2....+....3....+....4....+....5....+....6....+....7.. -*USRICALL 0.511232 -0.026532 0.0148447 0.0020532 2.3e-6 BBCOL_H -*USRICALL 0.59894741-0.0002269-0.0001212-0.00035714.8293E-10 BBCOL_H -* USRICALL 0.511232 -0.026532 0.0148447 0.0020532 2.3e-6 BBCOL_H -USRICALL {format_fluka_float(bb_int['theta2'])}{format_fluka_float(bb_int['xs'])}{format_fluka_float(bb_int['sigma_p_x2'])}{format_fluka_float(bb_int['sigma_p_y2'])}{Z:8}.0{A:8}.0BBEAMCOL +* W(1): Theta2-> Polar angle (rad) between b2 and -z direction. +* W(2): Azimuthal angle (deg) defining the crossing plane. +* W(3): Sigma_x for Gaussian sampling in b2, dir (x', y', 1) wrt to px +* W(4): Sigma_y for Gaussian sampling in b2, dir (x', y', 1) wrt to py +* W(5): Z b2 +* W(6): A b2 +USRICALL {format_fluka_float(bb_int['theta2'])}{format_fluka_float(bb_int['xs'])}{format_fluka_float(bb_int['sigma_p_x2'])}{format_fluka_float(bb_int['sigma_p_y2'])}{format_fluka_float(bb_int['Z'])}{format_fluka_float(bb_int['A'])}BBEAMCOL USRICALL {format_fluka_float(bb_int['betx'])}{format_fluka_float(bb_int['alfx'])}{format_fluka_float(bb_int['dx'])}{format_fluka_float(bb_int['dpx'])}{format_fluka_float(bb_int['rms_emx'])} BBCOL_H USRICALL {format_fluka_float(bb_int['bety'])}{format_fluka_float(bb_int['alfy'])}{format_fluka_float(bb_int['dy'])}{format_fluka_float(bb_int['dpy'])}{format_fluka_float(bb_int['rms_emy'])} BBCOL_V -* ..+....1....+....2....+....3....+....4....+....5....+....6....+....7.. -*USRICALL #OFFSET BBCOL_O +* * ..+....1....+....2....+....3....+....4....+....5....+....6....+....7.. +* *USRICALL #OFFSET BBCOL_O USRICALL 0.0 0.0 0.0 0.0 0.0 0.0BBCOL_O -*USRICALL #SIGDPP BBCOL_L +* *USRICALL #SIGDPP BBCOL_L USRICALL 0.0 0.0 0.0 0.0 0.0 0.0BBCOL_L * """ From 9d5a24458f88dbded453437c10eba5ac87c06391 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9?= Date: Fri, 27 Jun 2025 19:03:18 +0200 Subject: [PATCH 06/26] there is channeling, however there are some problems in the definition of the shape of the crystal. TO CORRECT --- xcoll/scattering_routines/fluka/generic_prototype.py | 2 +- xcoll/scattering_routines/fluka/includes.py | 3 ++- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/xcoll/scattering_routines/fluka/generic_prototype.py b/xcoll/scattering_routines/fluka/generic_prototype.py index 26b097ddd..d54fdcc3b 100644 --- a/xcoll/scattering_routines/fluka/generic_prototype.py +++ b/xcoll/scattering_routines/fluka/generic_prototype.py @@ -256,7 +256,7 @@ def _material_file(fedb_tag, material, **kwargs): def _crystal_body_file(fedb_tag, length, bending_radius, width, height, **kwargs): template_body = f"""\ -RPP {fedb_tag}_B 0.0 {width*(100+10)} -{height*(100+10)/2} {height*(100+10)/2} 0.00001 {length*(100+20)} +RPP {fedb_tag}_B 0.0 {width*(100+10)} -{height*(100+10)/2} {height*(100+10)/2} -{length*(100+20)} {length*(100+20)} ZCC {fedb_tag}Z1 {bending_radius*100} 0.0 {bending_radius*100} ZCC {fedb_tag}Z2 {bending_radius*100} 0.0 {bending_radius*100-width*100} PLA {fedb_tag}P1 1.0 0.0 {np.cos(length/bending_radius)/np.sin(length/bending_radius)} {bending_radius*100} 0.0 0.0 diff --git a/xcoll/scattering_routines/fluka/includes.py b/xcoll/scattering_routines/fluka/includes.py index 4248cb1b3..569b16a90 100644 --- a/xcoll/scattering_routines/fluka/includes.py +++ b/xcoll/scattering_routines/fluka/includes.py @@ -221,7 +221,8 @@ def _assignmat_include_file(): pos = crystal.dependant_assemblies[0].fluka_position pos1 = format_fluka_float(pos[3]-50.0+1e-4) pos2 = format_fluka_float(pos[4]) - pos3 = format_fluka_float(pos[5]+1e-4) + pos3 = format_fluka_float(pos[5]-(l)-0.04+1e-4) + #pos3 = format_fluka_float(pos[5]+1e-4) template += f"""\ * ..+....1....+....2....+....3....+....4....+....5....+....6....+....7.. CRYSTAL {name:>8}{format_fluka_float(bang)}{format_fluka_float(l)} 0.0 0.0 300.0 110 From 36371a17f14670d301024cc224d3bdcb3220dcec Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9?= Date: Tue, 1 Jul 2025 12:41:24 +0200 Subject: [PATCH 07/26] Crystal working.. TODO: Tilts, test, examples --- xcoll/scattering_routines/fluka/engine.py | 9 +++++++-- xcoll/scattering_routines/fluka/generic_prototype.py | 8 +++++--- xcoll/scattering_routines/fluka/includes.py | 9 ++++----- 3 files changed, 16 insertions(+), 10 deletions(-) diff --git a/xcoll/scattering_routines/fluka/engine.py b/xcoll/scattering_routines/fluka/engine.py index bbb1eaf08..dd545eb5c 100644 --- a/xcoll/scattering_routines/fluka/engine.py +++ b/xcoll/scattering_routines/fluka/engine.py @@ -246,8 +246,13 @@ def _match_input_file(self): continue self._assert_element(ee) ee.fluka_id = input_dict[name]['fluka_id'] - ee.length_front = (input_dict[name]['length'] - ee.length)/2 - ee.length_back = (input_dict[name]['length'] - ee.length)/2 + if ee.assembly.is_crystal: + ee.length_front = (input_dict[name]['length'])/2 + ee.length_back = (input_dict[name]['length'] - 2*ee.length)/2 + else: + ee.length_front = (input_dict[name]['length'] - 2*ee.length)/2 + ee.length_back = (input_dict[name]['length'] - 2*ee.length)/2 + jaw = input_dict[name]['jaw'] if jaw is not None and not hasattr(jaw, '__iter__'): jaw = [jaw, -jaw] diff --git a/xcoll/scattering_routines/fluka/generic_prototype.py b/xcoll/scattering_routines/fluka/generic_prototype.py index d54fdcc3b..36152c623 100644 --- a/xcoll/scattering_routines/fluka/generic_prototype.py +++ b/xcoll/scattering_routines/fluka/generic_prototype.py @@ -257,9 +257,10 @@ def _material_file(fedb_tag, material, **kwargs): def _crystal_body_file(fedb_tag, length, bending_radius, width, height, **kwargs): template_body = f"""\ RPP {fedb_tag}_B 0.0 {width*(100+10)} -{height*(100+10)/2} {height*(100+10)/2} -{length*(100+20)} {length*(100+20)} -ZCC {fedb_tag}Z1 {bending_radius*100} 0.0 {bending_radius*100} -ZCC {fedb_tag}Z2 {bending_radius*100} 0.0 {bending_radius*100-width*100} +YCC {fedb_tag}Z1 0.0 {bending_radius*100} {bending_radius*100} +YCC {fedb_tag}Z2 0.0 {bending_radius*100} {bending_radius*100-width*100} PLA {fedb_tag}P1 1.0 0.0 {np.cos(length/bending_radius)/np.sin(length/bending_radius)} {bending_radius*100} 0.0 0.0 +XYP {fedb_tag}P2 0.0 """ body_file = _write_file("bodies", f"generic_{fedb_tag}_B.bodies", template_body) @@ -274,10 +275,11 @@ def _crystal_body_file(fedb_tag, length, bending_radius, width, height, **kwargs def _crystal_region_file(fedb_tag, **kwargs): template_body_reg = f"""\ -{fedb_tag}_B 5 | +{fedb_tag}_B +{fedb_tag}Z1 -{fedb_tag}Z2 +{fedb_tag}P1 +{fedb_tag}_B 5 | +{fedb_tag}_B +{fedb_tag}Z1 -{fedb_tag}Z2 +{fedb_tag}P1 - {fedb_tag}P2 {fedb_tag}B2 5 | +{fedb_tag}_B +{fedb_tag}Z2 | +{fedb_tag}_B -{fedb_tag}Z1 | +{fedb_tag}_B -{fedb_tag}P1 + | +{fedb_tag}_B +{fedb_tag}P2 -{fedb_tag}Z2 """ _write_file("regions", f"generic_{fedb_tag}_B.regions", template_body_reg) diff --git a/xcoll/scattering_routines/fluka/includes.py b/xcoll/scattering_routines/fluka/includes.py index 569b16a90..cf1f831cb 100644 --- a/xcoll/scattering_routines/fluka/includes.py +++ b/xcoll/scattering_routines/fluka/includes.py @@ -217,12 +217,11 @@ def _assignmat_include_file(): raise ValueError(f"Crystal {name} has no position in the prototypes file.") elif len(crystal.dependant_assemblies) > 1: raise ValueError(f"Crystal {name} has multiple positions in the prototypes file. " - + "(too many dependant aseemblies).") + + "(too many dependant asemblies).") pos = crystal.dependant_assemblies[0].fluka_position - pos1 = format_fluka_float(pos[3]-50.0+1e-4) + pos1 = format_fluka_float(pos[3]-50.0+1e-4) # XXX Does it always work with the 50cm shift? pos2 = format_fluka_float(pos[4]) - pos3 = format_fluka_float(pos[5]-(l)-0.04+1e-4) - #pos3 = format_fluka_float(pos[5]+1e-4) + pos3 = format_fluka_float(pos[5]) template += f"""\ * ..+....1....+....2....+....3....+....4....+....5....+....6....+....7.. CRYSTAL {name:>8}{format_fluka_float(bang)}{format_fluka_float(l)} 0.0 0.0 300.0 110 @@ -360,7 +359,7 @@ def _physics_include_file(*, verbose, lower_momentum_cut, photon_lower_momentum_ PHYSICS 1.D+5 1.D+5 1.D+5 1.D+5 1.D+5 1.D+5PEATHRES PHYSICS 1. 0.005 0.15 2.0 2.0 3.0IONSPLIT PHYSICS 2. EM-DISSO -* beam-beam collisoins +* beam-beam collisions PHYSICS 8000.0 LIMITS *THRESHOL 0.0 0.0 8000.0 8000.0 0.0 0.0 """ From 89ef285d45aa79cb3466777842837022cfc1cd29 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9?= Date: Wed, 2 Jul 2025 10:25:57 +0200 Subject: [PATCH 08/26] missing capacity; Now example works --- examples/fluka_crystal.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/examples/fluka_crystal.py b/examples/fluka_crystal.py index e2d5b3cec..eeb6b6333 100644 --- a/examples/fluka_crystal.py +++ b/examples/fluka_crystal.py @@ -18,10 +18,10 @@ # context = xo.ContextCupy() # For CUDA GPUs # context = xo.ContextPyopencl() # For OpenCL GPUs -num_part = 100 +num_part = int(10_000) coll = xc.FlukaCrystal(length=0.002, material='si', bending_angle=149e-6, - width=0.002, height=0.05, side='+', jaw=1.e-12, # Hack to make it work for now; jaw should be almost zero + width=0.002, height=0.05, side='+', jaw=-1e-3, # Hack to make it work for now; jaw should be almost zero _context=context) # Connect to FLUKA @@ -34,7 +34,7 @@ px_init = np.random.uniform(low=-50.e-6, high=250.e-6, size=num_part) y_init = np.random.normal(loc=0., scale=1e-3, size=num_part) py_init = np.random.normal(loc=0., scale=5.e-6, size=num_part) -part = xp.build_particles(x=x_init, px=px_init, y=y_init, py=py_init, particle_ref=xc.fluka.engine.particle_ref, _context=context) +part = xp.build_particles(x=x_init, px=px_init, y=y_init, py=py_init, particle_ref=xc.fluka.engine.particle_ref, _context=context, _capacity=xc.fluka.engine.capacity) part_init = part.copy() From be074d15764186e4e836d7af8c89930b2a9ce04d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9?= Date: Wed, 2 Jul 2025 11:00:25 +0200 Subject: [PATCH 09/26] added width and height --- xcoll/beam_elements/fluka.py | 8 +++++++- xcoll/scattering_routines/fluka/fluka_input.py | 4 ++-- xcoll/scattering_routines/fluka/generic_prototype.py | 3 ++- 3 files changed, 11 insertions(+), 4 deletions(-) diff --git a/xcoll/beam_elements/fluka.py b/xcoll/beam_elements/fluka.py index f91397e17..5dcbdfe7e 100644 --- a/xcoll/beam_elements/fluka.py +++ b/xcoll/beam_elements/fluka.py @@ -226,6 +226,7 @@ def __init__(self, **kwargs): import xcoll as xc to_assign = {} generic = False + import pdb; pdb.set_trace() if '_xobject' not in kwargs: kwargs.setdefault('_tracking', True) kwargs.setdefault('_acc_ionisation_loss', -1.) @@ -252,6 +253,10 @@ def __init__(self, **kwargs): side = kwargs.pop('side', None) if side is None: raise ValueError('Need to provide side!') + if 'width' in kwargs: + width = kwargs.pop('width', None) + if 'height' in kwargs: + height = kwargs.pop('height', None) generic = True super().__init__(**kwargs) for key, val in to_assign.items(): @@ -262,7 +267,8 @@ def __init__(self, **kwargs): side = self._get_side_from_input(side) self.assembly = create_generic_assembly(is_crystal=True, material=material, side=side, length=self.length, - bending_radius=bending_radius) + bending_radius=bending_radius, + width=width, height=height) if not hasattr(self, '_equivalent_drift'): self._equivalent_drift = xt.Drift(length=self.length) diff --git a/xcoll/scattering_routines/fluka/fluka_input.py b/xcoll/scattering_routines/fluka/fluka_input.py index 143942beb..7e7e322fc 100644 --- a/xcoll/scattering_routines/fluka/fluka_input.py +++ b/xcoll/scattering_routines/fluka/fluka_input.py @@ -160,8 +160,8 @@ def _element_dict_to_fluka(element_dict, dump=False): offset = (ee._jaw_LU + ee._jaw_LD + ee._jaw_RU + ee._jaw_RD) / 4 tilt_1 = round(ee.tilt_L, 9) tilt_2 = round(ee.tilt_R, 9) - if abs(tilt_1) > 1.e-12 or abs(tilt_2) > 1.e-12: - raise NotImplementedError(f"Collimator {name}: Tilts are not (yet) supported in FLUKA-Xcoll!") + # if abs(tilt_1) > 1.e-12 or abs(tilt_2) > 1.e-12: + # raise NotImplementedError(f"Collimator {name}: Tilts are not (yet) supported in FLUKA-Xcoll!") if nsig is None: nsig = 1 diff --git a/xcoll/scattering_routines/fluka/generic_prototype.py b/xcoll/scattering_routines/fluka/generic_prototype.py index db3047c7c..414746c89 100644 --- a/xcoll/scattering_routines/fluka/generic_prototype.py +++ b/xcoll/scattering_routines/fluka/generic_prototype.py @@ -115,7 +115,8 @@ def _validate_kwargs(kwargs): if field not in kwargs: raise ValueError(f"Need to provide {field}!") for field, opt_value in _generic_crystal_optional_fields.items(): - kwargs.setdefault(field, opt_value) + if field not in kwargs: + kwargs.setdefault(field, opt_value) else: for field in _generic_required_fields: if field not in kwargs: From d51a611124a5af3a3aa59a0d47cb96d91fcc91ac Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9?= Date: Tue, 8 Jul 2025 09:59:46 +0200 Subject: [PATCH 10/26] bug in elements front and back length --- xcoll/scattering_routines/fluka/engine.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/xcoll/scattering_routines/fluka/engine.py b/xcoll/scattering_routines/fluka/engine.py index 77427bd9f..8fa24dda1 100644 --- a/xcoll/scattering_routines/fluka/engine.py +++ b/xcoll/scattering_routines/fluka/engine.py @@ -248,8 +248,8 @@ def _match_input_file(self): ee.length_front = (input_dict[name]['length'])/2 ee.length_back = (input_dict[name]['length'] - 2*ee.length)/2 else: - ee.length_front = (input_dict[name]['length'] - 2*ee.length)/2 - ee.length_back = (input_dict[name]['length'] - 2*ee.length)/2 + ee.length_front = (input_dict[name]['length'] - ee.length)/2 + ee.length_back = (input_dict[name]['length'] - ee.length)/2 jaw = input_dict[name]['jaw'] if jaw is not None and not hasattr(jaw, '__iter__'): From 2ec2d0825abcb9347c3a195bf65750fd070ef0d0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9?= Date: Tue, 8 Jul 2025 10:01:05 +0200 Subject: [PATCH 11/26] reincluded collimators shift when closed orbit is not None --- xcoll/scattering_routines/fluka/track.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/xcoll/scattering_routines/fluka/track.py b/xcoll/scattering_routines/fluka/track.py index 8a0d33f43..9c5938c21 100644 --- a/xcoll/scattering_routines/fluka/track.py +++ b/xcoll/scattering_routines/fluka/track.py @@ -53,7 +53,15 @@ def track(coll, particles): + "FLUKA.\nIn any case, please stop and restart the FlukaEngine now.") _drift(coll, particles, -coll.length_front) + if coll.co is not None: # FLUKA collimators are centered; need to shift + dx = coll.co[1][0] + dy = coll.co[1][1] + particles.x -= dx + particles.y -= dy track_core(coll, particles) + if coll.co is not None: + particles.x += dx + particles.y += dy _drift(coll, particles, -coll.length_back) From 92f2d1ccfc748ddb8451dc9088df32b16ca8f2dd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9?= Date: Mon, 14 Jul 2025 13:51:58 +0200 Subject: [PATCH 12/26] bug in tilts solved; offset to think of to deal you that and the closed orbit (removed for the moment) --- xcoll/beam_elements/fluka.py | 1 - xcoll/scattering_routines/fluka/fluka_input.py | 13 ++++++++----- 2 files changed, 8 insertions(+), 6 deletions(-) diff --git a/xcoll/beam_elements/fluka.py b/xcoll/beam_elements/fluka.py index 5dcbdfe7e..c9c7c7bce 100644 --- a/xcoll/beam_elements/fluka.py +++ b/xcoll/beam_elements/fluka.py @@ -226,7 +226,6 @@ def __init__(self, **kwargs): import xcoll as xc to_assign = {} generic = False - import pdb; pdb.set_trace() if '_xobject' not in kwargs: kwargs.setdefault('_tracking', True) kwargs.setdefault('_acc_ionisation_loss', -1.) diff --git a/xcoll/scattering_routines/fluka/fluka_input.py b/xcoll/scattering_routines/fluka/fluka_input.py index 7e7e322fc..caa0c2356 100644 --- a/xcoll/scattering_routines/fluka/fluka_input.py +++ b/xcoll/scattering_routines/fluka/fluka_input.py @@ -157,11 +157,13 @@ def _element_dict_to_fluka(element_dict, dump=False): elif ee.gap_R is not None: nsig = ee.gap_R half_gap = (ee._jaw_LU + ee._jaw_LD - ee._jaw_RU - ee._jaw_RD) / 4 - offset = (ee._jaw_LU + ee._jaw_LD + ee._jaw_RU + ee._jaw_RD) / 4 - tilt_1 = round(ee.tilt_L, 9) - tilt_2 = round(ee.tilt_R, 9) - # if abs(tilt_1) > 1.e-12 or abs(tilt_2) > 1.e-12: - # raise NotImplementedError(f"Collimator {name}: Tilts are not (yet) supported in FLUKA-Xcoll!") + offset = 0 # (ee._jaw_LU + ee._jaw_LD + ee._jaw_RU + ee._jaw_RD) / 4 + + tilt_1 = round(ee.tilt_L, 9) if ee.tilt_L is not None else 0.0 + tilt_2 = round(ee.tilt_R, 9) if ee.tilt_R is not None else 0.0 + + if abs(tilt_1) > 1.e-12 or abs(tilt_2) > 1.e-12: + raise NotImplementedError(f"Collimator {name}: Tilts are not (yet) supported in FLUKA-Xcoll!") if nsig is None: nsig = 1 @@ -211,6 +213,7 @@ def _fluka_builder(collimator_dict): with open('linebuilder.log', 'w') as f: with redirect_stdout(f): + # input_file, coll_dict = fb.fluka_builder(args_fb, auto_accept=True, verbose=True) input_file, coll_dict = fb.fluka_builder(args_fb, auto_accept=True) # Restore system state From 7b97ea82820f0bb8567f29b6bd27edd8386ef671 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9?= Date: Fri, 25 Jul 2025 14:46:29 +0200 Subject: [PATCH 13/26] to clean up, version with hit miss particles --- xcoll/scattering_routines/fluka/track.py | 296 +++++++++++++++++++++-- 1 file changed, 281 insertions(+), 15 deletions(-) diff --git a/xcoll/scattering_routines/fluka/track.py b/xcoll/scattering_routines/fluka/track.py index 9c5938c21..ef8c16ba7 100644 --- a/xcoll/scattering_routines/fluka/track.py +++ b/xcoll/scattering_routines/fluka/track.py @@ -21,6 +21,7 @@ def _drift(coll, particles, length): def track(coll, particles): import xcoll as xc + # Initialize ionisation loss accumulation variable if coll._acc_ionisation_loss < 0: coll._acc_ionisation_loss = 0. @@ -52,30 +53,214 @@ def track(coll, particles): + "with a value large enough to accommodate secondaries outside of " + "FLUKA.\nIn any case, please stop and restart the FlukaEngine now.") - _drift(coll, particles, -coll.length_front) + import time + + import xpart as xp + + start = time.time() + + coll_everest = xc.EverestCollimator(length=coll.length, material=xc.materials.MolybdenumGraphite) + coll_everest.jaw = coll.jaw + + part_everest = particles.copy() + coll_everest.track(part_everest) + + mask_hit_base = ( + (abs(particles.px - part_everest.px) > 1e-12) | + (abs(particles.py - part_everest.py) > 1e-12) | + (particles.pdg_id == -999999999) + ) + + end = time.time() + print(f"Tracking Everest collimator took {end - start:.4f} seconds") + + + start = time.time() + mask_lost = (particles.state != 1) & (particles.pdg_id != -999999999) + + # Final hit mask includes both + mask_hit = mask_hit_base # | mask_lost + + +# mask_hit = ((abs(particles.px - part_everest.px) > 1e-12) | (abs(particles.py - part_everest.py) >1e-12) | (particles.pdg_id == -999999999)) +# + # mask_lost = ((particles.state != 1) & (particles.pdg_id != -999999999)) + mask_miss_no_lost = ~((abs(particles.px - part_everest.px) > 1e-12) | (abs(particles.py - part_everest.py) > 1e-12)) & (particles.pdg_id != -999999999) + mask_miss = mask_miss_no_lost | mask_lost + + + hit_id = particles.particle_id[mask_hit] + miss_id = particles.particle_id[mask_miss] + is_hit = np.isin(particles.particle_id, hit_id) + is_miss = np.isin(particles.particle_id, miss_id) + + if mask_lost.any(): + particles_lost = particles.filter(mask_lost) + else: + particles_lost = None + + particles_hit = 0 + particles_miss = 0 + particles_hit = particles.filter(is_hit) + end = time.time() + print(f"Tracking Everest collimator took {end - start:.4f} seconds") + + print(f"Tracking {particles._num_active_particles} particles (FLUKA)... ", end='') + print(f"Hit: {particles_hit._num_active_particles}") + print(f"Collimator: {coll.name}") + + + # if ((particles_hit._num_active_particles == 0) & (particles_miss == 0)): + # # track_core(coll, particles) + # _drift(coll, particles, coll.length) + # return + + + # check if ~is_hit has any true value + if np.any(~is_hit): + particles_miss = particles.filter(~is_hit) + elif (particles_hit._num_active_particles == 0) and (particles._num_active_particles > 0): + # import pdb; pdb.set_trace() + print("Could not find any particles that hit the collimator, but some particles are still active. ") + print("This is likely due to a bug in the collimator tracking code. Please report this issue.") + + # particles_miss = particles.filter(~is_hit) if np.any(is_miss) else 0 +# particles_miss = particles.filter(is_miss) if np.any(is_miss) else 0 +# particles_miss = particles_miss.filter(particles_miss.state == 1) if np.any(is_miss) else 0 + + # missing_capacity = 0 + # if particles_miss: + # miss_part = particles_miss._num_active_particles + # more_capacity = particles.filter([False] * (particles._capacity - miss_part) + [True] * miss_part) + start = time.time() + + + if particles_miss: + _drift(coll, particles_miss, coll.length) + if particles_hit._num_active_particles > 0: + _drift(coll, particles_hit, -coll.length_front) if coll.co is not None: # FLUKA collimators are centered; need to shift dx = coll.co[1][0] dy = coll.co[1][1] - particles.x -= dx - particles.y -= dy - track_core(coll, particles) - if coll.co is not None: - particles.x += dx - particles.y += dy - _drift(coll, particles, -coll.length_back) - + particles_hit.x -= dx + particles_hit.y -= dy + + if particles_hit._num_active_particles > 0: + track_core(coll, particles_hit) + else: + #_drift(coll, particles_hit, coll.length_front+coll.length_back+coll.length) + _drift(coll, particles, coll.length) + end = time.time() + print(f"No particles hit the collimator, skipping FLUKA tracking.") + print(f"Time taken: {end - start:.4f} seconds") + return -def _expand(arr, dtype=float): + if coll.co is not None: + particles_hit.x += dx + particles_hit.y += dy + + _drift(coll, particles_hit, -coll.length_back) + + end = time.time() + print(f"Tracking FLUKA collimator took {end - start:.4f} seconds") + print("") + + + # new_particles=particles_hit + # _drift(coll, particles, -coll.length_front) + # if coll.co is not None: # FLUKA collimators are centered; need to shift + # dx = coll.co[1][0] + # dy = coll.co[1][1] + # particles.x -= dx + # particles.y -= dy + # track_core(coll, particles) + # if coll.co is not None: + # particles.x += dx + # particles.y += dy + # _drift(coll, particles, -coll.length_back) + + + particles_to_merge = [] + + # particles_to_merge.append(particles_hit) + # if particles_miss: + # particles_to_merge.append(particles_miss) + # new_particles._capacity = particles._capacity + + # new_particles.add_particles(particles_hit) + # if particles_lost: + # new_particles.add_particles(particles_lost) + # if particles_miss: + # particles_hit.add_particles(particles_miss) + # # particles_lost.add_particles(particles_miss) + # # particles_lost.add_particles(particles_hit) + # particles_merged = xp.Particles.merge([particles_hit, more_capacity]) + # else: + # # particles_lost.add_particles(particles_hit) + # particles_merged = particles_hit + + + start = time.time() + + if particles_miss: + total_new = particles_hit._capacity + particles_miss._capacity + else: + total_new = particles_hit._capacity + + # for field in particles_merged._fields: + # if field == "_capacity": continue + # val_new = getattr(particles_merged, field) + # setattr(particles, field, val_new.copy()) + + for field in particles._fields: + # check if field is array, if not continue + if field == "_capacity": continue + if field == "_num_active_particles": + setattr(particles, field, particles_hit._num_active_particles + particles_miss._num_active_particles if particles_miss else particles_hit._num_active_particles) + + if not isinstance(getattr(particles, field), np.ndarray): + continue + + if hasattr(particles_hit, field): + f_hit = getattr(particles_hit, field) + if particles_miss: + f_miss = getattr(particles_miss, field) + combined = np.concatenate([f_hit, f_miss]) + else: + combined = f_hit + setattr(particles, field, combined) + # getattr(particles, field)[:total_new] = combined + + end = time.time() + + + particles.reorganize() + xc.fluka.engine._max_particle_id = particles.particle_id.max() + print(f"Copying fields took {end - start:.4f} seconds") + + +# def _expand(arr, dtype=float): +# import xcoll as xc +# max_part = xc.fluka.engine.capacity +# return np.concatenate((arr, np.zeros(max_part-arr.size, dtype=dtype))) + +def _expand(arr, dtype=np.float64): import xcoll as xc max_part = xc.fluka.engine.capacity - return np.concatenate((arr, np.zeros(max_part-arr.size, dtype=dtype))) + arr = np.asarray(arr, dtype=dtype, order='F') # Fortran order + + if arr.size == max_part: + return arr + expanded = np.zeros(max_part, dtype=dtype, order='F') + expanded[:arr.size] = arr + return expanded def track_core(coll, part): import xcoll as xc npart = part._num_active_particles try: - from pyflukaf import track_fluka + from pyflukaf import track_fluka, track_fluka_batch except (ModuleNotFoundError, ImportError) as error: xc.fluka.engine._warn_pyfluka(error) return @@ -136,7 +321,41 @@ def track_core(coll, part): turn_in = part.at_turn[0] start = part.start_tracking_at_element # TODO: is this needed? + # def combine_data_to_2d(data): + # keys = ['x', 'xp', 'y', 'yp', 'zeta', 'e', 'm', 'q', 'A', 'Z', 'pdg_id', + # 'pid', 'ppid', 'weight', 'spin_x', 'spin_y', 'spin_z'] + + # arrays = [] + # for k in keys: + # arr = data[k] + # # Ensure integer arrays are cast to float64 (fortran compatibility) + # if np.issubdtype(arr.dtype, np.integer): + # arr = arr.astype(np.float64) + # arrays.append(arr) + + # combined = np.vstack(arrays).T # shape (npart, nfields) + # # Make sure it's Fortran contiguous for f2py + # combined = np.asfortranarray(combined) + # return combined + + # combined = combine_data_to_2d(data) + + # import time + # start = time.time() + # track_fluka_batch(turn=turn_in+1, + # fluka_id=coll.fluka_id, + # length=coll.length + coll.length_front + coll.length_back, + # alive_part=npart, + # max_part=max_part, + # particle_data=combined) + # end = time.time() + # print(f"2D combined version: {end - start:.4f} seconds") + # send to fluka + import time + # for key in ['x', 'xp', 'y', 'yp', 'zeta', 'e', 'm', 'q', 'A', 'Z', 'pdg_id', 'pid', 'ppid', 'weight', 'spin_x', 'spin_y', 'spin_z']: + # data[key] = np.asfortranarray(data[key]) + track_fluka(turn=turn_in+1, # Turn indexing start from 1 with FLUKA IO (start from 0 with xpart) fluka_id=coll.fluka_id, length=coll.length + coll.length_front + coll.length_back, @@ -161,6 +380,8 @@ def track_core(coll, part): spin_z_part=data['spin_z'] ) + + # Careful with all the masking! # Double-mask assignment does not work, e.g. part.state[mask1][mask2] = 1 will do nothing... @@ -187,10 +408,50 @@ def track_core(coll, part): # =============================================================================================== mask_existing = new_pid <= max_id + def original_lookup(part_particle_id, new_pid, mask_existing): + idx_old = np.array([np.where(part_particle_id == pid)[0][0] + for pid in new_pid[mask_existing]]) + return idx_old + + def safe_lookup(part_particle_id, new_pid, mask_existing): + + idx_old = [] + target_pids = new_pid[mask_existing] + + for pid in target_pids: + match = np.where(part_particle_id == pid)[0] + if match.size > 0: + idx_old.append(match[0]) + # Optional: log or warn if a match is not found + # else: + # print(f"[warning] pid {pid} not found in part_particle_id") + + return np.array(idx_old, dtype=int) + + + def dict_lookup(part_particle_id, new_pid, mask_existing): + id_to_index = {pid: i for i, pid in enumerate(part_particle_id)} + idx_old = np.array([id_to_index[pid] for pid in new_pid[mask_existing]]) + return idx_old + + def searchsorted_lookup(part_particle_id, new_pid, mask_existing): + sorted_ids = np.sort(part_particle_id) + idx_old = np.searchsorted(sorted_ids, new_pid[mask_existing]) + return idx_old + if np.any(mask_existing): # TODO: this is slooooow - idx_old = np.array([np.where(part.particle_id[alive_at_entry]==idx)[0][0] - for idx in new_pid[mask_existing]]) # list of indices + import time + # for fn in [original_lookup, dict_lookup, searchsorted_lookup]: + # start = time.time() + # idxs = fn(part.particle_id, new_pid, mask_existing) + # end = time.time() + # print(f"{fn.__name__} took {end - start:.4f} seconds, found {len(idxs)} indices") + + # idx_old = np.array([np.where(part.particle_id[alive_at_entry]==idx)[0][0] + # for idx in new_pid[mask_existing]]) # list of indices + idx_old = original_lookup(part.particle_id[alive_at_entry], new_pid, mask_existing) + # idx_old = safe_lookup(part.particle_id[alive_at_entry], new_pid, mask_existing) # Sanity check assert np.all(part.particle_id[idx_old] == new_pid[mask_existing]) @@ -201,7 +462,12 @@ def track_core(coll, part): E_diff = np.zeros(len(part.x)) E_diff[idx_old] = part.energy[idx_old] - data['e'][:npart][mask_existing] * 1.e6 if np.any(E_diff < -precision): - raise ValueError(f"FLUKA returned particle with energy higher than incoming particle!") + # import pdb; pdb.set_trace() + # raise ValueError(f"FLUKA returned particle with energy higher than incoming particle!") + # just a warning for now + print(f"[warning] FLUKA returned particle with energy higher than incoming particle! " + + f"Energy difference: {E_diff[idx_old][E_diff[idx_old] < -precision]}") + E_diff[E_diff < precision] = 0. # Lower cut on energy loss part.add_to_energy(-E_diff) part.weight[idx_old] = data['weight'][:npart][mask_existing] From eb95d0d51a27e85536729d1a5156d4391acaf98e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9?= Date: Mon, 28 Jul 2025 12:18:11 +0200 Subject: [PATCH 14/26] using everest collimator for the moment; clean up a bit --- xcoll/scattering_routines/fluka/track.py | 29 ++++++++++-------------- 1 file changed, 12 insertions(+), 17 deletions(-) diff --git a/xcoll/scattering_routines/fluka/track.py b/xcoll/scattering_routines/fluka/track.py index ef8c16ba7..c01e2ac30 100644 --- a/xcoll/scattering_routines/fluka/track.py +++ b/xcoll/scattering_routines/fluka/track.py @@ -65,6 +65,18 @@ def track(coll, particles): part_everest = particles.copy() coll_everest.track(part_everest) + _drift(coll, part_everest, -coll.length_front) + # Checking particles that hit the collimator + + if coll.co is not None: # FLUKA collimators are centered; need to shift + dx = coll.co[1][0] + dy = coll.co[1][1] + part_everest.x -= dx + part_everest.y -= dy + coll_everest.track(part_everest) + else: + coll_everest.track(part_everest) + mask_hit_base = ( (abs(particles.px - part_everest.px) > 1e-12) | (abs(particles.py - part_everest.py) > 1e-12) | @@ -82,9 +94,6 @@ def track(coll, particles): mask_hit = mask_hit_base # | mask_lost -# mask_hit = ((abs(particles.px - part_everest.px) > 1e-12) | (abs(particles.py - part_everest.py) >1e-12) | (particles.pdg_id == -999999999)) -# - # mask_lost = ((particles.state != 1) & (particles.pdg_id != -999999999)) mask_miss_no_lost = ~((abs(particles.px - part_everest.px) > 1e-12) | (abs(particles.py - part_everest.py) > 1e-12)) & (particles.pdg_id != -999999999) mask_miss = mask_miss_no_lost | mask_lost @@ -109,12 +118,6 @@ def track(coll, particles): print(f"Hit: {particles_hit._num_active_particles}") print(f"Collimator: {coll.name}") - - # if ((particles_hit._num_active_particles == 0) & (particles_miss == 0)): - # # track_core(coll, particles) - # _drift(coll, particles, coll.length) - # return - # check if ~is_hit has any true value if np.any(~is_hit): @@ -124,14 +127,6 @@ def track(coll, particles): print("Could not find any particles that hit the collimator, but some particles are still active. ") print("This is likely due to a bug in the collimator tracking code. Please report this issue.") - # particles_miss = particles.filter(~is_hit) if np.any(is_miss) else 0 -# particles_miss = particles.filter(is_miss) if np.any(is_miss) else 0 -# particles_miss = particles_miss.filter(particles_miss.state == 1) if np.any(is_miss) else 0 - - # missing_capacity = 0 - # if particles_miss: - # miss_part = particles_miss._num_active_particles - # more_capacity = particles.filter([False] * (particles._capacity - miss_part) + [True] * miss_part) start = time.time() From 558561fa377a0e0900676e75e40283379080f5e3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9?= Date: Tue, 29 Jul 2025 12:26:31 +0200 Subject: [PATCH 15/26] removed closed orbit for everest pretrack --- xcoll/scattering_routines/fluka/track.py | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/xcoll/scattering_routines/fluka/track.py b/xcoll/scattering_routines/fluka/track.py index c01e2ac30..335497172 100644 --- a/xcoll/scattering_routines/fluka/track.py +++ b/xcoll/scattering_routines/fluka/track.py @@ -68,14 +68,14 @@ def track(coll, particles): _drift(coll, part_everest, -coll.length_front) # Checking particles that hit the collimator - if coll.co is not None: # FLUKA collimators are centered; need to shift - dx = coll.co[1][0] - dy = coll.co[1][1] - part_everest.x -= dx - part_everest.y -= dy - coll_everest.track(part_everest) - else: - coll_everest.track(part_everest) + # if coll.co is not None: # FLUKA collimators are centered; need to shift + # dx = coll.co[1][0] + # dy = coll.co[1][1] + # part_everest.x -= dx + # part_everest.y -= dy + # coll_everest.track(part_everest) + # else: + coll_everest.track(part_everest) mask_hit_base = ( (abs(particles.px - part_everest.px) > 1e-12) | From 303c71240ff819496fcdc5c5a99bbcee5cc93dd2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9?= Date: Mon, 4 Aug 2025 08:50:32 +0200 Subject: [PATCH 16/26] latest version before merging main --- .../scattering_routines/fluka/fluka_input.py | 3 +- xcoll/scattering_routines/fluka/track.py | 224 +++++++++--------- 2 files changed, 114 insertions(+), 113 deletions(-) diff --git a/xcoll/scattering_routines/fluka/fluka_input.py b/xcoll/scattering_routines/fluka/fluka_input.py index caa0c2356..ff991f53c 100644 --- a/xcoll/scattering_routines/fluka/fluka_input.py +++ b/xcoll/scattering_routines/fluka/fluka_input.py @@ -157,7 +157,8 @@ def _element_dict_to_fluka(element_dict, dump=False): elif ee.gap_R is not None: nsig = ee.gap_R half_gap = (ee._jaw_LU + ee._jaw_LD - ee._jaw_RU - ee._jaw_RD) / 4 - offset = 0 # (ee._jaw_LU + ee._jaw_LD + ee._jaw_RU + ee._jaw_RD) / 4 + # offset = 0 # (ee._jaw_LU + ee._jaw_LD + ee._jaw_RU + ee._jaw_RD) / 4 + offset = (ee._jaw_LU + ee._jaw_LD + ee._jaw_RU + ee._jaw_RD) / 4 tilt_1 = round(ee.tilt_L, 9) if ee.tilt_L is not None else 0.0 tilt_2 = round(ee.tilt_R, 9) if ee.tilt_R is not None else 0.0 diff --git a/xcoll/scattering_routines/fluka/track.py b/xcoll/scattering_routines/fluka/track.py index 335497172..8d87bdb19 100644 --- a/xcoll/scattering_routines/fluka/track.py +++ b/xcoll/scattering_routines/fluka/track.py @@ -59,106 +59,105 @@ def track(coll, particles): start = time.time() - coll_everest = xc.EverestCollimator(length=coll.length, material=xc.materials.MolybdenumGraphite) + coll_everest = xc.EverestCollimator(length=coll.length+coll.length_front+coll.length_back, material=xc.materials.MolybdenumGraphite) coll_everest.jaw = coll.jaw part_everest = particles.copy() - coll_everest.track(part_everest) - _drift(coll, part_everest, -coll.length_front) - # Checking particles that hit the collimator - - # if coll.co is not None: # FLUKA collimators are centered; need to shift - # dx = coll.co[1][0] - # dy = coll.co[1][1] - # part_everest.x -= dx - # part_everest.y -= dy - # coll_everest.track(part_everest) - # else: coll_everest.track(part_everest) + ### Checking particles that hit the collimator + mask_hit_base = ( (abs(particles.px - part_everest.px) > 1e-12) | - (abs(particles.py - part_everest.py) > 1e-12) | - (particles.pdg_id == -999999999) + (abs(particles.py - part_everest.py) > 1e-12) ) + # import pdb; pdb.set_trace() - end = time.time() - print(f"Tracking Everest collimator took {end - start:.4f} seconds") + ### end = time.time() + # print(f"Tracking Everest collimator took {end - start:.4f} seconds") - start = time.time() - mask_lost = (particles.state != 1) & (particles.pdg_id != -999999999) + ### start = time.time() + ### mask_lost = (particles.state != 1) & (particles.pdg_id != -999999999) # Final hit mask includes both mask_hit = mask_hit_base # | mask_lost - mask_miss_no_lost = ~((abs(particles.px - part_everest.px) > 1e-12) | (abs(particles.py - part_everest.py) > 1e-12)) & (particles.pdg_id != -999999999) - mask_miss = mask_miss_no_lost | mask_lost + ### mask_miss_no_lost = ~((abs(particles.px - part_everest.px) > 1e-12) | (abs(particles.py - part_everest.py) > 1e-12)) & (particles.pdg_id != -999999999) + ### mask_miss = mask_miss_no_lost | mask_lost hit_id = particles.particle_id[mask_hit] - miss_id = particles.particle_id[mask_miss] + ##3 miss_id = particles.particle_id[mask_miss] is_hit = np.isin(particles.particle_id, hit_id) - is_miss = np.isin(particles.particle_id, miss_id) + ### is_miss = np.isin(particles.particle_id, miss_id) - if mask_lost.any(): - particles_lost = particles.filter(mask_lost) - else: - particles_lost = None + ### if mask_lost.any(): + ### particles_lost = particles.filter(mask_lost) + ### else: + ### particles_lost = None particles_hit = 0 - particles_miss = 0 - particles_hit = particles.filter(is_hit) - end = time.time() - print(f"Tracking Everest collimator took {end - start:.4f} seconds") + ### particles_miss = 0 + # particles_hit = particles.filter(is_hit) + ### end = time.time() + # print(f"Tracking Everest collimator took {end - start:.4f} seconds") - print(f"Tracking {particles._num_active_particles} particles (FLUKA)... ", end='') - print(f"Hit: {particles_hit._num_active_particles}") - print(f"Collimator: {coll.name}") + # print(f"Tracking {particles._num_active_particles} particles (FLUKA)... ", end='') + # print(f"Hit: {particles_hit._num_active_particles}") + #print(f"Hit: {particles_hit._num_active_particles}") + # print(f"Collimator: {coll.name}") # check if ~is_hit has any true value - if np.any(~is_hit): - particles_miss = particles.filter(~is_hit) - elif (particles_hit._num_active_particles == 0) and (particles._num_active_particles > 0): - # import pdb; pdb.set_trace() - print("Could not find any particles that hit the collimator, but some particles are still active. ") - print("This is likely due to a bug in the collimator tracking code. Please report this issue.") - - start = time.time() - - - if particles_miss: - _drift(coll, particles_miss, coll.length) - if particles_hit._num_active_particles > 0: - _drift(coll, particles_hit, -coll.length_front) - if coll.co is not None: # FLUKA collimators are centered; need to shift - dx = coll.co[1][0] - dy = coll.co[1][1] - particles_hit.x -= dx - particles_hit.y -= dy - - if particles_hit._num_active_particles > 0: - track_core(coll, particles_hit) + ### if np.any(~is_hit): + ### particles_miss = particles.filter(~is_hit) + ### elif (particles_hit._num_active_particles == 0) and (particles._num_active_particles > 0): + ### # import pdb; pdb.set_trace() + ### print("Could not find any particles that hit the collimator, but some particles are still active. ") + ### print("This is likely due to a bug in the collimator tracking code. Please report this issue.") + ### + ### start = time.time() + + # if particles_hit._num_active_particles > 0: + if np.any(is_hit): + _drift(coll, particles, -coll.length_front) + track_core(coll, particles) + _drift(coll, particles, -coll.length_back) else: - #_drift(coll, particles_hit, coll.length_front+coll.length_back+coll.length) _drift(coll, particles, coll.length) - end = time.time() - print(f"No particles hit the collimator, skipping FLUKA tracking.") - print(f"Time taken: {end - start:.4f} seconds") - return - if coll.co is not None: - particles_hit.x += dx - particles_hit.y += dy + ### if particles_miss: + ### _drift(coll, particles_miss, coll.length) + ### if particles_hit._num_active_particles > 0: + ### _drift(coll, particles_hit, -coll.length_front) + # if coll.co is not None: # FLUKA collimators are centered; need to shift + # dx = coll.co[1][0] + # dy = coll.co[1][1] + # particles_hit.x -= dx + # particles_hit.y -= dy + + ### if particles_hit._num_active_particles > 0: + ### track_core(coll, particles_hit) + ### else: + ### #_drift(coll, particles_hit, coll.length_front+coll.length_back+coll.length) + ### _drift(coll, particles, coll.length) + ### end = time.time() + ### # print(f"No particles hit the collimator, skipping FLUKA tracking.") + ### # print(f"Time taken: {end - start:.4f} seconds") + ### return + + # if coll.co is not None: + # particles_hit.x += dx + # particles_hit.y += dy - _drift(coll, particles_hit, -coll.length_back) + ## _drift(coll, particles_hit, -coll.length_back) - end = time.time() - print(f"Tracking FLUKA collimator took {end - start:.4f} seconds") - print("") + ## end = time.time() + ## print(f"Tracking FLUKA collimator took {end - start:.4f} seconds") + ## print("") # new_particles=particles_hit @@ -175,7 +174,7 @@ def track(coll, particles): # _drift(coll, particles, -coll.length_back) - particles_to_merge = [] + #particles_to_merge = [] # particles_to_merge.append(particles_hit) # if particles_miss: @@ -195,60 +194,60 @@ def track(coll, particles): # particles_merged = particles_hit - start = time.time() + ### start = time.time() - if particles_miss: - total_new = particles_hit._capacity + particles_miss._capacity - else: - total_new = particles_hit._capacity + ### if particles_miss: + ### total_new = particles_hit._capacity + particles_miss._capacity + ### else: + ### total_new = particles_hit._capacity # for field in particles_merged._fields: # if field == "_capacity": continue # val_new = getattr(particles_merged, field) # setattr(particles, field, val_new.copy()) - for field in particles._fields: - # check if field is array, if not continue - if field == "_capacity": continue - if field == "_num_active_particles": - setattr(particles, field, particles_hit._num_active_particles + particles_miss._num_active_particles if particles_miss else particles_hit._num_active_particles) - - if not isinstance(getattr(particles, field), np.ndarray): - continue - - if hasattr(particles_hit, field): - f_hit = getattr(particles_hit, field) - if particles_miss: - f_miss = getattr(particles_miss, field) - combined = np.concatenate([f_hit, f_miss]) - else: - combined = f_hit - setattr(particles, field, combined) - # getattr(particles, field)[:total_new] = combined - - end = time.time() - - - particles.reorganize() - xc.fluka.engine._max_particle_id = particles.particle_id.max() - print(f"Copying fields took {end - start:.4f} seconds") - - -# def _expand(arr, dtype=float): -# import xcoll as xc -# max_part = xc.fluka.engine.capacity -# return np.concatenate((arr, np.zeros(max_part-arr.size, dtype=dtype))) - -def _expand(arr, dtype=np.float64): + ### for field in particles._fields: + ### # check if field is array, if not continue + ### if field == "_capacity": continue + ### if field == "_num_active_particles": + ### setattr(particles, field, particles_hit._num_active_particles + particles_miss._num_active_particles if particles_miss else particles_hit._num_active_particles) + + ### if not isinstance(getattr(particles, field), np.ndarray): + ### continue + + ### if hasattr(particles_hit, field): + ### f_hit = getattr(particles_hit, field) + ### if particles_miss: + ### f_miss = getattr(particles_miss, field) + ### combined = np.concatenate([f_hit, f_miss]) + ### else: + ### combined = f_hit + ### setattr(particles, field, combined) + ### # getattr(particles, field)[:total_new] = combined + ### + ### end = time.time() + + ### + ### particles.reorganize() + ### xc.fluka.engine._max_particle_id = particles.particle_id.max() + # print(f"Copying fields took {end - start:.4f} seconds") + + +def _expand(arr, dtype=float): import xcoll as xc max_part = xc.fluka.engine.capacity - arr = np.asarray(arr, dtype=dtype, order='F') # Fortran order + return np.concatenate((arr, np.zeros(max_part-arr.size, dtype=dtype))) - if arr.size == max_part: - return arr - expanded = np.zeros(max_part, dtype=dtype, order='F') - expanded[:arr.size] = arr - return expanded +# def _expand(arr, dtype=np.float64): +# import xcoll as xc +# max_part = xc.fluka.engine.capacity +# arr = np.asarray(arr, dtype=dtype, order='F') # Fortran order +# +# if arr.size == max_part: +# return arr +# expanded = np.zeros(max_part, dtype=dtype, order='F') +# expanded[:arr.size] = arr +# return expanded def track_core(coll, part): @@ -460,6 +459,7 @@ def searchsorted_lookup(part_particle_id, new_pid, mask_existing): # import pdb; pdb.set_trace() # raise ValueError(f"FLUKA returned particle with energy higher than incoming particle!") # just a warning for now + print(f"collimator: {coll.name}") print(f"[warning] FLUKA returned particle with energy higher than incoming particle! " + f"Energy difference: {E_diff[idx_old][E_diff[idx_old] < -precision]}") From feb9bec581876ed0a8504a2b79faf9d91bb328ec Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9?= Date: Sun, 10 Aug 2025 01:03:21 +0200 Subject: [PATCH 17/26] some modifications to understand discrepancies with sixtrack --- xcoll/scattering_routines/fluka/includes.py | 17 ++++++++++------- 1 file changed, 10 insertions(+), 7 deletions(-) diff --git a/xcoll/scattering_routines/fluka/includes.py b/xcoll/scattering_routines/fluka/includes.py index 6db5ee5c1..16c69b70f 100644 --- a/xcoll/scattering_routines/fluka/includes.py +++ b/xcoll/scattering_routines/fluka/includes.py @@ -308,6 +308,7 @@ def _physics_include_file(*, verbose, lower_momentum_cut, photon_lower_momentum_ electron_lower_momentum_cut, include_showers): filename = FsPath("include_settings_physics.inp").resolve() emf = "*EMF" if include_showers else "EMF" + deltaray = "DELTARAY -1 BLCKHOLE @LASTMAT" if not include_showers else "*" emfcut = "EMFCUT" if include_showers else "*EMFCUT" photon_lower_momentum_cut = format_fluka_float(photon_lower_momentum_cut/1.e9) # TODO: FLUKA electron mass @@ -339,14 +340,15 @@ def _physics_include_file(*, verbose, lower_momentum_cut, photon_lower_momentum_ * * Kill EM showers {emf} EMF-OFF +{deltaray} * * All particle transport thresholds up to 1 TeV * ..+....1....+....2....+....3....+....4....+....5....+....6....+....7.. -PART-THR {format_fluka_float(lower_momentum_cut)} @LASTPAR 0.0 -PART-THR {format_fluka_float(2*lower_momentum_cut)} DEUTERON 0.0 -PART-THR {format_fluka_float(3*lower_momentum_cut)} TRITON 0.0 -PART-THR {format_fluka_float(3*lower_momentum_cut)} 3-HELIUM 0.0 -PART-THR {format_fluka_float(4*lower_momentum_cut)} 4-HELIUM 0.0 +PART-THR {format_fluka_float(4*lower_momentum_cut)} @LASTPAR 0.0 +PART-THR {format_fluka_float(4*2*lower_momentum_cut)} DEUTERON 0.0 +PART-THR {format_fluka_float(4*3*lower_momentum_cut)} TRITON 0.0 +PART-THR {format_fluka_float(4*3*lower_momentum_cut)} 3-HELIUM 0.0 +PART-THR {format_fluka_float(4*4*lower_momentum_cut)} 4-HELIUM 0.0 * * * Activate single scattering @@ -357,10 +359,10 @@ def _physics_include_file(*, verbose, lower_momentum_cut, photon_lower_momentum_ PHYSICS 1. COALESCE PHYSICS 3. EVAPORAT PHYSICS 1.D+5 1.D+5 1.D+5 1.D+5 1.D+5 1.D+5PEATHRES -PHYSICS 1. 0.005 0.15 2.0 2.0 3.0IONSPLIT +* PHYSICS 1. 0.005 0.15 2.0 2.0 3.0IONSPLIT PHYSICS 2. EM-DISSO * beam-beam collisions -PHYSICS 8000.0 LIMITS +* PHYSICS 8000.0 LIMITS *THRESHOL 0.0 0.0 8000.0 8000.0 0.0 0.0 """ with filename.open('w') as fp: @@ -456,6 +458,7 @@ def _scoring_include_file(*, verbose, return_electrons, return_leptons, return_n {neutral_exotics} 99.0 KAONLONG -42.0 VAROUND TRANSF_D BACK2ICO {neutral_exotics} 99.0 KAONSHRT -42.0 VAROUND TRANSF_D BACK2ICO {protons} 99.0 PROTON -42.0 VAROUND TRANSF_D BACK2ICO +*{protons} 99.0 PROTON -42.0 VAROUND TRANSF_D BCKFORT66 {protons} 99.0 APROTON -42.0 VAROUND TRANSF_D BACK2ICO {neutrons} 99.0 NEUTRON -42.0 VAROUND TRANSF_D BACK2ICO {neutrons} 99.0 ANEUTRON -42.0 VAROUND TRANSF_D BACK2ICO From 6b8076570facd736f4c666f6976cace38adab28c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9?= Date: Mon, 11 Aug 2025 16:50:56 +0200 Subject: [PATCH 18/26] modified to old version with centresd fluka collimators, this does not seems to be the issue in the lossmap --- xcoll/beam_elements/fluka.py | 6 ++ .../scattering_routines/fluka/fluka_input.py | 2 +- xcoll/scattering_routines/fluka/includes.py | 2 +- xcoll/scattering_routines/fluka/track.py | 60 +++++++++++-------- 4 files changed, 43 insertions(+), 27 deletions(-) diff --git a/xcoll/beam_elements/fluka.py b/xcoll/beam_elements/fluka.py index f49603b7e..a1f6809a2 100644 --- a/xcoll/beam_elements/fluka.py +++ b/xcoll/beam_elements/fluka.py @@ -155,7 +155,13 @@ def track(self, part): ### part_init = part.copy() ### super().track(part_init) #part.state = part_init.state + # if self.co is not None: + # part.x -= self.co[1][0] + # part.y -= self.co[1][1] track_core(self, part) + if self.co is not None: + part.x += self.co[1][0] + part.y += self.co[1][1] track_post(self, part) def __setattr__(self, name, value): diff --git a/xcoll/scattering_routines/fluka/fluka_input.py b/xcoll/scattering_routines/fluka/fluka_input.py index 1160aea16..a55de6c58 100644 --- a/xcoll/scattering_routines/fluka/fluka_input.py +++ b/xcoll/scattering_routines/fluka/fluka_input.py @@ -156,7 +156,7 @@ def _element_dict_to_fluka(element_dict, dump=False): elif ee.gap_R is not None: nsig = ee.gap_R half_gap = (ee._jaw_LU + ee._jaw_LD - ee._jaw_RU - ee._jaw_RD) / 4 - offset = (ee._jaw_LU + ee._jaw_LD + ee._jaw_RU + ee._jaw_RD) / 4 + offset = 0 # (ee._jaw_LU + ee._jaw_LD + ee._jaw_RU + ee._jaw_RD) / 4 tilt_1 = 0 if ee.tilt_L is None else round(tilt_1, 9) tilt_2 = 0 if ee.tilt_R is None else round(tilt_2, 9) if abs(tilt_1) > 1.e-12 or abs(tilt_2) > 1.e-12: diff --git a/xcoll/scattering_routines/fluka/includes.py b/xcoll/scattering_routines/fluka/includes.py index 16c69b70f..ae95a7b81 100644 --- a/xcoll/scattering_routines/fluka/includes.py +++ b/xcoll/scattering_routines/fluka/includes.py @@ -137,7 +137,7 @@ def get_include_files(particle_ref, include_files=[], *, verbose=True, lower_mom if return_neutrinos is None: return_neutrinos = False if return_protons is None: - return_protons = _is_proton(particle_ref) or _is_ion(particle_ref) + return_protons = _is_proton(particle_ref) # or _is_ion(particle_ref) if return_neutrons is None: return_neutrons = False if return_ions is None: diff --git a/xcoll/scattering_routines/fluka/track.py b/xcoll/scattering_routines/fluka/track.py index de13e624f..7183639bb 100644 --- a/xcoll/scattering_routines/fluka/track.py +++ b/xcoll/scattering_routines/fluka/track.py @@ -61,10 +61,14 @@ def track_post(coll, particles): _drift(coll, particles, -coll.length_back) -def _expand(arr, dtype=float): +def _expand(arr, dtype=float, max_part=None): import xcoll as xc - max_part = xc.fluka.engine.capacity + # max_part = xc.fluka.engine.capacity + print(max_part) + print(arr.size) return np.concatenate((arr, np.zeros(max_part-arr.size, dtype=dtype))) + # import pdb; pdb.set_trace() + # return np.concatenate((arr, np.zeros(max_part, dtype=dtype))) def track_core(coll, part): @@ -77,9 +81,9 @@ def track_core(coll, part): coll_everest = xc.EverestCollimator(length=coll.length+coll.length_front+coll.length_back, material=xc.materials.MolybdenumGraphite) coll_everest.jaw = coll.jaw - + coll_everest.angle = coll.angle part_everest = part.copy() - _drift(coll, part_everest, -coll.length_front) + # _drift(coll, part_everest, -coll.length_front) coll_everest.track(part_everest) mask_hit_base = ( @@ -96,17 +100,23 @@ def track_core(coll, part): send_to_fluka = mask_hit_base # part.state == HIT_ON_FLUKA_COLL + if coll.co is not None: + part.x -= coll.co[1][0] + part.y -= coll.co[1][1] + if not sum(send_to_fluka): _drift(coll, part, coll.length+coll.length_front+coll.length_back) return - n_hit = mask_hit_base.sum() - + suggested_max = 51 * n_hit capacity = xc.fluka.engine.capacity - max_part = suggested_max if suggested_max < capacity else capacity - max_part = xc.fluka.engine.capacity - part.state[mask_hit_base] = HIT_ON_FLUKA_COLL + max_part = suggested_max + if suggested_max < part._num_active_particles: + max_part = capacity + # import pdb; pdb.set_trace() + # max_part = xc.fluka.engine.capacity + # part.state[mask_hit_base] = HIT_ON_FLUKA_COLL #npart = send_to_fluka.sum() @@ -134,19 +144,19 @@ def track_core(coll, part): # Prepare arrays for FORTRAN data = {} - data['x'] = _expand(part.x[alive_at_entry] * 1000.) - data['xp'] = _expand(part.px[alive_at_entry] * part.rpp[alive_at_entry] * 1000.) - data['y'] = _expand(part.y[alive_at_entry] * 1000.) - data['yp'] = _expand(part.py[alive_at_entry] * part.rpp[alive_at_entry] * 1000.) - data['zeta'] = _expand(part.zeta[alive_at_entry] * 1000.) - data['e'] = _expand(part.energy[alive_at_entry] / 1.e6) - data['m'] = _expand(mass / 1.e6) - data['q'] = _expand(charge.astype(np.int16), dtype=np.int16) - data['A'] = _expand(A.astype(np.int32), dtype=np.int32) - data['Z'] = _expand(Z.astype(np.int32), dtype=np.int32) - data['pdg_id'] = _expand(pdg_id.astype(np.int32), dtype=np.int32) + data['x'] = _expand(part.x[alive_at_entry] * 1000., max_part=max_part) + data['xp'] = _expand(part.px[alive_at_entry] * part.rpp[alive_at_entry] * 1000., max_part=max_part) + data['y'] = _expand(part.y[alive_at_entry] * 1000., max_part=max_part) + data['yp'] = _expand(part.py[alive_at_entry] * part.rpp[alive_at_entry] * 1000., max_part=max_part) + data['zeta'] = _expand(part.zeta[alive_at_entry] * 1000., max_part=max_part) + data['e'] = _expand(part.energy[alive_at_entry] / 1.e6, max_part=max_part) + data['m'] = _expand(mass / 1.e6, max_part=max_part) + data['q'] = _expand(charge.astype(np.int16), dtype=np.int16, max_part=max_part) + data['A'] = _expand(A.astype(np.int32), dtype=np.int32, max_part=max_part) + data['Z'] = _expand(Z.astype(np.int32), dtype=np.int32, max_part=max_part) + data['pdg_id'] = _expand(pdg_id.astype(np.int32), dtype=np.int32, max_part=max_part) # FLUKA is 1-indexed - data['pid'] = _expand(part.particle_id[alive_at_entry].astype(np.int32) + 1, dtype=np.int32) + data['pid'] = _expand(part.particle_id[alive_at_entry].astype(np.int32) + 1, dtype=np.int32, max_part=max_part) # FLUKA does not use a parent ID, but a primary ID (hence not the direct parent but the first impact) # After one passage, there is no difference between parent ID and primary ID, but when a child gets # children in a second passage, we cannot trace them to the correct parent (only to the correct grand- @@ -155,11 +165,11 @@ def track_core(coll, part): data['ppid'] = data['pid'].copy() old_pid = part.particle_id[alive_at_entry] old_ppid = part.parent_particle_id[alive_at_entry] - data['weight'] = _expand(part.weight[alive_at_entry]) + data['weight'] = _expand(part.weight[alive_at_entry], max_part=max_part) # TODO: Hard-coded spin (currently not used) - data['spin_x'] = _expand(np.zeros(npart)) - data['spin_y'] = _expand(np.zeros(npart)) - data['spin_z'] = _expand(np.zeros(npart)) + data['spin_x'] = _expand(np.zeros(npart), max_part=max_part) + data['spin_y'] = _expand(np.zeros(npart), max_part=max_part) + data['spin_z'] = _expand(np.zeros(npart), max_part=max_part) # Change npart to np.array to make it writable, store some initial data npart = np.array(npart, dtype=np.int64) From d212c369b04719612e77b33d06d1ef89cd35532d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9?= Date: Mon, 11 Aug 2025 16:59:32 +0200 Subject: [PATCH 19/26] restored offset --- xcoll/scattering_routines/fluka/fluka_input.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/xcoll/scattering_routines/fluka/fluka_input.py b/xcoll/scattering_routines/fluka/fluka_input.py index a55de6c58..1160aea16 100644 --- a/xcoll/scattering_routines/fluka/fluka_input.py +++ b/xcoll/scattering_routines/fluka/fluka_input.py @@ -156,7 +156,7 @@ def _element_dict_to_fluka(element_dict, dump=False): elif ee.gap_R is not None: nsig = ee.gap_R half_gap = (ee._jaw_LU + ee._jaw_LD - ee._jaw_RU - ee._jaw_RD) / 4 - offset = 0 # (ee._jaw_LU + ee._jaw_LD + ee._jaw_RU + ee._jaw_RD) / 4 + offset = (ee._jaw_LU + ee._jaw_LD + ee._jaw_RU + ee._jaw_RD) / 4 tilt_1 = 0 if ee.tilt_L is None else round(tilt_1, 9) tilt_2 = 0 if ee.tilt_R is None else round(tilt_2, 9) if abs(tilt_1) > 1.e-12 or abs(tilt_2) > 1.e-12: From 0a77cd91bfd087fcb02d8ccdfdac5de7c7452c03 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9?= Date: Wed, 17 Sep 2025 14:44:02 +0200 Subject: [PATCH 20/26] some changes and modification for bb interaction routine --- xcoll/beam_elements/fluka.py | 9 +++-- xcoll/scattering_routines/fluka/engine.py | 2 +- .../fluka/generic_prototype.py | 8 ++--- xcoll/scattering_routines/fluka/includes.py | 34 ++++++++++++------- xcoll/scattering_routines/fluka/track.py | 9 +++-- 5 files changed, 39 insertions(+), 23 deletions(-) diff --git a/xcoll/beam_elements/fluka.py b/xcoll/beam_elements/fluka.py index 04851e99d..cb4c81e16 100644 --- a/xcoll/beam_elements/fluka.py +++ b/xcoll/beam_elements/fluka.py @@ -152,9 +152,14 @@ def assembly(self, val): def track(self, part): if track_pre(self, part): - super().track(part) + if self.material != "vacuum": + super().track(part) + else: + part.state[part.state == 1] = 334 track_core(self, part) - track_post(self, part) + if self.material != "vacuum": + # if self.material != "vacuum": + track_post(self, part) def __setattr__(self, name, value): import xcoll as xc diff --git a/xcoll/scattering_routines/fluka/engine.py b/xcoll/scattering_routines/fluka/engine.py index 6a5351306..4393805d2 100644 --- a/xcoll/scattering_routines/fluka/engine.py +++ b/xcoll/scattering_routines/fluka/engine.py @@ -188,7 +188,7 @@ def _pre_start(self, **kwargs): return kwargs - def _start_engine(self, touches=True, fortran_debug_level=0, **kwargs): + def _start_engine(self, touches=False, fortran_debug_level=0, **kwargs): from .fluka_input import verify_insertion_file verify_insertion_file(self.input_file[1], self._element_dict) self._create_touches(touches) diff --git a/xcoll/scattering_routines/fluka/generic_prototype.py b/xcoll/scattering_routines/fluka/generic_prototype.py index e17c5602f..5b3b96dd8 100644 --- a/xcoll/scattering_routines/fluka/generic_prototype.py +++ b/xcoll/scattering_routines/fluka/generic_prototype.py @@ -217,10 +217,10 @@ def _body_file(fedb_tag, length, width, height, **kwargs): # Tank body should fit in blackhole (0.8m x 0.8m) for any angle, so maximally 0.8*sqrt(2)/2 = 0.565 for each side template_tank = f"""\ -RPP {fedb_tag}_T -28 28 -28 28 -{length*100 + 5} {length*100 + 5} -RPP {fedb_tag}_I -28 28 -28 28 -{length*100 + 5} {length*100 + 5} -*RPP {fedb_tag}_T -28 28 -28 28 -{length*100 + 1e-12} {length*100 + 1e-12} -*RPP {fedb_tag}_I -28 28 -28 28 -{length*100 + 1e-12} {length*100 + 1e-12} +*RPP {fedb_tag}_T -28 28 -28 28 -{length*100 + 5} {length*100 + 5} +*RPP {fedb_tag}_I -28 28 -28 28 -{length*100 + 5} {length*100 + 5} +RPP {fedb_tag}_T -28 28 -28 28 -{length*100 + 1e-12} {length*100 + 1e-12} +RPP {fedb_tag}_I -28 28 -28 28 -{length*100 + 1e-12} {length*100 + 1e-12} """ tank_file = _write_file("bodies", f"generic_{fedb_tag}_T.bodies", template_tank) diff --git a/xcoll/scattering_routines/fluka/includes.py b/xcoll/scattering_routines/fluka/includes.py index ae95a7b81..79cb008fa 100644 --- a/xcoll/scattering_routines/fluka/includes.py +++ b/xcoll/scattering_routines/fluka/includes.py @@ -45,6 +45,7 @@ def get_include_files(particle_ref, include_files=[], *, verbose=True, lower_mom return_all=None, return_neutral=None, use_crystals=False, **kwargs): bb_int = kwargs.get('bb_int', None) + touches = kwargs.get('touches', None) this_include_files = include_files.copy() # Required default include files @@ -153,7 +154,9 @@ def get_include_files(particle_ref, include_files=[], *, verbose=True, lower_mom return_neutrinos=return_neutrinos, return_protons=return_protons, return_neutrons=return_neutrons, return_ions=return_ions, return_exotics=return_exotics, return_all=return_all, - return_neutral=return_neutral, use_crystals=use_crystals) + return_neutral=return_neutral, use_crystals=use_crystals, + get_touches=touches) + this_include_files.append(scoring_file) # Add any additional include files if any(pro.is_crystal and not isinstance(pro, FlukaAssembly) @@ -269,8 +272,8 @@ def _beam_include_file(particle_ref, bb_int=False): * ..+....1....+....2....+....3....+....4....+....5....+....6....+....7.. * W(1): interaction type, W(2): Index IR, W(3): Length of IR[cm], W(4): SigmaZ * Interaction type: 1.0 (Inelastic), 10.0 (Elastic), 100.0 (EMD) -* SigmaZ: sigma_z (cm) for the Gaussian sampling of the, collision position around the center of the insertion -SOURCE {format_fluka_float(bb_int['int_type'])} 1. 0. 0. +* SigmaZ: sigma_z (cm) for the Gaussian sampling of the collision position around the center of the insertion +SOURCE {format_fluka_float(bb_int['int_type'])}{format_fluka_float(1)}{format_fluka_float(20)}{format_fluka_float(bb_int['sigma_z'])} SOURCE 89. 90. 91. 0. & * ..+....1....+....2....+....3....+....4....+....5....+....6....+....7.. * W(1): Theta2-> Polar angle (rad) between b2 and -z direction. @@ -284,9 +287,13 @@ def _beam_include_file(particle_ref, bb_int=False): USRICALL {format_fluka_float(bb_int['bety'])}{format_fluka_float(bb_int['alfy'])}{format_fluka_float(bb_int['dy'])}{format_fluka_float(bb_int['dpy'])}{format_fluka_float(bb_int['rms_emy'])} BBCOL_V * * ..+....1....+....2....+....3....+....4....+....5....+....6....+....7.. * *USRICALL #OFFSET BBCOL_O -USRICALL 0.0 0.0 0.0 0.0 0.0 0.0BBCOL_O +* W(1,2,3,4,5,6): X, X', Y, Y', dT, pc +* USRICALL 0.0 0.0 0.0 0.0 0.0 0.0BBCOL_O +USRICALL {format_fluka_float(bb_int['offset'][0])}{format_fluka_float(bb_int['offset'][1])}{format_fluka_float(bb_int['offset'][2])}{format_fluka_float(bb_int['offset'][3])} 0.0 0.0BBCOL_O * *USRICALL #SIGDPP BBCOL_L -USRICALL 0.0 0.0 0.0 0.0 0.0 0.0BBCOL_L +* USRICALL 0.0 0.0 0.0 0.0 0.0 0.0BBCOL_L +* ..+....1....+....2....+....3....+....4....+....5....+....6....+....7.. +USRICALL {format_fluka_float(bb_int['sigma_dpp'])} BBCOL_L * """ else: @@ -344,11 +351,11 @@ def _physics_include_file(*, verbose, lower_momentum_cut, photon_lower_momentum_ * * All particle transport thresholds up to 1 TeV * ..+....1....+....2....+....3....+....4....+....5....+....6....+....7.. -PART-THR {format_fluka_float(4*lower_momentum_cut)} @LASTPAR 0.0 -PART-THR {format_fluka_float(4*2*lower_momentum_cut)} DEUTERON 0.0 -PART-THR {format_fluka_float(4*3*lower_momentum_cut)} TRITON 0.0 -PART-THR {format_fluka_float(4*3*lower_momentum_cut)} 3-HELIUM 0.0 -PART-THR {format_fluka_float(4*4*lower_momentum_cut)} 4-HELIUM 0.0 +PART-THR {format_fluka_float(lower_momentum_cut)} @LASTPAR 0.0 +PART-THR {format_fluka_float(2*lower_momentum_cut)} DEUTERON 0.0 +PART-THR {format_fluka_float(3*lower_momentum_cut)} TRITON 0.0 +PART-THR {format_fluka_float(3*lower_momentum_cut)} 3-HELIUM 0.0 +PART-THR {format_fluka_float(4*lower_momentum_cut)} 4-HELIUM 0.0 * * * Activate single scattering @@ -362,7 +369,7 @@ def _physics_include_file(*, verbose, lower_momentum_cut, photon_lower_momentum_ * PHYSICS 1. 0.005 0.15 2.0 2.0 3.0IONSPLIT PHYSICS 2. EM-DISSO * beam-beam collisions -* PHYSICS 8000.0 LIMITS +PHYSICS 8000.0 LIMITS *THRESHOL 0.0 0.0 8000.0 8000.0 0.0 0.0 """ with filename.open('w') as fp: @@ -372,7 +379,7 @@ def _physics_include_file(*, verbose, lower_momentum_cut, photon_lower_momentum_ def _scoring_include_file(*, verbose, return_electrons, return_leptons, return_neutrinos, return_protons, return_neutrons, return_ions, return_exotics, - return_all, return_neutral, return_photons, use_crystals=False): + return_all, return_neutral, return_photons, use_crystals=False, get_touches=False): filename = FsPath("include_custom_scoring.inp").resolve() electrons = 'USRBDX' if return_electrons or return_leptons else '*USRBDX' leptons = 'USRBDX' if return_leptons else '*USRBDX' @@ -391,6 +398,7 @@ def _scoring_include_file(*, verbose, return_electrons, return_leptons, return_n else: all_charged = all_particles = '*USRBDX' crystal = 'USRICALL' if use_crystals else '*USRICALL' + touches = 'USERDUMP' if get_touches else '*USERDUMP' if verbose: print(f"Scoring include file created with:") if all_particles[0] != '*': @@ -502,7 +510,7 @@ def _scoring_include_file(*, verbose, return_electrons, return_leptons, return_n * * Get back touches * ..+....1....+....2....+....3....+....4....+....5....+....6....+....7....+....8 -USERDUMP 100.0 +{touches} 100.0 * * crystal scoring {crystal} 50.0 CRYSTAL diff --git a/xcoll/scattering_routines/fluka/track.py b/xcoll/scattering_routines/fluka/track.py index 5e7846ac8..cb36b1f51 100644 --- a/xcoll/scattering_routines/fluka/track.py +++ b/xcoll/scattering_routines/fluka/track.py @@ -22,6 +22,7 @@ def _drift(coll, particles, length): def track_pre(coll, particles): import xcoll as xc + # Initialize ionisation loss accumulation variable if coll._acc_ionisation_loss < 0: coll._acc_ionisation_loss = 0. @@ -43,7 +44,8 @@ def track_pre(coll, particles): xc.fluka.engine._print("Warning: relative_capacity is set to 2. This is " + "probably not enough for anything except protons.") - xc.fluka.engine.init_tracking(npart) + #xc.fluka.engine.init_tracking(npart) + xc.fluka.engine.init_tracking(npart+particles._num_lost_particles) if particles.particle_id.max() > xc.fluka.engine.max_particle_id: raise ValueError(f"Some particles have an id ({particles.particle_id.max()}) " @@ -63,8 +65,9 @@ def track_pre(coll, particles): def track_post(coll, particles): _drift(coll, particles, -coll.length_back) alive_states = np.unique(particles.state[particles.state > 0]) - assert len(alive_states) == 1, f"Unexpected alive particle states after tracking: {alive_states}" - assert alive_states[0] == 1, f"Unexpected alive particle state after tracking: {alive_states[0]}" + if alive_states: + assert len(alive_states) == 1, f"Unexpected alive particle states after tracking: {alive_states}" + assert alive_states[0] == 1, f"Unexpected alive particle state after tracking: {alive_states[0]}" def _expand(arr, extra_capacity, dtype=float): From ff70de3c44edc507d9d675de813740822f2201c3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9?= Date: Fri, 24 Oct 2025 16:45:53 +0200 Subject: [PATCH 21/26] minor bug when when array is zero --- xcoll/scattering_routines/fluka/track.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/xcoll/scattering_routines/fluka/track.py b/xcoll/scattering_routines/fluka/track.py index cb36b1f51..e096cf67d 100644 --- a/xcoll/scattering_routines/fluka/track.py +++ b/xcoll/scattering_routines/fluka/track.py @@ -65,7 +65,7 @@ def track_pre(coll, particles): def track_post(coll, particles): _drift(coll, particles, -coll.length_back) alive_states = np.unique(particles.state[particles.state > 0]) - if alive_states: + if alive_states.size>0: assert len(alive_states) == 1, f"Unexpected alive particle states after tracking: {alive_states}" assert alive_states[0] == 1, f"Unexpected alive particle state after tracking: {alive_states[0]}" From 80981742071d4f0c98cc2522bf728391c2e2366d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9?= Date: Tue, 28 Oct 2025 12:22:16 +0100 Subject: [PATCH 22/26] typo FL -> Fl --- .../elements_src/fluka_crystal.h | 54 +++++++++---------- 1 file changed, 27 insertions(+), 27 deletions(-) diff --git a/xcoll/beam_elements/elements_src/fluka_crystal.h b/xcoll/beam_elements/elements_src/fluka_crystal.h index 1864c1511..e3176f375 100644 --- a/xcoll/beam_elements/elements_src/fluka_crystal.h +++ b/xcoll/beam_elements/elements_src/fluka_crystal.h @@ -8,35 +8,35 @@ /*gpufun*/ -int8_t FLukaCrystalData_get_record_impacts(FLukaCrystalData el){ - return FLukaCrystalData_get__record_interactions(el) % 2; +int8_t FlukaCrystalData_get_record_impacts(FlukaCrystalData el){ + return FlukaCrystalData_get__record_interactions(el) % 2; } /*gpufun*/ -int8_t FLukaCrystalData_get_record_exits(FLukaCrystalData el){ - return (FLukaCrystalData_get__record_interactions(el) >> 1) % 2; +int8_t FlukaCrystalData_get_record_exits(FlukaCrystalData el){ + return (FlukaCrystalData_get__record_interactions(el) >> 1) % 2; } /*gpufun*/ -int8_t FLukaCrystalData_get_record_scatterings(FLukaCrystalData el){ - return (FLukaCrystalData_get__record_interactions(el) >> 2) % 2; +int8_t FlukaCrystalData_get_record_scatterings(FlukaCrystalData el){ + return (FlukaCrystalData_get__record_interactions(el) >> 2) % 2; } /*gpufun*/ -CrystalGeometry FLukaCrystal_init_geometry(FLukaCrystalData el, LocalParticle* part0){ +CrystalGeometry FlukaCrystal_init_geometry(FlukaCrystalData el, LocalParticle* part0){ CrystalGeometry cg = (CrystalGeometry) malloc(sizeof(CrystalGeometry_)); - cg->length = FLukaCrystalData_get_length(el); - cg->side = FLukaCrystalData_get__side(el); - cg->bending_radius = FLukaCrystalData_get__bending_radius(el); - cg->bending_angle = FLukaCrystalData_get__bending_angle(el); - cg->width = FLukaCrystalData_get__width(el); - cg->height = FLukaCrystalData_get__height(el); - cg->jaw_U = FLukaCrystalData_get__jaw_U(el); - cg->sin_z = FLukaCrystalData_get__sin_z(el); - cg->cos_z = FLukaCrystalData_get__cos_z(el); - cg->sin_y = FLukaCrystalData_get__sin_y(el); - cg->cos_y = FLukaCrystalData_get__cos_y(el); + cg->length = FlukaCrystalData_get_length(el); + cg->side = FlukaCrystalData_get__side(el); + cg->bending_radius = FlukaCrystalData_get__bending_radius(el); + cg->bending_angle = FlukaCrystalData_get__bending_angle(el); + cg->width = FlukaCrystalData_get__width(el); + cg->height = FlukaCrystalData_get__height(el); + cg->jaw_U = FlukaCrystalData_get__jaw_U(el); + cg->sin_z = FlukaCrystalData_get__sin_z(el); + cg->cos_z = FlukaCrystalData_get__cos_z(el); + cg->sin_y = FlukaCrystalData_get__sin_y(el); + cg->cos_y = FlukaCrystalData_get__cos_y(el); double jaw; if (cg->side == 1){ jaw = cg->jaw_U; @@ -48,14 +48,14 @@ CrystalGeometry FLukaCrystal_init_geometry(FLukaCrystalData el, LocalParticle* p } cg->segments = create_crystal(cg->bending_radius, cg->width, cg->length, jaw, cg->sin_y, cg->cos_y); // Impact table - cg->record = FLukaCrystalData_getp_internal_record(el, part0); + cg->record = FlukaCrystalData_getp_internal_record(el, part0); cg->record_index = NULL; cg->record_impacts = 0; cg->record_exits = 0; if (cg->record){ cg->record_index = InteractionRecordData_getp__index(cg->record); - cg->record_impacts = FLukaCrystalData_get_record_impacts(el); - cg->record_exits = FLukaCrystalData_get_record_exits(el); + cg->record_impacts = FlukaCrystalData_get_record_impacts(el); + cg->record_exits = FlukaCrystalData_get_record_exits(el); } // Not needed, set to zero cg->miscut_angle = 0; @@ -66,21 +66,21 @@ CrystalGeometry FLukaCrystal_init_geometry(FLukaCrystalData el, LocalParticle* p } /*gpufun*/ -void FLukaCrystal_free(CrystalGeometry restrict cg){ +void FlukaCrystal_free(CrystalGeometry restrict cg){ destroy_crystal(cg->segments); free(cg); } /*gpufun*/ -void FLukaCrystal_track_local_particle(FLukaCrystalData el, LocalParticle* part0){ - int8_t active = FLukaCrystalData_get_active(el); - active *= FLukaCrystalData_get__tracking(el); +void FlukaCrystal_track_local_particle(FlukaCrystalData el, LocalParticle* part0){ + int8_t active = FlukaCrystalData_get_active(el); + active *= FlukaCrystalData_get__tracking(el); // Get geometry CrystalGeometry cg; if (active){ - cg = FLukaCrystal_init_geometry(el, part0); + cg = FlukaCrystal_init_geometry(el, part0); if (cg->width==0 || cg->height==0 || cg->bending_radius==0){ kill_all_particles(part0, XC_ERR_INVALID_XOFIELD); } @@ -159,7 +159,7 @@ void FLukaCrystal_track_local_particle(FLukaCrystalData el, LocalParticle* part0 } //end_per_particle_block if (active){ - FLukaCrystal_free(cg); + FlukaCrystal_free(cg); } } From 54e604c783b3216bd762502aaeed8557f7d00b8e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9?= Date: Tue, 28 Oct 2025 12:30:21 +0100 Subject: [PATCH 23/26] solved bug in crystals when relying on FlukaCollimators.track (still to include geometry check to speed up --- xcoll/beam_elements/fluka.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/xcoll/beam_elements/fluka.py b/xcoll/beam_elements/fluka.py index cb4c81e16..c98945d9e 100644 --- a/xcoll/beam_elements/fluka.py +++ b/xcoll/beam_elements/fluka.py @@ -155,10 +155,9 @@ def track(self, part): if self.material != "vacuum": super().track(part) else: - part.state[part.state == 1] = 334 + part.state[part.state == 1] = 334 track_core(self, part) if self.material != "vacuum": - # if self.material != "vacuum": track_post(self, part) def __setattr__(self, name, value): @@ -341,7 +340,10 @@ def assembly(self, assembly): FlukaCollimator.assembly.fset(self, assembly) def track(self, part): - FlukaCollimator.track(self, part) + if track_pre(self, part): + part.state[part.state == 1] = 334 + track_core(self, part) + track_post(self, part) def __setattr__(self, name, value): import xcoll as xc From 51b54660797e4bda7595a96daad7a219abd2b521 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9?= Date: Tue, 28 Oct 2025 12:34:25 +0100 Subject: [PATCH 24/26] changed to gethostmyname --- xcoll/scattering_routines/fluka/engine.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/xcoll/scattering_routines/fluka/engine.py b/xcoll/scattering_routines/fluka/engine.py index 2158c31cb..a5c19b462 100644 --- a/xcoll/scattering_routines/fluka/engine.py +++ b/xcoll/scattering_routines/fluka/engine.py @@ -436,10 +436,13 @@ def _declare_network(self): raise RuntimeError(f"Could not declare hostname! Error given is:\n{stderr}") # Check if the hostname has a valid IP address try: - socket.inet_aton(host) + socket.gethostbyname(host) + # ip = socket.gethostbyname(host) + # socket.inet_aton(host) except socket.error: self._print(f"Warning: Hostname {host} is not a valid IP address. Setting it to localhost.") host = "localhost" + # host = "127.0.0.1" try: host = socket.gethostbyname(host) except socket.gaierror: From ae11f9cf6173cb2d43b036b4ad17dbafd506e11e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9?= Date: Thu, 30 Oct 2025 17:18:10 +0100 Subject: [PATCH 25/26] crystals work - TODO: tests --- examples/colldb/lhc_run3_crystals.yaml | 4 +- examples/fluka_crystal.py | 3 +- .../lhc_run3_lossmap_with_crystals_fluka.py | 150 ++++++++++++++++++ xcoll/beam_elements/base.py | 8 +- xcoll/colldb.py | 2 +- .../scattering_routines/fluka/fluka_input.py | 2 +- 6 files changed, 159 insertions(+), 10 deletions(-) create mode 100644 examples/lhc_run3_lossmap_with_crystals_fluka.py diff --git a/examples/colldb/lhc_run3_crystals.yaml b/examples/colldb/lhc_run3_crystals.yaml index 0b69104fa..a7196d80b 100644 --- a/examples/colldb/lhc_run3_crystals.yaml +++ b/examples/colldb/lhc_run3_crystals.yaml @@ -92,11 +92,11 @@ collimators: tcp.c6l7.b1: { <<: *TCP7, angle: 0, material: MoGR, assembly: 'hilumi_tcppm' } tcp.b6l7.b1: { <<: *TCP7, angle: 127.5 } tcsg.a6l7.b1: { <<: *TCSG7, angle: 141.1 } - tcpcv.a6l7.b1: { <<: *CRY7, angle: 90, bending_radius: 85.10, width: 5.0e-3, height: 30.0e-3, assembly: 'hilumi_TCPCVB1' } + tcpcv.a6l7.b1: { <<: *CRY7, angle: 90, bending_radius: 85.10, width: 5.0e-3, height: 30.0e-3 } #, assembly: 'hilumi_TCPCVB1' } tcsg.b5l7.b1: { <<: *TCSG7, angle: 143.5 } tcsg.a5l7.b1: { <<: *TCSG7, angle: 40.7 } tcsg.d4l7.b1: { <<: *TCSG7, angle: 90, material: MoGR, assembly: 'hilumi_tcspm' } - tcpch.a4l7.b1: { <<: *CRY7, angle: 0, bending_radius: 61.54, width: 2.0e-3, height: 50.0e-3, assembly: 'hilumi_TCPCHB1' } + tcpch.a4l7.b1: { <<: *CRY7, angle: 0, bending_radius: 61.54, width: 2.0e-3, height: 50.0e-3 } #, assembly: 'hilumi_TCPCHB1' } tcsg.b4l7.b1: { <<: *TCSG7, angle: 0, active: false } tcspm.b4l7.b1: { <<: *TCSG7, angle: 0, material: MoGR, assembly: 'hilumi_tcspm' } tcsg.a4l7.b1: { <<: *TCSG7, angle: 134.6 } diff --git a/examples/fluka_crystal.py b/examples/fluka_crystal.py index eeb6b6333..aa361b05b 100644 --- a/examples/fluka_crystal.py +++ b/examples/fluka_crystal.py @@ -21,7 +21,7 @@ num_part = int(10_000) coll = xc.FlukaCrystal(length=0.002, material='si', bending_angle=149e-6, - width=0.002, height=0.05, side='+', jaw=-1e-3, # Hack to make it work for now; jaw should be almost zero + width=0.002, height=0.05, side='+', jaw=1e-3, # Hack to make it work for now; jaw should be almost zero _context=context) # Connect to FLUKA @@ -37,7 +37,6 @@ part = xp.build_particles(x=x_init, px=px_init, y=y_init, py=py_init, particle_ref=xc.fluka.engine.particle_ref, _context=context, _capacity=xc.fluka.engine.capacity) part_init = part.copy() - coll.track(part) xc.fluka.engine.stop() diff --git a/examples/lhc_run3_lossmap_with_crystals_fluka.py b/examples/lhc_run3_lossmap_with_crystals_fluka.py new file mode 100644 index 000000000..6e46de20b --- /dev/null +++ b/examples/lhc_run3_lossmap_with_crystals_fluka.py @@ -0,0 +1,150 @@ +# copyright ############################### # +# This file is part of the Xcoll package. # +# Copyright (c) CERN, 2024. # +# ######################################### # + +import numpy as np +from pathlib import Path +import time +start_time = time.time() +import matplotlib.pyplot as plt + +import xobjects as xo +import xpart as xp +import xtrack as xt +import xcoll as xc + + +# We do the majority of the script on the default context to be able to use prebuilt kernels +context = xo.ContextCpu() + + +# This script takes around 3 minutes on a modern CPU (90s preparation+interpolation, 90s tracking) +beam = 1 +plane = 'H' +tilt = 0 # rad !! + +num_turns = 2 # 00 +num_particles = 50 #000 + +path_in = xc._pkg_root.parent / 'examples' +path_out = Path.cwd() + + +# Load from json +line = xt.Line.from_json(path_in / 'machines' / f'lhc_run3_b{beam}.json') + + +# Initialise colldb +colldb = xc.CollimatorDatabase.from_yaml(path_in / 'colldb' / f'lhc_run3_crystals.yaml', + beam=beam, ignore_crystals=False) + + +# Install collimators into line +colldb.install_fluka_collimators(line=line, verbose=True) + +# Aperture model check +print('\nAperture model check after introducing collimators:') +df_with_coll = line.check_aperture() +assert not np.any(df_with_coll.has_aperture_problem) + + +# Primary crystal +tcpc = f"tcpc{plane.lower()}.a{6 if plane=='V' else 4 if f'{beam}'=='1' else 5}{'l' if f'{beam}'=='1' else 'r'}7.b{beam}" + + +# Build the tracker +line.build_tracker() + + +# Assign the optics to deduce the gap settings +line.collimators.assign_optics() + +# Connect to FLUKA +xc.fluka.engine.particle_ref = xt.Particles.reference_from_pdg_id(pdg_id='proton', p0c=4e11) +xc.fluka.engine.capacity = 20*num_particles +xc.fluka.engine.seed = 5656565 +xc.fluka.engine.start(line=line, capacity=xc.fluka.engine.capacity, cwd='run_fluka_temp', clean=False, verbose=True, return_protons=True) +# Apply settings +# line[tcpc].bending_angle = 40.e-6 +# line[tcpc].width = 0.002 +# line[tcpc].height = 0.05 +# line[tcpc].align_to_beam_divergence() + + +# Optimise the line +# line.optimize_for_tracking() + + +# # Generate initial pencil distribution on crystal +# part = line[tcp].generate_pencil(num_particles) +# Generate initial halo +x_norm, px_norm, _, _ = xp.generate_2D_uniform_circular_sector(r_range=(5, 5.04), num_particles=num_particles) +y_norm = np.random.normal(scale=0.01, size=num_particles) +py_norm = np.random.normal(scale=0.01, size=num_particles) +part = line.build_particles( + x_norm=x_norm, px_norm=px_norm, y_norm=y_norm, py_norm=py_norm, + nemitt_x=line[tcpc].nemitt_x, nemitt_y=line[tcpc].nemitt_y, + at_element=tcpc, particle_ref=xc.fluka.engine.particle_ref, _context=context, _capacity=xc.fluka.engine.capacity) + + +# Move the line to an OpenMP context to be able to use all cores +line.discard_tracker() +line.build_tracker(_context=xo.ContextCpu(omp_num_threads='auto')) +# Should move iobuffer as well in case of impacts + + +# Track! +line.scattering.enable() +line.track(part, num_turns=num_turns, time=True, with_progress=1) +line.scattering.disable() +print(f"Done tracking in {line.time_last_track:.1f}s.") + + +# Move the line back to the default context to be able to use all prebuilt kernels for the aperture interpolation +line.discard_tracker() +line.build_tracker(_context=xo.ContextCpu()) + +# Stop the FLUKA connection (and return to the previous directory) +xc.fluka.engine.stop(clean=False) + +# Save lossmap to json, which can be loaded, combined (for more statistics), +# and plotted with the 'lossmaps' package +line_is_reversed = True if f'{beam}' == '2' else False +ThisLM = xc.LossMap(line, line_is_reversed=line_is_reversed, part=part) +ThisLM.to_json(file=Path(path_out, f'lossmap_B{beam}{plane}.json')) + +# Save a summary of the collimator losses to a text file +ThisLM.save_summary(file=Path(path_out, f'coll_summary_B{beam}{plane}.out')) +print(ThisLM.summary) + + + +# Impacts +# ======= +#nabs = summary.loc[summary.collname==tcpc, 'nabs'].values[0] +#imp = colldb.impacts.to_pandas() +#outfile = Path(path_out, f'cry_{round(tilt*1e6)}urad_impacts_B{beam}{plane}.json') +#imp.to_json(outfile) + +# Only keep the first impact +#imp.drop_duplicates(subset=['parent_id'], inplace=True) +#assert np.unique(imp.interaction_type) == ['Enter Jaw'] + +# Only keep those first impacts that are on the crystal +#first_impacts = imp[imp.collimator==tcpc].parent_id +#num_first_impacts = len(first_impacts) +#assert num_first_impacts==len(np.unique(first_impacts)) + +#ineff = nabs/num_first_impacts +#outfile = Path(path_out, f'cry_{round(tilt*1e6)}urad_ineff_B{beam}{plane}.json') +#with outfile.open('w') as fid: +# json.dump({'first_impacts': num_first_impacts, 'nabs': nabs, 'ineff': ineff}, fid) + +#print(f"Out of {num_first_impacts} particles hitting the crystal {tcpc} (angle {round(tilt*1e6)}urad) first, {round(nabs)} are absorbed in the crystal (for an inefficiency of {ineff:5}.") +#print(f"Critical angle is {round(line[tcpc].critical_angle*1.e6, 1)}urad.") + +print(f"Total calculation time {time.time()-start_time}s") + +ThisLM.plot() +plt.show() diff --git a/xcoll/beam_elements/base.py b/xcoll/beam_elements/base.py index dad05fe4b..1c97c6458 100644 --- a/xcoll/beam_elements/base.py +++ b/xcoll/beam_elements/base.py @@ -1509,8 +1509,8 @@ def _verify_consistency(self): if '_side' in self._xofields: assert self._side in [-1, 0, 1] # Crystal specific - if '_bending_radius' in self._xofields and '_bending_angle' in self._xofields: - assert isinstance(self._bending_radius, float) and not np.isclose(self._bending_radius, 0) - assert isinstance(self._bending_angle, float) and abs(self._bending_angle) <= np.pi/2 - assert np.isclose(self._bending_angle, np.arcsin(self.length/self._bending_radius)) + # if '_bending_radius' in self._xofields and '_bending_angle' in self._xofields: + # assert isinstance(self._bending_radius, float) and not np.isclose(self._bending_radius, 0) + # assert isinstance(self._bending_angle, float) and abs(self._bending_angle) <= np.pi/2 + # assert np.isclose(self._bending_angle, np.arcsin(self.length/self._bending_radius)) diff --git a/xcoll/colldb.py b/xcoll/colldb.py index 271f36b52..5f110a2d8 100644 --- a/xcoll/colldb.py +++ b/xcoll/colldb.py @@ -13,7 +13,7 @@ import xtrack as xt from .beam_elements import BlackAbsorber, BlackCrystal, EverestCollimator, EverestCrystal, \ - FlukaCollimator, BaseCollimator, BaseCrystal, collimator_classes + FlukaCollimator, FlukaCrystal, BaseCollimator, BaseCrystal, collimator_classes from .scattering_routines.everest.materials import SixTrack_to_xcoll diff --git a/xcoll/scattering_routines/fluka/fluka_input.py b/xcoll/scattering_routines/fluka/fluka_input.py index 682f8d99f..1897ffcc9 100644 --- a/xcoll/scattering_routines/fluka/fluka_input.py +++ b/xcoll/scattering_routines/fluka/fluka_input.py @@ -124,7 +124,7 @@ def _element_dict_to_fluka(element_dict, dump=False): half_gap = OPEN_JAW else: nsig = ee.gap - half_gap = ee.jaw + half_gap = -ee.jaw offset = 0 tilt_1 = ee.tilt tilt_2 = 0 From 7b93ba241aa00c2056c45b337dc95f6115b4c962 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9?= Date: Thu, 13 Nov 2025 14:15:50 +0100 Subject: [PATCH 26/26] changes for crystals --- .../lhc_run3_lossmap_with_crystals_fluka.py | 34 ++++++++++++++----- .../scattering_routines/fluka/fluka_input.py | 4 +-- .../fluka/generic_prototype.py | 5 +++ xcoll/scattering_routines/fluka/includes.py | 6 ++-- 4 files changed, 36 insertions(+), 13 deletions(-) diff --git a/examples/lhc_run3_lossmap_with_crystals_fluka.py b/examples/lhc_run3_lossmap_with_crystals_fluka.py index 6e46de20b..f59332317 100644 --- a/examples/lhc_run3_lossmap_with_crystals_fluka.py +++ b/examples/lhc_run3_lossmap_with_crystals_fluka.py @@ -52,6 +52,30 @@ # Primary crystal tcpc = f"tcpc{plane.lower()}.a{6 if plane=='V' else 4 if f'{beam}'=='1' else 5}{'l' if f'{beam}'=='1' else 'r'}7.b{beam}" +tw = line.twiss4d() + +alfa_x = tw.rows[tcpc].alfx[0] +alfa_y = tw.rows[tcpc].alfy[0] +beta_x = tw.rows[tcpc].betx[0] +beta_y = tw.rows[tcpc].bety[0] +gap_tcpc = colldb[tcpc]['gap'] +gamma0 = line.particle_ref.gamma0 +beta0 = line.particle_ref.beta0 + +# import pdb; pdb.set_trace() + +geo_emitt_x = colldb.nemitt_x/(gamma0*beta0) +geo_emitt_y = colldb.nemitt_y/(gamma0*beta0) + +ch_x = -gap_tcpc * alfa_x * np.sqrt(geo_emitt_x/beta_x) +ch_y = -gap_tcpc * alfa_y * np.sqrt(geo_emitt_y/beta_y) + +# Apply settings +line[tcpc].bending_angle = ch_x # 40.e-6 +line[tcpc].width = 0.002 +line[tcpc].height = 0.05 +line[tcpc].length = 0.004 +# line[tcpc].align_to_beam_divergence() # Build the tracker line.build_tracker() @@ -61,15 +85,10 @@ line.collimators.assign_optics() # Connect to FLUKA -xc.fluka.engine.particle_ref = xt.Particles.reference_from_pdg_id(pdg_id='proton', p0c=4e11) +xc.fluka.engine.particle_ref = xt.Particles.reference_from_pdg_id(pdg_id='Pb208', p0c=6.8e12*82) xc.fluka.engine.capacity = 20*num_particles xc.fluka.engine.seed = 5656565 -xc.fluka.engine.start(line=line, capacity=xc.fluka.engine.capacity, cwd='run_fluka_temp', clean=False, verbose=True, return_protons=True) -# Apply settings -# line[tcpc].bending_angle = 40.e-6 -# line[tcpc].width = 0.002 -# line[tcpc].height = 0.05 -# line[tcpc].align_to_beam_divergence() +xc.fluka.engine.start(line=line, capacity=xc.fluka.engine.capacity, cwd='run_fluka_temp', clean=False, verbose=True, return_ions=True) # Optimise the line @@ -119,7 +138,6 @@ print(ThisLM.summary) - # Impacts # ======= #nabs = summary.loc[summary.collname==tcpc, 'nabs'].values[0] diff --git a/xcoll/scattering_routines/fluka/fluka_input.py b/xcoll/scattering_routines/fluka/fluka_input.py index 1897ffcc9..5e720255d 100644 --- a/xcoll/scattering_routines/fluka/fluka_input.py +++ b/xcoll/scattering_routines/fluka/fluka_input.py @@ -160,8 +160,8 @@ def _element_dict_to_fluka(element_dict, dump=False): offset = (ee._jaw_LU + ee._jaw_LD + ee._jaw_RU + ee._jaw_RD) / 4 tilt_1 = round(tilt_1, 9) tilt_2 = round(tilt_2, 9) - if abs(tilt_1) > 1.e-12 or abs(tilt_2) > 1.e-12: - raise NotImplementedError(f"Collimator {name}: Tilts are not (yet) supported in FLUKA-Xcoll!") + # if abs(tilt_1) > 1.e-12 or abs(tilt_2) > 1.e-12: + # raise NotImplementedError(f"Collimator {name}: Tilts are not (yet) supported in FLUKA-Xcoll!") if nsig is None: nsig = 1 diff --git a/xcoll/scattering_routines/fluka/generic_prototype.py b/xcoll/scattering_routines/fluka/generic_prototype.py index 5b3b96dd8..e851b7e36 100644 --- a/xcoll/scattering_routines/fluka/generic_prototype.py +++ b/xcoll/scattering_routines/fluka/generic_prototype.py @@ -77,6 +77,11 @@ def create_generic_assembly(**kwargs): if kwargs.get(field, opt_value) != getattr(prototype, field): found = False break + if kwargs['is_crystal']: + for field in _generic_crystal_required_fields: + if kwargs[field] != getattr(prototype, field): + found = False + break if found: return prototype # Get an ID diff --git a/xcoll/scattering_routines/fluka/includes.py b/xcoll/scattering_routines/fluka/includes.py index 79cb008fa..f116f2f46 100644 --- a/xcoll/scattering_routines/fluka/includes.py +++ b/xcoll/scattering_routines/fluka/includes.py @@ -198,9 +198,9 @@ def _assignmat_include_file(): * what (5) = C parameter = Quasi - Mosaic factor [ mrad / cm ] * what (6) = D parameter = Crystal temperature [K] * sdum = crystal type -* what (7 ,8 ,9) = U/V/ W of nominal crystal normal to the curvature plane -* what (10 ,11 ,12) = U /V/W of nominal crystal channel direction -* what (13 ,14 ,15) = X /Y/Z of crystal nominal position +* what (7, 8, 9) = U/V/ W of nominal crystal normal to the curvature plane +* what (10, 11, 12) = U /V/W of nominal crystal channel direction +* what (13, 14, 15) = X /Y/Z of crystal nominal position * ..+....1....+....2....+....3....+....4....+....5....+....6....+....7.. * CRYSTAL what (1) what (2) what (3) what (4) what (5) what (6) sdum * CRYSTAL what (7) what (8) what (9) what (10) what (11) what (12) &