| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // G4PhaseSpaceDecayChannel class implementation
27 //
28 // Author: H.Kurashige, 27 July 1996
29 // --------------------------------------------------------------------
30
31 #include "G4PhaseSpaceDecayChannel.hh"
32
33 #include "G4DecayProducts.hh"
34 #include "G4LorentzRotation.hh"
35 #include "G4LorentzVector.hh"
36 #include "G4ParticleDefinition.hh"
37 #include "G4PhysicalConstants.hh"
38 #include "G4SystemOfUnits.hh"
39 #include "G4VDecayChannel.hh"
40 #include "Randomize.hh"
41
42 G4PhaseSpaceDecayChannel::G4PhaseSpaceDecayChannel(G4int Verbose)
43 : G4VDecayChannel("Phase Space", Verbose)
44 {}
45
46 G4PhaseSpaceDecayChannel::G4PhaseSpaceDecayChannel(const G4String& theParentName, G4double theBR,
47 G4int theNumberOfDaughters,
48 const G4String& theDaughterName1,
49 const G4String& theDaughterName2,
50 const G4String& theDaughterName3,
51 const G4String& theDaughterName4,
52 const G4String& theDaughterName5)
53 : G4VDecayChannel("Phase Space", theParentName, theBR, theNumberOfDaughters, theDaughterName1,
54 theDaughterName2, theDaughterName3, theDaughterName4, theDaughterName5)
55 {}
56
57 G4DecayProducts* G4PhaseSpaceDecayChannel::DecayIt(G4double parentMass)
58 {
59 #ifdef G4VERBOSE
60 if (GetVerboseLevel() > 1) G4cout << "G4PhaseSpaceDecayChannel::DecayIt()" << G4endl;
61 #endif
62
63 G4DecayProducts* products = nullptr;
64
65 CheckAndFillParent();
66 CheckAndFillDaughters();
67
68 if (parentMass > 0.0)
69 current_parent_mass.Put(parentMass);
70 else
71 current_parent_mass.Put(G4MT_parent_mass);
72
73 switch (numberOfDaughters) {
74 case 0:
75 #ifdef G4VERBOSE
76 if (GetVerboseLevel() > 0) {
77 G4cout << "G4PhaseSpaceDecayChannel::DecayIt() -";
78 G4cout << " daughters not defined " << G4endl;
79 }
80 #endif
81 break;
82 case 1:
83 products = OneBodyDecayIt();
84 break;
85 case 2:
86 products = TwoBodyDecayIt();
87 break;
88 case 3:
89 products = ThreeBodyDecayIt();
90 break;
91 default:
92 products = ManyBodyDecayIt();
93 break;
94 }
95 #ifdef G4VERBOSE
96 if ((products == nullptr) && (GetVerboseLevel() > 0)) {
97 G4cout << "G4PhaseSpaceDecayChannel::DecayIt() - ";
98 G4cout << *parent_name << " cannot decay " << G4endl;
99 DumpInfo();
100 }
101 #endif
102 return products;
103 }
104
105 G4DecayProducts* G4PhaseSpaceDecayChannel::OneBodyDecayIt()
106 {
107 #ifdef G4VERBOSE
108 if (GetVerboseLevel() > 1) G4cout << "G4PhaseSpaceDecayChannel::OneBodyDecayIt()" << G4endl;
109 #endif
110 // parent mass
111 G4double parentmass = current_parent_mass.Get();
112
113 // create parent G4DynamicParticle at rest
114 G4ThreeVector dummy;
115 auto parentparticle = new G4DynamicParticle(G4MT_parent, dummy, 0.0, parentmass);
116 // create G4Decayproducts
117 auto products = new G4DecayProducts(*parentparticle);
118 delete parentparticle;
119
120 // create daughter G4DynamicParticle at rest
121 auto daughterparticle = new G4DynamicParticle(G4MT_daughters[0], dummy, 0.0);
122 if (useGivenDaughterMass) daughterparticle->SetMass(givenDaughterMasses[0]);
123 products->PushProducts(daughterparticle);
124
125 #ifdef G4VERBOSE
126 if (GetVerboseLevel() > 1) {
127 G4cout << "G4PhaseSpaceDecayChannel::OneBodyDecayIt() -";
128 G4cout << " create decay products in rest frame " << G4endl;
129 products->DumpInfo();
130 }
131 #endif
132 return products;
133 }
134
135 G4DecayProducts* G4PhaseSpaceDecayChannel::TwoBodyDecayIt()
136 {
137 #ifdef G4VERBOSE
138 if (GetVerboseLevel() > 1) G4cout << "G4PhaseSpaceDecayChannel::TwoBodyDecayIt()" << G4endl;
139 #endif
140 // parent mass
141 G4double parentmass = current_parent_mass.Get();
142
143 // daughters'mass, width
144 G4double daughtermass[2], daughterwidth[2];
145 daughtermass[0] = G4MT_daughters_mass[0];
146 daughtermass[1] = G4MT_daughters_mass[1];
147 daughterwidth[0] = G4MT_daughters_width[0];
148 daughterwidth[1] = G4MT_daughters_width[1];
149
150 // create parent G4DynamicParticle at rest
151 G4ThreeVector dummy;
152 auto parentparticle = new G4DynamicParticle(G4MT_parent, dummy, 0.0, parentmass);
153 // create G4Decayproducts
154 auto products = new G4DecayProducts(*parentparticle);
155 delete parentparticle;
156
157 if (!useGivenDaughterMass) {
158 G4bool withWidth = (daughterwidth[0] > 1.0e-3 * daughtermass[0])
159 || (daughterwidth[1] > 1.0e-3 * daughtermass[1]);
160 if (withWidth) {
161 G4double sumofdaughterwidthsq =
162 daughterwidth[0] * daughterwidth[0] + daughterwidth[1] * daughterwidth[1];
163 // check parent mass and daughter mass
164 G4double maxDev =
165 (parentmass - daughtermass[0] - daughtermass[1]) / std::sqrt(sumofdaughterwidthsq);
166 if (maxDev <= -1.0 * rangeMass) {
167 #ifdef G4VERBOSE
168 if (GetVerboseLevel() > 0) {
169 G4cout << "G4PhaseSpaceDecayChannel::TwoBodyDecayIt()" << G4endl
170 << "Sum of daughter mass is larger than parent mass!" << G4endl;
171 G4cout << "Parent :" << G4MT_parent->GetParticleName() << " "
172 << current_parent_mass.Get() / GeV << G4endl;
173 G4cout << "Daughter 1 :" << G4MT_daughters[0]->GetParticleName() << " "
174 << daughtermass[0] / GeV << G4endl;
175 G4cout << "Daughter 2:" << G4MT_daughters[1]->GetParticleName() << " "
176 << daughtermass[1] / GeV << G4endl;
177 }
178 #endif
179 G4Exception("G4PhaseSpaceDecayChannel::TwoBodyDecayIt()", "PART112", JustWarning,
180 "Cannot create decay products: sum of daughter mass is \
181 larger than parent mass!");
182 return products;
183 }
184 G4double dm1 = daughtermass[0];
185 if (daughterwidth[0] > 0.) dm1 = DynamicalMass(daughtermass[0], daughterwidth[0], maxDev);
186 G4double dm2 = daughtermass[1];
187 if (daughterwidth[1] > 0.) dm2 = DynamicalMass(daughtermass[1], daughterwidth[1], maxDev);
188 while (dm1 + dm2 > parentmass) // Loop checking, 09.08.2015, K.Kurashige
189 {
190 dm1 = DynamicalMass(daughtermass[0], daughterwidth[0], maxDev);
191 dm2 = DynamicalMass(daughtermass[1], daughterwidth[1], maxDev);
192 }
193 daughtermass[0] = dm1;
194 daughtermass[1] = dm2;
195 }
196 }
197 else {
198 // use given daughter mass;
199 daughtermass[0] = givenDaughterMasses[0];
200 daughtermass[1] = givenDaughterMasses[1];
201 }
202 if (parentmass < daughtermass[0] + daughtermass[1]) {
203 #ifdef G4VERBOSE
204 if (GetVerboseLevel() > 0) {
205 G4cout << "G4PhaseSpaceDecayChannel::TwoBodyDecayIt()" << G4endl
206 << "Sum of daughter mass is larger than parent mass!" << G4endl;
207 G4cout << "Parent :" << G4MT_parent->GetParticleName() << " "
208 << current_parent_mass.Get() / GeV << G4endl;
209 G4cout << "Daughter 1 :" << G4MT_daughters[0]->GetParticleName() << " "
210 << daughtermass[0] / GeV << G4endl;
211 G4cout << "Daughter 2:" << G4MT_daughters[1]->GetParticleName() << " "
212 << daughtermass[1] / GeV << G4endl;
213 if (useGivenDaughterMass) {
214 G4cout << "Daughter Mass is given." << G4endl;
215 }
216 }
217 #endif
218 G4Exception("G4PhaseSpaceDecayChannel::TwoBodyDecayIt()", "PART112", JustWarning,
219 "Cannot create decay products: sum of daughter mass is \
220 larger than parent mass!");
221 return products;
222 }
223
224 // calculate daughter momentum
225 G4double daughtermomentum = Pmx(parentmass, daughtermass[0], daughtermass[1]);
226
227 G4double costheta = 2. * G4UniformRand() - 1.0;
228 G4double sintheta = std::sqrt((1.0 - costheta) * (1.0 + costheta));
229 G4double phi = twopi * G4UniformRand() * rad;
230 G4ThreeVector direction(sintheta * std::cos(phi), sintheta * std::sin(phi), costheta);
231
232 // create daughter G4DynamicParticle
233 G4double Ekin = std::sqrt(daughtermomentum * daughtermomentum + daughtermass[0] * daughtermass[0])
234 - daughtermass[0];
235 auto daughterparticle =
236 new G4DynamicParticle(G4MT_daughters[0], direction, Ekin, daughtermass[0]);
237 products->PushProducts(daughterparticle);
238 Ekin = std::sqrt(daughtermomentum * daughtermomentum + daughtermass[1] * daughtermass[1])
239 - daughtermass[1];
240 daughterparticle =
241 new G4DynamicParticle(G4MT_daughters[1], -1.0 * direction, Ekin, daughtermass[1]);
242 products->PushProducts(daughterparticle);
243
244 #ifdef G4VERBOSE
245 if (GetVerboseLevel() > 1) {
246 G4cout << "G4PhaseSpaceDecayChannel::TwoBodyDecayIt() -";
247 G4cout << " Create decay products in rest frame " << G4endl;
248 products->DumpInfo();
249 }
250 #endif
251 return products;
252 }
253
254 G4DecayProducts* G4PhaseSpaceDecayChannel::ThreeBodyDecayIt()
255 {
256 // Algorithm of this code originally written in GDECA3 of GEANT3
257
258 #ifdef G4VERBOSE
259 if (GetVerboseLevel() > 1) G4cout << "G4PhaseSpaceDecayChannel::ThreeBodyDecayIt()" << G4endl;
260 #endif
261 // parent mass
262 G4double parentmass = current_parent_mass.Get();
263 // daughters'mass
264 G4double daughtermass[3], daughterwidth[3];
265 G4double sumofdaughtermass = 0.0;
266 G4double sumofdaughterwidthsq = 0.0;
267 G4bool withWidth = false;
268 for (G4int index = 0; index < 3; ++index) {
269 daughtermass[index] = G4MT_daughters_mass[index];
270 sumofdaughtermass += daughtermass[index];
271 daughterwidth[index] = G4MT_daughters_width[index];
272 sumofdaughterwidthsq += daughterwidth[index] * daughterwidth[index];
273 withWidth = withWidth || (daughterwidth[index] > 1.0e-3 * daughtermass[index]);
274 }
275
276 // create parent G4DynamicParticle at rest
277 G4ThreeVector dummy;
278 auto parentparticle = new G4DynamicParticle(G4MT_parent, dummy, 0.0, parentmass);
279
280 // create G4Decayproducts
281 auto products = new G4DecayProducts(*parentparticle);
282 delete parentparticle;
283
284 if (!useGivenDaughterMass) {
285 if (withWidth) {
286 G4double maxDev = (parentmass - sumofdaughtermass) / std::sqrt(sumofdaughterwidthsq);
287 if (maxDev <= -1.0 * rangeMass) {
288 #ifdef G4VERBOSE
289 if (GetVerboseLevel() > 0) {
290 G4cout << "G4PhaseSpaceDecayChannel::ThreeBodyDecayIt()" << G4endl
291 << "Sum of daughter mass is larger than parent mass!" << G4endl;
292 G4cout << "Parent :" << G4MT_parent->GetParticleName() << " "
293 << current_parent_mass.Get() / GeV << G4endl;
294 G4cout << "Daughter 1 :" << G4MT_daughters[0]->GetParticleName() << " "
295 << daughtermass[0] / GeV << G4endl;
296 G4cout << "Daughter 2:" << G4MT_daughters[1]->GetParticleName() << " "
297 << daughtermass[1] / GeV << G4endl;
298 G4cout << "Daughter 3:" << G4MT_daughters[2]->GetParticleName() << " "
299 << daughtermass[2] / GeV << G4endl;
300 }
301 #endif
302 G4Exception("G4PhaseSpaceDecayChannel::ThreeBodyDecayIt()", "PART112", JustWarning,
303 "Cannot create decay products: sum of daughter mass \
304 is larger than parent mass!");
305 return products;
306 }
307 G4double dm1 = daughtermass[0];
308 if (daughterwidth[0] > 0.) dm1 = DynamicalMass(daughtermass[0], daughterwidth[0], maxDev);
309 G4double dm2 = daughtermass[1];
310 if (daughterwidth[1] > 0.) dm2 = DynamicalMass(daughtermass[1], daughterwidth[1], maxDev);
311 G4double dm3 = daughtermass[2];
312 if (daughterwidth[2] > 0.) dm3 = DynamicalMass(daughtermass[2], daughterwidth[2], maxDev);
313 while (dm1 + dm2 + dm3 > parentmass) // Loop checking, 09.08.2015, K.Kurashige
314 {
315 dm1 = DynamicalMass(daughtermass[0], daughterwidth[0], maxDev);
316 dm2 = DynamicalMass(daughtermass[1], daughterwidth[1], maxDev);
317 dm3 = DynamicalMass(daughtermass[2], daughterwidth[2], maxDev);
318 }
319 daughtermass[0] = dm1;
320 daughtermass[1] = dm2;
321 daughtermass[2] = dm3;
322 sumofdaughtermass = dm1 + dm2 + dm3;
323 }
324 }
325 else {
326 // use given daughter mass;
327 daughtermass[0] = givenDaughterMasses[0];
328 daughtermass[1] = givenDaughterMasses[1];
329 daughtermass[2] = givenDaughterMasses[2];
330 sumofdaughtermass = daughtermass[0] + daughtermass[1] + daughtermass[2];
331 }
332
333 if (sumofdaughtermass > parentmass) {
334 #ifdef G4VERBOSE
335 if (GetVerboseLevel() > 0) {
336 G4cout << "G4PhaseSpaceDecayChannel::ThreeBodyDecayIt()" << G4endl
337 << "Sum of daughter mass is larger than parent mass!" << G4endl;
338 G4cout << "Parent :" << G4MT_parent->GetParticleName() << " "
339 << current_parent_mass.Get() / GeV << G4endl;
340 G4cout << "Daughter 1 :" << G4MT_daughters[0]->GetParticleName() << " "
341 << daughtermass[0] / GeV << G4endl;
342 G4cout << "Daughter 2:" << G4MT_daughters[1]->GetParticleName() << " "
343 << daughtermass[1] / GeV << G4endl;
344 G4cout << "Daughter 3:" << G4MT_daughters[2]->GetParticleName() << " "
345 << daughtermass[2] / GeV << G4endl;
346 }
347 if (useGivenDaughterMass) {
348 G4cout << "Daughter Mass is given." << G4endl;
349 }
350 #endif
351 G4Exception("G4PhaseSpaceDecayChannel::ThreeBodyDecayIt", "PART112", JustWarning,
352 "Can not create decay products: sum of daughter mass \
353 is larger than parent mass!");
354 return products;
355 }
356
357 // calculate daughter momentum
358 // Generate two
359 G4double rd1, rd2, rd;
360 G4double daughtermomentum[3];
361 G4double momentummax = 0.0, momentumsum = 0.0;
362 G4double energy;
363 const std::size_t MAX_LOOP = 10000;
364
365 for (std::size_t loop_counter = 0; loop_counter < MAX_LOOP; ++loop_counter) {
366 rd1 = G4UniformRand();
367 rd2 = G4UniformRand();
368 if (rd2 > rd1) {
369 rd = rd1;
370 rd1 = rd2;
371 rd2 = rd;
372 }
373 momentummax = 0.0;
374 momentumsum = 0.0;
375 // daughter 0
376 energy = rd2 * (parentmass - sumofdaughtermass);
377 daughtermomentum[0] = std::sqrt(energy * energy + 2.0 * energy * daughtermass[0]);
378 if (daughtermomentum[0] > momentummax) momentummax = daughtermomentum[0];
379 momentumsum += daughtermomentum[0];
380 // daughter 1
381 energy = (1. - rd1) * (parentmass - sumofdaughtermass);
382 daughtermomentum[1] = std::sqrt(energy * energy + 2.0 * energy * daughtermass[1]);
383 if (daughtermomentum[1] > momentummax) momentummax = daughtermomentum[1];
384 momentumsum += daughtermomentum[1];
385 // daughter 2
386 energy = (rd1 - rd2) * (parentmass - sumofdaughtermass);
387 daughtermomentum[2] = std::sqrt(energy * energy + 2.0 * energy * daughtermass[2]);
388 if (daughtermomentum[2] > momentummax) momentummax = daughtermomentum[2];
389 momentumsum += daughtermomentum[2];
390 if (momentummax <= momentumsum - momentummax) break;
391 }
392
393 // output message
394 #ifdef G4VERBOSE
395 if (GetVerboseLevel() > 1) {
396 G4cout << " daughter 0:" << daughtermomentum[0] / GeV << "[GeV/c]" << G4endl;
397 G4cout << " daughter 1:" << daughtermomentum[1] / GeV << "[GeV/c]" << G4endl;
398 G4cout << " daughter 2:" << daughtermomentum[2] / GeV << "[GeV/c]" << G4endl;
399 G4cout << " momentum sum:" << momentumsum / GeV << "[GeV/c]" << G4endl;
400 }
401 #endif
402
403 // create daughter G4DynamicParticle
404 G4double costheta, sintheta, phi, sinphi, cosphi;
405 G4double costhetan, sinthetan, phin, sinphin, cosphin;
406 costheta = 2. * G4UniformRand() - 1.0;
407 sintheta = std::sqrt((1.0 - costheta) * (1.0 + costheta));
408 phi = twopi * G4UniformRand() * rad;
409 sinphi = std::sin(phi);
410 cosphi = std::cos(phi);
411
412 G4ThreeVector direction0(sintheta * cosphi, sintheta * sinphi, costheta);
413 G4double Ekin =
414 std::sqrt(daughtermomentum[0] * daughtermomentum[0] + daughtermass[0] * daughtermass[0])
415 - daughtermass[0];
416 auto daughterparticle =
417 new G4DynamicParticle(G4MT_daughters[0], direction0, Ekin, daughtermass[0]);
418 products->PushProducts(daughterparticle);
419
420 costhetan = (daughtermomentum[1] * daughtermomentum[1] - daughtermomentum[2] * daughtermomentum[2]
421 - daughtermomentum[0] * daughtermomentum[0])
422 / (2.0 * daughtermomentum[2] * daughtermomentum[0]);
423 sinthetan = std::sqrt((1.0 - costhetan) * (1.0 + costhetan));
424 phin = twopi * G4UniformRand() * rad;
425 sinphin = std::sin(phin);
426 cosphin = std::cos(phin);
427 G4ThreeVector direction2;
428 direction2.setX(sinthetan * cosphin * costheta * cosphi - sinthetan * sinphin * sinphi
429 + costhetan * sintheta * cosphi);
430 direction2.setY(sinthetan * cosphin * costheta * sinphi + sinthetan * sinphin * cosphi
431 + costhetan * sintheta * sinphi);
432 direction2.setZ(-sinthetan * cosphin * sintheta + costhetan * costheta);
433 G4ThreeVector pmom = daughtermomentum[2] * direction2 / direction2.mag();
434 Ekin = std::sqrt(pmom.mag2() + daughtermass[2] * daughtermass[2]) - daughtermass[2];
435 daughterparticle =
436 new G4DynamicParticle(G4MT_daughters[2], pmom / pmom.mag(), Ekin, daughtermass[2]);
437 products->PushProducts(daughterparticle);
438
439 pmom = (direction0 * daughtermomentum[0] + direction2 * (daughtermomentum[2] / direction2.mag()))
440 * (-1.0);
441 Ekin = std::sqrt(pmom.mag2() + daughtermass[1] * daughtermass[1]) - daughtermass[1];
442 daughterparticle =
443 new G4DynamicParticle(G4MT_daughters[1], pmom / pmom.mag(), Ekin, daughtermass[1]);
444 products->PushProducts(daughterparticle);
445
446 #ifdef G4VERBOSE
447 if (GetVerboseLevel() > 1) {
448 G4cout << "G4PhaseSpaceDecayChannel::ThreeBodyDecayIt -";
449 G4cout << " create decay products in rest frame " << G4endl;
450 products->DumpInfo();
451 }
452 #endif
453 return products;
454 }
455
456 G4DecayProducts* G4PhaseSpaceDecayChannel::ManyBodyDecayIt()
457 {
458 // Algorithm of this code originally written in FORTRAN by M.Asai
459 // **************************************************************
460 // NBODY - N-body phase space Monte-Carlo generator
461 // Makoto Asai - Hiroshima Institute of Technology
462 // Revised release : 19/Apr/1995
463
464 G4int index, index2;
465
466 #ifdef G4VERBOSE
467 if (GetVerboseLevel() > 1) G4cout << "G4PhaseSpaceDecayChannel::ManyBodyDecayIt()" << G4endl;
468 #endif
469
470 // parent mass
471 G4double parentmass = current_parent_mass.Get();
472
473 // parent particle
474 G4ThreeVector dummy;
475 auto parentparticle = new G4DynamicParticle(G4MT_parent, dummy, 0.0, parentmass);
476
477 // create G4Decayproducts
478 auto products = new G4DecayProducts(*parentparticle);
479 delete parentparticle;
480
481 // daughters'mass
482 auto daughtermass = new G4double[numberOfDaughters];
483
484 G4double sumofdaughtermass = 0.0;
485 for (index = 0; index < numberOfDaughters; ++index) {
486 if (!useGivenDaughterMass) {
487 daughtermass[index] = G4MT_daughters_mass[index];
488 }
489 else {
490 daughtermass[index] = givenDaughterMasses[index];
491 }
492 sumofdaughtermass += daughtermass[index];
493 }
494 if (sumofdaughtermass > parentmass) {
495 #ifdef G4VERBOSE
496 if (GetVerboseLevel() > 0) {
497 G4cout << "G4PhaseSpaceDecayChannel::ManyBodyDecayIt()" << G4endl
498 << "Sum of daughter mass is larger than parent mass!" << G4endl;
499 G4cout << "Parent :" << G4MT_parent->GetParticleName() << " "
500 << current_parent_mass.Get() / GeV << G4endl;
501 G4cout << "Daughter 1 :" << G4MT_daughters[0]->GetParticleName() << " "
502 << daughtermass[0] / GeV << G4endl;
503 G4cout << "Daughter 2:" << G4MT_daughters[1]->GetParticleName() << " "
504 << daughtermass[1] / GeV << G4endl;
505 G4cout << "Daughter 3:" << G4MT_daughters[2]->GetParticleName() << " "
506 << daughtermass[2] / GeV << G4endl;
507 G4cout << "Daughter 4:" << G4MT_daughters[3]->GetParticleName() << " "
508 << daughtermass[3] / GeV << G4endl;
509 G4cout << "Daughter 5:" << G4MT_daughters[4]->GetParticleName() << " "
510 << daughtermass[4] / GeV << G4endl;
511 }
512 #endif
513 G4Exception("G4PhaseSpaceDecayChannel::ManyBodyDecayIt()", "PART112", JustWarning,
514 "Can not create decay products: sum of daughter mass \
515 is larger than parent mass!");
516 delete[] daughtermass;
517 return products;
518 }
519
520 // Calculate daughter momentum
521 auto daughtermomentum = new G4double[numberOfDaughters];
522 G4ThreeVector direction;
523 G4DynamicParticle** daughterparticle;
524 auto sm = new G4double[numberOfDaughters];
525 G4double tmas;
526 G4double weight = 1.0;
527 G4int numberOfTry = 0;
528 const std::size_t MAX_LOOP = 10000;
529
530 for (std::size_t loop_counter = 0; loop_counter < MAX_LOOP; ++loop_counter) {
531 // Generate rundom number in descending order
532 G4double temp;
533 auto rd = new G4double[numberOfDaughters];
534 rd[0] = 1.0;
535 for (index = 1; index < numberOfDaughters - 1; ++index) {
536 rd[index] = G4UniformRand();
537 }
538 rd[numberOfDaughters - 1] = 0.0;
539 for (index = 1; index < numberOfDaughters - 1; ++index) {
540 for (index2 = index + 1; index2 < numberOfDaughters; ++index2) {
541 if (rd[index] < rd[index2]) {
542 temp = rd[index];
543 rd[index] = rd[index2];
544 rd[index2] = temp;
545 }
546 }
547 }
548
549 // calculate virtual mass
550 tmas = parentmass - sumofdaughtermass;
551 temp = sumofdaughtermass;
552 for (index = 0; index < numberOfDaughters; ++index) {
553 sm[index] = rd[index] * tmas + temp;
554 temp -= daughtermass[index];
555 if (GetVerboseLevel() > 1) {
556 G4cout << index << " rundom number:" << rd[index];
557 G4cout << " virtual mass:" << sm[index] / GeV << "[GeV/c/c]" << G4endl;
558 }
559 }
560 delete[] rd;
561
562 // Calculate daughter momentum
563 weight = 1.0;
564 G4bool smOK = true;
565 for (index = 0; index < numberOfDaughters - 1 && smOK; ++index) {
566 smOK = (sm[index] - daughtermass[index] - sm[index + 1] >= 0.);
567 }
568 if (!smOK) continue;
569
570 index = numberOfDaughters - 1;
571 daughtermomentum[index] = Pmx(sm[index - 1], daughtermass[index - 1], sm[index]);
572 #ifdef G4VERBOSE
573 if (GetVerboseLevel() > 1) {
574 G4cout << " daughter " << index << ":" << *daughters_name[index];
575 G4cout << " momentum:" << daughtermomentum[index] / GeV << "[GeV/c]" << G4endl;
576 }
577 #endif
578 for (index = numberOfDaughters - 2; index >= 0; --index) {
579 // calculate
580 daughtermomentum[index] = Pmx(sm[index], daughtermass[index], sm[index + 1]);
581 if (daughtermomentum[index] < 0.0) {
582 // !!! illegal momentum !!!
583 #ifdef G4VERBOSE
584 if (GetVerboseLevel() > 0) {
585 G4cout << "G4PhaseSpaceDecayChannel::ManyBodyDecayIt()" << G4endl;
586 G4cout << " Cannot calculate daughter momentum " << G4endl;
587 G4cout << " parent:" << *parent_name;
588 G4cout << " mass:" << parentmass / GeV << "[GeV/c/c]" << G4endl;
589 G4cout << " daughter " << index << ":" << *daughters_name[index];
590 G4cout << " mass:" << daughtermass[index] / GeV << "[GeV/c/c]";
591 G4cout << " mass:" << daughtermomentum[index] / GeV << "[GeV/c]" << G4endl;
592 if (useGivenDaughterMass) {
593 G4cout << "Daughter Mass is given." << G4endl;
594 }
595 }
596 #endif
597 delete[] sm;
598 delete[] daughtermass;
599 delete[] daughtermomentum;
600 delete products;
601 G4Exception("G4PhaseSpaceDecayChannel::ManyBodyDecayIt()", "PART112", JustWarning,
602 "Can not create decay products: sum of daughter mass \
603 is larger than parent mass");
604 return nullptr; // Error detection
605 }
606
607 // calculate weight of this events
608 weight *= daughtermomentum[index] / sm[index];
609 #ifdef G4VERBOSE
610 if (GetVerboseLevel() > 1) {
611 G4cout << " daughter " << index << ":" << *daughters_name[index];
612 G4cout << " momentum:" << daughtermomentum[index] / GeV << "[GeV/c]" << G4endl;
613 }
614 #endif
615 }
616
617 #ifdef G4VERBOSE
618 if (GetVerboseLevel() > 1) {
619 G4cout << " weight: " << weight << G4endl;
620 }
621 #endif
622
623 // exit if number of Try exceeds 100
624 if (++numberOfTry > 100) {
625 #ifdef G4VERBOSE
626 if (GetVerboseLevel() > 0) {
627 G4cout << "G4PhaseSpaceDecayChannel::ManyBodyDecayIt()" << G4endl;
628 G4cout << "Cannot determine Decay Kinematics " << G4endl;
629 G4cout << "parent : " << *parent_name << G4endl;
630 G4cout << "daughters : ";
631 for (index = 0; index < numberOfDaughters; ++index) {
632 G4cout << *daughters_name[index] << " , ";
633 }
634 G4cout << G4endl;
635 }
636 #endif
637 G4Exception("G4PhaseSpaceDecayChannel::ManyBodyDecayIt()", "PART113", JustWarning,
638 "Cannot decay: Decay Kinematics cannot be calculated");
639
640 delete[] sm;
641 delete[] daughtermass;
642 delete[] daughtermomentum;
643 delete products;
644 return nullptr; // Error detection
645 }
646 if (weight < G4UniformRand()) break;
647 }
648
649 #ifdef G4VERBOSE
650 if (GetVerboseLevel() > 1) {
651 G4cout << "Start calculation of daughters momentum vector " << G4endl;
652 }
653 #endif
654
655 G4double costheta, sintheta, phi;
656 G4double beta;
657 daughterparticle = new G4DynamicParticle*[numberOfDaughters];
658
659 index = numberOfDaughters - 2;
660 costheta = 2. * G4UniformRand() - 1.0;
661 sintheta = std::sqrt((1.0 - costheta) * (1.0 + costheta));
662 phi = twopi * G4UniformRand() * rad;
663 direction.setZ(costheta);
664 direction.setY(sintheta * std::sin(phi));
665 direction.setX(sintheta * std::cos(phi));
666 daughterparticle[index] =
667 new G4DynamicParticle(G4MT_daughters[index], direction * daughtermomentum[index]);
668 daughterparticle[index + 1] =
669 new G4DynamicParticle(G4MT_daughters[index + 1], direction * (-1.0 * daughtermomentum[index]));
670
671 for (index = numberOfDaughters - 3; index >= 0; --index) {
672 // calculate momentum direction
673 costheta = 2. * G4UniformRand() - 1.0;
674 sintheta = std::sqrt((1.0 - costheta) * (1.0 + costheta));
675 phi = twopi * G4UniformRand() * rad;
676 direction.setZ(costheta);
677 direction.setY(sintheta * std::sin(phi));
678 direction.setX(sintheta * std::cos(phi));
679
680 // boost already created particles
681 beta = daughtermomentum[index];
682 beta /=
683 std::sqrt(daughtermomentum[index] * daughtermomentum[index] + sm[index + 1] * sm[index + 1]);
684 for (index2 = index + 1; index2 < numberOfDaughters; ++index2) {
685 G4LorentzVector p4;
686 // make G4LorentzVector for secondaries
687 p4 = daughterparticle[index2]->Get4Momentum();
688
689 // boost secondaries to new frame
690 p4.boost(direction.x() * beta, direction.y() * beta, direction.z() * beta);
691
692 // change energy/momentum
693 daughterparticle[index2]->Set4Momentum(p4);
694 }
695 // create daughter G4DynamicParticle
696 daughterparticle[index] =
697 new G4DynamicParticle(G4MT_daughters[index], direction * (-1.0 * daughtermomentum[index]));
698 }
699
700 // add daughters to G4Decayproducts
701 for (index = 0; index < numberOfDaughters; ++index) {
702 products->PushProducts(daughterparticle[index]);
703 }
704
705 #ifdef G4VERBOSE
706 if (GetVerboseLevel() > 1) {
707 G4cout << "G4PhaseSpaceDecayChannel::ManyBodyDecayIt() -";
708 G4cout << " create decay products in rest frame " << G4endl;
709 products->DumpInfo();
710 }
711 #endif
712
713 delete[] daughterparticle;
714 delete[] daughtermomentum;
715 delete[] daughtermass;
716 delete[] sm;
717
718 return products;
719 }
720
721 G4bool G4PhaseSpaceDecayChannel::SetDaughterMasses(G4double masses[])
722 {
723 for (G4int idx = 0; idx < numberOfDaughters; ++idx) {
724 givenDaughterMasses[idx] = masses[idx];
725 }
726 useGivenDaughterMass = true;
727 return useGivenDaughterMass;
728 }
729
730 G4bool G4PhaseSpaceDecayChannel::SampleDaughterMasses()
731 {
732 useGivenDaughterMass = false;
733 return useGivenDaughterMass;
734 }
735
736 G4bool G4PhaseSpaceDecayChannel::IsOKWithParentMass(G4double parentMass)
737 {
738 if (!useGivenDaughterMass) return G4VDecayChannel::IsOKWithParentMass(parentMass);
739
740 CheckAndFillParent();
741 CheckAndFillDaughters();
742
743 G4double sumOfDaughterMassMin = 0.0;
744 for (G4int index = 0; index < numberOfDaughters; ++index) {
745 sumOfDaughterMassMin += givenDaughterMasses[index];
746 }
747 return (parentMass >= sumOfDaughterMassMin);
748 }
749
750 G4double G4PhaseSpaceDecayChannel::Pmx(G4double e, G4double p1, G4double p2)
751 {
752 // calcurate momentum of daughter particles in two-body decay
753 G4double ppp = (e + p1 + p2) * (e + p1 - p2) * (e - p1 + p2) * (e - p1 - p2) / (4.0 * e * e);
754 if (ppp > 0) return std::sqrt(ppp);
755 return -1.;
756 }
757