Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-05-10 02:20:41 

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