#include <G4PAIPhotonModel.hh>
Inheritance diagram for G4PAIPhotonModel:
Definition at line 66 of file G4PAIPhotonModel.hh.
G4PAIPhotonModel::G4PAIPhotonModel | ( | const G4ParticleDefinition * | p = 0 , |
|
const G4String & | nam = "PAIPhoton" | |||
) |
Definition at line 75 of file G4PAIPhotonModel.cc.
References G4Electron::Electron(), and G4Positron::Positron().
00076 : G4VEmModel(nam),G4VEmFluctuationModel(nam), 00077 fLowestKineticEnergy(10.0*keV), 00078 fHighestKineticEnergy(100.*TeV), 00079 fTotBin(200), 00080 fMeanNumber(20), 00081 fParticle(0), 00082 fHighKinEnergy(100.*TeV), 00083 fLowKinEnergy(2.0*MeV), 00084 fTwoln10(2.0*log(10.0)), 00085 fBg2lim(0.0169), 00086 fTaulim(8.4146e-3) 00087 { 00088 fVerbose = 0; 00089 fElectron = G4Electron::Electron(); 00090 fPositron = G4Positron::Positron(); 00091 00092 fProtonEnergyVector = new G4PhysicsLogVector(fLowestKineticEnergy, 00093 fHighestKineticEnergy, 00094 fTotBin); 00095 fPAItransferTable = 0; 00096 fPAIphotonTable = 0; 00097 fPAIplasmonTable = 0; 00098 00099 fPAIdEdxTable = 0; 00100 fSandiaPhotoAbsCof = 0; 00101 fdEdxVector = 0; 00102 00103 fLambdaVector = 0; 00104 fdNdxCutVector = 0; 00105 fdNdxCutPhotonVector = 0; 00106 fdNdxCutPlasmonVector = 0; 00107 00108 fSandiaIntervalNumber = 0; 00109 fMatIndex = 0; 00110 00111 fParticleChange = 0; 00112 00113 if(p) { SetParticle(p); } 00114 else { SetParticle(fElectron); } 00115 00116 isInitialised = false; 00117 }
G4PAIPhotonModel::~G4PAIPhotonModel | ( | ) | [virtual] |
Definition at line 121 of file G4PAIPhotonModel.cc.
References G4PhysicsTable::clearAndDestroy().
00122 { 00123 if(fProtonEnergyVector) delete fProtonEnergyVector; 00124 if(fdEdxVector) delete fdEdxVector ; 00125 if ( fLambdaVector) delete fLambdaVector; 00126 if ( fdNdxCutVector) delete fdNdxCutVector; 00127 00128 if( fPAItransferTable ) 00129 { 00130 fPAItransferTable->clearAndDestroy(); 00131 delete fPAItransferTable ; 00132 } 00133 if( fPAIphotonTable ) 00134 { 00135 fPAIphotonTable->clearAndDestroy(); 00136 delete fPAIphotonTable ; 00137 } 00138 if( fPAIplasmonTable ) 00139 { 00140 fPAIplasmonTable->clearAndDestroy(); 00141 delete fPAIplasmonTable ; 00142 } 00143 if(fSandiaPhotoAbsCof) 00144 { 00145 for(G4int i=0;i<fSandiaIntervalNumber;i++) 00146 { 00147 delete[] fSandiaPhotoAbsCof[i]; 00148 } 00149 delete[] fSandiaPhotoAbsCof; 00150 } 00151 }
void G4PAIPhotonModel::BuildLambdaVector | ( | const G4MaterialCutsCouple * | matCutsCouple | ) |
Definition at line 401 of file G4PAIPhotonModel.cc.
References DBL_MAX, DBL_MIN, G4cout, G4endl, GetdNdxPhotonCut(), GetdNdxPlasmonCut(), G4GeometryTolerance::GetInstance(), G4ProductionCutsTable::GetMaterialCutsCouple(), G4ProductionCutsTable::GetProductionCutsTable(), G4GeometryTolerance::GetSurfaceTolerance(), G4ProductionCutsTable::GetTableSize(), idxG4ElectronCut, idxG4GammaCut, G4InuclParticleNames::lambda, and G4PhysicsVector::PutValue().
Referenced by Initialise().
00402 { 00403 G4int i ; 00404 G4double dNdxCut,dNdxPhotonCut,dNdxPlasmonCut, lambda; 00405 G4double kCarTolerance = G4GeometryTolerance::GetInstance() 00406 ->GetSurfaceTolerance(); 00407 00408 const G4ProductionCutsTable* theCoupleTable= 00409 G4ProductionCutsTable::GetProductionCutsTable(); 00410 00411 size_t numOfCouples = theCoupleTable->GetTableSize(); 00412 size_t jMatCC; 00413 00414 for (jMatCC = 0 ; jMatCC < numOfCouples ; jMatCC++ ) 00415 { 00416 if( matCutsCouple == theCoupleTable->GetMaterialCutsCouple(jMatCC) ) break; 00417 } 00418 if( jMatCC == numOfCouples && jMatCC > 0 ) jMatCC--; 00419 00420 const vector<G4double>* deltaCutInKineticEnergy = theCoupleTable-> 00421 GetEnergyCutsVector(idxG4ElectronCut); 00422 const vector<G4double>* photonCutInKineticEnergy = theCoupleTable-> 00423 GetEnergyCutsVector(idxG4GammaCut); 00424 00425 if (fLambdaVector) delete fLambdaVector; 00426 if (fdNdxCutVector) delete fdNdxCutVector; 00427 if (fdNdxCutPhotonVector) delete fdNdxCutPhotonVector; 00428 if (fdNdxCutPlasmonVector) delete fdNdxCutPlasmonVector; 00429 00430 fLambdaVector = new G4PhysicsLogVector( fLowestKineticEnergy, 00431 fHighestKineticEnergy, 00432 fTotBin ); 00433 fdNdxCutVector = new G4PhysicsLogVector( fLowestKineticEnergy, 00434 fHighestKineticEnergy, 00435 fTotBin ); 00436 fdNdxCutPhotonVector = new G4PhysicsLogVector( fLowestKineticEnergy, 00437 fHighestKineticEnergy, 00438 fTotBin ); 00439 fdNdxCutPlasmonVector = new G4PhysicsLogVector( fLowestKineticEnergy, 00440 fHighestKineticEnergy, 00441 fTotBin ); 00442 00443 G4double deltaCutInKineticEnergyNow = (*deltaCutInKineticEnergy)[jMatCC]; 00444 G4double photonCutInKineticEnergyNow = (*photonCutInKineticEnergy)[jMatCC]; 00445 00446 if(fVerbose > 0) 00447 { 00448 G4cout<<"PAIPhotonModel deltaCutInKineticEnergyNow = " 00449 <<deltaCutInKineticEnergyNow/keV<<" keV"<<G4endl; 00450 G4cout<<"PAIPhotonModel photonCutInKineticEnergyNow = " 00451 <<photonCutInKineticEnergyNow/keV<<" keV"<<G4endl; 00452 } 00453 for ( i = 0 ; i <= fTotBin ; i++ ) 00454 { 00455 dNdxPhotonCut = GetdNdxPhotonCut(i,photonCutInKineticEnergyNow); 00456 dNdxPlasmonCut = GetdNdxPlasmonCut(i,deltaCutInKineticEnergyNow); 00457 00458 dNdxCut = dNdxPhotonCut + dNdxPlasmonCut; 00459 lambda = dNdxCut <= DBL_MIN ? DBL_MAX: 1.0/dNdxCut; 00460 00461 if (lambda <= 1000*kCarTolerance) lambda = 1000*kCarTolerance; // Mmm ??? 00462 00463 fLambdaVector->PutValue(i, lambda); 00464 00465 fdNdxCutVector->PutValue(i, dNdxCut); 00466 fdNdxCutPhotonVector->PutValue(i, dNdxPhotonCut); 00467 fdNdxCutPlasmonVector->PutValue(i, dNdxPlasmonCut); 00468 } 00469 }
void G4PAIPhotonModel::BuildPAIonisationTable | ( | ) |
Definition at line 285 of file G4PAIPhotonModel.cc.
References DBL_MIN, G4PAIxSection::GetIntegralCerenkov(), G4PAIxSection::GetIntegralPAIdEdx(), G4PAIxSection::GetIntegralPAIxSection(), G4PAIxSection::GetIntegralPlasmon(), G4PhysicsVector::GetLowEdgeEnergy(), G4PAIxSection::GetMeanEnergyLoss(), G4PAIxSection::GetSplineEnergy(), G4PAIxSection::GetSplineSize(), G4PhysicsTable::insertAt(), MaxSecondaryEnergy(), G4PhysicsVector::PutValue(), and G4PhysicsFreeVector::PutValue().
Referenced by Initialise().
00286 { 00287 G4double LowEdgeEnergy , ionloss ; 00288 G4double /*massRatio,*/ tau, Tmax, Tmin, Tkin, deltaLow, /*gamma,*/ bg2 ; 00289 /* 00290 if( fPAItransferTable ) 00291 { 00292 fPAItransferTable->clearAndDestroy() ; 00293 delete fPAItransferTable ; 00294 } 00295 */ 00296 fPAItransferTable = new G4PhysicsTable(fTotBin); 00297 /* 00298 if( fPAIratioTable ) 00299 { 00300 fPAIratioTable->clearAndDestroy() ; 00301 delete fPAIratioTable ; 00302 } 00303 */ 00304 fPAIphotonTable = new G4PhysicsTable(fTotBin); 00305 fPAIplasmonTable = new G4PhysicsTable(fTotBin); 00306 /* 00307 if( fPAIdEdxTable ) 00308 { 00309 fPAIdEdxTable->clearAndDestroy() ; 00310 delete fPAIdEdxTable ; 00311 } 00312 */ 00313 fPAIdEdxTable = new G4PhysicsTable(fTotBin); 00314 00315 // if(fdEdxVector) delete fdEdxVector ; 00316 fdEdxVector = new G4PhysicsLogVector( fLowestKineticEnergy, 00317 fHighestKineticEnergy, 00318 fTotBin ) ; 00319 Tmin = fSandiaPhotoAbsCof[0][0] ; // low energy Sandia interval 00320 deltaLow = 100.*eV; // 0.5*eV ; 00321 00322 for (G4int i = 0 ; i <= fTotBin ; i++) //The loop for the kinetic energy 00323 { 00324 LowEdgeEnergy = fProtonEnergyVector->GetLowEdgeEnergy(i) ; 00325 tau = LowEdgeEnergy/proton_mass_c2 ; 00326 // if(tau < 0.01) tau = 0.01 ; 00327 //gamma = tau +1. ; 00328 // G4cout<<"gamma = "<<gamma<<endl ; 00329 bg2 = tau*(tau + 2. ) ; 00330 //massRatio = electron_mass_c2/proton_mass_c2 ; 00331 Tmax = MaxSecondaryEnergy(fParticle, LowEdgeEnergy); 00332 // G4cout<<"proton Tkin = "<<LowEdgeEnergy/MeV<<" MeV" 00333 // <<" Tmax = "<<Tmax/MeV<<" MeV"<<G4endl; 00334 // Tkin = DeltaCutInKineticEnergyNow ; 00335 00336 // if ( DeltaCutInKineticEnergyNow > Tmax) // was < 00337 Tkin = Tmax ; 00338 if ( Tkin < Tmin + deltaLow ) // low energy safety 00339 { 00340 Tkin = Tmin + deltaLow ; 00341 } 00342 G4PAIxSection protonPAI( fMatIndex, 00343 Tkin, 00344 bg2, 00345 fSandiaPhotoAbsCof, 00346 fSandiaIntervalNumber ) ; 00347 00348 00349 // G4cout<<"ionloss = "<<ionloss*cm/keV<<" keV/cm"<<endl ; 00350 // G4cout<<"n1 = "<<protonPAI.GetIntegralPAIxSection(1)*cm<<" 1/cm"<<endl ; 00351 // G4cout<<"protonPAI.GetSplineSize() = "<< 00352 // protonPAI.GetSplineSize()<<G4endl<<G4endl ; 00353 00354 G4PhysicsFreeVector* transferVector = new 00355 G4PhysicsFreeVector(protonPAI.GetSplineSize()) ; 00356 G4PhysicsFreeVector* photonVector = new 00357 G4PhysicsFreeVector(protonPAI.GetSplineSize()) ; 00358 G4PhysicsFreeVector* plasmonVector = new 00359 G4PhysicsFreeVector(protonPAI.GetSplineSize()) ; 00360 G4PhysicsFreeVector* dEdxVector = new 00361 G4PhysicsFreeVector(protonPAI.GetSplineSize()) ; 00362 00363 for( G4int k = 0 ; k < protonPAI.GetSplineSize() ; k++ ) 00364 { 00365 transferVector->PutValue( k , 00366 protonPAI.GetSplineEnergy(k+1), 00367 protonPAI.GetIntegralPAIxSection(k+1) ) ; 00368 photonVector->PutValue( k , 00369 protonPAI.GetSplineEnergy(k+1), 00370 protonPAI.GetIntegralCerenkov(k+1) ) ; 00371 plasmonVector->PutValue( k , 00372 protonPAI.GetSplineEnergy(k+1), 00373 protonPAI.GetIntegralPlasmon(k+1) ) ; 00374 dEdxVector->PutValue( k , 00375 protonPAI.GetSplineEnergy(k+1), 00376 protonPAI.GetIntegralPAIdEdx(k+1) ) ; 00377 } 00378 ionloss = protonPAI.GetMeanEnergyLoss() ; // total <dE/dx> 00379 if ( ionloss <= 0.) ionloss = DBL_MIN ; 00380 fdEdxVector->PutValue(i,ionloss) ; 00381 00382 fPAItransferTable->insertAt(i,transferVector) ; 00383 fPAIphotonTable->insertAt(i,photonVector) ; 00384 fPAIplasmonTable->insertAt(i,plasmonVector) ; 00385 fPAIdEdxTable->insertAt(i,dEdxVector) ; 00386 00387 } // end of Tkin loop 00388 // theLossTable->insert(fdEdxVector); 00389 // end of material loop 00390 // G4cout<<"G4PAIonisation::BuildPAIonisationTable() have been called"<<G4endl ; 00391 // G4cout<<"G4PAIonisation::BuildLossTable() have been called"<<G4endl ; 00392 }
G4double G4PAIPhotonModel::ComputeDEDXPerVolume | ( | const G4Material * | , | |
const G4ParticleDefinition * | , | |||
G4double | kineticEnergy, | |||
G4double | cutEnergy | |||
) | [virtual] |
Reimplemented from G4VEmModel.
Definition at line 650 of file G4PAIPhotonModel.cc.
References G4VEmModel::CurrentCouple(), GetdEdxCut(), G4ParticleDefinition::GetPDGCharge(), and G4ParticleDefinition::GetPDGMass().
00654 { 00655 G4int iTkin,iPlace; 00656 size_t jMat; 00657 00658 //G4double cut = std::min(MaxSecondaryEnergy(p, kineticEnergy), cutEnergy); 00659 G4double cut = cutEnergy; 00660 00661 G4double particleMass = p->GetPDGMass(); 00662 G4double scaledTkin = kineticEnergy*proton_mass_c2/particleMass; 00663 G4double charge = p->GetPDGCharge()/eplus; 00664 G4double charge2 = charge*charge; 00665 G4double dEdx = 0.; 00666 const G4MaterialCutsCouple* matCC = CurrentCouple(); 00667 00668 for( jMat = 0 ;jMat < fMaterialCutsCoupleVector.size() ; ++jMat ) 00669 { 00670 if( matCC == fMaterialCutsCoupleVector[jMat] ) break; 00671 } 00672 if(jMat == fMaterialCutsCoupleVector.size() && jMat > 0) jMat--; 00673 00674 fPAIdEdxTable = fPAIdEdxBank[jMat]; 00675 fdEdxVector = fdEdxTable[jMat]; 00676 for(iTkin = 0 ; iTkin <= fTotBin ; iTkin++) 00677 { 00678 if(scaledTkin < fProtonEnergyVector->GetLowEdgeEnergy(iTkin)) break ; 00679 } 00680 iPlace = iTkin - 1; 00681 if(iPlace < 0) iPlace = 0; 00682 dEdx = charge2*( (*fdEdxVector)(iPlace) - GetdEdxCut(iPlace,cut) ) ; 00683 00684 if( dEdx < 0.) dEdx = 0.; 00685 return dEdx; 00686 }
void G4PAIPhotonModel::ComputeSandiaPhotoAbsCof | ( | ) |
Definition at line 237 of file G4PAIPhotonModel.cc.
References G4Material::GetMaterialTable(), G4SandiaTable::GetPhotoAbsorpCof(), G4SandiaTable::SandiaIntervals(), and G4SandiaTable::SandiaMixing().
Referenced by Initialise().
00238 { 00239 G4int i, j, numberOfElements ; 00240 const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable(); 00241 00242 G4SandiaTable thisMaterialSandiaTable(fMatIndex) ; 00243 numberOfElements = (*theMaterialTable)[fMatIndex]-> 00244 GetNumberOfElements(); 00245 G4int* thisMaterialZ = new G4int[numberOfElements] ; 00246 00247 for(i=0;i<numberOfElements;i++) 00248 { 00249 thisMaterialZ[i] = 00250 (G4int)(*theMaterialTable)[fMatIndex]->GetElement(i)->GetZ() ; 00251 } 00252 fSandiaIntervalNumber = thisMaterialSandiaTable.SandiaIntervals 00253 (thisMaterialZ,numberOfElements) ; 00254 00255 fSandiaIntervalNumber = thisMaterialSandiaTable.SandiaMixing 00256 ( thisMaterialZ , 00257 (*theMaterialTable)[fMatIndex]->GetFractionVector() , 00258 numberOfElements,fSandiaIntervalNumber) ; 00259 00260 fSandiaPhotoAbsCof = new G4double*[fSandiaIntervalNumber] ; 00261 00262 for(i=0;i<fSandiaIntervalNumber;i++) fSandiaPhotoAbsCof[i] = new G4double[5] ; 00263 00264 for( i = 0 ; i < fSandiaIntervalNumber ; i++ ) 00265 { 00266 fSandiaPhotoAbsCof[i][0] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i+1,0) ; 00267 00268 for( j = 1; j < 5 ; j++ ) 00269 { 00270 fSandiaPhotoAbsCof[i][j] = thisMaterialSandiaTable. 00271 GetPhotoAbsorpCof(i+1,j)* 00272 (*theMaterialTable)[fMatIndex]->GetDensity() ; 00273 } 00274 } 00275 delete[] thisMaterialZ ; 00276 }
G4double G4PAIPhotonModel::CrossSectionPerVolume | ( | const G4Material * | , | |
const G4ParticleDefinition * | , | |||
G4double | kineticEnergy, | |||
G4double | cutEnergy, | |||
G4double | maxEnergy | |||
) | [virtual] |
Reimplemented from G4VEmModel.
Definition at line 690 of file G4PAIPhotonModel.cc.
References G4VEmModel::CurrentCouple(), GetdNdxPhotonCut(), GetdNdxPlasmonCut(), G4ProductionCutsTable::GetMaterialCutsCouple(), G4ParticleDefinition::GetPDGCharge(), G4ParticleDefinition::GetPDGMass(), G4ProductionCutsTable::GetProductionCutsTable(), G4ProductionCutsTable::GetTableSize(), idxG4GammaCut, and MaxSecondaryEnergy().
00695 { 00696 G4int iTkin,iPlace; 00697 size_t jMat, jMatCC; 00698 G4double tmax = std::min(MaxSecondaryEnergy(p, kineticEnergy), maxEnergy); 00699 if(cutEnergy >= tmax) return 0.0; 00700 G4double particleMass = p->GetPDGMass(); 00701 G4double scaledTkin = kineticEnergy*proton_mass_c2/particleMass; 00702 G4double charge = p->GetPDGCharge(); 00703 G4double charge2 = charge*charge, cross, cross1, cross2; 00704 G4double photon1, photon2, plasmon1, plasmon2; 00705 00706 const G4MaterialCutsCouple* matCC = CurrentCouple(); 00707 00708 const G4ProductionCutsTable* theCoupleTable= 00709 G4ProductionCutsTable::GetProductionCutsTable(); 00710 00711 size_t numOfCouples = theCoupleTable->GetTableSize(); 00712 00713 for (jMatCC = 0 ; jMatCC < numOfCouples ; jMatCC++ ) 00714 { 00715 if( matCC == theCoupleTable->GetMaterialCutsCouple(jMatCC) ) break; 00716 } 00717 if( jMatCC == numOfCouples && jMatCC > 0 ) jMatCC--; 00718 00719 const vector<G4double>* photonCutInKineticEnergy = theCoupleTable-> 00720 GetEnergyCutsVector(idxG4GammaCut); 00721 00722 G4double photonCut = (*photonCutInKineticEnergy)[jMatCC] ; 00723 00724 for( jMat = 0 ;jMat < fMaterialCutsCoupleVector.size() ; ++jMat ) 00725 { 00726 if( matCC == fMaterialCutsCoupleVector[jMat] ) break; 00727 } 00728 if(jMat == fMaterialCutsCoupleVector.size() && jMat > 0) jMat--; 00729 00730 fPAItransferTable = fPAIxscBank[jMat]; 00731 fPAIphotonTable = fPAIphotonBank[jMat]; 00732 fPAIplasmonTable = fPAIplasmonBank[jMat]; 00733 00734 for(iTkin = 0 ; iTkin <= fTotBin ; iTkin++) 00735 { 00736 if(scaledTkin < fProtonEnergyVector->GetLowEdgeEnergy(iTkin)) break ; 00737 } 00738 iPlace = iTkin - 1; 00739 if(iPlace < 0) iPlace = 0; 00740 00741 // G4cout<<"iPlace = "<<iPlace<<"; tmax = " 00742 // <<tmax<<"; cutEnergy = "<<cutEnergy<<G4endl; 00743 photon1 = GetdNdxPhotonCut(iPlace,tmax); 00744 photon2 = GetdNdxPhotonCut(iPlace,photonCut); 00745 00746 plasmon1 = GetdNdxPlasmonCut(iPlace,tmax); 00747 plasmon2 = GetdNdxPlasmonCut(iPlace,cutEnergy); 00748 00749 cross1 = photon1 + plasmon1; 00750 // G4cout<<"cross1 = "<<cross1<<G4endl; 00751 cross2 = photon2 + plasmon2; 00752 // G4cout<<"cross2 = "<<cross2<<G4endl; 00753 cross = (cross2 - cross1)*charge2; 00754 // G4cout<<"cross = "<<cross<<G4endl; 00755 00756 if( cross < 0. ) cross = 0.; 00757 return cross; 00758 }
void G4PAIPhotonModel::DefineForRegion | ( | const G4Region * | r | ) | [virtual] |
G4double G4PAIPhotonModel::Dispersion | ( | const G4Material * | , | |
const G4DynamicParticle * | , | |||
G4double & | , | |||
G4double & | ||||
) | [virtual] |
Implements G4VEmFluctuationModel.
Definition at line 1240 of file G4PAIPhotonModel.cc.
References SampleFluctuations().
01244 { 01245 G4double loss, sumLoss=0., sumLoss2=0., sigma2, meanLoss=0.; 01246 for(G4int i = 0 ; i < fMeanNumber; i++) 01247 { 01248 loss = SampleFluctuations(material,aParticle,tmax,step,meanLoss); 01249 sumLoss += loss; 01250 sumLoss2 += loss*loss; 01251 } 01252 meanLoss = sumLoss/fMeanNumber; 01253 sigma2 = meanLoss*meanLoss + (sumLoss2-2*sumLoss*meanLoss)/fMeanNumber; 01254 return sigma2; 01255 }
G4double G4PAIPhotonModel::GetAlongStepTransfer | ( | G4PhysicsTable * | , | |
G4PhysicsLogVector * | , | |||
G4int | iPlace, | |||
G4double | scaledTkin, | |||
G4double | step, | |||
G4double | cof | |||
) |
Definition at line 1094 of file G4PAIPhotonModel.cc.
References DBL_MAX, G4UniformRand, GetEnergyTransfer(), G4PhysicsVector::GetLowEdgeEnergy(), G4InuclParticleNames::lambda, and position.
Referenced by SampleFluctuations().
01098 { 01099 G4int iTkin = iPlace + 1, iTransfer; 01100 G4double loss = 0., position, E1, E2, W1, W2, W, dNdxCut1, dNdxCut2, meanNumber; 01101 G4double lambda, stepDelta, stepSum=0. ; 01102 G4long numOfCollisions=0; 01103 G4bool numb = true; 01104 01105 dNdxCut1 = (*pVector)(iPlace) ; 01106 01107 // G4cout<<"iPlace = "<<iPlace<<endl ; 01108 01109 if(iTkin == fTotBin) // Fermi plato, try from left 01110 { 01111 meanNumber = ((*(*pTable)(iPlace))(0) - dNdxCut1)*cof; 01112 if(meanNumber < 0.) meanNumber = 0. ; 01113 // numOfCollisions = RandPoisson::shoot(meanNumber) ; 01114 if( meanNumber > 0.) lambda = step/meanNumber; 01115 else lambda = DBL_MAX; 01116 while(numb) 01117 { 01118 stepDelta = CLHEP::RandExponential::shoot(lambda); 01119 stepSum += stepDelta; 01120 if(stepSum >= step) break; 01121 numOfCollisions++; 01122 } 01123 01124 // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl ; 01125 01126 while(numOfCollisions) 01127 { 01128 position = dNdxCut1+ 01129 ((*(*pTable)(iPlace))(0) - dNdxCut1)*G4UniformRand() ; 01130 01131 for( iTransfer = 0; 01132 iTransfer < G4int((*pTable)(iPlace)->GetVectorLength()); iTransfer++ ) 01133 { 01134 if(position >= (*(*pTable)(iPlace))(iTransfer)) break ; 01135 } 01136 loss += GetEnergyTransfer(pTable,iPlace,position,iTransfer); 01137 numOfCollisions-- ; 01138 } 01139 } 01140 else 01141 { 01142 dNdxCut2 = (*pVector)(iPlace+1) ; 01143 01144 if(iTkin == 0) // Tkin is too small, trying from right only 01145 { 01146 meanNumber = ((*(*pTable)(iPlace+1))(0) - dNdxCut2)*cof; 01147 if( meanNumber < 0. ) meanNumber = 0. ; 01148 // numOfCollisions = CLHEP::RandPoisson::shoot(meanNumber) ; 01149 if( meanNumber > 0.) lambda = step/meanNumber; 01150 else lambda = DBL_MAX; 01151 while(numb) 01152 { 01153 stepDelta = CLHEP::RandExponential::shoot(lambda); 01154 stepSum += stepDelta; 01155 if(stepSum >= step) break; 01156 numOfCollisions++; 01157 } 01158 01159 // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl ; 01160 01161 while(numOfCollisions) 01162 { 01163 position = dNdxCut2+ 01164 ((*(*pTable)(iPlace+1))(0) - dNdxCut2)*G4UniformRand(); 01165 01166 for( iTransfer = 0; 01167 iTransfer < G4int((*pTable)(iPlace+1)->GetVectorLength()); iTransfer++ ) 01168 { 01169 if(position >= (*(*pTable)(iPlace+1))(iTransfer)) break ; 01170 } 01171 loss += GetEnergyTransfer(pTable,iPlace+1,position,iTransfer); 01172 numOfCollisions-- ; 01173 } 01174 } 01175 else // general case: Tkin between two vectors of the material 01176 { 01177 E1 = fProtonEnergyVector->GetLowEdgeEnergy(iTkin - 1) ; 01178 E2 = fProtonEnergyVector->GetLowEdgeEnergy(iTkin) ; 01179 W = 1.0/(E2 - E1) ; 01180 W1 = (E2 - scaledTkin)*W ; 01181 W2 = (scaledTkin - E1)*W ; 01182 01183 // G4cout<<"(*(*pTable)(iPlace))(0) = "<< 01184 // (*(*pTable)(iPlace))(0)<<G4endl ; 01185 // G4cout<<"(*(*pTable)(iPlace+1))(0) = "<< 01186 // (*(*pTable)(iPlace+1))(0)<<G4endl ; 01187 01188 meanNumber=( ((*(*pTable)(iPlace))(0)-dNdxCut1)*W1 + 01189 ((*(*pTable)(iPlace+1))(0)-dNdxCut2)*W2 )*cof; 01190 if(meanNumber<0.0) meanNumber = 0.0; 01191 // numOfCollisions = CLHEP::RandPoisson::shoot(meanNumber) ; 01192 if( meanNumber > 0.) lambda = step/meanNumber; 01193 else lambda = DBL_MAX; 01194 while(numb) 01195 { 01196 stepDelta = CLHEP::RandExponential::shoot(lambda); 01197 stepSum += stepDelta; 01198 if(stepSum >= step) break; 01199 numOfCollisions++; 01200 } 01201 01202 // G4cout<<"numOfCollisions = "<<numOfCollisions<<endl ; 01203 01204 while(numOfCollisions) 01205 { 01206 position = dNdxCut1*W1 + dNdxCut2*W2 + 01207 ( ( (*(*pTable)(iPlace ))(0) - dNdxCut1)*W1 + 01208 01209 ( (*(*pTable)(iPlace+1))(0) - dNdxCut2)*W2 )*G4UniformRand(); 01210 01211 // G4cout<<position<<"\t" ; 01212 01213 for( iTransfer = 0; 01214 iTransfer < G4int((*pTable)(iPlace)->GetVectorLength()); iTransfer++ ) 01215 { 01216 if( position >= 01217 ( (*(*pTable)(iPlace))(iTransfer)*W1 + 01218 (*(*pTable)(iPlace+1))(iTransfer)*W2) ) 01219 { 01220 break ; 01221 } 01222 } 01223 // loss += (*pTable)(iPlace)->GetLowEdgeEnergy(iTransfer) ; 01224 loss += GetEnergyTransfer(pTable,iPlace,position,iTransfer); 01225 numOfCollisions-- ; 01226 } 01227 } 01228 } 01229 01230 return loss ; 01231 01232 }
Definition at line 609 of file G4PAIPhotonModel.cc.
Referenced by ComputeDEDXPerVolume().
00610 { 00611 G4int iTransfer; 00612 G4double x1, x2, y1, y2, dEdxCut; 00613 // G4cout<<"iPlace = "<<iPlace<<"; "<<"transferCut = "<<transferCut<<G4endl; 00614 // G4cout<<"size = "<<G4int((*fPAIdEdxTable)(iPlace)->GetVectorLength()) 00615 // <<G4endl; 00616 for( iTransfer = 0 ; 00617 iTransfer < G4int((*fPAIdEdxTable)(iPlace)->GetVectorLength()) ; 00618 iTransfer++) 00619 { 00620 if(transferCut <= (*fPAIdEdxTable)(iPlace)->GetLowEdgeEnergy(iTransfer)) 00621 { 00622 break ; 00623 } 00624 } 00625 if ( iTransfer >= G4int((*fPAIdEdxTable)(iPlace)->GetVectorLength()) ) 00626 { 00627 iTransfer = (*fPAIdEdxTable)(iPlace)->GetVectorLength() - 1 ; 00628 } 00629 y1 = (*(*fPAIdEdxTable)(iPlace))(iTransfer-1) ; 00630 y2 = (*(*fPAIdEdxTable)(iPlace))(iTransfer) ; 00631 // G4cout<<"y1 = "<<y1<<"; "<<"y2 = "<<y2<<G4endl; 00632 x1 = (*fPAIdEdxTable)(iPlace)->GetLowEdgeEnergy(iTransfer-1) ; 00633 x2 = (*fPAIdEdxTable)(iPlace)->GetLowEdgeEnergy(iTransfer) ; 00634 // G4cout<<"x1 = "<<x1<<"; "<<"x2 = "<<x2<<G4endl; 00635 00636 if ( y1 == y2 ) dEdxCut = y2 ; 00637 else 00638 { 00639 // if ( x1 == x2 ) dEdxCut = y1 + (y2 - y1)*G4UniformRand() ; 00640 // if ( std::abs(x1-x2) <= eV ) dEdxCut = y1 + (y2 - y1)*G4UniformRand() ; 00641 if ( std::abs(x1-x2) <= eV ) dEdxCut = y1 + (y2 - y1)*0.5 ; 00642 else dEdxCut = y1 + (transferCut - x1)*(y2 - y1)/(x2 - x1) ; 00643 } 00644 // G4cout<<""<<dEdxCut<<G4endl; 00645 return dEdxCut ; 00646 }
Definition at line 476 of file G4PAIPhotonModel.cc.
00477 { 00478 G4int iTransfer; 00479 G4double x1, x2, y1, y2, dNdxCut; 00480 // G4cout<<"iPlace = "<<iPlace<<"; "<<"transferCut = "<<transferCut<<G4endl; 00481 // G4cout<<"size = "<<G4int((*fPAItransferTable)(iPlace)->GetVectorLength()) 00482 // <<G4endl; 00483 for( iTransfer = 0 ; 00484 iTransfer < G4int((*fPAItransferTable)(iPlace)->GetVectorLength()) ; 00485 iTransfer++) 00486 { 00487 if(transferCut <= (*fPAItransferTable)(iPlace)->GetLowEdgeEnergy(iTransfer)) 00488 { 00489 break ; 00490 } 00491 } 00492 if ( iTransfer >= G4int((*fPAItransferTable)(iPlace)->GetVectorLength()) ) 00493 { 00494 iTransfer = (*fPAItransferTable)(iPlace)->GetVectorLength() - 1 ; 00495 } 00496 y1 = (*(*fPAItransferTable)(iPlace))(iTransfer-1) ; 00497 y2 = (*(*fPAItransferTable)(iPlace))(iTransfer) ; 00498 // G4cout<<"y1 = "<<y1<<"; "<<"y2 = "<<y2<<G4endl; 00499 x1 = (*fPAItransferTable)(iPlace)->GetLowEdgeEnergy(iTransfer-1) ; 00500 x2 = (*fPAItransferTable)(iPlace)->GetLowEdgeEnergy(iTransfer) ; 00501 // G4cout<<"x1 = "<<x1<<"; "<<"x2 = "<<x2<<G4endl; 00502 00503 if ( y1 == y2 ) dNdxCut = y2 ; 00504 else 00505 { 00506 // if ( x1 == x2 ) dNdxCut = y1 + (y2 - y1)*G4UniformRand() ; 00507 // if ( std::abs(x1-x2) <= eV ) dNdxCut = y1 + (y2 - y1)*G4UniformRand() ; 00508 if ( std::abs(x1-x2) <= eV ) dNdxCut = y1 + (y2 - y1)*0.5 ; 00509 else dNdxCut = y1 + (transferCut - x1)*(y2 - y1)/(x2 - x1) ; 00510 } 00511 // G4cout<<""<<dNdxCut<<G4endl; 00512 return dNdxCut ; 00513 }
Definition at line 520 of file G4PAIPhotonModel.cc.
Referenced by BuildLambdaVector(), and CrossSectionPerVolume().
00521 { 00522 G4int iTransfer; 00523 G4double x1, x2, y1, y2, dNdxCut; 00524 // G4cout<<"iPlace = "<<iPlace<<"; "<<"transferCut = "<<transferCut<<G4endl; 00525 // G4cout<<"size = "<<G4int((*fPAIphotonTable)(iPlace)->GetVectorLength()) 00526 // <<G4endl; 00527 for( iTransfer = 0 ; 00528 iTransfer < G4int((*fPAIphotonTable)(iPlace)->GetVectorLength()) ; 00529 iTransfer++) 00530 { 00531 if(transferCut <= (*fPAIphotonTable)(iPlace)->GetLowEdgeEnergy(iTransfer)) 00532 { 00533 break ; 00534 } 00535 } 00536 if ( iTransfer >= G4int((*fPAIphotonTable)(iPlace)->GetVectorLength()) ) 00537 { 00538 iTransfer = (*fPAIphotonTable)(iPlace)->GetVectorLength() - 1 ; 00539 } 00540 y1 = (*(*fPAIphotonTable)(iPlace))(iTransfer-1) ; 00541 y2 = (*(*fPAIphotonTable)(iPlace))(iTransfer) ; 00542 // G4cout<<"y1 = "<<y1<<"; "<<"y2 = "<<y2<<G4endl; 00543 x1 = (*fPAIphotonTable)(iPlace)->GetLowEdgeEnergy(iTransfer-1) ; 00544 x2 = (*fPAIphotonTable)(iPlace)->GetLowEdgeEnergy(iTransfer) ; 00545 // G4cout<<"x1 = "<<x1<<"; "<<"x2 = "<<x2<<G4endl; 00546 00547 if ( y1 == y2 ) dNdxCut = y2 ; 00548 else 00549 { 00550 // if ( x1 == x2 ) dNdxCut = y1 + (y2 - y1)*G4UniformRand() ; 00551 // if ( std::abs(x1-x2) <= eV ) dNdxCut = y1 + (y2 - y1)*G4UniformRand() ; 00552 if ( std::abs(x1-x2) <= eV ) dNdxCut = y1 + (y2 - y1)*0.5 ; 00553 else dNdxCut = y1 + (transferCut - x1)*(y2 - y1)/(x2 - x1) ; 00554 } 00555 // G4cout<<""<<dNdxPhotonCut<<G4endl; 00556 return dNdxCut ; 00557 }
Definition at line 564 of file G4PAIPhotonModel.cc.
Referenced by BuildLambdaVector(), and CrossSectionPerVolume().
00565 { 00566 G4int iTransfer; 00567 G4double x1, x2, y1, y2, dNdxCut; 00568 00569 // G4cout<<"iPlace = "<<iPlace<<"; "<<"transferCut = "<<transferCut<<G4endl; 00570 // G4cout<<"size = "<<G4int((*fPAIPlasmonTable)(iPlace)->GetVectorLength()) 00571 // <<G4endl; 00572 for( iTransfer = 0 ; 00573 iTransfer < G4int((*fPAIplasmonTable)(iPlace)->GetVectorLength()) ; 00574 iTransfer++) 00575 { 00576 if(transferCut <= (*fPAIplasmonTable)(iPlace)->GetLowEdgeEnergy(iTransfer)) 00577 { 00578 break ; 00579 } 00580 } 00581 if ( iTransfer >= G4int((*fPAIplasmonTable)(iPlace)->GetVectorLength()) ) 00582 { 00583 iTransfer = (*fPAIplasmonTable)(iPlace)->GetVectorLength() - 1 ; 00584 } 00585 y1 = (*(*fPAIplasmonTable)(iPlace))(iTransfer-1) ; 00586 y2 = (*(*fPAIplasmonTable)(iPlace))(iTransfer) ; 00587 // G4cout<<"y1 = "<<y1<<"; "<<"y2 = "<<y2<<G4endl; 00588 x1 = (*fPAIplasmonTable)(iPlace)->GetLowEdgeEnergy(iTransfer-1) ; 00589 x2 = (*fPAIplasmonTable)(iPlace)->GetLowEdgeEnergy(iTransfer) ; 00590 // G4cout<<"x1 = "<<x1<<"; "<<"x2 = "<<x2<<G4endl; 00591 00592 if ( y1 == y2 ) dNdxCut = y2 ; 00593 else 00594 { 00595 // if ( x1 == x2 ) dNdxCut = y1 + (y2 - y1)*G4UniformRand() ; 00596 // if ( std::abs(x1-x2) <= eV ) dNdxCut = y1 + (y2 - y1)*G4UniformRand() ; 00597 if ( std::abs(x1-x2) <= eV ) dNdxCut = y1 + (y2 - y1)*0.5 ; 00598 else dNdxCut = y1 + (transferCut - x1)*(y2 - y1)/(x2 - x1) ; 00599 } 00600 // G4cout<<""<<dNdxPlasmonCut<<G4endl; 00601 return dNdxCut ; 00602 }
G4double G4PAIPhotonModel::GetEnergyTransfer | ( | G4PhysicsTable * | , | |
G4int | iPlace, | |||
G4double | position, | |||
G4int | iTransfer | |||
) |
Definition at line 987 of file G4PAIPhotonModel.cc.
References G4UniformRand.
Referenced by GetAlongStepTransfer(), and GetPostStepTransfer().
00989 { 00990 G4int iTransferMax; 00991 G4double x1, x2, y1, y2, energyTransfer; 00992 00993 if(iTransfer == 0) 00994 { 00995 energyTransfer = (*pTable)(iPlace)->GetLowEdgeEnergy(iTransfer); 00996 } 00997 else 00998 { 00999 iTransferMax = G4int((*pTable)(iPlace)->GetVectorLength()); 01000 01001 if ( iTransfer >= iTransferMax) iTransfer = iTransferMax - 1; 01002 01003 y1 = (*(*pTable)(iPlace))(iTransfer-1); 01004 y2 = (*(*fPAItransferTable)(iPlace))(iTransfer); 01005 01006 x1 = (*pTable)(iPlace)->GetLowEdgeEnergy(iTransfer-1); 01007 x2 = (*pTable)(iPlace)->GetLowEdgeEnergy(iTransfer); 01008 01009 if ( x1 == x2 ) energyTransfer = x2; 01010 else 01011 { 01012 if ( y1 == y2 ) energyTransfer = x1 + (x2 - x1)*G4UniformRand(); 01013 else 01014 { 01015 energyTransfer = x1 + (position - y1)*(x2 - x1)/(y2 - y1); 01016 } 01017 } 01018 } 01019 return energyTransfer; 01020 }
G4double G4PAIPhotonModel::GetPostStepTransfer | ( | G4PhysicsTable * | , | |
G4PhysicsLogVector * | , | |||
G4int | iPlace, | |||
G4double | scaledTkin | |||
) |
Definition at line 909 of file G4PAIPhotonModel.cc.
References G4UniformRand, GetEnergyTransfer(), G4PhysicsVector::GetLowEdgeEnergy(), and position.
Referenced by SampleSecondaries().
00912 { 00913 // G4cout<<"G4PAIPhotonModel::GetPostStepTransfer"<<G4endl ; 00914 00915 G4int iTkin = iPlace+1, iTransfer; 00916 G4double transfer = 0.0, position, dNdxCut1, dNdxCut2, E1, E2, W1, W2, W ; 00917 00918 dNdxCut1 = (*pVector)(iPlace) ; 00919 00920 // G4cout<<"iPlace = "<<iPlace<<endl ; 00921 00922 if(iTkin == fTotBin) // Fermi plato, try from left 00923 { 00924 position = dNdxCut1*G4UniformRand() ; 00925 00926 for( iTransfer = 0; 00927 iTransfer < G4int((*pTable)(iPlace)->GetVectorLength()); iTransfer++ ) 00928 { 00929 if(position >= (*(*pTable)(iPlace))(iTransfer)) break ; 00930 } 00931 transfer = GetEnergyTransfer(pTable,iPlace,position,iTransfer); 00932 } 00933 else 00934 { 00935 dNdxCut2 = (*pVector)(iPlace+1) ; 00936 if(iTkin == 0) // Tkin is too small, trying from right only 00937 { 00938 position = dNdxCut2*G4UniformRand() ; 00939 00940 for( iTransfer = 0; 00941 iTransfer < G4int((*pTable)(iPlace+1)->GetVectorLength()); iTransfer++ ) 00942 { 00943 if(position >= (*(*pTable)(iPlace+1))(iTransfer)) break ; 00944 } 00945 transfer = GetEnergyTransfer(pTable,iPlace+1,position,iTransfer); 00946 } 00947 else // general case: Tkin between two vectors of the material 00948 { 00949 E1 = fProtonEnergyVector->GetLowEdgeEnergy(iTkin - 1) ; 00950 E2 = fProtonEnergyVector->GetLowEdgeEnergy(iTkin) ; 00951 W = 1.0/(E2 - E1) ; 00952 W1 = (E2 - scaledTkin)*W ; 00953 W2 = (scaledTkin - E1)*W ; 00954 00955 position = ( dNdxCut1*W1 + dNdxCut2*W2 )*G4UniformRand() ; 00956 00957 // G4cout<<position<<"\t" ; 00958 00959 G4int iTrMax1, iTrMax2, iTrMax; 00960 00961 iTrMax1 = G4int((*pTable)(iPlace)->GetVectorLength()); 00962 iTrMax2 = G4int((*pTable)(iPlace+1)->GetVectorLength()); 00963 00964 if (iTrMax1 >= iTrMax2) iTrMax = iTrMax2; 00965 else iTrMax = iTrMax1; 00966 00967 for( iTransfer = 0; iTransfer < iTrMax; iTransfer++ ) 00968 { 00969 if( position >= 00970 ( (*(*pTable)(iPlace))(iTransfer)*W1 + 00971 (*(*pTable)(iPlace+1))(iTransfer)*W2) ) break ; 00972 } 00973 transfer = GetEnergyTransfer(pTable, iPlace, position, iTransfer); 00974 } 00975 } 00976 // G4cout<<"PAIPhotonModel PostStepTransfer = "<<transfer/keV<<" keV"<<G4endl ; 00977 if(transfer < 0.0 ) transfer = 0.0 ; 00978 return transfer ; 00979 }
void G4PAIPhotonModel::Initialise | ( | const G4ParticleDefinition * | , | |
const G4DataVector & | ||||
) | [virtual] |
Implements G4VEmModel.
Definition at line 169 of file G4PAIPhotonModel.cc.
References BuildLambdaVector(), BuildPAIonisationTable(), ComputeSandiaPhotoAbsCof(), G4Material::GetMaterialTable(), G4Material::GetNumberOfMaterials(), G4VEmModel::GetParticleChangeForLoss(), and G4ProductionCutsTable::GetProductionCutsTable().
00171 { 00172 // G4cout<<"G4PAIPhotonModel::Initialise for "<<p->GetParticleName()<<G4endl; 00173 if(isInitialised) { return; } 00174 isInitialised = true; 00175 00176 if(!fParticle) SetParticle(p); 00177 00178 fParticleChange = GetParticleChangeForLoss(); 00179 00180 const G4ProductionCutsTable* theCoupleTable = 00181 G4ProductionCutsTable::GetProductionCutsTable(); 00182 00183 for(size_t iReg = 0; iReg < fPAIRegionVector.size();++iReg) // region loop 00184 { 00185 const G4Region* curReg = fPAIRegionVector[iReg]; 00186 00187 vector<G4Material*>::const_iterator matIter = curReg->GetMaterialIterator(); 00188 size_t jMat; 00189 size_t numOfMat = curReg->GetNumberOfMaterials(); 00190 00191 // for(size_t jMat = 0; jMat < curReg->GetNumberOfMaterials();++jMat){} 00192 const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable(); 00193 size_t numberOfMat = G4Material::GetNumberOfMaterials(); 00194 00195 for(jMat = 0 ; jMat < numOfMat; ++jMat) // region material loop 00196 { 00197 const G4MaterialCutsCouple* matCouple = theCoupleTable-> 00198 GetMaterialCutsCouple( *matIter, curReg->GetProductionCuts() ); 00199 fMaterialCutsCoupleVector.push_back(matCouple); 00200 00201 size_t iMatGlob; 00202 for(iMatGlob = 0 ; iMatGlob < numberOfMat ; iMatGlob++ ) 00203 { 00204 if( *matIter == (*theMaterialTable)[iMatGlob]) break ; 00205 } 00206 fMatIndex = iMatGlob; 00207 00208 ComputeSandiaPhotoAbsCof(); 00209 BuildPAIonisationTable(); 00210 00211 fPAIxscBank.push_back(fPAItransferTable); 00212 fPAIphotonBank.push_back(fPAIphotonTable); 00213 fPAIplasmonBank.push_back(fPAIplasmonTable); 00214 fPAIdEdxBank.push_back(fPAIdEdxTable); 00215 fdEdxTable.push_back(fdEdxVector); 00216 00217 BuildLambdaVector(matCouple); 00218 00219 fdNdxCutTable.push_back(fdNdxCutVector); 00220 fdNdxCutPhotonTable.push_back(fdNdxCutPhotonVector); 00221 fdNdxCutPlasmonTable.push_back(fdNdxCutPlasmonVector); 00222 fLambdaTable.push_back(fLambdaVector); 00223 00224 00225 matIter++; 00226 } 00227 } 00228 }
void G4PAIPhotonModel::InitialiseMe | ( | const G4ParticleDefinition * | ) | [virtual] |
G4double G4PAIPhotonModel::MaxSecondaryEnergy | ( | const G4ParticleDefinition * | , | |
G4double | kinEnergy | |||
) | [protected, virtual] |
Reimplemented from G4VEmModel.
Definition at line 1259 of file G4PAIPhotonModel.cc.
References G4ParticleDefinition::GetPDGMass().
Referenced by BuildPAIonisationTable(), and CrossSectionPerVolume().
01261 { 01262 G4double tmax = kinEnergy; 01263 if(p == fElectron) tmax *= 0.5; 01264 else if(p != fPositron) { 01265 G4double mass = p->GetPDGMass(); 01266 G4double ratio= electron_mass_c2/mass; 01267 G4double gamma= kinEnergy/mass + 1.0; 01268 tmax = 2.0*electron_mass_c2*(gamma*gamma - 1.) / 01269 (1. + 2.0*gamma*ratio + ratio*ratio); 01270 } 01271 return tmax; 01272 }
G4double G4PAIPhotonModel::SampleFluctuations | ( | const G4Material * | , | |
const G4DynamicParticle * | , | |||
G4double & | , | |||
G4double & | , | |||
G4double & | ||||
) | [virtual] |
Implements G4VEmFluctuationModel.
Definition at line 1029 of file G4PAIPhotonModel.cc.
References GetAlongStepTransfer(), G4DynamicParticle::GetDefinition(), G4DynamicParticle::GetKineticEnergy(), G4ParticleDefinition::GetPDGCharge(), and G4ParticleDefinition::GetPDGMass().
Referenced by Dispersion().
01034 { 01035 size_t jMat; 01036 for( jMat = 0 ;jMat < fMaterialCutsCoupleVector.size() ; ++jMat ) 01037 { 01038 if( material == fMaterialCutsCoupleVector[jMat]->GetMaterial() ) break; 01039 } 01040 if(jMat == fMaterialCutsCoupleVector.size() && jMat > 0) jMat--; 01041 01042 fPAItransferTable = fPAIxscBank[jMat]; 01043 fPAIphotonTable = fPAIphotonBank[jMat]; 01044 fPAIplasmonTable = fPAIplasmonBank[jMat]; 01045 01046 fdNdxCutVector = fdNdxCutTable[jMat]; 01047 fdNdxCutPhotonVector = fdNdxCutPhotonTable[jMat]; 01048 fdNdxCutPlasmonVector = fdNdxCutPlasmonTable[jMat]; 01049 01050 G4int iTkin, iPlace ; 01051 01052 // G4cout<<"G4PAIPhotonModel::SampleFluctuations"<<G4endl ; 01053 01054 G4double loss, photonLoss, plasmonLoss, charge2 ; 01055 01056 01057 G4double Tkin = aParticle->GetKineticEnergy() ; 01058 G4double MassRatio = proton_mass_c2/aParticle->GetDefinition()->GetPDGMass() ; 01059 G4double charge = aParticle->GetDefinition()->GetPDGCharge() ; 01060 charge2 = charge*charge ; 01061 G4double scaledTkin = Tkin*MassRatio ; 01062 G4double cof = step*charge2; 01063 01064 for( iTkin = 0; iTkin <= fTotBin; iTkin++) 01065 { 01066 if(scaledTkin < fProtonEnergyVector->GetLowEdgeEnergy(iTkin)) break ; 01067 } 01068 iPlace = iTkin - 1 ; 01069 if( iPlace < 0 ) iPlace = 0; 01070 01071 photonLoss = GetAlongStepTransfer(fPAIphotonTable,fdNdxCutPhotonVector, 01072 iPlace,scaledTkin,step,cof); 01073 01074 // G4cout<<"PAIPhotonModel AlongStepPhotonLoss = "<<photonLoss/keV<<" keV"<<G4endl ; 01075 01076 plasmonLoss = GetAlongStepTransfer(fPAIplasmonTable,fdNdxCutPlasmonVector, 01077 iPlace,scaledTkin,step,cof); 01078 01079 // G4cout<<"PAIPhotonModel AlongStepPlasmonLoss = "<<plasmonLoss/keV<<" keV"<<G4endl ; 01080 01081 loss = photonLoss + plasmonLoss; 01082 01083 // G4cout<<"PAIPhotonModel AlongStepLoss = "<<loss/keV<<" keV"<<G4endl ; 01084 01085 return loss; 01086 }
void G4PAIPhotonModel::SampleSecondaries | ( | std::vector< G4DynamicParticle * > * | , | |
const G4MaterialCutsCouple * | , | |||
const G4DynamicParticle * | , | |||
G4double | tmin, | |||
G4double | maxEnergy | |||
) | [virtual] |
Implements G4VEmModel.
Definition at line 766 of file G4PAIPhotonModel.cc.
References G4Electron::Electron(), G4cout, G4endl, G4UniformRand, G4Gamma::Gamma(), G4DynamicParticle::GetKineticEnergy(), G4DynamicParticle::GetMass(), G4DynamicParticle::GetMomentumDirection(), GetPostStepTransfer(), G4VEmModel::MaxSecondaryKinEnergy(), G4DynamicParticle::SetDefinition(), G4DynamicParticle::SetKineticEnergy(), G4DynamicParticle::SetMomentumDirection(), G4ParticleChangeForLoss::SetProposedKineticEnergy(), and G4ParticleChangeForLoss::SetProposedMomentumDirection().
00771 { 00772 size_t jMat; 00773 for( jMat = 0 ;jMat < fMaterialCutsCoupleVector.size() ; ++jMat ) 00774 { 00775 if( matCC == fMaterialCutsCoupleVector[jMat] ) break; 00776 } 00777 if( jMat == fMaterialCutsCoupleVector.size() && jMat > 0 ) jMat--; 00778 00779 fPAItransferTable = fPAIxscBank[jMat]; 00780 fPAIphotonTable = fPAIphotonBank[jMat]; 00781 fPAIplasmonTable = fPAIplasmonBank[jMat]; 00782 00783 fdNdxCutVector = fdNdxCutTable[jMat]; 00784 fdNdxCutPhotonVector = fdNdxCutPhotonTable[jMat]; 00785 fdNdxCutPlasmonVector = fdNdxCutPlasmonTable[jMat]; 00786 00787 G4double tmax = min(MaxSecondaryKinEnergy(dp), maxEnergy); 00788 if( tmin >= tmax && fVerbose > 0) 00789 { 00790 G4cout<<"G4PAIPhotonModel::SampleSecondary: tmin >= tmax "<<G4endl; 00791 } 00792 00793 G4ThreeVector direction = dp->GetMomentumDirection(); 00794 G4double particleMass = dp->GetMass(); 00795 G4double kineticEnergy = dp->GetKineticEnergy(); 00796 G4double scaledTkin = kineticEnergy*fMass/particleMass; 00797 G4double totalEnergy = kineticEnergy + particleMass; 00798 G4double pSquare = kineticEnergy*(totalEnergy+particleMass); 00799 00800 G4int iTkin; 00801 for(iTkin=0;iTkin<=fTotBin;iTkin++) 00802 { 00803 if(scaledTkin < fProtonEnergyVector->GetLowEdgeEnergy(iTkin)) break ; 00804 } 00805 G4int iPlace = iTkin - 1 ; 00806 if(iPlace < 0) iPlace = 0; 00807 00808 G4double dNdxPhotonCut = (*fdNdxCutPhotonVector)(iPlace) ; 00809 G4double dNdxPlasmonCut = (*fdNdxCutPlasmonVector)(iPlace) ; 00810 G4double dNdxCut = dNdxPhotonCut + dNdxPlasmonCut; 00811 00812 G4double ratio; 00813 if (dNdxCut > 0.) ratio = dNdxPhotonCut/dNdxCut; 00814 else ratio = 0.; 00815 00816 if(ratio < G4UniformRand() ) // secondary e- 00817 { 00818 G4double deltaTkin = GetPostStepTransfer(fPAIplasmonTable, fdNdxCutPlasmonVector, 00819 iPlace, scaledTkin); 00820 00821 // G4cout<<"PAIPhotonModel PlasmonPostStepTransfer = "<<deltaTkin/keV<<" keV"<<G4endl ; 00822 00823 if( deltaTkin <= 0. ) 00824 { 00825 G4cout<<"G4PAIPhotonModel::SampleSecondary e- deltaTkin = "<<deltaTkin<<G4endl; 00826 } 00827 if( deltaTkin <= 0.) return; 00828 00829 G4double deltaTotalMomentum = sqrt(deltaTkin*(deltaTkin + 2. * electron_mass_c2 )); 00830 G4double totalMomentum = sqrt(pSquare); 00831 G4double costheta = deltaTkin*(totalEnergy + electron_mass_c2) 00832 /(deltaTotalMomentum * totalMomentum); 00833 00834 if( costheta > 0.99999 ) costheta = 0.99999; 00835 G4double sintheta = 0.0; 00836 G4double sin2 = 1. - costheta*costheta; 00837 if( sin2 > 0.) sintheta = sqrt(sin2); 00838 00839 // direction of the delta electron 00840 00841 G4double phi = twopi*G4UniformRand(); 00842 G4double dirx = sintheta*cos(phi), diry = sintheta*sin(phi), dirz = costheta; 00843 00844 G4ThreeVector deltaDirection(dirx,diry,dirz); 00845 deltaDirection.rotateUz(direction); 00846 00847 // primary change 00848 00849 kineticEnergy -= deltaTkin; 00850 G4ThreeVector dir = totalMomentum*direction - deltaTotalMomentum*deltaDirection; 00851 direction = dir.unit(); 00852 fParticleChange->SetProposedMomentumDirection(direction); 00853 00854 // create G4DynamicParticle object for e- delta ray 00855 00856 G4DynamicParticle* deltaRay = new G4DynamicParticle; 00857 deltaRay->SetDefinition(G4Electron::Electron()); 00858 deltaRay->SetKineticEnergy( deltaTkin ); 00859 deltaRay->SetMomentumDirection(deltaDirection); 00860 vdp->push_back(deltaRay); 00861 00862 } 00863 else // secondary 'Cherenkov' photon 00864 { 00865 G4double deltaTkin = GetPostStepTransfer(fPAIphotonTable, fdNdxCutPhotonVector, 00866 iPlace,scaledTkin); 00867 00868 // G4cout<<"PAIPhotonModel PhotonPostStepTransfer = "<<deltaTkin/keV<<" keV"<<G4endl ; 00869 00870 if( deltaTkin <= 0. ) 00871 { 00872 G4cout<<"G4PAIPhotonModel::SampleSecondary gamma deltaTkin = "<<deltaTkin<<G4endl; 00873 } 00874 if( deltaTkin <= 0.) return; 00875 00876 G4double costheta = 0.; // G4UniformRand(); // VG: ??? for start only 00877 G4double sintheta = sqrt((1.+costheta)*(1.-costheta)); 00878 00879 // direction of the 'Cherenkov' photon 00880 G4double phi = twopi*G4UniformRand(); 00881 G4double dirx = sintheta*cos(phi), diry = sintheta*sin(phi), dirz = costheta; 00882 00883 G4ThreeVector deltaDirection(dirx,diry,dirz); 00884 deltaDirection.rotateUz(direction); 00885 00886 // primary change 00887 kineticEnergy -= deltaTkin; 00888 00889 // create G4DynamicParticle object for photon ray 00890 00891 G4DynamicParticle* photonRay = new G4DynamicParticle; 00892 photonRay->SetDefinition( G4Gamma::Gamma() ); 00893 photonRay->SetKineticEnergy( deltaTkin ); 00894 photonRay->SetMomentumDirection(deltaDirection); 00895 00896 vdp->push_back(photonRay); 00897 } 00898 00899 fParticleChange->SetProposedKineticEnergy(kineticEnergy); 00900 }