00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00030
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00063
00064 #include "G4ios.hh"
00065 #include "G4PhysicalConstants.hh"
00066 #include "G4SystemOfUnits.hh"
00067 #include "G4Poisson.hh"
00068 #include "G4EmProcessSubType.hh"
00069
00070 #include "G4LossTableManager.hh"
00071 #include "G4MaterialCutsCouple.hh"
00072 #include "G4ParticleDefinition.hh"
00073
00074 #include "G4Cerenkov.hh"
00075
00077
00079
00081
00083
00084
00085
00086
00087
00089
00091
00092 G4Cerenkov::G4Cerenkov(const G4String& processName, G4ProcessType type)
00093 : G4VProcess(processName, type)
00094 {
00095 SetProcessSubType(fCerenkov);
00096
00097 fTrackSecondariesFirst = false;
00098 fMaxBetaChange = 0.;
00099 fMaxPhotons = 0;
00100
00101 thePhysicsTable = NULL;
00102
00103 if (verboseLevel>0) {
00104 G4cout << GetProcessName() << " is created " << G4endl;
00105 }
00106
00107 BuildThePhysicsTable();
00108 }
00109
00110
00111
00112
00113
00115
00117
00118 G4Cerenkov::~G4Cerenkov()
00119 {
00120 if (thePhysicsTable != NULL) {
00121 thePhysicsTable->clearAndDestroy();
00122 delete thePhysicsTable;
00123 }
00124 }
00125
00127
00129
00130
00131
00132
00133 G4VParticleChange*
00134 G4Cerenkov::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
00135
00136
00137
00138
00139
00140
00141
00142
00143 {
00145
00147
00148 aParticleChange.Initialize(aTrack);
00149
00150 const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
00151 const G4Material* aMaterial = aTrack.GetMaterial();
00152
00153 G4StepPoint* pPreStepPoint = aStep.GetPreStepPoint();
00154 G4StepPoint* pPostStepPoint = aStep.GetPostStepPoint();
00155
00156 G4ThreeVector x0 = pPreStepPoint->GetPosition();
00157 G4ThreeVector p0 = aStep.GetDeltaPosition().unit();
00158 G4double t0 = pPreStepPoint->GetGlobalTime();
00159
00160 G4MaterialPropertiesTable* aMaterialPropertiesTable =
00161 aMaterial->GetMaterialPropertiesTable();
00162 if (!aMaterialPropertiesTable) return pParticleChange;
00163
00164 G4MaterialPropertyVector* Rindex =
00165 aMaterialPropertiesTable->GetProperty("RINDEX");
00166 if (!Rindex) return pParticleChange;
00167
00168
00169 const G4double charge = aParticle->GetDefinition()->GetPDGCharge();
00170
00171
00172 const G4double beta = (pPreStepPoint ->GetBeta() +
00173 pPostStepPoint->GetBeta())/2.;
00174
00175 G4double MeanNumberOfPhotons =
00176 GetAverageNumberOfPhotons(charge,beta,aMaterial,Rindex);
00177
00178 if (MeanNumberOfPhotons <= 0.0) {
00179
00180
00181
00182 aParticleChange.SetNumberOfSecondaries(0);
00183
00184 return pParticleChange;
00185
00186 }
00187
00188 G4double step_length;
00189 step_length = aStep.GetStepLength();
00190
00191 MeanNumberOfPhotons = MeanNumberOfPhotons * step_length;
00192
00193 G4int NumPhotons = (G4int) G4Poisson(MeanNumberOfPhotons);
00194
00195 if (NumPhotons <= 0) {
00196
00197
00198
00199 aParticleChange.SetNumberOfSecondaries(0);
00200
00201 return pParticleChange;
00202 }
00203
00205
00206 aParticleChange.SetNumberOfSecondaries(NumPhotons);
00207
00208 if (fTrackSecondariesFirst) {
00209 if (aTrack.GetTrackStatus() == fAlive )
00210 aParticleChange.ProposeTrackStatus(fSuspend);
00211 }
00212
00214
00215 G4double Pmin = Rindex->GetMinLowEdgeEnergy();
00216 G4double Pmax = Rindex->GetMaxLowEdgeEnergy();
00217 G4double dp = Pmax - Pmin;
00218
00219 G4double nMax = Rindex->GetMaxValue();
00220
00221 G4double BetaInverse = 1./beta;
00222
00223 G4double maxCos = BetaInverse / nMax;
00224 G4double maxSin2 = (1.0 - maxCos) * (1.0 + maxCos);
00225
00226 const G4double beta1 = pPreStepPoint ->GetBeta();
00227 const G4double beta2 = pPostStepPoint->GetBeta();
00228
00229 G4double MeanNumberOfPhotons1 =
00230 GetAverageNumberOfPhotons(charge,beta1,aMaterial,Rindex);
00231 G4double MeanNumberOfPhotons2 =
00232 GetAverageNumberOfPhotons(charge,beta2,aMaterial,Rindex);
00233
00234 for (G4int i = 0; i < NumPhotons; i++) {
00235
00236
00237
00238 G4double rand;
00239 G4double sampledEnergy, sampledRI;
00240 G4double cosTheta, sin2Theta;
00241
00242
00243
00244 do {
00245 rand = G4UniformRand();
00246 sampledEnergy = Pmin + rand * dp;
00247 sampledRI = Rindex->Value(sampledEnergy);
00248 cosTheta = BetaInverse / sampledRI;
00249
00250 sin2Theta = (1.0 - cosTheta)*(1.0 + cosTheta);
00251 rand = G4UniformRand();
00252
00253 } while (rand*maxSin2 > sin2Theta);
00254
00255
00256
00257
00258 rand = G4UniformRand();
00259
00260 G4double phi = twopi*rand;
00261 G4double sinPhi = std::sin(phi);
00262 G4double cosPhi = std::cos(phi);
00263
00264
00265
00266
00267
00268 G4double sinTheta = std::sqrt(sin2Theta);
00269 G4double px = sinTheta*cosPhi;
00270 G4double py = sinTheta*sinPhi;
00271 G4double pz = cosTheta;
00272
00273
00274
00275
00276
00277
00278 G4ParticleMomentum photonMomentum(px, py, pz);
00279
00280
00281
00282
00283 photonMomentum.rotateUz(p0);
00284
00285
00286
00287 G4double sx = cosTheta*cosPhi;
00288 G4double sy = cosTheta*sinPhi;
00289 G4double sz = -sinTheta;
00290
00291 G4ThreeVector photonPolarization(sx, sy, sz);
00292
00293
00294
00295 photonPolarization.rotateUz(p0);
00296
00297
00298
00299 G4DynamicParticle* aCerenkovPhoton =
00300 new G4DynamicParticle(G4OpticalPhoton::OpticalPhoton(),
00301 photonMomentum);
00302 aCerenkovPhoton->SetPolarization
00303 (photonPolarization.x(),
00304 photonPolarization.y(),
00305 photonPolarization.z());
00306
00307 aCerenkovPhoton->SetKineticEnergy(sampledEnergy);
00308
00309
00310
00311 G4double delta, NumberOfPhotons, N;
00312
00313 do {
00314 rand = G4UniformRand();
00315 delta = rand * aStep.GetStepLength();
00316 NumberOfPhotons = MeanNumberOfPhotons1 - delta *
00317 (MeanNumberOfPhotons1-MeanNumberOfPhotons2)/
00318 aStep.GetStepLength();
00319 N = G4UniformRand() *
00320 std::max(MeanNumberOfPhotons1,MeanNumberOfPhotons2);
00321 } while (N > NumberOfPhotons);
00322
00323 G4double deltaTime = delta /
00324 ((pPreStepPoint->GetVelocity()+
00325 pPostStepPoint->GetVelocity())/2.);
00326
00327 G4double aSecondaryTime = t0 + deltaTime;
00328
00329 G4ThreeVector aSecondaryPosition =
00330 x0 + rand * aStep.GetDeltaPosition();
00331
00332 G4Track* aSecondaryTrack =
00333 new G4Track(aCerenkovPhoton,aSecondaryTime,aSecondaryPosition);
00334
00335 aSecondaryTrack->SetTouchableHandle(
00336 aStep.GetPreStepPoint()->GetTouchableHandle());
00337
00338 aSecondaryTrack->SetParentID(aTrack.GetTrackID());
00339
00340 aParticleChange.AddSecondary(aSecondaryTrack);
00341 }
00342
00343 if (verboseLevel>0) {
00344 G4cout <<"\n Exiting from G4Cerenkov::DoIt -- NumberOfSecondaries = "
00345 << aParticleChange.GetNumberOfSecondaries() << G4endl;
00346 }
00347
00348 return pParticleChange;
00349 }
00350
00351
00352
00353
00354
00355 void G4Cerenkov::BuildThePhysicsTable()
00356 {
00357 if (thePhysicsTable) return;
00358
00359 const G4MaterialTable* theMaterialTable=
00360 G4Material::GetMaterialTable();
00361 G4int numOfMaterials = G4Material::GetNumberOfMaterials();
00362
00363
00364
00365 thePhysicsTable = new G4PhysicsTable(numOfMaterials);
00366
00367
00368
00369 for (G4int i=0 ; i < numOfMaterials; i++)
00370 {
00371 G4PhysicsOrderedFreeVector* aPhysicsOrderedFreeVector =
00372 new G4PhysicsOrderedFreeVector();
00373
00374
00375
00376
00377 G4Material* aMaterial = (*theMaterialTable)[i];
00378
00379 G4MaterialPropertiesTable* aMaterialPropertiesTable =
00380 aMaterial->GetMaterialPropertiesTable();
00381
00382 if (aMaterialPropertiesTable) {
00383
00384 G4MaterialPropertyVector* theRefractionIndexVector =
00385 aMaterialPropertiesTable->GetProperty("RINDEX");
00386
00387 if (theRefractionIndexVector) {
00388
00389
00390
00391
00392 G4double currentRI = (*theRefractionIndexVector)[0];
00393
00394 if (currentRI > 1.0) {
00395
00396
00397
00398
00399 G4double currentPM = theRefractionIndexVector->
00400 Energy(0);
00401 G4double currentCAI = 0.0;
00402
00403 aPhysicsOrderedFreeVector->
00404 InsertValues(currentPM , currentCAI);
00405
00406
00407
00408 G4double prevPM = currentPM;
00409 G4double prevCAI = currentCAI;
00410 G4double prevRI = currentRI;
00411
00412
00413
00414
00415 for (size_t ii = 1;
00416 ii < theRefractionIndexVector->GetVectorLength();
00417 ++ii)
00418 {
00419 currentRI = (*theRefractionIndexVector)[ii];
00420 currentPM = theRefractionIndexVector->Energy(ii);
00421
00422 currentCAI = 0.5*(1.0/(prevRI*prevRI) +
00423 1.0/(currentRI*currentRI));
00424
00425 currentCAI = prevCAI +
00426 (currentPM - prevPM) * currentCAI;
00427
00428 aPhysicsOrderedFreeVector->
00429 InsertValues(currentPM, currentCAI);
00430
00431 prevPM = currentPM;
00432 prevCAI = currentCAI;
00433 prevRI = currentRI;
00434 }
00435
00436 }
00437 }
00438 }
00439
00440
00441
00442
00443
00444
00445 thePhysicsTable->insertAt(i,aPhysicsOrderedFreeVector);
00446
00447 }
00448 }
00449
00450
00451
00452
00453
00454 G4double G4Cerenkov::GetMeanFreePath(const G4Track&,
00455 G4double,
00456 G4ForceCondition*)
00457 {
00458 return 1.;
00459 }
00460
00461 G4double G4Cerenkov::PostStepGetPhysicalInteractionLength(
00462 const G4Track& aTrack,
00463 G4double,
00464 G4ForceCondition* condition)
00465 {
00466 *condition = NotForced;
00467 G4double StepLimit = DBL_MAX;
00468
00469 const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
00470 const G4Material* aMaterial = aTrack.GetMaterial();
00471 const G4MaterialCutsCouple* couple = aTrack.GetMaterialCutsCouple();
00472
00473 G4double kineticEnergy = aParticle->GetKineticEnergy();
00474 const G4ParticleDefinition* particleType = aParticle->GetDefinition();
00475 G4double mass = particleType->GetPDGMass();
00476
00477
00478 G4double beta = aParticle->GetTotalMomentum() /
00479 aParticle->GetTotalEnergy();
00480
00481 G4double gamma = aParticle->GetTotalEnergy()/mass;
00482
00483 G4MaterialPropertiesTable* aMaterialPropertiesTable =
00484 aMaterial->GetMaterialPropertiesTable();
00485
00486 G4MaterialPropertyVector* Rindex = NULL;
00487
00488 if (aMaterialPropertiesTable)
00489 Rindex = aMaterialPropertiesTable->GetProperty("RINDEX");
00490
00491 G4double nMax;
00492 if (Rindex) {
00493 nMax = Rindex->GetMaxValue();
00494 } else {
00495 return StepLimit;
00496 }
00497
00498 G4double BetaMin = 1./nMax;
00499 if ( BetaMin >= 1. ) return StepLimit;
00500
00501 G4double GammaMin = 1./std::sqrt(1.-BetaMin*BetaMin);
00502
00503 if (gamma < GammaMin ) return StepLimit;
00504
00505 G4double kinEmin = mass*(GammaMin-1.);
00506
00507 G4double RangeMin = G4LossTableManager::Instance()->
00508 GetRange(particleType,
00509 kinEmin,
00510 couple);
00511 G4double Range = G4LossTableManager::Instance()->
00512 GetRange(particleType,
00513 kineticEnergy,
00514 couple);
00515
00516 G4double Step = Range - RangeMin;
00517 if (Step < 1.*um ) return StepLimit;
00518
00519 if (Step > 0. && Step < StepLimit) StepLimit = Step;
00520
00521
00522
00523
00524
00525 if (fMaxPhotons > 0) {
00526
00527
00528 const G4double charge = aParticle->
00529 GetDefinition()->GetPDGCharge();
00530
00531 G4double MeanNumberOfPhotons =
00532 GetAverageNumberOfPhotons(charge,beta,aMaterial,Rindex);
00533
00534 Step = 0.;
00535 if (MeanNumberOfPhotons > 0.0) Step = fMaxPhotons /
00536 MeanNumberOfPhotons;
00537
00538 if (Step > 0. && Step < StepLimit) StepLimit = Step;
00539 }
00540
00541
00542 if (fMaxBetaChange > 0.) {
00543
00544 G4double dedx = G4LossTableManager::Instance()->
00545 GetDEDX(particleType,
00546 kineticEnergy,
00547 couple);
00548
00549 G4double deltaGamma = gamma -
00550 1./std::sqrt(1.-beta*beta*
00551 (1.-fMaxBetaChange)*
00552 (1.-fMaxBetaChange));
00553
00554 Step = mass * deltaGamma / dedx;
00555
00556 if (Step > 0. && Step < StepLimit) StepLimit = Step;
00557
00558 }
00559
00560 *condition = StronglyForced;
00561 return StepLimit;
00562 }
00563
00564
00565
00566
00567
00568
00569
00570 G4double
00571 G4Cerenkov::GetAverageNumberOfPhotons(const G4double charge,
00572 const G4double beta,
00573 const G4Material* aMaterial,
00574 G4MaterialPropertyVector* Rindex) const
00575 {
00576 const G4double Rfact = 369.81/(eV * cm);
00577
00578 if(beta <= 0.0)return 0.0;
00579
00580 G4double BetaInverse = 1./beta;
00581
00582
00583
00584
00585
00586 G4int materialIndex = aMaterial->GetIndex();
00587
00588
00589
00590 G4PhysicsOrderedFreeVector* CerenkovAngleIntegrals =
00591 (G4PhysicsOrderedFreeVector*)((*thePhysicsTable)(materialIndex));
00592
00593 if(!(CerenkovAngleIntegrals->IsFilledVectorExist()))return 0.0;
00594
00595
00596 G4double Pmin = Rindex->GetMinLowEdgeEnergy();
00597 G4double Pmax = Rindex->GetMaxLowEdgeEnergy();
00598
00599
00600 G4double nMin = Rindex->GetMinValue();
00601 G4double nMax = Rindex->GetMaxValue();
00602
00603
00604 G4double CAImax = CerenkovAngleIntegrals->GetMaxValue();
00605
00606 G4double dp, ge;
00607
00608
00609
00610 if (nMax < BetaInverse) {
00611 dp = 0;
00612 ge = 0;
00613 }
00614
00615
00616
00617 else if (nMin > BetaInverse) {
00618 dp = Pmax - Pmin;
00619 ge = CAImax;
00620 }
00621
00622
00623
00624
00625
00626
00627
00628 else {
00629 Pmin = Rindex->GetEnergy(BetaInverse);
00630 dp = Pmax - Pmin;
00631
00632
00633
00634 G4bool isOutRange;
00635 G4double CAImin = CerenkovAngleIntegrals->
00636 GetValue(Pmin, isOutRange);
00637 ge = CAImax - CAImin;
00638
00639 if (verboseLevel>0) {
00640 G4cout << "CAImin = " << CAImin << G4endl;
00641 G4cout << "ge = " << ge << G4endl;
00642 }
00643 }
00644
00645
00646 G4double NumPhotons = Rfact * charge/eplus * charge/eplus *
00647 (dp - ge * BetaInverse*BetaInverse);
00648
00649 return NumPhotons;
00650 }