| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27 // -------------------------------------------
28 //
29 // GEANT4 Class
30 // File name: G4PAIModelData.cc
31 //
32 // Author: V. Ivanchenko based on V.Gri
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
53 //////////////////////////////////////////////
54
55 using namespace std;
56
57 G4PAIModelData::G4PAIModelData(G4double tmin,
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, lowes
66 fHighestKineticEnergy = tmax;
67 if(tmax < 10*fLowestKineticEnergy) {
68 fHighestKineticEnergy = 10*fLowestKineticE
69 } else if(tmax > highestTkin) {
70 fHighestKineticEnergy = std::max(highestTk
71 }
72 fTotBin = (G4int)(nPerDecade*
73 std::log10(fHighestKineticEnergy/fLowe
74
75 fParticleEnergyVector = new G4PhysicsLogVect
76 fHighestKineticEnergy,
77 fTotBin);
78 if(0 < ver) {
79 G4cout << "### G4PAIModelData: Nbins= " <<
80 << " Tlowest(keV)= " << lowestTkin/keV
81 << " Tmin(keV)= " << fLowestKineticEnergy
82 << " Tmax(GeV)= " << fHighestKineticEnerg
83 << G4endl;
84 }
85 }
86
87 //////////////////////////////////////////////
88
89 G4PAIModelData::~G4PAIModelData()
90 {
91 std::size_t n = fPAIxscBank.size();
92 if(0 < n) {
93 for(std::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
108 //////////////////////////////////////////////
109
110 void G4PAIModelData::Initialise(const G4Materi
111 G4PAIModel* mo
112 {
113 const G4Material* mat = couple->GetMaterial(
114 fSandia.Initialize(const_cast<G4Material*>(m
115
116 auto PAItransferTable = new G4PhysicsTable(f
117 auto PAIdEdxTable = new G4PhysicsTable(fTotB
118 auto dEdxMeanVector =
119 new G4PhysicsLogVector(fLowestKineticEnerg
120 fHighestKineticEnergy,
121 fTotBin);
122 // low energy Sandia interval
123 G4double Tmin = fSandia.GetSandiaMatTablePAI
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
131 G4double Tmax = model->ComputeMaxEnergy(ki
132 G4double tau = kinEnergy/proton_mass_c2;
133 G4double bg2 = tau*( tau + 2. );
134
135 if (Tmax < Tmin + deltaLow ) { Tmax = Tmin
136
137 fPAIySection.Initialize(mat, Tmax, bg2, &f
138
139 //G4cout << i << ". TransferMax(keV)= "<<
140 // << " E(MeV)= " << kinEnergy/MeV <<
141
142 G4int n = fPAIySection.GetSplineSize();
143 G4int kmin = 0;
144 for(G4int k = 0; k < n; ++k) {
145 if(fPAIySection.GetIntegralPAIySection(k
146 kmin = k;
147 } else {
148 break;
149 }
150 }
151 n -= kmin;
152
153 auto transferVector = new G4PhysicsFreeVec
154 auto dEdxVector = new G4PhysicsFreeVector(
155
156 //G4double tr0 = 0.0;
157 G4double tr = 0.0;
158 for(G4int k = kmin; k < n; ++k)
159 {
160 G4double t = fPAIySection.GetSplineEner
161 tr = fPAIySection.GetIntegralPAIySection
162 //if(tr >= tr0) { tr0 = tr; }
163 //else { G4cout << "G4PAIModelData::Init
164 // << t/MeV << " IntegralTransfer
165 // << " < " << tr0 << G4endl; }
166 transferVector->PutValue(k, t, t*tr);
167 dEdxVector->PutValue(k, t, fPAIySection.
168 }
169 //G4cout << "TransferVector:" << G4endl;
170 //G4cout << *transferVector << G4endl;
171 //G4cout << "DEDXVector:" << G4endl;
172 //G4cout << *dEdxVector << G4endl;
173
174 G4double ionloss = std::max(fPAIySection.G
175 dEdxMeanVector->PutValue(i,ionloss);
176
177 PAItransferTable->insertAt(i,transferVecto
178 PAIdEdxTable->insertAt(i,dEdxVector);
179
180 } // end of Tkin loop`
181 fPAIxscBank.push_back(PAItransferTable);
182 fPAIdEdxBank.push_back(PAIdEdxTable);
183 //G4cout << "dEdxMeanVector: " << G4endl;
184 //G4cout << *dEdxMeanVector << G4endl;
185 fdEdxTable.push_back(dEdxMeanVector);
186 }
187
188 //////////////////////////////////////////////
189
190 G4double G4PAIModelData::DEDXPerVolume(G4int c
191 G4double cut) const
192 {
193 // VI: iPlace is the low edge index of the b
194 // iPlace is in interval from 0 to (N-1)
195 std::size_t iPlace(0);
196 G4double dEdx = fdEdxTable[coupleIndex]->Val
197 std::size_t nPlace = fParticleEnergyVector->
198 /*
199 G4cout << "G4PAIModelData::DEDXPerVolume: co
200 << " Tscaled= " << scaledTkin << " cut= " <
201 << " iPlace= " << iPlace << " nPlace= " <<
202 */
203 G4bool one = true;
204 if(scaledTkin >= fParticleEnergyVector->Ener
205 else if(scaledTkin > fParticleEnergyVector->
206 one = false;
207 }
208
209 // VI: apply interpolation of the vector
210 G4double del = (*(fPAIdEdxBank[coupleIndex]
211 //G4cout << "dEdx= " << dEdx << " del= " <<
212 if(!one) {
213 G4double del2 = (*(fPAIdEdxBank[coupleInde
214 G4double E1 = fParticleEnergyVector->Energ
215 G4double E2 = fParticleEnergyVector->Energ
216 G4double W = 1.0/(E2 - E1);
217 G4double W1 = (E2 - scaledTkin)*W;
218 G4double W2 = (scaledTkin - E1)*W;
219 del *= W1;
220 del += W2*del2;
221 }
222 dEdx -= del;
223 //G4cout << "dEdx= " << dEdx << " del= " <<
224
225 dEdx = std::max(dEdx, 0.);
226 return dEdx;
227 }
228
229 //////////////////////////////////////////////
230
231 G4double
232 G4PAIModelData::CrossSectionPerVolume(G4int co
233 G4double scaledTkin,
234 G4double tcut, G4double tmax) co
235 {
236 G4double cross, cross1, cross2;
237
238 // iPlace is in interval from 0 to (N-1)
239 std::size_t iPlace = fParticleEnergyVector->
240 std::size_t nPlace = fParticleEnergyVector->
241
242 G4bool one = true;
243 if(scaledTkin >= fParticleEnergyVector->Ener
244 else if(scaledTkin > fParticleEnergyVector->
245 one = false;
246 }
247 G4PhysicsTable* table = fPAIxscBank[coupleIn
248
249 //G4cout<<"iPlace = "<<iPlace<<"; tmax = "
250 // <<tmax<<"; cutEnergy = "<<cutEnergy<<G4en
251 cross1 = (*table)(iPlace)->Value(tmax)/tmax;
252 //G4cout<<"cross1 = "<<cross1<<G4endl;
253 cross2 = (*table)(iPlace)->Value(tcut)/tcut;
254 //G4cout<<"cross2 = "<<cross2<<G4endl;
255 cross = (cross2-cross1);
256 //G4cout<<"cross = "<<cross<<G4endl;
257 if(!one) {
258 cross2 = (*table)(iPlace+1)->Value(tcut)/t
259 - (*table)(iPlace+1)->Value(tmax)/tmax;
260
261 G4double E1 = fParticleEnergyVector->Energ
262 G4double E2 = fParticleEnergyVector->Energ
263 G4double W = 1.0/(E2 - E1);
264 G4double W1 = (E2 - scaledTkin)*W;
265 G4double W2 = (scaledTkin - E1)*W;
266 cross *= W1;
267 cross += W2*cross2;
268 }
269
270 cross = std::max(cross, 0.0);
271 return cross;
272 }
273
274 //////////////////////////////////////////////
275
276 G4double G4PAIModelData::SampleAlongStepTransf
277
278 G4double scaledTkin,
279 G4double tmax,
280 G4double stepFactor) const
281 {
282 //G4cout << "=== G4PAIModelData::SampleAlong
283 G4double loss = 0.0;
284
285 std::size_t iPlace = fParticleEnergyVector->
286 std::size_t nPlace = fParticleEnergyVector->
287
288 G4bool one = true;
289 if(scaledTkin >= fParticleEnergyVector->Ener
290 else if(scaledTkin > fParticleEnergyVector->
291 one = false;
292 }
293
294 G4double meanNumber = 0.0;
295 G4double meanN11 = 0.0;
296 G4double meanN12 = 0.0;
297 G4double meanN21 = 0.0;
298 G4double meanN22 = 0.0;
299
300 G4PhysicsVector* v1 = (*(fPAIxscBank[coupleI
301 G4PhysicsVector* v2 = nullptr;
302
303 G4double e1 = v1->Energy(0);
304 G4double e2 = std::min(tmax, v1->GetMaxEnerg
305
306 if(e2 >= e1) {
307 meanN11 = (*v1)[0]/e1;
308 meanN12 = v1->Value(e2)/e2;
309 meanNumber = (meanN11 - meanN12)*stepFacto
310 }
311 //G4cout<<"iPlace = "<<iPlace<< " meanN11= "
312 // << " meanN12= " << meanN12 << G4endl;
313
314 G4double W1 = 1.0;
315 G4double W2 = 0.0;
316 if(!one) {
317 v2 = (*(fPAIxscBank[coupleIndex]))(iPlace+
318
319 e1 = v2->Energy(0);
320 e2 = std::min(tmax, v2->GetMaxEnergy());
321 if(e2 >= e1) {
322 meanN21 = (*v2)[0]/e1;
323 meanN22 = v2->Value(e2)/e2;
324 G4double E1 = fParticleEnergyVector->Ene
325 G4double E2 = fParticleEnergyVector->Ene
326 G4double W = 1.0/(E2 - E1);
327 W1 = (E2 - scaledTkin)*W;
328 W2 = (scaledTkin - E1)*W;
329 meanNumber *= W1;
330 meanNumber += (meanN21 - meanN22)*stepFa
331 }
332 }
333
334 if(meanNumber < 0.0) { return loss; }
335 G4int numOfCollisions = (G4int)G4Poisson(mea
336
337 //G4cout << "meanNumber= " << meanNumber <<
338
339 if(0 == numOfCollisions) { return loss; }
340
341 G4double position, omega, omega2;
342 for(G4int i=0; i< numOfCollisions; ++i) {
343 G4double rand = G4UniformRand();
344 position = meanN12 + (meanN11 - meanN12)*r
345 omega = GetEnergyTransfer(coupleIndex, iPl
346 //G4cout << "omega(keV)= " << omega/keV <<
347 if(!one) {
348 position = meanN22 + (meanN21 - meanN22)
349 omega2 = GetEnergyTransfer(coupleIndex,
350 omega *= W1;
351 omega += omega2*W2;
352 }
353 //G4cout << "omega(keV)= " << omega/keV <<
354
355 loss += omega;
356 if(loss > kinEnergy) { break; }
357 }
358
359 //G4cout<<"PAIModelData AlongStepLoss = "<<l
360 if(loss > kinEnergy) { loss = kinEnergy; }
361 else if(loss < 0.) { loss = 0.; }
362 return loss;
363 }
364
365 //////////////////////////////////////////////
366 //
367 // Returns post step PAI energy transfer > cut
368 // according to passed scaled kinetic energy o
369
370 G4double G4PAIModelData::SamplePostStepTransfe
371 G4double scaledTkin,
372 G4double tmin,
373 G4double tmax) const
374 {
375 //G4cout<<"=== G4PAIModelData::SamplePostSte
376 // << " Tkin= " << scaledTkin << " Tmax= "
377 G4double transfer = 0.0;
378 G4double rand = G4UniformRand();
379
380 std::size_t nPlace = fParticleEnergyVector->
381 std::size_t iPlace = fParticleEnergyVector->
382
383 G4bool one = true;
384 if(scaledTkin >= fParticleEnergyVector->Ener
385 else if(scaledTkin > fParticleEnergyVector->
386 one = false;
387 }
388 G4PhysicsTable* table = fPAIxscBank[coupleIn
389 G4PhysicsVector* v1 = (*table)[iPlace];
390
391 G4double emin = std::max(tmin, v1->Energy(0)
392 G4double emax = std::min(tmax, v1->GetMaxEne
393 if(emax < emin) { return transfer; }
394
395 G4double dNdx1 = v1->Value(emin)/emin;
396 G4double dNdx2 = v1->Value(emax)/emax;
397 /*
398 G4cout << "iPlace= " << iPlace << " nPlace=
399 << " emin= " << emin << " emax= " << emax
400 << " dNdx1= " << dNdx1 << " dNdx2= " << dNd
401 << " one: " << one << G4endl;
402 */
403 G4double position = dNdx2 + (dNdx1 - dNdx2)*
404 transfer = GetEnergyTransfer(coupleIndex, iP
405
406 //G4cout<<"PAImodel PostStepTransfer = "<<tr
407 // << " position= " << position << G4endl;
408
409 if(!one) {
410
411 G4PhysicsVector* v2 = (*table)[iPlace+1];
412 emin = std::max(tmin, v2->Energy(0));
413 emax = std::min(tmax, v2->GetMaxEnergy());
414 if(emin <= emax) {
415 dNdx1 = v2->Value(emin)/emin;
416 dNdx2 = v2->Value(emax)/emax;
417
418 //G4cout << " emax2= " << emax
419 // << " dNdx2= " << dNdx2 << " dNdx1
420
421 G4double E1 = fParticleEnergyVector->Ene
422 G4double E2 = fParticleEnergyVector->Ene
423 G4double W = 1.0/(E2 - E1);
424 G4double W1 = (E2 - scaledTkin)*W;
425 G4double W2 = (scaledTkin - E1)*W;
426
427 //G4cout<< "E1= " << E1 << " E2= " << E2
428 // << " W1= " << W1 << " W2= " << W2
429
430 position = dNdx2 + (dNdx1 - dNdx2)*rand;
431 G4double tr2 = GetEnergyTransfer(coupleI
432
433 //G4cout<<"PAImodel PostStepTransfer1 =
434 // << " position= " << position << G4
435 transfer *= W1;
436 transfer += tr2*W2;
437 }
438 }
439 //G4cout<<"PAImodel PostStepTransfer = "<<tr
440 // << " position= " << position << G4endl;
441 transfer = std::max(transfer, 0.0);
442 return transfer;
443 }
444
445 //////////////////////////////////////////////
446 //
447 // Returns PAI energy transfer according to pa
448 // indexes of particle kinetic enegry and rand
449
450 G4double G4PAIModelData::GetEnergyTransfer(G4i
451 std::size_t iPlace,
452 G4double position) const
453 {
454 G4PhysicsVector* v = (*(fPAIxscBank[coupleIn
455 if(position*v->Energy(0) >= (*v)[0]) { retur
456
457 std::size_t iTransferMax = v->GetVectorLengt
458
459 std::size_t iTransfer;
460 G4double x1(0.0), x2(0.0), y1(0.0), y2(0.0),
461
462 //G4cout << "iPlace= " << iPlace << " iTrans
463 for(iTransfer=1; iTransfer<=iTransferMax; ++
464 x2 = v->Energy(iTransfer);
465 y2 = (*v)[iTransfer]/x2;
466 if(position >= y2) { break; }
467 if(iTransfer == iTransferMax) { return v->
468 }
469
470 x1 = v->Energy(iTransfer-1);
471 y1 = (*v)[iTransfer-1]/x1;
472 /*
473 G4cout << "i= " << iTransfer << " imax= " <<
474 << " x1= " << x1 << " x2= " << x2
475 << " y1= " << y1 << " y2= " << y2 << G4en
476 */
477 energyTransfer = x1;
478 if ( x1 != x2 ) {
479 if ( y1 == y2 ) {
480 energyTransfer += (x2 - x1)*G4UniformRan
481 } else {
482 if(x1*1.1 < x2) {
483 const G4int nbins = 5;
484 G4double del = (x2 - x1)/G4int(nbins);
485 x2 = x1;
486 for(G4int i=1; i<=nbins; ++i) {
487 x2 += del;
488 y2 = v->Value(x2)/x2;
489 if(position >= y2) {
490 break;
491 }
492 x1 = x2;
493 y1 = y2;
494 }
495 }
496 //G4cout << "x1(keV)= " << x1/keV << " x
497 // << " y1= " << y1 << " y2= " << y2 <
498 energyTransfer = (y2 - y1)*x1*x2/(positi
499 }
500 }
501 //G4cout << "x1(keV)= " << x1/keV << " x2(ke
502 // << " y1= " << y1 << " y2= " << y2 << "
503 // << " E(keV)= " << energyTransfer/keV <<
504 return energyTransfer;
505 }
506
507 //////////////////////////////////////////////
508
509