Geant4.10
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Public Member Functions
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 94 of file G4PAIxSection.cc.

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

Definition at line 131 of file G4PAIxSection.cc.

References G4Material::GetDensity(), G4Material::GetIndex(), and G4MaterialCutsCouple::GetMaterial().

132 {
133  fDensity = matCC->GetMaterial()->GetDensity();
134  G4int matIndex = matCC->GetMaterial()->GetIndex();
135  fMaterialIndex = matIndex;
136  fSandia = new G4SandiaTable(matIndex);
137  fVerbose = 0;
138 
139  G4int i, j;
140  fMatSandiaMatrix = new G4OrderedTable();
141 
142  for (i = 0; i < fSandia->GetMaxInterval()-1; i++)
143  {
144  fMatSandiaMatrix->push_back(new G4DataVector(5,0.));
145  }
146  for (i = 0; i < fSandia->GetMaxInterval()-1; i++)
147  {
148  (*(*fMatSandiaMatrix)[i])[0] = fSandia->GetSandiaMatTable(i,0);
149 
150  for(j = 1; j < 5; j++)
151  {
152  (*(*fMatSandiaMatrix)[i])[j] = fSandia->GetSandiaMatTable(i,j)*fDensity;
153  }
154  }
156  // fEnergyInterval = fA1 = fA2 = fA3 = fA4 = 0;
157 }
size_t GetIndex() const
Definition: G4Material.hh:260
G4double GetDensity() const
Definition: G4Material.hh:178
G4int GetMaxInterval() const
void ComputeLowEnergyCof()
int G4int
Definition: G4Types.hh:78
G4double GetSandiaMatTable(G4int, G4int)
const G4Material * GetMaterial() const
G4PAIxSection::G4PAIxSection ( G4int  materialIndex,
G4double  maxEnergyTransfer 
)

Definition at line 161 of file G4PAIxSection.cc.

References G4Material::GetMaterialTable().

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

Definition at line 273 of file G4PAIxSection.cc.

References G4Material::GetMaterialTable().

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

Definition at line 434 of file G4PAIxSection.cc.

References G4Material::GetMaterialTable(), G4SandiaTable::GetPhotoAbsorpCof(), G4SandiaTable::SandiaIntervals(), and G4SandiaTable::SandiaMixing().

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

Definition at line 585 of file G4PAIxSection.cc.

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

Member Function Documentation

void G4PAIxSection::ComputeLowEnergyCof ( const G4Material material)

Definition at line 746 of file G4PAIxSection.cc.

References G4Material::GetElement(), G4Material::GetNumberOfElements(), and G4Element::GetZ().

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:200
int G4int
Definition: G4Types.hh:78
size_t GetNumberOfElements() const
Definition: G4Material.hh:184
double G4double
Definition: G4Types.hh:76
void G4PAIxSection::ComputeLowEnergyCof ( )

Definition at line 779 of file G4PAIxSection.cc.

References G4Material::GetMaterialTable().

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:564
std::vector< G4Material * > G4MaterialTable
int G4int
Definition: G4Types.hh:78
double G4double
Definition: G4Types.hh:76
G4double G4PAIxSection::DifPAIxSection ( G4int  intervalNumber,
G4double  betaGammaSq 
)

Definition at line 1224 of file G4PAIxSection.cc.

References python.hepunit::electron_mass_c2, python.hepunit::fine_structure_const, python.hepunit::hbarc, and python.hepunit::pi.

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

Definition at line 2319 of file G4PAIxSection.cc.

References G4UniformRand, and position.

2320 {
2321  G4int iTransfer ;
2322 
2323  G4double energyTransfer, position;
2324 
2325  position = fIntegralCerenkov[1]*G4UniformRand();
2326 
2327  for( iTransfer = 1; iTransfer <= fSplineNumber; iTransfer++ )
2328  {
2329  if( position >= fIntegralCerenkov[iTransfer] ) break;
2330  }
2331  if(iTransfer > fSplineNumber) iTransfer--;
2332 
2333  energyTransfer = fSplineEnergy[iTransfer];
2334 
2335  if(iTransfer > 1)
2336  {
2337  energyTransfer -= (fSplineEnergy[iTransfer]-fSplineEnergy[iTransfer-1])*G4UniformRand();
2338  }
2339  return energyTransfer;
2340 }
int G4int
Definition: G4Types.hh:78
#define G4UniformRand()
Definition: Randomize.hh:87
int position
Definition: filter.cc:7
double G4double
Definition: G4Types.hh:76
G4double G4PAIxSection::GetDifPAIxSection ( G4int  i)
inline

Definition at line 193 of file G4PAIxSection.hh.

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

Definition at line 1128 of file G4PAIxSection.cc.

References python.hepunit::cm2, energy(), g(), and python.hepunit::keV.

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

Definition at line 191 of file G4PAIxSection.hh.

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

Definition at line 2240 of file G4PAIxSection.cc.

References G4UniformRand, and position.

2241 {
2242  G4int iTransfer ;
2243 
2244  G4double energyTransfer, position;
2245 
2246  position = fIntegralPAIxSection[1]*G4UniformRand();
2247 
2248  for( iTransfer = 1; iTransfer <= fSplineNumber; iTransfer++ )
2249  {
2250  if( position >= fIntegralPAIxSection[iTransfer] ) break;
2251  }
2252  if(iTransfer > fSplineNumber) iTransfer--;
2253 
2254  energyTransfer = fSplineEnergy[iTransfer];
2255 
2256  if(iTransfer > 1)
2257  {
2258  energyTransfer -= (fSplineEnergy[iTransfer]-fSplineEnergy[iTransfer-1])*G4UniformRand();
2259  }
2260  return energyTransfer;
2261 }
int G4int
Definition: G4Types.hh:78
#define G4UniformRand()
Definition: Randomize.hh:87
int position
Definition: filter.cc:7
double G4double
Definition: G4Types.hh:76
G4double G4PAIxSection::GetIntegralCerenkov ( G4int  i) const
inline

Definition at line 324 of file G4PAIxSection.hh.

Referenced by G4PAIPhotonModel::BuildPAIonisationTable().

325 {
326  if(i < 1 || i > fSplineNumber) { CallError(i, "GetIntegralCerenkov"); }
327  return fIntegralCerenkov[i];
328 }
G4double G4PAIxSection::GetIntegralMM ( G4int  i) const
inline

Definition at line 330 of file G4PAIxSection.hh.

331 {
332  if(i < 1 || i > fSplineNumber) { CallError(i, "GetIntegralMM"); }
333  return fIntegralMM[i];
334 }
G4double G4PAIxSection::GetIntegralPAIdEdx ( G4int  i) const
inline

Definition at line 318 of file G4PAIxSection.hh.

Referenced by G4PAIPhotonModel::BuildPAIonisationTable().

319 {
320  if(i < 1 || i > fSplineNumber) { CallError(i, "GetIntegralPAIdEdx"); }
321  return fIntegralPAIdEdx[i];
322 }
G4double G4PAIxSection::GetIntegralPAIxSection ( G4int  i) const
inline

Definition at line 312 of file G4PAIxSection.hh.

Referenced by G4PAIPhotonModel::BuildPAIonisationTable().

313 {
314  if(i < 1 || i > fSplineNumber) { CallError(i, "GetIntegralPAIxSection"); }
315  return fIntegralPAIxSection[i];
316 }
G4double G4PAIxSection::GetIntegralPlasmon ( G4int  i) const
inline

Definition at line 336 of file G4PAIxSection.hh.

Referenced by G4PAIPhotonModel::BuildPAIonisationTable().

337 {
338  if(i < 1 || i > fSplineNumber) { CallError(i, "GetIntegralPlasmon"); }
339  return fIntegralPlasmon[i];
340 }
G4double G4PAIxSection::GetIntegralResonance ( G4int  i) const
inline

Definition at line 342 of file G4PAIxSection.hh.

343 {
344  if(i < 1 || i > fSplineNumber) { CallError(i, "GetIntegralResonance"); }
345  return fIntegralResonance[i];
346 }
G4int G4PAIxSection::GetIntervalNumber ( ) const
inline

Definition at line 189 of file G4PAIxSection.hh.

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

Definition at line 301 of file G4PAIxSection.hh.

302 {
303  return fLorentzFactor[j];
304 }
G4double G4PAIxSection::GetLowEnergyCof ( ) const
inline

Definition at line 207 of file G4PAIxSection.hh.

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

Definition at line 200 of file G4PAIxSection.hh.

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

Definition at line 199 of file G4PAIxSection.hh.

Referenced by G4PAIPhotonModel::BuildPAIonisationTable().

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

Definition at line 201 of file G4PAIxSection.hh.

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

Definition at line 202 of file G4PAIxSection.hh.

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

Definition at line 203 of file G4PAIxSection.hh.

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

Definition at line 2346 of file G4PAIxSection.cc.

References G4UniformRand, and position.

2347 {
2348  G4int iTransfer ;
2349 
2350  G4double energyTransfer, position;
2351 
2352  position = fIntegralMM[1]*G4UniformRand();
2353 
2354  for( iTransfer = 1; iTransfer <= fSplineNumber; iTransfer++ )
2355  {
2356  if( position >= fIntegralMM[iTransfer] ) break;
2357  }
2358  if(iTransfer > fSplineNumber) iTransfer--;
2359 
2360  energyTransfer = fSplineEnergy[iTransfer];
2361 
2362  if(iTransfer > 1)
2363  {
2364  energyTransfer -= (fSplineEnergy[iTransfer]-fSplineEnergy[iTransfer-1])*G4UniformRand();
2365  }
2366  return energyTransfer;
2367 }
int G4int
Definition: G4Types.hh:78
#define G4UniformRand()
Definition: Randomize.hh:87
int position
Definition: filter.cc:7
double G4double
Definition: G4Types.hh:76
G4double G4PAIxSection::GetNormalizationCof ( ) const
inline

Definition at line 205 of file G4PAIxSection.hh.

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

Definition at line 185 of file G4PAIxSection.hh.

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

Definition at line 194 of file G4PAIxSection.hh.

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

Definition at line 195 of file G4PAIxSection.hh.

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

Definition at line 196 of file G4PAIxSection.hh.

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

Definition at line 197 of file G4PAIxSection.hh.

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

Definition at line 296 of file G4PAIxSection.hh.

297 {
298  return fPAItable[i][j];
299 }
G4double G4PAIxSection::GetPhotonRange ( G4double  energy)

Definition at line 1095 of file G4PAIxSection.cc.

References DBL_MAX, DBL_MIN, and G4InuclParticleNames::lambda.

1096 {
1097  G4int i;
1098  G4double energy2, energy3, energy4, result, lambda;
1099 
1100  energy2 = energy1*energy1;
1101  energy3 = energy2*energy1;
1102  energy4 = energy3*energy1;
1103 
1104  // G4double* SandiaCof = fSandia->GetSandiaCofForMaterialPAI(energy1);
1105  // result = SandiaCof[0]/energy1+SandiaCof[1]/energy2+SandiaCof[2]/energy3+SandiaCof[3]/energy4;
1106  // result *= fDensity;
1107 
1108  for( i = 1; i <= fIntervalNumber; i++ )
1109  {
1110  if( energy1 < fEnergyInterval[i]) break;
1111  }
1112  i--;
1113  if(i == 0) i = 1;
1114 
1115  result = fA1[i]/energy1+fA2[i]/energy2+fA3[i]/energy3+fA4[i]/energy4;
1116 
1117  if( result > DBL_MIN ) lambda = 1./result;
1118  else lambda = DBL_MAX;
1119 
1120  return lambda;
1121 }
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 2399 of file G4PAIxSection.cc.

References G4UniformRand, and position.

2400 {
2401  G4int iTransfer ;
2402 
2403  G4double energyTransfer, position;
2404 
2405  position = fIntegralPlasmon[1]*G4UniformRand();
2406 
2407  for( iTransfer = 1; iTransfer <= fSplineNumber; iTransfer++ )
2408  {
2409  if( position >= fIntegralPlasmon[iTransfer] ) break;
2410  }
2411  if(iTransfer > fSplineNumber) iTransfer--;
2412 
2413  energyTransfer = fSplineEnergy[iTransfer];
2414 
2415  if(iTransfer > 1)
2416  {
2417  energyTransfer -= (fSplineEnergy[iTransfer]-fSplineEnergy[iTransfer-1])*G4UniformRand();
2418  }
2419  return energyTransfer;
2420 }
int G4int
Definition: G4Types.hh:78
#define G4UniformRand()
Definition: Randomize.hh:87
int position
Definition: filter.cc:7
double G4double
Definition: G4Types.hh:76
G4double G4PAIxSection::GetResonanceEnergyTransfer ( )

Definition at line 2453 of file G4PAIxSection.cc.

References G4UniformRand, and position.

2454 {
2455  G4int iTransfer ;
2456 
2457  G4double energyTransfer, position;
2458 
2459  position = fIntegralResonance[1]*G4UniformRand();
2460 
2461  for( iTransfer = 1; iTransfer <= fSplineNumber; iTransfer++ )
2462  {
2463  if( position >= fIntegralResonance[iTransfer] ) break;
2464  }
2465  if(iTransfer > fSplineNumber) iTransfer--;
2466 
2467  energyTransfer = fSplineEnergy[iTransfer];
2468 
2469  if(iTransfer > 1)
2470  {
2471  energyTransfer -= (fSplineEnergy[iTransfer]-fSplineEnergy[iTransfer-1])*G4UniformRand();
2472  }
2473  return energyTransfer;
2474 }
int G4int
Definition: G4Types.hh:78
#define G4UniformRand()
Definition: Randomize.hh:87
int position
Definition: filter.cc:7
double G4double
Definition: G4Types.hh:76
G4double G4PAIxSection::GetRutherfordEnergyTransfer ( )

Definition at line 2481 of file G4PAIxSection.cc.

References G4UniformRand, and position.

2482 {
2483  G4int iTransfer ;
2484 
2485  G4double energyTransfer, position;
2486 
2487  position = (fIntegralPlasmon[1]-fIntegralResonance[1])*G4UniformRand();
2488 
2489  for( iTransfer = 1; iTransfer <= fSplineNumber; iTransfer++ )
2490  {
2491  if( position >= (fIntegralPlasmon[iTransfer]-fIntegralResonance[iTransfer]) ) break;
2492  }
2493  if(iTransfer > fSplineNumber) iTransfer--;
2494 
2495  energyTransfer = fSplineEnergy[iTransfer];
2496 
2497  if(iTransfer > 1)
2498  {
2499  energyTransfer -= (fSplineEnergy[iTransfer]-fSplineEnergy[iTransfer-1])*G4UniformRand();
2500  }
2501  return energyTransfer;
2502 }
int G4int
Definition: G4Types.hh:78
#define G4UniformRand()
Definition: Randomize.hh:87
int position
Definition: filter.cc:7
double G4double
Definition: G4Types.hh:76
G4double G4PAIxSection::GetSplineEnergy ( G4int  i) const
inline

Definition at line 306 of file G4PAIxSection.hh.

Referenced by G4PAIPhotonModel::BuildPAIonisationTable().

307 {
308  if(i < 1 || i > fSplineNumber) { CallError(i, "GetSplineEnergy"); }
309  return fSplineEnergy[i];
310 }
G4int G4PAIxSection::GetSplineSize ( ) const
inline

Definition at line 187 of file G4PAIxSection.hh.

Referenced by G4PAIPhotonModel::BuildPAIonisationTable().

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

Definition at line 2267 of file G4PAIxSection.cc.

References G4Poisson().

2268 {
2269  G4long numOfCollisions;
2270  G4double meanNumber, loss = 0.0;
2271 
2272  // G4cout<<" G4PAIxSection::GetStepCerenkovLoss "<<G4endl;
2273 
2274  meanNumber = fIntegralCerenkov[1]*step;
2275  numOfCollisions = G4Poisson(meanNumber);
2276 
2277  // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl;
2278 
2279  while(numOfCollisions)
2280  {
2281  loss += GetCerenkovEnergyTransfer();
2282  numOfCollisions--;
2283  }
2284  // G4cout<<"PAI Cerenkov loss = "<<loss/keV<<" keV"<<G4endl;
2285 
2286  return loss;
2287 }
G4long G4Poisson(G4double mean)
Definition: G4Poisson.hh:51
long G4long
Definition: G4Types.hh:80
G4double GetCerenkovEnergyTransfer()
double G4double
Definition: G4Types.hh:76
G4double G4PAIxSection::GetStepEnergyLoss ( G4double  step)

Definition at line 2214 of file G4PAIxSection.cc.

References G4Poisson().

2215 {
2216  G4long numOfCollisions;
2217  G4double meanNumber, loss = 0.0;
2218 
2219  // G4cout<<" G4PAIxSection::GetStepEnergyLoss "<<G4endl;
2220 
2221  meanNumber = fIntegralPAIxSection[1]*step;
2222  numOfCollisions = G4Poisson(meanNumber);
2223 
2224  // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl;
2225 
2226  while(numOfCollisions)
2227  {
2228  loss += GetEnergyTransfer();
2229  numOfCollisions--;
2230  }
2231  // G4cout<<"PAI energy loss = "<<loss/keV<<" keV"<<G4endl;
2232 
2233  return loss;
2234 }
G4long G4Poisson(G4double mean)
Definition: G4Poisson.hh:51
long G4long
Definition: G4Types.hh:80
double G4double
Definition: G4Types.hh:76
G4double GetEnergyTransfer()
G4double G4PAIxSection::GetStepMMLoss ( G4double  step)

Definition at line 2293 of file G4PAIxSection.cc.

References G4Poisson().

2294 {
2295  G4long numOfCollisions;
2296  G4double meanNumber, loss = 0.0;
2297 
2298  // G4cout<<" G4PAIxSection::GetStepMMLoss "<<G4endl;
2299 
2300  meanNumber = fIntegralMM[1]*step;
2301  numOfCollisions = G4Poisson(meanNumber);
2302 
2303  // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl;
2304 
2305  while(numOfCollisions)
2306  {
2307  loss += GetMMEnergyTransfer();
2308  numOfCollisions--;
2309  }
2310  // G4cout<<"PAI MM-Cerenkov loss = "<<loss/keV<<" keV"<<G4endl;
2311 
2312  return loss;
2313 }
G4long G4Poisson(G4double mean)
Definition: G4Poisson.hh:51
long G4long
Definition: G4Types.hh:80
G4double GetMMEnergyTransfer()
double G4double
Definition: G4Types.hh:76
G4double G4PAIxSection::GetStepPlasmonLoss ( G4double  step)

Definition at line 2373 of file G4PAIxSection.cc.

References G4Poisson().

2374 {
2375  G4long numOfCollisions;
2376  G4double meanNumber, loss = 0.0;
2377 
2378  // G4cout<<" G4PAIxSection::GetStepPlasmonLoss "<<G4endl;
2379 
2380  meanNumber = fIntegralPlasmon[1]*step;
2381  numOfCollisions = G4Poisson(meanNumber);
2382 
2383  // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl;
2384 
2385  while(numOfCollisions)
2386  {
2387  loss += GetPlasmonEnergyTransfer();
2388  numOfCollisions--;
2389  }
2390  // G4cout<<"PAI Plasmon loss = "<<loss/keV<<" keV"<<G4endl;
2391 
2392  return loss;
2393 }
G4long G4Poisson(G4double mean)
Definition: G4Poisson.hh:51
G4double GetPlasmonEnergyTransfer()
long G4long
Definition: G4Types.hh:80
double G4double
Definition: G4Types.hh:76
G4double G4PAIxSection::GetStepResonanceLoss ( G4double  step)

Definition at line 2426 of file G4PAIxSection.cc.

References G4Poisson().

2427 {
2428  G4long numOfCollisions;
2429  G4double meanNumber, loss = 0.0;
2430 
2431  // G4cout<<" G4PAIxSection::GetStepCreLosnkovs "<<G4endl;
2432 
2433  meanNumber = fIntegralResonance[1]*step;
2434  numOfCollisions = G4Poisson(meanNumber);
2435 
2436  // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl;
2437 
2438  while(numOfCollisions)
2439  {
2440  loss += GetResonanceEnergyTransfer();
2441  numOfCollisions--;
2442  }
2443  // G4cout<<"PAI resonance loss = "<<loss/keV<<" keV"<<G4endl;
2444 
2445  return loss;
2446 }
G4long G4Poisson(G4double mean)
Definition: G4Poisson.hh:51
long G4long
Definition: G4Types.hh:80
G4double GetResonanceEnergyTransfer()
double G4double
Definition: G4Types.hh:76
G4double G4PAIxSection::ImPartDielectricConst ( G4int  intervalNumber,
G4double  energy 
)

Definition at line 1075 of file G4PAIxSection.cc.

References python.hepunit::hbarc.

1077 {
1078  G4double energy2,energy3,energy4,result;
1079 
1080  energy2 = energy1*energy1;
1081  energy3 = energy2*energy1;
1082  energy4 = energy3*energy1;
1083 
1084  result = fA1[k]/energy1+fA2[k]/energy2+fA3[k]/energy3+fA4[k]/energy4;
1085  result *=hbarc/energy1;
1086 
1087  return result;
1088 
1089 } // end of ImPartDielectricConst
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.

References python.hepunit::eV, G4cout, G4endl, G4Material::GetDensity(), G4Material::GetElectronDensity(), G4SandiaTable::GetLowerI1(), G4SandiaTable::GetMaxInterval(), G4SandiaTable::GetSandiaMatTablePAI(), and python.hepunit::keV.

Referenced by G4PAIPhotonModel::BuildPAIonisationTable().

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:178
void NormShift(G4double betaGammaSq)
G4int GetMaxInterval() const
G4double PAIdNdxPlasmon(G4int intervalNumber, G4double betaGammaSq)
G4double GetSandiaMatTablePAI(G4int, G4int)
void ComputeLowEnergyCof()
int G4int
Definition: G4Types.hh:78
G4GLOB_DLL std::ostream G4cout
G4double GetElectronDensity() const
Definition: G4Material.hh:215
G4bool GetLowerI1()
G4double PAIdNdxMM(G4int intervalNumber, G4double betaGammaSq)
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
G4double DifPAIxSection(G4int intervalNumber, G4double betaGammaSq)
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 1521 of file G4PAIxSection.cc.

1522 {
1523  G4int i, k;
1524  fIntegralCerenkov[fSplineNumber] = 0;
1525  fIntegralCerenkov[0] = 0;
1526  k = fIntervalNumber -1;
1527 
1528  for( i = fSplineNumber-1; i >= 1; i-- )
1529  {
1530  if(fSplineEnergy[i] >= fEnergyInterval[k])
1531  {
1532  fIntegralCerenkov[i] = fIntegralCerenkov[i+1] + SumOverInterCerenkov(i);
1533  // G4cout<<"int: i = "<<i<<"; sumC = "<<fIntegralCerenkov[i]<<G4endl;
1534  }
1535  else
1536  {
1537  fIntegralCerenkov[i] = fIntegralCerenkov[i+1] +
1538  SumOverBordCerenkov(i+1,fEnergyInterval[k]);
1539  k--;
1540  // G4cout<<"bord: i = "<<i<<"; sumC = "<<fIntegralCerenkov[i]<<G4endl;
1541  }
1542  }
1543 
1544 } // 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 1552 of file G4PAIxSection.cc.

1553 {
1554  G4int i, k;
1555  fIntegralMM[fSplineNumber] = 0;
1556  fIntegralMM[0] = 0;
1557  k = fIntervalNumber -1;
1558 
1559  for( i = fSplineNumber-1; i >= 1; i-- )
1560  {
1561  if(fSplineEnergy[i] >= fEnergyInterval[k])
1562  {
1563  fIntegralMM[i] = fIntegralMM[i+1] + SumOverInterMM(i);
1564  // G4cout<<"int: i = "<<i<<"; sumC = "<<fIntegralMM[i]<<G4endl;
1565  }
1566  else
1567  {
1568  fIntegralMM[i] = fIntegralMM[i+1] +
1569  SumOverBordMM(i+1,fEnergyInterval[k]);
1570  k--;
1571  // G4cout<<"bord: i = "<<i<<"; sumC = "<<fIntegralMM[i]<<G4endl;
1572  }
1573  }
1574 
1575 } // 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 1489 of file G4PAIxSection.cc.

References G4cout, and G4endl.

1490 {
1491  fIntegralPAIxSection[fSplineNumber] = 0;
1492  fIntegralPAIdEdx[fSplineNumber] = 0;
1493  fIntegralPAIxSection[0] = 0;
1494  G4int i, k = fIntervalNumber -1;
1495 
1496  for( i = fSplineNumber-1; i >= 1; i--)
1497  {
1498  if(fSplineEnergy[i] >= fEnergyInterval[k])
1499  {
1500  fIntegralPAIxSection[i] = fIntegralPAIxSection[i+1] + SumOverInterval(i);
1501  fIntegralPAIdEdx[i] = fIntegralPAIdEdx[i+1] + SumOverIntervaldEdx(i);
1502  }
1503  else
1504  {
1505  fIntegralPAIxSection[i] = fIntegralPAIxSection[i+1] +
1506  SumOverBorder(i+1,fEnergyInterval[k]);
1507  fIntegralPAIdEdx[i] = fIntegralPAIdEdx[i+1] +
1508  SumOverBorderdEdx(i+1,fEnergyInterval[k]);
1509  k--;
1510  }
1511  if(fVerbose>0) G4cout<<"i = "<<i<<"; k = "<<k<<"; intPAIxsc[i] = "<<fIntegralPAIxSection[i]<<G4endl;
1512  }
1513 } // 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 1583 of file G4PAIxSection.cc.

1584 {
1585  fIntegralPlasmon[fSplineNumber] = 0;
1586  fIntegralPlasmon[0] = 0;
1587  G4int k = fIntervalNumber -1;
1588  for(G4int i=fSplineNumber-1;i>=1;i--)
1589  {
1590  if(fSplineEnergy[i] >= fEnergyInterval[k])
1591  {
1592  fIntegralPlasmon[i] = fIntegralPlasmon[i+1] + SumOverInterPlasmon(i);
1593  }
1594  else
1595  {
1596  fIntegralPlasmon[i] = fIntegralPlasmon[i+1] +
1597  SumOverBordPlasmon(i+1,fEnergyInterval[k]);
1598  k--;
1599  }
1600  }
1601 
1602 } // 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 1610 of file G4PAIxSection.cc.

1611 {
1612  fIntegralResonance[fSplineNumber] = 0;
1613  fIntegralResonance[0] = 0;
1614  G4int k = fIntervalNumber -1;
1615  for(G4int i=fSplineNumber-1;i>=1;i--)
1616  {
1617  if(fSplineEnergy[i] >= fEnergyInterval[k])
1618  {
1619  fIntegralResonance[i] = fIntegralResonance[i+1] + SumOverInterResonance(i);
1620  }
1621  else
1622  {
1623  fIntegralResonance[i] = fIntegralResonance[i+1] +
1624  SumOverBordResonance(i+1,fEnergyInterval[k]);
1625  k--;
1626  }
1627  }
1628 
1629 } // 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.

References python.hepunit::electron_mass_c2, python.hepunit::fine_structure_const, G4cout, G4endl, python.hepunit::hbarc, python.hepunit::keV, python.hepunit::pi, and test::x.

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)
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
double G4double
Definition: G4Types.hh:76
G4double DifPAIxSection(G4int intervalNumber, G4double betaGammaSq)
G4double RutherfordIntegral(G4int intervalNumber, G4double limitLow, G4double limitHigh)
G4double G4PAIxSection::PAIdNdxCerenkov ( G4int  intervalNumber,
G4double  betaGammaSq 
)

Definition at line 1301 of file G4PAIxSection.cc.

References python.hepunit::fine_structure_const, python.hepunit::hbarc, and python.hepunit::pi.

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

Definition at line 1359 of file G4PAIxSection.cc.

References python.hepunit::fine_structure_const, python.hepunit::hbarc, and python.hepunit::pi.

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

Definition at line 1409 of file G4PAIxSection.cc.

References python.hepunit::electron_mass_c2, python.hepunit::fine_structure_const, python.hepunit::hbarc, and python.hepunit::pi.

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

Definition at line 1449 of file G4PAIxSection.cc.

References python.hepunit::electron_mass_c2, python.hepunit::fine_structure_const, python.hepunit::hbarc, and python.hepunit::pi.

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

Definition at line 1170 of file G4PAIxSection.cc.

References plottest35::c1, python.hepunit::hbarc, and python.hepunit::pi.

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

Definition at line 1054 of file G4PAIxSection.cc.

References plottest35::c1.

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

Definition at line 209 of file G4PAIxSection.hh.

References test::v.

209 {fVerbose=v;};
void G4PAIxSection::SplainPAI ( G4double  betaGammaSq)

Definition at line 950 of file G4PAIxSection.cc.

References test::a, test::b, G4cout, G4endl, and test::x.

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  } // close 'while'
1045 
1046 } // end of SplainPAI
G4double PAIdNdxResonance(G4int intervalNumber, G4double betaGammaSq)
G4double PAIdNdxCerenkov(G4int intervalNumber, G4double betaGammaSq)
G4double PAIdNdxPlasmon(G4int intervalNumber, G4double betaGammaSq)
G4double RePartDielectricConst(G4double energy)
int G4int
Definition: G4Types.hh:78
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 1981 of file G4PAIxSection.cc.

References test::a, test::b, and test::c.

1983 {
1984  G4double x0,x1,y0,yy1,a,b,e0,c,d,result;
1985 
1986  e0 = en0;
1987  x0 = fSplineEnergy[i];
1988  x1 = fSplineEnergy[i+1];
1989  y0 = fdNdxCerenkov[i];
1990  yy1 = fdNdxCerenkov[i+1];
1991 
1992  // G4cout<<G4endl;
1993  // G4cout<<"SumBordC, i = "<<i<<"; en0 = "<<en0<<"; x0 ="<<x0<<"; x1 = "<<x1
1994  // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
1995  c = x1/x0;
1996  d = e0/x0;
1997  a = log10(yy1/y0)/log10(c);
1998  // b0 = log10(y0) - a*log10(x0);
1999  b = y0/pow(x0,a); // pow(10.,b0);
2000 
2001  a += 1.0;
2002  if( a == 0 ) result = b*log(x0/e0);
2003  else result = y0*(x0 - e0*pow(d,a-1))/a;
2004  a += 1.0;
2005 
2006  if( a == 0 ) fIntegralCerenkov[0] += b*log(x0/e0);
2007  else fIntegralCerenkov[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
2008 
2009 // G4cout<<"a = "<<a<<"; b0 = "<<b0<<"; b = "<<b<<"; result = "<<result<<G4endl;
2010 
2011  x0 = fSplineEnergy[i - 1];
2012  x1 = fSplineEnergy[i - 2];
2013  y0 = fdNdxCerenkov[i - 1];
2014  yy1 = fdNdxCerenkov[i - 2];
2015 
2016  // G4cout<<"x0 ="<<x0<<"; x1 = "<<x1
2017  // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
2018 
2019  c = x1/x0;
2020  d = e0/x0;
2021  a = log10(yy1/y0)/log10(x1/x0);
2022  // b0 = log10(y0) - a*log10(x0);
2023  b = y0/pow(x0,a); // pow(10.,b0);
2024 
2025  a += 1.0;
2026  if( a == 0 ) result += b*log(e0/x0);
2027  else result += y0*(e0*pow(d,a-1) - x0 )/a;
2028  a += 1.0;
2029 
2030  if( a == 0 ) fIntegralCerenkov[0] += b*log(e0/x0);
2031  else fIntegralCerenkov[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
2032 
2033  // G4cout<<"a = "<<a<<"; b0 = "<<b0<<"; b = "
2034  // <<b<<"; result = "<<result<<G4endl;
2035 
2036  return result;
2037 
2038 }
double G4double
Definition: G4Types.hh:76
G4double G4PAIxSection::SumOverBorder ( G4int  intervalNumber,
G4double  energy 
)

Definition at line 1855 of file G4PAIxSection.cc.

References test::a, test::b, G4cout, and G4endl.

1857 {
1858  G4double x0,x1,y0,yy1,a,b,/*c,*/d,e0,result;
1859 
1860  e0 = en0;
1861  x0 = fSplineEnergy[i];
1862  x1 = fSplineEnergy[i+1];
1863  y0 = fDifPAIxSection[i];
1864  yy1 = fDifPAIxSection[i+1];
1865 
1866  //c = x1/x0;
1867  d = e0/x0;
1868  a = log10(yy1/y0)/log10(x1/x0);
1869 
1870  if(fVerbose>0) G4cout<<"SumOverBorder, a = "<<a<<G4endl;
1871 
1872  // b0 = log10(y0) - a*log10(x0);
1873  b = y0/pow(x0,a); // pow(10.,b);
1874 
1875  a += 1.;
1876  if( std::fabs(a) < 1.e-6 )
1877  {
1878  result = b*log(x0/e0);
1879  }
1880  else
1881  {
1882  result = y0*(x0 - e0*pow(d,a-1))/a;
1883  }
1884  a += 1.;
1885  if( std::fabs(a) < 1.e-6 )
1886  {
1887  fIntegralPAIxSection[0] += b*log(x0/e0);
1888  }
1889  else
1890  {
1891  fIntegralPAIxSection[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
1892  }
1893  x0 = fSplineEnergy[i - 1];
1894  x1 = fSplineEnergy[i - 2];
1895  y0 = fDifPAIxSection[i - 1];
1896  yy1 = fDifPAIxSection[i - 2];
1897 
1898  //c = x1/x0;
1899  d = e0/x0;
1900  a = log10(yy1/y0)/log10(x1/x0);
1901  // b0 = log10(y0) - a*log10(x0);
1902  b = y0/pow(x0,a);
1903  a += 1.;
1904  if( std::fabs(a) < 1.e-6 )
1905  {
1906  result += b*log(e0/x0);
1907  }
1908  else
1909  {
1910  result += y0*(e0*pow(d,a-1) - x0)/a;
1911  }
1912  a += 1.;
1913  if( std::fabs(a) < 1.e-6 )
1914  {
1915  fIntegralPAIxSection[0] += b*log(e0/x0);
1916  }
1917  else
1918  {
1919  fIntegralPAIxSection[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
1920  }
1921  return result;
1922 
1923 }
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 1927 of file G4PAIxSection.cc.

References test::a, and test::b.

1929 {
1930  G4double x0,x1,y0,yy1,a,b,/*c,*/d,e0,result;
1931 
1932  e0 = en0;
1933  x0 = fSplineEnergy[i];
1934  x1 = fSplineEnergy[i+1];
1935  y0 = fDifPAIxSection[i];
1936  yy1 = fDifPAIxSection[i+1];
1937 
1938  //c = x1/x0;
1939  d = e0/x0;
1940  a = log10(yy1/y0)/log10(x1/x0);
1941  // b0 = log10(y0) - a*log10(x0);
1942  b = y0/pow(x0,a); // pow(10.,b);
1943 
1944  a += 2;
1945  if(a == 0)
1946  {
1947  result = b*log(x0/e0);
1948  }
1949  else
1950  {
1951  result = y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
1952  }
1953  x0 = fSplineEnergy[i - 1];
1954  x1 = fSplineEnergy[i - 2];
1955  y0 = fDifPAIxSection[i - 1];
1956  yy1 = fDifPAIxSection[i - 2];
1957 
1958  // c = x1/x0;
1959  d = e0/x0;
1960  a = log10(yy1/y0)/log10(x1/x0);
1961  // b0 = log10(y0) - a*log10(x0);
1962  b = y0/pow(x0,a);
1963  a += 2;
1964  if(a == 0)
1965  {
1966  result += b*log(e0/x0);
1967  }
1968  else
1969  {
1970  result += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
1971  }
1972  return result;
1973 
1974 }
double G4double
Definition: G4Types.hh:76
G4double G4PAIxSection::SumOverBordMM ( G4int  intervalNumber,
G4double  energy 
)

Definition at line 2045 of file G4PAIxSection.cc.

References test::a, test::b, and test::c.

2047 {
2048  G4double x0,x1,y0,yy1,a,b,e0,c,d,result;
2049 
2050  e0 = en0;
2051  x0 = fSplineEnergy[i];
2052  x1 = fSplineEnergy[i+1];
2053  y0 = fdNdxMM[i];
2054  yy1 = fdNdxMM[i+1];
2055 
2056  // G4cout<<G4endl;
2057  // G4cout<<"SumBordC, i = "<<i<<"; en0 = "<<en0<<"; x0 ="<<x0<<"; x1 = "<<x1
2058  // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
2059  c = x1/x0;
2060  d = e0/x0;
2061  a = log10(yy1/y0)/log10(c);
2062  // b0 = log10(y0) - a*log10(x0);
2063  b = y0/pow(x0,a); // pow(10.,b0);
2064 
2065  a += 1.0;
2066  if( a == 0 ) result = b*log(x0/e0);
2067  else result = y0*(x0 - e0*pow(d,a-1))/a;
2068  a += 1.0;
2069 
2070  if( a == 0 ) fIntegralMM[0] += b*log(x0/e0);
2071  else fIntegralMM[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
2072 
2073 // G4cout<<"a = "<<a<<"; b0 = "<<b0<<"; b = "<<b<<"; result = "<<result<<G4endl;
2074 
2075  x0 = fSplineEnergy[i - 1];
2076  x1 = fSplineEnergy[i - 2];
2077  y0 = fdNdxMM[i - 1];
2078  yy1 = fdNdxMM[i - 2];
2079 
2080  // G4cout<<"x0 ="<<x0<<"; x1 = "<<x1
2081  // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
2082 
2083  c = x1/x0;
2084  d = e0/x0;
2085  a = log10(yy1/y0)/log10(x1/x0);
2086  // b0 = log10(y0) - a*log10(x0);
2087  b = y0/pow(x0,a); // pow(10.,b0);
2088 
2089  a += 1.0;
2090  if( a == 0 ) result += b*log(e0/x0);
2091  else result += y0*(e0*pow(d,a-1) - x0 )/a;
2092  a += 1.0;
2093 
2094  if( a == 0 ) fIntegralMM[0] += b*log(e0/x0);
2095  else fIntegralMM[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
2096 
2097  // G4cout<<"a = "<<a<<"; b0 = "<<b0<<"; b = "
2098  // <<b<<"; result = "<<result<<G4endl;
2099 
2100  return result;
2101 
2102 }
double G4double
Definition: G4Types.hh:76
G4double G4PAIxSection::SumOverBordPlasmon ( G4int  intervalNumber,
G4double  energy 
)

Definition at line 2109 of file G4PAIxSection.cc.

References test::a, test::b, and test::c.

2111 {
2112  G4double x0,x1,y0,yy1,a,b,c,d,e0,result;
2113 
2114  e0 = en0;
2115  x0 = fSplineEnergy[i];
2116  x1 = fSplineEnergy[i+1];
2117  y0 = fdNdxPlasmon[i];
2118  yy1 = fdNdxPlasmon[i+1];
2119 
2120  c = x1/x0;
2121  d = e0/x0;
2122  a = log10(yy1/y0)/log10(c);
2123  // b0 = log10(y0) - a*log10(x0);
2124  b = y0/pow(x0,a); //pow(10.,b);
2125 
2126  a += 1.0;
2127  if( a == 0 ) result = b*log(x0/e0);
2128  else result = y0*(x0 - e0*pow(d,a-1))/a;
2129  a += 1.0;
2130 
2131  if( a == 0 ) fIntegralPlasmon[0] += b*log(x0/e0);
2132  else fIntegralPlasmon[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
2133 
2134  x0 = fSplineEnergy[i - 1];
2135  x1 = fSplineEnergy[i - 2];
2136  y0 = fdNdxPlasmon[i - 1];
2137  yy1 = fdNdxPlasmon[i - 2];
2138 
2139  c = x1/x0;
2140  d = e0/x0;
2141  a = log10(yy1/y0)/log10(c);
2142  // b0 = log10(y0) - a*log10(x0);
2143  b = y0/pow(x0,a);// pow(10.,b0);
2144 
2145  a += 1.0;
2146  if( a == 0 ) result += b*log(e0/x0);
2147  else result += y0*(e0*pow(d,a-1) - x0)/a;
2148  a += 1.0;
2149 
2150  if( a == 0 ) fIntegralPlasmon[0] += b*log(e0/x0);
2151  else fIntegralPlasmon[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
2152 
2153  return result;
2154 
2155 }
double G4double
Definition: G4Types.hh:76
G4double G4PAIxSection::SumOverBordResonance ( G4int  intervalNumber,
G4double  energy 
)

Definition at line 2162 of file G4PAIxSection.cc.

References test::a, test::b, and test::c.

2164 {
2165  G4double x0,x1,y0,yy1,a,b,c,d,e0,result;
2166 
2167  e0 = en0;
2168  x0 = fSplineEnergy[i];
2169  x1 = fSplineEnergy[i+1];
2170  y0 = fdNdxResonance[i];
2171  yy1 = fdNdxResonance[i+1];
2172 
2173  c = x1/x0;
2174  d = e0/x0;
2175  a = log10(yy1/y0)/log10(c);
2176  // b0 = log10(y0) - a*log10(x0);
2177  b = y0/pow(x0,a); //pow(10.,b);
2178 
2179  a += 1.0;
2180  if( a == 0 ) result = b*log(x0/e0);
2181  else result = y0*(x0 - e0*pow(d,a-1))/a;
2182  a += 1.0;
2183 
2184  if( a == 0 ) fIntegralResonance[0] += b*log(x0/e0);
2185  else fIntegralResonance[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a;
2186 
2187  x0 = fSplineEnergy[i - 1];
2188  x1 = fSplineEnergy[i - 2];
2189  y0 = fdNdxResonance[i - 1];
2190  yy1 = fdNdxResonance[i - 2];
2191 
2192  c = x1/x0;
2193  d = e0/x0;
2194  a = log10(yy1/y0)/log10(c);
2195  // b0 = log10(y0) - a*log10(x0);
2196  b = y0/pow(x0,a);// pow(10.,b0);
2197 
2198  a += 1.0;
2199  if( a == 0 ) result += b*log(e0/x0);
2200  else result += y0*(e0*pow(d,a-1) - x0)/a;
2201  a += 1.0;
2202 
2203  if( a == 0 ) fIntegralResonance[0] += b*log(e0/x0);
2204  else fIntegralResonance[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a;
2205 
2206  return result;
2207 
2208 }
double G4double
Definition: G4Types.hh:76
G4double G4PAIxSection::SumOverInterCerenkov ( G4int  intervalNumber)

Definition at line 1716 of file G4PAIxSection.cc.

References test::a, test::b, and test::c.

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

Definition at line 1752 of file G4PAIxSection.cc.

References test::a, test::b, and test::c.

1753 {
1754  G4double x0,x1,y0,yy1,a,b,c,result;
1755 
1756  x0 = fSplineEnergy[i];
1757  x1 = fSplineEnergy[i+1];
1758 
1759  if( std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6) return 0.;
1760 
1761  y0 = fdNdxMM[i];
1762  yy1 = fdNdxMM[i+1];
1763  // G4cout<<"SumC, i = "<<i<<"; x0 ="<<x0<<"; x1 = "<<x1
1764  // <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
1765 
1766  c = x1/x0;
1767  a = log10(yy1/y0)/log10(c);
1768  b = y0/pow(x0,a);
1769 
1770  a += 1.0;
1771  if(a == 0) result = b*log(c);
1772  else result = y0*(x1*pow(c,a-1) - x0)/a;
1773  a += 1.0;
1774 
1775  if( a == 0 ) fIntegralMM[0] += b*log(x1/x0);
1776  else fIntegralMM[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
1777  // G4cout<<"a = "<<a<<"; b = "<<b<<"; result = "<<result<<G4endl;
1778  return result;
1779 
1780 } // end of SumOverInterMM
double G4double
Definition: G4Types.hh:76
G4double G4PAIxSection::SumOverInterPlasmon ( G4int  intervalNumber)

Definition at line 1788 of file G4PAIxSection.cc.

References test::a, test::b, and test::c.

1789 {
1790  G4double x0,x1,y0,yy1,a,b,c,result;
1791 
1792  x0 = fSplineEnergy[i];
1793  x1 = fSplineEnergy[i+1];
1794 
1795  if( std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6) return 0.;
1796 
1797  y0 = fdNdxPlasmon[i];
1798  yy1 = fdNdxPlasmon[i+1];
1799  c =x1/x0;
1800  a = log10(yy1/y0)/log10(c);
1801  // b = log10(y0) - a*log10(x0);
1802  b = y0/pow(x0,a);
1803 
1804  a += 1.0;
1805  if(a == 0) result = b*log(x1/x0);
1806  else result = y0*(x1*pow(c,a-1) - x0)/a;
1807  a += 1.0;
1808 
1809  if( a == 0 ) fIntegralPlasmon[0] += b*log(x1/x0);
1810  else fIntegralPlasmon[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
1811 
1812  return result;
1813 
1814 } // end of SumOverInterPlasmon
double G4double
Definition: G4Types.hh:76
G4double G4PAIxSection::SumOverInterResonance ( G4int  intervalNumber)

Definition at line 1822 of file G4PAIxSection.cc.

References test::a, test::b, and test::c.

1823 {
1824  G4double x0,x1,y0,yy1,a,b,c,result;
1825 
1826  x0 = fSplineEnergy[i];
1827  x1 = fSplineEnergy[i+1];
1828 
1829  if( std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6) return 0.;
1830 
1831  y0 = fdNdxResonance[i];
1832  yy1 = fdNdxResonance[i+1];
1833  c =x1/x0;
1834  a = log10(yy1/y0)/log10(c);
1835  // b = log10(y0) - a*log10(x0);
1836  b = y0/pow(x0,a);
1837 
1838  a += 1.0;
1839  if(a == 0) result = b*log(x1/x0);
1840  else result = y0*(x1*pow(c,a-1) - x0)/a;
1841  a += 1.0;
1842 
1843  if( a == 0 ) fIntegralResonance[0] += b*log(x1/x0);
1844  else fIntegralResonance[0] += y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
1845 
1846  return result;
1847 
1848 } // end of SumOverInterResonance
double G4double
Definition: G4Types.hh:76
G4double G4PAIxSection::SumOverInterval ( G4int  intervalNumber)

Definition at line 1637 of file G4PAIxSection.cc.

References test::a, test::b, test::c, G4cout, and G4endl.

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

Definition at line 1682 of file G4PAIxSection.cc.

References test::a, test::b, and test::c.

1683 {
1684  G4double x0,x1,y0,yy1,a,b,c,result;
1685 
1686  x0 = fSplineEnergy[i];
1687  x1 = fSplineEnergy[i+1];
1688 
1689  if( std::abs( 2.*(x1-x0)/(x1+x0) ) < 1.e-6) return 0.;
1690 
1691  y0 = fDifPAIxSection[i];
1692  yy1 = fDifPAIxSection[i+1];
1693  c = x1/x0;
1694  a = log10(yy1/y0)/log10(c);
1695  // b = log10(y0) - a*log10(x0);
1696  b = y0/pow(x0,a);
1697  a += 2;
1698  if(a == 0)
1699  {
1700  result = b*log(x1/x0);
1701  }
1702  else
1703  {
1704  result = y0*(x1*x1*pow(c,a-2) - x0*x0)/a;
1705  }
1706  return result;
1707 
1708 } // end of SumOverInterval
double G4double
Definition: G4Types.hh:76

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