Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
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//
27// $Id$
28//
29//
30// ------------------------------------------------------------
31// GEANT 4 class header file
32//
33// History: first implementation, based on object model of
34// 27 July 1996 H.Kurashige
35// 20 Oct 1996 H.Kurashige
36// 30 May 1997 H.Kurashige
37// ------------------------------------------------------------
38
41#include "G4SystemOfUnits.hh"
42#include "G4DecayProducts.hh"
43#include "G4VDecayChannel.hh"
45#include "Randomize.hh"
46#include "G4LorentzVector.hh"
47#include "G4LorentzRotation.hh"
48
49
51 :G4VDecayChannel("Phase Space", Verbose)
52{
53
54}
55
57 const G4String& theParentName,
58 G4double theBR,
59 G4int theNumberOfDaughters,
60 const G4String& theDaughterName1,
61 const G4String& theDaughterName2,
62 const G4String& theDaughterName3,
63 const G4String& theDaughterName4 )
64 :G4VDecayChannel("Phase Space",
65 theParentName,theBR,
66 theNumberOfDaughters,
67 theDaughterName1,
68 theDaughterName2,
69 theDaughterName3,
70 theDaughterName4),
71 current_parent_mass(0.0)
72{
73
74}
75
77{
78}
79
81{
82#ifdef G4VERBOSE
83 if (GetVerboseLevel()>1) G4cout << "G4PhaseSpaceDecayChannel::DecayIt ";
84#endif
85
86 G4DecayProducts * products = 0;
87
88 if (parent == 0) FillParent();
89 if (daughters == 0) FillDaughters();
90
91 if (parentMass >0.0) current_parent_mass = parentMass;
92 else current_parent_mass = parent_mass;
93
94 switch (numberOfDaughters){
95 case 0:
96#ifdef G4VERBOSE
97 if (GetVerboseLevel()>0) {
98 G4cout << "G4PhaseSpaceDecayChannel::DecayIt ";
99 G4cout << " daughters not defined " <<G4endl;
100 }
101#endif
102 break;
103 case 1:
104 products = OneBodyDecayIt();
105 break;
106 case 2:
107 products = TwoBodyDecayIt();
108 break;
109 case 3:
110 products = ThreeBodyDecayIt();
111 break;
112 default:
113 products = ManyBodyDecayIt();
114 break;
115 }
116#ifdef G4VERBOSE
117 if ((products == 0) && (GetVerboseLevel()>0)) {
118 G4cout << "G4PhaseSpaceDecayChannel::DecayIt ";
119 G4cout << *parent_name << " can not decay " << G4endl;
120 DumpInfo();
121 }
122#endif
123 return products;
124}
125
126G4DecayProducts *G4PhaseSpaceDecayChannel::OneBodyDecayIt()
127{
128#ifdef G4VERBOSE
129 if (GetVerboseLevel()>1) G4cout << "G4PhaseSpaceDecayChannel::OneBodyDecayIt()"<<G4endl;
130#endif
131
132 //create parent G4DynamicParticle at rest
133 G4ThreeVector dummy;
134 G4DynamicParticle * parentparticle = new G4DynamicParticle( parent, dummy, 0.0);
135 //create G4Decayproducts
136 G4DecayProducts *products = new G4DecayProducts(*parentparticle);
137 delete parentparticle;
138
139 //create daughter G4DynamicParticle at rest
140 G4DynamicParticle * daughterparticle = new G4DynamicParticle( daughters[0], dummy, 0.0);
141 products->PushProducts(daughterparticle);
142
143#ifdef G4VERBOSE
144 if (GetVerboseLevel()>1) {
145 G4cout << "G4PhaseSpaceDecayChannel::OneBodyDecayIt ";
146 G4cout << " create decay products in rest frame " <<G4endl;
147 products->DumpInfo();
148 }
149#endif
150 return products;
151}
152
153G4DecayProducts *G4PhaseSpaceDecayChannel::TwoBodyDecayIt()
154{
155#ifdef G4VERBOSE
156 if (GetVerboseLevel()>1) G4cout << "G4PhaseSpaceDecayChannel::TwoBodyDecayIt()"<<G4endl;
157#endif
158 // parent mass
159 G4double parentmass = current_parent_mass;
160
161 //daughters'mass
162 G4double daughtermass[2];
163 G4double daughtermomentum;
164 daughtermass[0] = daughters_mass[0];
165 daughtermass[1] = daughters_mass[1];
166
167 //create parent G4DynamicParticle at rest
168 G4ThreeVector dummy;
169 G4DynamicParticle * parentparticle = new G4DynamicParticle( parent, dummy, 0.0);
170 //create G4Decayproducts
171 G4DecayProducts *products = new G4DecayProducts(*parentparticle);
172 delete parentparticle;
173
174 //calculate daughter momentum
175 daughtermomentum = Pmx(parentmass,daughtermass[0],daughtermass[1]);
176 if (daughtermomentum <0.0) {
177#ifdef G4VERBOSE
178 if (GetVerboseLevel()>0) {
179 G4cerr << "G4PhaseSpaceDecayChannel::TwoBodyDecayIt "
180 << "sum of daughter mass is larger than parent mass" << G4endl;
181 G4cerr << "parent :" << parent->GetParticleName() << " " << current_parent_mass/GeV << G4endl;
182 G4cerr << "daughter 1 :" << daughters[0]->GetParticleName() << " " << daughtermass[0]/GeV << G4endl;
183 G4cerr << "daughter 2:" << daughters[1]->GetParticleName() << " " << daughtermass[1]/GeV << G4endl;
184 }
185#endif
186 G4Exception("G4PhaseSpaceDecayChannel::TwoBodyDecayIt",
187 "PART112", JustWarning,
188 "Can not create decay products: sum of daughter mass is larger than parent mass");
189 return products;
190 }
191
192 G4double costheta = 2.*G4UniformRand()-1.0;
193 G4double sintheta = std::sqrt((1.0 - costheta)*(1.0 + costheta));
194 G4double phi = twopi*G4UniformRand()*rad;
195 G4ThreeVector direction(sintheta*std::cos(phi),sintheta*std::sin(phi),costheta);
196
197 //create daughter G4DynamicParticle
198 G4DynamicParticle * daughterparticle = new G4DynamicParticle( daughters[0], direction*daughtermomentum);
199 products->PushProducts(daughterparticle);
200 daughterparticle = new G4DynamicParticle( daughters[1], direction*(-1.0*daughtermomentum));
201 products->PushProducts(daughterparticle);
202
203#ifdef G4VERBOSE
204 if (GetVerboseLevel()>1) {
205 G4cout << "G4PhaseSpaceDecayChannel::TwoBodyDecayIt ";
206 G4cout << " create decay products in rest frame " <<G4endl;
207 products->DumpInfo();
208 }
209#endif
210 return products;
211}
212
213G4DecayProducts *G4PhaseSpaceDecayChannel::ThreeBodyDecayIt()
214// algorism of this code is originally written in GDECA3 of GEANT3
215{
216#ifdef G4VERBOSE
217 if (GetVerboseLevel()>1) G4cout << "G4PhaseSpaceDecayChannel::ThreeBodyDecayIt()"<<G4endl;
218#endif
219 // parent mass
220 G4double parentmass = current_parent_mass;
221 //daughters'mass
222 G4double daughtermass[3];
223 G4double sumofdaughtermass = 0.0;
224 for (G4int index=0; index<3; index++){
225 daughtermass[index] = daughters_mass[index];
226 sumofdaughtermass += daughtermass[index];
227 }
228
229 //create parent G4DynamicParticle at rest
230 G4ThreeVector dummy;
231 G4DynamicParticle * parentparticle = new G4DynamicParticle( parent, dummy, 0.0);
232
233
234 //create G4Decayproducts
235 G4DecayProducts *products = new G4DecayProducts(*parentparticle);
236 delete parentparticle;
237
238 if (sumofdaughtermass >parentmass) {
239#ifdef G4VERBOSE
240 if (GetVerboseLevel()>0) {
241 G4cerr << "G4PhaseSpaceDecayChannel::ThreeBodyDecayIt "
242 << "sum of daughter mass is larger than parent mass" << G4endl;
243 G4cerr << "parent :" << parent->GetParticleName() << " " << current_parent_mass/GeV << G4endl;
244 G4cerr << "daughter 1 :" << daughters[0]->GetParticleName() << " " << daughtermass[0]/GeV << G4endl;
245 G4cerr << "daughter 2:" << daughters[1]->GetParticleName() << " " << daughtermass[1]/GeV << G4endl;
246 G4cerr << "daughter 3:" << daughters[2]->GetParticleName() << " " << daughtermass[2]/GeV << G4endl;
247 }
248#endif
249 G4Exception("G4PhaseSpaceDecayChannel::ThreeBodyDecayIt",
250 "PART112", JustWarning,
251 "Can not create decay products: sum of daughter mass is larger than parent mass");
252 return products;
253 }
254
255
256
257 //calculate daughter momentum
258 // Generate two
259 G4double rd1, rd2, rd;
260 G4double daughtermomentum[3];
261 G4double momentummax=0.0, momentumsum = 0.0;
262 G4double energy;
263
264 do {
265 rd1 = G4UniformRand();
266 rd2 = G4UniformRand();
267 if (rd2 > rd1) {
268 rd = rd1;
269 rd1 = rd2;
270 rd2 = rd;
271 }
272 momentummax = 0.0;
273 momentumsum = 0.0;
274 // daughter 0
275 energy = rd2*(parentmass - sumofdaughtermass);
276 daughtermomentum[0] = std::sqrt(energy*energy + 2.0*energy* daughtermass[0]);
277 if ( daughtermomentum[0] >momentummax )momentummax = daughtermomentum[0];
278 momentumsum += daughtermomentum[0];
279 // daughter 1
280 energy = (1.-rd1)*(parentmass - sumofdaughtermass);
281 daughtermomentum[1] = std::sqrt(energy*energy + 2.0*energy* daughtermass[1]);
282 if ( daughtermomentum[1] >momentummax )momentummax = daughtermomentum[1];
283 momentumsum += daughtermomentum[1];
284 // daughter 2
285 energy = (rd1-rd2)*(parentmass - sumofdaughtermass);
286 daughtermomentum[2] = std::sqrt(energy*energy + 2.0*energy* daughtermass[2]);
287 if ( daughtermomentum[2] >momentummax )momentummax = daughtermomentum[2];
288 momentumsum += daughtermomentum[2];
289 } while (momentummax > momentumsum - momentummax );
290 // output message
291#ifdef G4VERBOSE
292 if (GetVerboseLevel()>1) {
293 G4cout << " daughter 0:" << daughtermomentum[0]/GeV << "[GeV/c]" <<G4endl;
294 G4cout << " daughter 1:" << daughtermomentum[1]/GeV << "[GeV/c]" <<G4endl;
295 G4cout << " daughter 2:" << daughtermomentum[2]/GeV << "[GeV/c]" <<G4endl;
296 G4cout << " momentum sum:" << momentumsum/GeV << "[GeV/c]" <<G4endl;
297 }
298#endif
299
300 //create daughter G4DynamicParticle
301 G4double costheta, sintheta, phi, sinphi, cosphi;
302 G4double costhetan, sinthetan, phin, sinphin, cosphin;
303 costheta = 2.*G4UniformRand()-1.0;
304 sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
305 phi = twopi*G4UniformRand()*rad;
306 sinphi = std::sin(phi);
307 cosphi = std::cos(phi);
308 G4ThreeVector direction0(sintheta*cosphi,sintheta*sinphi,costheta);
309 G4DynamicParticle * daughterparticle
310 = new G4DynamicParticle( daughters[0], direction0*daughtermomentum[0]);
311 products->PushProducts(daughterparticle);
312
313 costhetan = (daughtermomentum[1]*daughtermomentum[1]-daughtermomentum[2]*daughtermomentum[2]-daughtermomentum[0]*daughtermomentum[0])/(2.0*daughtermomentum[2]*daughtermomentum[0]);
314 sinthetan = std::sqrt((1.0-costhetan)*(1.0+costhetan));
315 phin = twopi*G4UniformRand()*rad;
316 sinphin = std::sin(phin);
317 cosphin = std::cos(phin);
318 G4ThreeVector direction2;
319 direction2.setX( sinthetan*cosphin*costheta*cosphi - sinthetan*sinphin*sinphi + costhetan*sintheta*cosphi);
320 direction2.setY( sinthetan*cosphin*costheta*sinphi + sinthetan*sinphin*cosphi + costhetan*sintheta*sinphi);
321 direction2.setZ( -sinthetan*cosphin*sintheta + costhetan*costheta);
322 daughterparticle = new G4DynamicParticle( daughters[2], direction2*(daughtermomentum[2]/direction2.mag()));
323 products->PushProducts(daughterparticle);
324
325 daughterparticle =
327 daughters[1],
328 (direction0*daughtermomentum[0] + direction2*(daughtermomentum[2]/direction2.mag()))*(-1.0)
329 );
330 products->PushProducts(daughterparticle);
331
332#ifdef G4VERBOSE
333 if (GetVerboseLevel()>1) {
334 G4cout << "G4PhaseSpaceDecayChannel::ThreeBodyDecayIt ";
335 G4cout << " create decay products in rest frame " <<G4endl;
336 products->DumpInfo();
337 }
338#endif
339 return products;
340}
341
342G4DecayProducts *G4PhaseSpaceDecayChannel::ManyBodyDecayIt()
343// algorithm of this code is originally written in FORTRAN by M.Asai
344//******************************************************************
345// NBODY
346// N-body phase space Monte-Carlo generator
347// Makoto Asai
348// Hiroshima Institute of Technology
349// (asai@kekvax.kek.jp)
350// Revised release : 19/Apr/1995
351//
352{
353 G4int index, index2;
354
355 //return value
356 G4DecayProducts *products;
357
358#ifdef G4VERBOSE
359 if (GetVerboseLevel()>1) G4cout << "G4PhaseSpaceDecayChannel::ManyBodyDecayIt()"<<G4endl;
360#endif
361 // parent mass
362 G4double parentmass = current_parent_mass;
363 //daughters'mass
364 G4double *daughtermass = new G4double[numberOfDaughters];
365 G4double sumofdaughtermass = 0.0;
366 for (index=0; index<numberOfDaughters; index++){
367 daughtermass[index] = daughters_mass[index];
368 sumofdaughtermass += daughtermass[index];
369 }
370
371 //Calculate daughter momentum
372 G4double *daughtermomentum = new G4double[numberOfDaughters];
373 G4ThreeVector direction;
374 G4DynamicParticle **daughterparticle;
376 G4double tmas;
377 G4double weight = 1.0;
378 G4int numberOfTry = 0;
379
380 do {
381 //Generate rundom number in descending order
382 G4double temp;
384 rd[0] = 1.0;
385 for(index =1; index < numberOfDaughters -1; index++) rd[index] = G4UniformRand();
386 rd[ numberOfDaughters -1] = 0.0;
387 for(index =1; index < numberOfDaughters -1; index++) {
388 for(index2 = index+1; index2 < numberOfDaughters; index2++) {
389 if (rd[index] < rd[index2]){
390 temp = rd[index];
391 rd[index] = rd[index2];
392 rd[index2] = temp;
393 }
394 }
395 }
396
397 //calcurate virtual mass
398 tmas = parentmass - sumofdaughtermass;
399 temp = sumofdaughtermass;
400 for(index =0; index < numberOfDaughters; index++) {
401 sm[index] = rd[index]*tmas + temp;
402 temp -= daughtermass[index];
403 if (GetVerboseLevel()>1) {
404 G4cout << index << " rundom number:" << rd[index];
405 G4cout << " virtual mass:" << sm[index]/GeV << "[GeV/c/c]" <<G4endl;
406 }
407 }
408 delete [] rd;
409
410 //Calculate daughter momentum
411 weight = 1.0;
412 index =numberOfDaughters-1;
413 daughtermomentum[index]= Pmx( sm[index-1],daughtermass[index-1], sm[index]);
414#ifdef G4VERBOSE
415 if (GetVerboseLevel()>1) {
416 G4cout << " daughter " << index << ":" << *daughters_name[index];
417 G4cout << " momentum:" << daughtermomentum[index]/GeV << "[GeV/c]" <<G4endl;
418 }
419#endif
420 for(index =numberOfDaughters-2; index>=0; index--) {
421 // calculate
422 daughtermomentum[index]= Pmx( sm[index],daughtermass[index], sm[index +1]);
423 if(daughtermomentum[index] < 0.0) {
424 // !!! illegal momentum !!!
425#ifdef G4VERBOSE
426 if (GetVerboseLevel()>0) {
427 G4cout << "G4PhaseSpaceDecayChannel::ManyBodyDecayIt ";
428 G4cout << " can not calculate daughter momentum " <<G4endl;
429 G4cout << " parent:" << *parent_name;
430 G4cout << " mass:" << parentmass/GeV << "[GeV/c/c]" <<G4endl;
431 G4cout << " daughter " << index << ":" << *daughters_name[index];
432 G4cout << " mass:" << daughtermass[index]/GeV << "[GeV/c/c]" ;
433 G4cout << " mass:" << daughtermomentum[index]/GeV << "[GeV/c]" <<G4endl;
434 }
435#endif
436 delete [] sm;
437 delete [] daughtermass;
438 delete [] daughtermomentum;
439
440 G4Exception("G4PhaseSpaceDecayChannel::ManyBodyDecayIt",
441 "PART112", JustWarning,
442 "Can not create decay products: sum of daughter mass is larger than parent mass");
443
444 return 0; // Error detection
445
446 } else {
447 // calculate weight of this events
448 weight *= daughtermomentum[index]/sm[index];
449#ifdef G4VERBOSE
450 if (GetVerboseLevel()>1) {
451 G4cout << " daughter " << index << ":" << *daughters_name[index];
452 G4cout << " momentum:" << daughtermomentum[index]/GeV << "[GeV/c]" <<G4endl;
453 }
454#endif
455 }
456 }
457
458#ifdef G4VERBOSE
459 if (GetVerboseLevel()>1) {
460 G4cout << " weight: " << weight <<G4endl;
461 }
462#endif
463
464 // exit if number of Try exceeds 100
465 if (numberOfTry++ >100) {
466#ifdef G4VERBOSE
467 if (GetVerboseLevel()>0) {
468 G4cout << "G4PhaseSpaceDecayChannel::ManyBodyDecayIt: ";
469 G4cout << " can not determine Decay Kinematics " << G4endl;
470 G4cout << "parent : " << *parent_name << G4endl;
471 G4cout << "daughters : ";
472 for(index=0; index<numberOfDaughters; index++) {
473 G4cout << *daughters_name[index] << " , ";
474 }
475 G4cout << G4endl;
476 }
477#endif
478 G4Exception("G4PhaseSpaceDecayChannel::ManyBodyDecayIt: ",
479 "PART113", JustWarning,
480 " Cannot decay : Decay Kinematics cannot be calculated ");
481
482 delete [] sm;
483 delete [] daughtermass;
484 delete [] daughtermomentum;
485 return 0; // Error detection
486 }
487 } while ( weight > G4UniformRand());
488
489#ifdef G4VERBOSE
490 if (GetVerboseLevel()>1) {
491 G4cout << "Start calulation of daughters momentum vector "<<G4endl;
492 }
493#endif
494
495 G4double costheta, sintheta, phi;
496 G4double beta;
497 daughterparticle = new G4DynamicParticle*[numberOfDaughters];
498
499 index = numberOfDaughters -2;
500 costheta = 2.*G4UniformRand()-1.0;
501 sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
502 phi = twopi*G4UniformRand()*rad;
503 direction.setZ(costheta);
504 direction.setY(sintheta*std::sin(phi));
505 direction.setX(sintheta*std::cos(phi));
506 daughterparticle[index] = new G4DynamicParticle( daughters[index], direction*daughtermomentum[index] );
507 daughterparticle[index+1] = new G4DynamicParticle( daughters[index+1], direction*(-1.0*daughtermomentum[index]) );
508
509 for (index = numberOfDaughters -3; index >= 0; index--) {
510 //calculate momentum direction
511 costheta = 2.*G4UniformRand()-1.0;
512 sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
513 phi = twopi*G4UniformRand()*rad;
514 direction.setZ(costheta);
515 direction.setY(sintheta*std::sin(phi));
516 direction.setX(sintheta*std::cos(phi));
517
518 // boost already created particles
519 beta = daughtermomentum[index];
520 beta /= std::sqrt( daughtermomentum[index]*daughtermomentum[index] + sm[index+1]*sm[index+1] );
521 for (index2 = index+1; index2<numberOfDaughters; index2++) {
523 // make G4LorentzVector for secondaries
524 p4 = daughterparticle[index2]->Get4Momentum();
525
526 // boost secondaries to new frame
527 p4.boost( direction.x()*beta, direction.y()*beta, direction.z()*beta);
528
529 // change energy/momentum
530 daughterparticle[index2]->Set4Momentum(p4);
531 }
532 //create daughter G4DynamicParticle
533 daughterparticle[index]= new G4DynamicParticle( daughters[index], direction*(-1.0*daughtermomentum[index]));
534 }
535
536 //create G4Decayproducts
537 G4DynamicParticle *parentparticle;
538 direction.setX(1.0); direction.setY(0.0); direction.setZ(0.0);
539 parentparticle = new G4DynamicParticle( parent, direction, 0.0);
540 products = new G4DecayProducts(*parentparticle);
541 delete parentparticle;
542 for (index = 0; index<numberOfDaughters; index++) {
543 products->PushProducts(daughterparticle[index]);
544 }
545
546#ifdef G4VERBOSE
547 if (GetVerboseLevel()>1) {
548 G4cout << "G4PhaseSpaceDecayChannel::ManyBodyDecayIt ";
549 G4cout << " create decay products in rest frame " << G4endl;
550 products->DumpInfo();
551 }
552#endif
553
554 delete [] daughterparticle;
555 delete [] daughtermomentum;
556 delete [] daughtermass;
557 delete [] sm;
558
559 return products;
560}
561
562
564{
565 // calcurate momentum of daughter particles in two-body decay
566 G4double ppp = (e+p1+p2)*(e+p1-p2)*(e-p1+p2)*(e-p1-p2)/(4.0*e*e);
567 if (ppp>0) return std::sqrt(ppp);
568 else return -1.;
569}
570
571
572
@ JustWarning
double G4double
Definition: G4Types.hh:64
int G4int
Definition: G4Types.hh:66
#define G4endl
Definition: G4ios.hh:52
G4DLLIMPORT std::ostream G4cerr
G4DLLIMPORT std::ostream G4cout
#define G4UniformRand()
Definition: Randomize.hh:53
double z() const
double x() const
void setY(double)
double y() const
void setZ(double)
double mag() const
void setX(double)
HepLorentzVector & boost(double, double, double)
void DumpInfo() const
G4int PushProducts(G4DynamicParticle *aParticle)
G4LorentzVector Get4Momentum() const
void Set4Momentum(const G4LorentzVector &momentum)
const G4String & GetParticleName() const
static G4double Pmx(G4double e, G4double p1, G4double p2)
virtual G4DecayProducts * DecayIt(G4double)
G4String * parent_name
G4String ** daughters_name
G4int GetVerboseLevel() const
G4double * daughters_mass
G4ParticleDefinition * parent
G4ParticleDefinition ** daughters
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41