Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Public Member Functions | Protected Attributes | List of all members
G4DNABornIonisationModel1 Class Reference

#include <G4DNABornIonisationModel1.hh>

Inheritance diagram for G4DNABornIonisationModel1:
Inheritance graph
[legend]
Collaboration diagram for G4DNABornIonisationModel1:
Collaboration graph
[legend]

Public Member Functions

 G4DNABornIonisationModel1 (const G4ParticleDefinition *p=0, const G4String &nam="DNABornIonisationModel")
 
virtual ~G4DNABornIonisationModel1 ()
 
virtual void Initialise (const G4ParticleDefinition *, const G4DataVector &=*(new G4DataVector()))
 
virtual G4double CrossSectionPerVolume (const G4Material *material, const G4ParticleDefinition *p, G4double ekin, G4double emin, G4double emax)
 
virtual void SampleSecondaries (std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
 
virtual G4double GetPartialCrossSection (const G4Material *, G4int, const G4ParticleDefinition *, G4double)
 
double DifferentialCrossSection (G4ParticleDefinition *aParticleDefinition, G4double k, G4double energyTransfer, G4int shell)
 
G4double TransferedEnergy (G4ParticleDefinition *aParticleDefinition, G4double incomingParticleEnergy, G4int shell, G4double random)
 
void SelectFasterComputation (G4bool input)
 
void SelectStationary (G4bool input)
 
void SelectSPScaling (G4bool input)
 
- Public Member Functions inherited from G4VEmModel
 G4VEmModel (const G4String &nam)
 
virtual ~G4VEmModel ()
 
virtual void InitialiseLocal (const G4ParticleDefinition *, G4VEmModel *masterModel)
 
virtual void InitialiseForMaterial (const G4ParticleDefinition *, const G4Material *)
 
virtual void InitialiseForElement (const G4ParticleDefinition *, G4int Z)
 
virtual G4double ComputeDEDXPerVolume (const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=DBL_MAX)
 
virtual G4double ComputeCrossSectionPerAtom (const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0., G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
virtual G4double ComputeCrossSectionPerShell (const G4ParticleDefinition *, G4int Z, G4int shellIdx, G4double kinEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
virtual G4double ChargeSquareRatio (const G4Track &)
 
virtual G4double GetChargeSquareRatio (const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
 
virtual G4double GetParticleCharge (const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
 
virtual void StartTracking (G4Track *)
 
virtual void CorrectionsAlongStep (const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double &eloss, G4double &niel, G4double length)
 
virtual G4double Value (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy)
 
virtual G4double MinPrimaryEnergy (const G4Material *, const G4ParticleDefinition *, G4double cut=0.0)
 
virtual G4double MinEnergyCut (const G4ParticleDefinition *, const G4MaterialCutsCouple *)
 
virtual void SetupForMaterial (const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
 
virtual void DefineForRegion (const G4Region *)
 
virtual void ModelDescription (std::ostream &outFile) const
 
void InitialiseElementSelectors (const G4ParticleDefinition *, const G4DataVector &)
 
std::vector
< G4EmElementSelector * > * 		GetElementSelectors ()
 
void SetElementSelectors (std::vector< G4EmElementSelector * > *)
 
virtual G4double ComputeDEDX (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=DBL_MAX)
 
G4double CrossSection (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
G4double ComputeMeanFreePath (const G4ParticleDefinition *, G4double kineticEnergy, const G4Material *, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
G4double ComputeCrossSectionPerAtom (const G4ParticleDefinition *, const G4Element *, G4double kinEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
G4int SelectIsotopeNumber (const G4Element *)
 
const G4Element * SelectRandomAtom (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
const G4Element * SelectRandomAtom (const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
G4int SelectRandomAtomNumber (const G4Material *)
 
void SetParticleChange (G4VParticleChange *, G4VEmFluctuationModel *f=nullptr)
 
void SetCrossSectionTable (G4PhysicsTable *, G4bool isLocal)
 
G4ElementData * GetElementData ()
 
G4PhysicsTable * GetCrossSectionTable ()
 
G4VEmFluctuationModel * GetModelOfFluctuations ()
 
G4VEmAngularDistribution * GetAngularDistribution ()
 
void SetAngularDistribution (G4VEmAngularDistribution *)
 
G4double HighEnergyLimit () const
 
G4double LowEnergyLimit () const
 
G4double HighEnergyActivationLimit () const
 
G4double LowEnergyActivationLimit () const
 
G4double PolarAngleLimit () const
 
G4double SecondaryThreshold () const
 
G4bool LPMFlag () const
 
G4bool DeexcitationFlag () const
 
G4bool ForceBuildTableFlag () const
 
G4bool UseAngularGeneratorFlag () const
 
void SetAngularGeneratorFlag (G4bool)
 
void SetHighEnergyLimit (G4double)
 
void SetLowEnergyLimit (G4double)
 
void SetActivationHighEnergyLimit (G4double)
 
void SetActivationLowEnergyLimit (G4double)
 
G4bool IsActive (G4double kinEnergy)
 
void SetPolarAngleLimit (G4double)
 
void SetSecondaryThreshold (G4double)
 
void SetLPMFlag (G4bool val)
 
void SetDeexcitationFlag (G4bool val)
 
void SetForceBuildTable (G4bool val)
 
void SetFluctuationFlag (G4bool val)
 
void SetMasterThread (G4bool val)
 
G4bool IsMaster () const
 
G4double MaxSecondaryKinEnergy (const G4DynamicParticle *dynParticle)
 
const G4String & GetName () const
 
void SetCurrentCouple (const G4MaterialCutsCouple *)
 
const G4Element * GetCurrentElement () const
 
const G4Isotope * GetCurrentIsotope () const
 
G4bool IsLocked () const
 
void SetLocked (G4bool)
 

Protected Attributes

G4ParticleChangeForGamma * fParticleChangeForGamma
 
- Protected Attributes inherited from G4VEmModel
G4ElementData * fElementData
 
G4VParticleChange * pParticleChange
 
G4PhysicsTable * xSectionTable
 
const std::vector< G4double > * theDensityFactor
 
const std::vector< G4int > * theDensityIdx
 
size_t idxTable
 
G4bool lossFlucFlag
 

Additional Inherited Members

- Protected Member Functions inherited from G4VEmModel
G4ParticleChangeForLoss * GetParticleChangeForLoss ()
 
G4ParticleChangeForGamma * GetParticleChangeForGamma ()
 
virtual G4double MaxSecondaryEnergy (const G4ParticleDefinition *, G4double kineticEnergy)
 
const G4MaterialCutsCouple * CurrentCouple () const
 
void SetCurrentElement (const G4Element *)
 
- Static Protected Attributes inherited from G4VEmModel
static const G4double inveplus = 1.0/CLHEP::eplus
 

Detailed Description

Definition at line 48 of file G4DNABornIonisationModel1.hh.

Constructor & Destructor Documentation

G4DNABornIonisationModel1::G4DNABornIonisationModel1 ( const G4ParticleDefinition *  p = 0,
const G4String &  nam = "DNABornIonisationModel" 
)

Definition at line 46 of file G4DNABornIonisationModel1.cc.

47  :
48  G4VEmModel(nam), isInitialised(false)
49 {
50  verboseLevel = 0;
51  // Verbosity scale:
52  // 0 = nothing
53  // 1 = warning for energy non-conservation
54  // 2 = details of energy budget
55  // 3 = calculation of cross sections, file openings, sampling of atoms
56  // 4 = entering in methods
57 
58  if (verboseLevel > 0)
59  {
60  G4cout << "Born ionisation model is constructed " << G4endl;
61  }
62 
63  //Mark this model as "applicable" for atomic deexcitation
64  SetDeexcitationFlag(true);
65  fAtomDeexcitation = 0;
66  fParticleChangeForGamma = 0;
67  fpMolWaterDensity = 0;
68 
69  // define default angular generator
70  SetAngularDistribution(new G4DNABornAngle());
71 
72  // Selection of computation method
73 
74  fasterCode = false;
75 
76  // Selection of stationary mode
77 
78  statCode = false;
79 
80  // Selection of SP scaling
81 
82  spScaling = true;
83 }
G4VEmModel::G4VEmModel
G4VEmModel(const G4String &nam)
Definition: G4VEmModel.cc:68
G4DNABornIonisationModel1::fParticleChangeForGamma
G4ParticleChangeForGamma * fParticleChangeForGamma
Definition: G4DNABornIonisationModel1.hh:90
G4cout
G4GLOB_DLL std::ostream G4cout
G4VEmModel::SetAngularDistribution
void SetAngularDistribution(G4VEmAngularDistribution *)
Definition: G4VEmModel.hh:626
G4endl
#define G4endl
Definition: G4ios.hh:61
G4VEmModel::SetDeexcitationFlag
void SetDeexcitationFlag(G4bool val)
Definition: G4VEmModel.hh:788

Here is the call graph for this function:

G4DNABornIonisationModel1::~G4DNABornIonisationModel1 ( )
virtual

Definition at line 87 of file G4DNABornIonisationModel1.cc.

88 {
89  // Cross section
90 
91  std::map<G4String, G4DNACrossSectionDataSet*, std::less<G4String> >::iterator pos;
92  for (pos = tableData.begin(); pos != tableData.end(); ++pos)
93  {
94  G4DNACrossSectionDataSet* table = pos->second;
95  delete table;
96  }
97 
98  // Final state
99 
100  eVecm.clear();
101  pVecm.clear();
102 }
pos
static const G4double pos
Definition: G4ElectroNuclearCrossSection.cc:66

Member Function Documentation

G4double G4DNABornIonisationModel1::CrossSectionPerVolume ( const G4Material *  material,
const G4ParticleDefinition *  p,
G4double  ekin,
G4double  emin,
G4double  emax 
)
virtual

Reimplemented from G4VEmModel.

Definition at line 315 of file G4DNABornIonisationModel1.cc.

320 {
321  if (verboseLevel > 3)
322  {
323  G4cout << "Calling CrossSectionPerVolume() of G4DNABornIonisationModel1"
324  << G4endl;
325 
326  }
327 
328  if (
329  particleDefinition != G4Proton::ProtonDefinition()
330  &&
331  particleDefinition != G4Electron::ElectronDefinition()
332  )
333 
334  return 0;
335 
336  // Calculate total cross section for model
337 
338  G4double lowLim = 0;
339  G4double highLim = 0;
340  G4double sigma=0;
341 
342  G4double waterDensity = (*fpMolWaterDensity)[material->GetIndex()];
343 
344  if(waterDensity!= 0.0)
345  // if (material == nistwater || material->GetBaseMaterial() == nistwater)
346  {
347  const G4String& particleName = particleDefinition->GetParticleName();
348 
349  std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
350  pos1 = lowEnergyLimit.find(particleName);
351  if (pos1 != lowEnergyLimit.end())
352  {
353  lowLim = pos1->second;
354  }
355 
356  std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
357  pos2 = highEnergyLimit.find(particleName);
358  if (pos2 != highEnergyLimit.end())
359  {
360  highLim = pos2->second;
361  }
362 
363  if (ekin >= lowLim && ekin < highLim)
364  {
365  std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
366  pos = tableData.find(particleName);
367 
368  if (pos != tableData.end())
369  {
370  G4DNACrossSectionDataSet* table = pos->second;
371  if (table != 0)
372  {
373  sigma = table->FindValue(ekin);
374 
375  // ICRU49 electronic SP scaling - ZF, SI
376 
377  if (particleDefinition == G4Proton::ProtonDefinition() && ekin < 70*MeV && spScaling)
378  {
379  G4double A = 1.39241700556072800000E-009 ;
380  G4double B = -8.52610412942622630000E-002 ;
381  sigma = sigma * G4Exp(A*(ekin/eV)+B);
382  }
383  //
384 
385  }
386  }
387  else
388  {
389  G4Exception("G4DNABornIonisationModel1::CrossSectionPerVolume","em0002",
390  FatalException,"Model not applicable to particle type.");
391  }
392  }
393 
394  if (verboseLevel > 2)
395  {
396  G4cout << "__________________________________" << G4endl;
397  G4cout << "G4DNABornIonisationModel1 - XS INFO START" << G4endl;
398  G4cout << "Kinetic energy(eV)=" << ekin/eV << " particle : " << particleName << G4endl;
399  G4cout << "Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
400  G4cout << "Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) << G4endl;
401  G4cout << "G4DNABornIonisationModel1 - XS INFO END" << G4endl;
402  }
403  } // if (waterMaterial)
404 
405  return sigma*waterDensity;
406 }
G4Electron::ElectronDefinition
static G4Electron * ElectronDefinition()
Definition: G4Electron.cc:89
G4Material::GetIndex
size_t GetIndex() const
Definition: G4Material.hh:262
G4Proton::ProtonDefinition
static G4Proton * ProtonDefinition()
Definition: G4Proton.cc:88
B
double B(double temperature)
Definition: G4DNAElectronHoleRecombination.cc:69
G4cout
G4GLOB_DLL std::ostream G4cout
A
double A(double temperature)
Definition: G4DNAElectronHoleRecombination.cc:60
cm
static constexpr double cm
Definition: G4SIunits.hh:119
eV
static constexpr double eV
Definition: G4SIunits.hh:215
G4DNACrossSectionDataSet::FindValue
virtual G4double FindValue(G4double e, G4int componentId=0) const
Definition: G4DNACrossSectionDataSet.cc:484
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4Exp
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
FatalException
Definition: G4ExceptionSeverity.hh:60
G4endl
#define G4endl
Definition: G4ios.hh:61
MeV
static constexpr double MeV
Definition: G4SIunits.hh:214
G4double
double G4double
Definition: G4Types.hh:76
pos
static const G4double pos
Definition: G4ElectroNuclearCrossSection.cc:66

Here is the call graph for this function:

double G4DNABornIonisationModel1::DifferentialCrossSection ( G4ParticleDefinition *  aParticleDefinition,
G4double  k,
G4double  energyTransfer,
G4int  shell 
)

Definition at line 725 of file G4DNABornIonisationModel1.cc.

729 {
730  G4double sigma = 0.;
731 
732  if (energyTransfer >= waterStructure.IonisationEnergy(ionizationLevelIndex))
733  {
734  G4double valueT1 = 0;
735  G4double valueT2 = 0;
736  G4double valueE21 = 0;
737  G4double valueE22 = 0;
738  G4double valueE12 = 0;
739  G4double valueE11 = 0;
740 
741  G4double xs11 = 0;
742  G4double xs12 = 0;
743  G4double xs21 = 0;
744  G4double xs22 = 0;
745 
746  if (particleDefinition == G4Electron::ElectronDefinition())
747  {
748  // k should be in eV and energy transfer eV also
749 
750  std::vector<double>::iterator t2 = std::upper_bound(eTdummyVec.begin(),
751  eTdummyVec.end(),
752  k);
753 
754  std::vector<double>::iterator t1 = t2 - 1;
755 
756  // SI : the following condition avoids situations where energyTransfer >last vector element
757  if (energyTransfer <= eVecm[(*t1)].back()
758  && energyTransfer <= eVecm[(*t2)].back())
759  {
760  std::vector<double>::iterator e12 =
761  std::upper_bound(eVecm[(*t1)].begin(),
762  eVecm[(*t1)].end(),
763  energyTransfer);
764  std::vector<double>::iterator e11 = e12 - 1;
765 
766  std::vector<double>::iterator e22 =
767  std::upper_bound(eVecm[(*t2)].begin(),
768  eVecm[(*t2)].end(),
769  energyTransfer);
770  std::vector<double>::iterator e21 = e22 - 1;
771 
772  valueT1 = *t1;
773  valueT2 = *t2;
774  valueE21 = *e21;
775  valueE22 = *e22;
776  valueE12 = *e12;
777  valueE11 = *e11;
778 
779  xs11 = eDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE11];
780  xs12 = eDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE12];
781  xs21 = eDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE21];
782  xs22 = eDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE22];
783 
784  }
785 
786  }
787 
788  if (particleDefinition == G4Proton::ProtonDefinition())
789  {
790  // k should be in eV and energy transfer eV also
791  std::vector<double>::iterator t2 = std::upper_bound(pTdummyVec.begin(),
792  pTdummyVec.end(),
793  k);
794  std::vector<double>::iterator t1 = t2 - 1;
795 
796  std::vector<double>::iterator e12 = std::upper_bound(pVecm[(*t1)].begin(),
797  pVecm[(*t1)].end(),
798  energyTransfer);
799  std::vector<double>::iterator e11 = e12 - 1;
800 
801  std::vector<double>::iterator e22 = std::upper_bound(pVecm[(*t2)].begin(),
802  pVecm[(*t2)].end(),
803  energyTransfer);
804  std::vector<double>::iterator e21 = e22 - 1;
805 
806  valueT1 = *t1;
807  valueT2 = *t2;
808  valueE21 = *e21;
809  valueE22 = *e22;
810  valueE12 = *e12;
811  valueE11 = *e11;
812 
813  xs11 = pDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE11];
814  xs12 = pDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE12];
815  xs21 = pDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE21];
816  xs22 = pDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE22];
817 
818  }
819 
820  G4double xsProduct = xs11 * xs12 * xs21 * xs22;
821  if (xsProduct != 0.)
822  {
823  sigma = QuadInterpolator(valueE11,
824  valueE12,
825  valueE21,
826  valueE22,
827  xs11,
828  xs12,
829  xs21,
830  xs22,
831  valueT1,
832  valueT2,
833  k,
834  energyTransfer);
835  }
836  }
837 
838  return sigma;
839 }
G4Electron::ElectronDefinition
static G4Electron * ElectronDefinition()
Definition: G4Electron.cc:89
G4Proton::ProtonDefinition
static G4Proton * ProtonDefinition()
Definition: G4Proton.cc:88
plottest35.t1
tuple t1
Definition: plottest35.py:33
G4double
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

G4double G4DNABornIonisationModel1::GetPartialCrossSection ( const G4Material *  ,
G4int  level,
const G4ParticleDefinition *  particle,
G4double  kineticEnergy 
)
virtual

Reimplemented from G4VEmModel.

Definition at line 931 of file G4DNABornIonisationModel1.cc.

935 {
936  std::map<G4String, G4DNACrossSectionDataSet*, std::less<G4String> >::iterator pos;
937  pos = tableData.find(particle->GetParticleName());
938 
939  if (pos != tableData.end())
940  {
941  G4DNACrossSectionDataSet* table = pos->second;
942  return table->GetComponent(level)->FindValue(kineticEnergy);
943  }
944 
945  return 0;
946 }
G4VEMDataSet::FindValue
virtual G4double FindValue(G4double x, G4int componentId=0) const =0
G4DNACrossSectionDataSet::GetComponent
virtual const G4VEMDataSet * GetComponent(G4int componentId) const
Definition: G4DNACrossSectionDataSet.hh:69
G4ParticleDefinition::GetParticleName
const G4String & GetParticleName() const
Definition: G4ParticleDefinition.hh:120
pos
static const G4double pos
Definition: G4ElectroNuclearCrossSection.cc:66

Here is the call graph for this function:

void G4DNABornIonisationModel1::Initialise ( const G4ParticleDefinition *  particle,
const G4DataVector &  = *(new G4DataVector()) 
)
virtual

Implements G4VEmModel.

Definition at line 106 of file G4DNABornIonisationModel1.cc.

108 {
109 
110  if (verboseLevel > 3)
111  {
112  G4cout << "Calling G4DNABornIonisationModel1::Initialise()" << G4endl;
113  }
114 
115  // Energy limits
116 
117  G4String fileElectron("dna/sigma_ionisation_e_born");
118  G4String fileProton("dna/sigma_ionisation_p_born");
119 
120  G4ParticleDefinition* electronDef = G4Electron::ElectronDefinition();
121  G4ParticleDefinition* protonDef = G4Proton::ProtonDefinition();
122 
123  G4String electron;
124  G4String proton;
125 
126  G4double scaleFactor = (1.e-22 / 3.343) * m*m;
127 
128  char *path = getenv("G4LEDATA");
129 
130  // *** ELECTRON
131 
132  electron = electronDef->GetParticleName();
133 
134  tableFile[electron] = fileElectron;
135 
136  lowEnergyLimit[electron] = 11. * eV;
137  highEnergyLimit[electron] = 1. * MeV;
138 
139  // Cross section
140 
141  G4DNACrossSectionDataSet* tableE = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, eV,scaleFactor );
142  tableE->LoadData(fileElectron);
143 
144  tableData[electron] = tableE;
145 
146  // Final state
147 
148  std::ostringstream eFullFileName;
149 
150  if (fasterCode) eFullFileName << path << "/dna/sigmadiff_cumulated_ionisation_e_born_hp.dat";
151  if (!fasterCode) eFullFileName << path << "/dna/sigmadiff_ionisation_e_born.dat";
152 
153  std::ifstream eDiffCrossSection(eFullFileName.str().c_str());
154 
155  if (!eDiffCrossSection)
156  {
157  if (fasterCode) G4Exception("G4DNABornIonisationModel1::Initialise","em0003",
158  FatalException,"Missing data file:/dna/sigmadiff_cumulated_ionisation_e_born_hp.dat");
159 
160  if (!fasterCode) G4Exception("G4DNABornIonisationModel1::Initialise","em0003",
161  FatalException,"Missing data file:/dna/sigmadiff_ionisation_e_born.dat");
162  }
163 
164  //
165 
166  // Clear the arrays for re-initialization case (MT mode)
167  // March 25th, 2014 - Vaclav Stepan, Sebastien Incerti
168 
169  eTdummyVec.clear();
170  pTdummyVec.clear();
171 
172  eVecm.clear();
173  pVecm.clear();
174 
175  for (int j=0; j<5; j++)
176  {
177  eProbaShellMap[j].clear();
178  pProbaShellMap[j].clear();
179 
180  eDiffCrossSectionData[j].clear();
181  pDiffCrossSectionData[j].clear();
182 
183  eNrjTransfData[j].clear();
184  pNrjTransfData[j].clear();
185  }
186  //
187 
188  eTdummyVec.push_back(0.);
189  while(!eDiffCrossSection.eof())
190  {
191  double tDummy;
192  double eDummy;
193  eDiffCrossSection>>tDummy>>eDummy;
194  if (tDummy != eTdummyVec.back()) eTdummyVec.push_back(tDummy);
195 
196  double tmp;
197  for (int j=0; j<5; j++)
198  {
199  eDiffCrossSection>> tmp;
200 
201  eDiffCrossSectionData[j][tDummy][eDummy] = tmp;
202 
203  if (fasterCode)
204  {
205  eNrjTransfData[j][tDummy][eDiffCrossSectionData[j][tDummy][eDummy]]=eDummy;
206  eProbaShellMap[j][tDummy].push_back(eDiffCrossSectionData[j][tDummy][eDummy]);
207  }
208 
209  // SI - only if eof is not reached
210  if (!eDiffCrossSection.eof() && !fasterCode) eDiffCrossSectionData[j][tDummy][eDummy]*=scaleFactor;
211 
212  if (!fasterCode) eVecm[tDummy].push_back(eDummy);
213 
214  }
215  }
216 
217  // *** PROTON
218 
219  proton = protonDef->GetParticleName();
220 
221  tableFile[proton] = fileProton;
222 
223  lowEnergyLimit[proton] = 500. * keV;
224  highEnergyLimit[proton] = 100. * MeV;
225 
226  // Cross section
227 
228  G4DNACrossSectionDataSet* tableP = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, eV,scaleFactor );
229  tableP->LoadData(fileProton);
230 
231  tableData[proton] = tableP;
232 
233  // Final state
234 
235  std::ostringstream pFullFileName;
236 
237  if (fasterCode) pFullFileName << path << "/dna/sigmadiff_cumulated_ionisation_p_born_hp.dat";
238 
239  if (!fasterCode) pFullFileName << path << "/dna/sigmadiff_ionisation_p_born.dat";
240 
241  std::ifstream pDiffCrossSection(pFullFileName.str().c_str());
242 
243  if (!pDiffCrossSection)
244  {
245  if (fasterCode) G4Exception("G4DNABornIonisationModel1::Initialise","em0003",
246  FatalException,"Missing data file:/dna/sigmadiff_cumulated_ionisation_p_born_hp.dat");
247 
248  if (!fasterCode) G4Exception("G4DNABornIonisationModel1::Initialise","em0003",
249  FatalException,"Missing data file:/dna/sigmadiff_ionisation_p_born.dat");
250  }
251 
252  pTdummyVec.push_back(0.);
253  while(!pDiffCrossSection.eof())
254  {
255  double tDummy;
256  double eDummy;
257  pDiffCrossSection>>tDummy>>eDummy;
258  if (tDummy != pTdummyVec.back()) pTdummyVec.push_back(tDummy);
259  for (int j=0; j<5; j++)
260  {
261  pDiffCrossSection>>pDiffCrossSectionData[j][tDummy][eDummy];
262 
263  if (fasterCode)
264  {
265  pNrjTransfData[j][tDummy][pDiffCrossSectionData[j][tDummy][eDummy]]=eDummy;
266  pProbaShellMap[j][tDummy].push_back(pDiffCrossSectionData[j][tDummy][eDummy]);
267  }
268 
269  // SI - only if eof is not reached !
270  if (!pDiffCrossSection.eof() && !fasterCode) pDiffCrossSectionData[j][tDummy][eDummy]*=scaleFactor;
271 
272  if (!fasterCode) pVecm[tDummy].push_back(eDummy);
273  }
274  }
275 
276  //
277 
278  if (particle==electronDef)
279  {
280  SetLowEnergyLimit(lowEnergyLimit[electron]);
281  SetHighEnergyLimit(highEnergyLimit[electron]);
282  }
283 
284  if (particle==protonDef)
285  {
286  SetLowEnergyLimit(lowEnergyLimit[proton]);
287  SetHighEnergyLimit(highEnergyLimit[proton]);
288  }
289 
290  if( verboseLevel>0 )
291  {
292  G4cout << "Born ionisation model is initialized " << G4endl
293  << "Energy range: "
294  << LowEnergyLimit() / eV << " eV - "
295  << HighEnergyLimit() / keV << " keV for "
296  << particle->GetParticleName()
297  << G4endl;
298  }
299 
300  // Initialize water density pointer
301  fpMolWaterDensity = G4DNAMolecularMaterial::Instance()->
302  GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
303 
304  //
305  fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
306 
307  if (isInitialised)
308  { return;}
309  fParticleChangeForGamma = GetParticleChangeForGamma();
310  isInitialised = true;
311 }
G4Electron::ElectronDefinition
static G4Electron * ElectronDefinition()
Definition: G4Electron.cc:89
G4VEmModel::LowEnergyLimit
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:643
G4LossTableManager::Instance
static G4LossTableManager * Instance()
Definition: G4LossTableManager.cc:112
G4VEmModel::HighEnergyLimit
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:636
G4Material::GetMaterial
static G4Material * GetMaterial(const G4String &name, G4bool warning=true)
Definition: G4Material.cc:602
G4Proton::ProtonDefinition
static G4Proton * ProtonDefinition()
Definition: G4Proton.cc:88
G4DNACrossSectionDataSet::LoadData
virtual G4bool LoadData(const G4String &argFileName)
Definition: G4DNACrossSectionDataSet.cc:84
G4ParticleDefinition::GetParticleName
const G4String & GetParticleName() const
Definition: G4ParticleDefinition.hh:120
G4VEmModel::SetHighEnergyLimit
void SetHighEnergyLimit(G4double)
Definition: G4VEmModel.hh:732
G4DNABornIonisationModel1::fParticleChangeForGamma
G4ParticleChangeForGamma * fParticleChangeForGamma
Definition: G4DNABornIonisationModel1.hh:90
G4cout
G4GLOB_DLL std::ostream G4cout
m
static constexpr double m
Definition: G4SIunits.hh:129
eV
static constexpr double eV
Definition: G4SIunits.hh:215
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4DNAMolecularMaterial::Instance
static G4DNAMolecularMaterial * Instance()
Definition: G4DNAMolecularMaterial.cc:75
FatalException
Definition: G4ExceptionSeverity.hh:60
G4endl
#define G4endl
Definition: G4ios.hh:61
MeV
static constexpr double MeV
Definition: G4SIunits.hh:214
G4LossTableManager::AtomDeexcitation
G4VAtomDeexcitation * AtomDeexcitation()
Definition: G4LossTableManager.hh:444
G4double
double G4double
Definition: G4Types.hh:76
G4VEmModel::SetLowEnergyLimit
void SetLowEnergyLimit(G4double)
Definition: G4VEmModel.hh:739
keV
static constexpr double keV
Definition: G4SIunits.hh:216
G4VEmModel::GetParticleChangeForGamma
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:132

Here is the call graph for this function:

void G4DNABornIonisationModel1::SampleSecondaries ( std::vector< G4DynamicParticle * > *  fvect,
const G4MaterialCutsCouple *  couple,
const G4DynamicParticle *  particle,
G4double  tmin,
G4double  maxEnergy 
)
virtual

Implements G4VEmModel.

Definition at line 410 of file G4DNABornIonisationModel1.cc.

415 {
416 
417  if (verboseLevel > 3)
418  {
419  G4cout << "Calling SampleSecondaries() of G4DNABornIonisationModel1"
420  << G4endl;
421  }
422 
423  G4double lowLim = 0;
424  G4double highLim = 0;
425 
426  G4double k = particle->GetKineticEnergy();
427 
428  const G4String& particleName = particle->GetDefinition()->GetParticleName();
429 
430  std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
431  pos1 = lowEnergyLimit.find(particleName);
432 
433  if (pos1 != lowEnergyLimit.end())
434  {
435  lowLim = pos1->second;
436  }
437 
438  std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
439  pos2 = highEnergyLimit.find(particleName);
440 
441  if (pos2 != highEnergyLimit.end())
442  {
443  highLim = pos2->second;
444  }
445 
446  if (k >= lowLim && k < highLim)
447  {
448  G4ParticleMomentum primaryDirection = particle->GetMomentumDirection();
449  G4double particleMass = particle->GetDefinition()->GetPDGMass();
450  G4double totalEnergy = k + particleMass;
451  G4double pSquare = k * (totalEnergy + particleMass);
452  G4double totalMomentum = std::sqrt(pSquare);
453 
454  G4int ionizationShell = 0;
455 
456  if (!fasterCode) ionizationShell = RandomSelect(k,particleName);
457 
458  // SI: The following protection is necessary to avoid infinite loops :
459  // sigmadiff_ionisation_e_born.dat has non zero partial xs at 18 eV for shell 3 (ionizationShell ==2)
460  // sigmadiff_cumulated_ionisation_e_born.dat has zero cumulated partial xs at 18 eV for shell 3 (ionizationShell ==2)
461  // this is due to the fact that the max allowed transfered energy is (18+10.79)/2=17.025 eV and only transfered energies
462  // strictly above this value have non zero partial xs in sigmadiff_ionisation_e_born.dat (starting at trans = 17.12 eV)
463 
464  if (fasterCode)
465  do
466  {
467  ionizationShell = RandomSelect(k,particleName);
468  }while (k<19*eV && ionizationShell==2 && particle->GetDefinition()==G4Electron::ElectronDefinition());
469 
470  // AM: sample deexcitation
471  // here we assume that H_{2}O electronic levels are the same as Oxygen.
472  // this can be considered true with a rough 10% error in energy on K-shell,
473 
474  G4int secNumberInit = 0;// need to know at a certain point the energy of secondaries
475  G4int secNumberFinal = 0;// So I'll make the diference and then sum the energies
476 
477  G4double bindingEnergy = 0;
478  bindingEnergy = waterStructure.IonisationEnergy(ionizationShell);
479 
480  //SI: additional protection if tcs interpolation method is modified
481  if (k<bindingEnergy) return;
482  //
483 
484  G4int Z = 8;
485  if(fAtomDeexcitation)
486  {
487  G4AtomicShellEnumerator as = fKShell;
488 
489  if (ionizationShell <5 && ionizationShell >1)
490  {
491  as = G4AtomicShellEnumerator(4-ionizationShell);
492  }
493  else if (ionizationShell <2)
494  {
495  as = G4AtomicShellEnumerator(3);
496  }
497 
498  // FOR DEBUG ONLY
499  // if (ionizationShell == 4) {
500  //
501  // G4cout << "Z: " << Z << " as: " << as
502  // << " ionizationShell: " << ionizationShell << " bindingEnergy: "<< bindingEnergy/eV << G4endl;
503  // G4cout << "Press <Enter> key to continue..." << G4endl;
504  // G4cin.ignore();
505  // }
506 
507  const G4AtomicShell* shell = fAtomDeexcitation->GetAtomicShell(Z, as);
508  secNumberInit = fvect->size();
509  fAtomDeexcitation->GenerateParticles(fvect, shell, Z, 0, 0);
510  secNumberFinal = fvect->size();
511  }
512 
513  G4double secondaryKinetic=-1000*eV;
514 
515  if (fasterCode == false)
516  {
517  secondaryKinetic = RandomizeEjectedElectronEnergy(particle->GetDefinition(),k,ionizationShell);
518  }
519  // SI - 01/04/2014
520  else
521  {
522  secondaryKinetic = RandomizeEjectedElectronEnergyFromCumulatedDcs(particle->GetDefinition(),k,ionizationShell);
523  }
524  //
525 
526  G4ThreeVector deltaDirection =
527  GetAngularDistribution()->SampleDirectionForShell(particle, secondaryKinetic,
528  Z, ionizationShell,
529  couple->GetMaterial());
530 
531  if (particle->GetDefinition() == G4Electron::ElectronDefinition())
532  {
533  G4double deltaTotalMomentum = std::sqrt(secondaryKinetic*(secondaryKinetic + 2.*electron_mass_c2 ));
534 
535  G4double finalPx = totalMomentum*primaryDirection.x() - deltaTotalMomentum*deltaDirection.x();
536  G4double finalPy = totalMomentum*primaryDirection.y() - deltaTotalMomentum*deltaDirection.y();
537  G4double finalPz = totalMomentum*primaryDirection.z() - deltaTotalMomentum*deltaDirection.z();
538  G4double finalMomentum = std::sqrt(finalPx*finalPx + finalPy*finalPy + finalPz*finalPz);
539  finalPx /= finalMomentum;
540  finalPy /= finalMomentum;
541  finalPz /= finalMomentum;
542 
543  G4ThreeVector direction;
544  direction.set(finalPx,finalPy,finalPz);
545 
546  fParticleChangeForGamma->ProposeMomentumDirection(direction.unit());
547  }
548 
549  else fParticleChangeForGamma->ProposeMomentumDirection(primaryDirection);
550 
551  // note that secondaryKinetic is the energy of the delta ray, not of all secondaries.
552  G4double scatteredEnergy = k-bindingEnergy-secondaryKinetic;
553  G4double deexSecEnergy = 0;
554  for (G4int j=secNumberInit; j < secNumberFinal; j++)
555  {
556  deexSecEnergy = deexSecEnergy + (*fvect)[j]->GetKineticEnergy();
557  }
558 
559  if (!statCode)
560  {
561  fParticleChangeForGamma->SetProposedKineticEnergy(scatteredEnergy);
562  fParticleChangeForGamma->ProposeLocalEnergyDeposit(k-scatteredEnergy-secondaryKinetic-deexSecEnergy);
563  }
564  else
565  {
566  fParticleChangeForGamma->SetProposedKineticEnergy(k);
567  fParticleChangeForGamma->ProposeLocalEnergyDeposit(k-scatteredEnergy);
568  }
569 
570  // SI - 01/04/2014
571  if (secondaryKinetic>0)
572  {
573  G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(),deltaDirection,secondaryKinetic);
574  fvect->push_back(dp);
575  }
576  //
577 
578  const G4Track * theIncomingTrack = fParticleChangeForGamma->GetCurrentTrack();
579  G4DNAChemistryManager::Instance()->CreateWaterMolecule(eIonizedMolecule,
580  ionizationShell,
581  theIncomingTrack);
582  }
583 }
CLHEP::Hep3Vector::set
void set(double x, double y, double z)
G4Electron::ElectronDefinition
static G4Electron * ElectronDefinition()
Definition: G4Electron.cc:89
G4DynamicParticle::GetKineticEnergy
G4double GetKineticEnergy() const
CLHEP::Hep3Vector::x
double x() const
G4VEmModel::GetAngularDistribution
G4VEmAngularDistribution * GetAngularDistribution()
Definition: G4VEmModel.hh:619
G4DynamicParticle::GetDefinition
G4ParticleDefinition * GetDefinition() const
G4int
int G4int
Definition: G4Types.hh:78
G4ParticleChangeForGamma::ProposeMomentumDirection
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
Definition: G4ParticleChangeForGamma.hh:147
fKShell
Definition: G4AtomicShellEnumerator.hh:53
G4ParticleDefinition::GetParticleName
const G4String & GetParticleName() const
Definition: G4ParticleDefinition.hh:120
CLHEP::Hep3Vector::z
double z() const
G4VParticleChange::ProposeLocalEnergyDeposit
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
G4VEmAngularDistribution::SampleDirectionForShell
virtual G4ThreeVector & SampleDirectionForShell(const G4DynamicParticle *dp, G4double finalTotalEnergy, G4int Z, G4int shellID, const G4Material *)
Definition: G4VEmAngularDistribution.cc:60
G4VAtomDeexcitation::GetAtomicShell
virtual const G4AtomicShell * GetAtomicShell(G4int Z, G4AtomicShellEnumerator shell)=0
G4DNABornIonisationModel1::fParticleChangeForGamma
G4ParticleChangeForGamma * fParticleChangeForGamma
Definition: G4DNABornIonisationModel1.hh:90
G4cout
G4GLOB_DLL std::ostream G4cout
G4DynamicParticle::GetMomentumDirection
const G4ThreeVector & GetMomentumDirection() const
eV
static constexpr double eV
Definition: G4SIunits.hh:215
python.hepunit.electron_mass_c2
float electron_mass_c2
Definition: hepunit.py:274
G4DNAChemistryManager::Instance
static G4DNAChemistryManager * Instance()
Definition: G4DNAChemistryManager.cc:117
eIonizedMolecule
Definition: G4DNAChemistryManager.hh:70
G4ParticleDefinition::GetPDGMass
G4double GetPDGMass() const
Definition: G4ParticleDefinition.hh:122
G4DNAChemistryManager::CreateWaterMolecule
void CreateWaterMolecule(ElectronicModification, G4int, const G4Track *)
Definition: G4DNAChemistryManager.cc:572
CLHEP::Hep3Vector::unit
Hep3Vector unit() const
CLHEP::Hep3Vector::y
double y() const
G4ParticleChangeForGamma::GetCurrentTrack
const G4Track * GetCurrentTrack() const
Definition: G4ParticleChangeForGamma.hh:154
G4Electron::Electron
static G4Electron * Electron()
Definition: G4Electron.cc:94
G4ParticleChangeForGamma::SetProposedKineticEnergy
void SetProposedKineticEnergy(G4double proposedKinEnergy)
Definition: G4ParticleChangeForGamma.hh:129
G4endl
#define G4endl
Definition: G4ios.hh:61
G4double
double G4double
Definition: G4Types.hh:76
G4VAtomDeexcitation::GenerateParticles
void GenerateParticles(std::vector< G4DynamicParticle * > *secVect, const G4AtomicShell *, G4int Z, G4int coupleIndex)
Definition: G4VAtomDeexcitation.hh:274
G4InuclSpecialFunctions::bindingEnergy
G4double bindingEnergy(G4int A, G4int Z)
Definition: bindingEnergy.cc:37
G4AtomicShellEnumerator
G4AtomicShellEnumerator
Definition: G4AtomicShellEnumerator.hh:50
Z
G4int Z
Definition: G4ParticleHPJENDLHEData.cc:262
G4MaterialCutsCouple::GetMaterial
const G4Material * GetMaterial() const
Definition: G4MaterialCutsCouple.hh:150

Here is the call graph for this function:

void G4DNABornIonisationModel1::SelectFasterComputation ( G4bool  input)
inline

Definition at line 175 of file G4DNABornIonisationModel1.hh.

176 {
177  fasterCode = input;
178 }
void G4DNABornIonisationModel1::SelectSPScaling ( G4bool  input)
inline

Definition at line 189 of file G4DNABornIonisationModel1.hh.

190 {
191  spScaling = input;
192 }
void G4DNABornIonisationModel1::SelectStationary ( G4bool  input)
inline

Definition at line 182 of file G4DNABornIonisationModel1.hh.

183 {
184  statCode = input;
185 }
G4double G4DNABornIonisationModel1::TransferedEnergy ( G4ParticleDefinition *  aParticleDefinition,
G4double  incomingParticleEnergy,
G4int  shell,
G4double  random 
)

Definition at line 1036 of file G4DNABornIonisationModel1.cc.

1040 {
1041  G4double nrj = 0.;
1042 
1043  G4double valueK1 = 0;
1044  G4double valueK2 = 0;
1045  G4double valuePROB21 = 0;
1046  G4double valuePROB22 = 0;
1047  G4double valuePROB12 = 0;
1048  G4double valuePROB11 = 0;
1049 
1050  G4double nrjTransf11 = 0;
1051  G4double nrjTransf12 = 0;
1052  G4double nrjTransf21 = 0;
1053  G4double nrjTransf22 = 0;
1054 
1055  if (particleDefinition == G4Electron::ElectronDefinition())
1056  {
1057  // k should be in eV
1058  std::vector<double>::iterator k2 = std::upper_bound(eTdummyVec.begin(),
1059  eTdummyVec.end(),
1060  k);
1061  std::vector<double>::iterator k1 = k2 - 1;
1062 
1063  /*
1064  G4cout << "----> k=" << k
1065  << " " << *k1
1066  << " " << *k2
1067  << " " << random
1068  << " " << ionizationLevelIndex
1069  << " " << eProbaShellMap[ionizationLevelIndex][(*k1)].back()
1070  << " " << eProbaShellMap[ionizationLevelIndex][(*k2)].back()
1071  << G4endl;
1072  */
1073 
1074  // SI : the following condition avoids situations where random >last vector element
1075  if (random <= eProbaShellMap[ionizationLevelIndex][(*k1)].back()
1076  && random <= eProbaShellMap[ionizationLevelIndex][(*k2)].back())
1077  {
1078  std::vector<double>::iterator prob12 =
1079  std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k1)].begin(),
1080  eProbaShellMap[ionizationLevelIndex][(*k1)].end(),
1081  random);
1082 
1083  std::vector<double>::iterator prob11 = prob12 - 1;
1084 
1085  std::vector<double>::iterator prob22 =
1086  std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
1087  eProbaShellMap[ionizationLevelIndex][(*k2)].end(),
1088  random);
1089 
1090  std::vector<double>::iterator prob21 = prob22 - 1;
1091 
1092  valueK1 = *k1;
1093  valueK2 = *k2;
1094  valuePROB21 = *prob21;
1095  valuePROB22 = *prob22;
1096  valuePROB12 = *prob12;
1097  valuePROB11 = *prob11;
1098 
1099  /*
1100  G4cout << " " << random << " " << valuePROB11 << " "
1101  << valuePROB12 << " " << valuePROB21 << " " << valuePROB22 << G4endl;
1102  */
1103 
1104  nrjTransf11 = eNrjTransfData[ionizationLevelIndex][valueK1][valuePROB11];
1105  nrjTransf12 = eNrjTransfData[ionizationLevelIndex][valueK1][valuePROB12];
1106  nrjTransf21 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
1107  nrjTransf22 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
1108 
1109  /*
1110  G4cout << " " << ionizationLevelIndex << " "
1111  << random << " " <<valueK1 << " " << valueK2 << G4endl;
1112 
1113  G4cout << " " << random << " " << nrjTransf11 << " "
1114  << nrjTransf12 << " " << nrjTransf21 << " " <<nrjTransf22 << G4endl;
1115  */
1116 
1117  }
1118  // Avoids cases where cum xs is zero for k1 and is not for k2 (with always k1<k2)
1119  if (random > eProbaShellMap[ionizationLevelIndex][(*k1)].back())
1120  {
1121  std::vector<double>::iterator prob22 =
1122  std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
1123  eProbaShellMap[ionizationLevelIndex][(*k2)].end(),
1124  random);
1125 
1126  std::vector<double>::iterator prob21 = prob22 - 1;
1127 
1128  valueK1 = *k1;
1129  valueK2 = *k2;
1130  valuePROB21 = *prob21;
1131  valuePROB22 = *prob22;
1132 
1133  //G4cout << " " << random << " " << valuePROB21 << " " << valuePROB22 << G4endl;
1134 
1135  nrjTransf21 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
1136  nrjTransf22 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
1137 
1138  G4double interpolatedvalue2 = Interpolate(valuePROB21,
1139  valuePROB22,
1140  random,
1141  nrjTransf21,
1142  nrjTransf22);
1143 
1144  // zeros are explicitely set
1145 
1146  G4double value = Interpolate(valueK1, valueK2, k, 0., interpolatedvalue2);
1147 
1148  /*
1149  G4cout << " " << ionizationLevelIndex << " "
1150  << random << " " <<valueK1 << " " << valueK2 << G4endl;
1151 
1152  G4cout << " " << random << " " << nrjTransf11 << " "
1153  << nrjTransf12 << " " << nrjTransf21 << " " <<nrjTransf22 << G4endl;
1154 
1155  G4cout << "ici" << " " << value << G4endl;
1156  */
1157 
1158  return value;
1159  }
1160  }
1161  //
1162  else if (particleDefinition == G4Proton::ProtonDefinition())
1163  {
1164  // k should be in eV
1165 
1166  std::vector<double>::iterator k2 = std::upper_bound(pTdummyVec.begin(),
1167  pTdummyVec.end(),
1168  k);
1169 
1170  std::vector<double>::iterator k1 = k2 - 1;
1171 
1172  /*
1173  G4cout << "----> k=" << k
1174  << " " << *k1
1175  << " " << *k2
1176  << " " << random
1177  << " " << ionizationLevelIndex
1178  << " " << pProbaShellMap[ionizationLevelIndex][(*k1)].back()
1179  << " " << pProbaShellMap[ionizationLevelIndex][(*k2)].back()
1180  << G4endl;
1181  */
1182 
1183  // SI : the following condition avoids situations where random > last vector element,
1184  // for eg. when the last element is zero
1185  if (random <= pProbaShellMap[ionizationLevelIndex][(*k1)].back()
1186  && random <= pProbaShellMap[ionizationLevelIndex][(*k2)].back())
1187  {
1188  std::vector<double>::iterator prob12 =
1189  std::upper_bound(pProbaShellMap[ionizationLevelIndex][(*k1)].begin(),
1190  pProbaShellMap[ionizationLevelIndex][(*k1)].end(),
1191  random);
1192 
1193  std::vector<double>::iterator prob11 = prob12 - 1;
1194 
1195  std::vector<double>::iterator prob22 =
1196  std::upper_bound(pProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
1197  pProbaShellMap[ionizationLevelIndex][(*k2)].end(),
1198  random);
1199 
1200  std::vector<double>::iterator prob21 = prob22 - 1;
1201 
1202  valueK1 = *k1;
1203  valueK2 = *k2;
1204  valuePROB21 = *prob21;
1205  valuePROB22 = *prob22;
1206  valuePROB12 = *prob12;
1207  valuePROB11 = *prob11;
1208 
1209  /*
1210  G4cout << " " << random << " " << valuePROB11 << " "
1211  << valuePROB12 << " " << valuePROB21 << " " << valuePROB22 << G4endl;
1212  */
1213 
1214  nrjTransf11 = pNrjTransfData[ionizationLevelIndex][valueK1][valuePROB11];
1215  nrjTransf12 = pNrjTransfData[ionizationLevelIndex][valueK1][valuePROB12];
1216  nrjTransf21 = pNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
1217  nrjTransf22 = pNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
1218 
1219  /*
1220  G4cout << " " << ionizationLevelIndex << " "
1221  << random << " " <<valueK1 << " " << valueK2 << G4endl;
1222 
1223  G4cout << " " << random << " " << nrjTransf11 << " "
1224  << nrjTransf12 << " " << nrjTransf21 << " " <<nrjTransf22 << G4endl;
1225  */
1226  }
1227 
1228  // Avoids cases where cum xs is zero for k1 and is not for k2 (with always k1<k2)
1229 
1230  if (random > pProbaShellMap[ionizationLevelIndex][(*k1)].back())
1231  {
1232  std::vector<double>::iterator prob22 =
1233  std::upper_bound(pProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
1234  pProbaShellMap[ionizationLevelIndex][(*k2)].end(),
1235  random);
1236 
1237  std::vector<double>::iterator prob21 = prob22 - 1;
1238 
1239  valueK1 = *k1;
1240  valueK2 = *k2;
1241  valuePROB21 = *prob21;
1242  valuePROB22 = *prob22;
1243 
1244  //G4cout << " " << random << " " << valuePROB21 << " " << valuePROB22 << G4endl;
1245 
1246  nrjTransf21 = pNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
1247  nrjTransf22 = pNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
1248 
1249  G4double interpolatedvalue2 = Interpolate(valuePROB21,
1250  valuePROB22,
1251  random,
1252  nrjTransf21,
1253  nrjTransf22);
1254 
1255  // zeros are explicitely set
1256 
1257  G4double value = Interpolate(valueK1, valueK2, k, 0., interpolatedvalue2);
1258 
1259  /*
1260  G4cout << " " << ionizationLevelIndex << " "
1261  << random << " " <<valueK1 << " " << valueK2 << G4endl;
1262 
1263  G4cout << " " << random << " " << nrjTransf11 << " "
1264  << nrjTransf12 << " " << nrjTransf21 << " " <<nrjTransf22 << G4endl;
1265 
1266  G4cout << "ici" << " " << value << G4endl;
1267  */
1268 
1269  return value;
1270  }
1271  }
1272  // End electron and proton cases
1273 
1274  G4double nrjTransfProduct = nrjTransf11 * nrjTransf12 * nrjTransf21
1275  * nrjTransf22;
1276  //G4cout << "nrjTransfProduct=" << nrjTransfProduct << G4endl;
1277 
1278  if (nrjTransfProduct != 0.)
1279  {
1280  nrj = QuadInterpolator(valuePROB11,
1281  valuePROB12,
1282  valuePROB21,
1283  valuePROB22,
1284  nrjTransf11,
1285  nrjTransf12,
1286  nrjTransf21,
1287  nrjTransf22,
1288  valueK1,
1289  valueK2,
1290  k,
1291  random);
1292  }
1293  //G4cout << nrj << endl;
1294 
1295  return nrj;
1296 }
G4Electron::ElectronDefinition
static G4Electron * ElectronDefinition()
Definition: G4Electron.cc:89
G4Proton::ProtonDefinition
static G4Proton * ProtonDefinition()
Definition: G4Proton.cc:88
value
const XML_Char int const XML_Char * value
Definition: expat.h:331
G4double
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

Member Data Documentation

G4ParticleChangeForGamma* G4DNABornIonisationModel1::fParticleChangeForGamma
protected

Definition at line 90 of file G4DNABornIonisationModel1.hh.

The documentation for this class was generated from the following files: