Changeset 961 for trunk/source/processes/electromagnetic/xrays/src
- Timestamp:
- Apr 6, 2009, 12:21:12 PM (15 years ago)
- Location:
- trunk/source/processes/electromagnetic/xrays/src
- Files:
-
- 12 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/source/processes/electromagnetic/xrays/src/G4Cerenkov.cc
r819 r961 25 25 // 26 26 // 27 // $Id: G4Cerenkov.cc,v 1.2 3 2007/10/15 20:05:23gum Exp $28 // GEANT4 tag $Name: $27 // $Id: G4Cerenkov.cc,v 1.26 2008/11/14 20:16:51 gum Exp $ 28 // GEANT4 tag $Name: geant4-09-02-ref-02 $ 29 29 // 30 30 //////////////////////////////////////////////////////////////////////// … … 65 65 #include "G4ios.hh" 66 66 #include "G4Poisson.hh" 67 #include "G4EmProcessSubType.hh" 68 69 #include "G4LossTableManager.hh" 70 #include "G4MaterialCutsCouple.hh" 71 #include "G4ParticleDefinition.hh" 72 67 73 #include "G4Cerenkov.hh" 68 74 … … 86 92 87 93 G4Cerenkov::G4Cerenkov(const G4String& processName, G4ProcessType type) 88 : G4V DiscreteProcess(processName, type)94 : G4VProcess(processName, type) 89 95 { 90 96 G4cout << "G4Cerenkov::G4Cerenkov constructor" << G4endl; 91 G4cout << "NOTE: this is now a G4V DiscreteProcess!" << G4endl;97 G4cout << "NOTE: this is now a G4VProcess!" << G4endl; 92 98 G4cout << "Required change in UserPhysicsList: " << G4endl; 93 99 G4cout << "change: pmanager->AddContinuousProcess(theCerenkovProcess);" << G4endl; … … 95 101 G4cout << " pmanager->SetProcessOrdering(theCerenkovProcess,idxPostStep);" << G4endl; 96 102 103 SetProcessSubType(fCerenkov); 104 97 105 fTrackSecondariesFirst = false; 106 fMaxBetaChange = 0.; 98 107 fMaxPhotons = 0; 99 108 … … 141 150 142 151 { 143 144 152 ////////////////////////////////////////////////////// 145 153 // Should we ensure that the material is dispersive? … … 160 168 G4MaterialPropertiesTable* aMaterialPropertiesTable = 161 169 aMaterial->GetMaterialPropertiesTable(); 162 if (!aMaterialPropertiesTable) 163 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep); 170 if (!aMaterialPropertiesTable) return pParticleChange; 164 171 165 172 const G4MaterialPropertyVector* Rindex = 166 173 aMaterialPropertiesTable->GetProperty("RINDEX"); 167 if (!Rindex) 168 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep); 174 if (!Rindex) return pParticleChange; 169 175 170 176 // particle charge … … 184 190 aParticleChange.SetNumberOfSecondaries(0); 185 191 186 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);192 return pParticleChange; 187 193 188 194 } … … 201 207 aParticleChange.SetNumberOfSecondaries(0); 202 208 203 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);209 return pParticleChange; 204 210 } 205 211 … … 215 221 //////////////////////////////////////////////////////////////// 216 222 217 G4double Pmin = Rindex->GetMinPhoton Momentum();218 G4double Pmax = Rindex->GetMaxPhoton Momentum();223 G4double Pmin = Rindex->GetMinPhotonEnergy(); 224 G4double Pmax = Rindex->GetMaxPhotonEnergy(); 219 225 G4double dp = Pmax - Pmin; 220 226 … … 226 232 G4double maxSin2 = (1.0 - maxCos) * (1.0 + maxCos); 227 233 234 const G4double beta1 = pPreStepPoint ->GetBeta(); 235 const G4double beta2 = pPostStepPoint->GetBeta(); 236 237 G4double MeanNumberOfPhotons1 = 238 GetAverageNumberOfPhotons(charge,beta1,aMaterial,Rindex); 239 G4double MeanNumberOfPhotons2 = 240 GetAverageNumberOfPhotons(charge,beta2,aMaterial,Rindex); 241 228 242 for (G4int i = 0; i < NumPhotons; i++) { 229 243 230 // Determine photon momentum244 // Determine photon energy 231 245 232 246 G4double rand; 233 G4double sampled Momentum, sampledRI;247 G4double sampledEnergy, sampledRI; 234 248 G4double cosTheta, sin2Theta; 235 249 236 // sample a momentum250 // sample an energy 237 251 238 252 do { 239 253 rand = G4UniformRand(); 240 sampled Momentum= Pmin + rand * dp;241 sampledRI = Rindex->GetProperty(sampled Momentum);254 sampledEnergy = Pmin + rand * dp; 255 sampledRI = Rindex->GetProperty(sampledEnergy); 242 256 cosTheta = BetaInverse / sampledRI; 243 257 … … 256 270 G4double cosPhi = cos(phi); 257 271 258 // calculate x,y, and z components of photon momentum272 // calculate x,y, and z components of photon energy 259 273 // (in coord system with primary particle direction 260 274 // aligned with the z axis) … … 299 313 photonPolarization.z()); 300 314 301 aCerenkovPhoton->SetKineticEnergy(sampled Momentum);315 aCerenkovPhoton->SetKineticEnergy(sampledEnergy); 302 316 303 317 // Generate new G4Track object: 304 318 305 rand = G4UniformRand(); 306 307 G4double delta = rand * aStep.GetStepLength(); 319 G4double delta, NumberOfPhotons, N; 320 321 do { 322 rand = G4UniformRand(); 323 delta = rand * aStep.GetStepLength(); 324 NumberOfPhotons = MeanNumberOfPhotons1 - delta * 325 (MeanNumberOfPhotons1-MeanNumberOfPhotons2)/ 326 aStep.GetStepLength(); 327 N = G4UniformRand() * 328 std::max(MeanNumberOfPhotons1,MeanNumberOfPhotons2); 329 } while (N > NumberOfPhotons); 330 308 331 G4double deltaTime = delta / 309 332 ((pPreStepPoint->GetVelocity()+ … … 318 341 new G4Track(aCerenkovPhoton,aSecondaryTime,aSecondaryPosition); 319 342 320 aSecondaryTrack->SetTouchableHandle((G4VTouchable*)0); 343 aSecondaryTrack->SetTouchableHandle( 344 aStep.GetPreStepPoint()->GetTouchableHandle()); 321 345 322 346 aSecondaryTrack->SetParentID(aTrack.GetTrackID()); … … 330 354 } 331 355 332 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);356 return pParticleChange; 333 357 } 334 358 … … 372 396 373 397 // Retrieve the first refraction index in vector 374 // of (photon momentum, refraction index) pairs398 // of (photon energy, refraction index) pairs 375 399 376 400 theRefractionIndexVector->ResetIterator(); … … 382 406 if (currentRI > 1.0) { 383 407 384 // Create first (photon momentum, Cerenkov Integral)408 // Create first (photon energy, Cerenkov Integral) 385 409 // pair 386 410 387 411 G4double currentPM = theRefractionIndexVector-> 388 GetPhoton Momentum();412 GetPhotonEnergy(); 389 413 G4double currentCAI = 0.0; 390 414 … … 398 422 G4double prevRI = currentRI; 399 423 400 // loop over all (photon momentum, refraction index)424 // loop over all (photon energy, refraction index) 401 425 // pairs stored for this material 402 426 … … 407 431 408 432 currentPM = theRefractionIndexVector-> 409 GetPhoton Momentum();433 GetPhotonEnergy(); 410 434 411 435 currentCAI = 0.5*(1.0/(prevRI*prevRI) + … … 441 465 // 442 466 443 G4double G4Cerenkov::GetMeanFreePath(const G4Track& aTrack, 467 G4double G4Cerenkov::GetMeanFreePath(const G4Track&, 468 G4double, 469 G4ForceCondition*) 470 { 471 return 1.; 472 } 473 474 G4double G4Cerenkov::PostStepGetPhysicalInteractionLength( 475 const G4Track& aTrack, 444 476 G4double, 445 477 G4ForceCondition* condition) 446 478 { 447 *condition = StronglyForced; 448 449 // If user has defined an average maximum number of photons to 450 // be generated in a Step, then return the Step length for that 451 // number of photons. 452 453 if (fMaxPhotons <= 0) return DBL_MAX; 479 *condition = NotForced; 480 G4double StepLimit = DBL_MAX; 454 481 455 482 const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle(); 456 483 const G4Material* aMaterial = aTrack.GetMaterial(); 457 458 G4MaterialPropertiesTable* aMaterialPropertiesTable = 459 aMaterial->GetMaterialPropertiesTable(); 460 if (!aMaterialPropertiesTable) return DBL_MAX; 461 462 const G4MaterialPropertyVector* Rindex = 463 aMaterialPropertiesTable->GetProperty("RINDEX"); 464 if (!Rindex) return DBL_MAX; 465 466 // particle charge 467 const G4double charge = aParticle->GetDefinition()->GetPDGCharge(); 484 const G4MaterialCutsCouple* couple = aTrack.GetMaterialCutsCouple(); 485 486 const G4double kineticEnergy = aParticle->GetKineticEnergy(); 487 const G4ParticleDefinition* particleType = aParticle->GetDefinition(); 488 const G4double mass = particleType->GetPDGMass(); 468 489 469 490 // particle beta 470 491 const G4double beta = aParticle->GetTotalMomentum() / 471 492 aParticle->GetTotalEnergy(); 472 473 G4double MeanNumberOfPhotons = 474 GetAverageNumberOfPhotons(charge,beta,aMaterial,Rindex); 475 476 if(MeanNumberOfPhotons <= 0.0) return DBL_MAX; 477 478 G4double StepLimit = fMaxPhotons / MeanNumberOfPhotons; 479 480 return StepLimit; 493 // particle gamma 494 const G4double gamma = 1./std::sqrt(1.-beta*beta); 495 496 G4MaterialPropertiesTable* aMaterialPropertiesTable = 497 aMaterial->GetMaterialPropertiesTable(); 498 499 const G4MaterialPropertyVector* Rindex = NULL; 500 501 if (aMaterialPropertiesTable) 502 Rindex = aMaterialPropertiesTable->GetProperty("RINDEX"); 503 504 G4double nMax; 505 if (Rindex) { 506 nMax = Rindex->GetMaxProperty(); 507 } else { 508 return StepLimit; 509 } 510 511 G4double BetaMin = 1./nMax; 512 if ( BetaMin >= 1. ) return StepLimit; 513 514 G4double GammaMin = 1./std::sqrt(1.-BetaMin*BetaMin); 515 516 if (gamma < GammaMin ) return StepLimit; 517 518 G4double kinEmin = mass*(GammaMin-1.); 519 520 G4double RangeMin = G4LossTableManager::Instance()-> 521 GetRange(particleType, 522 kinEmin, 523 couple); 524 G4double Range = G4LossTableManager::Instance()-> 525 GetRange(particleType, 526 kineticEnergy, 527 couple); 528 529 G4double Step = Range - RangeMin; 530 if (Step < 1.*um ) return StepLimit; 531 532 if (Step > 0. && Step < StepLimit) StepLimit = Step; 533 534 // If user has defined an average maximum number of photons to 535 // be generated in a Step, then calculate the Step length for 536 // that number of photons. 537 538 if (fMaxPhotons > 0) { 539 540 // particle charge 541 const G4double charge = aParticle-> 542 GetDefinition()->GetPDGCharge(); 543 544 G4double MeanNumberOfPhotons = 545 GetAverageNumberOfPhotons(charge,beta,aMaterial,Rindex); 546 547 G4double Step = 0.; 548 if (MeanNumberOfPhotons > 0.0) Step = fMaxPhotons / 549 MeanNumberOfPhotons; 550 551 if (Step > 0. && Step < StepLimit) StepLimit = Step; 552 } 553 554 // If user has defined an maximum allowed change in beta per step 555 if (fMaxBetaChange > 0.) { 556 557 G4double dedx = G4LossTableManager::Instance()-> 558 GetDEDX(particleType, 559 kineticEnergy, 560 couple); 561 562 G4double deltaGamma = gamma - 563 1./std::sqrt(1.-beta*beta* 564 (1.-fMaxBetaChange)* 565 (1.-fMaxBetaChange)); 566 567 G4double Step = mass * deltaGamma / dedx; 568 569 if (Step > 0. && Step < StepLimit) StepLimit = Step; 570 571 } 572 573 *condition = StronglyForced; 574 return StepLimit; 481 575 } 482 576 … … 512 606 if(!(CerenkovAngleIntegrals->IsFilledVectorExist()))return 0.0; 513 607 514 // Min and Max photon momenta515 G4double Pmin = Rindex->GetMinPhoton Momentum();516 G4double Pmax = Rindex->GetMaxPhoton Momentum();608 // Min and Max photon energies 609 G4double Pmin = Rindex->GetMinPhotonEnergy(); 610 G4double Pmax = Rindex->GetMaxPhotonEnergy(); 517 611 518 612 // Min and Max Refraction Indices … … 541 635 // If n(Pmin) < 1/Beta, and n(Pmax) >= 1/Beta, then 542 636 // we need to find a P such that the value of n(P) == 1/Beta. 543 // Interpolation is performed by the GetPhoton Momentum() and637 // Interpolation is performed by the GetPhotonEnergy() and 544 638 // GetProperty() methods of the G4MaterialPropertiesTable and 545 639 // the GetValue() method of G4PhysicsVector. 546 640 547 641 else { 548 Pmin = Rindex->GetPhoton Momentum(BetaInverse);642 Pmin = Rindex->GetPhotonEnergy(BetaInverse); 549 643 dp = Pmax - Pmin; 550 644 -
trunk/source/processes/electromagnetic/xrays/src/G4ForwardXrayTR.cc
r819 r961 26 26 // 27 27 // $Id: G4ForwardXrayTR.cc,v 1.14 2007/05/11 14:23:04 gcosmo Exp $ 28 // GEANT4 tag $Name: $28 // GEANT4 tag $Name: geant4-09-02-ref-02 $ 29 29 // 30 30 // G4ForwardXrayTR class -- implementation file -
trunk/source/processes/electromagnetic/xrays/src/G4GammaXTRadiator.cc
r819 r961 26 26 // 27 27 // $Id: G4GammaXTRadiator.cc,v 1.5 2006/06/29 19:56:07 gunter Exp $ 28 // GEANT4 tag $Name: $28 // GEANT4 tag $Name: geant4-09-02-ref-02 $ 29 29 // 30 30 -
trunk/source/processes/electromagnetic/xrays/src/G4RegularXTRadiator.cc
r819 r961 26 26 // 27 27 // $Id: G4RegularXTRadiator.cc,v 1.9 2006/06/29 19:56:09 gunter Exp $ 28 // GEANT4 tag $Name: $28 // GEANT4 tag $Name: geant4-09-02-ref-02 $ 29 29 // 30 30 -
trunk/source/processes/electromagnetic/xrays/src/G4Scintillation.cc
r819 r961 25 25 // 26 26 // 27 // $Id: G4Scintillation.cc,v 1. 26 2006/06/29 19:56:11 gunterExp $28 // GEANT4 tag $Name: geant4-09-0 1-patch-02 $27 // $Id: G4Scintillation.cc,v 1.30 2008/10/22 01:19:11 gum Exp $ 28 // GEANT4 tag $Name: geant4-09-02-ref-02 $ 29 29 // 30 30 //////////////////////////////////////////////////////////////////////// … … 64 64 65 65 #include "G4ios.hh" 66 #include "G4EmProcessSubType.hh" 67 66 68 #include "G4Scintillation.hh" 67 69 … … 88 90 : G4VRestDiscreteProcess(processName, type) 89 91 { 92 SetProcessSubType(fScintillation); 93 90 94 fTrackSecondariesFirst = false; 91 95 … … 101 105 102 106 BuildThePhysicsTable(); 107 108 emSaturation = NULL; 103 109 } 104 110 … … 180 186 G4double ScintillationYield = aMaterialPropertiesTable-> 181 187 GetConstProperty("SCINTILLATIONYIELD"); 188 ScintillationYield *= YieldFactor; 189 182 190 G4double ResolutionScale = aMaterialPropertiesTable-> 183 191 GetConstProperty("RESOLUTIONSCALE"); 184 192 185 ScintillationYield = YieldFactor * ScintillationYield; 186 187 G4double MeanNumberOfPhotons = ScintillationYield * TotalEnergyDeposit; 193 // Birks law saturation: 194 195 G4double constBirks = 0.0; 196 197 constBirks = aMaterial->GetIonisation()->GetBirksConstant(); 198 199 G4double MeanNumberOfPhotons; 200 201 if (emSaturation) { 202 MeanNumberOfPhotons = ScintillationYield* 203 (emSaturation->VisibleEnergyDeposition(&aStep)); 204 } else { 205 MeanNumberOfPhotons = ScintillationYield*TotalEnergyDeposit; 206 } 188 207 189 208 G4int NumPhotons; 190 if (MeanNumberOfPhotons > 10.) { 209 210 if (MeanNumberOfPhotons > 10.) 211 { 191 212 G4double sigma = ResolutionScale * sqrt(MeanNumberOfPhotons); 192 213 NumPhotons = G4int(G4RandGauss::shoot(MeanNumberOfPhotons,sigma)+0.5); 193 214 } 194 else { 215 else 216 { 195 217 NumPhotons = G4int(G4Poisson(MeanNumberOfPhotons)); 196 218 } 197 219 198 if (NumPhotons <= 0) {199 220 if (NumPhotons <= 0) 221 { 200 222 // return unchanged particle and no secondaries 201 223 … … 274 296 for (G4int i = 0; i < Num; i++) { 275 297 276 // Determine photon momentum298 // Determine photon energy 277 299 278 300 G4double CIIvalue = G4UniformRand()*CIImax; 279 G4double sampled Momentum=301 G4double sampledEnergy = 280 302 ScintillationIntegral->GetEnergy(CIIvalue); 281 303 282 304 if (verboseLevel>1) { 283 G4cout << "sampled Momentum = " << sampledMomentum<< G4endl;305 G4cout << "sampledEnergy = " << sampledEnergy << G4endl; 284 306 G4cout << "CIIvalue = " << CIIvalue << G4endl; 285 307 } … … 330 352 photonPolarization.z()); 331 353 332 aScintillationPhoton->SetKineticEnergy(sampled Momentum);354 aScintillationPhoton->SetKineticEnergy(sampledEnergy); 333 355 334 356 // Generate new G4Track object: … … 358 380 new G4Track(aScintillationPhoton,aSecondaryTime,aSecondaryPosition); 359 381 360 aSecondaryTrack->SetTouchableHandle((G4VTouchable*)0); 382 aSecondaryTrack->SetTouchableHandle( 383 aStep.GetPreStepPoint()->GetTouchableHandle()); 384 // aSecondaryTrack->SetTouchableHandle((G4VTouchable*)0); 361 385 362 386 aSecondaryTrack->SetParentID(aTrack.GetTrackID()); … … 417 441 418 442 // Retrieve the first intensity point in vector 419 // of (photon momentum, intensity) pairs443 // of (photon energy, intensity) pairs 420 444 421 445 theFastLightVector->ResetIterator(); … … 427 451 if (currentIN >= 0.0) { 428 452 429 // Create first (photon momentum, Scintillation453 // Create first (photon energy, Scintillation 430 454 // Integral pair 431 455 432 456 G4double currentPM = theFastLightVector-> 433 GetPhoton Momentum();457 GetPhotonEnergy(); 434 458 435 459 G4double currentCII = 0.0; … … 444 468 G4double prevIN = currentIN; 445 469 446 // loop over all (photon momentum, intensity)470 // loop over all (photon energy, intensity) 447 471 // pairs stored for this material 448 472 … … 450 474 { 451 475 currentPM = theFastLightVector-> 452 GetPhoton Momentum();476 GetPhotonEnergy(); 453 477 454 478 currentIN=theFastLightVector-> … … 477 501 478 502 // Retrieve the first intensity point in vector 479 // of (photon momentum, intensity) pairs503 // of (photon energy, intensity) pairs 480 504 481 505 theSlowLightVector->ResetIterator(); … … 487 511 if (currentIN >= 0.0) { 488 512 489 // Create first (photon momentum, Scintillation513 // Create first (photon energy, Scintillation 490 514 // Integral pair 491 515 492 516 G4double currentPM = theSlowLightVector-> 493 GetPhoton Momentum();517 GetPhotonEnergy(); 494 518 495 519 G4double currentCII = 0.0; … … 504 528 G4double prevIN = currentIN; 505 529 506 // loop over all (photon momentum, intensity)530 // loop over all (photon energy, intensity) 507 531 // pairs stored for this material 508 532 … … 510 534 { 511 535 currentPM = theSlowLightVector-> 512 GetPhoton Momentum();536 GetPhotonEnergy(); 513 537 514 538 currentIN=theSlowLightVector-> -
trunk/source/processes/electromagnetic/xrays/src/G4StrawTubeXTRadiator.cc
r819 r961 26 26 // 27 27 // $Id: G4StrawTubeXTRadiator.cc,v 1.6 2007/09/29 17:49:34 vnivanch Exp $ 28 // GEANT4 tag $Name: $28 // GEANT4 tag $Name: geant4-09-02-ref-02 $ 29 29 // 30 30 -
trunk/source/processes/electromagnetic/xrays/src/G4SynchrotronRadiation.cc
r819 r961 26 26 // 27 27 // $Id: G4SynchrotronRadiation.cc,v 1.5 2006/06/29 19:56:15 gunter Exp $ 28 // GEANT4 tag $Name: $28 // GEANT4 tag $Name: geant4-09-02-ref-02 $ 29 29 // 30 30 // -------------------------------------------------------------- -
trunk/source/processes/electromagnetic/xrays/src/G4SynchrotronRadiationInMat.cc
r819 r961 26 26 // 27 27 // $Id: G4SynchrotronRadiationInMat.cc,v 1.2 2006/06/29 19:56:17 gunter Exp $ 28 // GEANT4 tag $Name: $28 // GEANT4 tag $Name: geant4-09-02-ref-02 $ 29 29 // 30 30 // -------------------------------------------------------------- -
trunk/source/processes/electromagnetic/xrays/src/G4TransitionRadiation.cc
r819 r961 26 26 // 27 27 // $Id: G4TransitionRadiation.cc,v 1.7 2006/06/29 19:56:19 gunter Exp $ 28 // GEANT4 tag $Name: $28 // GEANT4 tag $Name: geant4-09-02-ref-02 $ 29 29 // 30 30 // G4TransitionRadiation class -- implementation file -
trunk/source/processes/electromagnetic/xrays/src/G4TransparentRegXTRadiator.cc
r819 r961 26 26 // 27 27 // $Id: G4TransparentRegXTRadiator.cc,v 1.11 2007/09/29 17:49:34 vnivanch Exp $ 28 // GEANT4 tag $Name: $28 // GEANT4 tag $Name: geant4-09-02-ref-02 $ 29 29 // 30 30 -
trunk/source/processes/electromagnetic/xrays/src/G4VTransitionRadiation.cc
r819 r961 26 26 // 27 27 // $Id: G4VTransitionRadiation.cc,v 1.5 2006/06/29 19:56:23 gunter Exp $ 28 // GEANT4 tag $Name: $28 // GEANT4 tag $Name: geant4-09-02-ref-02 $ 29 29 // 30 30 // G4VTransitionRadiation class -- implementation file -
trunk/source/processes/electromagnetic/xrays/src/G4VXTRenergyLoss.cc
r819 r961 26 26 // 27 27 // $Id: G4VXTRenergyLoss.cc,v 1.44 2007/09/29 17:49:34 vnivanch Exp $ 28 // GEANT4 tag $Name: $28 // GEANT4 tag $Name: geant4-09-02-ref-02 $ 29 29 // 30 30 // History:
Note: See TracChangeset
for help on using the changeset viewer.