Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 // 25 // >> 26 // $Id: G4hCoulombScatteringModel.cc 76536 2013-11-12 15:17:41Z gcosmo $ >> 27 // 26 // ------------------------------------------- 28 // ------------------------------------------------------------------- 27 // 29 // 28 // GEANT4 Class file 30 // GEANT4 Class file 29 // 31 // 30 // 32 // 31 // File name: G4hCoulombScatteringModel 33 // File name: G4hCoulombScatteringModel 32 // 34 // 33 // Author: Vladimir Ivanchenko 35 // Author: Vladimir Ivanchenko 34 // 36 // 35 // Creation date: 08.06.2012 from G4eCoulombSc 37 // Creation date: 08.06.2012 from G4eCoulombScatteringModel 36 // 38 // 37 // Modifications: 39 // Modifications: 38 // 40 // 39 // 41 // 40 // Class Description: 42 // Class Description: 41 // 43 // 42 // ------------------------------------------- 44 // ------------------------------------------------------------------- 43 // 45 // 44 //....oooOO0OOooo........oooOO0OOooo........oo 46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 45 //....oooOO0OOooo........oooOO0OOooo........oo 47 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 46 48 47 #include "G4hCoulombScatteringModel.hh" 49 #include "G4hCoulombScatteringModel.hh" 48 #include "G4PhysicalConstants.hh" 50 #include "G4PhysicalConstants.hh" 49 #include "G4SystemOfUnits.hh" 51 #include "G4SystemOfUnits.hh" 50 #include "Randomize.hh" 52 #include "Randomize.hh" 51 #include "G4DataVector.hh" 53 #include "G4DataVector.hh" 52 #include "G4ElementTable.hh" 54 #include "G4ElementTable.hh" 53 #include "G4ParticleChangeForGamma.hh" 55 #include "G4ParticleChangeForGamma.hh" 54 #include "G4Proton.hh" 56 #include "G4Proton.hh" 55 #include "G4ParticleTable.hh" 57 #include "G4ParticleTable.hh" 56 #include "G4IonTable.hh" 58 #include "G4IonTable.hh" 57 #include "G4ProductionCutsTable.hh" 59 #include "G4ProductionCutsTable.hh" 58 #include "G4NucleiProperties.hh" 60 #include "G4NucleiProperties.hh" 59 #include "G4Pow.hh" 61 #include "G4Pow.hh" >> 62 #include "G4LossTableManager.hh" >> 63 #include "G4LossTableBuilder.hh" 60 #include "G4NistManager.hh" 64 #include "G4NistManager.hh" 61 65 62 //....oooOO0OOooo........oooOO0OOooo........oo 66 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 63 67 64 G4hCoulombScatteringModel::G4hCoulombScatterin << 68 using namespace std; 65 : G4VEmModel("hCoulombScattering"), << 69 >> 70 G4hCoulombScatteringModel::G4hCoulombScatteringModel(const G4String& nam) >> 71 : G4VEmModel(nam), 66 cosThetaMin(1.0), 72 cosThetaMin(1.0), 67 cosThetaMax(-1.0), 73 cosThetaMax(-1.0), 68 isCombined(combined) << 74 isInitialised(false) 69 { 75 { 70 fParticleChange = nullptr; << 76 fParticleChange = 0; 71 fNistManager = G4NistManager::Instance(); 77 fNistManager = G4NistManager::Instance(); 72 theIonTable = G4ParticleTable::GetParticleT << 78 theIonTable = G4ParticleTable::GetParticleTable()->GetIonTable(); 73 theProton = G4Proton::Proton(); << 79 theProton = G4Proton::Proton(); 74 currentMaterial = nullptr; << 80 currentMaterial = 0; 75 fixedCut = -1.0; << 76 81 77 pCuts = nullptr; << 82 pCuts = 0; 78 83 79 recoilThreshold = 0.0; // by default does no << 84 lowEnergyThreshold = 1*keV; // particle will be killed for lower energy >> 85 recoilThreshold = 0.*keV; // by default does not work 80 86 81 particle = nullptr; << 87 particle = 0; 82 currentCouple = nullptr; << 88 currentCouple = 0; 83 wokvi = new G4WentzelVIRelXSection(); 89 wokvi = new G4WentzelVIRelXSection(); 84 90 85 currentMaterialIndex = 0; 91 currentMaterialIndex = 0; 86 mass = CLHEP::proton_mass_c2; << 92 >> 93 cosTetMinNuc = 1.0; >> 94 cosTetMaxNuc = -1.0; 87 elecRatio = 0.0; 95 elecRatio = 0.0; >> 96 mass = proton_mass_c2; 88 } 97 } 89 98 90 //....oooOO0OOooo........oooOO0OOooo........oo 99 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 91 100 92 G4hCoulombScatteringModel::~G4hCoulombScatteri 101 G4hCoulombScatteringModel::~G4hCoulombScatteringModel() 93 { 102 { 94 delete wokvi; 103 delete wokvi; 95 } 104 } 96 105 97 //....oooOO0OOooo........oooOO0OOooo........oo 106 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 98 107 99 void G4hCoulombScatteringModel::Initialise(con << 108 void G4hCoulombScatteringModel::Initialise(const G4ParticleDefinition* p, 100 const G4DataVector& cuts) 109 const G4DataVector& cuts) 101 { 110 { 102 SetupParticle(part); << 111 SetupParticle(p); 103 currentCouple = nullptr; << 112 currentCouple = 0; 104 << 113 cosThetaMin = cos(PolarAngleLimit()); 105 // defined theta limit between single and mu << 114 wokvi->Initialise(p, cosThetaMin); 106 isCombined = true; << 115 /* 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 116 G4cout << "G4hCoulombScatteringModel: " << particle->GetParticleName() 121 << " 1-cos(ThetaLimit)= " << 1 - cos 117 << " 1-cos(ThetaLimit)= " << 1 - cosThetaMin 122 << " cos(thetaMax)= " << cosThetaMax 118 << " cos(thetaMax)= " << cosThetaMax 123 << G4endl; 119 << G4endl; 124 */ 120 */ 125 pCuts = &cuts; << 121 pCuts = G4ProductionCutsTable::GetProductionCutsTable()->GetEnergyCutsVector(3); 126 //G4ProductionCutsTable::GetProductionCutsTa << 122 //G4cout << "!!! G4hCoulombScatteringModel::Initialise for " 127 /* << 123 // << p->GetParticleName() << " cos(TetMin)= " << cosThetaMin 128 G4cout << "!!! G4hCoulombScatteringModel::In << 124 // << " cos(TetMax)= " << cosThetaMax <<G4endl; 129 << part->GetParticleName() << " cos(TetM << 125 // G4cout << "cut0= " << cuts[0] << " cut1= " << cuts[1] << G4endl; 130 << " cos(TetMax)= " << cosThetaMax <<G4e << 126 if(!isInitialised) { 131 G4cout << "cut= " << (*pCuts)[0] << " cut1= << 127 isInitialised = true; 132 */ << 133 if(!fParticleChange) { << 134 fParticleChange = GetParticleChangeForGamm 128 fParticleChange = GetParticleChangeForGamma(); 135 } 129 } 136 if(IsMaster() && mass < CLHEP::GeV && part-> << 130 if(mass < GeV) { 137 InitialiseElementSelectors(part, cuts); << 131 InitialiseElementSelectors(p,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 } 132 } 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 } 133 } 177 134 178 //....oooOO0OOooo........oooOO0OOooo........oo 135 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 179 136 180 G4double G4hCoulombScatteringModel::ComputeCro 137 G4double G4hCoulombScatteringModel::ComputeCrossSectionPerAtom( 181 const G4ParticleDefinition* p, 138 const G4ParticleDefinition* p, 182 G4double kinEnergy, 139 G4double kinEnergy, 183 G4double Z, G4double, 140 G4double Z, G4double, 184 G4double cutEnergy, G4double) 141 G4double cutEnergy, G4double) 185 { 142 { 186 //G4cout << "### G4hCoulombScatteringModel:: 143 //G4cout << "### G4hCoulombScatteringModel::ComputeCrossSectionPerAtom for " 187 //<< p->GetParticleName()<<" Z= "<<Z<<" e(Me << 144 // << p->GetParticleName()<<" Z= "<<Z<<" e(MeV)= "<< kinEnergy/MeV << G4endl; 188 G4double cross = 0.0; 145 G4double cross = 0.0; 189 elecRatio = 0.0; << 190 if(p != particle) { SetupParticle(p); } 146 if(p != particle) { SetupParticle(p); } 191 147 192 // cross section is set to zero to avoid pro 148 // cross section is set to zero to avoid problems in sample secondary 193 if(kinEnergy <= 0.0) { return cross; } 149 if(kinEnergy <= 0.0) { return cross; } 194 DefineMaterial(CurrentCouple()); 150 DefineMaterial(CurrentCouple()); 195 << 151 cosTetMinNuc = wokvi->SetupKinematic(kinEnergy, currentMaterial); 196 G4int iz = G4lrint(Z); << 152 if(cosThetaMax < cosTetMinNuc) { 197 G4double tmass = (1 == iz) ? proton_mass_c2 << 153 G4int iz = G4int(Z); 198 fNistManager->GetAtomicMassAmu(iz)*amu_c2; << 154 cosTetMinNuc = wokvi->SetupTarget(iz, cutEnergy); 199 wokvi->SetTargetMass(tmass); << 155 cosTetMaxNuc = cosThetaMax; 200 << 156 if(iz == 1 && cosTetMaxNuc < 0.0 && particle == theProton) { 201 G4double costmin = << 157 cosTetMaxNuc = 0.0; 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 } 158 } 214 /* << 159 cross = wokvi->ComputeNuclearCrossSection(cosTetMinNuc, cosTetMaxNuc); 215 if(p->GetParticleName() == "mu+") << 160 elecRatio = wokvi->ComputeElectronCrossSection(cosTetMinNuc, cosThetaMax); >> 161 cross += elecRatio; >> 162 if(cross > 0.0) { elecRatio /= cross; } >> 163 } >> 164 /* >> 165 if(p->GetParticleName() == "e-") 216 G4cout << "e(MeV)= " << kinEnergy/MeV << " c 166 G4cout << "e(MeV)= " << kinEnergy/MeV << " cross(b)= " << cross/barn 217 << " 1-costmin= " << 1-costmin << 167 << " 1-cosTetMinNuc= " << 1-cosTetMinNuc 218 << " 1-costmax= " << 1-costmax << 168 << " 1-cosTetMaxNuc= " << 1-cosTetMaxNuc 219 << " 1-cosThetaMax= " << 1-cosThetaMax << 220 << G4endl; 169 << G4endl; 221 */ << 170 */ 222 } << 223 return cross; 171 return cross; 224 } 172 } 225 173 226 //....oooOO0OOooo........oooOO0OOooo........oo 174 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 227 175 228 void G4hCoulombScatteringModel::SampleSecondar 176 void G4hCoulombScatteringModel::SampleSecondaries( 229 std::vector<G4DynamicParticle* 177 std::vector<G4DynamicParticle*>* fvect, 230 const G4MaterialCutsCouple* couple, 178 const G4MaterialCutsCouple* couple, 231 const G4DynamicParticle* dp, 179 const G4DynamicParticle* dp, 232 G4double cutEnergy, 180 G4double cutEnergy, 233 G4double) 181 G4double) 234 { 182 { 235 G4double kinEnergy = dp->GetKineticEnergy(); 183 G4double kinEnergy = dp->GetKineticEnergy(); >> 184 >> 185 // absorb particle below low-energy limit to avoid situation >> 186 // when a particle has no energy loss >> 187 if(kinEnergy < lowEnergyThreshold) { >> 188 fParticleChange->SetProposedKineticEnergy(0.0); >> 189 fParticleChange->ProposeLocalEnergyDeposit(kinEnergy); >> 190 fParticleChange->ProposeNonIonizingEnergyDeposit(kinEnergy); >> 191 return; >> 192 } 236 SetupParticle(dp->GetDefinition()); 193 SetupParticle(dp->GetDefinition()); 237 DefineMaterial(couple); 194 DefineMaterial(couple); 238 195 >> 196 //G4cout << "G4hCoulombScatteringModel::SampleSecondaries e(MeV)= " >> 197 // << kinEnergy << " " << particle->GetParticleName() >> 198 // << " cut= " << cutEnergy<< G4endl; >> 199 239 // Choose nucleus 200 // Choose nucleus 240 G4double cut = (0.0 < fixedCut) ? fixedCut : << 201 const G4Element* currentElement = >> 202 SelectRandomAtom(couple,particle,kinEnergy,cutEnergy,kinEnergy); 241 203 242 const G4Element* elm = SelectRandomAtom(coup << 204 G4double Z = currentElement->GetZ(); 243 kinEnergy,cut,kinEnergy); << 244 205 245 G4int iz = elm->GetZasInt(); << 206 if(ComputeCrossSectionPerAtom(particle,kinEnergy, Z, 246 G4int ia = SelectIsotopeNumber(elm); << 207 kinEnergy, cutEnergy, kinEnergy) == 0.0) 247 G4double mass2 = G4NucleiProperties::GetNucl << 208 { return; } 248 << 209 249 wokvi->SetTargetMass(mass2); << 210 G4int iz = G4int(Z); 250 wokvi->SetupKinematic(kinEnergy, currentMate << 211 G4int ia = SelectIsotopeNumber(currentElement); 251 G4double costmin = wokvi->SetupTarget(iz, cu << 212 G4double targetMass = G4NucleiProperties::GetNuclearMass(ia, iz); 252 G4double costmax = (1 == iz && particle == t << 213 wokvi->SetTargetMass(targetMass); 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 214 260 G4ThreeVector newDirection = 215 G4ThreeVector newDirection = 261 wokvi->SampleSingleScattering(costmin, cos << 216 wokvi->SampleSingleScattering(cosTetMinNuc, cosThetaMax, elecRatio); >> 217 G4double cost = newDirection.z(); >> 218 >> 219 G4ThreeVector direction = dp->GetMomentumDirection(); >> 220 newDirection.rotateUz(direction); 262 221 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 222 fParticleChange->ProposeMomentumDirection(newDirection); 288 << 223 289 // recoil << 224 // recoil sampling assuming a small recoil 290 v0 -= v1; << 225 // and first order correction to primary 4-momentum 291 G4double trec = std::max(v0.e() - mass2, 0.0 << 226 G4double mom2 = wokvi->GetMomentumSquare(); >> 227 G4double trec = mom2*(1.0 - cost) >> 228 /(targetMass + (mass + kinEnergy)*(1.0 - cost)); >> 229 >> 230 // the check likely not needed >> 231 if(trec > kinEnergy) { trec = kinEnergy; } >> 232 G4double finalT = kinEnergy - trec; 292 G4double edep = 0.0; 233 G4double edep = 0.0; >> 234 //G4cout<<"G4eCoulombScatteringModel: finalT= "<<finalT<<" Trec= " >> 235 // <<trec << " Z= " << iz << " A= " << ia<<G4endl; 293 236 294 G4double tcut = recoilThreshold; 237 G4double tcut = recoilThreshold; 295 if(pCuts) { tcut= std::max(tcut,(*pCuts)[cur 238 if(pCuts) { tcut= std::max(tcut,(*pCuts)[currentMaterialIndex]); } 296 << 239 297 if(trec > tcut) { 240 if(trec > tcut) { 298 G4ParticleDefinition* ion = theIonTable->G 241 G4ParticleDefinition* ion = theIonTable->GetIon(iz, ia, 0); 299 newDirection = v0.vect().unit(); << 242 G4ThreeVector dir = (direction*sqrt(mom2) - 300 newDirection.rotateUz(dir); << 243 newDirection*sqrt(finalT*(2*mass + finalT))).unit(); 301 auto newdp = new G4DynamicParticle(ion, ne << 244 G4DynamicParticle* newdp = new G4DynamicParticle(ion, dir, trec); 302 fvect->push_back(newdp); 245 fvect->push_back(newdp); 303 } else if(trec > 0.0) { << 246 } else { 304 edep = trec; 247 edep = trec; 305 fParticleChange->ProposeNonIonizingEnergyD 248 fParticleChange->ProposeNonIonizingEnergyDeposit(edep); 306 } 249 } 307 250 308 // finelize primary energy and energy balanc 251 // finelize primary energy and energy balance 309 G4double finalT = v1.e() - mass; << 252 // this threshold may be applied only because for low-enegry 310 if(finalT < 0.0) { << 253 // e+e- msc model is applied 311 edep += finalT; << 254 if(finalT <= lowEnergyThreshold) { >> 255 edep += finalT; 312 finalT = 0.0; 256 finalT = 0.0; 313 } 257 } 314 edep = std::max(edep, 0.0); << 315 fParticleChange->SetProposedKineticEnergy(fi 258 fParticleChange->SetProposedKineticEnergy(finalT); 316 fParticleChange->ProposeLocalEnergyDeposit(e 259 fParticleChange->ProposeLocalEnergyDeposit(edep); 317 } 260 } 318 261 319 //....oooOO0OOooo........oooOO0OOooo........oo 262 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... >> 263 >> 264 320 265