Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
Public Member Functions | Static Public Member Functions | Protected Member Functions | List of all members
G4GeneralPhaseSpaceDecay Class Reference

#include <G4GeneralPhaseSpaceDecay.hh>

+ Inheritance diagram for G4GeneralPhaseSpaceDecay:

Public Member Functions

 G4GeneralPhaseSpaceDecay (G4int Verbose=1)
 
 G4GeneralPhaseSpaceDecay (const G4String &theParentName, G4double theBR, G4int theNumberOfDaughters, const G4String &theDaughterName1, const G4String &theDaughterName2="", const G4String &theDaughterName3="")
 
 G4GeneralPhaseSpaceDecay (const G4String &theParentName, G4double theParentMass, G4double theBR, G4int theNumberOfDaughters, const G4String &theDaughterName1, const G4String &theDaughterName2="", const G4String &theDaughterName3="")
 
 G4GeneralPhaseSpaceDecay (const G4String &theParentName, G4double theParentMass, G4double theBR, G4int theNumberOfDaughters, const G4String &theDaughterName1, const G4String &theDaughterName2, const G4String &theDaughterName3, const G4double *masses)
 
virtual ~G4GeneralPhaseSpaceDecay ()
 
G4double GetParentMass () const
 
void SetParentMass (const G4double aParentMass)
 
virtual G4DecayProducts * DecayIt (G4double mass=0.0)
 
- Public Member Functions inherited from G4VDecayChannel
 G4VDecayChannel (const G4String &aName, G4int Verbose=1)
 
 G4VDecayChannel (const G4String &aName, const G4String &theParentName, G4double theBR, G4int theNumberOfDaughters, const G4String &theDaughterName1, const G4String &theDaughterName2="", const G4String &theDaughterName3="", const G4String &theDaughterName4="")
 
virtual ~G4VDecayChannel ()
 
G4int operator== (const G4VDecayChannel &right) const
 
G4int operator!= (const G4VDecayChannel &right) const
 
G4int operator< (const G4VDecayChannel &right) const
 
virtual G4DecayProducts * DecayIt (G4double parentMass=-1.0)=0
 
const G4String & GetKinematicsName () const
 
G4double GetBR () const
 
G4int GetNumberOfDaughters () const
 
G4ParticleDefinition * GetParent ()
 
G4ParticleDefinition * GetDaughter (G4int anIndex)
 
G4int GetAngularMomentum ()
 
const G4String & GetParentName () const
 
const G4String & GetDaughterName (G4int anIndex) const
 
G4double GetParentMass () const
 
G4double GetDaughterMass (G4int anIndex) const
 
void SetParent (const G4ParticleDefinition *particle_type)
 
void SetParent (const G4String &particle_name)
 
void SetBR (G4double value)
 
void SetNumberOfDaughters (G4int value)
 
void SetDaughter (G4int anIndex, const G4ParticleDefinition *particle_type)
 
void SetDaughter (G4int anIndex, const G4String &particle_name)
 
void SetVerboseLevel (G4int value)
 
G4int GetVerboseLevel () const
 
void DumpInfo ()
 

Static Public Member Functions

static G4double Pmx (G4double e, G4double p1, G4double p2)
 

Protected Member Functions

G4DecayProducts * OneBodyDecayIt ()
 
G4DecayProducts * TwoBodyDecayIt ()
 
G4DecayProducts * ThreeBodyDecayIt ()
 
G4DecayProducts * ManyBodyDecayIt ()
 
- Protected Member Functions inherited from G4VDecayChannel
 G4VDecayChannel ()
 
 G4VDecayChannel (const G4VDecayChannel &)
 
G4VDecayChannel & operator= (const G4VDecayChannel &)
 
void ClearDaughtersName ()
 
void FillDaughters ()
 
void FillParent ()
 

Additional Inherited Members

- Protected Attributes inherited from G4VDecayChannel
G4String kinematics_name
 
G4double rbranch
 
G4int numberOfDaughters
 
G4String * parent_name
 
G4String ** daughters_name
 
G4ParticleTable * particletable
 
G4ParticleDefinition * parent
 
G4ParticleDefinition ** daughters
 
G4double parent_mass
 
G4double * daughters_mass
 
G4int verboseLevel
 
- Static Protected Attributes inherited from G4VDecayChannel
static const G4String noName = " "
 

Detailed Description

Definition at line 45 of file G4GeneralPhaseSpaceDecay.hh.

Constructor & Destructor Documentation

◆ G4GeneralPhaseSpaceDecay() [1/4]

G4GeneralPhaseSpaceDecay::G4GeneralPhaseSpaceDecay ( G4int  Verbose = 1)

Definition at line 50 of file G4GeneralPhaseSpaceDecay.cc.

50 :
51 G4VDecayChannel("Phase Space", Verbose),
52 parentmass(0.), theDaughterMasses(0)
53{
54 if (GetVerboseLevel()>1) G4cout << "G4GeneralPhaseSpaceDecay:: constructor " << G4endl;
55}
G4endl
#define G4endl
Definition: G4ios.hh:52
G4cout
G4DLLIMPORT std::ostream G4cout
G4VDecayChannel::GetVerboseLevel
G4int GetVerboseLevel() const
Definition: G4VDecayChannel.hh:276

◆ G4GeneralPhaseSpaceDecay() [2/4]

G4GeneralPhaseSpaceDecay::G4GeneralPhaseSpaceDecay ( const G4String &  theParentName,
G4double  theBR,
G4int  theNumberOfDaughters,
const G4String &  theDaughterName1,
const G4String &  theDaughterName2 = "",
const G4String &  theDaughterName3 = "" 
)

Definition at line 57 of file G4GeneralPhaseSpaceDecay.cc.

62 :
63 G4VDecayChannel("Phase Space",
64 theParentName,theBR,
65 theNumberOfDaughters,
66 theDaughterName1,
67 theDaughterName2,
68 theDaughterName3),
69 theDaughterMasses(0)
70{
71 if (GetVerboseLevel()>1) G4cout << "G4GeneralPhaseSpaceDecay:: constructor " << G4endl;
72
73 // Set the parent particle (resonance) mass to the (default) PDG vale
74 if (parent != NULL)
75 {
76 parentmass = parent->GetPDGMass();
77 } else {
78 parentmass=0.;
79 }
80
81}
G4VDecayChannel::parent
G4ParticleDefinition * parent
Definition: G4VDecayChannel.hh:146

◆ G4GeneralPhaseSpaceDecay() [3/4]

G4GeneralPhaseSpaceDecay::G4GeneralPhaseSpaceDecay ( const G4String &  theParentName,
G4double  theParentMass,
G4double  theBR,
G4int  theNumberOfDaughters,
const G4String &  theDaughterName1,
const G4String &  theDaughterName2 = "",
const G4String &  theDaughterName3 = "" 
)

Definition at line 83 of file G4GeneralPhaseSpaceDecay.cc.

89 :
90 G4VDecayChannel("Phase Space",
91 theParentName,theBR,
92 theNumberOfDaughters,
93 theDaughterName1,
94 theDaughterName2,
95 theDaughterName3),
96 parentmass(theParentMass),
97 theDaughterMasses(0)
98{
99 if (GetVerboseLevel()>1) G4cout << "G4GeneralPhaseSpaceDecay:: constructor " << G4endl;
100}

◆ G4GeneralPhaseSpaceDecay() [4/4]

G4GeneralPhaseSpaceDecay::G4GeneralPhaseSpaceDecay ( const G4String &  theParentName,
G4double  theParentMass,
G4double  theBR,
G4int  theNumberOfDaughters,
const G4String &  theDaughterName1,
const G4String &  theDaughterName2,
const G4String &  theDaughterName3,
const G4double *  masses 
)

Definition at line 102 of file G4GeneralPhaseSpaceDecay.cc.

109 :
110 G4VDecayChannel("Phase Space",
111 theParentName,theBR,
112 theNumberOfDaughters,
113 theDaughterName1,
114 theDaughterName2,
115 theDaughterName3),
116 parentmass(theParentMass),
117 theDaughterMasses(masses)
118{
119 if (GetVerboseLevel()>1) G4cout << "G4GeneralPhaseSpaceDecay:: constructor " << G4endl;
120}

◆ ~G4GeneralPhaseSpaceDecay()

G4GeneralPhaseSpaceDecay::~G4GeneralPhaseSpaceDecay ( )
virtual

Definition at line 122 of file G4GeneralPhaseSpaceDecay.cc.

123{
124}

Member Function Documentation

◆ DecayIt()

G4DecayProducts * G4GeneralPhaseSpaceDecay::DecayIt ( G4double  mass = 0.0)
virtual

Implements G4VDecayChannel.

Reimplemented in G4NuclearDecayChannel.

Definition at line 126 of file G4GeneralPhaseSpaceDecay.cc.

127{
128 if (GetVerboseLevel()>1) G4cout << "G4GeneralPhaseSpaceDecay::DecayIt ";
129 G4DecayProducts * products = NULL;
130
131 if (parent == NULL) FillParent();
132 if (daughters == NULL) FillDaughters();
133
134 switch (numberOfDaughters){
135 case 0:
136 if (GetVerboseLevel()>0) {
137 G4cout << "G4GeneralPhaseSpaceDecay::DecayIt ";
138 G4cout << " daughters not defined " <<G4endl;
139 }
140 break;
141 case 1:
142 products = OneBodyDecayIt();
143 break;
144 case 2:
145 products = TwoBodyDecayIt();
146 break;
147 case 3:
148 products = ThreeBodyDecayIt();
149 break;
150 default:
151 products = ManyBodyDecayIt();
152 break;
153 }
154 if ((products == NULL) && (GetVerboseLevel()>0)) {
155 G4cout << "G4GeneralPhaseSpaceDecay::DecayIt ";
156 G4cout << *parent_name << " can not decay " << G4endl;
157 DumpInfo();
158 }
159 return products;
160}
G4VDecayChannel::parent_name
G4String * parent_name
Definition: G4VDecayChannel.hh:133
G4VDecayChannel::daughters
G4ParticleDefinition ** daughters
Definition: G4VDecayChannel.hh:147

Referenced by G4KineticTrack::Decay().

◆ GetParentMass()

G4double G4GeneralPhaseSpaceDecay::GetParentMass ( ) const
inline

Definition at line 98 of file G4GeneralPhaseSpaceDecay.hh.

99{
100 return parentmass;
101}

◆ ManyBodyDecayIt()

G4DecayProducts * G4GeneralPhaseSpaceDecay::ManyBodyDecayIt ( )
protected

Definition at line 353 of file G4GeneralPhaseSpaceDecay.cc.

363{
364 //return value
365 G4DecayProducts *products;
366
367 if (GetVerboseLevel()>1) G4cout << "G4GeneralPhaseSpaceDecay::ManyBodyDecayIt()"<<G4endl;
368
369 //daughters'mass
370 G4double *daughtermass = new G4double[numberOfDaughters];
371 G4double sumofdaughtermass = 0.0;
372 for (G4int index=0; index<numberOfDaughters; index++){
373 daughtermass[index] = daughters[index]->GetPDGMass();
374 sumofdaughtermass += daughtermass[index];
375 }
376
377 //Calculate daughter momentum
378 G4double *daughtermomentum = new G4double[numberOfDaughters];
379 G4ParticleMomentum direction;
380 G4DynamicParticle **daughterparticle;
381 G4double *sm = new G4double[numberOfDaughters];
382 G4double tmas;
383 G4double weight = 1.0;
384 G4int numberOfTry = 0;
385 G4int index1;
386
387 do {
388 //Generate rundom number in descending order
389 G4double temp;
390 G4double *rd = new G4double[numberOfDaughters];
391 rd[0] = 1.0;
392 for(index1 =1; index1 < numberOfDaughters -1; index1++)
393 rd[index1] = G4UniformRand();
394 rd[ numberOfDaughters -1] = 0.0;
395 for(index1 =1; index1 < numberOfDaughters -1; index1++) {
396 for(G4int index2 = index1+1; index2 < numberOfDaughters; index2++) {
397 if (rd[index1] < rd[index2]){
398 temp = rd[index1];
399 rd[index1] = rd[index2];
400 rd[index2] = temp;
401 }
402 }
403 }
404
405 //calcurate virtual mass
406 tmas = parentmass - sumofdaughtermass;
407 temp = sumofdaughtermass;
408 for(index1 =0; index1 < numberOfDaughters; index1++) {
409 sm[index1] = rd[index1]*tmas + temp;
410 temp -= daughtermass[index1];
411 if (GetVerboseLevel()>1) {
412 G4cout << index1 << " rundom number:" << rd[index1];
413 G4cout << " virtual mass:" << sm[index1]/GeV << "[GeV/c/c]" <<G4endl;
414 }
415 }
416 delete [] rd;
417
418 //Calculate daughter momentum
419 weight = 1.0;
420 index1 =numberOfDaughters-1;
421 daughtermomentum[index1]= Pmx( sm[index1-1],daughtermass[index1-1],sm[index1]);
422 if (GetVerboseLevel()>1) {
423 G4cout << " daughter " << index1 << ":" << *daughters_name[index1];
424 G4cout << " momentum:" << daughtermomentum[index1]/GeV << "[GeV/c]" <<G4endl;
425 }
426 for(index1 =numberOfDaughters-2; index1>=0; index1--) {
427 // calculate
428 daughtermomentum[index1]= Pmx( sm[index1],daughtermass[index1], sm[index1 +1]);
429 if(daughtermomentum[index1] < 0.0) {
430 // !!! illegal momentum !!!
431 if (GetVerboseLevel()>0) {
432 G4cout << "G4GeneralPhaseSpaceDecay::ManyBodyDecayIt ";
433 G4cout << " can not calculate daughter momentum " <<G4endl;
434 G4cout << " parent:" << *parent_name;
435 G4cout << " mass:" << parentmass/GeV << "[GeV/c/c]" <<G4endl;
436 G4cout << " daughter " << index1 << ":" << *daughters_name[index1];
437 G4cout << " mass:" << daughtermass[index1]/GeV << "[GeV/c/c]" ;
438 G4cout << " mass:" << daughtermomentum[index1]/GeV << "[GeV/c]" <<G4endl;
439 }
440 delete [] sm;
441 delete [] daughtermass;
442 delete [] daughtermomentum;
443 return NULL; // Error detection
444
445 } else {
446 // calculate weight of this events
447 weight *= daughtermomentum[index1]/sm[index1];
448 if (GetVerboseLevel()>1) {
449 G4cout << " daughter " << index1 << ":" << *daughters_name[index1];
450 G4cout << " momentum:" << daughtermomentum[index1]/GeV << "[GeV/c]" <<G4endl;
451 }
452 }
453 }
454 if (GetVerboseLevel()>1) {
455 G4cout << " weight: " << weight <<G4endl;
456 }
457
458 // exit if number of Try exceeds 100
459 if (numberOfTry++ >100) {
460 if (GetVerboseLevel()>0) {
461 G4cout << "G4GeneralPhaseSpaceDecay::ManyBodyDecayIt: ";
462 G4cout << " can not determine Decay Kinematics " << G4endl;
463 }
464 delete [] sm;
465 delete [] daughtermass;
466 delete [] daughtermomentum;
467 return NULL; // Error detection
468 }
469 } while ( weight > G4UniformRand());
470 if (GetVerboseLevel()>1) {
471 G4cout << "Start calulation of daughters momentum vector "<<G4endl;
472 }
473
474 G4double costheta, sintheta, phi;
475 G4double beta;
476 daughterparticle = new G4DynamicParticle*[numberOfDaughters];
477
478 index1 = numberOfDaughters -2;
479 costheta = 2.*G4UniformRand()-1.0;
480 sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
481 phi = twopi*G4UniformRand()*rad;
482 direction.setZ(costheta);
483 direction.setY(sintheta*std::sin(phi));
484 direction.setX(sintheta*std::cos(phi));
485 daughterparticle[index1] = new G4DynamicParticle( daughters[index1], direction*daughtermomentum[index1] );
486 daughterparticle[index1+1] = new G4DynamicParticle( daughters[index1+1], direction*(-1.0*daughtermomentum[index1]) );
487
488 for (index1 = numberOfDaughters -3; index1 >= 0; index1--) {
489 //calculate momentum direction
490 costheta = 2.*G4UniformRand()-1.0;
491 sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
492 phi = twopi*G4UniformRand()*rad;
493 direction.setZ(costheta);
494 direction.setY(sintheta*std::sin(phi));
495 direction.setX(sintheta*std::cos(phi));
496
497 // boost already created particles
498 beta = daughtermomentum[index1];
499 beta /= std::sqrt( daughtermomentum[index1]*daughtermomentum[index1] + sm[index1+1]*sm[index1+1] );
500 for (G4int index2 = index1+1; index2<numberOfDaughters; index2++) {
501 G4LorentzVector p4;
502 // make G4LorentzVector for secondaries
503 p4 = daughterparticle[index2]->Get4Momentum();
504
505 // boost secondaries to new frame
506 p4.boost( direction.x()*beta, direction.y()*beta, direction.z()*beta);
507
508 // change energy/momentum
509 daughterparticle[index2]->Set4Momentum(p4);
510 }
511 //create daughter G4DynamicParticle
512 daughterparticle[index1]= new G4DynamicParticle( daughters[index1], direction*(-1.0*daughtermomentum[index1]));
513 }
514
515 //create G4Decayproducts
516 G4DynamicParticle *parentparticle;
517 direction.setX(1.0); direction.setY(0.0); direction.setZ(0.0);
518 parentparticle = new G4DynamicParticle( parent, direction, 0.0);
519 products = new G4DecayProducts(*parentparticle);
520 delete parentparticle;
521 for (index1 = 0; index1<numberOfDaughters; index1++) {
522 products->PushProducts(daughterparticle[index1]);
523 }
524 if (GetVerboseLevel()>1) {
525 G4cout << "G4GeneralPhaseSpaceDecay::ManyBodyDecayIt ";
526 G4cout << " create decay products in rest frame " << G4endl;
527 products->DumpInfo();
528 }
529
530 delete [] daughterparticle;
531 delete [] daughtermomentum;
532 delete [] daughtermass;
533 delete [] sm;
534
535 return products;
536}
G4double
double G4double
Definition: G4Types.hh:64
G4int
int G4int
Definition: G4Types.hh:66
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:53
CLHEP::Hep3Vector::z
double z() const
CLHEP::Hep3Vector::x
double x() const
CLHEP::Hep3Vector::setY
void setY(double)
CLHEP::Hep3Vector::y
double y() const
CLHEP::Hep3Vector::setZ
void setZ(double)
CLHEP::Hep3Vector::setX
void setX(double)
CLHEP::HepLorentzVector::boost
HepLorentzVector & boost(double, double, double)
Definition: LorentzVector.cc:59
G4DecayProducts::DumpInfo
void DumpInfo() const
Definition: G4DecayProducts.cc:311
G4DecayProducts::PushProducts
G4int PushProducts(G4DynamicParticle *aParticle)
Definition: G4DecayProducts.cc:153
G4DynamicParticle::Get4Momentum
G4LorentzVector Get4Momentum() const
G4DynamicParticle::Set4Momentum
void Set4Momentum(const G4LorentzVector &momentum)
Definition: G4DynamicParticle.cc:351
G4GeneralPhaseSpaceDecay::Pmx
static G4double Pmx(G4double e, G4double p1, G4double p2)
Definition: G4GeneralPhaseSpaceDecay.hh:111
G4VDecayChannel::daughters_name
G4String ** daughters_name
Definition: G4VDecayChannel.hh:135

Referenced by DecayIt().

◆ OneBodyDecayIt()

G4DecayProducts * G4GeneralPhaseSpaceDecay::OneBodyDecayIt ( )
protected

Definition at line 162 of file G4GeneralPhaseSpaceDecay.cc.

163{
164 if (GetVerboseLevel()>1) G4cout << "G4GeneralPhaseSpaceDecay::OneBodyDecayIt()"<<G4endl;
165
166// G4double daughtermass = daughters[0]->GetPDGMass();
167
168 //create parent G4DynamicParticle at rest
169 G4ParticleMomentum dummy;
170 G4DynamicParticle * parentparticle = new G4DynamicParticle(parent, dummy, 0.0);
171
172 //create G4Decayproducts
173 G4DecayProducts *products = new G4DecayProducts(*parentparticle);
174 delete parentparticle;
175
176 //create daughter G4DynamicParticle at rest
177 G4DynamicParticle * daughterparticle = new G4DynamicParticle(daughters[0], dummy, 0.0);
178 products->PushProducts(daughterparticle);
179
180 if (GetVerboseLevel()>1)
181 {
182 G4cout << "G4GeneralPhaseSpaceDecay::OneBodyDecayIt ";
183 G4cout << " create decay products in rest frame " <<G4endl;
184 products->DumpInfo();
185 }
186 return products;
187}

Referenced by DecayIt(), and G4NuclearDecayChannel::DecayIt().

◆ Pmx()

G4double G4GeneralPhaseSpaceDecay::Pmx ( G4double  e,
G4double  p1,
G4double  p2 
)
inlinestatic

Definition at line 111 of file G4GeneralPhaseSpaceDecay.hh.

112{
113 // calculate momentum of daughter particles in two-body decay
114 if (e-p1-p2 < 0 )
115 {
116 G4HadronicException(__FILE__, __LINE__, "G4GeneralPhaseSpaceDecay::Pmx energy in cms > mass1+mass2");
117 }
118 G4double ppp = (e+p1+p2)*(e+p1-p2)*(e-p1+p2)*(e-p1-p2)/(4.0*e*e);
119 if (ppp>0) return std::sqrt(ppp);
120 else return -1.;
121}

Referenced by ManyBodyDecayIt(), and TwoBodyDecayIt().

◆ SetParentMass()

void G4GeneralPhaseSpaceDecay::SetParentMass ( const G4double  aParentMass)
inline

Definition at line 103 of file G4GeneralPhaseSpaceDecay.hh.

104{
105 parentmass = aParentMass;
106}

Referenced by G4NuclearDecayChannel::DecayIt().

◆ ThreeBodyDecayIt()

G4DecayProducts * G4GeneralPhaseSpaceDecay::ThreeBodyDecayIt ( )
protected

Definition at line 239 of file G4GeneralPhaseSpaceDecay.cc.

241{
242 if (GetVerboseLevel()>1) G4cout << "G4GeneralPhaseSpaceDecay::ThreeBodyDecayIt()"<<G4endl;
243
244 //daughters'mass
245 G4double daughtermass[3];
246 G4double sumofdaughtermass = 0.0;
247 for (G4int index=0; index<3; index++)
248 {
249 if ( theDaughterMasses )
250 {
251 daughtermass[index]= *(theDaughterMasses+index);
252 } else {
253 daughtermass[index] = daughters[index]->GetPDGMass();
254 }
255 sumofdaughtermass += daughtermass[index];
256 }
257
258 //create parent G4DynamicParticle at rest
259 G4ParticleMomentum dummy;
260 G4DynamicParticle * parentparticle = new G4DynamicParticle( parent, dummy, 0.0);
261
262 //create G4Decayproducts
263 G4DecayProducts *products = new G4DecayProducts(*parentparticle);
264 delete parentparticle;
265
266 //calculate daughter momentum
267 // Generate two
268 G4double rd1, rd2, rd;
269 G4double daughtermomentum[3];
270 G4double momentummax=0.0, momentumsum = 0.0;
271 G4double energy;
272
273 do
274 {
275 rd1 = G4UniformRand();
276 rd2 = G4UniformRand();
277 if (rd2 > rd1)
278 {
279 rd = rd1;
280 rd1 = rd2;
281 rd2 = rd;
282 }
283 momentummax = 0.0;
284 momentumsum = 0.0;
285 // daughter 0
286
287 energy = rd2*(parentmass - sumofdaughtermass);
288 daughtermomentum[0] = std::sqrt(energy*energy + 2.0*energy* daughtermass[0]);
289 if ( daughtermomentum[0] >momentummax )momentummax = daughtermomentum[0];
290 momentumsum += daughtermomentum[0];
291
292 // daughter 1
293 energy = (1.-rd1)*(parentmass - sumofdaughtermass);
294 daughtermomentum[1] = std::sqrt(energy*energy + 2.0*energy* daughtermass[1]);
295 if ( daughtermomentum[1] >momentummax )momentummax = daughtermomentum[1];
296 momentumsum += daughtermomentum[1];
297
298 // daughter 2
299 energy = (rd1-rd2)*(parentmass - sumofdaughtermass);
300 daughtermomentum[2] = std::sqrt(energy*energy + 2.0*energy* daughtermass[2]);
301 if ( daughtermomentum[2] >momentummax )momentummax = daughtermomentum[2];
302 momentumsum += daughtermomentum[2];
303 } while (momentummax > momentumsum - momentummax );
304
305 // output message
306 if (GetVerboseLevel()>1) {
307 G4cout << " daughter 0:" << daughtermomentum[0]/GeV << "[GeV/c]" <<G4endl;
308 G4cout << " daughter 1:" << daughtermomentum[1]/GeV << "[GeV/c]" <<G4endl;
309 G4cout << " daughter 2:" << daughtermomentum[2]/GeV << "[GeV/c]" <<G4endl;
310 G4cout << " momentum sum:" << momentumsum/GeV << "[GeV/c]" <<G4endl;
311 }
312
313 //create daughter G4DynamicParticle
314 G4double costheta, sintheta, phi, sinphi, cosphi;
315 G4double costhetan, sinthetan, phin, sinphin, cosphin;
316 costheta = 2.*G4UniformRand()-1.0;
317 sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
318 phi = twopi*G4UniformRand()*rad;
319 sinphi = std::sin(phi);
320 cosphi = std::cos(phi);
321 G4ParticleMomentum direction0(sintheta*cosphi,sintheta*sinphi,costheta);
322 G4double Etotal=std::sqrt( daughtermass[0]*daughtermass[0] + daughtermomentum[0]*daughtermomentum[0]);
323 G4DynamicParticle * daughterparticle
324 = new G4DynamicParticle( daughters[0], Etotal, direction0*daughtermomentum[0]);
325 products->PushProducts(daughterparticle);
326
327 costhetan = (daughtermomentum[1]*daughtermomentum[1]-daughtermomentum[2]*daughtermomentum[2]-daughtermomentum[0]*daughtermomentum[0])/(2.0*daughtermomentum[2]*daughtermomentum[0]);
328 sinthetan = std::sqrt((1.0-costhetan)*(1.0+costhetan));
329 phin = twopi*G4UniformRand()*rad;
330 sinphin = std::sin(phin);
331 cosphin = std::cos(phin);
332 G4ParticleMomentum direction2;
333 direction2.setX( sinthetan*cosphin*costheta*cosphi - sinthetan*sinphin*sinphi + costhetan*sintheta*cosphi);
334 direction2.setY( sinthetan*cosphin*costheta*sinphi + sinthetan*sinphin*cosphi + costhetan*sintheta*sinphi);
335 direction2.setZ( -sinthetan*cosphin*sintheta + costhetan*costheta);
336 Etotal=std::sqrt( daughtermass[2]*daughtermass[2] + daughtermomentum[2]*daughtermomentum[2]/direction2.mag2());
337 daughterparticle = new G4DynamicParticle( daughters[2],Etotal, direction2*(daughtermomentum[2]/direction2.mag()));
338 products->PushProducts(daughterparticle);
339 G4ThreeVector mom=(direction0*daughtermomentum[0] + direction2*(daughtermomentum[2]/direction2.mag()))*(-1.0);
340 Etotal= std::sqrt( daughtermass[1]*daughtermass[1] + mom.mag2() );
341 daughterparticle =
342 new G4DynamicParticle(daughters[1], Etotal, mom);
343 products->PushProducts(daughterparticle);
344
345 if (GetVerboseLevel()>1) {
346 G4cout << "G4GeneralPhaseSpaceDecay::ThreeBodyDecayIt ";
347 G4cout << " create decay products in rest frame " <<G4endl;
348 products->DumpInfo();
349 }
350 return products;
351}
CLHEP::Hep3Vector::mag2
double mag2() const
CLHEP::Hep3Vector::mag
double mag() const

Referenced by DecayIt().

◆ TwoBodyDecayIt()

G4DecayProducts * G4GeneralPhaseSpaceDecay::TwoBodyDecayIt ( )
protected

Definition at line 189 of file G4GeneralPhaseSpaceDecay.cc.

190{
191 if (GetVerboseLevel()>1) G4cout << "G4GeneralPhaseSpaceDecay::TwoBodyDecayIt()"<<G4endl;
192
193 //daughters'mass
194 G4double daughtermass[2];
195 G4double daughtermomentum;
196 if ( theDaughterMasses )
197 {
198 daughtermass[0]= *(theDaughterMasses);
199 daughtermass[1] = *(theDaughterMasses+1);
200 } else {
201 daughtermass[0] = daughters[0]->GetPDGMass();
202 daughtermass[1] = daughters[1]->GetPDGMass();
203 }
204
205// G4double sumofdaughtermass = daughtermass[0] + daughtermass[1];
206
207 //create parent G4DynamicParticle at rest
208 G4ParticleMomentum dummy;
209 G4DynamicParticle * parentparticle = new G4DynamicParticle( parent, dummy, 0.0);
210
211 //create G4Decayproducts @@GF why dummy parentparticle?
212 G4DecayProducts *products = new G4DecayProducts(*parentparticle);
213 delete parentparticle;
214
215 //calculate daughter momentum
216 daughtermomentum = Pmx(parentmass,daughtermass[0],daughtermass[1]);
217 G4double costheta = 2.*G4UniformRand()-1.0;
218 G4double sintheta = std::sqrt((1.0 - costheta)*(1.0 + costheta));
219 G4double phi = twopi*G4UniformRand()*rad;
220 G4ParticleMomentum direction(sintheta*std::cos(phi),sintheta*std::sin(phi),costheta);
221
222 //create daughter G4DynamicParticle
223 G4double Etotal= std::sqrt(daughtermass[0]*daughtermass[0] + daughtermomentum*daughtermomentum);
224 G4DynamicParticle * daughterparticle = new G4DynamicParticle( daughters[0],Etotal, direction*daughtermomentum);
225 products->PushProducts(daughterparticle);
226 Etotal= std::sqrt(daughtermass[1]*daughtermass[1] + daughtermomentum*daughtermomentum);
227 daughterparticle = new G4DynamicParticle( daughters[1],Etotal, direction*(-1.0*daughtermomentum));
228 products->PushProducts(daughterparticle);
229
230 if (GetVerboseLevel()>1)
231 {
232 G4cout << "G4GeneralPhaseSpaceDecay::TwoBodyDecayIt ";
233 G4cout << " create decay products in rest frame " <<G4endl;
234 products->DumpInfo();
235 }
236 return products;
237}

Referenced by DecayIt(), and G4NuclearDecayChannel::DecayIt().

The documentation for this class was generated from the following files: