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 // GEANT4 Class file 29 // 30 // 31 // File name: G4hCoulombScatteringModel 32 // 33 // Author: Vladimir Ivanchenko 34 // 35 // Creation date: 08.06.2012 from G4eCoulombSc 36 // 37 // Modifications: 38 // 39 // 40 // Class Description: 41 // 42 // ------------------------------------------- 43 // 44 //....oooOO0OOooo........oooOO0OOooo........oo 45 //....oooOO0OOooo........oooOO0OOooo........oo 46 47 #include "G4hCoulombScatteringModel.hh" 48 #include "G4PhysicalConstants.hh" 49 #include "G4SystemOfUnits.hh" 50 #include "Randomize.hh" 51 #include "G4DataVector.hh" 52 #include "G4ElementTable.hh" 53 #include "G4ParticleChangeForGamma.hh" 54 #include "G4Proton.hh" 55 #include "G4ParticleTable.hh" 56 #include "G4IonTable.hh" 57 #include "G4ProductionCutsTable.hh" 58 #include "G4NucleiProperties.hh" 59 #include "G4Pow.hh" 60 #include "G4NistManager.hh" 61 62 //....oooOO0OOooo........oooOO0OOooo........oo 63 64 G4hCoulombScatteringModel::G4hCoulombScatterin 65 : G4VEmModel("hCoulombScattering"), 66 cosThetaMin(1.0), 67 cosThetaMax(-1.0), 68 isCombined(combined) 69 { 70 fParticleChange = nullptr; 71 fNistManager = G4NistManager::Instance(); 72 theIonTable = G4ParticleTable::GetParticleT 73 theProton = G4Proton::Proton(); 74 currentMaterial = nullptr; 75 fixedCut = -1.0; 76 77 pCuts = nullptr; 78 79 recoilThreshold = 0.0; // by default does no 80 81 particle = nullptr; 82 currentCouple = nullptr; 83 wokvi = new G4WentzelVIRelXSection(); 84 85 currentMaterialIndex = 0; 86 mass = CLHEP::proton_mass_c2; 87 elecRatio = 0.0; 88 } 89 90 //....oooOO0OOooo........oooOO0OOooo........oo 91 92 G4hCoulombScatteringModel::~G4hCoulombScatteri 93 { 94 delete wokvi; 95 } 96 97 //....oooOO0OOooo........oooOO0OOooo........oo 98 99 void G4hCoulombScatteringModel::Initialise(con 100 const G4DataVector& cuts) 101 { 102 SetupParticle(part); 103 currentCouple = nullptr; 104 105 // defined theta limit between single and mu 106 isCombined = true; 107 G4double tet = PolarAngleLimit(); 108 109 if(tet <= 0.0) { 110 cosThetaMin = 1.0; 111 isCombined = false; 112 } else if(tet >= CLHEP::pi) { 113 cosThetaMin = -1.0; 114 } else { 115 cosThetaMin = std::cos(tet); 116 } 117 118 wokvi->Initialise(part, cosThetaMin); 119 /* 120 G4cout << "G4hCoulombScatteringModel: " << p 121 << " 1-cos(ThetaLimit)= " << 1 - cos 122 << " cos(thetaMax)= " << cosThetaMax 123 << G4endl; 124 */ 125 pCuts = &cuts; 126 //G4ProductionCutsTable::GetProductionCutsTa 127 /* 128 G4cout << "!!! G4hCoulombScatteringModel::In 129 << part->GetParticleName() << " cos(TetM 130 << " cos(TetMax)= " << cosThetaMax <<G4e 131 G4cout << "cut= " << (*pCuts)[0] << " cut1= 132 */ 133 if(!fParticleChange) { 134 fParticleChange = GetParticleChangeForGamm 135 } 136 if(IsMaster() && mass < CLHEP::GeV && part-> 137 InitialiseElementSelectors(part, cuts); 138 } 139 } 140 141 //....oooOO0OOooo........oooOO0OOooo........oo 142 143 void G4hCoulombScatteringModel::InitialiseLoca 144 G4VEmModel* masterModel) 145 { 146 SetElementSelectors(masterModel->GetElementS 147 } 148 149 //....oooOO0OOooo........oooOO0OOooo........oo 150 151 G4double 152 G4hCoulombScatteringModel::MinPrimaryEnergy(co 153 const G4ParticleDefinition* part 154 G4double) 155 { 156 SetupParticle(part); 157 158 // define cut using cuts for proton 159 G4double cut = 160 std::max(recoilThreshold, (*pCuts)[Current 161 162 // find out lightest element 163 const G4ElementVector* theElementVector = ma 164 std::size_t nelm = material->GetNumberOfElem 165 166 // select lightest element 167 G4int Z = 300; 168 for (std::size_t j=0; j<nelm; ++j) { 169 Z = std::min(Z,(*theElementVector)[j]->Get 170 } 171 G4int A = G4lrint(fNistManager->GetAtomicMas 172 G4double targetMass = G4NucleiProperties::Ge 173 G4double t = std::max(cut, 0.5*(cut + std::s 174 175 return t; 176 } 177 178 //....oooOO0OOooo........oooOO0OOooo........oo 179 180 G4double G4hCoulombScatteringModel::ComputeCro 181 const G4ParticleDefinition* p, 182 G4double kinEnergy, 183 G4double Z, G4double, 184 G4double cutEnergy, G4double) 185 { 186 //G4cout << "### G4hCoulombScatteringModel:: 187 //<< p->GetParticleName()<<" Z= "<<Z<<" e(Me 188 G4double cross = 0.0; 189 elecRatio = 0.0; 190 if(p != particle) { SetupParticle(p); } 191 192 // cross section is set to zero to avoid pro 193 if(kinEnergy <= 0.0) { return cross; } 194 DefineMaterial(CurrentCouple()); 195 196 G4int iz = G4lrint(Z); 197 G4double tmass = (1 == iz) ? proton_mass_c2 198 fNistManager->GetAtomicMassAmu(iz)*amu_c2; 199 wokvi->SetTargetMass(tmass); 200 201 G4double costmin = 202 wokvi->SetupKinematic(kinEnergy, currentMa 203 204 if(cosThetaMax < costmin) { 205 G4double cut = (0.0 < fixedCut) ? fixedCut 206 costmin = wokvi->SetupTarget(iz, cut); 207 G4double costmax = 208 (1 == iz && particle == theProton && cos 209 ? 0.0 : cosThetaMax; 210 if(costmin > costmax) { 211 cross = wokvi->ComputeNuclearCrossSectio 212 + wokvi->ComputeElectronCrossSection(costmin 213 } 214 /* 215 if(p->GetParticleName() == "mu+") 216 G4cout << "e(MeV)= " << kinEnergy/MeV << " c 217 << " 1-costmin= " << 1-costmin 218 << " 1-costmax= " << 1-costmax 219 << " 1-cosThetaMax= " << 1-cosThetaMax 220 << G4endl; 221 */ 222 } 223 return cross; 224 } 225 226 //....oooOO0OOooo........oooOO0OOooo........oo 227 228 void G4hCoulombScatteringModel::SampleSecondar 229 std::vector<G4DynamicParticle* 230 const G4MaterialCutsCouple* couple, 231 const G4DynamicParticle* dp, 232 G4double cutEnergy, 233 G4double) 234 { 235 G4double kinEnergy = dp->GetKineticEnergy(); 236 SetupParticle(dp->GetDefinition()); 237 DefineMaterial(couple); 238 239 // Choose nucleus 240 G4double cut = (0.0 < fixedCut) ? fixedCut : 241 242 const G4Element* elm = SelectRandomAtom(coup 243 kinEnergy,cut,kinEnergy); 244 245 G4int iz = elm->GetZasInt(); 246 G4int ia = SelectIsotopeNumber(elm); 247 G4double mass2 = G4NucleiProperties::GetNucl 248 249 wokvi->SetTargetMass(mass2); 250 wokvi->SetupKinematic(kinEnergy, currentMate 251 G4double costmin = wokvi->SetupTarget(iz, cu 252 G4double costmax = (1 == iz && particle == t 253 ? 0.0 : cosThetaMax; 254 if(costmin <= costmax) { return; } 255 256 G4double cross = wokvi->ComputeNuclearCrossS 257 G4double ecross = wokvi->ComputeElectronCros 258 G4double ratio = ecross/(cross + ecross); 259 260 G4ThreeVector newDirection = 261 wokvi->SampleSingleScattering(costmin, cos 262 263 // kinematics in the Lab system 264 G4double ptot = std::sqrt(kinEnergy*(kinEner 265 G4double e1 = mass + kinEnergy; 266 267 // Lab. system kinematics along projectile d 268 G4LorentzVector v0 = G4LorentzVector(0, 0, p 269 G4LorentzVector v1 = G4LorentzVector(0, 0, p 270 G4ThreeVector bst = v0.boostVector(); 271 v1.boost(-bst); 272 // CM projectile 273 G4double momCM = v1.pz(); 274 275 // Momentum after scattering of incident par 276 v1.setX(momCM*newDirection.x()); 277 v1.setY(momCM*newDirection.y()); 278 v1.setZ(momCM*newDirection.z()); 279 280 // CM--->Lab 281 v1.boost(bst); 282 283 G4ThreeVector dir = dp->GetMomentumDirection 284 newDirection = v1.vect().unit(); 285 newDirection.rotateUz(dir); 286 287 fParticleChange->ProposeMomentumDirection(ne 288 289 // recoil 290 v0 -= v1; 291 G4double trec = std::max(v0.e() - mass2, 0.0 292 G4double edep = 0.0; 293 294 G4double tcut = recoilThreshold; 295 if(pCuts) { tcut= std::max(tcut,(*pCuts)[cur 296 297 if(trec > tcut) { 298 G4ParticleDefinition* ion = theIonTable->G 299 newDirection = v0.vect().unit(); 300 newDirection.rotateUz(dir); 301 auto newdp = new G4DynamicParticle(ion, ne 302 fvect->push_back(newdp); 303 } else if(trec > 0.0) { 304 edep = trec; 305 fParticleChange->ProposeNonIonizingEnergyD 306 } 307 308 // finelize primary energy and energy balanc 309 G4double finalT = v1.e() - mass; 310 if(finalT < 0.0) { 311 edep += finalT; 312 finalT = 0.0; 313 } 314 edep = std::max(edep, 0.0); 315 fParticleChange->SetProposedKineticEnergy(fi 316 fParticleChange->ProposeLocalEnergyDeposit(e 317 } 318 319 //....oooOO0OOooo........oooOO0OOooo........oo 320