Geant4 11.1.1
Toolkit for the simulation of the passage of particles through matter
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
G4PhaseSpaceDecayChannel.cc
Go to the documentation of this file.
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 "G4ParticleDefinition.hh"
32#include "G4PhysicalConstants.hh"
33#include "G4SystemOfUnits.hh"
34#include "G4DecayProducts.hh"
35#include "G4VDecayChannel.hh"
36#include "G4PhaseSpaceDecayChannel.hh"
37#include "Randomize.hh"
38#include "G4LorentzVector.hh"
39#include "G4LorentzRotation.hh"
40
41// --------------------------------------------------------------------
42G4PhaseSpaceDecayChannel::G4PhaseSpaceDecayChannel(G4int Verbose)
43 : G4VDecayChannel("Phase Space", Verbose)
44{
45}
46
47G4PhaseSpaceDecayChannel::
48G4PhaseSpaceDecayChannel(const G4String& theParentName,
49 G4double theBR,
50 G4int theNumberOfDaughters,
51 const G4String& theDaughterName1,
52 const G4String& theDaughterName2,
53 const G4String& theDaughterName3,
54 const G4String& theDaughterName4,
55 const G4String& theDaughterName5 )
56 : G4VDecayChannel("Phase Space", theParentName,theBR, theNumberOfDaughters,
57 theDaughterName1, theDaughterName2,
58 theDaughterName3, theDaughterName4, theDaughterName5)
59{
60}
61
62// --------------------------------------------------------------------
63G4PhaseSpaceDecayChannel::~G4PhaseSpaceDecayChannel()
64{
65}
66
67// --------------------------------------------------------------------
68G4DecayProducts* G4PhaseSpaceDecayChannel::DecayIt(G4double parentMass)
69{
70#ifdef G4VERBOSE
71 if (GetVerboseLevel()>1)
72 G4cout << "G4PhaseSpaceDecayChannel::DecayIt()" << G4endl;
73#endif
74
75 G4DecayProducts* products = nullptr;
76
77 CheckAndFillParent();
78 CheckAndFillDaughters();
79
80 if (parentMass >0.0) current_parent_mass.Put( parentMass );
81 else current_parent_mass.Put(G4MT_parent_mass);
82
83 switch (numberOfDaughters)
84 {
85 case 0:
86#ifdef G4VERBOSE
87 if (GetVerboseLevel()>0)
88 {
89 G4cout << "G4PhaseSpaceDecayChannel::DecayIt() -";
90 G4cout << " daughters not defined " << G4endl;
91 }
92#endif
93 break;
94 case 1:
95 products = OneBodyDecayIt();
96 break;
97 case 2:
98 products = TwoBodyDecayIt();
99 break;
100 case 3:
101 products = ThreeBodyDecayIt();
102 break;
103 default:
104 products = ManyBodyDecayIt();
105 break;
106 }
107#ifdef G4VERBOSE
108 if ((products == nullptr) && (GetVerboseLevel()>0))
109 {
110 G4cout << "G4PhaseSpaceDecayChannel::DecayIt() - ";
111 G4cout << *parent_name << " cannot decay " << G4endl;
112 DumpInfo();
113 }
114#endif
115 return products;
116}
117
118// --------------------------------------------------------------------
119G4DecayProducts* G4PhaseSpaceDecayChannel::OneBodyDecayIt()
120{
121#ifdef G4VERBOSE
122 if (GetVerboseLevel()>1)
123 G4cout << "G4PhaseSpaceDecayChannel::OneBodyDecayIt()" << G4endl;
124#endif
125 // parent mass
126 G4double parentmass = current_parent_mass.Get();
127
128 // create parent G4DynamicParticle at rest
129 G4ThreeVector dummy;
130 G4DynamicParticle* parentparticle
131 = new G4DynamicParticle( G4MT_parent, dummy, 0.0, parentmass);
132 // create G4Decayproducts
133 G4DecayProducts *products = new G4DecayProducts(*parentparticle);
134 delete parentparticle;
135
136 // create daughter G4DynamicParticle at rest
137 G4DynamicParticle* daughterparticle
138 = new G4DynamicParticle( G4MT_daughters[0], dummy, 0.0);
139 if (useGivenDaughterMass) daughterparticle->SetMass(givenDaughterMasses[0]);
140 products->PushProducts(daughterparticle);
141
142#ifdef G4VERBOSE
143 if (GetVerboseLevel()>1)
144 {
145 G4cout << "G4PhaseSpaceDecayChannel::OneBodyDecayIt() -";
146 G4cout << " create decay products in rest frame " << G4endl;
147 products->DumpInfo();
148 }
149#endif
150 return products;
151}
152
153// --------------------------------------------------------------------
154G4DecayProducts* G4PhaseSpaceDecayChannel::TwoBodyDecayIt()
155{
156#ifdef G4VERBOSE
157 if (GetVerboseLevel()>1)
158 G4cout << "G4PhaseSpaceDecayChannel::TwoBodyDecayIt()" << G4endl;
159#endif
160 // parent mass
161 G4double parentmass = current_parent_mass.Get();
162
163 // daughters'mass, width
164 G4double daughtermass[2], daughterwidth[2];
165 daughtermass[0] = G4MT_daughters_mass[0];
166 daughtermass[1] = G4MT_daughters_mass[1];
167 daughterwidth[0] = G4MT_daughters_width[0];
168 daughterwidth[1] = G4MT_daughters_width[1];
169
170 // create parent G4DynamicParticle at rest
171 G4ThreeVector dummy;
172 G4DynamicParticle* parentparticle
173 = new G4DynamicParticle( G4MT_parent, dummy, 0.0, parentmass);
174 // create G4Decayproducts
175 G4DecayProducts* products = new G4DecayProducts(*parentparticle);
176 delete parentparticle;
177
178 if (!useGivenDaughterMass)
179 {
180 G4bool withWidth = (daughterwidth[0]>1.0e-3*daughtermass[0])
181 || (daughterwidth[1]>1.0e-3*daughtermass[1]);
182 if (withWidth)
183 {
184 G4double sumofdaughterwidthsq = daughterwidth[0]*daughterwidth[0]
185 + daughterwidth[1]*daughterwidth[1];
186 // check parent mass and daughter mass
187 G4double maxDev = (parentmass - daughtermass[0] - daughtermass[1] )
188 / std::sqrt(sumofdaughterwidthsq);
189 if (maxDev <= -1.0*rangeMass )
190 {
191#ifdef G4VERBOSE
192 if (GetVerboseLevel()>0)
193 {
194 G4cout << "G4PhaseSpaceDecayChannel::TwoBodyDecayIt()" << G4endl
195 << "Sum of daughter mass is larger than parent mass!"
196 << G4endl;
197 G4cout << "Parent :" << G4MT_parent->GetParticleName()
198 << " " << current_parent_mass.Get()/GeV << G4endl;
199 G4cout << "Daughter 1 :" << G4MT_daughters[0]->GetParticleName()
200 << " " << daughtermass[0]/GeV << G4endl;
201 G4cout << "Daughter 2:" << G4MT_daughters[1]->GetParticleName()
202 << " " << daughtermass[1]/GeV << G4endl;
203 }
204#endif
205 G4Exception("G4PhaseSpaceDecayChannel::TwoBodyDecayIt()",
206 "PART112", JustWarning,
207 "Cannot create decay products: sum of daughter mass is \
208 larger than parent mass!");
209 return products;
210 }
211 G4double dm1=daughtermass[0];
212 if (daughterwidth[0] > 0.)
213 dm1 = DynamicalMass(daughtermass[0],daughterwidth[0], maxDev);
214 G4double dm2= daughtermass[1];
215 if (daughterwidth[1] > 0.)
216 dm2 = DynamicalMass(daughtermass[1],daughterwidth[1], maxDev);
217 while (dm1+dm2>parentmass) // Loop checking, 09.08.2015, K.Kurashige
218 {
219 dm1 = DynamicalMass(daughtermass[0],daughterwidth[0], maxDev);
220 dm2 = DynamicalMass(daughtermass[1],daughterwidth[1], maxDev);
221 }
222 daughtermass[0] = dm1;
223 daughtermass[1] = dm2;
224 }
225 }
226 else
227 {
228 // use given daughter mass;
229 daughtermass[0] = givenDaughterMasses[0];
230 daughtermass[1] = givenDaughterMasses[1];
231 }
232 if (parentmass < daughtermass[0] + daughtermass[1] )
233 {
234#ifdef G4VERBOSE
235 if (GetVerboseLevel()>0)
236 {
237 G4cout << "G4PhaseSpaceDecayChannel::TwoBodyDecayIt()" << G4endl
238 << "Sum of daughter mass is larger than parent mass!" << G4endl;
239 G4cout << "Parent :" << G4MT_parent->GetParticleName()
240 << " " << current_parent_mass.Get()/GeV << G4endl;
241 G4cout << "Daughter 1 :" << G4MT_daughters[0]->GetParticleName()
242 << " " << daughtermass[0]/GeV << G4endl;
243 G4cout << "Daughter 2:" << G4MT_daughters[1]->GetParticleName()
244 << " " << daughtermass[1]/GeV << G4endl;
245 if (useGivenDaughterMass)
246 {
247 G4cout << "Daughter Mass is given." << G4endl;
248 }
249 }
250#endif
251 G4Exception("G4PhaseSpaceDecayChannel::TwoBodyDecayIt()",
252 "PART112", JustWarning,
253 "Cannot create decay products: sum of daughter mass is \
254 larger than parent mass!");
255 return products;
256 }
257
258 // calculate daughter momentum
259 G4double daughtermomentum = Pmx(parentmass,daughtermass[0],daughtermass[1]);
260
261 G4double costheta = 2.*G4UniformRand()-1.0;
262 G4double sintheta = std::sqrt((1.0 - costheta)*(1.0 + costheta));
263 G4double phi = twopi*G4UniformRand()*rad;
264 G4ThreeVector direction(sintheta*std::cos(phi),
265 sintheta*std::sin(phi), costheta);
266
267 // create daughter G4DynamicParticle
268 G4double Ekin = std::sqrt(daughtermomentum*daughtermomentum
269 + daughtermass[0]*daughtermass[0]) - daughtermass[0];
270 G4DynamicParticle* daughterparticle
271 = new G4DynamicParticle(G4MT_daughters[0],direction,Ekin,daughtermass[0]);
272 products->PushProducts(daughterparticle);
273 Ekin = std::sqrt(daughtermomentum*daughtermomentum
274 + daughtermass[1]*daughtermass[1]) - daughtermass[1];
275 daughterparticle
276 = new G4DynamicParticle(G4MT_daughters[1], -1.0*direction,
277 Ekin, daughtermass[1]);
278 products->PushProducts(daughterparticle);
279
280#ifdef G4VERBOSE
281 if (GetVerboseLevel()>1)
282 {
283 G4cout << "G4PhaseSpaceDecayChannel::TwoBodyDecayIt() -";
284 G4cout << " Create decay products in rest frame " << G4endl;
285 products->DumpInfo();
286 }
287#endif
288 return products;
289}
290
291// --------------------------------------------------------------------
292G4DecayProducts* G4PhaseSpaceDecayChannel::ThreeBodyDecayIt()
293{
294 // Algorithm of this code originally written in GDECA3 of GEANT3
295
296#ifdef G4VERBOSE
297 if (GetVerboseLevel()>1)
298 G4cout << "G4PhaseSpaceDecayChannel::ThreeBodyDecayIt()" << G4endl;
299#endif
300 // parent mass
301 G4double parentmass = current_parent_mass.Get();
302 // daughters'mass
303 G4double daughtermass[3], daughterwidth[3];
304 G4double sumofdaughtermass = 0.0;
305 G4double sumofdaughterwidthsq = 0.0;
306 G4bool withWidth = false;
307 for (G4int index=0; index<3; ++index)
308 {
309 daughtermass[index] = G4MT_daughters_mass[index];
310 sumofdaughtermass += daughtermass[index];
311 daughterwidth[index] = G4MT_daughters_width[index];
312 sumofdaughterwidthsq += daughterwidth[index]*daughterwidth[index];
313 withWidth = withWidth ||(daughterwidth[index]>1.0e-3*daughtermass[index]);
314 }
315
316 // create parent G4DynamicParticle at rest
317 G4ThreeVector dummy;
318 G4DynamicParticle* parentparticle
319 = new G4DynamicParticle( G4MT_parent, dummy, 0.0, parentmass);
320
321 // create G4Decayproducts
322 G4DecayProducts* products = new G4DecayProducts(*parentparticle);
323 delete parentparticle;
324
325 if (!useGivenDaughterMass)
326 {
327 if (withWidth)
328 {
329 G4double maxDev = (parentmass - sumofdaughtermass )
330 / std::sqrt(sumofdaughterwidthsq);
331 if (maxDev <= -1.0*rangeMass )
332 {
333#ifdef G4VERBOSE
334 if (GetVerboseLevel()>0)
335 {
336 G4cout << "G4PhaseSpaceDecayChannel::ThreeBodyDecayIt()" << G4endl
337 << "Sum of daughter mass is larger than parent mass!"
338 << G4endl;
339 G4cout << "Parent :" << G4MT_parent->GetParticleName()
340 << " " << current_parent_mass.Get()/GeV << G4endl;
341 G4cout << "Daughter 1 :" << G4MT_daughters[0]->GetParticleName()
342 << " " << daughtermass[0]/GeV << G4endl;
343 G4cout << "Daughter 2:" << G4MT_daughters[1]->GetParticleName()
344 << " " << daughtermass[1]/GeV << G4endl;
345 G4cout << "Daughter 3:" << G4MT_daughters[2]->GetParticleName()
346 << " " << daughtermass[2]/GeV << G4endl;
347 }
348#endif
349 G4Exception("G4PhaseSpaceDecayChannel::ThreeBodyDecayIt()",
350 "PART112", JustWarning,
351 "Cannot create decay products: sum of daughter mass \
352 is larger than parent mass!");
353 return products;
354 }
355 G4double dm1=daughtermass[0];
356 if (daughterwidth[0] > 0.)
357 dm1= DynamicalMass(daughtermass[0],daughterwidth[0], maxDev);
358 G4double dm2= daughtermass[1];
359 if (daughterwidth[1] > 0.)
360 dm2= DynamicalMass(daughtermass[1],daughterwidth[1], maxDev);
361 G4double dm3= daughtermass[2];
362 if (daughterwidth[2] > 0.)
363 dm3= DynamicalMass(daughtermass[2],daughterwidth[2], maxDev);
364 while (dm1+dm2+dm3>parentmass) // Loop checking, 09.08.2015, K.Kurashige
365 {
366 dm1= DynamicalMass(daughtermass[0],daughterwidth[0], maxDev);
367 dm2= DynamicalMass(daughtermass[1],daughterwidth[1], maxDev);
368 dm3= DynamicalMass(daughtermass[2],daughterwidth[2], maxDev);
369 }
370 daughtermass[0] = dm1;
371 daughtermass[1] = dm2;
372 daughtermass[2] = dm3;
373 sumofdaughtermass = dm1+dm2+dm3;
374 }
375 }
376 else
377 {
378 // use given daughter mass;
379 daughtermass[0] = givenDaughterMasses[0];
380 daughtermass[1] = givenDaughterMasses[1];
381 daughtermass[2] = givenDaughterMasses[2];
382 sumofdaughtermass = daughtermass[0] + daughtermass[1] + daughtermass[2];
383 }
384
385 if (sumofdaughtermass >parentmass)
386 {
387#ifdef G4VERBOSE
388 if (GetVerboseLevel()>0)
389 {
390 G4cout << "G4PhaseSpaceDecayChannel::ThreeBodyDecayIt()" << G4endl
391 << "Sum of daughter mass is larger than parent mass!" << G4endl;
392 G4cout << "Parent :" << G4MT_parent->GetParticleName()
393 << " " << current_parent_mass.Get()/GeV << G4endl;
394 G4cout << "Daughter 1 :" << G4MT_daughters[0]->GetParticleName()
395 << " " << daughtermass[0]/GeV << G4endl;
396 G4cout << "Daughter 2:" << G4MT_daughters[1]->GetParticleName()
397 << " " << daughtermass[1]/GeV << G4endl;
398 G4cout << "Daughter 3:" << G4MT_daughters[2]->GetParticleName()
399 << " " << daughtermass[2]/GeV << G4endl;
400 }
401 if (useGivenDaughterMass)
402 {
403 G4cout << "Daughter Mass is given." << G4endl;
404 }
405#endif
406 G4Exception("G4PhaseSpaceDecayChannel::ThreeBodyDecayIt",
407 "PART112", JustWarning,
408 "Can not create decay products: sum of daughter mass \
409 is larger than parent mass!");
410 return products;
411 }
412
413 // calculate daughter momentum
414 // Generate two
415 G4double rd1, rd2, rd;
416 G4double daughtermomentum[3];
417 G4double momentummax=0.0, momentumsum = 0.0;
418 G4double energy;
419 const std::size_t MAX_LOOP=10000;
420
421 for (std::size_t loop_counter=0; loop_counter<MAX_LOOP; ++loop_counter)
422 {
423 rd1 = G4UniformRand();
424 rd2 = G4UniformRand();
425 if (rd2 > rd1)
426 {
427 rd = rd1;
428 rd1 = rd2;
429 rd2 = rd;
430 }
431 momentummax = 0.0;
432 momentumsum = 0.0;
433 // daughter 0
434 energy = rd2*(parentmass - sumofdaughtermass);
435 daughtermomentum[0] = std::sqrt(energy*energy
436 + 2.0*energy* daughtermass[0]);
437 if ( daughtermomentum[0] > momentummax )
438 momentummax = daughtermomentum[0];
439 momentumsum += daughtermomentum[0];
440 // daughter 1
441 energy = (1.-rd1)*(parentmass - sumofdaughtermass);
442 daughtermomentum[1] = std::sqrt(energy*energy
443 + 2.0*energy*daughtermass[1]);
444 if ( daughtermomentum[1] >momentummax )momentummax = daughtermomentum[1];
445 momentumsum += daughtermomentum[1];
446 // daughter 2
447 energy = (rd1-rd2)*(parentmass - sumofdaughtermass);
448 daughtermomentum[2] = std::sqrt(energy*energy
449 + 2.0*energy* daughtermass[2]);
450 if ( daughtermomentum[2] > momentummax )
451 momentummax = daughtermomentum[2];
452 momentumsum += daughtermomentum[2];
453 if (momentummax <= momentumsum - momentummax ) break;
454 }
455
456 // output message
457#ifdef G4VERBOSE
458 if (GetVerboseLevel()>1)
459 {
460 G4cout << " daughter 0:" << daughtermomentum[0]/GeV << "[GeV/c]"
461 << G4endl;
462 G4cout << " daughter 1:" << daughtermomentum[1]/GeV << "[GeV/c]"
463 << G4endl;
464 G4cout << " daughter 2:" << daughtermomentum[2]/GeV << "[GeV/c]"
465 << G4endl;
466 G4cout << " momentum sum:" << momentumsum/GeV << "[GeV/c]"
467 << G4endl;
468 }
469#endif
470
471 // create daughter G4DynamicParticle
472 G4double costheta, sintheta, phi, sinphi, cosphi;
473 G4double costhetan, sinthetan, phin, sinphin, cosphin;
474 costheta = 2.*G4UniformRand()-1.0;
475 sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
476 phi = twopi*G4UniformRand()*rad;
477 sinphi = std::sin(phi);
478 cosphi = std::cos(phi);
479
480 G4ThreeVector direction0(sintheta*cosphi,sintheta*sinphi,costheta);
481 G4double Ekin = std::sqrt(daughtermomentum[0]*daughtermomentum[0]
482 + daughtermass[0]*daughtermass[0]) - daughtermass[0];
483 G4DynamicParticle* daughterparticle
484 = new G4DynamicParticle(G4MT_daughters[0],direction0,Ekin,daughtermass[0]);
485 products->PushProducts(daughterparticle);
486
487 costhetan = (daughtermomentum[1]*daughtermomentum[1]
488 - daughtermomentum[2]*daughtermomentum[2]
489 - daughtermomentum[0]*daughtermomentum[0])
490 / (2.0*daughtermomentum[2]*daughtermomentum[0]);
491 sinthetan = std::sqrt((1.0-costhetan)*(1.0+costhetan));
492 phin = twopi*G4UniformRand()*rad;
493 sinphin = std::sin(phin);
494 cosphin = std::cos(phin);
495 G4ThreeVector direction2;
496 direction2.setX( sinthetan*cosphin*costheta*cosphi
497 - sinthetan*sinphin*sinphi + costhetan*sintheta*cosphi);
498 direction2.setY( sinthetan*cosphin*costheta*sinphi
499 + sinthetan*sinphin*cosphi + costhetan*sintheta*sinphi);
500 direction2.setZ( -sinthetan*cosphin*sintheta + costhetan*costheta);
501 G4ThreeVector pmom = daughtermomentum[2]*direction2/direction2.mag();
502 Ekin = std::sqrt(pmom.mag2() + daughtermass[2]*daughtermass[2])
503 - daughtermass[2];
504 daughterparticle = new G4DynamicParticle( G4MT_daughters[2],
505 pmom/pmom.mag(),
506 Ekin, daughtermass[2] );
507 products->PushProducts(daughterparticle);
508
509 pmom = (direction0*daughtermomentum[0]
510 + direction2*(daughtermomentum[2]/direction2.mag()))*(-1.0);
511 Ekin = std::sqrt(pmom.mag2() + daughtermass[1]*daughtermass[1])
512 - daughtermass[1];
513 daughterparticle = new G4DynamicParticle( G4MT_daughters[1],
514 pmom/pmom.mag(),
515 Ekin, daughtermass[1] );
516 products->PushProducts(daughterparticle);
517
518#ifdef G4VERBOSE
519 if (GetVerboseLevel()>1)
520 {
521 G4cout << "G4PhaseSpaceDecayChannel::ThreeBodyDecayIt -";
522 G4cout << " create decay products in rest frame " << G4endl;
523 products->DumpInfo();
524 }
525#endif
526 return products;
527}
528
529// --------------------------------------------------------------------
530G4DecayProducts* G4PhaseSpaceDecayChannel::ManyBodyDecayIt()
531{
532 // Algorithm of this code originally written in FORTRAN by M.Asai
533 // **************************************************************
534 // NBODY - N-body phase space Monte-Carlo generator
535 // Makoto Asai - Hiroshima Institute of Technology
536 // Revised release : 19/Apr/1995
537
538 G4int index, index2;
539
540#ifdef G4VERBOSE
541 if (GetVerboseLevel()>1)
542 G4cout << "G4PhaseSpaceDecayChannel::ManyBodyDecayIt()" << G4endl;
543#endif
544
545 // parent mass
546 G4double parentmass = current_parent_mass.Get();
547
548 // parent particle
549 G4ThreeVector dummy;
550 G4DynamicParticle* parentparticle
551 = new G4DynamicParticle( G4MT_parent, dummy, 0.0, parentmass);
552
553 // create G4Decayproducts
554 G4DecayProducts* products = new G4DecayProducts(*parentparticle);
555 delete parentparticle;
556
557 // daughters'mass
558 G4double* daughtermass = new G4double[numberOfDaughters];
559
560 G4double sumofdaughtermass = 0.0;
561 for (index=0; index<numberOfDaughters; ++index)
562 {
563 if (!useGivenDaughterMass)
564 {
565 daughtermass[index] = G4MT_daughters_mass[index];
566 }
567 else
568 {
569 daughtermass[index] = givenDaughterMasses[index];
570 }
571 sumofdaughtermass += daughtermass[index];
572 }
573 if (sumofdaughtermass >parentmass)
574 {
575#ifdef G4VERBOSE
576 if (GetVerboseLevel()>0)
577 {
578 G4cout << "G4PhaseSpaceDecayChannel::ManyBodyDecayIt()" << G4endl
579 << "Sum of daughter mass is larger than parent mass!" << G4endl;
580 G4cout << "Parent :" << G4MT_parent->GetParticleName()
581 << " " << current_parent_mass.Get()/GeV << G4endl;
582 G4cout << "Daughter 1 :" << G4MT_daughters[0]->GetParticleName()
583 << " " << daughtermass[0]/GeV << G4endl;
584 G4cout << "Daughter 2:" << G4MT_daughters[1]->GetParticleName()
585 << " " << daughtermass[1]/GeV << G4endl;
586 G4cout << "Daughter 3:" << G4MT_daughters[2]->GetParticleName()
587 << " " << daughtermass[2]/GeV << G4endl;
588 G4cout << "Daughter 4:" << G4MT_daughters[3]->GetParticleName()
589 << " " << daughtermass[3]/GeV << G4endl;
590 G4cout << "Daughter 5:" << G4MT_daughters[4]->GetParticleName()
591 << " " << daughtermass[4]/GeV << G4endl;
592 }
593#endif
594 G4Exception("G4PhaseSpaceDecayChannel::ManyBodyDecayIt()",
595 "PART112", JustWarning,
596 "Can not create decay products: sum of daughter mass \
597 is larger than parent mass!");
598 delete [] daughtermass;
599 return products;
600 }
601
602 // Calculate daughter momentum
603 G4double* daughtermomentum = new G4double[numberOfDaughters];
604 G4ThreeVector direction;
605 G4DynamicParticle** daughterparticle;
606 G4double* sm = new G4double[numberOfDaughters];
607 G4double tmas;
608 G4double weight = 1.0;
609 G4int numberOfTry = 0;
610 const std::size_t MAX_LOOP=10000;
611
612 for (std::size_t loop_counter=0; loop_counter<MAX_LOOP; ++loop_counter)
613 {
614 // Generate rundom number in descending order
615 G4double temp;
616 G4double* rd = new G4double[numberOfDaughters];
617 rd[0] = 1.0;
618 for( index=1; index<numberOfDaughters-1; ++index )
619 {
620 rd[index] = G4UniformRand();
621 }
622 rd[ numberOfDaughters -1] = 0.0;
623 for( index=1; index<numberOfDaughters-1; ++index)
624 {
625 for( index2=index+1; index2<numberOfDaughters; ++index2)
626 {
627 if (rd[index] < rd[index2])
628 {
629 temp = rd[index];
630 rd[index] = rd[index2];
631 rd[index2] = temp;
632 }
633 }
634 }
635
636 // calculate virtual mass
637 tmas = parentmass - sumofdaughtermass;
638 temp = sumofdaughtermass;
639 for( index=0; index<numberOfDaughters; ++index )
640 {
641 sm[index] = rd[index]*tmas + temp;
642 temp -= daughtermass[index];
643 if (GetVerboseLevel()>1)
644 {
645 G4cout << index << " rundom number:" << rd[index];
646 G4cout << " virtual mass:" << sm[index]/GeV << "[GeV/c/c]" << G4endl;
647 }
648 }
649 delete [] rd;
650
651 // Calculate daughter momentum
652 weight = 1.0;
653 G4bool smOK=true;
654 for( index=0; index<numberOfDaughters-1 && smOK; ++index )
655 {
656 smOK = (sm[index]-daughtermass[index]- sm[index+1] >=0.);
657 }
658 if (!smOK) continue;
659
660 index = numberOfDaughters-1;
661 daughtermomentum[index]= Pmx(sm[index-1],daughtermass[index-1],sm[index]);
662#ifdef G4VERBOSE
663 if (GetVerboseLevel()>1)
664 {
665 G4cout << " daughter " << index << ":" << *daughters_name[index];
666 G4cout << " momentum:" << daughtermomentum[index]/GeV << "[GeV/c]"
667 << G4endl;
668 }
669#endif
670 for( index=numberOfDaughters-2; index>=0; --index)
671 {
672 // calculate
673 daughtermomentum[index]= Pmx( sm[index],daughtermass[index],sm[index +1]);
674 if(daughtermomentum[index] < 0.0)
675 {
676 // !!! illegal momentum !!!
677#ifdef G4VERBOSE
678 if (GetVerboseLevel()>0)
679 {
680 G4cout << "G4PhaseSpaceDecayChannel::ManyBodyDecayIt()" << G4endl;
681 G4cout << " Cannot calculate daughter momentum " << G4endl;
682 G4cout << " parent:" << *parent_name;
683 G4cout << " mass:" << parentmass/GeV << "[GeV/c/c]" << G4endl;
684 G4cout << " daughter " << index << ":" << *daughters_name[index];
685 G4cout << " mass:" << daughtermass[index]/GeV << "[GeV/c/c]";
686 G4cout << " mass:" << daughtermomentum[index]/GeV << "[GeV/c]"
687 << G4endl;
688 if (useGivenDaughterMass)
689 {
690 G4cout << "Daughter Mass is given." << G4endl;
691 }
692 }
693#endif
694 delete [] sm;
695 delete [] daughtermass;
696 delete [] daughtermomentum;
697 delete products;
698 G4Exception("G4PhaseSpaceDecayChannel::ManyBodyDecayIt()",
699 "PART112", JustWarning,
700 "Can not create decay products: sum of daughter mass \
701 is larger than parent mass");
702 return nullptr; // Error detection
703
704 }
705 else
706 {
707 // calculate weight of this events
708 weight *= daughtermomentum[index]/sm[index];
709#ifdef G4VERBOSE
710 if (GetVerboseLevel()>1)
711 {
712 G4cout << " daughter " << index << ":" << *daughters_name[index];
713 G4cout << " momentum:" << daughtermomentum[index]/GeV << "[GeV/c]"
714 << G4endl;
715 }
716#endif
717 }
718 }
719
720#ifdef G4VERBOSE
721 if (GetVerboseLevel()>1)
722 {
723 G4cout << " weight: " << weight << G4endl;
724 }
725#endif
726
727 // exit if number of Try exceeds 100
728 if (++numberOfTry > 100)
729 {
730#ifdef G4VERBOSE
731 if (GetVerboseLevel()>0)
732 {
733 G4cout << "G4PhaseSpaceDecayChannel::ManyBodyDecayIt()" << G4endl;
734 G4cout << "Cannot determine Decay Kinematics " << G4endl;
735 G4cout << "parent : " << *parent_name << G4endl;
736 G4cout << "daughters : ";
737 for(index=0; index<numberOfDaughters; ++index)
738 {
739 G4cout << *daughters_name[index] << " , ";
740 }
741 G4cout << G4endl;
742 }
743#endif
744 G4Exception("G4PhaseSpaceDecayChannel::ManyBodyDecayIt()",
745 "PART113", JustWarning,
746 "Cannot decay: Decay Kinematics cannot be calculated");
747
748 delete [] sm;
749 delete [] daughtermass;
750 delete [] daughtermomentum;
751 delete products;
752 return nullptr; // Error detection
753 }
754 if ( weight < G4UniformRand()) break;
755 }
756
757#ifdef G4VERBOSE
758 if (GetVerboseLevel()>1)
759 {
760 G4cout << "Start calculation of daughters momentum vector " << G4endl;
761 }
762#endif
763
764 G4double costheta, sintheta, phi;
765 G4double beta;
766 daughterparticle = new G4DynamicParticle*[numberOfDaughters];
767
768 index = numberOfDaughters -2;
769 costheta = 2.*G4UniformRand()-1.0;
770 sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
771 phi = twopi*G4UniformRand()*rad;
772 direction.setZ(costheta);
773 direction.setY(sintheta*std::sin(phi));
774 direction.setX(sintheta*std::cos(phi));
775 daughterparticle[index] = new G4DynamicParticle( G4MT_daughters[index],
776 direction*daughtermomentum[index] );
777 daughterparticle[index+1] = new G4DynamicParticle( G4MT_daughters[index+1],
778 direction*(-1.0*daughtermomentum[index]) );
779
780 for (index = numberOfDaughters-3; index >= 0; --index)
781 {
782 // calculate momentum direction
783 costheta = 2.*G4UniformRand()-1.0;
784 sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
785 phi = twopi*G4UniformRand()*rad;
786 direction.setZ(costheta);
787 direction.setY(sintheta*std::sin(phi));
788 direction.setX(sintheta*std::cos(phi));
789
790 // boost already created particles
791 beta = daughtermomentum[index];
792 beta /= std::sqrt( daughtermomentum[index]*daughtermomentum[index]
793 + sm[index+1]*sm[index+1] );
794 for (index2 = index+1; index2<numberOfDaughters; ++index2)
795 {
796 G4LorentzVector p4;
797 // make G4LorentzVector for secondaries
798 p4 = daughterparticle[index2]->Get4Momentum();
799
800 // boost secondaries to new frame
801 p4.boost( direction.x()*beta, direction.y()*beta, direction.z()*beta);
802
803 // change energy/momentum
804 daughterparticle[index2]->Set4Momentum(p4);
805 }
806 // create daughter G4DynamicParticle
807 daughterparticle[index] = new G4DynamicParticle( G4MT_daughters[index],
808 direction*(-1.0*daughtermomentum[index]));
809 }
810
811 // add daughters to G4Decayproducts
812 for (index = 0; index<numberOfDaughters; ++index)
813 {
814 products->PushProducts(daughterparticle[index]);
815 }
816
817#ifdef G4VERBOSE
818 if (GetVerboseLevel()>1)
819 {
820 G4cout << "G4PhaseSpaceDecayChannel::ManyBodyDecayIt() -";
821 G4cout << " create decay products in rest frame " << G4endl;
822 products->DumpInfo();
823 }
824#endif
825
826 delete [] daughterparticle;
827 delete [] daughtermomentum;
828 delete [] daughtermass;
829 delete [] sm;
830
831 return products;
832}
833
834// --------------------------------------------------------------------
835G4bool G4PhaseSpaceDecayChannel::SetDaughterMasses( G4double masses[] )
836{
837 for (G4int idx=0; idx<numberOfDaughters; ++idx)
838 {
839 givenDaughterMasses[idx] = masses[idx];
840 }
841 useGivenDaughterMass = true;
842 return useGivenDaughterMass;
843}
844
845// --------------------------------------------------------------------
846G4bool G4PhaseSpaceDecayChannel::SampleDaughterMasses()
847{
848 useGivenDaughterMass = false;
849 return useGivenDaughterMass;
850}
851
852// --------------------------------------------------------------------
853G4bool G4PhaseSpaceDecayChannel::IsOKWithParentMass(G4double parentMass)
854{
855 if (!useGivenDaughterMass)
856 return G4VDecayChannel::IsOKWithParentMass(parentMass);
857
858 CheckAndFillParent();
859 CheckAndFillDaughters();
860
861 G4double sumOfDaughterMassMin=0.0;
862 for (G4int index=0; index<numberOfDaughters; ++index)
863 {
864 sumOfDaughterMassMin += givenDaughterMasses[index];
865 }
866 return (parentMass >= sumOfDaughterMassMin);
867}
868
869// --------------------------------------------------------------------
870G4double G4PhaseSpaceDecayChannel::Pmx(G4double e, G4double p1, G4double p2)
871{
872 // calcurate momentum of daughter particles in two-body decay
873 G4double ppp = (e+p1+p2)*(e+p1-p2)*(e-p1+p2)*(e-p1-p2)/(4.0*e*e);
874 if (ppp>0) return std::sqrt(ppp);
875 else return -1.;
876}
JustWarning
@ JustWarning
Definition: G4ExceptionSeverity.hh:65
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.cc:59
G4double
double G4double
Definition: G4Types.hh:83
G4bool
bool G4bool
Definition: G4Types.hh:86
G4int
int G4int
Definition: G4Types.hh:85
G4endl
#define G4endl
Definition: G4ios.hh:57
G4cout
G4GLOB_DLL std::ostream G4cout
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:52
CLHEP::Hep3Vector::z
double z() const
CLHEP::Hep3Vector::x
double x() const
CLHEP::Hep3Vector::setY
void setY(double)
CLHEP::Hep3Vector::mag2
double mag2() const
CLHEP::Hep3Vector::y
double y() const
CLHEP::Hep3Vector::setZ
void setZ(double)
CLHEP::Hep3Vector::mag
double mag() const
CLHEP::Hep3Vector::setX
void setX(double)
CLHEP::HepLorentzVector::boost
HepLorentzVector & boost(double, double, double)
Definition: LorentzVector.cc:55
G4Cache::Get
value_type & Get() const
Definition: G4Cache.hh:315
G4Cache::Put
void Put(const value_type &val) const
Definition: G4Cache.hh:321
G4DecayProducts::DumpInfo
void DumpInfo() const
Definition: G4DecayProducts.cc:334
G4DecayProducts::PushProducts
G4int PushProducts(G4DynamicParticle *aParticle)
Definition: G4DecayProducts.cc:151
G4DynamicParticle::SetMass
void SetMass(G4double mass)
G4DynamicParticle::Get4Momentum
G4LorentzVector Get4Momentum() const
G4DynamicParticle::Set4Momentum
void Set4Momentum(const G4LorentzVector &momentum)
Definition: G4DynamicParticle.cc:410
G4ParticleDefinition::GetParticleName
const G4String & GetParticleName() const
Definition: G4ParticleDefinition.hh:109
G4PhaseSpaceDecayChannel::Pmx
static G4double Pmx(G4double e, G4double p1, G4double p2)
Definition: G4PhaseSpaceDecayChannel.cc:870
G4PhaseSpaceDecayChannel::DecayIt
virtual G4DecayProducts * DecayIt(G4double)
Definition: G4PhaseSpaceDecayChannel.cc:68
G4PhaseSpaceDecayChannel::IsOKWithParentMass
G4bool IsOKWithParentMass(G4double parentMass)
Definition: G4PhaseSpaceDecayChannel.cc:853
G4PhaseSpaceDecayChannel::SetDaughterMasses
G4bool SetDaughterMasses(G4double masses[])
Definition: G4PhaseSpaceDecayChannel.cc:835
G4VDecayChannel::G4MT_daughters
G4ParticleDefinition ** G4MT_daughters
Definition: G4VDecayChannel.hh:173
G4VDecayChannel::parent_name
G4String * parent_name
Definition: G4VDecayChannel.hh:156
G4VDecayChannel::daughters_name
G4String ** daughters_name
Definition: G4VDecayChannel.hh:158
G4VDecayChannel::DynamicalMass
G4double DynamicalMass(G4double massPDG, G4double width, G4double maxDev=1.0) const
Definition: G4VDecayChannel.cc:586
G4VDecayChannel::GetVerboseLevel
G4int GetVerboseLevel() const
Definition: G4VDecayChannel.hh:310
G4VDecayChannel::G4MT_parent_mass
G4double G4MT_parent_mass
Definition: G4VDecayChannel.hh:174
G4VDecayChannel::G4MT_parent
G4ParticleDefinition * G4MT_parent
Definition: G4VDecayChannel.hh:172
G4VDecayChannel::IsOKWithParentMass
virtual G4bool IsOKWithParentMass(G4double parentMass)
Definition: G4VDecayChannel.cc:607
G4VDecayChannel::CheckAndFillDaughters
void CheckAndFillDaughters()
Definition: G4VDecayChannel.hh:337
G4VDecayChannel::G4MT_daughters_mass
G4double * G4MT_daughters_mass
Definition: G4VDecayChannel.hh:175
G4VDecayChannel::G4MT_daughters_width
G4double * G4MT_daughters_width
Definition: G4VDecayChannel.hh:176
G4INCL::KinematicsUtils::energy
G4double energy(const ThreeVector &p, const G4double m)
Definition: G4INCLKinematicsUtils.cc:158