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 // 26 // 27 // ------------------------------------------- 27 // ------------------------------------------------------------------- 28 // 28 // 29 // GEANT4 Class 29 // GEANT4 Class 30 // File name: G4PAIPhotModel.cc 30 // File name: G4PAIPhotModel.cc 31 // 31 // 32 // Author: Vladimir.Grichine@cern.ch on base o 32 // Author: Vladimir.Grichine@cern.ch on base of G4PAIModel MT interface 33 // 33 // 34 // Creation date: 07.10.2013 34 // Creation date: 07.10.2013 35 // 35 // 36 // Modifications: 36 // Modifications: 37 // 37 // 38 // 38 // 39 39 40 #include "G4PAIPhotModel.hh" 40 #include "G4PAIPhotModel.hh" 41 41 42 #include "G4SystemOfUnits.hh" 42 #include "G4SystemOfUnits.hh" 43 #include "G4PhysicalConstants.hh" 43 #include "G4PhysicalConstants.hh" 44 #include "G4Region.hh" 44 #include "G4Region.hh" >> 45 #include "G4PhysicsLogVector.hh" >> 46 #include "G4PhysicsFreeVector.hh" >> 47 #include "G4PhysicsTable.hh" 45 #include "G4ProductionCutsTable.hh" 48 #include "G4ProductionCutsTable.hh" 46 #include "G4MaterialCutsCouple.hh" 49 #include "G4MaterialCutsCouple.hh" 47 #include "G4MaterialTable.hh" 50 #include "G4MaterialTable.hh" >> 51 #include "G4SandiaTable.hh" >> 52 #include "G4OrderedTable.hh" 48 #include "G4RegionStore.hh" 53 #include "G4RegionStore.hh" 49 54 50 #include "Randomize.hh" 55 #include "Randomize.hh" 51 #include "G4Electron.hh" 56 #include "G4Electron.hh" 52 #include "G4Positron.hh" 57 #include "G4Positron.hh" 53 #include "G4Gamma.hh" 58 #include "G4Gamma.hh" 54 #include "G4Poisson.hh" 59 #include "G4Poisson.hh" 55 #include "G4Step.hh" 60 #include "G4Step.hh" 56 #include "G4Material.hh" 61 #include "G4Material.hh" 57 #include "G4DynamicParticle.hh" 62 #include "G4DynamicParticle.hh" 58 #include "G4ParticleDefinition.hh" 63 #include "G4ParticleDefinition.hh" 59 #include "G4ParticleChangeForLoss.hh" 64 #include "G4ParticleChangeForLoss.hh" 60 #include "G4PAIPhotData.hh" 65 #include "G4PAIPhotData.hh" 61 #include "G4DeltaAngle.hh" 66 #include "G4DeltaAngle.hh" 62 67 63 ////////////////////////////////////////////// 68 //////////////////////////////////////////////////////////////////////// 64 69 65 using namespace std; 70 using namespace std; 66 71 67 G4PAIPhotModel::G4PAIPhotModel(const G4Particl 72 G4PAIPhotModel::G4PAIPhotModel(const G4ParticleDefinition* p, const G4String& nam) 68 : G4VEmModel(nam),G4VEmFluctuationModel(nam) 73 : G4VEmModel(nam),G4VEmFluctuationModel(nam), 69 fVerbose(0), 74 fVerbose(0), 70 fModelData(nullptr), 75 fModelData(nullptr), 71 fParticle(nullptr) 76 fParticle(nullptr) 72 { 77 { 73 fElectron = G4Electron::Electron(); 78 fElectron = G4Electron::Electron(); 74 fPositron = G4Positron::Positron(); 79 fPositron = G4Positron::Positron(); 75 80 76 fParticleChange = nullptr; 81 fParticleChange = nullptr; 77 82 78 if(p) { SetParticle(p); } 83 if(p) { SetParticle(p); } 79 else { SetParticle(fElectron); } 84 else { SetParticle(fElectron); } 80 85 81 // default generator 86 // default generator 82 SetAngularDistribution(new G4DeltaAngle()); 87 SetAngularDistribution(new G4DeltaAngle()); 83 fLowestTcut = 12.5*CLHEP::eV; 88 fLowestTcut = 12.5*CLHEP::eV; 84 } 89 } 85 90 86 ////////////////////////////////////////////// 91 //////////////////////////////////////////////////////////////////////////// 87 92 88 G4PAIPhotModel::~G4PAIPhotModel() 93 G4PAIPhotModel::~G4PAIPhotModel() 89 { 94 { >> 95 //G4cout << "G4PAIPhotModel::~G4PAIPhotModel() " << this << G4endl; 90 if(IsMaster()) { delete fModelData; fModelDa 96 if(IsMaster()) { delete fModelData; fModelData = nullptr; } 91 } 97 } 92 98 93 ////////////////////////////////////////////// 99 //////////////////////////////////////////////////////////////////////////// 94 100 95 void G4PAIPhotModel::Initialise(const G4Partic 101 void G4PAIPhotModel::Initialise(const G4ParticleDefinition* p, 96 const G4DataVector& cuts) 102 const G4DataVector& cuts) 97 { 103 { 98 if(fVerbose > 1) << 104 if(fVerbose > 0) 99 { 105 { 100 G4cout<<"G4PAIPhotModel::Initialise for "< 106 G4cout<<"G4PAIPhotModel::Initialise for "<<p->GetParticleName()<<G4endl; 101 } 107 } 102 SetParticle(p); 108 SetParticle(p); 103 fParticleChange = GetParticleChangeForLoss() 109 fParticleChange = GetParticleChangeForLoss(); 104 110 105 if( IsMaster() ) 111 if( IsMaster() ) 106 { 112 { >> 113 >> 114 InitialiseElementSelectors(p, cuts); >> 115 107 delete fModelData; 116 delete fModelData; 108 fMaterialCutsCoupleVector.clear(); 117 fMaterialCutsCoupleVector.clear(); 109 118 110 G4double tmin = LowEnergyLimit()*fRatio; 119 G4double tmin = LowEnergyLimit()*fRatio; 111 G4double tmax = HighEnergyLimit()*fRatio; 120 G4double tmax = HighEnergyLimit()*fRatio; 112 fModelData = new G4PAIPhotData(tmin, tmax, 121 fModelData = new G4PAIPhotData(tmin, tmax, fVerbose); 113 122 114 // Prepare initialization 123 // Prepare initialization 115 const G4MaterialTable* theMaterialTable = 124 const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable(); 116 size_t numOfMat = G4Material::GetNumberO 125 size_t numOfMat = G4Material::GetNumberOfMaterials(); 117 size_t numRegions = fPAIRegionVector.size( 126 size_t numRegions = fPAIRegionVector.size(); 118 127 119 // protect for unit tests 128 // protect for unit tests 120 if(0 == numRegions) { 129 if(0 == numRegions) { 121 G4Exception("G4PAIModel::Initialise()"," 130 G4Exception("G4PAIModel::Initialise()","em0106",JustWarning, 122 "no G4Regions are registered 131 "no G4Regions are registered for the PAI model - World is used"); 123 fPAIRegionVector.push_back(G4RegionStore 132 fPAIRegionVector.push_back(G4RegionStore::GetInstance() 124 ->GetRegion("DefaultRegionForTheWorld 133 ->GetRegion("DefaultRegionForTheWorld", false)); 125 numRegions = 1; 134 numRegions = 1; 126 } 135 } 127 136 128 for( size_t iReg = 0; iReg < numRegions; + 137 for( size_t iReg = 0; iReg < numRegions; ++iReg ) 129 { 138 { 130 const G4Region* curReg = fPAIRegionVecto 139 const G4Region* curReg = fPAIRegionVector[iReg]; 131 G4Region* reg = const_cast<G4Region*>(cu 140 G4Region* reg = const_cast<G4Region*>(curReg); 132 141 133 for(size_t jMat = 0; jMat < numOfMat; ++ 142 for(size_t jMat = 0; jMat < numOfMat; ++jMat) 134 { 143 { 135 G4Material* mat = (*theMaterialTable)[jMat]; 144 G4Material* mat = (*theMaterialTable)[jMat]; 136 const G4MaterialCutsCouple* cutCouple = reg- 145 const G4MaterialCutsCouple* cutCouple = reg->FindCouple(mat); 137 if(nullptr != cutCouple) << 146 //G4cout << "Couple <" << fCutCouple << G4endl; >> 147 if(cutCouple) 138 { 148 { 139 if(fVerbose > 1) << 149 if(fVerbose>0) 140 { 150 { 141 G4cout << "Reg <" <<curReg->GetName() << 151 G4cout << "Reg <" <<curReg->GetName() << "> mat <" 142 << mat->GetName() << "> fCouple= " 152 << mat->GetName() << "> fCouple= " 143 << cutCouple << ", idx= " << cutCouple-> 153 << cutCouple << ", idx= " << cutCouple->GetIndex() 144 <<" " << p->GetParticleName() 154 <<" " << p->GetParticleName() 145 <<", cuts.size() = " << cuts.size() << G 155 <<", cuts.size() = " << cuts.size() << G4endl; 146 } 156 } 147 // check if this couple is not already ini 157 // check if this couple is not already initialized 148 158 149 size_t n = fMaterialCutsCoupleVector.size( 159 size_t n = fMaterialCutsCoupleVector.size(); 150 160 151 G4bool isnew = true; 161 G4bool isnew = true; 152 if( 0 < n ) 162 if( 0 < n ) 153 { 163 { 154 for(size_t i=0; i<fMaterialCutsCoupleVec 164 for(size_t i=0; i<fMaterialCutsCoupleVector.size(); ++i) 155 { 165 { 156 if(cutCouple == fMaterialCutsCoupleVec 166 if(cutCouple == fMaterialCutsCoupleVector[i]) { 157 isnew = false; 167 isnew = false; 158 break; 168 break; 159 } 169 } 160 } 170 } 161 } 171 } 162 // initialise data banks 172 // initialise data banks 163 if(isnew) { 173 if(isnew) { 164 fMaterialCutsCoupleVector.push_back(cutC 174 fMaterialCutsCoupleVector.push_back(cutCouple); 165 G4double deltaCutInKinEnergy = cuts[cutC 175 G4double deltaCutInKinEnergy = cuts[cutCouple->GetIndex()]; 166 fModelData->Initialise(cutCouple, deltaC 176 fModelData->Initialise(cutCouple, deltaCutInKinEnergy, this); 167 } 177 } 168 } 178 } 169 } 179 } 170 } 180 } 171 InitialiseElementSelectors(p, cuts); << 172 } 181 } 173 } 182 } 174 183 175 ////////////////////////////////////////////// 184 ///////////////////////////////////////////////////////////////////////// 176 185 177 void G4PAIPhotModel::InitialiseLocal(const G4P 186 void G4PAIPhotModel::InitialiseLocal(const G4ParticleDefinition*, 178 G4VEmModel* masterModel) 187 G4VEmModel* masterModel) 179 { 188 { 180 fModelData = static_cast<G4PAIPhotModel*>(ma 189 fModelData = static_cast<G4PAIPhotModel*>(masterModel)->GetPAIPhotData(); 181 fMaterialCutsCoupleVector = static_cast<G4PA 190 fMaterialCutsCoupleVector = static_cast<G4PAIPhotModel*>(masterModel)->GetVectorOfCouples(); 182 SetElementSelectors( masterModel->GetElement 191 SetElementSelectors( masterModel->GetElementSelectors() ); 183 } 192 } 184 193 185 ////////////////////////////////////////////// 194 ////////////////////////////////////////////////////////////////////////////// 186 195 187 G4double G4PAIPhotModel::MinEnergyCut(const G4 196 G4double G4PAIPhotModel::MinEnergyCut(const G4ParticleDefinition*, 188 const G4MaterialCutsCouple*) 197 const G4MaterialCutsCouple*) 189 { 198 { 190 return fLowestTcut; 199 return fLowestTcut; 191 } 200 } 192 201 193 ////////////////////////////////////////////// 202 ////////////////////////////////////////////////////////////////////////////// 194 203 195 G4double G4PAIPhotModel::ComputeDEDXPerVolume( 204 G4double G4PAIPhotModel::ComputeDEDXPerVolume(const G4Material*, 196 const G4ParticleDefinition* p, 205 const G4ParticleDefinition* p, 197 G4double kineticEnergy, 206 G4double kineticEnergy, 198 G4double cutEnergy) 207 G4double cutEnergy) 199 { 208 { 200 G4int coupleIndex = FindCoupleIndex(CurrentC 209 G4int coupleIndex = FindCoupleIndex(CurrentCouple()); 201 if(0 > coupleIndex) { return 0.0; } 210 if(0 > coupleIndex) { return 0.0; } 202 211 203 G4double cut = std::min(MaxSecondaryEnergy(p 212 G4double cut = std::min(MaxSecondaryEnergy(p, kineticEnergy), cutEnergy); >> 213 204 G4double scaledTkin = kineticEnergy*fRatio; 214 G4double scaledTkin = kineticEnergy*fRatio; 205 G4double dedx = fChargeSquare*fModelData->DE << 215 206 return dedx; << 216 return fChargeSquare*fModelData->DEDXPerVolume(coupleIndex, scaledTkin, cut); 207 } 217 } 208 218 209 ////////////////////////////////////////////// 219 ///////////////////////////////////////////////////////////////////////// 210 220 211 G4double G4PAIPhotModel::CrossSectionPerVolume 221 G4double G4PAIPhotModel::CrossSectionPerVolume( const G4Material*, 212 const G4ParticleDefinition* p, 222 const G4ParticleDefinition* p, 213 G4double kineticEnergy, 223 G4double kineticEnergy, 214 G4double cutEnergy, 224 G4double cutEnergy, 215 G4double maxEnergy ) 225 G4double maxEnergy ) 216 { 226 { 217 G4int coupleIndex = FindCoupleIndex(CurrentC 227 G4int coupleIndex = FindCoupleIndex(CurrentCouple()); 218 if(0 > coupleIndex) { return 0.0; } << 228 >> 229 if(0 > coupleIndex) return 0.0; 219 230 220 G4double tmax = std::min(MaxSecondaryEnergy( 231 G4double tmax = std::min(MaxSecondaryEnergy(p, kineticEnergy), maxEnergy); 221 if(tmax <= cutEnergy) { return 0.0; } << 232 >> 233 if(tmax <= cutEnergy) return 0.0; 222 234 223 G4double scaledTkin = kineticEnergy*fRatio; 235 G4double scaledTkin = kineticEnergy*fRatio; 224 G4double xs = fChargeSquare*fModelData->Cros << 236 G4double xsc = fChargeSquare*fModelData->CrossSectionPerVolume(coupleIndex, 225 scal << 237 scaledTkin, 226 return xs; << 238 cutEnergy, tmax); >> 239 return xsc; 227 } 240 } 228 241 229 ////////////////////////////////////////////// 242 /////////////////////////////////////////////////////////////////////////// 230 // 243 // 231 // It is analog of PostStepDoIt in terms of se 244 // It is analog of PostStepDoIt in terms of secondary electron. 232 // 245 // 233 246 234 void G4PAIPhotModel::SampleSecondaries(std::ve 247 void G4PAIPhotModel::SampleSecondaries(std::vector<G4DynamicParticle*>* vdp, 235 const G4MaterialCutsCouple* matCC, 248 const G4MaterialCutsCouple* matCC, 236 const G4DynamicParticle* dp, 249 const G4DynamicParticle* dp, 237 G4double tmin, 250 G4double tmin, 238 G4double maxEnergy) 251 G4double maxEnergy) 239 { 252 { 240 G4int coupleIndex = FindCoupleIndex(matCC); 253 G4int coupleIndex = FindCoupleIndex(matCC); 241 if(0 > coupleIndex) { return; } 254 if(0 > coupleIndex) { return; } 242 255 243 SetParticle(dp->GetDefinition()); 256 SetParticle(dp->GetDefinition()); 244 257 245 G4double kineticEnergy = dp->GetKineticEnerg 258 G4double kineticEnergy = dp->GetKineticEnergy(); 246 259 247 G4double tmax = MaxSecondaryEnergy(fParticle 260 G4double tmax = MaxSecondaryEnergy(fParticle, kineticEnergy); 248 261 249 if( maxEnergy < tmax) tmax = maxEnergy; 262 if( maxEnergy < tmax) tmax = maxEnergy; 250 if( tmin >= tmax) return; 263 if( tmin >= tmax) return; 251 264 252 G4ThreeVector direction = dp->GetMomentumDir 265 G4ThreeVector direction = dp->GetMomentumDirection(); 253 G4double scaledTkin = kineticEnergy*fRat 266 G4double scaledTkin = kineticEnergy*fRatio; 254 G4double totalEnergy = kineticEnergy + fM 267 G4double totalEnergy = kineticEnergy + fMass; 255 G4double totalMomentum = sqrt(kineticEnergy 268 G4double totalMomentum = sqrt(kineticEnergy*(totalEnergy + fMass)); 256 G4double plRatio = fModelData->GetPla 269 G4double plRatio = fModelData->GetPlasmonRatio(coupleIndex, scaledTkin); 257 270 258 if( G4UniformRand() <= plRatio ) 271 if( G4UniformRand() <= plRatio ) 259 { 272 { 260 G4double deltaTkin = fModelData->SamplePos 273 G4double deltaTkin = fModelData->SamplePostStepPlasmonTransfer(coupleIndex, scaledTkin); 261 274 262 // G4cout<<"G4PAIPhotModel::SampleSecondar 275 // G4cout<<"G4PAIPhotModel::SampleSecondaries; dp "<<dp->GetParticleDefinition()->GetParticleName() 263 // <<"; Tkin = "<<kineticEnergy/keV<<" keV 276 // <<"; Tkin = "<<kineticEnergy/keV<<" keV; transfer = "<<deltaTkin/keV<<" keV "<<G4endl; 264 277 265 if( deltaTkin <= 0. && fVerbose > 0) 278 if( deltaTkin <= 0. && fVerbose > 0) 266 { 279 { 267 G4cout<<"G4PAIPhotModel::SampleSecondary << 280 G4cout<<"G4PAIPhotModel::SampleSecondary e- deltaTkin = "<<deltaTkin<<G4endl; 268 } 281 } 269 if( deltaTkin <= 0.) { return; } << 282 if( deltaTkin <= 0.) return; 270 283 271 if( deltaTkin > tmax) { deltaTkin = tmax; << 284 if( deltaTkin > tmax) deltaTkin = tmax; 272 285 273 const G4Element* anElement = SelectTargetA 286 const G4Element* anElement = SelectTargetAtom(matCC,fParticle,kineticEnergy, 274 287 dp->GetLogKineticEnergy()); 275 G4int Z = anElement->GetZasInt(); << 288 G4int Z = G4lrint(anElement->GetZ()); 276 289 277 auto deltaRay = new G4DynamicParticle(fEle << 290 G4DynamicParticle* deltaRay = new G4DynamicParticle(fElectron, 278 GetAngularDistribution()->SampleDirectio 291 GetAngularDistribution()->SampleDirection(dp, deltaTkin, 279 Z, matCC->GetMaterial()), 292 Z, matCC->GetMaterial()), 280 deltaTkin); 293 deltaTkin); 281 294 282 // primary change 295 // primary change 283 296 284 kineticEnergy -= deltaTkin; 297 kineticEnergy -= deltaTkin; 285 298 286 if( kineticEnergy <= 0. ) // kill primary 299 if( kineticEnergy <= 0. ) // kill primary as local? energy deposition 287 { 300 { 288 fParticleChange->SetProposedKineticEnerg 301 fParticleChange->SetProposedKineticEnergy(0.0); 289 fParticleChange->ProposeLocalEnergyDepos 302 fParticleChange->ProposeLocalEnergyDeposit(kineticEnergy+deltaTkin); >> 303 // fParticleChange->ProposeTrackStatus(fStopAndKill); 290 return; 304 return; 291 } 305 } 292 else 306 else 293 { 307 { 294 G4ThreeVector dir = totalMomentum*direct 308 G4ThreeVector dir = totalMomentum*direction - deltaRay->GetMomentum(); 295 direction = dir.unit(); 309 direction = dir.unit(); 296 fParticleChange->SetProposedKineticEnerg 310 fParticleChange->SetProposedKineticEnergy(kineticEnergy); 297 fParticleChange->SetProposedMomentumDire 311 fParticleChange->SetProposedMomentumDirection(direction); 298 vdp->push_back(deltaRay); 312 vdp->push_back(deltaRay); 299 } 313 } 300 } 314 } 301 else // secondary X-ray CR photon 315 else // secondary X-ray CR photon 302 { 316 { 303 G4double deltaTkin = fModelData->Sampl 317 G4double deltaTkin = fModelData->SamplePostStepPhotonTransfer(coupleIndex, scaledTkin); 304 318 305 // G4cout<<"PAIPhotonModel PhotonPostStep 319 // G4cout<<"PAIPhotonModel PhotonPostStepTransfer = "<<deltaTkin/keV<<" keV"<<G4endl; 306 320 307 if( deltaTkin <= 0. ) 321 if( deltaTkin <= 0. ) 308 { 322 { 309 G4cout<<"G4PAIPhotonModel::SampleSeconda 323 G4cout<<"G4PAIPhotonModel::SampleSecondary gamma deltaTkin = "<<deltaTkin<<G4endl; 310 } 324 } 311 if( deltaTkin <= 0.) return; 325 if( deltaTkin <= 0.) return; 312 326 313 if( deltaTkin >= kineticEnergy ) // stop p 327 if( deltaTkin >= kineticEnergy ) // stop primary 314 { 328 { 315 deltaTkin = kineticEnergy; 329 deltaTkin = kineticEnergy; 316 kineticEnergy = 0.0; 330 kineticEnergy = 0.0; 317 } 331 } 318 G4double costheta = 0.; // G4UniformRand() 332 G4double costheta = 0.; // G4UniformRand(); // VG: ??? for start only 319 G4double sintheta = sqrt((1.+costheta)*(1. 333 G4double sintheta = sqrt((1.+costheta)*(1.-costheta)); 320 334 321 // direction of the 'Cherenkov' photon 335 // direction of the 'Cherenkov' photon 322 G4double phi = twopi*G4UniformRand(); 336 G4double phi = twopi*G4UniformRand(); 323 G4double dirx = sintheta*cos(phi), diry = 337 G4double dirx = sintheta*cos(phi), diry = sintheta*sin(phi), dirz = costheta; 324 338 325 G4ThreeVector deltaDirection(dirx,diry,dir 339 G4ThreeVector deltaDirection(dirx,diry,dirz); 326 deltaDirection.rotateUz(direction); 340 deltaDirection.rotateUz(direction); 327 341 328 if( kineticEnergy > 0.) // primary change 342 if( kineticEnergy > 0.) // primary change 329 { 343 { 330 kineticEnergy -= deltaTkin; 344 kineticEnergy -= deltaTkin; 331 fParticleChange->SetProposedKineticEnerg 345 fParticleChange->SetProposedKineticEnergy(kineticEnergy); 332 } 346 } 333 else // stop primary, but pass X-ray CR 347 else // stop primary, but pass X-ray CR 334 { 348 { 335 // fParticleChange->ProposeLocalEnergyDe 349 // fParticleChange->ProposeLocalEnergyDeposit(deltaTkin); 336 fParticleChange->SetProposedKineticEnerg 350 fParticleChange->SetProposedKineticEnergy(0.0); 337 } 351 } 338 // create G4DynamicParticle object for pho 352 // create G4DynamicParticle object for photon ray 339 353 340 auto photonRay = new G4DynamicParticle; << 354 G4DynamicParticle* photonRay = new G4DynamicParticle; 341 photonRay->SetDefinition( G4Gamma::Gamma() 355 photonRay->SetDefinition( G4Gamma::Gamma() ); 342 photonRay->SetKineticEnergy( deltaTkin ); 356 photonRay->SetKineticEnergy( deltaTkin ); 343 photonRay->SetMomentumDirection(deltaDirec 357 photonRay->SetMomentumDirection(deltaDirection); 344 358 345 vdp->push_back(photonRay); 359 vdp->push_back(photonRay); 346 } 360 } 347 return; 361 return; 348 } 362 } 349 363 350 ////////////////////////////////////////////// 364 /////////////////////////////////////////////////////////////////////// 351 365 352 G4double G4PAIPhotModel::SampleFluctuations( 366 G4double G4PAIPhotModel::SampleFluctuations( 353 const G4MaterialCutsC 367 const G4MaterialCutsCouple* matCC, 354 const G4DynamicPartic 368 const G4DynamicParticle* aParticle, 355 const G4double, const 369 const G4double, const G4double, 356 const G4double step, 370 const G4double step, const G4double eloss) 357 { 371 { 358 // return 0.; 372 // return 0.; 359 G4int coupleIndex = FindCoupleIndex(matCC); 373 G4int coupleIndex = FindCoupleIndex(matCC); 360 if(0 > coupleIndex) { return eloss; } 374 if(0 > coupleIndex) { return eloss; } 361 375 362 SetParticle(aParticle->GetDefinition()); 376 SetParticle(aParticle->GetDefinition()); 363 377 >> 378 364 // G4cout << "G4PAIPhotModel::SampleFluctuat 379 // G4cout << "G4PAIPhotModel::SampleFluctuations step(mm)= "<< step/mm 365 // << " Eloss(keV)= " << eloss/keV << " in 380 // << " Eloss(keV)= " << eloss/keV << " in " 366 // << matCC->GetMaterial()->GetName() << G4e 381 // << matCC->GetMaterial()->GetName() << G4endl; >> 382 367 383 368 G4double Tkin = aParticle->GetKineticE 384 G4double Tkin = aParticle->GetKineticEnergy(); 369 G4double scaledTkin = Tkin*fRatio; 385 G4double scaledTkin = Tkin*fRatio; 370 386 371 G4double loss = fModelData->SampleAlongStepP 387 G4double loss = fModelData->SampleAlongStepPhotonTransfer(coupleIndex, Tkin, 372 scaledTkin, << 388 scaledTkin, 373 step*fChargeSquare); << 389 step*fChargeSquare); 374 loss += fModelData->SampleAlongStepPlasmonTr << 390 loss += fModelData->SampleAlongStepPlasmonTransfer(coupleIndex, Tkin, 375 << 391 scaledTkin, >> 392 step*fChargeSquare); >> 393 376 394 377 // G4cout<<" PAIPhotModel::SampleFluctuatio 395 // G4cout<<" PAIPhotModel::SampleFluctuations loss = "<<loss/keV<<" keV, on step = " 378 // <<step/mm<<" mm"<<G4endl; 396 // <<step/mm<<" mm"<<G4endl; 379 return loss; 397 return loss; 380 398 381 } 399 } 382 400 383 ////////////////////////////////////////////// 401 ////////////////////////////////////////////////////////////////////// 384 // 402 // 385 // Returns the statistical estimation of the e 403 // Returns the statistical estimation of the energy loss distribution variance 386 // 404 // 387 405 388 406 389 G4double G4PAIPhotModel::Dispersion(const G4Ma 407 G4double G4PAIPhotModel::Dispersion(const G4Material* material, 390 const G4Dy 408 const G4DynamicParticle* aParticle, 391 const G4double tcut, 409 const G4double tcut, 392 const G4double tmax, 410 const G4double tmax, 393 const G4double step) 411 const G4double step) 394 { 412 { 395 G4double particleMass = aParticle->GetMass( 413 G4double particleMass = aParticle->GetMass(); 396 G4double electronDensity = material->GetElec 414 G4double electronDensity = material->GetElectronDensity(); 397 G4double kineticEnergy = aParticle->GetKinet 415 G4double kineticEnergy = aParticle->GetKineticEnergy(); 398 G4double q = aParticle->GetCharge()/eplus; 416 G4double q = aParticle->GetCharge()/eplus; 399 G4double etot = kineticEnergy + particleMass 417 G4double etot = kineticEnergy + particleMass; 400 G4double beta2 = kineticEnergy*(kineticEnerg 418 G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*particleMass)/(etot*etot); 401 G4double siga = (tmax/beta2 - 0.5*tcut) * t 419 G4double siga = (tmax/beta2 - 0.5*tcut) * twopi_mc2_rcl2 * step 402 * electronDensity * q * q; 420 * electronDensity * q * q; 403 421 404 return siga; 422 return siga; 405 } 423 } 406 424 407 ////////////////////////////////////////////// 425 ///////////////////////////////////////////////////////////////////// 408 426 409 G4double G4PAIPhotModel::MaxSecondaryEnergy( c 427 G4double G4PAIPhotModel::MaxSecondaryEnergy( const G4ParticleDefinition* p, 410 G4double kinEnergy) 428 G4double kinEnergy) 411 { 429 { 412 SetParticle(p); 430 SetParticle(p); 413 G4double tmax = kinEnergy; 431 G4double tmax = kinEnergy; 414 if(p == fElectron) { tmax *= 0.5; } 432 if(p == fElectron) { tmax *= 0.5; } 415 else if(p != fPositron) { 433 else if(p != fPositron) { 416 G4double ratio= electron_mass_c2/fMass; 434 G4double ratio= electron_mass_c2/fMass; 417 G4double gamma= kinEnergy/fMass + 1.0; 435 G4double gamma= kinEnergy/fMass + 1.0; 418 tmax = 2.0*electron_mass_c2*(gamma*gamma - 436 tmax = 2.0*electron_mass_c2*(gamma*gamma - 1.) / 419 (1. + 2.0*gamma*ratio + rati 437 (1. + 2.0*gamma*ratio + ratio*ratio); 420 } 438 } 421 return tmax; 439 return tmax; 422 } 440 } 423 441 424 ////////////////////////////////////////////// 442 /////////////////////////////////////////////////////////////// 425 443 426 void G4PAIPhotModel::DefineForRegion(const G4R 444 void G4PAIPhotModel::DefineForRegion(const G4Region* r) 427 { 445 { 428 fPAIRegionVector.push_back(r); 446 fPAIRegionVector.push_back(r); 429 } 447 } 430 448 431 // 449 // 432 // 450 // 433 ////////////////////////////////////////////// 451 ///////////////////////////////////////////////// 434 452