CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
NuclearInteractionFTFSimulator.cc
Go to the documentation of this file.
1 // system headers
2 #include <mutex>
3 
4 // Framework Headers
5 #include "FWCore/MessageLogger/interface/MessageLogger.h"
6 #include "FWCore/Utilities/interface/thread_safety_macros.h"
7 
8 // Fast Sim headers
9 #include "FastSimulation/MaterialEffects/interface/NuclearInteractionFTFSimulator.h"
10 #include "FastSimulation/MaterialEffects/interface/CMSDummyDeexcitation.h"
11 #include "FastSimulation/Utilities/interface/RandomEngineAndDistribution.h"
12 #include "FastSimulation/ParticlePropagator/interface/ParticlePropagator.h"
13 #include "FastSimulation/Particle/interface/makeParticle.h"
14 
15 // Geant4 headers
16 #include "G4ParticleDefinition.hh"
17 #include "G4DynamicParticle.hh"
18 #include "G4Track.hh"
19 #include "G4Step.hh"
20 #include "G4StepPoint.hh"
21 #include "G4TheoFSGenerator.hh"
22 #include "G4FTFModel.hh"
23 #include "G4ExcitedStringDecay.hh"
24 #include "G4LundStringFragmentation.hh"
25 #include "G4GeneratorPrecompoundInterface.hh"
26 #include "G4CascadeInterface.hh"
27 #include "G4DiffuseElastic.hh"
28 
29 #include "G4Proton.hh"
30 #include "G4Neutron.hh"
31 #include "G4PionPlus.hh"
32 #include "G4PionMinus.hh"
33 #include "G4AntiProton.hh"
34 #include "G4AntiNeutron.hh"
35 #include "G4KaonPlus.hh"
36 #include "G4KaonMinus.hh"
37 #include "G4KaonZeroLong.hh"
38 #include "G4KaonZeroShort.hh"
39 #include "G4KaonZero.hh"
40 #include "G4AntiKaonZero.hh"
41 #include "G4GenericIon.hh"
42 
43 #include "G4Lambda.hh"
44 #include "G4OmegaMinus.hh"
45 #include "G4SigmaMinus.hh"
46 #include "G4SigmaPlus.hh"
47 #include "G4SigmaZero.hh"
48 #include "G4XiMinus.hh"
49 #include "G4XiZero.hh"
50 #include "G4AntiLambda.hh"
51 #include "G4AntiOmegaMinus.hh"
52 #include "G4AntiSigmaMinus.hh"
53 #include "G4AntiSigmaPlus.hh"
54 #include "G4AntiSigmaZero.hh"
55 #include "G4AntiXiMinus.hh"
56 #include "G4AntiXiZero.hh"
57 #include "G4AntiAlpha.hh"
58 #include "G4AntiDeuteron.hh"
59 #include "G4AntiTriton.hh"
60 #include "G4AntiHe3.hh"
61 
62 #include "G4Material.hh"
63 #include "G4DecayPhysics.hh"
64 #include "G4ParticleTable.hh"
65 #include "G4IonTable.hh"
66 #include "G4ProcessManager.hh"
67 #include "G4PhysicsLogVector.hh"
68 #include "G4SystemOfUnits.hh"
69 
70 static std::once_flag initializeOnce;
71 CMS_THREAD_GUARD(initializeOnce) const G4ParticleDefinition* NuclearInteractionFTFSimulator::theG4Hadron[] = {nullptr};
72 CMS_THREAD_GUARD(initializeOnce) int NuclearInteractionFTFSimulator::theId[] = {0};
73 
74 const double fact = 1.0 / CLHEP::GeV;
75 
76 // inelastic interaction length corrections per particle and energy
77 const double corrfactors_inel[numHadrons][npoints] = {
78  {1.0872, 1.1026, 1.111, 1.111, 1.0105, 0.97622, 0.9511, 0.9526, 0.97591, 0.99277, 1.0099,
79  1.015, 1.0217, 1.0305, 1.0391, 1.0438, 1.0397, 1.0328, 1.0232, 1.0123, 1.0},
80  {1.0416, 1.1044, 1.1467, 1.1273, 1.026, 0.99085, 0.96572, 0.96724, 0.99091, 1.008, 1.0247,
81  1.0306, 1.0378, 1.0427, 1.0448, 1.0438, 1.0397, 1.0328, 1.0232, 1.0123, 1.0},
82  {0.5308, 0.53589, 0.67059, 0.80253, 0.82341, 0.79083, 0.85967, 0.90248, 0.93792, 0.9673, 1.0034,
83  1.022, 1.0418, 1.0596, 1.0749, 1.079, 1.0704, 1.0576, 1.0408, 1.0214, 1.0},
84  {0.49107, 0.50571, 0.64149, 0.77209, 0.80472, 0.78166, 0.83509, 0.8971, 0.93234, 0.96154, 0.99744,
85  1.0159, 1.0355, 1.0533, 1.0685, 1.0732, 1.0675, 1.0485, 1.0355, 1.0191, 1.0},
86  {1.9746, 1.7887, 1.5645, 1.2817, 1.0187, 0.95216, 0.9998, 1.035, 1.0498, 1.0535, 1.0524,
87  1.0495, 1.0461, 1.0424, 1.0383, 1.0338, 1.0287, 1.0228, 1.0161, 1.0085, 1.0},
88  {0.46028, 0.59514, 0.70355, 0.70698, 0.62461, 0.65103, 0.71945, 0.77753, 0.83582, 0.88422, 0.92117,
89  0.94889, 0.96963, 0.98497, 0.99596, 1.0033, 1.0075, 1.0091, 1.0081, 1.005, 1.0},
90  {0.75016, 0.89607, 0.97185, 0.91083, 0.77425, 0.77412, 0.8374, 0.88848, 0.93104, 0.96174, 0.98262,
91  0.99684, 1.0065, 1.0129, 1.0168, 1.0184, 1.018, 1.0159, 1.0121, 1.0068, 1.0},
92  {0.75016, 0.89607, 0.97185, 0.91083, 0.77425, 0.77412, 0.8374, 0.88848, 0.93104, 0.96174, 0.98262,
93  0.99684, 1.0065, 1.0129, 1.0168, 1.0184, 1.018, 1.0159, 1.0121, 1.0068, 1.0},
94  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.0},
95  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.0},
96  {1.1006, 1.1332, 1.121, 1.1008, 1.086, 1.077, 1.0717, 1.0679, 1.0643, 1.0608, 1.057,
97  1.053, 1.0487, 1.0441, 1.0392, 1.0338, 1.028, 1.0218, 1.015, 1.0078, 1.0},
98  {1.1318, 1.1255, 1.1062, 1.0904, 1.0802, 1.0742, 1.0701, 1.0668, 1.0636, 1.0602, 1.0566,
99  1.0527, 1.0485, 1.044, 1.0391, 1.0337, 1.028, 1.0217, 1.015, 1.0078, 1.0},
100  {1.1094, 1.1332, 1.1184, 1.0988, 1.0848, 1.0765, 1.0714, 1.0677, 1.0642, 1.0607, 1.0569,
101  1.053, 1.0487, 1.0441, 1.0391, 1.0338, 1.028, 1.0218, 1.015, 1.0078, 1.0},
102  {1.1087, 1.1332, 1.1187, 1.099, 1.0849, 1.0765, 1.0715, 1.0677, 1.0642, 1.0607, 1.057,
103  1.053, 1.0487, 1.0441, 1.0391, 1.0338, 1.028, 1.0218, 1.015, 1.0078, 1.0},
104  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.0},
105  {1.1192, 1.132, 1.1147, 1.0961, 1.0834, 1.0758, 1.0711, 1.0674, 1.064, 1.0606, 1.0569,
106  1.0529, 1.0486, 1.0441, 1.0391, 1.0338, 1.028, 1.0218, 1.015, 1.0078, 1.0},
107  {1.1188, 1.1321, 1.1149, 1.0963, 1.0834, 1.0758, 1.0711, 1.0675, 1.0641, 1.0606, 1.0569,
108  1.0529, 1.0486, 1.0441, 1.0391, 1.0338, 1.028, 1.0218, 1.015, 1.0078, 1.0},
109  {0.50776, 0.5463, 0.5833, 0.61873, 0.65355, 0.68954, 0.72837, 0.7701, 0.81267, 0.85332, 0.89037,
110  0.92329, 0.95177, 0.97539, 0.99373, 1.0066, 1.014, 1.0164, 1.0144, 1.0087, 1.0},
111  {0.50787, 0.5464, 0.58338, 0.6188, 0.65361, 0.6896, 0.72841, 0.77013, 0.8127, 0.85333, 0.89038,
112  0.92329, 0.95178, 0.9754, 0.99373, 1.0066, 1.014, 1.0164, 1.0144, 1.0087, 1.0},
113  {1.1006, 1.1332, 1.121, 1.1008, 1.086, 1.077, 1.0717, 1.0679, 1.0643, 1.0608, 1.057,
114  1.053, 1.0487, 1.0441, 1.0392, 1.0338, 1.028, 1.0218, 1.015, 1.0078, 1.0},
115  {1.1318, 1.1255, 1.1062, 1.0904, 1.0802, 1.0742, 1.0701, 1.0668, 1.0636, 1.0602, 1.0566,
116  1.0527, 1.0485, 1.044, 1.0391, 1.0337, 1.028, 1.0217, 1.015, 1.0078, 1.0},
117  {1.1094, 1.1332, 1.1184, 1.0988, 1.0848, 1.0765, 1.0714, 1.0677, 1.0642, 1.0607, 1.0569,
118  1.053, 1.0487, 1.0441, 1.0391, 1.0338, 1.028, 1.0218, 1.015, 1.0078, 1.0},
119  {1.1087, 1.1332, 1.1187, 1.099, 1.0849, 1.0765, 1.0715, 1.0677, 1.0642, 1.0607, 1.057,
120  1.053, 1.0487, 1.0441, 1.0391, 1.0338, 1.028, 1.0218, 1.015, 1.0078, 1.0},
121  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.0},
122  {1.1192, 1.132, 1.1147, 1.0961, 1.0834, 1.0758, 1.0711, 1.0674, 1.064, 1.0606, 1.0569,
123  1.0529, 1.0486, 1.0441, 1.0391, 1.0338, 1.028, 1.0218, 1.015, 1.0078, 1.0},
124  {1.1188, 1.1321, 1.1149, 1.0963, 1.0834, 1.0758, 1.0711, 1.0675, 1.0641, 1.0606, 1.0569,
125  1.0529, 1.0486, 1.0441, 1.0391, 1.0338, 1.028, 1.0218, 1.015, 1.0078, 1.0},
126  {0.47677, 0.51941, 0.56129, 0.60176, 0.64014, 0.67589, 0.70891, 0.73991, 0.77025, 0.80104, 0.83222,
127  0.86236, 0.8901, 0.91518, 0.9377, 0.95733, 0.97351, 0.98584, 0.9942, 0.99879, 1.0},
128  {0.49361, 0.53221, 0.56976, 0.60563, 0.63954, 0.67193, 0.70411, 0.73777, 0.77378, 0.81114, 0.84754,
129  0.88109, 0.91113, 0.93745, 0.95974, 0.97762, 0.99081, 0.99929, 1.0033, 1.0034, 1.0},
130  {0.4873, 0.52744, 0.56669, 0.60443, 0.64007, 0.67337, 0.70482, 0.73572, 0.76755, 0.80086, 0.83456,
131  0.86665, 0.8959, 0.92208, 0.94503, 0.96437, 0.97967, 0.99072, 0.99756, 1.0005, 1.0},
132  {0.48729, 0.52742, 0.56668, 0.60442, 0.64006, 0.67336, 0.70482, 0.73571, 0.76754, 0.80086, 0.83455,
133  0.86665, 0.8959, 0.92208, 0.94503, 0.96437, 0.97967, 0.99072, 0.99756, 1.0005, 1.0},
134 };
135 
136 // elastic interaction length corrections per particle and energy
137 const double corrfactors_el[numHadrons][npoints] = {
138  {0.58834, 1.1238, 1.7896, 1.4409, 0.93175, 0.80587, 0.80937, 0.83954, 0.87453, 0.91082, 0.94713,
139  0.98195, 1.0134, 1.0397, 1.0593, 1.071, 1.0739, 1.0678, 1.053, 1.03, 1.0},
140  {0.40938, 0.92337, 1.3365, 1.1607, 1.008, 0.82206, 0.81163, 0.79489, 0.82919, 0.91812, 0.96688,
141  1.0225, 1.0734, 1.0833, 1.0874, 1.0854, 1.0773, 1.0637, 1.0448, 1.0235, 1.0},
142  {0.43699, 0.42165, 0.46594, 0.64917, 0.85314, 0.80782, 0.83204, 0.91162, 1.0155, 1.0665, 1.0967,
143  1.1125, 1.1275, 1.1376, 1.1464, 1.1477, 1.1312, 1.1067, 1.0751, 1.039, 1.0},
144  {0.3888, 0.39527, 0.43921, 0.62834, 0.8164, 0.79866, 0.82272, 0.90163, 1.0045, 1.055, 1.0849,
145  1.1005, 1.1153, 1.1253, 1.134, 1.1365, 1.1255, 1.0895, 1.0652, 1.0348, 1.0},
146  {0.32004, 0.31119, 0.30453, 0.30004, 0.31954, 0.40148, 0.5481, 0.74485, 0.99317, 1.1642, 1.2117,
147  1.2351, 1.2649, 1.3054, 1.375, 1.4992, 1.4098, 1.3191, 1.2232, 1.118, 1.0},
148  {0.10553, 0.14623, 0.20655, 0.26279, 0.19996, 0.40125, 0.5139, 0.71271, 0.89269, 1.0108, 1.1673,
149  1.3052, 1.4149, 1.429, 1.4521, 1.4886, 1.4006, 1.3116, 1.2177, 1.1151, 1.0},
150  {0.106, 0.14692, 0.20755, 0.26257, 0.20089, 0.40236, 0.51452, 0.71316, 0.89295, 1.0109, 1.1673,
151  1.3053, 1.4149, 1.429, 1.4521, 1.4886, 1.4006, 1.3116, 1.2177, 1.1151, 1.0},
152  {0.31991, 0.31111, 0.30445, 0.30004, 0.31995, 0.40221, 0.54884, 0.74534, 0.99364, 1.1644, 1.2117,
153  1.2351, 1.265, 1.3054, 1.375, 1.4992, 1.4098, 1.3191, 1.2232, 1.118, 1.0},
154  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.0},
155  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.0},
156  {0.37579, 0.39922, 0.37445, 0.32631, 0.39002, 0.42161, 0.54251, 0.69127, 0.90332, 1.0664, 1.1346,
157  1.1481, 1.1692, 1.2036, 1.2625, 1.3633, 1.2913, 1.2215, 1.1516, 1.0788, 1.0},
158  {0.31756, 0.33409, 0.25339, 0.35525, 0.52989, 0.63382, 0.7453, 0.93505, 1.1464, 1.2942, 1.3161,
159  1.328, 1.3393, 1.3525, 1.374, 1.4051, 1.3282, 1.2523, 1.1745, 1.0916, 1.0},
160  {0.38204, 0.39694, 0.36502, 0.33367, 0.39229, 0.43119, 0.54898, 0.70169, 0.91004, 1.0696, 1.1348,
161  1.1483, 1.1694, 1.2038, 1.2627, 1.3632, 1.2913, 1.2215, 1.1516, 1.0788, 1.0},
162  {0.38143, 0.39716, 0.36609, 0.33294, 0.39207, 0.43021, 0.54834, 0.70066, 0.90945, 1.0693, 1.1348,
163  1.1482, 1.1694, 1.2038, 1.2627, 1.3632, 1.2913, 1.2215, 1.1516, 1.0788, 1.0},
164  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.0},
165  {0.29564, 0.32645, 0.29986, 0.30611, 0.48808, 0.59902, 0.71207, 0.8832, 1.1164, 1.2817, 1.3154,
166  1.3273, 1.3389, 1.3521, 1.3736, 1.4056, 1.3285, 1.2524, 1.1746, 1.0916, 1.0},
167  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.0},
168  {0.3265, 0.3591, 0.39232, 0.42635, 0.46259, 0.50365, 0.55244, 0.61014, 0.67446, 0.74026, 0.80252,
169  0.85858, 0.90765, 0.94928, 0.9827, 1.0071, 1.0221, 1.0279, 1.0253, 1.0155, 1.0},
170  {0.13808, 0.15585, 0.17798, 0.2045, 0.22596, 0.25427, 0.33214, 0.44821, 0.5856, 0.74959, 0.89334,
171  1.0081, 1.0964, 1.1248, 1.173, 1.2548, 1.1952, 1.1406, 1.0903, 1.0437, 1.0},
172  {0.20585, 0.23253, 0.26371, 0.28117, 0.30433, 0.35417, 0.44902, 0.58211, 0.73486, 0.90579, 1.0395,
173  1.1488, 1.2211, 1.2341, 1.2553, 1.2877, 1.2245, 1.1654, 1.1093, 1.0547, 1.0},
174  {0.2852, 0.32363, 0.31419, 0.35164, 0.45463, 0.54331, 0.66908, 0.81735, 0.98253, 1.1557, 1.2557,
175  1.3702, 1.4186, 1.401, 1.374, 1.3325, 1.2644, 1.1991, 1.1348, 1.0694, 1.0},
176  {0.20928, 0.23671, 0.2664, 0.28392, 0.30584, 0.35929, 0.45725, 0.5893, 0.74047, 0.9101, 1.0407,
177  1.1503, 1.2212, 1.2342, 1.2554, 1.2876, 1.2245, 1.1654, 1.1093, 1.0547, 1.0},
178  {0.11897, 0.13611, 0.15796, 0.1797, 0.21335, 0.26367, 0.34705, 0.46115, 0.6016, 0.7759, 0.91619,
179  1.0523, 1.1484, 1.1714, 1.2098, 1.2721, 1.2106, 1.1537, 1.1004, 1.0496, 1.0},
180  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.0},
181  {0.26663, 0.30469, 0.32886, 0.33487, 0.41692, 0.51616, 0.63323, 0.78162, 0.95551, 1.1372, 1.2502,
182  1.3634, 1.4189, 1.4013, 1.3743, 1.3329, 1.2646, 1.1992, 1.1349, 1.0694, 1.0},
183  {0.16553, 0.19066, 0.21468, 0.23609, 0.30416, 0.38821, 0.49644, 0.63386, 0.80299, 0.99907, 1.1304,
184  1.2724, 1.3535, 1.3475, 1.3381, 1.3219, 1.2549, 1.191, 1.1287, 1.0659, 1.0},
185  {0.37736, 0.41414, 0.45135, 0.48843, 0.52473, 0.55973, 0.59348, 0.62696, 0.66202, 0.70042, 0.74241,
186  0.786, 0.82819, 0.86688, 0.90128, 0.93107, 0.95589, 0.97532, 0.98908, 0.99719, 1.0},
187  {0.34354, 0.37692, 0.4109, 0.44492, 0.47873, 0.51296, 0.54937, 0.59047, 0.63799, 0.69117, 0.74652,
188  0.7998, 0.84832, 0.89111, 0.92783, 0.95798, 0.98095, 0.99635, 1.0043, 1.0052, 1.0},
189  {0.36364, 0.39792, 0.43277, 0.4676, 0.50186, 0.53538, 0.56884, 0.604, 0.64308, 0.68729, 0.73544,
190  0.7842, 0.83019, 0.87156, 0.90777, 0.93854, 0.96346, 0.98209, 0.99421, 1, 1.0},
191  {0.36362, 0.39791, 0.43276, 0.46759, 0.50185, 0.53537, 0.56883, 0.604, 0.64307, 0.68728, 0.73544,
192  0.7842, 0.83019, 0.87156, 0.90777, 0.93854, 0.96346, 0.98209, 0.99421, 1, 1.0},
193 };
194 
195 // inelastic interaction length in Silicon at 1 TeV per particle
196 const double nuclInelLength[numHadrons] = {
197  4.5606, 4.4916, 5.7511, 5.7856, 6.797, 6.8373, 6.8171, 6.8171, 0, 0, 4.6926, 4.6926, 4.6926, 4.6926, 0,
198  4.6926, 4.6926, 4.3171, 4.3171, 4.6926, 4.6926, 4.6926, 4.6926, 0, 4.6926, 4.6926, 2.509, 2.9048, 2.5479, 2.5479};
199 
200 // elastic interaction length in Silicon at 1 TeV per particle
201 const double nuclElLength[numHadrons] = {
202  9.248, 9.451, 11.545, 11.671, 32.081, 34.373, 34.373, 32.081, 0, 0, 15.739, 20.348, 15.739, 15.739, 0,
203  20.349, 0, 9.7514, 12.864, 15.836, 20.516, 15.836, 15.744, 0, 20.517, 20.44, 4.129, 6.0904, 4.5204, 4.5204};
204 
205 NuclearInteractionFTFSimulator::NuclearInteractionFTFSimulator(unsigned int distAlgo,
206  double distCut,
207  double elimit,
208  double eth)
209  : curr4Mom(0., 0., 0., 0.),
210  vectProj(0., 0., 1.),
211  theBoost(0., 0., 0.),
212  theBertiniLimit(elimit),
213  theEnergyLimit(eth),
214  theDistCut(distCut),
215  distMin(1E99),
216  theDistAlgo(distAlgo) {
217  // FTF model
218  theHadronicModel = new G4TheoFSGenerator("FTF");
219  theStringModel = new G4FTFModel();
220  G4GeneratorPrecompoundInterface* cascade = new G4GeneratorPrecompoundInterface(new CMSDummyDeexcitation());
221  theLund = new G4LundStringFragmentation();
222  theStringDecay = new G4ExcitedStringDecay(theLund);
223  theStringModel->SetFragmentationModel(theStringDecay);
224 
225  theHadronicModel->SetTransport(cascade);
226  theHadronicModel->SetHighEnergyGenerator(theStringModel);
227  theHadronicModel->SetMinEnergy(theEnergyLimit);
228 
229  // Bertini Cascade
230  theBertiniCascade = new G4CascadeInterface();
231 
232  theDiffuseElastic = new G4DiffuseElastic();
233 
234  // Geant4 particles and cross sections
235  std::call_once(initializeOnce, []() {
236  theG4Hadron[0] = G4Proton::Proton();
237  theG4Hadron[1] = G4Neutron::Neutron();
238  theG4Hadron[2] = G4PionPlus::PionPlus();
239  theG4Hadron[3] = G4PionMinus::PionMinus();
240  theG4Hadron[4] = G4KaonPlus::KaonPlus();
241  theG4Hadron[5] = G4KaonMinus::KaonMinus();
242  theG4Hadron[6] = G4KaonZeroLong::KaonZeroLong();
243  theG4Hadron[7] = G4KaonZeroShort::KaonZeroShort();
244  theG4Hadron[8] = G4KaonZero::KaonZero();
245  theG4Hadron[9] = G4AntiKaonZero::AntiKaonZero();
246  theG4Hadron[10] = G4Lambda::Lambda();
247  theG4Hadron[11] = G4OmegaMinus::OmegaMinus();
248  theG4Hadron[12] = G4SigmaMinus::SigmaMinus();
249  theG4Hadron[13] = G4SigmaPlus::SigmaPlus();
250  theG4Hadron[14] = G4SigmaZero::SigmaZero();
251  theG4Hadron[15] = G4XiMinus::XiMinus();
252  theG4Hadron[16] = G4XiZero::XiZero();
253  theG4Hadron[17] = G4AntiProton::AntiProton();
254  theG4Hadron[18] = G4AntiNeutron::AntiNeutron();
255  theG4Hadron[19] = G4AntiLambda::AntiLambda();
256  theG4Hadron[20] = G4AntiOmegaMinus::AntiOmegaMinus();
257  theG4Hadron[21] = G4AntiSigmaMinus::AntiSigmaMinus();
258  theG4Hadron[22] = G4AntiSigmaPlus::AntiSigmaPlus();
259  theG4Hadron[23] = G4AntiSigmaZero::AntiSigmaZero();
260  theG4Hadron[24] = G4AntiXiMinus::AntiXiMinus();
261  theG4Hadron[25] = G4AntiXiZero::AntiXiZero();
262  theG4Hadron[26] = G4AntiAlpha::AntiAlpha();
263  theG4Hadron[27] = G4AntiDeuteron::AntiDeuteron();
264  theG4Hadron[28] = G4AntiTriton::AntiTriton();
265  theG4Hadron[29] = G4AntiHe3::AntiHe3();
266 
267  // other Geant4 particles
268  G4ParticleDefinition* ion = G4GenericIon::GenericIon();
269  ion->SetProcessManager(new G4ProcessManager(ion));
270  G4DecayPhysics decays;
271  decays.ConstructParticle();
272  G4ParticleTable* partTable = G4ParticleTable::GetParticleTable();
273  partTable->SetVerboseLevel(0);
274  partTable->SetReadiness();
275 
276  for (int i = 0; i < numHadrons; ++i) {
277  theId[i] = theG4Hadron[i]->GetPDGEncoding();
278  }
279  });
280 
281  // local objects
282  vect = new G4PhysicsLogVector(npoints - 1, 100 * MeV, TeV);
283  intLengthElastic = intLengthInelastic = 0.0;
284  currIdx = 0;
285  index = 0;
286  currTrack = nullptr;
287  currParticle = theG4Hadron[0];
288 
289  // fill projectile particle definitions
290  dummyStep = new G4Step();
291  dummyStep->SetPreStepPoint(new G4StepPoint());
292 
293  // target is always Silicon
294  targetNucleus.SetParameters(28, 14);
295 }
296 
297 NuclearInteractionFTFSimulator::~NuclearInteractionFTFSimulator() {
298  delete theStringDecay;
299  delete theStringModel;
300  delete theLund;
301  delete vect;
302 }
303 
304 void NuclearInteractionFTFSimulator::compute(ParticlePropagator& Particle, RandomEngineAndDistribution const* random) {
305  //std::cout << "#### Primary " << Particle.particle().pid() << " E(GeV)= "
306  // << Particle.particle().momentum().e() << std::endl;
307 
308  int thePid = Particle.particle().pid();
309  if (thePid != theId[currIdx]) {
310  currParticle = nullptr;
311  currIdx = 0;
312  for (; currIdx < numHadrons; ++currIdx) {
313  if (theId[currIdx] == thePid) {
314  currParticle = theG4Hadron[currIdx];
315  // neutral kaons
316  if (7 == currIdx || 8 == currIdx) {
317  currParticle = theG4Hadron[9];
318  if (random->flatShoot() > 0.5) {
319  currParticle = theG4Hadron[10];
320  }
321  }
322  break;
323  }
324  }
325  }
326  if (!currParticle) {
327  return;
328  }
329 
330  // fill projectile for Geant4
331  double mass = currParticle->GetPDGMass();
332  double ekin = CLHEP::GeV * Particle.particle().momentum().e() - mass;
333 
334  // check interaction length
335  intLengthElastic = nuclElLength[currIdx];
336  intLengthInelastic = nuclInelLength[currIdx];
337  if (0.0 == intLengthInelastic) {
338  return;
339  }
340 
341  // apply corrections
342  if (ekin <= vect->Energy(0)) {
343  intLengthElastic *= corrfactors_el[currIdx][0];
344  intLengthInelastic *= corrfactors_inel[currIdx][0];
345  } else if (ekin < vect->Energy(npoints - 1)) {
346  index = vect->FindBin(ekin, index);
347  double e1 = vect->Energy(index);
348  double e2 = vect->Energy(index + 1);
349  intLengthElastic *=
350  ((corrfactors_el[currIdx][index] * (e2 - ekin) + corrfactors_el[currIdx][index + 1] * (ekin - e1)) / (e2 - e1));
351  intLengthInelastic *=
352  ((corrfactors_inel[currIdx][index] * (e2 - ekin) + corrfactors_inel[currIdx][index + 1] * (ekin - e1)) /
353  (e2 - e1));
354  }
355  /*
356  std::cout << " Primary " << currParticle->GetParticleName()
357  << " E(GeV)= " << e*fact << std::endl;
358  */
359 
360  double currInteractionLength =
361  -G4Log(random->flatShoot()) * intLengthElastic * intLengthInelastic / (intLengthElastic + intLengthInelastic);
362  /*
363  std::cout << "*NuclearInteractionFTFSimulator::compute: R(X0)= " << radLengths
364  << " Rnuc(X0)= " << theNuclIntLength[currIdx] << " IntLength(X0)= "
365  << currInteractionLength << std::endl;
366  */
367  // Check position of nuclear interaction
368  if (currInteractionLength > radLengths) {
369  return;
370  }
371 
372  // fill projectile for Geant4
373  double px = Particle.particle().momentum().px();
374  double py = Particle.particle().momentum().py();
375  double pz = Particle.particle().momentum().pz();
376  double ptot = sqrt(px * px + py * py + pz * pz);
377  double norm = 1. / ptot;
378  G4ThreeVector dir(px * norm, py * norm, pz * norm);
379  /*
380  std::cout << " Primary " << currParticle->GetParticleName()
381  << " E(GeV)= " << e*fact << " P(GeV/c)= ("
382  << px << " " << py << " " << pz << ")" << std::endl;
383  */
384 
385  G4DynamicParticle* dynParticle = new G4DynamicParticle(theG4Hadron[currIdx], dir, ekin);
386  currTrack = new G4Track(dynParticle, 0.0, vectProj);
387  currTrack->SetStep(dummyStep);
388 
389  theProjectile.Initialise(*currTrack);
390  delete currTrack;
391 
392  G4HadFinalState* result;
393 
394  // elastic interaction
395  if (random->flatShoot() * (intLengthElastic + intLengthInelastic) > intLengthElastic) {
396  result = theDiffuseElastic->ApplyYourself(theProjectile, targetNucleus);
397  G4ThreeVector dirnew = result->GetMomentumChange().unit();
398  double cost = (dir * dirnew);
399  double sint = std::sqrt((1. - cost) * (1. + cost));
400 
401  curr4Mom.set(ptot * dirnew.x(), ptot * dirnew.y(), ptot * dirnew.z(), Particle.particle().momentum().e());
402 
403  // Always create a daughter if the kink is large engough
404  if (sint > theDistCut) {
405  saveDaughter(Particle, curr4Mom, thePid);
406  } else {
407  Particle.particle().setMomentum(curr4Mom.px(), curr4Mom.py(), curr4Mom.pz(), curr4Mom.e());
408  }
409 
410  // inelastic interaction
411  } else {
412  // Bertini cascade for low-energy hadrons (except light anti-nuclei)
413  // FTFP is applied above energy limit and for all anti-hyperons and anti-ions
414  if (ekin <= theBertiniLimit && currIdx < 17) {
415  result = theBertiniCascade->ApplyYourself(theProjectile, targetNucleus);
416  } else {
417  result = theHadronicModel->ApplyYourself(theProjectile, targetNucleus);
418  }
419  if (result) {
420  int nsec = result->GetNumberOfSecondaries();
421  if (0 < nsec) {
422  result->SetTrafoToLab(theProjectile.GetTrafoToLab());
423  _theUpdatedState.clear();
424 
425  //std::cout << " " << nsec << " secondaries" << std::endl;
426  // Generate angle
427  double phi = random->flatShoot() * CLHEP::twopi;
428  theClosestChargedDaughterId = -1;
429  distMin = 1e99;
430 
431  // rotate and store secondaries
432  for (int j = 0; j < nsec; ++j) {
433  const G4DynamicParticle* dp = result->GetSecondary(j)->GetParticle();
434  int thePid = dp->GetParticleDefinition()->GetPDGEncoding();
435 
436  // rotate around primary direction
437  curr4Mom = dp->Get4Momentum();
438  curr4Mom.rotate(phi, vectProj);
439  curr4Mom *= result->GetTrafoToLab();
440  /*
441  std::cout << j << ". " << dp->GetParticleDefinition()->GetParticleName()
442  << " " << thePid
443  << " " << curr4Mom*fact << std::endl;
444  */
445  // prompt 2-gamma decay for pi0, eta, eta'p
446  if (111 == thePid || 221 == thePid || 331 == thePid) {
447  theBoost = curr4Mom.boostVector();
448  double e = 0.5 * dp->GetParticleDefinition()->GetPDGMass();
449  double fi = random->flatShoot() * CLHEP::twopi;
450  double cth = 2 * random->flatShoot() - 1.0;
451  double sth = sqrt((1.0 - cth) * (1.0 + cth));
452  G4LorentzVector lv(e * sth * cos(fi), e * sth * sin(fi), e * cth, e);
453  lv.boost(theBoost);
454  saveDaughter(Particle, lv, 22);
455  curr4Mom -= lv;
456  if (curr4Mom.e() > theEnergyLimit) {
457  saveDaughter(Particle, curr4Mom, 22);
458  }
459  } else {
460  if (curr4Mom.e() > theEnergyLimit + dp->GetParticleDefinition()->GetPDGMass()) {
461  saveDaughter(Particle, curr4Mom, thePid);
462  }
463  }
464  }
465  result->Clear();
466  }
467  }
468  }
469 }
470 
471 void NuclearInteractionFTFSimulator::saveDaughter(ParticlePropagator& Particle, const G4LorentzVector& lv, int pdgid) {
472  unsigned int idx = _theUpdatedState.size();
473  _theUpdatedState.emplace_back(
474  makeParticle(Particle.particleDataTable(),
475  pdgid,
476  XYZTLorentzVector{lv.px() * fact, lv.py() * fact, lv.pz() * fact, lv.e() * fact},
477  Particle.particle().vertex()));
478 
479  // Store the closest daughter index (for later tracking purposes, so charged particles only)
480  double distance = distanceToPrimary(Particle.particle(), _theUpdatedState[idx]);
481  // Find the closest daughter, if closer than a given upper limit.
482  if (distance < distMin && distance < theDistCut) {
483  distMin = distance;
484  theClosestChargedDaughterId = idx;
485  }
486  // std::cout << _theUpdatedState[idx] << std::endl;
487 }
488 
489 double NuclearInteractionFTFSimulator::distanceToPrimary(const RawParticle& Particle,
490  const RawParticle& aDaughter) const {
491  double distance = 2E99;
492  // Compute the distance only for charged primaries
493  if (fabs(Particle.charge()) > 1E-6) {
494  // The secondary must have the same charge
495  double chargeDiff = fabs(aDaughter.charge() - Particle.charge());
496  if (fabs(chargeDiff) < 1E-6) {
497  // Here are two distance definitions * to be tuned *
498  switch (theDistAlgo) {
499  case 1:
500  // sin(theta12)
501  distance = (aDaughter.Vect().Unit().Cross(Particle.Vect().Unit())).R();
502  break;
503 
504  case 2:
505  // sin(theta12) * p1/p2
506  distance = (aDaughter.Vect().Cross(Particle.Vect())).R() / aDaughter.Vect().Mag2();
507  break;
508 
509  default:
510  // Should not happen
511  break;
512  }
513  }
514  }
515  return distance;
516 }
RawParticle::setMomentum
void setMomentum(const XYZTLorentzVector &vtx)
set the momentum
Definition: RawParticle.h:328
nuclInelLength
const double nuclInelLength[numHadrons]
Definition: NuclearInteractionFTFSimulator.cc:196
ParticlePropagator::particleDataTable
const HepPDT::ParticleDataTable * particleDataTable() const
Definition: ParticlePropagator.h:100
RawParticle::Vect
XYZVector Vect() const
the momentum threevector
Definition: RawParticle.h:323
GeV
const double GeV
Definition: MathUtil.h:16
mps_fire.i
i
Definition: mps_fire.py:428
HLT_FULL_cff.distance
tuple distance
Definition: HLT_FULL_cff.py:7955
RandomEngineAndDistribution::flatShoot
double flatShoot(double xmin=0.0, double xmax=1.0) const
Definition: RandomEngineAndDistribution.h:27
NuclearInteractionFTFSimulator::compute
void compute(ParticlePropagator &Particle, RandomEngineAndDistribution const *) override
Generate a nuclear interaction according to the probability that it happens.
Definition: NuclearInteractionFTFSimulator.cc:304
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
thread_safety_macros.h
RawParticle::pid
int pid() const
get the HEP particle ID number
Definition: RawParticle.h:277
BaseParticlePropagator::particle
RawParticle const & particle() const
The particle being propagated.
Definition: BaseParticlePropagator.h:164
corrfactors_inel
const double corrfactors_inel[numHadrons][npoints]
Definition: NuclearInteractionFTFSimulator.cc:77
NuclearInteractionFTFSimulator::~NuclearInteractionFTFSimulator
~NuclearInteractionFTFSimulator() override
Default Destructor.
Definition: NuclearInteractionFTFSimulator.cc:297
mps_fire.result
tuple result
Definition: mps_fire.py:311
makeParticle
RawParticle makeParticle(HepPDT::ParticleDataTable const *, int id, const math::XYZTLorentzVector &p)
Definition: makeParticle.cc:28
MeV
const double MeV
RawParticle::momentum
const XYZTLorentzVector & momentum() const
the momentum fourvector
Definition: RawParticle.h:321
NuclearInteractionFTFSimulator::saveDaughter
void saveDaughter(ParticlePropagator &Particle, const G4LorentzVector &lv, int pdgid)
Definition: NuclearInteractionFTFSimulator.cc:471
DeadROC_duringRun.dir
string dir
Definition: DeadROC_duringRun.py:23
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
initializeOnce
static std::once_flag initializeOnce
Definition: NuclearInteractionFTFSimulator.cc:70
HLT_FULL_cff.distCut
tuple distCut
Definition: HLT_FULL_cff.py:52288
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
CMS_THREAD_GUARD
#define CMS_THREAD_GUARD(_var_)
nuclElLength
const double nuclElLength[numHadrons]
Definition: NuclearInteractionFTFSimulator.cc:201
npoints
static const int npoints
Definition: NuclearInteractionFTFSimulator.h:38
RawParticle::charge
double charge() const
get the MEASURED charge
Definition: RawParticle.h:294
NuclearInteractionFTFSimulator::distanceToPrimary
double distanceToPrimary(const RawParticle &Particle, const RawParticle &aDaughter) const
Definition: NuclearInteractionFTFSimulator.cc:489
numHadrons
static const int numHadrons
Definition: NuclearInteractionFTFSimulator.h:37
corrfactors_el
const double corrfactors_el[numHadrons][npoints]
Definition: NuclearInteractionFTFSimulator.cc:137
NuclearInteractionFTFSimulator::NuclearInteractionFTFSimulator
NuclearInteractionFTFSimulator(unsigned int distAlgo, double distCut, double elimit, double eth)
Constructor.
Definition: NuclearInteractionFTFSimulator.cc:205
NuclearInteractionFTFSimulator::currParticle
const G4ParticleDefinition * currParticle
Definition: NuclearInteractionFTFSimulator.h:71
RawParticle::vertex
const XYZTLorentzVector & vertex() const
the vertex fourvector
Definition: RawParticle.h:320
fact
const double fact
Definition: NuclearInteractionFTFSimulator.cc:74
NuclearInteractionFTFSimulator::theG4Hadron
static const G4ParticleDefinition * theG4Hadron[numHadrons]
Definition: NuclearInteractionFTFSimulator.h:56
MaterialEffectsSimulator::_theUpdatedState
std::vector< RawParticle > _theUpdatedState
Definition: MaterialEffectsSimulator.h:82
XYZTLorentzVector
math::XYZTLorentzVector XYZTLorentzVector
Definition: RawParticle.h:25
Generated for CMSSW Reference Manual by doxygen 1.8.5