Geant4 10.7.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)
 
 G4GeneralPhaseSpaceDecay (const G4String &theParentName, G4double theParentMass, G4double theBR, G4int theNumberOfDaughters, const G4String &theDaughterName1, const G4String &theDaughterName2, const G4String &theDaughterName3, const G4String &theDaughterName4, 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 ()
 
G4bool operator== (const G4VDecayChannel &r) const
 
G4bool operator!= (const G4VDecayChannel &r) const
 
G4bool 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 ()
 
G4double GetRangeMass () const
 
void SetRangeMass (G4double val)
 
virtual G4bool IsOKWithParentMass (G4double parentMass)
 
void SetPolarization (const G4ThreeVector &)
 
const G4ThreeVector & GetPolarization () const
 

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
void ClearDaughtersName ()
 
void CheckAndFillDaughters ()
 
void CheckAndFillParent ()
 
G4double DynamicalMass (G4double massPDG, G4double width, G4double maxDev=1.0) const
 
 G4VDecayChannel ()
 
 G4VDecayChannel (const G4VDecayChannel &)
 
G4VDecayChannel & operator= (const G4VDecayChannel &)
 

Additional Inherited Members

- Protected Attributes inherited from G4VDecayChannel
G4String kinematics_name = ""
 
G4double rbranch = 0.0
 
G4String * parent_name = nullptr
 
G4String ** daughters_name = nullptr
 
G4double rangeMass = 2.5
 
G4ThreeVector parent_polarization
 
G4ParticleTable * particletable = nullptr
 
G4ParticleDefinition * G4MT_parent = nullptr
 
G4ParticleDefinition ** G4MT_daughters = nullptr
 
G4double G4MT_parent_mass = 0.0
 
G4double * G4MT_daughters_mass = nullptr
 
G4double * G4MT_daughters_width = nullptr
 
G4Mutex daughtersMutex
 
G4Mutex parentMutex
 
G4int numberOfDaughters = 0
 
G4int verboseLevel = 1
 
- Static Protected Attributes inherited from G4VDecayChannel
static const G4String noName = " "
 

Detailed Description

Definition at line 44 of file G4GeneralPhaseSpaceDecay.hh.

Constructor & Destructor Documentation

◆ G4GeneralPhaseSpaceDecay() [1/5]

G4GeneralPhaseSpaceDecay::G4GeneralPhaseSpaceDecay ( G4int  Verbose = 1)

Definition at line 48 of file G4GeneralPhaseSpaceDecay.cc.

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

◆ G4GeneralPhaseSpaceDecay() [2/5]

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

Definition at line 55 of file G4GeneralPhaseSpaceDecay.cc.

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

◆ G4GeneralPhaseSpaceDecay() [3/5]

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

Definition at line 81 of file G4GeneralPhaseSpaceDecay.cc.

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

◆ G4GeneralPhaseSpaceDecay() [4/5]

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 100 of file G4GeneralPhaseSpaceDecay.cc.

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

◆ G4GeneralPhaseSpaceDecay() [5/5]

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

Definition at line 120 of file G4GeneralPhaseSpaceDecay.cc.

128 :
129 G4VDecayChannel("Phase Space",
130 theParentName,theBR,
131 theNumberOfDaughters,
132 theDaughterName1,
133 theDaughterName2,
134 theDaughterName3,
135 theDaughterName4),
136 parentmass(theParentMass),
137 theDaughterMasses(masses)
138{
139 if (GetVerboseLevel()>1) G4cout << "G4GeneralPhaseSpaceDecay:: constructor " << G4endl;
140}

◆ ~G4GeneralPhaseSpaceDecay()

G4GeneralPhaseSpaceDecay::~G4GeneralPhaseSpaceDecay ( )
virtual

Definition at line 142 of file G4GeneralPhaseSpaceDecay.cc.

143{
144}

Member Function Documentation

◆ DecayIt()

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

Implements G4VDecayChannel.

Definition at line 146 of file G4GeneralPhaseSpaceDecay.cc.

147{
148 if (GetVerboseLevel()>1) G4cout << "G4GeneralPhaseSpaceDecay::DecayIt ";
149 G4DecayProducts * products = NULL;
150
151 CheckAndFillParent();
152 CheckAndFillDaughters();
153
154 switch (numberOfDaughters){
155 case 0:
156 if (GetVerboseLevel()>0) {
157 G4cout << "G4GeneralPhaseSpaceDecay::DecayIt ";
158 G4cout << " daughters not defined " <<G4endl;
159 }
160 break;
161 case 1:
162 products = OneBodyDecayIt();
163 break;
164 case 2:
165 products = TwoBodyDecayIt();
166 break;
167 case 3:
168 products = ThreeBodyDecayIt();
169 break;
170 default:
171 products = ManyBodyDecayIt();
172 break;
173 }
174 if ((products == NULL) && (GetVerboseLevel()>0)) {
175 G4cout << "G4GeneralPhaseSpaceDecay::DecayIt ";
176 G4cout << *parent_name << " can not decay " << G4endl;
177 DumpInfo();
178 }
179 return products;
180}
G4VDecayChannel::parent_name
G4String * parent_name
Definition: G4VDecayChannel.hh:155
G4VDecayChannel::CheckAndFillDaughters
void CheckAndFillDaughters()
Definition: G4VDecayChannel.hh:336

Referenced by G4KineticTrack::Decay().

◆ GetParentMass()

G4double G4GeneralPhaseSpaceDecay::GetParentMass ( ) const
inline

Definition at line 107 of file G4GeneralPhaseSpaceDecay.hh.

108{
109 return parentmass;
110}

◆ ManyBodyDecayIt()

G4DecayProducts * G4GeneralPhaseSpaceDecay::ManyBodyDecayIt ( )
protected

Definition at line 381 of file G4GeneralPhaseSpaceDecay.cc.

391{
392 //return value
393 G4DecayProducts *products;
394
395 if (GetVerboseLevel()>1) G4cout << "G4GeneralPhaseSpaceDecay::ManyBodyDecayIt()"<<G4endl;
396
397 //daughters'mass
398 G4double *daughtermass = new G4double[numberOfDaughters];
399 G4double sumofdaughtermass = 0.0;
400 for (G4int index=0; index<numberOfDaughters; index++){
401 daughtermass[index] = G4MT_daughters[index]->GetPDGMass();
402 sumofdaughtermass += daughtermass[index];
403 }
404
405 //Calculate daughter momentum
406 G4double *daughtermomentum = new G4double[numberOfDaughters];
407 G4ParticleMomentum direction;
408 G4DynamicParticle **daughterparticle;
409 G4double *sm = new G4double[numberOfDaughters];
410 G4double tmas;
411 G4double weight = 1.0;
412 G4int numberOfTry = 0;
413 G4int index1;
414
415 do {
416 //Generate rundom number in descending order
417 G4double temp;
418 G4double *rd = new G4double[numberOfDaughters];
419 rd[0] = 1.0;
420 for(index1 =1; index1 < numberOfDaughters -1; index1++)
421 rd[index1] = G4UniformRand();
422 rd[ numberOfDaughters -1] = 0.0;
423 for(index1 =1; index1 < numberOfDaughters -1; index1++) {
424 for(G4int index2 = index1+1; index2 < numberOfDaughters; index2++) {
425 if (rd[index1] < rd[index2]){
426 temp = rd[index1];
427 rd[index1] = rd[index2];
428 rd[index2] = temp;
429 }
430 }
431 }
432
433 //calcurate virtual mass
434 tmas = parentmass - sumofdaughtermass;
435 temp = sumofdaughtermass;
436 for(index1 =0; index1 < numberOfDaughters; index1++) {
437 sm[index1] = rd[index1]*tmas + temp;
438 temp -= daughtermass[index1];
439 if (GetVerboseLevel()>1) {
440 G4cout << index1 << " rundom number:" << rd[index1];
441 G4cout << " virtual mass:" << sm[index1]/GeV << "[GeV/c/c]" <<G4endl;
442 }
443 }
444 delete [] rd;
445
446 //Calculate daughter momentum
447 weight = 1.0;
448 index1 =numberOfDaughters-1;
449 daughtermomentum[index1]= Pmx( sm[index1-1],daughtermass[index1-1],sm[index1]);
450 if (GetVerboseLevel()>1) {
451 G4cout << " daughter " << index1 << ":" << *daughters_name[index1];
452 G4cout << " momentum:" << daughtermomentum[index1]/GeV << "[GeV/c]" <<G4endl;
453 }
454 for(index1 =numberOfDaughters-2; index1>=0; index1--) {
455 // calculate
456 daughtermomentum[index1]= Pmx( sm[index1],daughtermass[index1], sm[index1 +1]);
457 if(daughtermomentum[index1] < 0.0) {
458 // !!! illegal momentum !!!
459 if (GetVerboseLevel()>0) {
460 G4cout << "G4GeneralPhaseSpaceDecay::ManyBodyDecayIt ";
461 G4cout << " can not calculate daughter momentum " <<G4endl;
462 G4cout << " parent:" << *parent_name;
463 G4cout << " mass:" << parentmass/GeV << "[GeV/c/c]" <<G4endl;
464 G4cout << " daughter " << index1 << ":" << *daughters_name[index1];
465 G4cout << " mass:" << daughtermass[index1]/GeV << "[GeV/c/c]" ;
466 G4cout << " mass:" << daughtermomentum[index1]/GeV << "[GeV/c]" <<G4endl;
467 }
468 delete [] sm;
469 delete [] daughtermass;
470 delete [] daughtermomentum;
471 return NULL; // Error detection
472
473 } else {
474 // calculate weight of this events
475 weight *= daughtermomentum[index1]/sm[index1];
476 if (GetVerboseLevel()>1) {
477 G4cout << " daughter " << index1 << ":" << *daughters_name[index1];
478 G4cout << " momentum:" << daughtermomentum[index1]/GeV << "[GeV/c]" <<G4endl;
479 }
480 }
481 }
482 if (GetVerboseLevel()>1) {
483 G4cout << " weight: " << weight <<G4endl;
484 }
485
486 // exit if number of Try exceeds 100
487 if (numberOfTry++ >100) {
488 if (GetVerboseLevel()>0) {
489 G4cout << "G4GeneralPhaseSpaceDecay::ManyBodyDecayIt: ";
490 G4cout << " can not determine Decay Kinematics " << G4endl;
491 }
492 delete [] sm;
493 delete [] daughtermass;
494 delete [] daughtermomentum;
495 return NULL; // Error detection
496 }
497 } while ( weight > G4UniformRand()); /* Loop checking, 02.11.2015, A.Ribon */
498 if (GetVerboseLevel()>1) {
499 G4cout << "Start calculation of daughters momentum vector "<<G4endl;
500 }
501
502 G4double costheta, sintheta, phi;
503 G4double beta;
504 daughterparticle = new G4DynamicParticle*[numberOfDaughters];
505
506 index1 = numberOfDaughters -2;
507 costheta = 2.*G4UniformRand()-1.0;
508 sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
509 phi = twopi*G4UniformRand()*rad;
510 direction.setZ(costheta);
511 direction.setY(sintheta*std::sin(phi));
512 direction.setX(sintheta*std::cos(phi));
513 daughterparticle[index1] = new G4DynamicParticle( G4MT_daughters[index1], direction*daughtermomentum[index1] );
514 daughterparticle[index1+1] = new G4DynamicParticle( G4MT_daughters[index1+1], direction*(-1.0*daughtermomentum[index1]) );
515
516 for (index1 = numberOfDaughters -3; index1 >= 0; index1--) {
517 //calculate momentum direction
518 costheta = 2.*G4UniformRand()-1.0;
519 sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
520 phi = twopi*G4UniformRand()*rad;
521 direction.setZ(costheta);
522 direction.setY(sintheta*std::sin(phi));
523 direction.setX(sintheta*std::cos(phi));
524
525 // boost already created particles
526 beta = daughtermomentum[index1];
527 beta /= std::sqrt( daughtermomentum[index1]*daughtermomentum[index1] + sm[index1+1]*sm[index1+1] );
528 for (G4int index2 = index1+1; index2<numberOfDaughters; index2++) {
529 G4LorentzVector p4;
530 // make G4LorentzVector for secondaries
531 p4 = daughterparticle[index2]->Get4Momentum();
532
533 // boost secondaries to new frame
534 p4.boost( direction.x()*beta, direction.y()*beta, direction.z()*beta);
535
536 // change energy/momentum
537 daughterparticle[index2]->Set4Momentum(p4);
538 }
539 //create daughter G4DynamicParticle
540 daughterparticle[index1]= new G4DynamicParticle( G4MT_daughters[index1], direction*(-1.0*daughtermomentum[index1]));
541 }
542
543 //create G4Decayproducts
544 G4DynamicParticle *parentparticle;
545 direction.setX(1.0); direction.setY(0.0); direction.setZ(0.0);
546 parentparticle = new G4DynamicParticle( G4MT_parent, direction, 0.0);
547 products = new G4DecayProducts(*parentparticle);
548 delete parentparticle;
549 for (index1 = 0; index1<numberOfDaughters; index1++) {
550 products->PushProducts(daughterparticle[index1]);
551 }
552 if (GetVerboseLevel()>1) {
553 G4cout << "G4GeneralPhaseSpaceDecay::ManyBodyDecayIt ";
554 G4cout << " create decay products in rest frame " << G4endl;
555 products->DumpInfo();
556 }
557
558 delete [] daughterparticle;
559 delete [] daughtermomentum;
560 delete [] daughtermass;
561 delete [] sm;
562
563 return products;
564}
G4double
double G4double
Definition: G4Types.hh:83
G4int
int G4int
Definition: G4Types.hh:85
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::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:55
G4DecayProducts::DumpInfo
void DumpInfo() const
Definition: G4DecayProducts.cc:334
G4DecayProducts::PushProducts
G4int PushProducts(G4DynamicParticle *aParticle)
Definition: G4DecayProducts.cc:151
G4DynamicParticle::Get4Momentum
G4LorentzVector Get4Momentum() const
G4DynamicParticle::Set4Momentum
void Set4Momentum(const G4LorentzVector &momentum)
Definition: G4DynamicParticle.cc:410
G4GeneralPhaseSpaceDecay::Pmx
static G4double Pmx(G4double e, G4double p1, G4double p2)
Definition: G4GeneralPhaseSpaceDecay.hh:120
G4VDecayChannel::G4MT_daughters
G4ParticleDefinition ** G4MT_daughters
Definition: G4VDecayChannel.hh:172
G4VDecayChannel::daughters_name
G4String ** daughters_name
Definition: G4VDecayChannel.hh:157

Referenced by DecayIt().

◆ OneBodyDecayIt()

G4DecayProducts * G4GeneralPhaseSpaceDecay::OneBodyDecayIt ( )
protected

Definition at line 182 of file G4GeneralPhaseSpaceDecay.cc.

183{
184 if (GetVerboseLevel()>1) G4cout << "G4GeneralPhaseSpaceDecay::OneBodyDecayIt()"<<G4endl;
185
186// G4double daughtermass = daughters[0]->GetPDGMass();
187
188 //create parent G4DynamicParticle at rest
189 G4ParticleMomentum dummy;
190 G4DynamicParticle * parentparticle = new G4DynamicParticle(G4MT_parent, dummy, 0.0);
191
192 //create G4Decayproducts
193 G4DecayProducts *products = new G4DecayProducts(*parentparticle);
194 delete parentparticle;
195
196 //create daughter G4DynamicParticle at rest
197 G4DynamicParticle * daughterparticle = new G4DynamicParticle(G4MT_daughters[0], dummy, 0.0);
198 products->PushProducts(daughterparticle);
199
200 if (GetVerboseLevel()>1)
201 {
202 G4cout << "G4GeneralPhaseSpaceDecay::OneBodyDecayIt ";
203 G4cout << " create decay products in rest frame " <<G4endl;
204 products->DumpInfo();
205 }
206 return products;
207}

Referenced by DecayIt().

◆ Pmx()

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

Definition at line 120 of file G4GeneralPhaseSpaceDecay.hh.

121{
122 // calculate momentum of daughter particles in two-body decay
123 if (e-p1-p2 < 0 )
124 {
125 throw G4HadronicException(__FILE__, __LINE__, "G4GeneralPhaseSpaceDecay::Pmx energy in cms < mass1+mass2");
126 }
127 G4double ppp = (e+p1+p2)*(e+p1-p2)*(e-p1+p2)*(e-p1-p2)/(4.0*e*e);
128 if (ppp>0) return std::sqrt(ppp);
129 else return -1.;
130}

Referenced by ManyBodyDecayIt(), and TwoBodyDecayIt().

◆ SetParentMass()

void G4GeneralPhaseSpaceDecay::SetParentMass ( const G4double  aParentMass)
inline

Definition at line 112 of file G4GeneralPhaseSpaceDecay.hh.

113{
114 parentmass = aParentMass;
115}

◆ ThreeBodyDecayIt()

G4DecayProducts * G4GeneralPhaseSpaceDecay::ThreeBodyDecayIt ( )
protected

Definition at line 259 of file G4GeneralPhaseSpaceDecay.cc.

261{
262 if (GetVerboseLevel()>1) G4cout << "G4GeneralPhaseSpaceDecay::ThreeBodyDecayIt()"<<G4endl;
263
264 //daughters'mass
265 G4double daughtermass[3];
266 G4double sumofdaughtermass = 0.0;
267 for (G4int index=0; index<3; index++)
268 {
269 if ( theDaughterMasses )
270 {
271 daughtermass[index]= *(theDaughterMasses+index);
272 } else {
273 daughtermass[index] = G4MT_daughters[index]->GetPDGMass();
274 }
275 sumofdaughtermass += daughtermass[index];
276 }
277
278 //create parent G4DynamicParticle at rest
279 G4ParticleMomentum dummy;
280 G4DynamicParticle * parentparticle = new G4DynamicParticle( G4MT_parent, dummy, 0.0);
281
282 //create G4Decayproducts
283 G4DecayProducts *products = new G4DecayProducts(*parentparticle);
284 delete parentparticle;
285
286 //calculate daughter momentum
287 // Generate two
288 G4double rd1, rd2, rd;
289 G4double daughtermomentum[3];
290 G4double momentummax=0.0, momentumsum = 0.0;
291 G4double energy;
292 const G4int maxNumberOfLoops = 10000;
293 G4int loopCounter = 0;
294
295 do
296 {
297 rd1 = G4UniformRand();
298 rd2 = G4UniformRand();
299 if (rd2 > rd1)
300 {
301 rd = rd1;
302 rd1 = rd2;
303 rd2 = rd;
304 }
305 momentummax = 0.0;
306 momentumsum = 0.0;
307 // daughter 0
308
309 energy = rd2*(parentmass - sumofdaughtermass);
310 daughtermomentum[0] = std::sqrt(energy*energy + 2.0*energy* daughtermass[0]);
311 if ( daughtermomentum[0] >momentummax )momentummax = daughtermomentum[0];
312 momentumsum += daughtermomentum[0];
313
314 // daughter 1
315 energy = (1.-rd1)*(parentmass - sumofdaughtermass);
316 daughtermomentum[1] = std::sqrt(energy*energy + 2.0*energy* daughtermass[1]);
317 if ( daughtermomentum[1] >momentummax )momentummax = daughtermomentum[1];
318 momentumsum += daughtermomentum[1];
319
320 // daughter 2
321 energy = (rd1-rd2)*(parentmass - sumofdaughtermass);
322 daughtermomentum[2] = std::sqrt(energy*energy + 2.0*energy* daughtermass[2]);
323 if ( daughtermomentum[2] >momentummax )momentummax = daughtermomentum[2];
324 momentumsum += daughtermomentum[2];
325 } while ( ( momentummax > momentumsum - momentummax ) && /* Loop checking, 02.11.2015, A.Ribon */
326 ++loopCounter < maxNumberOfLoops );
327 if ( loopCounter >= maxNumberOfLoops ) {
328 G4ExceptionDescription ed;
329 ed << " Failed sampling after maxNumberOfLoops attempts : forced exit" << G4endl;
330 G4Exception( " G4GeneralPhaseSpaceDecay::ThreeBodyDecayIt ", "HAD_PHASESPACE_001", FatalException, ed );
331 }
332
333 // output message
334 if (GetVerboseLevel()>1) {
335 G4cout << " daughter 0:" << daughtermomentum[0]/GeV << "[GeV/c]" <<G4endl;
336 G4cout << " daughter 1:" << daughtermomentum[1]/GeV << "[GeV/c]" <<G4endl;
337 G4cout << " daughter 2:" << daughtermomentum[2]/GeV << "[GeV/c]" <<G4endl;
338 G4cout << " momentum sum:" << momentumsum/GeV << "[GeV/c]" <<G4endl;
339 }
340
341 //create daughter G4DynamicParticle
342 G4double costheta, sintheta, phi, sinphi, cosphi;
343 G4double costhetan, sinthetan, phin, sinphin, cosphin;
344 costheta = 2.*G4UniformRand()-1.0;
345 sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
346 phi = twopi*G4UniformRand()*rad;
347 sinphi = std::sin(phi);
348 cosphi = std::cos(phi);
349 G4ParticleMomentum direction0(sintheta*cosphi,sintheta*sinphi,costheta);
350 G4double Etotal=std::sqrt( daughtermass[0]*daughtermass[0] + daughtermomentum[0]*daughtermomentum[0]);
351 G4DynamicParticle * daughterparticle
352 = new G4DynamicParticle( G4MT_daughters[0], Etotal, direction0*daughtermomentum[0]);
353 products->PushProducts(daughterparticle);
354
355 costhetan = (daughtermomentum[1]*daughtermomentum[1]-daughtermomentum[2]*daughtermomentum[2]-daughtermomentum[0]*daughtermomentum[0])/(2.0*daughtermomentum[2]*daughtermomentum[0]);
356 sinthetan = std::sqrt((1.0-costhetan)*(1.0+costhetan));
357 phin = twopi*G4UniformRand()*rad;
358 sinphin = std::sin(phin);
359 cosphin = std::cos(phin);
360 G4ParticleMomentum direction2;
361 direction2.setX( sinthetan*cosphin*costheta*cosphi - sinthetan*sinphin*sinphi + costhetan*sintheta*cosphi);
362 direction2.setY( sinthetan*cosphin*costheta*sinphi + sinthetan*sinphin*cosphi + costhetan*sintheta*sinphi);
363 direction2.setZ( -sinthetan*cosphin*sintheta + costhetan*costheta);
364 Etotal=std::sqrt( daughtermass[2]*daughtermass[2] + daughtermomentum[2]*daughtermomentum[2]/direction2.mag2());
365 daughterparticle = new G4DynamicParticle( G4MT_daughters[2],Etotal, direction2*(daughtermomentum[2]/direction2.mag()));
366 products->PushProducts(daughterparticle);
367 G4ThreeVector mom=(direction0*daughtermomentum[0] + direction2*(daughtermomentum[2]/direction2.mag()))*(-1.0);
368 Etotal= std::sqrt( daughtermass[1]*daughtermass[1] + mom.mag2() );
369 daughterparticle =
370 new G4DynamicParticle(G4MT_daughters[1], Etotal, mom);
371 products->PushProducts(daughterparticle);
372
373 if (GetVerboseLevel()>1) {
374 G4cout << "G4GeneralPhaseSpaceDecay::ThreeBodyDecayIt ";
375 G4cout << " create decay products in rest frame " <<G4endl;
376 products->DumpInfo();
377 }
378 return products;
379}
FatalException
@ FatalException
Definition: G4ExceptionSeverity.hh:61
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.cc:35
G4ExceptionDescription
std::ostringstream G4ExceptionDescription
Definition: G4Exception.hh:40
CLHEP::Hep3Vector::mag2
double mag2() const
CLHEP::Hep3Vector::mag
double mag() const
G4INCL::KinematicsUtils::energy
G4double energy(const ThreeVector &p, const G4double m)
Definition: G4INCLKinematicsUtils.cc:158

Referenced by DecayIt().

◆ TwoBodyDecayIt()

G4DecayProducts * G4GeneralPhaseSpaceDecay::TwoBodyDecayIt ( )
protected

Definition at line 209 of file G4GeneralPhaseSpaceDecay.cc.

210{
211 if (GetVerboseLevel()>1) G4cout << "G4GeneralPhaseSpaceDecay::TwoBodyDecayIt()"<<G4endl;
212
213 //daughters'mass
214 G4double daughtermass[2];
215 G4double daughtermomentum;
216 if ( theDaughterMasses )
217 {
218 daughtermass[0]= *(theDaughterMasses);
219 daughtermass[1] = *(theDaughterMasses+1);
220 } else {
221 daughtermass[0] = G4MT_daughters[0]->GetPDGMass();
222 daughtermass[1] = G4MT_daughters[1]->GetPDGMass();
223 }
224
225// G4double sumofdaughtermass = daughtermass[0] + daughtermass[1];
226
227 //create parent G4DynamicParticle at rest
228 G4ParticleMomentum dummy;
229 G4DynamicParticle * parentparticle = new G4DynamicParticle( G4MT_parent, dummy, 0.0);
230
231 //create G4Decayproducts @@GF why dummy parentparticle?
232 G4DecayProducts *products = new G4DecayProducts(*parentparticle);
233 delete parentparticle;
234
235 //calculate daughter momentum
236 daughtermomentum = Pmx(parentmass,daughtermass[0],daughtermass[1]);
237 G4double costheta = 2.*G4UniformRand()-1.0;
238 G4double sintheta = std::sqrt((1.0 - costheta)*(1.0 + costheta));
239 G4double phi = twopi*G4UniformRand()*rad;
240 G4ParticleMomentum direction(sintheta*std::cos(phi),sintheta*std::sin(phi),costheta);
241
242 //create daughter G4DynamicParticle
243 G4double Etotal= std::sqrt(daughtermass[0]*daughtermass[0] + daughtermomentum*daughtermomentum);
244 G4DynamicParticle * daughterparticle = new G4DynamicParticle( G4MT_daughters[0],Etotal, direction*daughtermomentum);
245 products->PushProducts(daughterparticle);
246 Etotal= std::sqrt(daughtermass[1]*daughtermass[1] + daughtermomentum*daughtermomentum);
247 daughterparticle = new G4DynamicParticle( G4MT_daughters[1],Etotal, direction*(-1.0*daughtermomentum));
248 products->PushProducts(daughterparticle);
249
250 if (GetVerboseLevel()>1)
251 {
252 G4cout << "G4GeneralPhaseSpaceDecay::TwoBodyDecayIt ";
253 G4cout << " create decay products in rest frame " <<G4endl;
254 products->DumpInfo();
255 }
256 return products;
257}

Referenced by DecayIt().

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