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