Geant4.10
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4PAIModelData.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // $Id: G4PAIModelData.cc 72008 2013-07-03 08:46:39Z vnivanch $
27 //
28 // -------------------------------------------------------------------
29 //
30 // GEANT4 Class
31 // File name: G4PAIModelData.cc
32 //
33 // Author: V. Ivanchenko based on V.Grichine code of G4PAIModel
34 //
35 // Creation date: 16.08.2013
36 //
37 // Modifications:
38 //
39 
40 #include "G4PAIModelData.hh"
41 #include "G4PAIModel.hh"
42 
43 #include "G4SystemOfUnits.hh"
44 #include "G4PhysicalConstants.hh"
45 
46 #include "G4PhysicsLogVector.hh"
47 #include "G4PhysicsFreeVector.hh"
48 #include "G4PhysicsTable.hh"
49 #include "G4MaterialCutsCouple.hh"
50 #include "G4SandiaTable.hh"
51 #include "Randomize.hh"
52 #include "G4Poisson.hh"
53 
54 ////////////////////////////////////////////////////////////////////////
55 
56 using namespace std;
57 
59 {
60  const G4int nPerDecade = 10;
61  const G4double lowestTkin = 50*keV;
62  const G4double highestTkin = 10*TeV;
63 
64  fPAIySection.SetVerbose(ver);
65 
66  fLowestKineticEnergy = std::max(tmin, lowestTkin);
67  fHighestKineticEnergy = tmax;
68  if(tmax < 10*fLowestKineticEnergy) {
69  fHighestKineticEnergy = 10*fLowestKineticEnergy;
70  } else if(tmax > highestTkin) {
71  fHighestKineticEnergy = std::max(highestTkin, 10*fLowestKineticEnergy);
72  }
73  fTotBin = (G4int)(nPerDecade*
74  std::log10(fHighestKineticEnergy/fLowestKineticEnergy));
75 
76  fParticleEnergyVector = new G4PhysicsLogVector(fLowestKineticEnergy,
77  fHighestKineticEnergy,
78  fTotBin);
79  if(0 < ver) {
80  G4cout << "### G4PAIModelData: Nbins= " << fTotBin
81  << " Tmin(MeV)= " << fLowestKineticEnergy/MeV
82  << " Tmax(GeV)= " << fHighestKineticEnergy/GeV
83  << " tmin(keV)= " << tmin/keV << G4endl;
84  }
85 }
86 
87 ////////////////////////////////////////////////////////////////////////////
88 
90 {
91  size_t n = fPAIxscBank.size();
92  if(0 < n) {
93  for(size_t i=0; i<n; ++i) {
94  if(fPAIxscBank[i]) {
95  delete fPAIxscBank[i];
96  }
97  if(fPAIdEdxBank[i]) {
98  delete fPAIdEdxBank[i];
99  }
100  }
101  }
102 }
103 
104 ///////////////////////////////////////////////////////////////////////////////
105 
107  G4double cut, G4PAIModel* model)
108 {
109  const G4Material* mat = couple->GetMaterial();
110  fSandia.Initialize(const_cast<G4Material*>(mat));
111 
112  G4PhysicsTable* PAItransferTable = new G4PhysicsTable(fTotBin+1);
113  G4PhysicsTable* PAIdEdxTable = new G4PhysicsTable(fTotBin+1);
114 
115  G4PhysicsLogVector* dEdxCutVector =
116  new G4PhysicsLogVector(fLowestKineticEnergy,
117  fHighestKineticEnergy,
118  fTotBin);
119 
120  G4PhysicsLogVector* dEdxMeanVector =
121  new G4PhysicsLogVector(fLowestKineticEnergy,
122  fHighestKineticEnergy,
123  fTotBin);
124 
125  G4PhysicsLogVector* dNdxCutVector =
126  new G4PhysicsLogVector(fLowestKineticEnergy,
127  fHighestKineticEnergy,
128  fTotBin);
129 
130  // low energy Sandia interval
131  G4double Tmin = fSandia.GetSandiaMatTablePAI(0,0);
132 
133  // energy safety
134  const G4double deltaLow = 100.*eV;
135 
136  for (G4int i = 0; i <= fTotBin; ++i) {
137 
138  G4double kinEnergy = fParticleEnergyVector->Energy(i);
139  G4double Tmax = model->ComputeMaxEnergy(kinEnergy);
140  G4double tau = kinEnergy/proton_mass_c2;
141  G4double bg2 = tau*( tau + 2. );
142 
143  if (Tmax < Tmin + deltaLow ) { Tmax = Tmin + deltaLow; }
144 
145  fPAIySection.Initialize(mat, Tmax, bg2, &fSandia);
146 
147  //G4cout << i << ". TransferMax(keV)= "<< Tmax/keV << " cut(keV)= "
148  // << cut/keV << " E(MeV)= " << kinEnergy/MeV << G4endl;
149 
150  G4int n = fPAIySection.GetSplineSize();
151  G4PhysicsFreeVector* transferVector = new G4PhysicsFreeVector(n);
152  G4PhysicsFreeVector* dEdxVector = new G4PhysicsFreeVector(n);
153 
154  for( G4int k = 0; k < n; k++ )
155  {
156  G4double t = fPAIySection.GetSplineEnergy(k+1);
157  transferVector->PutValue(k , t,
158  t*fPAIySection.GetIntegralPAIySection(k+1));
159  dEdxVector->PutValue(k, t, fPAIySection.GetIntegralPAIdEdx(k+1));
160  }
161  // G4cout << *transferVector << G4endl;
162 
163  G4double ionloss = fPAIySection.GetMeanEnergyLoss();// total <dE/dx>
164 
165  if(ionloss < 0.0) ionloss = 0.0;
166 
167  dEdxMeanVector->PutValue(i,ionloss);
168 
169  G4double dNdxCut = transferVector->Value(cut)/cut;
170  G4double dEdxCut = dEdxVector->Value(cut)/cut;
171  //G4cout << "i= " << i << " x= " << dNdxCut << G4endl;
172  if(dNdxCut < 0.0) { dNdxCut = 0.0; }
173  dNdxCutVector->PutValue(i, dNdxCut);
174  dEdxCutVector->PutValue(i, dEdxCut);
175 
176  PAItransferTable->insertAt(i,transferVector);
177  PAIdEdxTable->insertAt(i,dEdxVector);
178 
179  } // end of Tkin loop
180 
181  fPAIxscBank.push_back(PAItransferTable);
182  fPAIdEdxBank.push_back(PAIdEdxTable);
183  fdEdxTable.push_back(dEdxMeanVector);
184 
185  fdNdxCutTable.push_back(dNdxCutVector);
186  fdEdxCutTable.push_back(dEdxCutVector);
187 }
188 
189 //////////////////////////////////////////////////////////////////////////////
190 
192  G4double cut) const
193 {
194  // VI: iPlace is the low edge index of the bin
195  // iPlace is in interval from 0 to (N-1)
196  size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
197  size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
198 
199  G4bool one = true;
200  if(scaledTkin >= fParticleEnergyVector->Energy(nPlace)) { iPlace = nPlace; }
201  else if(scaledTkin > fParticleEnergyVector->Energy(0)) {
202  one = false;
203  }
204 
205  // VI: apply interpolation of the vector
206  G4double dEdx = fdEdxTable[coupleIndex]->Value(scaledTkin);
207  G4double del = (*(fPAIdEdxBank[coupleIndex]))(iPlace)->Value(cut);
208  if(!one) {
209  G4double del2 = (*(fPAIdEdxBank[coupleIndex]))(iPlace+1)->Value(cut);
210  G4double E1 = fParticleEnergyVector->Energy(iPlace);
211  G4double E2 = fParticleEnergyVector->Energy(iPlace+1);
212  G4double W = 1.0/(E2 - E1);
213  G4double W1 = (E2 - scaledTkin)*W;
214  G4double W2 = (scaledTkin - E1)*W;
215  del *= W1;
216  del += W2*del2;
217  }
218  dEdx -= del;
219 
220  if( dEdx < 0.) { dEdx = 0.; }
221  return dEdx;
222 }
223 
224 /////////////////////////////////////////////////////////////////////////
225 
226 G4double
228  G4double scaledTkin,
229  G4double tcut, G4double tmax) const
230 {
231  G4double cross, cross1, cross2;
232 
233  // iPlace is in interval from 0 to (N-1)
234  size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
235  size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
236 
237  G4bool one = true;
238  if(scaledTkin >= fParticleEnergyVector->Energy(nPlace)) { iPlace = nPlace; }
239  else if(scaledTkin > fParticleEnergyVector->Energy(0)) {
240  one = false;
241  }
242  G4PhysicsTable* table = fPAIxscBank[coupleIndex];
243 
244  //G4cout<<"iPlace = "<<iPlace<<"; tmax = "
245  // <<tmax<<"; cutEnergy = "<<cutEnergy<<G4endl;
246  cross1 = (*table)(iPlace)->Value(tmax)/tmax;
247  //G4cout<<"cross1 = "<<cross1<<G4endl;
248  cross2 = (*table)(iPlace)->Value(tcut)/tcut;
249  //G4cout<<"cross2 = "<<cross2<<G4endl;
250  cross = (cross2-cross1);
251  //G4cout<<"cross = "<<cross<<G4endl;
252  if(!one) {
253  cross2 = (*table)(iPlace+1)->Value(tcut)/tcut
254  - (*table)(iPlace+1)->Value(tmax)/tmax;
255 
256  G4double E1 = fParticleEnergyVector->Energy(iPlace);
257  G4double E2 = fParticleEnergyVector->Energy(iPlace+1);
258  G4double W = 1.0/(E2 - E1);
259  G4double W1 = (E2 - scaledTkin)*W;
260  G4double W2 = (scaledTkin - E1)*W;
261  cross *= W1;
262  cross += W2*cross2;
263  }
264 
265  if( cross < 0.0) { cross = 0.0; }
266  return cross;
267 }
268 
269 ///////////////////////////////////////////////////////////////////////
270 
272  G4double kinEnergy,
273  G4double scaledTkin,
274  G4double stepFactor) const
275 {
276  G4double loss = 0.0;
277  G4double omega;
278  G4double position, E1, E2, W1, W2, W, dNdxCut1, dNdxCut2, meanNumber;
279 
280  size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
281  size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
282 
283  G4bool one = true;
284  if(scaledTkin >= fParticleEnergyVector->Energy(nPlace)) { iPlace = nPlace; }
285  else if(scaledTkin > fParticleEnergyVector->Energy(0)) {
286  one = false;
287  }
288  G4PhysicsLogVector* vcut = fdNdxCutTable[coupleIndex];
289  G4PhysicsVector* v1 = (*(fPAIxscBank[coupleIndex]))(iPlace);
290  G4PhysicsVector* v2 = 0;
291 
292  dNdxCut1 = (*vcut)[iPlace];
293  G4double e1 = v1->Energy(0);
294  G4double e2 = e1;
295 
296  G4double meanN1 = ((*v1)[0]/e1 - dNdxCut1)*stepFactor;
297  meanNumber = meanN1;
298 
299  //G4cout<<"iPlace = "<<iPlace<< " meanN1= " << meanN1
300  // << " dNdxCut1= " << dNdxCut1 << G4endl;
301 
302  dNdxCut2 = dNdxCut1;
303  W1 = 1.0;
304  W2 = 0.0;
305  if(!one) {
306  v2 = (*(fPAIxscBank[coupleIndex]))(iPlace+1);
307  dNdxCut2 = (*vcut)[iPlace+1];
308  e2 = v2->Energy(0);
309  G4double meanN2 = ((*v2)[0]/e2 - dNdxCut2)*stepFactor;
310  E1 = fParticleEnergyVector->Energy(iPlace);
311  E2 = fParticleEnergyVector->Energy(iPlace+1);
312  W = 1.0/(E2 - E1);
313  W1 = (E2 - scaledTkin)*W;
314  W2 = (scaledTkin - E1)*W;
315  meanNumber = W1*meanN1 + W2*meanN2;
316  //G4cout<<"meanN= " << meanNumber << " meanN2= " << meanN2
317  // << " dNdxCut2= " << dNdxCut2 << G4endl;
318  }
319  if(meanNumber < 0.0) { return 0.0; }
320 
321  G4int numOfCollisions = G4Poisson(meanNumber);
322 
323  //G4cout << "N= " << numOfCollisions << G4endl;
324 
325  if(0 == numOfCollisions) { return 0.0; }
326 
327  for(G4int i=0; i< numOfCollisions; ++i) {
328  G4double rand = G4UniformRand();
329  position = dNdxCut1 + ((*v1)[0]/e1 - dNdxCut1)*rand;
330  omega = GetEnergyTransfer(coupleIndex, iPlace, position);
331  //G4cout << "omega(keV)= " << omega/keV << G4endl;
332  if(!one) {
333  position = dNdxCut2 + ((*v2)[0]/e2 - dNdxCut2)*rand;
334  G4double omega2 = GetEnergyTransfer(coupleIndex, iPlace+1, position);
335  omega = omega*W1 + omega2*W2;
336  }
337  //G4cout << "omega(keV)= " << omega/keV << G4endl;
338 
339  loss += omega;
340  if(loss > kinEnergy) { break; }
341  }
342 
343  // G4cout<<"PAIModelData AlongStepLoss = "<<loss/keV<<" keV, on step = "
344  //<<step/mm<<" mm"<<G4endl;
345  if(loss > kinEnergy) { loss = kinEnergy; }
346  else if(loss < 0.) { loss = 0.; }
347  return loss;
348 }
349 
350 ///////////////////////////////////////////////////////////////////////
351 //
352 // Returns post step PAI energy transfer > cut electron energy
353 // according to passed scaled kinetic energy of particle
354 
356  G4double scaledTkin) const
357 {
358  //G4cout<<"G4PAIModelData::GetPostStepTransfer"<<G4endl;
359  G4double transfer = 0.0;
360  G4double rand = G4UniformRand();
361 
362  size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
363 
364  // size_t iTransfer, iTr1, iTr2;
365  G4double position, dNdxCut1, dNdxCut2, E1, E2, W1, W2, W;
366 
367  G4PhysicsVector* cutv = fdNdxCutTable[coupleIndex];
368 
369  // Fermi plato, try from left
370  if(scaledTkin >= fParticleEnergyVector->GetMaxEnergy())
371  {
372  position = (*cutv)[nPlace]*rand;
373  transfer = GetEnergyTransfer(coupleIndex, nPlace, position);
374  }
375  else if(scaledTkin <= fParticleEnergyVector->Energy(0))
376  {
377  position = (*cutv)[0]*rand;
378  transfer = GetEnergyTransfer(coupleIndex, 0, position);
379  }
380  else
381  {
382  size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
383 
384  dNdxCut1 = (*cutv)[iPlace];
385  dNdxCut2 = (*cutv)[iPlace+1];
386 
387  E1 = fParticleEnergyVector->Energy(iPlace);
388  E2 = fParticleEnergyVector->Energy(iPlace+1);
389  W = 1.0/(E2 - E1);
390  W1 = (E2 - scaledTkin)*W;
391  W2 = (scaledTkin - E1)*W;
392 
393  //G4cout<<"iPlace= " << " dNdxCut1 = "<<dNdxCut1
394  // <<" dNdxCut2 = "<<dNdxCut2<< " W1= " << W1 << " W2= " << W2 <<G4endl;
395 
396  position = dNdxCut1*rand;
397  G4double tr1 = GetEnergyTransfer(coupleIndex, iPlace, position);
398 
399  position = dNdxCut2*rand;
400  G4double tr2 = GetEnergyTransfer(coupleIndex, iPlace+1, position);
401 
402  transfer = tr1*W1 + tr2*W2;
403  }
404  //G4cout<<"PAImodel PostStepTransfer = "<<transfer/keV<<" keV"<<G4endl;
405  if(transfer < 0.0 ) { transfer = 0.0; }
406  return transfer;
407 }
408 
409 ///////////////////////////////////////////////////////////////////////
410 //
411 // Returns PAI energy transfer according to passed
412 // indexes of particle kinetic enegry and random x-section
413 
414 G4double G4PAIModelData::GetEnergyTransfer(G4int coupleIndex,
415  size_t iPlace,
416  G4double position) const
417 {
418  G4PhysicsVector* v = (*(fPAIxscBank[coupleIndex]))(iPlace);
419  if(position*v->Energy(0) >= (*v)[0]) { return v->Energy(0); }
420 
421  size_t iTransferMax = v->GetVectorLength() - 1;
422 
423  size_t iTransfer;
424  G4double x1(0.0), x2(0.0), y1(0.0), y2(0.0), energyTransfer;
425 
426  for(iTransfer=1; iTransfer<=iTransferMax; ++iTransfer) {
427  x2 = v->Energy(iTransfer);
428  y2 = (*v)[iTransfer]/x2;
429  if(position >= y2) { break; }
430  }
431 
432  x1 = v->Energy(iTransfer-1);
433  y1 = (*v)[iTransfer-1]/x1;
434  //G4cout << "i= " << iTransfer << " imax= " << iTransferMax
435  // << " x1= " << x1 << " x2= " << x2 << G4endl;
436 
437  energyTransfer = x1;
438  if ( x1 != x2 ) {
439  if ( y1 == y2 ) {
440  energyTransfer += (x2 - x1)*G4UniformRand();
441  } else {
442  if(x1*1.1 < x2) {
443  const G4int nbins = 5;
444  G4double del = (x2 - x1)/G4int(nbins);
445  x2 = x1;
446  for(G4int i=1; i<=nbins; ++i) {
447  x2 += del;
448  y2 = v->Value(x2)/x2;
449  if(position >= y2) { break; }
450  x1 = x2;
451  y1 = y2;
452  }
453  }
454  energyTransfer = (y2 - y1)*x1*x2/(position*(x1 - x2) - y1*x1 + y2*x2);
455  }
456  }
457  // G4cout << "x1(keV)= " << x1/keV << " x2(keV)= " << x2/keV
458  // << " y1= " << y1 << " y2= " << y2 << " pos= " << position
459  // << " E(keV)= " << energyTransfer/keV << G4endl;
460  return energyTransfer;
461 }
462 
463 //////////////////////////////////////////////////////////////////////
464 
G4long G4Poisson(G4double mean)
Definition: G4Poisson.hh:51
void PutValue(size_t binNumber, G4double binValue, G4double dataValue)
size_t GetVectorLength() const
int G4int
Definition: G4Types.hh:78
G4double SampleAlongStepTransfer(G4int coupleIndex, G4double kinEnergy, G4double scaledTkin, G4double stepFactor) const
G4double DEDXPerVolume(G4int coupleIndex, G4double scaledTkin, G4double cut) const
#define G4UniformRand()
Definition: Randomize.hh:87
G4GLOB_DLL std::ostream G4cout
bool G4bool
Definition: G4Types.hh:79
void PutValue(size_t index, G4double theValue)
const XML_Char XML_Content * model
float proton_mass_c2
Definition: hepunit.py:275
const G4int n
G4double Energy(size_t index) const
G4double Value(G4double theEnergy, size_t &lastidx) const
G4double CrossSectionPerVolume(G4int coupleIndex, G4double scaledTkin, G4double tcut, G4double tmax) const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
int position
Definition: filter.cc:7
void Initialise(const G4MaterialCutsCouple *, G4double cut, G4PAIModel *)
void insertAt(size_t, G4PhysicsVector *)
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
G4double SamplePostStepTransfer(G4int coupleIndex, G4double scaledTkin) const
G4PAIModelData(G4double tmin, G4double tmax, G4int verbose)
G4double ComputeMaxEnergy(G4double scaledEnergy)
Definition: G4PAIModel.hh:158
const G4Material * GetMaterial() const