120 for (
G4int index=0; index<3; index++){
122 sumofdaughtermass += daughtermass[index];
132 delete parentparticle;
146 G4double W_mue = (EMMU*EMMU+EMASS*EMASS)/(2.*EMMU);
159 const size_t MAX_LOOP=10000;
160 for (
size_t loop_count =0; loop_count <MAX_LOOP; ++loop_count){
166 x = x0 + rndm*(1.-x0);
172 F_IS = 1./6.*(-2.*x_squared+3.*x-x0_squared);
173 F_AS = 1./6.*std::sqrt(x_squared-x0_squared)*(2.*x-2.+std::sqrt(1.-x0_squared));
175 G_IS = 2./9.*(michel_rho-0.75)*(4.*x_squared-3.*x-x0_squared);
176 G_IS = G_IS + michel_eta*(1.-
x)*x0;
178 G_AS = 3.*(michel_xsi-1.)*(1.-
x);
179 G_AS = G_AS+2.*(michel_xsi*michel_delta-0.75)*(4.*x-4.+std::sqrt(1.-x0_squared));
180 G_AS = 1./9.*std::sqrt(x_squared-x0_squared)*G_AS;
186 const G4double omega = std::log(EMMU/EMASS);
189 G4double F = 6.*F_IS + R_IS/std::sqrt(x_squared-x0_squared);
193 G4double R_AS = F_theta(x,x0,omega);
199 G4double G = 6.*F_AS - R_AS/std::sqrt(x_squared-x0_squared);
201 FG = std::sqrt(x_squared-x0_squared)*F*(1.+(G/F)*ctheta);
204 G4cout<<
"***Problem in Muon Decay *** : FG > FG_max"<<
G4endl;
210 if (FG >= rndm*FG_max)
break;
219 if(energy < EMASS) energy = EMASS;
224 daughtermomentum[0] = std::sqrt(energy*energy - EMASS*EMASS);
226 G4double stheta = std::sqrt(1.-ctheta*ctheta);
248 G4double energy2 = parentmass*(1.0 - x/2.0);
249 G4double vmass = std::sqrt((energy2-daughtermomentum[0])*(energy2+daughtermomentum[0]));
250 G4double beta = -1.0*daughtermomentum[0]/energy2;
252 G4double sinthetan = std::sqrt((1.0-costhetan)*(1.0+costhetan));
257 G4ThreeVector direction1(sinthetan*cosphin,sinthetan*sinphin,costhetan);
266 p4.
boost( direction0.x()*beta, direction0.y()*beta, direction0.z()*beta);
269 p4.
boost( direction0.x()*beta, direction0.y()*beta, direction0.z()*beta);
279 G4cout <<
"G4MuonDecayChannelWithSpin::DecayIt ";
280 G4cout <<
" create decay products in rest frame " <<
G4endl;
void CheckAndFillDaughters()
G4int PushProducts(G4DynamicParticle *aParticle)
G4ParticleDefinition * G4MT_parent
G4ParticleDefinition ** G4MT_daughters
static constexpr double rad
static constexpr double twopi
G4double GetTotalMomentum() const
G4GLOB_DLL std::ostream G4cout
HepLorentzVector & boost(double, double, double)
G4ThreeVector parent_polarization
G4LorentzVector Get4Momentum() const
void Set4Momentum(const G4LorentzVector &momentum)
G4int GetVerboseLevel() const
G4double GetPDGMass() const
G4double energy(const ThreeVector &p, const G4double m)
void CheckAndFillParent()