Geant4 11.1.1
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
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