-
Notifications
You must be signed in to change notification settings - Fork 7
/
pumas.c
15607 lines (14322 loc) · 626 KB
/
pumas.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
/*
* Copyright (C) 2017 Université Clermont Auvergne, CNRS/IN2P3, LPC
* Author: Valentin NIESS ([email protected])
*
* This software is a C library whose purpose is to transport high energy
* muons or taus in various media.
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU Lesser General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>
*/
/* The PUMAS API. */
#include "pumas.h"
#ifdef _WIN32
/* For rand_s on Windows */
#define _CRT_RAND_S
#endif
/* The C standard library. */
#include <assert.h>
#include <ctype.h>
#include <errno.h>
#include <float.h>
#include <math.h>
#include <stdarg.h>
#include <stddef.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
/* For debugging with gdb, on linux. */
#ifndef GDB_MODE
#define GDB_MODE 0
#endif
#if (GDB_MODE)
#ifndef __USE_GNU
#define __USE_GNU
#endif
#include <fenv.h>
#endif
/* Some tuning factors as macros. */
/**
* Number of schemes to tabulate for the computation of the energy loss.
*
* The NO_LOSS scheme has no tabulation and the detailed one uses the same
* tabulation than the hybrid scheme.
*/
#define N_SCHEMES 2
/**
* Order of expansion for the computation of the magnetic deflection,
* when using CSDA in a uniform medium of infinite extension.
*/
#define N_LARMOR_ORDERS 8
/**
* Number of inelastic processes for which Discrete Energy Losses (DELs) are
* simulated.
*/
#define N_DEL_PROCESSES 4
/**
* Default cutoff between Continuous Energy Loss (CEL) and DELs.
*/
#define DEFAULT_CUTOFF 5E-02
/**
* Default ratio for the elastic scattering.
*/
#define DEFAULT_ELASTIC_RATIO 5E-02
/**
* Exponent of the differential cross section approximation in Backward
* Monte-Carlo (BMC).
*/
#define BMC_ALPHA 2.0
/**
* Maximum path length for Elastic Hard Scattering (EHS) events, in kg/m^-2.
*/
#define EHS_PATH_MAX 1E+09
/**
* Default accuracy (relative step length limit) of the Monte Carlo integration.
*/
#define DEFAULT_ACCURACY 1E-02
/**
* Maximum deflection angle for a soft scattering event, in degrees.
*/
#define MAX_SOFT_ANGLE 1E+00
/**
* Grammage ratio for small steps, in CSDA mode.
*
* Below this ratio, the step energy loss is computed from the stopping power,
* using finite differences. Above, the interpolated grammage integral is used.
*/
#define EPSILON_X 3E-03
/**
* Minimum step size.
*/
#define STEP_MIN 1E-07
/**
* Tuning parameters for the tabulation of the DCS
*/
/* Minimum number of nodes on the low side */
#define DCS_MODEL_N_MIN 10
/* Nodes density per decade on the low side */
#define DCS_MODEL_PER_DECADE 15
/* Constant number of nodes on the high side */
#define DCS_MODEL_N_REVERSE 20
/* Cut value between low and high sides */
#define DCS_MODEL_X_REVERSE 0.6
/* Resolution for the high side sampling */
#define DCS_MODEL_DXMIN 1E-03
/* Upper cut for the tabulation */
#define DCS_MODEL_MAX_FRACTION 0.95
/**
* Minimum kinetic energy for using the DCS tabulation.
*/
#define DCS_MODEL_MIN_KINETIC 10.
/* Some constants, as macros. */
/**
* Fine-structure constant
*/
#define ALPHA_EM 7.2973525693E-03
/**
* Planck constant in GeV m.
*/
#define HBAR_C 1.973269804E-16
/**
* Bohr radius in m.
*/
#define BOHR_RADIUS 0.529177210903E-10
/**
* The muon decay length in m.
*/
#define MUON_C_TAU 658.654
/**
* The tau decay length in m.
*/
#define TAU_C_TAU 87.03E-06
/**
* Larmor magnetic factor in m^-1 GeV/c T^-1.
*/
#define LARMOR_FACTOR 0.299792458
/**
* The electron mass in GeV/c^2.
*/
#define ELECTRON_MASS 0.510998910E-03
/**
* The electron classical radius, in m.
*/
#define ELECTRON_RADIUS 2.817940285E-15
/**
* The muon mass in GeV/c^2.
*/
#define MUON_MASS 0.10565839
/**
* The tau mass in GeV/c^2.
*/
#define TAU_MASS 1.77682
/**
* The proton mass in GeV/c^2.
*/
#define PROTON_MASS 0.938272
/**
* The neutron mass in GeV/c^2.
*/
#define NEUTRON_MASS 0.939565
/**
* The pion mass in GeV/c^2.
*/
#define PION_MASS 0.13957018
#ifndef M_PI
/**
* Define pi, if unknown.
*/
#define M_PI 3.14159265358979323846
#endif
/**
* Avogadro's number
*/
#define AVOGADRO_NUMBER 6.02214076E+23
/**
* Default Bremsstrahlung model
*/
#define DEFAULT_BREMSSTRAHLUNG "SSR"
/**
* Default pair production model
*/
#define DEFAULT_PAIR_PRODUCTION "SSR"
/**
* Default photonuclear model
*/
#define DEFAULT_PHOTONUCLEAR "DRSS"
/* Helper macros for managing errors. */
#define ERROR_INITIALISE(caller) \
struct error_context error_data = {.code = PUMAS_RETURN_SUCCESS, \
.function = (pumas_function_t *)caller }; \
struct error_context * error_ = &error_data;
#define ERROR_MESSAGE(rc, message) \
error_format(error_, rc, __FILE__, __LINE__, message), \
error_raise(error_)
#define ERROR_FORMAT(rc, format, ...) \
error_format(error_, rc, __FILE__, __LINE__, format, __VA_ARGS__), \
error_raise(error_)
#define ERROR_REGISTER(rc, message) \
error_format(error_, rc, __FILE__, __LINE__, message)
#define ERROR_VREGISTER(rc, format, ...) \
error_format(error_, rc, __FILE__, __LINE__, format, __VA_ARGS__)
#define ERROR_RAISE() error_raise(error_)
#define ERROR_NULL_PHYSICS() \
ERROR_MESSAGE(PUMAS_RETURN_PHYSICS_ERROR, \
"a NULL physics pointer was provided")
#define ERROR_NOT_INITIALISED() \
ERROR_MESSAGE(PUMAS_RETURN_PHYSICS_ERROR, \
"the Physics has not been initialised")
#define ERROR_INVALID_SCHEME(scheme) \
ERROR_FORMAT(PUMAS_RETURN_INDEX_ERROR, \
"invalid energy loss scheme [%d]", scheme)
#define ERROR_INVALID_ELEMENT(element) \
ERROR_FORMAT( \
PUMAS_RETURN_INDEX_ERROR, "invalid element index [%d]", element)
#define ERROR_INVALID_MATERIAL(material) \
ERROR_FORMAT( \
PUMAS_RETURN_INDEX_ERROR, "invalid material index [%d]", material)
#define ERROR_REGISTER_MEMORY() \
ERROR_REGISTER(PUMAS_RETURN_MEMORY_ERROR, "could not allocate memory")
#define ERROR_REGISTER_EOF(path) \
ERROR_VREGISTER(PUMAS_RETURN_END_OF_FILE, \
"abnormal end of file when parsing `%s'", path)
#define ERROR_REGISTER_UNEXPECTED_TAG(tag, path, line) \
ERROR_VREGISTER(PUMAS_RETURN_FORMAT_ERROR, \
"unexpected XML tag `%s' [@%s:%d]", tag, path, line)
#define ERROR_REGISTER_TOO_LONG(path, line) \
ERROR_VREGISTER(PUMAS_RETURN_TOO_LONG, \
"XML node is too long [@%s:%d]", path, line)
#define ERROR_REGISTER_INVALID_XML_VALUE(value, path, line) \
ERROR_VREGISTER(PUMAS_RETURN_FORMAT_ERROR, \
"invalid XML value `%s' [@%s:%d]", value, path, line)
#define ERROR_REGISTER_INVALID_XML_ATTRIBUTE(attribute, node, path, line) \
ERROR_VREGISTER(PUMAS_RETURN_FORMAT_ERROR, \
"invalid XML attribute `%s' for element %s [@%s:%d]", attribute, \
node, path, line)
#define ERROR_REGISTER_NEGATIVE_DENSITY(material) \
ERROR_VREGISTER( \
PUMAS_RETURN_DENSITY_ERROR, "negative density for `%s'", material)
/* Function prototypes for the DCS implementations. */
/**
* Handle for a DCS computation function.
*/
struct atomic_element;
typedef double(dcs_function_t)(const struct pumas_physics * physics,
const struct atomic_element * element, double K, double q);
/**
* Handle for a polar angle sampling function.
*/
typedef double(polar_function_t)(const struct pumas_physics * physics,
struct pumas_context * context, const struct atomic_element * element,
double ki, double kf);
/**
* Generic doubly differential cross-section.
*/
typedef double ddcs_t(
double Z, double A, double m, double K, double q, double Q2);
/* A collection of low level flags. */
/**
* Node keys for the MDF files.
*/
enum mdf_key {
/** An atomic component. */
MDF_KEY_ATOMIC_COMPONENT,
/** A composite material. */
MDF_KEY_COMPOSITE,
/** A component of a composite material. */
MDF_KEY_COMPOSITE_COMPONENT,
/** The root PUMAS node. */
MDF_KEY_PUMAS,
/** An atomic element. */
MDF_KEY_ELEMENT,
/** A base material. */
MDF_KEY_MATERIAL,
/** Any other node. */
MDF_KEY_OTHER
};
/**
* Nodes hierarchy in MDFs.
*/
enum mdf_depth {
/** Outside of the root PUMAS node. */
MDF_DEPTH_EXTERN = 0,
/** Inside the root pumas node. */
MDF_DEPTH_ROOT,
/** Inside an element node. */
MDF_DEPTH_ELEMENT,
/** Inside a material node. */
MDF_DEPTH_MATERIAL,
/** Inside a composite material node. */
MDF_DEPTH_COMPOSITE
};
/**
* Tags for operations relative to the parsing of materials in MDFs.
*/
enum mdf_settings_operation {
/** Free the names table. */
MDF_INDEX_FREE = -1,
/** Initialise the names table. */
MDF_INDEX_INITIALISE,
/** Add a new base material to the table. */
MDF_INDEX_WRITE_MATERIAL,
/** Finalise a base material info. */
MDF_INDEX_FINALISE_MATERIAL,
/** Add a new composite material to the table. */
MDF_INDEX_INITIALISE_COMPOSITE,
/** Update a composite material with a new component. */
MDF_INDEX_UPDATE_COMPOSITE,
/** Finalise the composite material. */
MDF_INDEX_FINALISE_COMPOSITE
};
/* Low level data structures. */
/**
* Temporary data for the computation of the Coulomb scattering.
*
* The DCS is given by a product of 3 Wentzel distributions with 1 atomic
* screening parameter and 2 nuclear screening ones. A pole reduction allows
* the analytical computation of the various momenta: total cross-section, first
* transport path length, ...
*/
struct coulomb_data {
/** The partial cross section for hard events. */
double cs_hard;
/** The normalisation of the macroscopic cross-section. */
double normalisation;
/** The number of screening parameters */
int n_parameters;
/** The amplitudes of the atomic screening terms. */
double amplitude[3];
/** The atomic and nuclear screening parameters. */
double screening[4];
/** The 1st order coefficients of the pole reduction of the DCS. */
long double a[3];
/** The 2nd order coefficients of the pole reduction of the DCS. */
long double b[3];
/** The nuclear coefficients of the pole reduction of the DCS. */
long double c[4];
/** The spin correction factor. */
double fspin;
/**
* The parameters of the relativistic transform from Center of Mass (CM)
* frame to the Laboratory one.
*/
double fCM[2];
};
/**
* Temporary data for the rendering of a Coulomb DEL. The event is randomised
* using the inverse transform method. A root bracketing algorithm is used,
* seeded with an approximate solution.
*/
struct coulomb_workspace {
/** The index of the propagation material. */
int material;
/** The index of the hard scatterer element. */
int ihard;
/** The targeted cross section. */
double cs_h;
/** Placeholder for the table of Coulomb scattering data. */
struct coulomb_data data[];
};
/**
* Low level container for the local properties of a propagation medium.
*/
struct medium_locals {
/** The public API properties exposed to the end user. */
struct pumas_locals api;
/** A flag telling if the material has a magnetic field or not. */
int magnetized;
/** The local's physics */
const struct pumas_physics * physics;
};
/**
* Data for the default per context PRNG
*/
struct pumas_random_data {
/*
* Version tag for the random data format. Increment whenever the
* structure changes.
*/
#define RANDOM_BINARY_DUMP_TAG 0
/** The initial seed */
unsigned long seed;
/** Index in the PRNG buffer */
int index;
#define MT_PERIOD 624
/** PRNG buffer (Mersenne Twister) */
unsigned long buffer[MT_PERIOD];
};
/**
* The local data managed by a simulation context.
*/
struct simulation_context {
/** The public API settings exposed to the end user. */
struct pumas_context api;
/** Handle for Physics tables. */
const struct pumas_physics * physics;
/** Lifetime limit for the decay process. */
double lifetime;
/**
* Last kinetic energy indices used in the tables.
*
* We keep a memory of the last indices accessed in tables. The memory
* has a depth of 2 which should be a good compromise since most use
* cases consider differences between an initial and final state.
*/
int index_K_last[2];
/** Last grammage indices used in the tables. */
int index_X_last[2];
/** Flag for the first step, for integration of various quantities. */
int step_first;
/** Tracking of stepping events. */
enum pumas_event step_event;
/** The expected next event during the stepping. */
enum pumas_event step_foreseen;
/** The kinetic limit converted to grammage. */
double step_X_limit;
/** The scaterring 1st transport path length of the previous step. */
double step_invlb1;
/** Data for the default PRNG. */
struct pumas_random_data * random_data;
/** Flag for the parity check of the Gaussian random generator.
*
* Gaussian variates are generated in pair using the Box-Muller
* transform.
*/
int randn_done;
/** The next Gaussian variate. */
double randn_next;
/**
* Pointer to the worspace for the temporary storage of intermediary
* computations.
*/
struct coulomb_workspace * workspace;
/** Size of the user extended memory. */
int extra_memory;
/**
* Placeholder for variable data storage with -double- memory alignment.
*
* Extra bytes are allocated for the workspace, the error stack and
* any extended memory for end user usage.
*/
long double data[];
};
/**
* Handle for a stack of recorded frames.
*
* The frames are allocated in bunches. The bunches are managed as a chained
* list.
*/
struct frame_stack {
/** The memory size left, in bytes. */
int size;
/** Pointer to the next memory segment. */
struct frame_stack * next;
/** Pointer to the first frame in the stack. */
struct pumas_frame * frame;
/** Placeholder for frames */
struct pumas_frame frames[];
};
/**
* Low level container for a frame recorder.
*/
struct frame_recorder {
/** The public API data exposed to the end user. */
struct pumas_recorder api;
/** Link to the last record. */
struct pumas_frame * last;
/** Link to the 1st entry of the chained list of stacks. */
struct frame_stack * stack;
/** Placeholder for extra data. */
double data[];
};
/**
* Data relative to an atomic element.
*/
struct atomic_element {
/** The element atomic number. */
double Z;
/** The element atomic mass, in g/mol. */
double A;
/** The element Mean Excitation Energy (MEE), in eV. */
double I;
/** The element tabulation index */
int index;
/** The element name. */
char * name;
/** Placeholder for user data with -double- memory alignment. */
double data[];
};
/**
* An atomic component of a material.
*/
struct material_component {
/** The atomic element index. */
int element;
/** The element mass fraction in the material. */
double fraction;
};
/**
* A component of a composite material.
*/
struct composite_component {
/** The constituent base material index. */
int material;
/** The component mass fraction in the composite. */
double fraction;
};
/**
* Handle for a composite material.
*/
struct composite_material {
/** The number of sub components. */
int n_components;
/** Placeholder for the sub components' data. */
struct composite_component component[];
};
/**
* Temporary data for the parsing of a MDF.
*/
struct mdf_buffer {
/** Flag for the dry initialisation mode. */
int dry_mode;
/** Handle to the MDF. */
FILE * fid;
/** Current position in the read buffer. */
char * pos;
/** The number of bytes left to be read. */
int left;
/** The total size of the read buffer. */
int size;
/** The current line number in the MDF file. */
int line;
/** Current node depth during the MDF parsing. */
enum mdf_depth depth;
/** Counter for the number of base materials in a composite. */
int materials_in;
/** Counter for the number of elements in a material. */
int elements_in;
/** Pointer to the current MDF. */
const char * mdf_path;
/** The number of kinetic energy rows in a dE/dX file. */
int n_energies;
/** The total number of materials, base and composites. */
int n_materials;
/** The number of composite materials. */
int n_composites;
/** The number of atomic elements. */
int n_elements;
/** The header length of dE/dX files. */
int n_energy_loss_header;
/** The total number of atomic element components. */
int n_components;
/** The maximum number of atomic elements in a single material. */
int max_components;
/** The total byte size for the storage of composite materials. */
int size_composite;
/** The total size of the path to the dE/dX files. */
int size_dedx_path;
/** The total size of elements names. */
int size_elements_names;
/** The total size of materials names. */
int size_materials_names;
/** Placeholder for the read buffer. */
char data[];
};
/*!
* Pointers to the data fields of a node in a MDF file.
*
* The pointers refer to an mdf_buffer object. If the buffer is refilled the
* links are no mode valid.
*/
struct mdf_node {
/** The node key. */
enum mdf_key key;
/** Flag telling if this is a head node. */
int head;
/** Flag telling if this is a tail node. */
int tail;
/** The first attribute names. */
union attribute1 {
/** The node name. */
char * name;
} at1;
/** The second attribute names. */
union attribute2 {
/** The energy loss file name. */
char * file;
/** The atomic number. */
char * Z;
/** The mass fraction. */
char * fraction;
} at2;
/** The third attribute names. */
union attribute3 {
/** The atomic mass. */
char * A;
/** The composite component's density. */
char * density;
} at3;
/** The fourth attribute names. */
union attribute4 {
/** The mean excitation energy (MEE). */
char * I;
} at4;
};
/**
* Temporary data for a DEL event.
*/
struct del_info {
union {
/** The energy transfer. */
double Q;
/** The Monte-Carlo weight factor. */
double weight;
}
/** Data specific to reverse Monte-Carlo. */
reverse;
/** The index of the inelastic sub-process. */
int process;
/** The index of the target element. */
int element;
};
/**
* Global data shared by all simulation contexts.
*/
struct pumas_physics {
/*
* Version tag for the physics data format. Increment whenever the
* structure changes.
*/
#define PHYSICS_BINARY_DUMP_TAG 13
/** The total byte size of the shared data. */
int size;
/** The number of kinetic energy values in the dE/dX tables. */
int n_energies;
/** The total number of materials, basic and composites. */
int n_materials;
/** The total number of composite materials. */
int n_composites;
/** The total number of declared atomic_elements. */
int n_elements;
/** The total number of atomic components in materials. */
int n_components;
/** The maximum number of atomic components in a single material. */
int max_components;
/** The number of header lines in a dE/dX table. */
int n_energy_loss_header;
/**
* Offset for the tabulation of DCS and their polynomial
* approximations.
*/
int dcs_model_offset;
/** The number of tabulated dcs_values */
int n_table_dcs;
/** The transported particle type. */
enum pumas_particle particle;
/** The transported particle decay length, in m. */
double ctau;
/** The transported particle rest mass, in GeV. */
double mass;
/** The relative cutoff between CEL and DELs. */
double cutoff;
/** Ratio of EHS path length w.r.t. the first transport path length. */
double elastic_ratio;
/** Path to the current MDF. */
char * mdf_path;
/** Path where the dE/dX files are stored. */
char * dedx_path;
/** Names of the dE/dX files. */
char ** dedx_filename;
/** The tabulated values of the kinetic energy. */
double * table_K;
double * table_K_dX;
double * table_K_dNI_in;
double * table_K_dNI_el;
/** The tabulated values of the total grammage (CSDA range). */
double * table_X;
double * table_X_dK;
double * table_X_dT;
/** The tabulated values of the total proper time. */
double * table_T;
double * table_T_dK;
/** The tabulated values of the average energy loss. */
double * table_dE;
double * table_dE_dK;
/** The tabulated values of the energy straggling. */
double * table_Omega;
double * table_Omega_dK;
/** The tabulated values of the EHS number of interaction lengths. */
double * table_NI_el;
double * table_NI_el_dK;
/**
* The tabulated values of the number of interaction lengths for
* inelastic DELs.
*/
double * table_NI_in;
double * table_NI_in_dK;
/** The tabulated cross section values. */
double * table_CS;
double * table_CS_dK;
/** The tabulated cross section fractions. */
double * table_CSf;
double * table_CSf_dK;
/** The tabulated cross section normalisation. */
double * table_CSn;
double * table_CSn_dK;
/* The tabulated data for DCSs **/
float * table_DCS;
float * table_DCS_x;
float * table_DCS_envelope;
/** The element wise fractional threshold for DELs. */
double * table_Xt;
/** The total kinetic threshold for DELs. */
double * table_Kt;
/**
* The tabulated values of the magnetic deflection momenta, within
* the CSDA.
*/
double * table_Li;
double * table_Li_dK;
/** The last tabulated value of the ionisation energy loss. */
double * table_a_max;
/** The last tabulated value of the radiative energy loss. */
double * table_b_max;
/**
* The tabulated angular cutoff values for the splitting of Coulomb
* scattering.
*/
double * table_Mu0;
double * table_Mu0_dK;
/** The tabulated interaction lengths for DEL Coulomb events. */
double * table_Lb;
double * table_Lb_dK;
/** The tabulated multiple scattering 1st moment. */
double * table_Ms1;
double * table_Ms1_dK;
/** The number of elements in a material. */
int * elements_in;
/** The reference density of a material. */
double * material_density;
/** The relative electronic density of a material. */
double * material_ZoA;
/** The mean excitation energy of a base material. */
double * material_I;
/** The Sternheimer scaling ratio of a material. */
double * material_aS;
/** The properties of an atomic element . */
struct atomic_element ** element;
/** The composition of a base material. */
struct material_component ** composition;
/** The composition of a composite material. */
struct composite_material ** composite;
/** The material names. */
char ** material_name;
/** The Bremsstrahlung model. */
char * model_bremsstrahlung;
/** The pair_production model. */
char * model_pair_production;
/** The photonuclear model. */
char * model_photonuclear;
/** The Bremsstrahlung DCS. */
pumas_dcs_t * dcs_bremsstrahlung;
/** The pair_production DCS. */
pumas_dcs_t * dcs_pair_production;
/** The photonuclear DCS. */
pumas_dcs_t * dcs_photonuclear;
/**
* Placeholder for shared data storage with -double- memory alignment.
*/
double data[];
};
struct error_context {
enum pumas_return code;
pumas_function_t * function;
#define ERROR_MSG_LENGTH 1024
char message[ERROR_MSG_LENGTH];
};
/**
* Handle for an atomic element within a material.
*
* This structure is a proxy exposing some data of an atomic element within
* the material last processed by `pumas_tabulation_tabulate`.
*/
struct physics_element {
/** Linked list pointer to the previous element. */
struct physics_element * prev;
/** Linked list pointer to the next element. */
struct physics_element * next;
/** The element index. */
int index;
/** The mass fraction of the element in the current material. */
double fraction;
};
/**
* Handle for tabulation data.
*
* This structure gathers data related to the tabulation of the energy loss of
* materials with the `physics_tabulate` function. **Note** that the two last
* parameters: *path* and *elements* should not be set directly. They are filled
* in (updated) by the `physics_tabulate` function. Note also that if no energy
* grid is provided then a default one is set.
*
* **Warning**: the energy grid should not be changed between successive calls
* to `physics_tabulate`. If a new energy grid is needed then a new
* `physics_tabulation_data` object must be created.
*/
struct physics_tabulation_data {
/** The number of kinetic energy values to tabulate. Providing a value
* of zero or less results in a default energy grid being set.
*/
int n_energies;
/** Array of kinetic energy values to tabulate. Providing a `NULL`
* value results in a default energy grid being set.
*/
double * energy;
/** Flag to enable overwriting an existing energy loss file. */
int overwrite;
/** Path to a directory where the tabulation should be written. */
char * outdir;
/** Index of the material to tabulate */
int material;
/** Path to the energy loss file of the last tabulated material. */
char * path;
/** List of atomic elements contained in the tabulated material(s). */
struct physics_element * elements;
/** Temporary work data */
double * work;
};
/**
* Default error handler.
*/
static void default_error_handler(
enum pumas_return rc, pumas_function_t * caller, const char * message)
{
/* Dump the error summary */
fputs("pumas: library error. See details below\n", stderr);
fprintf(stderr, "error: %s\n", message);
/* Exit to the OS */
exit(EXIT_FAILURE);
}
/**
* Shared data for the error handling.
*/
static struct {
pumas_handler_cb * handler;
int catch;
struct error_context catch_error;
} s_error = { &default_error_handler, 0 };
/* Prototypes of low level static functions. */
/**
* Encapsulations of the tabulated CEL and DEL properties.
*/
static double cel_grammage(const struct pumas_physics * physics,
struct pumas_context * context, enum pumas_mode scheme, int material,
double kinetic);
static double cel_grammage_as_time(const struct pumas_physics * physics,
struct pumas_context * context, enum pumas_mode scheme, int material,
double time);
static double cel_proper_time(const struct pumas_physics * physics,
struct pumas_context * context, enum pumas_mode scheme, int material,
double kinetic);
static double cel_kinetic_energy(const struct pumas_physics * physics,
struct pumas_context * context, enum pumas_mode scheme, int material,
double grammage);
static double cel_energy_loss(const struct pumas_physics * physics,
struct pumas_context * context, enum pumas_mode scheme, int material,
double kinetic);
static double cel_straggling(const struct pumas_physics * physics,
struct pumas_context * context, int material, double kinetic);
static double cel_magnetic_rotation(const struct pumas_physics * physics,
struct pumas_context * context, int material, double kinetic);
static double del_cross_section(const struct pumas_physics * physics,
struct pumas_context * context, int material, double kinetic);
static double del_interaction_length(const struct pumas_physics * physics,
struct pumas_context * context, int material, double kinetic);
static double del_kinetic_from_interaction_length(
const struct pumas_physics * physics, struct pumas_context * context,
int material, double nI);
static double ehs_interaction_length(const struct pumas_physics * physics,
struct pumas_context * context, enum pumas_mode scheme, int material,
double kinetic);
static double ehs_kinetic_from_interaction_length(
const struct pumas_physics * physics, struct pumas_context * context,
enum pumas_mode scheme, int material, double nI);
/**
* Routines related to DCS: implementation and handling.
*/
static inline dcs_function_t * dcs_get(int process);
static inline int dcs_get_index(dcs_function_t * dcs_func);
static inline void dcs_get_range(int process, double Z, double mass,
double energy, double * qmin, double * qmax);
static double dcs_bremsstrahlung(const struct pumas_physics * physics,
const struct atomic_element * element, double K, double q);
static double dcs_bremsstrahlung_transport(const struct pumas_physics * physics,
const struct atomic_element * element, double K, double cutoff);
static double dcs_pair_production(const struct pumas_physics * physics,
const struct atomic_element * element, double K, double q);
static double dcs_pair_production_transport(
const struct pumas_physics * physics, const struct atomic_element * element,
double K, double cutoff);
static inline void dcs_pair_production_range(
double Z, double m, double K, double * qmin, double * qmax);
static double dcs_photonuclear(const struct pumas_physics * physics,
const struct atomic_element * element, double K, double q);
static inline void dcs_photonuclear_range(
double m, double K, double * qmin, double * qmax);
static inline int dcs_photonuclear_check(double m, double K, double q);
static inline ddcs_t * dcs_photonuclear_ddcs(
const struct pumas_physics * physics, pumas_dcs_t ** dcs);
static double dcs_photonuclear_transport(const struct pumas_physics * physics,
const struct atomic_element * element, double K, double cutoff);
static double dcs_ionisation(const struct pumas_physics * physics,
const struct atomic_element * element, double K, double q);
static double dcs_ionisation_integrate(const struct pumas_physics * physics,
int mode, const struct atomic_element * element, double K, double xlow,
double xhigh);
static double dcs_ionisation_randomise(const struct pumas_physics * physics,
struct pumas_context * context, const struct atomic_element * element,
double K, double xlow);
static double dcs_evaluate(const struct pumas_physics * physics,
struct pumas_context * context, dcs_function_t * dcs_func,
const struct atomic_element * element, double K, double q);
/**
* Implementations of polar angle distributions and accessor.
*/
static inline polar_function_t * polar_get(int process);
static double polar_bremsstrahlung(const struct pumas_physics * physics,
struct pumas_context * context, const struct atomic_element * element,
double ki, double kf);
static double polar_pair_production(const struct pumas_physics * physics,
struct pumas_context * context, const struct atomic_element * element,
double ki, double kf);
static double polar_photonuclear(const struct pumas_physics * physics,
struct pumas_context * context, const struct atomic_element * element,
double ki, double kf);
static double polar_ionisation(const struct pumas_physics * physics,
struct pumas_context * context, const struct atomic_element * element,
double ki, double kf);
/**
* Low level routines for the propagation in matter.
*/
static enum pumas_event transport_with_csda(
const struct pumas_physics * physics, struct pumas_context * context,
struct pumas_state * state, struct pumas_medium * medium,
struct medium_locals * locals, struct error_context * error_);
static enum pumas_return transport_csda_deflect(
const struct pumas_physics * physics, struct pumas_context * context,
struct pumas_state * state, struct pumas_medium * medium,
struct medium_locals * locals, double ki, double distance,
struct error_context * error_);
static enum pumas_return csda_magnetic_transport(
const struct pumas_physics * physics, struct pumas_context * context,
int material, double density, double magnet, double charge, double kinetic,
double phase, double * x, double * y, double * z,
struct error_context * error_);
static enum pumas_event transport_with_stepping(
const struct pumas_physics * physics, struct pumas_context * context,