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