Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4PAIPhotonModel.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: G4PAIPhotonModel.cc
32 //
33 // Author: Vladimir.Grichine@cern.ch based on G4PAIModel class
34 //
35 // Creation date: 20.05.2004
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 // 11.04.05 Major optimisation of internal interfaces (V.Ivantchenko)
42 //
43 
44 #include "G4PAIPhotonModel.hh"
45 
46 #include "G4PhysicalConstants.hh"
47 #include "G4SystemOfUnits.hh"
48 
49 #include "G4Region.hh"
50 #include "G4PhysicsLogVector.hh"
51 #include "G4PhysicsFreeVector.hh"
52 #include "G4PhysicsTable.hh"
53 #include "G4ProductionCutsTable.hh"
54 #include "G4MaterialCutsCouple.hh"
55 #include "G4MaterialTable.hh"
56 #include "G4SandiaTable.hh"
57 #include "G4PAIxSection.hh"
58 
59 #include "Randomize.hh"
60 #include "G4Electron.hh"
61 #include "G4Positron.hh"
62 #include "G4Gamma.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 #include "G4GeometryTolerance.hh"
70 
72 
73 using namespace std;
74 
75 G4PAIPhotonModel::G4PAIPhotonModel(const G4ParticleDefinition* p, const G4String& nam)
76  : G4VEmModel(nam),G4VEmFluctuationModel(nam),
77  fLowestKineticEnergy(10.0*keV),
78  fHighestKineticEnergy(100.*TeV),
79  fTotBin(200),
80  fMeanNumber(20),
81  fParticle(0),
82  fHighKinEnergy(100.*TeV),
83  fLowKinEnergy(2.0*MeV),
84  fTwoln10(2.0*log(10.0)),
85  fBg2lim(0.0169),
86  fTaulim(8.4146e-3)
87 {
88  fVerbose = 0;
89  fElectron = G4Electron::Electron();
90  fPositron = G4Positron::Positron();
91 
92  fProtonEnergyVector = new G4PhysicsLogVector(fLowestKineticEnergy,
93  fHighestKineticEnergy,
94  fTotBin);
95  fPAItransferTable = 0;
96  fPAIphotonTable = 0;
97  fPAIplasmonTable = 0;
98 
99  fPAIdEdxTable = 0;
100  fSandiaPhotoAbsCof = 0;
101  fdEdxVector = 0;
102 
103  fLambdaVector = 0;
104  fdNdxCutVector = 0;
105  fdNdxCutPhotonVector = 0;
106  fdNdxCutPlasmonVector = 0;
107 
108  fSandiaIntervalNumber = 0;
109  fMatIndex = 0;
110 
111  fParticleChange = 0;
112 
113  if(p) { SetParticle(p); }
114  else { SetParticle(fElectron); }
115 
116  isInitialised = false;
117 }
118 
120 
121 G4PAIPhotonModel::~G4PAIPhotonModel()
122 {
123  if(fProtonEnergyVector) delete fProtonEnergyVector;
124  if(fdEdxVector) delete fdEdxVector ;
125  if ( fLambdaVector) delete fLambdaVector;
126  if ( fdNdxCutVector) delete fdNdxCutVector;
127 
128  if( fPAItransferTable )
129  {
130  fPAItransferTable->clearAndDestroy();
131  delete fPAItransferTable ;
132  }
133  if( fPAIphotonTable )
134  {
135  fPAIphotonTable->clearAndDestroy();
136  delete fPAIphotonTable ;
137  }
138  if( fPAIplasmonTable )
139  {
140  fPAIplasmonTable->clearAndDestroy();
141  delete fPAIplasmonTable ;
142  }
143  if(fSandiaPhotoAbsCof)
144  {
145  for(G4int i=0;i<fSandiaIntervalNumber;i++)
146  {
147  delete[] fSandiaPhotoAbsCof[i];
148  }
149  delete[] fSandiaPhotoAbsCof;
150  }
151 }
152 
154 
155 void G4PAIPhotonModel::SetParticle(const G4ParticleDefinition* p)
156 {
157  fParticle = p;
158  fMass = fParticle->GetPDGMass();
159  fSpin = fParticle->GetPDGSpin();
160  G4double q = fParticle->GetPDGCharge()/eplus;
161  fChargeSquare = q*q;
162  fLowKinEnergy *= fMass/proton_mass_c2;
163  fRatio = electron_mass_c2/fMass;
164  fQc = fMass/fRatio;
165 }
166 
168 
169 void G4PAIPhotonModel::Initialise(const G4ParticleDefinition* p,
170  const G4DataVector&)
171 {
172  // G4cout<<"G4PAIPhotonModel::Initialise for "<<p->GetParticleName()<<G4endl;
173  if(isInitialised) { return; }
174  isInitialised = true;
175 
176  if(!fParticle) SetParticle(p);
177 
178  fParticleChange = GetParticleChangeForLoss();
179 
180  const G4ProductionCutsTable* theCoupleTable =
181  G4ProductionCutsTable::GetProductionCutsTable();
182 
183  for(size_t iReg = 0; iReg < fPAIRegionVector.size();++iReg) // region loop
184  {
185  const G4Region* curReg = fPAIRegionVector[iReg];
186 
187  vector<G4Material*>::const_iterator matIter = curReg->GetMaterialIterator();
188  size_t jMat;
189  size_t numOfMat = curReg->GetNumberOfMaterials();
190 
191  // for(size_t jMat = 0; jMat < curReg->GetNumberOfMaterials();++jMat){}
192  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
193  size_t numberOfMat = G4Material::GetNumberOfMaterials();
194 
195  for(jMat = 0 ; jMat < numOfMat; ++jMat) // region material loop
196  {
197  const G4MaterialCutsCouple* matCouple = theCoupleTable->
198  GetMaterialCutsCouple( *matIter, curReg->GetProductionCuts() );
199  fMaterialCutsCoupleVector.push_back(matCouple);
200 
201  size_t iMatGlob;
202  for(iMatGlob = 0 ; iMatGlob < numberOfMat ; iMatGlob++ )
203  {
204  if( *matIter == (*theMaterialTable)[iMatGlob]) break ;
205  }
206  fMatIndex = iMatGlob;
207 
208  ComputeSandiaPhotoAbsCof();
209  BuildPAIonisationTable();
210 
211  fPAIxscBank.push_back(fPAItransferTable);
212  fPAIphotonBank.push_back(fPAIphotonTable);
213  fPAIplasmonBank.push_back(fPAIplasmonTable);
214  fPAIdEdxBank.push_back(fPAIdEdxTable);
215  fdEdxTable.push_back(fdEdxVector);
216 
217  BuildLambdaVector(matCouple);
218 
219  fdNdxCutTable.push_back(fdNdxCutVector);
220  fdNdxCutPhotonTable.push_back(fdNdxCutPhotonVector);
221  fdNdxCutPlasmonTable.push_back(fdNdxCutPlasmonVector);
222  fLambdaTable.push_back(fLambdaVector);
223 
224 
225  matIter++;
226  }
227  }
228 }
229 
231 
232 void G4PAIPhotonModel::InitialiseMe(const G4ParticleDefinition*)
233 {}
234 
236 
237 void G4PAIPhotonModel::ComputeSandiaPhotoAbsCof()
238 {
239  G4int i, j, numberOfElements ;
240  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
241 
242  G4SandiaTable thisMaterialSandiaTable(fMatIndex) ;
243  numberOfElements = (*theMaterialTable)[fMatIndex]->
244  GetNumberOfElements();
245  G4int* thisMaterialZ = new G4int[numberOfElements] ;
246 
247  for(i=0;i<numberOfElements;i++)
248  {
249  thisMaterialZ[i] =
250  (G4int)(*theMaterialTable)[fMatIndex]->GetElement(i)->GetZ() ;
251  }
252  fSandiaIntervalNumber = thisMaterialSandiaTable.SandiaIntervals
253  (thisMaterialZ,numberOfElements) ;
254 
255  fSandiaIntervalNumber = thisMaterialSandiaTable.SandiaMixing
256  ( thisMaterialZ ,
257  (*theMaterialTable)[fMatIndex]->GetFractionVector() ,
258  numberOfElements,fSandiaIntervalNumber) ;
259 
260  fSandiaPhotoAbsCof = new G4double*[fSandiaIntervalNumber] ;
261 
262  for(i=0;i<fSandiaIntervalNumber;i++) fSandiaPhotoAbsCof[i] = new G4double[5] ;
263 
264  for( i = 0 ; i < fSandiaIntervalNumber ; i++ )
265  {
266  fSandiaPhotoAbsCof[i][0] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i+1,0) ;
267 
268  for( j = 1; j < 5 ; j++ )
269  {
270  fSandiaPhotoAbsCof[i][j] = thisMaterialSandiaTable.
271  GetPhotoAbsorpCof(i+1,j)*
272  (*theMaterialTable)[fMatIndex]->GetDensity() ;
273  }
274  }
275  delete[] thisMaterialZ ;
276 }
277 
279 //
280 // Build tables for the ionization energy loss
281 // the tables are built for MATERIALS
282 // *********
283 
284 void
285 G4PAIPhotonModel::BuildPAIonisationTable()
286 {
287  G4double LowEdgeEnergy , ionloss ;
288  G4double /*massRatio,*/ tau, Tmax, Tmin, Tkin, deltaLow, /*gamma,*/ bg2 ;
289  /*
290  if( fPAItransferTable )
291  {
292  fPAItransferTable->clearAndDestroy() ;
293  delete fPAItransferTable ;
294  }
295  */
296  fPAItransferTable = new G4PhysicsTable(fTotBin);
297  /*
298  if( fPAIratioTable )
299  {
300  fPAIratioTable->clearAndDestroy() ;
301  delete fPAIratioTable ;
302  }
303  */
304  fPAIphotonTable = new G4PhysicsTable(fTotBin);
305  fPAIplasmonTable = new G4PhysicsTable(fTotBin);
306  /*
307  if( fPAIdEdxTable )
308  {
309  fPAIdEdxTable->clearAndDestroy() ;
310  delete fPAIdEdxTable ;
311  }
312  */
313  fPAIdEdxTable = new G4PhysicsTable(fTotBin);
314 
315  // if(fdEdxVector) delete fdEdxVector ;
316  fdEdxVector = new G4PhysicsLogVector( fLowestKineticEnergy,
317  fHighestKineticEnergy,
318  fTotBin ) ;
319  Tmin = fSandiaPhotoAbsCof[0][0] ; // low energy Sandia interval
320  deltaLow = 100.*eV; // 0.5*eV ;
321 
322  for (G4int i = 0 ; i <= fTotBin ; i++) //The loop for the kinetic energy
323  {
324  LowEdgeEnergy = fProtonEnergyVector->GetLowEdgeEnergy(i) ;
325  tau = LowEdgeEnergy/proton_mass_c2 ;
326  // if(tau < 0.01) tau = 0.01 ;
327  //gamma = tau +1. ;
328  // G4cout<<"gamma = "<<gamma<<endl ;
329  bg2 = tau*(tau + 2. ) ;
330  //massRatio = electron_mass_c2/proton_mass_c2 ;
331  Tmax = MaxSecondaryEnergy(fParticle, LowEdgeEnergy);
332  // G4cout<<"proton Tkin = "<<LowEdgeEnergy/MeV<<" MeV"
333  // <<" Tmax = "<<Tmax/MeV<<" MeV"<<G4endl;
334  // Tkin = DeltaCutInKineticEnergyNow ;
335 
336  // if ( DeltaCutInKineticEnergyNow > Tmax) // was <
337  Tkin = Tmax ;
338  if ( Tkin < Tmin + deltaLow ) // low energy safety
339  {
340  Tkin = Tmin + deltaLow ;
341  }
342  G4PAIxSection protonPAI( fMatIndex,
343  Tkin,
344  bg2,
345  fSandiaPhotoAbsCof,
346  fSandiaIntervalNumber ) ;
347 
348 
349  // G4cout<<"ionloss = "<<ionloss*cm/keV<<" keV/cm"<<endl ;
350  // G4cout<<"n1 = "<<protonPAI.GetIntegralPAIxSection(1)*cm<<" 1/cm"<<endl ;
351  // G4cout<<"protonPAI.GetSplineSize() = "<<
352  // protonPAI.GetSplineSize()<<G4endl<<G4endl ;
353 
354  G4PhysicsFreeVector* transferVector = new
355  G4PhysicsFreeVector(protonPAI.GetSplineSize()) ;
356  G4PhysicsFreeVector* photonVector = new
357  G4PhysicsFreeVector(protonPAI.GetSplineSize()) ;
358  G4PhysicsFreeVector* plasmonVector = new
359  G4PhysicsFreeVector(protonPAI.GetSplineSize()) ;
360  G4PhysicsFreeVector* dEdxVector = new
361  G4PhysicsFreeVector(protonPAI.GetSplineSize()) ;
362 
363  for( G4int k = 0 ; k < protonPAI.GetSplineSize() ; k++ )
364  {
365  transferVector->PutValue( k ,
366  protonPAI.GetSplineEnergy(k+1),
367  protonPAI.GetIntegralPAIxSection(k+1) ) ;
368  photonVector->PutValue( k ,
369  protonPAI.GetSplineEnergy(k+1),
370  protonPAI.GetIntegralCerenkov(k+1) ) ;
371  plasmonVector->PutValue( k ,
372  protonPAI.GetSplineEnergy(k+1),
373  protonPAI.GetIntegralPlasmon(k+1) ) ;
374  dEdxVector->PutValue( k ,
375  protonPAI.GetSplineEnergy(k+1),
376  protonPAI.GetIntegralPAIdEdx(k+1) ) ;
377  }
378  ionloss = protonPAI.GetMeanEnergyLoss() ; // total <dE/dx>
379  if ( ionloss <= 0.) ionloss = DBL_MIN ;
380  fdEdxVector->PutValue(i,ionloss) ;
381 
382  fPAItransferTable->insertAt(i,transferVector) ;
383  fPAIphotonTable->insertAt(i,photonVector) ;
384  fPAIplasmonTable->insertAt(i,plasmonVector) ;
385  fPAIdEdxTable->insertAt(i,dEdxVector) ;
386 
387  } // end of Tkin loop
388  // theLossTable->insert(fdEdxVector);
389  // end of material loop
390  // G4cout<<"G4PAIonisation::BuildPAIonisationTable() have been called"<<G4endl ;
391  // G4cout<<"G4PAIonisation::BuildLossTable() have been called"<<G4endl ;
392 }
393 
395 //
396 // Build mean free path tables for the delta ray production process
397 // tables are built for MATERIALS
398 //
399 
400 void
401 G4PAIPhotonModel::BuildLambdaVector(const G4MaterialCutsCouple* matCutsCouple)
402 {
403  G4int i ;
404  G4double dNdxCut,dNdxPhotonCut,dNdxPlasmonCut, lambda;
405  G4double kCarTolerance = G4GeometryTolerance::GetInstance()
406  ->GetSurfaceTolerance();
407 
408  const G4ProductionCutsTable* theCoupleTable=
409  G4ProductionCutsTable::GetProductionCutsTable();
410 
411  size_t numOfCouples = theCoupleTable->GetTableSize();
412  size_t jMatCC;
413 
414  for (jMatCC = 0 ; jMatCC < numOfCouples ; jMatCC++ )
415  {
416  if( matCutsCouple == theCoupleTable->GetMaterialCutsCouple(jMatCC) ) break;
417  }
418  if( jMatCC == numOfCouples && jMatCC > 0 ) jMatCC--;
419 
420  const vector<G4double>* deltaCutInKineticEnergy = theCoupleTable->
421  GetEnergyCutsVector(idxG4ElectronCut);
422  const vector<G4double>* photonCutInKineticEnergy = theCoupleTable->
423  GetEnergyCutsVector(idxG4GammaCut);
424 
425  if (fLambdaVector) delete fLambdaVector;
426  if (fdNdxCutVector) delete fdNdxCutVector;
427  if (fdNdxCutPhotonVector) delete fdNdxCutPhotonVector;
428  if (fdNdxCutPlasmonVector) delete fdNdxCutPlasmonVector;
429 
430  fLambdaVector = new G4PhysicsLogVector( fLowestKineticEnergy,
431  fHighestKineticEnergy,
432  fTotBin );
433  fdNdxCutVector = new G4PhysicsLogVector( fLowestKineticEnergy,
434  fHighestKineticEnergy,
435  fTotBin );
436  fdNdxCutPhotonVector = new G4PhysicsLogVector( fLowestKineticEnergy,
437  fHighestKineticEnergy,
438  fTotBin );
439  fdNdxCutPlasmonVector = new G4PhysicsLogVector( fLowestKineticEnergy,
440  fHighestKineticEnergy,
441  fTotBin );
442 
443  G4double deltaCutInKineticEnergyNow = (*deltaCutInKineticEnergy)[jMatCC];
444  G4double photonCutInKineticEnergyNow = (*photonCutInKineticEnergy)[jMatCC];
445 
446  if(fVerbose > 0)
447  {
448  G4cout<<"PAIPhotonModel deltaCutInKineticEnergyNow = "
449  <<deltaCutInKineticEnergyNow/keV<<" keV"<<G4endl;
450  G4cout<<"PAIPhotonModel photonCutInKineticEnergyNow = "
451  <<photonCutInKineticEnergyNow/keV<<" keV"<<G4endl;
452  }
453  for ( i = 0 ; i <= fTotBin ; i++ )
454  {
455  dNdxPhotonCut = GetdNdxPhotonCut(i,photonCutInKineticEnergyNow);
456  dNdxPlasmonCut = GetdNdxPlasmonCut(i,deltaCutInKineticEnergyNow);
457 
458  dNdxCut = dNdxPhotonCut + dNdxPlasmonCut;
459  lambda = dNdxCut <= DBL_MIN ? DBL_MAX: 1.0/dNdxCut;
460 
461  if (lambda <= 1000*kCarTolerance) lambda = 1000*kCarTolerance; // Mmm ???
462 
463  fLambdaVector->PutValue(i, lambda);
464 
465  fdNdxCutVector->PutValue(i, dNdxCut);
466  fdNdxCutPhotonVector->PutValue(i, dNdxPhotonCut);
467  fdNdxCutPlasmonVector->PutValue(i, dNdxPlasmonCut);
468  }
469 }
470 
472 //
473 // Returns integral PAI cross section for energy transfers >= transferCut
474 
475 G4double
476 G4PAIPhotonModel::GetdNdxCut( G4int iPlace, G4double transferCut)
477 {
478  G4int iTransfer;
479  G4double x1, x2, y1, y2, dNdxCut;
480  // G4cout<<"iPlace = "<<iPlace<<"; "<<"transferCut = "<<transferCut<<G4endl;
481  // G4cout<<"size = "<<G4int((*fPAItransferTable)(iPlace)->GetVectorLength())
482  // <<G4endl;
483  for( iTransfer = 0 ;
484  iTransfer < G4int((*fPAItransferTable)(iPlace)->GetVectorLength()) ;
485  iTransfer++)
486  {
487  if(transferCut <= (*fPAItransferTable)(iPlace)->GetLowEdgeEnergy(iTransfer))
488  {
489  break ;
490  }
491  }
492  if ( iTransfer >= G4int((*fPAItransferTable)(iPlace)->GetVectorLength()) )
493  {
494  iTransfer = (*fPAItransferTable)(iPlace)->GetVectorLength() - 1 ;
495  }
496  y1 = (*(*fPAItransferTable)(iPlace))(iTransfer-1) ;
497  y2 = (*(*fPAItransferTable)(iPlace))(iTransfer) ;
498  // G4cout<<"y1 = "<<y1<<"; "<<"y2 = "<<y2<<G4endl;
499  x1 = (*fPAItransferTable)(iPlace)->GetLowEdgeEnergy(iTransfer-1) ;
500  x2 = (*fPAItransferTable)(iPlace)->GetLowEdgeEnergy(iTransfer) ;
501  // G4cout<<"x1 = "<<x1<<"; "<<"x2 = "<<x2<<G4endl;
502 
503  if ( y1 == y2 ) dNdxCut = y2 ;
504  else
505  {
506  // if ( x1 == x2 ) dNdxCut = y1 + (y2 - y1)*G4UniformRand() ;
507  // if ( std::abs(x1-x2) <= eV ) dNdxCut = y1 + (y2 - y1)*G4UniformRand() ;
508  if ( std::abs(x1-x2) <= eV ) dNdxCut = y1 + (y2 - y1)*0.5 ;
509  else dNdxCut = y1 + (transferCut - x1)*(y2 - y1)/(x2 - x1) ;
510  }
511  // G4cout<<""<<dNdxCut<<G4endl;
512  return dNdxCut ;
513 }
514 
516 //
517 // Returns integral PAI cherenkovcross section for energy transfers >= transferCut
518 
519 G4double
520 G4PAIPhotonModel::GetdNdxPhotonCut( G4int iPlace, G4double transferCut)
521 {
522  G4int iTransfer;
523  G4double x1, x2, y1, y2, dNdxCut;
524  // G4cout<<"iPlace = "<<iPlace<<"; "<<"transferCut = "<<transferCut<<G4endl;
525  // G4cout<<"size = "<<G4int((*fPAIphotonTable)(iPlace)->GetVectorLength())
526  // <<G4endl;
527  for( iTransfer = 0 ;
528  iTransfer < G4int((*fPAIphotonTable)(iPlace)->GetVectorLength()) ;
529  iTransfer++)
530  {
531  if(transferCut <= (*fPAIphotonTable)(iPlace)->GetLowEdgeEnergy(iTransfer))
532  {
533  break ;
534  }
535  }
536  if ( iTransfer >= G4int((*fPAIphotonTable)(iPlace)->GetVectorLength()) )
537  {
538  iTransfer = (*fPAIphotonTable)(iPlace)->GetVectorLength() - 1 ;
539  }
540  y1 = (*(*fPAIphotonTable)(iPlace))(iTransfer-1) ;
541  y2 = (*(*fPAIphotonTable)(iPlace))(iTransfer) ;
542  // G4cout<<"y1 = "<<y1<<"; "<<"y2 = "<<y2<<G4endl;
543  x1 = (*fPAIphotonTable)(iPlace)->GetLowEdgeEnergy(iTransfer-1) ;
544  x2 = (*fPAIphotonTable)(iPlace)->GetLowEdgeEnergy(iTransfer) ;
545  // G4cout<<"x1 = "<<x1<<"; "<<"x2 = "<<x2<<G4endl;
546 
547  if ( y1 == y2 ) dNdxCut = y2 ;
548  else
549  {
550  // if ( x1 == x2 ) dNdxCut = y1 + (y2 - y1)*G4UniformRand() ;
551  // if ( std::abs(x1-x2) <= eV ) dNdxCut = y1 + (y2 - y1)*G4UniformRand() ;
552  if ( std::abs(x1-x2) <= eV ) dNdxCut = y1 + (y2 - y1)*0.5 ;
553  else dNdxCut = y1 + (transferCut - x1)*(y2 - y1)/(x2 - x1) ;
554  }
555  // G4cout<<""<<dNdxPhotonCut<<G4endl;
556  return dNdxCut ;
557 }
558 
560 //
561 // Returns integral PAI cross section for energy transfers >= transferCut
562 
563 G4double
564 G4PAIPhotonModel::GetdNdxPlasmonCut( G4int iPlace, G4double transferCut)
565 {
566  G4int iTransfer;
567  G4double x1, x2, y1, y2, dNdxCut;
568 
569  // G4cout<<"iPlace = "<<iPlace<<"; "<<"transferCut = "<<transferCut<<G4endl;
570  // G4cout<<"size = "<<G4int((*fPAIPlasmonTable)(iPlace)->GetVectorLength())
571  // <<G4endl;
572  for( iTransfer = 0 ;
573  iTransfer < G4int((*fPAIplasmonTable)(iPlace)->GetVectorLength()) ;
574  iTransfer++)
575  {
576  if(transferCut <= (*fPAIplasmonTable)(iPlace)->GetLowEdgeEnergy(iTransfer))
577  {
578  break ;
579  }
580  }
581  if ( iTransfer >= G4int((*fPAIplasmonTable)(iPlace)->GetVectorLength()) )
582  {
583  iTransfer = (*fPAIplasmonTable)(iPlace)->GetVectorLength() - 1 ;
584  }
585  y1 = (*(*fPAIplasmonTable)(iPlace))(iTransfer-1) ;
586  y2 = (*(*fPAIplasmonTable)(iPlace))(iTransfer) ;
587  // G4cout<<"y1 = "<<y1<<"; "<<"y2 = "<<y2<<G4endl;
588  x1 = (*fPAIplasmonTable)(iPlace)->GetLowEdgeEnergy(iTransfer-1) ;
589  x2 = (*fPAIplasmonTable)(iPlace)->GetLowEdgeEnergy(iTransfer) ;
590  // G4cout<<"x1 = "<<x1<<"; "<<"x2 = "<<x2<<G4endl;
591 
592  if ( y1 == y2 ) dNdxCut = y2 ;
593  else
594  {
595  // if ( x1 == x2 ) dNdxCut = y1 + (y2 - y1)*G4UniformRand() ;
596  // if ( std::abs(x1-x2) <= eV ) dNdxCut = y1 + (y2 - y1)*G4UniformRand() ;
597  if ( std::abs(x1-x2) <= eV ) dNdxCut = y1 + (y2 - y1)*0.5 ;
598  else dNdxCut = y1 + (transferCut - x1)*(y2 - y1)/(x2 - x1) ;
599  }
600  // G4cout<<""<<dNdxPlasmonCut<<G4endl;
601  return dNdxCut ;
602 }
603 
605 //
606 // Returns integral dEdx for energy transfers >= transferCut
607 
608 G4double
609 G4PAIPhotonModel::GetdEdxCut( G4int iPlace, G4double transferCut)
610 {
611  G4int iTransfer;
612  G4double x1, x2, y1, y2, dEdxCut;
613  // G4cout<<"iPlace = "<<iPlace<<"; "<<"transferCut = "<<transferCut<<G4endl;
614  // G4cout<<"size = "<<G4int((*fPAIdEdxTable)(iPlace)->GetVectorLength())
615  // <<G4endl;
616  for( iTransfer = 0 ;
617  iTransfer < G4int((*fPAIdEdxTable)(iPlace)->GetVectorLength()) ;
618  iTransfer++)
619  {
620  if(transferCut <= (*fPAIdEdxTable)(iPlace)->GetLowEdgeEnergy(iTransfer))
621  {
622  break ;
623  }
624  }
625  if ( iTransfer >= G4int((*fPAIdEdxTable)(iPlace)->GetVectorLength()) )
626  {
627  iTransfer = (*fPAIdEdxTable)(iPlace)->GetVectorLength() - 1 ;
628  }
629  y1 = (*(*fPAIdEdxTable)(iPlace))(iTransfer-1) ;
630  y2 = (*(*fPAIdEdxTable)(iPlace))(iTransfer) ;
631  // G4cout<<"y1 = "<<y1<<"; "<<"y2 = "<<y2<<G4endl;
632  x1 = (*fPAIdEdxTable)(iPlace)->GetLowEdgeEnergy(iTransfer-1) ;
633  x2 = (*fPAIdEdxTable)(iPlace)->GetLowEdgeEnergy(iTransfer) ;
634  // G4cout<<"x1 = "<<x1<<"; "<<"x2 = "<<x2<<G4endl;
635 
636  if ( y1 == y2 ) dEdxCut = y2 ;
637  else
638  {
639  // if ( x1 == x2 ) dEdxCut = y1 + (y2 - y1)*G4UniformRand() ;
640  // if ( std::abs(x1-x2) <= eV ) dEdxCut = y1 + (y2 - y1)*G4UniformRand() ;
641  if ( std::abs(x1-x2) <= eV ) dEdxCut = y1 + (y2 - y1)*0.5 ;
642  else dEdxCut = y1 + (transferCut - x1)*(y2 - y1)/(x2 - x1) ;
643  }
644  // G4cout<<""<<dEdxCut<<G4endl;
645  return dEdxCut ;
646 }
647 
649 
650 G4double G4PAIPhotonModel::ComputeDEDXPerVolume(const G4Material*,
651  const G4ParticleDefinition* p,
652  G4double kineticEnergy,
653  G4double cutEnergy)
654 {
655  G4int iTkin,iPlace;
656  size_t jMat;
657 
658  //G4double cut = std::min(MaxSecondaryEnergy(p, kineticEnergy), cutEnergy);
659  G4double cut = cutEnergy;
660 
661  G4double particleMass = p->GetPDGMass();
662  G4double scaledTkin = kineticEnergy*proton_mass_c2/particleMass;
663  G4double charge = p->GetPDGCharge()/eplus;
664  G4double charge2 = charge*charge;
665  G4double dEdx = 0.;
666  const G4MaterialCutsCouple* matCC = CurrentCouple();
667 
668  for( jMat = 0 ;jMat < fMaterialCutsCoupleVector.size() ; ++jMat )
669  {
670  if( matCC == fMaterialCutsCoupleVector[jMat] ) break;
671  }
672  if(jMat == fMaterialCutsCoupleVector.size() && jMat > 0) jMat--;
673 
674  fPAIdEdxTable = fPAIdEdxBank[jMat];
675  fdEdxVector = fdEdxTable[jMat];
676  for(iTkin = 0 ; iTkin <= fTotBin ; iTkin++)
677  {
678  if(scaledTkin < fProtonEnergyVector->GetLowEdgeEnergy(iTkin)) break ;
679  }
680  iPlace = iTkin - 1;
681  if(iPlace < 0) iPlace = 0;
682  dEdx = charge2*( (*fdEdxVector)(iPlace) - GetdEdxCut(iPlace,cut) ) ;
683 
684  if( dEdx < 0.) dEdx = 0.;
685  return dEdx;
686 }
687 
689 
690 G4double G4PAIPhotonModel::CrossSectionPerVolume( const G4Material*,
691  const G4ParticleDefinition* p,
692  G4double kineticEnergy,
693  G4double cutEnergy,
694  G4double maxEnergy )
695 {
696  G4int iTkin,iPlace;
697  size_t jMat, jMatCC;
698  G4double tmax = std::min(MaxSecondaryEnergy(p, kineticEnergy), maxEnergy);
699  if(cutEnergy >= tmax) return 0.0;
700  G4double particleMass = p->GetPDGMass();
701  G4double scaledTkin = kineticEnergy*proton_mass_c2/particleMass;
702  G4double charge = p->GetPDGCharge();
703  G4double charge2 = charge*charge, cross, cross1, cross2;
704  G4double photon1, photon2, plasmon1, plasmon2;
705 
706  const G4MaterialCutsCouple* matCC = CurrentCouple();
707 
708  const G4ProductionCutsTable* theCoupleTable=
709  G4ProductionCutsTable::GetProductionCutsTable();
710 
711  size_t numOfCouples = theCoupleTable->GetTableSize();
712 
713  for (jMatCC = 0 ; jMatCC < numOfCouples ; jMatCC++ )
714  {
715  if( matCC == theCoupleTable->GetMaterialCutsCouple(jMatCC) ) break;
716  }
717  if( jMatCC == numOfCouples && jMatCC > 0 ) jMatCC--;
718 
719  const vector<G4double>* photonCutInKineticEnergy = theCoupleTable->
720  GetEnergyCutsVector(idxG4GammaCut);
721 
722  G4double photonCut = (*photonCutInKineticEnergy)[jMatCC] ;
723 
724  for( jMat = 0 ;jMat < fMaterialCutsCoupleVector.size() ; ++jMat )
725  {
726  if( matCC == fMaterialCutsCoupleVector[jMat] ) break;
727  }
728  if(jMat == fMaterialCutsCoupleVector.size() && jMat > 0) jMat--;
729 
730  fPAItransferTable = fPAIxscBank[jMat];
731  fPAIphotonTable = fPAIphotonBank[jMat];
732  fPAIplasmonTable = fPAIplasmonBank[jMat];
733 
734  for(iTkin = 0 ; iTkin <= fTotBin ; iTkin++)
735  {
736  if(scaledTkin < fProtonEnergyVector->GetLowEdgeEnergy(iTkin)) break ;
737  }
738  iPlace = iTkin - 1;
739  if(iPlace < 0) iPlace = 0;
740 
741  // G4cout<<"iPlace = "<<iPlace<<"; tmax = "
742  // <<tmax<<"; cutEnergy = "<<cutEnergy<<G4endl;
743  photon1 = GetdNdxPhotonCut(iPlace,tmax);
744  photon2 = GetdNdxPhotonCut(iPlace,photonCut);
745 
746  plasmon1 = GetdNdxPlasmonCut(iPlace,tmax);
747  plasmon2 = GetdNdxPlasmonCut(iPlace,cutEnergy);
748 
749  cross1 = photon1 + plasmon1;
750  // G4cout<<"cross1 = "<<cross1<<G4endl;
751  cross2 = photon2 + plasmon2;
752  // G4cout<<"cross2 = "<<cross2<<G4endl;
753  cross = (cross2 - cross1)*charge2;
754  // G4cout<<"cross = "<<cross<<G4endl;
755 
756  if( cross < 0. ) cross = 0.;
757  return cross;
758 }
759 
761 //
762 // It is analog of PostStepDoIt in terms of secondary electron or photon to
763 // be returned as G4Dynamicparticle*.
764 //
765 
766 void G4PAIPhotonModel::SampleSecondaries(std::vector<G4DynamicParticle*>* vdp,
767  const G4MaterialCutsCouple* matCC,
768  const G4DynamicParticle* dp,
769  G4double tmin,
770  G4double maxEnergy)
771 {
772  size_t jMat;
773  for( jMat = 0 ;jMat < fMaterialCutsCoupleVector.size() ; ++jMat )
774  {
775  if( matCC == fMaterialCutsCoupleVector[jMat] ) break;
776  }
777  if( jMat == fMaterialCutsCoupleVector.size() && jMat > 0 ) jMat--;
778 
779  fPAItransferTable = fPAIxscBank[jMat];
780  fPAIphotonTable = fPAIphotonBank[jMat];
781  fPAIplasmonTable = fPAIplasmonBank[jMat];
782 
783  fdNdxCutVector = fdNdxCutTable[jMat];
784  fdNdxCutPhotonVector = fdNdxCutPhotonTable[jMat];
785  fdNdxCutPlasmonVector = fdNdxCutPlasmonTable[jMat];
786 
787  G4double tmax = min(MaxSecondaryKinEnergy(dp), maxEnergy);
788  if( tmin >= tmax && fVerbose > 0)
789  {
790  G4cout<<"G4PAIPhotonModel::SampleSecondary: tmin >= tmax "<<G4endl;
791  }
792 
793  G4ThreeVector direction = dp->GetMomentumDirection();
794  G4double particleMass = dp->GetMass();
795  G4double kineticEnergy = dp->GetKineticEnergy();
796  G4double scaledTkin = kineticEnergy*fMass/particleMass;
797  G4double totalEnergy = kineticEnergy + particleMass;
798  G4double pSquare = kineticEnergy*(totalEnergy+particleMass);
799 
800  G4int iTkin;
801  for(iTkin=0;iTkin<=fTotBin;iTkin++)
802  {
803  if(scaledTkin < fProtonEnergyVector->GetLowEdgeEnergy(iTkin)) break ;
804  }
805  G4int iPlace = iTkin - 1 ;
806  if(iPlace < 0) iPlace = 0;
807 
808  G4double dNdxPhotonCut = (*fdNdxCutPhotonVector)(iPlace) ;
809  G4double dNdxPlasmonCut = (*fdNdxCutPlasmonVector)(iPlace) ;
810  G4double dNdxCut = dNdxPhotonCut + dNdxPlasmonCut;
811 
812  G4double ratio;
813  if (dNdxCut > 0.) ratio = dNdxPhotonCut/dNdxCut;
814  else ratio = 0.;
815 
816  if(ratio < G4UniformRand() ) // secondary e-
817  {
818  G4double deltaTkin = GetPostStepTransfer(fPAIplasmonTable, fdNdxCutPlasmonVector,
819  iPlace, scaledTkin);
820 
821 // G4cout<<"PAIPhotonModel PlasmonPostStepTransfer = "<<deltaTkin/keV<<" keV"<<G4endl ;
822 
823  if( deltaTkin <= 0. )
824  {
825  G4cout<<"G4PAIPhotonModel::SampleSecondary e- deltaTkin = "<<deltaTkin<<G4endl;
826  }
827  if( deltaTkin <= 0.) return;
828 
829  G4double deltaTotalMomentum = sqrt(deltaTkin*(deltaTkin + 2. * electron_mass_c2 ));
830  G4double totalMomentum = sqrt(pSquare);
831  G4double costheta = deltaTkin*(totalEnergy + electron_mass_c2)
832  /(deltaTotalMomentum * totalMomentum);
833 
834  if( costheta > 0.99999 ) costheta = 0.99999;
835  G4double sintheta = 0.0;
836  G4double sin2 = 1. - costheta*costheta;
837  if( sin2 > 0.) sintheta = sqrt(sin2);
838 
839  // direction of the delta electron
840 
841  G4double phi = twopi*G4UniformRand();
842  G4double dirx = sintheta*cos(phi), diry = sintheta*sin(phi), dirz = costheta;
843 
844  G4ThreeVector deltaDirection(dirx,diry,dirz);
845  deltaDirection.rotateUz(direction);
846 
847  // primary change
848 
849  kineticEnergy -= deltaTkin;
850  G4ThreeVector dir = totalMomentum*direction - deltaTotalMomentum*deltaDirection;
851  direction = dir.unit();
852  fParticleChange->SetProposedMomentumDirection(direction);
853 
854  // create G4DynamicParticle object for e- delta ray
855 
856  G4DynamicParticle* deltaRay = new G4DynamicParticle;
857  deltaRay->SetDefinition(G4Electron::Electron());
858  deltaRay->SetKineticEnergy( deltaTkin );
859  deltaRay->SetMomentumDirection(deltaDirection);
860  vdp->push_back(deltaRay);
861 
862  }
863  else // secondary 'Cherenkov' photon
864  {
865  G4double deltaTkin = GetPostStepTransfer(fPAIphotonTable, fdNdxCutPhotonVector,
866  iPlace,scaledTkin);
867 
868  // G4cout<<"PAIPhotonModel PhotonPostStepTransfer = "<<deltaTkin/keV<<" keV"<<G4endl ;
869 
870  if( deltaTkin <= 0. )
871  {
872  G4cout<<"G4PAIPhotonModel::SampleSecondary gamma deltaTkin = "<<deltaTkin<<G4endl;
873  }
874  if( deltaTkin <= 0.) return;
875 
876  G4double costheta = 0.; // G4UniformRand(); // VG: ??? for start only
877  G4double sintheta = sqrt((1.+costheta)*(1.-costheta));
878 
879  // direction of the 'Cherenkov' photon
880  G4double phi = twopi*G4UniformRand();
881  G4double dirx = sintheta*cos(phi), diry = sintheta*sin(phi), dirz = costheta;
882 
883  G4ThreeVector deltaDirection(dirx,diry,dirz);
884  deltaDirection.rotateUz(direction);
885 
886  // primary change
887  kineticEnergy -= deltaTkin;
888 
889  // create G4DynamicParticle object for photon ray
890 
891  G4DynamicParticle* photonRay = new G4DynamicParticle;
892  photonRay->SetDefinition( G4Gamma::Gamma() );
893  photonRay->SetKineticEnergy( deltaTkin );
894  photonRay->SetMomentumDirection(deltaDirection);
895 
896  vdp->push_back(photonRay);
897  }
898 
899  fParticleChange->SetProposedKineticEnergy(kineticEnergy);
900 }
901 
902 
904 //
905 // Returns post step PAI energy transfer > cut electron/photon energy according to passed
906 // scaled kinetic energy of particle
907 
908 G4double
909 G4PAIPhotonModel::GetPostStepTransfer( G4PhysicsTable* pTable,
910  G4PhysicsLogVector* pVector,
911  G4int iPlace, G4double scaledTkin )
912 {
913  // G4cout<<"G4PAIPhotonModel::GetPostStepTransfer"<<G4endl ;
914 
915  G4int iTkin = iPlace+1, iTransfer;
916  G4double transfer = 0.0, position, dNdxCut1, dNdxCut2, E1, E2, W1, W2, W ;
917 
918  dNdxCut1 = (*pVector)(iPlace) ;
919 
920  // G4cout<<"iPlace = "<<iPlace<<endl ;
921 
922  if(iTkin == fTotBin) // Fermi plato, try from left
923  {
924  position = dNdxCut1*G4UniformRand() ;
925 
926  for( iTransfer = 0;
927  iTransfer < G4int((*pTable)(iPlace)->GetVectorLength()); iTransfer++ )
928  {
929  if(position >= (*(*pTable)(iPlace))(iTransfer)) break ;
930  }
931  transfer = GetEnergyTransfer(pTable,iPlace,position,iTransfer);
932  }
933  else
934  {
935  dNdxCut2 = (*pVector)(iPlace+1) ;
936  if(iTkin == 0) // Tkin is too small, trying from right only
937  {
938  position = dNdxCut2*G4UniformRand() ;
939 
940  for( iTransfer = 0;
941  iTransfer < G4int((*pTable)(iPlace+1)->GetVectorLength()); iTransfer++ )
942  {
943  if(position >= (*(*pTable)(iPlace+1))(iTransfer)) break ;
944  }
945  transfer = GetEnergyTransfer(pTable,iPlace+1,position,iTransfer);
946  }
947  else // general case: Tkin between two vectors of the material
948  {
949  E1 = fProtonEnergyVector->GetLowEdgeEnergy(iTkin - 1) ;
950  E2 = fProtonEnergyVector->GetLowEdgeEnergy(iTkin) ;
951  W = 1.0/(E2 - E1) ;
952  W1 = (E2 - scaledTkin)*W ;
953  W2 = (scaledTkin - E1)*W ;
954 
955  position = ( dNdxCut1*W1 + dNdxCut2*W2 )*G4UniformRand() ;
956 
957  // G4cout<<position<<"\t" ;
958 
959  G4int iTrMax1, iTrMax2, iTrMax;
960 
961  iTrMax1 = G4int((*pTable)(iPlace)->GetVectorLength());
962  iTrMax2 = G4int((*pTable)(iPlace+1)->GetVectorLength());
963 
964  if (iTrMax1 >= iTrMax2) iTrMax = iTrMax2;
965  else iTrMax = iTrMax1;
966 
967  for( iTransfer = 0; iTransfer < iTrMax; iTransfer++ )
968  {
969  if( position >=
970  ( (*(*pTable)(iPlace))(iTransfer)*W1 +
971  (*(*pTable)(iPlace+1))(iTransfer)*W2) ) break ;
972  }
973  transfer = GetEnergyTransfer(pTable, iPlace, position, iTransfer);
974  }
975  }
976  // G4cout<<"PAIPhotonModel PostStepTransfer = "<<transfer/keV<<" keV"<<G4endl ;
977  if(transfer < 0.0 ) transfer = 0.0 ;
978  return transfer ;
979 }
980 
982 //
983 // Returns random PAI energy transfer according to passed
984 // indexes of particle
985 
986 G4double
987 G4PAIPhotonModel::GetEnergyTransfer( G4PhysicsTable* pTable, G4int iPlace,
988  G4double position, G4int iTransfer )
989 {
990  G4int iTransferMax;
991  G4double x1, x2, y1, y2, energyTransfer;
992 
993  if(iTransfer == 0)
994  {
995  energyTransfer = (*pTable)(iPlace)->GetLowEdgeEnergy(iTransfer);
996  }
997  else
998  {
999  iTransferMax = G4int((*pTable)(iPlace)->GetVectorLength());
1000 
1001  if ( iTransfer >= iTransferMax) iTransfer = iTransferMax - 1;
1002 
1003  y1 = (*(*pTable)(iPlace))(iTransfer-1);
1004  y2 = (*(*fPAItransferTable)(iPlace))(iTransfer);
1005 
1006  x1 = (*pTable)(iPlace)->GetLowEdgeEnergy(iTransfer-1);
1007  x2 = (*pTable)(iPlace)->GetLowEdgeEnergy(iTransfer);
1008 
1009  if ( x1 == x2 ) energyTransfer = x2;
1010  else
1011  {
1012  if ( y1 == y2 ) energyTransfer = x1 + (x2 - x1)*G4UniformRand();
1013  else
1014  {
1015  energyTransfer = x1 + (position - y1)*(x2 - x1)/(y2 - y1);
1016  }
1017  }
1018  }
1019  return energyTransfer;
1020 }
1021 
1023 //
1024 // Works like AlongStepDoIt method of process family
1025 
1026 
1027 
1028 
1029 G4double G4PAIPhotonModel::SampleFluctuations( const G4Material* material,
1030  const G4DynamicParticle* aParticle,
1031  G4double&,
1032  G4double& step,
1033  G4double&)
1034 {
1035  size_t jMat;
1036  for( jMat = 0 ;jMat < fMaterialCutsCoupleVector.size() ; ++jMat )
1037  {
1038  if( material == fMaterialCutsCoupleVector[jMat]->GetMaterial() ) break;
1039  }
1040  if(jMat == fMaterialCutsCoupleVector.size() && jMat > 0) jMat--;
1041 
1042  fPAItransferTable = fPAIxscBank[jMat];
1043  fPAIphotonTable = fPAIphotonBank[jMat];
1044  fPAIplasmonTable = fPAIplasmonBank[jMat];
1045 
1046  fdNdxCutVector = fdNdxCutTable[jMat];
1047  fdNdxCutPhotonVector = fdNdxCutPhotonTable[jMat];
1048  fdNdxCutPlasmonVector = fdNdxCutPlasmonTable[jMat];
1049 
1050  G4int iTkin, iPlace ;
1051 
1052  // G4cout<<"G4PAIPhotonModel::SampleFluctuations"<<G4endl ;
1053 
1054  G4double loss, photonLoss, plasmonLoss, charge2 ;
1055 
1056 
1057  G4double Tkin = aParticle->GetKineticEnergy() ;
1058  G4double MassRatio = proton_mass_c2/aParticle->GetDefinition()->GetPDGMass() ;
1059  G4double charge = aParticle->GetDefinition()->GetPDGCharge() ;
1060  charge2 = charge*charge ;
1061  G4double scaledTkin = Tkin*MassRatio ;
1062  G4double cof = step*charge2;
1063 
1064  for( iTkin = 0; iTkin <= fTotBin; iTkin++)
1065  {
1066  if(scaledTkin < fProtonEnergyVector->GetLowEdgeEnergy(iTkin)) break ;
1067  }
1068  iPlace = iTkin - 1 ;
1069  if( iPlace < 0 ) iPlace = 0;
1070 
1071  photonLoss = GetAlongStepTransfer(fPAIphotonTable,fdNdxCutPhotonVector,
1072 iPlace,scaledTkin,step,cof);
1073 
1074  // G4cout<<"PAIPhotonModel AlongStepPhotonLoss = "<<photonLoss/keV<<" keV"<<G4endl ;
1075 
1076  plasmonLoss = GetAlongStepTransfer(fPAIplasmonTable,fdNdxCutPlasmonVector,
1077 iPlace,scaledTkin,step,cof);
1078 
1079  // G4cout<<"PAIPhotonModel AlongStepPlasmonLoss = "<<plasmonLoss/keV<<" keV"<<G4endl ;
1080 
1081  loss = photonLoss + plasmonLoss;
1082 
1083  // G4cout<<"PAIPhotonModel AlongStepLoss = "<<loss/keV<<" keV"<<G4endl ;
1084 
1085  return loss;
1086 }
1087 
1089 //
1090 // Returns along step PAI energy transfer < cut electron/photon energy according to passed
1091 // scaled kinetic energy of particle and cof = step*charge*charge
1092 
1093 G4double
1094 G4PAIPhotonModel::GetAlongStepTransfer( G4PhysicsTable* pTable,
1095  G4PhysicsLogVector* pVector,
1096  G4int iPlace, G4double scaledTkin,G4double step,
1097  G4double cof )
1098 {
1099  G4int iTkin = iPlace + 1, iTransfer;
1100  G4double loss = 0., position, E1, E2, W1, W2, W, dNdxCut1, dNdxCut2, meanNumber;
1101  G4double lambda, stepDelta, stepSum=0. ;
1102  G4long numOfCollisions=0;
1103  G4bool numb = true;
1104 
1105  dNdxCut1 = (*pVector)(iPlace) ;
1106 
1107  // G4cout<<"iPlace = "<<iPlace<<endl ;
1108 
1109  if(iTkin == fTotBin) // Fermi plato, try from left
1110  {
1111  meanNumber = ((*(*pTable)(iPlace))(0) - dNdxCut1)*cof;
1112  if(meanNumber < 0.) meanNumber = 0. ;
1113  // numOfCollisions = RandPoisson::shoot(meanNumber) ;
1114  if( meanNumber > 0.) lambda = step/meanNumber;
1115  else lambda = DBL_MAX;
1116  while(numb)
1117  {
1118  stepDelta = CLHEP::RandExponential::shoot(lambda);
1119  stepSum += stepDelta;
1120  if(stepSum >= step) break;
1121  numOfCollisions++;
1122  }
1123 
1124  // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl ;
1125 
1126  while(numOfCollisions)
1127  {
1128  position = dNdxCut1+
1129  ((*(*pTable)(iPlace))(0) - dNdxCut1)*G4UniformRand() ;
1130 
1131  for( iTransfer = 0;
1132  iTransfer < G4int((*pTable)(iPlace)->GetVectorLength()); iTransfer++ )
1133  {
1134  if(position >= (*(*pTable)(iPlace))(iTransfer)) break ;
1135  }
1136  loss += GetEnergyTransfer(pTable,iPlace,position,iTransfer);
1137  numOfCollisions-- ;
1138  }
1139  }
1140  else
1141  {
1142  dNdxCut2 = (*pVector)(iPlace+1) ;
1143 
1144  if(iTkin == 0) // Tkin is too small, trying from right only
1145  {
1146  meanNumber = ((*(*pTable)(iPlace+1))(0) - dNdxCut2)*cof;
1147  if( meanNumber < 0. ) meanNumber = 0. ;
1148  // numOfCollisions = CLHEP::RandPoisson::shoot(meanNumber) ;
1149  if( meanNumber > 0.) lambda = step/meanNumber;
1150  else lambda = DBL_MAX;
1151  while(numb)
1152  {
1153  stepDelta = CLHEP::RandExponential::shoot(lambda);
1154  stepSum += stepDelta;
1155  if(stepSum >= step) break;
1156  numOfCollisions++;
1157  }
1158 
1159  // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl ;
1160 
1161  while(numOfCollisions)
1162  {
1163  position = dNdxCut2+
1164  ((*(*pTable)(iPlace+1))(0) - dNdxCut2)*G4UniformRand();
1165 
1166  for( iTransfer = 0;
1167  iTransfer < G4int((*pTable)(iPlace+1)->GetVectorLength()); iTransfer++ )
1168  {
1169  if(position >= (*(*pTable)(iPlace+1))(iTransfer)) break ;
1170  }
1171  loss += GetEnergyTransfer(pTable,iPlace+1,position,iTransfer);
1172  numOfCollisions-- ;
1173  }
1174  }
1175  else // general case: Tkin between two vectors of the material
1176  {
1177  E1 = fProtonEnergyVector->GetLowEdgeEnergy(iTkin - 1) ;
1178  E2 = fProtonEnergyVector->GetLowEdgeEnergy(iTkin) ;
1179  W = 1.0/(E2 - E1) ;
1180  W1 = (E2 - scaledTkin)*W ;
1181  W2 = (scaledTkin - E1)*W ;
1182 
1183  // G4cout<<"(*(*pTable)(iPlace))(0) = "<<
1184  // (*(*pTable)(iPlace))(0)<<G4endl ;
1185  // G4cout<<"(*(*pTable)(iPlace+1))(0) = "<<
1186  // (*(*pTable)(iPlace+1))(0)<<G4endl ;
1187 
1188  meanNumber=( ((*(*pTable)(iPlace))(0)-dNdxCut1)*W1 +
1189  ((*(*pTable)(iPlace+1))(0)-dNdxCut2)*W2 )*cof;
1190  if(meanNumber<0.0) meanNumber = 0.0;
1191  // numOfCollisions = CLHEP::RandPoisson::shoot(meanNumber) ;
1192  if( meanNumber > 0.) lambda = step/meanNumber;
1193  else lambda = DBL_MAX;
1194  while(numb)
1195  {
1196  stepDelta = CLHEP::RandExponential::shoot(lambda);
1197  stepSum += stepDelta;
1198  if(stepSum >= step) break;
1199  numOfCollisions++;
1200  }
1201 
1202  // G4cout<<"numOfCollisions = "<<numOfCollisions<<endl ;
1203 
1204  while(numOfCollisions)
1205  {
1206  position = dNdxCut1*W1 + dNdxCut2*W2 +
1207  ( ( (*(*pTable)(iPlace ))(0) - dNdxCut1)*W1 +
1208 
1209  ( (*(*pTable)(iPlace+1))(0) - dNdxCut2)*W2 )*G4UniformRand();
1210 
1211  // G4cout<<position<<"\t" ;
1212 
1213  for( iTransfer = 0;
1214  iTransfer < G4int((*pTable)(iPlace)->GetVectorLength()); iTransfer++ )
1215  {
1216  if( position >=
1217  ( (*(*pTable)(iPlace))(iTransfer)*W1 +
1218  (*(*pTable)(iPlace+1))(iTransfer)*W2) )
1219  {
1220  break ;
1221  }
1222  }
1223  // loss += (*pTable)(iPlace)->GetLowEdgeEnergy(iTransfer) ;
1224  loss += GetEnergyTransfer(pTable,iPlace,position,iTransfer);
1225  numOfCollisions-- ;
1226  }
1227  }
1228  }
1229 
1230  return loss ;
1231 
1232 }
1233 
1235 //
1236 // Returns the statistical estimation of the energy loss distribution variance
1237 //
1238 
1239 
1240 G4double G4PAIPhotonModel::Dispersion( const G4Material* material,
1241  const G4DynamicParticle* aParticle,
1242  G4double& tmax,
1243  G4double& step )
1244 {
1245  G4double loss, sumLoss=0., sumLoss2=0., sigma2, meanLoss=0.;
1246  for(G4int i = 0 ; i < fMeanNumber; i++)
1247  {
1248  loss = SampleFluctuations(material,aParticle,tmax,step,meanLoss);
1249  sumLoss += loss;
1250  sumLoss2 += loss*loss;
1251  }
1252  meanLoss = sumLoss/fMeanNumber;
1253  sigma2 = meanLoss*meanLoss + (sumLoss2-2*sumLoss*meanLoss)/fMeanNumber;
1254  return sigma2;
1255 }
1256 
1258 
1259 G4double G4PAIPhotonModel::MaxSecondaryEnergy( const G4ParticleDefinition* p,
1260  G4double kinEnergy)
1261 {
1262  G4double tmax = kinEnergy;
1263  if(p == fElectron) tmax *= 0.5;
1264  else if(p != fPositron) {
1265  G4double mass = p->GetPDGMass();
1266  G4double ratio= electron_mass_c2/mass;
1267  G4double gamma= kinEnergy/mass + 1.0;
1268  tmax = 2.0*electron_mass_c2*(gamma*gamma - 1.) /
1269  (1. + 2.0*gamma*ratio + ratio*ratio);
1270  }
1271  return tmax;
1272 }
1273 
1275 
1276 void G4PAIPhotonModel::DefineForRegion(const G4Region* r)
1277 {
1278  fPAIRegionVector.push_back(r);
1279 }
1280 
1281 
1282 //
1283 //
1285 
1286 
1287 
1288 
1289 
1290 
Generated on Sat May 25 2013 14:33:37 for Geant4 by   doxygen 1.8.4