source: trunk/source/processes/electromagnetic/lowenergy/src/G4hLowEnergyIonisation.cc @ 1196

Last change on this file since 1196 was 1192, checked in by garnier, 15 years ago

update par rapport a CVS
File size: 66.5 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//      GEANT 4 class implementation file
29//
30//      History: based on object model of
31//      2nd December 1995, G.Cosmo
32//      ---------- G4hLowEnergyIonisation physics process -------
33//                by Vladimir Ivanchenko, 14 July 1999
34//                was made on the base of G4hIonisation class
35//                developed by Laszlo Urban
36// ************************************************************
37// It is the extention of the ionisation process for the slow
38// charged hadrons.
39// ************************************************************
40// 28 July   1999 V.Ivanchenko cleen up
41// 17 August 1999 G.Mancinelli added ICRU parametrisations for protons
42// 20 August 1999 G.Mancinelli added ICRU tables for alpha
43// 31 August 1999 V.Ivanchenko update and cleen up
44// 30 Sept.  1999 V.Ivanchenko minor upgrade
45// 12 Dec.   1999 S. Chauvie added Barkas correction
46// 19 Jan.   2000 V.Ivanchenko minor changing in Barkas corrections
47// 02 April  2000 S. Chauvie linearization of Barkas effect
48// 03 April  2000 V.Ivanchenko Nuclear Stopping power for antiprotons
49// 23 May    2000 MG Pia  Clean up for QAO model
50// 24 May    2000 MG Pia  Code properly indented to improve legibility
51// 17 July   2000 V.Ivanchenko Bug in scaling AlongStepDoIt method
52// 25 July   2000 V.Ivanchenko New design iteration
53// 17 August 2000 V.Ivanchenko Add ion fluctuation models
54// 18 August 2000 V.Ivanchenko Bug fixed in GetConstrain
55// 22 August 2000 V.Ivanchenko Insert paramStepLimit and
56//                reorganise access to Barkas and Bloch terms 
57// 04 Sept.  2000 V.Ivanchenko rename fluctuations
58// 05 Sept.  2000 V.Ivanchenko clean up
59// 03 Oct.   2000 V.Ivanchenko CodeWizard clean up
60// 03 Nov.   2000 V.Ivanchenko MinKineticEnergy=LowestKineticEnergy=10eV
61// 05 Nov.   2000 MG Pia - Removed const cast previously introduced to get
62//                the code compiled (const G4Material* now introduced in
63//                electromagnetic/utils utils-V02-00-03 tag)
64//                (this is going back and forth, to cope with Michel's
65//                utils tag not being accepted yet by system testing)
66// 21 Nov.  2000 V.Ivanchenko Fix a problem in fluctuations
67// 23 Nov.  2000 V.Ivanchenko Ion type fluctuations only for charge>0
68// 10 May   2001 V.Ivanchenko Clean up againist Linux compilation with -Wall
69// 23 May   2001 V.Ivanchenko Minor fix in PostStepDoIt
70// 07 June  2001 V.Ivanchenko Clean up AntiProtonDEDX + add print out
71// 18 June  2001 V.Ivanchenko Cleanup print out
72// 18 Oct.  2001 V.Ivanchenko Add fluorescence
73// 30 Oct.  2001 V.Ivanchenko Add minGammaEnergy and minElectronEnergy
74// 07 Dec   2001 V.Ivanchenko Add SetFluorescence method
75// 15 Feb   2002 V.Ivanchenko Fix problem of Generic Ions
76// 25 Mar   2002 V.Ivanchenko Fix problem of fluorescence below threshold
77// 28 Mar   2002 V.Ivanchenko Set fluorescence off by default
78// 09 Apr   2002 V.Ivanchenko Fix table problem of GenericIons
79// 28 May   2002 V.Ivanchenko Remove flag fStopAndKill
80// 31 May   2002 V.Ivanchenko Add path of Fluo + Auger cuts to
81//                            AtomicDeexcitation
82// 03 Jun   2002 MGP          Restore fStopAndKill
83// 10 Jun   2002 V.Ivanchenko Restore fStopButAlive
84// 12 Jun   2002 V.Ivanchenko Fix in fluctuations - if tmax<2*Ipot Gaussian
85//                            fluctuations enables
86// 20 Sept  2002 V.Ivanchenko Clean up energy ranges for models
87// 07 Oct   2002 V.Ivanchenko Clean up initialisation of fluorescence
88// 28 Oct   2002 V.Ivanchenko Optimal binning for dE/dx
89// 10 Dec   2002 V.Ivanchenko antiProtonLowEnergy -> 25 keV, QEG model below
90// 21 Jan   2003 V.Ivanchenko Cut per region
91// 10 Mar   2003 V.Ivanchenko Use SubTypes for ions
92// 12 Apr   2003 V.Ivanchenko Cut per region for fluo AlongStep
93// 18 Apr   2003 V.Ivanchenko finalRange redefinition
94// 26 Apr   2003 V.Ivanchenko fix for stepLimit
95// 30 Mar   2004 S.Saliceti add shellCS data member and expFlag variable,
96//                          atom total cross section for the Empiric Model
97// 28 May   2004 V.Ivanchenko fix for ionisation of antiprotons in complex materials
98// 30 Aug   2004 V.Ivanchenko use energy limit for parameterisation from model
99// 03 Oct   2005 V.Ivanchenko change logic of definition of high energy limit for
100//               parametrised proton model: min(user value, model limit)
101// 26 Jan   2005 S. Chauvie added PrintInfoDefinition() for antiproton
102// 30 Sep   2009 ALF Removed dependencies to old shell Ionisation XS models
103// 11 Nov   2009 ALF Code cleaning for the Dec release
104//
105// -----------------------------------------------------------------------
106
107//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
108//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
109
110#include "G4hLowEnergyIonisation.hh"
111#include "globals.hh"
112#include "G4ios.hh"
113#include "Randomize.hh"
114#include "G4Poisson.hh"
115#include "G4UnitsTable.hh"
116#include "G4EnergyLossTables.hh"
117#include "G4Material.hh"
118#include "G4DynamicParticle.hh"
119#include "G4ParticleDefinition.hh"
120#include "G4AtomicDeexcitation.hh"
121#include "G4AtomicTransitionManager.hh"
122#include "G4ShellVacancy.hh"
123#include "G4VhShellCrossSection.hh"
124#include "G4VEMDataSet.hh"
125#include "G4EMDataSet.hh"
126#include "G4CompositeEMDataSet.hh"
127#include "G4Gamma.hh"
128#include "G4LogLogInterpolation.hh"
129#include "G4SemiLogInterpolation.hh"
130#include "G4ProcessManager.hh"
131#include "G4ProductionCutsTable.hh"
132#include "G4teoCrossSection.hh"
133#include "G4empCrossSection.hh"
134//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
135
136G4hLowEnergyIonisation::G4hLowEnergyIonisation(const G4String& processName)
137  : G4hLowEnergyLoss(processName),
138    theBetheBlochModel(0),
139    theProtonModel(0),
140    theAntiProtonModel(0),
141    theIonEffChargeModel(0),
142    theNuclearStoppingModel(0),
143    theIonChuFluctuationModel(0),
144    theIonYangFluctuationModel(0),
145    theProtonTable("ICRU_R49p"),
146    theAntiProtonTable("ICRU_R49p"),
147    theNuclearTable("ICRU_R49"),
148    nStopping(true),
149    theBarkas(true),
150    theMeanFreePathTable(0),
151    paramStepLimit (0.005),
152    shellVacancy(0),
153    shellCS(0),
154    theFluo(false)
155{ 
156  InitializeMe();
157}
158
159//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
160
161void G4hLowEnergyIonisation::InitializeMe()
162{
163  LowestKineticEnergy  = 10.0*eV ;
164  HighestKineticEnergy = 100.0*GeV ;
165  MinKineticEnergy     = 10.0*eV ; 
166  TotBin               = 360 ;
167  protonLowEnergy      = 1.*keV ;
168  protonHighEnergy     = 100.*MeV ;
169  antiProtonLowEnergy  = 25.*keV ;
170  antiProtonHighEnergy = 2.*MeV ;
171  minGammaEnergy       = 25.*keV;
172  minElectronEnergy    = 25.*keV;
173  verboseLevel         = 0;
174
175  shellCS = new G4teoCrossSection("analytical");
176
177}
178
179//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
180
181G4hLowEnergyIonisation::~G4hLowEnergyIonisation()
182{
183  if (theMeanFreePathTable) {
184    theMeanFreePathTable->clearAndDestroy();
185    delete theMeanFreePathTable;
186  }
187  if(theBetheBlochModel)delete theBetheBlochModel;
188  if(theProtonModel)delete theProtonModel;
189  if(theAntiProtonModel)delete theAntiProtonModel;
190  if(theNuclearStoppingModel)delete theNuclearStoppingModel;
191  if(theIonEffChargeModel)delete theIonEffChargeModel;
192  if(theIonChuFluctuationModel)delete theIonChuFluctuationModel;
193  if(theIonYangFluctuationModel)delete theIonYangFluctuationModel;
194  if(shellVacancy) delete shellVacancy;
195  if(shellCS) delete shellCS;
196  cutForDelta.clear();
197  G4int length = zFluoDataVector.size();
198  if(length) {
199    for(G4int i=0; i<length; i++) {
200      delete zFluoDataVector[i];
201    }
202    zFluoDataVector.clear();
203  }
204}
205
206//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
207
208void G4hLowEnergyIonisation::SetElectronicStoppingPowerModel(
209                             const G4ParticleDefinition* aParticle,
210                             const G4String& dedxTable)
211  // This method defines the ionisation parametrisation method via its name
212{
213  if(0 < aParticle->GetPDGCharge()) {
214    SetProtonElectronicStoppingPowerModel(dedxTable) ;
215  } else {
216    SetAntiProtonElectronicStoppingPowerModel(dedxTable) ;
217  }
218}
219
220//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
221
222void G4hLowEnergyIonisation::InitializeParametrisation() 
223
224{
225  // Define models for parametrisation of electronic energy losses
226  theBetheBlochModel = new G4hBetheBlochModel("Bethe-Bloch") ;
227  theProtonModel = new G4hParametrisedLossModel(theProtonTable) ;
228  protonHighEnergy = std::min(protonHighEnergy,theProtonModel->HighEnergyLimit(0, 0));
229  theAntiProtonModel = new G4QAOLowEnergyLoss(theAntiProtonTable) ;
230  theNuclearStoppingModel = new G4hNuclearStoppingModel(theNuclearTable) ;
231  theIonEffChargeModel = new G4hIonEffChargeSquare("Ziegler1988") ;
232  theIonChuFluctuationModel = new G4IonChuFluctuationModel("Chu") ;
233  theIonYangFluctuationModel = new G4IonYangFluctuationModel("Yang") ;
234}
235
236//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
237
238void G4hLowEnergyIonisation::BuildPhysicsTable(
239                       const G4ParticleDefinition& aParticleType)
240
241  //  just call BuildLossTable+BuildLambdaTable
242{
243  if(verboseLevel > 0) {
244    G4cout << "G4hLowEnergyIonisation::BuildPhysicsTable for "
245           << aParticleType.GetParticleName()
246           << " mass(MeV)= " << aParticleType.GetPDGMass()/MeV
247           << " charge= " << aParticleType.GetPDGCharge()/eplus
248           << " type= " << aParticleType.GetParticleType()
249           << G4endl;
250
251    if(verboseLevel > 1) {
252      G4ProcessVector* pv = aParticleType.GetProcessManager()->GetProcessList();
253
254      G4cout << " 0: " << (*pv)[0]->GetProcessName() << " " << (*pv)[0]
255             << " 1: " << (*pv)[1]->GetProcessName() << " " << (*pv)[1]
256        //        << " 2: " << (*pv)[2]->GetProcessName() << " " << (*pv)[2]
257             << G4endl;
258      G4cout << "ionModel= " << theIonEffChargeModel
259             << " MFPtable= " << theMeanFreePathTable
260             << " iniMass= " << initialMass
261             << G4endl;
262    }
263  }
264
265  if(aParticleType.GetParticleType() == "nucleus" &&
266     aParticleType.GetParticleName() != "GenericIon" &&
267     aParticleType.GetParticleSubType() == "generic")
268  {
269
270     G4EnergyLossTables::Register(&aParticleType,
271              theDEDXpTable,
272              theRangepTable,
273              theInverseRangepTable,
274              theLabTimepTable,
275              theProperTimepTable,
276              LowestKineticEnergy, HighestKineticEnergy,
277              proton_mass_c2/aParticleType.GetPDGMass(),
278              TotBin);
279
280     return;
281  }
282
283  if( !CutsWhereModified() && theLossTable) return;
284
285  InitializeParametrisation() ;
286  G4Proton* theProton = G4Proton::Proton();
287  G4AntiProton* theAntiProton = G4AntiProton::AntiProton();
288
289  charge = aParticleType.GetPDGCharge()/eplus;
290  chargeSquare = charge*charge ;
291
292  const G4ProductionCutsTable* theCoupleTable=
293        G4ProductionCutsTable::GetProductionCutsTable();
294  size_t numOfCouples = theCoupleTable->GetTableSize();
295
296  cutForDelta.clear();
297  cutForGamma.clear();
298
299  for (size_t j=0; j<numOfCouples; j++) {
300
301    // get material parameters needed for the energy loss calculation
302    const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(j);
303    const G4Material* material= couple->GetMaterial();
304
305    // the cut cannot be below lowest limit
306    G4double tCut = (*(theCoupleTable->GetEnergyCutsVector(1)))[j];
307    if(tCut > HighestKineticEnergy) tCut = HighestKineticEnergy;
308
309    G4double excEnergy = material->GetIonisation()->GetMeanExcitationEnergy();
310
311    tCut = std::max(tCut,excEnergy);
312    cutForDelta.push_back(tCut);
313
314    // the cut cannot be below lowest limit
315    tCut = (*(theCoupleTable->GetEnergyCutsVector(0)))[j];
316    if(tCut > HighestKineticEnergy) tCut = HighestKineticEnergy;
317    tCut = std::max(tCut,minGammaEnergy);
318    cutForGamma.push_back(tCut);
319  }
320
321  if(verboseLevel > 0) {
322    G4cout << "Cuts are defined " << G4endl;
323  }
324
325  if(0.0 < charge)
326    {
327      {
328        BuildLossTable(*theProton) ;
329
330//      The following vector has a fixed dimension (see src/G4hLowEnergyLoss.cc for more details)       
331//      It happended in the past that caused memory corruption errors. The problem is still pending, even if temporary solved
332//        G4cout << "[NOTE]: __LINE__=" << __LINE__ << ", aParticleType=" << aParticleType.GetParticleName() << ", theProton=" << theProton << ", theLossTable=" << theLossTable << ", CounterOfpProcess=" << CounterOfpProcess << G4endl;
333       
334        RecorderOfpProcess[CounterOfpProcess] = theLossTable ;
335        CounterOfpProcess++;
336      }
337  } else {
338      {
339        BuildLossTable(*theAntiProton) ;
340       
341//      The following vector has a fixed dimension (see src/G4hLowEnergyLoss.cc for more details)       
342//      It happended in the past that caused memory corruption errors. The problem is still pending, even if temporary solved
343//        G4cout << "[NOTE]: __LINE__=" << __LINE__ << ", aParticleType=" << aParticleType.GetParticleName() << ", theAntiProton=" << theAntiProton << ", theLossTable=" << theLossTable << ", CounterOfpbarProcess=" << CounterOfpbarProcess << G4endl;
344       
345        RecorderOfpbarProcess[CounterOfpbarProcess] = theLossTable ;
346        CounterOfpbarProcess++;
347      }
348  }
349
350  if(verboseLevel > 0) {
351    G4cout << "G4hLowEnergyIonisation::BuildPhysicsTable: "
352           << "Loss table is built "
353//         << theLossTable
354           << G4endl;
355  }
356
357  BuildLambdaTable(aParticleType) ;
358  BuildDataForFluorescence(aParticleType);
359
360  if(verboseLevel > 1) {
361    G4cout << (*theMeanFreePathTable) << G4endl;
362  }
363
364  if(verboseLevel > 0) {
365    G4cout << "G4hLowEnergyIonisation::BuildPhysicsTable: "
366           << "DEDX table will be built "
367//         << theDEDXpTable << " " << theDEDXpbarTable
368//         << " " << theRangepTable << " " << theRangepbarTable
369           << G4endl;
370  }
371
372  BuildDEDXTable(aParticleType) ;
373
374  if(verboseLevel > 1) {
375    G4cout << (*theDEDXpTable) << G4endl;
376  }
377
378  if((&aParticleType == theProton) ||  (&aParticleType == theAntiProton)) PrintInfoDefinition() ;
379
380  if(verboseLevel > 0) {
381    G4cout << "G4hLowEnergyIonisation::BuildPhysicsTable: end for "
382           << aParticleType.GetParticleName() << G4endl;
383  }
384}
385
386//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
387
388void G4hLowEnergyIonisation::BuildLossTable(
389                             const G4ParticleDefinition& aParticleType)
390{
391
392  // Initialisation
393  G4double lowEdgeEnergy , ionloss, ionlossBB, paramB ;
394  G4double lowEnergy, highEnergy;
395  G4Proton* theProton = G4Proton::Proton();
396
397  if(aParticleType == *theProton) {
398    lowEnergy = protonLowEnergy ;
399    highEnergy = protonHighEnergy ;
400    charge = 1.0 ;
401  } else {
402    lowEnergy = antiProtonLowEnergy ;
403    highEnergy = antiProtonHighEnergy ;
404    charge = -1.0 ;
405  }
406  chargeSquare = 1.0 ;
407
408  const G4ProductionCutsTable* theCoupleTable=
409        G4ProductionCutsTable::GetProductionCutsTable();
410  size_t numOfCouples = theCoupleTable->GetTableSize();
411
412  if ( theLossTable) {
413    theLossTable->clearAndDestroy();
414    delete theLossTable;
415  }
416
417  theLossTable = new G4PhysicsTable(numOfCouples);
418
419  //  loop for materials
420  for (size_t j=0; j<numOfCouples; j++) {
421
422    // create physics vector and fill it
423    G4PhysicsLogVector* aVector = new G4PhysicsLogVector(LowestKineticEnergy,
424                                                         HighestKineticEnergy,
425                                                         TotBin);
426
427    // get material parameters needed for the energy loss calculation
428    const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(j);
429    const G4Material* material= couple->GetMaterial();
430
431    if ( charge > 0.0 ) {
432      ionloss = ProtonParametrisedDEDX(couple,highEnergy) ;
433    } else {
434      ionloss = AntiProtonParametrisedDEDX(couple,highEnergy) ;
435    }
436
437    ionlossBB = theBetheBlochModel->TheValue(&aParticleType,material,highEnergy) ;
438    ionlossBB -= DeltaRaysEnergy(couple,highEnergy,proton_mass_c2) ;
439
440
441    paramB =  ionloss/ionlossBB - 1.0 ;
442
443    // now comes the loop for the kinetic energy values
444    for (G4int i = 0 ; i < TotBin ; i++) {
445      lowEdgeEnergy = aVector->GetLowEdgeEnergy(i) ;
446
447      // low energy part for this material, parametrised energy loss formulae
448      if ( lowEdgeEnergy < highEnergy ) {
449
450        if ( charge > 0.0 ) {
451          ionloss = ProtonParametrisedDEDX(couple,lowEdgeEnergy) ;
452        } else {
453          ionloss = AntiProtonParametrisedDEDX(couple,lowEdgeEnergy) ;
454        }
455
456      } else {
457
458        // high energy part for this material, Bethe-Bloch formula
459        ionloss = theBetheBlochModel->TheValue(theProton,material,
460                                               lowEdgeEnergy) ;
461
462        ionloss -= DeltaRaysEnergy(couple,lowEdgeEnergy,proton_mass_c2) ;
463
464        ionloss *= (1.0 + paramB*highEnergy/lowEdgeEnergy) ;
465      }
466
467      // now put the loss into the vector
468      if(verboseLevel > 1) {
469        G4cout << "E(MeV)= " << lowEdgeEnergy/MeV
470               << "  dE/dx(MeV/mm)= " << ionloss*mm/MeV
471               << " in " << material->GetName() << G4endl;
472      }
473      aVector->PutValue(i,ionloss) ;
474    }
475    // Insert vector for this material into the table
476    theLossTable->insert(aVector) ;
477  }
478}
479
480//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
481
482void G4hLowEnergyIonisation::BuildDataForFluorescence(
483                         const G4ParticleDefinition& aParticleType)
484{
485
486  if(verboseLevel > 1) {
487    G4cout << "G4hLowEnergyIonisation::BuildDataForFluorescence for "
488           << aParticleType.GetParticleName() << " is started" << G4endl;
489  }
490
491  // fill data for fluorescence
492
493  deexcitationManager.SetCutForSecondaryPhotons(minGammaEnergy);
494  deexcitationManager.SetCutForAugerElectrons(minElectronEnergy);
495
496  G4double mass = aParticleType.GetPDGMass();
497  const G4ProductionCutsTable* theCoupleTable=
498        G4ProductionCutsTable::GetProductionCutsTable();
499  size_t numOfCouples = theCoupleTable->GetTableSize();
500
501  if (shellVacancy != 0) delete shellVacancy;
502  shellVacancy = new G4ShellVacancy();
503  G4DataVector* ksi = 0;
504  G4DataVector* ksi1 = 0;
505  G4DataVector* energy = 0;
506  G4DataVector* energy1 = 0;
507  size_t binForFluo = TotBin/10;
508  G4int length = zFluoDataVector.size();
509  if(length > 0) {
510    for(G4int i=0; i<length; i++) {
511      G4VEMDataSet* x = zFluoDataVector[i];
512      delete x;
513    }
514    zFluoDataVector.clear();
515  }
516
517  G4PhysicsLogVector* bVector = new G4PhysicsLogVector(LowestKineticEnergy,
518                                                       HighestKineticEnergy,
519                                                       binForFluo);
520  const G4AtomicTransitionManager* transitionManager =
521                             G4AtomicTransitionManager::Instance();
522
523  G4double bindingEnergy;
524  //  G4double x;
525  //  G4double y;
526
527  //  loop for materials
528  for (size_t j=0; j<numOfCouples; j++) {
529
530    // get material parameters needed for the energy loss calculation
531    const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(j);
532    const G4Material* material= couple->GetMaterial();
533
534    const G4ElementVector* theElementVector = material->GetElementVector();
535    size_t NumberOfElements = material->GetNumberOfElements() ;
536    const G4double* theAtomicNumDensityVector =
537                    material->GetAtomicNumDensityVector();
538    G4VDataSetAlgorithm* interp = new G4SemiLogInterpolation();
539    G4VEMDataSet* xsis = new G4CompositeEMDataSet(interp, 1., 1.);
540    G4VDataSetAlgorithm* interp1 = new G4SemiLogInterpolation();
541    G4VEMDataSet* xsis1 = new G4CompositeEMDataSet(interp1, 1., 1.);
542
543    G4double tCut = cutForDelta[j];
544    G4double elDensity = 1.;
545
546    for (size_t iel=0; iel<NumberOfElements; iel++ ) {
547
548      G4int Z = (G4int)((*theElementVector)[iel]->GetZ());
549      G4int nShells = transitionManager->NumberOfShells(Z);
550      energy = new G4DataVector();
551      ksi    = new G4DataVector();
552      energy1= new G4DataVector();
553      ksi1   = new G4DataVector();
554      //if(NumberOfElements > 1)
555      elDensity = theAtomicNumDensityVector[iel]/((G4double)nShells);
556
557      for (size_t j = 0; j<binForFluo; j++) {
558
559        G4double tkin  = bVector->GetLowEdgeEnergy(j);
560        G4double gamma = tkin/mass + 1.;
561        G4double beta2 = 1.0 - 1.0/(gamma*gamma);
562        G4double r     = electron_mass_c2/mass;
563        G4double tmax  = 2.*electron_mass_c2*(gamma*gamma - 1.)/(1. + 2.*gamma*r + r*r);
564        G4double cross   = 0.;
565        G4double cross1  = 0.;
566        G4double eAverage= 0.;
567        G4double tmin = std::min(tCut,tmax);
568        G4double rel;
569
570        for (G4int n=0; n<nShells; n++) {
571
572          bindingEnergy = transitionManager->Shell(Z, n)->BindingEnergy();
573          if (tmin > bindingEnergy) {
574            rel = std::log(tmin/bindingEnergy);
575            eAverage   += rel - beta2*(tmin - bindingEnergy)/tmax;
576            cross      += 1.0/bindingEnergy - 1.0/tmin - beta2*rel/tmax;
577          }
578          if (tmax > tmin) {
579            cross1     += 1.0/tmin - 1.0/tmax - beta2*std::log(tmax/tmin)/tmax;
580          }
581        }
582
583        cross1 *= elDensity;
584        energy1->push_back(tkin);
585        ksi1->push_back(cross1);
586
587        if(eAverage > 0.) cross /= eAverage;
588        else              cross  = 0.;
589
590        energy->push_back(tkin);
591        ksi->push_back(cross);
592      }
593      G4VDataSetAlgorithm* algo = interp->Clone();
594      G4VEMDataSet* set = new G4EMDataSet(Z,energy,ksi,algo,1.,1.);
595      xsis->AddComponent(set);
596      G4VDataSetAlgorithm* algo1 = interp1->Clone();
597      G4VEMDataSet* set1 = new G4EMDataSet(Z,energy1,ksi1,algo1,1.,1.);
598      xsis1->AddComponent(set1);
599    }
600    if(verboseLevel > 1) {
601      G4cout << "### Shell inverse cross sections for "
602             << material->GetName() << G4endl;
603      xsis->PrintData();
604      G4cout << "### Atom cross sections for "
605             << material->GetName() << G4endl;
606      xsis1->PrintData();
607    }
608    shellVacancy->AddXsiTable(xsis);
609    zFluoDataVector.push_back(xsis1);
610  }
611  delete bVector;
612}
613
614//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
615
616void G4hLowEnergyIonisation::BuildLambdaTable(
617                       const G4ParticleDefinition& aParticleType)
618
619{
620  // Build mean free path tables for the delta ray production process
621  //     tables are built for MATERIALS
622
623  if(verboseLevel > 1) {
624    G4cout << "G4hLowEnergyIonisation::BuildLambdaTable for "
625           << aParticleType.GetParticleName() << " is started" << G4endl;
626  }
627
628
629  G4double lowEdgeEnergy, value;
630  charge = aParticleType.GetPDGCharge()/eplus ;
631  chargeSquare = charge*charge ;
632  initialMass = aParticleType.GetPDGMass();
633
634  const G4ProductionCutsTable* theCoupleTable=
635        G4ProductionCutsTable::GetProductionCutsTable();
636  size_t numOfCouples = theCoupleTable->GetTableSize();
637
638
639  if (theMeanFreePathTable) {
640    theMeanFreePathTable->clearAndDestroy();
641    delete theMeanFreePathTable;
642  }
643
644  theMeanFreePathTable = new G4PhysicsTable(numOfCouples);
645
646  // loop for materials
647
648  for (size_t J=0 ; J < numOfCouples; J++) {
649
650    //create physics vector then fill it ....
651    G4PhysicsLogVector* aVector = new G4PhysicsLogVector(LowestKineticEnergy,
652                                                         HighestKineticEnergy,
653                                                         TotBin);
654
655    // compute the (macroscopic) cross section first
656    const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(J);
657    const G4Material* material= couple->GetMaterial();
658
659    const G4ElementVector* theElementVector =
660                           material->GetElementVector() ;
661    const G4double* theAtomicNumDensityVector =
662                           material->GetAtomicNumDensityVector();
663    const G4int NumberOfElements = material->GetNumberOfElements() ;
664
665      // get the electron kinetic energy cut for the actual material,
666      //  it will be used in ComputeMicroscopicCrossSection
667      // ( it is the SAME for ALL the ELEMENTS in THIS MATERIAL )
668      //   ------------------------------------------------------
669
670    G4double deltaCut = cutForDelta[J];
671
672    for ( G4int i = 0 ; i < TotBin ; i++ ) {
673      lowEdgeEnergy = aVector->GetLowEdgeEnergy(i) ;
674      G4double sigma = 0.0 ;
675      G4int Z;
676     
677      for (G4int iel=0; iel<NumberOfElements; iel++ ) {
678        Z = (G4int) (*theElementVector)[iel]->GetZ();
679        totalCrossSectionMap [Z] = ComputeMicroscopicCrossSection(
680                                                                  aParticleType,
681                                                                  lowEdgeEnergy,
682                                                                  Z,
683                                                                  deltaCut ) ; 
684        sigma += theAtomicNumDensityVector[iel]*ComputeMicroscopicCrossSection(
685                                                                               aParticleType,
686                                                                               lowEdgeEnergy,
687                                                                               Z,
688                                                                               deltaCut ) ; 
689       
690      }
691
692      // mean free path = 1./macroscopic cross section
693
694      value = sigma<=0 ? DBL_MAX : 1./sigma ;
695
696      aVector->PutValue(i, value) ;
697    }
698
699    theMeanFreePathTable->insert(aVector);
700  }
701
702}
703
704
705//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
706
707G4double G4hLowEnergyIonisation::ComputeMicroscopicCrossSection(
708                           const G4ParticleDefinition& aParticleType,
709                                 G4double kineticEnergy,
710                                 G4double atomicNumber,
711                                 G4double deltaCutInEnergy) const
712{
713  //******************************************************************
714  // cross section formula is OK for spin=0, 1/2, 1 only !
715  // *****************************************************************
716
717  // calculates the microscopic cross section in GEANT4 internal units
718  //    ( it is called for elements , AtomicNumber = z )
719
720  G4double energy, gamma, beta2, tmax, var;
721  G4double totalCrossSection = 0.0 ;
722
723  G4double particleMass = initialMass;
724
725  // get particle data ...................................
726
727  energy = kineticEnergy + particleMass;
728
729  // some kinematics......................
730
731  gamma  = energy/particleMass;
732  beta2  = 1.0 - 1.0/(gamma*gamma);
733  var    = electron_mass_c2/particleMass;
734  tmax   = 2.*electron_mass_c2*(gamma*gamma - 1.)/(1. + 2.*gamma*var + var*var);
735
736  // now you can calculate the total cross section
737
738  if( tmax > deltaCutInEnergy ) {
739
740    var=deltaCutInEnergy/tmax;
741    totalCrossSection = (1.0 - var*(1.0 - beta2*std::log(var))) / deltaCutInEnergy ;
742    G4double spin = aParticleType.GetPDGSpin() ;
743
744    // +term for spin=1/2 particle
745    if( 0.5 == spin )
746      totalCrossSection +=  0.5 * (tmax - deltaCutInEnergy) / (energy*energy);
747
748    // +term for spin=1 particle
749    else if( 0.9 < spin )
750      totalCrossSection += -std::log(var)/(3.0*deltaCutInEnergy) +
751        (tmax - deltaCutInEnergy) * ( (5.0+ 1.0/var)*0.25 / (energy*energy) -
752         beta2 / (tmax * deltaCutInEnergy) ) / 3.0 ;
753
754    totalCrossSection *= twopi_mc2_rcl2 * atomicNumber / beta2 ;
755  }
756
757  return totalCrossSection ;
758}
759
760//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
761
762G4double G4hLowEnergyIonisation::GetMeanFreePath(const G4Track& trackData,
763                                                 G4double, // previousStepSize
764                                                 enum G4ForceCondition* condition)
765{
766   const G4DynamicParticle* aParticle = trackData.GetDynamicParticle();
767   const G4MaterialCutsCouple* couple = trackData.GetMaterialCutsCouple();
768   const G4Material* material = couple->GetMaterial();
769   G4double meanFreePath;
770   G4bool isOutRange ;
771
772   *condition = NotForced ;
773
774   G4double kineticEnergy = (aParticle->GetKineticEnergy())*initialMass/(aParticle->GetMass());
775   charge = aParticle->GetCharge()/eplus;
776   chargeSquare = theIonEffChargeModel->TheValue(aParticle, material);
777
778   if(kineticEnergy < LowestKineticEnergy) meanFreePath = DBL_MAX;
779
780   else {
781     if(kineticEnergy > HighestKineticEnergy)
782                    kineticEnergy = HighestKineticEnergy;
783     meanFreePath = (((*theMeanFreePathTable)(couple->GetIndex()))->
784                    GetValue(kineticEnergy,isOutRange))/chargeSquare;
785     }
786
787   return meanFreePath ;
788}
789
790//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
791
792G4double G4hLowEnergyIonisation::GetConstraints(
793                                 const G4DynamicParticle* particle,
794                                 const G4MaterialCutsCouple* couple)
795{
796  // returns the Step limit
797  // dEdx is calculated as well as the range
798  // based on Effective Charge Approach
799
800  const G4Material* material = couple->GetMaterial();
801  G4Proton* theProton = G4Proton::Proton();
802  G4AntiProton* theAntiProton = G4AntiProton::AntiProton();
803
804  G4double stepLimit = 0.0 ;
805  G4double dx, highEnergy;
806
807  G4double massRatio = proton_mass_c2/(particle->GetMass()) ;
808  G4double kineticEnergy = particle->GetKineticEnergy() ;
809
810  // Scale the kinetic energy
811
812  G4double tscaled = kineticEnergy*massRatio ;
813  fBarkas = 0.0;
814
815  if(charge > 0.0) {
816
817    highEnergy = protonHighEnergy ;
818
819    fRangeNow = G4EnergyLossTables::GetRange(theProton, tscaled, couple);
820    dx = G4EnergyLossTables::GetRange(theProton, highEnergy, couple);
821    fdEdx = G4EnergyLossTables::GetDEDX(theProton, tscaled, couple)
822          * chargeSquare ;
823
824        // Correction for positive ions
825    if(theBarkas && tscaled > highEnergy) { 
826        fBarkas = BarkasTerm(material,tscaled)*std::sqrt(chargeSquare)*chargeSquare
827                + BlochTerm(material,tscaled,chargeSquare);
828    }
829    // Antiprotons and negative hadrons
830  } else {
831
832    highEnergy = antiProtonHighEnergy ;
833    fRangeNow = G4EnergyLossTables::GetRange(theAntiProton, tscaled, couple);
834    dx = G4EnergyLossTables::GetRange(theAntiProton, highEnergy, couple);
835    fdEdx = G4EnergyLossTables::GetDEDX(theAntiProton, tscaled, couple)
836          * chargeSquare ;
837
838    if(theBarkas && tscaled > highEnergy) { 
839        fBarkas = -BarkasTerm(material,tscaled)*std::sqrt(chargeSquare)*chargeSquare
840                + BlochTerm(material,tscaled,chargeSquare);
841    }
842  }
843  /*
844  const G4Material* mat = couple->GetMaterial();
845  G4double fac = gram/(MeV*cm2*mat->GetDensity());
846  G4cout << particle->GetDefinition()->GetParticleName()
847         << " in " << mat->GetName()
848         << " E(MeV)= " << kineticEnergy/MeV
849         << " dedx(MeV*cm^2/g)= " << fdEdx*fac
850         << " barcas(MeV*cm^2/gram)= " << fBarkas*fac
851         << " Q^2= " << chargeSquare
852         << G4endl;
853  */
854  // scaling back
855  fRangeNow /= (chargeSquare*massRatio) ;
856  dx        /= (chargeSquare*massRatio) ;
857
858  stepLimit  = fRangeNow ;
859  G4double r = std::min(finalRange, couple->GetProductionCuts()
860                 ->GetProductionCut(idxG4ElectronCut));
861
862  if (fRangeNow > r) {
863    stepLimit = dRoverRange*fRangeNow + r*(1.0 - dRoverRange)*(2.0 - r/fRangeNow);
864    if (stepLimit > fRangeNow) stepLimit = fRangeNow;
865  }
866  // compute the (random) Step limit in standard energy range
867  if(tscaled > highEnergy ) {
868
869    // add Barkas correction directly to dedx
870    fdEdx  += fBarkas;
871 
872    if(stepLimit > fRangeNow - dx*0.9) stepLimit = fRangeNow - dx*0.9 ;
873
874  // Step limit in low energy range
875  } else {
876    G4double x = dx*paramStepLimit;
877    if (stepLimit > x) stepLimit = x;
878  }
879  return stepLimit ;
880}
881
882//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
883
884G4VParticleChange* G4hLowEnergyIonisation::AlongStepDoIt(
885                                           const G4Track& trackData,
886                                           const G4Step& stepData)
887{
888  // compute the energy loss after a step
889  G4Proton* theProton = G4Proton::Proton();
890  G4AntiProton* theAntiProton = G4AntiProton::AntiProton();
891  G4double finalT = 0.0 ;
892
893  aParticleChange.Initialize(trackData) ;
894
895  const G4MaterialCutsCouple* couple = trackData.GetMaterialCutsCouple();
896  const G4Material* material = couple->GetMaterial();
897
898  // get the actual (true) Step length from stepData
899  const G4double step = stepData.GetStepLength() ;
900
901  const G4DynamicParticle* particle = trackData.GetDynamicParticle() ;
902
903  G4double kineticEnergy = particle->GetKineticEnergy() ;
904  G4double massRatio = proton_mass_c2/(particle->GetMass()) ;
905  G4double tscaled= kineticEnergy*massRatio ;
906  G4double eloss = 0.0 ;
907  G4double nloss = 0.0 ;
908
909
910    // very small particle energy
911  if(kineticEnergy < MinKineticEnergy) {
912
913    eloss = kineticEnergy ;
914
915    // particle energy outside tabulated energy range
916  } else if( kineticEnergy > HighestKineticEnergy) {
917    eloss = step*fdEdx ;
918
919    // big step
920  } else if(step >= fRangeNow ) {
921    eloss = kineticEnergy ;
922
923    // tabulated range
924  } else {
925
926    // step longer than linear step limit
927    if(step > linLossLimit*fRangeNow) {
928
929      G4double rscaled= fRangeNow*massRatio*chargeSquare ;
930      G4double sscaled=   step   *massRatio*chargeSquare ;
931
932      if(charge > 0.0) {
933        eloss = G4EnergyLossTables::GetPreciseEnergyFromRange(
934                                    theProton,rscaled, couple) -
935                G4EnergyLossTables::GetPreciseEnergyFromRange(
936                                    theProton,rscaled-sscaled,couple) ;
937
938      } else {
939        eloss = G4EnergyLossTables::GetPreciseEnergyFromRange(
940                                    theAntiProton,rscaled,couple) -
941                G4EnergyLossTables::GetPreciseEnergyFromRange(
942                                    theAntiProton,rscaled-sscaled,couple) ;
943      }
944      eloss /= massRatio ;
945
946      // Barkas correction at big step     
947      eloss += fBarkas*step;
948
949    // step shorter than linear step limit
950    } else {
951      eloss = step*fdEdx ;
952    }
953    if(nStopping && tscaled < protonHighEnergy) {
954      nloss = (theNuclearStoppingModel->TheValue(particle, material))*step;
955    }
956  }
957
958  if(eloss < 0.0) eloss = 0.0;
959
960  finalT = kineticEnergy - eloss - nloss;
961
962  if( EnlossFlucFlag && 0.0 < eloss && finalT > MinKineticEnergy) {
963
964    //  now the electron loss with fluctuation
965    eloss = ElectronicLossFluctuation(particle, couple, eloss, step) ;
966    if(eloss < 0.0) eloss = 0.0;
967    finalT = kineticEnergy - eloss - nloss;
968  }
969
970  //  stop particle if the kinetic energy <= MinKineticEnergy
971  if (finalT*massRatio <= MinKineticEnergy ) {
972
973     finalT = 0.0;
974      if(!particle->GetDefinition()->GetProcessManager()->
975                     GetAtRestProcessVector()->size())
976        aParticleChange.ProposeTrackStatus(fStopAndKill);
977      else
978        aParticleChange.ProposeTrackStatus(fStopButAlive);
979  }
980
981  aParticleChange.ProposeEnergy( finalT );
982  eloss = kineticEnergy-finalT;
983
984  // Deexcitation only of ionised atoms
985  G4double hMass = particle->GetMass();
986  std::vector<G4DynamicParticle*>* newpart = 0;
987  G4DynamicParticle* part = 0;
988
989  if(theFluo) newpart = DeexciteAtom(couple, kineticEnergy, hMass, eloss);
990
991  if(newpart != 0) {
992
993    size_t nSecondaries = newpart->size();
994    aParticleChange.SetNumberOfSecondaries(nSecondaries);
995    G4Track* newtrack = 0;
996    const G4StepPoint* preStep = stepData.GetPreStepPoint();
997    const G4StepPoint* postStep = stepData.GetPostStepPoint();
998    G4ThreeVector r = preStep->GetPosition();
999    G4ThreeVector deltaR = postStep->GetPosition();
1000    deltaR -= r;
1001    G4double t = preStep->GetGlobalTime();
1002    G4double deltaT = postStep->GetGlobalTime();
1003    deltaT -= t;
1004    G4double time, q, e;
1005    G4ThreeVector position;
1006
1007    for(size_t i=0; i<nSecondaries; i++) {
1008
1009      part = (*newpart)[i];
1010      if(part) {
1011
1012        e = part->GetKineticEnergy();
1013        if(e <= eloss) {
1014
1015          eloss -= e;
1016          q = G4UniformRand();
1017          time = deltaT*q + t;
1018          position  = deltaR*q;
1019          position += r;
1020          newtrack = new G4Track(part, time, position);
1021          aParticleChange.AddSecondary(newtrack);
1022
1023        } else {
1024
1025          delete part;
1026
1027        }
1028      }
1029    }
1030    delete newpart;
1031  }
1032
1033  aParticleChange.ProposeLocalEnergyDeposit(eloss);
1034  return &aParticleChange ;
1035}
1036
1037//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1038
1039G4double G4hLowEnergyIonisation::ProtonParametrisedDEDX(
1040                                 const G4MaterialCutsCouple* couple,
1041                                       G4double kineticEnergy) const
1042{
1043  const G4Material* material = couple->GetMaterial();
1044  G4Proton* theProton = G4Proton::Proton();
1045  G4double eloss = 0.0;
1046
1047    // Free Electron Gas Model
1048  if(kineticEnergy < protonLowEnergy) {
1049    eloss = (theProtonModel->TheValue(theProton, material, protonLowEnergy))
1050          * std::sqrt(kineticEnergy/protonLowEnergy) ;
1051
1052    // Parametrisation
1053  } else {
1054    eloss = theProtonModel->TheValue(theProton, material, kineticEnergy) ;
1055  }
1056
1057  // Delta rays energy
1058  eloss -= DeltaRaysEnergy(couple,kineticEnergy,proton_mass_c2) ;
1059
1060  if(verboseLevel > 2) {
1061    G4cout << "p E(MeV)= " << kineticEnergy/MeV
1062           << " dE/dx(MeV/mm)= " << eloss*mm/MeV
1063           << " for " << material->GetName()
1064           << " model: " << theProtonModel << G4endl;
1065  }
1066
1067  if(eloss < 0.0) eloss = 0.0 ;
1068
1069  return eloss ;
1070}
1071
1072//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1073
1074G4double G4hLowEnergyIonisation::AntiProtonParametrisedDEDX(
1075                                 const G4MaterialCutsCouple* couple,
1076                                       G4double kineticEnergy) const
1077{
1078  const G4Material* material = couple->GetMaterial();
1079  G4AntiProton* theAntiProton = G4AntiProton::AntiProton();
1080  G4double eloss = 0.0 ;
1081
1082  // Antiproton model is used
1083  if(theAntiProtonModel->IsInCharge(theAntiProton,material)) {
1084    if(kineticEnergy < antiProtonLowEnergy) {
1085      eloss = theAntiProtonModel->TheValue(theAntiProton,material,antiProtonLowEnergy)
1086            * std::sqrt(kineticEnergy/antiProtonLowEnergy) ;
1087
1088    // Parametrisation
1089    } else {
1090      eloss = theAntiProtonModel->TheValue(theAntiProton,material,
1091                                           kineticEnergy);
1092    }
1093
1094  // The proton model is used + Barkas correction
1095  } else {
1096    if(kineticEnergy < protonLowEnergy) {
1097      eloss = theProtonModel->TheValue(G4Proton::Proton(),material,protonLowEnergy)
1098          * std::sqrt(kineticEnergy/protonLowEnergy) ;
1099
1100    // Parametrisation
1101    } else {
1102      eloss = theProtonModel->TheValue(G4Proton::Proton(),material,
1103                                       kineticEnergy);
1104    }
1105    //if(theBarkas) eloss -= 2.0*BarkasTerm(material, kineticEnergy);
1106  }
1107
1108  // Delta rays energy
1109  eloss -= DeltaRaysEnergy(couple,kineticEnergy,proton_mass_c2) ;
1110
1111  if(verboseLevel > 2) {
1112    G4cout << "pbar E(MeV)= " << kineticEnergy/MeV
1113           << " dE/dx(MeV/mm)= " << eloss*mm/MeV
1114           << " for " << material->GetName()
1115           << " model: " << theProtonModel << G4endl;
1116  }
1117
1118  if(eloss < 0.0) eloss = 0.0 ;
1119
1120  return eloss ;
1121}
1122
1123//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1124
1125G4double G4hLowEnergyIonisation::DeltaRaysEnergy(
1126                                 const G4MaterialCutsCouple* couple,
1127                                       G4double kineticEnergy,
1128                                       G4double particleMass) const
1129{
1130  G4double dloss = 0.0 ;
1131
1132  G4double deltaCutNow = cutForDelta[(couple->GetIndex())] ;
1133  const G4Material* material = couple->GetMaterial();
1134  G4double electronDensity = material->GetElectronDensity();
1135  G4double eexc = material->GetIonisation()->GetMeanExcitationEnergy();
1136
1137  G4double tau = kineticEnergy/particleMass ;
1138  G4double rateMass = electron_mass_c2/particleMass ;
1139
1140  // some local variables
1141
1142  G4double gamma,bg2,beta2,tmax,x ;
1143
1144  gamma = tau + 1.0 ;
1145  bg2 = tau*(tau+2.0) ;
1146  beta2 = bg2/(gamma*gamma) ;
1147  tmax = 2.*electron_mass_c2*bg2/(1.0+2.0*gamma*rateMass+rateMass*rateMass) ;
1148
1149  // Validity range for delta electron cross section
1150  G4double deltaCut = std::max(deltaCutNow, eexc);
1151
1152  if ( deltaCut < tmax) {
1153    x = deltaCut / tmax ;
1154    dloss = ( beta2 * (x - 1.0) - std::log(x) ) * twopi_mc2_rcl2
1155          * electronDensity / beta2 ;
1156  }
1157  return dloss ;
1158}
1159
1160//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1161
1162G4VParticleChange* G4hLowEnergyIonisation::PostStepDoIt(
1163                                           const G4Track& trackData,
1164                                           const G4Step& stepData)
1165{
1166  // Units are expressed in GEANT4 internal units.
1167
1168  G4double KineticEnergy,TotalEnergy,TotalMomentum,betasquare,
1169           DeltaKineticEnergy,DeltaTotalMomentum,costheta,sintheta,phi,
1170           dirx,diry,dirz,finalKineticEnergy,finalPx,finalPy,finalPz,
1171           x,xc,grej,Psquare,Esquare,rate,finalMomentum ;
1172
1173  aParticleChange.Initialize(trackData) ;
1174  const G4MaterialCutsCouple* couple = trackData.GetMaterialCutsCouple();
1175
1176  const G4DynamicParticle* aParticle = trackData.GetDynamicParticle() ;
1177
1178  // some kinematics
1179
1180  ParticleMass=aParticle->GetDefinition()->GetPDGMass();
1181  KineticEnergy=aParticle->GetKineticEnergy();
1182  TotalEnergy=KineticEnergy + ParticleMass ;
1183  Psquare=KineticEnergy*(TotalEnergy+ParticleMass) ;
1184  Esquare=TotalEnergy*TotalEnergy;
1185  betasquare=Psquare/Esquare;
1186  G4ThreeVector ParticleDirection = aParticle->GetMomentumDirection() ;
1187
1188  G4double gamma= KineticEnergy/ParticleMass + 1.;
1189  G4double r    = electron_mass_c2/ParticleMass;
1190  G4double tmax = 2.*electron_mass_c2*(gamma*gamma - 1.)/(1. + 2.*gamma*r + r*r);
1191
1192  // Validity range for delta electron cross section
1193  G4double DeltaCut = cutForDelta[couple->GetIndex()];
1194
1195  // This should not be a case
1196  if(DeltaCut >= tmax)
1197       return G4VContinuousDiscreteProcess::PostStepDoIt(trackData,stepData);
1198
1199  xc   = DeltaCut / tmax;
1200  rate = tmax / TotalEnergy;
1201  rate = rate*rate ;
1202  G4double spin = aParticle->GetDefinition()->GetPDGSpin() ;
1203
1204  // sampling follows ...
1205      do {
1206        x=xc/(1.-(1.-xc)*G4UniformRand());
1207
1208        if(0.0 == spin) {
1209          grej = 1.0 - betasquare * x ;
1210
1211        } else if (0.5 == spin) {
1212          grej = (1.0 - betasquare * x + 0.5*x*x*rate) / (1.0 + 0.5 * rate) ;
1213
1214        } else {
1215          grej = (1.0 - betasquare * x ) * (1.0 + x/ (3.0*xc)) +
1216            x * x * rate * (1.0 + 0.5 * x / xc) / 3.0 /
1217            (1.0 + 1.0/(3.0*xc) + rate *(1.0+ 0.5/xc) /3.0) ;
1218        }
1219
1220      } while( G4UniformRand() > grej );
1221
1222
1223  DeltaKineticEnergy = x * tmax;
1224
1225  DeltaTotalMomentum = std::sqrt(DeltaKineticEnergy * (DeltaKineticEnergy +
1226                                                  2. * electron_mass_c2 )) ;
1227  TotalMomentum = std::sqrt(Psquare) ;
1228  costheta = DeltaKineticEnergy * (TotalEnergy + electron_mass_c2)
1229    /(DeltaTotalMomentum * TotalMomentum) ;
1230
1231  //  protection against costheta > 1 or < -1   ---------------
1232  if ( costheta < -1. )
1233    costheta = -1. ;
1234  if ( costheta > +1. )
1235    costheta = +1. ;
1236
1237  //  direction of the delta electron  ........
1238  phi = twopi * G4UniformRand() ;
1239  sintheta = std::sqrt(1. - costheta*costheta);
1240  dirx = sintheta * std::cos(phi) ;
1241  diry = sintheta * std::sin(phi) ;
1242  dirz = costheta ;
1243
1244  G4ThreeVector DeltaDirection(dirx,diry,dirz) ;
1245  DeltaDirection.rotateUz(ParticleDirection) ;
1246
1247  // create G4DynamicParticle object for delta ray
1248  G4DynamicParticle *theDeltaRay = new G4DynamicParticle;
1249  theDeltaRay->SetKineticEnergy( DeltaKineticEnergy );
1250  theDeltaRay->SetMomentumDirection(DeltaDirection.x(),
1251                                    DeltaDirection.y(),
1252                                    DeltaDirection.z());
1253  theDeltaRay->SetDefinition(G4Electron::Electron());
1254
1255  // fill aParticleChange
1256  finalKineticEnergy = KineticEnergy - DeltaKineticEnergy ;
1257
1258  // Generation of Fluorescence and Auger
1259  size_t nSecondaries = 0;
1260  size_t totalNumber  = 1;
1261  std::vector<G4DynamicParticle*>* secondaryVector = 0;
1262  G4DynamicParticle* aSecondary = 0;
1263  G4ParticleDefinition* type = 0;
1264
1265  // Select atom and shell
1266  G4int Z = SelectRandomAtom(couple, KineticEnergy);
1267
1268  //   G4cout << "Fluorescence is switched :" << theFluo << G4endl;
1269
1270  // Fluorescence data start from element 6
1271  if(theFluo && Z > 5) {
1272
1273
1274
1275    // Atom total cross section     
1276    shellCS->SetTotalCS(totalCrossSectionMap[Z]);   
1277
1278    G4int shell = shellCS->SelectRandomShell(Z, KineticEnergy,ParticleMass,DeltaKineticEnergy);
1279
1280    if (shell!=-1) {       
1281     
1282      const G4AtomicShell* atomicShell =
1283        (G4AtomicTransitionManager::Instance())->Shell(Z, shell);
1284      G4double bindingEnergy = atomicShell->BindingEnergy();
1285     
1286      if(verboseLevel > 1) {
1287        G4cout << "PostStep Z= " << Z << " shell= " << shell
1288               << " bindingE(keV)= " << bindingEnergy/keV
1289               << " finalE(keV)= " << finalKineticEnergy/keV
1290               << G4endl;
1291      }
1292     
1293     
1294     
1295      if (finalKineticEnergy >= bindingEnergy
1296          && (bindingEnergy >= minGammaEnergy
1297              ||  bindingEnergy >= minElectronEnergy) ) {
1298       
1299        G4int shellId = atomicShell->ShellId();
1300        secondaryVector = deexcitationManager.GenerateParticles(Z, shellId);
1301       
1302        if (secondaryVector != 0) {
1303         
1304          nSecondaries = secondaryVector->size();
1305          for (size_t i = 0; i<nSecondaries; i++) {
1306           
1307            aSecondary = (*secondaryVector)[i];
1308            if (aSecondary) {
1309             
1310              G4double e = aSecondary->GetKineticEnergy();
1311              type = aSecondary->GetDefinition();
1312              if (e < finalKineticEnergy &&
1313                  ((type == G4Gamma::Gamma() && e > minGammaEnergy ) ||
1314                   (type == G4Electron::Electron() && e > minElectronEnergy ))) {
1315               
1316                finalKineticEnergy -= e;
1317                totalNumber++;
1318               
1319              } else {
1320               
1321                delete aSecondary;
1322                (*secondaryVector)[i] = 0;
1323              }
1324            }
1325          }
1326        }
1327      }
1328    }
1329  }
1330 
1331  // Save delta-electrons
1332
1333  G4double edep = 0.0;
1334
1335  if (finalKineticEnergy > MinKineticEnergy)
1336    {
1337      finalPx = TotalMomentum*ParticleDirection.x()
1338        - DeltaTotalMomentum*DeltaDirection.x();
1339      finalPy = TotalMomentum*ParticleDirection.y()
1340        - DeltaTotalMomentum*DeltaDirection.y();
1341      finalPz = TotalMomentum*ParticleDirection.z()
1342        - DeltaTotalMomentum*DeltaDirection.z();
1343      finalMomentum =
1344        std::sqrt(finalPx*finalPx+finalPy*finalPy+finalPz*finalPz) ;
1345      finalPx /= finalMomentum ;
1346      finalPy /= finalMomentum ;
1347      finalPz /= finalMomentum ;
1348
1349      aParticleChange.ProposeMomentumDirection( finalPx,finalPy,finalPz );
1350    }
1351  else
1352    {
1353      edep = finalKineticEnergy;
1354      finalKineticEnergy = 0.;
1355      aParticleChange.ProposeMomentumDirection(ParticleDirection.x(),
1356                      ParticleDirection.y(),ParticleDirection.z());
1357      if(!aParticle->GetDefinition()->GetProcessManager()->
1358                     GetAtRestProcessVector()->size())
1359        aParticleChange.ProposeTrackStatus(fStopAndKill);
1360      else
1361        aParticleChange.ProposeTrackStatus(fStopButAlive);
1362    }
1363
1364  aParticleChange.ProposeEnergy( finalKineticEnergy );
1365  aParticleChange.ProposeLocalEnergyDeposit (edep);
1366  aParticleChange.SetNumberOfSecondaries(totalNumber);
1367  aParticleChange.AddSecondary(theDeltaRay);
1368
1369  // Save Fluorescence and Auger
1370
1371  if (secondaryVector) {
1372
1373    for (size_t l = 0; l < nSecondaries; l++) {
1374
1375      aSecondary = (*secondaryVector)[l];
1376      if(aSecondary) {
1377        aParticleChange.AddSecondary(aSecondary);
1378      }
1379    }
1380    delete secondaryVector;
1381  }
1382
1383  return G4VContinuousDiscreteProcess::PostStepDoIt(trackData,stepData);
1384}
1385
1386//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1387
1388
1389
1390void G4hLowEnergyIonisation::SelectShellIonisationCS(G4String val) {
1391
1392  if (shellCS) delete shellCS;
1393
1394  if (val == "empirical") {
1395    shellCS = new G4empCrossSection();
1396  }
1397 else {
1398    shellCS = new G4teoCrossSection(val);
1399  }
1400}
1401
1402
1403
1404//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1405
1406
1407std::vector<G4DynamicParticle*>* G4hLowEnergyIonisation::DeexciteAtom(const G4MaterialCutsCouple* couple,
1408                                           G4double incidentEnergy,
1409                                           G4double hMass,
1410                                           G4double eLoss)
1411{
1412
1413  if (verboseLevel > 1) {
1414        G4cout << "DeexciteAtom: cutForPhotons(keV)= " << minGammaEnergy/keV
1415               << "  cutForElectrons(keV)= " << minElectronEnergy/keV
1416               << "  eLoss(MeV)= " << eLoss
1417               << G4endl;
1418  }
1419
1420
1421
1422  if(eLoss < minGammaEnergy && eLoss < minElectronEnergy) return 0;
1423
1424  const G4Material* material = couple->GetMaterial();
1425  G4int index    = couple->GetIndex();
1426  //  G4double eexc  = material->GetIonisation()->GetMeanExcitationEnergy();
1427  G4double gamma = incidentEnergy/hMass + 1;
1428  G4double beta2 = 1.0 - 1.0/(gamma*gamma);
1429  G4double r     = electron_mass_c2/hMass;
1430  G4double tmax  = 2.*electron_mass_c2*(gamma*gamma - 1.)/(1. + 2.*gamma*r + r*r);
1431  G4double tcut  = std::min(tmax,cutForDelta[index]);
1432  const G4AtomicTransitionManager* transitionManager =
1433                             G4AtomicTransitionManager::Instance();
1434
1435  size_t nElements = material->GetNumberOfElements();
1436  const G4ElementVector* theElementVector = material->GetElementVector();
1437  G4bool stop = true;
1438
1439  for (size_t j=0; j<nElements; j++) {
1440
1441    G4int Z = (G4int)((*theElementVector)[j]->GetZ());
1442    G4double maxE = transitionManager->Shell(Z, 0)->BindingEnergy();
1443
1444    if (Z > 5 && maxE < tcut && (maxE > minGammaEnergy || maxE > minElectronEnergy) ) {
1445      stop = false;
1446      break;
1447    }
1448  }
1449
1450  if(stop) return 0;
1451
1452  // create vector of tracks of secondary particles
1453
1454  std::vector<G4DynamicParticle*>* partVector =
1455         new std::vector<G4DynamicParticle*>;
1456  std::vector<G4DynamicParticle*>* secVector = 0;
1457  G4DynamicParticle* aSecondary = 0;
1458  G4ParticleDefinition* type = 0;
1459  G4double e, tkin, grej;
1460  G4ThreeVector position;
1461  G4int shell, shellId;
1462
1463  // sample secondaries
1464
1465  G4double etot = 0.0;
1466  std::vector<G4int> n = shellVacancy->GenerateNumberOfIonisations(couple,
1467                                         incidentEnergy, eLoss);
1468
1469  for (size_t i=0; i<nElements; i++) {
1470
1471    size_t nVacancies = n[i];
1472    G4int Z = (G4int)((*theElementVector)[i]->GetZ());
1473    G4double maxE = transitionManager->Shell(Z, 0)->BindingEnergy();
1474
1475    if (nVacancies && Z  > 5 && maxE < tcut && (maxE > minGammaEnergy || maxE > minElectronEnergy)) {
1476      for(size_t j=0; j<nVacancies; j++) {
1477
1478        // sampling follows
1479        do {
1480          tkin = tcut/(1.0 + (tcut/maxE - 1.0)*G4UniformRand());
1481          grej = 1.0 - beta2 * tkin/tmax;
1482
1483        } while( G4UniformRand() > grej );
1484
1485        shell = shellCS->SelectRandomShell(Z,incidentEnergy,hMass,tkin);
1486
1487        shellId = transitionManager->Shell(Z, shell)->ShellId();
1488        G4double maxE = transitionManager->Shell(Z, shell)->BindingEnergy();
1489
1490        if (maxE>minGammaEnergy || maxE>minElectronEnergy ) {
1491          secVector = deexcitationManager.GenerateParticles(Z, shellId);
1492        } else {
1493          secVector = 0;
1494        }
1495
1496        if (secVector) {
1497
1498          for (size_t l = 0; l<secVector->size(); l++) {
1499
1500            aSecondary = (*secVector)[l];
1501            if(aSecondary) {
1502
1503              e = aSecondary->GetKineticEnergy();
1504              type = aSecondary->GetDefinition();
1505              if ( etot + e <= eLoss &&
1506                   ( (type == G4Gamma::Gamma() && e > minGammaEnergy ) ||
1507                   (type == G4Electron::Electron() && e > minElectronEnergy) ) ) {
1508
1509                     etot += e;
1510                     partVector->push_back(aSecondary);
1511
1512              } else {
1513                     delete aSecondary;
1514              }
1515            }
1516          }
1517          delete secVector;
1518        }
1519      }
1520    }
1521  }
1522
1523  if(partVector->empty()) {
1524    delete partVector;
1525    return 0;
1526  }
1527
1528  return partVector;
1529}
1530
1531//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1532
1533G4int G4hLowEnergyIonisation::SelectRandomAtom(const G4MaterialCutsCouple* couple,
1534                                                     G4double kineticEnergy) const
1535{
1536  const G4Material* material = couple->GetMaterial();
1537  G4int nElements = material->GetNumberOfElements();
1538  G4int Z = 0;
1539
1540  if(nElements == 1) {
1541    Z = (G4int)(material->GetZ());
1542    return Z;
1543  }
1544
1545  const G4ElementVector* theElementVector = material->GetElementVector();
1546  std::vector<G4double> p;
1547  G4int index = couple->GetIndex();
1548
1549  G4double norm = 0.0;
1550  for (G4int j=0; j<nElements; j++) {
1551
1552    const G4VEMDataSet* set = (zFluoDataVector[index])->GetComponent(j);
1553    G4double cross    = set->FindValue(kineticEnergy);
1554
1555    p.push_back(cross);
1556    norm += cross;
1557  }
1558
1559  if(norm == 0.0) return 0;
1560
1561  G4double q = norm*G4UniformRand();
1562
1563  for (G4int i=0; i<nElements; i++) {
1564
1565    if(p[i] > q) {
1566       Z = (G4int)((*theElementVector)[i]->GetZ());
1567       break;
1568    }
1569    q -= p[i];
1570  }
1571
1572  return Z;
1573}
1574
1575//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1576
1577G4double G4hLowEnergyIonisation::ComputeDEDX(
1578                                 const G4ParticleDefinition* aParticle,
1579                                 const G4MaterialCutsCouple* couple,
1580                                 G4double kineticEnergy)
1581{
1582  const G4Material* material = couple->GetMaterial();
1583  G4Proton* theProton = G4Proton::Proton();
1584  G4AntiProton* theAntiProton = G4AntiProton::AntiProton();
1585  G4double dedx = 0.0 ;
1586
1587  G4double tscaled = kineticEnergy*proton_mass_c2/(aParticle->GetPDGMass()) ;
1588  charge  = aParticle->GetPDGCharge() ;
1589
1590  if(charge>0.0) {
1591    if(tscaled > protonHighEnergy) {
1592      dedx=G4EnergyLossTables::GetDEDX(theProton,tscaled,couple) ;
1593
1594    } else {
1595      dedx=ProtonParametrisedDEDX(couple,tscaled) ;
1596    }
1597
1598  } else {
1599    if(tscaled > antiProtonHighEnergy) {
1600      dedx=G4EnergyLossTables::GetDEDX(theAntiProton,tscaled,couple);
1601
1602    } else {
1603      dedx=AntiProtonParametrisedDEDX(couple,tscaled) ;
1604    }
1605  }
1606  dedx *= theIonEffChargeModel->TheValue(aParticle, material, kineticEnergy) ;
1607
1608  return dedx ;
1609}
1610
1611//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1612
1613G4double G4hLowEnergyIonisation::BarkasTerm(const G4Material* material,
1614                                                  G4double kineticEnergy) const
1615//Function to compute the Barkas term for protons:
1616//
1617//Ref. Z_1^3 effect in the stopping power of matter for charged particles
1618//     J.C Ashley and R.H.Ritchie
1619//     Physical review B Vol.5 No.7 1 April 1972 pagg. 2393-2397
1620//
1621{
1622  static double FTable[47][2] = {
1623   { 0.02, 21.5},
1624   { 0.03, 20.0},
1625   { 0.04, 18.0},
1626   { 0.05, 15.6},
1627   { 0.06, 15.0},
1628   { 0.07, 14.0},
1629   { 0.08, 13.5},
1630   { 0.09, 13.},
1631   { 0.1,  12.2},
1632   { 0.2,  9.25},
1633   { 0.3,  7.0},
1634   { 0.4,  6.0},
1635   { 0.5,  4.5},
1636   { 0.6,  3.5},
1637   { 0.7,  3.0},
1638   { 0.8,  2.5},
1639   { 0.9,  2.0},
1640   { 1.0,  1.7},
1641   { 1.2,  1.2},
1642   { 1.3,  1.0},
1643   { 1.4,  0.86},
1644   { 1.5,  0.7},
1645   { 1.6,  0.61},
1646   { 1.7,  0.52},
1647   { 1.8,  0.5},
1648   { 1.9,  0.43},
1649   { 2.0,  0.42},
1650   { 2.1,  0.3},
1651   { 2.4,  0.2},
1652   { 3.0,  0.13},
1653   { 3.08, 0.1},
1654   { 3.1,  0.09},
1655   { 3.3,  0.08},
1656   { 3.5,  0.07},
1657   { 3.8,  0.06},
1658   { 4.0,  0.051},
1659   { 4.1,  0.04},
1660   { 4.8,  0.03},
1661   { 5.0,  0.024},
1662   { 5.1,  0.02},
1663   { 6.0,  0.013},
1664   { 6.5,  0.01},
1665   { 7.0,  0.009},
1666   { 7.1,  0.008},
1667   { 8.0,  0.006},
1668   { 9.0,  0.0032},
1669   { 10.0, 0.0025} };
1670
1671  // Information on particle and material
1672  G4double kinE  = kineticEnergy ;
1673  if(0.5*MeV > kinE) kinE = 0.5*MeV ;
1674  G4double gamma = 1.0 + kinE / proton_mass_c2 ;
1675  G4double beta2 = 1.0 - 1.0/(gamma*gamma) ;
1676  if(0.0 >= beta2) return 0.0;
1677
1678  G4double BarkasTerm = 0.0;
1679  G4double AMaterial = 0.0;
1680  G4double ZMaterial = 0.0;
1681  const G4ElementVector* theElementVector = material->GetElementVector();
1682  G4int numberOfElements = material->GetNumberOfElements();
1683
1684  for (G4int i = 0; i<numberOfElements; i++) {
1685
1686    AMaterial = (*theElementVector)[i]->GetA()*mole/g;
1687    ZMaterial = (*theElementVector)[i]->GetZ();
1688
1689    G4double X = 137.0 * 137.0 * beta2 / ZMaterial;
1690
1691    // Variables to compute L_1
1692    G4double Eta0Chi = 0.8;
1693    G4double EtaChi = Eta0Chi * ( 1.0 + 6.02*std::pow( ZMaterial,-1.19 ) );
1694    G4double W = ( EtaChi * std::pow( ZMaterial,1.0/6.0 ) ) / std::sqrt(X);
1695    G4double FunctionOfW = FTable[46][1]*FTable[46][0]/W ;
1696
1697    for(G4int j=0; j<47; j++) {
1698
1699      if( W < FTable[j][0] ) {
1700
1701        if(0 == j) {
1702          FunctionOfW = FTable[0][1] ;
1703
1704        } else {
1705          FunctionOfW = (FTable[j][1] - FTable[j-1][1]) * (W - FTable[j-1][0])
1706                      / (FTable[j][0] - FTable[j-1][0])
1707                      +  FTable[j-1][1] ;
1708        }
1709
1710        break;
1711      }
1712
1713    }
1714
1715    BarkasTerm += FunctionOfW /( std::sqrt(ZMaterial * X) * X);
1716  }
1717
1718  BarkasTerm *= twopi_mc2_rcl2 * (material->GetElectronDensity()) / beta2 ;
1719
1720  return BarkasTerm;
1721}
1722
1723//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1724
1725G4double G4hLowEnergyIonisation::BlochTerm(const G4Material* material,
1726                                                 G4double kineticEnergy,
1727                                                 G4double cSquare) const
1728//Function to compute the Bloch term for protons:
1729//
1730//Ref. Z_1^3 effect in the stopping power of matter for charged particles
1731//     J.C Ashley and R.H.Ritchie
1732//     Physical review B Vol.5 No.7 1 April 1972 pagg. 2393-2397
1733//
1734{
1735  G4double eloss = 0.0 ;
1736  G4double gamma = 1.0 + kineticEnergy / proton_mass_c2 ;
1737  G4double beta2 = 1.0 - 1.0/(gamma*gamma) ;
1738  G4double y = cSquare / (137.0*137.0*beta2) ;
1739
1740  if(y < 0.05) {
1741    eloss = 1.202 ;
1742
1743  } else {
1744    eloss = 1.0 / (1.0 + y) ;
1745    G4double de = eloss ;
1746
1747    for(G4int i=2; de>eloss*0.01; i++) {
1748      de = 1.0/( i * (i*i + y)) ;
1749      eloss += de ;
1750    }
1751  }
1752  eloss *= -1.0 * y * cSquare * twopi_mc2_rcl2 *
1753            (material->GetElectronDensity()) / beta2 ;
1754
1755  return eloss;
1756}
1757
1758//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1759
1760G4double G4hLowEnergyIonisation::ElectronicLossFluctuation(
1761                                 const G4DynamicParticle* particle,
1762                                 const G4MaterialCutsCouple* couple,
1763                                       G4double meanLoss,
1764                                       G4double step) const
1765//  calculate actual loss from the mean loss
1766//  The model used to get the fluctuation is essentially the same
1767// as in Glandz in Geant3.
1768{
1769  // data members to speed up the fluctuation calculation
1770  //  G4int imat ;
1771  //  G4double f1Fluct,f2Fluct,e1Fluct,e2Fluct,rateFluct,ipotFluct;
1772  //  G4double e1LogFluct,e2LogFluct,ipotLogFluct;
1773
1774  static const G4double minLoss = 1.*eV ;
1775  static const G4double kappa = 10. ;
1776  static const G4double theBohrBeta2 = 50.0 * keV/proton_mass_c2 ;
1777
1778  const G4Material* material = couple->GetMaterial();
1779  G4int    imaterial   = couple->GetIndex() ;
1780  G4double ipotFluct   = material->GetIonisation()->GetMeanExcitationEnergy() ;
1781  G4double electronDensity = material->GetElectronDensity() ;
1782  G4double zeff = electronDensity/(material->GetTotNbOfAtomsPerVolume()) ;
1783
1784  // get particle data
1785  G4double tkin   = particle->GetKineticEnergy();
1786  G4double particleMass = particle->GetMass() ;
1787  G4double deltaCutInKineticEnergyNow = cutForDelta[imaterial];
1788
1789  // shortcut for very very small loss
1790  if(meanLoss < minLoss) return meanLoss ;
1791
1792  // Validity range for delta electron cross section
1793  G4double threshold = std::max(deltaCutInKineticEnergyNow,ipotFluct);
1794  G4double loss, siga;
1795
1796  G4double rmass = electron_mass_c2/particleMass;
1797  G4double tau   = tkin/particleMass;
1798  G4double tau1 = tau+1.0;
1799  G4double tau2 = tau*(tau+2.);
1800  G4double tmax    = 2.*electron_mass_c2*tau2/(1.+2.*tau1*rmass+rmass*rmass);
1801
1802
1803  if(tmax > threshold) tmax = threshold;
1804  G4double beta2 = tau2/(tau1*tau1);
1805
1806  // Gaussian fluctuation
1807  if(meanLoss > kappa*tmax || tmax < kappa*ipotFluct )
1808  {
1809    siga = tmax * (1.0-0.5*beta2) * step * twopi_mc2_rcl2
1810         * electronDensity / beta2 ;
1811
1812    // High velocity or negatively charged particle
1813    if( beta2 > 3.0*theBohrBeta2*zeff || charge < 0.0) {
1814      siga = std::sqrt( siga * chargeSquare ) ;
1815
1816    // Low velocity - additional ion charge fluctuations according to
1817    // Q.Yang et al., NIM B61(1991)149-155.
1818    } else {
1819      G4double chu = theIonChuFluctuationModel->TheValue(particle, material);
1820      G4double yang = theIonYangFluctuationModel->TheValue(particle, material);
1821      siga = std::sqrt( siga * (chargeSquare * chu + yang)) ;
1822    }
1823
1824    do {
1825        loss = G4RandGauss::shoot(meanLoss,siga);
1826    } while (loss < 0. || loss > 2.0*meanLoss);
1827    return loss;
1828  }
1829
1830  // Non Gaussian fluctuation
1831  static const G4double probLim = 0.01 ;
1832  static const G4double sumaLim = -std::log(probLim) ;
1833  static const G4double alim = 10.;
1834
1835  G4double suma,w1,w2,C,e0,lossc,w;
1836  G4double a1,a2,a3;
1837  G4int p1,p2,p3;
1838  G4int nb;
1839  G4double corrfac, na,alfa,rfac,namean,sa,alfa1,ea,sea;
1840  G4double dp3;
1841
1842  G4double f1Fluct     = material->GetIonisation()->GetF1fluct();
1843  G4double f2Fluct     = material->GetIonisation()->GetF2fluct();
1844  G4double e1Fluct     = material->GetIonisation()->GetEnergy1fluct();
1845  G4double e2Fluct     = material->GetIonisation()->GetEnergy2fluct();
1846  G4double e1LogFluct  = material->GetIonisation()->GetLogEnergy1fluct();
1847  G4double e2LogFluct  = material->GetIonisation()->GetLogEnergy2fluct();
1848  G4double rateFluct   = material->GetIonisation()->GetRateionexcfluct();
1849  G4double ipotLogFluct= material->GetIonisation()->GetLogMeanExcEnergy();
1850
1851  w1 = tmax/ipotFluct;
1852  w2 = std::log(2.*electron_mass_c2*tau2);
1853
1854  C = meanLoss*(1.-rateFluct)/(w2-ipotLogFluct-beta2);
1855
1856  a1 = C*f1Fluct*(w2-e1LogFluct-beta2)/e1Fluct;
1857  a2 = C*f2Fluct*(w2-e2LogFluct-beta2)/e2Fluct;
1858  a3 = rateFluct*meanLoss*(tmax-ipotFluct)/(ipotFluct*tmax*std::log(w1));
1859  if(a1 < 0.0) a1 = 0.0;
1860  if(a2 < 0.0) a2 = 0.0;
1861  if(a3 < 0.0) a3 = 0.0;
1862
1863  suma = a1+a2+a3;
1864
1865  loss = 0.;
1866
1867
1868  if(suma < sumaLim)             // very small Step
1869    {
1870      e0 = material->GetIonisation()->GetEnergy0fluct();
1871
1872      if(tmax == ipotFluct)
1873      {
1874        a3 = meanLoss/e0;
1875
1876        if(a3>alim)
1877        {
1878          siga=std::sqrt(a3) ;
1879          p3 = std::max(0,G4int(G4RandGauss::shoot(a3,siga)+0.5));
1880        }
1881        else
1882          p3 = G4Poisson(a3);
1883
1884        loss = p3*e0 ;
1885
1886        if(p3 > 0)
1887          loss += (1.-2.*G4UniformRand())*e0 ;
1888
1889      }
1890      else
1891      {
1892        tmax = tmax-ipotFluct+e0 ;
1893        a3 = meanLoss*(tmax-e0)/(tmax*e0*std::log(tmax/e0));
1894
1895        if(a3>alim)
1896        {
1897          siga=std::sqrt(a3) ;
1898          p3 = std::max(0,int(G4RandGauss::shoot(a3,siga)+0.5));
1899        }
1900        else
1901          p3 = G4Poisson(a3);
1902
1903        if(p3 > 0)
1904        {
1905          w = (tmax-e0)/tmax ;
1906          if(p3 > nmaxCont2)
1907          {
1908            dp3 = G4float(p3) ;
1909            corrfac = dp3/G4float(nmaxCont2) ;
1910            p3 = nmaxCont2 ;
1911          }
1912          else
1913            corrfac = 1. ;
1914
1915          for(G4int i=0; i<p3; i++) loss += 1./(1.-w*G4UniformRand()) ;
1916          loss *= e0*corrfac ;
1917        }
1918      }
1919    }
1920
1921  else                              // not so small Step
1922    {
1923      // excitation type 1
1924      if(a1>alim)
1925      {
1926        siga=std::sqrt(a1) ;
1927        p1 = std::max(0,G4int(G4RandGauss::shoot(a1,siga)+0.5));
1928      }
1929      else
1930       p1 = G4Poisson(a1);
1931
1932      // excitation type 2
1933      if(a2>alim)
1934      {
1935        siga=std::sqrt(a2) ;
1936        p2 = std::max(0,G4int(G4RandGauss::shoot(a2,siga)+0.5));
1937      }
1938      else
1939        p2 = G4Poisson(a2);
1940
1941      loss = p1*e1Fluct+p2*e2Fluct;
1942
1943      // smearing to avoid unphysical peaks
1944      if(p2 > 0)
1945        loss += (1.-2.*G4UniformRand())*e2Fluct;
1946      else if (loss>0.)
1947        loss += (1.-2.*G4UniformRand())*e1Fluct;
1948
1949      // ionisation .......................................
1950     if(a3 > 0.)
1951     {
1952      if(a3>alim)
1953      {
1954        siga=std::sqrt(a3) ;
1955        p3 = std::max(0,G4int(G4RandGauss::shoot(a3,siga)+0.5));
1956      }
1957      else
1958        p3 = G4Poisson(a3);
1959
1960      lossc = 0.;
1961      if(p3 > 0)
1962      {
1963        na = 0.;
1964        alfa = 1.;
1965        if (p3 > nmaxCont2)
1966        {
1967          dp3        = G4float(p3);
1968          rfac       = dp3/(G4float(nmaxCont2)+dp3);
1969          namean     = G4float(p3)*rfac;
1970          sa         = G4float(nmaxCont1)*rfac;
1971          na         = G4RandGauss::shoot(namean,sa);
1972          if (na > 0.)
1973          {
1974            alfa   = w1*G4float(nmaxCont2+p3)/
1975                    (w1*G4float(nmaxCont2)+G4float(p3));
1976            alfa1  = alfa*std::log(alfa)/(alfa-1.);
1977            ea     = na*ipotFluct*alfa1;
1978            sea    = ipotFluct*std::sqrt(na*(alfa-alfa1*alfa1));
1979            lossc += G4RandGauss::shoot(ea,sea);
1980          }
1981        }
1982
1983        nb = G4int(G4float(p3)-na);
1984        if (nb > 0)
1985        {
1986          w2 = alfa*ipotFluct;
1987          w  = (tmax-w2)/tmax;
1988          for (G4int k=0; k<nb; k++) lossc += w2/(1.-w*G4UniformRand());
1989        }
1990      }
1991      loss += lossc;
1992     }
1993    }
1994
1995  return loss ;
1996}
1997
1998//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1999
2000void G4hLowEnergyIonisation::SetCutForSecondaryPhotons(G4double cut)
2001{
2002  minGammaEnergy = cut;
2003}
2004
2005//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
2006
2007void G4hLowEnergyIonisation::SetCutForAugerElectrons(G4double cut)
2008{
2009  minElectronEnergy = cut;
2010}
2011
2012//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
2013
2014void G4hLowEnergyIonisation::ActivateAugerElectronProduction(G4bool val)
2015{
2016  deexcitationManager.ActivateAugerElectronProduction(val);
2017}
2018
2019//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
2020
2021void G4hLowEnergyIonisation::PrintInfoDefinition() const
2022{
2023  G4String comments = "  Knock-on electron cross sections . ";
2024  comments += "\n        Good description above the mean excitation energy.\n";
2025  comments += "        Delta ray energy sampled from  differential Xsection.";
2026
2027  G4cout << G4endl << GetProcessName() << ":  " << comments
2028         << "\n        PhysicsTables from " << LowestKineticEnergy / eV << " eV "
2029         << " to " << HighestKineticEnergy / TeV << " TeV "
2030         << " in " << TotBin << " bins."
2031 << "\n        Electronic stopping power model is  "
2032 << theProtonTable
2033         << "\n        from " << protonLowEnergy / keV << " keV "
2034         << " to " << protonHighEnergy / MeV << " MeV " << "." << G4endl ;
2035  G4cout << "\n        Parametrisation model for antiprotons is  "
2036         << theAntiProtonTable
2037         << "\n        from " << antiProtonLowEnergy / keV << " keV "
2038         << " to " << antiProtonHighEnergy / MeV << " MeV " << "." << G4endl ;
2039  if(theBarkas){
2040  G4cout << "        Parametrization of the Barkas effect is switched on."
2041         << G4endl ;
2042  }
2043  if(nStopping) {
2044  G4cout << "        Nuclear stopping power model is " << theNuclearTable
2045         << G4endl ;
2046  }
2047
2048  G4bool printHead = true;
2049
2050  const G4ProductionCutsTable* theCoupleTable=
2051        G4ProductionCutsTable::GetProductionCutsTable();
2052  size_t numOfCouples = theCoupleTable->GetTableSize();
2053
2054  // loop for materials
2055
2056  for (size_t j=0 ; j < numOfCouples; j++) {
2057
2058    const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(j);
2059    const G4Material* material= couple->GetMaterial();
2060    G4double deltaCutNow = cutForDelta[(couple->GetIndex())] ;
2061    G4double eexc = material->GetIonisation()->GetMeanExcitationEnergy();
2062
2063    if(eexc > deltaCutNow) {
2064      if(printHead) {
2065        printHead = false ;
2066
2067        G4cout << "           material       min.delta energy(keV) " << G4endl;
2068        G4cout << G4endl;
2069      }
2070
2071      G4cout << std::setw(20) << material->GetName()
2072             << std::setw(15) << eexc/keV << G4endl;
2073    }
2074  }
2075}
Note: See TracBrowser for help on using the repository browser.