64 theBraggDirectEMModel =
new G4BraggModel(projectileDefinition);
74 DefineProjectileProperty();
85 G4bool IsScatProjToProjCase,
106 adjointPrimKinEnergy,
108 IsScatProjToProjCase);
117 G4double projectileTotalEnergy = projectileM0+projectileKinEnergy;
118 G4double projectileP2 = projectileTotalEnergy*projectileTotalEnergy - projectileM0*projectileM0;
125 if (IsScatProjToProjCase) {
128 G4double companionTotalEnergy =companionM0+ projectileKinEnergy-adjointPrimKinEnergy;
129 G4double companionP2 = companionTotalEnergy*companionTotalEnergy - companionM0*companionM0;
134 G4double P_parallel = (adjointPrimP*adjointPrimP + projectileP2 - companionP2)/(2.*adjointPrimP);
135 G4double P_perp = std::sqrt( projectileP2 - P_parallel*P_parallel);
139 projectileMomentum.
rotateUz(dir_parallel);
143 if (!IsScatProjToProjCase ){
161 G4bool IsScatProjToProjCase,
179 if (!IsScatProjToProjCase){
181 eEnergy=adjointPrimKinEnergy;
184 if (Emin>=Emax)
return;
187 newCS=newCS*(b-
a)/eEnergy;
195 if (Emin>=Emax)
return;
196 G4double diff1=Emin-adjointPrimKinEnergy;
197 G4double diff2=Emax-adjointPrimKinEnergy;
199 G4double t1=adjointPrimKinEnergy*(1./diff1-1./diff2);
200 G4double t2=adjointPrimKinEnergy*(1./Emin-1./Emax);
205 newCS=newCS*sum_t/adjointPrimKinEnergy/adjointPrimKinEnergy;
209 projectileKinEnergy =adjointPrimKinEnergy +1./(1./diff1-q);
214 projectileKinEnergy =1./(1./Emin-q);
217 projectileKinEnergy=Emin*std::pow(Emax/Emin,
G4UniformRand());
220 eEnergy=projectileKinEnergy-adjointPrimKinEnergy;
227 G4double diffCS_perAtom_Used=
twopi_mc2_rcl2*mass*adjointPrimKinEnergy/projectileKinEnergy/projectileKinEnergy/eEnergy/eEnergy;
243 w_corr*=diffCS/diffCS_perAtom_Used;
260 G4double projectileTotalEnergy = projectileM0+projectileKinEnergy;
261 G4double projectileP2 = projectileTotalEnergy*projectileTotalEnergy - projectileM0*projectileM0;
268 if (IsScatProjToProjCase) {
271 G4double companionTotalEnergy =companionM0+ projectileKinEnergy-adjointPrimKinEnergy;
272 G4double companionP2 = companionTotalEnergy*companionTotalEnergy - companionM0*companionM0;
277 G4double P_parallel = (adjointPrimP*adjointPrimP + projectileP2 - companionP2)/(2.*adjointPrimP);
278 G4double P_perp = std::sqrt( projectileP2 - P_parallel*P_parallel);
282 projectileMomentum.
rotateUz(dir_parallel);
286 if (!IsScatProjToProjCase ){
320 if (kinEnergyProj>Emin_proj && kinEnergyProj<=Emax_proj){
327 if (kinEnergyProj >2.*
MeV){
337 dSigmadEprod=(sigma1-sigma2)/dE;
338 if (dSigmadEprod>1.) {
341 G4cout<<
"dsigma "<<kinEnergyProj/
MeV<<
'\t'<<kinEnergyProd/
MeV<<
'\t'<<dSigmadEprod<<
G4endl;
354 G4double deltaKinEnergy = kinEnergyProd;
367 G4double totEnergy = kinEnergyProj + mass;
368 G4double etot2 = totEnergy*totEnergy;
369 G4double beta2 = kinEnergyProj*(kinEnergyProj + 2.0*mass)/etot2;
372 f = 1.0 - beta2*deltaKinEnergy/Tmax;
374 f1 = 0.5*deltaKinEnergy*deltaKinEnergy/etot2;
381 gg *= (1.0 + magMoment2*(x2 - f1/
f)/(1.0 + x2));
384 G4cout <<
"### G4BetheBlochModel in Adjoint Sim WARNING: g= " <<
g
401 void G4AdjointhIonisationModel::DefineProjectileProperty()
407 pname !=
"deuteron" && pname !=
"triton") {
417 ratio2 = ratio*ratio;
418 one_plus_ratio_2=(1+ratio)*(1+ratio);
419 one_minus_ratio_2=(1-ratio)*(1-ratio);
422 magMoment2 = magmom*magmom - 1.0;
426 if(spin == 0.0 && mass <
GeV) {x = 0.736*
GeV;}
427 else if(mass >
GeV) {
432 tlimit = 2.0/formfact;
440 G4bool IsScatProjToProjCase)
448 if (!IsScatProjToProjCase ){
452 Cross*=(1./Emin_proj -1./Emax_proj)/primEnergy;
465 G4double diff1=Emin_proj-primEnergy;
466 G4double diff2=Emax_proj-primEnergy;
467 G4double t1=(1./diff1+1./Emin_proj-1./diff2-1./Emax_proj)/primEnergy;
469 G4double t2=2.*std::log(Emax_proj/Emin_proj)/primEnergy/primEnergy;
481 G4double Tmax=PrimAdjEnergy*one_plus_ratio_2/(one_minus_ratio_2-2.*ratio*PrimAdjEnergy/mass);
487 {
return PrimAdjEnergy+Tcut;
497 {
G4double Tmin= (2*PrimAdjEnergy-4*mass + std::sqrt(4.*PrimAdjEnergy*PrimAdjEnergy +16.*mass*mass + 8.*PrimAdjEnergy*mass*(1/ratio +ratio)))/4.;
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
G4double GetPostStepWeightCorrection()
const G4DynamicParticle * GetDynamicParticle() const
G4double SampleAdjSecEnergyFromCSMatrix(size_t MatrixIndex, G4double prim_energy, G4bool IsScatProjToProjCase)
G4ParticleDefinition * theDirectPrimaryPartDef
void ProposeParentWeight(G4double finalWeight)
G4double GetZ13(G4double Z)
const G4MaterialCutsCouple * GetMaterialCutsCouple() const
virtual G4double GetSecondAdjEnergyMaxForScatProjToProjCase(G4double PrimAdjEnergy)
void RapidSampleSecondaries(const G4Track &aTrack, G4bool IsScatProjToProjCase, G4ParticleChange *fParticleChange)
static G4AdjointElectron * AdjointElectron()
G4double CS_biasing_factor
G4Material * currentMaterial
G4ParticleDefinition * theAdjEquivOfDirectPrimPartDef
static G4NistManager * Instance()
const G4String & GetParticleName() const
virtual G4double GetSecondAdjEnergyMinForScatProjToProjCase(G4double PrimAdjEnergy, G4double Tcut=0)
virtual ~G4AdjointhIonisationModel()
G4double GetTotalMomentum() const
void SetSecondaryWeightByProcess(G4bool)
void SetParentWeightByProcess(G4bool)
G4VEmModel * theDirectEMModel
G4GLOB_DLL std::ostream G4cout
G4double GetElectronDensity() const
const G4ThreeVector & GetMomentumDirection() const
Hep3Vector & rotateUz(const Hep3Vector &)
static G4Proton * Proton()
const G4String & GetParticleType() const
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0., G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
void DefineCurrentMaterial(const G4MaterialCutsCouple *couple)
G4bool UseOnlyOneMatrixForAllElements
virtual G4double AdjointCrossSection(const G4MaterialCutsCouple *aCouple, G4double primEnergy, G4bool IsScatProjToProjCase)
G4double GetPDGMass() const
G4double mass_ratio_projectile
static G4AdjointProton * AdjointProton()
virtual void SampleSecondaries(const G4Track &aTrack, G4bool IsScatProjToProjCase, G4ParticleChange *fParticleChange)
void ProposeEnergy(G4double finalEnergy)
void AddSecondary(G4Track *aSecondary)
virtual G4double AdjointCrossSection(const G4MaterialCutsCouple *aCouple, G4double primEnergy, G4bool IsScatProjToProjCase)
G4double GetWeight() const
G4double GetPDGSpin() const
G4ParticleDefinition * theAdjEquivOfDirectSecondPartDef
G4double currentTcutForDirectSecond
G4bool UseMatrixPerElement
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
G4double GetPDGMagneticMoment() const
void ProposeTrackStatus(G4TrackStatus status)
virtual G4double GetSecondAdjEnergyMinForProdToProjCase(G4double PrimAdjEnergy)
G4double GetPDGCharge() const
static G4AdjointCSManager * GetAdjointCSManager()
G4int GetLeptonNumber() const
virtual void CorrectPostStepWeight(G4ParticleChange *fParticleChange, G4double old_weight, G4double adjointPrimKinEnergy, G4double projectileKinEnergy, G4bool IsScatProjToProjCase)
virtual G4double DiffCrossSectionPerAtomPrimToSecond(G4double kinEnergyProj, G4double kinEnergyProd, G4double Z, G4double A=0.)
G4AdjointhIonisationModel(G4ParticleDefinition *projectileDefinition)
virtual G4double GetSecondAdjEnergyMaxForProdToProjCase(G4double PrimAdjEnergy)
G4bool second_part_of_same_type