107{
108
109
110
111
112#ifdef G4VERBOSE
114#endif
115
118
119
121
122
125 for (
G4int index = 0; index < 3; ++index) {
127 sumofdaughtermass += daughtermass[index];
128 }
129 G4double xmax = parentmass - sumofdaughtermass;
130
131
134
135
137 delete parentparticle;
138
139
141
142
149 const std::size_t MAX_LOOP = 10000;
150 for (std::size_t loop_counter = 0; loop_counter < MAX_LOOP; ++loop_counter) {
152 p = std::sqrt(x * (x + 2.0 * dm));
154 r = p * (x + dm) * (xmax - x) * (xmax - x) * (1.0 +
aENuCorr * p / (x + dm) * w);
156 if (r > r0) break;
157 }
158
159
160
161
163 G4double theta = std::acos(costheta) * rad;
168
169
170 daughtermomentum[0] = p;
172 direction0 = rm * direction0;
173 auto daughterparticle0 =
175 products->PushProducts(daughterparticle0);
176
177
179 eNu = (parentmass - daughtermass[2]) * (parentmass + daughtermass[2]) + (dm * dm)
180 - 2. * parentmass * (x + dm);
181 eNu /= 2. * (parentmass + p * w - (x + dm));
184 G4double sinn = std::sqrt((1.0 - cosn) * (1.0 + cosn));
185
186 G4ThreeVector direction1(sinn * std::cos(phin), sinn * std::sin(phin), cosn);
187 direction1 = rm * direction1;
189 products->PushProducts(daughterparticle1);
190
191
193 eP = parentmass - eNu - (x + dm) - daughtermass[2];
196 G4double pP = std::sqrt(eP * (eP + 2. * daughtermass[2]));
197 G4ThreeVector direction2(pPx / pP * std::cos(phin), pPx / pP * std::sin(phin), pPz / pP);
198 direction2 = rm * direction2;
200 products->PushProducts(daughterparticle2);
201
202
203#ifdef G4VERBOSE
205 G4cout <<
"G4NeutronBetaDecayChannel::DecayIt ";
206 G4cout <<
" create decay products in rest frame " <<
G4endl;
207 products->DumpInfo();
208 }
209#endif
210 return products;
211}
HepRotation & rotateZ(double delta)
HepRotation & rotateY(double delta)
G4double GetPDGMass() const
G4ParticleDefinition ** G4MT_daughters
void CheckAndFillParent()
G4ParticleDefinition * G4MT_parent
void CheckAndFillDaughters()