Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4PAIModel.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$
27 //
28 // -------------------------------------------------------------------
29 //
30 // GEANT4 Class
31 // File name: G4PAIModel.cc
32 //
33 // Author: Vladimir.Grichine@cern.ch on base of Vladimir Ivanchenko model interface
34 //
35 // Creation date: 05.10.2003
36 //
37 // Modifications:
38 //
39 // 17.08.04 V.Grichine, bug fixed for Tkin<=0 in SampleSecondary
40 // 16.08.04 V.Grichine, bug fixed in massRatio for DEDX, CrossSection, SampleSecondary
41 // 08.04.05 Major optimisation of internal interfaces (V.Ivantchenko)
42 // 26.07.09 Fixed logic to work with several materials (V.Ivantchenko)
43 // 21.11.10 V. Grichine verbose flag for protons and G4PAYySection to check sandia table
44 //
45 
46 #include "G4PAIModel.hh"
47 
48 #include "G4PhysicalConstants.hh"
49 #include "G4SystemOfUnits.hh"
50 #include "G4Region.hh"
51 #include "G4PhysicsLogVector.hh"
52 #include "G4PhysicsFreeVector.hh"
53 #include "G4PhysicsTable.hh"
54 #include "G4ProductionCutsTable.hh"
55 #include "G4MaterialCutsCouple.hh"
56 #include "G4MaterialTable.hh"
57 #include "G4SandiaTable.hh"
58 #include "G4OrderedTable.hh"
59 
60 #include "Randomize.hh"
61 #include "G4Electron.hh"
62 #include "G4Positron.hh"
63 #include "G4Poisson.hh"
64 #include "G4Step.hh"
65 #include "G4Material.hh"
66 #include "G4DynamicParticle.hh"
67 #include "G4ParticleDefinition.hh"
68 #include "G4ParticleChangeForLoss.hh"
69 
71 
72 using namespace std;
73 
74 G4PAIModel::G4PAIModel(const G4ParticleDefinition* p, const G4String& nam)
75  : G4VEmModel(nam),G4VEmFluctuationModel(nam),
76  fVerbose(0),
77  fLowestGamma(1.005),
78  fHighestGamma(10000.),
79  fTotBin(200),
80  fMeanNumber(20),
81  fParticle(0),
82  fHighKinEnergy(100.*TeV),
83  fTwoln10(2.0*log(10.0)),
84  fBg2lim(0.0169),
85  fTaulim(8.4146e-3)
86 {
87  fElectron = G4Electron::Electron();
88  fPositron = G4Positron::Positron();
89 
90  fPAItransferTable = 0;
91  fPAIdEdxTable = 0;
92  fSandiaPhotoAbsCof = 0;
93  fdEdxVector = 0;
94  fLambdaVector = 0;
95  fdNdxCutVector = 0;
96  fParticleEnergyVector = 0;
97  fSandiaIntervalNumber = 0;
98  fMatIndex = 0;
99  fDeltaCutInKinEnergy = 0.0;
100  fParticleChange = 0;
101  fMaterial = 0;
102  fCutCouple = 0;
103 
104  if(p) { SetParticle(p); }
105  else { SetParticle(fElectron); }
106 
107  isInitialised = false;
108 }
109 
111 
112 G4PAIModel::~G4PAIModel()
113 {
114  // G4cout << "PAI: start destruction" << G4endl;
115  if(fParticleEnergyVector) delete fParticleEnergyVector;
116  if(fdEdxVector) delete fdEdxVector ;
117  if(fLambdaVector) delete fLambdaVector;
118  if(fdNdxCutVector) delete fdNdxCutVector;
119 
120  if( fPAItransferTable )
121  {
122  fPAItransferTable->clearAndDestroy();
123  delete fPAItransferTable ;
124  }
125  if( fPAIdEdxTable )
126  {
127  fPAIdEdxTable->clearAndDestroy();
128  delete fPAIdEdxTable ;
129  }
130  if(fSandiaPhotoAbsCof)
131  {
132  for(G4int i=0;i<fSandiaIntervalNumber;i++)
133  {
134  delete[] fSandiaPhotoAbsCof[i];
135  }
136  delete[] fSandiaPhotoAbsCof;
137  }
138  //G4cout << "PAI: end destruction" << G4endl;
139 }
140 
142 
143 void G4PAIModel::SetParticle(const G4ParticleDefinition* p)
144 {
145  if(fParticle == p) { return; }
146  fParticle = p;
147  fMass = fParticle->GetPDGMass();
148  fSpin = fParticle->GetPDGSpin();
149  G4double q = fParticle->GetPDGCharge()/eplus;
150  fChargeSquare = q*q;
151  fLowKinEnergy = 0.2*MeV*fMass/proton_mass_c2;
152  fRatio = electron_mass_c2/fMass;
153  fQc = fMass/fRatio;
154  fLowestKineticEnergy = fMass*(fLowestGamma - 1.0);
155  fHighestKineticEnergy = fMass*(fHighestGamma - 1.0);
156 }
157 
159 
160 void G4PAIModel::Initialise(const G4ParticleDefinition* p,
161  const G4DataVector&)
162 {
163  if( fVerbose > 0 && p->GetParticleName()=="proton")
164  {
165  G4cout<<"G4PAIModel::Initialise for "<<p->GetParticleName()<<G4endl;
166  fPAIySection.SetVerbose(1);
167  }
168  else fPAIySection.SetVerbose(0);
169 
170  if(isInitialised) { return; }
171  isInitialised = true;
172 
173  SetParticle(p);
174 
175  fParticleEnergyVector = new G4PhysicsLogVector(fLowestKineticEnergy,
176  fHighestKineticEnergy,
177  fTotBin);
178 
179  fParticleChange = GetParticleChangeForLoss();
180 
181  // Prepare initialization
182  const G4ProductionCutsTable* theCoupleTable =
183  G4ProductionCutsTable::GetProductionCutsTable();
184  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
185  size_t numOfMat = G4Material::GetNumberOfMaterials();
186  size_t numRegions = fPAIRegionVector.size();
187 
188  for(size_t iReg = 0; iReg < numRegions; ++iReg) // region loop
189  {
190  const G4Region* curReg = fPAIRegionVector[iReg];
191 
192  for(size_t jMat = 0; jMat < numOfMat; ++jMat) // region material loop
193  {
194  fMaterial = (*theMaterialTable)[jMat];
195  fCutCouple = theCoupleTable->GetMaterialCutsCouple( fMaterial,
196  curReg->GetProductionCuts() );
197  //G4cout << "Reg <" <<curReg->GetName() << "> mat <"
198  // << fMaterial->GetName() << "> fCouple= "
199  // << fCutCouple<<" " << p->GetParticleName() <<G4endl;
200  if( fCutCouple )
201  {
202  fMaterialCutsCoupleVector.push_back(fCutCouple);
203 
204  fPAItransferTable = new G4PhysicsTable(fTotBin+1);
205  fPAIdEdxTable = new G4PhysicsTable(fTotBin+1);
206 
207  fDeltaCutInKinEnergy =
208  (*theCoupleTable->GetEnergyCutsVector(1))[fCutCouple->GetIndex()];
209 
210  //ComputeSandiaPhotoAbsCof();
211  BuildPAIonisationTable();
212 
213  fPAIxscBank.push_back(fPAItransferTable);
214  fPAIdEdxBank.push_back(fPAIdEdxTable);
215  fdEdxTable.push_back(fdEdxVector);
216 
217  BuildLambdaVector();
218  fdNdxCutTable.push_back(fdNdxCutVector);
219  fLambdaTable.push_back(fLambdaVector);
220  }
221  }
222  }
223 }
224 
226 
227 void G4PAIModel::InitialiseMe(const G4ParticleDefinition*)
228 {}
229 
231 
232 void G4PAIModel::ComputeSandiaPhotoAbsCof()
233 {
234  G4int i, j;
235  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
236 
237  G4SandiaTable thisMaterialSandiaTable(fMatIndex) ;
238  G4int numberOfElements =
239  (*theMaterialTable)[fMatIndex]->GetNumberOfElements();
240 
241  G4int* thisMaterialZ = new G4int[numberOfElements] ;
242 
243  for(i=0;i<numberOfElements;i++)
244  {
245  thisMaterialZ[i] =
246  (G4int)(*theMaterialTable)[fMatIndex]->GetElement(i)->GetZ() ;
247  }
248  fSandiaIntervalNumber = thisMaterialSandiaTable.SandiaIntervals
249  (thisMaterialZ,numberOfElements) ;
250 
251  fSandiaIntervalNumber = thisMaterialSandiaTable.SandiaMixing
252  ( thisMaterialZ ,
253  (*theMaterialTable)[fMatIndex]->GetFractionVector() ,
254  numberOfElements,fSandiaIntervalNumber) ;
255 
256  fSandiaPhotoAbsCof = new G4double*[fSandiaIntervalNumber] ;
257 
258  for(i=0; i<fSandiaIntervalNumber; i++)
259  {
260  fSandiaPhotoAbsCof[i] = new G4double[5];
261  }
262 
263  for( i = 0 ; i < fSandiaIntervalNumber ; i++ )
264  {
265  fSandiaPhotoAbsCof[i][0] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i+1,0);
266 
267  for( j = 1; j < 5 ; j++ )
268  {
269  fSandiaPhotoAbsCof[i][j] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i+1,j)*
270  (*theMaterialTable)[fMatIndex]->GetDensity() ;
271  }
272  }
273  delete[] thisMaterialZ;
274 }
275 
277 //
278 // Build tables for the ionization energy loss
279 // the tables are built for MATERIALS
280 // *********
281 
282 void G4PAIModel::BuildPAIonisationTable()
283 {
284  G4double LowEdgeEnergy , ionloss ;
285  G4double tau, Tmax, Tmin, Tkin, deltaLow, /*gamma,*/ bg2 ;
286 
287  fdEdxVector = new G4PhysicsLogVector( fLowestKineticEnergy,
288  fHighestKineticEnergy,
289  fTotBin);
290  G4SandiaTable* sandia = fMaterial->GetSandiaTable();
291 
292  Tmin = sandia->GetSandiaCofForMaterialPAI(0,0)*keV;
293  deltaLow = 100.*eV; // 0.5*eV ;
294 
295  for (G4int i = 0 ; i <= fTotBin ; i++) //The loop for the kinetic energy
296  {
297  LowEdgeEnergy = fParticleEnergyVector->GetLowEdgeEnergy(i) ;
298  tau = LowEdgeEnergy/fMass ;
299  //gamma = tau +1. ;
300  // G4cout<<"gamma = "<<gamma<<endl ;
301  bg2 = tau*( tau + 2. );
302  Tmax = MaxSecondaryEnergy(fParticle, LowEdgeEnergy);
303  // Tmax = std::min(fDeltaCutInKinEnergy, Tmax);
304  Tkin = Tmax ;
305 
306  if ( Tmax < Tmin + deltaLow ) // low energy safety
307  Tkin = Tmin + deltaLow ;
308 
309  fPAIySection.Initialize(fMaterial, Tkin, bg2);
310 
311  // G4cout<<"ionloss = "<<ionloss*cm/keV<<" keV/cm"<<endl ;
312  // G4cout<<"n1 = "<<protonPAI.GetIntegralPAIxSection(1)*cm<<" 1/cm"<<endl ;
313  // G4cout<<"protonPAI.GetSplineSize() = "<<
314  // protonPAI.GetSplineSize()<<G4endl<<G4endl ;
315 
316  G4int n = fPAIySection.GetSplineSize();
317  G4PhysicsFreeVector* transferVector = new G4PhysicsFreeVector(n) ;
318  G4PhysicsFreeVector* dEdxVector = new G4PhysicsFreeVector(n);
319 
320  for( G4int k = 0 ; k < n; k++ )
321  {
322  transferVector->PutValue( k ,
323  fPAIySection.GetSplineEnergy(k+1),
324  fPAIySection.GetIntegralPAIySection(k+1) ) ;
325  dEdxVector->PutValue( k ,
326  fPAIySection.GetSplineEnergy(k+1),
327  fPAIySection.GetIntegralPAIdEdx(k+1) ) ;
328  }
329  ionloss = fPAIySection.GetMeanEnergyLoss() ; // total <dE/dx>
330 
331  if ( ionloss < DBL_MIN) { ionloss = 0.0; }
332  fdEdxVector->PutValue(i,ionloss) ;
333 
334  fPAItransferTable->insertAt(i,transferVector) ;
335  fPAIdEdxTable->insertAt(i,dEdxVector) ;
336 
337  } // end of Tkin loop
338 }
339 
341 //
342 // Build mean free path tables for the delta ray production process
343 // tables are built for MATERIALS
344 //
345 
346 void G4PAIModel::BuildLambdaVector()
347 {
348  fLambdaVector = new G4PhysicsLogVector( fLowestKineticEnergy,
349  fHighestKineticEnergy,
350  fTotBin ) ;
351  fdNdxCutVector = new G4PhysicsLogVector( fLowestKineticEnergy,
352  fHighestKineticEnergy,
353  fTotBin ) ;
354  if(fVerbose > 1)
355  {
356  G4cout<<"PAIModel DeltaCutInKineticEnergyNow = "
357  <<fDeltaCutInKinEnergy/keV<<" keV"<<G4endl;
358  }
359  for (G4int i = 0 ; i <= fTotBin ; i++ )
360  {
361  G4double dNdxCut = GetdNdxCut(i,fDeltaCutInKinEnergy) ;
362  //G4cout << "i= " << i << " x= " << dNdxCut << G4endl;
363  if(dNdxCut < 0.0) dNdxCut = 0.0;
364  // G4double lambda = dNdxCut <= DBL_MIN ? DBL_MAX: 1.0/dNdxCut ;
365  G4double lambda = DBL_MAX;
366  if(dNdxCut > 0.0) lambda = 1.0/dNdxCut;
367 
368  fLambdaVector->PutValue(i, lambda) ;
369  fdNdxCutVector->PutValue(i, dNdxCut) ;
370  }
371 }
372 
374 //
375 // Returns integral PAI cross section for energy transfers >= transferCut
376 
377 G4double
378 G4PAIModel::GetdNdxCut( G4int iPlace, G4double transferCut)
379 {
380  G4int iTransfer;
381  G4double x1, x2, y1, y2, dNdxCut;
382  // G4cout<<"iPlace = "<<iPlace<<"; "<<"transferCut = "<<transferCut<<G4endl;
383  // G4cout<<"size = "<<G4int((*fPAItransferTable)(iPlace)->GetVectorLength())
384  // <<G4endl;
385  for( iTransfer = 0 ;
386  iTransfer < G4int((*fPAItransferTable)(iPlace)->GetVectorLength()) ;
387  iTransfer++)
388  {
389  if(transferCut <= (*fPAItransferTable)(iPlace)->GetLowEdgeEnergy(iTransfer))
390  {
391  break ;
392  }
393  }
394  if ( iTransfer >= G4int((*fPAItransferTable)(iPlace)->GetVectorLength()) )
395  {
396  iTransfer = (*fPAItransferTable)(iPlace)->GetVectorLength() - 1 ;
397  }
398  y1 = (*(*fPAItransferTable)(iPlace))(iTransfer-1) ;
399  y2 = (*(*fPAItransferTable)(iPlace))(iTransfer) ;
400  if(y1 < 0.0) y1 = 0.0;
401  if(y2 < 0.0) y2 = 0.0;
402  // G4cout<<"y1 = "<<y1<<"; "<<"y2 = "<<y2<<G4endl;
403  x1 = (*fPAItransferTable)(iPlace)->GetLowEdgeEnergy(iTransfer-1) ;
404  x2 = (*fPAItransferTable)(iPlace)->GetLowEdgeEnergy(iTransfer) ;
405  // G4cout<<"x1 = "<<x1<<"; "<<"x2 = "<<x2<<G4endl;
406 
407  if ( y1 == y2 ) dNdxCut = y2 ;
408  else
409  {
410  // if ( x1 == x2 ) dNdxCut = y1 + (y2 - y1)*G4UniformRand() ;
411  // if ( std::abs(x1-x2) <= eV ) dNdxCut = y1 + (y2 - y1)*G4UniformRand() ;
412  if ( std::abs(x1-x2) <= eV ) dNdxCut = y1 + (y2 - y1)*0.5 ;
413  else dNdxCut = y1 + (transferCut - x1)*(y2 - y1)/(x2 - x1) ;
414  }
415  // G4cout<<""<<dNdxCut<<G4endl;
416  if(dNdxCut < 0.0) dNdxCut = 0.0;
417  return dNdxCut ;
418 }
420 //
421 // Returns integral dEdx for energy transfers >= transferCut
422 
423 G4double
424 G4PAIModel::GetdEdxCut( G4int iPlace, G4double transferCut)
425 {
426  G4int iTransfer;
427  G4double x1, x2, y1, y2, dEdxCut;
428  //G4cout<<"iPlace = "<<iPlace<<"; "<<"transferCut = "<<transferCut<<G4endl;
429  // G4cout<<"size = "<<G4int((*fPAIdEdxTable)(iPlace)->GetVectorLength())
430  // <<G4endl;
431  for( iTransfer = 0 ;
432  iTransfer < G4int((*fPAIdEdxTable)(iPlace)->GetVectorLength()) ;
433  iTransfer++)
434  {
435  if(transferCut <= (*fPAIdEdxTable)(iPlace)->GetLowEdgeEnergy(iTransfer))
436  {
437  break ;
438  }
439  }
440  if ( iTransfer >= G4int((*fPAIdEdxTable)(iPlace)->GetVectorLength()) )
441  {
442  iTransfer = (*fPAIdEdxTable)(iPlace)->GetVectorLength() - 1 ;
443  }
444  y1 = (*(*fPAIdEdxTable)(iPlace))(iTransfer-1) ;
445  y2 = (*(*fPAIdEdxTable)(iPlace))(iTransfer) ;
446  if(y1 < 0.0) y1 = 0.0;
447  if(y2 < 0.0) y2 = 0.0;
448  //G4cout<<"y1 = "<<y1<<"; "<<"y2 = "<<y2<<G4endl;
449  x1 = (*fPAIdEdxTable)(iPlace)->GetLowEdgeEnergy(iTransfer-1) ;
450  x2 = (*fPAIdEdxTable)(iPlace)->GetLowEdgeEnergy(iTransfer) ;
451  //G4cout<<"x1 = "<<x1<<"; "<<"x2 = "<<x2<<G4endl;
452 
453  if ( y1 == y2 ) dEdxCut = y2 ;
454  else
455  {
456  // if ( x1 == x2 ) dEdxCut = y1 + (y2 - y1)*G4UniformRand() ;
457  // if ( std::abs(x1-x2) <= eV ) dEdxCut = y1 + (y2 - y1)*G4UniformRand() ;
458  if ( std::abs(x1-x2) <= eV ) dEdxCut = y1 + (y2 - y1)*0.5 ;
459  else dEdxCut = y1 + (transferCut - x1)*(y2 - y1)/(x2 - x1) ;
460  }
461  //G4cout<<""<<dEdxCut<<G4endl;
462  if(dEdxCut < 0.0) dEdxCut = 0.0;
463  return dEdxCut ;
464 }
465 
467 
468 G4double G4PAIModel::ComputeDEDXPerVolume(const G4Material*,
469  const G4ParticleDefinition* p,
470  G4double kineticEnergy,
471  G4double cutEnergy)
472 {
473  G4int iTkin,iPlace;
474 
475  //G4double cut = std::min(MaxSecondaryEnergy(p, kineticEnergy), cutEnergy);
476  G4double cut = cutEnergy;
477 
478  G4double massRatio = fMass/p->GetPDGMass();
479  G4double scaledTkin = kineticEnergy*massRatio;
480  G4double charge = p->GetPDGCharge();
481  G4double charge2 = charge*charge;
482  const G4MaterialCutsCouple* matCC = CurrentCouple();
483 
484  size_t jMat = 0;
485  for(;jMat < fMaterialCutsCoupleVector.size(); ++jMat )
486  {
487  if( matCC == fMaterialCutsCoupleVector[jMat] ) break;
488  }
489  //G4cout << "jMat= " << jMat
490  // << " jMax= " << fMaterialCutsCoupleVector.size()
491  // << " matCC: " << matCC;
492  //if(matCC) G4cout << " mat: " << matCC->GetMaterial()->GetName();
493  // G4cout << G4endl;
494  if(jMat == fMaterialCutsCoupleVector.size()) return 0.0;
495 
496  fPAIdEdxTable = fPAIdEdxBank[jMat];
497  fdEdxVector = fdEdxTable[jMat];
498  for(iTkin = 0 ; iTkin <= fTotBin ; iTkin++)
499  {
500  if(scaledTkin < fParticleEnergyVector->GetLowEdgeEnergy(iTkin)) break ;
501  }
502  //G4cout << "E= " << scaledTkin << " iTkin= " << iTkin
503  // << " " << fdEdxVector->GetVectorLength()<<G4endl;
504  iPlace = iTkin - 1;
505  if(iPlace < 0) iPlace = 0;
506  else if(iPlace >= fTotBin) iPlace = fTotBin-1;
507  G4double dEdx = charge2*( (*fdEdxVector)(iPlace) - GetdEdxCut(iPlace,cut) );
508  if( dEdx < 0.) dEdx = 0.;
509  return dEdx;
510 }
511 
513 
514 G4double G4PAIModel::CrossSectionPerVolume( const G4Material*,
515  const G4ParticleDefinition* p,
516  G4double kineticEnergy,
517  G4double cutEnergy,
518  G4double maxEnergy )
519 {
520  G4int iTkin,iPlace;
521  G4double tmax = std::min(MaxSecondaryEnergy(p, kineticEnergy), maxEnergy);
522  if(tmax <= cutEnergy) return 0.0;
523  G4double massRatio = fMass/p->GetPDGMass();
524  G4double scaledTkin = kineticEnergy*massRatio;
525  G4double charge = p->GetPDGCharge();
526  G4double charge2 = charge*charge, cross, cross1, cross2;
527  const G4MaterialCutsCouple* matCC = CurrentCouple();
528 
529  size_t jMat = 0;
530  for(;jMat < fMaterialCutsCoupleVector.size(); ++jMat )
531  {
532  if( matCC == fMaterialCutsCoupleVector[jMat] ) break;
533  }
534  if(jMat == fMaterialCutsCoupleVector.size()) return 0.0;
535 
536  fPAItransferTable = fPAIxscBank[jMat];
537 
538  for(iTkin = 0 ; iTkin <= fTotBin ; iTkin++)
539  {
540  if(scaledTkin < fParticleEnergyVector->GetLowEdgeEnergy(iTkin)) break ;
541  }
542  iPlace = iTkin - 1;
543  if(iPlace < 0) iPlace = 0;
544  else if(iPlace >= fTotBin) iPlace = fTotBin-1;
545 
546  //G4cout<<"iPlace = "<<iPlace<<"; tmax = "
547  // <<tmax<<"; cutEnergy = "<<cutEnergy<<G4endl;
548  cross1 = GetdNdxCut(iPlace,tmax) ;
549  //G4cout<<"cross1 = "<<cross1<<G4endl;
550  cross2 = GetdNdxCut(iPlace,cutEnergy) ;
551  //G4cout<<"cross2 = "<<cross2<<G4endl;
552  cross = (cross2-cross1)*charge2;
553  //G4cout<<"cross = "<<cross<<G4endl;
554  if( cross < DBL_MIN) cross = 0.0;
555  // if( cross2 < DBL_MIN) cross2 = DBL_MIN;
556 
557  // return cross2;
558  return cross;
559 }
560 
562 //
563 // It is analog of PostStepDoIt in terms of secondary electron.
564 //
565 
566 void G4PAIModel::SampleSecondaries(std::vector<G4DynamicParticle*>* vdp,
567  const G4MaterialCutsCouple* matCC,
568  const G4DynamicParticle* dp,
569  G4double tmin,
570  G4double maxEnergy)
571 {
572  size_t jMat = 0;
573  for(;jMat < fMaterialCutsCoupleVector.size(); ++jMat )
574  {
575  if( matCC == fMaterialCutsCoupleVector[jMat] ) break;
576  }
577  if(jMat == fMaterialCutsCoupleVector.size()) return;
578 
579  fPAItransferTable = fPAIxscBank[jMat];
580  fdNdxCutVector = fdNdxCutTable[jMat];
581 
582  G4double tmax = std::min(MaxSecondaryKinEnergy(dp), maxEnergy);
583  if( tmin >= tmax && fVerbose > 0)
584  {
585  G4cout<<"G4PAIModel::SampleSecondary: tmin >= tmax "<<G4endl;
586  }
587  G4ThreeVector direction= dp->GetMomentumDirection();
588  G4double particleMass = dp->GetMass();
589  G4double kineticEnergy = dp->GetKineticEnergy();
590 
591  G4double massRatio = fMass/particleMass;
592  G4double scaledTkin = kineticEnergy*massRatio;
593  G4double totalEnergy = kineticEnergy + particleMass;
594  G4double pSquare = kineticEnergy*(totalEnergy+particleMass);
595 
596  G4double deltaTkin = GetPostStepTransfer(scaledTkin);
597 
598  // G4cout<<"G4PAIModel::SampleSecondaries; deltaKIn = "<<deltaTkin/keV<<" keV "<<G4endl;
599 
600  if( deltaTkin <= 0. && fVerbose > 0)
601  {
602  G4cout<<"G4PAIModel::SampleSecondary e- deltaTkin = "<<deltaTkin<<G4endl;
603  }
604  if( deltaTkin <= 0.) return;
605 
606  if( deltaTkin > tmax) deltaTkin = tmax;
607 
608  G4double deltaTotalMomentum = sqrt(deltaTkin*(deltaTkin + 2. * electron_mass_c2 ));
609  G4double totalMomentum = sqrt(pSquare);
610  G4double costheta = deltaTkin*(totalEnergy + electron_mass_c2)
611  /(deltaTotalMomentum * totalMomentum);
612 
613  if( costheta > 0.99999 ) costheta = 0.99999;
614  G4double sintheta = 0.0;
615  G4double sin2 = 1. - costheta*costheta;
616  if( sin2 > 0.) sintheta = sqrt(sin2);
617 
618  // direction of the delta electron
619  G4double phi = twopi*G4UniformRand();
620  G4double dirx = sintheta*cos(phi), diry = sintheta*sin(phi), dirz = costheta;
621 
622  G4ThreeVector deltaDirection(dirx,diry,dirz);
623  deltaDirection.rotateUz(direction);
624  deltaDirection.unit();
625 
626  // primary change
627  kineticEnergy -= deltaTkin;
628  G4ThreeVector dir = totalMomentum*direction - deltaTotalMomentum*deltaDirection;
629  direction = dir.unit();
630  fParticleChange->SetProposedKineticEnergy(kineticEnergy);
631  fParticleChange->SetProposedMomentumDirection(direction);
632 
633  // create G4DynamicParticle object for e- delta ray
634  G4DynamicParticle* deltaRay = new G4DynamicParticle;
635  deltaRay->SetDefinition(G4Electron::Electron());
636  deltaRay->SetKineticEnergy( deltaTkin ); // !!! trick for last steps /2.0 ???
637  deltaRay->SetMomentumDirection(deltaDirection);
638 
639  vdp->push_back(deltaRay);
640 }
641 
642 
644 //
645 // Returns post step PAI energy transfer > cut electron energy according to passed
646 // scaled kinetic energy of particle
647 
648 G4double
649 G4PAIModel::GetPostStepTransfer( G4double scaledTkin )
650 {
651  // G4cout<<"G4PAIModel::GetPostStepTransfer"<<G4endl ;
652 
653  G4int iTkin, iTransfer, iPlace ;
654  G4double transfer = 0.0, position, dNdxCut1, dNdxCut2, E1, E2, W1, W2, W ;
655 
656  for(iTkin=0;iTkin<=fTotBin;iTkin++)
657  {
658  if(scaledTkin < fParticleEnergyVector->GetLowEdgeEnergy(iTkin)) break ;
659  }
660  iPlace = iTkin - 1 ;
661  // G4cout<<"from search, iPlace = "<<iPlace<<G4endl ;
662  if(iPlace < 0) iPlace = 0;
663  else if(iPlace >= fTotBin) iPlace = fTotBin-1;
664  dNdxCut1 = (*fdNdxCutVector)(iPlace) ;
665  // G4cout<<"dNdxCut1 = "<<dNdxCut1<<G4endl ;
666 
667 
668  if(iTkin == fTotBin) // Fermi plato, try from left
669  {
670  position = dNdxCut1*G4UniformRand() ;
671 
672  for( iTransfer = 0;
673  iTransfer < G4int((*fPAItransferTable)(iPlace)->GetVectorLength()); iTransfer++ )
674  {
675  if(position >= (*(*fPAItransferTable)(iPlace))(iTransfer)) break ;
676  }
677  transfer = GetEnergyTransfer(iPlace,position,iTransfer);
678  }
679  else
680  {
681  dNdxCut2 = (*fdNdxCutVector)(iPlace+1) ;
682  // G4cout<<"dNdxCut2 = "<<dNdxCut2<<G4endl ;
683  if(iTkin == 0) // Tkin is too small, trying from right only
684  {
685  position = dNdxCut2*G4UniformRand() ;
686 
687  for( iTransfer = 0;
688  iTransfer < G4int((*fPAItransferTable)(iPlace+1)->GetVectorLength()); iTransfer++ )
689  {
690  if(position >= (*(*fPAItransferTable)(iPlace+1))(iTransfer)) break ;
691  }
692  transfer = GetEnergyTransfer(iPlace+1,position,iTransfer);
693  }
694  else // general case: Tkin between two vectors of the material
695  {
696  E1 = fParticleEnergyVector->GetLowEdgeEnergy(iTkin - 1) ;
697  E2 = fParticleEnergyVector->GetLowEdgeEnergy(iTkin) ;
698  W = 1.0/(E2 - E1) ;
699  W1 = (E2 - scaledTkin)*W ;
700  W2 = (scaledTkin - E1)*W ;
701 
702  position = ( dNdxCut1*W1 + dNdxCut2*W2 )*G4UniformRand() ;
703 
704  // G4cout<<position<<"\t" ;
705 
706  G4int iTrMax1, iTrMax2, iTrMax;
707 
708  iTrMax1 = G4int((*fPAItransferTable)(iPlace)->GetVectorLength());
709  iTrMax2 = G4int((*fPAItransferTable)(iPlace+1)->GetVectorLength());
710 
711  if (iTrMax1 >= iTrMax2) iTrMax = iTrMax2;
712  else iTrMax = iTrMax1;
713 
714 
715  for( iTransfer = 0; iTransfer < iTrMax; iTransfer++ )
716  {
717  if( position >=
718  ( (*(*fPAItransferTable)(iPlace))(iTransfer)*W1 +
719  (*(*fPAItransferTable)(iPlace+1))(iTransfer)*W2) ) break ;
720  }
721  transfer = GetEnergyTransfer(iPlace,position,iTransfer);
722  }
723  }
724  // G4cout<<"PAImodel PostStepTransfer = "<<transfer/keV<<" keV"<<G4endl ;
725  if(transfer < 0.0 ) transfer = 0.0 ;
726  // if(transfer < DBL_MIN ) transfer = DBL_MIN;
727 
728  return transfer ;
729 }
730 
732 //
733 // Returns random PAI energy transfer according to passed
734 // indexes of particle kinetic
735 
736 G4double
737 G4PAIModel::GetEnergyTransfer( G4int iPlace, G4double position, G4int iTransfer )
738 {
739  G4int iTransferMax;
740  G4double x1, x2, y1, y2, energyTransfer;
741 
742  if(iTransfer == 0)
743  {
744  energyTransfer = (*fPAItransferTable)(iPlace)->GetLowEdgeEnergy(iTransfer);
745  }
746  else
747  {
748  iTransferMax = G4int((*fPAItransferTable)(iPlace)->GetVectorLength());
749 
750  if ( iTransfer >= iTransferMax ) iTransfer = iTransferMax - 1;
751 
752  y1 = (*(*fPAItransferTable)(iPlace))(iTransfer-1);
753  y2 = (*(*fPAItransferTable)(iPlace))(iTransfer);
754 
755  x1 = (*fPAItransferTable)(iPlace)->GetLowEdgeEnergy(iTransfer-1);
756  x2 = (*fPAItransferTable)(iPlace)->GetLowEdgeEnergy(iTransfer);
757 
758  if ( x1 == x2 ) energyTransfer = x2;
759  else
760  {
761  if ( y1 == y2 ) energyTransfer = x1 + (x2 - x1)*G4UniformRand();
762  else
763  {
764  energyTransfer = x1 + (position - y1)*(x2 - x1)/(y2 - y1);
765  }
766  }
767  }
768  return energyTransfer;
769 }
770 
772 
773 G4double G4PAIModel::SampleFluctuations( const G4Material* material,
774  const G4DynamicParticle* aParticle,
775  G4double&,
776  G4double& step,
777  G4double&)
778 {
779  size_t jMat = 0;
780  for(;jMat < fMaterialCutsCoupleVector.size(); ++jMat )
781  {
782  if( material == fMaterialCutsCoupleVector[jMat]->GetMaterial() ) break;
783  }
784  if(jMat == fMaterialCutsCoupleVector.size()) return 0.0;
785 
786  fPAItransferTable = fPAIxscBank[jMat];
787  fdNdxCutVector = fdNdxCutTable[jMat];
788 
789  G4int iTkin, iTransfer, iPlace;
790  G4long numOfCollisions=0;
791 
792  //G4cout<<"G4PAIModel::SampleFluctuations"<<G4endl ;
793  //G4cout<<"in: "<<fMaterialCutsCoupleVector[jMat]->GetMaterial()->GetName()<<G4endl;
794  G4double loss = 0.0, charge2 ;
795  G4double stepSum = 0., stepDelta, lambda, omega;
796  G4double position, E1, E2, W1, W2, W, dNdxCut1, dNdxCut2, meanNumber;
797  G4bool numb = true;
798  G4double Tkin = aParticle->GetKineticEnergy() ;
799  G4double MassRatio = fMass/aParticle->GetDefinition()->GetPDGMass() ;
800  G4double charge = aParticle->GetDefinition()->GetPDGCharge() ;
801  charge2 = charge*charge ;
802  G4double TkinScaled = Tkin*MassRatio ;
803 
804  for(iTkin=0;iTkin<=fTotBin;iTkin++)
805  {
806  if(TkinScaled < fParticleEnergyVector->GetLowEdgeEnergy(iTkin)) break ;
807  }
808  iPlace = iTkin - 1 ;
809  if(iPlace < 0) iPlace = 0;
810  else if(iPlace >= fTotBin) iPlace = fTotBin - 1;
811  //G4cout<<"from search, iPlace = "<<iPlace<<G4endl ;
812  dNdxCut1 = (*fdNdxCutVector)(iPlace) ;
813  //G4cout<<"dNdxCut1 = "<<dNdxCut1<<G4endl ;
814 
815  if(iTkin == fTotBin) // Fermi plato, try from left
816  {
817  meanNumber =((*(*fPAItransferTable)(iPlace))(0)-dNdxCut1)*step*charge2;
818  if(meanNumber < 0.) meanNumber = 0. ;
819  // numOfCollisions = RandPoisson::shoot(meanNumber) ;
820  // numOfCollisions = G4Poisson(meanNumber) ;
821  if( meanNumber > 0.) lambda = step/meanNumber;
822  else lambda = DBL_MAX;
823  while(numb)
824  {
825  stepDelta = CLHEP::RandExponential::shoot(lambda);
826  stepSum += stepDelta;
827  if(stepSum >= step) break;
828  numOfCollisions++;
829  }
830  //G4cout<<"##1 numOfCollisions = "<<numOfCollisions<<G4endl ;
831 
832  while(numOfCollisions)
833  {
834  position = dNdxCut1+
835  ((*(*fPAItransferTable)(iPlace))(0)-dNdxCut1)*G4UniformRand() ;
836 
837  for( iTransfer = 0;
838  iTransfer < G4int((*fPAItransferTable)(iPlace)->GetVectorLength()); iTransfer++ )
839  {
840  if(position >= (*(*fPAItransferTable)(iPlace))(iTransfer)) break ;
841  }
842  omega = GetEnergyTransfer(iPlace,position,iTransfer);
843  //G4cout<<"G4PAIModel::SampleFluctuations, omega = "<<omega/keV<<" keV; "<<"\t";
844  loss += omega;
845  numOfCollisions-- ;
846  }
847  }
848  else
849  {
850  dNdxCut2 = (*fdNdxCutVector)(iPlace+1) ;
851  //G4cout<<"dNdxCut2 = "<<dNdxCut2<< " iTkin= "<<iTkin<<" iPlace= "<<iPlace<<G4endl;
852 
853  if(iTkin == 0) // Tkin is too small, trying from right only
854  {
855  meanNumber =((*(*fPAItransferTable)(iPlace+1))(0)-dNdxCut2)*step*charge2;
856  if( meanNumber < 0. ) meanNumber = 0. ;
857  // numOfCollisions = CLHEP::RandPoisson::shoot(meanNumber) ;
858  // numOfCollisions = G4Poisson(meanNumber) ;
859  if( meanNumber > 0.) lambda = step/meanNumber;
860  else lambda = DBL_MAX;
861  while(numb)
862  {
863  stepDelta = CLHEP::RandExponential::shoot(lambda);
864  stepSum += stepDelta;
865  if(stepSum >= step) break;
866  numOfCollisions++;
867  }
868 
869  //G4cout<<"##2 numOfCollisions = "<<numOfCollisions<<G4endl ;
870 
871  while(numOfCollisions)
872  {
873  position = dNdxCut2+
874  ((*(*fPAItransferTable)(iPlace+1))(0)-dNdxCut2)*G4UniformRand();
875 
876  for( iTransfer = 0;
877  iTransfer < G4int((*fPAItransferTable)(iPlace+1)->GetVectorLength()); iTransfer++ )
878  {
879  if(position >= (*(*fPAItransferTable)(iPlace+1))(iTransfer)) break ;
880  }
881  omega = GetEnergyTransfer(iPlace,position,iTransfer);
882  //G4cout<<omega/keV<<"\t";
883  loss += omega;
884  numOfCollisions-- ;
885  }
886  }
887  else // general case: Tkin between two vectors of the material
888  {
889  E1 = fParticleEnergyVector->GetLowEdgeEnergy(iTkin - 1) ;
890  E2 = fParticleEnergyVector->GetLowEdgeEnergy(iTkin) ;
891  W = 1.0/(E2 - E1) ;
892  W1 = (E2 - TkinScaled)*W ;
893  W2 = (TkinScaled - E1)*W ;
894 
895  //G4cout << fPAItransferTable->size() << G4endl;
896  //G4cout<<"(*(*fPAItransferTable)(iPlace))(0) = "<<
897  // (*(*fPAItransferTable)(iPlace))(0)<<G4endl ;
898  //G4cout<<"(*(*fPAItransferTable)(iPlace+1))(0) = "<<
899  // (*(*fPAItransferTable)(iPlace+1))(0)<<G4endl ;
900 
901  meanNumber=( ((*(*fPAItransferTable)(iPlace))(0)-dNdxCut1)*W1 +
902  ((*(*fPAItransferTable)(iPlace+1))(0)-dNdxCut2)*W2 )*step*charge2;
903  if(meanNumber<0.0) meanNumber = 0.0;
904  // numOfCollisions = RandPoisson::shoot(meanNumber) ;
905  // numOfCollisions = G4Poisson(meanNumber) ;
906  if( meanNumber > 0.) lambda = step/meanNumber;
907  else lambda = DBL_MAX;
908  while(numb)
909  {
910  stepDelta = CLHEP::RandExponential::shoot(lambda);
911  stepSum += stepDelta;
912  if(stepSum >= step) break;
913  numOfCollisions++;
914  }
915 
916  //G4cout<<"##3 numOfCollisions = "<<numOfCollisions<<endl ;
917 
918  while(numOfCollisions)
919  {
920  position = dNdxCut1*W1 + dNdxCut2*W2 +
921  ( ( (*(*fPAItransferTable)(iPlace))(0)-dNdxCut1 )*W1 +
922  dNdxCut2+
923  ( (*(*fPAItransferTable)(iPlace+1))(0)-dNdxCut2 )*W2 )*G4UniformRand();
924 
925  // G4cout<<position<<"\t" ;
926 
927  for( iTransfer = 0;
928  iTransfer < G4int((*fPAItransferTable)(iPlace)->GetVectorLength()); iTransfer++ )
929  {
930  if( position >=
931  ( (*(*fPAItransferTable)(iPlace))(iTransfer)*W1 +
932  (*(*fPAItransferTable)(iPlace+1))(iTransfer)*W2) )
933  {
934  break ;
935  }
936  }
937  omega = GetEnergyTransfer(iPlace,position,iTransfer);
938  // G4cout<<omega/keV<<"\t";
939  loss += omega;
940  numOfCollisions-- ;
941  }
942  }
943  }
944  //G4cout<<"PAIModel AlongStepLoss = "<<loss/keV<<" keV, on step = "
945  // <<step/mm<<" mm"<<G4endl ;
946  if(loss > Tkin) loss=Tkin;
947  if(loss < 0. ) loss = 0.;
948  return loss ;
949 
950 }
951 
953 //
954 // Returns the statistical estimation of the energy loss distribution variance
955 //
956 
957 
958 G4double G4PAIModel::Dispersion( const G4Material* material,
959  const G4DynamicParticle* aParticle,
960  G4double& tmax,
961  G4double& step )
962 {
963  G4double loss, sumLoss=0., sumLoss2=0., sigma2, meanLoss=0.;
964  for(G4int i = 0 ; i < fMeanNumber; i++)
965  {
966  loss = SampleFluctuations(material,aParticle,tmax,step,meanLoss);
967  sumLoss += loss;
968  sumLoss2 += loss*loss;
969  }
970  meanLoss = sumLoss/fMeanNumber;
971  sigma2 = meanLoss*meanLoss + (sumLoss2-2*sumLoss*meanLoss)/fMeanNumber;
972  return sigma2;
973 }
974 
976 
977 G4double G4PAIModel::MaxSecondaryEnergy( const G4ParticleDefinition* p,
978  G4double kinEnergy)
979 {
980  G4double tmax = kinEnergy;
981  if(p == fElectron) tmax *= 0.5;
982  else if(p != fPositron) {
983  G4double mass = p->GetPDGMass();
984  G4double ratio= electron_mass_c2/mass;
985  G4double gamma= kinEnergy/mass + 1.0;
986  tmax = 2.0*electron_mass_c2*(gamma*gamma - 1.) /
987  (1. + 2.0*gamma*ratio + ratio*ratio);
988  }
989  return tmax;
990 }
991 
993 
994 void G4PAIModel::DefineForRegion(const G4Region* r)
995 {
996  fPAIRegionVector.push_back(r);
997 }
998 
999 //
1000 //
1002 
1003 
1004 
1005 
1006 
1007 
Generated on Sat May 25 2013 14:33:37 for Geant4 by   doxygen 1.8.4