Geant4.10
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Public Member Functions | Protected Member Functions
G4PAIPhotonModel Class Reference

#include <G4PAIPhotonModel.hh>

Inheritance diagram for G4PAIPhotonModel:
G4VEmModel G4VEmFluctuationModel

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 void InitTest (const G4ParticleDefinition *, G4MaterialCutsCouple *, G4double, G4double)
 
virtual G4double ComputeDEDXPerVolume (const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy)
 
virtual G4double CrossSectionPerVolume (const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy, G4double maxEnergy)
 
G4double GetXscPerVolume (const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double photonCut, G4double cutEnergy, G4double maxEnergy)
 
virtual void SampleSecondaries (std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
 
G4double TestSecondaries (G4MaterialCutsCouple *, G4DynamicParticle *, G4double tmin, G4double maxEnergy)
 
virtual G4double SampleFluctuations (const G4MaterialCutsCouple *, 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)
 
void BuildLambdaVector (const G4MaterialCutsCouple *matCutsCouple, G4double, G4double)
 
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 InitialiseLocal (const G4ParticleDefinition *, G4VEmModel *masterModel)
 
virtual void InitialiseForMaterial (const G4ParticleDefinition *, const G4Material *)
 
virtual void InitialiseForElement (const G4ParticleDefinition *, G4int Z)
 
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 *, G4double cut=0.0)
 
virtual G4double MinEnergyCut (const G4ParticleDefinition *, const G4MaterialCutsCouple *)
 
virtual void SetupForMaterial (const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
 
void InitialiseElementSelectors (const G4ParticleDefinition *, const G4DataVector &)
 
std::vector
< G4EmElementSelector * > * 		GetElementSelectors ()
 
void SetElementSelectors (std::vector< G4EmElementSelector * > *)
 
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 G4Element * SelectRandomAtom (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
const G4Element * SelectRandomAtom (const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
G4int SelectRandomAtomNumber (const G4Material *)
 
void SetParticleChange (G4VParticleChange *, G4VEmFluctuationModel *f=0)
 
void SetCrossSectionTable (G4PhysicsTable *, G4bool isLocal)
 
G4ElementData * GetElementData ()
 
G4PhysicsTable * GetCrossSectionTable ()
 
G4VEmFluctuationModel * GetModelOfFluctuations ()
 
G4VEmAngularDistribution * GetAngularDistribution ()
 
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
 
G4bool UseAngularGeneratorFlag () const
 
void SetAngularGeneratorFlag (G4bool)
 
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 SetForceBuildTable (G4bool val)
 
void SetMasterThread (G4bool val)
 
G4bool IsMaster () const
 
G4double MaxSecondaryKinEnergy (const G4DynamicParticle *dynParticle)
 
const G4String & GetName () const
 
void SetCurrentCouple (const G4MaterialCutsCouple *)
 
const G4Element * GetCurrentElement () const
 
- Public Member Functions inherited from G4VEmFluctuationModel
 G4VEmFluctuationModel (const G4String &nam)
 
virtual ~G4VEmFluctuationModel ()
 
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
G4ParticleChangeForLoss * GetParticleChangeForLoss ()
 
G4ParticleChangeForGamma * GetParticleChangeForGamma ()
 
const G4MaterialCutsCouple * CurrentCouple () const
 
void SetCurrentElement (const G4Element *)
 

Additional Inherited Members

- Protected Attributes inherited from G4VEmModel
G4ElementData * fElementData
 
G4VParticleChange * pParticleChange
 
G4PhysicsTable * xSectionTable
 
const std::vector< G4double > * theDensityFactor
 
const std::vector< G4int > * theDensityIdx
 
size_t idxTable
 

Detailed Description

Definition at line 69 of file G4PAIPhotonModel.hh.

Constructor & Destructor Documentation

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

Definition at line 80 of file G4PAIPhotonModel.cc.

References G4Electron::Electron(), G4Gamma::Gamma(), and G4Positron::Positron().

81  : G4VEmModel(nam),G4VEmFluctuationModel(nam),
82  fLowestKineticEnergy(10.0*keV),
83  fHighestKineticEnergy(100.*TeV),
84  fTotBin(200),
85  fMeanNumber(20),
86  fParticle(0),
87  fHighKinEnergy(100.*TeV),
88  fLowKinEnergy(2.0*MeV),
89  fTwoln10(2.0*log(10.0)),
90  fBg2lim(0.0169),
91  fTaulim(8.4146e-3)
92 {
93  fVerbose = 0;
94  fGamma = G4Gamma::Gamma();
95  fElectron = G4Electron::Electron();
96  fPositron = G4Positron::Positron();
97 
98  fProtonEnergyVector = new G4PhysicsLogVector(fLowestKineticEnergy,
99  fHighestKineticEnergy,
100  fTotBin);
101  fPAItransferTable = 0;
102  fPAIphotonTable = 0;
103  fPAIplasmonTable = 0;
104 
105  fPAIdEdxTable = 0;
106  fSandiaPhotoAbsCof = 0;
107  fdEdxVector = 0;
108 
109  fLambdaVector = 0;
110  fdNdxCutVector = 0;
111  fdNdxCutPhotonVector = 0;
112  fdNdxCutPlasmonVector = 0;
113 
114  fSandiaIntervalNumber = 0;
115  fMatIndex = 0;
116  fCutCouple = 0;
117  fMaterial = 0;
118 
119  fParticleChange = 0;
120 
121  if(p) { SetParticle(p); }
122  else { SetParticle(fElectron); }
123 
124  isInitialised = false;
125 }
G4VEmModel::G4VEmModel
G4VEmModel(const G4String &nam)
Definition: G4VEmModel.cc:65
G4Gamma::Gamma
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
G4Positron::Positron
static G4Positron * Positron()
Definition: G4Positron.cc:94
G4Electron::Electron
static G4Electron * Electron()
Definition: G4Electron.cc:94
G4VEmFluctuationModel::G4VEmFluctuationModel
G4VEmFluctuationModel(const G4String &nam)
Definition: G4VEmFluctuationModel.cc:55
G4PAIPhotonModel::~G4PAIPhotonModel ( )
virtual

Definition at line 129 of file G4PAIPhotonModel.cc.

130 {
131  // if(fdEdxVector) delete fdEdxVector;
132  // if ( fLambdaVector) delete fLambdaVector;
133  // if ( fdNdxCutVector) delete fdNdxCutVector;
134 
135  if( fPAItransferTable )
136  {
137  delete fPAItransferTable;
138  }
139  if( fPAIphotonTable )
140  {
141  delete fPAIphotonTable;
142  }
143  if( fPAIplasmonTable )
144  {
145  delete fPAIplasmonTable;
146  }
147 }

Member Function Documentation

void G4PAIPhotonModel::BuildLambdaVector ( const G4MaterialCutsCouple *  matCutsCouple)

Definition at line 470 of file G4PAIPhotonModel.cc.

References DBL_MAX, DBL_MIN, G4cout, G4endl, GetdNdxPhotonCut(), GetdNdxPlasmonCut(), G4ProductionCutsTable::GetEnergyCutsVector(), G4GeometryTolerance::GetInstance(), G4ProductionCutsTable::GetMaterialCutsCouple(), G4ProductionCutsTable::GetProductionCutsTable(), G4GeometryTolerance::GetSurfaceTolerance(), G4ProductionCutsTable::GetTableSize(), idxG4ElectronCut, idxG4GammaCut, python.hepunit::keV, G4InuclParticleNames::lambda, and G4PhysicsVector::PutValue().

Referenced by Initialise(), and InitTest().

471 {
472  G4int i;
473  G4double dNdxCut,dNdxPhotonCut,dNdxPlasmonCut, lambda, deltaCutInKineticEnergyNow, photonCutInKineticEnergyNow;
474  G4double kCarTolerance = G4GeometryTolerance::GetInstance()
475  ->GetSurfaceTolerance();
476 
477  G4ProductionCutsTable* theCoupleTable=
478  G4ProductionCutsTable::GetProductionCutsTable();
479 
480  // G4EmCalculator converter;
481 
482  // const G4Material* material = matCutsCouple->GetMaterial();
483 
484  // G4double rangeGamma = matCutsCouple->GetProductionCuts()->GetProductionCut(0);
485  // G4double rangeElectron = matCutsCouple->GetProductionCuts()->GetProductionCut(1);
486 
487  size_t numOfCouples = theCoupleTable->GetTableSize();
488  size_t jMatCC;
489 
490  for (jMatCC = 0; jMatCC < numOfCouples; jMatCC++ )
491  {
492  if( matCutsCouple == theCoupleTable->GetMaterialCutsCouple(jMatCC) ) break;
493  }
494  if( jMatCC == numOfCouples && jMatCC > 0 ) jMatCC--;
495 
496  const vector<G4double>* deltaCutInKineticEnergy = theCoupleTable->GetEnergyCutsVector(idxG4ElectronCut);
497  const vector<G4double>* photonCutInKineticEnergy = theCoupleTable->GetEnergyCutsVector(idxG4GammaCut);
498 
499  if (fLambdaVector) delete fLambdaVector;
500  if (fdNdxCutVector) delete fdNdxCutVector;
501  if (fdNdxCutPhotonVector) delete fdNdxCutPhotonVector;
502  if (fdNdxCutPlasmonVector) delete fdNdxCutPlasmonVector;
503 
504  fLambdaVector = new G4PhysicsLogVector( fLowestKineticEnergy,
505  fHighestKineticEnergy,
506  fTotBin );
507  fdNdxCutVector = new G4PhysicsLogVector( fLowestKineticEnergy,
508  fHighestKineticEnergy,
509  fTotBin );
510  fdNdxCutPhotonVector = new G4PhysicsLogVector( fLowestKineticEnergy,
511  fHighestKineticEnergy,
512  fTotBin );
513  fdNdxCutPlasmonVector = new G4PhysicsLogVector( fLowestKineticEnergy,
514  fHighestKineticEnergy,
515  fTotBin );
516 
517  deltaCutInKineticEnergyNow = (*deltaCutInKineticEnergy)[jMatCC];
518  photonCutInKineticEnergyNow = (*photonCutInKineticEnergy)[jMatCC];
519 
520  // deltaCutInKineticEnergyNow = theCoupleTable->ConvertRangeToEnergy(fElectron,material,rangeElectron);
521  // photonCutInKineticEnergyNow = theCoupleTable->ConvertRangeToEnergy(fGamma,material,rangeGamma);
522 
523  // photonCutInKineticEnergyNow = converter.GetKinEnergy(rangeGamma, fGamma,material);
524  // deltaCutInKineticEnergyNow = converter.GetKinEnergy(rangeElectron, fElectron,material);
525 
526  if(fVerbose > 0)
527  {
528  G4cout<<"PAIPhotonModel deltaCutInKineticEnergyNow = "
529  <<deltaCutInKineticEnergyNow/keV<<" keV"<<G4endl;
530  G4cout<<"PAIPhotonModel photonCutInKineticEnergyNow = "
531  <<photonCutInKineticEnergyNow/keV<<" keV"<<G4endl;
532  }
533  for ( i = 0; i <= fTotBin; i++ )
534  {
535  dNdxPhotonCut = GetdNdxPhotonCut(i,photonCutInKineticEnergyNow);
536  dNdxPlasmonCut = GetdNdxPlasmonCut(i,deltaCutInKineticEnergyNow);
537 
538  dNdxCut = dNdxPhotonCut + dNdxPlasmonCut;
539  lambda = dNdxCut <= DBL_MIN ? DBL_MAX: 1.0/dNdxCut;
540 
541  if (lambda <= 1000*kCarTolerance) lambda = 1000*kCarTolerance; // Mmm ???
542 
543  fLambdaVector->PutValue(i, lambda);
544 
545  fdNdxCutVector->PutValue(i, dNdxCut);
546  fdNdxCutPhotonVector->PutValue(i, dNdxPhotonCut);
547  fdNdxCutPlasmonVector->PutValue(i, dNdxPlasmonCut);
548  }
549 }
G4ProductionCutsTable::GetEnergyCutsVector
const std::vector< G4double > * GetEnergyCutsVector(size_t pcIdx) const
Definition: G4ProductionCutsTable.hh:250
G4GeometryTolerance::GetSurfaceTolerance
G4double GetSurfaceTolerance() const
Definition: G4GeometryTolerance.cc:87
idxG4GammaCut
Definition: G4ProductionCuts.hh:53
G4PAIPhotonModel::GetdNdxPlasmonCut
G4double GetdNdxPlasmonCut(G4int iPlace, G4double transferCut)
Definition: G4PAIPhotonModel.cc:723
G4int
int G4int
Definition: G4Types.hh:78
idxG4ElectronCut
Definition: G4ProductionCuts.hh:54
G4cout
G4GLOB_DLL std::ostream G4cout
G4PAIPhotonModel::GetdNdxPhotonCut
G4double GetdNdxPhotonCut(G4int iPlace, G4double transferCut)
Definition: G4PAIPhotonModel.cc:678
G4PhysicsVector::PutValue
void PutValue(size_t index, G4double theValue)
G4ProductionCutsTable::GetProductionCutsTable
static G4ProductionCutsTable * GetProductionCutsTable()
Definition: G4ProductionCutsTable.cc:63
G4ProductionCutsTable::GetMaterialCutsCouple
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
Definition: G4ProductionCutsTable.hh:262
DBL_MIN
#define DBL_MIN
Definition: templates.hh:75
G4endl
#define G4endl
Definition: G4ios.hh:61
G4double
double G4double
Definition: G4Types.hh:76
DBL_MAX
#define DBL_MAX
Definition: templates.hh:83
G4GeometryTolerance::GetInstance
static G4GeometryTolerance * GetInstance()
Definition: G4GeometryTolerance.cc:74
void G4PAIPhotonModel::BuildLambdaVector ( const G4MaterialCutsCouple *  matCutsCouple,
G4double  photEnergy,
G4double  eTkin 
)

Definition at line 558 of file G4PAIPhotonModel.cc.

References DBL_MAX, DBL_MIN, G4cout, G4endl, GetdNdxPhotonCut(), GetdNdxPlasmonCut(), G4GeometryTolerance::GetInstance(), G4ProductionCutsTable::GetMaterialCutsCouple(), G4ProductionCutsTable::GetProductionCutsTable(), G4GeometryTolerance::GetSurfaceTolerance(), G4ProductionCutsTable::GetTableSize(), python.hepunit::keV, G4InuclParticleNames::lambda, and G4PhysicsVector::PutValue().

559 {
560  G4int i;
561  G4double dNdxCut,dNdxPhotonCut,dNdxPlasmonCut, lambda, deltaCutInKineticEnergyNow, photonCutInKineticEnergyNow;
562  G4double kCarTolerance = G4GeometryTolerance::GetInstance()
563  ->GetSurfaceTolerance();
564 
565  G4ProductionCutsTable* theCoupleTable=
566  G4ProductionCutsTable::GetProductionCutsTable();
567 
568 
569  // const G4Material* material = matCutsCouple->GetMaterial();
570 
571 
572  size_t numOfCouples = theCoupleTable->GetTableSize();
573  size_t jMatCC;
574 
575  for (jMatCC = 0; jMatCC < numOfCouples; jMatCC++ )
576  {
577  if( matCutsCouple == theCoupleTable->GetMaterialCutsCouple(jMatCC) ) break;
578  }
579  if( jMatCC == numOfCouples && jMatCC > 0 ) jMatCC--;
580 
581  if (fLambdaVector) delete fLambdaVector;
582  if (fdNdxCutVector) delete fdNdxCutVector;
583  if (fdNdxCutPhotonVector) delete fdNdxCutPhotonVector;
584  if (fdNdxCutPlasmonVector) delete fdNdxCutPlasmonVector;
585 
586  fLambdaVector = new G4PhysicsLogVector( fLowestKineticEnergy,
587  fHighestKineticEnergy,
588  fTotBin );
589  fdNdxCutVector = new G4PhysicsLogVector( fLowestKineticEnergy,
590  fHighestKineticEnergy,
591  fTotBin );
592  fdNdxCutPhotonVector = new G4PhysicsLogVector( fLowestKineticEnergy,
593  fHighestKineticEnergy,
594  fTotBin );
595  fdNdxCutPlasmonVector = new G4PhysicsLogVector( fLowestKineticEnergy,
596  fHighestKineticEnergy,
597  fTotBin );
598 
599 
600  photonCutInKineticEnergyNow = photEnergy;
601  deltaCutInKineticEnergyNow = eTkin;
602 
603  if(fVerbose > 0)
604  {
605  G4cout<<"PAIPhotonModel deltaCutInKineticEnergyNow = "
606  <<deltaCutInKineticEnergyNow/keV<<" keV"<<G4endl;
607  G4cout<<"PAIPhotonModel photonCutInKineticEnergyNow = "
608  <<photonCutInKineticEnergyNow/keV<<" keV"<<G4endl;
609  }
610  for ( i = 0; i <= fTotBin; i++ )
611  {
612  dNdxPhotonCut = GetdNdxPhotonCut(i,photonCutInKineticEnergyNow);
613  dNdxPlasmonCut = GetdNdxPlasmonCut(i,deltaCutInKineticEnergyNow);
614 
615  dNdxCut = dNdxPhotonCut + dNdxPlasmonCut;
616  lambda = dNdxCut <= DBL_MIN ? DBL_MAX: 1.0/dNdxCut;
617 
618  if (lambda <= 1000*kCarTolerance) lambda = 1000*kCarTolerance; // Mmm ???
619 
620  fLambdaVector->PutValue(i, lambda);
621 
622  fdNdxCutVector->PutValue(i, dNdxCut);
623  fdNdxCutPhotonVector->PutValue(i, dNdxPhotonCut);
624  fdNdxCutPlasmonVector->PutValue(i, dNdxPlasmonCut);
625  }
626 }
G4GeometryTolerance::GetSurfaceTolerance
G4double GetSurfaceTolerance() const
Definition: G4GeometryTolerance.cc:87
G4PAIPhotonModel::GetdNdxPlasmonCut
G4double GetdNdxPlasmonCut(G4int iPlace, G4double transferCut)
Definition: G4PAIPhotonModel.cc:723
G4int
int G4int
Definition: G4Types.hh:78
G4cout
G4GLOB_DLL std::ostream G4cout
G4PAIPhotonModel::GetdNdxPhotonCut
G4double GetdNdxPhotonCut(G4int iPlace, G4double transferCut)
Definition: G4PAIPhotonModel.cc:678
G4PhysicsVector::PutValue
void PutValue(size_t index, G4double theValue)
G4ProductionCutsTable::GetProductionCutsTable
static G4ProductionCutsTable * GetProductionCutsTable()
Definition: G4ProductionCutsTable.cc:63
G4ProductionCutsTable::GetMaterialCutsCouple
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
Definition: G4ProductionCutsTable.hh:262
DBL_MIN
#define DBL_MIN
Definition: templates.hh:75
G4endl
#define G4endl
Definition: G4ios.hh:61
G4double
double G4double
Definition: G4Types.hh:76
DBL_MAX
#define DBL_MAX
Definition: templates.hh:83
G4GeometryTolerance::GetInstance
static G4GeometryTolerance * GetInstance()
Definition: G4GeometryTolerance.cc:74
void G4PAIPhotonModel::BuildPAIonisationTable ( )

Definition at line 350 of file G4PAIPhotonModel.cc.

References DBL_MIN, python.hepunit::eV, G4PAIxSection::GetIntegralCerenkov(), G4PAIxSection::GetIntegralPAIdEdx(), G4PAIxSection::GetIntegralPAIxSection(), G4PAIxSection::GetIntegralPlasmon(), G4PhysicsVector::GetLowEdgeEnergy(), G4PAIxSection::GetMeanEnergyLoss(), G4SandiaTable::GetSandiaMatTablePAI(), G4PAIxSection::GetSplineEnergy(), G4PAIxSection::GetSplineSize(), G4PAIxSection::Initialize(), G4PhysicsTable::insertAt(), MaxSecondaryEnergy(), python.hepunit::proton_mass_c2, G4PhysicsFreeVector::PutValue(), and G4PhysicsVector::PutValue().

Referenced by Initialise(), and InitTest().

351 {
352  G4double LowEdgeEnergy , ionloss;
353  G4double /*massRatio,*/ tau, Tmax, Tmin, Tkin, deltaLow, /*gamma,*/ bg2;
354  /*
355  if( fPAItransferTable )
356  {
357  fPAItransferTable->clearAndDestroy();
358  delete fPAItransferTable;
359  }
360  */
361  fPAItransferTable = new G4PhysicsTable(fTotBin);
362  /*
363  if( fPAIratioTable )
364  {
365  fPAIratioTable->clearAndDestroy();
366  delete fPAIratioTable;
367  }
368  */
369  fPAIphotonTable = new G4PhysicsTable(fTotBin);
370  fPAIplasmonTable = new G4PhysicsTable(fTotBin);
371  /*
372  if( fPAIdEdxTable )
373  {
374  fPAIdEdxTable->clearAndDestroy();
375  delete fPAIdEdxTable;
376  }
377  */
378  fPAIdEdxTable = new G4PhysicsTable(fTotBin);
379 
380  // if(fdEdxVector) delete fdEdxVector;
381  fdEdxVector = new G4PhysicsLogVector( fLowestKineticEnergy,
382  fHighestKineticEnergy,
383  fTotBin );
384  // Tmin = fSandiaPhotoAbsCof[0][0]; // low energy Sandia interval
385  Tmin = fSandia.GetSandiaMatTablePAI(0,0); // low energy Sandia interval
386  deltaLow = 100.*eV; // 0.5*eV;
387 
388  for (G4int i = 0; i <= fTotBin; i++) //The loop for the kinetic energy
389  {
390  LowEdgeEnergy = fProtonEnergyVector->GetLowEdgeEnergy(i);
391  tau = LowEdgeEnergy/proton_mass_c2;
392  // if(tau < 0.01) tau = 0.01;
393  //gamma = tau +1.;
394  // G4cout<<"gamma = "<<gamma<<endl;
395  bg2 = tau*(tau + 2. );
396  //massRatio = electron_mass_c2/proton_mass_c2;
397  Tmax = MaxSecondaryEnergy(fParticle, LowEdgeEnergy);
398  // G4cout<<"proton Tkin = "<<LowEdgeEnergy/MeV<<" MeV"
399  // <<" Tmax = "<<Tmax/MeV<<" MeV"<<G4endl;
400  // Tkin = DeltaCutInKineticEnergyNow;
401 
402  // if ( DeltaCutInKineticEnergyNow > Tmax) // was <
403  Tkin = Tmax;
404  if ( Tkin < Tmin + deltaLow ) // low energy safety
405  {
406  Tkin = Tmin + deltaLow;
407  }
408  fPAIxSection.Initialize(fMaterial, Tkin, bg2,
409  &fSandia);
410  /*
411  G4PAIxSection protonPAI( fMatIndex,
412  Tkin,
413  bg2,
414  fSandiaPhotoAbsCof,
415  fSandiaIntervalNumber );
416 
417  */
418  // G4cout<<"ionloss = "<<ionloss*cm/keV<<" keV/cm"<<endl;
419  // G4cout<<"n1 = "<<fPAIxSection.GetIntegralPAIxSection(1)*cm<<" 1/cm"<<endl;
420  // G4cout<<"fPAIxSection.GetSplineSize() = "<<
421  // fPAIxSection.GetSplineSize()<<G4endl<<G4endl;
422 
423  G4PhysicsFreeVector* transferVector = new
424  G4PhysicsFreeVector(fPAIxSection.GetSplineSize());
425  G4PhysicsFreeVector* photonVector = new
426  G4PhysicsFreeVector(fPAIxSection.GetSplineSize());
427  G4PhysicsFreeVector* plasmonVector = new
428  G4PhysicsFreeVector(fPAIxSection.GetSplineSize());
429  G4PhysicsFreeVector* dEdxVector = new
430  G4PhysicsFreeVector(fPAIxSection.GetSplineSize());
431 
432  for( G4int k = 0; k < fPAIxSection.GetSplineSize(); k++ )
433  {
434  transferVector->PutValue( k ,
435  fPAIxSection.GetSplineEnergy(k+1),
436  fPAIxSection.GetIntegralPAIxSection(k+1) );
437  photonVector->PutValue( k ,
438  fPAIxSection.GetSplineEnergy(k+1),
439  fPAIxSection.GetIntegralCerenkov(k+1) );
440  plasmonVector->PutValue( k ,
441  fPAIxSection.GetSplineEnergy(k+1),
442  fPAIxSection.GetIntegralPlasmon(k+1) );
443  dEdxVector->PutValue( k ,
444  fPAIxSection.GetSplineEnergy(k+1),
445  fPAIxSection.GetIntegralPAIdEdx(k+1) );
446  }
447  ionloss = fPAIxSection.GetMeanEnergyLoss(); // total <dE/dx>
448  if ( ionloss <= 0.) ionloss = DBL_MIN;
449  fdEdxVector->PutValue(i,ionloss);
450 
451  fPAItransferTable->insertAt(i,transferVector);
452  fPAIphotonTable->insertAt(i,photonVector);
453  fPAIplasmonTable->insertAt(i,plasmonVector);
454  fPAIdEdxTable->insertAt(i,dEdxVector);
455 
456  } // end of Tkin loop
457  // theLossTable->insert(fdEdxVector);
458  // end of material loop
459  // G4cout<<"G4PAIonisation::BuildPAIonisationTable() have been called"<<G4endl;
460  // G4cout<<"G4PAIonisation::BuildLossTable() have been called"<<G4endl;
461 }
G4PhysicsFreeVector::PutValue
void PutValue(size_t binNumber, G4double binValue, G4double dataValue)
Definition: G4PhysicsFreeVector.cc:94
G4PAIxSection::GetIntegralPlasmon
G4double GetIntegralPlasmon(G4int i) const
Definition: G4PAIxSection.hh:336
G4PAIxSection::GetIntegralPAIxSection
G4double GetIntegralPAIxSection(G4int i) const
Definition: G4PAIxSection.hh:312
G4SandiaTable::GetSandiaMatTablePAI
G4double GetSandiaMatTablePAI(G4int, G4int)
Definition: G4SandiaTable.cc:986
G4PAIxSection::GetIntegralPAIdEdx
G4double GetIntegralPAIdEdx(G4int i) const
Definition: G4PAIxSection.hh:318
G4PhysicsVector::GetLowEdgeEnergy
G4double GetLowEdgeEnergy(size_t binNumber) const
Definition: G4PhysicsVector.cc:157
G4int
int G4int
Definition: G4Types.hh:78
G4PAIxSection::GetSplineSize
G4int GetSplineSize() const
Definition: G4PAIxSection.hh:187
G4PAIPhotonModel::MaxSecondaryEnergy
G4double MaxSecondaryEnergy(const G4ParticleDefinition *, G4double kinEnergy)
Definition: G4PAIPhotonModel.cc:1674
G4PAIxSection::GetIntegralCerenkov
G4double GetIntegralCerenkov(G4int i) const
Definition: G4PAIxSection.hh:324
G4PhysicsVector::PutValue
void PutValue(size_t index, G4double theValue)
python.hepunit.proton_mass_c2
float proton_mass_c2
Definition: hepunit.py:275
G4PAIxSection::Initialize
void Initialize(const G4Material *material, G4double maxEnergyTransfer, G4double betaGammaSq, G4SandiaTable *)
Definition: G4PAIxSection.cc:606
DBL_MIN
#define DBL_MIN
Definition: templates.hh:75
G4PhysicsTable::insertAt
void insertAt(size_t, G4PhysicsVector *)
G4double
double G4double
Definition: G4Types.hh:76
G4PAIxSection::GetMeanEnergyLoss
G4double GetMeanEnergyLoss() const
Definition: G4PAIxSection.hh:199
G4PAIxSection::GetSplineEnergy
G4double GetSplineEnergy(G4int i) const
Definition: G4PAIxSection.hh:306
G4double G4PAIPhotonModel::ComputeDEDXPerVolume ( const G4Material *  ,
const G4ParticleDefinition *  p,
G4double  kineticEnergy,
G4double  cutEnergy 
)
virtual

Reimplemented from G4VEmModel.

Definition at line 811 of file G4PAIPhotonModel.cc.

References G4VEmModel::CurrentCouple(), python.hepunit::eplus, GetdEdxCut(), G4ParticleDefinition::GetPDGCharge(), G4ParticleDefinition::GetPDGMass(), and python.hepunit::proton_mass_c2.

815 {
816  G4int iTkin,iPlace;
817  size_t jMat;
818 
819  //G4double cut = std::min(MaxSecondaryEnergy(p, kineticEnergy), cutEnergy);
820  G4double cut = cutEnergy;
821 
822  G4double particleMass = p->GetPDGMass();
823  G4double scaledTkin = kineticEnergy*proton_mass_c2/particleMass;
824  G4double charge = p->GetPDGCharge()/eplus;
825  G4double charge2 = charge*charge;
826  G4double dEdx = 0.;
827  const G4MaterialCutsCouple* matCC = CurrentCouple();
828 
829  for( jMat = 0;jMat < fMaterialCutsCoupleVector.size(); ++jMat )
830  {
831  if( matCC == fMaterialCutsCoupleVector[jMat] ) break;
832  }
833  if(jMat == fMaterialCutsCoupleVector.size() && jMat > 0) jMat--;
834  /*
835  G4cout << "G4PAIPhotonModel::ComputeDEDXPerVolume: jMat= " << jMat
836  << " jMax= " << fMaterialCutsCoupleVector.size()
837  << " matCC: " << matCC;
838  if(matCC) G4cout << " mat: " << matCC->GetMaterial()->GetName();
839  G4cout << G4endl;
840  G4cout << fPAIdEdxTable << " " << fdEdxVector << " "
841  << fProtonEnergyVector << G4endl;
842  */
843  fPAIdEdxTable = fPAIdEdxBank[jMat];
844  fdEdxVector = fdEdxTable[jMat];
845  for(iTkin = 0; iTkin <= fTotBin; iTkin++)
846  {
847  if(scaledTkin < fProtonEnergyVector->GetLowEdgeEnergy(iTkin)) break;
848  }
849  iPlace = iTkin - 1;
850  if(iPlace < 0) iPlace = 0;
851  else if(iPlace > fTotBin) iPlace = fTotBin;
852  dEdx = charge2*( (*fdEdxVector)(iPlace) - GetdEdxCut(iPlace,cut) );
853 
854  if( dEdx < 0.) dEdx = 0.;
855  return dEdx;
856 }
G4PAIPhotonModel::GetdEdxCut
G4double GetdEdxCut(G4int iPlace, G4double transferCut)
Definition: G4PAIPhotonModel.cc:769
G4int
int G4int
Definition: G4Types.hh:78
G4VEmModel::CurrentCouple
const G4MaterialCutsCouple * CurrentCouple() const
Definition: G4VEmModel.hh:426
python.hepunit.proton_mass_c2
float proton_mass_c2
Definition: hepunit.py:275
G4ParticleDefinition::GetPDGMass
G4double GetPDGMass() const
Definition: G4ParticleDefinition.hh:161
G4double
double G4double
Definition: G4Types.hh:76
G4ParticleDefinition::GetPDGCharge
G4double GetPDGCharge() const
Definition: G4ParticleDefinition.hh:163
void G4PAIPhotonModel::ComputeSandiaPhotoAbsCof ( )

Definition at line 302 of file G4PAIPhotonModel.cc.

References G4Material::GetMaterialTable(), G4SandiaTable::GetPhotoAbsorpCof(), G4SandiaTable::SandiaIntervals(), and G4SandiaTable::SandiaMixing().

303 {
304  G4int i, j, numberOfElements;
305  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
306 
307  G4SandiaTable thisMaterialSandiaTable(fMatIndex);
308  numberOfElements = (*theMaterialTable)[fMatIndex]->
309  GetNumberOfElements();
310  G4int* thisMaterialZ = new G4int[numberOfElements];
311 
312  for(i=0;i<numberOfElements;i++)
313  {
314  thisMaterialZ[i] =
315  (G4int)(*theMaterialTable)[fMatIndex]->GetElement(i)->GetZ();
316  }
317  fSandiaIntervalNumber = thisMaterialSandiaTable.SandiaIntervals
318  (thisMaterialZ,numberOfElements);
319 
320  fSandiaIntervalNumber = thisMaterialSandiaTable.SandiaMixing
321  ( thisMaterialZ ,
322  (*theMaterialTable)[fMatIndex]->GetFractionVector() ,
323  numberOfElements,fSandiaIntervalNumber);
324 
325  fSandiaPhotoAbsCof = new G4double*[fSandiaIntervalNumber];
326 
327  for(i=0;i<fSandiaIntervalNumber;i++) fSandiaPhotoAbsCof[i] = new G4double[5];
328 
329  for( i = 0; i < fSandiaIntervalNumber; i++ )
330  {
331  fSandiaPhotoAbsCof[i][0] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i+1,0);
332 
333  for( j = 1; j < 5; j++ )
334  {
335  fSandiaPhotoAbsCof[i][j] = thisMaterialSandiaTable.
336  GetPhotoAbsorpCof(i+1,j)*
337  (*theMaterialTable)[fMatIndex]->GetDensity();
338  }
339  }
340  delete[] thisMaterialZ;
341 }
G4Material::GetMaterialTable
static G4MaterialTable * GetMaterialTable()
Definition: G4Material.cc:564
G4MaterialTable
std::vector< G4Material * > G4MaterialTable
Definition: G4MaterialTable.hh:41
G4int
int G4int
Definition: G4Types.hh:78
G4double
double G4double
Definition: G4Types.hh:76
G4double G4PAIPhotonModel::CrossSectionPerVolume ( const G4Material *  ,
const G4ParticleDefinition *  p,
G4double  kineticEnergy,
G4double  cutEnergy,
G4double  maxEnergy 
)
virtual

Reimplemented from G4VEmModel.

Definition at line 862 of file G4PAIPhotonModel.cc.

References G4VEmModel::CurrentCouple(), GetdNdxPhotonCut(), GetdNdxPlasmonCut(), G4ProductionCutsTable::GetMaterialCutsCouple(), G4ParticleDefinition::GetPDGCharge(), G4ParticleDefinition::GetPDGMass(), G4ProductionCutsTable::GetProductionCutsTable(), G4ProductionCutsTable::GetTableSize(), idxG4GammaCut, MaxSecondaryEnergy(), G4INCL::Math::min(), and python.hepunit::proton_mass_c2.

867 {
868  G4int iTkin,iPlace;
869  size_t jMat, jMatCC;
870  G4double tmax = std::min(MaxSecondaryEnergy(p, kineticEnergy), maxEnergy);
871  if(cutEnergy >= tmax) return 0.0;
872  G4double particleMass = p->GetPDGMass();
873  G4double scaledTkin = kineticEnergy*proton_mass_c2/particleMass;
874  G4double charge = p->GetPDGCharge();
875  G4double charge2 = charge*charge, cross, cross1, cross2;
876  G4double photon1, photon2, plasmon1, plasmon2;
877 
878  const G4MaterialCutsCouple* matCC = CurrentCouple();
879 
880  const G4ProductionCutsTable* theCoupleTable=
881  G4ProductionCutsTable::GetProductionCutsTable();
882 
883  size_t numOfCouples = theCoupleTable->GetTableSize();
884 
885  for (jMatCC = 0; jMatCC < numOfCouples; jMatCC++ )
886  {
887  if( matCC == theCoupleTable->GetMaterialCutsCouple(jMatCC) ) break;
888  }
889  if( jMatCC == numOfCouples && jMatCC > 0 ) jMatCC--;
890 
891  const vector<G4double>* photonCutInKineticEnergy = theCoupleTable->
892  GetEnergyCutsVector(idxG4GammaCut);
893 
894  G4double photonCut = (*photonCutInKineticEnergy)[jMatCC];
895 
896  for( jMat = 0;jMat < fMaterialCutsCoupleVector.size(); ++jMat )
897  {
898  if( matCC == fMaterialCutsCoupleVector[jMat] ) break;
899  }
900  if(jMat == fMaterialCutsCoupleVector.size() && jMat > 0) jMat--;
901 
902  fPAItransferTable = fPAIxscBank[jMat];
903  fPAIphotonTable = fPAIphotonBank[jMat];
904  fPAIplasmonTable = fPAIplasmonBank[jMat];
905 
906  for(iTkin = 0; iTkin <= fTotBin; iTkin++)
907  {
908  if(scaledTkin < fProtonEnergyVector->GetLowEdgeEnergy(iTkin)) break;
909  }
910  iPlace = iTkin - 1;
911  if(iPlace < 0) iPlace = 0;
912 
913  // G4cout<<"iPlace = "<<iPlace<<"; tmax = "
914  // <<tmax<<"; cutEnergy = "<<cutEnergy<<G4endl;
915  photon1 = GetdNdxPhotonCut(iPlace,tmax);
916  photon2 = GetdNdxPhotonCut(iPlace,photonCut);
917 
918  plasmon1 = GetdNdxPlasmonCut(iPlace,tmax);
919  plasmon2 = GetdNdxPlasmonCut(iPlace,cutEnergy);
920 
921  cross1 = photon1 + plasmon1;
922  // G4cout<<"cross1 = "<<cross1<<G4endl;
923  cross2 = photon2 + plasmon2;
924  // G4cout<<"cross2 = "<<cross2<<G4endl;
925  cross = (cross2 - cross1)*charge2;
926  // G4cout<<"cross = "<<cross<<G4endl;
927 
928  if( cross < 0. ) cross = 0.;
929  return cross;
930 }
idxG4GammaCut
Definition: G4ProductionCuts.hh:53
G4PAIPhotonModel::GetdNdxPlasmonCut
G4double GetdNdxPlasmonCut(G4int iPlace, G4double transferCut)
Definition: G4PAIPhotonModel.cc:723
G4int
int G4int
Definition: G4Types.hh:78
G4VEmModel::CurrentCouple
const G4MaterialCutsCouple * CurrentCouple() const
Definition: G4VEmModel.hh:426
G4PAIPhotonModel::GetdNdxPhotonCut
G4double GetdNdxPhotonCut(G4int iPlace, G4double transferCut)
Definition: G4PAIPhotonModel.cc:678
G4PAIPhotonModel::MaxSecondaryEnergy
G4double MaxSecondaryEnergy(const G4ParticleDefinition *, G4double kinEnergy)
Definition: G4PAIPhotonModel.cc:1674
python.hepunit.proton_mass_c2
float proton_mass_c2
Definition: hepunit.py:275
G4ProductionCutsTable::GetProductionCutsTable
static G4ProductionCutsTable * GetProductionCutsTable()
Definition: G4ProductionCutsTable.cc:63
G4ParticleDefinition::GetPDGMass
G4double GetPDGMass() const
Definition: G4ParticleDefinition.hh:161
G4ProductionCutsTable::GetMaterialCutsCouple
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
Definition: G4ProductionCutsTable.hh:262
G4INCL::Math::min
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
Definition: G4INCLGlobals.hh:121
G4double
double G4double
Definition: G4Types.hh:76
G4ParticleDefinition::GetPDGCharge
G4double GetPDGCharge() const
Definition: G4ParticleDefinition.hh:163
void G4PAIPhotonModel::DefineForRegion ( const G4Region *  r)
virtual

Reimplemented from G4VEmModel.

Definition at line 1692 of file G4PAIPhotonModel.cc.

1693 {
1694  fPAIRegionVector.push_back(r);
1695 }
G4double G4PAIPhotonModel::Dispersion ( const G4Material *  material,
const G4DynamicParticle *  aParticle,
G4double  tmax,
G4double  step 
)
virtual

Implements G4VEmFluctuationModel.

Definition at line 1642 of file G4PAIPhotonModel.cc.

References python.hepunit::eplus, G4DynamicParticle::GetCharge(), G4Material::GetElectronDensity(), G4DynamicParticle::GetKineticEnergy(), G4DynamicParticle::GetMass(), and python.hepunit::twopi_mc2_rcl2.

1646 {
1647  G4double particleMass = aParticle->GetMass();
1648  G4double electronDensity = material->GetElectronDensity();
1649  G4double kineticEnergy = aParticle->GetKineticEnergy();
1650  G4double q = aParticle->GetCharge()/eplus;
1651  G4double etot = kineticEnergy + particleMass;
1652  G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*particleMass)/(etot*etot);
1653  G4double siga = (1.0/beta2 - 0.5) * twopi_mc2_rcl2 * tmax * step
1654  * electronDensity * q * q;
1655 
1656  return siga;
1657 
1658  /*
1659  G4double loss, sumLoss=0., sumLoss2=0., sigma2, meanLoss=0.;
1660  for(G4int i = 0; i < fMeanNumber; i++)
1661  {
1662  loss = SampleFluctuations(material,aParticle,tmax,step,meanLoss);
1663  sumLoss += loss;
1664  sumLoss2 += loss*loss;
1665  }
1666  meanLoss = sumLoss/fMeanNumber;
1667  sigma2 = meanLoss*meanLoss + (sumLoss2-2*sumLoss*meanLoss)/fMeanNumber;
1668  return sigma2;
1669  */
1670 }
G4DynamicParticle::GetKineticEnergy
G4double GetKineticEnergy() const
G4Material::GetElectronDensity
G4double GetElectronDensity() const
Definition: G4Material.hh:215
G4DynamicParticle::GetMass
G4double GetMass() const
G4DynamicParticle::GetCharge
G4double GetCharge() const
G4double
double G4double
Definition: G4Types.hh:76
G4double G4PAIPhotonModel::GetAlongStepTransfer ( G4PhysicsTable *  pTable,
G4PhysicsLogVector *  pVector,
G4int  iPlace,
G4double  scaledTkin,
G4double  step,
G4double  cof 
)

Definition at line 1496 of file G4PAIPhotonModel.cc.

References DBL_MAX, G4UniformRand, GetEnergyTransfer(), G4PhysicsVector::GetLowEdgeEnergy(), G4InuclParticleNames::lambda, position, and G4INCL::DeJongSpin::shoot().

Referenced by SampleFluctuations().

1500 {
1501  G4int iTkin = iPlace + 1, iTransfer;
1502  G4double loss = 0., position, E1, E2, W1, W2, W, dNdxCut1, dNdxCut2, meanNumber;
1503  G4double lambda, stepDelta, stepSum=0.;
1504  G4long numOfCollisions=0;
1505  G4bool numb = true;
1506 
1507  dNdxCut1 = (*pVector)(iPlace);
1508 
1509  // G4cout<<"iPlace = "<<iPlace<<endl;
1510 
1511  if(iTkin == fTotBin) // Fermi plato, try from left
1512  {
1513  meanNumber = ((*(*pTable)(iPlace))(0) - dNdxCut1)*cof;
1514  if(meanNumber < 0.) meanNumber = 0.;
1515  // numOfCollisions = RandPoisson::shoot(meanNumber);
1516  if( meanNumber > 0.) lambda = step/meanNumber;
1517  else lambda = DBL_MAX;
1518  while(numb)
1519  {
1520  stepDelta = G4RandExponential::shoot(lambda);
1521  stepSum += stepDelta;
1522  if(stepSum >= step) break;
1523  numOfCollisions++;
1524  }
1525 
1526  // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl;
1527 
1528  while(numOfCollisions)
1529  {
1530  position = dNdxCut1+
1531  ((*(*pTable)(iPlace))(0) - dNdxCut1)*G4UniformRand();
1532 
1533  for( iTransfer = 0;
1534  iTransfer < G4int((*pTable)(iPlace)->GetVectorLength()); iTransfer++ )
1535  {
1536  if(position >= (*(*pTable)(iPlace))(iTransfer)) break;
1537  }
1538  loss += GetEnergyTransfer(pTable,iPlace,position,iTransfer);
1539  numOfCollisions--;
1540  }
1541  }
1542  else
1543  {
1544  dNdxCut2 = (*pVector)(iPlace+1);
1545 
1546  if(iTkin == 0) // Tkin is too small, trying from right only
1547  {
1548  meanNumber = ((*(*pTable)(iPlace+1))(0) - dNdxCut2)*cof;
1549  if( meanNumber < 0. ) meanNumber = 0.;
1550  // numOfCollisions = G4RandPoisson::shoot(meanNumber);
1551  if( meanNumber > 0.) lambda = step/meanNumber;
1552  else lambda = DBL_MAX;
1553  while(numb)
1554  {
1555  stepDelta = G4RandExponential::shoot(lambda);
1556  stepSum += stepDelta;
1557  if(stepSum >= step) break;
1558  numOfCollisions++;
1559  }
1560 
1561  // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl;
1562 
1563  while(numOfCollisions)
1564  {
1565  position = dNdxCut2+
1566  ((*(*pTable)(iPlace+1))(0) - dNdxCut2)*G4UniformRand();
1567 
1568  for( iTransfer = 0;
1569  iTransfer < G4int((*pTable)(iPlace+1)->GetVectorLength()); iTransfer++ )
1570  {
1571  if(position >= (*(*pTable)(iPlace+1))(iTransfer)) break;
1572  }
1573  loss += GetEnergyTransfer(pTable,iPlace+1,position,iTransfer);
1574  numOfCollisions--;
1575  }
1576  }
1577  else // general case: Tkin between two vectors of the material
1578  {
1579  E1 = fProtonEnergyVector->GetLowEdgeEnergy(iTkin - 1);
1580  E2 = fProtonEnergyVector->GetLowEdgeEnergy(iTkin) ;
1581  W = 1.0/(E2 - E1);
1582  W1 = (E2 - scaledTkin)*W;
1583  W2 = (scaledTkin - E1)*W;
1584 
1585  // G4cout<<"(*(*pTable)(iPlace))(0) = "<<
1586  // (*(*pTable)(iPlace))(0)<<G4endl;
1587  // G4cout<<"(*(*pTable)(iPlace+1))(0) = "<<
1588  // (*(*pTable)(iPlace+1))(0)<<G4endl;
1589 
1590  meanNumber=( ((*(*pTable)(iPlace))(0)-dNdxCut1)*W1 +
1591  ((*(*pTable)(iPlace+1))(0)-dNdxCut2)*W2 )*cof;
1592  if(meanNumber<0.0) meanNumber = 0.0;
1593  // numOfCollisions = G4RandPoisson::shoot(meanNumber);
1594  if( meanNumber > 0.) lambda = step/meanNumber;
1595  else lambda = DBL_MAX;
1596  while(numb)
1597  {
1598  stepDelta = G4RandExponential::shoot(lambda);
1599  stepSum += stepDelta;
1600  if(stepSum >= step) break;
1601  numOfCollisions++;
1602  }
1603 
1604  // G4cout<<"numOfCollisions = "<<numOfCollisions<<endl;
1605 
1606  while(numOfCollisions)
1607  {
1608  position = dNdxCut1*W1 + dNdxCut2*W2 +
1609  ( ( (*(*pTable)(iPlace ))(0) - dNdxCut1)*W1 +
1610 
1611  ( (*(*pTable)(iPlace+1))(0) - dNdxCut2)*W2 )*G4UniformRand();
1612 
1613  // G4cout<<position<<"\t";
1614 
1615  for( iTransfer = 0;
1616  iTransfer < G4int((*pTable)(iPlace)->GetVectorLength()); iTransfer++ )
1617  {
1618  if( position >=
1619  ( (*(*pTable)(iPlace))(iTransfer)*W1 +
1620  (*(*pTable)(iPlace+1))(iTransfer)*W2) )
1621  {
1622  break;
1623  }
1624  }
1625  // loss += (*pTable)(iPlace)->GetLowEdgeEnergy(iTransfer);
1626  loss += GetEnergyTransfer(pTable,iPlace,position,iTransfer);
1627  numOfCollisions--;
1628  }
1629  }
1630  }
1631 
1632  return loss;
1633 
1634 }
G4INCL::DeJongSpin::shoot
ThreeVector shoot(const G4int Ap, const G4int Af)
Definition: G4INCLDeJongSpin.cc:67
G4PAIPhotonModel::GetEnergyTransfer
G4double GetEnergyTransfer(G4PhysicsTable *, G4int iPlace, G4double position, G4int iTransfer)
Definition: G4PAIPhotonModel.cc:1391
G4long
long G4long
Definition: G4Types.hh:80
G4PhysicsVector::GetLowEdgeEnergy
G4double GetLowEdgeEnergy(size_t binNumber) const
Definition: G4PhysicsVector.cc:157
G4int
int G4int
Definition: G4Types.hh:78
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:87
G4bool
bool G4bool
Definition: G4Types.hh:79
position
int position
Definition: filter.cc:7
G4double
double G4double
Definition: G4Types.hh:76
DBL_MAX
#define DBL_MAX
Definition: templates.hh:83
G4double G4PAIPhotonModel::GetdEdxCut ( G4int  iPlace,
G4double  transferCut 
)

Definition at line 769 of file G4PAIPhotonModel.cc.

References python.hepunit::eV.

Referenced by ComputeDEDXPerVolume().

770 {
771  G4int iTransfer;
772  G4double x1, x2, y1, y2, dEdxCut;
773  // G4cout<<"iPlace = "<<iPlace<<"; "<<"transferCut = "<<transferCut<<G4endl;
774  // G4cout<<"size = "<<G4int((*fPAIdEdxTable)(iPlace)->GetVectorLength())
775  // <<G4endl;
776  for( iTransfer = 0;
777  iTransfer < G4int((*fPAIdEdxTable)(iPlace)->GetVectorLength());
778  iTransfer++)
779  {
780  if(transferCut <= (*fPAIdEdxTable)(iPlace)->GetLowEdgeEnergy(iTransfer))
781  {
782  break;
783  }
784  }
785  if ( iTransfer >= G4int((*fPAIdEdxTable)(iPlace)->GetVectorLength()) )
786  {
787  iTransfer = (*fPAIdEdxTable)(iPlace)->GetVectorLength() - 1;
788  }
789  if (iTransfer == 0) return (*(*fPAIdEdxTable)(iPlace))(iTransfer);
790  y1 = (*(*fPAIdEdxTable)(iPlace))(iTransfer-1);
791  y2 = (*(*fPAIdEdxTable)(iPlace))(iTransfer);
792  // G4cout<<"y1 = "<<y1<<"; "<<"y2 = "<<y2<<G4endl;
793  x1 = (*fPAIdEdxTable)(iPlace)->GetLowEdgeEnergy(iTransfer-1);
794  x2 = (*fPAIdEdxTable)(iPlace)->GetLowEdgeEnergy(iTransfer);
795  // G4cout<<"x1 = "<<x1<<"; "<<"x2 = "<<x2<<G4endl;
796 
797  if ( y1 == y2 ) dEdxCut = y2;
798  else
799  {
800  // if ( x1 == x2 ) dEdxCut = y1 + (y2 - y1)*G4UniformRand();
801  // if ( std::abs(x1-x2) <= eV ) dEdxCut = y1 + (y2 - y1)*G4UniformRand();
802  if ( std::abs(x1-x2) <= eV ) dEdxCut = y1 + (y2 - y1)*0.5;
803  else dEdxCut = y1 + (transferCut - x1)*(y2 - y1)/(x2 - x1);
804  }
805  // G4cout<<""<<dEdxCut<<G4endl;
806  return dEdxCut;
807 }
G4int
int G4int
Definition: G4Types.hh:78
G4double
double G4double
Definition: G4Types.hh:76
G4double G4PAIPhotonModel::GetdNdxCut ( G4int  iPlace,
G4double  transferCut 
)

Definition at line 633 of file G4PAIPhotonModel.cc.

References python.hepunit::eV.

634 {
635  G4int iTransfer;
636  G4double x1, x2, y1, y2, dNdxCut;
637  // G4cout<<"iPlace = "<<iPlace<<"; "<<"transferCut = "<<transferCut<<G4endl;
638  // G4cout<<"size = "<<G4int((*fPAItransferTable)(iPlace)->GetVectorLength())
639  // <<G4endl;
640  for( iTransfer = 0;
641  iTransfer < G4int((*fPAItransferTable)(iPlace)->GetVectorLength());
642  iTransfer++)
643  {
644  if(transferCut <= (*fPAItransferTable)(iPlace)->GetLowEdgeEnergy(iTransfer))
645  {
646  break;
647  }
648  }
649  if ( iTransfer >= G4int((*fPAItransferTable)(iPlace)->GetVectorLength()) )
650  {
651  iTransfer = (*fPAItransferTable)(iPlace)->GetVectorLength() - 1;
652  }
653  if (iTransfer == 0) return (*(*fPAItransferTable)(iPlace))(iTransfer);
654  y1 = (*(*fPAItransferTable)(iPlace))(iTransfer-1);
655  y2 = (*(*fPAItransferTable)(iPlace))(iTransfer);
656  // G4cout<<"y1 = "<<y1<<"; "<<"y2 = "<<y2<<G4endl;
657  x1 = (*fPAItransferTable)(iPlace)->GetLowEdgeEnergy(iTransfer-1);
658  x2 = (*fPAItransferTable)(iPlace)->GetLowEdgeEnergy(iTransfer);
659  // G4cout<<"x1 = "<<x1<<"; "<<"x2 = "<<x2<<G4endl;
660 
661  if ( y1 == y2 ) dNdxCut = y2;
662  else
663  {
664  // if ( x1 == x2 ) dNdxCut = y1 + (y2 - y1)*G4UniformRand();
665  // if ( std::abs(x1-x2) <= eV ) dNdxCut = y1 + (y2 - y1)*G4UniformRand();
666  if ( std::abs(x1-x2) <= eV ) dNdxCut = y1 + (y2 - y1)*0.5;
667  else dNdxCut = y1 + (transferCut - x1)*(y2 - y1)/(x2 - x1);
668  }
669  // G4cout<<""<<dNdxCut<<G4endl;
670  return dNdxCut;
671 }
G4int
int G4int
Definition: G4Types.hh:78
G4double
double G4double
Definition: G4Types.hh:76
G4double G4PAIPhotonModel::GetdNdxPhotonCut ( G4int  iPlace,
G4double  transferCut 
)

Definition at line 678 of file G4PAIPhotonModel.cc.

References python.hepunit::eV.

Referenced by BuildLambdaVector(), CrossSectionPerVolume(), and GetXscPerVolume().

679 {
680  G4int iTransfer;
681  G4double x1, x2, y1, y2, dNdxCut;
682  // G4cout<<"iPlace = "<<iPlace<<"; "<<"transferCut = "<<transferCut<<G4endl;
683  // G4cout<<"size = "<<G4int((*fPAIphotonTable)(iPlace)->GetVectorLength())
684  // <<G4endl;
685  for( iTransfer = 0;
686  iTransfer < G4int((*fPAIphotonTable)(iPlace)->GetVectorLength());
687  iTransfer++)
688  {
689  if(transferCut <= (*fPAIphotonTable)(iPlace)->GetLowEdgeEnergy(iTransfer))
690  {
691  break;
692  }
693  }
694  if ( iTransfer >= G4int((*fPAIphotonTable)(iPlace)->GetVectorLength()) )
695  {
696  iTransfer = (*fPAIphotonTable)(iPlace)->GetVectorLength() - 1;
697  }
698  if (iTransfer == 0) return (*(*fPAIphotonTable)(iPlace))(iTransfer);
699  y1 = (*(*fPAIphotonTable)(iPlace))(iTransfer-1);
700  y2 = (*(*fPAIphotonTable)(iPlace))(iTransfer);
701  // G4cout<<"y1 = "<<y1<<"; "<<"y2 = "<<y2<<G4endl;
702  x1 = (*fPAIphotonTable)(iPlace)->GetLowEdgeEnergy(iTransfer-1);
703  x2 = (*fPAIphotonTable)(iPlace)->GetLowEdgeEnergy(iTransfer);
704  // G4cout<<"x1 = "<<x1<<"; "<<"x2 = "<<x2<<G4endl;
705 
706  if ( y1 == y2 ) dNdxCut = y2;
707  else
708  {
709  // if ( x1 == x2 ) dNdxCut = y1 + (y2 - y1)*G4UniformRand();
710  // if ( std::abs(x1-x2) <= eV ) dNdxCut = y1 + (y2 - y1)*G4UniformRand();
711  if ( std::abs(x1-x2) <= eV ) dNdxCut = y1 + (y2 - y1)*0.5;
712  else dNdxCut = y1 + (transferCut - x1)*(y2 - y1)/(x2 - x1);
713  }
714  // G4cout<<""<<dNdxPhotonCut<<G4endl;
715  return dNdxCut;
716 }
G4int
int G4int
Definition: G4Types.hh:78
G4double
double G4double
Definition: G4Types.hh:76
G4double G4PAIPhotonModel::GetdNdxPlasmonCut ( G4int  iPlace,
G4double  transferCut 
)

Definition at line 723 of file G4PAIPhotonModel.cc.

References python.hepunit::eV.

Referenced by BuildLambdaVector(), CrossSectionPerVolume(), and GetXscPerVolume().

724 {
725  G4int iTransfer;
726  G4double x1, x2, y1, y2, dNdxCut;
727 
728  // G4cout<<"iPlace = "<<iPlace<<"; "<<"transferCut = "<<transferCut<<G4endl;
729  // G4cout<<"size = "<<G4int((*fPAIPlasmonTable)(iPlace)->GetVectorLength())
730  // <<G4endl;
731  for( iTransfer = 0;
732  iTransfer < G4int((*fPAIplasmonTable)(iPlace)->GetVectorLength());
733  iTransfer++)
734  {
735  if(transferCut <= (*fPAIplasmonTable)(iPlace)->GetLowEdgeEnergy(iTransfer))
736  {
737  break;
738  }
739  }
740  if ( iTransfer >= G4int((*fPAIplasmonTable)(iPlace)->GetVectorLength()) )
741  {
742  iTransfer = (*fPAIplasmonTable)(iPlace)->GetVectorLength() - 1;
743  }
744  if (iTransfer == 0) return (*(*fPAIplasmonTable)(iPlace))(iTransfer);
745  y1 = (*(*fPAIplasmonTable)(iPlace))(iTransfer-1);
746  y2 = (*(*fPAIplasmonTable)(iPlace))(iTransfer);
747  // G4cout<<"y1 = "<<y1<<"; "<<"y2 = "<<y2<<G4endl;
748  x1 = (*fPAIplasmonTable)(iPlace)->GetLowEdgeEnergy(iTransfer-1);
749  x2 = (*fPAIplasmonTable)(iPlace)->GetLowEdgeEnergy(iTransfer);
750  // G4cout<<"x1 = "<<x1<<"; "<<"x2 = "<<x2<<G4endl;
751 
752  if ( y1 == y2 ) dNdxCut = y2;
753  else
754  {
755  // if ( x1 == x2 ) dNdxCut = y1 + (y2 - y1)*G4UniformRand();
756  // if ( std::abs(x1-x2) <= eV ) dNdxCut = y1 + (y2 - y1)*G4UniformRand();
757  if ( std::abs(x1-x2) <= eV ) dNdxCut = y1 + (y2 - y1)*0.5;
758  else dNdxCut = y1 + (transferCut - x1)*(y2 - y1)/(x2 - x1);
759  }
760  // G4cout<<""<<dNdxPlasmonCut<<G4endl;
761  return dNdxCut;
762 }
G4int
int G4int
Definition: G4Types.hh:78
G4double
double G4double
Definition: G4Types.hh:76
G4double G4PAIPhotonModel::GetEnergyTransfer ( G4PhysicsTable *  pTable,
G4int  iPlace,
G4double  position,
G4int  iTransfer 
)

Definition at line 1391 of file G4PAIPhotonModel.cc.

References G4UniformRand.

Referenced by GetAlongStepTransfer(), and GetPostStepTransfer().

1393 {
1394  G4int iTransferMax;
1395  G4double x1, x2, y1, y2, energyTransfer;
1396 
1397  if(iTransfer == 0)
1398  {
1399  energyTransfer = (*pTable)(iPlace)->GetLowEdgeEnergy(iTransfer);
1400  }
1401  else
1402  {
1403  iTransferMax = G4int((*pTable)(iPlace)->GetVectorLength());
1404 
1405  if ( iTransfer >= iTransferMax) iTransfer = iTransferMax - 1;
1406 
1407  y1 = (*(*pTable)(iPlace))(iTransfer-1);
1408  y2 = (*(*pTable)(iPlace))(iTransfer);
1409 
1410  x1 = (*pTable)(iPlace)->GetLowEdgeEnergy(iTransfer-1);
1411  x2 = (*pTable)(iPlace)->GetLowEdgeEnergy(iTransfer);
1412 
1413  if ( x1 == x2 ) energyTransfer = x2;
1414  else
1415  {
1416  if ( y1 == y2 ) energyTransfer = x1 + (x2 - x1)*G4UniformRand();
1417  else
1418  {
1419  energyTransfer = x1 + (position - y1)*(x2 - x1)/(y2 - y1);
1420  }
1421  }
1422  }
1423  return energyTransfer;
1424 }
G4int
int G4int
Definition: G4Types.hh:78
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:87
G4double
double G4double
Definition: G4Types.hh:76
G4double G4PAIPhotonModel::GetPostStepTransfer ( G4PhysicsTable *  pTable,
G4PhysicsLogVector *  pVector,
G4int  iPlace,
G4double  scaledTkin 
)

Definition at line 1313 of file G4PAIPhotonModel.cc.

References G4UniformRand, GetEnergyTransfer(), G4PhysicsVector::GetLowEdgeEnergy(), and position.

Referenced by SampleSecondaries(), and TestSecondaries().

1316 {
1317  // G4cout<<"G4PAIPhotonModel::GetPostStepTransfer"<<G4endl;
1318 
1319  G4int iTkin = iPlace+1, iTransfer;
1320  G4double transfer = 0.0, position, dNdxCut1, dNdxCut2, E1, E2, W1, W2, W;
1321 
1322  dNdxCut1 = (*pVector)(iPlace);
1323 
1324  // G4cout<<"iPlace = "<<iPlace<<endl;
1325 
1326  if(iTkin == fTotBin) // Fermi plato, try from left
1327  {
1328  position = dNdxCut1*G4UniformRand();
1329 
1330  for( iTransfer = 0;
1331  iTransfer < G4int((*pTable)(iPlace)->GetVectorLength()); iTransfer++ )
1332  {
1333  if(position >= (*(*pTable)(iPlace))(iTransfer)) break;
1334  }
1335  transfer = GetEnergyTransfer(pTable,iPlace,position,iTransfer);
1336  }
1337  else
1338  {
1339  dNdxCut2 = (*pVector)(iPlace+1);
1340  if(iTkin == 0) // Tkin is too small, trying from right only
1341  {
1342  position = dNdxCut2*G4UniformRand();
1343 
1344  for( iTransfer = 0;
1345  iTransfer < G4int((*pTable)(iPlace+1)->GetVectorLength()); iTransfer++ )
1346  {
1347  if(position >= (*(*pTable)(iPlace+1))(iTransfer)) break;
1348  }
1349  transfer = GetEnergyTransfer(pTable,iPlace+1,position,iTransfer);
1350  }
1351  else // general case: Tkin between two vectors of the material
1352  {
1353  E1 = fProtonEnergyVector->GetLowEdgeEnergy(iTkin - 1);
1354  E2 = fProtonEnergyVector->GetLowEdgeEnergy(iTkin) ;
1355  W = 1.0/(E2 - E1);
1356  W1 = (E2 - scaledTkin)*W;
1357  W2 = (scaledTkin - E1)*W;
1358 
1359  position = ( dNdxCut1*W1 + dNdxCut2*W2 )*G4UniformRand();
1360 
1361  // G4cout<<position<<"\t";
1362 
1363  G4int iTrMax1, iTrMax2, iTrMax;
1364 
1365  iTrMax1 = G4int((*pTable)(iPlace)->GetVectorLength());
1366  iTrMax2 = G4int((*pTable)(iPlace+1)->GetVectorLength());
1367 
1368  if (iTrMax1 >= iTrMax2) iTrMax = iTrMax2;
1369  else iTrMax = iTrMax1;
1370 
1371  for( iTransfer = 0; iTransfer < iTrMax; iTransfer++ )
1372  {
1373  if( position >=
1374  ( (*(*pTable)(iPlace))(iTransfer)*W1 +
1375  (*(*pTable)(iPlace+1))(iTransfer)*W2) ) break;
1376  }
1377  transfer = GetEnergyTransfer(pTable, iPlace, position, iTransfer);
1378  }
1379  }
1380  // G4cout<<"PAIPhotonModel PostStepTransfer = "<<transfer/keV<<" keV"<<G4endl;
1381  if( transfer < 0.0 ) transfer = 0.0;
1382  return transfer;
1383 }
G4PAIPhotonModel::GetEnergyTransfer
G4double GetEnergyTransfer(G4PhysicsTable *, G4int iPlace, G4double position, G4int iTransfer)
Definition: G4PAIPhotonModel.cc:1391
G4PhysicsVector::GetLowEdgeEnergy
G4double GetLowEdgeEnergy(size_t binNumber) const
Definition: G4PhysicsVector.cc:157
G4int
int G4int
Definition: G4Types.hh:78
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:87
position
int position
Definition: filter.cc:7
G4double
double G4double
Definition: G4Types.hh:76
G4double G4PAIPhotonModel::GetXscPerVolume ( const G4Material *  ,
const G4ParticleDefinition *  p,
G4double  kineticEnergy,
G4double  photonCut,
G4double  cutEnergy,
G4double  maxEnergy 
)

Definition at line 935 of file G4PAIPhotonModel.cc.

References G4VEmModel::CurrentCouple(), GetdNdxPhotonCut(), GetdNdxPlasmonCut(), G4ProductionCutsTable::GetMaterialCutsCouple(), G4ParticleDefinition::GetPDGCharge(), G4ParticleDefinition::GetPDGMass(), G4ProductionCutsTable::GetProductionCutsTable(), G4ProductionCutsTable::GetTableSize(), MaxSecondaryEnergy(), G4INCL::Math::min(), and python.hepunit::proton_mass_c2.

941 {
942  G4int iTkin,iPlace;
943  size_t jMat, jMatCC;
944  G4double tmax = std::min(MaxSecondaryEnergy(p, kineticEnergy), maxEnergy);
945  if(cutEnergy >= tmax) return 0.0;
946  G4double particleMass = p->GetPDGMass();
947  G4double scaledTkin = kineticEnergy*proton_mass_c2/particleMass;
948  G4double charge = p->GetPDGCharge();
949  G4double charge2 = charge*charge, cross, cross1, cross2;
950  G4double photon1, photon2, plasmon1, plasmon2;
951 
952  const G4MaterialCutsCouple* matCC = CurrentCouple();
953 
954  const G4ProductionCutsTable* theCoupleTable=
955  G4ProductionCutsTable::GetProductionCutsTable();
956 
957  size_t numOfCouples = theCoupleTable->GetTableSize();
958 
959  for (jMatCC = 0; jMatCC < numOfCouples; jMatCC++ )
960  {
961  if( matCC == theCoupleTable->GetMaterialCutsCouple(jMatCC) ) break;
962  }
963  if( jMatCC == numOfCouples && jMatCC > 0 ) jMatCC--;
964 
965  // const vector<G4double>* photonCutInKineticEnergy = theCoupleTable->
966  // GetEnergyCutsVector(idxG4GammaCut);
967  // G4double photonCut = (*photonCutInKineticEnergy)[jMatCC];
968 
969  for( jMat = 0;jMat < fMaterialCutsCoupleVector.size(); ++jMat )
970  {
971  if( matCC == fMaterialCutsCoupleVector[jMat] ) break;
972  }
973  if(jMat == fMaterialCutsCoupleVector.size() && jMat > 0) jMat--;
974 
975  fPAItransferTable = fPAIxscBank[jMat];
976  fPAIphotonTable = fPAIphotonBank[jMat];
977  fPAIplasmonTable = fPAIplasmonBank[jMat];
978 
979  for(iTkin = 0; iTkin <= fTotBin; iTkin++)
980  {
981  if(scaledTkin < fProtonEnergyVector->GetLowEdgeEnergy(iTkin)) break;
982  }
983  iPlace = iTkin - 1;
984  if(iPlace < 0) iPlace = 0;
985 
986  // G4cout<<"iPlace = "<<iPlace<<"; tmax = "
987  // <<tmax<<"; cutEnergy = "<<cutEnergy<<G4endl;
988  photon1 = GetdNdxPhotonCut(iPlace,tmax);
989  photon2 = GetdNdxPhotonCut(iPlace,photonCut);
990 
991  plasmon1 = GetdNdxPlasmonCut(iPlace,tmax);
992  plasmon2 = GetdNdxPlasmonCut(iPlace,cutEnergy);
993 
994  cross1 = photon1 + plasmon1;
995  // G4cout<<"cross1 = "<<cross1<<G4endl;
996  cross2 = photon2 + plasmon2;
997  // G4cout<<"cross2 = "<<cross2<<G4endl;
998  cross = (cross2 - cross1)*charge2;
999  // G4cout<<"cross = "<<cross<<G4endl;
1000 
1001  if( cross < 0. ) cross = 0.;
1002  return cross;
1003 }
G4PAIPhotonModel::GetdNdxPlasmonCut
G4double GetdNdxPlasmonCut(G4int iPlace, G4double transferCut)
Definition: G4PAIPhotonModel.cc:723
G4int
int G4int
Definition: G4Types.hh:78
G4VEmModel::CurrentCouple
const G4MaterialCutsCouple * CurrentCouple() const
Definition: G4VEmModel.hh:426
G4PAIPhotonModel::GetdNdxPhotonCut
G4double GetdNdxPhotonCut(G4int iPlace, G4double transferCut)
Definition: G4PAIPhotonModel.cc:678
G4PAIPhotonModel::MaxSecondaryEnergy
G4double MaxSecondaryEnergy(const G4ParticleDefinition *, G4double kinEnergy)
Definition: G4PAIPhotonModel.cc:1674
python.hepunit.proton_mass_c2
float proton_mass_c2
Definition: hepunit.py:275
G4ProductionCutsTable::GetProductionCutsTable
static G4ProductionCutsTable * GetProductionCutsTable()
Definition: G4ProductionCutsTable.cc:63
G4ParticleDefinition::GetPDGMass
G4double GetPDGMass() const
Definition: G4ParticleDefinition.hh:161
G4ProductionCutsTable::GetMaterialCutsCouple
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
Definition: G4ProductionCutsTable.hh:262
G4INCL::Math::min
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
Definition: G4INCLGlobals.hh:121
G4double
double G4double
Definition: G4Types.hh:76
G4ParticleDefinition::GetPDGCharge
G4double GetPDGCharge() const
Definition: G4ParticleDefinition.hh:163
void G4PAIPhotonModel::Initialise ( const G4ParticleDefinition *  p,
const G4DataVector &   
)
virtual

Implements G4VEmModel.

Definition at line 165 of file G4PAIPhotonModel.cc.

References BuildLambdaVector(), BuildPAIonisationTable(), G4Region::FindCouple(), G4Material::GetMaterialTable(), G4Material::GetNumberOfMaterials(), G4VEmModel::GetParticleChangeForLoss(), G4SandiaTable::Initialize(), and eplot::material.

167 {
168  // G4cout<<"G4PAIPhotonModel::Initialise for "<<p->GetParticleName()<<G4endl;
169  if(isInitialised) { return; }
170  isInitialised = true;
171 
172  if(!fParticle) SetParticle(p);
173 
174  fParticleChange = GetParticleChangeForLoss();
175 
176  //const G4ProductionCutsTable* theCoupleTable =
177  // G4ProductionCutsTable::GetProductionCutsTable();
178  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
179  size_t numOfMat = G4Material::GetNumberOfMaterials();
180  size_t numRegions = fPAIRegionVector.size();
181 
182  for(size_t iReg = 0; iReg < numRegions; ++iReg) // region loop
183  {
184  const G4Region* curReg = fPAIRegionVector[iReg];
185  G4Region* reg = const_cast<G4Region*>(curReg);
186 
187  for(size_t jMat = 0; jMat < numOfMat; ++jMat) // material loop
188  {
189  G4Material* material = (*theMaterialTable)[jMat];
190  const G4MaterialCutsCouple* couple = reg->FindCouple(material);
191 
192  // theCoupleTable->GetMaterialCutsCouple( material,
193  // curReg->GetProductionCuts() );
194 
195  if(couple) {
196  //G4cout << "Reg <" <<curReg->GetName() << "> mat <"
197  // << material->GetName() << "> fCouple= "
198  // << couple<<" " << p->GetParticleName() <<G4endl;
199 
200  fMaterialCutsCoupleVector.push_back(couple);
201 
202  fMatIndex = jMat;
203  fMaterial = material;
204 
205  // ComputeSandiaPhotoAbsCof();
206  fSandia.Initialize(material);
207  BuildPAIonisationTable();
208  /*
209  if( fSandiaPhotoAbsCof ) // delete SANDIA cofs've been used for pai-xsc
210  {
211  for( G4int i = 0;i < fSandiaIntervalNumber;i++)
212  {
213  delete[] fSandiaPhotoAbsCof[i];
214  }
215  delete[] fSandiaPhotoAbsCof;
216  fSandiaPhotoAbsCof = 0;
217  }
218  */
219  fPAIxscBank.push_back(fPAItransferTable);
220  fPAIphotonBank.push_back(fPAIphotonTable);
221  fPAIplasmonBank.push_back(fPAIplasmonTable);
222  fPAIdEdxBank.push_back(fPAIdEdxTable);
223  fdEdxTable.push_back(fdEdxVector);
224 
225  BuildLambdaVector(couple);
226 
227  fdNdxCutTable.push_back(fdNdxCutVector);
228  fdNdxCutPhotonTable.push_back(fdNdxCutPhotonVector);
229  fdNdxCutPlasmonTable.push_back(fdNdxCutPlasmonVector);
230  fLambdaTable.push_back(fLambdaVector);
231  }
232  }
233  }
234 }
G4VEmModel::GetParticleChangeForLoss
G4ParticleChangeForLoss * GetParticleChangeForLoss()
Definition: G4VEmModel.cc:107
G4Material::GetMaterialTable
static G4MaterialTable * GetMaterialTable()
Definition: G4Material.cc:564
G4MaterialTable
std::vector< G4Material * > G4MaterialTable
Definition: G4MaterialTable.hh:41
G4PAIPhotonModel::BuildLambdaVector
void BuildLambdaVector(const G4MaterialCutsCouple *matCutsCouple)
Definition: G4PAIPhotonModel.cc:470
eplot.material
string material
Definition: eplot.py:19
G4Region::FindCouple
G4MaterialCutsCouple * FindCouple(G4Material *mat)
G4Material::GetNumberOfMaterials
static size_t GetNumberOfMaterials()
Definition: G4Material.cc:571
G4SandiaTable::Initialize
void Initialize(G4Material *)
Definition: G4SandiaTable.cc:625
virtual void G4PAIPhotonModel::InitialiseMe ( const G4ParticleDefinition *  )
inlinevirtual

Reimplemented from G4VEmFluctuationModel.

Definition at line 80 of file G4PAIPhotonModel.hh.

80 {};
void G4PAIPhotonModel::InitTest ( const G4ParticleDefinition *  p,
G4MaterialCutsCouple *  couple,
G4double  phE,
G4double  eTkin 
)
virtual

Definition at line 238 of file G4PAIPhotonModel.cc.

References BuildLambdaVector(), BuildPAIonisationTable(), G4MaterialCutsCouple::GetMaterial(), G4Material::GetMaterialTable(), G4Material::GetName(), G4Material::GetNumberOfMaterials(), G4VEmModel::GetParticleChangeForLoss(), G4SandiaTable::Initialize(), and eplot::material.

240 {
241  // G4cout<<"G4PAIPhotonModel::InitTest for "<<p->GetParticleName()<<G4endl;
242  if(isInitialised) { return; }
243  isInitialised = true;
244 
245  if( !fParticle ) SetParticle(p);
246 
247  fParticleChange = GetParticleChangeForLoss();
248 
249  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
250 
251  size_t jMat, numOfMat = G4Material::GetNumberOfMaterials();
252 
253 
254  // const G4MaterialCutsCouple* couple = new G4MaterialCutsCouple(material,cuts);
255 
256  if( couple )
257  {
258  const G4Material* material = couple->GetMaterial();
259 
260  fMaterialCutsCoupleVector.push_back(couple);
261 
262  for( jMat = 0; jMat < numOfMat; ++jMat ) // material loop
263  {
264  if( material->GetName() == (*theMaterialTable)[jMat]->GetName() ) break;
265  }
266  fMatIndex = jMat;
267  G4Material* mat = (*theMaterialTable)[jMat];
268  fMaterial = mat;
269 
270 
271  // ComputeSandiaPhotoAbsCof();
272  fSandia.Initialize(mat);
273  BuildPAIonisationTable();
274  /*
275  if( fSandiaPhotoAbsCof ) // delete SANDIA cofs've been used for pai-xsc
276  {
277  for( G4int i = 0;i < fSandiaIntervalNumber;i++)
278  {
279  delete[] fSandiaPhotoAbsCof[i];
280  }
281  delete[] fSandiaPhotoAbsCof;
282  fSandiaPhotoAbsCof = 0;
283  }
284  */
285  fPAIxscBank.push_back(fPAItransferTable);
286  fPAIphotonBank.push_back(fPAIphotonTable);
287  fPAIplasmonBank.push_back(fPAIplasmonTable);
288  fPAIdEdxBank.push_back(fPAIdEdxTable);
289  fdEdxTable.push_back(fdEdxVector);
290 
291  BuildLambdaVector(couple,phE,eTkin);
292 
293  fdNdxCutTable.push_back(fdNdxCutVector);
294  fdNdxCutPhotonTable.push_back(fdNdxCutPhotonVector);
295  fdNdxCutPlasmonTable.push_back(fdNdxCutPlasmonVector);
296  fLambdaTable.push_back(fLambdaVector);
297  }
298 }
G4VEmModel::GetParticleChangeForLoss
G4ParticleChangeForLoss * GetParticleChangeForLoss()
Definition: G4VEmModel.cc:107
G4Material::GetName
const G4String & GetName() const
Definition: G4Material.hh:176
G4Material::GetMaterialTable
static G4MaterialTable * GetMaterialTable()
Definition: G4Material.cc:564
G4MaterialTable
std::vector< G4Material * > G4MaterialTable
Definition: G4MaterialTable.hh:41
G4PAIPhotonModel::BuildLambdaVector
void BuildLambdaVector(const G4MaterialCutsCouple *matCutsCouple)
Definition: G4PAIPhotonModel.cc:470
eplot.material
string material
Definition: eplot.py:19
G4Material::GetNumberOfMaterials
static size_t GetNumberOfMaterials()
Definition: G4Material.cc:571
G4SandiaTable::Initialize
void Initialize(G4Material *)
Definition: G4SandiaTable.cc:625
G4MaterialCutsCouple::GetMaterial
const G4Material * GetMaterial() const
Definition: G4MaterialCutsCouple.hh:150
G4double G4PAIPhotonModel::MaxSecondaryEnergy ( const G4ParticleDefinition *  p,
G4double  kinEnergy 
)
protectedvirtual

Reimplemented from G4VEmModel.

Definition at line 1674 of file G4PAIPhotonModel.cc.

References python.hepunit::electron_mass_c2, and G4ParticleDefinition::GetPDGMass().

Referenced by BuildPAIonisationTable(), CrossSectionPerVolume(), GetXscPerVolume(), SampleSecondaries(), and TestSecondaries().

1676 {
1677  G4double tmax = kinEnergy;
1678  if(p == fElectron) tmax *= 0.5;
1679  else if(p != fPositron)
1680  {
1681  G4double mass = p->GetPDGMass();
1682  G4double ratio= electron_mass_c2/mass;
1683  G4double gamma= kinEnergy/mass + 1.0;
1684  tmax = 2.0*electron_mass_c2*(gamma*gamma - 1.) /
1685  (1. + 2.0*gamma*ratio + ratio*ratio);
1686  }
1687  return tmax;
1688 }
python.hepunit.electron_mass_c2
float electron_mass_c2
Definition: hepunit.py:274
G4ParticleDefinition::GetPDGMass
G4double GetPDGMass() const
Definition: G4ParticleDefinition.hh:161
G4double
double G4double
Definition: G4Types.hh:76
G4double G4PAIPhotonModel::SampleFluctuations ( const G4MaterialCutsCouple *  matCC,
const G4DynamicParticle *  aParticle,
G4double  ,
G4double  step,
G4double  eloss 
)
virtual

Implements G4VEmFluctuationModel.

Definition at line 1431 of file G4PAIPhotonModel.cc.

References GetAlongStepTransfer(), G4DynamicParticle::GetDefinition(), G4DynamicParticle::GetKineticEnergy(), G4ParticleDefinition::GetPDGCharge(), G4ParticleDefinition::GetPDGMass(), and python.hepunit::proton_mass_c2.

1436 {
1437  size_t jMat = 0;
1438  for(;jMat < fMaterialCutsCoupleVector.size(); ++jMat )
1439  {
1440  if( matCC == fMaterialCutsCoupleVector[jMat] ) break;
1441  }
1442  if(jMat == fMaterialCutsCoupleVector.size()) { return eloss; }
1443 
1444  fPAItransferTable = fPAIxscBank[jMat];
1445  fPAIphotonTable = fPAIphotonBank[jMat];
1446  fPAIplasmonTable = fPAIplasmonBank[jMat];
1447 
1448  fdNdxCutVector = fdNdxCutTable[jMat];
1449  fdNdxCutPhotonVector = fdNdxCutPhotonTable[jMat];
1450  fdNdxCutPlasmonVector = fdNdxCutPlasmonTable[jMat];
1451 
1452  G4int iTkin, iPlace ;
1453 
1454  // G4cout<<"G4PAIPhotonModel::SampleFluctuations"<<G4endl;
1455 
1456  G4double loss, photonLoss, plasmonLoss, charge2;
1457 
1458 
1459  G4double Tkin = aParticle->GetKineticEnergy();
1460  G4double MassRatio = proton_mass_c2/aParticle->GetDefinition()->GetPDGMass();
1461  G4double charge = aParticle->GetDefinition()->GetPDGCharge();
1462  charge2 = charge*charge;
1463  G4double scaledTkin = Tkin*MassRatio;
1464  G4double cof = step*charge2;
1465 
1466  for( iTkin = 0; iTkin <= fTotBin; iTkin++)
1467  {
1468  if(scaledTkin < fProtonEnergyVector->GetLowEdgeEnergy(iTkin)) break;
1469  }
1470  iPlace = iTkin - 1;
1471  if( iPlace < 0 ) iPlace = 0;
1472 
1473  photonLoss = GetAlongStepTransfer(fPAIphotonTable,fdNdxCutPhotonVector,
1474 iPlace,scaledTkin,step,cof);
1475 
1476  // G4cout<<"PAIPhotonModel AlongStepPhotonLoss = "<<photonLoss/keV<<" keV"<<G4endl;
1477 
1478  plasmonLoss = GetAlongStepTransfer(fPAIplasmonTable,fdNdxCutPlasmonVector,
1479 iPlace,scaledTkin,step,cof);
1480 
1481  // G4cout<<"PAIPhotonModel AlongStepPlasmonLoss = "<<plasmonLoss/keV<<" keV"<<G4endl;
1482 
1483  loss = photonLoss + plasmonLoss;
1484 
1485  // G4cout<<"PAIPhotonModel AlongStepLoss = "<<loss/keV<<" keV"<<G4endl;
1486 
1487  return loss;
1488 }
G4DynamicParticle::GetKineticEnergy
G4double GetKineticEnergy() const
G4DynamicParticle::GetDefinition
G4ParticleDefinition * GetDefinition() const
G4int
int G4int
Definition: G4Types.hh:78
python.hepunit.proton_mass_c2
float proton_mass_c2
Definition: hepunit.py:275
G4PAIPhotonModel::GetAlongStepTransfer
G4double GetAlongStepTransfer(G4PhysicsTable *, G4PhysicsLogVector *, G4int iPlace, G4double scaledTkin, G4double step, G4double cof)
Definition: G4PAIPhotonModel.cc:1496
G4ParticleDefinition::GetPDGMass
G4double GetPDGMass() const
Definition: G4ParticleDefinition.hh:161
G4double
double G4double
Definition: G4Types.hh:76
G4ParticleDefinition::GetPDGCharge
G4double GetPDGCharge() const
Definition: G4ParticleDefinition.hh:163
void G4PAIPhotonModel::SampleSecondaries ( std::vector< G4DynamicParticle * > *  vdp,
const G4MaterialCutsCouple *  matCC,
const G4DynamicParticle *  dp,
G4double  tmin,
G4double  maxEnergy 
)
virtual

Implements G4VEmModel.

Definition at line 1011 of file G4PAIPhotonModel.cc.

References G4Electron::Electron(), python.hepunit::electron_mass_c2, fStopAndKill, G4cout, G4endl, G4UniformRand, G4Gamma::Gamma(), G4DynamicParticle::GetDefinition(), G4DynamicParticle::GetKineticEnergy(), G4DynamicParticle::GetMass(), G4DynamicParticle::GetMomentumDirection(), GetPostStepTransfer(), MaxSecondaryEnergy(), G4INCL::Math::min(), G4VParticleChange::ProposeTrackStatus(), python.hepunit::proton_mass_c2, CLHEP::Hep3Vector::rotateUz(), G4DynamicParticle::SetDefinition(), G4DynamicParticle::SetKineticEnergy(), G4DynamicParticle::SetMomentumDirection(), G4ParticleChangeForLoss::SetProposedKineticEnergy(), G4ParticleChangeForLoss::SetProposedMomentumDirection(), python.hepunit::twopi, and CLHEP::Hep3Vector::unit().

1016 {
1017  size_t jMat;
1018  for( jMat = 0;jMat < fMaterialCutsCoupleVector.size(); ++jMat )
1019  {
1020  if( matCC == fMaterialCutsCoupleVector[jMat] ) break;
1021  }
1022  if( jMat == fMaterialCutsCoupleVector.size() && jMat > 0 ) jMat--;
1023 
1024  fPAItransferTable = fPAIxscBank[jMat];
1025  fPAIphotonTable = fPAIphotonBank[jMat];
1026  fPAIplasmonTable = fPAIplasmonBank[jMat];
1027 
1028  fdNdxCutVector = fdNdxCutTable[jMat];
1029  fdNdxCutPhotonVector = fdNdxCutPhotonTable[jMat];
1030  fdNdxCutPlasmonVector = fdNdxCutPlasmonTable[jMat];
1031 
1032  G4double tmax = std::min(MaxSecondaryEnergy(dp->GetDefinition(),dp->GetKineticEnergy()), maxEnergy);
1033  if( tmin >= tmax && fVerbose > 0)
1034  {
1035  G4cout<<"G4PAIPhotonModel::SampleSecondary: tmin >= tmax "<<G4endl;
1036  }
1037 
1038  G4ThreeVector direction = dp->GetMomentumDirection();
1039  G4double particleMass = dp->GetMass();
1040  G4double kineticEnergy = dp->GetKineticEnergy();
1041  G4double scaledTkin = kineticEnergy*proton_mass_c2/particleMass; // fMass
1042  G4double totalEnergy = kineticEnergy + particleMass;
1043  G4double pSquare = kineticEnergy*(totalEnergy+particleMass);
1044 
1045  G4int iTkin;
1046  for(iTkin=0;iTkin<=fTotBin;iTkin++)
1047  {
1048  if(scaledTkin < fProtonEnergyVector->GetLowEdgeEnergy(iTkin)) break;
1049  }
1050  G4int iPlace = iTkin - 1;
1051  if(iPlace < 0) iPlace = 0;
1052 
1053  G4double dNdxPhotonCut = (*fdNdxCutPhotonVector)(iPlace);
1054  G4double dNdxPlasmonCut = (*fdNdxCutPlasmonVector)(iPlace);
1055  G4double dNdxCut = dNdxPhotonCut + dNdxPlasmonCut;
1056 
1057  G4double ratio;
1058  if (dNdxCut > 0.) ratio = dNdxPhotonCut/dNdxCut;
1059  else return; // ratio = 0.;
1060 
1061  if(ratio < G4UniformRand() ) // secondary e-
1062  {
1063  G4double deltaTkin = GetPostStepTransfer(fPAIplasmonTable, fdNdxCutPlasmonVector,
1064  iPlace, scaledTkin);
1065 
1066 // G4cout<<"PAIPhotonModel PlasmonPostStepTransfer = "<<deltaTkin/keV<<" keV"<<G4endl;
1067 
1068  if( deltaTkin <= 0. )
1069  {
1070  G4cout<<"G4PAIPhotonModel::SampleSecondary e- deltaTkin = "<<deltaTkin<<G4endl;
1071  }
1072  if( deltaTkin <= 0.) return;
1073 
1074  if( deltaTkin >= kineticEnergy ) // stop primary
1075  {
1076  deltaTkin = kineticEnergy;
1077  kineticEnergy = 0.0;
1078  }
1079  G4double deltaTotalMomentum = sqrt(deltaTkin*(deltaTkin + 2. * electron_mass_c2 ));
1080  G4double totalMomentum = sqrt(pSquare);
1081  G4double costheta = deltaTkin*(totalEnergy + electron_mass_c2)
1082  /(deltaTotalMomentum * totalMomentum);
1083 
1084  if( costheta > 0.99999 ) costheta = 0.99999;
1085  G4double sintheta = 0.0;
1086  G4double sin2 = 1. - costheta*costheta;
1087  if( sin2 > 0.) sintheta = sqrt(sin2);
1088 
1089  // direction of the delta electron
1090 
1091  G4double phi = twopi*G4UniformRand();
1092  G4double dirx = sintheta*cos(phi), diry = sintheta*sin(phi), dirz = costheta;
1093 
1094  G4ThreeVector deltaDirection(dirx,diry,dirz);
1095  deltaDirection.rotateUz(direction);
1096 
1097  if( kineticEnergy > 0.) // primary change
1098  {
1099  kineticEnergy -= deltaTkin;
1100  G4ThreeVector dir = totalMomentum*direction - deltaTotalMomentum*deltaDirection;
1101  direction = dir.unit();
1102  fParticleChange->SetProposedKineticEnergy(kineticEnergy);
1103  fParticleChange->SetProposedMomentumDirection(direction);
1104  }
1105  else // stop primary
1106  {
1107  fParticleChange->ProposeTrackStatus(fStopAndKill);
1108  fParticleChange->SetProposedKineticEnergy(0.0);
1109  }
1110 
1111  // create G4DynamicParticle object for e- delta ray
1112 
1113  G4DynamicParticle* deltaRay = new G4DynamicParticle;
1114  deltaRay->SetDefinition(G4Electron::Electron());
1115  deltaRay->SetKineticEnergy( deltaTkin );
1116  deltaRay->SetMomentumDirection(deltaDirection);
1117  vdp->push_back(deltaRay);
1118 
1119  }
1120  else // secondary 'Cherenkov' photon
1121  {
1122  G4double deltaTkin = GetPostStepTransfer(fPAIphotonTable, fdNdxCutPhotonVector,
1123  iPlace,scaledTkin);
1124 
1125  // G4cout<<"PAIPhotonModel PhotonPostStepTransfer = "<<deltaTkin/keV<<" keV"<<G4endl;
1126 
1127  if( deltaTkin <= 0. )
1128  {
1129  G4cout<<"G4PAIPhotonModel::SampleSecondary gamma deltaTkin = "<<deltaTkin<<G4endl;
1130  }
1131  if( deltaTkin <= 0.) return;
1132 
1133  if( deltaTkin >= kineticEnergy ) // stop primary
1134  {
1135  deltaTkin = kineticEnergy;
1136  kineticEnergy = 0.0;
1137  }
1138  G4double costheta = 0.; // G4UniformRand(); // VG: ??? for start only
1139  G4double sintheta = sqrt((1.+costheta)*(1.-costheta));
1140 
1141  // direction of the 'Cherenkov' photon
1142  G4double phi = twopi*G4UniformRand();
1143  G4double dirx = sintheta*cos(phi), diry = sintheta*sin(phi), dirz = costheta;
1144 
1145  G4ThreeVector deltaDirection(dirx,diry,dirz);
1146  deltaDirection.rotateUz(direction);
1147 
1148  if( kineticEnergy > 0.) // primary change
1149  {
1150  kineticEnergy -= deltaTkin;
1151  fParticleChange->SetProposedKineticEnergy(kineticEnergy);
1152  }
1153  else // stop primary
1154  {
1155  fParticleChange->ProposeTrackStatus(fStopAndKill);
1156  fParticleChange->SetProposedKineticEnergy(0.0);
1157  }
1158  // create G4DynamicParticle object for photon ray
1159 
1160  G4DynamicParticle* photonRay = new G4DynamicParticle;
1161  photonRay->SetDefinition( G4Gamma::Gamma() );
1162  photonRay->SetKineticEnergy( deltaTkin );
1163  photonRay->SetMomentumDirection(deltaDirection);
1164 
1165  vdp->push_back(photonRay);
1166  }
1167 }
G4DynamicParticle::GetKineticEnergy
G4double GetKineticEnergy() const
G4DynamicParticle::SetMomentumDirection
void SetMomentumDirection(const G4ThreeVector &aDirection)
G4DynamicParticle::GetDefinition
G4ParticleDefinition * GetDefinition() const
G4int
int G4int
Definition: G4Types.hh:78
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:87
G4cout
G4GLOB_DLL std::ostream G4cout
G4DynamicParticle::GetMass
G4double GetMass() const
G4DynamicParticle::GetMomentumDirection
const G4ThreeVector & GetMomentumDirection() const
G4PAIPhotonModel::MaxSecondaryEnergy
G4double MaxSecondaryEnergy(const G4ParticleDefinition *, G4double kinEnergy)
Definition: G4PAIPhotonModel.cc:1674
G4ParticleChangeForLoss::SetProposedKineticEnergy
void SetProposedKineticEnergy(G4double proposedKinEnergy)
Definition: G4ParticleChangeForLoss.hh:150
python.hepunit.proton_mass_c2
float proton_mass_c2
Definition: hepunit.py:275
python.hepunit.electron_mass_c2
float electron_mass_c2
Definition: hepunit.py:274
G4Gamma::Gamma
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
G4ParticleChangeForLoss::SetProposedMomentumDirection
void SetProposedMomentumDirection(const G4ThreeVector &dir)
Definition: G4ParticleChangeForLoss.hh:194
G4DynamicParticle::SetKineticEnergy
void SetKineticEnergy(G4double aEnergy)
CLHEP::Hep3Vector::unit
Hep3Vector unit() const
G4INCL::Math::min
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
Definition: G4INCLGlobals.hh:121
fStopAndKill
Definition: G4TrackStatus.hh:56
G4Electron::Electron
static G4Electron * Electron()
Definition: G4Electron.cc:94
G4endl
#define G4endl
Definition: G4ios.hh:61
G4double
double G4double
Definition: G4Types.hh:76
G4VParticleChange::ProposeTrackStatus
void ProposeTrackStatus(G4TrackStatus status)
G4DynamicParticle::SetDefinition
void SetDefinition(const G4ParticleDefinition *aParticleDefinition)
Definition: G4DynamicParticle.cc:273
G4PAIPhotonModel::GetPostStepTransfer
G4double GetPostStepTransfer(G4PhysicsTable *, G4PhysicsLogVector *, G4int iPlace, G4double scaledTkin)
Definition: G4PAIPhotonModel.cc:1313
G4double G4PAIPhotonModel::TestSecondaries ( G4MaterialCutsCouple *  matCC,
G4DynamicParticle *  dp,
G4double  tmin,
G4double  maxEnergy 
)

Definition at line 1173 of file G4PAIPhotonModel.cc.

References G4Electron::Electron(), python.hepunit::electron_mass_c2, G4cout, G4endl, G4UniformRand, G4Gamma::Gamma(), G4DynamicParticle::GetDefinition(), G4DynamicParticle::GetKineticEnergy(), G4DynamicParticle::GetMass(), G4MaterialCutsCouple::GetMaterial(), G4DynamicParticle::GetMomentumDirection(), G4Material::GetName(), GetPostStepTransfer(), MaxSecondaryEnergy(), G4INCL::Math::min(), python.hepunit::proton_mass_c2, CLHEP::Hep3Vector::rotateUz(), G4DynamicParticle::SetDefinition(), G4DynamicParticle::SetKineticEnergy(), G4DynamicParticle::SetMomentumDirection(), G4ParticleChangeForLoss::SetProposedKineticEnergy(), G4ParticleChangeForLoss::SetProposedMomentumDirection(), python.hepunit::twopi, and CLHEP::Hep3Vector::unit().

1176 {
1177  size_t jMat;
1178  for( jMat = 0;jMat < fMaterialCutsCoupleVector.size(); ++jMat )
1179  {
1180  if( matCC->GetMaterial()->GetName() == fMaterialCutsCoupleVector[jMat]->GetMaterial()->GetName() ) break;
1181  }
1182  if( jMat == fMaterialCutsCoupleVector.size() && jMat > 0 ) jMat--;
1183 
1184  fPAItransferTable = fPAIxscBank[jMat];
1185  fPAIphotonTable = fPAIphotonBank[jMat];
1186  fPAIplasmonTable = fPAIplasmonBank[jMat];
1187 
1188  fdNdxCutVector = fdNdxCutTable[jMat];
1189  fdNdxCutPhotonVector = fdNdxCutPhotonTable[jMat];
1190  fdNdxCutPlasmonVector = fdNdxCutPlasmonTable[jMat];
1191 
1192  G4double tmax = std::min(MaxSecondaryEnergy(dp->GetDefinition(),dp->GetKineticEnergy()), maxEnergy);
1193  if( tmin >= tmax && fVerbose > 0)
1194  {
1195  G4cout<<"G4PAIPhotonModel::TestSecondaries: tmin >= tmax "<<G4endl;
1196  }
1197 
1198  G4ThreeVector direction = dp->GetMomentumDirection();
1199  G4double particleMass = dp->GetMass();
1200  G4double kineticEnergy = dp->GetKineticEnergy();
1201  G4double scaledTkin = kineticEnergy*proton_mass_c2/particleMass; // fMass
1202  G4double totalEnergy = kineticEnergy + particleMass;
1203  G4double pSquare = kineticEnergy*(totalEnergy+particleMass);
1204 
1205  G4int iTkin;
1206  for(iTkin=0;iTkin<=fTotBin;iTkin++)
1207  {
1208  if(scaledTkin < fProtonEnergyVector->GetLowEdgeEnergy(iTkin)) break;
1209  }
1210  G4int iPlace = iTkin - 1;
1211  if(iPlace < 0) iPlace = 0;
1212 
1213  G4double dNdxPhotonCut = (*fdNdxCutPhotonVector)(iPlace);
1214  G4double dNdxPlasmonCut = (*fdNdxCutPlasmonVector)(iPlace);
1215  G4double dNdxCut = dNdxPhotonCut + dNdxPlasmonCut;
1216 
1217  G4double ratio, deltaTkin;
1218  if (dNdxCut > 0.) ratio = dNdxPhotonCut/dNdxCut;
1219  else ratio = 0.;
1220 
1221  if(ratio < G4UniformRand() ) // secondary e-
1222  {
1223  deltaTkin = GetPostStepTransfer(fPAIplasmonTable, fdNdxCutPlasmonVector,
1224  iPlace, scaledTkin);
1225 
1226 // G4cout<<"PAIPhotonModel PlasmonPostStepTransfer = "<<deltaTkin/keV<<" keV"<<G4endl;
1227 
1228  if( deltaTkin <= 0. )
1229  {
1230  G4cout<<"G4PAIPhotonModel::SampleSecondary e- deltaTkin = "<<deltaTkin<<G4endl;
1231  }
1232  if( deltaTkin <= 0.) return 0.;
1233 
1234  G4double deltaTotalMomentum = sqrt(deltaTkin*(deltaTkin + 2. * electron_mass_c2 ));
1235  G4double totalMomentum = sqrt(pSquare);
1236  G4double costheta = deltaTkin*(totalEnergy + electron_mass_c2)
1237  /(deltaTotalMomentum * totalMomentum);
1238 
1239  if( costheta > 0.99999 ) costheta = 0.99999;
1240  G4double sintheta = 0.0;
1241  G4double sin2 = 1. - costheta*costheta;
1242  if( sin2 > 0.) sintheta = sqrt(sin2);
1243 
1244  // direction of the delta electron
1245 
1246  G4double phi = twopi*G4UniformRand();
1247  G4double dirx = sintheta*cos(phi), diry = sintheta*sin(phi), dirz = costheta;
1248 
1249  G4ThreeVector deltaDirection(dirx,diry,dirz);
1250  deltaDirection.rotateUz(direction);
1251 
1252  // primary change
1253 
1254  kineticEnergy -= deltaTkin;
1255  G4ThreeVector dir = totalMomentum*direction - deltaTotalMomentum*deltaDirection;
1256  direction = dir.unit();
1257  fParticleChange->SetProposedMomentumDirection(direction);
1258 
1259  // create G4DynamicParticle object for e- delta ray
1260 
1261  G4DynamicParticle* deltaRay = new G4DynamicParticle;
1262  deltaRay->SetDefinition(G4Electron::Electron());
1263  deltaRay->SetKineticEnergy( deltaTkin );
1264  deltaRay->SetMomentumDirection(deltaDirection);
1265 
1266  }
1267  else // secondary 'Cherenkov' photon
1268  {
1269  deltaTkin = GetPostStepTransfer(fPAIphotonTable, fdNdxCutPhotonVector,
1270  iPlace,scaledTkin);
1271 
1272  // G4cout<<"PAIPhotonModel PhotonPostStepTransfer = "<<deltaTkin/keV<<" keV"<<G4endl;
1273 
1274  if( deltaTkin <= 0. )
1275  {
1276  G4cout<<"G4PAIPhotonModel::SampleSecondary gamma deltaTkin = "<<deltaTkin<<G4endl;
1277  }
1278  if( deltaTkin <= 0.) return 0.;
1279 
1280  G4double costheta = 0.; // G4UniformRand(); // VG: ??? for start only
1281  G4double sintheta = sqrt((1.+costheta)*(1.-costheta));
1282 
1283  // direction of the 'Cherenkov' photon
1284  G4double phi = twopi*G4UniformRand();
1285  G4double dirx = sintheta*cos(phi), diry = sintheta*sin(phi), dirz = costheta;
1286 
1287  G4ThreeVector deltaDirection(dirx,diry,dirz);
1288  deltaDirection.rotateUz(direction);
1289 
1290  // primary change
1291  kineticEnergy -= deltaTkin;
1292 
1293  // create G4DynamicParticle object for photon ray
1294 
1295  G4DynamicParticle* photonRay = new G4DynamicParticle;
1296  photonRay->SetDefinition( G4Gamma::Gamma() );
1297  photonRay->SetKineticEnergy( deltaTkin );
1298  photonRay->SetMomentumDirection(deltaDirection);
1299 
1300  }
1301  fParticleChange->SetProposedKineticEnergy(kineticEnergy);
1302 
1303  return deltaTkin;
1304 }
G4DynamicParticle::GetKineticEnergy
G4double GetKineticEnergy() const
G4Material::GetName
const G4String & GetName() const
Definition: G4Material.hh:176
G4DynamicParticle::SetMomentumDirection
void SetMomentumDirection(const G4ThreeVector &aDirection)
G4DynamicParticle::GetDefinition
G4ParticleDefinition * GetDefinition() const
G4int
int G4int
Definition: G4Types.hh:78
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:87
G4cout
G4GLOB_DLL std::ostream G4cout
G4DynamicParticle::GetMass
G4double GetMass() const
G4DynamicParticle::GetMomentumDirection
const G4ThreeVector & GetMomentumDirection() const
G4PAIPhotonModel::MaxSecondaryEnergy
G4double MaxSecondaryEnergy(const G4ParticleDefinition *, G4double kinEnergy)
Definition: G4PAIPhotonModel.cc:1674
G4ParticleChangeForLoss::SetProposedKineticEnergy
void SetProposedKineticEnergy(G4double proposedKinEnergy)
Definition: G4ParticleChangeForLoss.hh:150
python.hepunit.proton_mass_c2
float proton_mass_c2
Definition: hepunit.py:275
python.hepunit.electron_mass_c2
float electron_mass_c2
Definition: hepunit.py:274
G4Gamma::Gamma
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
G4ParticleChangeForLoss::SetProposedMomentumDirection
void SetProposedMomentumDirection(const G4ThreeVector &dir)
Definition: G4ParticleChangeForLoss.hh:194
G4DynamicParticle::SetKineticEnergy
void SetKineticEnergy(G4double aEnergy)
CLHEP::Hep3Vector::unit
Hep3Vector unit() const
G4INCL::Math::min
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
Definition: G4INCLGlobals.hh:121
G4Electron::Electron
static G4Electron * Electron()
Definition: G4Electron.cc:94
G4endl
#define G4endl
Definition: G4ios.hh:61
G4double
double G4double
Definition: G4Types.hh:76
G4DynamicParticle::SetDefinition
void SetDefinition(const G4ParticleDefinition *aParticleDefinition)
Definition: G4DynamicParticle.cc:273
G4PAIPhotonModel::GetPostStepTransfer
G4double GetPostStepTransfer(G4PhysicsTable *, G4PhysicsLogVector *, G4int iPlace, G4double scaledTkin)
Definition: G4PAIPhotonModel.cc:1313
G4MaterialCutsCouple::GetMaterial
const G4Material * GetMaterial() const
Definition: G4MaterialCutsCouple.hh:150
The documentation for this class was generated from the following files:
Generated on Wed Apr 30 2014 15:56:51 for Geant4.10 by   doxygen 1.8.7