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