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