127{
128
129
130
131
132#ifdef G4VERBOSE
134#endif
135
138
139
141 const G4int N_DAUGHTER=3;
142
143
146 for (
G4int index=0; index<N_DAUGHTER; ++index)
147 {
149 sumofdaughtermass += daughtermass[index];
150 }
151
152
156
158 delete parentparticle;
159
160
161 G4double daughtermomentum[N_DAUGHTER];
162
163 G4double xmax = (1.0+daughtermass[0]*daughtermass[0]/parentmass/parentmass);
165
167
169 G4double EMax=parentmass/2-daughtermass[0];
170
171 const std::size_t MAX_LOOP=1000;
172
173 for (std::size_t loop1=0; loop1<MAX_LOOP; ++loop1)
174 {
176 for (std::size_t loop2 =0; loop2<MAX_LOOP; ++loop2)
177 {
180 if (gam <= x*(1.-x)) break;
181 x = xmax;
182 }
183 Ene=x;
184 if ( Ene >= (1.-Ee)) break;
185 Ene = 1.-Ee;
186 }
188
189
190
191 G4double costheta,sintheta,rphi,rtheta,rpsi;
192 costheta= 1.-2./Ee-2./Ene+2./Ene/Ee;
193 sintheta=std::sqrt(1.-costheta*costheta);
194
195
199
201 rot.
set(rphi,rtheta,rpsi);
202
203
204 daughtermomentum[0]=std::sqrt(Ee*Ee*EMax*EMax+2.0*Ee*EMax * daughtermass[0]);
206
207 direction0 *= rot;
208
211
213
214
215
216 daughtermomentum[1] = std::sqrt(Ene*Ene*EMax*EMax
217 +2.0*Ene*EMax*daughtermass[1]);
219
220 direction1 *= rot;
221
225
226
227
228 daughtermomentum[2] = std::sqrt(Enm*Enm*EMax*EMax
229 +2.0*Enm*EMax*daughtermass[2]);
230 G4ThreeVector direction2(-Ene/Enm*sintheta,0,-Ee/Enm-Ene/Enm*costheta);
231
232 direction2 *= rot;
233
237
238
239#ifdef G4VERBOSE
241 {
242 G4cout <<
"G4MuonDecayChannel::DecayIt()";
243 G4cout <<
" create decay products in rest frame " <<
G4endl;
245 }
246#endif
247 return products;
248}
HepRotation & set(const Hep3Vector &axis, double delta)
G4int PushProducts(G4DynamicParticle *aParticle)
G4double GetPDGMass() const
G4ParticleDefinition ** G4MT_daughters
void CheckAndFillParent()
G4ParticleDefinition * G4MT_parent
void CheckAndFillDaughters()