Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/adjoint/src/G4AdjointIonIonisationModel.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  1 //
  2 // ********************************************************************
  3 // * License and Disclaimer                                           *
  4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.                             *
 10 // *                                                                  *
 11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                                                  *
 18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // ********************************************************************
 25 //
 26 
 27 #include "G4AdjointIonIonisationModel.hh"
 28 
 29 #include "G4AdjointCSManager.hh"
 30 #include "G4AdjointElectron.hh"
 31 #include "G4BetheBlochModel.hh"
 32 #include "G4BraggIonModel.hh"
 33 #include "G4GenericIon.hh"
 34 #include "G4NistManager.hh"
 35 #include "G4ParticleChange.hh"
 36 #include "G4PhysicalConstants.hh"
 37 #include "G4SystemOfUnits.hh"
 38 #include "G4TrackStatus.hh"
 39 
 40 ////////////////////////////////////////////////////////////////////////////////
 41 G4AdjointIonIonisationModel::G4AdjointIonIonisationModel()
 42   : G4VEmAdjointModel("Adjoint_IonIonisation")
 43 {
 44   fUseMatrix               = true;
 45   fUseMatrixPerElement     = true;
 46   fApplyCutInRange         = true;
 47   fOneMatrixForAllElements = true;
 48   fSecondPartSameType      = false;
 49 
 50   // The direct EM Model is taken as BetheBloch. It is only used for the
 51   // computation of the differential cross section.
 52   // The Bragg model could be used as an alternative as it offers the same
 53   // differential cross section
 54 
 55   fBetheBlochDirectEMModel  = new G4BetheBlochModel(G4GenericIon::GenericIon());
 56   fBraggIonDirectEMModel    = new G4BraggIonModel(G4GenericIon::GenericIon());
 57   fAdjEquivDirectSecondPart = G4AdjointElectron::AdjointElectron();
 58   fDirectPrimaryPart        = nullptr;
 59 }
 60 
 61 ////////////////////////////////////////////////////////////////////////////////
 62 G4AdjointIonIonisationModel::~G4AdjointIonIonisationModel() {}
 63 
 64 ////////////////////////////////////////////////////////////////////////////////
 65 void G4AdjointIonIonisationModel::SampleSecondaries(
 66   const G4Track& aTrack, G4bool isScatProjToProj,
 67   G4ParticleChange* fParticleChange)
 68 {
 69   const G4DynamicParticle* theAdjointPrimary = aTrack.GetDynamicParticle();
 70 
 71   // Elastic inverse scattering
 72   G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
 73   G4double adjointPrimP         = theAdjointPrimary->GetTotalMomentum();
 74 
 75   if(adjointPrimKinEnergy > GetHighEnergyLimit() * 0.999)
 76   {
 77     return;
 78   }
 79 
 80   // Sample secondary energy
 81   G4double projectileKinEnergy =
 82     SampleAdjSecEnergyFromCSMatrix(adjointPrimKinEnergy, isScatProjToProj);
 83   // Caution !!!this weight correction should be always applied
 84   CorrectPostStepWeight(fParticleChange, aTrack.GetWeight(),
 85                         adjointPrimKinEnergy, projectileKinEnergy,
 86                         isScatProjToProj);
 87 
 88   // Kinematics:
 89   // we consider a two body elastic scattering for the forward processes where
 90   // the projectile knock on an e- at rest and gives it part of its  energy
 91   G4double projectileM0          = fAdjEquivDirectPrimPart->GetPDGMass();
 92   G4double projectileTotalEnergy = projectileM0 + projectileKinEnergy;
 93   G4double projectileP2 =
 94     projectileTotalEnergy * projectileTotalEnergy - projectileM0 * projectileM0;
 95 
 96   // Companion
 97   G4double companionM0 = fAdjEquivDirectPrimPart->GetPDGMass();
 98   if(isScatProjToProj)
 99   {
100     companionM0 = fAdjEquivDirectSecondPart->GetPDGMass();
101   }
102   G4double companionTotalEnergy =
103     companionM0 + projectileKinEnergy - adjointPrimKinEnergy;
104   G4double companionP2 =
105     companionTotalEnergy * companionTotalEnergy - companionM0 * companionM0;
106 
107   // Projectile momentum
108   G4double P_parallel =
109     (adjointPrimP * adjointPrimP + projectileP2 - companionP2) /
110     (2. * adjointPrimP);
111   G4double P_perp = std::sqrt(projectileP2 - P_parallel * P_parallel);
112   G4ThreeVector dir_parallel = theAdjointPrimary->GetMomentumDirection();
113   G4double phi               = G4UniformRand() * twopi;
114   G4ThreeVector projectileMomentum =
115     G4ThreeVector(P_perp * std::cos(phi), P_perp * std::sin(phi), P_parallel);
116   projectileMomentum.rotateUz(dir_parallel);
117 
118   if(!isScatProjToProj)
119   {  // kill the primary and add a secondary
120     fParticleChange->ProposeTrackStatus(fStopAndKill);
121     fParticleChange->AddSecondary(
122       new G4DynamicParticle(fAdjEquivDirectPrimPart, projectileMomentum));
123   }
124   else
125   {
126     fParticleChange->ProposeEnergy(projectileKinEnergy);
127     fParticleChange->ProposeMomentumDirection(projectileMomentum.unit());
128   }
129 }
130 
131 ////////////////////////////////////////////////////////////////////////////////
132 G4double G4AdjointIonIonisationModel::DiffCrossSectionPerAtomPrimToSecond(
133   G4double kinEnergyProj, G4double kinEnergyProd, G4double Z, G4double A)
134 {
135   // Probably that here the Bragg Model should be also used for
136   // kinEnergyProj/nuc<2MeV
137   G4double dSigmadEprod = 0.;
138   G4double Emax_proj    = GetSecondAdjEnergyMaxForProdToProj(kinEnergyProd);
139   G4double Emin_proj    = GetSecondAdjEnergyMinForProdToProj(kinEnergyProd);
140 
141   G4double kinEnergyProjScaled = fMassRatio * kinEnergyProj;
142 
143   // the produced particle should have a kinetic energy smaller than the
144   // projectile
145   if(kinEnergyProj > Emin_proj && kinEnergyProj <= Emax_proj)
146   {
147     G4double Tmax = kinEnergyProj;
148 
149     G4double E1 = kinEnergyProd;
150     G4double E2 = kinEnergyProd * 1.000001;
151     G4double dE = (E2 - E1);
152     G4double sigma1, sigma2;
153     fDirectModel = fBraggIonDirectEMModel;
154     if(kinEnergyProjScaled > 2. * MeV && !fUseOnlyBragg)
155       fDirectModel = fBetheBlochDirectEMModel;
156     sigma1 = fDirectModel->ComputeCrossSectionPerAtom(
157       fDirectPrimaryPart, kinEnergyProj, Z, A, E1, 1.e20);
158     sigma2 = fDirectModel->ComputeCrossSectionPerAtom(
159       fDirectPrimaryPart, kinEnergyProj, Z, A, E2, 1.e20);
160 
161     dSigmadEprod = (sigma1 - sigma2) / dE;
162 
163     if(dSigmadEprod > 1.)
164     {
165       G4cout << "sigma1 " << kinEnergyProj / MeV << '\t' << kinEnergyProd / MeV
166              << '\t' << sigma1 << G4endl;
167       G4cout << "sigma2 " << kinEnergyProj / MeV << '\t' << kinEnergyProd / MeV
168              << '\t' << sigma2 << G4endl;
169       G4cout << "dsigma " << kinEnergyProj / MeV << '\t' << kinEnergyProd / MeV
170              << '\t' << dSigmadEprod << G4endl;
171     }
172 
173     if(fDirectModel == fBetheBlochDirectEMModel)
174     {
175       // correction of differential cross section at high energy to correct for
176       // the suppression of particle at secondary at high energy used in the
177       // Bethe Bloch Model. This correction consist to multiply by g the
178       // probability function used to test the rejection of a secondary Source
179       // code taken from G4BetheBlochModel::SampleSecondaries
180 
181       G4double deltaKinEnergy = kinEnergyProd;
182 
183       G4double x = fFormFact * deltaKinEnergy;
184       if(x > 1.e-6)
185       {
186         G4double totEnergy = kinEnergyProj + fMass;
187         G4double etot2     = totEnergy * totEnergy;
188         G4double beta2 = kinEnergyProj * (kinEnergyProj + 2.0 * fMass) / etot2;
189         G4double f1    = 0.0;
190         G4double f     = 1.0 - beta2 * deltaKinEnergy / Tmax;
191         if(0.5 == fSpin)
192         {
193           f1 = 0.5 * deltaKinEnergy * deltaKinEnergy / etot2;
194           f += f1;
195         }
196         G4double x1 = 1.0 + x;
197         G4double gg = 1.0 / (x1 * x1);
198         if(0.5 == fSpin)
199         {
200           G4double x2 =
201             0.5 * electron_mass_c2 * deltaKinEnergy / (fMass * fMass);
202           gg *= (1.0 + fMagMoment2 * (x2 - f1 / f) / (1.0 + x2));
203         }
204         if(gg > 1.0)
205         {
206           G4cout << "### G4BetheBlochModel in Adjoint Sim WARNING: gg= " << gg
207                  << G4endl;
208           gg = 1.;
209         }
210         dSigmadEprod *= gg;
211       }
212     }
213   }
214   return dSigmadEprod;
215 }
216 
217 ////////////////////////////////////////////////////////////////////////////////
218 void G4AdjointIonIonisationModel::SetIon(G4ParticleDefinition* adj_ion,
219                                          G4ParticleDefinition* fwd_ion)
220 {
221   fDirectPrimaryPart      = fwd_ion;
222   fAdjEquivDirectPrimPart = adj_ion;
223 
224   DefineProjectileProperty();
225 }
226 
227 ////////////////////////////////////////////////////////////////////////////////
228 void G4AdjointIonIonisationModel::CorrectPostStepWeight(
229   G4ParticleChange* fParticleChange, G4double old_weight,
230   G4double adjointPrimKinEnergy, G4double projectileKinEnergy, G4bool)
231 {
232   // It is needed because the direct cross section used to compute the
233   // differential cross section is not the one used in
234   // the direct model where the GenericIon stuff is considered with correction
235   // of effective charge.  In the direct model the samnepl of secondaries does
236   // not reflect the integral cross section. The integral fwd cross section that
237   // we used to compute the differential CS match the sample of secondaries in
238   // the forward case despite the fact that its is not the same total CS than in
239   // the FWD case. For this reason an extra weight correction is needed at the
240   // end.
241 
242   G4double new_weight = old_weight;
243 
244   // the correction of CS due to the problem explained above
245   G4double kinEnergyProjScaled = fMassRatio * projectileKinEnergy;
246   fDirectModel                 = fBraggIonDirectEMModel;
247   if(kinEnergyProjScaled > 2. * MeV && !fUseOnlyBragg)
248     fDirectModel = fBetheBlochDirectEMModel;
249   G4double UsedFwdCS = fDirectModel->ComputeCrossSectionPerAtom(
250     fDirectPrimaryPart, projectileKinEnergy, 1, 1, fTcutSecond, 1.e20);
251   G4double chargeSqRatio = 1.;
252   if(fChargeSquare > 1.)
253     chargeSqRatio = fDirectModel->GetChargeSquareRatio(
254       fDirectPrimaryPart, fCurrentMaterial, projectileKinEnergy);
255   G4double CorrectFwdCS =
256     chargeSqRatio * fDirectModel->ComputeCrossSectionPerAtom(
257                       G4GenericIon::GenericIon(), kinEnergyProjScaled, 1, 1,
258                       fTcutSecond, 1.e20);
259   // May be some check is needed if UsedFwdCS ==0 probably that then we should
260   // avoid a secondary to be produced,
261   if(UsedFwdCS > 0.)
262     new_weight *= CorrectFwdCS / UsedFwdCS;
263 
264   // additional CS correction needed for cross section biasing in general.
265   // May be wrong for ions. Most of the time not used.
266   new_weight *=
267     G4AdjointCSManager::GetAdjointCSManager()->GetPostStepWeightCorrection() /
268     fCsBiasingFactor;
269 
270   new_weight *= projectileKinEnergy / adjointPrimKinEnergy;
271 
272   fParticleChange->SetParentWeightByProcess(false);
273   fParticleChange->SetSecondaryWeightByProcess(false);
274   fParticleChange->ProposeParentWeight(new_weight);
275 }
276 
277 ////////////////////////////////////////////////////////////////////////////////
278 void G4AdjointIonIonisationModel::DefineProjectileProperty()
279 {
280   // Slightly modified code taken from G4BetheBlochModel::SetParticle
281   G4String pname = fDirectPrimaryPart->GetParticleName();
282 
283   fMass           = fDirectPrimaryPart->GetPDGMass();
284   fMassRatio      = G4GenericIon::GenericIon()->GetPDGMass() / fMass;
285   fSpin           = fDirectPrimaryPart->GetPDGSpin();
286   G4double q      = fDirectPrimaryPart->GetPDGCharge() / eplus;
287   fChargeSquare   = q * q;
288   fRatio          = electron_mass_c2 / fMass;
289   fOnePlusRatio2  = (1. + fRatio) * (1. + fRatio);
290   fOneMinusRatio2 = (1. - fRatio) * (1. - fRatio);
291   G4double magmom = fDirectPrimaryPart->GetPDGMagneticMoment() * fMass /
292                     (0.5 * eplus * hbar_Planck * c_squared);
293   fMagMoment2 = magmom * magmom - 1.0;
294   if(fDirectPrimaryPart->GetLeptonNumber() == 0)
295   {
296     G4double x = 0.8426 * GeV;
297     if(fSpin == 0.0 && fMass < GeV)
298     {
299       x = 0.736 * GeV;
300     }
301     else if(fMass > GeV)
302     {
303       x /= G4NistManager::Instance()->GetZ13(fMass / proton_mass_c2);
304     }
305     fFormFact = 2.0 * electron_mass_c2 / (x * x);
306   }
307 }
308 
309 //////////////////////////////////////////////////////////////////////////////
310 G4double G4AdjointIonIonisationModel::GetSecondAdjEnergyMaxForScatProjToProj(
311   G4double primAdjEnergy)
312 {
313   return primAdjEnergy * fOnePlusRatio2 /
314          (fOneMinusRatio2 - 2. * fRatio * primAdjEnergy / fMass);
315 }
316 
317 //////////////////////////////////////////////////////////////////////////////
318 G4double G4AdjointIonIonisationModel::GetSecondAdjEnergyMinForScatProjToProj(
319   G4double primAdjEnergy, G4double tcut)
320 {
321   return primAdjEnergy + tcut;
322 }
323 
324 //////////////////////////////////////////////////////////////////////////////
325 G4double G4AdjointIonIonisationModel::GetSecondAdjEnergyMaxForProdToProj(
326   G4double)
327 {
328   return GetHighEnergyLimit();
329 }
330 
331 //////////////////////////////////////////////////////////////////////////////
332 G4double G4AdjointIonIonisationModel::GetSecondAdjEnergyMinForProdToProj(
333   G4double primAdjEnergy)
334 {
335   return (2. * primAdjEnergy - 4. * fMass +
336           std::sqrt(4. * primAdjEnergy * primAdjEnergy + 16. * fMass * fMass +
337                     8. * primAdjEnergy * fMass * (1. / fRatio + fRatio))) /
338          4.;
339 }
340