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 G4DecayProductsDecayIt (G4double mass=0.0)

Static Public Member Functions

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

Protected Member Functions

G4DecayProductsOneBodyDecayIt ()
G4DecayProductsTwoBodyDecayIt ()
G4DecayProductsThreeBodyDecayIt ()
G4DecayProductsManyBodyDecayIt ()

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