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 #include "globals.hh"
00038 #include <iostream>
00039 #include "G4SampleResonance.hh"
00040 #include "G4DecayTable.hh"
00041 #include "Randomize.hh"
00042 #include "G4HadronicException.hh"
00043
00044 G4SampleResonance::minMassMapType G4SampleResonance::minMassCache;
00045
00046 G4double G4SampleResonance::GetMinimumMass(const G4ParticleDefinition* p) const
00047 {
00048
00049 G4double minResonanceMass = DBL_MAX;
00050
00051 if ( p->IsShortLived() )
00052 {
00053 minMassMapIterator iter = minMassCache.find(p);
00054 if ( iter!=minMassCache.end() )
00055 {
00056 minResonanceMass = (*iter).second;
00057 }
00058 else
00059 {
00060
00061
00062 const G4DecayTable* theDecays = const_cast<G4ParticleDefinition*>(p)->GetDecayTable();
00063 const G4int nDecays = theDecays->entries();
00064
00065 for (G4int i=0; i<nDecays; i++)
00066 {
00067 const G4VDecayChannel* aDecay = theDecays->GetDecayChannel(i);
00068 const G4int nDaughters = aDecay->GetNumberOfDaughters();
00069
00070 G4double minChannelMass = 0;
00071
00072 for (G4int j=0; j<nDaughters; j++)
00073 {
00074 const G4ParticleDefinition* aDaughter = const_cast<G4VDecayChannel*>(aDecay)->GetDaughter(j);
00075 G4double minMass = GetMinimumMass(aDaughter);
00076 if (!minMass) minMass = DBL_MAX;
00077 minChannelMass+=minMass;
00078 }
00079
00080 if (minChannelMass < minResonanceMass) minResonanceMass = minChannelMass;
00081
00082 }
00083
00084 G4SampleResonance* self = const_cast<G4SampleResonance*>(this);
00085 (self->minMassCache)[p] = minResonanceMass;
00086
00087 }
00088 }
00089 else
00090 {
00091
00092 minResonanceMass = p->GetPDGMass();
00093
00094 }
00095
00096
00097 return minResonanceMass;
00098 }
00099
00100
00101
00102 G4double G4SampleResonance::SampleMass(const G4ParticleDefinition* p, const G4double maxMass) const
00103 {
00104 return SampleMass(p->GetPDGMass(), p->GetPDGWidth(), GetMinimumMass(p), maxMass);
00105 }
00106
00107
00108 G4double G4SampleResonance::SampleMass(const G4double poleMass,
00109 const G4double gamma,
00110 const G4double minMass,
00111 const G4double maxMass) const
00112 {
00113
00114
00115
00116
00117 if ( minMass > maxMass )
00118 {
00119 throw G4HadronicException(__FILE__, __LINE__,
00120 "SampleResonanceMass: mass range negative (minMass>maxMass)");
00121 }
00122
00123 G4double returnMass;
00124
00125 if ( gamma < DBL_EPSILON )
00126 {
00127 returnMass = std::max(minMass, std::min(maxMass, poleMass));
00128 }
00129 else
00130 {
00131 double fmin = BrWigInt0(minMass, gamma, poleMass);
00132 double fmax = BrWigInt0(maxMass, gamma, poleMass);
00133 double f = fmin + (fmax-fmin)*G4UniformRand();
00134 returnMass = BrWigInv(f, gamma, poleMass);
00135 }
00136
00137 return returnMass;
00138 }
00139
00140
00141
00142