diff --git a/examples/colldb/lhc_run3_crystals.yaml b/examples/colldb/lhc_run3_crystals.yaml index 406d3945..ef7a292a 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 e2d5b3ce..aa361b05 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,10 +34,9 @@ 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() - 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 00000000..f5933231 --- /dev/null +++ b/examples/lhc_run3_lossmap_with_crystals_fluka.py @@ -0,0 +1,168 @@ +# 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}" + +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() + + +# 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='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_ions=True) + + +# 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 2a747410..201c95d9 100644 --- a/xcoll/beam_elements/base.py +++ b/xcoll/beam_elements/base.py @@ -1502,8 +1502,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/beam_elements/elements_src/fluka_crystal.h b/xcoll/beam_elements/elements_src/fluka_crystal.h index 1864c151..e3176f37 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); } } diff --git a/xcoll/beam_elements/fluka.py b/xcoll/beam_elements/fluka.py index 04851e99..c98945d9 100644 --- a/xcoll/beam_elements/fluka.py +++ b/xcoll/beam_elements/fluka.py @@ -152,9 +152,13 @@ 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": + track_post(self, part) def __setattr__(self, name, value): import xcoll as xc @@ -336,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 diff --git a/xcoll/colldb.py b/xcoll/colldb.py index 3524f598..46493c05 100644 --- a/xcoll/colldb.py +++ b/xcoll/colldb.py @@ -14,7 +14,7 @@ import xtrack as xt from .beam_elements import (BlackAbsorber, BlackCrystal, EverestCollimator, EverestCrystal, - FlukaCollimator, Geant4Collimator, collimator_classes) + FlukaCollimator, FlukaCrystal, Geant4Collimator, collimator_classes) def _initialise_None(dct): diff --git a/xcoll/scattering_routines/fluka/engine.py b/xcoll/scattering_routines/fluka/engine.py index 3cf0e94e..3433aca8 100644 --- a/xcoll/scattering_routines/fluka/engine.py +++ b/xcoll/scattering_routines/fluka/engine.py @@ -168,7 +168,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) @@ -250,8 +250,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'] - 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__'): jaw = [jaw, -jaw] @@ -409,10 +414,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: diff --git a/xcoll/scattering_routines/fluka/fluka_input.py b/xcoll/scattering_routines/fluka/fluka_input.py index d42e702c..5e720255 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 @@ -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 @@ -211,6 +211,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 diff --git a/xcoll/scattering_routines/fluka/generic_prototype.py b/xcoll/scattering_routines/fluka/generic_prototype.py index 4f9f770d..e851b7e3 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 @@ -217,8 +222,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 + 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) @@ -260,10 +267,11 @@ 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)} -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} +RPP {fedb_tag}_B 0.0 {width*(100+10)} -{height*(100+10)/2} {height*(100+10)/2} -{length*(100+20)} {length*(100+20)} +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) @@ -278,10 +286,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 c4513962..f116f2f4 100644 --- a/xcoll/scattering_routines/fluka/includes.py +++ b/xcoll/scattering_routines/fluka/includes.py @@ -43,10 +43,14 @@ 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) + touches = kwargs.get('touches', 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 " @@ -134,7 +138,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: @@ -150,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) @@ -192,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) & @@ -214,9 +220,9 @@ 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]) + 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]) template += f"""\ @@ -232,10 +238,11 @@ def _assignmat_include_file(): return filename -def _beam_include_file(particle_ref): +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 + hi_prope = "*" if pdg_id in fluka_names: name = fluka_names[pdg_id] @@ -248,7 +255,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 * @@ -260,12 +266,46 @@ def _beam_include_file(particle_ref): * BEAMPOS * -* Only asking for loss map and touches map as in Sixtrack +""" + 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 +* 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'])}{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. +* 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 +* 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 +* ..+....1....+....2....+....3....+....4....+....5....+....6....+....7.. +USRICALL {format_fluka_float(bb_int['sigma_dpp'])} BBCOL_L +* +""" + else: + template += f"""\ +* 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&& """ + + with filename.open('w') as fp: fp.write(template) return filename @@ -275,6 +315,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 @@ -306,6 +347,7 @@ 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.. @@ -324,8 +366,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 *THRESHOL 0.0 0.0 8000.0 8000.0 0.0 0.0 """ with filename.open('w') as fp: @@ -335,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' @@ -354,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] != '*': @@ -421,6 +466,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 @@ -464,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 1baf93dd..a51a761b 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.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]}" def _expand(arr, available_capacity, dtype=float): @@ -208,7 +211,6 @@ def track_core(coll, part): idx_old = np.array([np.where(part.particle_id==idx)[0][0] for idx in new_pid[mask_existing]]) # list of indices - # import ipdb; ipdb.set_trace() # Sanity check assert np.all(part.particle_id[idx_old] == new_pid[mask_existing]) assert np.all(part.parent_particle_id[idx_old] == new_ppid[mask_existing])