G4GeneralPhaseSpaceDecay Class Reference

#include <G4GeneralPhaseSpaceDecay.hh>

Inheritance diagram for G4GeneralPhaseSpaceDecay:

G4VDecayChannel G4NuclearDecayChannel G4AlphaDecayChannel G4BetaMinusDecayChannel G4BetaPlusDecayChannel G4ITDecayChannel G4KshellECDecayChannel G4LshellECDecayChannel G4MshellECDecayChannel

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)

Static Public Member Functions

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

Protected Member Functions

G4DecayProducts * OneBodyDecayIt ()
G4DecayProducts * TwoBodyDecayIt ()
G4DecayProducts * ThreeBodyDecayIt ()
G4DecayProducts * ManyBodyDecayIt ()

Detailed Description

Definition at line 45 of file G4GeneralPhaseSpaceDecay.hh.


Constructor & Destructor Documentation

G4GeneralPhaseSpaceDecay::G4GeneralPhaseSpaceDecay ( G4int  Verbose = 1  ) 

Definition at line 50 of file G4GeneralPhaseSpaceDecay.cc.

References G4cout, G4endl, and G4VDecayChannel::GetVerboseLevel().

00050                                                                 : 
00051                           G4VDecayChannel("Phase Space", Verbose),
00052                           parentmass(0.), theDaughterMasses(0)
00053 {
00054   if (GetVerboseLevel()>1) G4cout << "G4GeneralPhaseSpaceDecay:: constructor " << G4endl;
00055 }

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.

References G4cout, G4endl, G4ParticleDefinition::GetPDGMass(), G4VDecayChannel::GetVerboseLevel(), and G4VDecayChannel::parent.

00062                                                                                      :
00063                                    G4VDecayChannel("Phase Space",
00064                                                    theParentName,theBR,
00065                                                    theNumberOfDaughters,
00066                                                    theDaughterName1,
00067                                                    theDaughterName2,
00068                                                    theDaughterName3),
00069                                    theDaughterMasses(0)
00070 {
00071   if (GetVerboseLevel()>1) G4cout << "G4GeneralPhaseSpaceDecay:: constructor " << G4endl;
00072   
00073   //   Set the parent particle (resonance) mass to the (default) PDG vale
00074   if (parent != NULL)
00075   {
00076       parentmass = parent->GetPDGMass();
00077   } else {
00078           parentmass=0.;
00079   }
00080 
00081 }

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.

References G4cout, G4endl, and G4VDecayChannel::GetVerboseLevel().

00089                                                                                      :
00090                                    G4VDecayChannel("Phase Space",
00091                                                    theParentName,theBR,
00092                                                    theNumberOfDaughters,
00093                                                    theDaughterName1,
00094                                                    theDaughterName2,
00095                                                    theDaughterName3),
00096                                    parentmass(theParentMass),
00097                                    theDaughterMasses(0)
00098 {
00099   if (GetVerboseLevel()>1) G4cout << "G4GeneralPhaseSpaceDecay:: constructor " << G4endl;
00100 }

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.

References G4cout, G4endl, and G4VDecayChannel::GetVerboseLevel().

00109                                                                            :
00110                                    G4VDecayChannel("Phase Space",
00111                                                    theParentName,theBR,
00112                                                    theNumberOfDaughters,
00113                                                    theDaughterName1,
00114                                                    theDaughterName2,
00115                                                    theDaughterName3),
00116                                    parentmass(theParentMass),
00117                                    theDaughterMasses(masses)
00118 {
00119   if (GetVerboseLevel()>1) G4cout << "G4GeneralPhaseSpaceDecay:: constructor " << G4endl;
00120 }

G4GeneralPhaseSpaceDecay::~G4GeneralPhaseSpaceDecay (  )  [virtual]

Definition at line 122 of file G4GeneralPhaseSpaceDecay.cc.

00123 {
00124 }


Member Function Documentation

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

Implements G4VDecayChannel.

Reimplemented in G4NuclearDecayChannel.

Definition at line 126 of file G4GeneralPhaseSpaceDecay.cc.

References G4VDecayChannel::daughters, G4VDecayChannel::DumpInfo(), G4VDecayChannel::FillDaughters(), G4VDecayChannel::FillParent(), G4cout, G4endl, G4VDecayChannel::GetVerboseLevel(), ManyBodyDecayIt(), G4VDecayChannel::numberOfDaughters, OneBodyDecayIt(), G4VDecayChannel::parent, G4VDecayChannel::parent_name, ThreeBodyDecayIt(), and TwoBodyDecayIt().

Referenced by G4KineticTrack::Decay().

00127 {
00128   if (GetVerboseLevel()>1) G4cout << "G4GeneralPhaseSpaceDecay::DecayIt ";
00129   G4DecayProducts * products = NULL;
00130  
00131   if (parent == NULL) FillParent();  
00132   if (daughters == NULL) FillDaughters();
00133   
00134   switch (numberOfDaughters){
00135   case 0:
00136     if (GetVerboseLevel()>0) {
00137       G4cout << "G4GeneralPhaseSpaceDecay::DecayIt ";
00138       G4cout << " daughters not defined " <<G4endl;
00139     }
00140     break;
00141   case 1:
00142     products =  OneBodyDecayIt();
00143     break;
00144   case 2:
00145     products =  TwoBodyDecayIt();
00146     break;
00147   case 3:
00148     products =  ThreeBodyDecayIt();
00149     break;
00150   default:
00151     products =  ManyBodyDecayIt();
00152     break;
00153   }
00154   if ((products == NULL) && (GetVerboseLevel()>0)) {
00155     G4cout << "G4GeneralPhaseSpaceDecay::DecayIt ";
00156     G4cout << *parent_name << " can not decay " << G4endl;
00157     DumpInfo();
00158   }
00159   return products;
00160 }

G4double G4GeneralPhaseSpaceDecay::GetParentMass (  )  const [inline]

Reimplemented from G4VDecayChannel.

Definition at line 98 of file G4GeneralPhaseSpaceDecay.hh.

00099 {
00100   return parentmass;
00101 }

G4DecayProducts * G4GeneralPhaseSpaceDecay::ManyBodyDecayIt (  )  [protected]

Definition at line 353 of file G4GeneralPhaseSpaceDecay.cc.

References G4VDecayChannel::daughters, G4VDecayChannel::daughters_name, G4DecayProducts::DumpInfo(), G4cout, G4endl, G4UniformRand, G4DynamicParticle::Get4Momentum(), G4VDecayChannel::GetVerboseLevel(), G4VDecayChannel::numberOfDaughters, G4VDecayChannel::parent, G4VDecayChannel::parent_name, Pmx(), G4DecayProducts::PushProducts(), G4DynamicParticle::Set4Momentum(), and G4InuclParticleNames::sm.

Referenced by DecayIt().

00361                     : 19/Apr/1995
00362 //
00363 {
00364   //return value
00365   G4DecayProducts *products;
00366 
00367   if (GetVerboseLevel()>1) G4cout << "G4GeneralPhaseSpaceDecay::ManyBodyDecayIt()"<<G4endl;
00368 
00369   //daughters'mass
00370   G4double *daughtermass = new G4double[numberOfDaughters]; 
00371   G4double sumofdaughtermass = 0.0;
00372   for (G4int index=0; index<numberOfDaughters; index++){
00373     daughtermass[index] = daughters[index]->GetPDGMass();
00374     sumofdaughtermass += daughtermass[index];
00375   }
00376   
00377   //Calculate daughter momentum
00378   G4double *daughtermomentum = new G4double[numberOfDaughters];
00379   G4ParticleMomentum direction;  
00380   G4DynamicParticle **daughterparticle;
00381   G4double *sm = new G4double[numberOfDaughters];
00382   G4double tmas;
00383   G4double weight = 1.0;
00384   G4int    numberOfTry = 0; 
00385   G4int index1;
00386 
00387   do {
00388     //Generate rundom number in descending order 
00389     G4double temp;
00390     G4double *rd = new G4double[numberOfDaughters];
00391     rd[0] = 1.0;
00392     for(index1 =1; index1 < numberOfDaughters -1; index1++)
00393       rd[index1] = G4UniformRand(); 
00394     rd[ numberOfDaughters -1] = 0.0;
00395     for(index1 =1; index1 < numberOfDaughters -1; index1++) {
00396       for(G4int index2 = index1+1; index2 < numberOfDaughters; index2++) {
00397         if (rd[index1] < rd[index2]){
00398           temp         = rd[index1];
00399           rd[index1]   = rd[index2];
00400           rd[index2]   = temp;
00401         }
00402       }
00403     }
00404     
00405     //calcurate virtual mass 
00406     tmas = parentmass -  sumofdaughtermass;
00407     temp =  sumofdaughtermass; 
00408     for(index1 =0; index1 < numberOfDaughters; index1++) {
00409       sm[index1] = rd[index1]*tmas + temp;
00410       temp -= daughtermass[index1];
00411       if (GetVerboseLevel()>1) {
00412         G4cout << index1 << "  rundom number:" << rd[index1]; 
00413         G4cout << "   virtual mass:" << sm[index1]/GeV << "[GeV/c/c]" <<G4endl; 
00414       }
00415     }
00416     delete [] rd;
00417 
00418     //Calculate daughter momentum
00419     weight = 1.0;
00420     index1 =numberOfDaughters-1;
00421     daughtermomentum[index1]= Pmx( sm[index1-1],daughtermass[index1-1],sm[index1]);
00422     if (GetVerboseLevel()>1) {
00423       G4cout << "     daughter " << index1 << ":" << *daughters_name[index1];
00424       G4cout << " momentum:" << daughtermomentum[index1]/GeV << "[GeV/c]" <<G4endl;
00425     }
00426     for(index1 =numberOfDaughters-2; index1>=0; index1--) {
00427       // calculate 
00428       daughtermomentum[index1]= Pmx( sm[index1],daughtermass[index1], sm[index1 +1]);
00429       if(daughtermomentum[index1] < 0.0) {
00430         // !!! illegal momentum !!!
00431         if (GetVerboseLevel()>0) {
00432           G4cout << "G4GeneralPhaseSpaceDecay::ManyBodyDecayIt ";
00433           G4cout << "     can not calculate daughter momentum " <<G4endl;
00434           G4cout << "     parent:" << *parent_name;
00435           G4cout << " mass:" << parentmass/GeV << "[GeV/c/c]" <<G4endl;
00436           G4cout << "     daughter " << index1 << ":" << *daughters_name[index1];
00437           G4cout << " mass:" << daughtermass[index1]/GeV << "[GeV/c/c]" ;
00438           G4cout << " mass:" << daughtermomentum[index1]/GeV << "[GeV/c]" <<G4endl;
00439         }
00440         delete [] sm;
00441         delete [] daughtermass;
00442         delete [] daughtermomentum;
00443         return NULL;   // Error detection
00444 
00445       } else {
00446         // calculate weight of this events
00447         weight *=  daughtermomentum[index1]/sm[index1];
00448         if (GetVerboseLevel()>1) {
00449           G4cout << "     daughter " << index1 << ":" << *daughters_name[index1];
00450           G4cout << " momentum:" << daughtermomentum[index1]/GeV << "[GeV/c]" <<G4endl;
00451         }
00452       }
00453     }
00454     if (GetVerboseLevel()>1) {
00455       G4cout << "    weight: " << weight <<G4endl;
00456     }
00457     
00458     // exit if number of Try exceeds 100
00459     if (numberOfTry++ >100) {
00460       if (GetVerboseLevel()>0) {
00461         G4cout << "G4GeneralPhaseSpaceDecay::ManyBodyDecayIt: ";
00462         G4cout << " can not determine Decay Kinematics " << G4endl;
00463       }
00464       delete [] sm;
00465       delete [] daughtermass;
00466       delete [] daughtermomentum;
00467       return NULL;  // Error detection
00468     }
00469   } while ( weight > G4UniformRand());
00470   if (GetVerboseLevel()>1) {
00471       G4cout << "Start calulation of daughters momentum vector "<<G4endl;
00472   }
00473   
00474   G4double costheta, sintheta, phi;
00475   G4double beta;
00476   daughterparticle = new G4DynamicParticle*[numberOfDaughters];
00477 
00478   index1 = numberOfDaughters -2;
00479   costheta = 2.*G4UniformRand()-1.0;
00480   sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
00481   phi  = twopi*G4UniformRand()*rad;
00482   direction.setZ(costheta);
00483   direction.setY(sintheta*std::sin(phi));
00484   direction.setX(sintheta*std::cos(phi));
00485   daughterparticle[index1] = new G4DynamicParticle( daughters[index1], direction*daughtermomentum[index1] );
00486   daughterparticle[index1+1] = new G4DynamicParticle( daughters[index1+1], direction*(-1.0*daughtermomentum[index1]) );
00487 
00488   for (index1 = numberOfDaughters -3;  index1 >= 0; index1--) {
00489     //calculate momentum direction
00490     costheta = 2.*G4UniformRand()-1.0;
00491     sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
00492     phi  = twopi*G4UniformRand()*rad;
00493     direction.setZ(costheta);
00494     direction.setY(sintheta*std::sin(phi));
00495     direction.setX(sintheta*std::cos(phi));
00496 
00497     // boost already created particles 
00498     beta = daughtermomentum[index1];
00499     beta /= std::sqrt( daughtermomentum[index1]*daughtermomentum[index1] + sm[index1+1]*sm[index1+1] );
00500     for (G4int index2 = index1+1; index2<numberOfDaughters; index2++) {
00501       G4LorentzVector p4;
00502       // make G4LorentzVector for secondaries
00503       p4 = daughterparticle[index2]->Get4Momentum();
00504 
00505       // boost secondaries to  new frame 
00506       p4.boost( direction.x()*beta, direction.y()*beta, direction.z()*beta);
00507 
00508       // change energy/momentum
00509       daughterparticle[index2]->Set4Momentum(p4);
00510     }
00511     //create daughter G4DynamicParticle 
00512     daughterparticle[index1]= new G4DynamicParticle( daughters[index1], direction*(-1.0*daughtermomentum[index1]));
00513   }
00514 
00515   //create G4Decayproducts
00516   G4DynamicParticle *parentparticle; 
00517   direction.setX(1.0);  direction.setY(0.0); direction.setZ(0.0);
00518   parentparticle = new G4DynamicParticle( parent, direction, 0.0);
00519   products = new G4DecayProducts(*parentparticle);
00520   delete parentparticle;
00521   for (index1 = 0; index1<numberOfDaughters; index1++) {
00522     products->PushProducts(daughterparticle[index1]);
00523   }
00524   if (GetVerboseLevel()>1) { 
00525     G4cout << "G4GeneralPhaseSpaceDecay::ManyBodyDecayIt ";
00526     G4cout << "  create decay products in rest frame " << G4endl;
00527     products->DumpInfo();
00528   }
00529 
00530   delete [] daughterparticle;
00531   delete [] daughtermomentum;
00532   delete [] daughtermass;
00533   delete [] sm;
00534   
00535   return products;
00536 }

G4DecayProducts * G4GeneralPhaseSpaceDecay::OneBodyDecayIt (  )  [protected]

Definition at line 162 of file G4GeneralPhaseSpaceDecay.cc.

References G4VDecayChannel::daughters, G4DecayProducts::DumpInfo(), G4cout, G4endl, G4VDecayChannel::GetVerboseLevel(), G4VDecayChannel::parent, and G4DecayProducts::PushProducts().

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

00163 {
00164   if (GetVerboseLevel()>1) G4cout << "G4GeneralPhaseSpaceDecay::OneBodyDecayIt()"<<G4endl;
00165 
00166 //  G4double daughtermass = daughters[0]->GetPDGMass();
00167 
00168   //create parent G4DynamicParticle at rest
00169   G4ParticleMomentum dummy;
00170   G4DynamicParticle * parentparticle = new G4DynamicParticle(parent, dummy, 0.0);
00171 
00172   //create G4Decayproducts
00173   G4DecayProducts *products = new G4DecayProducts(*parentparticle);
00174   delete parentparticle;
00175 
00176   //create daughter G4DynamicParticle at rest
00177   G4DynamicParticle * daughterparticle = new G4DynamicParticle(daughters[0], dummy, 0.0);
00178   products->PushProducts(daughterparticle);
00179 
00180   if (GetVerboseLevel()>1) 
00181     {
00182      G4cout << "G4GeneralPhaseSpaceDecay::OneBodyDecayIt ";
00183      G4cout << "  create decay products in rest frame " <<G4endl;
00184      products->DumpInfo();
00185     }
00186   return products;
00187 }

G4double G4GeneralPhaseSpaceDecay::Pmx ( G4double  e,
G4double  p1,
G4double  p2 
) [inline, static]

Definition at line 111 of file G4GeneralPhaseSpaceDecay.hh.

Referenced by ManyBodyDecayIt(), and TwoBodyDecayIt().

00112 {
00113    // calculate momentum of daughter particles in two-body decay
00114    if (e-p1-p2 < 0 )
00115    {  
00116         G4HadronicException(__FILE__, __LINE__, "G4GeneralPhaseSpaceDecay::Pmx energy in cms > mass1+mass2");   
00117    }
00118    G4double ppp = (e+p1+p2)*(e+p1-p2)*(e-p1+p2)*(e-p1-p2)/(4.0*e*e);
00119    if (ppp>0) return std::sqrt(ppp);
00120    else       return -1.;
00121 }

void G4GeneralPhaseSpaceDecay::SetParentMass ( const G4double  aParentMass  )  [inline]

Definition at line 103 of file G4GeneralPhaseSpaceDecay.hh.

Referenced by G4NuclearDecayChannel::DecayIt().

00104 {
00105   parentmass = aParentMass;
00106 }

G4DecayProducts * G4GeneralPhaseSpaceDecay::ThreeBodyDecayIt (  )  [protected]

Definition at line 239 of file G4GeneralPhaseSpaceDecay.cc.

References G4VDecayChannel::daughters, G4DecayProducts::DumpInfo(), G4cout, G4endl, G4UniformRand, G4ParticleDefinition::GetPDGMass(), G4VDecayChannel::GetVerboseLevel(), G4VDecayChannel::parent, and G4DecayProducts::PushProducts().

Referenced by DecayIt().

00241 {
00242   if (GetVerboseLevel()>1) G4cout << "G4GeneralPhaseSpaceDecay::ThreeBodyDecayIt()"<<G4endl;
00243 
00244   //daughters'mass
00245   G4double daughtermass[3]; 
00246   G4double sumofdaughtermass = 0.0;
00247   for (G4int index=0; index<3; index++)
00248     {
00249      if ( theDaughterMasses )
00250      { 
00251          daughtermass[index]= *(theDaughterMasses+index);
00252      } else {   
00253          daughtermass[index] = daughters[index]->GetPDGMass();
00254      }   
00255      sumofdaughtermass += daughtermass[index];
00256     }
00257   
00258   //create parent G4DynamicParticle at rest
00259   G4ParticleMomentum dummy;
00260   G4DynamicParticle * parentparticle = new G4DynamicParticle( parent, dummy, 0.0);
00261 
00262   //create G4Decayproducts
00263   G4DecayProducts *products = new G4DecayProducts(*parentparticle);
00264   delete parentparticle;
00265 
00266   //calculate daughter momentum
00267   //  Generate two 
00268   G4double rd1, rd2, rd;
00269   G4double daughtermomentum[3];
00270   G4double momentummax=0.0, momentumsum = 0.0;
00271   G4double energy;
00272 
00273   do 
00274     {
00275      rd1 = G4UniformRand();
00276      rd2 = G4UniformRand();
00277      if (rd2 > rd1) 
00278        {
00279         rd  = rd1;
00280         rd1 = rd2;
00281         rd2 = rd;
00282        } 
00283      momentummax = 0.0;
00284      momentumsum = 0.0;
00285      // daughter 0
00286 
00287      energy = rd2*(parentmass - sumofdaughtermass);
00288      daughtermomentum[0] = std::sqrt(energy*energy + 2.0*energy* daughtermass[0]);
00289      if ( daughtermomentum[0] >momentummax )momentummax =  daughtermomentum[0];
00290      momentumsum  +=  daughtermomentum[0];
00291 
00292      // daughter 1
00293      energy = (1.-rd1)*(parentmass - sumofdaughtermass);
00294      daughtermomentum[1] = std::sqrt(energy*energy + 2.0*energy* daughtermass[1]);
00295      if ( daughtermomentum[1] >momentummax )momentummax =  daughtermomentum[1];
00296      momentumsum  +=  daughtermomentum[1];
00297 
00298      // daughter 2
00299      energy = (rd1-rd2)*(parentmass - sumofdaughtermass);
00300      daughtermomentum[2] = std::sqrt(energy*energy + 2.0*energy* daughtermass[2]);
00301      if ( daughtermomentum[2] >momentummax )momentummax =  daughtermomentum[2];
00302      momentumsum  +=  daughtermomentum[2];
00303     } while (momentummax >  momentumsum - momentummax );
00304 
00305   // output message
00306   if (GetVerboseLevel()>1) {
00307     G4cout << "     daughter 0:" << daughtermomentum[0]/GeV << "[GeV/c]" <<G4endl;
00308     G4cout << "     daughter 1:" << daughtermomentum[1]/GeV << "[GeV/c]" <<G4endl;
00309     G4cout << "     daughter 2:" << daughtermomentum[2]/GeV << "[GeV/c]" <<G4endl;
00310     G4cout << "   momentum sum:" << momentumsum/GeV << "[GeV/c]" <<G4endl;
00311   }
00312 
00313   //create daughter G4DynamicParticle 
00314   G4double costheta, sintheta, phi, sinphi, cosphi; 
00315   G4double costhetan, sinthetan, phin, sinphin, cosphin; 
00316   costheta = 2.*G4UniformRand()-1.0;
00317   sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
00318   phi  = twopi*G4UniformRand()*rad;
00319   sinphi = std::sin(phi);
00320   cosphi = std::cos(phi);
00321   G4ParticleMomentum direction0(sintheta*cosphi,sintheta*sinphi,costheta);
00322   G4double Etotal=std::sqrt( daughtermass[0]*daughtermass[0] + daughtermomentum[0]*daughtermomentum[0]);
00323   G4DynamicParticle * daughterparticle 
00324          = new G4DynamicParticle( daughters[0], Etotal, direction0*daughtermomentum[0]);
00325   products->PushProducts(daughterparticle);
00326 
00327   costhetan = (daughtermomentum[1]*daughtermomentum[1]-daughtermomentum[2]*daughtermomentum[2]-daughtermomentum[0]*daughtermomentum[0])/(2.0*daughtermomentum[2]*daughtermomentum[0]);
00328   sinthetan = std::sqrt((1.0-costhetan)*(1.0+costhetan));
00329   phin  = twopi*G4UniformRand()*rad;
00330   sinphin = std::sin(phin);
00331   cosphin = std::cos(phin);
00332   G4ParticleMomentum direction2;
00333   direction2.setX( sinthetan*cosphin*costheta*cosphi - sinthetan*sinphin*sinphi + costhetan*sintheta*cosphi); 
00334   direction2.setY( sinthetan*cosphin*costheta*sinphi + sinthetan*sinphin*cosphi + costhetan*sintheta*sinphi); 
00335   direction2.setZ( -sinthetan*cosphin*sintheta + costhetan*costheta);
00336   Etotal=std::sqrt( daughtermass[2]*daughtermass[2] + daughtermomentum[2]*daughtermomentum[2]/direction2.mag2());
00337   daughterparticle = new G4DynamicParticle( daughters[2],Etotal, direction2*(daughtermomentum[2]/direction2.mag()));
00338   products->PushProducts(daughterparticle);
00339   G4ThreeVector mom=(direction0*daughtermomentum[0] + direction2*(daughtermomentum[2]/direction2.mag()))*(-1.0);
00340   Etotal= std::sqrt( daughtermass[1]*daughtermass[1] + mom.mag2() );
00341   daughterparticle = 
00342        new G4DynamicParticle(daughters[1], Etotal, mom);
00343   products->PushProducts(daughterparticle);
00344 
00345   if (GetVerboseLevel()>1) {
00346      G4cout << "G4GeneralPhaseSpaceDecay::ThreeBodyDecayIt ";
00347      G4cout << "  create decay products in rest frame " <<G4endl;
00348      products->DumpInfo();
00349   }
00350   return products;
00351 }

G4DecayProducts * G4GeneralPhaseSpaceDecay::TwoBodyDecayIt (  )  [protected]

Definition at line 189 of file G4GeneralPhaseSpaceDecay.cc.

References G4VDecayChannel::daughters, G4DecayProducts::DumpInfo(), G4cout, G4endl, G4UniformRand, G4ParticleDefinition::GetPDGMass(), G4VDecayChannel::GetVerboseLevel(), G4VDecayChannel::parent, Pmx(), and G4DecayProducts::PushProducts().

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

00190 {
00191   if (GetVerboseLevel()>1) G4cout << "G4GeneralPhaseSpaceDecay::TwoBodyDecayIt()"<<G4endl;
00192   
00193   //daughters'mass
00194   G4double daughtermass[2]; 
00195   G4double daughtermomentum;
00196   if ( theDaughterMasses )
00197   { 
00198      daughtermass[0]= *(theDaughterMasses);
00199      daughtermass[1] = *(theDaughterMasses+1);
00200   } else {   
00201      daughtermass[0] = daughters[0]->GetPDGMass();
00202      daughtermass[1] = daughters[1]->GetPDGMass();
00203   }
00204   
00205 //  G4double sumofdaughtermass =  daughtermass[0] + daughtermass[1];
00206 
00207   //create parent G4DynamicParticle at rest
00208   G4ParticleMomentum dummy;
00209   G4DynamicParticle * parentparticle = new G4DynamicParticle( parent, dummy, 0.0);
00210 
00211   //create G4Decayproducts  @@GF why dummy parentparticle?
00212   G4DecayProducts *products = new G4DecayProducts(*parentparticle);
00213   delete parentparticle;
00214 
00215   //calculate daughter momentum
00216   daughtermomentum = Pmx(parentmass,daughtermass[0],daughtermass[1]);
00217   G4double costheta = 2.*G4UniformRand()-1.0;
00218   G4double sintheta = std::sqrt((1.0 - costheta)*(1.0 + costheta));
00219   G4double phi  = twopi*G4UniformRand()*rad;
00220   G4ParticleMomentum direction(sintheta*std::cos(phi),sintheta*std::sin(phi),costheta);
00221 
00222   //create daughter G4DynamicParticle
00223   G4double Etotal= std::sqrt(daughtermass[0]*daughtermass[0] + daughtermomentum*daughtermomentum); 
00224   G4DynamicParticle * daughterparticle = new G4DynamicParticle( daughters[0],Etotal, direction*daughtermomentum);
00225   products->PushProducts(daughterparticle);
00226   Etotal= std::sqrt(daughtermass[1]*daughtermass[1] + daughtermomentum*daughtermomentum);
00227   daughterparticle = new G4DynamicParticle( daughters[1],Etotal, direction*(-1.0*daughtermomentum));
00228   products->PushProducts(daughterparticle);
00229 
00230   if (GetVerboseLevel()>1) 
00231     {
00232      G4cout << "G4GeneralPhaseSpaceDecay::TwoBodyDecayIt ";
00233      G4cout << "  create decay products in rest frame " <<G4endl;
00234      products->DumpInfo();
00235     }
00236   return products;
00237 }


The documentation for this class was generated from the following files:
Generated on Mon May 27 17:52:04 2013 for Geant4 by  doxygen 1.4.7