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