CMS 3D CMS Logo
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 }
NuclearInteractionFTFSimulator::theProjectile
G4HadProjectile theProjectile
Definition: NuclearInteractionFTFSimulator.h:74
npoints
static const int npoints
Definition: NuclearInteractionFTFSimulator.h:38
NuclearInteractionFTFSimulator::theEnergyLimit
double theEnergyLimit
Definition: NuclearInteractionFTFSimulator.h:80
RawParticle
Definition: RawParticle.h:37
bJpsiMuMuTrigSettings_cff.decays
decays
Definition: bJpsiMuMuTrigSettings_cff.py:15
mps_fire.i
i
Definition: mps_fire.py:428
NuclearInteractionFTFSimulator::vect
G4PhysicsLogVector * vect
Definition: NuclearInteractionFTFSimulator.h:58
MessageLogger.h
EcalCondDBWriter_cfi.Energy
Energy
Definition: EcalCondDBWriter_cfi.py:152
makeParticle
RawParticle makeParticle(HepPDT::ParticleDataTable const *, int id, const math::XYZTLorentzVector &p)
Definition: makeParticle.cc:28
numHadrons
static const int numHadrons
Definition: NuclearInteractionFTFSimulator.h:37
multPhiCorr_741_25nsDY_cfi.py
py
Definition: multPhiCorr_741_25nsDY_cfi.py:12
MaterialEffectsSimulator::_theUpdatedState
std::vector< RawParticle > _theUpdatedState
Definition: MaterialEffectsSimulator.h:82
nuclInelLength
const double nuclInelLength[numHadrons]
Definition: NuclearInteractionFTFSimulator.cc:196
NuclearInteractionFTFSimulator::targetNucleus
G4Nucleus targetNucleus
Definition: NuclearInteractionFTFSimulator.h:73
MeV
const double MeV
NuclearInteractionFTFSimulator::NuclearInteractionFTFSimulator
NuclearInteractionFTFSimulator(unsigned int distAlgo, double distCut, double elimit, double eth)
Constructor.
Definition: NuclearInteractionFTFSimulator.cc:205
corrfactors_inel
const double corrfactors_inel[numHadrons][npoints]
Definition: NuclearInteractionFTFSimulator.cc:77
XYZTLorentzVector
math::XYZTLorentzVector XYZTLorentzVector
Definition: RawParticle.h:25
MaterialEffectsSimulator::theClosestChargedDaughterId
int theClosestChargedDaughterId
Definition: MaterialEffectsSimulator.h:94
NuclearInteractionFTFSimulator::theBoost
G4ThreeVector theBoost
Definition: NuclearInteractionFTFSimulator.h:77
RawParticle::charge
double charge() const
get the MEASURED charge
Definition: RawParticle.h:294
NuclearInteractionFTFSimulator::dummyStep
G4Step * dummyStep
Definition: NuclearInteractionFTFSimulator.h:69
RandomEngineAndDistribution.h
NuclearInteractionFTFSimulator.h
BPhysicsValidation_cfi.Lambda
Lambda
Definition: BPhysicsValidation_cfi.py:63
heavyIonCSV_trainingSettings.idx
idx
Definition: heavyIonCSV_trainingSettings.py:5
NuclearInteractionFTFSimulator::theStringDecay
G4ExcitedStringDecay * theStringDecay
Definition: NuclearInteractionFTFSimulator.h:62
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
NuclearInteractionFTFSimulator::distanceToPrimary
double distanceToPrimary(const RawParticle &Particle, const RawParticle &aDaughter) const
Definition: NuclearInteractionFTFSimulator.cc:489
NuclearInteractionFTFSimulator::theBertiniCascade
G4CascadeInterface * theBertiniCascade
Definition: NuclearInteractionFTFSimulator.h:66
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
phase2TrackerDigitizer_cfi.SigmaZero
SigmaZero
Definition: phase2TrackerDigitizer_cfi.py:28
Calorimetry_cff.dp
dp
Definition: Calorimetry_cff.py:157
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
NuclearInteractionFTFSimulator::theG4Hadron
static const G4ParticleDefinition * theG4Hadron[numHadrons]
Definition: NuclearInteractionFTFSimulator.h:56
NuclearInteractionFTFSimulator::theDistAlgo
unsigned int theDistAlgo
Definition: NuclearInteractionFTFSimulator.h:90
ParticlePropagator.h
RawParticle::Vect
XYZVector Vect() const
the momentum threevector
Definition: RawParticle.h:323
corrfactors_el
const double corrfactors_el[numHadrons][npoints]
Definition: NuclearInteractionFTFSimulator.cc:137
NuclearInteractionFTFSimulator::currIdx
int currIdx
Definition: NuclearInteractionFTFSimulator.h:88
MaterialEffectsSimulator::radLengths
double radLengths
Definition: MaterialEffectsSimulator.h:84
GeV
const double GeV
Definition: MathUtil.h:16
thread_safety_macros.h
StorageManager_cfg.e1
e1
Definition: StorageManager_cfg.py:16
fact
const double fact
Definition: NuclearInteractionFTFSimulator.cc:74
NuclearInteractionFTFSimulator::theDiffuseElastic
G4DiffuseElastic * theDiffuseElastic
Definition: NuclearInteractionFTFSimulator.h:67
ParticlePropagator
Definition: ParticlePropagator.h:28
NuclearInteractionFTFSimulator::compute
void compute(ParticlePropagator &Particle, RandomEngineAndDistribution const *) override
Generate a nuclear interaction according to the probability that it happens.
Definition: NuclearInteractionFTFSimulator.cc:304
initializeOnce
static std::once_flag initializeOnce
Definition: NuclearInteractionFTFSimulator.cc:70
NuclearInteractionFTFSimulator::intLengthElastic
double intLengthElastic
Definition: NuclearInteractionFTFSimulator.h:85
NuclearInteractionFTFSimulator::theBertiniLimit
double theBertiniLimit
Definition: NuclearInteractionFTFSimulator.h:79
makeParticle.h
NuclearInteractionFTFSimulator::theStringModel
G4FTFModel * theStringModel
Definition: NuclearInteractionFTFSimulator.h:61
NuclearInteractionFTFSimulator::currTrack
G4Track * currTrack
Definition: NuclearInteractionFTFSimulator.h:70
NuclearInteractionFTFSimulator::theId
static int theId[numHadrons]
Definition: NuclearInteractionFTFSimulator.h:57
RandomEngineAndDistribution::flatShoot
double flatShoot(double xmin=0.0, double xmax=1.0) const
Definition: RandomEngineAndDistribution.h:27
NuclearInteractionFTFSimulator::index
size_t index
Definition: NuclearInteractionFTFSimulator.h:89
NuclearInteractionFTFSimulator::distMin
double distMin
Definition: NuclearInteractionFTFSimulator.h:83
nuclElLength
const double nuclElLength[numHadrons]
Definition: NuclearInteractionFTFSimulator.cc:201
DDAxes::phi
NuclearInteractionFTFSimulator::vectProj
G4ThreeVector vectProj
Definition: NuclearInteractionFTFSimulator.h:76
multPhiCorr_741_25nsDY_cfi.px
px
Definition: multPhiCorr_741_25nsDY_cfi.py:10
NuclearInteractionFTFSimulator::currParticle
const G4ParticleDefinition * currParticle
Definition: NuclearInteractionFTFSimulator.h:71
NuclearInteractionFTFSimulator::~NuclearInteractionFTFSimulator
~NuclearInteractionFTFSimulator() override
Default Destructor.
Definition: NuclearInteractionFTFSimulator.cc:297
muonTagProbeFilters_cff.distMin
distMin
Definition: muonTagProbeFilters_cff.py:61
NuclearInteractionFTFSimulator::theHadronicModel
G4TheoFSGenerator * theHadronicModel
Definition: NuclearInteractionFTFSimulator.h:60
NuclearInteractionFTFSimulator::theLund
G4LundStringFragmentation * theLund
Definition: NuclearInteractionFTFSimulator.h:63
MaterialEffects_cfi.distCut
distCut
Definition: MaterialEffects_cfi.py:128
EgHLTOffHistBins_cfi.mass
mass
Definition: EgHLTOffHistBins_cfi.py:34
Particle
Definition: Particle.py:1
mps_fire.result
result
Definition: mps_fire.py:311
MaterialEffects_cfi.distAlgo
distAlgo
Definition: MaterialEffects_cfi.py:127
dqmiolumiharvest.j
j
Definition: dqmiolumiharvest.py:66
EgammaValidation_cff.pdgid
pdgid
Definition: EgammaValidation_cff.py:30
HLT_FULL_cff.distance
distance
Definition: HLT_FULL_cff.py:7733
CMS_THREAD_GUARD
#define CMS_THREAD_GUARD(_var_)
Definition: thread_safety_macros.h:6
CMSDummyDeexcitation
Definition: CMSDummyDeexcitation.h:20
NuclearInteractionFTFSimulator::saveDaughter
void saveDaughter(ParticlePropagator &Particle, const G4LorentzVector &lv, int pdgid)
Definition: NuclearInteractionFTFSimulator.cc:471
CMSDummyDeexcitation.h
NuclearInteractionFTFSimulator::theDistCut
double theDistCut
Definition: NuclearInteractionFTFSimulator.h:82
NuclearInteractionFTFSimulator::intLengthInelastic
double intLengthInelastic
Definition: NuclearInteractionFTFSimulator.h:86
DeadROC_duringRun.dir
dir
Definition: DeadROC_duringRun.py:23
NuclearInteractionFTFSimulator::curr4Mom
G4LorentzVector curr4Mom
Definition: NuclearInteractionFTFSimulator.h:75
MillePedeFileConverter_cfg.e
e
Definition: MillePedeFileConverter_cfg.py:37
RandomEngineAndDistribution
Definition: RandomEngineAndDistribution.h:18
Generated for CMSSW Reference Manual by doxygen 1.8.16