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 // 27 // ------------------------------------------- 28 // 29 // GEANT4 Class file 30 // 31 // 32 // File name: G4DynamicParticleIonisation 33 // 34 // Author: Vladimir Ivanchenko 35 // 36 // Creation date: 17.08.2024 37 // 38 // ------------------------------------------- 39 // 40 //....oooOO0OOooo........oooOO0OOooo........oo 41 //....oooOO0OOooo........oooOO0OOooo........oo 42 43 #include "G4DynamicParticleIonisation.hh" 44 #include "G4PhysicalConstants.hh" 45 #include "G4SystemOfUnits.hh" 46 #include "G4DynamicParticleFluctuation.hh" 47 #include "G4EmSecondaryParticleType.hh" 48 #include "G4Electron.hh" 49 #include "G4EmParameters.hh" 50 #include "G4EmProcessSubType.hh" 51 #include "G4LossTableManager.hh" 52 #include "G4MaterialCutsCouple.hh" 53 #include "G4ProductionCutsTable.hh" 54 #include "G4Material.hh" 55 #include "G4Step.hh" 56 #include "G4Track.hh" 57 #include "G4Log.hh" 58 59 namespace 60 { 61 constexpr G4double ekinLimit = 0.2*CLHEP::Me 62 const G4double twoln10 = 2*G4Log(10.0); 63 } 64 65 //....oooOO0OOooo........oooOO0OOooo........oo 66 67 G4DynamicParticleIonisation::G4DynamicParticle 68 : G4VContinuousDiscreteProcess("dynPartIoni" 69 { 70 SetVerboseLevel(1); 71 SetProcessSubType(fDynamicIonisation); 72 theElectron = G4Electron::Electron(); 73 74 lManager = G4LossTableManager::Instance(); 75 lManager->Register(this); 76 77 fUrban = new G4DynamicParticleFluctuation(); 78 79 // define these flags only once 80 auto param = G4EmParameters::Instance(); 81 fFluct = param->LossFluctuation(); 82 fLinLimit = 5*param->LinearLossLimit(); 83 } 84 85 //....oooOO0OOooo........oooOO0OOooo........oo 86 87 G4DynamicParticleIonisation::~G4DynamicParticl 88 { 89 lManager->DeRegister(this); 90 } 91 92 //....oooOO0OOooo........oooOO0OOooo........oo 93 94 void 95 G4DynamicParticleIonisation::BuildPhysicsTable 96 { 97 auto theCoupleTable = G4ProductionCutsTable: 98 fCuts = theCoupleTable->GetEnergyCutsVector( 99 } 100 101 //....oooOO0OOooo........oooOO0OOooo........oo 102 103 void G4DynamicParticleIonisation::PreStepIniti 104 { 105 fCouple = track.GetMaterialCutsCouple(); 106 fMaterial = fCouple->GetMaterial(); 107 auto dpart = track.GetDynamicParticle(); 108 fEkinPreStep = dpart->GetKineticEnergy(); 109 fMass = std::max(dpart->GetMass(), CLHEP::el 110 fCharge = dpart->GetCharge()/CLHEP::eplus; 111 fRatio = fMass/CLHEP::proton_mass_c2; 112 fLowestEkin = ekinLimit*fRatio; 113 G4double tau = fEkinPreStep/fMass; 114 G4double ratio = CLHEP::electron_mass_c2/fMa 115 fTmax = 2.0*CLHEP::electron_mass_c2*tau*(tau 116 (1. + 2.0*(tau + 1.)*ratio + ratio*ratio); 117 fCut = (*fCuts)[fCouple->GetIndex()]; 118 fCut = std::max(fCut, fMaterial->GetIonisati 119 fCut = std::min(fCut, fTmax); 120 } 121 122 //....oooOO0OOooo........oooOO0OOooo........oo 123 124 G4double G4DynamicParticleIonisation::AlongSte 125 const G4Trac 126 G4GPILSelect 127 { 128 *selection = CandidateForSelection; 129 130 // no step limit for the time being 131 return DBL_MAX; 132 } 133 134 //....oooOO0OOooo........oooOO0OOooo........oo 135 136 G4double G4DynamicParticleIonisation::PostStep 137 const G4Trac 138 G4ForceCondi 139 { 140 *condition = NotForced; 141 G4double x = DBL_MAX; 142 G4double xsec = 0.0; 143 144 PreStepInitialisation(track); 145 146 if (fCharge != 0.0) { 147 xsec = ComputeCrossSection(fEkinPreStep); 148 } 149 150 if (xsec <= 0.0) { 151 theNumberOfInteractionLengthLeft = -1.0; 152 currentInteractionLength = DBL_MAX; 153 154 } else { 155 if (theNumberOfInteractionLengthLeft < 0.0 156 157 theNumberOfInteractionLengthLeft = -G4Lo 158 theInitialNumberOfInteractionLength = th 159 160 } else if(currentInteractionLength < DBL_M 161 162 // subtract NumberOfInteractionLengthLef 163 theNumberOfInteractionLengthLeft -= 164 previousStepSize/currentInteractionLen 165 166 theNumberOfInteractionLengthLeft = 167 std::max(theNumberOfInteractionLengthL 168 } 169 currentInteractionLength = 1.0/xsec; 170 x = theNumberOfInteractionLengthLeft * cur 171 } 172 #ifdef G4VERBOSE 173 if (verboseLevel>2) { 174 G4cout << "G4DynamicParticleIonisation::Po 175 G4cout << " Process: " << GetProcessName( 176 << " for unknown particle Mass(GeV) 177 << " charge=" << fCharge 178 << " Material " << fMaterial->GetN 179 << " Ekin(MeV)=" << fEkinPreStep/C 180 << " MFP(cm)=" << currentInteracti 181 << " ProposedLength(cm)=" << x/CLH 182 } 183 #endif 184 return x; 185 } 186 187 //....oooOO0OOooo........oooOO0OOooo........oo 188 189 G4VParticleChange* 190 G4DynamicParticleIonisation::AlongStepDoIt(con 191 const G4Step& step) 192 { 193 fParticleChange.InitializeForAlongStep(track 194 195 // no energy loss 196 if (fCharge == 0.0) { return &fParticleChang 197 198 // stop low-energy object 199 if (fEkinPreStep <= fLowestEkin) { 200 fParticleChange.SetProposedKineticEnergy(0 201 fParticleChange.ProposeLocalEnergyDeposit( 202 return &fParticleChange; 203 } 204 205 G4double length = step.GetStepLength(); 206 G4double dedxPre = ComputeDEDX(fEkinPreStep) 207 G4double eloss = dedxPre*length; 208 G4double ekinPostStep = fEkinPreStep - eloss 209 210 // correction for large step if it is not th 211 if (fEkinPreStep*fLinLimit < eloss && ekinPo 212 G4double dedxPost = ComputeDEDX(ekinPostSt 213 eloss = (eloss + dedxPost*length)*0.5; 214 } 215 216 // do not sample fluctuations at the last st 217 if (fFluct && fEkinPreStep > eloss) { 218 eloss = fUrban->SampleFluctuations(fCouple 219 fCut, fTmax, length, eloss); 220 } 221 222 ekinPostStep = fEkinPreStep - eloss; 223 224 // stop low-energy object 225 if (ekinPostStep <= fLowestEkin) { 226 fParticleChange.SetProposedKineticEnergy(0 227 fParticleChange.ProposeLocalEnergyDeposit( 228 } else { 229 fParticleChange.SetProposedKineticEnergy(e 230 fParticleChange.ProposeLocalEnergyDeposit( 231 } 232 return &fParticleChange; 233 } 234 235 //....oooOO0OOooo........oooOO0OOooo........oo 236 237 G4VParticleChange* 238 G4DynamicParticleIonisation::PostStepDoIt(cons 239 { 240 theNumberOfInteractionLengthLeft = -1.0; 241 fParticleChange.InitializeForPostStep(track) 242 243 auto dp = track.GetDynamicParticle(); 244 G4double kinEnergy = dp->GetKineticEnergy(); 245 const G4double totEnergy = kinEnergy + fMass 246 const G4double beta2 = kinEnergy*(kinEnergy 247 248 G4double deltaKinEnergy, f; 249 250 CLHEP::HepRandomEngine* rndmEngineMod = G4Ra 251 G4double rndm[2]; 252 253 // sampling without nuclear size effect 254 do { 255 rndmEngineMod->flatArray(2, rndm); 256 deltaKinEnergy = fCut*fTmax/(fCut*(1.0 - r 257 f = 1.0 - beta2*deltaKinEnergy/fTmax; 258 // Loop checking, 14-Aug-2024, Vladimir Iv 259 } while( rndm[1] > f); 260 261 G4double deltaMomentum = 262 std::sqrt(deltaKinEnergy * (deltaKinEnergy 263 G4double cost = deltaKinEnergy * (totEnergy 264 (deltaMomentum * dp->GetTotalMomentum()) 265 cost = std::min(cost, 1.0); 266 const G4double sint = std::sqrt((1.0 - cost) 267 const G4double phi = CLHEP::twopi*rndmEngine 268 269 G4ThreeVector deltaDirection(sint*std::cos(p 270 deltaDirection.rotateUz(dp->GetMomentumDirec 271 272 // create G4DynamicParticle object for delta 273 auto delta = new G4DynamicParticle(theElectr 274 auto t = new G4Track(delta, track.GetGlobalT 275 t->SetTouchableHandle(track.GetTouchableHand 276 t->SetCreatorModelID(fSecID); 277 fParticleChange.AddSecondary(t); 278 279 // Change kinematics of primary particle 280 kinEnergy -= deltaKinEnergy; 281 G4ThreeVector finalP = dp->GetMomentum() - d 282 finalP = finalP.unit(); 283 284 fParticleChange.SetProposedKineticEnergy(kin 285 fParticleChange.SetProposedMomentumDirection 286 return &fParticleChange; 287 } 288 289 //....oooOO0OOooo........oooOO0OOooo........oo 290 291 G4double G4DynamicParticleIonisation::ComputeD 292 { 293 G4double tau = ekin/fMass; 294 G4double gam = tau + 1.0; 295 G4double bg2 = tau * (tau + 2.0); 296 G4double beta2 = bg2/(gam*gam); 297 G4double xc = fCut/fTmax; 298 299 G4double exc = fMaterial->GetIonisation()-> 300 G4double exc2 = exc*exc; 301 302 // general Bethe-Bloch formula 303 G4double dedx = G4Log(2.0*CLHEP::electron_ma 304 305 // density correction 306 G4double x = G4Log(bg2)/twoln10; 307 dedx -= fMaterial->GetIonisation()->DensityC 308 309 // now compute the total ionization loss per 310 dedx *= CLHEP::twopi_mc2_rcl2*fCharge*fCharg 311 dedx = std::max(dedx, 0.0); 312 return dedx; 313 } 314 315 //....oooOO0OOooo........oooOO0OOooo........oo 316 317 G4double G4DynamicParticleIonisation::ComputeC 318 { 319 G4double cross = 0.0; 320 if (fCut < fTmax) { 321 322 G4double totEnergy = ekin + fMass; 323 G4double energy2 = totEnergy*totEnergy; 324 G4double beta2 = ekin*(ekin + 2.0*fMass)/e 325 326 cross = (fTmax - fCut)/(fCut*fTmax*beta2) 327 cross *= CLHEP::twopi_mc2_rcl2*fCharge*fCh 328 } 329 return cross; 330 } 331 332 //....oooOO0OOooo........oooOO0OOooo........oo 333 334 G4double G4DynamicParticleIonisation::GetMeanF 335 336 { 337 // Note: this method is not used at run-time 338 // It might be eventually refined late 339 *condition = NotForced; 340 return DBL_MAX; 341 } 342 343 //....oooOO0OOooo........oooOO0OOooo........oo 344 345 G4double G4DynamicParticleIonisation::GetConti 346 G4double, G4double&) 347 { 348 return DBL_MAX; 349 } 350 351 //....oooOO0OOooo........oooOO0OOooo........oo 352 353 void G4DynamicParticleIonisation::ProcessDescr 354 { 355 out << "G4DynamicParticleIonisation: dynamic 356 } 357 358 //....oooOO0OOooo........oooOO0OOooo........oo 359