#include <G4AdjointhIonisationModel.hh>
Inheritance diagram for G4AdjointhIonisationModel:
Public Member Functions | |
G4AdjointhIonisationModel (G4ParticleDefinition *projectileDefinition) | |
virtual | ~G4AdjointhIonisationModel () |
virtual void | SampleSecondaries (const G4Track &aTrack, G4bool IsScatProjToProjCase, G4ParticleChange *fParticleChange) |
void | RapidSampleSecondaries (const G4Track &aTrack, G4bool IsScatProjToProjCase, G4ParticleChange *fParticleChange) |
virtual G4double | DiffCrossSectionPerAtomPrimToSecond (G4double kinEnergyProj, G4double kinEnergyProd, G4double Z, G4double A=0.) |
virtual G4double | AdjointCrossSection (const G4MaterialCutsCouple *aCouple, G4double primEnergy, G4bool IsScatProjToProjCase) |
virtual G4double | GetSecondAdjEnergyMaxForScatProjToProjCase (G4double PrimAdjEnergy) |
virtual G4double | GetSecondAdjEnergyMinForScatProjToProjCase (G4double PrimAdjEnergy, G4double Tcut=0) |
virtual G4double | GetSecondAdjEnergyMaxForProdToProjCase (G4double PrimAdjEnergy) |
virtual G4double | GetSecondAdjEnergyMinForProdToProjCase (G4double PrimAdjEnergy) |
Definition at line 72 of file G4AdjointhIonisationModel.hh.
G4AdjointhIonisationModel::G4AdjointhIonisationModel | ( | G4ParticleDefinition * | projectileDefinition | ) |
Definition at line 46 of file G4AdjointhIonisationModel.cc.
References G4AdjointElectron::AdjointElectron(), G4AdjointProton::AdjointProton(), G4VEmAdjointModel::ApplyCutInRange, G4VEmAdjointModel::CS_biasing_factor, G4Proton::Proton(), G4VEmAdjointModel::second_part_of_same_type, G4VEmAdjointModel::theAdjEquivOfDirectPrimPartDef, G4VEmAdjointModel::theAdjEquivOfDirectSecondPartDef, G4VEmAdjointModel::theDirectEMModel, G4VEmAdjointModel::theDirectPrimaryPartDef, G4VEmAdjointModel::UseMatrix, G4VEmAdjointModel::UseMatrixPerElement, and G4VEmAdjointModel::UseOnlyOneMatrixForAllElements.
00046 : 00047 G4VEmAdjointModel("Adjoint_hIonisation") 00048 { 00049 00050 00051 00052 UseMatrix =true; 00053 UseMatrixPerElement = true; 00054 ApplyCutInRange = true; 00055 UseOnlyOneMatrixForAllElements = true; 00056 CS_biasing_factor =1.; 00057 second_part_of_same_type =false; 00058 00059 //The direct EM Modfel is taken has BetheBloch it is only used for the computation 00060 // of the differential cross section. 00061 //The Bragg model could be used as an alternative as it offers the same differential cross section 00062 00063 theDirectEMModel = new G4BetheBlochModel(projectileDefinition); 00064 theBraggDirectEMModel = new G4BraggModel(projectileDefinition); 00065 theAdjEquivOfDirectSecondPartDef=G4AdjointElectron::AdjointElectron(); 00066 00067 theDirectPrimaryPartDef = projectileDefinition; 00068 theAdjEquivOfDirectPrimPartDef = 0; 00069 if (projectileDefinition == G4Proton::Proton()) { 00070 theAdjEquivOfDirectPrimPartDef = G4AdjointProton::AdjointProton(); 00071 00072 } 00073 00074 DefineProjectileProperty(); 00075 }
G4AdjointhIonisationModel::~G4AdjointhIonisationModel | ( | ) | [virtual] |
G4double G4AdjointhIonisationModel::AdjointCrossSection | ( | const G4MaterialCutsCouple * | aCouple, | |
G4double | primEnergy, | |||
G4bool | IsScatProjToProjCase | |||
) | [virtual] |
Reimplemented from G4VEmAdjointModel.
Definition at line 438 of file G4AdjointhIonisationModel.cc.
References G4VEmAdjointModel::AdjointCrossSection(), G4VEmAdjointModel::currentMaterial, G4VEmAdjointModel::currentTcutForDirectSecond, G4VEmAdjointModel::DefineCurrentMaterial(), G4Material::GetElectronDensity(), GetSecondAdjEnergyMaxForProdToProjCase(), GetSecondAdjEnergyMaxForScatProjToProjCase(), GetSecondAdjEnergyMinForProdToProjCase(), GetSecondAdjEnergyMinForScatProjToProjCase(), G4VEmAdjointModel::lastCS, and G4VEmAdjointModel::UseMatrix.
00441 { 00442 if (UseMatrix) return G4VEmAdjointModel::AdjointCrossSection(aCouple,primEnergy,IsScatProjToProjCase); 00443 DefineCurrentMaterial(aCouple); 00444 00445 00446 G4double Cross=currentMaterial->GetElectronDensity()*twopi_mc2_rcl2*mass; 00447 00448 if (!IsScatProjToProjCase ){ 00449 G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(primEnergy); 00450 G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(primEnergy); 00451 if (Emax_proj>Emin_proj && primEnergy > currentTcutForDirectSecond) { 00452 Cross*=(1./Emin_proj -1./Emax_proj)/primEnergy; 00453 } 00454 else Cross=0.; 00455 00456 00457 00458 00459 00460 00461 } 00462 else { 00463 G4double Emax_proj = GetSecondAdjEnergyMaxForScatProjToProjCase(primEnergy); 00464 G4double Emin_proj = GetSecondAdjEnergyMinForScatProjToProjCase(primEnergy,currentTcutForDirectSecond); 00465 G4double diff1=Emin_proj-primEnergy; 00466 G4double diff2=Emax_proj-primEnergy; 00467 G4double t1=(1./diff1+1./Emin_proj-1./diff2-1./Emax_proj)/primEnergy; 00468 //G4double t2=2.*std::log(diff2*Emin_proj/Emax_proj/diff1)/primEnergy/primEnergy; 00469 G4double t2=2.*std::log(Emax_proj/Emin_proj)/primEnergy/primEnergy; 00470 Cross*=(t1+t2); 00471 00472 00473 } 00474 lastCS =Cross; 00475 return Cross; 00476 }
G4double G4AdjointhIonisationModel::DiffCrossSectionPerAtomPrimToSecond | ( | G4double | kinEnergyProj, | |
G4double | kinEnergyProd, | |||
G4double | Z, | |||
G4double | A = 0. | |||
) | [virtual] |
Reimplemented from G4VEmAdjointModel.
Definition at line 306 of file G4AdjointhIonisationModel.cc.
References G4VEmModel::ComputeCrossSectionPerAtom(), G4cout, G4endl, GetSecondAdjEnergyMaxForProdToProjCase(), GetSecondAdjEnergyMinForProdToProjCase(), G4VEmAdjointModel::theDirectEMModel, and G4VEmAdjointModel::theDirectPrimaryPartDef.
Referenced by RapidSampleSecondaries().
00311 {//Probably that here the Bragg Model should be also used for kinEnergyProj/nuc<2MeV 00312 00313 00314 00315 G4double dSigmadEprod=0; 00316 G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd); 00317 G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd); 00318 00319 00320 if (kinEnergyProj>Emin_proj && kinEnergyProj<=Emax_proj){ //the produced particle should have a kinetic energy smaller than the projectile 00321 G4double Tmax=kinEnergyProj; 00322 00323 G4double E1=kinEnergyProd; 00324 G4double E2=kinEnergyProd*1.000001; 00325 G4double dE=(E2-E1); 00326 G4double sigma1,sigma2; 00327 if (kinEnergyProj >2.*MeV){ 00328 sigma1=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E1,1.e20); 00329 sigma2=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E2,1.e20); 00330 } 00331 else { 00332 sigma1=theBraggDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E1,1.e20); 00333 sigma2=theBraggDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E2,1.e20); 00334 } 00335 00336 00337 dSigmadEprod=(sigma1-sigma2)/dE; 00338 if (dSigmadEprod>1.) { 00339 G4cout<<"sigma1 "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<sigma1<<G4endl; 00340 G4cout<<"sigma2 "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<sigma2<<G4endl; 00341 G4cout<<"dsigma "<<kinEnergyProj/MeV<<'\t'<<kinEnergyProd/MeV<<'\t'<<dSigmadEprod<<G4endl; 00342 00343 } 00344 00345 00346 00347 //correction of differential cross section at high energy to correct for the suppression of particle at secondary at high 00348 //energy used in the Bethe Bloch Model. This correction consist to multiply by g the probability function used 00349 //to test the rejection of a secondary 00350 //------------------------- 00351 00352 //Source code taken from G4BetheBlochModel::SampleSecondaries 00353 00354 G4double deltaKinEnergy = kinEnergyProd; 00355 00356 //Part of the taken code 00357 //---------------------- 00358 00359 00360 00361 // projectile formfactor - suppresion of high energy 00362 // delta-electron production at high energy 00363 G4double x = formfact*deltaKinEnergy; 00364 if(x > 1.e-6) { 00365 00366 00367 G4double totEnergy = kinEnergyProj + mass; 00368 G4double etot2 = totEnergy*totEnergy; 00369 G4double beta2 = kinEnergyProj*(kinEnergyProj + 2.0*mass)/etot2; 00370 G4double f; 00371 G4double f1 = 0.0; 00372 f = 1.0 - beta2*deltaKinEnergy/Tmax; 00373 if( 0.5 == spin ) { 00374 f1 = 0.5*deltaKinEnergy*deltaKinEnergy/etot2; 00375 f += f1; 00376 } 00377 G4double x1 = 1.0 + x; 00378 G4double gg = 1.0/(x1*x1); 00379 if( 0.5 == spin ) { 00380 G4double x2 = 0.5*electron_mass_c2*deltaKinEnergy/(mass*mass); 00381 gg *= (1.0 + magMoment2*(x2 - f1/f)/(1.0 + x2)); 00382 } 00383 if(gg > 1.0) { 00384 G4cout << "### G4BetheBlochModel in Adjoint Sim WARNING: g= " << g 00385 << G4endl; 00386 gg=1.; 00387 } 00388 //G4cout<<"gg"<<gg<<G4endl; 00389 dSigmadEprod*=gg; 00390 } 00391 00392 } 00393 00394 return dSigmadEprod; 00395 }
G4double G4AdjointhIonisationModel::GetSecondAdjEnergyMaxForProdToProjCase | ( | G4double | PrimAdjEnergy | ) | [virtual] |
Reimplemented from G4VEmAdjointModel.
Definition at line 491 of file G4AdjointhIonisationModel.cc.
References G4VEmAdjointModel::HighEnergyLimit.
Referenced by AdjointCrossSection(), DiffCrossSectionPerAtomPrimToSecond(), and RapidSampleSecondaries().
00492 { return HighEnergyLimit; 00493 }
G4double G4AdjointhIonisationModel::GetSecondAdjEnergyMaxForScatProjToProjCase | ( | G4double | PrimAdjEnergy | ) | [virtual] |
Reimplemented from G4VEmAdjointModel.
Definition at line 479 of file G4AdjointhIonisationModel.cc.
Referenced by AdjointCrossSection(), and RapidSampleSecondaries().
00480 { 00481 G4double Tmax=PrimAdjEnergy*one_plus_ratio_2/(one_minus_ratio_2-2.*ratio*PrimAdjEnergy/mass); 00482 return Tmax; 00483 }
G4double G4AdjointhIonisationModel::GetSecondAdjEnergyMinForProdToProjCase | ( | G4double | PrimAdjEnergy | ) | [virtual] |
Reimplemented from G4VEmAdjointModel.
Definition at line 496 of file G4AdjointhIonisationModel.cc.
Referenced by AdjointCrossSection(), DiffCrossSectionPerAtomPrimToSecond(), and RapidSampleSecondaries().
00497 { G4double Tmin= (2*PrimAdjEnergy-4*mass + std::sqrt(4.*PrimAdjEnergy*PrimAdjEnergy +16.*mass*mass + 8.*PrimAdjEnergy*mass*(1/ratio +ratio)))/4.; 00498 return Tmin; 00499 }
G4double G4AdjointhIonisationModel::GetSecondAdjEnergyMinForScatProjToProjCase | ( | G4double | PrimAdjEnergy, | |
G4double | Tcut = 0 | |||
) | [virtual] |
Reimplemented from G4VEmAdjointModel.
Definition at line 486 of file G4AdjointhIonisationModel.cc.
Referenced by AdjointCrossSection(), and RapidSampleSecondaries().
void G4AdjointhIonisationModel::RapidSampleSecondaries | ( | const G4Track & | aTrack, | |
G4bool | IsScatProjToProjCase, | |||
G4ParticleChange * | fParticleChange | |||
) |
Definition at line 160 of file G4AdjointhIonisationModel.cc.
References G4ParticleChange::AddSecondary(), G4VEmAdjointModel::currentMaterial, G4VEmAdjointModel::currentTcutForDirectSecond, G4VEmAdjointModel::DefineCurrentMaterial(), DiffCrossSectionPerAtomPrimToSecond(), fStopAndKill, G4UniformRand, G4AdjointCSManager::GetAdjointCSManager(), G4Track::GetDynamicParticle(), G4Material::GetElectronDensity(), G4DynamicParticle::GetKineticEnergy(), G4Track::GetMaterialCutsCouple(), G4DynamicParticle::GetMomentumDirection(), G4ParticleDefinition::GetPDGMass(), G4AdjointCSManager::GetPostStepWeightCorrection(), GetSecondAdjEnergyMaxForProdToProjCase(), GetSecondAdjEnergyMaxForScatProjToProjCase(), GetSecondAdjEnergyMinForProdToProjCase(), GetSecondAdjEnergyMinForScatProjToProjCase(), G4DynamicParticle::GetTotalMomentum(), G4Track::GetWeight(), G4VEmAdjointModel::HighEnergyLimit, G4VEmAdjointModel::lastCS, G4ParticleChange::ProposeEnergy(), G4ParticleChange::ProposeMomentumDirection(), G4VParticleChange::ProposeParentWeight(), G4VParticleChange::ProposeTrackStatus(), G4VParticleChange::SetParentWeightByProcess(), G4VParticleChange::SetSecondaryWeightByProcess(), G4VEmAdjointModel::theAdjEquivOfDirectPrimPartDef, and G4VEmAdjointModel::theAdjEquivOfDirectSecondPartDef.
Referenced by SampleSecondaries().
00163 { 00164 00165 const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle(); 00166 DefineCurrentMaterial(aTrack.GetMaterialCutsCouple()); 00167 00168 00169 G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy(); 00170 G4double adjointPrimP =theAdjointPrimary->GetTotalMomentum(); 00171 00172 if (adjointPrimKinEnergy>HighEnergyLimit*0.999){ 00173 return; 00174 } 00175 00176 G4double projectileKinEnergy =0.; 00177 G4double eEnergy=0.; 00178 G4double newCS=currentMaterial->GetElectronDensity()*twopi_mc2_rcl2*mass; 00179 if (!IsScatProjToProjCase){//1/E^2 distribution 00180 00181 eEnergy=adjointPrimKinEnergy; 00182 G4double Emax = GetSecondAdjEnergyMaxForProdToProjCase(adjointPrimKinEnergy); 00183 G4double Emin= GetSecondAdjEnergyMinForProdToProjCase(adjointPrimKinEnergy); 00184 if (Emin>=Emax) return; 00185 G4double a=1./Emax; 00186 G4double b=1./Emin; 00187 newCS=newCS*(b-a)/eEnergy; 00188 00189 projectileKinEnergy =1./(b- (b-a)*G4UniformRand()); 00190 00191 00192 } 00193 else { G4double Emax = GetSecondAdjEnergyMaxForScatProjToProjCase(adjointPrimKinEnergy); 00194 G4double Emin = GetSecondAdjEnergyMinForScatProjToProjCase(adjointPrimKinEnergy,currentTcutForDirectSecond); 00195 if (Emin>=Emax) return; 00196 G4double diff1=Emin-adjointPrimKinEnergy; 00197 G4double diff2=Emax-adjointPrimKinEnergy; 00198 00199 G4double t1=adjointPrimKinEnergy*(1./diff1-1./diff2); 00200 G4double t2=adjointPrimKinEnergy*(1./Emin-1./Emax); 00201 /*G4double f31=diff1/Emin; 00202 G4double f32=diff2/Emax/f31;*/ 00203 G4double t3=2.*std::log(Emax/Emin); 00204 G4double sum_t=t1+t2+t3; 00205 newCS=newCS*sum_t/adjointPrimKinEnergy/adjointPrimKinEnergy; 00206 G4double t=G4UniformRand()*sum_t; 00207 if (t <=t1 ){ 00208 G4double q= G4UniformRand()*t1/adjointPrimKinEnergy ; 00209 projectileKinEnergy =adjointPrimKinEnergy +1./(1./diff1-q); 00210 00211 } 00212 else if (t <=t2 ) { 00213 G4double q= G4UniformRand()*t2/adjointPrimKinEnergy; 00214 projectileKinEnergy =1./(1./Emin-q); 00215 } 00216 else { 00217 projectileKinEnergy=Emin*std::pow(Emax/Emin,G4UniformRand()); 00218 00219 } 00220 eEnergy=projectileKinEnergy-adjointPrimKinEnergy; 00221 00222 00223 } 00224 00225 00226 00227 G4double diffCS_perAtom_Used=twopi_mc2_rcl2*mass*adjointPrimKinEnergy/projectileKinEnergy/projectileKinEnergy/eEnergy/eEnergy; 00228 00229 00230 00231 //Weight correction 00232 //----------------------- 00233 //First w_corr is set to the ratio between adjoint total CS and fwd total CS 00234 G4double w_corr=G4AdjointCSManager::GetAdjointCSManager()->GetPostStepWeightCorrection(); 00235 00236 //G4cout<<w_corr<<G4endl; 00237 w_corr*=newCS/lastCS; 00238 //G4cout<<w_corr<<G4endl; 00239 //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 00240 //Here we consider the true diffCS as the one obtained by the numerical differentiation over Tcut of the direct CS 00241 00242 G4double diffCS = DiffCrossSectionPerAtomPrimToSecond(projectileKinEnergy, eEnergy,1,1); 00243 w_corr*=diffCS/diffCS_perAtom_Used; 00244 //G4cout<<w_corr<<G4endl; 00245 00246 G4double new_weight = aTrack.GetWeight()*w_corr; 00247 fParticleChange->SetParentWeightByProcess(false); 00248 fParticleChange->SetSecondaryWeightByProcess(false); 00249 fParticleChange->ProposeParentWeight(new_weight); 00250 00251 00252 00253 00254 //Kinematic: 00255 //we consider a two body elastic scattering for the forward processes where the projectile knock on an e- at rest and gives 00256 // him part of its energy 00257 //---------------------------------------------------------------------------------------- 00258 00259 G4double projectileM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass(); 00260 G4double projectileTotalEnergy = projectileM0+projectileKinEnergy; 00261 G4double projectileP2 = projectileTotalEnergy*projectileTotalEnergy - projectileM0*projectileM0; 00262 00263 00264 00265 //Companion 00266 //----------- 00267 G4double companionM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass(); 00268 if (IsScatProjToProjCase) { 00269 companionM0=theAdjEquivOfDirectSecondPartDef->GetPDGMass(); 00270 } 00271 G4double companionTotalEnergy =companionM0+ projectileKinEnergy-adjointPrimKinEnergy; 00272 G4double companionP2 = companionTotalEnergy*companionTotalEnergy - companionM0*companionM0; 00273 00274 00275 //Projectile momentum 00276 //-------------------- 00277 G4double P_parallel = (adjointPrimP*adjointPrimP + projectileP2 - companionP2)/(2.*adjointPrimP); 00278 G4double P_perp = std::sqrt( projectileP2 - P_parallel*P_parallel); 00279 G4ThreeVector dir_parallel=theAdjointPrimary->GetMomentumDirection(); 00280 G4double phi =G4UniformRand()*2.*3.1415926; 00281 G4ThreeVector projectileMomentum = G4ThreeVector(P_perp*std::cos(phi),P_perp*std::sin(phi),P_parallel); 00282 projectileMomentum.rotateUz(dir_parallel); 00283 00284 00285 00286 if (!IsScatProjToProjCase ){ //kill the primary and add a secondary 00287 fParticleChange->ProposeTrackStatus(fStopAndKill); 00288 fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,projectileMomentum)); 00289 //G4cout<<"projectileMomentum "<<projectileMomentum<<G4endl; 00290 } 00291 else { 00292 fParticleChange->ProposeEnergy(projectileKinEnergy); 00293 fParticleChange->ProposeMomentumDirection(projectileMomentum.unit()); 00294 } 00295 00296 00297 00298 00299 00300 00301 00302 }
void G4AdjointhIonisationModel::SampleSecondaries | ( | const G4Track & | aTrack, | |
G4bool | IsScatProjToProjCase, | |||
G4ParticleChange * | fParticleChange | |||
) | [virtual] |
Implements G4VEmAdjointModel.
Definition at line 84 of file G4AdjointhIonisationModel.cc.
References G4ParticleChange::AddSecondary(), G4VEmAdjointModel::CorrectPostStepWeight(), fStopAndKill, G4UniformRand, G4DynamicParticle::GetKineticEnergy(), G4DynamicParticle::GetMomentumDirection(), G4ParticleDefinition::GetPDGMass(), G4DynamicParticle::GetTotalMomentum(), G4VParticleChange::GetWeight(), G4VEmAdjointModel::HighEnergyLimit, G4ParticleChange::ProposeEnergy(), G4ParticleChange::ProposeMomentumDirection(), G4VParticleChange::ProposeTrackStatus(), RapidSampleSecondaries(), G4VEmAdjointModel::SampleAdjSecEnergyFromCSMatrix(), G4VEmAdjointModel::theAdjEquivOfDirectPrimPartDef, G4VEmAdjointModel::theAdjEquivOfDirectSecondPartDef, and G4VEmAdjointModel::UseMatrix.
00087 { 00088 if (!UseMatrix) return RapidSampleSecondaries(aTrack,IsScatProjToProjCase,fParticleChange); 00089 00090 const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle(); 00091 00092 //Elastic inverse scattering 00093 //--------------------------------------------------------- 00094 G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy(); 00095 G4double adjointPrimP =theAdjointPrimary->GetTotalMomentum(); 00096 00097 if (adjointPrimKinEnergy>HighEnergyLimit*0.999){ 00098 return; 00099 } 00100 00101 //Sample secondary energy 00102 //----------------------- 00103 G4double projectileKinEnergy = SampleAdjSecEnergyFromCSMatrix(adjointPrimKinEnergy, IsScatProjToProjCase); 00104 CorrectPostStepWeight(fParticleChange, 00105 aTrack.GetWeight(), 00106 adjointPrimKinEnergy, 00107 projectileKinEnergy, 00108 IsScatProjToProjCase); //Caution !!!this weight correction should be always applied 00109 00110 00111 //Kinematic: 00112 //we consider a two body elastic scattering for the forward processes where the projectile knock on an e- at rest and gives 00113 // him part of its energy 00114 //---------------------------------------------------------------------------------------- 00115 00116 G4double projectileM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass(); 00117 G4double projectileTotalEnergy = projectileM0+projectileKinEnergy; 00118 G4double projectileP2 = projectileTotalEnergy*projectileTotalEnergy - projectileM0*projectileM0; 00119 00120 00121 00122 //Companion 00123 //----------- 00124 G4double companionM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass(); 00125 if (IsScatProjToProjCase) { 00126 companionM0=theAdjEquivOfDirectSecondPartDef->GetPDGMass(); 00127 } 00128 G4double companionTotalEnergy =companionM0+ projectileKinEnergy-adjointPrimKinEnergy; 00129 G4double companionP2 = companionTotalEnergy*companionTotalEnergy - companionM0*companionM0; 00130 00131 00132 //Projectile momentum 00133 //-------------------- 00134 G4double P_parallel = (adjointPrimP*adjointPrimP + projectileP2 - companionP2)/(2.*adjointPrimP); 00135 G4double P_perp = std::sqrt( projectileP2 - P_parallel*P_parallel); 00136 G4ThreeVector dir_parallel=theAdjointPrimary->GetMomentumDirection(); 00137 G4double phi =G4UniformRand()*2.*3.1415926; 00138 G4ThreeVector projectileMomentum = G4ThreeVector(P_perp*std::cos(phi),P_perp*std::sin(phi),P_parallel); 00139 projectileMomentum.rotateUz(dir_parallel); 00140 00141 00142 00143 if (!IsScatProjToProjCase ){ //kill the primary and add a secondary 00144 fParticleChange->ProposeTrackStatus(fStopAndKill); 00145 fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,projectileMomentum)); 00146 //G4cout<<"projectileMomentum "<<projectileMomentum<<G4endl; 00147 } 00148 else { 00149 fParticleChange->ProposeEnergy(projectileKinEnergy); 00150 fParticleChange->ProposeMomentumDirection(projectileMomentum.unit()); 00151 } 00152 00153 00154 00155 00156 }