source: trunk/source/processes/electromagnetic/pii/src/G4hImpactIonisation.cc @ 1350

Last change on this file since 1350 was 1350, checked in by garnier, 13 years ago

update to last version 4.9.4
File size: 53.6 KB
Line 
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// ------------------------------------------------------------
29// G4RDHadronIonisation
30//
31// $Id: G4hImpactIonisation.cc,v 1.4 2010/11/25 19:49:43 pia Exp $
32// GEANT4 tag $Name: geant4-09-04-ref-00 $
33//
34// Author: Maria Grazia Pia (MariaGrazia.Pia@ge.infn.it)
35//
36// 08 Sep 2008 - MGP - Created (initially based on G4hLowEnergyIonisation)
37//                     Added PIXE capabilities
38//                     Partial clean-up of the implementation (more needed)
39//                     Calculation of MicroscopicCrossSection delegated to specialised cla// Documentation available in:
40// M.G. Pia et al., PIXE Simulation With Geant4,
41// IEEE Trans. Nucl. Sci., vol. 56, no. 6, pp. 3614-3649, Dec. 2009.
42
43//
44// ------------------------------------------------------------
45 
46#include "G4hImpactIonisation.hh"
47#include "globals.hh"
48#include "G4ios.hh"
49#include "Randomize.hh"
50#include "G4Poisson.hh"
51#include "G4UnitsTable.hh"
52#include "G4EnergyLossTables.hh"
53#include "G4Material.hh"
54#include "G4DynamicParticle.hh"
55#include "G4ParticleDefinition.hh"
56#include "G4AtomicDeexcitation.hh"
57#include "G4AtomicTransitionManager.hh"
58#include "G4PixeCrossSectionHandler.hh"
59#include "G4IInterpolator.hh"
60#include "G4LogLogInterpolator.hh"
61#include "G4Gamma.hh"
62#include "G4Electron.hh"
63#include "G4Proton.hh"               
64#include "G4ProcessManager.hh"
65#include "G4ProductionCutsTable.hh"
66#include "G4PhysicsLogVector.hh"       
67#include "G4PhysicsLinearVector.hh"   
68
69#include "G4VLowEnergyModel.hh"
70#include "G4hNuclearStoppingModel.hh" 
71#include "G4hBetheBlochModel.hh"       
72#include "G4hParametrisedLossModel.hh"
73#include "G4QAOLowEnergyLoss.hh"       
74#include "G4hIonEffChargeSquare.hh"   
75#include "G4IonChuFluctuationModel.hh"
76#include "G4IonYangFluctuationModel.hh"
77
78#include "G4MaterialCutsCouple.hh"
79#include "G4Track.hh"
80#include "G4Step.hh"
81
82G4hImpactIonisation::G4hImpactIonisation(const G4String& processName)
83  : G4hRDEnergyLoss(processName),
84    betheBlochModel(0),
85    protonModel(0),
86    antiprotonModel(0),
87    theIonEffChargeModel(0),
88    theNuclearStoppingModel(0),
89    theIonChuFluctuationModel(0),
90    theIonYangFluctuationModel(0),
91    protonTable("ICRU_R49p"),
92    antiprotonTable("ICRU_R49p"),
93    theNuclearTable("ICRU_R49"),
94    nStopping(true),
95    theBarkas(true),
96    theMeanFreePathTable(0),
97    paramStepLimit (0.005),
98    pixeCrossSectionHandler(0)
99{ 
100  InitializeMe();
101}
102
103
104
105void G4hImpactIonisation::InitializeMe()
106{
107  LowestKineticEnergy  = 10.0*eV ;
108  HighestKineticEnergy = 100.0*GeV ;
109  MinKineticEnergy     = 10.0*eV ; 
110  TotBin               = 360 ;
111  protonLowEnergy      = 1.*keV ;
112  protonHighEnergy     = 100.*MeV ;
113  antiprotonLowEnergy  = 25.*keV ;
114  antiprotonHighEnergy = 2.*MeV ;
115  minGammaEnergy       = 100 * eV;
116  minElectronEnergy    = 250.* eV;
117  verboseLevel         = 0;
118 
119  // Min and max energy of incident particle for the calculation of shell cross sections
120  // for PIXE generation
121  eMinPixe = 1.* keV;
122  eMaxPixe = 200. * MeV;
123 
124  G4String defaultPixeModel("ecpssr"); 
125  modelK = defaultPixeModel;
126  modelL = defaultPixeModel;
127  modelM = defaultPixeModel;
128}
129
130
131G4hImpactIonisation::~G4hImpactIonisation()
132{
133  if (theMeanFreePathTable) 
134    {
135      theMeanFreePathTable->clearAndDestroy();
136      delete theMeanFreePathTable;
137    }
138
139  if (betheBlochModel) delete betheBlochModel;
140  if (protonModel) delete protonModel;
141  if (antiprotonModel) delete antiprotonModel;
142  if (theNuclearStoppingModel) delete theNuclearStoppingModel;
143  if (theIonEffChargeModel) delete theIonEffChargeModel;
144  if (theIonChuFluctuationModel) delete theIonChuFluctuationModel;
145  if (theIonYangFluctuationModel) delete theIonYangFluctuationModel;
146
147  delete pixeCrossSectionHandler;
148
149  // ---- MGP ---- The following is to be checked
150  //  if (shellVacancy) delete shellVacancy;
151
152  cutForDelta.clear();
153}
154
155// --------------------------------------------------------------------
156void G4hImpactIonisation::SetElectronicStoppingPowerModel(const G4ParticleDefinition* particle,
157                                                          const G4String& dedxTable)
158// This method defines the ionisation parametrisation method via its name
159{
160  if (particle->GetPDGCharge() > 0 ) 
161    {
162      // Positive charge
163      SetProtonElectronicStoppingPowerModel(dedxTable) ;
164    } 
165  else 
166    {
167      // Antiprotons
168      SetAntiProtonElectronicStoppingPowerModel(dedxTable) ;
169    }
170}
171
172
173// --------------------------------------------------------------------
174void G4hImpactIonisation::InitializeParametrisation() 
175
176{
177  // Define models for parametrisation of electronic energy losses
178  betheBlochModel = new G4hBetheBlochModel("Bethe-Bloch") ;
179  protonModel = new G4hParametrisedLossModel(protonTable) ;
180  protonHighEnergy = std::min(protonHighEnergy,protonModel->HighEnergyLimit(0, 0));
181  antiprotonModel = new G4QAOLowEnergyLoss(antiprotonTable) ;
182  theNuclearStoppingModel = new G4hNuclearStoppingModel(theNuclearTable) ;
183  theIonEffChargeModel = new G4hIonEffChargeSquare("Ziegler1988") ;
184  theIonChuFluctuationModel = new G4IonChuFluctuationModel("Chu") ;
185  theIonYangFluctuationModel = new G4IonYangFluctuationModel("Yang") ;
186}
187
188
189// --------------------------------------------------------------------
190void G4hImpactIonisation::BuildPhysicsTable(const G4ParticleDefinition& particleDef)
191
192//  just call BuildLossTable+BuildLambdaTable
193{
194
195  // Verbose print-out
196  if(verboseLevel > 0) 
197    {
198      G4cout << "G4hImpactIonisation::BuildPhysicsTable for "
199             << particleDef.GetParticleName()
200             << " mass(MeV)= " << particleDef.GetPDGMass()/MeV
201             << " charge= " << particleDef.GetPDGCharge()/eplus
202             << " type= " << particleDef.GetParticleType()
203             << G4endl;
204     
205      if(verboseLevel > 1) 
206        {
207          G4ProcessVector* pv = particleDef.GetProcessManager()->GetProcessList();
208         
209          G4cout << " 0: " << (*pv)[0]->GetProcessName() << " " << (*pv)[0]
210                 << " 1: " << (*pv)[1]->GetProcessName() << " " << (*pv)[1]
211            //        << " 2: " << (*pv)[2]->GetProcessName() << " " << (*pv)[2]
212                 << G4endl;
213          G4cout << "ionModel= " << theIonEffChargeModel
214                 << " MFPtable= " << theMeanFreePathTable
215                 << " iniMass= " << initialMass
216                 << G4endl;
217        }
218    }
219  // End of verbose print-out
220
221  if (particleDef.GetParticleType() == "nucleus" &&
222      particleDef.GetParticleName() != "GenericIon" &&
223      particleDef.GetParticleSubType() == "generic")
224    {
225
226      G4EnergyLossTables::Register(&particleDef,
227                                   theDEDXpTable,
228                                   theRangepTable,
229                                   theInverseRangepTable,
230                                   theLabTimepTable,
231                                   theProperTimepTable,
232                                   LowestKineticEnergy, HighestKineticEnergy,
233                                   proton_mass_c2/particleDef.GetPDGMass(),
234                                   TotBin);
235
236      return;
237    }
238
239  if( !CutsWhereModified() && theLossTable) return;
240
241  InitializeParametrisation() ;
242  G4Proton* proton = G4Proton::Proton();
243  G4AntiProton* antiproton = G4AntiProton::AntiProton();
244
245  charge = particleDef.GetPDGCharge() / eplus;
246  chargeSquare = charge*charge ;
247
248  const G4ProductionCutsTable* theCoupleTable=
249    G4ProductionCutsTable::GetProductionCutsTable();
250  size_t numOfCouples = theCoupleTable->GetTableSize();
251
252  cutForDelta.clear();
253  cutForGamma.clear();
254
255  for (size_t j=0; j<numOfCouples; j++) {
256
257    // get material parameters needed for the energy loss calculation
258    const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(j);
259    const G4Material* material= couple->GetMaterial();
260
261    // the cut cannot be below lowest limit
262    G4double tCut = (*(theCoupleTable->GetEnergyCutsVector(1)))[j];
263    if(tCut > HighestKineticEnergy) tCut = HighestKineticEnergy;
264
265    G4double excEnergy = material->GetIonisation()->GetMeanExcitationEnergy();
266
267    tCut = std::max(tCut,excEnergy);
268    cutForDelta.push_back(tCut);
269
270    // the cut cannot be below lowest limit
271    tCut = (*(theCoupleTable->GetEnergyCutsVector(0)))[j];
272    if(tCut > HighestKineticEnergy) tCut = HighestKineticEnergy;
273    tCut = std::max(tCut,minGammaEnergy);
274    cutForGamma.push_back(tCut);
275  }
276
277  if(verboseLevel > 0) {
278    G4cout << "Cuts are defined " << G4endl;
279  }
280
281  if(0.0 < charge)
282    {
283      {
284        BuildLossTable(*proton) ;
285
286        //      The following vector has a fixed dimension (see src/G4hImpactLoss.cc for more details)       
287        //      It happended in the past that caused memory corruption errors. The problem is still pending, even if temporary solved
288        //        G4cout << "[NOTE]: __LINE__=" << __LINE__ << ", particleDef=" << particleDef.GetParticleName() << ", proton=" << proton << ", theLossTable=" << theLossTable << ", CounterOfpProcess=" << CounterOfpProcess << G4endl;
289       
290        RecorderOfpProcess[CounterOfpProcess] = theLossTable ;
291        CounterOfpProcess++;
292      }
293    } else {
294    {
295      BuildLossTable(*antiproton) ;
296       
297      //      The following vector has a fixed dimension (see src/G4hImpactLoss.cc for more details)       
298      //      It happended in the past that caused memory corruption errors. The problem is still pending, even if temporary solved
299      //        G4cout << "[NOTE]: __LINE__=" << __LINE__ << ", particleDef=" << particleDef.GetParticleName() << ", antiproton=" << antiproton << ", theLossTable=" << theLossTable << ", CounterOfpbarProcess=" << CounterOfpbarProcess << G4endl;
300       
301      RecorderOfpbarProcess[CounterOfpbarProcess] = theLossTable ;
302      CounterOfpbarProcess++;
303    }
304  }
305
306  if(verboseLevel > 0) {
307    G4cout << "G4hImpactIonisation::BuildPhysicsTable: "
308           << "Loss table is built "
309      //           << theLossTable
310           << G4endl;
311  }
312
313  BuildLambdaTable(particleDef) ;
314  //  BuildDataForFluorescence(particleDef);
315
316  if(verboseLevel > 1) {
317    G4cout << (*theMeanFreePathTable) << G4endl;
318  }
319
320  if(verboseLevel > 0) {
321    G4cout << "G4hImpactIonisation::BuildPhysicsTable: "
322           << "DEDX table will be built "
323      //           << theDEDXpTable << " " << theDEDXpbarTable
324      //           << " " << theRangepTable << " " << theRangepbarTable
325           << G4endl;
326  }
327
328  BuildDEDXTable(particleDef) ;
329
330  if(verboseLevel > 1) {
331    G4cout << (*theDEDXpTable) << G4endl;
332  }
333
334  if((&particleDef == proton) ||  (&particleDef == antiproton)) PrintInfoDefinition() ;
335
336  if(verboseLevel > 0) {
337    G4cout << "G4hImpactIonisation::BuildPhysicsTable: end for "
338           << particleDef.GetParticleName() << G4endl;
339  }
340}
341
342
343
344
345
346// --------------------------------------------------------------------
347void G4hImpactIonisation::BuildLossTable(const G4ParticleDefinition& particleDef)
348{
349
350  // Initialisation
351  G4double lowEdgeEnergy , ionloss, ionlossBB, paramB ;
352  G4double lowEnergy, highEnergy;
353  G4Proton* proton = G4Proton::Proton();
354
355  if (particleDef == *proton) 
356    {
357      lowEnergy = protonLowEnergy ;
358      highEnergy = protonHighEnergy ;
359      charge = 1. ;
360    } 
361  else 
362    {
363      lowEnergy = antiprotonLowEnergy ;
364      highEnergy = antiprotonHighEnergy ;
365      charge = -1. ;
366    }
367  chargeSquare = 1. ;
368
369  const G4ProductionCutsTable* theCoupleTable=
370    G4ProductionCutsTable::GetProductionCutsTable();
371  size_t numOfCouples = theCoupleTable->GetTableSize();
372
373  if ( theLossTable) 
374    {
375      theLossTable->clearAndDestroy();
376      delete theLossTable;
377    }
378
379  theLossTable = new G4PhysicsTable(numOfCouples);
380
381  //  loop for materials
382  for (size_t j=0; j<numOfCouples; j++) {
383
384    // create physics vector and fill it
385    G4PhysicsLogVector* aVector = new G4PhysicsLogVector(LowestKineticEnergy,
386                                                         HighestKineticEnergy,
387                                                         TotBin);
388
389    // get material parameters needed for the energy loss calculation
390    const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(j);
391    const G4Material* material= couple->GetMaterial();
392
393    if ( charge > 0.0 ) {
394      ionloss = ProtonParametrisedDEDX(couple,highEnergy) ;
395    } else {
396      ionloss = AntiProtonParametrisedDEDX(couple,highEnergy) ;
397    }
398
399    ionlossBB = betheBlochModel->TheValue(&particleDef,material,highEnergy) ;
400    ionlossBB -= DeltaRaysEnergy(couple,highEnergy,proton_mass_c2) ;
401
402
403    paramB =  ionloss/ionlossBB - 1.0 ;
404
405    // now comes the loop for the kinetic energy values
406    for (G4int i = 0 ; i < TotBin ; i++) {
407      lowEdgeEnergy = aVector->GetLowEdgeEnergy(i) ;
408
409      // low energy part for this material, parametrised energy loss formulae
410      if ( lowEdgeEnergy < highEnergy ) {
411
412        if ( charge > 0.0 ) {
413          ionloss = ProtonParametrisedDEDX(couple,lowEdgeEnergy) ;
414        } else {
415          ionloss = AntiProtonParametrisedDEDX(couple,lowEdgeEnergy) ;
416        }
417
418      } else {
419
420        // high energy part for this material, Bethe-Bloch formula
421        ionloss = betheBlochModel->TheValue(proton,material,
422                                            lowEdgeEnergy) ;
423
424        ionloss -= DeltaRaysEnergy(couple,lowEdgeEnergy,proton_mass_c2) ;
425
426        ionloss *= (1.0 + paramB*highEnergy/lowEdgeEnergy) ;
427      }
428
429      // now put the loss into the vector
430      if(verboseLevel > 1) {
431        G4cout << "E(MeV)= " << lowEdgeEnergy/MeV
432               << "  dE/dx(MeV/mm)= " << ionloss*mm/MeV
433               << " in " << material->GetName() << G4endl;
434      }
435      aVector->PutValue(i,ionloss) ;
436    }
437    // Insert vector for this material into the table
438    theLossTable->insert(aVector) ;
439  }
440}
441
442
443
444// --------------------------------------------------------------------
445void G4hImpactIonisation::BuildLambdaTable(const G4ParticleDefinition& particleDef)
446
447{
448  // Build mean free path tables for the delta ray production process
449  //     tables are built for MATERIALS
450
451  if(verboseLevel > 1) {
452    G4cout << "G4hImpactIonisation::BuildLambdaTable for "
453           << particleDef.GetParticleName() << " is started" << G4endl;
454  }
455
456
457  G4double lowEdgeEnergy, value;
458  charge = particleDef.GetPDGCharge()/eplus ;
459  chargeSquare = charge*charge ;
460  initialMass = particleDef.GetPDGMass();
461
462  const G4ProductionCutsTable* theCoupleTable=
463    G4ProductionCutsTable::GetProductionCutsTable();
464  size_t numOfCouples = theCoupleTable->GetTableSize();
465
466
467  if (theMeanFreePathTable) {
468    theMeanFreePathTable->clearAndDestroy();
469    delete theMeanFreePathTable;
470  }
471
472  theMeanFreePathTable = new G4PhysicsTable(numOfCouples);
473
474  // loop for materials
475
476  for (size_t j=0 ; j < numOfCouples; j++) {
477
478    //create physics vector then fill it ....
479    G4PhysicsLogVector* aVector = new G4PhysicsLogVector(LowestKineticEnergy,
480                                                         HighestKineticEnergy,
481                                                         TotBin);
482
483    // compute the (macroscopic) cross section first
484    const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(j);
485    const G4Material* material= couple->GetMaterial();
486
487    const G4ElementVector* theElementVector =  material->GetElementVector() ;
488    const G4double* theAtomicNumDensityVector = material->GetAtomicNumDensityVector();
489    const G4int numberOfElements = material->GetNumberOfElements() ;
490
491    // get the electron kinetic energy cut for the actual material,
492    //  it will be used in ComputeMicroscopicCrossSection
493    // ( it is the SAME for ALL the ELEMENTS in THIS MATERIAL )
494    //   ------------------------------------------------------
495
496    G4double deltaCut = cutForDelta[j];
497
498    for ( G4int i = 0 ; i < TotBin ; i++ ) {
499      lowEdgeEnergy = aVector->GetLowEdgeEnergy(i) ;
500      G4double sigma = 0.0 ;
501      G4int Z;
502     
503      for (G4int iel=0; iel<numberOfElements; iel++ ) 
504        {
505          Z = (G4int) (*theElementVector)[iel]->GetZ();
506          // ---- MGP --- Corrected duplicated cross section calculation here from
507          // G4hLowEnergyIonisation original
508          G4double microCross = MicroscopicCrossSection( particleDef,
509                                                         lowEdgeEnergy,
510                                                         Z,
511                                                         deltaCut ) ;   
512          //totalCrossSectionMap [Z] = microCross;
513          sigma += theAtomicNumDensityVector[iel] * microCross ; 
514        }
515
516      // mean free path = 1./macroscopic cross section
517
518      value = sigma<=0 ? DBL_MAX : 1./sigma ;
519
520      aVector->PutValue(i, value) ;
521    }
522
523    theMeanFreePathTable->insert(aVector);
524  }
525
526}
527
528
529// --------------------------------------------------------------------
530G4double G4hImpactIonisation::MicroscopicCrossSection(const G4ParticleDefinition& particleDef,
531                                                      G4double kineticEnergy,
532                                                      G4double atomicNumber,
533                                                      G4double deltaCutInEnergy) const
534{
535  //******************************************************************
536    // cross section formula is OK for spin=0, 1/2, 1 only !
537    // *****************************************************************
538
539    // Calculates the microscopic cross section in GEANT4 internal units
540    // Formula documented in Geant4 Phys. Ref. Manual
541    // ( it is called for elements, AtomicNumber = z )
542
543    G4double totalCrossSection = 0.;
544
545  // Particle mass and energy
546  G4double particleMass = initialMass;
547  G4double energy = kineticEnergy + particleMass;
548
549  // Some kinematics
550  G4double gamma = energy / particleMass;
551  G4double beta2 = 1. - 1. / (gamma * gamma);
552  G4double var = electron_mass_c2 / particleMass;
553  G4double tMax   = 2. * electron_mass_c2 * (gamma*gamma - 1.) / (1. + 2.* gamma*var + var*var);
554
555  // Calculate the total cross section
556
557  if ( tMax > deltaCutInEnergy ) 
558    {
559      var = deltaCutInEnergy / tMax;
560      totalCrossSection = (1. - var * (1. - beta2 * std::log(var))) / deltaCutInEnergy ;
561     
562      G4double spin = particleDef.GetPDGSpin() ;
563     
564      // +term for spin=1/2 particle
565      if (spin == 0.5) 
566        {
567          totalCrossSection +=  0.5 * (tMax - deltaCutInEnergy) / (energy*energy);
568        }
569      // +term for spin=1 particle
570      else if (spin > 0.9 )
571        {
572          totalCrossSection += -std::log(var) / 
573            (3. * deltaCutInEnergy) + (tMax - deltaCutInEnergy) * ( (5. + 1. /var)*0.25 / (energy*energy) -
574                                                                    beta2 / (tMax * deltaCutInEnergy) ) / 3. ;
575        }
576      totalCrossSection *= twopi_mc2_rcl2 * atomicNumber / beta2 ;
577    }
578
579  //std::cout << "Microscopic = " << totalCrossSection/barn
580  //    << ", e = " << kineticEnergy/MeV <<std:: endl;
581
582  return totalCrossSection ;
583}
584
585
586
587// --------------------------------------------------------------------
588G4double G4hImpactIonisation::GetMeanFreePath(const G4Track& track,
589                                              G4double, // previousStepSize
590                                              enum G4ForceCondition* condition)
591{
592  const G4DynamicParticle* dynamicParticle = track.GetDynamicParticle();
593  const G4MaterialCutsCouple* couple = track.GetMaterialCutsCouple();
594  const G4Material* material = couple->GetMaterial();
595
596  G4double meanFreePath = DBL_MAX;
597  // ---- MGP ---- What is the meaning of the local variable isOutOfRange?
598  G4bool isOutRange = false;
599
600  *condition = NotForced ;
601
602  G4double kineticEnergy = (dynamicParticle->GetKineticEnergy())*initialMass/(dynamicParticle->GetMass());
603  charge = dynamicParticle->GetCharge()/eplus;
604  chargeSquare = theIonEffChargeModel->TheValue(dynamicParticle, material);
605
606  if (kineticEnergy < LowestKineticEnergy) 
607    {
608      meanFreePath = DBL_MAX;
609    }
610  else 
611    {
612      if (kineticEnergy > HighestKineticEnergy)  kineticEnergy = HighestKineticEnergy;
613      meanFreePath = (((*theMeanFreePathTable)(couple->GetIndex()))->
614                      GetValue(kineticEnergy,isOutRange))/chargeSquare;
615    }
616
617  return meanFreePath ;
618}
619
620
621// --------------------------------------------------------------------
622G4double G4hImpactIonisation::GetConstraints(const G4DynamicParticle* particle,
623                                             const G4MaterialCutsCouple* couple)
624{
625  // returns the Step limit
626  // dEdx is calculated as well as the range
627  // based on Effective Charge Approach
628
629  const G4Material* material = couple->GetMaterial();
630  G4Proton* proton = G4Proton::Proton();
631  G4AntiProton* antiproton = G4AntiProton::AntiProton();
632
633  G4double stepLimit = 0.;
634  G4double dx, highEnergy;
635
636  G4double massRatio = proton_mass_c2/(particle->GetMass()) ;
637  G4double kineticEnergy = particle->GetKineticEnergy() ;
638
639  // Scale the kinetic energy
640
641  G4double tScaled = kineticEnergy*massRatio ;
642  fBarkas = 0.;
643
644  if (charge > 0.) 
645    {
646      highEnergy = protonHighEnergy ;
647      fRangeNow = G4EnergyLossTables::GetRange(proton, tScaled, couple);
648      dx = G4EnergyLossTables::GetRange(proton, highEnergy, couple);
649      fdEdx = G4EnergyLossTables::GetDEDX(proton, tScaled, couple)
650        * chargeSquare ;
651     
652      // Correction for positive ions
653      if (theBarkas && tScaled > highEnergy) 
654        { 
655          fBarkas = BarkasTerm(material,tScaled)*std::sqrt(chargeSquare)*chargeSquare
656            + BlochTerm(material,tScaled,chargeSquare);
657        }
658      // Antiprotons and negative hadrons
659    } 
660  else 
661    {
662      highEnergy = antiprotonHighEnergy ;
663      fRangeNow = G4EnergyLossTables::GetRange(antiproton, tScaled, couple);
664      dx = G4EnergyLossTables::GetRange(antiproton, highEnergy, couple);
665      fdEdx = G4EnergyLossTables::GetDEDX(antiproton, tScaled, couple) * chargeSquare ;
666     
667      if (theBarkas && tScaled > highEnergy) 
668        { 
669          fBarkas = -BarkasTerm(material,tScaled)*std::sqrt(chargeSquare)*chargeSquare
670            + BlochTerm(material,tScaled,chargeSquare);
671        } 
672    } 
673  /*
674    const G4Material* mat = couple->GetMaterial();
675    G4double fac = gram/(MeV*cm2*mat->GetDensity());
676    G4cout << particle->GetDefinition()->GetParticleName()
677    << " in " << mat->GetName()
678    << " E(MeV)= " << kineticEnergy/MeV
679    << " dedx(MeV*cm^2/g)= " << fdEdx*fac
680    << " barcas(MeV*cm^2/gram)= " << fBarkas*fac
681    << " Q^2= " << chargeSquare
682    << G4endl;
683  */
684  // scaling back
685  fRangeNow /= (chargeSquare*massRatio) ;
686  dx        /= (chargeSquare*massRatio) ;
687
688  stepLimit  = fRangeNow ;
689  G4double r = std::min(finalRange, couple->GetProductionCuts()
690                        ->GetProductionCut(idxG4ElectronCut));
691
692  if (fRangeNow > r) 
693    { 
694      stepLimit = dRoverRange*fRangeNow + r*(1.0 - dRoverRange)*(2.0 - r/fRangeNow);
695      if (stepLimit > fRangeNow) stepLimit = fRangeNow;
696    } 
697  //   compute the (random) Step limit in standard energy range
698  if(tScaled > highEnergy ) 
699    {   
700      // add Barkas correction directly to dedx
701      fdEdx  += fBarkas;
702     
703      if(stepLimit > fRangeNow - dx*0.9) stepLimit = fRangeNow - dx*0.9 ;
704     
705      // Step limit in low energy range
706    } 
707  else 
708    {
709      G4double x = dx*paramStepLimit;
710      if (stepLimit > x) stepLimit = x;
711    }
712  return stepLimit;
713}
714
715
716// --------------------------------------------------------------------
717G4VParticleChange* G4hImpactIonisation::AlongStepDoIt(const G4Track& track,
718                                                      const G4Step& step)
719{
720  // compute the energy loss after a step
721  G4Proton* proton = G4Proton::Proton();
722  G4AntiProton* antiproton = G4AntiProton::AntiProton();
723  G4double finalT = 0.;
724
725  aParticleChange.Initialize(track) ;
726
727  const G4MaterialCutsCouple* couple = track.GetMaterialCutsCouple();
728  const G4Material* material = couple->GetMaterial();
729
730  // get the actual (true) Step length from step
731  const G4double stepLength = step.GetStepLength() ;
732
733  const G4DynamicParticle* particle = track.GetDynamicParticle() ;
734
735  G4double kineticEnergy = particle->GetKineticEnergy() ;
736  G4double massRatio = proton_mass_c2/(particle->GetMass()) ;
737  G4double tScaled = kineticEnergy * massRatio ;
738  G4double eLoss = 0.0 ;
739  G4double nLoss = 0.0 ;
740
741
742  // very small particle energy
743  if (kineticEnergy < MinKineticEnergy) 
744    {
745      eLoss = kineticEnergy ;
746      // particle energy outside tabulated energy range
747    } 
748 
749  else if( kineticEnergy > HighestKineticEnergy) 
750    {
751      eLoss = stepLength * fdEdx ;
752      // big step
753    } 
754  else if (stepLength >= fRangeNow ) 
755    {
756      eLoss = kineticEnergy ;
757     
758      // tabulated range
759    } 
760  else 
761    {
762      // step longer than linear step limit
763      if(stepLength > linLossLimit * fRangeNow) 
764        {
765          G4double rScaled = fRangeNow * massRatio * chargeSquare ;
766          G4double sScaled = stepLength * massRatio * chargeSquare ;
767         
768          if(charge > 0.0) 
769            {
770              eLoss = G4EnergyLossTables::GetPreciseEnergyFromRange(proton,rScaled, couple) -
771                G4EnergyLossTables::GetPreciseEnergyFromRange(proton,rScaled-sScaled,couple) ;
772           
773            } 
774          else 
775            {
776              // Antiproton
777              eLoss = G4EnergyLossTables::GetPreciseEnergyFromRange(antiproton,rScaled,couple) -
778                G4EnergyLossTables::GetPreciseEnergyFromRange(antiproton,rScaled-sScaled,couple) ;
779            }
780          eLoss /= massRatio ;
781         
782          // Barkas correction at big step     
783          eLoss += fBarkas * stepLength;
784         
785          // step shorter than linear step limit
786        } 
787      else 
788        {
789          eLoss = stepLength *fdEdx  ;
790        }
791      if (nStopping && tScaled < protonHighEnergy) 
792        {
793          nLoss = (theNuclearStoppingModel->TheValue(particle, material)) * stepLength;
794        }
795    }
796 
797  if (eLoss < 0.0) eLoss = 0.0;
798
799  finalT = kineticEnergy - eLoss - nLoss;
800
801  if ( EnlossFlucFlag && 0.0 < eLoss && finalT > MinKineticEnergy) 
802    {
803     
804      //  now the electron loss with fluctuation
805      eLoss = ElectronicLossFluctuation(particle, couple, eLoss, stepLength) ;
806      if (eLoss < 0.0) eLoss = 0.0;
807      finalT = kineticEnergy - eLoss - nLoss;
808    }
809
810  //  stop particle if the kinetic energy <= MinKineticEnergy
811  if (finalT*massRatio <= MinKineticEnergy ) 
812    {
813     
814      finalT = 0.0;
815      if (!particle->GetDefinition()->GetProcessManager()->GetAtRestProcessVector()->size())
816        aParticleChange.ProposeTrackStatus(fStopAndKill);
817      else
818        aParticleChange.ProposeTrackStatus(fStopButAlive);
819    }
820
821  aParticleChange.ProposeEnergy( finalT );
822  eLoss = kineticEnergy-finalT;
823
824  aParticleChange.ProposeLocalEnergyDeposit(eLoss);
825  return &aParticleChange ;
826}
827
828
829
830// --------------------------------------------------------------------
831G4double G4hImpactIonisation::ProtonParametrisedDEDX(const G4MaterialCutsCouple* couple,
832                                                     G4double kineticEnergy) const
833{
834  const G4Material* material = couple->GetMaterial();
835  G4Proton* proton = G4Proton::Proton();
836  G4double eLoss = 0.;
837
838  // Free Electron Gas Model
839  if(kineticEnergy < protonLowEnergy) {
840    eLoss = (protonModel->TheValue(proton, material, protonLowEnergy))
841      * std::sqrt(kineticEnergy/protonLowEnergy) ;
842
843    // Parametrisation
844  } else {
845    eLoss = protonModel->TheValue(proton, material, kineticEnergy) ;
846  }
847
848  // Delta rays energy
849  eLoss -= DeltaRaysEnergy(couple,kineticEnergy,proton_mass_c2) ;
850
851  if(verboseLevel > 2) {
852    G4cout << "p E(MeV)= " << kineticEnergy/MeV
853           << " dE/dx(MeV/mm)= " << eLoss*mm/MeV
854           << " for " << material->GetName()
855           << " model: " << protonModel << G4endl;
856  }
857
858  if(eLoss < 0.0) eLoss = 0.0 ;
859
860  return eLoss ;
861}
862
863
864
865// --------------------------------------------------------------------
866G4double G4hImpactIonisation::AntiProtonParametrisedDEDX(const G4MaterialCutsCouple* couple,
867                                                         G4double kineticEnergy) const
868{
869  const G4Material* material = couple->GetMaterial();
870  G4AntiProton* antiproton = G4AntiProton::AntiProton();
871  G4double eLoss = 0.0 ;
872
873  // Antiproton model is used
874  if(antiprotonModel->IsInCharge(antiproton,material)) {
875    if(kineticEnergy < antiprotonLowEnergy) {
876      eLoss = antiprotonModel->TheValue(antiproton,material,antiprotonLowEnergy)
877        * std::sqrt(kineticEnergy/antiprotonLowEnergy) ;
878
879      // Parametrisation
880    } else {
881      eLoss = antiprotonModel->TheValue(antiproton,material,
882                                        kineticEnergy);
883    }
884
885    // The proton model is used + Barkas correction
886  } else {
887    if(kineticEnergy < protonLowEnergy) {
888      eLoss = protonModel->TheValue(G4Proton::Proton(),material,protonLowEnergy)
889        * std::sqrt(kineticEnergy/protonLowEnergy) ;
890
891      // Parametrisation
892    } else {
893      eLoss = protonModel->TheValue(G4Proton::Proton(),material,
894                                    kineticEnergy);
895    }
896    //if(theBarkas) eLoss -= 2.0*BarkasTerm(material, kineticEnergy);
897  }
898
899  // Delta rays energy
900  eLoss -= DeltaRaysEnergy(couple,kineticEnergy,proton_mass_c2) ;
901
902  if(verboseLevel > 2) {
903    G4cout << "pbar E(MeV)= " << kineticEnergy/MeV
904           << " dE/dx(MeV/mm)= " << eLoss*mm/MeV
905           << " for " << material->GetName()
906           << " model: " << protonModel << G4endl;
907  }
908
909  if(eLoss < 0.0) eLoss = 0.0 ;
910
911  return eLoss ;
912}
913
914
915// --------------------------------------------------------------------
916G4double G4hImpactIonisation::DeltaRaysEnergy(const G4MaterialCutsCouple* couple,
917                                              G4double kineticEnergy,
918                                              G4double particleMass) const
919{
920  G4double dLoss = 0.;
921
922  G4double deltaCutNow = cutForDelta[(couple->GetIndex())] ;
923  const G4Material* material = couple->GetMaterial();
924  G4double electronDensity = material->GetElectronDensity();
925  G4double excitationEnergy = material->GetIonisation()->GetMeanExcitationEnergy();
926
927  G4double tau = kineticEnergy / particleMass ;
928  G4double rateMass = electron_mass_c2/particleMass ;
929
930  // some local variables
931
932  G4double gamma = tau + 1.0 ;
933  G4double bg2 = tau*(tau+2.0) ;
934  G4double beta2 = bg2/(gamma*gamma) ;
935  G4double tMax = 2.*electron_mass_c2*bg2/(1.0+2.0*gamma*rateMass+rateMass*rateMass) ;
936
937  // Validity range for delta electron cross section
938  G4double deltaCut = std::max(deltaCutNow, excitationEnergy);
939
940  if ( deltaCut < tMax) 
941    {
942      G4double x = deltaCut / tMax ;
943      dLoss = ( beta2 * (x-1.) - std::log(x) ) * twopi_mc2_rcl2 * electronDensity / beta2 ;
944    }
945  return dLoss ;
946}
947
948
949// -------------------------------------------------------------------------
950
951G4VParticleChange* G4hImpactIonisation::PostStepDoIt(const G4Track& track,
952                                                     const G4Step& step)
953{
954  // Units are expressed in GEANT4 internal units.
955
956  //  std::cout << "----- Calling PostStepDoIt ----- " << std::endl;
957
958  aParticleChange.Initialize(track) ;
959  const G4MaterialCutsCouple* couple = track.GetMaterialCutsCouple(); 
960  const G4DynamicParticle* aParticle = track.GetDynamicParticle() ;
961 
962  // Some kinematics
963
964  G4ParticleDefinition* definition = track.GetDefinition();
965  G4double mass = definition->GetPDGMass();
966  G4double kineticEnergy = aParticle->GetKineticEnergy();
967  G4double totalEnergy = kineticEnergy + mass ;
968  G4double pSquare = kineticEnergy *( totalEnergy + mass) ;
969  G4double eSquare = totalEnergy * totalEnergy;
970  G4double betaSquare = pSquare / eSquare;
971  G4ThreeVector particleDirection = aParticle->GetMomentumDirection() ;
972
973  G4double gamma = kineticEnergy / mass + 1.;
974  G4double r = electron_mass_c2 / mass;
975  G4double tMax = 2. * electron_mass_c2 *(gamma*gamma - 1.) / (1. + 2.*gamma*r + r*r);
976
977  // Validity range for delta electron cross section
978  G4double deltaCut = cutForDelta[couple->GetIndex()];
979
980  // This should not be a case
981  if (deltaCut >= tMax)
982    return G4VContinuousDiscreteProcess::PostStepDoIt(track,step);
983
984  G4double xc = deltaCut / tMax;
985  G4double rate = tMax / totalEnergy;
986  rate = rate*rate ;
987  G4double spin = aParticle->GetDefinition()->GetPDGSpin() ;
988
989  // Sampling follows ...
990  G4double x = 0.;
991  G4double gRej = 0.;
992
993  do {
994    x = xc / (1. - (1. - xc) * G4UniformRand());
995   
996    if (0.0 == spin) 
997      {
998        gRej = 1.0 - betaSquare * x ;
999      }
1000    else if (0.5 == spin) 
1001      {
1002        gRej = (1. - betaSquare * x + 0.5 * x*x * rate) / (1. + 0.5 * rate) ;
1003      } 
1004    else 
1005      {
1006        gRej = (1. - betaSquare * x ) * (1. + x/(3.*xc)) +
1007          x*x * rate * (1. + 0.5*x/xc) / 3.0 /
1008          (1. + 1./(3.*xc) + rate *(1.+ 0.5/xc) / 3.);
1009      }
1010   
1011  } while( G4UniformRand() > gRej );
1012
1013  G4double deltaKineticEnergy = x * tMax;
1014  G4double deltaTotalMomentum = std::sqrt(deltaKineticEnergy * 
1015                                          (deltaKineticEnergy + 2. * electron_mass_c2 ));
1016  G4double totalMomentum = std::sqrt(pSquare) ;
1017  G4double cosTheta = deltaKineticEnergy * (totalEnergy + electron_mass_c2) / (deltaTotalMomentum*totalMomentum);
1018
1019  //  protection against cosTheta > 1 or < -1
1020  if ( cosTheta < -1. ) cosTheta = -1.;
1021  if ( cosTheta > 1. ) cosTheta = 1.;
1022
1023  //  direction of the delta electron
1024  G4double phi = twopi * G4UniformRand();
1025  G4double sinTheta = std::sqrt(1. - cosTheta*cosTheta);
1026  G4double dirX = sinTheta * std::cos(phi);
1027  G4double dirY = sinTheta * std::sin(phi);
1028  G4double dirZ = cosTheta;
1029
1030  G4ThreeVector deltaDirection(dirX,dirY,dirZ);
1031  deltaDirection.rotateUz(particleDirection);
1032
1033  // create G4DynamicParticle object for delta ray
1034  G4DynamicParticle* deltaRay = new G4DynamicParticle;
1035  deltaRay->SetKineticEnergy( deltaKineticEnergy );
1036  deltaRay->SetMomentumDirection(deltaDirection.x(),
1037                                 deltaDirection.y(),
1038                                 deltaDirection.z());
1039  deltaRay->SetDefinition(G4Electron::Electron());
1040
1041  // fill aParticleChange
1042  G4double finalKineticEnergy = kineticEnergy - deltaKineticEnergy;
1043  size_t totalNumber = 1;
1044
1045  // Atomic relaxation
1046
1047  // ---- MGP ---- Temporary limitation: currently PIXE only for incident protons
1048
1049  size_t nSecondaries = 0;
1050  std::vector<G4DynamicParticle*>* secondaryVector = 0;
1051
1052  if (definition == G4Proton::ProtonDefinition())
1053    {
1054      const G4Material* material = couple->GetMaterial();
1055
1056      // Lazy initialization of pixeCrossSectionHandler
1057      if (pixeCrossSectionHandler == 0)
1058        {
1059          // Instantiate pixeCrossSectionHandler with selected shell cross section models
1060          // Ownership of interpolation is transferred to pixeCrossSectionHandler
1061          G4IInterpolator* interpolation = new G4LogLogInterpolator();
1062          pixeCrossSectionHandler = 
1063            new G4PixeCrossSectionHandler(interpolation,modelK,modelL,modelM,eMinPixe,eMaxPixe);
1064          G4String fileName("proton");
1065          pixeCrossSectionHandler->LoadShellData(fileName);
1066          //      pixeCrossSectionHandler->PrintData();
1067        }
1068     
1069      // Select an atom in the current material based on the total shell cross sections
1070      G4int Z = pixeCrossSectionHandler->SelectRandomAtom(material,kineticEnergy);
1071      //      std::cout << "G4hImpactIonisation::PostStepDoIt - Z = " << Z << std::endl;
1072
1073      //      G4double microscopicCross = MicroscopicCrossSection(*definition,
1074      //                                          kineticEnergy,
1075      //                                          Z, deltaCut);
1076  //    G4double crossFromShells = pixeCrossSectionHandler->FindValue(Z,kineticEnergy);
1077
1078      //std::cout << "G4hImpactIonisation: Z= "
1079      //                << Z
1080      //                << ", energy = "
1081      //                << kineticEnergy/MeV
1082      //                <<" MeV, microscopic = "
1083      //                << microscopicCross/barn
1084      //                << " barn, from shells = "
1085      //                << crossFromShells/barn
1086      //                << " barn"
1087      //                << std::endl;
1088
1089      // Select a shell in the target atom based on the individual shell cross sections
1090      G4int shellIndex = pixeCrossSectionHandler->SelectRandomShell(Z,kineticEnergy);
1091   
1092      G4AtomicTransitionManager* transitionManager = G4AtomicTransitionManager::Instance();
1093      const G4AtomicShell* atomicShell = transitionManager->Shell(Z,shellIndex);
1094      G4double bindingEnergy = atomicShell->BindingEnergy();
1095     
1096      //     if (verboseLevel > 1)
1097      //       {
1098      //         G4cout << "G4hImpactIonisation::PostStepDoIt - Z = "
1099      //         << Z
1100      //         << ", shell = "
1101      //         << shellIndex
1102      //         << ", bindingE (keV) = "
1103      //         << bindingEnergy/keV
1104      //         << G4endl;
1105      //       }
1106     
1107      // Generate PIXE if binding energy larger than cut for photons or electrons
1108     
1109      G4ParticleDefinition* type = 0;
1110     
1111      if (finalKineticEnergy >= bindingEnergy)
1112        //        && (bindingEnergy >= minGammaEnergy ||  bindingEnergy >= minElectronEnergy) )
1113        {
1114          // Vacancy in subshell shellIndex; shellId is the subshell identifier in EADL jargon
1115          G4int shellId = atomicShell->ShellId();
1116          // Atomic relaxation: generate secondaries
1117          secondaryVector = atomicDeexcitation.GenerateParticles(Z, shellId);
1118
1119          // ---- Debug ----
1120          //std::cout << "ShellId = "
1121          //    <<shellId << " ---- Atomic relaxation secondaries: ---- "
1122          //    << secondaryVector->size()
1123          //    << std::endl;
1124
1125          // ---- End debug ---
1126
1127          if (secondaryVector != 0) 
1128            {
1129              nSecondaries = secondaryVector->size();
1130              for (size_t i = 0; i<nSecondaries; i++) 
1131                {
1132                  G4DynamicParticle* aSecondary = (*secondaryVector)[i];
1133                  if (aSecondary) 
1134                    {
1135                      G4double e = aSecondary->GetKineticEnergy();
1136                      type = aSecondary->GetDefinition();
1137
1138                      // ---- Debug ----
1139                      //if (type == G4Gamma::GammaDefinition())
1140                      //        {                       
1141                      //          std::cout << "Z = " << Z
1142                      //                    << ", shell: " << shellId
1143                      //                    << ", PIXE photon energy (keV) = " << e/keV
1144                      //                    << std::endl;
1145                      //        }
1146                      // ---- End debug ---
1147
1148                      if (e < finalKineticEnergy &&
1149                          ((type == G4Gamma::Gamma() && e > minGammaEnergy ) ||
1150                           (type == G4Electron::Electron() && e > minElectronEnergy ))) 
1151                        {
1152                          // Subtract the energy of the emitted secondary from the primary
1153                          finalKineticEnergy -= e;
1154                          totalNumber++;
1155                          // ---- Debug ----
1156                          //if (type == G4Gamma::GammaDefinition())
1157                          //    {                       
1158                          //      std::cout << "Z = " << Z
1159                          //                << ", shell: " << shellId
1160                          //                << ", PIXE photon energy (keV) = " << e/keV
1161                          //                << std::endl;
1162                          //    }
1163                          // ---- End debug ---
1164                        } 
1165                      else 
1166                        {
1167                          // The atomic relaxation product has energy below the cut
1168                          // ---- Debug ----
1169                          // if (type == G4Gamma::GammaDefinition())
1170                          //    {                       
1171                          //      std::cout << "Z = " << Z
1172                          //
1173                          //                << ", PIXE photon energy = " << e/keV
1174                          //                << " keV below threshold " << minGammaEnergy/keV << " keV"
1175                          //                << std::endl;
1176                          //    }   
1177                          // ---- End debug ---
1178     
1179                          delete aSecondary;
1180                          (*secondaryVector)[i] = 0;
1181                        }
1182                    }
1183                }
1184            }
1185        }
1186    }
1187
1188
1189  // Save delta-electrons
1190
1191  G4double eDeposit = 0.;
1192
1193  if (finalKineticEnergy > MinKineticEnergy)
1194    {
1195      G4double finalPx = totalMomentum*particleDirection.x() - deltaTotalMomentum*deltaDirection.x();
1196      G4double finalPy = totalMomentum*particleDirection.y() - deltaTotalMomentum*deltaDirection.y();
1197      G4double finalPz = totalMomentum*particleDirection.z() - deltaTotalMomentum*deltaDirection.z();
1198      G4double finalMomentum = std::sqrt(finalPx*finalPx + finalPy*finalPy + finalPz*finalPz) ;
1199      finalPx /= finalMomentum;
1200      finalPy /= finalMomentum;
1201      finalPz /= finalMomentum;
1202
1203      aParticleChange.ProposeMomentumDirection( finalPx,finalPy,finalPz );
1204    }
1205  else
1206    {
1207      eDeposit = finalKineticEnergy;
1208      finalKineticEnergy = 0.;
1209      aParticleChange.ProposeMomentumDirection(particleDirection.x(),
1210                                               particleDirection.y(),
1211                                               particleDirection.z());
1212      if(!aParticle->GetDefinition()->GetProcessManager()->
1213         GetAtRestProcessVector()->size())
1214        aParticleChange.ProposeTrackStatus(fStopAndKill);
1215      else
1216        aParticleChange.ProposeTrackStatus(fStopButAlive);
1217    }
1218
1219  aParticleChange.ProposeEnergy(finalKineticEnergy);
1220  aParticleChange.ProposeLocalEnergyDeposit (eDeposit);
1221  aParticleChange.SetNumberOfSecondaries(totalNumber);
1222  aParticleChange.AddSecondary(deltaRay);
1223
1224  // ---- Debug ----
1225  //  std::cout << "RDHadronIonisation - finalKineticEnergy (MeV) = "
1226  //        << finalKineticEnergy/MeV
1227  //        << ", delta KineticEnergy (keV) = "
1228  //        << deltaKineticEnergy/keV
1229  //        << ", energy deposit (MeV) = "
1230  //        << eDeposit/MeV
1231  //        << std::endl;
1232  // ---- End debug ---
1233 
1234  // Save Fluorescence and Auger
1235
1236  if (secondaryVector != 0) 
1237    { 
1238      for (size_t l = 0; l < nSecondaries; l++) 
1239        { 
1240          G4DynamicParticle* secondary = (*secondaryVector)[l];
1241          if (secondary) aParticleChange.AddSecondary(secondary);
1242
1243          // ---- Debug ----
1244          //if (secondary != 0)
1245          // {
1246          //   if (secondary->GetDefinition() == G4Gamma::GammaDefinition())
1247          //    {
1248          //      G4double eX = secondary->GetKineticEnergy();                 
1249          //      std::cout << " PIXE photon of energy " << eX/keV
1250          //                << " keV added to ParticleChange; total number of secondaries is " << totalNumber
1251          //                << std::endl;
1252          //    }
1253          //}
1254          // ---- End debug ---   
1255
1256        }
1257      delete secondaryVector;
1258    }
1259
1260  return G4VContinuousDiscreteProcess::PostStepDoIt(track,step);
1261}
1262
1263// -------------------------------------------------------------------------
1264
1265G4double G4hImpactIonisation::ComputeDEDX(const G4ParticleDefinition* aParticle,
1266                                          const G4MaterialCutsCouple* couple,
1267                                          G4double kineticEnergy)
1268{
1269  const G4Material* material = couple->GetMaterial();
1270  G4Proton* proton = G4Proton::Proton();
1271  G4AntiProton* antiproton = G4AntiProton::AntiProton();
1272  G4double dedx = 0.;
1273
1274  G4double tScaled = kineticEnergy * proton_mass_c2 / (aParticle->GetPDGMass()) ;
1275  charge  = aParticle->GetPDGCharge() ;
1276
1277  if( charge > 0.) 
1278    {
1279      if (tScaled > protonHighEnergy) 
1280        {
1281          dedx = G4EnergyLossTables::GetDEDX(proton,tScaled,couple) ; 
1282        }
1283      else 
1284        {
1285          dedx = ProtonParametrisedDEDX(couple,tScaled) ;
1286        }
1287    } 
1288  else 
1289    {
1290      if (tScaled > antiprotonHighEnergy) 
1291        {
1292          dedx = G4EnergyLossTables::GetDEDX(antiproton,tScaled,couple);
1293        } 
1294      else
1295        {
1296          dedx = AntiProtonParametrisedDEDX(couple,tScaled) ;
1297        }
1298    }
1299  dedx *= theIonEffChargeModel->TheValue(aParticle, material, kineticEnergy) ;
1300 
1301  return dedx ;
1302}
1303
1304
1305G4double G4hImpactIonisation::BarkasTerm(const G4Material* material,
1306                                         G4double kineticEnergy) const
1307//Function to compute the Barkas term for protons:
1308//
1309//Ref. Z_1^3 effect in the stopping power of matter for charged particles
1310//     J.C Ashley and R.H.Ritchie
1311//     Physical review B Vol.5 No.7 1 April 1972 pagg. 2393-2397
1312//
1313{
1314  static double FTable[47][2] = {
1315    { 0.02, 21.5},
1316    { 0.03, 20.0},
1317    { 0.04, 18.0},
1318    { 0.05, 15.6},
1319    { 0.06, 15.0},
1320    { 0.07, 14.0},
1321    { 0.08, 13.5},
1322    { 0.09, 13.},
1323    { 0.1,  12.2},
1324    { 0.2,  9.25},
1325    { 0.3,  7.0},
1326    { 0.4,  6.0},
1327    { 0.5,  4.5},
1328    { 0.6,  3.5},
1329    { 0.7,  3.0},
1330    { 0.8,  2.5},
1331    { 0.9,  2.0},
1332    { 1.0,  1.7},
1333    { 1.2,  1.2},
1334    { 1.3,  1.0},
1335    { 1.4,  0.86},
1336    { 1.5,  0.7},
1337    { 1.6,  0.61},
1338    { 1.7,  0.52},
1339    { 1.8,  0.5},
1340    { 1.9,  0.43},
1341    { 2.0,  0.42},
1342    { 2.1,  0.3},
1343    { 2.4,  0.2},
1344    { 3.0,  0.13},
1345    { 3.08, 0.1},
1346    { 3.1,  0.09},
1347    { 3.3,  0.08},
1348    { 3.5,  0.07},
1349    { 3.8,  0.06},
1350    { 4.0,  0.051},
1351    { 4.1,  0.04},
1352    { 4.8,  0.03},
1353    { 5.0,  0.024},
1354    { 5.1,  0.02},
1355    { 6.0,  0.013},
1356    { 6.5,  0.01},
1357    { 7.0,  0.009},
1358    { 7.1,  0.008},
1359    { 8.0,  0.006},
1360    { 9.0,  0.0032},
1361    { 10.0, 0.0025} };
1362
1363  // Information on particle and material
1364  G4double kinE  = kineticEnergy ;
1365  if(0.5*MeV > kinE) kinE = 0.5*MeV ;
1366  G4double gamma = 1.0 + kinE / proton_mass_c2 ;
1367  G4double beta2 = 1.0 - 1.0/(gamma*gamma) ;
1368  if(0.0 >= beta2) return 0.0;
1369
1370  G4double BarkasTerm = 0.0;
1371  G4double AMaterial = 0.0;
1372  G4double ZMaterial = 0.0;
1373  const G4ElementVector* theElementVector = material->GetElementVector();
1374  G4int numberOfElements = material->GetNumberOfElements();
1375
1376  for (G4int i = 0; i<numberOfElements; i++) {
1377
1378    AMaterial = (*theElementVector)[i]->GetA()*mole/g;
1379    ZMaterial = (*theElementVector)[i]->GetZ();
1380
1381    G4double X = 137.0 * 137.0 * beta2 / ZMaterial;
1382
1383    // Variables to compute L_1
1384    G4double Eta0Chi = 0.8;
1385    G4double EtaChi = Eta0Chi * ( 1.0 + 6.02*std::pow( ZMaterial,-1.19 ) );
1386    G4double W = ( EtaChi * std::pow( ZMaterial,1.0/6.0 ) ) / std::sqrt(X);
1387    G4double FunctionOfW = FTable[46][1]*FTable[46][0]/W ;
1388
1389    for(G4int j=0; j<47; j++) {
1390
1391      if( W < FTable[j][0] ) {
1392
1393        if(0 == j) {
1394          FunctionOfW = FTable[0][1] ;
1395
1396        } else {
1397          FunctionOfW = (FTable[j][1] - FTable[j-1][1]) * (W - FTable[j-1][0])
1398            / (FTable[j][0] - FTable[j-1][0])
1399            +  FTable[j-1][1] ;
1400        }
1401
1402        break;
1403      }
1404
1405    }
1406
1407    BarkasTerm += FunctionOfW /( std::sqrt(ZMaterial * X) * X);
1408  }
1409
1410  BarkasTerm *= twopi_mc2_rcl2 * (material->GetElectronDensity()) / beta2 ;
1411
1412  return BarkasTerm;
1413}
1414
1415
1416
1417G4double G4hImpactIonisation::BlochTerm(const G4Material* material,
1418                                        G4double kineticEnergy,
1419                                        G4double cSquare) const
1420//Function to compute the Bloch term for protons:
1421//
1422//Ref. Z_1^3 effect in the stopping power of matter for charged particles
1423//     J.C Ashley and R.H.Ritchie
1424//     Physical review B Vol.5 No.7 1 April 1972 pagg. 2393-2397
1425//
1426{
1427  G4double eLoss = 0.0 ;
1428  G4double gamma = 1.0 + kineticEnergy / proton_mass_c2 ;
1429  G4double beta2 = 1.0 - 1.0/(gamma*gamma) ;
1430  G4double y = cSquare / (137.0*137.0*beta2) ;
1431
1432  if(y < 0.05) {
1433    eLoss = 1.202 ;
1434
1435  } else {
1436    eLoss = 1.0 / (1.0 + y) ;
1437    G4double de = eLoss ;
1438
1439    for(G4int i=2; de>eLoss*0.01; i++) {
1440      de = 1.0/( i * (i*i + y)) ;
1441      eLoss += de ;
1442    }
1443  }
1444  eLoss *= -1.0 * y * cSquare * twopi_mc2_rcl2 *
1445    (material->GetElectronDensity()) / beta2 ;
1446
1447  return eLoss;
1448}
1449
1450
1451
1452G4double G4hImpactIonisation::ElectronicLossFluctuation(
1453                                                        const G4DynamicParticle* particle,
1454                                                        const G4MaterialCutsCouple* couple,
1455                                                        G4double meanLoss,
1456                                                        G4double step) const
1457//  calculate actual loss from the mean loss
1458//  The model used to get the fluctuation is essentially the same
1459// as in Glandz in Geant3.
1460{
1461  // data members to speed up the fluctuation calculation
1462  //  G4int imat ;
1463  //  G4double f1Fluct,f2Fluct,e1Fluct,e2Fluct,rateFluct,ipotFluct;
1464  //  G4double e1LogFluct,e2LogFluct,ipotLogFluct;
1465
1466  static const G4double minLoss = 1.*eV ;
1467  static const G4double kappa = 10. ;
1468  static const G4double theBohrBeta2 = 50.0 * keV/proton_mass_c2 ;
1469
1470  const G4Material* material = couple->GetMaterial();
1471  G4int    imaterial   = couple->GetIndex() ;
1472  G4double ipotFluct   = material->GetIonisation()->GetMeanExcitationEnergy() ;
1473  G4double electronDensity = material->GetElectronDensity() ;
1474  G4double zeff = electronDensity/(material->GetTotNbOfAtomsPerVolume()) ;
1475
1476  // get particle data
1477  G4double tkin   = particle->GetKineticEnergy();
1478  G4double particleMass = particle->GetMass() ;
1479  G4double deltaCutInKineticEnergyNow = cutForDelta[imaterial];
1480
1481  // shortcut for very very small loss
1482  if(meanLoss < minLoss) return meanLoss ;
1483
1484  // Validity range for delta electron cross section
1485  G4double threshold = std::max(deltaCutInKineticEnergyNow,ipotFluct);
1486  G4double loss, siga;
1487
1488  G4double rmass = electron_mass_c2/particleMass;
1489  G4double tau   = tkin/particleMass;
1490  G4double tau1 = tau+1.0;
1491  G4double tau2 = tau*(tau+2.);
1492  G4double tMax    = 2.*electron_mass_c2*tau2/(1.+2.*tau1*rmass+rmass*rmass);
1493
1494
1495  if(tMax > threshold) tMax = threshold;
1496  G4double beta2 = tau2/(tau1*tau1);
1497
1498  // Gaussian fluctuation
1499  if(meanLoss > kappa*tMax || tMax < kappa*ipotFluct )
1500    {
1501      siga = tMax * (1.0-0.5*beta2) * step * twopi_mc2_rcl2
1502        * electronDensity / beta2 ;
1503
1504      // High velocity or negatively charged particle
1505      if( beta2 > 3.0*theBohrBeta2*zeff || charge < 0.0) {
1506        siga = std::sqrt( siga * chargeSquare ) ;
1507
1508        // Low velocity - additional ion charge fluctuations according to
1509        // Q.Yang et al., NIM B61(1991)149-155.
1510      } else {
1511        G4double chu = theIonChuFluctuationModel->TheValue(particle, material);
1512        G4double yang = theIonYangFluctuationModel->TheValue(particle, material);
1513        siga = std::sqrt( siga * (chargeSquare * chu + yang)) ;
1514      }
1515
1516      do {
1517        loss = G4RandGauss::shoot(meanLoss,siga);
1518      } while (loss < 0. || loss > 2.0*meanLoss);
1519      return loss;
1520    }
1521
1522  // Non Gaussian fluctuation
1523  static const G4double probLim = 0.01 ;
1524  static const G4double sumaLim = -std::log(probLim) ;
1525  static const G4double alim = 10.;
1526
1527  G4double suma,w1,w2,C,e0,lossc,w;
1528  G4double a1,a2,a3;
1529  G4int p1,p2,p3;
1530  G4int nb;
1531  G4double corrfac, na,alfa,rfac,namean,sa,alfa1,ea,sea;
1532  G4double dp3;
1533
1534  G4double f1Fluct     = material->GetIonisation()->GetF1fluct();
1535  G4double f2Fluct     = material->GetIonisation()->GetF2fluct();
1536  G4double e1Fluct     = material->GetIonisation()->GetEnergy1fluct();
1537  G4double e2Fluct     = material->GetIonisation()->GetEnergy2fluct();
1538  G4double e1LogFluct  = material->GetIonisation()->GetLogEnergy1fluct();
1539  G4double e2LogFluct  = material->GetIonisation()->GetLogEnergy2fluct();
1540  G4double rateFluct   = material->GetIonisation()->GetRateionexcfluct();
1541  G4double ipotLogFluct= material->GetIonisation()->GetLogMeanExcEnergy();
1542
1543  w1 = tMax/ipotFluct;
1544  w2 = std::log(2.*electron_mass_c2*tau2);
1545
1546  C = meanLoss*(1.-rateFluct)/(w2-ipotLogFluct-beta2);
1547
1548  a1 = C*f1Fluct*(w2-e1LogFluct-beta2)/e1Fluct;
1549  a2 = C*f2Fluct*(w2-e2LogFluct-beta2)/e2Fluct;
1550  a3 = rateFluct*meanLoss*(tMax-ipotFluct)/(ipotFluct*tMax*std::log(w1));
1551  if(a1 < 0.0) a1 = 0.0;
1552  if(a2 < 0.0) a2 = 0.0;
1553  if(a3 < 0.0) a3 = 0.0;
1554
1555  suma = a1+a2+a3;
1556
1557  loss = 0.;
1558
1559
1560  if(suma < sumaLim)             // very small Step
1561    {
1562      e0 = material->GetIonisation()->GetEnergy0fluct();
1563
1564      if(tMax == ipotFluct)
1565        {
1566          a3 = meanLoss/e0;
1567
1568          if(a3>alim)
1569            {
1570              siga=std::sqrt(a3) ;
1571              p3 = std::max(0,G4int(G4RandGauss::shoot(a3,siga)+0.5));
1572            }
1573          else
1574            p3 = G4Poisson(a3);
1575
1576          loss = p3*e0 ;
1577
1578          if(p3 > 0)
1579            loss += (1.-2.*G4UniformRand())*e0 ;
1580
1581        }
1582      else
1583        {
1584          tMax = tMax-ipotFluct+e0 ;
1585          a3 = meanLoss*(tMax-e0)/(tMax*e0*std::log(tMax/e0));
1586
1587          if(a3>alim)
1588            {
1589              siga=std::sqrt(a3) ;
1590              p3 = std::max(0,int(G4RandGauss::shoot(a3,siga)+0.5));
1591            }
1592          else
1593            p3 = G4Poisson(a3);
1594
1595          if(p3 > 0)
1596            {
1597              w = (tMax-e0)/tMax ;
1598              if(p3 > nmaxCont2)
1599                {
1600                  dp3 = G4float(p3) ;
1601                  corrfac = dp3/G4float(nmaxCont2) ;
1602                  p3 = nmaxCont2 ;
1603                }
1604              else
1605                corrfac = 1. ;
1606
1607              for(G4int i=0; i<p3; i++) loss += 1./(1.-w*G4UniformRand()) ;
1608              loss *= e0*corrfac ;
1609            }
1610        }
1611    }
1612
1613  else                              // not so small Step
1614    {
1615      // excitation type 1
1616      if(a1>alim)
1617        {
1618          siga=std::sqrt(a1) ;
1619          p1 = std::max(0,G4int(G4RandGauss::shoot(a1,siga)+0.5));
1620        }
1621      else
1622        p1 = G4Poisson(a1);
1623
1624      // excitation type 2
1625      if(a2>alim)
1626        {
1627          siga=std::sqrt(a2) ;
1628          p2 = std::max(0,G4int(G4RandGauss::shoot(a2,siga)+0.5));
1629        }
1630      else
1631        p2 = G4Poisson(a2);
1632
1633      loss = p1*e1Fluct+p2*e2Fluct;
1634
1635      // smearing to avoid unphysical peaks
1636      if(p2 > 0)
1637        loss += (1.-2.*G4UniformRand())*e2Fluct;
1638      else if (loss>0.)
1639        loss += (1.-2.*G4UniformRand())*e1Fluct;
1640
1641      // ionisation .......................................
1642      if(a3 > 0.)
1643        {
1644          if(a3>alim)
1645            {
1646              siga=std::sqrt(a3) ;
1647              p3 = std::max(0,G4int(G4RandGauss::shoot(a3,siga)+0.5));
1648            }
1649          else
1650            p3 = G4Poisson(a3);
1651
1652          lossc = 0.;
1653          if(p3 > 0)
1654            {
1655              na = 0.;
1656              alfa = 1.;
1657              if (p3 > nmaxCont2)
1658                {
1659                  dp3        = G4float(p3);
1660                  rfac       = dp3/(G4float(nmaxCont2)+dp3);
1661                  namean     = G4float(p3)*rfac;
1662                  sa         = G4float(nmaxCont1)*rfac;
1663                  na         = G4RandGauss::shoot(namean,sa);
1664                  if (na > 0.)
1665                    {
1666                      alfa   = w1*G4float(nmaxCont2+p3)/
1667                        (w1*G4float(nmaxCont2)+G4float(p3));
1668                      alfa1  = alfa*std::log(alfa)/(alfa-1.);
1669                      ea     = na*ipotFluct*alfa1;
1670                      sea    = ipotFluct*std::sqrt(na*(alfa-alfa1*alfa1));
1671                      lossc += G4RandGauss::shoot(ea,sea);
1672                    }
1673                }
1674
1675              nb = G4int(G4float(p3)-na);
1676              if (nb > 0)
1677                {
1678                  w2 = alfa*ipotFluct;
1679                  w  = (tMax-w2)/tMax;
1680                  for (G4int k=0; k<nb; k++) lossc += w2/(1.-w*G4UniformRand());
1681                }
1682            }
1683          loss += lossc;
1684        }
1685    }
1686
1687  return loss ;
1688}
1689
1690
1691
1692void G4hImpactIonisation::SetCutForSecondaryPhotons(G4double cut)
1693{
1694  minGammaEnergy = cut;
1695}
1696
1697
1698
1699void G4hImpactIonisation::SetCutForAugerElectrons(G4double cut)
1700{
1701  minElectronEnergy = cut;
1702}
1703
1704
1705
1706void G4hImpactIonisation::ActivateAugerElectronProduction(G4bool val)
1707{
1708  atomicDeexcitation.ActivateAugerElectronProduction(val);
1709}
1710
1711
1712
1713void G4hImpactIonisation::PrintInfoDefinition() const
1714{
1715  G4String comments = "  Knock-on electron cross sections . ";
1716  comments += "\n        Good description above the mean excitation energy.\n";
1717  comments += "        Delta ray energy sampled from  differential Xsection.";
1718
1719  G4cout << G4endl << GetProcessName() << ":  " << comments
1720         << "\n        PhysicsTables from " << LowestKineticEnergy / eV << " eV "
1721         << " to " << HighestKineticEnergy / TeV << " TeV "
1722         << " in " << TotBin << " bins."
1723         << "\n        Electronic stopping power model is  "
1724         << protonTable
1725         << "\n        from " << protonLowEnergy / keV << " keV "
1726         << " to " << protonHighEnergy / MeV << " MeV " << "." << G4endl ;
1727  G4cout << "\n        Parametrisation model for antiprotons is  "
1728         << antiprotonTable
1729         << "\n        from " << antiprotonLowEnergy / keV << " keV "
1730         << " to " << antiprotonHighEnergy / MeV << " MeV " << "." << G4endl ;
1731  if(theBarkas){
1732    G4cout << "        Parametrization of the Barkas effect is switched on."
1733           << G4endl ;
1734  }
1735  if(nStopping) {
1736    G4cout << "        Nuclear stopping power model is " << theNuclearTable
1737           << G4endl ;
1738  }
1739
1740  G4bool printHead = true;
1741
1742  const G4ProductionCutsTable* theCoupleTable=
1743    G4ProductionCutsTable::GetProductionCutsTable();
1744  size_t numOfCouples = theCoupleTable->GetTableSize();
1745
1746  // loop for materials
1747
1748  for (size_t j=0 ; j < numOfCouples; j++) {
1749
1750    const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(j);
1751    const G4Material* material= couple->GetMaterial();
1752    G4double deltaCutNow = cutForDelta[(couple->GetIndex())] ;
1753    G4double excitationEnergy = material->GetIonisation()->GetMeanExcitationEnergy();
1754
1755    if(excitationEnergy > deltaCutNow) {
1756      if(printHead) {
1757        printHead = false ;
1758
1759        G4cout << "           material       min.delta energy(keV) " << G4endl;
1760        G4cout << G4endl;
1761      }
1762
1763      G4cout << std::setw(20) << material->GetName()
1764             << std::setw(15) << excitationEnergy/keV << G4endl;
1765    }
1766  }
1767}
Note: See TracBrowser for help on using the repository browser.