Geant4  10.02.p03
G4DNAIonElasticModel Class Reference

#include <G4DNAIonElasticModel.hh>

Inheritance diagram for G4DNAIonElasticModel:
Collaboration diagram for G4DNAIonElasticModel:

Public Member Functions

 G4DNAIonElasticModel (const G4ParticleDefinition *p=0, const G4String &nam="DNAIonElasticModel")
 
virtual ~G4DNAIonElasticModel ()
 
virtual void Initialise (const G4ParticleDefinition *particuleDefinition, const 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)
 
void SetKillBelowThreshold (G4double threshold)
 
G4double GetKillBelowThreshold ()
 
void SelectStationary (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 GetPartialCrossSection (const G4Material *, G4int, const G4ParticleDefinition *, G4double)
 
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 *> *)
 
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 G4ElementSelectRandomAtom (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
const G4ElementSelectRandomAtom (const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
G4int SelectRandomAtomNumber (const G4Material *)
 
void SetParticleChange (G4VParticleChange *, G4VEmFluctuationModel *f=0)
 
void SetCrossSectionTable (G4PhysicsTable *, G4bool isLocal)
 
G4ElementDataGetElementData ()
 
G4PhysicsTableGetCrossSectionTable ()
 
G4VEmFluctuationModelGetModelOfFluctuations ()
 
G4VEmAngularDistributionGetAngularDistribution ()
 
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 SetMasterThread (G4bool val)
 
G4bool IsMaster () const
 
G4double MaxSecondaryKinEnergy (const G4DynamicParticle *dynParticle)
 
const G4StringGetName () const
 
void SetCurrentCouple (const G4MaterialCutsCouple *)
 
const G4ElementGetCurrentElement () const
 
const G4IsotopeGetCurrentIsotope () const
 
G4bool IsLocked () const
 
void SetLocked (G4bool)
 

Protected Attributes

G4ParticleChangeForGamma * fParticleChangeForGamma
 
- Protected Attributes inherited from G4VEmModel
G4ElementDatafElementData
 
G4VParticleChange * pParticleChange
 
G4PhysicsTablexSectionTable
 
const std::vector< G4double > * theDensityFactor
 
const std::vector< G4int > * theDensityIdx
 
size_t idxTable
 

Private Types

typedef std::map< double, std::map< double, double > > TriDimensionMap
 
typedef std::map< double, std::vector< double > > VecMap
 

Private Member Functions

G4double Theta (G4ParticleDefinition *aParticleDefinition, G4double k, G4double integrDiff)
 
G4double LinLogInterpolate (G4double e1, G4double e2, G4double e, G4double xs1, G4double xs2)
 
G4double LogLogInterpolate (G4double e1, G4double e2, G4double e, G4double xs1, G4double xs2)
 
G4double QuadInterpolator (G4double e11, G4double e12, G4double e21, G4double e22, G4double x11, G4double x12, G4double x21, G4double x22, G4double t1, G4double t2, G4double t, G4double e)
 
G4double LinLinInterpolate (G4double e1, G4double e2, G4double e, G4double xs1, G4double xs2)
 
G4double RandomizeThetaCM (G4double k, G4ParticleDefinition *aParticleDefinition)
 
G4DNAIonElasticModeloperator= (const G4DNAIonElasticModel &right)
 
 G4DNAIonElasticModel (const G4DNAIonElasticModel &)
 

Private Attributes

G4bool statCode
 
const std::vector< G4double > * fpMolWaterDensity
 
G4double killBelowEnergy
 
G4double lowEnergyLimit
 
G4double highEnergyLimit
 
G4bool isInitialised
 
G4int verboseLevel
 
G4double fParticle_Mass
 
G4DNACrossSectionDataSetfpTableData
 
TriDimensionMap fDiffCrossSectionData
 
std::vector< double > eTdummyVec
 
VecMap eVecm
 

Additional Inherited Members

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

Detailed Description

Definition at line 46 of file G4DNAIonElasticModel.hh.

Member Typedef Documentation

◆ TriDimensionMap

typedef std::map<double, std::map<double, double> > G4DNAIonElasticModel::TriDimensionMap
private

Definition at line 127 of file G4DNAIonElasticModel.hh.

◆ VecMap

typedef std::map<double, std::vector<double> > G4DNAIonElasticModel::VecMap
private

Definition at line 132 of file G4DNAIonElasticModel.hh.

Constructor & Destructor Documentation

◆ G4DNAIonElasticModel() [1/2]

G4DNAIonElasticModel::G4DNAIonElasticModel ( const G4ParticleDefinition p = 0,
const G4String nam = "DNAIonElasticModel" 
)

Definition at line 45 of file G4DNAIonElasticModel.cc.

46  :
47  G4VEmModel(nam), isInitialised(false)
48 {
49  //nistwater = G4NistManager::Instance()->FindOrBuildMaterial("G4_WATER");
50 
51  killBelowEnergy = 100 * eV;
52  lowEnergyLimit = 0 * eV;
53  highEnergyLimit = 1 * MeV;
56 
57  verboseLevel = 0;
58  // Verbosity scale:
59  // 0 = nothing
60  // 1 = warning for energy non-conservation
61  // 2 = details of energy budget
62  // 3 = calculation of cross sections, file openings, sampling of atoms
63  // 4 = entering in methods
64 
65  if(verboseLevel > 0)
66  {
67  G4cout << "Ion elastic model is constructed " << G4endl<< "Energy range: "
68  << lowEnergyLimit / eV << " eV - "
69  << highEnergyLimit / MeV << " MeV"
70  << G4endl;
71  }
74  fpTableData = 0;
75  fParticle_Mass = -1;
76 
77  // Selection of stationary mode
78 
79  statCode = false;
80 }
static const double MeV
Definition: G4SIunits.hh:211
G4VEmModel(const G4String &nam)
Definition: G4VEmModel.cc:69
void SetHighEnergyLimit(G4double)
Definition: G4VEmModel.hh:725
G4GLOB_DLL std::ostream G4cout
G4DNACrossSectionDataSet * fpTableData
G4ParticleChangeForGamma * fParticleChangeForGamma
const std::vector< G4double > * fpMolWaterDensity
static const double eV
Definition: G4SIunits.hh:212
#define G4endl
Definition: G4ios.hh:61
void SetLowEnergyLimit(G4double)
Definition: G4VEmModel.hh:732
Here is the call graph for this function:

◆ ~G4DNAIonElasticModel()

G4DNAIonElasticModel::~G4DNAIonElasticModel ( )
virtual

Definition at line 84 of file G4DNAIonElasticModel.cc.

85 {
86  // For total cross section
87  if(fpTableData) delete fpTableData;
88 }
G4DNACrossSectionDataSet * fpTableData

◆ G4DNAIonElasticModel() [2/2]

G4DNAIonElasticModel::G4DNAIonElasticModel ( const G4DNAIonElasticModel )
private

Member Function Documentation

◆ CrossSectionPerVolume()

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

Reimplemented from G4VEmModel.

Definition at line 247 of file G4DNAIonElasticModel.cc.

250 {
251  if(verboseLevel > 3)
252  {
253  G4cout << "Calling CrossSectionPerVolume() of G4DNAIonElasticModel"
254  << G4endl;
255  }
256 
257  // Calculate total cross section for model
258 
259  G4double sigma=0;
260 
261  G4double waterDensity = (*fpMolWaterDensity)[material->GetIndex()];
262 
263  if(waterDensity!= 0.0)
264  {
265  const G4String& particleName = p->GetParticleName();
266 
267  if (ekin < highEnergyLimit)
268  {
269  //SI : XS must not be zero otherwise sampling of secondaries method ignored
270  if (ekin < killBelowEnergy) return DBL_MAX;
271  //
272 
273  if (fpTableData != 0) // MK: is this check necessary?
274  {
275  sigma = fpTableData->FindValue(ekin);
276  }
277  else
278  {
279  G4Exception("G4DNAIonElasticModel::ComputeCrossSectionPerVolume","em0002",
280  FatalException,"Model not applicable to particle type.");
281  }
282  }
283 
284  if (verboseLevel > 2)
285  {
286  G4cout << "__________________________________" << G4endl;
287  G4cout << "G4DNAIonElasticModel - XS INFO START" << G4endl;
288  G4cout << "Kinetic energy(eV)=" << ekin/eV << " particle : " << particleName << G4endl;
289  G4cout << "Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
290  G4cout << "Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) << G4endl;
291  G4cout << "G4DNAIonElasticModel - XS INFO END" << G4endl;
292  }
293 
294  }
295 
296  return sigma*waterDensity;
297 }
static const double cm
Definition: G4SIunits.hh:118
size_t GetIndex() const
Definition: G4Material.hh:262
const G4String & GetParticleName() const
virtual G4double FindValue(G4double e, G4int componentId=0) const
G4GLOB_DLL std::ostream G4cout
G4DNACrossSectionDataSet * fpTableData
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
static const double eV
Definition: G4SIunits.hh:212
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
#define DBL_MAX
Definition: templates.hh:83
Here is the call graph for this function:

◆ GetKillBelowThreshold()

G4double G4DNAIonElasticModel::GetKillBelowThreshold ( )
inline

Definition at line 75 of file G4DNAIonElasticModel.hh.

76  {
77  return killBelowEnergy;
78  }
Here is the call graph for this function:

◆ Initialise()

void G4DNAIonElasticModel::Initialise ( const G4ParticleDefinition particuleDefinition,
const G4DataVector  
)
virtual

Implements G4VEmModel.

Definition at line 93 of file G4DNAIonElasticModel.cc.

96 {
97 
98  if(verboseLevel > 3)
99  {
100  G4cout << "Calling G4DNAIonElasticModel::Initialise()" << G4endl;
101  }
102 
103  // Energy limits
104 
106  {
107  G4cout << "G4DNAIonElasticModel: low energy limit increased from " <<
108  LowEnergyLimit()/eV << " eV to " << lowEnergyLimit/eV << " eV" << G4endl;
110  }
111 
113  {
114  G4cout << "G4DNAIonElasticModel: high energy limit decreased from " <<
115  HighEnergyLimit()/MeV << " MeV to " << highEnergyLimit/MeV << " MeV" << G4endl;
117  }
118 
119  // Reading of data files
120 
121  G4double scaleFactor = 1e-16*cm*cm;
122 
123  char *path = getenv("G4LEDATA");
124 
125  if (!path)
126  {
127  G4Exception("G4IonElasticModel::Initialise","em0006",
128  FatalException,"G4LEDATA environment variable not set.");
129  return;
130  }
131 
132  G4String totalXSFile;
133  std::ostringstream fullFileName;
134 
137  G4ParticleDefinition* protonDef =
139  G4ParticleDefinition* hydrogenDef = instance->GetIon("hydrogen");
140  G4ParticleDefinition* heliumDef = instance->GetIon("helium");
141  G4ParticleDefinition* alphaplusDef = instance->GetIon("alpha+");
142  G4ParticleDefinition* alphaplusplusDef = instance->GetIon("alpha++");
143  G4String proton, hydrogen, helium, alphaplus, alphaplusplus;
144 
145  if (
146  (particleDefinition == protonDef && protonDef != 0)
147  ||
148  (particleDefinition == hydrogenDef && hydrogenDef != 0)
149  )
150  {
151  // For total cross section of p,h
152  fParticle_Mass = 1.;
153  totalXSFile = "dna/sigma_elastic_proton_HTran";
154 
155  // For final state
156  fullFileName << path << "/dna/sigmadiff_cumulated_elastic_proton_HTran.dat";
157  }
158 
159  if (
160  (particleDefinition == instance->GetIon("helium") && heliumDef)
161  ||
162  (particleDefinition == instance->GetIon("alpha+") && alphaplusDef)
163  ||
164  (particleDefinition == instance->GetIon("alpha++") && alphaplusplusDef)
165  )
166  {
167  fParticle_Mass = 4.;
168  // For total cross section of he,he+,he++
169  totalXSFile = "dna/sigma_elastic_alpha_HTran";
170 
171  // For final state
172  fullFileName << path << "/dna/sigmadiff_cumulated_elastic_alpha_HTran.dat";
173  }
174 
176  fpTableData->LoadData(totalXSFile);
177  std::ifstream diffCrossSection(fullFileName.str().c_str());
178 
179  if (!diffCrossSection)
180  {
181  G4ExceptionDescription description;
182  description << "Missing data file:"
183  <<fullFileName.str().c_str()<< G4endl;
184  G4Exception("G4IonElasticModel::Initialise","em0003",
186  description);
187  }
188 
189  // Added clear for MT
190 
191  eTdummyVec.clear();
192  eVecm.clear();
193  fDiffCrossSectionData.clear();
194 
195  //
196 
197  eTdummyVec.push_back(0.);
198 
199  while(!diffCrossSection.eof())
200  {
201  double tDummy;
202  double eDummy;
203  diffCrossSection>>tDummy>>eDummy;
204 
205  // SI : mandatory eVecm initialization
206 
207  if (tDummy != eTdummyVec.back())
208  {
209  eTdummyVec.push_back(tDummy);
210  eVecm[tDummy].push_back(0.);
211  }
212 
213  diffCrossSection>>fDiffCrossSectionData[tDummy][eDummy];
214 
215  if (eDummy != eVecm[tDummy].back()) eVecm[tDummy].push_back(eDummy);
216 
217  }
218 
219  // End final state
220  if( verboseLevel>0 )
221  {
222  if (verboseLevel > 2)
223  {
224  G4cout << "Loaded cross section files for ion elastic model" << G4endl;
225  }
226  G4cout << "Ion elastic model is initialized " << G4endl
227  << "Energy range: "
228  << LowEnergyLimit() / eV << " eV - "
229  << HighEnergyLimit() / MeV << " MeV"
230  << G4endl;
231  }
232 
233  // Initialize water density pointer
236  GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
237 
238  if (isInitialised)
239  { return;}
241  isInitialised = true;
242 }
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:641
static const double cm
Definition: G4SIunits.hh:118
static const double MeV
Definition: G4SIunits.hh:211
G4ParticleDefinition * FindParticle(G4int PDGEncoding)
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
static G4Material * GetMaterial(const G4String &name, G4bool warning=true)
Definition: G4Material.cc:604
TriDimensionMap fDiffCrossSectionData
std::vector< double > eTdummyVec
virtual G4bool LoadData(const G4String &argFileName)
void SetHighEnergyLimit(G4double)
Definition: G4VEmModel.hh:725
G4GLOB_DLL std::ostream G4cout
G4DNACrossSectionDataSet * fpTableData
G4ParticleChangeForGamma * fParticleChangeForGamma
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:634
const std::vector< G4double > * fpMolWaterDensity
static G4DNAGenericIonsManager * Instance(void)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
static const double eV
Definition: G4SIunits.hh:212
static G4DNAMolecularMaterial * Instance()
static G4ParticleTable * GetParticleTable()
#define G4endl
Definition: G4ios.hh:61
static MCTruthManager * instance
double G4double
Definition: G4Types.hh:76
void SetLowEnergyLimit(G4double)
Definition: G4VEmModel.hh:732
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:134
G4ParticleDefinition * GetIon(const G4String &name)
Here is the call graph for this function:

◆ LinLinInterpolate()

G4double G4DNAIonElasticModel::LinLinInterpolate ( G4double  e1,
G4double  e2,
G4double  e,
G4double  xs1,
G4double  xs2 
)
private

Definition at line 426 of file G4DNAIonElasticModel.cc.

428 {
429  G4double d1 = xs1;
430  G4double d2 = xs2;
431  G4double value = (d1 + (d2 - d1) * (e - e1) / (e2 - e1));
432  return value;
433 }
static const G4double d1
static const G4double e2
static const G4double e1
double G4double
Definition: G4Types.hh:76
static const G4double d2
Here is the caller graph for this function:

◆ LinLogInterpolate()

G4double G4DNAIonElasticModel::LinLogInterpolate ( G4double  e1,
G4double  e2,
G4double  e,
G4double  xs1,
G4double  xs2 
)
private

Definition at line 438 of file G4DNAIonElasticModel.cc.

440 {
441  G4double d1 = std::log(xs1);
442  G4double d2 = std::log(xs2);
443  G4double value = std::exp(d1 + (d2 - d1) * (e - e1) / (e2 - e1));
444  return value;
445 }
static const G4double d1
static const G4double e2
static const G4double e1
double G4double
Definition: G4Types.hh:76
static const G4double d2

◆ LogLogInterpolate()

G4double G4DNAIonElasticModel::LogLogInterpolate ( G4double  e1,
G4double  e2,
G4double  e,
G4double  xs1,
G4double  xs2 
)
private

Definition at line 450 of file G4DNAIonElasticModel.cc.

452 {
453  G4double a = (std::log10(xs2) - std::log10(xs1))
454  / (std::log10(e2) - std::log10(e1));
455  G4double b = std::log10(xs2) - a * std::log10(e2);
456  G4double sigma = a * std::log10(e) + b;
457  G4double value = (std::pow(10., sigma));
458  return value;
459 }
static const G4double e2
static const G4double e1
double G4double
Definition: G4Types.hh:76

◆ operator=()

G4DNAIonElasticModel& G4DNAIonElasticModel::operator= ( const G4DNAIonElasticModel right)
private

◆ QuadInterpolator()

G4double G4DNAIonElasticModel::QuadInterpolator ( G4double  e11,
G4double  e12,
G4double  e21,
G4double  e22,
G4double  x11,
G4double  x12,
G4double  x21,
G4double  x22,
G4double  t1,
G4double  t2,
G4double  t,
G4double  e 
)
private

Definition at line 464 of file G4DNAIonElasticModel.cc.

470 {
471  // Log-Log
472  /*
473  G4double interpolatedvalue1 = LogLogInterpolate(e11, e12, e, xs11, xs12);
474  G4double interpolatedvalue2 = LogLogInterpolate(e21, e22, e, xs21, xs22);
475  G4double value = LogLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
476  */
477 
478  // Lin-Log
479  /*
480  G4double interpolatedvalue1 = LinLogInterpolate(e11, e12, e, xs11, xs12);
481  G4double interpolatedvalue2 = LinLogInterpolate(e21, e22, e, xs21, xs22);
482  G4double value = LinLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
483  */
484 
485 // Lin-Lin
486  G4double interpolatedvalue1 = LinLinInterpolate(e11, e12, e, xs11, xs12);
487  G4double interpolatedvalue2 = LinLinInterpolate(e21, e22, e, xs21, xs22);
488  G4double value = LinLinInterpolate(t1, t2, t, interpolatedvalue1,
489  interpolatedvalue2);
490 
491  return value;
492 }
TTree * t1
Definition: plottest35.C:26
G4double LinLinInterpolate(G4double e1, G4double e2, G4double e, G4double xs1, G4double xs2)
TTree * t2
Definition: plottest35.C:36
double G4double
Definition: G4Types.hh:76
Here is the call graph for this function:
Here is the caller graph for this function:

◆ RandomizeThetaCM()

G4double G4DNAIonElasticModel::RandomizeThetaCM ( G4double  k,
G4ParticleDefinition aParticleDefinition 
)
private

Definition at line 497 of file G4DNAIonElasticModel.cc.

499 {
500  G4double integrdiff = G4UniformRand();
501  return Theta(particleDefinition, k / eV, integrdiff);
502 }
#define G4UniformRand()
Definition: Randomize.hh:97
static const double eV
Definition: G4SIunits.hh:212
double G4double
Definition: G4Types.hh:76
G4double Theta(G4ParticleDefinition *aParticleDefinition, G4double k, G4double integrDiff)
Here is the call graph for this function:
Here is the caller graph for this function:

◆ SampleSecondaries()

void G4DNAIonElasticModel::SampleSecondaries ( std::vector< G4DynamicParticle *> *  ,
const G4MaterialCutsCouple ,
const G4DynamicParticle aDynamicParticle,
G4double  tmin,
G4double  maxEnergy 
)
virtual

Implements G4VEmModel.

Definition at line 302 of file G4DNAIonElasticModel.cc.

306 {
307 
308  if(verboseLevel > 3)
309  {
310  G4cout << "Calling SampleSecondaries() of G4DNAIonElasticModel" << G4endl;
311  }
312 
313  G4double particleEnergy0 = aDynamicParticle->GetKineticEnergy();
314 
315  if (particleEnergy0 < killBelowEnergy)
316  {
317  fParticleChangeForGamma->SetProposedKineticEnergy(0.);
318  fParticleChangeForGamma->ProposeTrackStatus(fStopAndKill);
319  fParticleChangeForGamma->ProposeLocalEnergyDeposit(particleEnergy0);
320  return;
321  }
322 
323  if (particleEnergy0>= killBelowEnergy && particleEnergy0 < highEnergyLimit)
324  {
325  G4double water_mass = 18.;
326 
327  G4double thetaCM = RandomizeThetaCM(particleEnergy0, aDynamicParticle->GetDefinition());
328 
329  //HT:convert to laboratory system
330 
331  G4double theta = std::atan(std::sin(thetaCM*CLHEP::pi/180)
332  /(fParticle_Mass/water_mass+std::cos(thetaCM*CLHEP::pi/180)));
333 
334  G4double cosTheta= std::cos(theta);
335 
336  //
337 
338  G4double phi = 2. * CLHEP::pi * G4UniformRand();
339 
340  G4ThreeVector zVers = aDynamicParticle->GetMomentumDirection();
341  G4ThreeVector xVers = zVers.orthogonal();
342  G4ThreeVector yVers = zVers.cross(xVers);
343 
344  G4double xDir = std::sqrt(1. - cosTheta*cosTheta);
345  G4double yDir = xDir;
346  xDir *= std::cos(phi);
347  yDir *= std::sin(phi);
348 
349  G4ThreeVector zPrimeVers((xDir*xVers + yDir*yVers + cosTheta*zVers));
350 
351  fParticleChangeForGamma->ProposeMomentumDirection(zPrimeVers.unit());
352 
353  G4double depositEnergyCM = 0;
354 
355  //HT: deposited energy
356  depositEnergyCM = 4. * particleEnergy0 * fParticle_Mass * water_mass *
357  (1-std::cos(thetaCM*CLHEP::pi/180))
358  / (2 * std::pow((fParticle_Mass+water_mass),2));
359 
360  //SI: added protection particleEnergy0 >= depositEnergyCM
361  if (!statCode && (particleEnergy0 >= depositEnergyCM) )
362 
363  fParticleChangeForGamma->SetProposedKineticEnergy(particleEnergy0 - depositEnergyCM);
364 
365  else fParticleChangeForGamma->SetProposedKineticEnergy(particleEnergy0);
366 
367  fParticleChangeForGamma->ProposeLocalEnergyDeposit(depositEnergyCM);
368  }
369 }
G4double RandomizeThetaCM(G4double k, G4ParticleDefinition *aParticleDefinition)
G4double GetKineticEnergy() const
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
Hep3Vector cross(const Hep3Vector &) const
G4ParticleChangeForGamma * fParticleChangeForGamma
static const double pi
Definition: SystemOfUnits.h:53
Hep3Vector orthogonal() const
const G4ThreeVector & GetMomentumDirection() const
G4ParticleDefinition * GetDefinition() const
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
Here is the call graph for this function:

◆ SelectStationary()

void G4DNAIonElasticModel::SelectStationary ( G4bool  input)
inline

Definition at line 148 of file G4DNAIonElasticModel.hh.

149 {
150  statCode = input;
151 }
Here is the caller graph for this function:

◆ SetKillBelowThreshold()

void G4DNAIonElasticModel::SetKillBelowThreshold ( G4double  threshold)

Definition at line 507 of file G4DNAIonElasticModel.cc.

508 {
509  killBelowEnergy = threshold;
510 
511  if(killBelowEnergy < 100 * eV)
512  {
513  G4cout << "*** WARNING : the G4DNAIonElasticModel class is not "
514  "activated below 100 eV !"
515  << G4endl;
516  }
517 }
G4GLOB_DLL std::ostream G4cout
static const double eV
Definition: G4SIunits.hh:212
#define G4endl
Definition: G4ios.hh:61

◆ Theta()

G4double G4DNAIonElasticModel::Theta ( G4ParticleDefinition aParticleDefinition,
G4double  k,
G4double  integrDiff 
)
private

Definition at line 374 of file G4DNAIonElasticModel.cc.

376 {
377  G4double theta = 0.;
378  G4double valueT1 = 0;
379  G4double valueT2 = 0;
380  G4double valueE21 = 0;
381  G4double valueE22 = 0;
382  G4double valueE12 = 0;
383  G4double valueE11 = 0;
384  G4double xs11 = 0;
385  G4double xs12 = 0;
386  G4double xs21 = 0;
387  G4double xs22 = 0;
388 
389  std::vector<double>::iterator t2 = std::upper_bound(eTdummyVec.begin(),
390  eTdummyVec.end(), k);
391  std::vector<double>::iterator t1 = t2 - 1;
392 
393  std::vector<double>::iterator e12 = std::upper_bound(eVecm[(*t1)].begin(),
394  eVecm[(*t1)].end(),
395  integrDiff);
396  std::vector<double>::iterator e11 = e12 - 1;
397 
398  std::vector<double>::iterator e22 = std::upper_bound(eVecm[(*t2)].begin(),
399  eVecm[(*t2)].end(),
400  integrDiff);
401  std::vector<double>::iterator e21 = e22 - 1;
402 
403  valueT1 = *t1;
404  valueT2 = *t2;
405  valueE21 = *e21;
406  valueE22 = *e22;
407  valueE12 = *e12;
408  valueE11 = *e11;
409 
410  xs11 = fDiffCrossSectionData[valueT1][valueE11];
411  xs12 = fDiffCrossSectionData[valueT1][valueE12];
412  xs21 = fDiffCrossSectionData[valueT2][valueE21];
413  xs22 = fDiffCrossSectionData[valueT2][valueE22];
414 
415  if(xs11 == 0 && xs12 == 0 && xs21 == 0 && xs22 == 0) return (0.);
416 
417  theta = QuadInterpolator(valueE11, valueE12, valueE21, valueE22, xs11, xs12,
418  xs21, xs22, valueT1, valueT2, k, integrDiff);
419 
420  return theta;
421 }
TTree * t1
Definition: plottest35.C:26
TriDimensionMap fDiffCrossSectionData
std::vector< double > eTdummyVec
TTree * t2
Definition: plottest35.C:36
G4double QuadInterpolator(G4double e11, G4double e12, G4double e21, G4double e22, G4double x11, G4double x12, G4double x21, G4double x22, G4double t1, G4double t2, G4double t, G4double e)
double G4double
Definition: G4Types.hh:76
Here is the call graph for this function:
Here is the caller graph for this function:

Member Data Documentation

◆ eTdummyVec

std::vector<double> G4DNAIonElasticModel::eTdummyVec
private

Definition at line 130 of file G4DNAIonElasticModel.hh.

◆ eVecm

VecMap G4DNAIonElasticModel::eVecm
private

Definition at line 133 of file G4DNAIonElasticModel.hh.

◆ fDiffCrossSectionData

TriDimensionMap G4DNAIonElasticModel::fDiffCrossSectionData
private

Definition at line 128 of file G4DNAIonElasticModel.hh.

◆ fParticle_Mass

G4double G4DNAIonElasticModel::fParticle_Mass
private

Definition at line 99 of file G4DNAIonElasticModel.hh.

◆ fParticleChangeForGamma

G4ParticleChangeForGamma* G4DNAIonElasticModel::fParticleChangeForGamma
protected

Definition at line 84 of file G4DNAIonElasticModel.hh.

◆ fpMolWaterDensity

const std::vector<G4double>* G4DNAIonElasticModel::fpMolWaterDensity
private

Definition at line 91 of file G4DNAIonElasticModel.hh.

◆ fpTableData

G4DNACrossSectionDataSet* G4DNAIonElasticModel::fpTableData
private

Definition at line 102 of file G4DNAIonElasticModel.hh.

◆ highEnergyLimit

G4double G4DNAIonElasticModel::highEnergyLimit
private

Definition at line 95 of file G4DNAIonElasticModel.hh.

◆ isInitialised

G4bool G4DNAIonElasticModel::isInitialised
private

Definition at line 96 of file G4DNAIonElasticModel.hh.

◆ killBelowEnergy

G4double G4DNAIonElasticModel::killBelowEnergy
private

Definition at line 93 of file G4DNAIonElasticModel.hh.

◆ lowEnergyLimit

G4double G4DNAIonElasticModel::lowEnergyLimit
private

Definition at line 94 of file G4DNAIonElasticModel.hh.

◆ statCode

G4bool G4DNAIonElasticModel::statCode
private

Definition at line 88 of file G4DNAIonElasticModel.hh.

◆ verboseLevel

G4int G4DNAIonElasticModel::verboseLevel
private

Definition at line 97 of file G4DNAIonElasticModel.hh.


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