Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
G4PAIPhotData.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 // $Id: G4PAIPhotData.cc 72008 2013-07-03 08:46:39Z vnivanch $
27 //
28 // -------------------------------------------------------------------
29 //
30 // GEANT4 Class
31 // File name: G4PAIPhotData.cc
32 //
33 // Author: V.Grichine based on G4PAIModelData MT
34 //
35 // Creation date: 07.10.2013
36 //
37 // Modifications:
38 //
39 
40 #include "G4PAIPhotData.hh"
41 #include "G4PAIPhotModel.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 "G4ProductionCutsTable.hh"
51 #include "G4SandiaTable.hh"
52 #include "Randomize.hh"
53 #include "G4Poisson.hh"
54 
56 
57 using namespace std;
58 
60 {
61  const G4int nPerDecade = 10;
62  const G4double lowestTkin = 50*keV;
63  const G4double highestTkin = 10*TeV;
64 
65  // fPAIxSection.SetVerbose(ver);
66 
67  fLowestKineticEnergy = std::max(tmin, lowestTkin);
68  fHighestKineticEnergy = tmax;
69 
70  if(tmax < 10*fLowestKineticEnergy)
71  {
72  fHighestKineticEnergy = 10*fLowestKineticEnergy;
73  }
74  else if(tmax > highestTkin)
75  {
76  fHighestKineticEnergy = std::max(highestTkin, 10*fLowestKineticEnergy);
77  }
78  fTotBin = (G4int)(nPerDecade*
79  std::log10(fHighestKineticEnergy/fLowestKineticEnergy));
80 
81  fParticleEnergyVector = new G4PhysicsLogVector(fLowestKineticEnergy,
82  fHighestKineticEnergy,
83  fTotBin);
84  if(0 < ver) {
85  G4cout << "### G4PAIPhotData: Nbins= " << fTotBin
86  << " Tmin(MeV)= " << fLowestKineticEnergy/MeV
87  << " Tmax(GeV)= " << fHighestKineticEnergy/GeV
88  << " tmin(keV)= " << tmin/keV << G4endl;
89  }
90 }
91 
93 
95 {
96  //G4cout << "G4PAIPhotData::~G4PAIPhotData() " << this << G4endl;
97  size_t n = fPAIxscBank.size();
98  if(0 < n)
99  {
100  for(size_t i=0; i<n; ++i)
101  {
102  if(fPAIxscBank[i])
103  {
104  fPAIxscBank[i]->clearAndDestroy();
105  delete fPAIxscBank[i];
106  fPAIxscBank[i] = 0;
107  }
108  if(fPAIdEdxBank[i])
109  {
110  fPAIdEdxBank[i]->clearAndDestroy();
111  delete fPAIdEdxBank[i];
112  fPAIdEdxBank[i]= 0;
113  }
114  delete fdEdxTable[i];
115  delete fdNdxCutTable[i];
116  fdEdxTable[i] = 0;
117  fdNdxCutTable[i] = 0;
118  }
119  }
120  delete fParticleEnergyVector;
121  fParticleEnergyVector = 0;
122  //G4cout << "G4PAIPhotData::~G4PAIPhotData() done for " << this << G4endl;
123 }
124 
126 
129 {
130  G4ProductionCutsTable* theCoupleTable=
132  size_t numOfCouples = theCoupleTable->GetTableSize();
133  size_t jMatCC;
134 
135  for (jMatCC = 0; jMatCC < numOfCouples; jMatCC++ )
136  {
137  if( couple == theCoupleTable->GetMaterialCutsCouple(jMatCC) ) break;
138  }
139  if( jMatCC == numOfCouples && jMatCC > 0 ) jMatCC--;
140 
141  const vector<G4double>* deltaCutInKineticEnergy = theCoupleTable->GetEnergyCutsVector(idxG4ElectronCut);
142  const vector<G4double>* photonCutInKineticEnergy = theCoupleTable->GetEnergyCutsVector(idxG4GammaCut);
143  G4double deltaCutInKineticEnergyNow = (*deltaCutInKineticEnergy)[jMatCC];
144  G4double photonCutInKineticEnergyNow = (*photonCutInKineticEnergy)[jMatCC];
145 
146  G4cout<<"G4PAIPhotData::Initialise: "<<"cut = "<<cut/keV<<" keV; cutEl = "
147  <<deltaCutInKineticEnergyNow/keV<<" keV; cutPh = "
148  <<photonCutInKineticEnergyNow/keV<<" keV"<<G4endl;
149 
150  // if( deltaCutInKineticEnergyNow != cut ) deltaCutInKineticEnergyNow = cut; // exception??
151 
152  G4PhysicsLogVector* dEdxCutVector =
153  new G4PhysicsLogVector(fLowestKineticEnergy,
154  fHighestKineticEnergy,
155  fTotBin);
156 
157  G4PhysicsLogVector* dNdxCutVector =
158  new G4PhysicsLogVector(fLowestKineticEnergy,
159  fHighestKineticEnergy,
160  fTotBin);
161  G4PhysicsLogVector* dNdxCutPhotonVector =
162  new G4PhysicsLogVector(fLowestKineticEnergy,
163  fHighestKineticEnergy,
164  fTotBin);
165  G4PhysicsLogVector* dNdxCutPlasmonVector =
166  new G4PhysicsLogVector(fLowestKineticEnergy,
167  fHighestKineticEnergy,
168  fTotBin);
169 
170  const G4Material* mat = couple->GetMaterial();
171  fSandia.Initialize(const_cast<G4Material*>(mat));
172 
173  G4PhysicsTable* PAItransferTable = new G4PhysicsTable(fTotBin+1);
174  G4PhysicsTable* PAIphotonTable = new G4PhysicsTable(fTotBin+1);
175  G4PhysicsTable* PAIplasmonTable = new G4PhysicsTable(fTotBin+1);
176 
177  G4PhysicsTable* PAIdEdxTable = new G4PhysicsTable(fTotBin+1);
178  G4PhysicsLogVector* dEdxMeanVector =
179  new G4PhysicsLogVector(fLowestKineticEnergy,
180  fHighestKineticEnergy,
181  fTotBin);
182 
183  // low energy Sandia interval
184  G4double Tmin = fSandia.GetSandiaMatTablePAI(0,0);
185 
186  // energy safety
187  const G4double deltaLow = 100.*eV;
188 
189  for (G4int i = 0; i <= fTotBin; ++i)
190  {
191  G4double kinEnergy = fParticleEnergyVector->Energy(i);
192  G4double Tmax = model->ComputeMaxEnergy(kinEnergy);
193  G4double tau = kinEnergy/proton_mass_c2;
194  G4double bg2 = tau*( tau + 2. );
195 
196  if ( Tmax < Tmin + deltaLow ) Tmax = Tmin + deltaLow;
197 
198  fPAIxSection.Initialize( mat, Tmax, bg2, &fSandia);
199 
200  //G4cout << i << ". TransferMax(keV)= "<< Tmax/keV << " cut(keV)= "
201  // << cut/keV << " E(MeV)= " << kinEnergy/MeV << G4endl;
202 
203  G4int n = fPAIxSection.GetSplineSize();
204 
205  G4PhysicsFreeVector* transferVector = new G4PhysicsFreeVector(n);
206  G4PhysicsFreeVector* photonVector = new G4PhysicsFreeVector(n);
207  G4PhysicsFreeVector* plasmonVector = new G4PhysicsFreeVector(n);
208 
209  G4PhysicsFreeVector* dEdxVector = new G4PhysicsFreeVector(n);
210 
211  for( G4int k = 0; k < n; k++ )
212  {
213  G4double t = fPAIxSection.GetSplineEnergy(k+1);
214 
215  transferVector->PutValue(k , t,
216  t*fPAIxSection.GetIntegralPAIxSection(k+1));
217  photonVector->PutValue(k , t,
218  t*fPAIxSection.GetIntegralCerenkov(k+1));
219  plasmonVector->PutValue(k , t,
220  t*fPAIxSection.GetIntegralPlasmon(k+1));
221 
222  dEdxVector->PutValue(k, t, fPAIxSection.GetIntegralPAIdEdx(k+1));
223  }
224  // G4cout << *transferVector << G4endl;
225 
226  G4double ionloss = fPAIxSection.GetMeanEnergyLoss();// total <dE/dx>
227 
228  if(ionloss < 0.0) ionloss = 0.0;
229 
230  dEdxMeanVector->PutValue(i,ionloss);
231 
232  G4double dNdxCut = transferVector->Value(deltaCutInKineticEnergyNow)/deltaCutInKineticEnergyNow;
233  G4double dNdxCutPhoton = photonVector->Value(photonCutInKineticEnergyNow)/photonCutInKineticEnergyNow;
234  G4double dNdxCutPlasmon = plasmonVector->Value(deltaCutInKineticEnergyNow)/deltaCutInKineticEnergyNow;
235 
236  G4double dEdxCut = dEdxVector->Value(cut)/cut;
237  //G4cout << "i= " << i << " x= " << dNdxCut << G4endl;
238 
239  if(dNdxCut < 0.0) { dNdxCut = 0.0; }
240  if(dNdxCutPhoton < 0.0) { dNdxCutPhoton = 0.0; }
241  if(dNdxCutPlasmon < 0.0) { dNdxCutPlasmon = 0.0; }
242 
243  dNdxCutVector->PutValue(i, dNdxCut);
244  dNdxCutPhotonVector->PutValue(i, dNdxCutPhoton);
245  dNdxCutPlasmonVector->PutValue(i, dNdxCutPlasmon);
246 
247  dEdxCutVector->PutValue(i, dEdxCut);
248 
249  PAItransferTable->insertAt(i,transferVector);
250  PAIphotonTable->insertAt(i,photonVector);
251  PAIplasmonTable->insertAt(i,plasmonVector);
252  PAIdEdxTable->insertAt(i,dEdxVector);
253 
254  } // end of Tkin loop
255 
256  fPAIxscBank.push_back(PAItransferTable);
257  fPAIphotonBank.push_back(PAIphotonTable);
258  fPAIplasmonBank.push_back(PAIplasmonTable);
259 
260  fPAIdEdxBank.push_back(PAIdEdxTable);
261  fdEdxTable.push_back(dEdxMeanVector);
262 
263  fdNdxCutTable.push_back(dNdxCutVector);
264  fdNdxCutPhotonTable.push_back(dNdxCutPhotonVector);
265  fdNdxCutPlasmonTable.push_back(dNdxCutPlasmonVector);
266 
267  fdEdxCutTable.push_back(dEdxCutVector);
268 }
269 
271 
273  G4double cut) const
274 {
275  // VI: iPlace is the low edge index of the bin
276  // iPlace is in interval from 0 to (N-1)
277  size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
278  size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
279 
280  G4bool one = true;
281  if(scaledTkin >= fParticleEnergyVector->Energy(nPlace)) { iPlace = nPlace; }
282  else if(scaledTkin > fParticleEnergyVector->Energy(0)) {
283  one = false;
284  }
285 
286  // VI: apply interpolation of the vector
287  G4double dEdx = fdEdxTable[coupleIndex]->Value(scaledTkin);
288  G4double del = (*(fPAIdEdxBank[coupleIndex]))(iPlace)->Value(cut);
289  if(!one) {
290  G4double del2 = (*(fPAIdEdxBank[coupleIndex]))(iPlace+1)->Value(cut);
291  G4double E1 = fParticleEnergyVector->Energy(iPlace);
292  G4double E2 = fParticleEnergyVector->Energy(iPlace+1);
293  G4double W = 1.0/(E2 - E1);
294  G4double W1 = (E2 - scaledTkin)*W;
295  G4double W2 = (scaledTkin - E1)*W;
296  del *= W1;
297  del += W2*del2;
298  }
299  dEdx -= del;
300 
301  if( dEdx < 0.) { dEdx = 0.; }
302  return dEdx;
303 }
304 
306 
307 G4double
309  G4double scaledTkin,
310  G4double tcut, G4double tmax) const
311 {
312  G4double cross, xscEl, xscEl2, xscPh, xscPh2;
313 
314  cross=tcut+tmax;
315 
316  // iPlace is in interval from 0 to (N-1)
317 
318  size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
319  size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
320 
321  G4bool one = true;
322 
323  if( scaledTkin >= fParticleEnergyVector->Energy(nPlace)) iPlace = nPlace;
324  else if( scaledTkin > fParticleEnergyVector->Energy(0) ) one = false;
325 
326 
327  xscEl2 = (*fdNdxCutPlasmonTable[coupleIndex])(iPlace);
328  xscPh2 = (*fdNdxCutPhotonTable[coupleIndex])(iPlace);
329 
330  xscPh = xscPh2;
331  xscEl = xscEl2;
332 
333  cross = xscPh + xscEl;
334 
335  if( !one )
336  {
337  xscEl2 = (*fdNdxCutPlasmonTable[coupleIndex])(iPlace+1);
338 
339  G4double E1 = fParticleEnergyVector->Energy(iPlace);
340  G4double E2 = fParticleEnergyVector->Energy(iPlace+1);
341 
342  G4double W = 1.0/(E2 - E1);
343  G4double W1 = (E2 - scaledTkin)*W;
344  G4double W2 = (scaledTkin - E1)*W;
345 
346  xscEl *= W1;
347  xscEl += W2*xscEl2;
348 
349  xscPh2 = (*fdNdxCutPhotonTable[coupleIndex])(iPlace+1);
350 
351  E1 = fParticleEnergyVector->Energy(iPlace);
352  E2 = fParticleEnergyVector->Energy(iPlace+1);
353 
354  W = 1.0/(E2 - E1);
355  W1 = (E2 - scaledTkin)*W;
356  W2 = (scaledTkin - E1)*W;
357 
358  xscPh *= W1;
359  xscPh += W2*xscPh2;
360 
361  cross = xscEl + xscPh;
362  }
363  if( cross < 0.0) cross = 0.0;
364 
365  return cross;
366 }
367 
369 
370 G4double
371 G4PAIPhotData::GetPlasmonRatio(G4int coupleIndex, G4double scaledTkin) const
372 {
373  G4double cross, xscEl, xscEl2, xscPh, xscPh2, plRatio;
374  // iPlace is in interval from 0 to (N-1)
375 
376  size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
377  size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
378 
379  G4bool one = true;
380 
381  if( scaledTkin >= fParticleEnergyVector->Energy(nPlace)) iPlace = nPlace;
382  else if( scaledTkin > fParticleEnergyVector->Energy(0) ) one = false;
383 
384 
385  xscEl2 = (*fdNdxCutPlasmonTable[coupleIndex])(iPlace);
386  xscPh2 = (*fdNdxCutPhotonTable[coupleIndex])(iPlace);
387 
388  xscPh = xscPh2;
389  xscEl = xscEl2;
390 
391  cross = xscPh + xscEl;
392 
393  if( !one )
394  {
395  xscEl2 = (*fdNdxCutPlasmonTable[coupleIndex])(iPlace+1);
396 
397  G4double E1 = fParticleEnergyVector->Energy(iPlace);
398  G4double E2 = fParticleEnergyVector->Energy(iPlace+1);
399 
400  G4double W = 1.0/(E2 - E1);
401  G4double W1 = (E2 - scaledTkin)*W;
402  G4double W2 = (scaledTkin - E1)*W;
403 
404  xscEl *= W1;
405  xscEl += W2*xscEl2;
406 
407  xscPh2 = (*fdNdxCutPhotonTable[coupleIndex])(iPlace+1);
408 
409  E1 = fParticleEnergyVector->Energy(iPlace);
410  E2 = fParticleEnergyVector->Energy(iPlace+1);
411 
412  W = 1.0/(E2 - E1);
413  W1 = (E2 - scaledTkin)*W;
414  W2 = (scaledTkin - E1)*W;
415 
416  xscPh *= W1;
417  xscPh += W2*xscPh2;
418 
419  cross = xscEl + xscPh;
420  }
421  if( cross <= 0.0)
422  {
423  plRatio = 2.0;
424  }
425  else
426  {
427  plRatio = xscEl/cross;
428 
429  if( plRatio > 1. || plRatio < 0.) plRatio = 2.0;
430  }
431  return plRatio;
432 }
433 
435 
437  G4double kinEnergy,
438  G4double scaledTkin,
439  G4double stepFactor) const
440 {
441  G4double loss = 0.0;
442  G4double omega;
443  G4double position, E1, E2, W1, W2, W, dNdxCut1, dNdxCut2, meanNumber;
444 
445  size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
446  size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
447 
448  G4bool one = true;
449 
450  if (scaledTkin >= fParticleEnergyVector->Energy(nPlace)) iPlace = nPlace;
451  else if(scaledTkin > fParticleEnergyVector->Energy(0)) one = false;
452 
453  G4PhysicsLogVector* vcut = fdNdxCutTable[coupleIndex];
454  G4PhysicsVector* v1 = (*(fPAIxscBank[coupleIndex]))(iPlace);
455  G4PhysicsVector* v2 = 0;
456 
457  dNdxCut1 = (*vcut)[iPlace];
458  G4double e1 = v1->Energy(0);
459  G4double e2 = e1;
460 
461  G4double meanN1 = ((*v1)[0]/e1 - dNdxCut1)*stepFactor;
462 
463  meanNumber = meanN1;
464 
465  // G4cout<<"iPlace = "<<iPlace<< " meanN1= " << meanN1
466  // <<" (*v1)[0]/e1 = "<<(*v1)[0]/e1<< " dNdxCut1= " << dNdxCut1 << G4endl;
467 
468  dNdxCut2 = dNdxCut1;
469  W1 = 1.0;
470  W2 = 0.0;
471  if(!one)
472  {
473  v2 = (*(fPAIxscBank[coupleIndex]))(iPlace+1);
474  dNdxCut2 = (*vcut)[iPlace+1];
475  e2 = v2->Energy(0);
476 
477  G4double meanN2 = ((*v2)[0]/e2 - dNdxCut2)*stepFactor;
478 
479  E1 = fParticleEnergyVector->Energy(iPlace);
480  E2 = fParticleEnergyVector->Energy(iPlace+1);
481  W = 1.0/(E2 - E1);
482  W1 = (E2 - scaledTkin)*W;
483  W2 = (scaledTkin - E1)*W;
484  meanNumber = W1*meanN1 + W2*meanN2;
485 
486  //G4cout<<"meanN= " << meanNumber << " meanN2= " << meanN2
487  // << " dNdxCut2= " << dNdxCut2 << G4endl;
488  }
489  if( meanNumber <= 0.0) return 0.0;
490 
491  G4int numOfCollisions = G4Poisson(meanNumber);
492 
493  //G4cout << "N= " << numOfCollisions << G4endl;
494 
495  if( 0 == numOfCollisions) return 0.0;
496 
497  for(G4int i=0; i< numOfCollisions; ++i)
498  {
499  G4double rand = G4UniformRand();
500  position = dNdxCut1 + ((*v1)[0]/e1 - dNdxCut1)*rand;
501  omega = GetEnergyTransfer(coupleIndex, iPlace, position);
502 
503  //G4cout << "omega(keV)= " << omega/keV << G4endl;
504 
505  if(!one)
506  {
507  position = dNdxCut2 + ((*v2)[0]/e2 - dNdxCut2)*rand;
508  G4double omega2 = GetEnergyTransfer(coupleIndex, iPlace+1, position);
509  omega = omega*W1 + omega2*W2;
510  }
511  //G4cout << "omega(keV)= " << omega/keV << G4endl;
512 
513  loss += omega;
514  if( loss > kinEnergy) { break; }
515  }
516 
517  // G4cout<<"PAIPhotData AlongStepLoss = "<<loss/keV<<" keV, on step = "
518  //<<step/mm<<" mm"<<G4endl;
519 
520  if ( loss > kinEnergy) loss = kinEnergy;
521  else if( loss < 0.) loss = 0.;
522 
523  return loss;
524 }
525 
527 
529  G4double kinEnergy,
530  G4double scaledTkin,
531  G4double stepFactor) const
532 {
533  G4double loss = 0.0;
534  G4double omega;
535  G4double position, E1, E2, W1, W2, W, dNdxCut1, dNdxCut2, meanNumber;
536 
537  size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
538  size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
539 
540  G4bool one = true;
541 
542  if (scaledTkin >= fParticleEnergyVector->Energy(nPlace)) iPlace = nPlace;
543  else if(scaledTkin > fParticleEnergyVector->Energy(0)) one = false;
544 
545  G4PhysicsLogVector* vcut = fdNdxCutPhotonTable[coupleIndex];
546  G4PhysicsVector* v1 = (*(fPAIphotonBank[coupleIndex]))(iPlace);
547  G4PhysicsVector* v2 = 0;
548 
549  dNdxCut1 = (*vcut)[iPlace];
550  G4double e1 = v1->Energy(0);
551  G4double e2 = e1;
552 
553  G4double meanN1 = ((*v1)[0]/e1 - dNdxCut1)*stepFactor;
554 
555  meanNumber = meanN1;
556 
557  // G4cout<<"iPlace = "<<iPlace<< " meanN1= " << meanN1
558  // <<" (*v1)[0]/e1 = "<<(*v1)[0]/e1<< " dNdxCut1= " << dNdxCut1 << G4endl;
559 
560  dNdxCut2 = dNdxCut1;
561  W1 = 1.0;
562  W2 = 0.0;
563  if(!one)
564  {
565  v2 = (*(fPAIphotonBank[coupleIndex]))(iPlace+1);
566  dNdxCut2 = (*vcut)[iPlace+1];
567  e2 = v2->Energy(0);
568 
569  G4double meanN2 = ((*v2)[0]/e2 - dNdxCut2)*stepFactor;
570 
571  E1 = fParticleEnergyVector->Energy(iPlace);
572  E2 = fParticleEnergyVector->Energy(iPlace+1);
573  W = 1.0/(E2 - E1);
574  W1 = (E2 - scaledTkin)*W;
575  W2 = (scaledTkin - E1)*W;
576  meanNumber = W1*meanN1 + W2*meanN2;
577 
578  //G4cout<<"meanN= " << meanNumber << " meanN2= " << meanN2
579  // << " dNdxCut2= " << dNdxCut2 << G4endl;
580  }
581  if( meanNumber <= 0.0) return 0.0;
582 
583  G4int numOfCollisions = G4Poisson(meanNumber);
584 
585  //G4cout << "N= " << numOfCollisions << G4endl;
586 
587  if( 0 == numOfCollisions) return 0.0;
588 
589  for(G4int i=0; i< numOfCollisions; ++i)
590  {
591  G4double rand = G4UniformRand();
592  position = dNdxCut1 + ((*v1)[0]/e1 - dNdxCut1)*rand;
593  omega = GetEnergyPhotonTransfer(coupleIndex, iPlace, position);
594 
595  //G4cout << "omega(keV)= " << omega/keV << G4endl;
596 
597  if(!one)
598  {
599  position = dNdxCut2 + ((*v2)[0]/e2 - dNdxCut2)*rand;
600  G4double omega2 = GetEnergyPhotonTransfer(coupleIndex, iPlace+1, position);
601  omega = omega*W1 + omega2*W2;
602  }
603  //G4cout << "omega(keV)= " << omega/keV << G4endl;
604 
605  loss += omega;
606  if( loss > kinEnergy) { break; }
607  }
608 
609  // G4cout<<"PAIPhotData AlongStepLoss = "<<loss/keV<<" keV, on step = "
610  //<<step/mm<<" mm"<<G4endl;
611 
612  if ( loss > kinEnergy) loss = kinEnergy;
613  else if( loss < 0.) loss = 0.;
614 
615  return loss;
616 }
617 
619 
621  G4double kinEnergy,
622  G4double scaledTkin,
623  G4double stepFactor) const
624 {
625  G4double loss = 0.0;
626  G4double omega;
627  G4double position, E1, E2, W1, W2, W, dNdxCut1, dNdxCut2, meanNumber;
628 
629  size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
630  size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
631 
632  G4bool one = true;
633 
634  if (scaledTkin >= fParticleEnergyVector->Energy(nPlace)) iPlace = nPlace;
635  else if(scaledTkin > fParticleEnergyVector->Energy(0)) one = false;
636 
637  G4PhysicsLogVector* vcut = fdNdxCutPlasmonTable[coupleIndex];
638  G4PhysicsVector* v1 = (*(fPAIplasmonBank[coupleIndex]))(iPlace);
639  G4PhysicsVector* v2 = 0;
640 
641  dNdxCut1 = (*vcut)[iPlace];
642  G4double e1 = v1->Energy(0);
643  G4double e2 = e1;
644 
645  G4double meanN1 = ((*v1)[0]/e1 - dNdxCut1)*stepFactor;
646 
647  meanNumber = meanN1;
648 
649  // G4cout<<"iPlace = "<<iPlace<< " meanN1= " << meanN1
650  // <<" (*v1)[0]/e1 = "<<(*v1)[0]/e1<< " dNdxCut1= " << dNdxCut1 << G4endl;
651 
652  dNdxCut2 = dNdxCut1;
653  W1 = 1.0;
654  W2 = 0.0;
655  if(!one)
656  {
657  v2 = (*(fPAIplasmonBank[coupleIndex]))(iPlace+1);
658  dNdxCut2 = (*vcut)[iPlace+1];
659  e2 = v2->Energy(0);
660 
661  G4double meanN2 = ((*v2)[0]/e2 - dNdxCut2)*stepFactor;
662 
663  E1 = fParticleEnergyVector->Energy(iPlace);
664  E2 = fParticleEnergyVector->Energy(iPlace+1);
665  W = 1.0/(E2 - E1);
666  W1 = (E2 - scaledTkin)*W;
667  W2 = (scaledTkin - E1)*W;
668  meanNumber = W1*meanN1 + W2*meanN2;
669 
670  //G4cout<<"meanN= " << meanNumber << " meanN2= " << meanN2
671  // << " dNdxCut2= " << dNdxCut2 << G4endl;
672  }
673  if( meanNumber <= 0.0) return 0.0;
674 
675  G4int numOfCollisions = G4Poisson(meanNumber);
676 
677  //G4cout << "N= " << numOfCollisions << G4endl;
678 
679  if( 0 == numOfCollisions) return 0.0;
680 
681  for(G4int i=0; i< numOfCollisions; ++i)
682  {
683  G4double rand = G4UniformRand();
684  position = dNdxCut1 + ((*v1)[0]/e1 - dNdxCut1)*rand;
685  omega = GetEnergyPlasmonTransfer(coupleIndex, iPlace, position);
686 
687  //G4cout << "omega(keV)= " << omega/keV << G4endl;
688 
689  if(!one)
690  {
691  position = dNdxCut2 + ((*v2)[0]/e2 - dNdxCut2)*rand;
692  G4double omega2 = GetEnergyPlasmonTransfer(coupleIndex, iPlace+1, position);
693  omega = omega*W1 + omega2*W2;
694  }
695  //G4cout << "omega(keV)= " << omega/keV << G4endl;
696 
697  loss += omega;
698  if( loss > kinEnergy) { break; }
699  }
700 
701  // G4cout<<"PAIPhotData AlongStepLoss = "<<loss/keV<<" keV, on step = "
702  //<<step/mm<<" mm"<<G4endl;
703 
704  if ( loss > kinEnergy) loss = kinEnergy;
705  else if( loss < 0.) loss = 0.;
706 
707  return loss;
708 }
709 
711 //
712 // Returns post step PAI energy transfer > cut electron energy
713 // according to passed scaled kinetic energy of particle
714 
716  G4double scaledTkin) const
717 {
718  //G4cout<<"G4PAIPhotData::GetPostStepTransfer"<<G4endl;
719  G4double transfer = 0.0;
720  G4double rand = G4UniformRand();
721 
722  size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
723 
724  // size_t iTransfer, iTr1, iTr2;
725  G4double position, dNdxCut1, dNdxCut2, E1, E2, W1, W2, W;
726 
727  G4PhysicsVector* cutv = fdNdxCutTable[coupleIndex];
728 
729  // Fermi plato, try from left
730  if( scaledTkin >= fParticleEnergyVector->GetMaxEnergy())
731  {
732  position = (*cutv)[nPlace]*rand;
733  transfer = GetEnergyTransfer(coupleIndex, nPlace, position);
734  }
735  else if( scaledTkin <= fParticleEnergyVector->Energy(0) )
736  {
737  position = (*cutv)[0]*rand;
738  transfer = GetEnergyTransfer(coupleIndex, 0, position);
739  }
740  else
741  {
742  size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
743 
744  dNdxCut1 = (*cutv)[iPlace];
745  dNdxCut2 = (*cutv)[iPlace+1];
746 
747  E1 = fParticleEnergyVector->Energy(iPlace);
748  E2 = fParticleEnergyVector->Energy(iPlace+1);
749  W = 1.0/(E2 - E1);
750  W1 = (E2 - scaledTkin)*W;
751  W2 = (scaledTkin - E1)*W;
752 
753  //G4cout<<"iPlace= " << " dNdxCut1 = "<<dNdxCut1
754  // <<" dNdxCut2 = "<<dNdxCut2<< " W1= " << W1 << " W2= " << W2 <<G4endl;
755 
756  position = dNdxCut1*rand;
757  G4double tr1 = GetEnergyTransfer(coupleIndex, iPlace, position);
758 
759  position = dNdxCut2*rand;
760  G4double tr2 = GetEnergyTransfer(coupleIndex, iPlace+1, position);
761 
762  transfer = tr1*W1 + tr2*W2;
763  }
764  //G4cout<<"PAImodel PostStepTransfer = "<<transfer/keV<<" keV"<<G4endl;
765  if(transfer < 0.0 ) { transfer = 0.0; }
766  return transfer;
767 }
768 
770 
772  G4double scaledTkin) const
773 {
774  //G4cout<<"G4PAIPhotData::GetPostStepTransfer"<<G4endl;
775  G4double transfer = 0.0;
776  G4double rand = G4UniformRand();
777 
778  size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
779 
780  // size_t iTransfer, iTr1, iTr2;
781  G4double position, dNdxCut1, dNdxCut2, E1, E2, W1, W2, W;
782 
783  G4PhysicsVector* cutv = fdNdxCutPhotonTable[coupleIndex];
784 
785  // Fermi plato, try from left
786 
787  if( scaledTkin >= fParticleEnergyVector->GetMaxEnergy())
788  {
789  position = (*cutv)[nPlace]*rand;
790  transfer = GetEnergyPhotonTransfer(coupleIndex, nPlace, position);
791  }
792  else if( scaledTkin <= fParticleEnergyVector->Energy(0) )
793  {
794  position = (*cutv)[0]*rand;
795  transfer = GetEnergyPhotonTransfer(coupleIndex, 0, position);
796  }
797  else
798  {
799  size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
800 
801  dNdxCut1 = (*cutv)[iPlace];
802  dNdxCut2 = (*cutv)[iPlace+1];
803 
804  E1 = fParticleEnergyVector->Energy(iPlace);
805  E2 = fParticleEnergyVector->Energy(iPlace+1);
806  W = 1.0/(E2 - E1);
807  W1 = (E2 - scaledTkin)*W;
808  W2 = (scaledTkin - E1)*W;
809 
810  //G4cout<<"iPlace= " << " dNdxCut1 = "<<dNdxCut1
811  // <<" dNdxCut2 = "<<dNdxCut2<< " W1= " << W1 << " W2= " << W2 <<G4endl;
812 
813  position = dNdxCut1*rand;
814 
815  G4double tr1 = GetEnergyPhotonTransfer(coupleIndex, iPlace, position);
816 
817  position = dNdxCut2*rand;
818  G4double tr2 = GetEnergyPhotonTransfer(coupleIndex, iPlace+1, position);
819 
820  transfer = tr1*W1 + tr2*W2;
821  }
822  //G4cout<<"PAImodel PostStepTransfer = "<<transfer/keV<<" keV"<<G4endl;
823  if(transfer < 0.0 ) { transfer = 0.0; }
824  return transfer;
825 }
826 
828 
830  G4double scaledTkin) const
831 {
832  //G4cout<<"G4PAIPhotData::GetPostStepTransfer"<<G4endl;
833  G4double transfer = 0.0;
834  G4double rand = G4UniformRand();
835 
836  size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
837 
838  // size_t iTransfer, iTr1, iTr2;
839  G4double position, dNdxCut1, dNdxCut2, E1, E2, W1, W2, W;
840 
841  G4PhysicsVector* cutv = fdNdxCutPlasmonTable[coupleIndex];
842 
843  // Fermi plato, try from left
844  if( scaledTkin >= fParticleEnergyVector->GetMaxEnergy())
845  {
846  position = (*cutv)[nPlace]*rand;
847  transfer = GetEnergyPlasmonTransfer(coupleIndex, nPlace, position);
848  }
849  else if( scaledTkin <= fParticleEnergyVector->Energy(0) )
850  {
851  position = (*cutv)[0]*rand;
852  transfer = GetEnergyPlasmonTransfer(coupleIndex, 0, position);
853  }
854  else
855  {
856  size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
857 
858  dNdxCut1 = (*cutv)[iPlace];
859  dNdxCut2 = (*cutv)[iPlace+1];
860 
861  E1 = fParticleEnergyVector->Energy(iPlace);
862  E2 = fParticleEnergyVector->Energy(iPlace+1);
863  W = 1.0/(E2 - E1);
864  W1 = (E2 - scaledTkin)*W;
865  W2 = (scaledTkin - E1)*W;
866 
867  //G4cout<<"iPlace= " << " dNdxCut1 = "<<dNdxCut1
868  // <<" dNdxCut2 = "<<dNdxCut2<< " W1= " << W1 << " W2= " << W2 <<G4endl;
869 
870  position = dNdxCut1*rand;
871  G4double tr1 = GetEnergyPlasmonTransfer(coupleIndex, iPlace, position);
872 
873  position = dNdxCut2*rand;
874  G4double tr2 = GetEnergyPlasmonTransfer(coupleIndex, iPlace+1, position);
875 
876  transfer = tr1*W1 + tr2*W2;
877  }
878  //G4cout<<"PAImodel PostStepPlasmonTransfer = "<<transfer/keV<<" keV"<<G4endl;
879 
880  if(transfer < 0.0 ) transfer = 0.0;
881 
882  return transfer;
883 }
884 
886 //
887 // Returns PAI energy transfer according to passed
888 // indexes of particle kinetic enegry and random x-section
889 
891  size_t iPlace,
892  G4double position) const
893 {
894  G4PhysicsVector* v = (*(fPAIxscBank[coupleIndex]))(iPlace);
895  if(position*v->Energy(0) >= (*v)[0]) { return v->Energy(0); }
896 
897  size_t iTransferMax = v->GetVectorLength() - 1;
898 
899  size_t iTransfer;
900  G4double x1(0.0), x2(0.0), y1(0.0), y2(0.0), energyTransfer;
901 
902  for(iTransfer=1; iTransfer<=iTransferMax; ++iTransfer) {
903  x2 = v->Energy(iTransfer);
904  y2 = (*v)[iTransfer]/x2;
905  if(position >= y2) { break; }
906  }
907 
908  x1 = v->Energy(iTransfer-1);
909  y1 = (*v)[iTransfer-1]/x1;
910  //G4cout << "i= " << iTransfer << " imax= " << iTransferMax
911  // << " x1= " << x1 << " x2= " << x2 << G4endl;
912 
913  energyTransfer = x1;
914  if ( x1 != x2 ) {
915  if ( y1 == y2 ) {
916  energyTransfer += (x2 - x1)*G4UniformRand();
917  } else {
918  if(x1*1.1 < x2) {
919  const G4int nbins = 5;
920  G4double del = (x2 - x1)/G4int(nbins);
921  x2 = x1;
922  for(G4int i=1; i<=nbins; ++i) {
923  x2 += del;
924  y2 = v->Value(x2)/x2;
925  if(position >= y2) { break; }
926  x1 = x2;
927  y1 = y2;
928  }
929  }
930  energyTransfer = (y2 - y1)*x1*x2/(position*(x1 - x2) - y1*x1 + y2*x2);
931  }
932  }
933  // G4cout << "x1(keV)= " << x1/keV << " x2(keV)= " << x2/keV
934  // << " y1= " << y1 << " y2= " << y2 << " pos= " << position
935  // << " E(keV)= " << energyTransfer/keV << G4endl;
936  return energyTransfer;
937 }
938 
940 
942  size_t iPlace,
943  G4double position) const
944 {
945  G4PhysicsVector* v = (*(fPAIphotonBank[coupleIndex]))(iPlace);
946  if(position*v->Energy(0) >= (*v)[0]) return v->Energy(0);
947 
948  size_t iTransferMax = v->GetVectorLength() - 1;
949 
950  size_t iTransfer;
951  G4double x1(0.0), x2(0.0), y1(0.0), y2(0.0), energyTransfer;
952 
953  for(iTransfer=1; iTransfer<=iTransferMax; ++iTransfer)
954  {
955  x2 = v->Energy(iTransfer);
956  y2 = (*v)[iTransfer]/x2;
957  if(position >= y2) break;
958  }
959  x1 = v->Energy(iTransfer-1);
960  y1 = (*v)[iTransfer-1]/x1;
961 
962  //G4cout << "i= " << iTransfer << " imax= " << iTransferMax
963  // << " x1= " << x1 << " x2= " << x2 << G4endl;
964 
965  energyTransfer = x1;
966 
967  if ( x1 != x2 )
968  {
969  if ( y1 == y2 )
970  {
971  energyTransfer += (x2 - x1)*G4UniformRand();
972  }
973  else
974  {
975  if( x1*1.1 < x2 )
976  {
977  const G4int nbins = 5;
978  G4double del = (x2 - x1)/G4int(nbins);
979  x2 = x1;
980 
981  for(G4int i=1; i<=nbins; ++i)
982  {
983  x2 += del;
984  y2 = v->Value(x2)/x2;
985  if(position >= y2) { break; }
986  x1 = x2;
987  y1 = y2;
988  }
989  }
990  energyTransfer = (y2 - y1)*x1*x2/(position*(x1 - x2) - y1*x1 + y2*x2);
991  }
992  }
993  // G4cout << "x1(keV)= " << x1/keV << " x2(keV)= " << x2/keV
994  // << " y1= " << y1 << " y2= " << y2 << " pos= " << position
995  // << " E(keV)= " << energyTransfer/keV << G4endl;
996 
997  return energyTransfer;
998 }
999 
1001 
1003  size_t iPlace,
1004  G4double position) const
1005 {
1006  G4PhysicsVector* v = (*(fPAIplasmonBank[coupleIndex]))(iPlace);
1007 
1008  if( position*v->Energy(0) >= (*v)[0] ) return v->Energy(0);
1009 
1010  size_t iTransferMax = v->GetVectorLength() - 1;
1011 
1012  size_t iTransfer;
1013  G4double x1(0.0), x2(0.0), y1(0.0), y2(0.0), energyTransfer;
1014 
1015  for(iTransfer = 1; iTransfer <= iTransferMax; ++iTransfer)
1016  {
1017  x2 = v->Energy(iTransfer);
1018  y2 = (*v)[iTransfer]/x2;
1019  if(position >= y2) break;
1020  }
1021  x1 = v->Energy(iTransfer-1);
1022  y1 = (*v)[iTransfer-1]/x1;
1023 
1024  //G4cout << "i= " << iTransfer << " imax= " << iTransferMax
1025  // << " x1= " << x1 << " x2= " << x2 << G4endl;
1026 
1027  energyTransfer = x1;
1028 
1029  if ( x1 != x2 )
1030  {
1031  if ( y1 == y2 )
1032  {
1033  energyTransfer += (x2 - x1)*G4UniformRand();
1034  }
1035  else
1036  {
1037  if(x1*1.1 < x2)
1038  {
1039  const G4int nbins = 5;
1040  G4double del = (x2 - x1)/G4int(nbins);
1041  x2 = x1;
1042 
1043  for( G4int i = 1; i <= nbins; ++i )
1044  {
1045  x2 += del;
1046  y2 = v->Value(x2)/x2;
1047 
1048  if(position >= y2) break;
1049 
1050  x1 = x2;
1051  y1 = y2;
1052  }
1053  }
1054  energyTransfer = (y2 - y1)*x1*x2/(position*(x1 - x2) - y1*x1 + y2*x2);
1055  }
1056  }
1057  // G4cout << "x1(keV)= " << x1/keV << " x2(keV)= " << x2/keV
1058  // << " y1= " << y1 << " y2= " << y2 << " pos= " << position
1059  // << " E(keV)= " << energyTransfer/keV << G4endl;
1060 
1061  return energyTransfer;
1062 }
1063 
1064 //
1065 //
1067 
G4double Energy(size_t index) const
G4double SamplePostStepPlasmonTransfer(G4int coupleIndex, G4double scaledTkin) const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4double GetEnergyPlasmonTransfer(G4int coupleIndex, size_t iPlace, G4double position) const
static constexpr double MeV
Definition: G4SIunits.hh:214
Float_t y1[n_points_granero]
Definition: compare.C:5
static constexpr double keV
Definition: G4SIunits.hh:216
G4double SamplePostStepTransfer(G4int coupleIndex, G4double scaledTkin) const
Float_t x1[n_points_granero]
Definition: compare.C:5
#define G4endl
Definition: G4ios.hh:61
void insertAt(size_t, G4PhysicsVector *)
G4double Value(G4double theEnergy, size_t &lastidx) const
G4double GetPlasmonRatio(G4int coupleIndex, G4double scaledTkin) const
const std::vector< G4double > * GetEnergyCutsVector(size_t pcIdx) const
Float_t y2[n_points_geant4]
Definition: compare.C:26
static constexpr double proton_mass_c2
static constexpr double TeV
Definition: G4SIunits.hh:218
double G4double
Definition: G4Types.hh:76
bool G4bool
Definition: G4Types.hh:79
#define position
Definition: xmlparse.cc:622
#define G4UniformRand()
Definition: Randomize.hh:53
G4PAIPhotData(G4double tmin, G4double tmax, G4int verbose)
static constexpr double eV
Definition: G4SIunits.hh:215
Float_t mat
G4double SampleAlongStepTransfer(G4int coupleIndex, G4double kinEnergy, G4double scaledTkin, G4double stepFactor) const
void Initialise(const G4MaterialCutsCouple *, G4double cut, G4PAIPhotModel *)
static G4ProductionCutsTable * GetProductionCutsTable()
G4double SamplePostStepPhotonTransfer(G4int coupleIndex, G4double scaledTkin) const
G4double SampleAlongStepPlasmonTransfer(G4int coupleIndex, G4double kinEnergy, G4double scaledTkin, G4double stepFactor) const
int G4int
Definition: G4Types.hh:78
G4double SampleAlongStepPhotonTransfer(G4int coupleIndex, G4double kinEnergy, G4double scaledTkin, G4double stepFactor) const
G4double CrossSectionPerVolume(G4int coupleIndex, G4double scaledTkin, G4double tcut, G4double tmax) const
G4GLOB_DLL std::ostream G4cout
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
G4long G4Poisson(G4double mean)
Definition: G4Poisson.hh:51
Char_t n[5]
G4double ComputeMaxEnergy(G4double scaledEnergy)
const G4Material * GetMaterial() const
void PutValue(size_t index, G4double energy, G4double dataValue)
Float_t x2[n_points_geant4]
Definition: compare.C:26
const XML_Char XML_Content * model
Definition: expat.h:151
static constexpr double GeV
Definition: G4SIunits.hh:217
void PutValue(size_t index, G4double theValue)
G4double DEDXPerVolume(G4int coupleIndex, G4double scaledTkin, G4double cut) const
G4double GetEnergyPhotonTransfer(G4int coupleIndex, size_t iPlace, G4double position) const
size_t GetVectorLength() const
G4double GetEnergyTransfer(G4int coupleIndex, size_t iPlace, G4double position) const