00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077 #ifndef G4VEmModel_h
00078 #define G4VEmModel_h 1
00079
00080 #include "globals.hh"
00081 #include "G4DynamicParticle.hh"
00082 #include "G4ParticleDefinition.hh"
00083 #include "G4MaterialCutsCouple.hh"
00084 #include "G4Material.hh"
00085 #include "G4Element.hh"
00086 #include "G4ElementVector.hh"
00087 #include "G4DataVector.hh"
00088 #include "G4VEmFluctuationModel.hh"
00089 #include "G4VEmAngularDistribution.hh"
00090 #include "G4EmElementSelector.hh"
00091 #include "Randomize.hh"
00092 #include <vector>
00093
00094 class G4PhysicsTable;
00095 class G4Region;
00096 class G4VParticleChange;
00097 class G4ParticleChangeForLoss;
00098 class G4ParticleChangeForGamma;
00099 class G4Track;
00100
00101 class G4VEmModel
00102 {
00103
00104 public:
00105
00106 G4VEmModel(const G4String& nam);
00107
00108 virtual ~G4VEmModel();
00109
00110
00111
00112
00113
00114 virtual void Initialise(const G4ParticleDefinition*,
00115 const G4DataVector&) = 0;
00116
00117 virtual void SampleSecondaries(std::vector<G4DynamicParticle*>*,
00118 const G4MaterialCutsCouple*,
00119 const G4DynamicParticle*,
00120 G4double tmin = 0.0,
00121 G4double tmax = DBL_MAX) = 0;
00122
00123
00124
00125
00126
00127
00128 virtual G4double ComputeDEDXPerVolume(const G4Material*,
00129 const G4ParticleDefinition*,
00130 G4double kineticEnergy,
00131 G4double cutEnergy = DBL_MAX);
00132
00133
00134 virtual G4double CrossSectionPerVolume(const G4Material*,
00135 const G4ParticleDefinition*,
00136 G4double kineticEnergy,
00137 G4double cutEnergy = 0.0,
00138 G4double maxEnergy = DBL_MAX);
00139
00140
00141 virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
00142 G4double kinEnergy,
00143 G4double Z,
00144 G4double A = 0.,
00145 G4double cutEnergy = 0.0,
00146 G4double maxEnergy = DBL_MAX);
00147
00148
00149 virtual G4double ChargeSquareRatio(const G4Track&);
00150
00151
00152 virtual G4double GetChargeSquareRatio(const G4ParticleDefinition*,
00153 const G4Material*,
00154 G4double kineticEnergy);
00155
00156
00157 virtual G4double GetParticleCharge(const G4ParticleDefinition*,
00158 const G4Material*,
00159 G4double kineticEnergy);
00160
00161
00162 virtual void StartTracking(G4Track*);
00163
00164
00165 virtual void CorrectionsAlongStep(const G4MaterialCutsCouple*,
00166 const G4DynamicParticle*,
00167 G4double& eloss,
00168 G4double& niel,
00169 G4double length);
00170
00171
00172 virtual G4double Value(const G4MaterialCutsCouple*,
00173 const G4ParticleDefinition*,
00174 G4double kineticEnergy);
00175
00176
00177 virtual G4double MinPrimaryEnergy(const G4Material*,
00178 const G4ParticleDefinition*);
00179
00180
00181 virtual void SetupForMaterial(const G4ParticleDefinition*,
00182 const G4Material*,
00183 G4double kineticEnergy);
00184
00185
00186 virtual void DefineForRegion(const G4Region*);
00187
00188 protected:
00189
00190
00191 G4ParticleChangeForLoss* GetParticleChangeForLoss();
00192
00193
00194 G4ParticleChangeForGamma* GetParticleChangeForGamma();
00195
00196
00197 virtual G4double MaxSecondaryEnergy(const G4ParticleDefinition*,
00198 G4double kineticEnergy);
00199
00200 public:
00201
00202
00203
00204
00205
00206
00207 void InitialiseElementSelectors(const G4ParticleDefinition*,
00208 const G4DataVector&);
00209
00210
00211 inline G4double ComputeDEDX(const G4MaterialCutsCouple*,
00212 const G4ParticleDefinition*,
00213 G4double kineticEnergy,
00214 G4double cutEnergy = DBL_MAX);
00215
00216
00217 inline G4double CrossSection(const G4MaterialCutsCouple*,
00218 const G4ParticleDefinition*,
00219 G4double kineticEnergy,
00220 G4double cutEnergy = 0.0,
00221 G4double maxEnergy = DBL_MAX);
00222
00223
00224 inline G4double ComputeMeanFreePath(const G4ParticleDefinition*,
00225 G4double kineticEnergy,
00226 const G4Material*,
00227 G4double cutEnergy = 0.0,
00228 G4double maxEnergy = DBL_MAX);
00229
00230
00231 inline G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
00232 const G4Element*,
00233 G4double kinEnergy,
00234 G4double cutEnergy = 0.0,
00235 G4double maxEnergy = DBL_MAX);
00236
00237
00238 inline G4int SelectIsotopeNumber(const G4Element*);
00239
00240
00241 inline const G4Element* SelectRandomAtom(const G4MaterialCutsCouple*,
00242 const G4ParticleDefinition*,
00243 G4double kineticEnergy,
00244 G4double cutEnergy = 0.0,
00245 G4double maxEnergy = DBL_MAX);
00246
00247
00248 const G4Element* SelectRandomAtom(const G4Material*,
00249 const G4ParticleDefinition*,
00250 G4double kineticEnergy,
00251 G4double cutEnergy = 0.0,
00252 G4double maxEnergy = DBL_MAX);
00253
00254
00255
00256
00257
00258 void SetParticleChange(G4VParticleChange*, G4VEmFluctuationModel* f=0);
00259
00260 void SetCrossSectionTable(G4PhysicsTable*);
00261
00262 inline G4PhysicsTable* GetCrossSectionTable();
00263
00264 inline G4VEmFluctuationModel* GetModelOfFluctuations();
00265
00266 inline G4VEmAngularDistribution* GetAngularDistribution();
00267
00268 inline void SetAngularDistribution(G4VEmAngularDistribution*);
00269
00270 inline G4double HighEnergyLimit() const;
00271
00272 inline G4double LowEnergyLimit() const;
00273
00274 inline G4double HighEnergyActivationLimit() const;
00275
00276 inline G4double LowEnergyActivationLimit() const;
00277
00278 inline G4double PolarAngleLimit() const;
00279
00280 inline G4double SecondaryThreshold() const;
00281
00282 inline G4bool LPMFlag() const;
00283
00284 inline G4bool DeexcitationFlag() const;
00285
00286 inline G4bool ForceBuildTableFlag() const;
00287
00288 inline void SetHighEnergyLimit(G4double);
00289
00290 inline void SetLowEnergyLimit(G4double);
00291
00292 inline void SetActivationHighEnergyLimit(G4double);
00293
00294 inline void SetActivationLowEnergyLimit(G4double);
00295
00296 inline G4bool IsActive(G4double kinEnergy);
00297
00298 inline void SetPolarAngleLimit(G4double);
00299
00300 inline void SetSecondaryThreshold(G4double);
00301
00302 inline void SetLPMFlag(G4bool val);
00303
00304 inline void SetDeexcitationFlag(G4bool val);
00305
00306 inline void ForceBuildTable(G4bool val);
00307
00308 inline G4double MaxSecondaryKinEnergy(const G4DynamicParticle* dynParticle);
00309
00310 inline const G4String& GetName() const;
00311
00312 inline void SetCurrentCouple(const G4MaterialCutsCouple*);
00313
00314 inline const G4Element* GetCurrentElement() const;
00315
00316 protected:
00317
00318 inline const G4MaterialCutsCouple* CurrentCouple() const;
00319
00320 inline void SetCurrentElement(const G4Element*);
00321
00322 private:
00323
00324
00325 G4VEmModel & operator=(const G4VEmModel &right);
00326 G4VEmModel(const G4VEmModel&);
00327
00328
00329
00330 G4VEmFluctuationModel* flucModel;
00331 G4VEmAngularDistribution* anglModel;
00332 const G4String name;
00333
00334
00335
00336 G4double lowLimit;
00337 G4double highLimit;
00338 G4double eMinActive;
00339 G4double eMaxActive;
00340 G4double polarAngleLimit;
00341 G4double secondaryThreshold;
00342 G4bool theLPMflag;
00343 G4bool flagDeexcitation;
00344 G4bool flagForceBuildTable;
00345
00346 G4int nSelectors;
00347 std::vector<G4EmElementSelector*> elmSelectors;
00348
00349 protected:
00350
00351 G4VParticleChange* pParticleChange;
00352 G4PhysicsTable* xSectionTable;
00353 const std::vector<G4double>* theDensityFactor;
00354 const std::vector<G4int>* theDensityIdx;
00355
00356
00357
00358 private:
00359
00360 const G4MaterialCutsCouple* fCurrentCouple;
00361 const G4Element* fCurrentElement;
00362
00363 G4int nsec;
00364 std::vector<G4double> xsec;
00365
00366 };
00367
00368
00369
00370 inline void G4VEmModel::SetCurrentCouple(const G4MaterialCutsCouple* p)
00371 {
00372 fCurrentCouple = p;
00373 }
00374
00375
00376
00377 inline const G4MaterialCutsCouple* G4VEmModel::CurrentCouple() const
00378 {
00379 return fCurrentCouple;
00380 }
00381
00382
00383
00384 inline void G4VEmModel::SetCurrentElement(const G4Element* elm)
00385 {
00386 fCurrentElement = elm;
00387 }
00388
00389
00390
00391 inline const G4Element* G4VEmModel::GetCurrentElement() const
00392 {
00393 return fCurrentElement;
00394 }
00395
00396
00397
00398 inline
00399 G4double G4VEmModel::MaxSecondaryKinEnergy(const G4DynamicParticle* dynPart)
00400 {
00401 return MaxSecondaryEnergy(dynPart->GetParticleDefinition(),
00402 dynPart->GetKineticEnergy());
00403 }
00404
00405
00406
00407 inline G4double G4VEmModel::ComputeDEDX(const G4MaterialCutsCouple* c,
00408 const G4ParticleDefinition* p,
00409 G4double kinEnergy,
00410 G4double cutEnergy)
00411 {
00412 fCurrentCouple = c;
00413 return ComputeDEDXPerVolume(c->GetMaterial(),p,kinEnergy,cutEnergy);
00414 }
00415
00416
00417
00418 inline G4double G4VEmModel::CrossSection(const G4MaterialCutsCouple* c,
00419 const G4ParticleDefinition* p,
00420 G4double kinEnergy,
00421 G4double cutEnergy,
00422 G4double maxEnergy)
00423 {
00424 fCurrentCouple = c;
00425 return CrossSectionPerVolume(c->GetMaterial(),p,kinEnergy,cutEnergy,maxEnergy);
00426 }
00427
00428
00429
00430 inline G4double G4VEmModel::ComputeMeanFreePath(const G4ParticleDefinition* p,
00431 G4double ekin,
00432 const G4Material* material,
00433 G4double emin,
00434 G4double emax)
00435 {
00436 G4double mfp = DBL_MAX;
00437 G4double cross = CrossSectionPerVolume(material,p,ekin,emin,emax);
00438 if (cross > DBL_MIN) { mfp = 1./cross; }
00439 return mfp;
00440 }
00441
00442
00443
00444 inline G4double G4VEmModel::ComputeCrossSectionPerAtom(
00445 const G4ParticleDefinition* part,
00446 const G4Element* elm,
00447 G4double kinEnergy,
00448 G4double cutEnergy,
00449 G4double maxEnergy)
00450 {
00451 fCurrentElement = elm;
00452 return ComputeCrossSectionPerAtom(part,kinEnergy,elm->GetZ(),elm->GetN(),
00453 cutEnergy,maxEnergy);
00454 }
00455
00456
00457
00458 inline const G4Element*
00459 G4VEmModel::SelectRandomAtom(const G4MaterialCutsCouple* couple,
00460 const G4ParticleDefinition* p,
00461 G4double kinEnergy,
00462 G4double cutEnergy,
00463 G4double maxEnergy)
00464 {
00465 fCurrentCouple = couple;
00466 if(nSelectors > 0) {
00467 fCurrentElement =
00468 elmSelectors[couple->GetIndex()]->SelectRandomAtom(kinEnergy);
00469 } else {
00470 fCurrentElement = SelectRandomAtom(couple->GetMaterial(),p,kinEnergy,
00471 cutEnergy,maxEnergy);
00472 }
00473 return fCurrentElement;
00474 }
00475
00476
00477
00478 inline G4int G4VEmModel::SelectIsotopeNumber(const G4Element* elm)
00479 {
00480 fCurrentElement = elm;
00481 G4int N = G4int(elm->GetN() + 0.5);
00482 G4int ni = elm->GetNumberOfIsotopes();
00483 if(ni > 0) {
00484 G4int idx = 0;
00485 if(ni > 1) {
00486 G4double* ab = elm->GetRelativeAbundanceVector();
00487 G4double x = G4UniformRand();
00488 for(; idx<ni; ++idx) {
00489 x -= ab[idx];
00490 if (x <= 0.0) { break; }
00491 }
00492 if(idx >= ni) { idx = ni - 1; }
00493 }
00494 N = elm->GetIsotope(idx)->GetN();
00495 }
00496 return N;
00497 }
00498
00499
00500
00501 inline G4VEmFluctuationModel* G4VEmModel::GetModelOfFluctuations()
00502 {
00503 return flucModel;
00504 }
00505
00506
00507
00508 inline G4VEmAngularDistribution* G4VEmModel::GetAngularDistribution()
00509 {
00510 return anglModel;
00511 }
00512
00513
00514
00515 inline void G4VEmModel::SetAngularDistribution(G4VEmAngularDistribution* p)
00516 {
00517 anglModel = p;
00518 }
00519
00520
00521
00522 inline G4double G4VEmModel::HighEnergyLimit() const
00523 {
00524 return highLimit;
00525 }
00526
00527
00528
00529 inline G4double G4VEmModel::LowEnergyLimit() const
00530 {
00531 return lowLimit;
00532 }
00533
00534
00535
00536 inline G4double G4VEmModel::HighEnergyActivationLimit() const
00537 {
00538 return eMaxActive;
00539 }
00540
00541
00542
00543 inline G4double G4VEmModel::LowEnergyActivationLimit() const
00544 {
00545 return eMinActive;
00546 }
00547
00548
00549
00550 inline G4double G4VEmModel::PolarAngleLimit() const
00551 {
00552 return polarAngleLimit;
00553 }
00554
00555
00556
00557 inline G4double G4VEmModel::SecondaryThreshold() const
00558 {
00559 return secondaryThreshold;
00560 }
00561
00562
00563
00564 inline G4bool G4VEmModel::LPMFlag() const
00565 {
00566 return theLPMflag;
00567 }
00568
00569
00570
00571 inline G4bool G4VEmModel::DeexcitationFlag() const
00572 {
00573 return flagDeexcitation;
00574 }
00575
00576
00577
00578 inline G4bool G4VEmModel::ForceBuildTableFlag() const
00579 {
00580 return flagForceBuildTable;
00581 }
00582
00583
00584
00585 inline void G4VEmModel::SetHighEnergyLimit(G4double val)
00586 {
00587 highLimit = val;
00588 }
00589
00590
00591
00592 inline void G4VEmModel::SetLowEnergyLimit(G4double val)
00593 {
00594 lowLimit = val;
00595 }
00596
00597
00598
00599 inline void G4VEmModel::SetActivationHighEnergyLimit(G4double val)
00600 {
00601 eMaxActive = val;
00602 }
00603
00604
00605
00606 inline void G4VEmModel::SetActivationLowEnergyLimit(G4double val)
00607 {
00608 eMinActive = val;
00609 }
00610
00611
00612
00613 inline G4bool G4VEmModel::IsActive(G4double kinEnergy)
00614 {
00615 return (kinEnergy >= eMinActive && kinEnergy <= eMaxActive);
00616 }
00617
00618
00619
00620 inline void G4VEmModel::SetPolarAngleLimit(G4double val)
00621 {
00622 polarAngleLimit = val;
00623 }
00624
00625
00626
00627 inline void G4VEmModel::SetSecondaryThreshold(G4double val)
00628 {
00629 secondaryThreshold = val;
00630 }
00631
00632
00633
00634 inline void G4VEmModel::SetLPMFlag(G4bool val)
00635 {
00636 theLPMflag = val;
00637 }
00638
00639
00640
00641 inline void G4VEmModel::SetDeexcitationFlag(G4bool val)
00642 {
00643 flagDeexcitation = val;
00644 }
00645
00646
00647
00648 inline void G4VEmModel::ForceBuildTable(G4bool val)
00649 {
00650 flagForceBuildTable = val;
00651 }
00652
00653
00654
00655 inline const G4String& G4VEmModel::GetName() const
00656 {
00657 return name;
00658 }
00659
00660 inline G4PhysicsTable* G4VEmModel::GetCrossSectionTable()
00661 {
00662 return xSectionTable;
00663 }
00664
00665
00666
00667 #endif
00668