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

Last change on this file since 1315 was 1315, checked in by garnier, 14 years ago

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