ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4PAIModelData.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4PAIModelData.cc
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 //
27 // -------------------------------------------------------------------
28 //
29 // GEANT4 Class
30 // File name: G4PAIModelData.cc
31 //
32 // Author: V. Ivanchenko based on V.Grichine code of G4PAIModel
33 //
34 // Creation date: 16.08.2013
35 //
36 // Modifications:
37 //
38 
39 #include "G4PAIModelData.hh"
40 #include "G4PAIModel.hh"
41 
42 #include "G4SystemOfUnits.hh"
43 #include "G4PhysicalConstants.hh"
44 
45 #include "G4PhysicsLogVector.hh"
46 #include "G4PhysicsFreeVector.hh"
47 #include "G4PhysicsTable.hh"
48 #include "G4MaterialCutsCouple.hh"
49 #include "G4SandiaTable.hh"
50 #include "Randomize.hh"
51 #include "G4Poisson.hh"
52 
54 
55 using namespace std;
56 
57 G4PAIModelData::G4PAIModelData(G4double tmin, G4double tmax, G4int ver)
58 {
59  const G4int nPerDecade = 10;
60  const G4double lowestTkin = 50*keV;
61  const G4double highestTkin = 10*TeV;
62 
63  fPAIySection.SetVerbose(ver);
64 
65  fLowestKineticEnergy = std::max(tmin, lowestTkin);
66  fHighestKineticEnergy = tmax;
67  if(tmax < 10*fLowestKineticEnergy) {
68  fHighestKineticEnergy = 10*fLowestKineticEnergy;
69  } else if(tmax > highestTkin) {
70  fHighestKineticEnergy = std::max(highestTkin, 10*fLowestKineticEnergy);
71  }
72  fTotBin = (G4int)(nPerDecade*
73  std::log10(fHighestKineticEnergy/fLowestKineticEnergy));
74 
75  fParticleEnergyVector = new G4PhysicsLogVector(fLowestKineticEnergy,
76  fHighestKineticEnergy,
77  fTotBin);
78  if(0 < ver) {
79  G4cout << "### G4PAIModelData: Nbins= " << fTotBin
80  << " Tlowest(keV)= " << lowestTkin/keV
81  << " Tmin(keV)= " << fLowestKineticEnergy/keV
82  << " Tmax(GeV)= " << fHighestKineticEnergy/GeV
83  << G4endl;
84  }
85 }
86 
88 
89 G4PAIModelData::~G4PAIModelData()
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  fPAIxscBank[i]->clearAndDestroy();
96  delete fPAIxscBank[i];
97  }
98  if(fPAIdEdxBank[i]) {
99  fPAIdEdxBank[i]->clearAndDestroy();
100  delete fPAIdEdxBank[i];
101  }
102  delete fdEdxTable[i];
103  }
104  }
105  delete fParticleEnergyVector;
106 }
107 
109 
110 void G4PAIModelData::Initialise(const G4MaterialCutsCouple* couple,
111  G4PAIModel* model)
112 {
113  const G4Material* mat = couple->GetMaterial();
114  fSandia.Initialize(const_cast<G4Material*>(mat));
115 
116  G4PhysicsTable* PAItransferTable = new G4PhysicsTable(fTotBin+1);
117  G4PhysicsTable* PAIdEdxTable = new G4PhysicsTable(fTotBin+1);
118  G4PhysicsLogVector* dEdxMeanVector =
119  new G4PhysicsLogVector(fLowestKineticEnergy,
120  fHighestKineticEnergy,
121  fTotBin);
122  // low energy Sandia interval
123  G4double Tmin = fSandia.GetSandiaMatTablePAI(0,0);
124 
125  // energy safety
126  static const G4double deltaLow = 100.*eV;
127 
128  for (G4int i = 0; i <= fTotBin; ++i) {
129 
130  G4double kinEnergy = fParticleEnergyVector->Energy(i);
131  G4double Tmax = model->ComputeMaxEnergy(kinEnergy);
132  G4double tau = kinEnergy/proton_mass_c2;
133  G4double bg2 = tau*( tau + 2. );
134 
135  if (Tmax < Tmin + deltaLow ) { Tmax = Tmin + deltaLow; }
136 
137  fPAIySection.Initialize(mat, Tmax, bg2, &fSandia);
138 
139  //G4cout << i << ". TransferMax(keV)= "<< Tmax/keV
140  // << " E(MeV)= " << kinEnergy/MeV << G4endl;
141 
142  G4int n = fPAIySection.GetSplineSize();
143  G4int kmin = 0;
144  for(G4int k = 0; k < n; ++k) {
145  if(fPAIySection.GetIntegralPAIySection(k+1) <= 0.0) {
146  kmin = k;
147  } else {
148  break;
149  }
150  }
151  n -= kmin;
152 
153  G4PhysicsFreeVector* transferVector = new G4PhysicsFreeVector(n);
154  G4PhysicsFreeVector* dEdxVector = new G4PhysicsFreeVector(n);
155 
156  //G4double tr0 = 0.0;
157  G4double tr = 0.0;
158  for(G4int k = kmin; k < n; ++k)
159  {
160  G4double t = fPAIySection.GetSplineEnergy(k+1);
161  tr = fPAIySection.GetIntegralPAIySection(k+1);
162  //if(tr >= tr0) { tr0 = tr; }
163  //else { G4cout << "G4PAIModelData::Initialise Warning: Ekin(MeV)= "
164  // << t/MeV << " IntegralTransfer= " << tr
165  // << " < " << tr0 << G4endl; }
166  transferVector->PutValue(k, t, t*tr);
167  dEdxVector->PutValue(k, t, fPAIySection.GetIntegralPAIdEdx(k+1));
168  }
169  //G4cout << "TransferVector:" << G4endl;
170  //G4cout << *transferVector << G4endl;
171  //G4cout << "DEDXVector:" << G4endl;
172  //G4cout << *dEdxVector << G4endl;
173 
174  G4double ionloss = fPAIySection.GetMeanEnergyLoss();// total <dE/dx>
175 
176  if(ionloss < 0.0) ionloss = 0.0;
177 
178  dEdxMeanVector->PutValue(i,ionloss);
179 
180  PAItransferTable->insertAt(i,transferVector);
181  PAIdEdxTable->insertAt(i,dEdxVector);
182 
183  //transferVector->SetSpline(true);
184  //transferVector->FillSecondDerivatives();
185  //dEdxVector->SetSpline(true);
186  //dEdxVector->FillSecondDerivatives();
187 
188  } // end of Tkin loop`
189  fPAIxscBank.push_back(PAItransferTable);
190  fPAIdEdxBank.push_back(PAIdEdxTable);
191  //G4cout << "dEdxMeanVector: " << G4endl;
192  //G4cout << *dEdxMeanVector << G4endl;
193  /*
194  dEdxMeanVector->SetSpline(true);
195  dEdxMeanVector->FillSecondDerivatives();
196  */
197  fdEdxTable.push_back(dEdxMeanVector);
198 }
199 
201 
202 G4double G4PAIModelData::DEDXPerVolume(G4int coupleIndex, G4double scaledTkin,
203  G4double cut) const
204 {
205  // VI: iPlace is the low edge index of the bin
206  // iPlace is in interval from 0 to (N-1)
207  size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
208  size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
209  /*
210  G4cout << "G4PAIModelData::DEDXPerVolume: coupleIdx= " << coupleIndex
211  << " Tscaled= " << scaledTkin << " cut= " << cut
212  << " iPlace= " << iPlace << " nPlace= " << nPlace << G4endl;
213  */
214  G4bool one = true;
215  if(scaledTkin >= fParticleEnergyVector->Energy(nPlace)) { iPlace = nPlace; }
216  else if(scaledTkin > fParticleEnergyVector->Energy(0)) {
217  one = false;
218  }
219 
220  // VI: apply interpolation of the vector
221  G4double dEdx = fdEdxTable[coupleIndex]->Value(scaledTkin);
222  G4double del = (*(fPAIdEdxBank[coupleIndex]))(iPlace)->Value(cut);
223  //G4cout << "dEdx= " << dEdx << " del= " << del << G4endl;
224  if(!one) {
225  G4double del2 = (*(fPAIdEdxBank[coupleIndex]))(iPlace+1)->Value(cut);
226  G4double E1 = fParticleEnergyVector->Energy(iPlace);
227  G4double E2 = fParticleEnergyVector->Energy(iPlace+1);
228  G4double W = 1.0/(E2 - E1);
229  G4double W1 = (E2 - scaledTkin)*W;
230  G4double W2 = (scaledTkin - E1)*W;
231  del *= W1;
232  del += W2*del2;
233  }
234  dEdx -= del;
235  //G4cout << "dEdx= " << dEdx << " del= " << del << G4endl;
236 
237  dEdx = std::max(dEdx, 0.);
238  return dEdx;
239 }
240 
242 
243 G4double
244 G4PAIModelData::CrossSectionPerVolume(G4int coupleIndex,
245  G4double scaledTkin,
246  G4double tcut, G4double tmax) const
247 {
248  G4double cross, cross1, cross2;
249 
250  // iPlace is in interval from 0 to (N-1)
251  size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
252  size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
253 
254  G4bool one = true;
255  if(scaledTkin >= fParticleEnergyVector->Energy(nPlace)) { iPlace = nPlace; }
256  else if(scaledTkin > fParticleEnergyVector->Energy(0)) {
257  one = false;
258  }
259  G4PhysicsTable* table = fPAIxscBank[coupleIndex];
260 
261  //G4cout<<"iPlace = "<<iPlace<<"; tmax = "
262  // <<tmax<<"; cutEnergy = "<<cutEnergy<<G4endl;
263  cross1 = (*table)(iPlace)->Value(tmax)/tmax;
264  //G4cout<<"cross1 = "<<cross1<<G4endl;
265  cross2 = (*table)(iPlace)->Value(tcut)/tcut;
266  //G4cout<<"cross2 = "<<cross2<<G4endl;
267  cross = (cross2-cross1);
268  //G4cout<<"cross = "<<cross<<G4endl;
269  if(!one) {
270  cross2 = (*table)(iPlace+1)->Value(tcut)/tcut
271  - (*table)(iPlace+1)->Value(tmax)/tmax;
272 
273  G4double E1 = fParticleEnergyVector->Energy(iPlace);
274  G4double E2 = fParticleEnergyVector->Energy(iPlace+1);
275  G4double W = 1.0/(E2 - E1);
276  G4double W1 = (E2 - scaledTkin)*W;
277  G4double W2 = (scaledTkin - E1)*W;
278  cross *= W1;
279  cross += W2*cross2;
280  }
281 
282  cross = std::max(cross, 0.0);
283  return cross;
284 }
285 
287 
288 G4double G4PAIModelData::SampleAlongStepTransfer(G4int coupleIndex,
289  G4double kinEnergy,
290  G4double scaledTkin,
291  G4double tmax,
292  G4double stepFactor) const
293 {
294  //G4cout << "=== G4PAIModelData::SampleAlongStepTransfer" << G4endl;
295  G4double loss = 0.0;
296 
297  size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
298  size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
299 
300  G4bool one = true;
301  if(scaledTkin >= fParticleEnergyVector->Energy(nPlace)) { iPlace = nPlace; }
302  else if(scaledTkin > fParticleEnergyVector->Energy(0)) {
303  one = false;
304  }
305 
306  G4double meanNumber = 0.0;
307  G4double meanN11 = 0.0;
308  G4double meanN12 = 0.0;
309  G4double meanN21 = 0.0;
310  G4double meanN22 = 0.0;
311 
312  G4PhysicsVector* v1 = (*(fPAIxscBank[coupleIndex]))(iPlace);
313  G4PhysicsVector* v2 = 0;
314 
315  G4double e1 = v1->Energy(0);
316  G4double e2 = std::min(tmax, v1->GetMaxEnergy());
317 
318  if(e2 >= e1) {
319  meanN11 = (*v1)[0]/e1;
320  meanN12 = v1->Value(e2)/e2;
321  meanNumber = (meanN11 - meanN12)*stepFactor;
322  }
323  //G4cout<<"iPlace = "<<iPlace<< " meanN11= " << meanN11
324  // << " meanN12= " << meanN12 << G4endl;
325 
326  G4double W1 = 1.0;
327  G4double W2 = 0.0;
328  if(!one) {
329  v2 = (*(fPAIxscBank[coupleIndex]))(iPlace+1);
330 
331  e1 = v2->Energy(0);
332  e2 = std::min(tmax, v2->GetMaxEnergy());
333  if(e2 >= e1) {
334  meanN21 = (*v2)[0]/e1;
335  meanN22 = v2->Value(e2)/e2;
336  G4double E1 = fParticleEnergyVector->Energy(iPlace);
337  G4double E2 = fParticleEnergyVector->Energy(iPlace+1);
338  G4double W = 1.0/(E2 - E1);
339  W1 = (E2 - scaledTkin)*W;
340  W2 = (scaledTkin - E1)*W;
341  meanNumber *= W1;
342  meanNumber += (meanN21 - meanN22)*stepFactor*W2;
343  }
344  }
345 
346  if(meanNumber < 0.0) { return loss; }
347  G4int numOfCollisions = G4Poisson(meanNumber);
348 
349  //G4cout << "meanNumber= " << meanNumber << " N= " << numOfCollisions << G4endl;
350 
351  if(0 == numOfCollisions) { return loss; }
352 
353  G4double position, omega, omega2;
354  for(G4int i=0; i< numOfCollisions; ++i) {
355  G4double rand = G4UniformRand();
356  position = meanN12 + (meanN11 - meanN12)*rand;
357  omega = GetEnergyTransfer(coupleIndex, iPlace, position);
358  //G4cout << "omega(keV)= " << omega/keV << G4endl;
359  if(!one) {
360  position = meanN22 + (meanN21 - meanN22)*rand;
361  omega2 = GetEnergyTransfer(coupleIndex, iPlace+1, position);
362  omega *= W1;
363  omega += omega2*W2;
364  }
365  //G4cout << "omega(keV)= " << omega/keV << G4endl;
366 
367  loss += omega;
368  if(loss > kinEnergy) { break; }
369  }
370 
371  //G4cout<<"PAIModelData AlongStepLoss = "<<loss/keV<<" keV"<<G4endl;
372  if(loss > kinEnergy) { loss = kinEnergy; }
373  else if(loss < 0.) { loss = 0.; }
374  return loss;
375 }
376 
378 //
379 // Returns post step PAI energy transfer > cut electron energy
380 // according to passed scaled kinetic energy of particle
381 
382 G4double G4PAIModelData::SamplePostStepTransfer(G4int coupleIndex,
383  G4double scaledTkin,
384  G4double tmin,
385  G4double tmax) const
386 {
387  //G4cout<<"=== G4PAIModelData::SamplePostStepTransfer idx= "<< coupleIndex
388  // << " Tkin= " << scaledTkin << " Tmax= " << tmax << G4endl;
389  G4double transfer = 0.0;
390  G4double rand = G4UniformRand();
391 
392  size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
393  size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
394 
395  G4bool one = true;
396  if(scaledTkin >= fParticleEnergyVector->Energy(nPlace)) { iPlace = nPlace; }
397  else if(scaledTkin > fParticleEnergyVector->Energy(0)) {
398  one = false;
399  }
400  G4PhysicsTable* table = fPAIxscBank[coupleIndex];
401  G4PhysicsVector* v1 = (*table)[iPlace];
402 
403  G4double emin = std::max(tmin, v1->Energy(0));
404  G4double emax = std::min(tmax, v1->GetMaxEnergy());
405  if(emax < emin) { return transfer; }
406 
407  G4double dNdx1 = v1->Value(emin)/emin;
408  G4double dNdx2 = v1->Value(emax)/emax;
409  /*
410  G4cout << "iPlace= " << iPlace << " nPlace= " << nPlace
411  << " emin= " << emin << " emax= " << emax
412  << " dNdx1= " << dNdx1 << " dNdx2= " << dNdx2
413  << " one: " << one << G4endl;
414  */
415  G4double position = dNdx2 + (dNdx1 - dNdx2)*rand;
416  transfer = GetEnergyTransfer(coupleIndex, iPlace, position);
417 
418  //G4cout<<"PAImodel PostStepTransfer = "<<transfer/keV<<" keV"
419  // << " position= " << position << G4endl;
420 
421  if(!one) {
422 
423  G4PhysicsVector* v2 = (*table)[iPlace+1];
424  emin = std::max(tmin, v2->Energy(0));
425  emax = std::min(tmax, v2->GetMaxEnergy());
426  if(emin <= emax) {
427  dNdx1 = v2->Value(emin)/emin;
428  dNdx2 = v2->Value(emax)/emax;
429 
430  //G4cout << " emax2= " << emax
431  // << " dNdx2= " << dNdx2 << " dNdx1= " << dNdx1 << G4endl;
432 
433  G4double E1 = fParticleEnergyVector->Energy(iPlace);
434  G4double E2 = fParticleEnergyVector->Energy(iPlace+1);
435  G4double W = 1.0/(E2 - E1);
436  G4double W1 = (E2 - scaledTkin)*W;
437  G4double W2 = (scaledTkin - E1)*W;
438 
439  //G4cout<< "E1= " << E1 << " E2= " << E2 <<" iPlace= " << iPlace
440  // << " W1= " << W1 << " W2= " << W2 <<G4endl;
441 
442  position = dNdx2 + (dNdx1 - dNdx2)*rand;
443  G4double tr2 = GetEnergyTransfer(coupleIndex, iPlace+1, position);
444 
445  //G4cout<<"PAImodel PostStepTransfer1 = "<<tr2/keV<<" keV"
446  // << " position= " << position << G4endl;
447  transfer *= W1;
448  transfer += tr2*W2;
449  }
450  }
451  //G4cout<<"PAImodel PostStepTransfer = "<<transfer/keV<<" keV"
452  // << " position= " << position << G4endl;
453  transfer = std::max(transfer, 0.0);
454  return transfer;
455 }
456 
458 //
459 // Returns PAI energy transfer according to passed
460 // indexes of particle kinetic enegry and random x-section
461 
462 G4double G4PAIModelData::GetEnergyTransfer(G4int coupleIndex,
463  size_t iPlace,
464  G4double position) const
465 {
466  G4PhysicsVector* v = (*(fPAIxscBank[coupleIndex]))(iPlace);
467  if(position*v->Energy(0) >= (*v)[0]) { return v->Energy(0); }
468 
469  size_t iTransferMax = v->GetVectorLength() - 1;
470 
471  size_t iTransfer;
472  G4double x1(0.0), x2(0.0), y1(0.0), y2(0.0), energyTransfer;
473 
474  //G4cout << "iPlace= " << iPlace << " iTransferMax= " << iTransferMax << G4endl;
475  for(iTransfer=1; iTransfer<=iTransferMax; ++iTransfer) {
476  x2 = v->Energy(iTransfer);
477  y2 = (*v)[iTransfer]/x2;
478  if(position >= y2) { break; }
479  if(iTransfer == iTransferMax) { return v->GetMaxEnergy(); }
480  }
481 
482  x1 = v->Energy(iTransfer-1);
483  y1 = (*v)[iTransfer-1]/x1;
484  /*
485  G4cout << "i= " << iTransfer << " imax= " << iTransferMax
486  << " x1= " << x1 << " x2= " << x2
487  << " y1= " << y1 << " y2= " << y2 << G4endl;
488  */
489  energyTransfer = x1;
490  if ( x1 != x2 ) {
491  if ( y1 == y2 ) {
492  energyTransfer += (x2 - x1)*G4UniformRand();
493  } else {
494  if(x1*1.1 < x2) {
495  const G4int nbins = 5;
496  G4double del = (x2 - x1)/G4int(nbins);
497  x2 = x1;
498  for(G4int i=1; i<=nbins; ++i) {
499  x2 += del;
500  y2 = v->Value(x2)/x2;
501  if(position >= y2) {
502  break;
503  }
504  x1 = x2;
505  y1 = y2;
506  }
507  }
508  //G4cout << "x1(keV)= " << x1/keV << " x2(keV)= " << x2/keV
509  // << " y1= " << y1 << " y2= " << y2 << " pos= " << position << G4endl;
510  energyTransfer = (y2 - y1)*x1*x2/(position*(x1 - x2) - y1*x1 + y2*x2);
511  }
512  }
513  //G4cout << "x1(keV)= " << x1/keV << " x2(keV)= " << x2/keV
514  // << " y1= " << y1 << " y2= " << y2 << " pos= " << position
515  // << " E(keV)= " << energyTransfer/keV << G4endl;
516  return energyTransfer;
517 }
518 
520