Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/standard/src/G4LindhardSorensenIonModel.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 //
 28 // GEANT4 Class header file
 29 //
 30 //
 31 // File name:     G4LindhardSorensenIonModel
 32 //
 33 // Author:        Alexander Bagulya & Vladimir Ivanchenko
 34 //
 35 // Creation date: 16.04.2018
 36 //
 37 //
 38 // -------------------------------------------------------------------
 39 //
 40 
 41 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 42 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 43 
 44 #include "G4LindhardSorensenIonModel.hh"
 45 #include "Randomize.hh"
 46 #include "G4PhysicalConstants.hh"
 47 #include "G4SystemOfUnits.hh"
 48 #include "G4Electron.hh"
 49 #include "G4LossTableManager.hh"
 50 #include "G4EmCorrections.hh"
 51 #include "G4ParticleChangeForLoss.hh"
 52 #include "G4Log.hh"
 53 #include "G4DeltaAngle.hh"
 54 #include "G4LindhardSorensenData.hh"
 55 #include "G4BraggModel.hh"
 56 #include "G4BetheBlochModel.hh"
 57 #include "G4IonICRU73Data.hh"
 58 #include "G4AutoLock.hh"
 59 
 60 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 61 
 62 using namespace std;
 63 
 64 G4LindhardSorensenData* G4LindhardSorensenIonModel::lsdata = nullptr;
 65 G4IonICRU73Data* G4LindhardSorensenIonModel::fIonData = nullptr;
 66 
 67 namespace
 68 {
 69   G4Mutex ionXSMutex = G4MUTEX_INITIALIZER;
 70 }
 71 
 72 G4LindhardSorensenIonModel::G4LindhardSorensenIonModel(const G4ParticleDefinition*, 
 73                                                        const G4String& nam)
 74   : G4VEmModel(nam),
 75     twoln10(2.0*G4Log(10.0))
 76 {
 77   theElectron = G4Electron::Electron();
 78   corr = G4LossTableManager::Instance()->EmCorrections();  
 79   nist = G4NistManager::Instance();
 80   fBraggModel = new G4BraggModel();
 81   fBBModel = new G4BetheBlochModel();
 82   fElimit = 2.0*CLHEP::MeV;
 83 }
 84 
 85 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 86 
 87 G4LindhardSorensenIonModel::~G4LindhardSorensenIonModel() {
 88   if(isFirst) {
 89     delete lsdata;
 90     delete fIonData;
 91     lsdata = nullptr;
 92     fIonData = nullptr;
 93   }
 94 }
 95 
 96 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 97 
 98 void G4LindhardSorensenIonModel::Initialise(const G4ParticleDefinition* p,
 99                                             const G4DataVector& ptr)
100 {
101   fBraggModel->Initialise(p, ptr);
102   fBBModel->Initialise(p, ptr);
103   SetParticle(p);
104 
105   // always false before the run
106   SetDeexcitationFlag(false);
107 
108   if(nullptr == fParticleChange) {
109     fParticleChange = GetParticleChangeForLoss();
110     if(UseAngularGeneratorFlag() && nullptr == GetAngularDistribution()) {
111       SetAngularDistribution(new G4DeltaAngle());
112     }
113   }
114   if(nullptr == lsdata) { 
115     G4AutoLock l(&ionXSMutex);
116     if(nullptr == lsdata) { 
117       isFirst = true;
118       lsdata = new G4LindhardSorensenData();
119       fIonData = new G4IonICRU73Data();
120     } 
121     l.unlock();
122   }
123   if(isFirst) {
124     fIonData->Initialise();
125   }
126 }
127 
128 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
129 
130 G4double 
131 G4LindhardSorensenIonModel::GetChargeSquareRatio(const G4ParticleDefinition* p,
132                                                  const G4Material* mat,
133                                                  G4double kinEnergy)
134 {
135   chargeSquare = corr->EffectiveChargeSquareRatio(p, mat, kinEnergy);
136   return chargeSquare;
137 }
138 
139 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
140 
141 G4double 
142 G4LindhardSorensenIonModel::GetParticleCharge(const G4ParticleDefinition* p,
143                                               const G4Material* mat,
144                                               G4double kinEnergy)
145 {
146   return corr->GetParticleCharge(p,mat,kinEnergy);
147 }
148 
149 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
150 
151 void G4LindhardSorensenIonModel::SetupParameters()
152 {
153   mass = particle->GetPDGMass();
154   spin = particle->GetPDGSpin();
155   charge = particle->GetPDGCharge()*inveplus;
156   Zin = G4lrint(std::abs(charge));
157   chargeSquare = charge*charge;
158   eRatio = CLHEP::electron_mass_c2/mass;
159   pRatio = CLHEP::proton_mass_c2/mass;
160   const G4double aMag = 
161     1./(0.5*CLHEP::eplus*CLHEP::hbar_Planck*CLHEP::c_squared);
162   G4double magmom = particle->GetPDGMagneticMoment()*mass*aMag;
163   magMoment2 = magmom*magmom - 1.0;
164   G4double x = 0.8426*CLHEP::GeV;
165   if(spin == 0.0 && mass < GeV) { x = 0.736*CLHEP::GeV; }
166   else if (Zin > 1) { x /= nist->GetA27(Zin); }  
167 
168   formfact = 2.0*CLHEP::electron_mass_c2/(x*x);
169   tlimit = 2.0/formfact;
170 }
171 
172 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
173 
174 G4double G4LindhardSorensenIonModel::MinEnergyCut(const G4ParticleDefinition*,
175                                          const G4MaterialCutsCouple* couple)
176 {
177   return couple->GetMaterial()->GetIonisation()->GetMeanExcitationEnergy();
178 }
179 
180 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
181 
182 G4double 
183 G4LindhardSorensenIonModel::ComputeCrossSectionPerElectron(
184                                                   const G4ParticleDefinition* p,
185                                                   G4double kinEnergy,
186                                                   G4double cutEnergy,
187                                                   G4double maxKinEnergy)        
188 {
189   // take into account formfactor
190   G4double tmax = MaxSecondaryEnergy(p, kinEnergy);
191   G4double emax = std::min(tmax, maxKinEnergy);
192   G4double escaled = kinEnergy*pRatio;
193   G4double cross = (escaled <= fElimit)
194     ? fBraggModel->ComputeCrossSectionPerElectron(p,kinEnergy,cutEnergy,emax)
195     : fBBModel->ComputeCrossSectionPerElectron(p,kinEnergy,cutEnergy,emax);
196   // G4cout << "LS: e= " << kinEnergy << " tmin= " << cutEnergy 
197   //        << " tmax= " << maxEnergy << " cross= " << cross << G4endl; 
198   return cross;
199 }
200 
201 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
202 
203 G4double G4LindhardSorensenIonModel::ComputeCrossSectionPerAtom(
204                                            const G4ParticleDefinition* p,
205                                                  G4double kineticEnergy,
206                                                  G4double Z, G4double,
207                                                  G4double cutEnergy,
208                                                  G4double maxEnergy)
209 {
210   return Z*ComputeCrossSectionPerElectron(p,kineticEnergy,cutEnergy,maxEnergy);
211 }
212 
213 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
214 
215 G4double G4LindhardSorensenIonModel::CrossSectionPerVolume(
216                                            const G4Material* material,
217                                            const G4ParticleDefinition* p,
218                                                  G4double kineticEnergy,
219                                                  G4double cutEnergy,
220                                                  G4double maxEnergy)
221 {
222   return material->GetElectronDensity()
223     *ComputeCrossSectionPerElectron(p,kineticEnergy,cutEnergy,maxEnergy);
224 }
225 
226 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
227 
228 G4double 
229 G4LindhardSorensenIonModel::ComputeDEDXPerVolume(const G4Material* mat,
230                                                  const G4ParticleDefinition* p,
231                                                  G4double kinEnergy,
232                                                  G4double cut)
233 {
234   // formfactor is taken into account in CorrectionsAlongStep(..)
235   G4double tmax      = MaxSecondaryEnergy(p, kinEnergy);
236   G4double cutEnergy = std::min(std::min(cut,tmax), tlimit);
237 
238   G4double escaled = kinEnergy*pRatio;
239   G4double dedx = (escaled <= fElimit)
240     ? fBraggModel->ComputeDEDXPerVolume(mat, p, kinEnergy, cutEnergy)
241     : fBBModel->ComputeDEDXPerVolume(mat, p, kinEnergy, cutEnergy);
242 
243   //G4cout << "E(MeV)=" << kinEnergy/MeV << " dedx=" << dedx 
244   //  << "  " << mat->GetName() << " Ecut(MeV)=" << cutEnergy << G4endl;
245   return dedx;
246 }
247 
248 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
249 
250 void G4LindhardSorensenIonModel::CorrectionsAlongStep(
251                                  const G4MaterialCutsCouple* couple,
252                                  const G4DynamicParticle* dp,
253                                  const G4double& length,
254                                        G4double& eloss)
255 {
256   // no correction at the last step
257   const G4double preKinEnergy = dp->GetKineticEnergy();
258   if(eloss >= preKinEnergy) { return; }
259 
260   const G4ParticleDefinition* p = dp->GetDefinition();
261   SetParticle(p);
262   const G4Material* mat = couple->GetMaterial();
263   const G4double eDensity = mat->GetElectronDensity();
264   const G4double e = std::max(preKinEnergy - eloss*0.5, preKinEnergy*0.5);
265   const G4double tmax = MaxSecondaryEnergy(p, e);
266   const G4double escaled = e*pRatio;
267   const G4double tau = e/mass;
268   const G4double q2 = corr->EffectiveChargeSquareRatio(p, mat, e);
269   const G4int Z = p->GetAtomicNumber();
270 
271   G4double res = 0.0;
272   if(escaled <= fElimit) {
273     // data from ICRU73 or ICRU90
274     if(Z > 2 && Z <= 80) {
275       res = fIonData->GetDEDX(mat, Z, escaled, G4Log(escaled));
276       /*
277   G4cout << "  GetDEDX for Z=" << Z << " in " << mat->GetName()
278   << " Escaled=" << escaled << " E=" 
279   << e << " dEdx=" << res << G4endl;
280       */
281     }
282     if(res > 0.0) {
283       auto pcuts = couple->GetProductionCuts();
284       G4double cut = (nullptr == pcuts) ? tmax : pcuts->GetProductionCut(1);
285       if(cut < tmax) {
286   const G4double x = cut/tmax;
287   res += (G4Log(x)*(tau + 1.)*(tau + 1.)/(tau * (tau + 2.0)) + 1.0 - x) 
288     *q2*CLHEP::twopi_mc2_rcl2*eDensity;
289       }
290       res *= length;
291     } else {
292       // simplified correction
293       res = eloss*q2/chargeSquare;
294     }
295   } else {
296     // Lindhard-Sorensen model
297     const G4double gam = tau + 1.0;
298     const G4double beta2 = tau * (tau+2.0)/(gam*gam);
299     G4double deltaL0 = 2.0*corr->BarkasCorrection(p, mat, e)*(charge-1.)/charge;
300     G4double deltaL = lsdata->GetDeltaL(Zin, gam); 
301 
302     res = eloss + 
303       CLHEP::twopi_mc2_rcl2*q2*eDensity*(deltaL+deltaL0)*length/beta2;
304     /*   
305     G4cout << "  E(GeV)=" << preKinEnergy/GeV << " eloss(MeV)=" << eloss
306      << " L= " << eloss*beta2/(twopi_mc2_rcl2*q2*eDensity*length)
307      << " dL0= " << deltaL0
308      << " dL= " << deltaL << " dE(MeV)=" << res - eloss << G4endl;
309     */
310   }
311   if(res > preKinEnergy || 2*res < eloss) { res = eloss; }
312   /*
313   G4cout << "  G4LindhardSorensenIonModel::CorrectionsAlongStep: E(GeV)=" 
314          << preKinEnergy/GeV << " eloss(MeV)=" << eloss 
315    << " res(MeV)=" << res << G4endl;
316   */
317   eloss = res;
318 }
319 
320 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
321 
322 void G4LindhardSorensenIonModel::SampleSecondaries(
323             std::vector<G4DynamicParticle*>* vdp,
324                                           const G4MaterialCutsCouple* couple,
325                                           const G4DynamicParticle* dp,
326                                           G4double cut,
327                                           G4double maxEnergy)
328 {
329   G4double kineticEnergy = dp->GetKineticEnergy();
330   // take into account formfactor
331   const G4double tmax = MaxSecondaryEnergy(dp->GetDefinition(), kineticEnergy);
332   const G4double minKinEnergy = std::min(cut, tmax);
333   const G4double maxKinEnergy = std::min(maxEnergy, tmax);
334   if(minKinEnergy >= maxKinEnergy) { return; }
335 
336   //G4cout << "G4LindhardSorensenIonModel::SampleSecondaries Emin= " 
337   // << minKinEnergy << " Emax= " << maxKinEnergy << G4endl;
338 
339   G4double totEnergy = kineticEnergy + mass;
340   G4double etot2 = totEnergy*totEnergy;
341   G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*mass)/etot2;
342 
343   G4double deltaKinEnergy, f; 
344   G4double f1 = 0.0;
345   G4double fmax = 1.0;
346   if( 0.0 < spin ) { fmax += 0.5*maxKinEnergy*maxKinEnergy/etot2; }
347 
348   CLHEP::HepRandomEngine* rndmEngineMod = G4Random::getTheEngine();
349   G4double rndm[2];
350 
351   // sampling without nuclear size effect
352   do {
353     rndmEngineMod->flatArray(2, rndm);
354     deltaKinEnergy = minKinEnergy*maxKinEnergy
355                     /(minKinEnergy*(1.0 - rndm[0]) + maxKinEnergy*rndm[0]);
356 
357     f = 1.0 - beta2*deltaKinEnergy/tmax;
358     if( 0.0 < spin ) {
359       f1 = 0.5*deltaKinEnergy*deltaKinEnergy/etot2;
360       f += f1;
361     }
362 
363     // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
364   } while( fmax*rndm[1] > f);
365 
366   // projectile formfactor - suppresion of high energy
367   // delta-electron production at high energy
368   
369   G4double x = formfact*deltaKinEnergy;
370   if(x > 1.e-6) {
371 
372     G4double x1 = 1.0 + x;
373     G4double grej  = 1.0/(x1*x1);
374     if( 0.0 < spin ) {
375       G4double x2 = 0.5*electron_mass_c2*deltaKinEnergy/(mass*mass);
376       grej *= (1.0 + magMoment2*(x2 - f1/f)/(1.0 + x2));
377     }
378     if(grej > 1.1) {
379       G4cout << "### G4LindhardSorensenIonModel WARNING: grej= " << grej
380              << "  " << dp->GetDefinition()->GetParticleName()
381              << " Ekin(MeV)= " <<  kineticEnergy
382              << " delEkin(MeV)= " << deltaKinEnergy
383              << G4endl;
384     }
385     if(rndmEngineMod->flat() > grej) { return; }
386   }
387 
388   G4ThreeVector deltaDirection;
389 
390   if(UseAngularGeneratorFlag()) {
391 
392     const G4Material* mat =  couple->GetMaterial();
393     G4int Z = SelectRandomAtomNumber(mat);
394 
395     deltaDirection = 
396       GetAngularDistribution()->SampleDirection(dp, deltaKinEnergy, Z, mat);
397 
398   } else {
399  
400     G4double deltaMomentum =
401       std::sqrt(deltaKinEnergy * (deltaKinEnergy + 2.0*electron_mass_c2));
402     G4double cost = deltaKinEnergy * (totEnergy + electron_mass_c2) /
403       (deltaMomentum * dp->GetTotalMomentum());
404     cost = std::min(cost, 1.0);
405     G4double sint = std::sqrt((1.0 - cost)*(1.0 + cost));
406 
407     G4double phi = CLHEP::twopi*rndmEngineMod->flat();
408 
409     deltaDirection.set(sint*std::cos(phi),sint*std::sin(phi), cost) ;
410     deltaDirection.rotateUz(dp->GetMomentumDirection());
411   }  
412   /*
413     G4cout << "### G4LindhardSorensenIonModel " 
414            << dp->GetDefinition()->GetParticleName()
415            << " Ekin(MeV)= " <<  kineticEnergy
416            << " delEkin(MeV)= " << deltaKinEnergy
417            << " tmin(MeV)= " << minKinEnergy
418            << " tmax(MeV)= " << maxKinEnergy
419            << " dir= " << dp->GetMomentumDirection()
420            << " dirDelta= " << deltaDirection
421            << G4endl;
422   */
423   // create G4DynamicParticle object for delta ray
424   auto delta = new G4DynamicParticle(theElectron,deltaDirection,deltaKinEnergy);
425 
426   vdp->push_back(delta);
427 
428   // Change kinematics of primary particle
429   kineticEnergy -= deltaKinEnergy;
430   G4ThreeVector finalP = dp->GetMomentum() - delta->GetMomentum();
431   finalP = finalP.unit();
432   
433   fParticleChange->SetProposedKineticEnergy(kineticEnergy);
434   fParticleChange->SetProposedMomentumDirection(finalP);
435 }
436 
437 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
438 
439 G4double 
440 G4LindhardSorensenIonModel::MaxSecondaryEnergy(const G4ParticleDefinition* pd,
441                                                G4double kinEnergy) 
442 {
443   // here particle type is checked for any method
444   SetParticle(pd);
445   G4double tau  = kinEnergy/mass;
446   return 2.0*CLHEP::electron_mass_c2*tau*(tau + 2.) /
447     (1. + 2.0*(tau + 1.)*eRatio + eRatio*eRatio);
448 }
449 
450 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
451