Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4PAIPhotonModel Class Reference

#include <G4PAIPhotonModel.hh>

+ Inheritance diagram for G4PAIPhotonModel:

Public Member Functions

 G4PAIPhotonModel (const G4ParticleDefinition *p=0, const G4String &nam="PAIPhoton")
 
virtual ~G4PAIPhotonModel ()
 
virtual void Initialise (const G4ParticleDefinition *, const G4DataVector &)
 
virtual void InitialiseMe (const G4ParticleDefinition *)
 
virtual G4double ComputeDEDXPerVolume (const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy)
 
virtual G4double CrossSectionPerVolume (const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy, G4double maxEnergy)
 
virtual void SampleSecondaries (std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
 
virtual G4double SampleFluctuations (const G4Material *, const G4DynamicParticle *, G4double &, G4double &, G4double &)
 
virtual G4double Dispersion (const G4Material *, const G4DynamicParticle *, G4double &, G4double &)
 
void DefineForRegion (const G4Region *r)
 
void ComputeSandiaPhotoAbsCof ()
 
void BuildPAIonisationTable ()
 
void BuildLambdaVector (const G4MaterialCutsCouple *matCutsCouple)
 
G4double GetdNdxCut (G4int iPlace, G4double transferCut)
 
G4double GetdNdxPhotonCut (G4int iPlace, G4double transferCut)
 
G4double GetdNdxPlasmonCut (G4int iPlace, G4double transferCut)
 
G4double GetdEdxCut (G4int iPlace, G4double transferCut)
 
G4double GetPostStepTransfer (G4PhysicsTable *, G4PhysicsLogVector *, G4int iPlace, G4double scaledTkin)
 
G4double GetAlongStepTransfer (G4PhysicsTable *, G4PhysicsLogVector *, G4int iPlace, G4double scaledTkin, G4double step, G4double cof)
 
G4double GetEnergyTransfer (G4PhysicsTable *, G4int iPlace, G4double position, G4int iTransfer)
 
- Public Member Functions inherited from G4VEmModel
 G4VEmModel (const G4String &nam)
 
virtual ~G4VEmModel ()
 
virtual void Initialise (const G4ParticleDefinition *, const G4DataVector &)=0
 
virtual void SampleSecondaries (std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin=0.0, G4double tmax=DBL_MAX)=0
 
virtual G4double ComputeDEDXPerVolume (const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=DBL_MAX)
 
virtual G4double CrossSectionPerVolume (const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
virtual G4double ComputeCrossSectionPerAtom (const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0., G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
virtual G4double ChargeSquareRatio (const G4Track &)
 
virtual G4double GetChargeSquareRatio (const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
 
virtual G4double GetParticleCharge (const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
 
virtual void StartTracking (G4Track *)
 
virtual void CorrectionsAlongStep (const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double &eloss, G4double &niel, G4double length)
 
virtual G4double Value (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy)
 
virtual G4double MinPrimaryEnergy (const G4Material *, const G4ParticleDefinition *)
 
virtual void SetupForMaterial (const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
 
virtual void DefineForRegion (const G4Region *)
 
void InitialiseElementSelectors (const G4ParticleDefinition *, const G4DataVector &)
 
G4double ComputeDEDX (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=DBL_MAX)
 
G4double CrossSection (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
G4double ComputeMeanFreePath (const G4ParticleDefinition *, G4double kineticEnergy, const G4Material *, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
G4double ComputeCrossSectionPerAtom (const G4ParticleDefinition *, const G4Element *, G4double kinEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
G4int SelectIsotopeNumber (const G4Element *)
 
const G4ElementSelectRandomAtom (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
const G4ElementSelectRandomAtom (const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
void SetParticleChange (G4VParticleChange *, G4VEmFluctuationModel *f=0)
 
void SetCrossSectionTable (G4PhysicsTable *)
 
G4PhysicsTableGetCrossSectionTable ()
 
G4VEmFluctuationModelGetModelOfFluctuations ()
 
G4VEmAngularDistributionGetAngularDistribution ()
 
void SetAngularDistribution (G4VEmAngularDistribution *)
 
G4double HighEnergyLimit () const
 
G4double LowEnergyLimit () const
 
G4double HighEnergyActivationLimit () const
 
G4double LowEnergyActivationLimit () const
 
G4double PolarAngleLimit () const
 
G4double SecondaryThreshold () const
 
G4bool LPMFlag () const
 
G4bool DeexcitationFlag () const
 
G4bool ForceBuildTableFlag () const
 
void SetHighEnergyLimit (G4double)
 
void SetLowEnergyLimit (G4double)
 
void SetActivationHighEnergyLimit (G4double)
 
void SetActivationLowEnergyLimit (G4double)
 
G4bool IsActive (G4double kinEnergy)
 
void SetPolarAngleLimit (G4double)
 
void SetSecondaryThreshold (G4double)
 
void SetLPMFlag (G4bool val)
 
void SetDeexcitationFlag (G4bool val)
 
void ForceBuildTable (G4bool val)
 
G4double MaxSecondaryKinEnergy (const G4DynamicParticle *dynParticle)
 
const G4StringGetName () const
 
void SetCurrentCouple (const G4MaterialCutsCouple *)
 
const G4ElementGetCurrentElement () const
 
- Public Member Functions inherited from G4VEmFluctuationModel
 G4VEmFluctuationModel (const G4String &nam)
 
virtual ~G4VEmFluctuationModel ()
 
virtual G4double SampleFluctuations (const G4Material *, const G4DynamicParticle *, G4double &tmax, G4double &length, G4double &meanLoss)=0
 
virtual G4double Dispersion (const G4Material *, const G4DynamicParticle *, G4double &tmax, G4double &length)=0
 
virtual void InitialiseMe (const G4ParticleDefinition *)
 
virtual void SetParticleAndCharge (const G4ParticleDefinition *, G4double q2)
 
G4String GetName () const
 

Protected Member Functions

G4double MaxSecondaryEnergy (const G4ParticleDefinition *, G4double kinEnergy)
 
- Protected Member Functions inherited from G4VEmModel
G4ParticleChangeForLossGetParticleChangeForLoss ()
 
G4ParticleChangeForGammaGetParticleChangeForGamma ()
 
virtual G4double MaxSecondaryEnergy (const G4ParticleDefinition *, G4double kineticEnergy)
 
const G4MaterialCutsCoupleCurrentCouple () const
 
void SetCurrentElement (const G4Element *)
 

Additional Inherited Members

- Protected Attributes inherited from G4VEmModel
G4VParticleChangepParticleChange
 
G4PhysicsTablexSectionTable
 
const std::vector< G4double > * theDensityFactor
 
const std::vector< G4int > * theDensityIdx
 

Detailed Description

Definition at line 66 of file G4PAIPhotonModel.hh.

Constructor & Destructor Documentation

◆ G4PAIPhotonModel()

G4PAIPhotonModel::G4PAIPhotonModel ( const G4ParticleDefinition p = 0,
const G4String nam = "PAIPhoton" 
)

Definition at line 75 of file G4PAIPhotonModel.cc.

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}
static G4Electron * Electron()
Definition: G4Electron.cc:94
static G4Positron * Positron()
Definition: G4Positron.cc:94

◆ ~G4PAIPhotonModel()

G4PAIPhotonModel::~G4PAIPhotonModel ( )
virtual

Definition at line 121 of file G4PAIPhotonModel.cc.

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}
int G4int
Definition: G4Types.hh:66
void clearAndDestroy()

Member Function Documentation

◆ BuildLambdaVector()

void G4PAIPhotonModel::BuildLambdaVector ( const G4MaterialCutsCouple matCutsCouple)

Definition at line 401 of file G4PAIPhotonModel.cc.

402{
403 G4int i ;
404 G4double dNdxCut,dNdxPhotonCut,dNdxPlasmonCut, lambda;
407
408 const G4ProductionCutsTable* theCoupleTable=
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}
@ idxG4ElectronCut
@ idxG4GammaCut
double G4double
Definition: G4Types.hh:64
#define G4endl
Definition: G4ios.hh:52
G4DLLIMPORT std::ostream G4cout
G4double GetSurfaceTolerance() const
static G4GeometryTolerance * GetInstance()
G4double GetdNdxPlasmonCut(G4int iPlace, G4double transferCut)
G4double GetdNdxPhotonCut(G4int iPlace, G4double transferCut)
void PutValue(size_t index, G4double theValue)
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
static G4ProductionCutsTable * GetProductionCutsTable()
#define DBL_MIN
Definition: templates.hh:75
#define DBL_MAX
Definition: templates.hh:83

Referenced by Initialise().

◆ BuildPAIonisationTable()

void G4PAIPhotonModel::BuildPAIonisationTable ( )

Definition at line 285 of file G4PAIPhotonModel.cc.

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}
G4double MaxSecondaryEnergy(const G4ParticleDefinition *, G4double kinEnergy)
void PutValue(size_t binNumber, G4double binValue, G4double dataValue)
void insertAt(size_t, G4PhysicsVector *)
virtual G4double GetLowEdgeEnergy(size_t binNumber) const

Referenced by Initialise().

◆ ComputeDEDXPerVolume()

G4double G4PAIPhotonModel::ComputeDEDXPerVolume ( const G4Material ,
const G4ParticleDefinition p,
G4double  kineticEnergy,
G4double  cutEnergy 
)
virtual

Reimplemented from G4VEmModel.

Definition at line 650 of file G4PAIPhotonModel.cc.

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}
G4double GetdEdxCut(G4int iPlace, G4double transferCut)
G4double GetPDGCharge() const
const G4MaterialCutsCouple * CurrentCouple() const
Definition: G4VEmModel.hh:377

◆ ComputeSandiaPhotoAbsCof()

void G4PAIPhotonModel::ComputeSandiaPhotoAbsCof ( )

Definition at line 237 of file G4PAIPhotonModel.cc.

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}
std::vector< G4Material * > G4MaterialTable
static const G4MaterialTable * GetMaterialTable()
Definition: G4Material.cc:562

Referenced by Initialise().

◆ CrossSectionPerVolume()

G4double G4PAIPhotonModel::CrossSectionPerVolume ( const G4Material ,
const G4ParticleDefinition p,
G4double  kineticEnergy,
G4double  cutEnergy,
G4double  maxEnergy 
)
virtual

Reimplemented from G4VEmModel.

Definition at line 690 of file G4PAIPhotonModel.cc.

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=
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}

◆ DefineForRegion()

void G4PAIPhotonModel::DefineForRegion ( const G4Region r)
virtual

Reimplemented from G4VEmModel.

Definition at line 1276 of file G4PAIPhotonModel.cc.

1277{
1278 fPAIRegionVector.push_back(r);
1279}

◆ Dispersion()

G4double G4PAIPhotonModel::Dispersion ( const G4Material material,
const G4DynamicParticle aParticle,
G4double tmax,
G4double step 
)
virtual

Implements G4VEmFluctuationModel.

Definition at line 1240 of file G4PAIPhotonModel.cc.

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}
virtual G4double SampleFluctuations(const G4Material *, const G4DynamicParticle *, G4double &, G4double &, G4double &)

◆ GetAlongStepTransfer()

G4double G4PAIPhotonModel::GetAlongStepTransfer ( G4PhysicsTable pTable,
G4PhysicsLogVector pVector,
G4int  iPlace,
G4double  scaledTkin,
G4double  step,
G4double  cof 
)

Definition at line 1094 of file G4PAIPhotonModel.cc.

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}
long G4long
Definition: G4Types.hh:68
bool G4bool
Definition: G4Types.hh:67
#define G4UniformRand()
Definition: Randomize.hh:53
G4double GetEnergyTransfer(G4PhysicsTable *, G4int iPlace, G4double position, G4int iTransfer)

Referenced by SampleFluctuations().

◆ GetdEdxCut()

G4double G4PAIPhotonModel::GetdEdxCut ( G4int  iPlace,
G4double  transferCut 
)

Definition at line 609 of file G4PAIPhotonModel.cc.

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}

Referenced by ComputeDEDXPerVolume().

◆ GetdNdxCut()

G4double G4PAIPhotonModel::GetdNdxCut ( G4int  iPlace,
G4double  transferCut 
)

Definition at line 476 of file G4PAIPhotonModel.cc.

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}

◆ GetdNdxPhotonCut()

G4double G4PAIPhotonModel::GetdNdxPhotonCut ( G4int  iPlace,
G4double  transferCut 
)

Definition at line 520 of file G4PAIPhotonModel.cc.

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}

Referenced by BuildLambdaVector(), and CrossSectionPerVolume().

◆ GetdNdxPlasmonCut()

G4double G4PAIPhotonModel::GetdNdxPlasmonCut ( G4int  iPlace,
G4double  transferCut 
)

Definition at line 564 of file G4PAIPhotonModel.cc.

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}

Referenced by BuildLambdaVector(), and CrossSectionPerVolume().

◆ GetEnergyTransfer()

G4double G4PAIPhotonModel::GetEnergyTransfer ( G4PhysicsTable pTable,
G4int  iPlace,
G4double  position,
G4int  iTransfer 
)

Definition at line 987 of file G4PAIPhotonModel.cc.

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}

Referenced by GetAlongStepTransfer(), and GetPostStepTransfer().

◆ GetPostStepTransfer()

G4double G4PAIPhotonModel::GetPostStepTransfer ( G4PhysicsTable pTable,
G4PhysicsLogVector pVector,
G4int  iPlace,
G4double  scaledTkin 
)

Definition at line 909 of file G4PAIPhotonModel.cc.

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}

Referenced by SampleSecondaries().

◆ Initialise()

void G4PAIPhotonModel::Initialise ( const G4ParticleDefinition p,
const G4DataVector  
)
virtual

Implements G4VEmModel.

Definition at line 169 of file G4PAIPhotonModel.cc.

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 =
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
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}
static size_t GetNumberOfMaterials()
Definition: G4Material.cc:569
void BuildLambdaVector(const G4MaterialCutsCouple *matCutsCouple)
G4ProductionCuts * GetProductionCuts() const
size_t GetNumberOfMaterials() const
std::vector< G4Material * >::const_iterator GetMaterialIterator() const
G4ParticleChangeForLoss * GetParticleChangeForLoss()
Definition: G4VEmModel.cc:95

◆ InitialiseMe()

void G4PAIPhotonModel::InitialiseMe ( const G4ParticleDefinition )
virtual

Reimplemented from G4VEmFluctuationModel.

Definition at line 232 of file G4PAIPhotonModel.cc.

233{}

◆ MaxSecondaryEnergy()

G4double G4PAIPhotonModel::MaxSecondaryEnergy ( const G4ParticleDefinition p,
G4double  kinEnergy 
)
protectedvirtual

Reimplemented from G4VEmModel.

Definition at line 1259 of file G4PAIPhotonModel.cc.

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}

Referenced by BuildPAIonisationTable(), and CrossSectionPerVolume().

◆ SampleFluctuations()

G4double G4PAIPhotonModel::SampleFluctuations ( const G4Material material,
const G4DynamicParticle aParticle,
G4double ,
G4double step,
G4double  
)
virtual

Implements G4VEmFluctuationModel.

Definition at line 1029 of file G4PAIPhotonModel.cc.

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,
1072iPlace,scaledTkin,step,cof);
1073
1074 // G4cout<<"PAIPhotonModel AlongStepPhotonLoss = "<<photonLoss/keV<<" keV"<<G4endl ;
1075
1076 plasmonLoss = GetAlongStepTransfer(fPAIplasmonTable,fdNdxCutPlasmonVector,
1077iPlace,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}
G4ParticleDefinition * GetDefinition() const
G4double GetKineticEnergy() const
G4double GetAlongStepTransfer(G4PhysicsTable *, G4PhysicsLogVector *, G4int iPlace, G4double scaledTkin, G4double step, G4double cof)

Referenced by Dispersion().

◆ SampleSecondaries()

void G4PAIPhotonModel::SampleSecondaries ( std::vector< G4DynamicParticle * > *  vdp,
const G4MaterialCutsCouple matCC,
const G4DynamicParticle dp,
G4double  tmin,
G4double  maxEnergy 
)
virtual

Implements G4VEmModel.

Definition at line 766 of file G4PAIPhotonModel.cc.

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;
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}
Hep3Vector unit() const
G4double GetMass() const
void SetMomentumDirection(const G4ThreeVector &aDirection)
const G4ThreeVector & GetMomentumDirection() const
void SetDefinition(const G4ParticleDefinition *aParticleDefinition)
void SetKineticEnergy(G4double aEnergy)
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
G4double GetPostStepTransfer(G4PhysicsTable *, G4PhysicsLogVector *, G4int iPlace, G4double scaledTkin)
void SetProposedKineticEnergy(G4double proposedKinEnergy)
void SetProposedMomentumDirection(const G4ThreeVector &dir)
G4double MaxSecondaryKinEnergy(const G4DynamicParticle *dynParticle)
Definition: G4VEmModel.hh:399

The documentation for this class was generated from the following files: