Geant4  10.03.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4GeneralPhaseSpaceDecay Class Reference

#include <G4GeneralPhaseSpaceDecay.hh>

Inheritance diagram for G4GeneralPhaseSpaceDecay:
Collaboration 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 G4DecayProductsDecayIt (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
 
const G4StringGetKinematicsName () const
 
G4double GetBR () const
 
G4int GetNumberOfDaughters () const
 
G4ParticleDefinitionGetParent ()
 
G4ParticleDefinitionGetDaughter (G4int anIndex)
 
G4int GetAngularMomentum ()
 
const G4StringGetParentName () const
 
const G4StringGetDaughterName (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 G4ThreeVectorGetPolarization () const
 

Static Public Member Functions

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

Protected Member Functions

G4DecayProductsOneBodyDecayIt ()
 
G4DecayProductsTwoBodyDecayIt ()
 
G4DecayProductsThreeBodyDecayIt ()
 
G4DecayProductsManyBodyDecayIt ()
 
- Protected Member Functions inherited from G4VDecayChannel
void ClearDaughtersName ()
 
void CheckAndFillDaughters ()
 
void CheckAndFillParent ()
 
 G4VDecayChannel ()
 
 G4VDecayChannel (const G4VDecayChannel &)
 
G4VDecayChanneloperator= (const G4VDecayChannel &)
 
G4double DynamicalMass (G4double massPDG, G4double width, G4double maxDev=+1.) const
 

Additional Inherited Members

- Protected Attributes inherited from G4VDecayChannel
G4String kinematics_name
 
G4double rbranch
 
G4int numberOfDaughters
 
G4Stringparent_name
 
G4String ** daughters_name
 
G4double rangeMass
 
G4ThreeVector parent_polarization
 
G4ParticleTableparticletable
 
G4int verboseLevel
 
G4ParticleDefinitionG4MT_parent
 
G4ParticleDefinition ** G4MT_daughters
 
G4double G4MT_parent_mass
 
G4doubleG4MT_daughters_mass
 
G4doubleG4MT_daughters_width
 
G4Mutex daughtersMutex
 
G4Mutex parentMutex
 
- Static Protected Attributes inherited from G4VDecayChannel
static const G4String noName = " "
 

Detailed Description

Definition at line 45 of file G4GeneralPhaseSpaceDecay.hh.

Constructor & Destructor Documentation

G4GeneralPhaseSpaceDecay::G4GeneralPhaseSpaceDecay ( G4int  Verbose = 1)

Definition at line 49 of file G4GeneralPhaseSpaceDecay.cc.

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

Here is the call graph for this function:

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

Definition at line 56 of file G4GeneralPhaseSpaceDecay.cc.

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

Here is the call graph for this function:

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

Definition at line 82 of file G4GeneralPhaseSpaceDecay.cc.

88  :
89  G4VDecayChannel("Phase Space",
90  theParentName,theBR,
91  theNumberOfDaughters,
92  theDaughterName1,
93  theDaughterName2,
94  theDaughterName3),
95  parentmass(theParentMass),
96  theDaughterMasses(0)
97 {
98  if (GetVerboseLevel()>1) G4cout << "G4GeneralPhaseSpaceDecay:: constructor " << G4endl;
99 }
G4GLOB_DLL std::ostream G4cout
G4int GetVerboseLevel() const
#define G4endl
Definition: G4ios.hh:61

Here is the call graph for this function:

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

108  :
109  G4VDecayChannel("Phase Space",
110  theParentName,theBR,
111  theNumberOfDaughters,
112  theDaughterName1,
113  theDaughterName2,
114  theDaughterName3),
115  parentmass(theParentMass),
116  theDaughterMasses(masses)
117 {
118  if (GetVerboseLevel()>1) G4cout << "G4GeneralPhaseSpaceDecay:: constructor " << G4endl;
119 }
G4GLOB_DLL std::ostream G4cout
G4int GetVerboseLevel() const
#define G4endl
Definition: G4ios.hh:61

Here is the call graph for this function:

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

129  :
130  G4VDecayChannel("Phase Space",
131  theParentName,theBR,
132  theNumberOfDaughters,
133  theDaughterName1,
134  theDaughterName2,
135  theDaughterName3,
136  theDaughterName4),
137  parentmass(theParentMass),
138  theDaughterMasses(masses)
139 {
140  if (GetVerboseLevel()>1) G4cout << "G4GeneralPhaseSpaceDecay:: constructor " << G4endl;
141 }
G4GLOB_DLL std::ostream G4cout
G4int GetVerboseLevel() const
#define G4endl
Definition: G4ios.hh:61

Here is the call graph for this function:

G4GeneralPhaseSpaceDecay::~G4GeneralPhaseSpaceDecay ( )
virtual

Definition at line 143 of file G4GeneralPhaseSpaceDecay.cc.

144 {
145 }

Member Function Documentation

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

Implements G4VDecayChannel.

Definition at line 147 of file G4GeneralPhaseSpaceDecay.cc.

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

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4GeneralPhaseSpaceDecay::GetParentMass ( ) const
inline

Definition at line 108 of file G4GeneralPhaseSpaceDecay.hh.

109 {
110  return parentmass;
111 }
G4DecayProducts * G4GeneralPhaseSpaceDecay::ManyBodyDecayIt ( )
protected

Definition at line 382 of file G4GeneralPhaseSpaceDecay.cc.

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

Here is the call graph for this function:

Here is the caller graph for this function:

G4DecayProducts * G4GeneralPhaseSpaceDecay::OneBodyDecayIt ( )
protected

Definition at line 183 of file G4GeneralPhaseSpaceDecay.cc.

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

Here is the call graph for this function:

Here is the caller graph for this function:

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

Definition at line 121 of file G4GeneralPhaseSpaceDecay.hh.

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

Here is the caller graph for this function:

void G4GeneralPhaseSpaceDecay::SetParentMass ( const G4double  aParentMass)
inline

Definition at line 113 of file G4GeneralPhaseSpaceDecay.hh.

114 {
115  parentmass = aParentMass;
116 }
G4DecayProducts * G4GeneralPhaseSpaceDecay::ThreeBodyDecayIt ( )
protected

Definition at line 260 of file G4GeneralPhaseSpaceDecay.cc.

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

Here is the call graph for this function:

Here is the caller graph for this function:

G4DecayProducts * G4GeneralPhaseSpaceDecay::TwoBodyDecayIt ( )
protected

Definition at line 210 of file G4GeneralPhaseSpaceDecay.cc.

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

Here is the call graph for this function:

Here is the caller graph for this function:


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