#include <G4AdjointBremsstrahlungModel.hh>
Inheritance diagram for G4AdjointBremsstrahlungModel:
Definition at line 63 of file G4AdjointBremsstrahlungModel.hh.
G4AdjointBremsstrahlungModel::G4AdjointBremsstrahlungModel | ( | G4VEmModel * | aModel | ) |
Definition at line 46 of file G4AdjointBremsstrahlungModel.cc.
References G4EmModelManager::AddEmModel(), G4AdjointElectron::AdjointElectron(), G4AdjointGamma::AdjointGamma(), G4Electron::Electron(), G4VEmAdjointModel::second_part_of_same_type, G4VEmAdjointModel::SetApplyCutInRange(), G4VEmAdjointModel::SetUseMatrix(), G4VEmAdjointModel::SetUseMatrixPerElement(), G4VEmAdjointModel::theAdjEquivOfDirectPrimPartDef, G4VEmAdjointModel::theAdjEquivOfDirectSecondPartDef, G4VEmAdjointModel::theDirectEMModel, and G4VEmAdjointModel::theDirectPrimaryPartDef.
00046 : 00047 G4VEmAdjointModel("AdjointeBremModel") 00048 { 00049 SetUseMatrix(false); 00050 SetUseMatrixPerElement(false); 00051 00052 theDirectStdBremModel = aModel; 00053 theDirectEMModel=theDirectStdBremModel; 00054 theEmModelManagerForFwdModels = new G4EmModelManager(); 00055 isDirectModelInitialised = false; 00056 G4VEmFluctuationModel* f=0; 00057 G4Region* r=0; 00058 theEmModelManagerForFwdModels->AddEmModel(1, theDirectStdBremModel, f, r); 00059 00060 SetApplyCutInRange(true); 00061 highKinEnergy= 100.*TeV; 00062 lowKinEnergy = 1.0*keV; 00063 00064 lastCZ =0.; 00065 00066 00067 theAdjEquivOfDirectPrimPartDef =G4AdjointElectron::AdjointElectron(); 00068 theAdjEquivOfDirectSecondPartDef=G4AdjointGamma::AdjointGamma(); 00069 theDirectPrimaryPartDef=G4Electron::Electron(); 00070 second_part_of_same_type=false; 00071 00072 /*UsePenelopeModel=false; 00073 if (UsePenelopeModel) { 00074 G4PenelopeBremsstrahlungModel* thePenelopeModel = new G4PenelopeBremsstrahlungModel(G4Electron::Electron(),"PenelopeBrem"); 00075 theEmModelManagerForFwdModels = new G4EmModelManager(); 00076 isPenelopeModelInitialised = false; 00077 G4VEmFluctuationModel* f=0; 00078 G4Region* r=0; 00079 theDirectEMModel=thePenelopeModel; 00080 theEmModelManagerForFwdModels->AddEmModel(1, thePenelopeModel, f, r); 00081 } 00082 */ 00083 00084 00085 00086 }
G4AdjointBremsstrahlungModel::G4AdjointBremsstrahlungModel | ( | ) |
Definition at line 89 of file G4AdjointBremsstrahlungModel.cc.
References G4EmModelManager::AddEmModel(), G4AdjointElectron::AdjointElectron(), G4AdjointGamma::AdjointGamma(), G4Electron::Electron(), G4VEmAdjointModel::second_part_of_same_type, G4VEmAdjointModel::SetApplyCutInRange(), G4VEmAdjointModel::SetUseMatrix(), G4VEmAdjointModel::SetUseMatrixPerElement(), G4VEmAdjointModel::theAdjEquivOfDirectPrimPartDef, G4VEmAdjointModel::theAdjEquivOfDirectSecondPartDef, G4VEmAdjointModel::theDirectEMModel, and G4VEmAdjointModel::theDirectPrimaryPartDef.
00089 : 00090 G4VEmAdjointModel("AdjointeBremModel") 00091 { 00092 SetUseMatrix(false); 00093 SetUseMatrixPerElement(false); 00094 00095 theDirectStdBremModel = new G4SeltzerBergerModel(); 00096 theDirectEMModel=theDirectStdBremModel; 00097 theEmModelManagerForFwdModels = new G4EmModelManager(); 00098 isDirectModelInitialised = false; 00099 G4VEmFluctuationModel* f=0; 00100 G4Region* r=0; 00101 theEmModelManagerForFwdModels->AddEmModel(1, theDirectStdBremModel, f, r); 00102 // theDirectPenelopeBremModel =0; 00103 SetApplyCutInRange(true); 00104 highKinEnergy= 1.*GeV; 00105 lowKinEnergy = 1.0*keV; 00106 lastCZ =0.; 00107 theAdjEquivOfDirectPrimPartDef =G4AdjointElectron::AdjointElectron(); 00108 theAdjEquivOfDirectSecondPartDef=G4AdjointGamma::AdjointGamma(); 00109 theDirectPrimaryPartDef=G4Electron::Electron(); 00110 second_part_of_same_type=false; 00111 00112 }
G4AdjointBremsstrahlungModel::~G4AdjointBremsstrahlungModel | ( | ) |
Definition at line 115 of file G4AdjointBremsstrahlungModel.cc.
00116 {if (theDirectStdBremModel) delete theDirectStdBremModel; 00117 if (theEmModelManagerForFwdModels) delete theEmModelManagerForFwdModels; 00118 }
G4double G4AdjointBremsstrahlungModel::AdjointCrossSection | ( | const G4MaterialCutsCouple * | aCouple, | |
G4double | primEnergy, | |||
G4bool | IsScatProjToProjCase | |||
) | [virtual] |
Reimplemented from G4VEmAdjointModel.
Definition at line 399 of file G4AdjointBremsstrahlungModel.cc.
References G4VEmAdjointModel::AdjointCrossSection(), G4VEmModel::CrossSectionPerVolume(), G4VEmAdjointModel::currentTcutForDirectSecond, G4VEmAdjointModel::DefineCurrentMaterial(), G4Electron::Electron(), G4Gamma::Gamma(), G4VEmAdjointModel::GetSecondAdjEnergyMaxForProdToProjCase(), G4VEmAdjointModel::GetSecondAdjEnergyMaxForScatProjToProjCase(), G4VEmAdjointModel::GetSecondAdjEnergyMinForProdToProjCase(), G4VEmAdjointModel::GetSecondAdjEnergyMinForScatProjToProjCase(), G4EmModelManager::Initialise(), G4VEmAdjointModel::theDirectEMModel, G4VEmAdjointModel::theDirectPrimaryPartDef, and G4VEmAdjointModel::UseMatrix.
Referenced by GetAdjointCrossSection().
00402 { if (!isDirectModelInitialised) { 00403 theEmModelManagerForFwdModels->Initialise(G4Electron::Electron(),G4Gamma::Gamma(),1.,0); 00404 isDirectModelInitialised =true; 00405 } 00406 if (UseMatrix) return G4VEmAdjointModel::AdjointCrossSection(aCouple,primEnergy,IsScatProjToProjCase); 00407 DefineCurrentMaterial(aCouple); 00408 G4double Cross=0.; 00409 lastCZ=theDirectEMModel->CrossSectionPerVolume(aCouple->GetMaterial(),theDirectPrimaryPartDef,100.*MeV,100.*MeV/std::exp(1.));//this give the constant above 00410 00411 if (!IsScatProjToProjCase ){ 00412 G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(primEnergy); 00413 G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(primEnergy); 00414 if (Emax_proj>Emin_proj && primEnergy > currentTcutForDirectSecond) Cross= lastCZ*std::log(Emax_proj/Emin_proj); 00415 } 00416 else { 00417 G4double Emax_proj = GetSecondAdjEnergyMaxForScatProjToProjCase(primEnergy); 00418 G4double Emin_proj = GetSecondAdjEnergyMinForScatProjToProjCase(primEnergy,currentTcutForDirectSecond); 00419 if (Emax_proj>Emin_proj) Cross= lastCZ*std::log((Emax_proj-primEnergy)*Emin_proj/Emax_proj/(Emin_proj-primEnergy)); 00420 00421 } 00422 return Cross; 00423 }
G4double G4AdjointBremsstrahlungModel::DiffCrossSectionPerVolumePrimToSecond | ( | const G4Material * | aMaterial, | |
G4double | kinEnergyProj, | |||
G4double | kinEnergyProd | |||
) | [virtual] |
Reimplemented from G4VEmAdjointModel.
Definition at line 310 of file G4AdjointBremsstrahlungModel.cc.
References DiffCrossSectionPerVolumePrimToSecondApproximated2(), G4Electron::Electron(), G4Gamma::Gamma(), and G4EmModelManager::Initialise().
Referenced by RapidSampleSecondaries().
00314 {if (!isDirectModelInitialised) { 00315 theEmModelManagerForFwdModels->Initialise(G4Electron::Electron(),G4Gamma::Gamma(),1.,0); 00316 isDirectModelInitialised =true; 00317 } 00318 00319 return DiffCrossSectionPerVolumePrimToSecondApproximated2(aMaterial, 00320 kinEnergyProj, 00321 kinEnergyProd); 00322 /*return G4VEmAdjointModel::DiffCrossSectionPerVolumePrimToSecond(aMaterial, 00323 kinEnergyProj, 00324 kinEnergyProd);*/ 00325 }
G4double G4AdjointBremsstrahlungModel::DiffCrossSectionPerVolumePrimToSecondApproximated1 | ( | const G4Material * | aMaterial, | |
G4double | kinEnergyProj, | |||
G4double | kinEnergyProd | |||
) |
Definition at line 329 of file G4AdjointBremsstrahlungModel.cc.
References G4VEmModel::CrossSectionPerVolume(), G4VEmAdjointModel::GetSecondAdjEnergyMaxForProdToProjCase(), G4VEmAdjointModel::GetSecondAdjEnergyMinForProdToProjCase(), G4VEmAdjointModel::theDirectEMModel, and G4VEmAdjointModel::theDirectPrimaryPartDef.
00334 { 00335 G4double dCrossEprod=0.; 00336 G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd); 00337 G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd); 00338 00339 00340 //In this approximation we consider that the secondary gammas are sampled with 1/Egamma energy distribution 00341 //This is what is applied in the discrete standard model before the rejection test that make a correction 00342 //The application of the same rejection function is not possible here. 00343 //The differentiation of the CS over Ecut does not produce neither a good differential CS. That is due to the 00344 // fact that in the discrete model the differential CS and the integrated CS are both fitted but separatly and 00345 // therefore do not allow a correct numerical differentiation of the integrated CS to get the differential one. 00346 // In the future we plan to use the brem secondary spectra from the G4Penelope implementation 00347 00348 if (kinEnergyProj>Emin_proj && kinEnergyProj<=Emax_proj){ 00349 G4double sigma=theDirectEMModel->CrossSectionPerVolume(aMaterial,theDirectPrimaryPartDef,kinEnergyProj,1.*keV); 00350 dCrossEprod=sigma/kinEnergyProd/std::log(kinEnergyProj/keV); 00351 } 00352 return dCrossEprod; 00353 00354 }
G4double G4AdjointBremsstrahlungModel::DiffCrossSectionPerVolumePrimToSecondApproximated2 | ( | const G4Material * | aMaterial, | |
G4double | kinEnergyProj, | |||
G4double | kinEnergyProd | |||
) |
Definition at line 358 of file G4AdjointBremsstrahlungModel.cc.
References C1, G4VEmModel::ComputeCrossSectionPerAtom(), G4Material::GetAtomicNumDensityVector(), G4Material::GetElementVector(), G4Material::GetNumberOfElements(), G4VEmAdjointModel::theDirectEMModel, and G4VEmAdjointModel::theDirectPrimaryPartDef.
Referenced by DiffCrossSectionPerVolumePrimToSecond().
00363 { 00364 //In this approximation we derive the direct cross section over Tcut=gamma energy, en after apply the Migdla correction factor 00365 //used in the direct model 00366 00367 G4double dCrossEprod=0.; 00368 00369 const G4ElementVector* theElementVector = material->GetElementVector(); 00370 const double* theAtomNumDensityVector = material->GetAtomicNumDensityVector(); 00371 G4double dum=0.; 00372 G4double E1=kinEnergyProd,E2=kinEnergyProd*1.001; 00373 G4double dE=E2-E1; 00374 for (size_t i=0; i<material->GetNumberOfElements(); i++) { 00375 G4double C1=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,(*theElementVector)[i]->GetZ(),dum ,E1); 00376 G4double C2=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,(*theElementVector)[i]->GetZ(),dum,E2); 00377 dCrossEprod += theAtomNumDensityVector[i] * (C1-C2)/dE; 00378 00379 } 00380 00381 //Now the Migdal correction 00382 /* 00383 G4double totalEnergy = kinEnergyProj+electron_mass_c2 ; 00384 G4double kp2 = MigdalConstant*totalEnergy*totalEnergy 00385 *(material->GetElectronDensity()); 00386 00387 00388 G4double MigdalFactor = 1./(1.+kp2/(kinEnergyProd*kinEnergyProd)); // its seems that the factor used in the CS compuation i the direct 00389 //model is different than the one used in the secondary sampling by a 00390 //factor (1.+kp2) To be checked! 00391 00392 dCrossEprod*=MigdalFactor; 00393 */ 00394 return dCrossEprod; 00395 00396 }
G4double G4AdjointBremsstrahlungModel::GetAdjointCrossSection | ( | const G4MaterialCutsCouple * | aCouple, | |
G4double | primEnergy, | |||
G4bool | IsScatProjToProjCase | |||
) | [virtual] |
Reimplemented from G4VEmAdjointModel.
Definition at line 425 of file G4AdjointBremsstrahlungModel.cc.
References AdjointCrossSection(), G4VEmModel::CrossSectionPerVolume(), G4VEmAdjointModel::GetAdjointCrossSection(), G4MaterialCutsCouple::GetMaterial(), G4VEmAdjointModel::theDirectEMModel, and G4VEmAdjointModel::theDirectPrimaryPartDef.
00428 { 00429 return AdjointCrossSection(aCouple, primEnergy,IsScatProjToProjCase); 00430 lastCZ=theDirectEMModel->CrossSectionPerVolume(aCouple->GetMaterial(),theDirectPrimaryPartDef,100.*MeV,100.*MeV/std::exp(1.));//this give the constant above 00431 return G4VEmAdjointModel::GetAdjointCrossSection(aCouple, primEnergy,IsScatProjToProjCase); 00432 00433 }
void G4AdjointBremsstrahlungModel::RapidSampleSecondaries | ( | const G4Track & | aTrack, | |
G4bool | IsScatProjToProjCase, | |||
G4ParticleChange * | fParticleChange | |||
) |
Definition at line 200 of file G4AdjointBremsstrahlungModel.cc.
References G4ParticleChange::AddSecondary(), G4VEmAdjointModel::currentMaterial, G4VEmAdjointModel::currentTcutForDirectSecond, G4VEmAdjointModel::DefineCurrentMaterial(), DiffCrossSectionPerVolumePrimToSecond(), fStopAndKill, G4UniformRand, G4AdjointCSManager::GetAdjointCSManager(), G4Track::GetDynamicParticle(), G4DynamicParticle::GetKineticEnergy(), G4Track::GetMaterialCutsCouple(), G4DynamicParticle::GetMomentumDirection(), G4ParticleDefinition::GetPDGMass(), G4AdjointCSManager::GetPostStepWeightCorrection(), G4VEmAdjointModel::GetSecondAdjEnergyMaxForProdToProjCase(), G4VEmAdjointModel::GetSecondAdjEnergyMaxForScatProjToProjCase(), G4VEmAdjointModel::GetSecondAdjEnergyMinForProdToProjCase(), G4VEmAdjointModel::GetSecondAdjEnergyMinForScatProjToProjCase(), G4DynamicParticle::GetTotalEnergy(), G4Track::GetWeight(), G4VEmAdjointModel::HighEnergyLimit, G4ParticleChange::ProposeEnergy(), G4ParticleChange::ProposeMomentumDirection(), G4VParticleChange::ProposeParentWeight(), G4VParticleChange::ProposeTrackStatus(), G4VParticleChange::SetParentWeightByProcess(), G4VParticleChange::SetSecondaryWeightByProcess(), and G4VEmAdjointModel::theAdjEquivOfDirectPrimPartDef.
Referenced by SampleSecondaries().
00203 { 00204 00205 const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle(); 00206 DefineCurrentMaterial(aTrack.GetMaterialCutsCouple()); 00207 00208 00209 G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy(); 00210 G4double adjointPrimTotalEnergy = theAdjointPrimary->GetTotalEnergy(); 00211 00212 if (adjointPrimKinEnergy>HighEnergyLimit*0.999){ 00213 return; 00214 } 00215 00216 G4double projectileKinEnergy =0.; 00217 G4double gammaEnergy=0.; 00218 G4double diffCSUsed=0.; 00219 if (!IsScatProjToProjCase){ 00220 gammaEnergy=adjointPrimKinEnergy; 00221 G4double Emax = GetSecondAdjEnergyMaxForProdToProjCase(adjointPrimKinEnergy); 00222 G4double Emin= GetSecondAdjEnergyMinForProdToProjCase(adjointPrimKinEnergy);; 00223 if (Emin>=Emax) return; 00224 projectileKinEnergy=Emin*std::pow(Emax/Emin,G4UniformRand()); 00225 diffCSUsed=lastCZ/projectileKinEnergy; 00226 00227 } 00228 else { G4double Emax = GetSecondAdjEnergyMaxForScatProjToProjCase(adjointPrimKinEnergy); 00229 G4double Emin = GetSecondAdjEnergyMinForScatProjToProjCase(adjointPrimKinEnergy,currentTcutForDirectSecond); 00230 if (Emin>=Emax) return; 00231 G4double f1=(Emin-adjointPrimKinEnergy)/Emin; 00232 G4double f2=(Emax-adjointPrimKinEnergy)/Emax/f1; 00233 //G4cout<<"f1 and f2 "<<f1<<'\t'<<f2<<G4endl; 00234 projectileKinEnergy=adjointPrimKinEnergy/(1.-f1*std::pow(f2,G4UniformRand())); 00235 gammaEnergy=projectileKinEnergy-adjointPrimKinEnergy; 00236 diffCSUsed=lastCZ*adjointPrimKinEnergy/projectileKinEnergy/gammaEnergy; 00237 00238 } 00239 00240 00241 00242 00243 //Weight correction 00244 //----------------------- 00245 //First w_corr is set to the ratio between adjoint total CS and fwd total CS 00246 G4double w_corr=G4AdjointCSManager::GetAdjointCSManager()->GetPostStepWeightCorrection(); 00247 00248 //Then another correction is needed due to the fact that a biaised differential CS has been used rather than the one consistent with the direct model 00249 //Here we consider the true diffCS as the one obtained by the numericla differentiation over Tcut of the direct CS, corrected by the Migdal term. 00250 //Basically any other differential CS diffCS could be used here (example Penelope). 00251 00252 G4double diffCS = DiffCrossSectionPerVolumePrimToSecond(currentMaterial, projectileKinEnergy, gammaEnergy); 00253 w_corr*=diffCS/diffCSUsed; 00254 00255 G4double new_weight = aTrack.GetWeight()*w_corr; 00256 fParticleChange->SetParentWeightByProcess(false); 00257 fParticleChange->SetSecondaryWeightByProcess(false); 00258 fParticleChange->ProposeParentWeight(new_weight); 00259 00260 //Kinematic 00261 //--------- 00262 G4double projectileM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass(); 00263 G4double projectileTotalEnergy = projectileM0+projectileKinEnergy; 00264 G4double projectileP2 = projectileTotalEnergy*projectileTotalEnergy - projectileM0*projectileM0; 00265 G4double projectileP = std::sqrt(projectileP2); 00266 00267 00268 //Angle of the gamma direction with the projectile taken from G4eBremsstrahlungModel 00269 //------------------------------------------------ 00270 G4double u; 00271 const G4double a1 = 0.625 , a2 = 3.*a1 , d = 27. ; 00272 00273 if (9./(9.+d) > G4UniformRand()) u = - std::log(G4UniformRand()*G4UniformRand())/a1; 00274 else u = - std::log(G4UniformRand()*G4UniformRand())/a2; 00275 00276 G4double theta = u*electron_mass_c2/projectileTotalEnergy; 00277 00278 G4double sint = std::sin(theta); 00279 G4double cost = std::cos(theta); 00280 00281 G4double phi = twopi * G4UniformRand() ; 00282 00283 G4ThreeVector projectileMomentum; 00284 projectileMomentum=G4ThreeVector(std::cos(phi)*sint,std::sin(phi)*sint,cost)*projectileP; //gamma frame 00285 if (IsScatProjToProjCase) {//the adjoint primary is the scattered e- 00286 G4ThreeVector gammaMomentum = (projectileTotalEnergy-adjointPrimTotalEnergy)*G4ThreeVector(0.,0.,1.); 00287 G4ThreeVector dirProd=projectileMomentum-gammaMomentum; 00288 G4double cost1 = std::cos(dirProd.angle(projectileMomentum)); 00289 G4double sint1 = std::sqrt(1.-cost1*cost1); 00290 projectileMomentum=G4ThreeVector(std::cos(phi)*sint1,std::sin(phi)*sint1,cost1)*projectileP; 00291 00292 } 00293 00294 projectileMomentum.rotateUz(theAdjointPrimary->GetMomentumDirection()); 00295 00296 00297 00298 if (!IsScatProjToProjCase ){ //kill the primary and add a secondary 00299 fParticleChange->ProposeTrackStatus(fStopAndKill); 00300 fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,projectileMomentum)); 00301 } 00302 else { 00303 fParticleChange->ProposeEnergy(projectileKinEnergy); 00304 fParticleChange->ProposeMomentumDirection(projectileMomentum.unit()); 00305 00306 } 00307 }
void G4AdjointBremsstrahlungModel::SampleSecondaries | ( | const G4Track & | aTrack, | |
G4bool | IsScatProjToProjCase, | |||
G4ParticleChange * | fParticleChange | |||
) | [virtual] |
Implements G4VEmAdjointModel.
Definition at line 122 of file G4AdjointBremsstrahlungModel.cc.
References G4ParticleChange::AddSecondary(), G4VEmAdjointModel::CorrectPostStepWeight(), G4VEmAdjointModel::DefineCurrentMaterial(), fStopAndKill, G4UniformRand, G4DynamicParticle::GetKineticEnergy(), G4DynamicParticle::GetMomentumDirection(), G4ParticleDefinition::GetPDGMass(), G4DynamicParticle::GetTotalEnergy(), G4VParticleChange::GetWeight(), G4VEmAdjointModel::HighEnergyLimit, G4ParticleChange::ProposeEnergy(), G4ParticleChange::ProposeMomentumDirection(), G4VParticleChange::ProposeTrackStatus(), RapidSampleSecondaries(), G4VEmAdjointModel::SampleAdjSecEnergyFromCSMatrix(), G4VEmAdjointModel::theAdjEquivOfDirectPrimPartDef, and G4VEmAdjointModel::UseMatrix.
00125 { 00126 if (!UseMatrix) return RapidSampleSecondaries(aTrack,IsScatProjToProjCase,fParticleChange); 00127 00128 const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle(); 00129 DefineCurrentMaterial(aTrack.GetMaterialCutsCouple()); 00130 00131 00132 G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy(); 00133 G4double adjointPrimTotalEnergy = theAdjointPrimary->GetTotalEnergy(); 00134 00135 if (adjointPrimKinEnergy>HighEnergyLimit*0.999){ 00136 return; 00137 } 00138 00139 G4double projectileKinEnergy = SampleAdjSecEnergyFromCSMatrix(adjointPrimKinEnergy, 00140 IsScatProjToProjCase); 00141 //Weight correction 00142 //----------------------- 00143 CorrectPostStepWeight(fParticleChange, 00144 aTrack.GetWeight(), 00145 adjointPrimKinEnergy, 00146 projectileKinEnergy, 00147 IsScatProjToProjCase); 00148 00149 00150 //Kinematic 00151 //--------- 00152 G4double projectileM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass(); 00153 G4double projectileTotalEnergy = projectileM0+projectileKinEnergy; 00154 G4double projectileP2 = projectileTotalEnergy*projectileTotalEnergy - projectileM0*projectileM0; 00155 G4double projectileP = std::sqrt(projectileP2); 00156 00157 00158 //Angle of the gamma direction with the projectile taken from G4eBremsstrahlungModel 00159 //------------------------------------------------ 00160 G4double u; 00161 const G4double a1 = 0.625 , a2 = 3.*a1 , d = 27. ; 00162 00163 if (9./(9.+d) > G4UniformRand()) u = - std::log(G4UniformRand()*G4UniformRand())/a1; 00164 else u = - std::log(G4UniformRand()*G4UniformRand())/a2; 00165 00166 G4double theta = u*electron_mass_c2/projectileTotalEnergy; 00167 00168 G4double sint = std::sin(theta); 00169 G4double cost = std::cos(theta); 00170 00171 G4double phi = twopi * G4UniformRand() ; 00172 00173 G4ThreeVector projectileMomentum; 00174 projectileMomentum=G4ThreeVector(std::cos(phi)*sint,std::sin(phi)*sint,cost)*projectileP; //gamma frame 00175 if (IsScatProjToProjCase) {//the adjoint primary is the scattered e- 00176 G4ThreeVector gammaMomentum = (projectileTotalEnergy-adjointPrimTotalEnergy)*G4ThreeVector(0.,0.,1.); 00177 G4ThreeVector dirProd=projectileMomentum-gammaMomentum; 00178 G4double cost1 = std::cos(dirProd.angle(projectileMomentum)); 00179 G4double sint1 = std::sqrt(1.-cost1*cost1); 00180 projectileMomentum=G4ThreeVector(std::cos(phi)*sint1,std::sin(phi)*sint1,cost1)*projectileP; 00181 00182 } 00183 00184 projectileMomentum.rotateUz(theAdjointPrimary->GetMomentumDirection()); 00185 00186 00187 00188 if (!IsScatProjToProjCase ){ //kill the primary and add a secondary 00189 fParticleChange->ProposeTrackStatus(fStopAndKill); 00190 fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,projectileMomentum)); 00191 } 00192 else { 00193 fParticleChange->ProposeEnergy(projectileKinEnergy); 00194 fParticleChange->ProposeMomentumDirection(projectileMomentum.unit()); 00195 00196 } 00197 }