File indexing completed on 2025-01-18 09:59:20
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046
0047
0048
0049
0050
0051
0052
0053
0054
0055
0056
0057
0058
0059
0060
0061
0062
0063
0064
0065
0066
0067
0068
0069
0070
0071
0072
0073
0074
0075 #ifndef G4VEmModel_h
0076 #define G4VEmModel_h 1
0077
0078 #include "globals.hh"
0079 #include "G4DynamicParticle.hh"
0080 #include "G4ParticleDefinition.hh"
0081 #include "G4MaterialCutsCouple.hh"
0082 #include "G4Material.hh"
0083 #include "G4Element.hh"
0084 #include "G4ElementVector.hh"
0085 #include "G4Isotope.hh"
0086 #include "G4DataVector.hh"
0087 #include "G4VEmFluctuationModel.hh"
0088 #include "G4VEmAngularDistribution.hh"
0089 #include "G4EmElementSelector.hh"
0090 #include <CLHEP/Random/RandomEngine.h>
0091 #include <vector>
0092
0093 class G4ElementData;
0094 class G4PhysicsTable;
0095 class G4Region;
0096 class G4VParticleChange;
0097 class G4ParticleChangeForLoss;
0098 class G4ParticleChangeForGamma;
0099 class G4Track;
0100 class G4LossTableManager;
0101
0102 class G4VEmModel
0103 {
0104
0105 public:
0106
0107 explicit G4VEmModel(const G4String& nam);
0108
0109 virtual ~G4VEmModel();
0110
0111
0112
0113
0114
0115 virtual void Initialise(const G4ParticleDefinition*, const G4DataVector&) = 0;
0116
0117 virtual void SampleSecondaries(std::vector<G4DynamicParticle*>*,
0118 const G4MaterialCutsCouple*,
0119 const G4DynamicParticle*,
0120 G4double tmin = 0.0,
0121 G4double tmax = DBL_MAX) = 0;
0122
0123
0124
0125
0126
0127
0128 virtual void InitialiseLocal(const G4ParticleDefinition*,
0129 G4VEmModel* masterModel);
0130
0131
0132 virtual void InitialiseForMaterial(const G4ParticleDefinition*,
0133 const G4Material*);
0134
0135
0136 virtual void InitialiseForElement(const G4ParticleDefinition*,
0137 G4int Z);
0138
0139
0140
0141
0142
0143
0144 virtual G4double ComputeDEDXPerVolume(const G4Material*,
0145 const G4ParticleDefinition*,
0146 G4double kineticEnergy,
0147 G4double cutEnergy = DBL_MAX);
0148
0149
0150 virtual G4double CrossSectionPerVolume(const G4Material*,
0151 const G4ParticleDefinition*,
0152 G4double kineticEnergy,
0153 G4double cutEnergy = 0.0,
0154 G4double maxEnergy = DBL_MAX);
0155
0156
0157 virtual G4double GetPartialCrossSection(const G4Material*,
0158 G4int level,
0159 const G4ParticleDefinition*,
0160 G4double kineticEnergy);
0161
0162
0163 virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
0164 G4double kinEnergy,
0165 G4double Z,
0166 G4double A = 0.,
0167 G4double cutEnergy = 0.0,
0168 G4double maxEnergy = DBL_MAX);
0169
0170
0171 virtual G4double ComputeCrossSectionPerShell(const G4ParticleDefinition*,
0172 G4int Z, G4int shellIdx,
0173 G4double kinEnergy,
0174 G4double cutEnergy = 0.0,
0175 G4double maxEnergy = DBL_MAX);
0176
0177
0178 virtual G4double ChargeSquareRatio(const G4Track&);
0179
0180
0181 virtual G4double GetChargeSquareRatio(const G4ParticleDefinition*,
0182 const G4Material*,
0183 G4double kineticEnergy);
0184
0185
0186 virtual G4double GetParticleCharge(const G4ParticleDefinition*,
0187 const G4Material*,
0188 G4double kineticEnergy);
0189
0190
0191 virtual void StartTracking(G4Track*);
0192
0193
0194 virtual void CorrectionsAlongStep(const G4MaterialCutsCouple*,
0195 const G4DynamicParticle*,
0196 const G4double& length,
0197 G4double& eloss);
0198
0199
0200 virtual G4double Value(const G4MaterialCutsCouple*,
0201 const G4ParticleDefinition*,
0202 G4double kineticEnergy);
0203
0204
0205 virtual G4double MinPrimaryEnergy(const G4Material*,
0206 const G4ParticleDefinition*,
0207 G4double cut = 0.0);
0208
0209
0210 virtual G4double MinEnergyCut(const G4ParticleDefinition*,
0211 const G4MaterialCutsCouple*);
0212
0213
0214 virtual void SetupForMaterial(const G4ParticleDefinition*,
0215 const G4Material*,
0216 G4double kineticEnergy);
0217
0218
0219 virtual void DefineForRegion(const G4Region*);
0220
0221
0222 virtual void FillNumberOfSecondaries(G4int& numberOfTriplets,
0223 G4int& numberOfRecoil);
0224
0225
0226 virtual void ModelDescription(std::ostream& outFile) const;
0227
0228 protected:
0229
0230
0231 G4ParticleChangeForLoss* GetParticleChangeForLoss();
0232
0233
0234 G4ParticleChangeForGamma* GetParticleChangeForGamma();
0235
0236
0237 virtual G4double MaxSecondaryEnergy(const G4ParticleDefinition*,
0238 G4double kineticEnergy);
0239
0240 public:
0241
0242
0243
0244
0245
0246
0247 void InitialiseElementSelectors(const G4ParticleDefinition*,
0248 const G4DataVector&);
0249
0250
0251 inline std::vector<G4EmElementSelector*>* GetElementSelectors();
0252
0253
0254 inline void SetElementSelectors(std::vector<G4EmElementSelector*>*);
0255
0256
0257 inline G4double ComputeDEDX( const G4MaterialCutsCouple*,
0258 const G4ParticleDefinition*,
0259 G4double kineticEnergy,
0260 G4double cutEnergy = DBL_MAX);
0261
0262
0263 inline G4double CrossSection(const G4MaterialCutsCouple*,
0264 const G4ParticleDefinition*,
0265 G4double kineticEnergy,
0266 G4double cutEnergy = 0.0,
0267 G4double maxEnergy = DBL_MAX);
0268
0269
0270 inline G4double ComputeMeanFreePath(const G4ParticleDefinition*,
0271 G4double kineticEnergy,
0272 const G4Material*,
0273 G4double cutEnergy = 0.0,
0274 G4double maxEnergy = DBL_MAX);
0275
0276
0277 inline G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
0278 const G4Element*,
0279 G4double kinEnergy,
0280 G4double cutEnergy = 0.0,
0281 G4double maxEnergy = DBL_MAX);
0282
0283
0284 inline const G4Element* SelectRandomAtom(const G4MaterialCutsCouple*,
0285 const G4ParticleDefinition*,
0286 G4double kineticEnergy,
0287 G4double cutEnergy = 0.0,
0288 G4double maxEnergy = DBL_MAX);
0289
0290 inline const G4Element* SelectTargetAtom(const G4MaterialCutsCouple*,
0291 const G4ParticleDefinition*,
0292 G4double kineticEnergy,
0293 G4double logKineticEnergy,
0294 G4double cutEnergy = 0.0,
0295 G4double maxEnergy = DBL_MAX);
0296
0297
0298 const G4Element* SelectRandomAtom(const G4Material*,
0299 const G4ParticleDefinition*,
0300 G4double kineticEnergy,
0301 G4double cutEnergy = 0.0,
0302 G4double maxEnergy = DBL_MAX);
0303
0304
0305 const G4Element* GetCurrentElement(const G4Material* mat = nullptr) const;
0306 G4int SelectRandomAtomNumber(const G4Material*) const;
0307
0308
0309 const G4Isotope* GetCurrentIsotope(const G4Element* elm = nullptr) const;
0310 G4int SelectIsotopeNumber(const G4Element*) const;
0311
0312
0313
0314
0315
0316 void SetParticleChange(G4VParticleChange*, G4VEmFluctuationModel* f=nullptr);
0317
0318 void SetCrossSectionTable(G4PhysicsTable*, G4bool isLocal);
0319
0320 inline G4ElementData* GetElementData();
0321
0322 inline G4PhysicsTable* GetCrossSectionTable();
0323
0324 inline G4VEmFluctuationModel* GetModelOfFluctuations();
0325
0326 inline G4VEmAngularDistribution* GetAngularDistribution();
0327
0328 inline G4VEmModel* GetTripletModel();
0329
0330 inline void SetTripletModel(G4VEmModel*);
0331
0332 inline void SetAngularDistribution(G4VEmAngularDistribution*);
0333
0334 inline G4double HighEnergyLimit() const;
0335
0336 inline G4double LowEnergyLimit() const;
0337
0338 inline G4double HighEnergyActivationLimit() const;
0339
0340 inline G4double LowEnergyActivationLimit() const;
0341
0342 inline G4double PolarAngleLimit() const;
0343
0344 inline G4double SecondaryThreshold() const;
0345
0346 inline G4bool DeexcitationFlag() const;
0347
0348 inline G4bool ForceBuildTableFlag() const;
0349
0350 inline G4bool UseAngularGeneratorFlag() const;
0351
0352 inline void SetAngularGeneratorFlag(G4bool);
0353
0354 inline void SetHighEnergyLimit(G4double);
0355
0356 inline void SetLowEnergyLimit(G4double);
0357
0358 inline void SetActivationHighEnergyLimit(G4double);
0359
0360 inline void SetActivationLowEnergyLimit(G4double);
0361
0362 inline G4bool IsActive(G4double kinEnergy) const;
0363
0364 inline void SetPolarAngleLimit(G4double);
0365
0366 inline void SetSecondaryThreshold(G4double);
0367
0368 inline void SetDeexcitationFlag(G4bool val);
0369
0370 inline void SetForceBuildTable(G4bool val);
0371
0372 inline void SetFluctuationFlag(G4bool val);
0373
0374 inline void SetMasterThread(G4bool val);
0375
0376 inline G4bool IsMaster() const;
0377
0378 inline void SetUseBaseMaterials(G4bool val);
0379
0380 inline G4bool UseBaseMaterials() const;
0381
0382 inline G4double MaxSecondaryKinEnergy(const G4DynamicParticle* dynParticle);
0383
0384 inline const G4String& GetName() const;
0385
0386 inline void SetCurrentCouple(const G4MaterialCutsCouple*);
0387
0388 inline G4bool IsLocked() const;
0389
0390 inline void SetLocked(G4bool);
0391
0392
0393 [[deprecated("Use G4EmParameters::Instance()->SetLPM instead")]]
0394 void SetLPMFlag(G4bool);
0395
0396
0397 G4VEmModel & operator=(const G4VEmModel &right) = delete;
0398 G4VEmModel(const G4VEmModel&) = delete;
0399
0400 protected:
0401
0402 inline const G4MaterialCutsCouple* CurrentCouple() const;
0403
0404 inline void SetCurrentElement(const G4Element*);
0405
0406 private:
0407
0408
0409
0410 G4VEmFluctuationModel* flucModel = nullptr;
0411 G4VEmAngularDistribution* anglModel = nullptr;
0412 G4VEmModel* fTripletModel = nullptr;
0413 const G4MaterialCutsCouple* fCurrentCouple = nullptr;
0414 const G4Element* fCurrentElement = nullptr;
0415 std::vector<G4EmElementSelector*>* elmSelectors = nullptr;
0416 G4LossTableManager* fEmManager;
0417
0418 protected:
0419
0420 G4ElementData* fElementData = nullptr;
0421 G4VParticleChange* pParticleChange = nullptr;
0422 G4PhysicsTable* xSectionTable = nullptr;
0423 const G4Material* pBaseMaterial = nullptr;
0424 const std::vector<G4double>* theDensityFactor = nullptr;
0425 const std::vector<G4int>* theDensityIdx = nullptr;
0426
0427 G4double inveplus;
0428 G4double pFactor = 1.0;
0429
0430 private:
0431
0432 G4double lowLimit;
0433 G4double highLimit;
0434 G4double eMinActive = 0.0;
0435 G4double eMaxActive = DBL_MAX;
0436 G4double secondaryThreshold = DBL_MAX;
0437 G4double polarAngleLimit;
0438
0439 G4int nSelectors = 0;
0440 G4int nsec = 5;
0441
0442 protected:
0443
0444 std::size_t currentCoupleIndex = 0;
0445 std::size_t basedCoupleIndex = 0;
0446 G4bool lossFlucFlag = true;
0447
0448 private:
0449
0450 G4bool flagDeexcitation = false;
0451 G4bool flagForceBuildTable = false;
0452 G4bool isMaster = true;
0453
0454 G4bool localTable = true;
0455 G4bool localElmSelectors = true;
0456 G4bool useAngularGenerator = false;
0457 G4bool useBaseMaterials = false;
0458 G4bool isLocked = false;
0459
0460 const G4String name;
0461 std::vector<G4double> xsec;
0462
0463 };
0464
0465
0466
0467 inline void G4VEmModel::SetCurrentCouple(const G4MaterialCutsCouple* ptr)
0468 {
0469 if(fCurrentCouple != ptr) {
0470 fCurrentCouple = ptr;
0471 basedCoupleIndex = currentCoupleIndex = ptr->GetIndex();
0472 pBaseMaterial = ptr->GetMaterial();
0473 pFactor = 1.0;
0474 if(useBaseMaterials) {
0475 basedCoupleIndex = (*theDensityIdx)[currentCoupleIndex];
0476 if(nullptr != pBaseMaterial->GetBaseMaterial())
0477 pBaseMaterial = pBaseMaterial->GetBaseMaterial();
0478 pFactor = (*theDensityFactor)[currentCoupleIndex];
0479 }
0480 }
0481 }
0482
0483
0484
0485 inline const G4MaterialCutsCouple* G4VEmModel::CurrentCouple() const
0486 {
0487 return fCurrentCouple;
0488 }
0489
0490
0491
0492 inline void G4VEmModel::SetCurrentElement(const G4Element* elm)
0493 {
0494 fCurrentElement = elm;
0495 }
0496
0497
0498
0499 inline
0500 G4double G4VEmModel::MaxSecondaryKinEnergy(const G4DynamicParticle* dynPart)
0501 {
0502 return MaxSecondaryEnergy(dynPart->GetParticleDefinition(),
0503 dynPart->GetKineticEnergy());
0504 }
0505
0506
0507
0508 inline G4double G4VEmModel::ComputeDEDX(const G4MaterialCutsCouple* couple,
0509 const G4ParticleDefinition* part,
0510 G4double kinEnergy,
0511 G4double cutEnergy)
0512 {
0513 SetCurrentCouple(couple);
0514 return pFactor*ComputeDEDXPerVolume(pBaseMaterial,part,kinEnergy,cutEnergy);
0515 }
0516
0517
0518
0519 inline G4double G4VEmModel::CrossSection(const G4MaterialCutsCouple* couple,
0520 const G4ParticleDefinition* part,
0521 G4double kinEnergy,
0522 G4double cutEnergy,
0523 G4double maxEnergy)
0524 {
0525 SetCurrentCouple(couple);
0526 return pFactor*CrossSectionPerVolume(pBaseMaterial,part,kinEnergy,
0527 cutEnergy,maxEnergy);
0528 }
0529
0530
0531
0532 inline
0533 G4double G4VEmModel::ComputeMeanFreePath(const G4ParticleDefinition* part,
0534 G4double ekin,
0535 const G4Material* material,
0536 G4double emin,
0537 G4double emax)
0538 {
0539 G4double cross = CrossSectionPerVolume(material,part,ekin,emin,emax);
0540 return (cross > 0.0) ? 1./cross : DBL_MAX;
0541 }
0542
0543
0544
0545 inline G4double
0546 G4VEmModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition* part,
0547 const G4Element* elm,
0548 G4double kinEnergy,
0549 G4double cutEnergy,
0550 G4double maxEnergy)
0551 {
0552 fCurrentElement = elm;
0553 return ComputeCrossSectionPerAtom(part,kinEnergy,elm->GetZ(),elm->GetN(),
0554 cutEnergy,maxEnergy);
0555 }
0556
0557
0558
0559 inline const G4Element*
0560 G4VEmModel::SelectRandomAtom(const G4MaterialCutsCouple* couple,
0561 const G4ParticleDefinition* part,
0562 G4double kinEnergy,
0563 G4double cutEnergy,
0564 G4double maxEnergy)
0565 {
0566 SetCurrentCouple(couple);
0567 fCurrentElement = (nSelectors > 0) ?
0568 ((*elmSelectors)[couple->GetIndex()])->SelectRandomAtom(kinEnergy) :
0569 SelectRandomAtom(pBaseMaterial,part,kinEnergy,cutEnergy,maxEnergy);
0570 return fCurrentElement;
0571 }
0572
0573
0574
0575 inline const G4Element*
0576 G4VEmModel::SelectTargetAtom(const G4MaterialCutsCouple* couple,
0577 const G4ParticleDefinition* part,
0578 G4double kinEnergy,
0579 G4double logKinE,
0580 G4double cutEnergy,
0581 G4double maxEnergy)
0582 {
0583 SetCurrentCouple(couple);
0584 fCurrentElement = (nSelectors > 0)
0585 ? ((*elmSelectors)[couple->GetIndex()])->SelectRandomAtom(kinEnergy,logKinE)
0586 : SelectRandomAtom(pBaseMaterial,part,kinEnergy,cutEnergy,maxEnergy);
0587 return fCurrentElement;
0588 }
0589
0590
0591
0592 inline G4VEmFluctuationModel* G4VEmModel::GetModelOfFluctuations()
0593 {
0594 return flucModel;
0595 }
0596
0597
0598
0599 inline G4VEmAngularDistribution* G4VEmModel::GetAngularDistribution()
0600 {
0601 return anglModel;
0602 }
0603
0604
0605
0606 inline void G4VEmModel::SetAngularDistribution(G4VEmAngularDistribution* p)
0607 {
0608 if(p != anglModel) {
0609 delete anglModel;
0610 anglModel = p;
0611 }
0612 }
0613
0614
0615
0616 inline G4VEmModel* G4VEmModel::GetTripletModel()
0617 {
0618 return fTripletModel;
0619 }
0620
0621
0622
0623 inline void G4VEmModel::SetTripletModel(G4VEmModel* p)
0624 {
0625 if(p != fTripletModel) {
0626 delete fTripletModel;
0627 fTripletModel = p;
0628 }
0629 }
0630
0631
0632
0633 inline G4double G4VEmModel::HighEnergyLimit() const
0634 {
0635 return highLimit;
0636 }
0637
0638
0639
0640 inline G4double G4VEmModel::LowEnergyLimit() const
0641 {
0642 return lowLimit;
0643 }
0644
0645
0646
0647 inline G4double G4VEmModel::HighEnergyActivationLimit() const
0648 {
0649 return eMaxActive;
0650 }
0651
0652
0653
0654 inline G4double G4VEmModel::LowEnergyActivationLimit() const
0655 {
0656 return eMinActive;
0657 }
0658
0659
0660
0661 inline G4double G4VEmModel::PolarAngleLimit() const
0662 {
0663 return polarAngleLimit;
0664 }
0665
0666
0667
0668 inline G4double G4VEmModel::SecondaryThreshold() const
0669 {
0670 return secondaryThreshold;
0671 }
0672
0673
0674
0675 inline G4bool G4VEmModel::DeexcitationFlag() const
0676 {
0677 return flagDeexcitation;
0678 }
0679
0680
0681
0682 inline G4bool G4VEmModel::ForceBuildTableFlag() const
0683 {
0684 return flagForceBuildTable;
0685 }
0686
0687
0688
0689 inline G4bool G4VEmModel::UseAngularGeneratorFlag() const
0690 {
0691 return useAngularGenerator;
0692 }
0693
0694
0695
0696 inline void G4VEmModel::SetAngularGeneratorFlag(G4bool val)
0697 {
0698 useAngularGenerator = val;
0699 }
0700
0701
0702
0703 inline void G4VEmModel::SetFluctuationFlag(G4bool val)
0704 {
0705 lossFlucFlag = val;
0706 }
0707
0708
0709
0710 inline void G4VEmModel::SetMasterThread(G4bool val)
0711 {
0712 isMaster = val;
0713 }
0714
0715
0716
0717 inline G4bool G4VEmModel::IsMaster() const
0718 {
0719 return isMaster;
0720 }
0721
0722
0723
0724 inline void G4VEmModel::SetUseBaseMaterials(G4bool val)
0725 {
0726 useBaseMaterials = val;
0727 }
0728
0729
0730
0731 inline G4bool G4VEmModel::UseBaseMaterials() const
0732 {
0733 return useBaseMaterials;
0734 }
0735
0736
0737
0738 inline void G4VEmModel::SetHighEnergyLimit(G4double val)
0739 {
0740 highLimit = val;
0741 }
0742
0743
0744
0745 inline void G4VEmModel::SetLowEnergyLimit(G4double val)
0746 {
0747 lowLimit = val;
0748 }
0749
0750
0751
0752 inline void G4VEmModel::SetActivationHighEnergyLimit(G4double val)
0753 {
0754 eMaxActive = val;
0755 }
0756
0757
0758
0759 inline void G4VEmModel::SetActivationLowEnergyLimit(G4double val)
0760 {
0761 eMinActive = val;
0762 }
0763
0764
0765
0766 inline G4bool G4VEmModel::IsActive(G4double kinEnergy) const
0767 {
0768 return (kinEnergy >= eMinActive && kinEnergy <= eMaxActive);
0769 }
0770
0771
0772
0773 inline void G4VEmModel::SetPolarAngleLimit(G4double val)
0774 {
0775 if(!isLocked) { polarAngleLimit = val; }
0776 }
0777
0778
0779
0780 inline void G4VEmModel::SetSecondaryThreshold(G4double val)
0781 {
0782 secondaryThreshold = val;
0783 }
0784
0785
0786
0787 inline void G4VEmModel::SetDeexcitationFlag(G4bool val)
0788 {
0789 flagDeexcitation = val;
0790 }
0791
0792
0793
0794 inline void G4VEmModel::SetForceBuildTable(G4bool val)
0795 {
0796 flagForceBuildTable = val;
0797 }
0798
0799
0800
0801 inline const G4String& G4VEmModel::GetName() const
0802 {
0803 return name;
0804 }
0805
0806
0807
0808 inline std::vector<G4EmElementSelector*>* G4VEmModel::GetElementSelectors()
0809 {
0810 return elmSelectors;
0811 }
0812
0813
0814
0815 inline void
0816 G4VEmModel::SetElementSelectors(std::vector<G4EmElementSelector*>* p)
0817 {
0818 if(p != elmSelectors) {
0819 elmSelectors = p;
0820 nSelectors = (nullptr != elmSelectors) ? G4int(elmSelectors->size()) : 0;
0821 localElmSelectors = false;
0822 }
0823 }
0824
0825
0826
0827 inline G4ElementData* G4VEmModel::GetElementData()
0828 {
0829 return fElementData;
0830 }
0831
0832
0833
0834 inline G4PhysicsTable* G4VEmModel::GetCrossSectionTable()
0835 {
0836 return xSectionTable;
0837 }
0838
0839
0840
0841 inline G4bool G4VEmModel::IsLocked() const
0842 {
0843 return isLocked;
0844 }
0845
0846
0847
0848 inline void G4VEmModel::SetLocked(G4bool val)
0849 {
0850 isLocked = val;
0851 }
0852
0853
0854
0855 #endif