Geant4  10.03.p03
 All Classes 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 
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  << " Tlowest(keV)= " << lowestTkin/keV
82  << " Tmin(keV)= " << fLowestKineticEnergy/keV
83  << " Tmax(GeV)= " << fHighestKineticEnergy/GeV
84  << G4endl;
85  }
86 }
87 
89 
91 {
92  size_t n = fPAIxscBank.size();
93  if(0 < n) {
94  for(size_t i=0; i<n; ++i) {
95  if(fPAIxscBank[i]) {
96  fPAIxscBank[i]->clearAndDestroy();
97  delete fPAIxscBank[i];
98  }
99  if(fPAIdEdxBank[i]) {
100  fPAIdEdxBank[i]->clearAndDestroy();
101  delete fPAIdEdxBank[i];
102  }
103  delete fdEdxTable[i];
104  }
105  }
106  delete fParticleEnergyVector;
107 }
108 
110 
112  G4PAIModel* model)
113 {
114  const G4Material* mat = couple->GetMaterial();
115  fSandia.Initialize(const_cast<G4Material*>(mat));
116 
117  G4PhysicsTable* PAItransferTable = new G4PhysicsTable(fTotBin+1);
118  G4PhysicsTable* PAIdEdxTable = new G4PhysicsTable(fTotBin+1);
119  G4PhysicsLogVector* dEdxMeanVector =
120  new G4PhysicsLogVector(fLowestKineticEnergy,
121  fHighestKineticEnergy,
122  fTotBin);
123  // low energy Sandia interval
124  G4double Tmin = fSandia.GetSandiaMatTablePAI(0,0);
125 
126  // energy safety
127  static const G4double deltaLow = 100.*eV;
128 
129  for (G4int i = 0; i <= fTotBin; ++i) {
130 
131  G4double kinEnergy = fParticleEnergyVector->Energy(i);
132  G4double Tmax = model->ComputeMaxEnergy(kinEnergy);
133  G4double tau = kinEnergy/proton_mass_c2;
134  G4double bg2 = tau*( tau + 2. );
135 
136  if (Tmax < Tmin + deltaLow ) { Tmax = Tmin + deltaLow; }
137 
138  fPAIySection.Initialize(mat, Tmax, bg2, &fSandia);
139 
140  //G4cout << i << ". TransferMax(keV)= "<< Tmax/keV
141  // << " E(MeV)= " << kinEnergy/MeV << G4endl;
142 
143  G4int n = fPAIySection.GetSplineSize();
144  G4int kmin = 0;
145  for(G4int k = 0; k < n; ++k) {
146  if(fPAIySection.GetIntegralPAIySection(k+1) <= 0.0) {
147  kmin = k;
148  } else {
149  break;
150  }
151  }
152  n -= kmin;
153 
154  G4PhysicsFreeVector* transferVector = new G4PhysicsFreeVector(n);
155  G4PhysicsFreeVector* dEdxVector = new G4PhysicsFreeVector(n);
156 
157  //G4double tr0 = 0.0;
158  G4double tr = 0.0;
159  for(G4int k = kmin; k < n; ++k)
160  {
161  G4double t = fPAIySection.GetSplineEnergy(k+1);
162  tr = fPAIySection.GetIntegralPAIySection(k+1);
163  //if(tr >= tr0) { tr0 = tr; }
164  //else { G4cout << "G4PAIModelData::Initialise Warning: Ekin(MeV)= "
165  // << t/MeV << " IntegralTransfer= " << tr
166  // << " < " << tr0 << G4endl; }
167  transferVector->PutValue(k, t, t*tr);
168  dEdxVector->PutValue(k, t, fPAIySection.GetIntegralPAIdEdx(k+1));
169  }
170  //G4cout << "TransferVector:" << G4endl;
171  //G4cout << *transferVector << G4endl;
172  //G4cout << "DEDXVector:" << G4endl;
173  //G4cout << *dEdxVector << G4endl;
174 
175  G4double ionloss = fPAIySection.GetMeanEnergyLoss();// total <dE/dx>
176 
177  if(ionloss < 0.0) ionloss = 0.0;
178 
179  dEdxMeanVector->PutValue(i,ionloss);
180 
181  PAItransferTable->insertAt(i,transferVector);
182  PAIdEdxTable->insertAt(i,dEdxVector);
183 
184  //transferVector->SetSpline(true);
185  //transferVector->FillSecondDerivatives();
186  //dEdxVector->SetSpline(true);
187  //dEdxVector->FillSecondDerivatives();
188 
189  } // end of Tkin loop`
190  fPAIxscBank.push_back(PAItransferTable);
191  fPAIdEdxBank.push_back(PAIdEdxTable);
192  //G4cout << "dEdxMeanVector: " << G4endl;
193  //G4cout << *dEdxMeanVector << G4endl;
194  /*
195  dEdxMeanVector->SetSpline(true);
196  dEdxMeanVector->FillSecondDerivatives();
197  */
198  fdEdxTable.push_back(dEdxMeanVector);
199 }
200 
202 
204  G4double cut) const
205 {
206  // VI: iPlace is the low edge index of the bin
207  // iPlace is in interval from 0 to (N-1)
208  size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
209  size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
210  /*
211  G4cout << "G4PAIModelData::DEDXPerVolume: coupleIdx= " << coupleIndex
212  << " Tscaled= " << scaledTkin << " cut= " << cut
213  << " iPlace= " << iPlace << " nPlace= " << nPlace << G4endl;
214  */
215  G4bool one = true;
216  if(scaledTkin >= fParticleEnergyVector->Energy(nPlace)) { iPlace = nPlace; }
217  else if(scaledTkin > fParticleEnergyVector->Energy(0)) {
218  one = false;
219  }
220 
221  // VI: apply interpolation of the vector
222  G4double dEdx = fdEdxTable[coupleIndex]->Value(scaledTkin);
223  G4double del = (*(fPAIdEdxBank[coupleIndex]))(iPlace)->Value(cut);
224  //G4cout << "dEdx= " << dEdx << " del= " << del << G4endl;
225  if(!one) {
226  G4double del2 = (*(fPAIdEdxBank[coupleIndex]))(iPlace+1)->Value(cut);
227  G4double E1 = fParticleEnergyVector->Energy(iPlace);
228  G4double E2 = fParticleEnergyVector->Energy(iPlace+1);
229  G4double W = 1.0/(E2 - E1);
230  G4double W1 = (E2 - scaledTkin)*W;
231  G4double W2 = (scaledTkin - E1)*W;
232  del *= W1;
233  del += W2*del2;
234  }
235  dEdx -= del;
236  //G4cout << "dEdx= " << dEdx << " del= " << del << G4endl;
237 
238  dEdx = std::max(dEdx, 0.);
239  return dEdx;
240 }
241 
243 
244 G4double
246  G4double scaledTkin,
247  G4double tcut, G4double tmax) const
248 {
249  G4double cross, cross1, cross2;
250 
251  // iPlace is in interval from 0 to (N-1)
252  size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
253  size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
254 
255  G4bool one = true;
256  if(scaledTkin >= fParticleEnergyVector->Energy(nPlace)) { iPlace = nPlace; }
257  else if(scaledTkin > fParticleEnergyVector->Energy(0)) {
258  one = false;
259  }
260  G4PhysicsTable* table = fPAIxscBank[coupleIndex];
261 
262  //G4cout<<"iPlace = "<<iPlace<<"; tmax = "
263  // <<tmax<<"; cutEnergy = "<<cutEnergy<<G4endl;
264  cross1 = (*table)(iPlace)->Value(tmax)/tmax;
265  //G4cout<<"cross1 = "<<cross1<<G4endl;
266  cross2 = (*table)(iPlace)->Value(tcut)/tcut;
267  //G4cout<<"cross2 = "<<cross2<<G4endl;
268  cross = (cross2-cross1);
269  //G4cout<<"cross = "<<cross<<G4endl;
270  if(!one) {
271  cross2 = (*table)(iPlace+1)->Value(tcut)/tcut
272  - (*table)(iPlace+1)->Value(tmax)/tmax;
273 
274  G4double E1 = fParticleEnergyVector->Energy(iPlace);
275  G4double E2 = fParticleEnergyVector->Energy(iPlace+1);
276  G4double W = 1.0/(E2 - E1);
277  G4double W1 = (E2 - scaledTkin)*W;
278  G4double W2 = (scaledTkin - E1)*W;
279  cross *= W1;
280  cross += W2*cross2;
281  }
282 
283  cross = std::max(cross, 0.0);
284  return cross;
285 }
286 
288 
290  G4double kinEnergy,
291  G4double scaledTkin,
292  G4double tmax,
293  G4double stepFactor) const
294 {
295  //G4cout << "=== G4PAIModelData::SampleAlongStepTransfer" << G4endl;
296  G4double loss = 0.0;
297 
298  size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
299  size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
300 
301  G4bool one = true;
302  if(scaledTkin >= fParticleEnergyVector->Energy(nPlace)) { iPlace = nPlace; }
303  else if(scaledTkin > fParticleEnergyVector->Energy(0)) {
304  one = false;
305  }
306 
307  G4double meanNumber = 0.0;
308  G4double meanN11 = 0.0;
309  G4double meanN12 = 0.0;
310  G4double meanN21 = 0.0;
311  G4double meanN22 = 0.0;
312 
313  G4PhysicsVector* v1 = (*(fPAIxscBank[coupleIndex]))(iPlace);
314  G4PhysicsVector* v2 = 0;
315 
316  G4double e1 = v1->Energy(0);
317  G4double e2 = std::min(tmax, v1->GetMaxEnergy());
318 
319  if(e2 >= e1) {
320  meanN11 = (*v1)[0]/e1;
321  meanN12 = v1->Value(e2)/e2;
322  meanNumber = (meanN11 - meanN12)*stepFactor;
323  }
324  //G4cout<<"iPlace = "<<iPlace<< " meanN11= " << meanN11
325  // << " meanN12= " << meanN12 << G4endl;
326 
327  G4double W1 = 1.0;
328  G4double W2 = 0.0;
329  if(!one) {
330  v2 = (*(fPAIxscBank[coupleIndex]))(iPlace+1);
331 
332  e1 = v2->Energy(0);
333  e2 = std::min(tmax, v2->GetMaxEnergy());
334  if(e2 >= e1) {
335  meanN21 = (*v2)[0]/e1;
336  meanN22 = v2->Value(e2)/e2;
337  G4double E1 = fParticleEnergyVector->Energy(iPlace);
338  G4double E2 = fParticleEnergyVector->Energy(iPlace+1);
339  G4double W = 1.0/(E2 - E1);
340  W1 = (E2 - scaledTkin)*W;
341  W2 = (scaledTkin - E1)*W;
342  meanNumber *= W1;
343  meanNumber += (meanN21 - meanN22)*stepFactor*W2;
344  }
345  }
346 
347  if(meanNumber < 0.0) { return loss; }
348  G4int numOfCollisions = G4Poisson(meanNumber);
349 
350  //G4cout << "meanNumber= " << meanNumber << " N= " << numOfCollisions << G4endl;
351 
352  if(0 == numOfCollisions) { return loss; }
353 
354  G4double position, omega, omega2;
355  for(G4int i=0; i< numOfCollisions; ++i) {
356  G4double rand = G4UniformRand();
357  position = meanN12 + (meanN11 - meanN12)*rand;
358  omega = GetEnergyTransfer(coupleIndex, iPlace, position);
359  //G4cout << "omega(keV)= " << omega/keV << G4endl;
360  if(!one) {
361  position = meanN22 + (meanN21 - meanN22)*rand;
362  omega2 = GetEnergyTransfer(coupleIndex, iPlace+1, position);
363  omega *= W1;
364  omega += omega2*W2;
365  }
366  //G4cout << "omega(keV)= " << omega/keV << G4endl;
367 
368  loss += omega;
369  if(loss > kinEnergy) { break; }
370  }
371 
372  //G4cout<<"PAIModelData AlongStepLoss = "<<loss/keV<<" keV"<<G4endl;
373  if(loss > kinEnergy) { loss = kinEnergy; }
374  else if(loss < 0.) { loss = 0.; }
375  return loss;
376 }
377 
379 //
380 // Returns post step PAI energy transfer > cut electron energy
381 // according to passed scaled kinetic energy of particle
382 
384  G4double scaledTkin,
385  G4double tmin,
386  G4double tmax) const
387 {
388  //G4cout<<"=== G4PAIModelData::SamplePostStepTransfer idx= "<< coupleIndex
389  // << " Tkin= " << scaledTkin << " Tmax= " << tmax << G4endl;
390  G4double transfer = 0.0;
391  G4double rand = G4UniformRand();
392 
393  size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
394  size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
395 
396  G4bool one = true;
397  if(scaledTkin >= fParticleEnergyVector->Energy(nPlace)) { iPlace = nPlace; }
398  else if(scaledTkin > fParticleEnergyVector->Energy(0)) {
399  one = false;
400  }
401  G4PhysicsTable* table = fPAIxscBank[coupleIndex];
402  G4PhysicsVector* v1 = (*table)[iPlace];
403 
404  G4double emin = std::max(tmin, v1->Energy(0));
405  G4double emax = std::min(tmax, v1->GetMaxEnergy());
406  if(emax < emin) { return transfer; }
407 
408  G4double dNdx1 = v1->Value(emin)/emin;
409  G4double dNdx2 = v1->Value(emax)/emax;
410  /*
411  G4cout << "iPlace= " << iPlace << " nPlace= " << nPlace
412  << " emin= " << emin << " emax= " << emax
413  << " dNdx1= " << dNdx1 << " dNdx2= " << dNdx2
414  << " one: " << one << G4endl;
415  */
416  G4double position = dNdx2 + (dNdx1 - dNdx2)*rand;
417  transfer = GetEnergyTransfer(coupleIndex, iPlace, position);
418 
419  //G4cout<<"PAImodel PostStepTransfer = "<<transfer/keV<<" keV"
420  // << " position= " << position << G4endl;
421 
422  if(!one) {
423 
424  G4PhysicsVector* v2 = (*table)[iPlace+1];
425  emin = std::max(tmin, v2->Energy(0));
426  emax = std::min(tmax, v2->GetMaxEnergy());
427  if(emin <= emax) {
428  dNdx1 = v2->Value(emin)/emin;
429  dNdx2 = v2->Value(emax)/emax;
430 
431  //G4cout << " emax2= " << emax
432  // << " dNdx2= " << dNdx2 << " dNdx1= " << dNdx1 << G4endl;
433 
434  G4double E1 = fParticleEnergyVector->Energy(iPlace);
435  G4double E2 = fParticleEnergyVector->Energy(iPlace+1);
436  G4double W = 1.0/(E2 - E1);
437  G4double W1 = (E2 - scaledTkin)*W;
438  G4double W2 = (scaledTkin - E1)*W;
439 
440  //G4cout<< "E1= " << E1 << " E2= " << E2 <<" iPlace= " << iPlace
441  // << " W1= " << W1 << " W2= " << W2 <<G4endl;
442 
443  position = dNdx2 + (dNdx1 - dNdx2)*rand;
444  G4double tr2 = GetEnergyTransfer(coupleIndex, iPlace+1, position);
445 
446  //G4cout<<"PAImodel PostStepTransfer1 = "<<tr2/keV<<" keV"
447  // << " position= " << position << G4endl;
448  transfer *= W1;
449  transfer += tr2*W2;
450  }
451  }
452  //G4cout<<"PAImodel PostStepTransfer = "<<transfer/keV<<" keV"
453  // << " position= " << position << G4endl;
454  transfer = std::max(transfer, 0.0);
455  return transfer;
456 }
457 
459 //
460 // Returns PAI energy transfer according to passed
461 // indexes of particle kinetic enegry and random x-section
462 
463 G4double G4PAIModelData::GetEnergyTransfer(G4int coupleIndex,
464  size_t iPlace,
465  G4double position) const
466 {
467  G4PhysicsVector* v = (*(fPAIxscBank[coupleIndex]))(iPlace);
468  if(position*v->Energy(0) >= (*v)[0]) { return v->Energy(0); }
469 
470  size_t iTransferMax = v->GetVectorLength() - 1;
471 
472  size_t iTransfer;
473  G4double x1(0.0), x2(0.0), y1(0.0), y2(0.0), energyTransfer;
474 
475  //G4cout << "iPlace= " << iPlace << " iTransferMax= " << iTransferMax << G4endl;
476  for(iTransfer=1; iTransfer<=iTransferMax; ++iTransfer) {
477  x2 = v->Energy(iTransfer);
478  y2 = (*v)[iTransfer]/x2;
479  if(position >= y2) { break; }
480  if(iTransfer == iTransferMax) { return v->GetMaxEnergy(); }
481  }
482 
483  x1 = v->Energy(iTransfer-1);
484  y1 = (*v)[iTransfer-1]/x1;
485  /*
486  G4cout << "i= " << iTransfer << " imax= " << iTransferMax
487  << " x1= " << x1 << " x2= " << x2
488  << " y1= " << y1 << " y2= " << y2 << G4endl;
489  */
490  energyTransfer = x1;
491  if ( x1 != x2 ) {
492  if ( y1 == y2 ) {
493  energyTransfer += (x2 - x1)*G4UniformRand();
494  } else {
495  if(x1*1.1 < x2) {
496  const G4int nbins = 5;
497  G4double del = (x2 - x1)/G4int(nbins);
498  x2 = x1;
499  for(G4int i=1; i<=nbins; ++i) {
500  x2 += del;
501  y2 = v->Value(x2)/x2;
502  if(position >= y2) {
503  break;
504  }
505  x1 = x2;
506  y1 = y2;
507  }
508  }
509  //G4cout << "x1(keV)= " << x1/keV << " x2(keV)= " << x2/keV
510  // << " y1= " << y1 << " y2= " << y2 << " pos= " << position << G4endl;
511  energyTransfer = (y2 - y1)*x1*x2/(position*(x1 - x2) - y1*x1 + y2*x2);
512  }
513  }
514  //G4cout << "x1(keV)= " << x1/keV << " x2(keV)= " << x2/keV
515  // << " y1= " << y1 << " y2= " << y2 << " pos= " << position
516  // << " E(keV)= " << energyTransfer/keV << G4endl;
517  return energyTransfer;
518 }
519 
521 
G4long G4Poisson(G4double mean)
Definition: G4Poisson.hh:51
G4double GetMaxEnergy() const
void PutValue(size_t index, G4double energy, G4double dataValue)
size_t GetVectorLength() const
int G4int
Definition: G4Types.hh:78
G4double SampleAlongStepTransfer(G4int coupleIndex, G4double kinEnergy, G4double scaledTkin, G4double tmax, G4double stepFactor) const
static constexpr double TeV
Definition: G4SIunits.hh:218
G4double DEDXPerVolume(G4int coupleIndex, G4double scaledTkin, G4double cut) const
#define position
Definition: xmlparse.cc:622
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
bool G4bool
Definition: G4Types.hh:79
void PutValue(size_t index, G4double theValue)
static constexpr double eV
Definition: G4SIunits.hh:215
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
tuple v
Definition: test.py:18
static const G4double emax
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
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
static constexpr double GeV
Definition: G4SIunits.hh:217
void insertAt(size_t, G4PhysicsVector *)
#define G4endl
Definition: G4ios.hh:61
G4double SamplePostStepTransfer(G4int coupleIndex, G4double scaledTkin, G4double tmin, G4double tmax) const
double G4double
Definition: G4Types.hh:76
void Initialise(const G4MaterialCutsCouple *, G4PAIModel *)
static constexpr double keV
Definition: G4SIunits.hh:216
const XML_Char XML_Content * model
Definition: expat.h:151
G4PAIModelData(G4double tmin, G4double tmax, G4int verbose)
G4double ComputeMaxEnergy(G4double scaledEnergy)
Definition: G4PAIModel.hh:166
const G4Material * GetMaterial() const