Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-03-17 11:00:54 

0001 // system headers
0002 #include <mutex>
0003 #include <cmath>
0004 #include <memory>
0005 
0006 // Framework Headers
0007 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0008 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0009 
0010 // Fast Sim headers
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 // Math
0022 #include "DataFormats/Math/interface/LorentzVector.h"
0023 
0024 // Geant4 headers
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 // Author: Vladimir Ivanchenko
0086 // Date: 20-Jan-2015
0087 //
0088 // Revision: Class structure modified to match SimplifiedGeometryPropagator
0089 //           Fixed a bug in which units where not properly converted from G4 to FastSim standard
0090 //           S. Kurz, 29 May 2017
0091 //////////////////////////////////////////////////////////
0092 
0093 typedef math::XYZTLorentzVector XYZTLorentzVector;
0094 
0095 namespace fastsim {
0096   //! Implementation of nuclear interactions of hadrons in the tracker layers (based on FTF model of Geant4).
0097   /*!
0098         Computes the probability for hadrons to interact with a nucleon of the tracker material (inelastically) and then sample nuclear interaction using FTF model of Geant4
0099         Also, another implementation of nuclear interactions can be used that is based on FullSim events (NuclearInteraction).
0100     */
0101   class NuclearInteractionFTF : public InteractionModel {
0102   public:
0103     //! Constructor.
0104     NuclearInteractionFTF(const std::string& name, const edm::ParameterSet& cfg);
0105 
0106     //! Default destructor.
0107     ~NuclearInteractionFTF() override;
0108 
0109     //! Perform the interaction.
0110     /*!
0111         \param particle The particle that interacts with the matter.
0112         \param layer The detector layer that interacts with the particle.
0113         \param secondaries Particles that are produced in the interaction (if any).
0114         \param random The Random Engine.
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;  //!< Number of G4 hadrons
0123     static const int npoints = 21;     //!< Number of steps to cover range of particle energies
0124 
0125     static const G4ParticleDefinition* theG4Hadron[numHadrons];  //!< The Geant4 particles
0126     static int theId[numHadrons];                                //!< The pdgIDs of the Geant4 particles
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;  //!< Bertini cascade for low-energy hadrons
0149     double theEnergyLimit;   //!< Minimum energy of interacting particle
0150 
0151     double theDistCut;  //!< Cut on deltaR for ClosestChargedDaughter algorithm (FastSim tracking)
0152 
0153     double intLengthElastic;    //!< Elastic interaction length
0154     double intLengthInelastic;  //!< Inelastic interaction length
0155 
0156     int currIdx;   //!< Index of particle in vector of Geant4 particles
0157     size_t index;  //!< Index for energy of particle (vectors parametrized as a function of energy)
0158 
0159     //! inelastic interaction length corrections per particle and energy
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     //! elastic interaction length corrections per particle and energy
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     //! inelastic interaction length in Silicon at 1 TeV per particle
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     //! elastic interaction length in Silicon at 1 TeV per particle
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   // Is this correct?
0292   // Thread safety
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 }  // namespace fastsim
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   // Angular distance of particle before/after interaction. If it too large, a daughter is created instead
0301   theDistCut = cfg.getParameter<double>("distCut");
0302   // Set energy limits for processes
0303   theBertiniLimit = cfg.getParameter<double>("bertiniLimit") * CLHEP::GeV;
0304   theEnergyLimit = cfg.getParameter<double>("energyLimit") * CLHEP::GeV;
0305 
0306   // FTF model
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   // Bertini Cascade
0319   theBertiniCascade = new G4CascadeInterface();
0320 
0321   theDiffuseElastic = new G4DiffuseElastic();
0322 
0323   // Geant4 particles and cross sections
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     // other Geant4 particles
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   // local objects
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   // fill projectile particle definitions
0379   dummyStep = new G4Step();
0380   dummyStep->SetPreStepPoint(new G4StepPoint());
0381 
0382   // target is always Silicon
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   // no valid PDGid
0400   //
0401   if (abs(thePid) <= 100 || abs(thePid) >= 1000000) {
0402     return;
0403   }
0404 
0405   // particle lost all its enrgy in previous interaction
0406   if (particle.momentum().E() < 1E-10) {
0407     return;
0408   }
0409 
0410   double radLengths = layer.getThickness(particle.position(), particle.momentum());
0411   // TEC layers have some fudge factor for nuclear interactions
0412   radLengths *= layer.getNuclearInteractionThicknessFactor();
0413   //
0414   // no material
0415   //
0416   if (radLengths < 1E-10) {
0417     return;
0418   }
0419 
0420   // get the G4 hadron
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         // neutral kaons
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   // no valid G4 hadron found
0440   if (!currParticle) {
0441     return;
0442   }
0443 
0444   // fill projectile for Geant4
0445   double mass = currParticle->GetPDGMass();
0446   double ekin = CLHEP::GeV * particle.momentum().e() - mass;
0447 
0448   // check interaction length
0449   intLengthElastic = nuclElLength[currIdx];
0450   intLengthInelastic = nuclInelLength[currIdx];
0451   if (0.0 == intLengthInelastic) {
0452     return;
0453   }
0454 
0455   // apply corrections
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   // Check position of nuclear interaction
0474   if (currInteractionLength > radLengths) {
0475     return;
0476   }
0477 
0478   // fill projectile for Geant4
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   // elastic interaction
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     // This vector is already in units of GeV (FastSim standard even if it is a G4Vector)
0503     curr4Mom.set(ptot * dirnew.x(), ptot * dirnew.y(), ptot * dirnew.z(), particle.momentum().e());
0504 
0505     // Always create a daughter if the kink is large enough
0506     if (sint > theDistCut) {
0507       // Create secondary
0508       secondaries.emplace_back(new fastsim::Particle(
0509           thePid, particle.position(), XYZTLorentzVector(curr4Mom.px(), curr4Mom.py(), curr4Mom.pz(), curr4Mom.e())));
0510 
0511       // ClosestChargedDaughter thing for FastSim tracking
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       // The particle is destroyed
0519       particle.momentum().SetXYZT(0., 0., 0., 0.);
0520     } else {
0521       // ... or just rotate particle
0522       particle.momentum().SetXYZT(curr4Mom.px(), curr4Mom.py(), curr4Mom.pz(), curr4Mom.e());
0523     }
0524 
0525     // inelastic interaction
0526   } else {
0527     // Bertini cascade for low-energy hadrons (except light anti-nuclei)
0528     // FTFP is applied above energy limit and for all anti-hyperons and anti-ions
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         // Generate angle
0541         double phi = random.flatShoot() * CLHEP::twopi;
0542 
0543         // rotate and store secondaries
0544         for (int j = 0; j < nsec; ++j) {
0545           const G4DynamicParticle* dp = result->GetSecondary(j)->GetParticle();
0546           int daughterID = dp->GetParticleDefinition()->GetPDGEncoding();
0547 
0548           // rotate around primary direction
0549           curr4Mom = dp->Get4Momentum();
0550           curr4Mom.rotate(phi, vectProj);
0551           curr4Mom *= result->GetTrafoToLab();
0552 
0553           // prompt 2-gamma decay for pi0, eta, eta'p
0554           // don't have a charge so the closest daughter thing does not have to be done
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             // create secondaries
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             // The mother particle is destroyed
0582             particle.momentum().SetXYZT(0., 0., 0., 0.);
0583           } else {
0584             // Copied from the original code, however I am not sure if all that is correct!
0585             // Energy of particles increases during the interaction!
0586             // std::cout << particle.momentum().e() << "(" << dynParticle->GetParticleDefinition()->GetPDGMass()/CLHEP::GeV << ") -> " << curr4Mom.e()/CLHEP::GeV << "(" << dp->GetParticleDefinition()->GetPDGMass()/CLHEP::GeV << ")" << std::endl;
0587 
0588             if (curr4Mom.e() > theEnergyLimit + dp->GetParticleDefinition()->GetPDGMass()) {
0589               // Create secondary
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               // ClosestChargedDaughter thing for FastSim tracking
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               // The mother particle is destroyed
0607               particle.momentum().SetXYZT(0., 0., 0., 0.);
0608             }
0609           }
0610         }
0611       }
0612 
0613       // Clear the final state
0614       result->Clear();
0615     }
0616   }
0617 }
0618 
0619 DEFINE_EDM_PLUGIN(fastsim::InteractionModelFactory, fastsim::NuclearInteractionFTF, "fastsim::NuclearInteractionFTF");
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!