#include <G4PAIModel.hh>
Inheritance diagram for G4PAIModel:
Definition at line 67 of file G4PAIModel.hh.
G4PAIModel::G4PAIModel | ( | const G4ParticleDefinition * | p = 0 , |
|
const G4String & | nam = "PAI" | |||
) |
Definition at line 74 of file G4PAIModel.cc.
References G4Electron::Electron(), and G4Positron::Positron().
00075 : G4VEmModel(nam),G4VEmFluctuationModel(nam), 00076 fVerbose(0), 00077 fLowestGamma(1.005), 00078 fHighestGamma(10000.), 00079 fTotBin(200), 00080 fMeanNumber(20), 00081 fParticle(0), 00082 fHighKinEnergy(100.*TeV), 00083 fTwoln10(2.0*log(10.0)), 00084 fBg2lim(0.0169), 00085 fTaulim(8.4146e-3) 00086 { 00087 fElectron = G4Electron::Electron(); 00088 fPositron = G4Positron::Positron(); 00089 00090 fPAItransferTable = 0; 00091 fPAIdEdxTable = 0; 00092 fSandiaPhotoAbsCof = 0; 00093 fdEdxVector = 0; 00094 fLambdaVector = 0; 00095 fdNdxCutVector = 0; 00096 fParticleEnergyVector = 0; 00097 fSandiaIntervalNumber = 0; 00098 fMatIndex = 0; 00099 fDeltaCutInKinEnergy = 0.0; 00100 fParticleChange = 0; 00101 fMaterial = 0; 00102 fCutCouple = 0; 00103 00104 if(p) { SetParticle(p); } 00105 else { SetParticle(fElectron); } 00106 00107 isInitialised = false; 00108 }
G4PAIModel::~G4PAIModel | ( | ) | [virtual] |
Definition at line 112 of file G4PAIModel.cc.
References G4PhysicsTable::clearAndDestroy().
00113 { 00114 // G4cout << "PAI: start destruction" << G4endl; 00115 if(fParticleEnergyVector) delete fParticleEnergyVector; 00116 if(fdEdxVector) delete fdEdxVector ; 00117 if(fLambdaVector) delete fLambdaVector; 00118 if(fdNdxCutVector) delete fdNdxCutVector; 00119 00120 if( fPAItransferTable ) 00121 { 00122 fPAItransferTable->clearAndDestroy(); 00123 delete fPAItransferTable ; 00124 } 00125 if( fPAIdEdxTable ) 00126 { 00127 fPAIdEdxTable->clearAndDestroy(); 00128 delete fPAIdEdxTable ; 00129 } 00130 if(fSandiaPhotoAbsCof) 00131 { 00132 for(G4int i=0;i<fSandiaIntervalNumber;i++) 00133 { 00134 delete[] fSandiaPhotoAbsCof[i]; 00135 } 00136 delete[] fSandiaPhotoAbsCof; 00137 } 00138 //G4cout << "PAI: end destruction" << G4endl; 00139 }
void G4PAIModel::BuildLambdaVector | ( | ) |
Definition at line 346 of file G4PAIModel.cc.
References DBL_MAX, G4cout, G4endl, GetdNdxCut(), G4InuclParticleNames::lambda, and G4PhysicsVector::PutValue().
Referenced by Initialise().
00347 { 00348 fLambdaVector = new G4PhysicsLogVector( fLowestKineticEnergy, 00349 fHighestKineticEnergy, 00350 fTotBin ) ; 00351 fdNdxCutVector = new G4PhysicsLogVector( fLowestKineticEnergy, 00352 fHighestKineticEnergy, 00353 fTotBin ) ; 00354 if(fVerbose > 1) 00355 { 00356 G4cout<<"PAIModel DeltaCutInKineticEnergyNow = " 00357 <<fDeltaCutInKinEnergy/keV<<" keV"<<G4endl; 00358 } 00359 for (G4int i = 0 ; i <= fTotBin ; i++ ) 00360 { 00361 G4double dNdxCut = GetdNdxCut(i,fDeltaCutInKinEnergy) ; 00362 //G4cout << "i= " << i << " x= " << dNdxCut << G4endl; 00363 if(dNdxCut < 0.0) dNdxCut = 0.0; 00364 // G4double lambda = dNdxCut <= DBL_MIN ? DBL_MAX: 1.0/dNdxCut ; 00365 G4double lambda = DBL_MAX; 00366 if(dNdxCut > 0.0) lambda = 1.0/dNdxCut; 00367 00368 fLambdaVector->PutValue(i, lambda) ; 00369 fdNdxCutVector->PutValue(i, dNdxCut) ; 00370 } 00371 }
void G4PAIModel::BuildPAIonisationTable | ( | ) |
Definition at line 282 of file G4PAIModel.cc.
References DBL_MIN, G4PAIySection::GetIntegralPAIdEdx(), G4PAIySection::GetIntegralPAIySection(), G4PhysicsVector::GetLowEdgeEnergy(), G4PAIySection::GetMeanEnergyLoss(), G4SandiaTable::GetSandiaCofForMaterialPAI(), G4Material::GetSandiaTable(), G4PAIySection::GetSplineEnergy(), G4PAIySection::GetSplineSize(), G4PAIySection::Initialize(), G4PhysicsTable::insertAt(), MaxSecondaryEnergy(), CLHEP::detail::n, G4PhysicsVector::PutValue(), and G4PhysicsFreeVector::PutValue().
Referenced by Initialise().
00283 { 00284 G4double LowEdgeEnergy , ionloss ; 00285 G4double tau, Tmax, Tmin, Tkin, deltaLow, /*gamma,*/ bg2 ; 00286 00287 fdEdxVector = new G4PhysicsLogVector( fLowestKineticEnergy, 00288 fHighestKineticEnergy, 00289 fTotBin); 00290 G4SandiaTable* sandia = fMaterial->GetSandiaTable(); 00291 00292 Tmin = sandia->GetSandiaCofForMaterialPAI(0,0)*keV; 00293 deltaLow = 100.*eV; // 0.5*eV ; 00294 00295 for (G4int i = 0 ; i <= fTotBin ; i++) //The loop for the kinetic energy 00296 { 00297 LowEdgeEnergy = fParticleEnergyVector->GetLowEdgeEnergy(i) ; 00298 tau = LowEdgeEnergy/fMass ; 00299 //gamma = tau +1. ; 00300 // G4cout<<"gamma = "<<gamma<<endl ; 00301 bg2 = tau*( tau + 2. ); 00302 Tmax = MaxSecondaryEnergy(fParticle, LowEdgeEnergy); 00303 // Tmax = std::min(fDeltaCutInKinEnergy, Tmax); 00304 Tkin = Tmax ; 00305 00306 if ( Tmax < Tmin + deltaLow ) // low energy safety 00307 Tkin = Tmin + deltaLow ; 00308 00309 fPAIySection.Initialize(fMaterial, Tkin, bg2); 00310 00311 // G4cout<<"ionloss = "<<ionloss*cm/keV<<" keV/cm"<<endl ; 00312 // G4cout<<"n1 = "<<protonPAI.GetIntegralPAIxSection(1)*cm<<" 1/cm"<<endl ; 00313 // G4cout<<"protonPAI.GetSplineSize() = "<< 00314 // protonPAI.GetSplineSize()<<G4endl<<G4endl ; 00315 00316 G4int n = fPAIySection.GetSplineSize(); 00317 G4PhysicsFreeVector* transferVector = new G4PhysicsFreeVector(n) ; 00318 G4PhysicsFreeVector* dEdxVector = new G4PhysicsFreeVector(n); 00319 00320 for( G4int k = 0 ; k < n; k++ ) 00321 { 00322 transferVector->PutValue( k , 00323 fPAIySection.GetSplineEnergy(k+1), 00324 fPAIySection.GetIntegralPAIySection(k+1) ) ; 00325 dEdxVector->PutValue( k , 00326 fPAIySection.GetSplineEnergy(k+1), 00327 fPAIySection.GetIntegralPAIdEdx(k+1) ) ; 00328 } 00329 ionloss = fPAIySection.GetMeanEnergyLoss() ; // total <dE/dx> 00330 00331 if ( ionloss < DBL_MIN) { ionloss = 0.0; } 00332 fdEdxVector->PutValue(i,ionloss) ; 00333 00334 fPAItransferTable->insertAt(i,transferVector) ; 00335 fPAIdEdxTable->insertAt(i,dEdxVector) ; 00336 00337 } // end of Tkin loop 00338 }
G4double G4PAIModel::ComputeDEDXPerVolume | ( | const G4Material * | , | |
const G4ParticleDefinition * | , | |||
G4double | kineticEnergy, | |||
G4double | cutEnergy | |||
) | [virtual] |
Reimplemented from G4VEmModel.
Definition at line 468 of file G4PAIModel.cc.
References G4VEmModel::CurrentCouple(), GetdEdxCut(), G4ParticleDefinition::GetPDGCharge(), and G4ParticleDefinition::GetPDGMass().
00472 { 00473 G4int iTkin,iPlace; 00474 00475 //G4double cut = std::min(MaxSecondaryEnergy(p, kineticEnergy), cutEnergy); 00476 G4double cut = cutEnergy; 00477 00478 G4double massRatio = fMass/p->GetPDGMass(); 00479 G4double scaledTkin = kineticEnergy*massRatio; 00480 G4double charge = p->GetPDGCharge(); 00481 G4double charge2 = charge*charge; 00482 const G4MaterialCutsCouple* matCC = CurrentCouple(); 00483 00484 size_t jMat = 0; 00485 for(;jMat < fMaterialCutsCoupleVector.size(); ++jMat ) 00486 { 00487 if( matCC == fMaterialCutsCoupleVector[jMat] ) break; 00488 } 00489 //G4cout << "jMat= " << jMat 00490 // << " jMax= " << fMaterialCutsCoupleVector.size() 00491 // << " matCC: " << matCC; 00492 //if(matCC) G4cout << " mat: " << matCC->GetMaterial()->GetName(); 00493 // G4cout << G4endl; 00494 if(jMat == fMaterialCutsCoupleVector.size()) return 0.0; 00495 00496 fPAIdEdxTable = fPAIdEdxBank[jMat]; 00497 fdEdxVector = fdEdxTable[jMat]; 00498 for(iTkin = 0 ; iTkin <= fTotBin ; iTkin++) 00499 { 00500 if(scaledTkin < fParticleEnergyVector->GetLowEdgeEnergy(iTkin)) break ; 00501 } 00502 //G4cout << "E= " << scaledTkin << " iTkin= " << iTkin 00503 // << " " << fdEdxVector->GetVectorLength()<<G4endl; 00504 iPlace = iTkin - 1; 00505 if(iPlace < 0) iPlace = 0; 00506 else if(iPlace >= fTotBin) iPlace = fTotBin-1; 00507 G4double dEdx = charge2*( (*fdEdxVector)(iPlace) - GetdEdxCut(iPlace,cut) ); 00508 if( dEdx < 0.) dEdx = 0.; 00509 return dEdx; 00510 }
void G4PAIModel::ComputeSandiaPhotoAbsCof | ( | ) |
Definition at line 232 of file G4PAIModel.cc.
References G4Material::GetMaterialTable(), G4SandiaTable::GetPhotoAbsorpCof(), G4SandiaTable::SandiaIntervals(), and G4SandiaTable::SandiaMixing().
00233 { 00234 G4int i, j; 00235 const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable(); 00236 00237 G4SandiaTable thisMaterialSandiaTable(fMatIndex) ; 00238 G4int numberOfElements = 00239 (*theMaterialTable)[fMatIndex]->GetNumberOfElements(); 00240 00241 G4int* thisMaterialZ = new G4int[numberOfElements] ; 00242 00243 for(i=0;i<numberOfElements;i++) 00244 { 00245 thisMaterialZ[i] = 00246 (G4int)(*theMaterialTable)[fMatIndex]->GetElement(i)->GetZ() ; 00247 } 00248 fSandiaIntervalNumber = thisMaterialSandiaTable.SandiaIntervals 00249 (thisMaterialZ,numberOfElements) ; 00250 00251 fSandiaIntervalNumber = thisMaterialSandiaTable.SandiaMixing 00252 ( thisMaterialZ , 00253 (*theMaterialTable)[fMatIndex]->GetFractionVector() , 00254 numberOfElements,fSandiaIntervalNumber) ; 00255 00256 fSandiaPhotoAbsCof = new G4double*[fSandiaIntervalNumber] ; 00257 00258 for(i=0; i<fSandiaIntervalNumber; i++) 00259 { 00260 fSandiaPhotoAbsCof[i] = new G4double[5]; 00261 } 00262 00263 for( i = 0 ; i < fSandiaIntervalNumber ; i++ ) 00264 { 00265 fSandiaPhotoAbsCof[i][0] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i+1,0); 00266 00267 for( j = 1; j < 5 ; j++ ) 00268 { 00269 fSandiaPhotoAbsCof[i][j] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i+1,j)* 00270 (*theMaterialTable)[fMatIndex]->GetDensity() ; 00271 } 00272 } 00273 delete[] thisMaterialZ; 00274 }
G4double G4PAIModel::CrossSectionPerVolume | ( | const G4Material * | , | |
const G4ParticleDefinition * | , | |||
G4double | kineticEnergy, | |||
G4double | cutEnergy, | |||
G4double | maxEnergy | |||
) | [virtual] |
Reimplemented from G4VEmModel.
Definition at line 514 of file G4PAIModel.cc.
References G4VEmModel::CurrentCouple(), DBL_MIN, GetdNdxCut(), G4ParticleDefinition::GetPDGCharge(), G4ParticleDefinition::GetPDGMass(), and MaxSecondaryEnergy().
00519 { 00520 G4int iTkin,iPlace; 00521 G4double tmax = std::min(MaxSecondaryEnergy(p, kineticEnergy), maxEnergy); 00522 if(tmax <= cutEnergy) return 0.0; 00523 G4double massRatio = fMass/p->GetPDGMass(); 00524 G4double scaledTkin = kineticEnergy*massRatio; 00525 G4double charge = p->GetPDGCharge(); 00526 G4double charge2 = charge*charge, cross, cross1, cross2; 00527 const G4MaterialCutsCouple* matCC = CurrentCouple(); 00528 00529 size_t jMat = 0; 00530 for(;jMat < fMaterialCutsCoupleVector.size(); ++jMat ) 00531 { 00532 if( matCC == fMaterialCutsCoupleVector[jMat] ) break; 00533 } 00534 if(jMat == fMaterialCutsCoupleVector.size()) return 0.0; 00535 00536 fPAItransferTable = fPAIxscBank[jMat]; 00537 00538 for(iTkin = 0 ; iTkin <= fTotBin ; iTkin++) 00539 { 00540 if(scaledTkin < fParticleEnergyVector->GetLowEdgeEnergy(iTkin)) break ; 00541 } 00542 iPlace = iTkin - 1; 00543 if(iPlace < 0) iPlace = 0; 00544 else if(iPlace >= fTotBin) iPlace = fTotBin-1; 00545 00546 //G4cout<<"iPlace = "<<iPlace<<"; tmax = " 00547 // <<tmax<<"; cutEnergy = "<<cutEnergy<<G4endl; 00548 cross1 = GetdNdxCut(iPlace,tmax) ; 00549 //G4cout<<"cross1 = "<<cross1<<G4endl; 00550 cross2 = GetdNdxCut(iPlace,cutEnergy) ; 00551 //G4cout<<"cross2 = "<<cross2<<G4endl; 00552 cross = (cross2-cross1)*charge2; 00553 //G4cout<<"cross = "<<cross<<G4endl; 00554 if( cross < DBL_MIN) cross = 0.0; 00555 // if( cross2 < DBL_MIN) cross2 = DBL_MIN; 00556 00557 // return cross2; 00558 return cross; 00559 }
void G4PAIModel::DefineForRegion | ( | const G4Region * | r | ) | [virtual] |
G4double G4PAIModel::Dispersion | ( | const G4Material * | , | |
const G4DynamicParticle * | , | |||
G4double & | , | |||
G4double & | ||||
) | [virtual] |
Implements G4VEmFluctuationModel.
Definition at line 958 of file G4PAIModel.cc.
References SampleFluctuations().
00962 { 00963 G4double loss, sumLoss=0., sumLoss2=0., sigma2, meanLoss=0.; 00964 for(G4int i = 0 ; i < fMeanNumber; i++) 00965 { 00966 loss = SampleFluctuations(material,aParticle,tmax,step,meanLoss); 00967 sumLoss += loss; 00968 sumLoss2 += loss*loss; 00969 } 00970 meanLoss = sumLoss/fMeanNumber; 00971 sigma2 = meanLoss*meanLoss + (sumLoss2-2*sumLoss*meanLoss)/fMeanNumber; 00972 return sigma2; 00973 }
Definition at line 424 of file G4PAIModel.cc.
Referenced by ComputeDEDXPerVolume().
00425 { 00426 G4int iTransfer; 00427 G4double x1, x2, y1, y2, dEdxCut; 00428 //G4cout<<"iPlace = "<<iPlace<<"; "<<"transferCut = "<<transferCut<<G4endl; 00429 // G4cout<<"size = "<<G4int((*fPAIdEdxTable)(iPlace)->GetVectorLength()) 00430 // <<G4endl; 00431 for( iTransfer = 0 ; 00432 iTransfer < G4int((*fPAIdEdxTable)(iPlace)->GetVectorLength()) ; 00433 iTransfer++) 00434 { 00435 if(transferCut <= (*fPAIdEdxTable)(iPlace)->GetLowEdgeEnergy(iTransfer)) 00436 { 00437 break ; 00438 } 00439 } 00440 if ( iTransfer >= G4int((*fPAIdEdxTable)(iPlace)->GetVectorLength()) ) 00441 { 00442 iTransfer = (*fPAIdEdxTable)(iPlace)->GetVectorLength() - 1 ; 00443 } 00444 y1 = (*(*fPAIdEdxTable)(iPlace))(iTransfer-1) ; 00445 y2 = (*(*fPAIdEdxTable)(iPlace))(iTransfer) ; 00446 if(y1 < 0.0) y1 = 0.0; 00447 if(y2 < 0.0) y2 = 0.0; 00448 //G4cout<<"y1 = "<<y1<<"; "<<"y2 = "<<y2<<G4endl; 00449 x1 = (*fPAIdEdxTable)(iPlace)->GetLowEdgeEnergy(iTransfer-1) ; 00450 x2 = (*fPAIdEdxTable)(iPlace)->GetLowEdgeEnergy(iTransfer) ; 00451 //G4cout<<"x1 = "<<x1<<"; "<<"x2 = "<<x2<<G4endl; 00452 00453 if ( y1 == y2 ) dEdxCut = y2 ; 00454 else 00455 { 00456 // if ( x1 == x2 ) dEdxCut = y1 + (y2 - y1)*G4UniformRand() ; 00457 // if ( std::abs(x1-x2) <= eV ) dEdxCut = y1 + (y2 - y1)*G4UniformRand() ; 00458 if ( std::abs(x1-x2) <= eV ) dEdxCut = y1 + (y2 - y1)*0.5 ; 00459 else dEdxCut = y1 + (transferCut - x1)*(y2 - y1)/(x2 - x1) ; 00460 } 00461 //G4cout<<""<<dEdxCut<<G4endl; 00462 if(dEdxCut < 0.0) dEdxCut = 0.0; 00463 return dEdxCut ; 00464 }
Definition at line 378 of file G4PAIModel.cc.
Referenced by BuildLambdaVector(), and CrossSectionPerVolume().
00379 { 00380 G4int iTransfer; 00381 G4double x1, x2, y1, y2, dNdxCut; 00382 // G4cout<<"iPlace = "<<iPlace<<"; "<<"transferCut = "<<transferCut<<G4endl; 00383 // G4cout<<"size = "<<G4int((*fPAItransferTable)(iPlace)->GetVectorLength()) 00384 // <<G4endl; 00385 for( iTransfer = 0 ; 00386 iTransfer < G4int((*fPAItransferTable)(iPlace)->GetVectorLength()) ; 00387 iTransfer++) 00388 { 00389 if(transferCut <= (*fPAItransferTable)(iPlace)->GetLowEdgeEnergy(iTransfer)) 00390 { 00391 break ; 00392 } 00393 } 00394 if ( iTransfer >= G4int((*fPAItransferTable)(iPlace)->GetVectorLength()) ) 00395 { 00396 iTransfer = (*fPAItransferTable)(iPlace)->GetVectorLength() - 1 ; 00397 } 00398 y1 = (*(*fPAItransferTable)(iPlace))(iTransfer-1) ; 00399 y2 = (*(*fPAItransferTable)(iPlace))(iTransfer) ; 00400 if(y1 < 0.0) y1 = 0.0; 00401 if(y2 < 0.0) y2 = 0.0; 00402 // G4cout<<"y1 = "<<y1<<"; "<<"y2 = "<<y2<<G4endl; 00403 x1 = (*fPAItransferTable)(iPlace)->GetLowEdgeEnergy(iTransfer-1) ; 00404 x2 = (*fPAItransferTable)(iPlace)->GetLowEdgeEnergy(iTransfer) ; 00405 // G4cout<<"x1 = "<<x1<<"; "<<"x2 = "<<x2<<G4endl; 00406 00407 if ( y1 == y2 ) dNdxCut = y2 ; 00408 else 00409 { 00410 // if ( x1 == x2 ) dNdxCut = y1 + (y2 - y1)*G4UniformRand() ; 00411 // if ( std::abs(x1-x2) <= eV ) dNdxCut = y1 + (y2 - y1)*G4UniformRand() ; 00412 if ( std::abs(x1-x2) <= eV ) dNdxCut = y1 + (y2 - y1)*0.5 ; 00413 else dNdxCut = y1 + (transferCut - x1)*(y2 - y1)/(x2 - x1) ; 00414 } 00415 // G4cout<<""<<dNdxCut<<G4endl; 00416 if(dNdxCut < 0.0) dNdxCut = 0.0; 00417 return dNdxCut ; 00418 }
Definition at line 737 of file G4PAIModel.cc.
References G4UniformRand.
Referenced by GetPostStepTransfer(), and SampleFluctuations().
00738 { 00739 G4int iTransferMax; 00740 G4double x1, x2, y1, y2, energyTransfer; 00741 00742 if(iTransfer == 0) 00743 { 00744 energyTransfer = (*fPAItransferTable)(iPlace)->GetLowEdgeEnergy(iTransfer); 00745 } 00746 else 00747 { 00748 iTransferMax = G4int((*fPAItransferTable)(iPlace)->GetVectorLength()); 00749 00750 if ( iTransfer >= iTransferMax ) iTransfer = iTransferMax - 1; 00751 00752 y1 = (*(*fPAItransferTable)(iPlace))(iTransfer-1); 00753 y2 = (*(*fPAItransferTable)(iPlace))(iTransfer); 00754 00755 x1 = (*fPAItransferTable)(iPlace)->GetLowEdgeEnergy(iTransfer-1); 00756 x2 = (*fPAItransferTable)(iPlace)->GetLowEdgeEnergy(iTransfer); 00757 00758 if ( x1 == x2 ) energyTransfer = x2; 00759 else 00760 { 00761 if ( y1 == y2 ) energyTransfer = x1 + (x2 - x1)*G4UniformRand(); 00762 else 00763 { 00764 energyTransfer = x1 + (position - y1)*(x2 - x1)/(y2 - y1); 00765 } 00766 } 00767 } 00768 return energyTransfer; 00769 }
Definition at line 649 of file G4PAIModel.cc.
References G4UniformRand, GetEnergyTransfer(), G4PhysicsVector::GetLowEdgeEnergy(), and position.
Referenced by SampleSecondaries().
00650 { 00651 // G4cout<<"G4PAIModel::GetPostStepTransfer"<<G4endl ; 00652 00653 G4int iTkin, iTransfer, iPlace ; 00654 G4double transfer = 0.0, position, dNdxCut1, dNdxCut2, E1, E2, W1, W2, W ; 00655 00656 for(iTkin=0;iTkin<=fTotBin;iTkin++) 00657 { 00658 if(scaledTkin < fParticleEnergyVector->GetLowEdgeEnergy(iTkin)) break ; 00659 } 00660 iPlace = iTkin - 1 ; 00661 // G4cout<<"from search, iPlace = "<<iPlace<<G4endl ; 00662 if(iPlace < 0) iPlace = 0; 00663 else if(iPlace >= fTotBin) iPlace = fTotBin-1; 00664 dNdxCut1 = (*fdNdxCutVector)(iPlace) ; 00665 // G4cout<<"dNdxCut1 = "<<dNdxCut1<<G4endl ; 00666 00667 00668 if(iTkin == fTotBin) // Fermi plato, try from left 00669 { 00670 position = dNdxCut1*G4UniformRand() ; 00671 00672 for( iTransfer = 0; 00673 iTransfer < G4int((*fPAItransferTable)(iPlace)->GetVectorLength()); iTransfer++ ) 00674 { 00675 if(position >= (*(*fPAItransferTable)(iPlace))(iTransfer)) break ; 00676 } 00677 transfer = GetEnergyTransfer(iPlace,position,iTransfer); 00678 } 00679 else 00680 { 00681 dNdxCut2 = (*fdNdxCutVector)(iPlace+1) ; 00682 // G4cout<<"dNdxCut2 = "<<dNdxCut2<<G4endl ; 00683 if(iTkin == 0) // Tkin is too small, trying from right only 00684 { 00685 position = dNdxCut2*G4UniformRand() ; 00686 00687 for( iTransfer = 0; 00688 iTransfer < G4int((*fPAItransferTable)(iPlace+1)->GetVectorLength()); iTransfer++ ) 00689 { 00690 if(position >= (*(*fPAItransferTable)(iPlace+1))(iTransfer)) break ; 00691 } 00692 transfer = GetEnergyTransfer(iPlace+1,position,iTransfer); 00693 } 00694 else // general case: Tkin between two vectors of the material 00695 { 00696 E1 = fParticleEnergyVector->GetLowEdgeEnergy(iTkin - 1) ; 00697 E2 = fParticleEnergyVector->GetLowEdgeEnergy(iTkin) ; 00698 W = 1.0/(E2 - E1) ; 00699 W1 = (E2 - scaledTkin)*W ; 00700 W2 = (scaledTkin - E1)*W ; 00701 00702 position = ( dNdxCut1*W1 + dNdxCut2*W2 )*G4UniformRand() ; 00703 00704 // G4cout<<position<<"\t" ; 00705 00706 G4int iTrMax1, iTrMax2, iTrMax; 00707 00708 iTrMax1 = G4int((*fPAItransferTable)(iPlace)->GetVectorLength()); 00709 iTrMax2 = G4int((*fPAItransferTable)(iPlace+1)->GetVectorLength()); 00710 00711 if (iTrMax1 >= iTrMax2) iTrMax = iTrMax2; 00712 else iTrMax = iTrMax1; 00713 00714 00715 for( iTransfer = 0; iTransfer < iTrMax; iTransfer++ ) 00716 { 00717 if( position >= 00718 ( (*(*fPAItransferTable)(iPlace))(iTransfer)*W1 + 00719 (*(*fPAItransferTable)(iPlace+1))(iTransfer)*W2) ) break ; 00720 } 00721 transfer = GetEnergyTransfer(iPlace,position,iTransfer); 00722 } 00723 } 00724 // G4cout<<"PAImodel PostStepTransfer = "<<transfer/keV<<" keV"<<G4endl ; 00725 if(transfer < 0.0 ) transfer = 0.0 ; 00726 // if(transfer < DBL_MIN ) transfer = DBL_MIN; 00727 00728 return transfer ; 00729 }
void G4PAIModel::Initialise | ( | const G4ParticleDefinition * | , | |
const G4DataVector & | ||||
) | [virtual] |
Implements G4VEmModel.
Definition at line 160 of file G4PAIModel.cc.
References BuildLambdaVector(), BuildPAIonisationTable(), G4cout, G4endl, G4ProductionCutsTable::GetEnergyCutsVector(), G4ProductionCutsTable::GetMaterialCutsCouple(), G4Material::GetMaterialTable(), G4Material::GetNumberOfMaterials(), G4VEmModel::GetParticleChangeForLoss(), G4ParticleDefinition::GetParticleName(), G4ProductionCutsTable::GetProductionCutsTable(), and G4PAIySection::SetVerbose().
00162 { 00163 if( fVerbose > 0 && p->GetParticleName()=="proton") 00164 { 00165 G4cout<<"G4PAIModel::Initialise for "<<p->GetParticleName()<<G4endl; 00166 fPAIySection.SetVerbose(1); 00167 } 00168 else fPAIySection.SetVerbose(0); 00169 00170 if(isInitialised) { return; } 00171 isInitialised = true; 00172 00173 SetParticle(p); 00174 00175 fParticleEnergyVector = new G4PhysicsLogVector(fLowestKineticEnergy, 00176 fHighestKineticEnergy, 00177 fTotBin); 00178 00179 fParticleChange = GetParticleChangeForLoss(); 00180 00181 // Prepare initialization 00182 const G4ProductionCutsTable* theCoupleTable = 00183 G4ProductionCutsTable::GetProductionCutsTable(); 00184 const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable(); 00185 size_t numOfMat = G4Material::GetNumberOfMaterials(); 00186 size_t numRegions = fPAIRegionVector.size(); 00187 00188 for(size_t iReg = 0; iReg < numRegions; ++iReg) // region loop 00189 { 00190 const G4Region* curReg = fPAIRegionVector[iReg]; 00191 00192 for(size_t jMat = 0; jMat < numOfMat; ++jMat) // region material loop 00193 { 00194 fMaterial = (*theMaterialTable)[jMat]; 00195 fCutCouple = theCoupleTable->GetMaterialCutsCouple( fMaterial, 00196 curReg->GetProductionCuts() ); 00197 //G4cout << "Reg <" <<curReg->GetName() << "> mat <" 00198 // << fMaterial->GetName() << "> fCouple= " 00199 // << fCutCouple<<" " << p->GetParticleName() <<G4endl; 00200 if( fCutCouple ) 00201 { 00202 fMaterialCutsCoupleVector.push_back(fCutCouple); 00203 00204 fPAItransferTable = new G4PhysicsTable(fTotBin+1); 00205 fPAIdEdxTable = new G4PhysicsTable(fTotBin+1); 00206 00207 fDeltaCutInKinEnergy = 00208 (*theCoupleTable->GetEnergyCutsVector(1))[fCutCouple->GetIndex()]; 00209 00210 //ComputeSandiaPhotoAbsCof(); 00211 BuildPAIonisationTable(); 00212 00213 fPAIxscBank.push_back(fPAItransferTable); 00214 fPAIdEdxBank.push_back(fPAIdEdxTable); 00215 fdEdxTable.push_back(fdEdxVector); 00216 00217 BuildLambdaVector(); 00218 fdNdxCutTable.push_back(fdNdxCutVector); 00219 fLambdaTable.push_back(fLambdaVector); 00220 } 00221 } 00222 } 00223 }
void G4PAIModel::InitialiseMe | ( | const G4ParticleDefinition * | ) | [virtual] |
G4double G4PAIModel::MaxSecondaryEnergy | ( | const G4ParticleDefinition * | , | |
G4double | kinEnergy | |||
) | [protected, virtual] |
Reimplemented from G4VEmModel.
Definition at line 977 of file G4PAIModel.cc.
References G4ParticleDefinition::GetPDGMass().
Referenced by BuildPAIonisationTable(), and CrossSectionPerVolume().
00979 { 00980 G4double tmax = kinEnergy; 00981 if(p == fElectron) tmax *= 0.5; 00982 else if(p != fPositron) { 00983 G4double mass = p->GetPDGMass(); 00984 G4double ratio= electron_mass_c2/mass; 00985 G4double gamma= kinEnergy/mass + 1.0; 00986 tmax = 2.0*electron_mass_c2*(gamma*gamma - 1.) / 00987 (1. + 2.0*gamma*ratio + ratio*ratio); 00988 } 00989 return tmax; 00990 }
G4double G4PAIModel::SampleFluctuations | ( | const G4Material * | , | |
const G4DynamicParticle * | , | |||
G4double & | , | |||
G4double & | , | |||
G4double & | ||||
) | [virtual] |
Implements G4VEmFluctuationModel.
Definition at line 773 of file G4PAIModel.cc.
References DBL_MAX, G4UniformRand, G4DynamicParticle::GetDefinition(), GetEnergyTransfer(), G4DynamicParticle::GetKineticEnergy(), G4PhysicsVector::GetLowEdgeEnergy(), G4ParticleDefinition::GetPDGCharge(), G4ParticleDefinition::GetPDGMass(), G4InuclParticleNames::lambda, and position.
Referenced by Dispersion().
00778 { 00779 size_t jMat = 0; 00780 for(;jMat < fMaterialCutsCoupleVector.size(); ++jMat ) 00781 { 00782 if( material == fMaterialCutsCoupleVector[jMat]->GetMaterial() ) break; 00783 } 00784 if(jMat == fMaterialCutsCoupleVector.size()) return 0.0; 00785 00786 fPAItransferTable = fPAIxscBank[jMat]; 00787 fdNdxCutVector = fdNdxCutTable[jMat]; 00788 00789 G4int iTkin, iTransfer, iPlace; 00790 G4long numOfCollisions=0; 00791 00792 //G4cout<<"G4PAIModel::SampleFluctuations"<<G4endl ; 00793 //G4cout<<"in: "<<fMaterialCutsCoupleVector[jMat]->GetMaterial()->GetName()<<G4endl; 00794 G4double loss = 0.0, charge2 ; 00795 G4double stepSum = 0., stepDelta, lambda, omega; 00796 G4double position, E1, E2, W1, W2, W, dNdxCut1, dNdxCut2, meanNumber; 00797 G4bool numb = true; 00798 G4double Tkin = aParticle->GetKineticEnergy() ; 00799 G4double MassRatio = fMass/aParticle->GetDefinition()->GetPDGMass() ; 00800 G4double charge = aParticle->GetDefinition()->GetPDGCharge() ; 00801 charge2 = charge*charge ; 00802 G4double TkinScaled = Tkin*MassRatio ; 00803 00804 for(iTkin=0;iTkin<=fTotBin;iTkin++) 00805 { 00806 if(TkinScaled < fParticleEnergyVector->GetLowEdgeEnergy(iTkin)) break ; 00807 } 00808 iPlace = iTkin - 1 ; 00809 if(iPlace < 0) iPlace = 0; 00810 else if(iPlace >= fTotBin) iPlace = fTotBin - 1; 00811 //G4cout<<"from search, iPlace = "<<iPlace<<G4endl ; 00812 dNdxCut1 = (*fdNdxCutVector)(iPlace) ; 00813 //G4cout<<"dNdxCut1 = "<<dNdxCut1<<G4endl ; 00814 00815 if(iTkin == fTotBin) // Fermi plato, try from left 00816 { 00817 meanNumber =((*(*fPAItransferTable)(iPlace))(0)-dNdxCut1)*step*charge2; 00818 if(meanNumber < 0.) meanNumber = 0. ; 00819 // numOfCollisions = RandPoisson::shoot(meanNumber) ; 00820 // numOfCollisions = G4Poisson(meanNumber) ; 00821 if( meanNumber > 0.) lambda = step/meanNumber; 00822 else lambda = DBL_MAX; 00823 while(numb) 00824 { 00825 stepDelta = CLHEP::RandExponential::shoot(lambda); 00826 stepSum += stepDelta; 00827 if(stepSum >= step) break; 00828 numOfCollisions++; 00829 } 00830 //G4cout<<"##1 numOfCollisions = "<<numOfCollisions<<G4endl ; 00831 00832 while(numOfCollisions) 00833 { 00834 position = dNdxCut1+ 00835 ((*(*fPAItransferTable)(iPlace))(0)-dNdxCut1)*G4UniformRand() ; 00836 00837 for( iTransfer = 0; 00838 iTransfer < G4int((*fPAItransferTable)(iPlace)->GetVectorLength()); iTransfer++ ) 00839 { 00840 if(position >= (*(*fPAItransferTable)(iPlace))(iTransfer)) break ; 00841 } 00842 omega = GetEnergyTransfer(iPlace,position,iTransfer); 00843 //G4cout<<"G4PAIModel::SampleFluctuations, omega = "<<omega/keV<<" keV; "<<"\t"; 00844 loss += omega; 00845 numOfCollisions-- ; 00846 } 00847 } 00848 else 00849 { 00850 dNdxCut2 = (*fdNdxCutVector)(iPlace+1) ; 00851 //G4cout<<"dNdxCut2 = "<<dNdxCut2<< " iTkin= "<<iTkin<<" iPlace= "<<iPlace<<G4endl; 00852 00853 if(iTkin == 0) // Tkin is too small, trying from right only 00854 { 00855 meanNumber =((*(*fPAItransferTable)(iPlace+1))(0)-dNdxCut2)*step*charge2; 00856 if( meanNumber < 0. ) meanNumber = 0. ; 00857 // numOfCollisions = CLHEP::RandPoisson::shoot(meanNumber) ; 00858 // numOfCollisions = G4Poisson(meanNumber) ; 00859 if( meanNumber > 0.) lambda = step/meanNumber; 00860 else lambda = DBL_MAX; 00861 while(numb) 00862 { 00863 stepDelta = CLHEP::RandExponential::shoot(lambda); 00864 stepSum += stepDelta; 00865 if(stepSum >= step) break; 00866 numOfCollisions++; 00867 } 00868 00869 //G4cout<<"##2 numOfCollisions = "<<numOfCollisions<<G4endl ; 00870 00871 while(numOfCollisions) 00872 { 00873 position = dNdxCut2+ 00874 ((*(*fPAItransferTable)(iPlace+1))(0)-dNdxCut2)*G4UniformRand(); 00875 00876 for( iTransfer = 0; 00877 iTransfer < G4int((*fPAItransferTable)(iPlace+1)->GetVectorLength()); iTransfer++ ) 00878 { 00879 if(position >= (*(*fPAItransferTable)(iPlace+1))(iTransfer)) break ; 00880 } 00881 omega = GetEnergyTransfer(iPlace,position,iTransfer); 00882 //G4cout<<omega/keV<<"\t"; 00883 loss += omega; 00884 numOfCollisions-- ; 00885 } 00886 } 00887 else // general case: Tkin between two vectors of the material 00888 { 00889 E1 = fParticleEnergyVector->GetLowEdgeEnergy(iTkin - 1) ; 00890 E2 = fParticleEnergyVector->GetLowEdgeEnergy(iTkin) ; 00891 W = 1.0/(E2 - E1) ; 00892 W1 = (E2 - TkinScaled)*W ; 00893 W2 = (TkinScaled - E1)*W ; 00894 00895 //G4cout << fPAItransferTable->size() << G4endl; 00896 //G4cout<<"(*(*fPAItransferTable)(iPlace))(0) = "<< 00897 // (*(*fPAItransferTable)(iPlace))(0)<<G4endl ; 00898 //G4cout<<"(*(*fPAItransferTable)(iPlace+1))(0) = "<< 00899 // (*(*fPAItransferTable)(iPlace+1))(0)<<G4endl ; 00900 00901 meanNumber=( ((*(*fPAItransferTable)(iPlace))(0)-dNdxCut1)*W1 + 00902 ((*(*fPAItransferTable)(iPlace+1))(0)-dNdxCut2)*W2 )*step*charge2; 00903 if(meanNumber<0.0) meanNumber = 0.0; 00904 // numOfCollisions = RandPoisson::shoot(meanNumber) ; 00905 // numOfCollisions = G4Poisson(meanNumber) ; 00906 if( meanNumber > 0.) lambda = step/meanNumber; 00907 else lambda = DBL_MAX; 00908 while(numb) 00909 { 00910 stepDelta = CLHEP::RandExponential::shoot(lambda); 00911 stepSum += stepDelta; 00912 if(stepSum >= step) break; 00913 numOfCollisions++; 00914 } 00915 00916 //G4cout<<"##3 numOfCollisions = "<<numOfCollisions<<endl ; 00917 00918 while(numOfCollisions) 00919 { 00920 position = dNdxCut1*W1 + dNdxCut2*W2 + 00921 ( ( (*(*fPAItransferTable)(iPlace))(0)-dNdxCut1 )*W1 + 00922 dNdxCut2+ 00923 ( (*(*fPAItransferTable)(iPlace+1))(0)-dNdxCut2 )*W2 )*G4UniformRand(); 00924 00925 // G4cout<<position<<"\t" ; 00926 00927 for( iTransfer = 0; 00928 iTransfer < G4int((*fPAItransferTable)(iPlace)->GetVectorLength()); iTransfer++ ) 00929 { 00930 if( position >= 00931 ( (*(*fPAItransferTable)(iPlace))(iTransfer)*W1 + 00932 (*(*fPAItransferTable)(iPlace+1))(iTransfer)*W2) ) 00933 { 00934 break ; 00935 } 00936 } 00937 omega = GetEnergyTransfer(iPlace,position,iTransfer); 00938 // G4cout<<omega/keV<<"\t"; 00939 loss += omega; 00940 numOfCollisions-- ; 00941 } 00942 } 00943 } 00944 //G4cout<<"PAIModel AlongStepLoss = "<<loss/keV<<" keV, on step = " 00945 // <<step/mm<<" mm"<<G4endl ; 00946 if(loss > Tkin) loss=Tkin; 00947 if(loss < 0. ) loss = 0.; 00948 return loss ; 00949 00950 }
void G4PAIModel::SampleSecondaries | ( | std::vector< G4DynamicParticle * > * | , | |
const G4MaterialCutsCouple * | , | |||
const G4DynamicParticle * | , | |||
G4double | tmin, | |||
G4double | maxEnergy | |||
) | [virtual] |
Implements G4VEmModel.
Definition at line 566 of file G4PAIModel.cc.
References G4Electron::Electron(), G4cout, G4endl, G4UniformRand, G4DynamicParticle::GetKineticEnergy(), G4DynamicParticle::GetMass(), G4DynamicParticle::GetMomentumDirection(), GetPostStepTransfer(), G4VEmModel::MaxSecondaryKinEnergy(), G4DynamicParticle::SetDefinition(), G4DynamicParticle::SetKineticEnergy(), G4DynamicParticle::SetMomentumDirection(), G4ParticleChangeForLoss::SetProposedKineticEnergy(), and G4ParticleChangeForLoss::SetProposedMomentumDirection().
00571 { 00572 size_t jMat = 0; 00573 for(;jMat < fMaterialCutsCoupleVector.size(); ++jMat ) 00574 { 00575 if( matCC == fMaterialCutsCoupleVector[jMat] ) break; 00576 } 00577 if(jMat == fMaterialCutsCoupleVector.size()) return; 00578 00579 fPAItransferTable = fPAIxscBank[jMat]; 00580 fdNdxCutVector = fdNdxCutTable[jMat]; 00581 00582 G4double tmax = std::min(MaxSecondaryKinEnergy(dp), maxEnergy); 00583 if( tmin >= tmax && fVerbose > 0) 00584 { 00585 G4cout<<"G4PAIModel::SampleSecondary: tmin >= tmax "<<G4endl; 00586 } 00587 G4ThreeVector direction= dp->GetMomentumDirection(); 00588 G4double particleMass = dp->GetMass(); 00589 G4double kineticEnergy = dp->GetKineticEnergy(); 00590 00591 G4double massRatio = fMass/particleMass; 00592 G4double scaledTkin = kineticEnergy*massRatio; 00593 G4double totalEnergy = kineticEnergy + particleMass; 00594 G4double pSquare = kineticEnergy*(totalEnergy+particleMass); 00595 00596 G4double deltaTkin = GetPostStepTransfer(scaledTkin); 00597 00598 // G4cout<<"G4PAIModel::SampleSecondaries; deltaKIn = "<<deltaTkin/keV<<" keV "<<G4endl; 00599 00600 if( deltaTkin <= 0. && fVerbose > 0) 00601 { 00602 G4cout<<"G4PAIModel::SampleSecondary e- deltaTkin = "<<deltaTkin<<G4endl; 00603 } 00604 if( deltaTkin <= 0.) return; 00605 00606 if( deltaTkin > tmax) deltaTkin = tmax; 00607 00608 G4double deltaTotalMomentum = sqrt(deltaTkin*(deltaTkin + 2. * electron_mass_c2 )); 00609 G4double totalMomentum = sqrt(pSquare); 00610 G4double costheta = deltaTkin*(totalEnergy + electron_mass_c2) 00611 /(deltaTotalMomentum * totalMomentum); 00612 00613 if( costheta > 0.99999 ) costheta = 0.99999; 00614 G4double sintheta = 0.0; 00615 G4double sin2 = 1. - costheta*costheta; 00616 if( sin2 > 0.) sintheta = sqrt(sin2); 00617 00618 // direction of the delta electron 00619 G4double phi = twopi*G4UniformRand(); 00620 G4double dirx = sintheta*cos(phi), diry = sintheta*sin(phi), dirz = costheta; 00621 00622 G4ThreeVector deltaDirection(dirx,diry,dirz); 00623 deltaDirection.rotateUz(direction); 00624 deltaDirection.unit(); 00625 00626 // primary change 00627 kineticEnergy -= deltaTkin; 00628 G4ThreeVector dir = totalMomentum*direction - deltaTotalMomentum*deltaDirection; 00629 direction = dir.unit(); 00630 fParticleChange->SetProposedKineticEnergy(kineticEnergy); 00631 fParticleChange->SetProposedMomentumDirection(direction); 00632 00633 // create G4DynamicParticle object for e- delta ray 00634 G4DynamicParticle* deltaRay = new G4DynamicParticle; 00635 deltaRay->SetDefinition(G4Electron::Electron()); 00636 deltaRay->SetKineticEnergy( deltaTkin ); // !!! trick for last steps /2.0 ??? 00637 deltaRay->SetMomentumDirection(deltaDirection); 00638 00639 vdp->push_back(deltaRay); 00640 }
void G4PAIModel::SetVerboseLevel | ( | G4int | verbose | ) | [inline] |