MaterialEffects/src/NuclearInteractionFTFSimulator.cc

0001 // system headers
0002 #include <mutex>
0003
0004 // Framework Headers
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006 #include "FWCore/Utilities/interface/thread_safety_macros.h"
0007
0008 // Fast Sim headers
0009 #include "FastSimulation/MaterialEffects/interface/NuclearInteractionFTFSimulator.h"
0010 #include "FastSimulation/MaterialEffects/interface/CMSDummyDeexcitation.h"
0011 #include "FastSimulation/Utilities/interface/RandomEngineAndDistribution.h"
0012 #include "FastSimulation/ParticlePropagator/interface/ParticlePropagator.h"
0013 #include "FastSimulation/Particle/interface/makeParticle.h"
0014
0015 // Geant4 headers
0016 #include "G4ParticleDefinition.hh"
0017 #include "G4DynamicParticle.hh"
0018 #include "G4Track.hh"
0019 #include "G4Step.hh"
0020 #include "G4StepPoint.hh"
0021 #include "G4TheoFSGenerator.hh"
0022 #include "G4FTFModel.hh"
0023 #include "G4ExcitedStringDecay.hh"
0024 #include "G4LundStringFragmentation.hh"
0025 #include "G4GeneratorPrecompoundInterface.hh"
0026 #include "G4CascadeInterface.hh"
0027 #include "G4DiffuseElastic.hh"
0028
0029 #include "G4Proton.hh"
0030 #include "G4Neutron.hh"
0031 #include "G4PionPlus.hh"
0032 #include "G4PionMinus.hh"
0033 #include "G4AntiProton.hh"
0034 #include "G4AntiNeutron.hh"
0035 #include "G4KaonPlus.hh"
0036 #include "G4KaonMinus.hh"
0037 #include "G4KaonZeroLong.hh"
0038 #include "G4KaonZeroShort.hh"
0039 #include "G4KaonZero.hh"
0040 #include "G4AntiKaonZero.hh"
0041 #include "G4GenericIon.hh"
0042
0043 #include "G4Lambda.hh"
0044 #include "G4OmegaMinus.hh"
0045 #include "G4SigmaMinus.hh"
0046 #include "G4SigmaPlus.hh"
0047 #include "G4SigmaZero.hh"
0048 #include "G4XiMinus.hh"
0049 #include "G4XiZero.hh"
0050 #include "G4AntiLambda.hh"
0051 #include "G4AntiOmegaMinus.hh"
0052 #include "G4AntiSigmaMinus.hh"
0053 #include "G4AntiSigmaPlus.hh"
0054 #include "G4AntiSigmaZero.hh"
0055 #include "G4AntiXiMinus.hh"
0056 #include "G4AntiXiZero.hh"
0057 #include "G4AntiAlpha.hh"
0058 #include "G4AntiDeuteron.hh"
0059 #include "G4AntiTriton.hh"
0060 #include "G4AntiHe3.hh"
0061
0062 #include "G4Material.hh"
0063 #include "G4DecayPhysics.hh"
0064 #include "G4ParticleTable.hh"
0065 #include "G4IonTable.hh"
0066 #include "G4ProcessManager.hh"
0067 #include "G4PhysicsLogVector.hh"
0068 #include <CLHEP/Units/SystemOfUnits.h>
0069
0070 static std::once_flag initializeOnce;
0071 CMS_THREAD_GUARD(initializeOnce) const G4ParticleDefinition* NuclearInteractionFTFSimulator::theG4Hadron[] = {nullptr};
0072 CMS_THREAD_GUARD(initializeOnce) int NuclearInteractionFTFSimulator::theId[] = {0};
0073
0074 const double fact = 1.0 / CLHEP::GeV;
0075
0076 // inelastic interaction length corrections per particle and energy
0077 const double corrfactors_inel[numHadrons][npoints] = {
0078     {1.0872, 1.1026, 1.111,  1.111,  1.0105, 0.97622, 0.9511, 0.9526, 0.97591, 0.99277, 1.0099,
0079      1.015,  1.0217, 1.0305, 1.0391, 1.0438, 1.0397,  1.0328, 1.0232, 1.0123,  1.0},
0080     {1.0416, 1.1044, 1.1467, 1.1273, 1.026,  0.99085, 0.96572, 0.96724, 0.99091, 1.008, 1.0247,
0081      1.0306, 1.0378, 1.0427, 1.0448, 1.0438, 1.0397,  1.0328,  1.0232,  1.0123,  1.0},
0082     {0.5308, 0.53589, 0.67059, 0.80253, 0.82341, 0.79083, 0.85967, 0.90248, 0.93792, 0.9673, 1.0034,
0083      1.022,  1.0418,  1.0596,  1.0749,  1.079,   1.0704,  1.0576,  1.0408,  1.0214,  1.0},
0084     {0.49107, 0.50571, 0.64149, 0.77209, 0.80472, 0.78166, 0.83509, 0.8971, 0.93234, 0.96154, 0.99744,
0085      1.0159,  1.0355,  1.0533,  1.0685,  1.0732,  1.0675,  1.0485,  1.0355, 1.0191,  1.0},
0086     {1.9746, 1.7887, 1.5645, 1.2817, 1.0187, 0.95216, 0.9998, 1.035,  1.0498, 1.0535, 1.0524,
0087      1.0495, 1.0461, 1.0424, 1.0383, 1.0338, 1.0287,  1.0228, 1.0161, 1.0085, 1.0},
0088     {0.46028, 0.59514, 0.70355, 0.70698, 0.62461, 0.65103, 0.71945, 0.77753, 0.83582, 0.88422, 0.92117,
0089      0.94889, 0.96963, 0.98497, 0.99596, 1.0033,  1.0075,  1.0091,  1.0081,  1.005,   1.0},
0090     {0.75016, 0.89607, 0.97185, 0.91083, 0.77425, 0.77412, 0.8374, 0.88848, 0.93104, 0.96174, 0.98262,
0091      0.99684, 1.0065,  1.0129,  1.0168,  1.0184,  1.018,   1.0159, 1.0121,  1.0068,  1.0},
0092     {0.75016, 0.89607, 0.97185, 0.91083, 0.77425, 0.77412, 0.8374, 0.88848, 0.93104, 0.96174, 0.98262,
0093      0.99684, 1.0065,  1.0129,  1.0168,  1.0184,  1.018,   1.0159, 1.0121,  1.0068,  1.0},
0094     {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.0},
0095     {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.0},
0096     {1.1006, 1.1332, 1.121,  1.1008, 1.086,  1.077, 1.0717, 1.0679, 1.0643, 1.0608, 1.057,
0097      1.053,  1.0487, 1.0441, 1.0392, 1.0338, 1.028, 1.0218, 1.015,  1.0078, 1.0},
0098     {1.1318, 1.1255, 1.1062, 1.0904, 1.0802, 1.0742, 1.0701, 1.0668, 1.0636, 1.0602, 1.0566,
0099      1.0527, 1.0485, 1.044,  1.0391, 1.0337, 1.028,  1.0217, 1.015,  1.0078, 1.0},
0100     {1.1094, 1.1332, 1.1184, 1.0988, 1.0848, 1.0765, 1.0714, 1.0677, 1.0642, 1.0607, 1.0569,
0101      1.053,  1.0487, 1.0441, 1.0391, 1.0338, 1.028,  1.0218, 1.015,  1.0078, 1.0},
0102     {1.1087, 1.1332, 1.1187, 1.099,  1.0849, 1.0765, 1.0715, 1.0677, 1.0642, 1.0607, 1.057,
0103      1.053,  1.0487, 1.0441, 1.0391, 1.0338, 1.028,  1.0218, 1.015,  1.0078, 1.0},
0104     {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.0},
0105     {1.1192, 1.132,  1.1147, 1.0961, 1.0834, 1.0758, 1.0711, 1.0674, 1.064,  1.0606, 1.0569,
0106      1.0529, 1.0486, 1.0441, 1.0391, 1.0338, 1.028,  1.0218, 1.015,  1.0078, 1.0},
0107     {1.1188, 1.1321, 1.1149, 1.0963, 1.0834, 1.0758, 1.0711, 1.0675, 1.0641, 1.0606, 1.0569,
0108      1.0529, 1.0486, 1.0441, 1.0391, 1.0338, 1.028,  1.0218, 1.015,  1.0078, 1.0},
0109     {0.50776, 0.5463,  0.5833,  0.61873, 0.65355, 0.68954, 0.72837, 0.7701, 0.81267, 0.85332, 0.89037,
0110      0.92329, 0.95177, 0.97539, 0.99373, 1.0066,  1.014,   1.0164,  1.0144, 1.0087,  1.0},
0111     {0.50787, 0.5464,  0.58338, 0.6188,  0.65361, 0.6896, 0.72841, 0.77013, 0.8127, 0.85333, 0.89038,
0112      0.92329, 0.95178, 0.9754,  0.99373, 1.0066,  1.014,  1.0164,  1.0144,  1.0087, 1.0},
0113     {1.1006, 1.1332, 1.121,  1.1008, 1.086,  1.077, 1.0717, 1.0679, 1.0643, 1.0608, 1.057,
0114      1.053,  1.0487, 1.0441, 1.0392, 1.0338, 1.028, 1.0218, 1.015,  1.0078, 1.0},
0115     {1.1318, 1.1255, 1.1062, 1.0904, 1.0802, 1.0742, 1.0701, 1.0668, 1.0636, 1.0602, 1.0566,
0116      1.0527, 1.0485, 1.044,  1.0391, 1.0337, 1.028,  1.0217, 1.015,  1.0078, 1.0},
0117     {1.1094, 1.1332, 1.1184, 1.0988, 1.0848, 1.0765, 1.0714, 1.0677, 1.0642, 1.0607, 1.0569,
0118      1.053,  1.0487, 1.0441, 1.0391, 1.0338, 1.028,  1.0218, 1.015,  1.0078, 1.0},
0119     {1.1087, 1.1332, 1.1187, 1.099,  1.0849, 1.0765, 1.0715, 1.0677, 1.0642, 1.0607, 1.057,
0120      1.053,  1.0487, 1.0441, 1.0391, 1.0338, 1.028,  1.0218, 1.015,  1.0078, 1.0},
0121     {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.0},
0122     {1.1192, 1.132,  1.1147, 1.0961, 1.0834, 1.0758, 1.0711, 1.0674, 1.064,  1.0606, 1.0569,
0123      1.0529, 1.0486, 1.0441, 1.0391, 1.0338, 1.028,  1.0218, 1.015,  1.0078, 1.0},
0124     {1.1188, 1.1321, 1.1149, 1.0963, 1.0834, 1.0758, 1.0711, 1.0675, 1.0641, 1.0606, 1.0569,
0125      1.0529, 1.0486, 1.0441, 1.0391, 1.0338, 1.028,  1.0218, 1.015,  1.0078, 1.0},
0126     {0.47677, 0.51941, 0.56129, 0.60176, 0.64014, 0.67589, 0.70891, 0.73991, 0.77025, 0.80104, 0.83222,
0127      0.86236, 0.8901,  0.91518, 0.9377,  0.95733, 0.97351, 0.98584, 0.9942,  0.99879, 1.0},
0128     {0.49361, 0.53221, 0.56976, 0.60563, 0.63954, 0.67193, 0.70411, 0.73777, 0.77378, 0.81114, 0.84754,
0129      0.88109, 0.91113, 0.93745, 0.95974, 0.97762, 0.99081, 0.99929, 1.0033,  1.0034,  1.0},
0130     {0.4873,  0.52744, 0.56669, 0.60443, 0.64007, 0.67337, 0.70482, 0.73572, 0.76755, 0.80086, 0.83456,
0131      0.86665, 0.8959,  0.92208, 0.94503, 0.96437, 0.97967, 0.99072, 0.99756, 1.0005,  1.0},
0132     {0.48729, 0.52742, 0.56668, 0.60442, 0.64006, 0.67336, 0.70482, 0.73571, 0.76754, 0.80086, 0.83455,
0133      0.86665, 0.8959,  0.92208, 0.94503, 0.96437, 0.97967, 0.99072, 0.99756, 1.0005,  1.0},
0134 };
0135
0136 // elastic interaction length corrections per particle and energy
0137 const double corrfactors_el[numHadrons][npoints] = {
0138     {0.58834, 1.1238, 1.7896, 1.4409, 0.93175, 0.80587, 0.80937, 0.83954, 0.87453, 0.91082, 0.94713,
0139      0.98195, 1.0134, 1.0397, 1.0593, 1.071,   1.0739,  1.0678,  1.053,   1.03,    1.0},
0140     {0.40938, 0.92337, 1.3365, 1.1607, 1.008,  0.82206, 0.81163, 0.79489, 0.82919, 0.91812, 0.96688,
0141      1.0225,  1.0734,  1.0833, 1.0874, 1.0854, 1.0773,  1.0637,  1.0448,  1.0235,  1.0},
0142     {0.43699, 0.42165, 0.46594, 0.64917, 0.85314, 0.80782, 0.83204, 0.91162, 1.0155, 1.0665, 1.0967,
0143      1.1125,  1.1275,  1.1376,  1.1464,  1.1477,  1.1312,  1.1067,  1.0751,  1.039,  1.0},
0144     {0.3888, 0.39527, 0.43921, 0.62834, 0.8164, 0.79866, 0.82272, 0.90163, 1.0045, 1.055, 1.0849,
0145      1.1005, 1.1153,  1.1253,  1.134,   1.1365, 1.1255,  1.0895,  1.0652,  1.0348, 1.0},
0146     {0.32004, 0.31119, 0.30453, 0.30004, 0.31954, 0.40148, 0.5481, 0.74485, 0.99317, 1.1642, 1.2117,
0147      1.2351,  1.2649,  1.3054,  1.375,   1.4992,  1.4098,  1.3191, 1.2232,  1.118,   1.0},
0148     {0.10553, 0.14623, 0.20655, 0.26279, 0.19996, 0.40125, 0.5139, 0.71271, 0.89269, 1.0108, 1.1673,
0149      1.3052,  1.4149,  1.429,   1.4521,  1.4886,  1.4006,  1.3116, 1.2177,  1.1151,  1.0},
0150     {0.106,  0.14692, 0.20755, 0.26257, 0.20089, 0.40236, 0.51452, 0.71316, 0.89295, 1.0109, 1.1673,
0151      1.3053, 1.4149,  1.429,   1.4521,  1.4886,  1.4006,  1.3116,  1.2177,  1.1151,  1.0},
0152     {0.31991, 0.31111, 0.30445, 0.30004, 0.31995, 0.40221, 0.54884, 0.74534, 0.99364, 1.1644, 1.2117,
0153      1.2351,  1.265,   1.3054,  1.375,   1.4992,  1.4098,  1.3191,  1.2232,  1.118,   1.0},
0154     {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.0},
0155     {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.0},
0156     {0.37579, 0.39922, 0.37445, 0.32631, 0.39002, 0.42161, 0.54251, 0.69127, 0.90332, 1.0664, 1.1346,
0157      1.1481,  1.1692,  1.2036,  1.2625,  1.3633,  1.2913,  1.2215,  1.1516,  1.0788,  1.0},
0158     {0.31756, 0.33409, 0.25339, 0.35525, 0.52989, 0.63382, 0.7453, 0.93505, 1.1464, 1.2942, 1.3161,
0159      1.328,   1.3393,  1.3525,  1.374,   1.4051,  1.3282,  1.2523, 1.1745,  1.0916, 1.0},
0160     {0.38204, 0.39694, 0.36502, 0.33367, 0.39229, 0.43119, 0.54898, 0.70169, 0.91004, 1.0696, 1.1348,
0161      1.1483,  1.1694,  1.2038,  1.2627,  1.3632,  1.2913,  1.2215,  1.1516,  1.0788,  1.0},
0162     {0.38143, 0.39716, 0.36609, 0.33294, 0.39207, 0.43021, 0.54834, 0.70066, 0.90945, 1.0693, 1.1348,
0163      1.1482,  1.1694,  1.2038,  1.2627,  1.3632,  1.2913,  1.2215,  1.1516,  1.0788,  1.0},
0164     {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.0},
0165     {0.29564, 0.32645, 0.29986, 0.30611, 0.48808, 0.59902, 0.71207, 0.8832, 1.1164, 1.2817, 1.3154,
0166      1.3273,  1.3389,  1.3521,  1.3736,  1.4056,  1.3285,  1.2524,  1.1746, 1.0916, 1.0},
0167     {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.0},
0168     {0.3265,  0.3591,  0.39232, 0.42635, 0.46259, 0.50365, 0.55244, 0.61014, 0.67446, 0.74026, 0.80252,
0169      0.85858, 0.90765, 0.94928, 0.9827,  1.0071,  1.0221,  1.0279,  1.0253,  1.0155,  1.0},
0170     {0.13808, 0.15585, 0.17798, 0.2045, 0.22596, 0.25427, 0.33214, 0.44821, 0.5856, 0.74959, 0.89334,
0171      1.0081,  1.0964,  1.1248,  1.173,  1.2548,  1.1952,  1.1406,  1.0903,  1.0437, 1.0},
0172     {0.20585, 0.23253, 0.26371, 0.28117, 0.30433, 0.35417, 0.44902, 0.58211, 0.73486, 0.90579, 1.0395,
0173      1.1488,  1.2211,  1.2341,  1.2553,  1.2877,  1.2245,  1.1654,  1.1093,  1.0547,  1.0},
0174     {0.2852, 0.32363, 0.31419, 0.35164, 0.45463, 0.54331, 0.66908, 0.81735, 0.98253, 1.1557, 1.2557,
0175      1.3702, 1.4186,  1.401,   1.374,   1.3325,  1.2644,  1.1991,  1.1348,  1.0694,  1.0},
0176     {0.20928, 0.23671, 0.2664, 0.28392, 0.30584, 0.35929, 0.45725, 0.5893, 0.74047, 0.9101, 1.0407,
0177      1.1503,  1.2212,  1.2342, 1.2554,  1.2876,  1.2245,  1.1654,  1.1093, 1.0547,  1.0},
0178     {0.11897, 0.13611, 0.15796, 0.1797, 0.21335, 0.26367, 0.34705, 0.46115, 0.6016, 0.7759, 0.91619,
0179      1.0523,  1.1484,  1.1714,  1.2098, 1.2721,  1.2106,  1.1537,  1.1004,  1.0496, 1.0},
0180     {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.0},
0181     {0.26663, 0.30469, 0.32886, 0.33487, 0.41692, 0.51616, 0.63323, 0.78162, 0.95551, 1.1372, 1.2502,
0182      1.3634,  1.4189,  1.4013,  1.3743,  1.3329,  1.2646,  1.1992,  1.1349,  1.0694,  1.0},
0183     {0.16553, 0.19066, 0.21468, 0.23609, 0.30416, 0.38821, 0.49644, 0.63386, 0.80299, 0.99907, 1.1304,
0184      1.2724,  1.3535,  1.3475,  1.3381,  1.3219,  1.2549,  1.191,   1.1287,  1.0659,  1.0},
0185     {0.37736, 0.41414, 0.45135, 0.48843, 0.52473, 0.55973, 0.59348, 0.62696, 0.66202, 0.70042, 0.74241,
0186      0.786,   0.82819, 0.86688, 0.90128, 0.93107, 0.95589, 0.97532, 0.98908, 0.99719, 1.0},
0187     {0.34354, 0.37692, 0.4109,  0.44492, 0.47873, 0.51296, 0.54937, 0.59047, 0.63799, 0.69117, 0.74652,
0188      0.7998,  0.84832, 0.89111, 0.92783, 0.95798, 0.98095, 0.99635, 1.0043,  1.0052,  1.0},
0189     {0.36364, 0.39792, 0.43277, 0.4676,  0.50186, 0.53538, 0.56884, 0.604,   0.64308, 0.68729, 0.73544,
0190      0.7842,  0.83019, 0.87156, 0.90777, 0.93854, 0.96346, 0.98209, 0.99421, 1,       1.0},
0191     {0.36362, 0.39791, 0.43276, 0.46759, 0.50185, 0.53537, 0.56883, 0.604,   0.64307, 0.68728, 0.73544,
0192      0.7842,  0.83019, 0.87156, 0.90777, 0.93854, 0.96346, 0.98209, 0.99421, 1,       1.0},
0193 };
0194
0195 // inelastic interaction length in Silicon at 1 TeV per particle
0196 const double nuclInelLength[numHadrons] = {
0197     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,
0198     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};
0199
0200 // elastic interaction length in Silicon at 1 TeV per particle
0201 const double nuclElLength[numHadrons] = {
0202     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,
0203     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};
0204
0205 NuclearInteractionFTFSimulator::NuclearInteractionFTFSimulator(unsigned int distAlgo,
0206                                                                double distCut,
0207                                                                double elimit,
0208                                                                double eth)
0209     : curr4Mom(0., 0., 0., 0.),
0210       vectProj(0., 0., 1.),
0211       theBoost(0., 0., 0.),
0212       theBertiniLimit(elimit),
0213       theEnergyLimit(eth),
0214       theDistCut(distCut),
0215       distMin(1E99),
0216       theDistAlgo(distAlgo) {
0217   // FTF model
0218   theHadronicModel = new G4TheoFSGenerator("FTF");
0219   theStringModel = new G4FTFModel();
0220   G4GeneratorPrecompoundInterface* cascade = new G4GeneratorPrecompoundInterface(new CMSDummyDeexcitation());
0221   theLund = new G4LundStringFragmentation();
0222   theStringDecay = new G4ExcitedStringDecay(theLund);
0223   theStringModel->SetFragmentationModel(theStringDecay);
0224
0225   theHadronicModel->SetTransport(cascade);
0226   theHadronicModel->SetHighEnergyGenerator(theStringModel);
0227   theHadronicModel->SetMinEnergy(theEnergyLimit);
0228
0229   // Bertini Cascade
0230   theBertiniCascade = new G4CascadeInterface();
0231
0232   theDiffuseElastic = new G4DiffuseElastic();
0233
0234   // Geant4 particles and cross sections
0235   std::call_once(initializeOnce, []() {
0236     theG4Hadron[0] = G4Proton::Proton();
0237     theG4Hadron[1] = G4Neutron::Neutron();
0238     theG4Hadron[2] = G4PionPlus::PionPlus();
0239     theG4Hadron[3] = G4PionMinus::PionMinus();
0240     theG4Hadron[4] = G4KaonPlus::KaonPlus();
0241     theG4Hadron[5] = G4KaonMinus::KaonMinus();
0242     theG4Hadron[6] = G4KaonZeroLong::KaonZeroLong();
0243     theG4Hadron[7] = G4KaonZeroShort::KaonZeroShort();
0244     theG4Hadron[8] = G4KaonZero::KaonZero();
0245     theG4Hadron[9] = G4AntiKaonZero::AntiKaonZero();
0246     theG4Hadron[10] = G4Lambda::Lambda();
0247     theG4Hadron[11] = G4OmegaMinus::OmegaMinus();
0248     theG4Hadron[12] = G4SigmaMinus::SigmaMinus();
0249     theG4Hadron[13] = G4SigmaPlus::SigmaPlus();
0250     theG4Hadron[14] = G4SigmaZero::SigmaZero();
0251     theG4Hadron[15] = G4XiMinus::XiMinus();
0252     theG4Hadron[16] = G4XiZero::XiZero();
0253     theG4Hadron[17] = G4AntiProton::AntiProton();
0254     theG4Hadron[18] = G4AntiNeutron::AntiNeutron();
0255     theG4Hadron[19] = G4AntiLambda::AntiLambda();
0256     theG4Hadron[20] = G4AntiOmegaMinus::AntiOmegaMinus();
0257     theG4Hadron[21] = G4AntiSigmaMinus::AntiSigmaMinus();
0258     theG4Hadron[22] = G4AntiSigmaPlus::AntiSigmaPlus();
0259     theG4Hadron[23] = G4AntiSigmaZero::AntiSigmaZero();
0260     theG4Hadron[24] = G4AntiXiMinus::AntiXiMinus();
0261     theG4Hadron[25] = G4AntiXiZero::AntiXiZero();
0262     theG4Hadron[26] = G4AntiAlpha::AntiAlpha();
0263     theG4Hadron[27] = G4AntiDeuteron::AntiDeuteron();
0264     theG4Hadron[28] = G4AntiTriton::AntiTriton();
0265     theG4Hadron[29] = G4AntiHe3::AntiHe3();
0266
0267     // other Geant4 particles
0268     G4ParticleDefinition* ion = G4GenericIon::GenericIon();
0269     ion->SetProcessManager(new G4ProcessManager(ion));
0270     G4DecayPhysics decays;
0271     decays.ConstructParticle();
0272     G4ParticleTable* partTable = G4ParticleTable::GetParticleTable();
0273     partTable->SetVerboseLevel(0);
0274     partTable->SetReadiness();
0275
0276     for (int i = 0; i < numHadrons; ++i) {
0277       theId[i] = theG4Hadron[i]->GetPDGEncoding();
0278     }
0279   });
0280
0281   // local objects
0282   vect = new G4PhysicsLogVector(npoints - 1, 100 * CLHEP::MeV, CLHEP::TeV);
0283   intLengthElastic = intLengthInelastic = 0.0;
0284   currIdx = 0;
0285   index = 0;
0286   currTrack = nullptr;
0287   currParticle = theG4Hadron[0];
0288
0289   // fill projectile particle definitions
0290   dummyStep = new G4Step();
0291   dummyStep->SetPreStepPoint(new G4StepPoint());
0292
0293   // target is always Silicon
0294   targetNucleus.SetParameters(28, 14);
0295 }
0296
0297 NuclearInteractionFTFSimulator::~NuclearInteractionFTFSimulator() {
0298   delete theStringDecay;
0299   delete theStringModel;
0300   delete theLund;
0301   delete vect;
0302 }
0303
0304 void NuclearInteractionFTFSimulator::compute(ParticlePropagator& Particle, RandomEngineAndDistribution const* random) {
0305   //std::cout << "#### Primary " << Particle.particle().pid() << " E(GeV)= "
0306   //        << Particle.particle().momentum().e() << std::endl;
0307
0308   int thePid = Particle.particle().pid();
0309   if (thePid != theId[currIdx]) {
0310     currParticle = nullptr;
0311     currIdx = 0;
0312     for (; currIdx < numHadrons; ++currIdx) {
0313       if (theId[currIdx] == thePid) {
0314         currParticle = theG4Hadron[currIdx];
0315         // neutral kaons
0316         if (7 == currIdx || 8 == currIdx) {
0317           currParticle = theG4Hadron[9];
0318           if (random->flatShoot() > 0.5) {
0319             currParticle = theG4Hadron[10];
0320           }
0321         }
0322         break;
0323       }
0324     }
0325   }
0326   if (!currParticle) {
0327     return;
0328   }
0329
0330   // fill projectile for Geant4
0331   double mass = currParticle->GetPDGMass();
0332   double ekin = CLHEP::GeV * Particle.particle().momentum().e() - mass;
0333
0334   // check interaction length
0335   intLengthElastic = nuclElLength[currIdx];
0336   intLengthInelastic = nuclInelLength[currIdx];
0337   if (0.0 == intLengthInelastic) {
0338     return;
0339   }
0340
0341   // apply corrections
0342   if (ekin <= vect->Energy(0)) {
0343     intLengthElastic *= corrfactors_el[currIdx][0];
0344     intLengthInelastic *= corrfactors_inel[currIdx][0];
0345   } else if (ekin < vect->Energy(npoints - 1)) {
0346     index = vect->FindBin(ekin, index);
0347     double e1 = vect->Energy(index);
0348     double e2 = vect->Energy(index + 1);
0349     intLengthElastic *=
0350         ((corrfactors_el[currIdx][index] * (e2 - ekin) + corrfactors_el[currIdx][index + 1] * (ekin - e1)) / (e2 - e1));
0351     intLengthInelastic *=
0352         ((corrfactors_inel[currIdx][index] * (e2 - ekin) + corrfactors_inel[currIdx][index + 1] * (ekin - e1)) /
0353          (e2 - e1));
0354   }
0355   /*
0356   std::cout << " Primary " <<  currParticle->GetParticleName()
0357         << "  E(GeV)= " << e*fact << std::endl;
0358   */
0359
0360   double currInteractionLength =
0361       -G4Log(random->flatShoot()) * intLengthElastic * intLengthInelastic / (intLengthElastic + intLengthInelastic);
0362   /*
0363   std::cout << "*NuclearInteractionFTFSimulator::compute: R(X0)= " << radLengths
0364         << " Rnuc(X0)= " << theNuclIntLength[currIdx] << "  IntLength(X0)= "
0365             << currInteractionLength << std::endl;
0366   */
0367   // Check position of nuclear interaction
0368   if (currInteractionLength > radLengths) {
0369     return;
0370   }
0371
0372   // fill projectile for Geant4
0373   double px = Particle.particle().momentum().px();
0374   double py = Particle.particle().momentum().py();
0375   double pz = Particle.particle().momentum().pz();
0376   double ptot = sqrt(px * px + py * py + pz * pz);
0377   double norm = 1. / ptot;
0378   G4ThreeVector dir(px * norm, py * norm, pz * norm);
0379   /*
0380   std::cout << " Primary " <<  currParticle->GetParticleName()
0381         << "  E(GeV)= " << e*fact << "  P(GeV/c)= ("
0382         << px << " " << py << " " << pz << ")" << std::endl;
0383   */
0384
0385   G4DynamicParticle* dynParticle = new G4DynamicParticle(theG4Hadron[currIdx], dir, ekin);
0386   currTrack = new G4Track(dynParticle, 0.0, vectProj);
0387   currTrack->SetStep(dummyStep);
0388
0389   theProjectile.Initialise(*currTrack);
0390   delete currTrack;
0391
0392   G4HadFinalState* result;
0393
0394   // elastic interaction
0395   if (random->flatShoot() * (intLengthElastic + intLengthInelastic) > intLengthElastic) {
0396     result = theDiffuseElastic->ApplyYourself(theProjectile, targetNucleus);
0397     G4ThreeVector dirnew = result->GetMomentumChange().unit();
0398     double cost = (dir * dirnew);
0399     double sint = std::sqrt((1. - cost) * (1. + cost));
0400
0401     curr4Mom.set(ptot * dirnew.x(), ptot * dirnew.y(), ptot * dirnew.z(), Particle.particle().momentum().e());
0402
0403     // Always create a daughter if the kink is large engough
0404     if (sint > theDistCut) {
0405       saveDaughter(Particle, curr4Mom, thePid);
0406     } else {
0407       Particle.particle().setMomentum(curr4Mom.px(), curr4Mom.py(), curr4Mom.pz(), curr4Mom.e());
0408     }
0409
0410     // inelastic interaction
0411   } else {
0412     // Bertini cascade for low-energy hadrons (except light anti-nuclei)
0413     // FTFP is applied above energy limit and for all anti-hyperons and anti-ions
0414     if (ekin <= theBertiniLimit && currIdx < 17) {
0415       result = theBertiniCascade->ApplyYourself(theProjectile, targetNucleus);
0416     } else {
0417       result = theHadronicModel->ApplyYourself(theProjectile, targetNucleus);
0418     }
0419     if (result) {
0420       int nsec = result->GetNumberOfSecondaries();
0421       if (0 < nsec) {
0422         result->SetTrafoToLab(theProjectile.GetTrafoToLab());
0423         _theUpdatedState.clear();
0424
0425         //std::cout << "   " << nsec << " secondaries" << std::endl;
0426         // Generate angle
0427         double phi = random->flatShoot() * CLHEP::twopi;
0428         theClosestChargedDaughterId = -1;
0429         distMin = 1e99;
0430
0431         // rotate and store secondaries
0432         for (int j = 0; j < nsec; ++j) {
0433           const G4DynamicParticle* dp = result->GetSecondary(j)->GetParticle();
0434           int thePid = dp->GetParticleDefinition()->GetPDGEncoding();
0435
0436           // rotate around primary direction
0437           curr4Mom = dp->Get4Momentum();
0438           curr4Mom.rotate(phi, vectProj);
0439           curr4Mom *= result->GetTrafoToLab();
0440           /*
0441         std::cout << j << ". " << dp->GetParticleDefinition()->GetParticleName()
0442         << "  " << thePid
0443         << "  " << curr4Mom*fact << std::endl;
0444       */
0445           // prompt 2-gamma decay for pi0, eta, eta'p
0446           if (111 == thePid || 221 == thePid || 331 == thePid) {
0447             theBoost = curr4Mom.boostVector();
0448             double e = 0.5 * dp->GetParticleDefinition()->GetPDGMass();
0449             double fi = random->flatShoot() * CLHEP::twopi;
0450             double cth = 2 * random->flatShoot() - 1.0;
0451             double sth = sqrt((1.0 - cth) * (1.0 + cth));
0452             G4LorentzVector lv(e * sth * cos(fi), e * sth * sin(fi), e * cth, e);
0453             lv.boost(theBoost);
0454             saveDaughter(Particle, lv, 22);
0455             curr4Mom -= lv;
0456             if (curr4Mom.e() > theEnergyLimit) {
0457               saveDaughter(Particle, curr4Mom, 22);
0458             }
0459           } else {
0460             if (curr4Mom.e() > theEnergyLimit + dp->GetParticleDefinition()->GetPDGMass()) {
0461               saveDaughter(Particle, curr4Mom, thePid);
0462             }
0463           }
0464         }
0465         result->Clear();
0466       }
0467     }
0468   }
0469 }
0470
0471 void NuclearInteractionFTFSimulator::saveDaughter(ParticlePropagator& Particle, const G4LorentzVector& lv, int pdgid) {
0472   unsigned int idx = _theUpdatedState.size();
0473   _theUpdatedState.emplace_back(
0474       makeParticle(Particle.particleDataTable(),
0475                    pdgid,
0476                    XYZTLorentzVector{lv.px() * fact, lv.py() * fact, lv.pz() * fact, lv.e() * fact},
0477                    Particle.particle().vertex()));
0478
0479   // Store the closest daughter index (for later tracking purposes, so charged particles only)
0480   double distance = distanceToPrimary(Particle.particle(), _theUpdatedState[idx]);
0481   // Find the closest daughter, if closer than a given upper limit.
0482   if (distance < distMin && distance < theDistCut) {
0483     distMin = distance;
0484     theClosestChargedDaughterId = idx;
0485   }
0486   // std::cout << _theUpdatedState[idx] << std::endl;
0487 }
0488
0489 double NuclearInteractionFTFSimulator::distanceToPrimary(const RawParticle& Particle,
0490                                                          const RawParticle& aDaughter) const {
0491   double distance = 2E99;
0492   // Compute the distance only for charged primaries
0493   if (fabs(Particle.charge()) > 1E-6) {
0494     // The secondary must have the same charge
0495     double chargeDiff = fabs(aDaughter.charge() - Particle.charge());
0496     if (fabs(chargeDiff) < 1E-6) {
0497       // Here are two distance definitions * to be tuned *
0498       switch (theDistAlgo) {
0499         case 1:
0500           // sin(theta12)
0501           distance = (aDaughter.Vect().Unit().Cross(Particle.Vect().Unit())).R();
0502           break;
0503
0504         case 2:
0505           // sin(theta12) * p1/p2
0506           distance = (aDaughter.Vect().Cross(Particle.Vect())).R() / aDaughter.Vect().Mag2();
0507           break;
0508
0509         default:
0510           // Should not happen
0511           break;
0512       }
0513     }
0514   }
0515   return distance;
0516 }