Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4PAIxSection Class Reference

#include <G4PAIxSection.hh>

Public Member Functions

 G4PAIxSection ()
 
 G4PAIxSection (G4MaterialCutsCouple *matCC)
 
 G4PAIxSection (G4int materialIndex, G4double maxEnergyTransfer)
 
 G4PAIxSection (G4int materialIndex, G4double maxEnergyTransfer, G4double betaGammaSq, G4double **photoAbsCof, G4int intNumber)
 
 G4PAIxSection (G4int materialIndex, G4double maxEnergyTransfer, G4double betaGammaSq)
 
 ~G4PAIxSection ()
 
void Initialize (const G4Material *material, G4double maxEnergyTransfer, G4double betaGammaSq, G4SandiaTable *)
 
void ComputeLowEnergyCof (const G4Material *material)
 
void ComputeLowEnergyCof ()
 
void InitPAI ()
 
void NormShift (G4double betaGammaSq)
 
void SplainPAI (G4double betaGammaSq)
 
G4double RutherfordIntegral (G4int intervalNumber, G4double limitLow, G4double limitHigh)
 
G4double ImPartDielectricConst (G4int intervalNumber, G4double energy)
 
G4double GetPhotonRange (G4double energy)
 
G4double GetElectronRange (G4double energy)
 
G4double RePartDielectricConst (G4double energy)
 
G4double DifPAIxSection (G4int intervalNumber, G4double betaGammaSq)
 
G4double PAIdNdxCerenkov (G4int intervalNumber, G4double betaGammaSq)
 
G4double PAIdNdxMM (G4int intervalNumber, G4double betaGammaSq)
 
G4double PAIdNdxPlasmon (G4int intervalNumber, G4double betaGammaSq)
 
G4double PAIdNdxResonance (G4int intervalNumber, G4double betaGammaSq)
 
void IntegralPAIxSection ()
 
void IntegralCerenkov ()
 
void IntegralMM ()
 
void IntegralPlasmon ()
 
void IntegralResonance ()
 
G4double SumOverInterval (G4int intervalNumber)
 
G4double SumOverIntervaldEdx (G4int intervalNumber)
 
G4double SumOverInterCerenkov (G4int intervalNumber)
 
G4double SumOverInterMM (G4int intervalNumber)
 
G4double SumOverInterPlasmon (G4int intervalNumber)
 
G4double SumOverInterResonance (G4int intervalNumber)
 
G4double SumOverBorder (G4int intervalNumber, G4double energy)
 
G4double SumOverBorderdEdx (G4int intervalNumber, G4double energy)
 
G4double SumOverBordCerenkov (G4int intervalNumber, G4double energy)
 
G4double SumOverBordMM (G4int intervalNumber, G4double energy)
 
G4double SumOverBordPlasmon (G4int intervalNumber, G4double energy)
 
G4double SumOverBordResonance (G4int intervalNumber, G4double energy)
 
G4double GetStepEnergyLoss (G4double step)
 
G4double GetStepCerenkovLoss (G4double step)
 
G4double GetStepMMLoss (G4double step)
 
G4double GetStepPlasmonLoss (G4double step)
 
G4double GetStepResonanceLoss (G4double step)
 
G4double GetEnergyTransfer ()
 
G4double GetCerenkovEnergyTransfer ()
 
G4double GetMMEnergyTransfer ()
 
G4double GetPlasmonEnergyTransfer ()
 
G4double GetResonanceEnergyTransfer ()
 
G4double GetRutherfordEnergyTransfer ()
 
G4int GetNumberOfGammas () const
 
G4int GetSplineSize () const
 
G4int GetIntervalNumber () const
 
G4double GetEnergyInterval (G4int i)
 
G4double GetDifPAIxSection (G4int i)
 
G4double GetPAIdNdxCerenkov (G4int i)
 
G4double GetPAIdNdxMM (G4int i)
 
G4double GetPAIdNdxPlasmon (G4int i)
 
G4double GetPAIdNdxResonance (G4int i)
 
G4double GetMeanEnergyLoss () const
 
G4double GetMeanCerenkovLoss () const
 
G4double GetMeanMMLoss () const
 
G4double GetMeanPlasmonLoss () const
 
G4double GetMeanResonanceLoss () const
 
G4double GetNormalizationCof () const
 
G4double GetLowEnergyCof () const
 
void SetVerbose (G4int v)
 
G4double GetPAItable (G4int i, G4int j) const
 
G4double GetLorentzFactor (G4int i) const
 
G4double GetSplineEnergy (G4int i) const
 
G4double GetIntegralPAIxSection (G4int i) const
 
G4double GetIntegralPAIdEdx (G4int i) const
 
G4double GetIntegralCerenkov (G4int i) const
 
G4double GetIntegralMM (G4int i) const
 
G4double GetIntegralPlasmon (G4int i) const
 
G4double GetIntegralResonance (G4int i) const
 

Detailed Description

Definition at line 68 of file G4PAIxSection.hh.

Constructor & Destructor Documentation

G4PAIxSection::G4PAIxSection ( )

Definition at line 92 of file G4PAIxSection.cc.

93 {
94  fSandia = 0;
95  fMatSandiaMatrix = 0;
96  fDensity = fElectronDensity = fNormalizationCof = fLowEnergyCof = 0.0;
97  fIntervalNumber = fSplineNumber = 0;
98  fVerbose = 0;
99 
100  fSplineEnergy = G4DataVector(fMaxSplineSize,0.0);
101  fRePartDielectricConst = G4DataVector(fMaxSplineSize,0.0);
102  fImPartDielectricConst = G4DataVector(fMaxSplineSize,0.0);
103  fIntegralTerm = G4DataVector(fMaxSplineSize,0.0);
104  fDifPAIxSection = G4DataVector(fMaxSplineSize,0.0);
105  fdNdxCerenkov = G4DataVector(fMaxSplineSize,0.0);
106  fdNdxPlasmon = G4DataVector(fMaxSplineSize,0.0);
107  fdNdxMM = G4DataVector(fMaxSplineSize,0.0);
108  fdNdxResonance = G4DataVector(fMaxSplineSize,0.0);
109  fIntegralPAIxSection = G4DataVector(fMaxSplineSize,0.0);
110  fIntegralPAIdEdx = G4DataVector(fMaxSplineSize,0.0);
111  fIntegralCerenkov = G4DataVector(fMaxSplineSize,0.0);
112  fIntegralPlasmon = G4DataVector(fMaxSplineSize,0.0);
113  fIntegralMM = G4DataVector(fMaxSplineSize,0.0);
114  fIntegralResonance = G4DataVector(fMaxSplineSize,0.0);
115 
116  fMaterialIndex = 0;
117 
118  for( G4int i = 0; i < 500; ++i )
119  {
120  for( G4int j = 0; j < 112; ++j ) fPAItable[i][j] = 0.0;
121  }
122 }
int G4int
Definition: G4Types.hh:78
G4PAIxSection::G4PAIxSection ( G4MaterialCutsCouple matCC)

Definition at line 129 of file G4PAIxSection.cc.

130 {
131  fDensity = matCC->GetMaterial()->GetDensity();
132  G4int matIndex = matCC->GetMaterial()->GetIndex();
133  fMaterialIndex = matIndex;
134 
135  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
136  fSandia = (*theMaterialTable)[matIndex]->GetSandiaTable();
137 
138  fVerbose = 0;
139 
140  G4int i, j;
141  fMatSandiaMatrix = new G4OrderedTable();
142 
143  for (i = 0; i < fSandia->GetMaxInterval()-1; i++)
144  {
145  fMatSandiaMatrix->push_back(new G4DataVector(5,0.));
146  }
147  for (i = 0; i < fSandia->GetMaxInterval()-1; i++)
148  {
149  (*(*fMatSandiaMatrix)[i])[0] = fSandia->GetSandiaMatTable(i,0);
150 
151  for(j = 1; j < 5; j++)
152  {
153  (*(*fMatSandiaMatrix)[i])[j] = fSandia->GetSandiaMatTable(i,j)*fDensity;
154  }
155  }
157  // fEnergyInterval = fA1 = fA2 = fA3 = fA4 = 0;
158 }
size_t GetIndex() const
Definition: G4Material.hh:262
static G4MaterialTable * GetMaterialTable()
Definition: G4Material.cc:587
std::vector< G4Material * > G4MaterialTable
G4double GetDensity() const
Definition: G4Material.hh:180
G4int GetMaxInterval() const
void ComputeLowEnergyCof()
int G4int
Definition: G4Types.hh:78
G4double GetSandiaMatTable(G4int, G4int) const
const G4Material * GetMaterial() const

Here is the call graph for this function:

G4PAIxSection::G4PAIxSection ( G4int  materialIndex,
G4double  maxEnergyTransfer 
)

Definition at line 162 of file G4PAIxSection.cc.

164 {
165  fSandia = 0;
166  fMatSandiaMatrix = 0;
167  fVerbose = 0;
168  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
169  G4int i, j;
170 
171  fMaterialIndex = materialIndex;
172  fDensity = (*theMaterialTable)[materialIndex]->GetDensity();
173  fElectronDensity = (*theMaterialTable)[materialIndex]->
174  GetElectronDensity();
175  fIntervalNumber = (*theMaterialTable)[materialIndex]->
176  GetSandiaTable()->GetMatNbOfIntervals();
177  fIntervalNumber--;
178  // G4cout<<fDensity<<"\t"<<fElectronDensity<<"\t"<<fIntervalNumber<<G4endl;
179 
180  fEnergyInterval = G4DataVector(fIntervalNumber+2,0.0);
181  fA1 = G4DataVector(fIntervalNumber+2,0.0);
182  fA2 = G4DataVector(fIntervalNumber+2,0.0);
183  fA3 = G4DataVector(fIntervalNumber+2,0.0);
184  fA4 = G4DataVector(fIntervalNumber+2,0.0);
185 
186  for(i = 1; i <= fIntervalNumber; i++ )
187  {
188  if(((*theMaterialTable)[materialIndex]->
189  GetSandiaTable()->GetSandiaCofForMaterial(i-1,0) >= maxEnergyTransfer) ||
190  i > fIntervalNumber )
191  {
192  fEnergyInterval[i] = maxEnergyTransfer;
193  fIntervalNumber = i;
194  break;
195  }
196  fEnergyInterval[i] = (*theMaterialTable)[materialIndex]->
197  GetSandiaTable()->GetSandiaCofForMaterial(i-1,0);
198  fA1[i] = (*theMaterialTable)[materialIndex]->
199  GetSandiaTable()->GetSandiaCofForMaterial(i-1,1);
200  fA2[i] = (*theMaterialTable)[materialIndex]->
201  GetSandiaTable()->GetSandiaCofForMaterial(i-1,2);
202  fA3[i] = (*theMaterialTable)[materialIndex]->
203  GetSandiaTable()->GetSandiaCofForMaterial(i-1,3);
204  fA4[i] = (*theMaterialTable)[materialIndex]->
205  GetSandiaTable()->GetSandiaCofForMaterial(i-1,4);
206  // G4cout<<i<<"\t"<<fEnergyInterval[i]<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
207  // <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl;
208  }
209  if(fEnergyInterval[fIntervalNumber] != maxEnergyTransfer)
210  {
211  fIntervalNumber++;
212  fEnergyInterval[fIntervalNumber] = maxEnergyTransfer;
213  }
214 
215  // Now checking, if two borders are too close together
216 
217  for(i=1;i<fIntervalNumber;i++)
218  {
219  if(fEnergyInterval[i+1]-fEnergyInterval[i] >
220  1.5*fDelta*(fEnergyInterval[i+1]+fEnergyInterval[i]))
221  {
222  continue;
223  }
224  else
225  {
226  for(j=i;j<fIntervalNumber;j++)
227  {
228  fEnergyInterval[j] = fEnergyInterval[j+1];
229  fA1[j] = fA1[j+1];
230  fA2[j] = fA2[j+1];
231  fA3[j] = fA3[j+1];
232  fA4[j] = fA4[j+1];
233  }
234  fIntervalNumber--;
235  i--;
236  }
237  }
238 
239 
240  /* *********************************
241 
242  fSplineEnergy = new G4double[fMaxSplineSize];
243  fRePartDielectricConst = new G4double[fMaxSplineSize];
244  fImPartDielectricConst = new G4double[fMaxSplineSize];
245  fIntegralTerm = new G4double[fMaxSplineSize];
246  fDifPAIxSection = new G4double[fMaxSplineSize];
247  fIntegralPAIxSection = new G4double[fMaxSplineSize];
248 
249  for(i=0;i<fMaxSplineSize;i++)
250  {
251  fSplineEnergy[i] = 0.0;
252  fRePartDielectricConst[i] = 0.0;
253  fImPartDielectricConst[i] = 0.0;
254  fIntegralTerm[i] = 0.0;
255  fDifPAIxSection[i] = 0.0;
256  fIntegralPAIxSection[i] = 0.0;
257  }
258  ************************************************** */
260  InitPAI(); // create arrays allocated above
261  /*
262  delete[] fEnergyInterval;
263  delete[] fA1;
264  delete[] fA2;
265  delete[] fA3;
266  delete[] fA4;
267  */
268 }
static G4MaterialTable * GetMaterialTable()
Definition: G4Material.cc:587
std::vector< G4Material * > G4MaterialTable
void ComputeLowEnergyCof()
int G4int
Definition: G4Types.hh:78

Here is the call graph for this function:

G4PAIxSection::G4PAIxSection ( G4int  materialIndex,
G4double  maxEnergyTransfer,
G4double  betaGammaSq,
G4double **  photoAbsCof,
G4int  intNumber 
)

Definition at line 274 of file G4PAIxSection.cc.

279 {
280  fSandia = 0;
281  fDensity = fElectronDensity = fNormalizationCof = fLowEnergyCof = 0.0;
282  fIntervalNumber = fSplineNumber = 0;
283  fVerbose = 0;
284 
285  fSplineEnergy = G4DataVector(500,0.0);
286  fRePartDielectricConst = G4DataVector(500,0.0);
287  fImPartDielectricConst = G4DataVector(500,0.0);
288  fIntegralTerm = G4DataVector(500,0.0);
289  fDifPAIxSection = G4DataVector(500,0.0);
290  fdNdxCerenkov = G4DataVector(500,0.0);
291  fdNdxPlasmon = G4DataVector(500,0.0);
292  fdNdxMM = G4DataVector(500,0.0);
293  fdNdxResonance = G4DataVector(500,0.0);
294  fIntegralPAIxSection = G4DataVector(500,0.0);
295  fIntegralPAIdEdx = G4DataVector(500,0.0);
296  fIntegralCerenkov = G4DataVector(500,0.0);
297  fIntegralPlasmon = G4DataVector(500,0.0);
298  fIntegralMM = G4DataVector(500,0.0);
299  fIntegralResonance = G4DataVector(500,0.0);
300 
301  for( G4int i = 0; i < 500; ++i )
302  {
303  for( G4int j = 0; j < 112; ++j ) fPAItable[i][j] = 0.0;
304  }
305 
306  fSandia = 0;
307  fMatSandiaMatrix = 0;
308  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
309  G4int i, j;
310 
311  fMaterialIndex = materialIndex;
312  fDensity = (*theMaterialTable)[materialIndex]->GetDensity();
313  fElectronDensity = (*theMaterialTable)[materialIndex]->GetElectronDensity();
314 
315  fIntervalNumber = intNumber;
316  fIntervalNumber--;
317  // G4cout<<fDensity<<"\t"<<fElectronDensity<<"\t"<<fIntervalNumber<<G4endl;
318 
319  fEnergyInterval = G4DataVector(fIntervalNumber+2,0.0);
320  fA1 = G4DataVector(fIntervalNumber+2,0.0);
321  fA2 = G4DataVector(fIntervalNumber+2,0.0);
322  fA3 = G4DataVector(fIntervalNumber+2,0.0);
323  fA4 = G4DataVector(fIntervalNumber+2,0.0);
324 
325 
326  /*
327  fEnergyInterval = new G4double[fIntervalNumber+2];
328  fA1 = new G4double[fIntervalNumber+2];
329  fA2 = new G4double[fIntervalNumber+2];
330  fA3 = new G4double[fIntervalNumber+2];
331  fA4 = new G4double[fIntervalNumber+2];
332  */
333  for( i = 1; i <= fIntervalNumber; i++ )
334  {
335  if( ( photoAbsCof[i-1][0] >= maxEnergyTransfer ) ||
336  i > fIntervalNumber )
337  {
338  fEnergyInterval[i] = maxEnergyTransfer;
339  fIntervalNumber = i;
340  break;
341  }
342  fEnergyInterval[i] = photoAbsCof[i-1][0];
343  fA1[i] = photoAbsCof[i-1][1];
344  fA2[i] = photoAbsCof[i-1][2];
345  fA3[i] = photoAbsCof[i-1][3];
346  fA4[i] = photoAbsCof[i-1][4];
347  // G4cout<<i<<"\t"<<fEnergyInterval[i]/keV<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
348  // <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl;
349  }
350  // G4cout<<"i last = "<<i<<"; "<<"fIntervalNumber = "<<fIntervalNumber<<G4endl;
351 
352  if(fEnergyInterval[fIntervalNumber] != maxEnergyTransfer)
353  {
354  fIntervalNumber++;
355  fEnergyInterval[fIntervalNumber] = maxEnergyTransfer;
356  }
357  // G4cout<<"after check of max energy transfer"<<G4endl;
358 
359  for( i = 1; i <= fIntervalNumber; i++ )
360  {
361  // G4cout<<i<<"\t"<<fEnergyInterval[i]/keV<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
362  // <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl;
363  }
364  // Now checking, if two borders are too close together
365 
366  for( i = 1; i < fIntervalNumber; i++ )
367  {
368  if(fEnergyInterval[i+1]-fEnergyInterval[i] >
369  1.5*fDelta*(fEnergyInterval[i+1]+fEnergyInterval[i]))
370  {
371  continue;
372  }
373  else
374  {
375  for(j=i;j<fIntervalNumber;j++)
376  {
377  fEnergyInterval[j] = fEnergyInterval[j+1];
378  fA1[j] = fA1[j+1];
379  fA2[j] = fA2[j+1];
380  fA3[j] = fA3[j+1];
381  fA4[j] = fA4[j+1];
382  }
383  fIntervalNumber--;
384  i--;
385  }
386  }
387  // G4cout<<"after check of close borders"<<G4endl;
388 
389  for( i = 1; i <= fIntervalNumber; i++ )
390  {
391  // G4cout<<i<<"\t"<<fEnergyInterval[i]/keV<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
392  // <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl;
393  }
394 
395  // Preparation of fSplineEnergy array corresponding to min ionisation, G~4
396 
398  G4double betaGammaSqRef =
399  fLorentzFactor[fRefGammaNumber]*fLorentzFactor[fRefGammaNumber] - 1;
400 
401  NormShift(betaGammaSqRef);
402  SplainPAI(betaGammaSqRef);
403 
404  // Preparation of integral PAI cross section for input betaGammaSq
405 
406  for(i = 1; i <= fSplineNumber; i++)
407  {
408  fdNdxCerenkov[i] = PAIdNdxCerenkov(i,betaGammaSq);
409  fdNdxMM[i] = PAIdNdxMM(i,betaGammaSq);
410  fdNdxPlasmon[i] = PAIdNdxPlasmon(i,betaGammaSq);
411  fdNdxResonance[i] = PAIdNdxResonance(i,betaGammaSq);
412  fDifPAIxSection[i] = DifPAIxSection(i,betaGammaSq);
413 
414  // G4cout<<i<<"; dNdxC = "<<fdNdxCerenkov[i]<<"; dNdxP = "<<fdNdxPlasmon[i]
415  // <<"; dNdxPAI = "<<fDifPAIxSection[i]<<G4endl;
416  }
418  IntegralMM();
419  IntegralPlasmon();
422  /*
423  delete[] fEnergyInterval;
424  delete[] fA1;
425  delete[] fA2;
426  delete[] fA3;
427  delete[] fA4;
428  */
429 }
void IntegralPAIxSection()
void SplainPAI(G4double betaGammaSq)
G4double PAIdNdxResonance(G4int intervalNumber, G4double betaGammaSq)
void IntegralCerenkov()
void IntegralPlasmon()
void IntegralResonance()
G4double PAIdNdxCerenkov(G4int intervalNumber, G4double betaGammaSq)
static G4MaterialTable * GetMaterialTable()
Definition: G4Material.cc:587
std::vector< G4Material * > G4MaterialTable
void NormShift(G4double betaGammaSq)
G4double PAIdNdxPlasmon(G4int intervalNumber, G4double betaGammaSq)
void ComputeLowEnergyCof()
int G4int
Definition: G4Types.hh:78
G4double PAIdNdxMM(G4int intervalNumber, G4double betaGammaSq)
double G4double
Definition: G4Types.hh:76
G4double DifPAIxSection(G4int intervalNumber, G4double betaGammaSq)

Here is the call graph for this function:

G4PAIxSection::G4PAIxSection ( G4int  materialIndex,
G4double  maxEnergyTransfer,
G4double  betaGammaSq 
)

Definition at line 435 of file G4PAIxSection.cc.

438 {
439  fSandia = 0;
440  fMatSandiaMatrix = 0;
441  fVerbose = 0;
442  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
443 
444  G4int i, j, numberOfElements;
445 
446  fMaterialIndex = materialIndex;
447  fDensity = (*theMaterialTable)[materialIndex]->GetDensity();
448  fElectronDensity = (*theMaterialTable)[materialIndex]->GetElectronDensity();
449  numberOfElements = (*theMaterialTable)[materialIndex]->GetNumberOfElements();
450 
451  G4int* thisMaterialZ = new G4int[numberOfElements];
452 
453  for( i = 0; i < numberOfElements; i++ )
454  {
455  thisMaterialZ[i] = (G4int)(*theMaterialTable)[materialIndex]->
456  GetElement(i)->GetZ();
457  }
458  // fSandia = new G4SandiaTable(materialIndex);
459  fSandia = (*theMaterialTable)[materialIndex]->GetSandiaTable();
460  G4SandiaTable thisMaterialSandiaTable(materialIndex);
461  fIntervalNumber = thisMaterialSandiaTable.SandiaIntervals(thisMaterialZ,
462  numberOfElements);
463  fIntervalNumber = thisMaterialSandiaTable.SandiaMixing
464  ( thisMaterialZ ,
465  (*theMaterialTable)[materialIndex]->GetFractionVector() ,
466  numberOfElements,fIntervalNumber);
467 
468  fIntervalNumber--;
469 
470  fEnergyInterval = G4DataVector(fIntervalNumber+2,0.0);
471  fA1 = G4DataVector(fIntervalNumber+2,0.0);
472  fA2 = G4DataVector(fIntervalNumber+2,0.0);
473  fA3 = G4DataVector(fIntervalNumber+2,0.0);
474  fA4 = G4DataVector(fIntervalNumber+2,0.0);
475 
476  /*
477  fEnergyInterval = new G4double[fIntervalNumber+2];
478  fA1 = new G4double[fIntervalNumber+2];
479  fA2 = new G4double[fIntervalNumber+2];
480  fA3 = new G4double[fIntervalNumber+2];
481  fA4 = new G4double[fIntervalNumber+2];
482  */
483  for( i = 1; i <= fIntervalNumber; i++ )
484  {
485  if((thisMaterialSandiaTable.GetPhotoAbsorpCof(i,0) >= maxEnergyTransfer) ||
486  i > fIntervalNumber)
487  {
488  fEnergyInterval[i] = maxEnergyTransfer;
489  fIntervalNumber = i;
490  break;
491  }
492  fEnergyInterval[i] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i,0);
493  fA1[i] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i,1)*fDensity;
494  fA2[i] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i,2)*fDensity;
495  fA3[i] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i,3)*fDensity;
496  fA4[i] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i,4)*fDensity;
497 
498  }
499  if(fEnergyInterval[fIntervalNumber] != maxEnergyTransfer)
500  {
501  fIntervalNumber++;
502  fEnergyInterval[fIntervalNumber] = maxEnergyTransfer;
503  fA1[fIntervalNumber] = fA1[fIntervalNumber-1];
504  fA2[fIntervalNumber] = fA2[fIntervalNumber-1];
505  fA3[fIntervalNumber] = fA3[fIntervalNumber-1];
506  fA4[fIntervalNumber] = fA4[fIntervalNumber-1];
507  }
508  for(i=1;i<=fIntervalNumber;i++)
509  {
510  // G4cout<<fEnergyInterval[i]<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
511  // <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl;
512  }
513  // Now checking, if two borders are too close together
514 
515  for( i = 1; i < fIntervalNumber; i++ )
516  {
517  if(fEnergyInterval[i+1]-fEnergyInterval[i] >
518  1.5*fDelta*(fEnergyInterval[i+1]+fEnergyInterval[i]))
519  {
520  continue;
521  }
522  else
523  {
524  for( j = i; j < fIntervalNumber; j++ )
525  {
526  fEnergyInterval[j] = fEnergyInterval[j+1];
527  fA1[j] = fA1[j+1];
528  fA2[j] = fA2[j+1];
529  fA3[j] = fA3[j+1];
530  fA4[j] = fA4[j+1];
531  }
532  fIntervalNumber--;
533  i--;
534  }
535  }
536 
537  /* *********************************
538  fSplineEnergy = new G4double[fMaxSplineSize];
539  fRePartDielectricConst = new G4double[fMaxSplineSize];
540  fImPartDielectricConst = new G4double[fMaxSplineSize];
541  fIntegralTerm = new G4double[fMaxSplineSize];
542  fDifPAIxSection = new G4double[fMaxSplineSize];
543  fIntegralPAIxSection = new G4double[fMaxSplineSize];
544 
545  for(i=0;i<fMaxSplineSize;i++)
546  {
547  fSplineEnergy[i] = 0.0;
548  fRePartDielectricConst[i] = 0.0;
549  fImPartDielectricConst[i] = 0.0;
550  fIntegralTerm[i] = 0.0;
551  fDifPAIxSection[i] = 0.0;
552  fIntegralPAIxSection[i] = 0.0;
553  }
554  */
555 
556  // Preparation of fSplineEnergy array corresponding to min ionisation, G~4
557 
559  G4double betaGammaSqRef =
560  fLorentzFactor[fRefGammaNumber]*fLorentzFactor[fRefGammaNumber] - 1;
561 
562  NormShift(betaGammaSqRef);
563  SplainPAI(betaGammaSqRef);
564 
565  // Preparation of integral PAI cross section for input betaGammaSq
566 
567  for(i = 1; i <= fSplineNumber; i++)
568  {
569  fDifPAIxSection[i] = DifPAIxSection(i,betaGammaSq);
570  fdNdxCerenkov[i] = PAIdNdxCerenkov(i,betaGammaSq);
571  fdNdxMM[i] = PAIdNdxMM(i,betaGammaSq);
572  fdNdxPlasmon[i] = PAIdNdxPlasmon(i,betaGammaSq);
573  fdNdxResonance[i] = PAIdNdxResonance(i,betaGammaSq);
574  }
577  IntegralMM();
578  IntegralPlasmon();
580 }
void IntegralPAIxSection()
void SplainPAI(G4double betaGammaSq)
G4double PAIdNdxResonance(G4int intervalNumber, G4double betaGammaSq)
void IntegralCerenkov()
void IntegralPlasmon()
void IntegralResonance()
G4double PAIdNdxCerenkov(G4int intervalNumber, G4double betaGammaSq)
static G4MaterialTable * GetMaterialTable()
Definition: G4Material.cc:587
std::vector< G4Material * > G4MaterialTable
void NormShift(G4double betaGammaSq)
G4double PAIdNdxPlasmon(G4int intervalNumber, G4double betaGammaSq)
void ComputeLowEnergyCof()
int G4int
Definition: G4Types.hh:78
G4double PAIdNdxMM(G4int intervalNumber, G4double betaGammaSq)
double G4double
Definition: G4Types.hh:76
G4double DifPAIxSection(G4int intervalNumber, G4double betaGammaSq)

Here is the call graph for this function:

G4PAIxSection::~G4PAIxSection ( )

Definition at line 586 of file G4PAIxSection.cc.

587 {
588  /* ************************
589  delete[] fSplineEnergy ;
590  delete[] fRePartDielectricConst;
591  delete[] fImPartDielectricConst;
592  delete[] fIntegralTerm ;
593  delete[] fDifPAIxSection ;
594  delete[] fIntegralPAIxSection ;
595  */
596  delete fMatSandiaMatrix;
597 }

Member Function Documentation

void G4PAIxSection::ComputeLowEnergyCof ( const G4Material material)

Definition at line 746 of file G4PAIxSection.cc.

747 {
748  G4int i, numberOfElements = material->GetNumberOfElements();
749  G4double sumZ = 0., sumCof = 0.;
750 
751  static const G4double p0 = 1.20923e+00;
752  static const G4double p1 = 3.53256e-01;
753  static const G4double p2 = -1.45052e-03;
754 
755  G4double* thisMaterialZ = new G4double[numberOfElements];
756  G4double* thisMaterialCof = new G4double[numberOfElements];
757 
758  for( i = 0; i < numberOfElements; i++ )
759  {
760  thisMaterialZ[i] = material->GetElement(i)->GetZ();
761  sumZ += thisMaterialZ[i];
762  thisMaterialCof[i] = p0+p1*thisMaterialZ[i]+p2*thisMaterialZ[i]*thisMaterialZ[i];
763  }
764  for( i = 0; i < numberOfElements; i++ )
765  {
766  sumCof += thisMaterialCof[i]*thisMaterialZ[i]/sumZ;
767  }
768  fLowEnergyCof = sumCof;
769  delete [] thisMaterialZ;
770  delete [] thisMaterialCof;
771  // G4cout<<"fLowEnergyCof = "<<fLowEnergyCof<<G4endl;
772 }
G4double GetZ() const
Definition: G4Element.hh:131
const G4Element * GetElement(G4int iel) const
Definition: G4Material.hh:202
int G4int
Definition: G4Types.hh:78
size_t GetNumberOfElements() const
Definition: G4Material.hh:186
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

void G4PAIxSection::ComputeLowEnergyCof ( )

Definition at line 779 of file G4PAIxSection.cc.

780 {
781  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
782  G4int i, numberOfElements = (*theMaterialTable)[fMaterialIndex]->GetNumberOfElements();
783  G4double sumZ = 0., sumCof = 0.;
784 
785  const G4double p0 = 1.20923e+00;
786  const G4double p1 = 3.53256e-01;
787  const G4double p2 = -1.45052e-03;
788 
789  G4double* thisMaterialZ = new G4double[numberOfElements];
790  G4double* thisMaterialCof = new G4double[numberOfElements];
791 
792  for( i = 0; i < numberOfElements; i++ )
793  {
794  thisMaterialZ[i] = (*theMaterialTable)[fMaterialIndex]->GetElement(i)->GetZ();
795  sumZ += thisMaterialZ[i];
796  thisMaterialCof[i] = p0+p1*thisMaterialZ[i]+p2*thisMaterialZ[i]*thisMaterialZ[i];
797  }
798  for( i = 0; i < numberOfElements; i++ )
799  {
800  sumCof += thisMaterialCof[i]*thisMaterialZ[i]/sumZ;
801  }
802  fLowEnergyCof = sumCof;
803  // G4cout<<"fLowEnergyCof = "<<fLowEnergyCof<<G4endl;
804  delete [] thisMaterialZ;
805  delete [] thisMaterialCof;
806 }
static G4MaterialTable * GetMaterialTable()
Definition: G4Material.cc:587
std::vector< G4Material * > G4MaterialTable
int G4int
Definition: G4Types.hh:78
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

G4double G4PAIxSection::DifPAIxSection ( G4int  intervalNumber,
G4double  betaGammaSq 
)

Definition at line 1225 of file G4PAIxSection.cc.

1227 {
1228  G4double cof,x1,x2,x3,x4,x5,x6,x7,x8,result;
1229 
1230  G4double betaBohr = fine_structure_const;
1231  // G4double betaBohr2 = fine_structure_const*fine_structure_const;
1232  // G4double betaBohr3 = betaBohr*betaBohr2; // *4.0;
1233 
1234  G4double be2 = betaGammaSq/(1 + betaGammaSq);
1235  G4double beta = sqrt(be2);
1236  // G4double be3 = beta*be2;
1237 
1238  cof = 1.;
1239  x1 = log(2*electron_mass_c2/fSplineEnergy[i]);
1240 
1241  if( betaGammaSq < 0.01 ) x2 = log(be2);
1242  else
1243  {
1244  x2 = -log( (1/betaGammaSq - fRePartDielectricConst[i])*
1245  (1/betaGammaSq - fRePartDielectricConst[i]) +
1246  fImPartDielectricConst[i]*fImPartDielectricConst[i] )/2;
1247  }
1248  if( fImPartDielectricConst[i] == 0.0 ||betaGammaSq < 0.01 )
1249  {
1250  x6 = 0.;
1251  }
1252  else
1253  {
1254  x3 = -fRePartDielectricConst[i] + 1/betaGammaSq;
1255  x5 = -1 - fRePartDielectricConst[i] +
1256  be2*((1 +fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) +
1257  fImPartDielectricConst[i]*fImPartDielectricConst[i]);
1258 
1259  x7 = atan2(fImPartDielectricConst[i],x3);
1260  x6 = x5 * x7;
1261  }
1262  // if(fImPartDielectricConst[i] == 0) x6 = 0.;
1263 
1264  x4 = ((x1 + x2)*fImPartDielectricConst[i] + x6)/hbarc;
1265 
1266  // if( x4 < 0.0 ) x4 = 0.0;
1267 
1268  x8 = (1 + fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) +
1269  fImPartDielectricConst[i]*fImPartDielectricConst[i];
1270 
1271  result = (x4 + cof*fIntegralTerm[i]/fSplineEnergy[i]/fSplineEnergy[i]);
1272 
1273  if( result < 1.0e-8 ) result = 1.0e-8;
1274 
1275  result *= fine_structure_const/be2/pi;
1276 
1277  // low energy correction
1278 
1279  G4double lowCof = fLowEnergyCof; // 6.0 ; // Ar ~ 4.; -> fLowCof as f(Z1,Z2)?
1280 
1281  result *= (1 - exp(-beta/betaBohr/lowCof));
1282 
1283 
1284  // result *= (1 - exp(-be2/betaBohr2/lowCof));
1285 
1286  // result *= (1 - exp(-be3/betaBohr3/lowCof)); // ~ be for be<<betaBohr
1287 
1288  // result *= (1 - exp(-be4/betaBohr4/lowCof));
1289 
1290  if(fDensity >= 0.1)
1291  {
1292  result /= x8;
1293  }
1294  return result;
1295 
1296 } // end of DifPAIxSection
G4double G4ParticleHPJENDLHEData::G4double result
float electron_mass_c2
Definition: hepunit.py:274
static constexpr double pi
Definition: G4SIunits.hh:75
double G4double
Definition: G4Types.hh:76
G4double G4PAIxSection::GetCerenkovEnergyTransfer ( )

Definition at line 2335 of file G4PAIxSection.cc.

2336 {
2337  G4int iTransfer ;
2338 
2339  G4double energyTransfer, position;
2340 
2341  position = fIntegralCerenkov[1]*G4UniformRand();
2342 
2343  for( iTransfer = 1; iTransfer <= fSplineNumber; iTransfer++ )
2344  {
2345  if( position >= fIntegralCerenkov[iTransfer] ) break;
2346  }
2347  if(iTransfer > fSplineNumber) iTransfer--;
2348 
2349  energyTransfer = fSplineEnergy[iTransfer];
2350 
2351  if(iTransfer > 1)
2352  {
2353  energyTransfer -= (fSplineEnergy[iTransfer]-fSplineEnergy[iTransfer-1])*G4UniformRand();
2354  }
2355  return energyTransfer;
2356 }
int G4int
Definition: G4Types.hh:78
#define position
Definition: xmlparse.cc:622
#define G4UniformRand()
Definition: Randomize.hh:97
double G4double
Definition: G4Types.hh:76
G4double G4PAIxSection::GetDifPAIxSection ( G4int  i)
inline

Definition at line 181 of file G4PAIxSection.hh.

181 { return fDifPAIxSection[i]; }
G4double G4PAIxSection::GetElectronRange ( G4double  energy)

Definition at line 1129 of file G4PAIxSection.cc.

1130 {
1131  G4double range;
1132  /*
1133  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
1134 
1135  G4double Z = (*theMaterialTable)[fMaterialIndex]->GetIonisation()->GetZeffective();
1136  G4double A = (*theMaterialTable)[fMaterialIndex]->GetA();
1137 
1138  energy /= keV; // energy in keV in parametrised formula
1139 
1140  if (energy < 10.)
1141  {
1142  range = 3.872e-3*A/Z;
1143  range *= pow(energy,1.492);
1144  }
1145  else
1146  {
1147  range = 6.97e-3*pow(energy,1.6);
1148  }
1149  */
1150  // Blum&Rolandi Particle Detection with Drift Chambers, p. 7
1151 
1152  G4double cofA = 5.37e-4*g/cm2/keV;
1153  G4double cofB = 0.9815;
1154  G4double cofC = 3.123e-3/keV;
1155  // energy /= keV;
1156 
1157  range = cofA*energy*( 1 - cofB/(1 + cofC*energy) );
1158 
1159  // range *= g/cm2;
1160  range /= fDensity;
1161 
1162  return range;
1163 }
static constexpr double cm2
Definition: G4SIunits.hh:120
function g(Y1, Y2, PT2)
Definition: hijing1.383.f:5205
const G4ParticleDefinition const G4Material *G4double range
G4double energy(const ThreeVector &p, const G4double m)
double G4double
Definition: G4Types.hh:76
static constexpr double keV
Definition: G4SIunits.hh:216

Here is the call graph for this function:

G4double G4PAIxSection::GetEnergyInterval ( G4int  i)
inline

Definition at line 179 of file G4PAIxSection.hh.

179 { return fEnergyInterval[i]; }
G4double G4PAIxSection::GetEnergyTransfer ( )

Definition at line 2254 of file G4PAIxSection.cc.

2255 {
2256  G4int iTransfer ;
2257 
2258  G4double energyTransfer, position;
2259 
2260  position = fIntegralPAIxSection[1]*G4UniformRand();
2261 
2262  for( iTransfer = 1; iTransfer <= fSplineNumber; iTransfer++ )
2263  {
2264  if( position >= fIntegralPAIxSection[iTransfer] ) break;
2265  }
2266  if(iTransfer > fSplineNumber) iTransfer--;
2267 
2268  energyTransfer = fSplineEnergy[iTransfer];
2269 
2270  if(iTransfer > 1)
2271  {
2272  energyTransfer -= (fSplineEnergy[iTransfer]-fSplineEnergy[iTransfer-1])*G4UniformRand();
2273  }
2274  return energyTransfer;
2275 }
int G4int
Definition: G4Types.hh:78
#define position
Definition: xmlparse.cc:622
#define G4UniformRand()
Definition: Randomize.hh:97
double G4double
Definition: G4Types.hh:76
G4double G4PAIxSection::GetIntegralCerenkov ( G4int  i) const
inline

Definition at line 309 of file G4PAIxSection.hh.

310 {
311  if(i < 1 || i > fSplineNumber) { CallError(i, "GetIntegralCerenkov"); }
312  return fIntegralCerenkov[i];
313 }
G4double G4PAIxSection::GetIntegralMM ( G4int  i) const
inline

Definition at line 315 of file G4PAIxSection.hh.

316 {
317  if(i < 1 || i > fSplineNumber) { CallError(i, "GetIntegralMM"); }
318  return fIntegralMM[i];
319 }
G4double G4PAIxSection::GetIntegralPAIdEdx ( G4int  i) const
inline

Definition at line 303 of file G4PAIxSection.hh.

304 {
305  if(i < 1 || i > fSplineNumber) { CallError(i, "GetIntegralPAIdEdx"); }
306  return fIntegralPAIdEdx[i];
307 }
G4double G4PAIxSection::GetIntegralPAIxSection ( G4int  i) const
inline

Definition at line 297 of file G4PAIxSection.hh.

298 {
299  if(i < 1 || i > fSplineNumber) { CallError(i, "GetIntegralPAIxSection"); }
300  return fIntegralPAIxSection[i];
301 }
G4double G4PAIxSection::GetIntegralPlasmon ( G4int  i) const
inline

Definition at line 321 of file G4PAIxSection.hh.

322 {
323  if(i < 1 || i > fSplineNumber) { CallError(i, "GetIntegralPlasmon"); }
324  return fIntegralPlasmon[i];
325 }
G4double G4PAIxSection::GetIntegralResonance ( G4int  i) const
inline

Definition at line 327 of file G4PAIxSection.hh.

328 {
329  if(i < 1 || i > fSplineNumber) { CallError(i, "GetIntegralResonance"); }
330  return fIntegralResonance[i];
331 }
G4int G4PAIxSection::GetIntervalNumber ( ) const
inline

Definition at line 177 of file G4PAIxSection.hh.

177 { return fIntervalNumber; }
G4double G4PAIxSection::GetLorentzFactor ( G4int  i) const
inline

Definition at line 286 of file G4PAIxSection.hh.

287 {
288  return fLorentzFactor[j];
289 }
G4double G4PAIxSection::GetLowEnergyCof ( ) const
inline

Definition at line 195 of file G4PAIxSection.hh.

195 { return fLowEnergyCof; }
G4double G4PAIxSection::GetMeanCerenkovLoss ( ) const
inline

Definition at line 188 of file G4PAIxSection.hh.

188 {return fIntegralCerenkov[0]; }
G4double G4PAIxSection::GetMeanEnergyLoss ( ) const
inline

Definition at line 187 of file G4PAIxSection.hh.

187 {return fIntegralPAIxSection[0]; }
G4double G4PAIxSection::GetMeanMMLoss ( ) const
inline

Definition at line 189 of file G4PAIxSection.hh.

189 {return fIntegralMM[0]; }
G4double G4PAIxSection::GetMeanPlasmonLoss ( ) const
inline

Definition at line 190 of file G4PAIxSection.hh.

190 {return fIntegralPlasmon[0]; }
G4double G4PAIxSection::GetMeanResonanceLoss ( ) const
inline

Definition at line 191 of file G4PAIxSection.hh.

191 {return fIntegralResonance[0]; }
G4double G4PAIxSection::GetMMEnergyTransfer ( )

Definition at line 2362 of file G4PAIxSection.cc.

2363 {
2364  G4int iTransfer ;
2365 
2366  G4double energyTransfer, position;
2367 
2368  position = fIntegralMM[1]*G4UniformRand();
2369 
2370  for( iTransfer = 1; iTransfer <= fSplineNumber; iTransfer++ )
2371  {
2372  if( position >= fIntegralMM[iTransfer] ) break;
2373  }
2374  if(iTransfer > fSplineNumber) iTransfer--;
2375 
2376  energyTransfer = fSplineEnergy[iTransfer];
2377 
2378  if(iTransfer > 1)
2379  {
2380  energyTransfer -= (fSplineEnergy[iTransfer]-fSplineEnergy[iTransfer-1])*G4UniformRand();
2381  }
2382  return energyTransfer;
2383 }
int G4int
Definition: G4Types.hh:78
#define position
Definition: xmlparse.cc:622
#define G4UniformRand()
Definition: Randomize.hh:97
double G4double
Definition: G4Types.hh:76
G4double G4PAIxSection::GetNormalizationCof ( ) const
inline

Definition at line 193 of file G4PAIxSection.hh.

193 { return fNormalizationCof; }
G4int G4PAIxSection::GetNumberOfGammas ( ) const
inline

Definition at line 173 of file G4PAIxSection.hh.

173 { return fNumberOfGammas; }
G4double G4PAIxSection::GetPAIdNdxCerenkov ( G4int  i)
inline

Definition at line 182 of file G4PAIxSection.hh.

182 { return fdNdxCerenkov[i]; }
G4double G4PAIxSection::GetPAIdNdxMM ( G4int  i)
inline

Definition at line 183 of file G4PAIxSection.hh.

183 { return fdNdxMM[i]; }
G4double G4PAIxSection::GetPAIdNdxPlasmon ( G4int  i)
inline

Definition at line 184 of file G4PAIxSection.hh.

184 { return fdNdxPlasmon[i]; }
G4double G4PAIxSection::GetPAIdNdxResonance ( G4int  i)
inline

Definition at line 185 of file G4PAIxSection.hh.

185 { return fdNdxResonance[i]; }
G4double G4PAIxSection::GetPAItable ( G4int  i,
G4int  j 
) const
inline

Definition at line 281 of file G4PAIxSection.hh.

282 {
283  return fPAItable[i][j];
284 }
G4double G4PAIxSection::GetPhotonRange ( G4double  energy)

Definition at line 1096 of file G4PAIxSection.cc.

1097 {
1098  G4int i;
1099  G4double energy2, energy3, energy4, result, lambda;
1100 
1101  energy2 = energy1*energy1;
1102  energy3 = energy2*energy1;
1103  energy4 = energy3*energy1;
1104 
1105  // G4double* SandiaCof = fSandia->GetSandiaCofForMaterialPAI(energy1);
1106  // result = SandiaCof[0]/energy1+SandiaCof[1]/energy2+SandiaCof[2]/energy3+SandiaCof[3]/energy4;
1107  // result *= fDensity;
1108 
1109  for( i = 1; i <= fIntervalNumber; i++ )
1110  {
1111  if( energy1 < fEnergyInterval[i]) break;
1112  }
1113  i--;
1114  if(i == 0) i = 1;
1115 
1116  result = fA1[i]/energy1+fA2[i]/energy2+fA3[i]/energy3+fA4[i]/energy4;
1117 
1118  if( result > DBL_MIN ) lambda = 1./result;
1119  else lambda = DBL_MAX;
1120 
1121  return lambda;
1122 }
G4double G4ParticleHPJENDLHEData::G4double result
int G4int
Definition: G4Types.hh:78
#define DBL_MIN
Definition: templates.hh:75
double G4double
Definition: G4Types.hh:76
#define DBL_MAX
Definition: templates.hh:83
G4double G4PAIxSection::GetPlasmonEnergyTransfer ( )

Definition at line 2416 of file G4PAIxSection.cc.

2417 {
2418  G4int iTransfer ;
2419 
2420  G4double energyTransfer, position;
2421 
2422  position = fIntegralPlasmon[1]*G4UniformRand();
2423 
2424  for( iTransfer = 1; iTransfer <= fSplineNumber; iTransfer++ )
2425  {
2426  if( position >= fIntegralPlasmon[iTransfer] ) break;
2427  }
2428  if(iTransfer > fSplineNumber) iTransfer--;
2429 
2430  energyTransfer = fSplineEnergy[iTransfer];
2431 
2432  if(iTransfer > 1)
2433  {
2434  energyTransfer -= (fSplineEnergy[iTransfer]-fSplineEnergy[iTransfer-1])*G4UniformRand();
2435  }
2436  return energyTransfer;
2437 }
int G4int
Definition: G4Types.hh:78
#define position
Definition: xmlparse.cc:622
#define G4UniformRand()
Definition: Randomize.hh:97
double G4double
Definition: G4Types.hh:76
G4double G4PAIxSection::GetResonanceEnergyTransfer ( )

Definition at line 2471 of file G4PAIxSection.cc.

2472 {
2473  G4int iTransfer ;
2474 
2475  G4double energyTransfer, position;
2476 
2477  position = fIntegralResonance[1]*G4UniformRand();
2478 
2479  for( iTransfer = 1; iTransfer <= fSplineNumber; iTransfer++ )
2480  {
2481  if( position >= fIntegralResonance[iTransfer] ) break;
2482  }
2483  if(iTransfer > fSplineNumber) iTransfer--;
2484 
2485  energyTransfer = fSplineEnergy[iTransfer];
2486 
2487  if(iTransfer > 1)
2488  {
2489  energyTransfer -= (fSplineEnergy[iTransfer]-fSplineEnergy[iTransfer-1])*G4UniformRand();
2490  }
2491  return energyTransfer;
2492 }
int G4int
Definition: G4Types.hh:78
#define position
Definition: xmlparse.cc:622
#define G4UniformRand()
Definition: Randomize.hh:97
double G4double
Definition: G4Types.hh:76
G4double G4PAIxSection::GetRutherfordEnergyTransfer ( )

Definition at line 2499 of file G4PAIxSection.cc.

2500 {
2501  G4int iTransfer ;
2502 
2503  G4double energyTransfer, position;
2504 
2505  position = (fIntegralPlasmon[1]-fIntegralResonance[1])*G4UniformRand();
2506 
2507  for( iTransfer = 1; iTransfer <= fSplineNumber; iTransfer++ )
2508  {
2509  if( position >= (fIntegralPlasmon[iTransfer]-fIntegralResonance[iTransfer]) ) break;
2510  }
2511  if(iTransfer > fSplineNumber) iTransfer--;
2512 
2513  energyTransfer = fSplineEnergy[iTransfer];
2514 
2515  if(iTransfer > 1)
2516  {
2517  energyTransfer -= (fSplineEnergy[iTransfer]-fSplineEnergy[iTransfer-1])*G4UniformRand();
2518  }
2519  return energyTransfer;
2520 }
int G4int
Definition: G4Types.hh:78
#define position
Definition: xmlparse.cc:622
#define G4UniformRand()
Definition: Randomize.hh:97
double G4double
Definition: G4Types.hh:76
G4double G4PAIxSection::GetSplineEnergy ( G4int  i) const
inline

Definition at line 291 of file G4PAIxSection.hh.

292 {
293  if(i < 1 || i > fSplineNumber) { CallError(i, "GetSplineEnergy"); }
294  return fSplineEnergy[i];
295 }
G4int G4PAIxSection::GetSplineSize ( ) const
inline

Definition at line 175 of file G4PAIxSection.hh.

175 { return fSplineNumber; }
G4double G4PAIxSection::GetStepCerenkovLoss ( G4double  step)

Definition at line 2281 of file G4PAIxSection.cc.

2282 {
2283  G4long numOfCollisions;
2284  G4double meanNumber, loss = 0.0;
2285 
2286  // G4cout<<" G4PAIxSection::GetStepCerenkovLoss "<<G4endl;
2287 
2288  meanNumber = fIntegralCerenkov[1]*step;
2289  numOfCollisions = G4Poisson(meanNumber);
2290 
2291  // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl;
2292 
2293  while(numOfCollisions)
2294  {
2295  loss += GetCerenkovEnergyTransfer();
2296  numOfCollisions--;
2297  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
2298  }
2299  // G4cout<<"PAI Cerenkov loss = "<<loss/keV<<" keV"<<G4endl;
2300 
2301  return loss;
2302 }
G4long G4Poisson(G4double mean)
Definition: G4Poisson.hh:51
long G4long
Definition: G4Types.hh:80
G4double GetCerenkovEnergyTransfer()
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

G4double G4PAIxSection::GetStepEnergyLoss ( G4double  step)

Definition at line 2227 of file G4PAIxSection.cc.

2228 {
2229  G4long numOfCollisions;
2230  G4double meanNumber, loss = 0.0;
2231 
2232  // G4cout<<" G4PAIxSection::GetStepEnergyLoss "<<G4endl;
2233 
2234  meanNumber = fIntegralPAIxSection[1]*step;
2235  numOfCollisions = G4Poisson(meanNumber);
2236 
2237  // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl;
2238 
2239  while(numOfCollisions)
2240  {
2241  loss += GetEnergyTransfer();
2242  numOfCollisions--;
2243  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
2244  }
2245  // G4cout<<"PAI energy loss = "<<loss/keV<<" keV"<<G4endl;
2246 
2247  return loss;
2248 }
G4long G4Poisson(G4double mean)
Definition: G4Poisson.hh:51
long G4long
Definition: G4Types.hh:80
double G4double
Definition: G4Types.hh:76
G4double GetEnergyTransfer()

Here is the call graph for this function:

G4double G4PAIxSection::GetStepMMLoss ( G4double  step)

Definition at line 2308 of file G4PAIxSection.cc.

2309 {
2310  G4long numOfCollisions;
2311  G4double meanNumber, loss = 0.0;
2312 
2313  // G4cout<<" G4PAIxSection::GetStepMMLoss "<<G4endl;
2314 
2315  meanNumber = fIntegralMM[1]*step;
2316  numOfCollisions = G4Poisson(meanNumber);
2317 
2318  // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl;
2319 
2320  while(numOfCollisions)
2321  {
2322  loss += GetMMEnergyTransfer();
2323  numOfCollisions--;
2324  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
2325  }
2326  // G4cout<<"PAI MM-Cerenkov loss = "<<loss/keV<<" keV"<<G4endl;
2327 
2328  return loss;
2329 }
G4long G4Poisson(G4double mean)
Definition: G4Poisson.hh:51
long G4long
Definition: G4Types.hh:80
G4double GetMMEnergyTransfer()
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

G4double G4PAIxSection::GetStepPlasmonLoss ( G4double  step)

Definition at line 2389 of file G4PAIxSection.cc.

2390 {
2391  G4long numOfCollisions;
2392  G4double meanNumber, loss = 0.0;
2393 
2394  // G4cout<<" G4PAIxSection::GetStepPlasmonLoss "<<G4endl;
2395 
2396  meanNumber = fIntegralPlasmon[1]*step;
2397  numOfCollisions = G4Poisson(meanNumber);
2398 
2399  // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl;
2400 
2401  while(numOfCollisions)
2402  {
2403  loss += GetPlasmonEnergyTransfer();
2404  numOfCollisions--;
2405  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
2406  }
2407  // G4cout<<"PAI Plasmon loss = "<<loss/keV<<" keV"<<G4endl;
2408 
2409  return loss;
2410 }
G4long G4Poisson(G4double mean)
Definition: G4Poisson.hh:51
G4double GetPlasmonEnergyTransfer()
long G4long
Definition: G4Types.hh:80
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

G4double G4PAIxSection::GetStepResonanceLoss ( G4double  step)

Definition at line 2443 of file G4PAIxSection.cc.

2444 {
2445  G4long numOfCollisions;
2446  G4double meanNumber, loss = 0.0;
2447 
2448  // G4cout<<" G4PAIxSection::GetStepCreLosnkovs "<<G4endl;
2449 
2450  meanNumber = fIntegralResonance[1]*step;
2451  numOfCollisions = G4Poisson(meanNumber);
2452 
2453  // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl;
2454 
2455  while(numOfCollisions)
2456  {
2457  loss += GetResonanceEnergyTransfer();
2458  numOfCollisions--;
2459  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
2460  }
2461  // G4cout<<"PAI resonance loss = "<<loss/keV<<" keV"<<G4endl;
2462 
2463  return loss;
2464 }
G4long G4Poisson(G4double mean)
Definition: G4Poisson.hh:51
long G4long
Definition: G4Types.hh:80
G4double GetResonanceEnergyTransfer()
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

G4double G4PAIxSection::ImPartDielectricConst ( G4int  intervalNumber,
G4double  energy 
)

Definition at line 1076 of file G4PAIxSection.cc.

1078 {
1079  G4double energy2,energy3,energy4,result;
1080 
1081  energy2 = energy1*energy1;
1082  energy3 = energy2*energy1;
1083  energy4 = energy3*energy1;
1084 
1085  result = fA1[k]/energy1+fA2[k]/energy2+fA3[k]/energy3+fA4[k]/energy4;
1086  result *=hbarc/energy1;
1087 
1088  return result;
1089 
1090 } // end of ImPartDielectricConst
G4double G4ParticleHPJENDLHEData::G4double result
double G4double
Definition: G4Types.hh:76
void G4PAIxSection::Initialize ( const G4Material material,
G4double  maxEnergyTransfer,
G4double  betaGammaSq,
G4SandiaTable sandia 
)

Definition at line 606 of file G4PAIxSection.cc.

610 {
611  if(fVerbose > 0)
612  {
613  G4cout<<G4endl;
614  G4cout<<"G4PAIxSection::Initialize(...,G4SandiaTable* sandia)"<<G4endl;
615  G4cout<<G4endl;
616  }
617  G4int i, j;
618 
619  fSandia = sandia;
620  fIntervalNumber = sandia->GetMaxInterval();
621  fDensity = material->GetDensity();
622  fElectronDensity = material->GetElectronDensity();
623 
624  // fIntervalNumber--;
625 
626  if( fVerbose > 0 )
627  {
628  G4cout<<"fDensity = "<<fDensity<<"\t"<<fElectronDensity<<"\t fIntervalNumber = "<<fIntervalNumber<<G4endl;
629  }
630  fEnergyInterval = G4DataVector(fIntervalNumber+2,0.0);
631  fA1 = G4DataVector(fIntervalNumber+2,0.0);
632  fA2 = G4DataVector(fIntervalNumber+2,0.0);
633  fA3 = G4DataVector(fIntervalNumber+2,0.0);
634  fA4 = G4DataVector(fIntervalNumber+2,0.0);
635 
636  for( i = 1; i <= fIntervalNumber; i++ )
637  {
638  if ( sandia->GetSandiaMatTablePAI(i-1,0) < 1.*eV && sandia->GetLowerI1() == false)
639  {
640  fIntervalNumber--;
641  continue;
642  }
643  if( ( sandia->GetSandiaMatTablePAI(i-1,0) >= maxEnergyTransfer ) || i >= fIntervalNumber )
644  {
645  fEnergyInterval[i] = maxEnergyTransfer;
646  fIntervalNumber = i;
647  break;
648  }
649  fEnergyInterval[i] = sandia->GetSandiaMatTablePAI(i-1,0);
650  fA1[i] = sandia->GetSandiaMatTablePAI(i-1,1);
651  fA2[i] = sandia->GetSandiaMatTablePAI(i-1,2);
652  fA3[i] = sandia->GetSandiaMatTablePAI(i-1,3);
653  fA4[i] = sandia->GetSandiaMatTablePAI(i-1,4);
654 
655  if( fVerbose > 0 )
656  {
657  G4cout<<i<<"\t"<<fEnergyInterval[i]/keV<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
658  <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl;
659  }
660  }
661  if( fVerbose > 0 ) G4cout<<"last i = "<<i<<"; "<<"fIntervalNumber = "<<fIntervalNumber<<G4endl;
662 
663  if( fEnergyInterval[fIntervalNumber] != maxEnergyTransfer )
664  {
665  fIntervalNumber++;
666  fEnergyInterval[fIntervalNumber] = maxEnergyTransfer;
667  }
668  if( fVerbose > 0 )
669  {
670  for( i = 1; i <= fIntervalNumber; i++ )
671  {
672  G4cout<<i<<"\t"<<fEnergyInterval[i]/keV<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
673  <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl;
674  }
675  }
676  if( fVerbose > 0 ) G4cout<<"Now checking, if two borders are too close together"<<G4endl;
677 
678  for( i = 1; i < fIntervalNumber; i++ )
679  {
680  if( fEnergyInterval[i+1]-fEnergyInterval[i] >
681  1.5*fDelta*(fEnergyInterval[i+1]+fEnergyInterval[i]) && fEnergyInterval[i] > 0.) continue;
682  else
683  {
684  if( fVerbose > 0 ) G4cout<<i<<"\t"<<fEnergyInterval[i]/keV<<"\t"<<fEnergyInterval[i+1]/keV;
685 
686  for( j = i; j < fIntervalNumber; j++ )
687  {
688  fEnergyInterval[j] = fEnergyInterval[j+1];
689  fA1[j] = fA1[j+1];
690  fA2[j] = fA2[j+1];
691  fA3[j] = fA3[j+1];
692  fA4[j] = fA4[j+1];
693  }
694  fIntervalNumber--;
695  i--;
696  }
697  }
698  if( fVerbose > 0 )
699  {
700  for( i = 1; i <= fIntervalNumber; i++ )
701  {
702  G4cout<<i<<"\t"<<fEnergyInterval[i]/keV<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
703  <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl;
704  }
705  }
706  // Preparation of fSplineEnergy array corresponding to min ionisation, G~4
707 
708  ComputeLowEnergyCof(material);
709 
710  G4double betaGammaSqRef =
711  fLorentzFactor[fRefGammaNumber]*fLorentzFactor[fRefGammaNumber] - 1;
712 
713  NormShift(betaGammaSqRef);
714  SplainPAI(betaGammaSqRef);
715 
716  // Preparation of integral PAI cross section for input betaGammaSq
717 
718  for( i = 1; i <= fSplineNumber; i++ )
719  {
720  fDifPAIxSection[i] = DifPAIxSection(i,betaGammaSq);
721 
722 
723  fdNdxCerenkov[i] = PAIdNdxCerenkov(i,betaGammaSq);
724  fdNdxMM[i] = PAIdNdxMM(i,betaGammaSq);
725  fdNdxPlasmon[i] = PAIdNdxPlasmon(i,betaGammaSq);
726  fdNdxResonance[i] = PAIdNdxResonance(i,betaGammaSq);
727  }
730  IntegralMM();
731  IntegralPlasmon();
733 
734  for( i = 1; i <= fSplineNumber; i++ )
735  {
736  if(fVerbose>0) G4cout<<i<<"; w = "<<fSplineEnergy[i]/keV<<" keV; dN/dx_>w = "<<fIntegralPAIxSection[i]<<" 1/mm"<<G4endl;
737  }
738 }
void IntegralPAIxSection()
void SplainPAI(G4double betaGammaSq)
G4double PAIdNdxResonance(G4int intervalNumber, G4double betaGammaSq)
void IntegralCerenkov()
void IntegralPlasmon()
void IntegralResonance()
G4double PAIdNdxCerenkov(G4int intervalNumber, G4double betaGammaSq)
G4double GetDensity() const
Definition: G4Material.hh:180
void NormShift(G4double betaGammaSq)
G4int GetMaxInterval() const
G4double PAIdNdxPlasmon(G4int intervalNumber, G4double betaGammaSq)
G4double GetSandiaMatTablePAI(G4int, G4int) const
void ComputeLowEnergyCof()
int G4int
Definition: G4Types.hh:78
G4GLOB_DLL std::ostream G4cout
G4double GetElectronDensity() const
Definition: G4Material.hh:217
G4bool GetLowerI1()
static constexpr double eV
Definition: G4SIunits.hh:215
G4double PAIdNdxMM(G4int intervalNumber, G4double betaGammaSq)
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
static constexpr double keV
Definition: G4SIunits.hh:216
G4double DifPAIxSection(G4int intervalNumber, G4double betaGammaSq)

Here is the call graph for this function:

void G4PAIxSection::InitPAI ( )

Definition at line 813 of file G4PAIxSection.cc.

814 {
815  G4int i;
816  G4double betaGammaSq = fLorentzFactor[fRefGammaNumber]*
817  fLorentzFactor[fRefGammaNumber] - 1;
818 
819  // Preparation of integral PAI cross section for reference gamma
820 
821  NormShift(betaGammaSq);
822  SplainPAI(betaGammaSq);
823 
826  IntegralMM();
827  IntegralPlasmon();
829 
830  for(i = 0; i<= fSplineNumber; i++)
831  {
832  fPAItable[i][fRefGammaNumber] = fIntegralPAIxSection[i];
833  if(i != 0)
834  {
835  fPAItable[i][0] = fSplineEnergy[i];
836  }
837  }
838  fPAItable[0][0] = fSplineNumber;
839 
840  for(G4int j = 1; j < 112; j++) // for other gammas
841  {
842  if( j == fRefGammaNumber ) continue;
843 
844  betaGammaSq = fLorentzFactor[j]*fLorentzFactor[j] - 1;
845 
846  for(i = 1; i <= fSplineNumber; i++)
847  {
848  fDifPAIxSection[i] = DifPAIxSection(i,betaGammaSq);
849  fdNdxCerenkov[i] = PAIdNdxCerenkov(i,betaGammaSq);
850  fdNdxMM[i] = PAIdNdxMM(i,betaGammaSq);
851  fdNdxPlasmon[i] = PAIdNdxPlasmon(i,betaGammaSq);
852  fdNdxResonance[i] = PAIdNdxResonance(i,betaGammaSq);
853  }
856  IntegralMM();
857  IntegralPlasmon();
859 
860  for(i = 0; i <= fSplineNumber; i++)
861  {
862  fPAItable[i][j] = fIntegralPAIxSection[i];
863  }
864  }
865 
866 }
void IntegralPAIxSection()
void SplainPAI(G4double betaGammaSq)
G4double PAIdNdxResonance(G4int intervalNumber, G4double betaGammaSq)
void IntegralCerenkov()
void IntegralPlasmon()
void IntegralResonance()
G4double PAIdNdxCerenkov(G4int intervalNumber, G4double betaGammaSq)
void NormShift(G4double betaGammaSq)
G4double PAIdNdxPlasmon(G4int intervalNumber, G4double betaGammaSq)
int G4int
Definition: G4Types.hh:78
G4double PAIdNdxMM(G4int intervalNumber, G4double betaGammaSq)
double G4double
Definition: G4Types.hh:76
G4double DifPAIxSection(G4int intervalNumber, G4double betaGammaSq)
void G4PAIxSection::IntegralCerenkov ( )

Definition at line 1522 of file G4PAIxSection.cc.

1523 {
1524  G4int i, k;
1525  fIntegralCerenkov[fSplineNumber] = 0;
1526  fIntegralCerenkov[0] = 0;
1527  k = fIntervalNumber -1;
1528 
1529  for( i = fSplineNumber-1; i >= 1; i-- )
1530  {
1531  if(fSplineEnergy[i] >= fEnergyInterval[k])
1532  {
1533  fIntegralCerenkov[i] = fIntegralCerenkov[i+1] + SumOverInterCerenkov(i);
1534  // G4cout<<"int: i = "<<i<<"; sumC = "<<fIntegralCerenkov[i]<<G4endl;
1535  }
1536  else
1537  {
1538  fIntegralCerenkov[i] = fIntegralCerenkov[i+1] +
1539  SumOverBordCerenkov(i+1,fEnergyInterval[k]);
1540  k--;
1541  // G4cout<<"bord: i = "<<i<<"; sumC = "<<fIntegralCerenkov[i]<<G4endl;
1542  }
1543  }
1544 
1545 } // end of IntegralCerenkov
int G4int
Definition: G4Types.hh:78
G4double SumOverBordCerenkov(G4int intervalNumber, G4double energy)
G4double SumOverInterCerenkov(G4int intervalNumber)
void G4PAIxSection::IntegralMM ( )

Definition at line 1553 of file G4PAIxSection.cc.

1554 {
1555  G4int i, k;
1556  fIntegralMM[fSplineNumber] = 0;
1557  fIntegralMM[0] = 0;
1558  k = fIntervalNumber -1;
1559 
1560  for( i = fSplineNumber-1; i >= 1; i-- )
1561  {
1562  if(fSplineEnergy[i] >= fEnergyInterval[k])
1563  {
1564  fIntegralMM[i] = fIntegralMM[i+1] + SumOverInterMM(i);
1565  // G4cout<<"int: i = "<<i<<"; sumC = "<<fIntegralMM[i]<<G4endl;
1566  }
1567  else
1568  {
1569  fIntegralMM[i] = fIntegralMM[i+1] +
1570  SumOverBordMM(i+1,fEnergyInterval[k]);
1571  k--;
1572  // G4cout<<"bord: i = "<<i<<"; sumC = "<<fIntegralMM[i]<<G4endl;
1573  }
1574  }
1575 
1576 } // end of IntegralMM
G4double SumOverInterMM(G4int intervalNumber)
int G4int
Definition: G4Types.hh:78
G4double SumOverBordMM(G4int intervalNumber, G4double energy)
void G4PAIxSection::IntegralPAIxSection ( )

Definition at line 1490 of file G4PAIxSection.cc.

1491 {
1492  fIntegralPAIxSection[fSplineNumber] = 0;
1493  fIntegralPAIdEdx[fSplineNumber] = 0;
1494  fIntegralPAIxSection[0] = 0;
1495  G4int i, k = fIntervalNumber -1;
1496 
1497  for( i = fSplineNumber-1; i >= 1; i--)
1498  {
1499  if(fSplineEnergy[i] >= fEnergyInterval[k])
1500  {
1501  fIntegralPAIxSection[i] = fIntegralPAIxSection[i+1] + SumOverInterval(i);
1502  fIntegralPAIdEdx[i] = fIntegralPAIdEdx[i+1] + SumOverIntervaldEdx(i);
1503  }
1504  else
1505  {
1506  fIntegralPAIxSection[i] = fIntegralPAIxSection[i+1] +
1507  SumOverBorder(i+1,fEnergyInterval[k]);
1508  fIntegralPAIdEdx[i] = fIntegralPAIdEdx[i+1] +
1509  SumOverBorderdEdx(i+1,fEnergyInterval[k]);
1510  k--;
1511  }
1512  if(fVerbose>0) G4cout<<"i = "<<i<<"; k = "<<k<<"; intPAIxsc[i] = "<<fIntegralPAIxSection[i]<<G4endl;
1513  }
1514 } // end of IntegralPAIxSection
G4double SumOverInterval(G4int intervalNumber)
int G4int
Definition: G4Types.hh:78
G4GLOB_DLL std::ostream G4cout
G4double SumOverBorder(G4int intervalNumber, G4double energy)
#define G4endl
Definition: G4ios.hh:61
G4double SumOverIntervaldEdx(G4int intervalNumber)
G4double SumOverBorderdEdx(G4int intervalNumber, G4double energy)
void G4PAIxSection::IntegralPlasmon ( )

Definition at line 1584 of file G4PAIxSection.cc.

1585 {
1586  fIntegralPlasmon[fSplineNumber] = 0;
1587  fIntegralPlasmon[0] = 0;
1588  G4int k = fIntervalNumber -1;
1589  for(G4int i=fSplineNumber-1;i>=1;i--)
1590  {
1591  if(fSplineEnergy[i] >= fEnergyInterval[k])
1592  {
1593  fIntegralPlasmon[i] = fIntegralPlasmon[i+1] + SumOverInterPlasmon(i);
1594  }
1595  else
1596  {
1597  fIntegralPlasmon[i] = fIntegralPlasmon[i+1] +
1598  SumOverBordPlasmon(i+1,fEnergyInterval[k]);
1599  k--;
1600  }
1601  }
1602 
1603 } // end of IntegralPlasmon
G4double SumOverBordPlasmon(G4int intervalNumber, G4double energy)
G4double SumOverInterPlasmon(G4int intervalNumber)
int G4int
Definition: G4Types.hh:78
void G4PAIxSection::IntegralResonance ( )

Definition at line 1611 of file G4PAIxSection.cc.

1612 {
1613  fIntegralResonance[fSplineNumber] = 0;
1614  fIntegralResonance[0] = 0;
1615  G4int k = fIntervalNumber -1;
1616  for(G4int i=fSplineNumber-1;i>=1;i--)
1617  {
1618  if(fSplineEnergy[i] >= fEnergyInterval[k])
1619  {
1620  fIntegralResonance[i] = fIntegralResonance[i+1] + SumOverInterResonance(i);
1621  }
1622  else
1623  {
1624  fIntegralResonance[i] = fIntegralResonance[i+1] +
1625  SumOverBordResonance(i+1,fEnergyInterval[k]);
1626  k--;
1627  }
1628  }
1629 
1630 } // end of IntegralResonance
G4double SumOverBordResonance(G4int intervalNumber, G4double energy)
int G4int
Definition: G4Types.hh:78
G4double SumOverInterResonance(G4int intervalNumber)
void G4PAIxSection::NormShift ( G4double  betaGammaSq)

Definition at line 873 of file G4PAIxSection.cc.

874 {
875  G4int i, j;
876 
877  if(fVerbose>0) G4cout<<" G4PAIxSection::NormShift call "<<G4endl;
878 
879 
880  for( i = 1; i <= fIntervalNumber-1; i++ )
881  {
882  for( j = 1; j <= 2; j++ )
883  {
884  fSplineNumber = (i-1)*2 + j;
885 
886  if( j == 1 ) fSplineEnergy[fSplineNumber] = fEnergyInterval[i ]*(1+fDelta);
887  else fSplineEnergy[fSplineNumber] = fEnergyInterval[i+1]*(1-fDelta);
888  if(fVerbose>0) G4cout<<"cn = "<<fSplineNumber<<"; "<<"w = "<<fSplineEnergy[fSplineNumber]/keV<<" keV"<<G4endl;
889  }
890  }
891  fIntegralTerm[1]=RutherfordIntegral(1,fEnergyInterval[1],fSplineEnergy[1]);
892 
893  j = 1;
894 
895  for( i = 2; i <= fSplineNumber; i++ )
896  {
897  if( fSplineEnergy[i]<fEnergyInterval[j+1] )
898  {
899  fIntegralTerm[i] = fIntegralTerm[i-1] +
900  RutherfordIntegral(j,fSplineEnergy[i-1],
901  fSplineEnergy[i] );
902  }
903  else
904  {
905  G4double x = RutherfordIntegral(j,fSplineEnergy[i-1],
906  fEnergyInterval[j+1] );
907  j++;
908  fIntegralTerm[i] = fIntegralTerm[i-1] + x +
909  RutherfordIntegral(j,fEnergyInterval[j],
910  fSplineEnergy[i] );
911  }
912  if(fVerbose>0) G4cout<<i<<" Shift: w = "<<fSplineEnergy[i]/keV<<" keV \t"<<fIntegralTerm[i]<<"\n"<<G4endl;
913  }
914  fNormalizationCof = 2*pi*pi*hbarc*hbarc*fine_structure_const/electron_mass_c2;
915  fNormalizationCof *= fElectronDensity/fIntegralTerm[fSplineNumber];
916 
917  // G4cout<<"fNormalizationCof = "<<fNormalizationCof<<G4endl;
918 
919  // Calculation of PAI differrential cross-section (1/(keV*cm))
920  // in the energy points near borders of energy intervals
921 
922  for(G4int k = 1; k <= fIntervalNumber-1; k++ )
923  {
924  for( j = 1; j <= 2; j++ )
925  {
926  i = (k-1)*2 + j;
927  fImPartDielectricConst[i] = fNormalizationCof*
928  ImPartDielectricConst(k,fSplineEnergy[i]);
929  fRePartDielectricConst[i] = fNormalizationCof*
930  RePartDielectricConst(fSplineEnergy[i]);
931  fIntegralTerm[i] *= fNormalizationCof;
932 
933  fDifPAIxSection[i] = DifPAIxSection(i,betaGammaSq);
934  fdNdxCerenkov[i] = PAIdNdxCerenkov(i,betaGammaSq);
935  fdNdxMM[i] = PAIdNdxMM(i,betaGammaSq);
936  fdNdxPlasmon[i] = PAIdNdxPlasmon(i,betaGammaSq);
937  fdNdxResonance[i] = PAIdNdxResonance(i,betaGammaSq);
938  if(fVerbose>0) G4cout<<i<<" Shift: w = "<<fSplineEnergy[i]/keV<<" keV, xsc = "<<fDifPAIxSection[i]<<"\n"<<G4endl;
939  }
940  }
941 
942 } // end of NormShift
G4double PAIdNdxResonance(G4int intervalNumber, G4double betaGammaSq)
G4double PAIdNdxCerenkov(G4int intervalNumber, G4double betaGammaSq)
G4double PAIdNdxPlasmon(G4int intervalNumber, G4double betaGammaSq)
tuple x
Definition: test.py:50
G4double RePartDielectricConst(G4double energy)
int G4int
Definition: G4Types.hh:78
G4GLOB_DLL std::ostream G4cout
G4double ImPartDielectricConst(G4int intervalNumber, G4double energy)
float electron_mass_c2
Definition: hepunit.py:274
G4double PAIdNdxMM(G4int intervalNumber, G4double betaGammaSq)
#define G4endl
Definition: G4ios.hh:61
static constexpr double pi
Definition: G4SIunits.hh:75
double G4double
Definition: G4Types.hh:76
static constexpr double keV
Definition: G4SIunits.hh:216
G4double DifPAIxSection(G4int intervalNumber, G4double betaGammaSq)
G4double RutherfordIntegral(G4int intervalNumber, G4double limitLow, G4double limitHigh)
G4double G4PAIxSection::PAIdNdxCerenkov ( G4int  intervalNumber,
G4double  betaGammaSq 
)

Definition at line 1302 of file G4PAIxSection.cc.

1304 {
1305  G4double logarithm, x3, x5, argument, modul2, dNdxC;
1306  G4double be2, betaBohr2, cofBetaBohr;
1307 
1308  cofBetaBohr = 4.0;
1310  G4double betaBohr4 = betaBohr2*betaBohr2*cofBetaBohr;
1311 
1312  be2 = betaGammaSq/(1 + betaGammaSq);
1313  G4double be4 = be2*be2;
1314 
1315  if( betaGammaSq < 0.01 ) logarithm = log(1.0+betaGammaSq); // 0.0;
1316  else
1317  {
1318  logarithm = -log( (1/betaGammaSq - fRePartDielectricConst[i])*
1319  (1/betaGammaSq - fRePartDielectricConst[i]) +
1320  fImPartDielectricConst[i]*fImPartDielectricConst[i] )*0.5;
1321  logarithm += log(1+1.0/betaGammaSq);
1322  }
1323 
1324  if( fImPartDielectricConst[i] == 0.0 || betaGammaSq < 0.01 )
1325  {
1326  argument = 0.0;
1327  }
1328  else
1329  {
1330  x3 = -fRePartDielectricConst[i] + 1.0/betaGammaSq;
1331  x5 = -1.0 - fRePartDielectricConst[i] +
1332  be2*((1.0 +fRePartDielectricConst[i])*(1.0 + fRePartDielectricConst[i]) +
1333  fImPartDielectricConst[i]*fImPartDielectricConst[i]);
1334  if( x3 == 0.0 ) argument = 0.5*pi;
1335  else argument = atan2(fImPartDielectricConst[i],x3);
1336  argument *= x5 ;
1337  }
1338  dNdxC = ( logarithm*fImPartDielectricConst[i] + argument )/hbarc;
1339 
1340  if(dNdxC < 1.0e-8) dNdxC = 1.0e-8;
1341 
1342  dNdxC *= fine_structure_const/be2/pi;
1343 
1344  dNdxC *= (1-exp(-be4/betaBohr4));
1345 
1346  if(fDensity >= 0.1)
1347  {
1348  modul2 = (1.0 + fRePartDielectricConst[i])*(1.0 + fRePartDielectricConst[i]) +
1349  fImPartDielectricConst[i]*fImPartDielectricConst[i];
1350  dNdxC /= modul2;
1351  }
1352  return dNdxC;
1353 
1354 } // end of PAIdNdxCerenkov
static constexpr double pi
Definition: G4SIunits.hh:75
double G4double
Definition: G4Types.hh:76
G4double G4PAIxSection::PAIdNdxMM ( G4int  intervalNumber,
G4double  betaGammaSq 
)

Definition at line 1360 of file G4PAIxSection.cc.

1362 {
1363  G4double logarithm, x3, x5, argument, dNdxC;
1364  G4double be2, be4, betaBohr2,betaBohr4,cofBetaBohr;
1365 
1366  cofBetaBohr = 4.0;
1368  betaBohr4 = betaBohr2*betaBohr2*cofBetaBohr;
1369 
1370  be2 = betaGammaSq/(1 + betaGammaSq);
1371  be4 = be2*be2;
1372 
1373  if( betaGammaSq < 0.01 ) logarithm = log(1.0+betaGammaSq); // 0.0;
1374  else
1375  {
1376  logarithm = -log( (1/betaGammaSq - fRePartDielectricConst[i])*
1377  (1/betaGammaSq - fRePartDielectricConst[i]) +
1378  fImPartDielectricConst[i]*fImPartDielectricConst[i] )*0.5;
1379  logarithm += log(1+1.0/betaGammaSq);
1380  }
1381 
1382  if( fImPartDielectricConst[i] == 0.0 || betaGammaSq < 0.01 )
1383  {
1384  argument = 0.0;
1385  }
1386  else
1387  {
1388  x3 = -fRePartDielectricConst[i] + 1.0/betaGammaSq;
1389  x5 = be2*( 1.0 + fRePartDielectricConst[i] ) - 1.0;
1390  if( x3 == 0.0 ) argument = 0.5*pi;
1391  else argument = atan2(fImPartDielectricConst[i],x3);
1392  argument *= x5 ;
1393  }
1394  dNdxC = ( logarithm*fImPartDielectricConst[i]*be2 + argument )/hbarc;
1395 
1396  if(dNdxC < 1.0e-8) dNdxC = 1.0e-8;
1397 
1398  dNdxC *= fine_structure_const/be2/pi;
1399 
1400  dNdxC *= (1-exp(-be4/betaBohr4));
1401  return dNdxC;
1402 
1403 } // end of PAIdNdxMM
static constexpr double pi
Definition: G4SIunits.hh:75
double G4double
Definition: G4Types.hh:76
G4double G4PAIxSection::PAIdNdxPlasmon ( G4int  intervalNumber,
G4double  betaGammaSq 
)

Definition at line 1410 of file G4PAIxSection.cc.

1412 {
1413  G4double resonance, modul2, dNdxP, cof = 1.;
1414  G4double be2, betaBohr;
1415 
1416  betaBohr = fine_structure_const;
1417  be2 = betaGammaSq/(1 + betaGammaSq);
1418 
1419  G4double beta = sqrt(be2);
1420 
1421  resonance = log(2*electron_mass_c2*be2/fSplineEnergy[i]);
1422  resonance *= fImPartDielectricConst[i]/hbarc;
1423 
1424 
1425  dNdxP = ( resonance + cof*fIntegralTerm[i]/fSplineEnergy[i]/fSplineEnergy[i] );
1426 
1427  if( dNdxP < 1.0e-8 ) dNdxP = 1.0e-8;
1428 
1429  dNdxP *= fine_structure_const/be2/pi;
1430 
1431  dNdxP *= (1 - exp(-beta/betaBohr/fLowEnergyCof));
1432 
1433  // dNdxP *= (1-exp(-be4/betaBohr4));
1434 
1435  if( fDensity >= 0.1 )
1436  {
1437  modul2 = (1 + fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) +
1438  fImPartDielectricConst[i]*fImPartDielectricConst[i];
1439  dNdxP /= modul2;
1440  }
1441  return dNdxP;
1442 
1443 } // end of PAIdNdxPlasmon
float electron_mass_c2
Definition: hepunit.py:274
static constexpr double pi
Definition: G4SIunits.hh:75
double G4double
Definition: G4Types.hh:76
G4double G4PAIxSection::PAIdNdxResonance ( G4int  intervalNumber,
G4double  betaGammaSq 
)

Definition at line 1450 of file G4PAIxSection.cc.

1452 {
1453  G4double resonance, modul2, dNdxP;
1454  G4double be2, be4, betaBohr2, betaBohr4, cofBetaBohr;
1455 
1456  cofBetaBohr = 4.0;
1458  betaBohr4 = betaBohr2*betaBohr2*cofBetaBohr;
1459 
1460  be2 = betaGammaSq/(1 + betaGammaSq);
1461  be4 = be2*be2;
1462 
1463  resonance = log(2*electron_mass_c2*be2/fSplineEnergy[i]);
1464  resonance *= fImPartDielectricConst[i]/hbarc;
1465 
1466 
1467  dNdxP = resonance;
1468 
1469  if( dNdxP < 1.0e-8 ) dNdxP = 1.0e-8;
1470 
1471  dNdxP *= fine_structure_const/be2/pi;
1472  dNdxP *= (1-exp(-be4/betaBohr4));
1473 
1474  if( fDensity >= 0.1 )
1475  {
1476  modul2 = (1 + fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) +
1477  fImPartDielectricConst[i]*fImPartDielectricConst[i];
1478  dNdxP /= modul2;
1479  }
1480  return dNdxP;
1481 
1482 } // end of PAIdNdxResonance
float electron_mass_c2
Definition: hepunit.py:274
static constexpr double pi
Definition: G4SIunits.hh:75
double G4double
Definition: G4Types.hh:76
G4double G4PAIxSection::RePartDielectricConst ( G4double  energy)

Definition at line 1171 of file G4PAIxSection.cc.

1172 {
1173  G4double x0, x02, x03, x04, x05, x1, x2, xx1 ,xx2 , xx12,
1174  c1, c2, c3, cof1, cof2, xln1, xln2, xln3, result;
1175 
1176  x0 = enb;
1177  result = 0;
1178 
1179  for(G4int i=1;i<=fIntervalNumber-1;i++)
1180  {
1181  x1 = fEnergyInterval[i];
1182  x2 = fEnergyInterval[i+1];
1183  xx1 = x1 - x0;
1184  xx2 = x2 - x0;
1185  xx12 = xx2/xx1;
1186 
1187  if(xx12<0)
1188  {
1189  xx12 = -xx12;
1190  }
1191  xln1 = log(x2/x1);
1192  xln2 = log(xx12);
1193  xln3 = log((x2 + x0)/(x1 + x0));
1194  x02 = x0*x0;
1195  x03 = x02*x0;
1196  x04 = x03*x0;
1197  x05 = x04*x0;
1198  c1 = (x2 - x1)/x1/x2;
1199  c2 = (x2 - x1)*(x2 +x1)/x1/x1/x2/x2;
1200  c3 = (x2 -x1)*(x1*x1 + x1*x2 + x2*x2)/x1/x1/x1/x2/x2/x2;
1201 
1202  result -= (fA1[i]/x02 + fA3[i]/x04)*xln1;
1203  result -= (fA2[i]/x02 + fA4[i]/x04)*c1;
1204  result -= fA3[i]*c2/2/x02;
1205  result -= fA4[i]*c3/3/x02;
1206 
1207  cof1 = fA1[i]/x02 + fA3[i]/x04;
1208  cof2 = fA2[i]/x03 + fA4[i]/x05;
1209 
1210  result += 0.5*(cof1 +cof2)*xln2;
1211  result += 0.5*(cof1 - cof2)*xln3;
1212  }
1213  result *= 2*hbarc/pi;
1214 
1215  return result;
1216 
1217 } // end of RePartDielectricConst
G4double G4ParticleHPJENDLHEData::G4double result
static c2_factory< G4double > c2
int G4int
Definition: G4Types.hh:78
static constexpr double pi
Definition: G4SIunits.hh:75
double G4double
Definition: G4Types.hh:76
tuple c1
Definition: plottest35.py:14
G4double G4PAIxSection::RutherfordIntegral ( G4int  intervalNumber,
G4double  limitLow,
G4double  limitHigh 
)

Definition at line 1055 of file G4PAIxSection.cc.

1058 {
1059  G4double c1, c2, c3;
1060  // G4cout<<"RI: x1 = "<<x1<<"; "<<"x2 = "<<x2<<G4endl;
1061  c1 = (x2 - x1)/x1/x2;
1062  c2 = (x2 - x1)*(x2 + x1)/x1/x1/x2/x2;
1063  c3 = (x2 - x1)*(x1*x1 + x1*x2 + x2*x2)/x1/x1/x1/x2/x2/x2;
1064  // G4cout<<" RI: c1 = "<<c1<<"; "<<"c2 = "<<c2<<"; "<<"c3 = "<<c3<<G4endl;
1065 
1066  return fA1[k]*log(x2/x1) + fA2[k]*c1 + fA3[k]*c2/2 + fA4[k]*c3/3;
1067 
1068 } // end of RutherfordIntegral
static c2_factory< G4double > c2
double G4double
Definition: G4Types.hh:76
tuple c1
Definition: plottest35.py:14
void G4PAIxSection::SetVerbose ( G4int  v)
inline

Definition at line 197 of file G4PAIxSection.hh.

197 {fVerbose=v;};
tuple v
Definition: test.py:18
void G4PAIxSection::SplainPAI ( G4double  betaGammaSq)

Definition at line 950 of file G4PAIxSection.cc.

951 {
952  G4int j, k = 1, i = 1;
953 
954  if(fVerbose>0) G4cout<<" G4PAIxSection::SplainPAI call "<<G4endl;
955 
956  while ( (i < fSplineNumber) && (fSplineNumber < fMaxSplineSize-1) )
957  {
958  // if( std::abs(fSplineEnergy[i+1]-fEnergyInterval[k+1]) > (fSplineEnergy[i+1]+fEnergyInterval[k+1])*5.e-7 )
959  if( fSplineEnergy[i+1] > fEnergyInterval[k+1] )
960  {
961  k++; // Here next energy point is in next energy interval
962  i++;
963  if(fVerbose>0) G4cout<<" in if: i = "<<i<<"; k = "<<k<<G4endl;
964  continue;
965  }
966  if(fVerbose>0) G4cout<<" out if: i = "<<i<<"; k = "<<k<<G4endl;
967 
968  // Shifting of arrayes for inserting the geometrical
969  // average of 'i' and 'i+1' energy points to 'i+1' place
970  fSplineNumber++;
971 
972  for( j = fSplineNumber; j >= i+2; j-- )
973  {
974  fSplineEnergy[j] = fSplineEnergy[j-1];
975  fImPartDielectricConst[j] = fImPartDielectricConst[j-1];
976  fRePartDielectricConst[j] = fRePartDielectricConst[j-1];
977  fIntegralTerm[j] = fIntegralTerm[j-1];
978 
979  fDifPAIxSection[j] = fDifPAIxSection[j-1];
980  fdNdxCerenkov[j] = fdNdxCerenkov[j-1];
981  fdNdxMM[j] = fdNdxMM[j-1];
982  fdNdxPlasmon[j] = fdNdxPlasmon[j-1];
983  fdNdxResonance[j] = fdNdxResonance[j-1];
984  }
985  G4double x1 = fSplineEnergy[i];
986  G4double x2 = fSplineEnergy[i+1];
987  G4double yy1 = fDifPAIxSection[i];
988  G4double y2 = fDifPAIxSection[i+1];
989 
990  if(fVerbose>0) G4cout<<"Spline: x1 = "<<x1<<"; x2 = "<<x2<<", yy1 = "<<yy1<<"; y2 = "<<y2<<G4endl;
991 
992 
993  G4double en1 = sqrt(x1*x2);
994  // G4double en1 = 0.5*(x1 + x2);
995 
996 
997  fSplineEnergy[i+1] = en1;
998 
999  // Calculation of logarithmic linear approximation
1000  // in this (enr) energy point, which number is 'i+1' now
1001 
1002  G4double a = log10(y2/yy1)/log10(x2/x1);
1003  G4double b = log10(yy1) - a*log10(x1);
1004  G4double y = a*log10(en1) + b;
1005 
1006  y = pow(10.,y);
1007 
1008  // Calculation of the PAI dif. cross-section at this point
1009 
1010  fImPartDielectricConst[i+1] = fNormalizationCof*
1011  ImPartDielectricConst(k,fSplineEnergy[i+1]);
1012  fRePartDielectricConst[i+1] = fNormalizationCof*
1013  RePartDielectricConst(fSplineEnergy[i+1]);
1014  fIntegralTerm[i+1] = fIntegralTerm[i] + fNormalizationCof*
1015  RutherfordIntegral(k,fSplineEnergy[i],
1016  fSplineEnergy[i+1]);
1017 
1018  fDifPAIxSection[i+1] = DifPAIxSection(i+1,betaGammaSq);
1019  fdNdxCerenkov[i+1] = PAIdNdxCerenkov(i+1,betaGammaSq);
1020  fdNdxMM[i+1] = PAIdNdxMM(i+1,betaGammaSq);
1021  fdNdxPlasmon[i+1] = PAIdNdxPlasmon(i+1,betaGammaSq);
1022  fdNdxResonance[i+1] = PAIdNdxResonance(i+1,betaGammaSq);
1023 
1024  // Condition for next division of this segment or to pass
1025 
1026  if(fVerbose>0) G4cout<<"Spline, a = "<<a<<"; b = "<<b<<"; new xsc = "<<y<<"; compxsc = "<<fDifPAIxSection[i+1]<<G4endl;
1027 
1028  // to higher energies
1029 
1030  G4double x = 2*(fDifPAIxSection[i+1] - y)/(fDifPAIxSection[i+1] + y);
1031 
1032  G4double delta = 2.*(fSplineEnergy[i+1]-fSplineEnergy[i])/(fSplineEnergy[i+1]+fSplineEnergy[i]);
1033 
1034  if( x < 0 )
1035  {
1036  x = -x;
1037  }
1038  if( x > fError && fSplineNumber < fMaxSplineSize-1 && delta > 2.*fDelta )
1039  {
1040  continue; // next division
1041  }
1042  i += 2; // pass to next segment
1043 
1044  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
1045  } // close 'while'
1046 
1047 } // end of SplainPAI
G4double PAIdNdxResonance(G4int intervalNumber, G4double betaGammaSq)
G4double PAIdNdxCerenkov(G4int intervalNumber, G4double betaGammaSq)
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
G4double PAIdNdxPlasmon(G4int intervalNumber, G4double betaGammaSq)
tuple x
Definition: test.py:50
G4double RePartDielectricConst(G4double energy)
int G4int
Definition: G4Types.hh:78
tuple b
Definition: test.py:12
G4GLOB_DLL std::ostream G4cout
G4double ImPartDielectricConst(G4int intervalNumber, G4double energy)
G4double PAIdNdxMM(G4int intervalNumber, G4double betaGammaSq)
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
G4double DifPAIxSection(G4int intervalNumber, G4double betaGammaSq)
G4double RutherfordIntegral(G4int intervalNumber, G4double limitLow, G4double limitHigh)
G4double G4PAIxSection::SumOverBordCerenkov ( G4int  intervalNumber,
G4double  energy 
)

Definition at line 1990 of file G4PAIxSection.cc.

1992 {
1993  G4double x0,x1,y0,yy1,a,b,e0,c,d,result;
1994 
1995  e0 = en0;
1996  x0 = fSplineEnergy[i];
1997  x1 = fSplineEnergy[i+1];
1998  y0 = fdNdxCerenkov[i];
1999  yy1 = fdNdxCerenkov[i+1];
2000 
2001  // G4cout<<G4endl;
2002  // G4cout<<"SumBordC, i = "<<i<<"; en0 = "<<en0<<"; x0 ="<<x0<<"; x1 = "<<x1
2003  // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
2004  c = x1/x0;
2005  d = e0/x0;
2006  a = log10(yy1/y0)/log10(c);
2007  if(a > 10.0) return 0.;
2008  // b0 = log10(y0) - a*log10(x0);
2009  b = y0/pow(x0,a); // pow(10.,b0);
2010 
2011  a += 1.0;
2012  if( a == 0 ) result = b*log(x0/e0);
2013  else result = y0*(x0 - e0*pow(d,a-1))/a;
2014  a += 1.0;
2015 
2016  if( a == 0 ) fIntegralCerenkov[0] += b*log(x0/e0);
2017  else fIntegralCerenkov[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
2018 
2019 // G4cout<<"a = "<<a<<"; b0 = "<<b0<<"; b = "<<b<<"; result = "<<result<<G4endl;
2020 
2021  x0 = fSplineEnergy[i - 1];
2022  x1 = fSplineEnergy[i - 2];
2023  y0 = fdNdxCerenkov[i - 1];
2024  yy1 = fdNdxCerenkov[i - 2];
2025 
2026  // G4cout<<"x0 ="<<x0<<"; x1 = "<<x1
2027  // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
2028 
2029  c = x1/x0;
2030  d = e0/x0;
2031  a = log10(yy1/y0)/log10(x1/x0);
2032  // b0 = log10(y0) - a*log10(x0);
2033  b = y0/pow(x0,a); // pow(10.,b0);
2034 
2035  a += 1.0;
2036  if( a == 0 ) result += b*log(e0/x0);
2037  else result += y0*(e0*pow(d,a-1) - x0 )/a;
2038  a += 1.0;
2039 
2040  if( a == 0 ) fIntegralCerenkov[0] += b*log(e0/x0);
2041  else fIntegralCerenkov[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
2042 
2043  // G4cout<<"a = "<<a<<"; b0 = "<<b0<<"; b = "
2044  // <<b<<"; result = "<<result<<G4endl;
2045 
2046  return result;
2047 
2048 }
G4double G4ParticleHPJENDLHEData::G4double result
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
tuple b
Definition: test.py:12
double G4double
Definition: G4Types.hh:76
tuple c
Definition: test.py:13
G4double G4PAIxSection::SumOverBorder ( G4int  intervalNumber,
G4double  energy 
)

Definition at line 1862 of file G4PAIxSection.cc.

1864 {
1865  G4double x0,x1,y0,yy1,a,b,/*c,*/d,e0,result;
1866 
1867  e0 = en0;
1868  x0 = fSplineEnergy[i];
1869  x1 = fSplineEnergy[i+1];
1870  y0 = fDifPAIxSection[i];
1871  yy1 = fDifPAIxSection[i+1];
1872 
1873  //c = x1/x0;
1874  d = e0/x0;
1875  a = log10(yy1/y0)/log10(x1/x0);
1876  if(a > 10.0) return 0.;
1877 
1878  if(fVerbose>0) G4cout<<"SumOverBorder, a = "<<a<<G4endl;
1879 
1880  // b0 = log10(y0) - a*log10(x0);
1881  b = y0/pow(x0,a); // pow(10.,b);
1882 
1883  a += 1.;
1884  if( std::fabs(a) < 1.e-6 )
1885  {
1886  result = b*log(x0/e0);
1887  }
1888  else
1889  {
1890  result = y0*(x0 - e0*pow(d,a-1))/a;
1891  }
1892  a += 1.;
1893  if( std::fabs(a) < 1.e-6 )
1894  {
1895  fIntegralPAIxSection[0] += b*log(x0/e0);
1896  }
1897  else
1898  {
1899  fIntegralPAIxSection[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
1900  }
1901  x0 = fSplineEnergy[i - 1];
1902  x1 = fSplineEnergy[i - 2];
1903  y0 = fDifPAIxSection[i - 1];
1904  yy1 = fDifPAIxSection[i - 2];
1905 
1906  //c = x1/x0;
1907  d = e0/x0;
1908  a = log10(yy1/y0)/log10(x1/x0);
1909  // b0 = log10(y0) - a*log10(x0);
1910  b = y0/pow(x0,a);
1911  a += 1.;
1912  if( std::fabs(a) < 1.e-6 )
1913  {
1914  result += b*log(e0/x0);
1915  }
1916  else
1917  {
1918  result += y0*(e0*pow(d,a-1) - x0)/a;
1919  }
1920  a += 1.;
1921  if( std::fabs(a) < 1.e-6 )
1922  {
1923  fIntegralPAIxSection[0] += b*log(e0/x0);
1924  }
1925  else
1926  {
1927  fIntegralPAIxSection[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
1928  }
1929  return result;
1930 
1931 }
G4double G4ParticleHPJENDLHEData::G4double result
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
tuple b
Definition: test.py:12
G4GLOB_DLL std::ostream G4cout
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
G4double G4PAIxSection::SumOverBorderdEdx ( G4int  intervalNumber,
G4double  energy 
)

Definition at line 1935 of file G4PAIxSection.cc.

1937 {
1938  G4double x0,x1,y0,yy1,a,b,/*c,*/d,e0,result;
1939 
1940  e0 = en0;
1941  x0 = fSplineEnergy[i];
1942  x1 = fSplineEnergy[i+1];
1943  y0 = fDifPAIxSection[i];
1944  yy1 = fDifPAIxSection[i+1];
1945 
1946  //c = x1/x0;
1947  d = e0/x0;
1948  a = log10(yy1/y0)/log10(x1/x0);
1949  if(a > 10.0) return 0.;
1950  // b0 = log10(y0) - a*log10(x0);
1951  b = y0/pow(x0,a); // pow(10.,b);
1952 
1953  a += 2;
1954  if(a == 0)
1955  {
1956  result = b*log(x0/e0);
1957  }
1958  else
1959  {
1960  result = y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
1961  }
1962  x0 = fSplineEnergy[i - 1];
1963  x1 = fSplineEnergy[i - 2];
1964  y0 = fDifPAIxSection[i - 1];
1965  yy1 = fDifPAIxSection[i - 2];
1966 
1967  // c = x1/x0;
1968  d = e0/x0;
1969  a = log10(yy1/y0)/log10(x1/x0);
1970  // b0 = log10(y0) - a*log10(x0);
1971  b = y0/pow(x0,a);
1972  a += 2;
1973  if(a == 0)
1974  {
1975  result += b*log(e0/x0);
1976  }
1977  else
1978  {
1979  result += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
1980  }
1981  return result;
1982 
1983 }
G4double G4ParticleHPJENDLHEData::G4double result
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
tuple b
Definition: test.py:12
double G4double
Definition: G4Types.hh:76
G4double G4PAIxSection::SumOverBordMM ( G4int  intervalNumber,
G4double  energy 
)

Definition at line 2055 of file G4PAIxSection.cc.

2057 {
2058  G4double x0,x1,y0,yy1,a,b,e0,c,d,result;
2059 
2060  e0 = en0;
2061  x0 = fSplineEnergy[i];
2062  x1 = fSplineEnergy[i+1];
2063  y0 = fdNdxMM[i];
2064  yy1 = fdNdxMM[i+1];
2065 
2066  // G4cout<<G4endl;
2067  // G4cout<<"SumBordC, i = "<<i<<"; en0 = "<<en0<<"; x0 ="<<x0<<"; x1 = "<<x1
2068  // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
2069  c = x1/x0;
2070  d = e0/x0;
2071  a = log10(yy1/y0)/log10(c);
2072  if(a > 10.0) return 0.;
2073  // b0 = log10(y0) - a*log10(x0);
2074  b = y0/pow(x0,a); // pow(10.,b0);
2075 
2076  a += 1.0;
2077  if( a == 0 ) result = b*log(x0/e0);
2078  else result = y0*(x0 - e0*pow(d,a-1))/a;
2079  a += 1.0;
2080 
2081  if( a == 0 ) fIntegralMM[0] += b*log(x0/e0);
2082  else fIntegralMM[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
2083 
2084 // G4cout<<"a = "<<a<<"; b0 = "<<b0<<"; b = "<<b<<"; result = "<<result<<G4endl;
2085 
2086  x0 = fSplineEnergy[i - 1];
2087  x1 = fSplineEnergy[i - 2];
2088  y0 = fdNdxMM[i - 1];
2089  yy1 = fdNdxMM[i - 2];
2090 
2091  // G4cout<<"x0 ="<<x0<<"; x1 = "<<x1
2092  // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
2093 
2094  c = x1/x0;
2095  d = e0/x0;
2096  a = log10(yy1/y0)/log10(x1/x0);
2097  // b0 = log10(y0) - a*log10(x0);
2098  b = y0/pow(x0,a); // pow(10.,b0);
2099 
2100  a += 1.0;
2101  if( a == 0 ) result += b*log(e0/x0);
2102  else result += y0*(e0*pow(d,a-1) - x0 )/a;
2103  a += 1.0;
2104 
2105  if( a == 0 ) fIntegralMM[0] += b*log(e0/x0);
2106  else fIntegralMM[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
2107 
2108  // G4cout<<"a = "<<a<<"; b0 = "<<b0<<"; b = "
2109  // <<b<<"; result = "<<result<<G4endl;
2110 
2111  return result;
2112 
2113 }
G4double G4ParticleHPJENDLHEData::G4double result
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
tuple b
Definition: test.py:12
double G4double
Definition: G4Types.hh:76
tuple c
Definition: test.py:13
G4double G4PAIxSection::SumOverBordPlasmon ( G4int  intervalNumber,
G4double  energy 
)

Definition at line 2120 of file G4PAIxSection.cc.

2122 {
2123  G4double x0,x1,y0,yy1,a,b,c,d,e0,result;
2124 
2125  e0 = en0;
2126  x0 = fSplineEnergy[i];
2127  x1 = fSplineEnergy[i+1];
2128  y0 = fdNdxPlasmon[i];
2129  yy1 = fdNdxPlasmon[i+1];
2130 
2131  c = x1/x0;
2132  d = e0/x0;
2133  a = log10(yy1/y0)/log10(c);
2134  if(a > 10.0) return 0.;
2135  // b0 = log10(y0) - a*log10(x0);
2136  b = y0/pow(x0,a); //pow(10.,b);
2137 
2138  a += 1.0;
2139  if( a == 0 ) result = b*log(x0/e0);
2140  else result = y0*(x0 - e0*pow(d,a-1))/a;
2141  a += 1.0;
2142 
2143  if( a == 0 ) fIntegralPlasmon[0] += b*log(x0/e0);
2144  else fIntegralPlasmon[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
2145 
2146  x0 = fSplineEnergy[i - 1];
2147  x1 = fSplineEnergy[i - 2];
2148  y0 = fdNdxPlasmon[i - 1];
2149  yy1 = fdNdxPlasmon[i - 2];
2150 
2151  c = x1/x0;
2152  d = e0/x0;
2153  a = log10(yy1/y0)/log10(c);
2154  // b0 = log10(y0) - a*log10(x0);
2155  b = y0/pow(x0,a);// pow(10.,b0);
2156 
2157  a += 1.0;
2158  if( a == 0 ) result += b*log(e0/x0);
2159  else result += y0*(e0*pow(d,a-1) - x0)/a;
2160  a += 1.0;
2161 
2162  if( a == 0 ) fIntegralPlasmon[0] += b*log(e0/x0);
2163  else fIntegralPlasmon[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
2164 
2165  return result;
2166 
2167 }
G4double G4ParticleHPJENDLHEData::G4double result
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
tuple b
Definition: test.py:12
double G4double
Definition: G4Types.hh:76
tuple c
Definition: test.py:13
G4double G4PAIxSection::SumOverBordResonance ( G4int  intervalNumber,
G4double  energy 
)

Definition at line 2174 of file G4PAIxSection.cc.

2176 {
2177  G4double x0,x1,y0,yy1,a,b,c,d,e0,result;
2178 
2179  e0 = en0;
2180  x0 = fSplineEnergy[i];
2181  x1 = fSplineEnergy[i+1];
2182  y0 = fdNdxResonance[i];
2183  yy1 = fdNdxResonance[i+1];
2184 
2185  c = x1/x0;
2186  d = e0/x0;
2187  a = log10(yy1/y0)/log10(c);
2188  if(a > 10.0) return 0.;
2189  // b0 = log10(y0) - a*log10(x0);
2190  b = y0/pow(x0,a); //pow(10.,b);
2191 
2192  a += 1.0;
2193  if( a == 0 ) result = b*log(x0/e0);
2194  else result = y0*(x0 - e0*pow(d,a-1))/a;
2195  a += 1.0;
2196 
2197  if( a == 0 ) fIntegralResonance[0] += b*log(x0/e0);
2198  else fIntegralResonance[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
2199 
2200  x0 = fSplineEnergy[i - 1];
2201  x1 = fSplineEnergy[i - 2];
2202  y0 = fdNdxResonance[i - 1];
2203  yy1 = fdNdxResonance[i - 2];
2204 
2205  c = x1/x0;
2206  d = e0/x0;
2207  a = log10(yy1/y0)/log10(c);
2208  // b0 = log10(y0) - a*log10(x0);
2209  b = y0/pow(x0,a);// pow(10.,b0);
2210 
2211  a += 1.0;
2212  if( a == 0 ) result += b*log(e0/x0);
2213  else result += y0*(e0*pow(d,a-1) - x0)/a;
2214  a += 1.0;
2215 
2216  if( a == 0 ) fIntegralResonance[0] += b*log(e0/x0);
2217  else fIntegralResonance[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
2218 
2219  return result;
2220 
2221 }
G4double G4ParticleHPJENDLHEData::G4double result
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
tuple b
Definition: test.py:12
double G4double
Definition: G4Types.hh:76
tuple c
Definition: test.py:13
G4double G4PAIxSection::SumOverInterCerenkov ( G4int  intervalNumber)

Definition at line 1719 of file G4PAIxSection.cc.

1720 {
1721  G4double x0,x1,y0,yy1,a,b,c,result;
1722 
1723  x0 = fSplineEnergy[i];
1724  x1 = fSplineEnergy[i+1];
1725 
1726  if(x1+x0 <= 0.0 || std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6) return 0.;
1727 
1728  y0 = fdNdxCerenkov[i];
1729  yy1 = fdNdxCerenkov[i+1];
1730  // G4cout<<"SumC, i = "<<i<<"; x0 ="<<x0<<"; x1 = "<<x1
1731  // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
1732 
1733  c = x1/x0;
1734  a = log10(yy1/y0)/log10(c);
1735  b = y0/pow(x0,a);
1736 
1737  a += 1.0;
1738  if(a == 0) result = b*log(c);
1739  else result = y0*(x1*pow(c,a-1) - x0)/a;
1740  a += 1.0;
1741 
1742  if( a == 0 ) fIntegralCerenkov[0] += b*log(x1/x0);
1743  else fIntegralCerenkov[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
1744  // G4cout<<"a = "<<a<<"; b = "<<b<<"; result = "<<result<<G4endl;
1745  return result;
1746 
1747 } // end of SumOverInterCerenkov
G4double G4ParticleHPJENDLHEData::G4double result
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
tuple b
Definition: test.py:12
double G4double
Definition: G4Types.hh:76
tuple c
Definition: test.py:13
G4double G4PAIxSection::SumOverInterMM ( G4int  intervalNumber)

Definition at line 1755 of file G4PAIxSection.cc.

1756 {
1757  G4double x0,x1,y0,yy1,a,b,c,result;
1758 
1759  x0 = fSplineEnergy[i];
1760  x1 = fSplineEnergy[i+1];
1761 
1762  if(x1+x0 <= 0.0 || std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6) return 0.;
1763 
1764  y0 = fdNdxMM[i];
1765  yy1 = fdNdxMM[i+1];
1766  //G4cout<<"SumC, i = "<<i<<"; x0 ="<<x0<<"; x1 = "<<x1
1767  // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
1768 
1769  c = x1/x0;
1770  //G4cout<<" c = "<<c<< " yy1/y0= " << yy1/y0 <<G4endl;
1771  a = log10(yy1/y0)/log10(c);
1772  if(a > 10.0) return 0.;
1773  b = y0/pow(x0,a);
1774 
1775  a += 1.0;
1776  if(a == 0) result = b*log(c);
1777  else result = y0*(x1*pow(c,a-1) - x0)/a;
1778  a += 1.0;
1779 
1780  if( a == 0 ) fIntegralMM[0] += b*log(c);
1781  else fIntegralMM[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
1782  //G4cout<<"a = "<<a<<"; b = "<<b<<"; result = "<<result<<G4endl;
1783  return result;
1784 
1785 } // end of SumOverInterMM
G4double G4ParticleHPJENDLHEData::G4double result
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
tuple b
Definition: test.py:12
double G4double
Definition: G4Types.hh:76
tuple c
Definition: test.py:13
G4double G4PAIxSection::SumOverInterPlasmon ( G4int  intervalNumber)

Definition at line 1793 of file G4PAIxSection.cc.

1794 {
1795  G4double x0,x1,y0,yy1,a,b,c,result;
1796 
1797  x0 = fSplineEnergy[i];
1798  x1 = fSplineEnergy[i+1];
1799 
1800  if(x1+x0 <= 0.0 || std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6) return 0.;
1801 
1802  y0 = fdNdxPlasmon[i];
1803  yy1 = fdNdxPlasmon[i+1];
1804  c =x1/x0;
1805  a = log10(yy1/y0)/log10(c);
1806  if(a > 10.0) return 0.;
1807  // b = log10(y0) - a*log10(x0);
1808  b = y0/pow(x0,a);
1809 
1810  a += 1.0;
1811  if(a == 0) result = b*log(x1/x0);
1812  else result = y0*(x1*pow(c,a-1) - x0)/a;
1813  a += 1.0;
1814 
1815  if( a == 0 ) fIntegralPlasmon[0] += b*log(x1/x0);
1816  else fIntegralPlasmon[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
1817 
1818  return result;
1819 
1820 } // end of SumOverInterPlasmon
G4double G4ParticleHPJENDLHEData::G4double result
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
tuple b
Definition: test.py:12
double G4double
Definition: G4Types.hh:76
tuple c
Definition: test.py:13
G4double G4PAIxSection::SumOverInterResonance ( G4int  intervalNumber)

Definition at line 1828 of file G4PAIxSection.cc.

1829 {
1830  G4double x0,x1,y0,yy1,a,b,c,result;
1831 
1832  x0 = fSplineEnergy[i];
1833  x1 = fSplineEnergy[i+1];
1834 
1835  if(x1+x0 <= 0.0 || std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6) return 0.;
1836 
1837  y0 = fdNdxResonance[i];
1838  yy1 = fdNdxResonance[i+1];
1839  c =x1/x0;
1840  a = log10(yy1/y0)/log10(c);
1841  if(a > 10.0) return 0.;
1842  // b = log10(y0) - a*log10(x0);
1843  b = y0/pow(x0,a);
1844 
1845  a += 1.0;
1846  if(a == 0) result = b*log(x1/x0);
1847  else result = y0*(x1*pow(c,a-1) - x0)/a;
1848  a += 1.0;
1849 
1850  if( a == 0 ) fIntegralResonance[0] += b*log(x1/x0);
1851  else fIntegralResonance[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
1852 
1853  return result;
1854 
1855 } // end of SumOverInterResonance
G4double G4ParticleHPJENDLHEData::G4double result
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
tuple b
Definition: test.py:12
double G4double
Definition: G4Types.hh:76
tuple c
Definition: test.py:13
G4double G4PAIxSection::SumOverInterval ( G4int  intervalNumber)

Definition at line 1638 of file G4PAIxSection.cc.

1639 {
1640  G4double x0,x1,y0,yy1,a,b,c,result;
1641 
1642  x0 = fSplineEnergy[i];
1643  x1 = fSplineEnergy[i+1];
1644  if(fVerbose>0) G4cout<<"SumOverInterval i= " << i << " x0 = "<<x0<<"; x1 = "<<x1<<G4endl;
1645 
1646  if( x1+x0 <= 0.0 || std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6) return 0.;
1647 
1648  y0 = fDifPAIxSection[i];
1649  yy1 = fDifPAIxSection[i+1];
1650 
1651  if(fVerbose>0) G4cout<<"x0 = "<<x0<<"; x1 = "<<x1<<", y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
1652 
1653  c = x1/x0;
1654  a = log10(yy1/y0)/log10(c);
1655 
1656  if(fVerbose>0) G4cout<<"SumOverInterval, a = "<<a<<"; c = "<<c<<G4endl;
1657 
1658  // b = log10(y0) - a*log10(x0);
1659  b = y0/pow(x0,a);
1660  a += 1.;
1661  if( std::fabs(a) < 1.e-6 )
1662  {
1663  result = b*log(x1/x0);
1664  }
1665  else
1666  {
1667  result = y0*(x1*pow(c,a-1) - x0)/a;
1668  }
1669  a += 1.;
1670  if( std::fabs(a) < 1.e-6 )
1671  {
1672  fIntegralPAIxSection[0] += b*log(x1/x0);
1673  }
1674  else
1675  {
1676  fIntegralPAIxSection[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
1677  }
1678  if(fVerbose>0) G4cout<<"SumOverInterval, result = "<<result<<G4endl;
1679  return result;
1680 
1681 } // end of SumOverInterval
G4double G4ParticleHPJENDLHEData::G4double result
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
tuple b
Definition: test.py:12
G4GLOB_DLL std::ostream G4cout
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
tuple c
Definition: test.py:13
G4double G4PAIxSection::SumOverIntervaldEdx ( G4int  intervalNumber)

Definition at line 1685 of file G4PAIxSection.cc.

1686 {
1687  G4double x0,x1,y0,yy1,a,b,c,result;
1688 
1689  x0 = fSplineEnergy[i];
1690  x1 = fSplineEnergy[i+1];
1691 
1692  if(x1+x0 <= 0.0 || std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6) return 0.;
1693 
1694  y0 = fDifPAIxSection[i];
1695  yy1 = fDifPAIxSection[i+1];
1696  c = x1/x0;
1697  a = log10(yy1/y0)/log10(c);
1698  // b = log10(y0) - a*log10(x0);
1699  b = y0/pow(x0,a);
1700  a += 2;
1701  if(a == 0)
1702  {
1703  result = b*log(x1/x0);
1704  }
1705  else
1706  {
1707  result = y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
1708  }
1709  return result;
1710 
1711 } // end of SumOverInterval
G4double G4ParticleHPJENDLHEData::G4double result
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
tuple b
Definition: test.py:12
double G4double
Definition: G4Types.hh:76
tuple c
Definition: test.py:13

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