Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4PAIPhotData.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // $Id: 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 
890 G4double G4PAIPhotData::GetEnergyTransfer(G4int coupleIndex,
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 
941 G4double G4PAIPhotData::GetEnergyPhotonTransfer(G4int coupleIndex,
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 
1002 G4double G4PAIPhotData::GetEnergyPlasmonTransfer(G4int coupleIndex,
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 SamplePostStepTransfer(G4int coupleIndex, G4double scaledTkin) const
G4double SamplePostStepPhotonTransfer(G4int coupleIndex, G4double scaledTkin) const
G4long G4Poisson(G4double mean)
Definition: G4Poisson.hh:51
const std::vector< G4double > * GetEnergyCutsVector(size_t pcIdx) const
G4double GetPlasmonRatio(G4int coupleIndex, G4double scaledTkin) const
void PutValue(size_t index, G4double energy, G4double dataValue)
G4double SampleAlongStepPhotonTransfer(G4int coupleIndex, G4double kinEnergy, G4double scaledTkin, G4double stepFactor) const
G4PAIPhotData(G4double tmin, G4double tmax, G4int verbose)
size_t GetVectorLength() const
int G4int
Definition: G4Types.hh:78
void Initialise(const G4MaterialCutsCouple *, G4double cut, G4PAIPhotModel *)
static constexpr double TeV
Definition: G4SIunits.hh:218
#define position
Definition: xmlparse.cc:622
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
G4double ComputeMaxEnergy(G4double scaledEnergy)
bool G4bool
Definition: G4Types.hh:79
G4double DEDXPerVolume(G4int coupleIndex, G4double scaledTkin, G4double cut) const
void PutValue(size_t index, G4double theValue)
static constexpr double eV
Definition: G4SIunits.hh:215
float proton_mass_c2
Definition: hepunit.py:275
const G4int n
G4double Energy(size_t index) const
G4double Value(G4double theEnergy, size_t &lastidx) const
tuple v
Definition: test.py:18
static G4ProductionCutsTable * GetProductionCutsTable()
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4double SamplePostStepPlasmonTransfer(G4int coupleIndex, G4double scaledTkin) const
G4double CrossSectionPerVolume(G4int coupleIndex, G4double scaledTkin, G4double tcut, G4double tmax) const
G4double SampleAlongStepPlasmonTransfer(G4int coupleIndex, G4double kinEnergy, G4double scaledTkin, G4double stepFactor) const
static constexpr double GeV
Definition: G4SIunits.hh:217
void insertAt(size_t, G4PhysicsVector *)
#define G4endl
Definition: G4ios.hh:61
static constexpr double MeV
Definition: G4SIunits.hh:214
double G4double
Definition: G4Types.hh:76
G4double SampleAlongStepTransfer(G4int coupleIndex, G4double kinEnergy, G4double scaledTkin, G4double stepFactor) const
static constexpr double keV
Definition: G4SIunits.hh:216
const XML_Char XML_Content * model
Definition: expat.h:151
const G4Material * GetMaterial() const