| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27 //
28 // -------------------------------------------
29 // GEANT 4 class header file
30 //
31 // ---------------- G4SampleResonance ---
32 // by Henning Weber, March 2001.
33 // helper class for sampling resonance ma
34 // -------------------------------------------
35
36
37 #include "globals.hh"
38 #include <iostream>
39 #include "G4SampleResonance.hh"
40 #include "G4DecayTable.hh"
41 #include "Randomize.hh"
42 #include "G4HadronicException.hh"
43
44 G4ThreadLocal G4SampleResonance::minMassMapTyp
45
46 G4double G4SampleResonance::GetMinimumMass(con
47 { ;;; if (!minMassCache_G4MT_TLS_) minMassC
48
49 G4double minResonanceMass = DBL_MAX;
50
51 if ( p->IsShortLived() )
52 {
53 minMassMapIterator iter = minMassCache.fin
54 if ( iter!=minMassCache.end() )
55 {
56 minResonanceMass = (*iter).second;
57 }
58 else
59 {
60 // G4cout << "--> request for " << p->Ge
61
62 const G4DecayTable* theDecays = p->GetDe
63 const G4int nDecays = theDecays->entries
64
65 // To find the minimum mass of the reson
66 // decay channels whose branching ratio
67 // This is needed to avoid that rare and
68 // (e.g. e+ e-) can set a very small min
69 // In the case that no channel with bran
70 // threshold has been found, consider th
71 // branching ratio (whatever its values)
72 // Note that this solution works also wh
73 // enhanced if both of the f
74 // 1. The enhanced rar
75 // 2. The decay with t
76 // i.e. not a rare decay which has be
77 const G4double thresholdChannelProbabili
78 G4double foundChannelAboveThresholdProba
79 G4double minMassMostProbableChannel = 0.
80 G4double highestChannelProbability = 0.0
81 for (G4int i=0; i<nDecays; i++)
82 {
83 const G4VDecayChannel* aDecay = theDec
84 G4double decayBr = aDecay->GetBR();
85 if (decayBr > std::min(highestChannelP
86 {
87 const G4int nDaughters = aDecay->Get
88 G4double minChannelMass = 0;
89 for (G4int j=0; j<nDaughters; j++)
90 {
91 const G4ParticleDefinition* aDaughte
92 G4double minMass = GetMinimumMass(aD
93 if (!minMass) minMass = DBL_MAX; //
94 minChannelMass+=minMass;
95 }
96 if (decayBr > highestChannelProbabil
97 {
98 highestChannelProbability = decayB
99 minMassMostProbableChannel = minCh
100 }
101 if (decayBr > thresholdChannelProbab
102 {
103 foundChannelAboveThresholdProbabil
104 if (minChannelMass < minResonanceM
105 }
106 }
107 }
108 if ( ! foundChannelAboveThresholdProbabi
109 minResonanceMass = minMassMostProbable
110 }
111 // replace this as soon as the compiler
112 G4SampleResonance* self = const_cast<G4S
113 //Andrea Dotti (13Jan2013): Change neede
114 //(self->minMassCache)[p] = minResonance
115 self->minMassCache_G4MT_TLS_->operator[]
116
117 }
118 }
119 else
120 {
121
122 minResonanceMass = p->GetPDGMass();
123
124 }
125 // G4cout << "minimal mass for " << p->GetPa
126
127 return minResonanceMass;
128 }
129
130
131
132 G4double G4SampleResonance::SampleMass(const G
133 { if (!minMassCache_G4MT_TLS_) minMassCache_G4
134 return SampleMass(p->GetPDGMass(), p->GetPDG
135 }
136
137
138 G4double G4SampleResonance::SampleMass(const G
139 const G4double gamma,
140 const G4double minMass,
141 const G4double maxMass) const
142 { if (!minMassCache_G4MT_TLS_) minMassCache_G4
143 // Chooses a mass randomly between minMass a
144 // according to a Breit-Wigner function
145 // width gamma and pole poleMass
146
147
148 //AR-14Nov2017 : protection for rare c
149 // dynamic mass, decays
150 // then that for the dau
151 // do not crash, but sim
152 // the sampling of the d
153 G4double protectedMinMass = minMass;
154 if ( minMass > maxMass )
155 {
156 //throw G4HadronicException(__FILE__, __LI
157 // "SampleResonanceMass: mass range
158 protectedMinMass = maxMass;
159 }
160
161 G4double returnMass;
162
163 if ( gamma < DBL_EPSILON )
164 {
165 returnMass = std::max(minMass, std::min(ma
166 }
167 else
168 {
169 //double fmin = BrWigInt0(minMass, gamma,
170 double fmin = BrWigInt0(protectedMinMass,
171 double fmax = BrWigInt0(maxMass, gamma, po
172 double f = fmin + (fmax-fmin)*G4UniformRan
173 returnMass = BrWigInv(f, gamma, poleMass);
174 }
175
176 return returnMass;
177 }
178
179
180
181
182