Geant4 Cross Reference |
1 // 1 2 // ******************************************* 3 // * License and Disclaimer 4 // * 5 // * The Geant4 software is copyright of th 6 // * the Geant4 Collaboration. It is provided 7 // * conditions of the Geant4 Software License 8 // * LICENSE and available at http://cern.ch/ 9 // * include a list of copyright holders. 10 // * 11 // * Neither the authors of this software syst 12 // * institutes,nor the agencies providing fin 13 // * work make any representation or warran 14 // * regarding this software system or assum 15 // * use. Please see the license in the file 16 // * for the full disclaimer and the limitatio 17 // * 18 // * This code implementation is the result 19 // * technical work of the GEANT4 collaboratio 20 // * By using, copying, modifying or distri 21 // * any work based on the software) you ag 22 // * use in resulting scientific publicati 23 // * acceptance of all terms of the Geant4 Sof 24 // ******************************************* 25 // 26 #include "Par03EMShowerModel.hh" 27 28 #include "Par03EMShowerMessenger.hh" 29 30 #include "G4Electron.hh" 31 #include "G4FastHit.hh" 32 #include "G4FastSimHitMaker.hh" 33 #include "G4Gamma.hh" 34 #include "G4Positron.hh" 35 #include "G4SystemOfUnits.hh" 36 #include "G4UnitsTable.hh" 37 #include "Randomize.hh" 38 39 Par03EMShowerModel::Par03EMShowerModel(G4Strin 40 : G4VFastSimulationModel(aModelName, aEnvelo 41 fMessenger(new Par03EMShowerMessenger(this 42 fHitMaker(new G4FastSimHitMaker) 43 {} 44 45 //....oooOO0OOooo........oooOO0OOooo........oo 46 47 Par03EMShowerModel::Par03EMShowerModel(G4Strin 48 : G4VFastSimulationModel(aModelName), 49 fMessenger(new Par03EMShowerMessenger(this 50 fHitMaker(new G4FastSimHitMaker) 51 {} 52 53 //....oooOO0OOooo........oooOO0OOooo........oo 54 55 Par03EMShowerModel::~Par03EMShowerModel() = de 56 57 //....oooOO0OOooo........oooOO0OOooo........oo 58 59 G4bool Par03EMShowerModel::IsApplicable(const 60 { 61 return &aParticleType == G4Electron::Electro 62 || &aParticleType == G4Positron::Posi 63 || &aParticleType == G4Gamma::GammaDe 64 } 65 66 //....oooOO0OOooo........oooOO0OOooo........oo 67 68 G4bool Par03EMShowerModel::ModelTrigger(const 69 { 70 // Check energy 71 if (aFastTrack.GetPrimaryTrack()->GetKinetic 72 return false; 73 } 74 // Check length of detector 75 // Calculate depth of the detector along sho 76 // will fit inside. Required max shower dept 77 // can be changed with UI command `/Par03/fa 78 G4double X0 = aFastTrack.GetPrimaryTrack()-> 79 auto particleDirection = aFastTrack.GetPrima 80 auto particlePosition = aFastTrack.GetPrimar 81 G4double detectorDepthInMM = 82 aFastTrack.GetEnvelopeSolid()->DistanceToO 83 G4double detectorDepthInX0 = detectorDepthIn 84 // check if detector depth is sufficient to 85 if (detectorDepthInX0 < fLongMaxDepth) { 86 return false; 87 } 88 return true; 89 } 90 91 //....oooOO0OOooo........oooOO0OOooo........oo 92 93 void Par03EMShowerModel::DoIt(const G4FastTrac 94 { 95 // Remove particle from further processing b 96 aFastStep.KillPrimaryTrack(); 97 aFastStep.ProposePrimaryTrackPathLength(0.0) 98 G4double energy = aFastTrack.GetPrimaryTrack 99 // No need to create any deposit, it will be 100 // G4FastSimHitMaker that will call the sens 101 aFastStep.ProposeTotalEnergyDeposited(0); 102 auto particlePosition = aFastTrack.GetPrimar 103 auto particleDirection = aFastTrack.GetPrima 104 105 // Calculate how to create energy deposits 106 // Following PDG 33.5 chapter 107 // material calculation assumes homogeneous 108 auto material = aFastTrack.GetPrimaryTrack() 109 G4double materialX0 = material->GetRadlen(); 110 G4double materialZ = material->GetZ(); 111 // EC estimation follows PDG fit to solids i 112 G4double materialEc = 610 * MeV / (materialZ 113 // RM estimation follows PDG Eq. (33.37) (rm 114 G4double materialRM = 21.2052 * MeV * materi 115 G4double particleY = energy / materialEc; 116 // Estimate shower maximum and alpha paramet 117 // that describes the longitudinal profile ( 118 // unless alpha is specified by UI command 119 if (fAlpha < 0) { 120 // from PDG Eq. (33.36) 121 G4double particleTmax = std::log(particleY 122 if (aFastTrack.GetPrimaryTrack()->GetParti 123 particleTmax += 0.5; 124 } 125 else { 126 particleTmax -= 0.5; 127 } 128 fAlpha = particleTmax * fBeta + 1; 129 } 130 // Unless sigma of Gaussian distribution des 131 // is specified by UI command, use value cal 132 if (fSigma < 0) { 133 // 90% of shower is contained within 1 * R 134 // 1.645 * std dev of Gaussian contains 90 135 fSigma = materialRM / 1.645; 136 } 137 138 // Calculate rotation matrix along the parti 139 // It will rotate the shower axes to match t 140 G4RotationMatrix rotMatrix = G4RotationMatri 141 double particleTheta = particleDirection.the 142 double particlePhi = particleDirection.phi() 143 double epsilon = 1e-3; 144 rotMatrix.rotateY(particleTheta); 145 // do not use (random) phi if x==y==0 146 if (!(std::fabs(particleDirection.x()) < eps 147 rotMatrix.rotateZ(particlePhi); 148 149 // Create hits 150 // First use rejecton sampling to sample fro 151 // then get random numbers from uniform dist 152 // from Gaussian for radius 153 G4ThreeVector position; 154 G4double gammaMax = Gamma((fAlpha - 1) / fBe 155 G4int generatedHits = 0; 156 while (generatedHits < fNbOfHits) { 157 G4double random1 = G4UniformRand() * fLong 158 G4double random2 = G4UniformRand() * gamma 159 if (Gamma(random1, fAlpha, fBeta) >= rando 160 // Generate corresponding rho (phi) from 161 G4double phiPosition = G4UniformRand() * 162 G4double rhoPosition = G4RandGauss::shoo 163 position = particlePosition 164 + rotMatrix 165 * G4ThreeVector(rhoPositi 166 rhoPositi 167 // Create energy deposit in the detector 168 // This will call appropriate sensitive 169 fHitMaker->make(G4FastHit(position, ener 170 generatedHits++; 171 } 172 } 173 } 174 175 //....oooOO0OOooo........oooOO0OOooo........oo 176 177 void Par03EMShowerModel::Print() const 178 { 179 G4cout << "Par03EMShowerModel: " << G4endl; 180 G4cout << "Gaussian distribution (transverse 181 << G4BestUnit(fSigma, "Length") << G4 182 if (fSigma < 0) 183 G4cout << "Negative sigma value means that 184 "from the value of the Moliere r 185 "taking into account that 90% of 186 "distribution (from mu - 1.645 s 187 "corresponds to area within 1 Mo 188 << G4endl; 189 G4cout << "Gamma distribution (along shower 190 << ", max depth = " << fLongMaxDepth 191 if (fAlpha < 0) 192 G4cout << "Negative alpha value means that 193 "from the critical energy of the 194 "type, and beta parameter.\n alp 195 "ln(E/E_C) + C\n where E is part 196 "energy in the material, and con 197 "0.5 for photons (Eq. (33.36) fr 198 << G4endl; 199 G4cout << "Number of created energy deposits 200 } 201