G4PhaseSpaceDecayChannel.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 //
00027 // $Id$
00028 //
00029 // 
00030 // ------------------------------------------------------------
00031 //      GEANT 4 class header file
00032 //
00033 //      History: first implementation, based on object model of
00034 //      27 July 1996 H.Kurashige
00035 //      20 Oct 1996 H.Kurashige
00036 //      30 May 1997 H.Kurashige
00037 // ------------------------------------------------------------
00038 
00039 #include "G4ParticleDefinition.hh"
00040 #include "G4PhysicalConstants.hh"
00041 #include "G4SystemOfUnits.hh"
00042 #include "G4DecayProducts.hh"
00043 #include "G4VDecayChannel.hh"
00044 #include "G4PhaseSpaceDecayChannel.hh"
00045 #include "Randomize.hh"
00046 #include "G4LorentzVector.hh"
00047 #include "G4LorentzRotation.hh"
00048 
00049 
00050 G4PhaseSpaceDecayChannel::G4PhaseSpaceDecayChannel(G4int Verbose)
00051  :G4VDecayChannel("Phase Space", Verbose)
00052 {
00053 
00054 }
00055 
00056 G4PhaseSpaceDecayChannel::G4PhaseSpaceDecayChannel(
00057                                const G4String& theParentName,
00058                                G4double        theBR,
00059                                G4int           theNumberOfDaughters,
00060                                const G4String& theDaughterName1,
00061                                const G4String& theDaughterName2,
00062                                const G4String& theDaughterName3,
00063                                const G4String& theDaughterName4    )
00064                          :G4VDecayChannel("Phase Space",
00065                                           theParentName,theBR,
00066                                           theNumberOfDaughters,
00067                                           theDaughterName1,
00068                                           theDaughterName2,
00069                                           theDaughterName3,
00070                                           theDaughterName4),
00071                           current_parent_mass(0.0)
00072 {
00073 
00074 }
00075 
00076 G4PhaseSpaceDecayChannel::~G4PhaseSpaceDecayChannel()
00077 {
00078 }
00079 
00080 G4DecayProducts *G4PhaseSpaceDecayChannel::DecayIt(G4double parentMass)
00081 {
00082 #ifdef G4VERBOSE
00083   if (GetVerboseLevel()>1) G4cout << "G4PhaseSpaceDecayChannel::DecayIt ";
00084 #endif
00085   
00086   G4DecayProducts * products = 0;
00087  
00088   if (parent == 0) FillParent();  
00089   if (daughters == 0) FillDaughters();
00090 
00091   if (parentMass >0.0) current_parent_mass = parentMass;
00092   else                 current_parent_mass = parent_mass;
00093 
00094   switch (numberOfDaughters){
00095   case 0:
00096 #ifdef G4VERBOSE
00097     if (GetVerboseLevel()>0) {
00098       G4cout << "G4PhaseSpaceDecayChannel::DecayIt ";
00099       G4cout << " daughters not defined " <<G4endl;
00100     }
00101 #endif
00102     break;
00103   case 1:
00104     products =  OneBodyDecayIt();
00105     break;
00106   case 2:
00107     products =  TwoBodyDecayIt();
00108     break;
00109   case 3:
00110     products =  ThreeBodyDecayIt();
00111     break;
00112   default:
00113     products =  ManyBodyDecayIt();
00114     break;
00115   }
00116 #ifdef G4VERBOSE
00117   if ((products == 0) && (GetVerboseLevel()>0)) {
00118     G4cout << "G4PhaseSpaceDecayChannel::DecayIt ";
00119     G4cout << *parent_name << " can not decay " << G4endl;
00120     DumpInfo();
00121   }
00122 #endif
00123   return products;
00124 }
00125 
00126 G4DecayProducts *G4PhaseSpaceDecayChannel::OneBodyDecayIt()
00127 {
00128 #ifdef G4VERBOSE
00129   if (GetVerboseLevel()>1) G4cout << "G4PhaseSpaceDecayChannel::OneBodyDecayIt()"<<G4endl;
00130 #endif
00131 
00132   //create parent G4DynamicParticle at rest
00133   G4ThreeVector dummy;
00134   G4DynamicParticle * parentparticle = new G4DynamicParticle( parent, dummy, 0.0);
00135   //create G4Decayproducts
00136   G4DecayProducts *products = new G4DecayProducts(*parentparticle);
00137   delete parentparticle;
00138 
00139   //create daughter G4DynamicParticle at rest
00140   G4DynamicParticle * daughterparticle = new G4DynamicParticle( daughters[0], dummy, 0.0);
00141   products->PushProducts(daughterparticle);
00142 
00143 #ifdef G4VERBOSE
00144   if (GetVerboseLevel()>1) {
00145      G4cout << "G4PhaseSpaceDecayChannel::OneBodyDecayIt ";
00146      G4cout << "  create decay products in rest frame " <<G4endl;
00147      products->DumpInfo();
00148   }
00149 #endif
00150   return products;
00151 }
00152 
00153 G4DecayProducts *G4PhaseSpaceDecayChannel::TwoBodyDecayIt()
00154 {
00155 #ifdef G4VERBOSE
00156   if (GetVerboseLevel()>1) G4cout << "G4PhaseSpaceDecayChannel::TwoBodyDecayIt()"<<G4endl;
00157 #endif
00158   // parent mass
00159   G4double parentmass = current_parent_mass;
00160   
00161   //daughters'mass
00162   G4double daughtermass[2]; 
00163   G4double daughtermomentum;
00164   daughtermass[0] = daughters_mass[0];
00165   daughtermass[1] = daughters_mass[1];
00166 
00167   //create parent G4DynamicParticle at rest
00168   G4ThreeVector dummy;
00169   G4DynamicParticle * parentparticle = new G4DynamicParticle( parent, dummy, 0.0);
00170   //create G4Decayproducts
00171   G4DecayProducts *products = new G4DecayProducts(*parentparticle);
00172   delete parentparticle;
00173 
00174   //calculate daughter momentum
00175   daughtermomentum = Pmx(parentmass,daughtermass[0],daughtermass[1]);
00176   if (daughtermomentum <0.0) {
00177 #ifdef G4VERBOSE
00178     if (GetVerboseLevel()>0) {
00179       G4cerr << "G4PhaseSpaceDecayChannel::TwoBodyDecayIt " 
00180              << "sum of daughter mass is larger than parent mass" << G4endl;
00181       G4cerr << "parent :" << parent->GetParticleName() << "  " << current_parent_mass/GeV << G4endl;
00182       G4cerr << "daughter 1 :" << daughters[0]->GetParticleName() << "  " << daughtermass[0]/GeV << G4endl;
00183       G4cerr << "daughter 2:" << daughters[1]->GetParticleName() << "  " << daughtermass[1]/GeV << G4endl;
00184     }
00185 #endif
00186     G4Exception("G4PhaseSpaceDecayChannel::TwoBodyDecayIt",
00187                 "PART112", JustWarning,
00188                 "Can not create decay products: sum of daughter mass is larger than parent mass");
00189     return products;
00190   }
00191 
00192   G4double costheta = 2.*G4UniformRand()-1.0;
00193   G4double sintheta = std::sqrt((1.0 - costheta)*(1.0 + costheta));
00194   G4double phi  = twopi*G4UniformRand()*rad;
00195   G4ThreeVector direction(sintheta*std::cos(phi),sintheta*std::sin(phi),costheta);
00196 
00197   //create daughter G4DynamicParticle 
00198   G4DynamicParticle * daughterparticle = new G4DynamicParticle( daughters[0], direction*daughtermomentum);
00199   products->PushProducts(daughterparticle);
00200   daughterparticle = new G4DynamicParticle( daughters[1], direction*(-1.0*daughtermomentum));
00201   products->PushProducts(daughterparticle);
00202 
00203 #ifdef G4VERBOSE
00204   if (GetVerboseLevel()>1) {
00205      G4cout << "G4PhaseSpaceDecayChannel::TwoBodyDecayIt ";
00206      G4cout << "  create decay products in rest frame " <<G4endl;
00207      products->DumpInfo();
00208   }
00209 #endif
00210   return products;
00211 }
00212 
00213 G4DecayProducts *G4PhaseSpaceDecayChannel::ThreeBodyDecayIt()
00214 // algorism of this code is originally written in GDECA3 of GEANT3
00215 {
00216 #ifdef G4VERBOSE
00217   if (GetVerboseLevel()>1) G4cout << "G4PhaseSpaceDecayChannel::ThreeBodyDecayIt()"<<G4endl;
00218 #endif
00219   // parent mass
00220   G4double parentmass = current_parent_mass;
00221   //daughters'mass
00222   G4double daughtermass[3]; 
00223   G4double sumofdaughtermass = 0.0;
00224   for (G4int index=0; index<3; index++){
00225     daughtermass[index] = daughters_mass[index];
00226     sumofdaughtermass += daughtermass[index];
00227   }
00228   
00229    //create parent G4DynamicParticle at rest
00230   G4ThreeVector dummy;
00231   G4DynamicParticle * parentparticle = new G4DynamicParticle( parent, dummy, 0.0);
00232 
00233 
00234   //create G4Decayproducts
00235   G4DecayProducts *products = new G4DecayProducts(*parentparticle);
00236   delete parentparticle;
00237 
00238   if (sumofdaughtermass >parentmass) {
00239 #ifdef G4VERBOSE
00240     if (GetVerboseLevel()>0) {
00241       G4cerr << "G4PhaseSpaceDecayChannel::ThreeBodyDecayIt " 
00242              << "sum of daughter mass is larger than parent mass" << G4endl;
00243       G4cerr << "parent :" << parent->GetParticleName() << "  " << current_parent_mass/GeV << G4endl;
00244       G4cerr << "daughter 1 :" << daughters[0]->GetParticleName() << "  " << daughtermass[0]/GeV << G4endl;
00245       G4cerr << "daughter 2:" << daughters[1]->GetParticleName() << "  " << daughtermass[1]/GeV << G4endl;
00246       G4cerr << "daughter 3:" << daughters[2]->GetParticleName() << "  " << daughtermass[2]/GeV << G4endl;
00247     }
00248 #endif
00249     G4Exception("G4PhaseSpaceDecayChannel::ThreeBodyDecayIt",
00250                 "PART112", JustWarning,
00251                 "Can not create decay products: sum of daughter mass is larger than parent mass");
00252     return products;
00253   }
00254 
00255 
00256 
00257   //calculate daughter momentum
00258   //  Generate two 
00259   G4double rd1, rd2, rd;
00260   G4double daughtermomentum[3];
00261   G4double momentummax=0.0, momentumsum = 0.0;
00262   G4double energy;
00263 
00264   do {
00265     rd1 = G4UniformRand();
00266     rd2 = G4UniformRand();
00267     if (rd2 > rd1) {
00268       rd  = rd1;
00269       rd1 = rd2;
00270       rd2 = rd;
00271     } 
00272     momentummax = 0.0;
00273     momentumsum = 0.0;
00274     // daughter 0
00275     energy = rd2*(parentmass - sumofdaughtermass);
00276     daughtermomentum[0] = std::sqrt(energy*energy + 2.0*energy* daughtermass[0]);
00277     if ( daughtermomentum[0] >momentummax )momentummax =  daughtermomentum[0];
00278     momentumsum  +=  daughtermomentum[0];
00279     // daughter 1
00280     energy = (1.-rd1)*(parentmass - sumofdaughtermass);
00281     daughtermomentum[1] = std::sqrt(energy*energy + 2.0*energy* daughtermass[1]);
00282     if ( daughtermomentum[1] >momentummax )momentummax =  daughtermomentum[1];
00283     momentumsum  +=  daughtermomentum[1];
00284     // daughter 2
00285     energy = (rd1-rd2)*(parentmass - sumofdaughtermass);
00286     daughtermomentum[2] = std::sqrt(energy*energy + 2.0*energy* daughtermass[2]);
00287     if ( daughtermomentum[2] >momentummax )momentummax =  daughtermomentum[2];
00288     momentumsum  +=  daughtermomentum[2];
00289   } while (momentummax >  momentumsum - momentummax );
00290   // output message
00291 #ifdef G4VERBOSE
00292   if (GetVerboseLevel()>1) {
00293     G4cout << "     daughter 0:" << daughtermomentum[0]/GeV << "[GeV/c]" <<G4endl;
00294     G4cout << "     daughter 1:" << daughtermomentum[1]/GeV << "[GeV/c]" <<G4endl;
00295     G4cout << "     daughter 2:" << daughtermomentum[2]/GeV << "[GeV/c]" <<G4endl;
00296     G4cout << "   momentum sum:" << momentumsum/GeV << "[GeV/c]" <<G4endl;
00297   }
00298 #endif
00299 
00300   //create daughter G4DynamicParticle 
00301   G4double costheta, sintheta, phi, sinphi, cosphi; 
00302   G4double costhetan, sinthetan, phin, sinphin, cosphin; 
00303   costheta = 2.*G4UniformRand()-1.0;
00304   sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
00305   phi  = twopi*G4UniformRand()*rad;
00306   sinphi = std::sin(phi);
00307   cosphi = std::cos(phi);
00308   G4ThreeVector direction0(sintheta*cosphi,sintheta*sinphi,costheta);
00309   G4DynamicParticle * daughterparticle 
00310          = new G4DynamicParticle( daughters[0], direction0*daughtermomentum[0]);
00311   products->PushProducts(daughterparticle);
00312 
00313   costhetan = (daughtermomentum[1]*daughtermomentum[1]-daughtermomentum[2]*daughtermomentum[2]-daughtermomentum[0]*daughtermomentum[0])/(2.0*daughtermomentum[2]*daughtermomentum[0]);
00314   sinthetan = std::sqrt((1.0-costhetan)*(1.0+costhetan));
00315   phin  = twopi*G4UniformRand()*rad;
00316   sinphin = std::sin(phin);
00317   cosphin = std::cos(phin);
00318   G4ThreeVector direction2;
00319   direction2.setX( sinthetan*cosphin*costheta*cosphi - sinthetan*sinphin*sinphi + costhetan*sintheta*cosphi); 
00320   direction2.setY( sinthetan*cosphin*costheta*sinphi + sinthetan*sinphin*cosphi + costhetan*sintheta*sinphi); 
00321   direction2.setZ( -sinthetan*cosphin*sintheta + costhetan*costheta);
00322   daughterparticle = new G4DynamicParticle( daughters[2], direction2*(daughtermomentum[2]/direction2.mag()));
00323   products->PushProducts(daughterparticle);
00324 
00325   daughterparticle = 
00326        new G4DynamicParticle( 
00327                 daughters[1],
00328                (direction0*daughtermomentum[0] + direction2*(daughtermomentum[2]/direction2.mag()))*(-1.0)   
00329                 );
00330   products->PushProducts(daughterparticle);
00331   
00332 #ifdef G4VERBOSE
00333   if (GetVerboseLevel()>1) {
00334      G4cout << "G4PhaseSpaceDecayChannel::ThreeBodyDecayIt ";
00335      G4cout << "  create decay products in rest frame " <<G4endl;
00336      products->DumpInfo();
00337   }
00338 #endif
00339   return products;
00340 }
00341 
00342 G4DecayProducts *G4PhaseSpaceDecayChannel::ManyBodyDecayIt()
00343 // algorithm of this code is originally written in FORTRAN by M.Asai
00344 //******************************************************************
00345 //  NBODY
00346 //   N-body phase space Monte-Carlo generator
00347 //  Makoto Asai 
00348 //   Hiroshima Institute of Technology
00349 //    (asai@kekvax.kek.jp)
00350 //  Revised release : 19/Apr/1995
00351 //
00352 {
00353   G4int index, index2;
00354 
00355   //return value
00356   G4DecayProducts *products;
00357 
00358 #ifdef G4VERBOSE
00359   if (GetVerboseLevel()>1) G4cout << "G4PhaseSpaceDecayChannel::ManyBodyDecayIt()"<<G4endl;
00360 #endif
00361   // parent mass
00362   G4double parentmass = current_parent_mass;
00363   //daughters'mass
00364   G4double *daughtermass = new G4double[numberOfDaughters]; 
00365   G4double sumofdaughtermass = 0.0;
00366   for (index=0; index<numberOfDaughters; index++){
00367     daughtermass[index] = daughters_mass[index];
00368     sumofdaughtermass += daughtermass[index];
00369   }
00370   
00371   //Calculate daughter momentum
00372   G4double *daughtermomentum = new G4double[numberOfDaughters];
00373   G4ThreeVector direction;  
00374   G4DynamicParticle **daughterparticle;
00375   G4double *sm = new G4double[numberOfDaughters];
00376   G4double tmas;
00377   G4double weight = 1.0;
00378   G4int    numberOfTry = 0;
00379  
00380   do {
00381     //Generate rundom number in descending order 
00382     G4double temp;
00383     G4double *rd = new G4double[numberOfDaughters];
00384     rd[0] = 1.0;
00385     for(index =1; index < numberOfDaughters -1; index++) rd[index] = G4UniformRand(); 
00386     rd[ numberOfDaughters -1] = 0.0;
00387     for(index =1; index < numberOfDaughters -1; index++) {
00388       for(index2 = index+1; index2 < numberOfDaughters; index2++) {
00389         if (rd[index] < rd[index2]){
00390           temp        = rd[index];
00391           rd[index]  = rd[index2];
00392           rd[index2] = temp;
00393         }
00394       }
00395     }
00396     
00397     //calcurate virtual mass 
00398     tmas = parentmass -  sumofdaughtermass;
00399     temp =  sumofdaughtermass; 
00400     for(index =0; index < numberOfDaughters; index++) {
00401       sm[index] = rd[index]*tmas + temp;
00402       temp -= daughtermass[index];
00403       if (GetVerboseLevel()>1) {
00404         G4cout << index << "  rundom number:" << rd[index]; 
00405         G4cout << "   virtual mass:" << sm[index]/GeV << "[GeV/c/c]" <<G4endl; 
00406       }
00407     }
00408     delete [] rd;
00409 
00410     //Calculate daughter momentum
00411     weight = 1.0;
00412     index =numberOfDaughters-1;
00413     daughtermomentum[index]= Pmx( sm[index-1],daughtermass[index-1], sm[index]);
00414 #ifdef G4VERBOSE
00415     if (GetVerboseLevel()>1) {
00416       G4cout << "     daughter " << index << ":" << *daughters_name[index];
00417       G4cout << " momentum:" << daughtermomentum[index]/GeV << "[GeV/c]" <<G4endl;
00418     }
00419 #endif
00420     for(index =numberOfDaughters-2; index>=0; index--) {
00421       // calculate 
00422       daughtermomentum[index]= Pmx( sm[index],daughtermass[index], sm[index +1]);
00423       if(daughtermomentum[index] < 0.0) {
00424         // !!! illegal momentum !!!
00425 #ifdef G4VERBOSE
00426         if (GetVerboseLevel()>0) {
00427           G4cout << "G4PhaseSpaceDecayChannel::ManyBodyDecayIt ";
00428           G4cout << "     can not calculate daughter momentum " <<G4endl;
00429           G4cout << "     parent:" << *parent_name;
00430           G4cout << " mass:" << parentmass/GeV << "[GeV/c/c]" <<G4endl;
00431           G4cout << "     daughter " << index << ":" << *daughters_name[index];
00432           G4cout << " mass:" << daughtermass[index]/GeV << "[GeV/c/c]" ;
00433           G4cout << " mass:" << daughtermomentum[index]/GeV << "[GeV/c]" <<G4endl;
00434         }
00435 #endif
00436         delete [] sm;
00437         delete [] daughtermass;
00438         delete [] daughtermomentum;
00439 
00440         G4Exception("G4PhaseSpaceDecayChannel::ManyBodyDecayIt",
00441                 "PART112", JustWarning,
00442                 "Can not create decay products: sum of daughter mass is larger than parent mass");
00443 
00444         return 0;   // Error detection
00445 
00446       } else {
00447         // calculate weight of this events
00448         weight *=  daughtermomentum[index]/sm[index];
00449 #ifdef G4VERBOSE
00450         if (GetVerboseLevel()>1) {
00451           G4cout << "     daughter " << index << ":" << *daughters_name[index];
00452           G4cout << " momentum:" << daughtermomentum[index]/GeV << "[GeV/c]" <<G4endl;
00453         }
00454 #endif
00455       }
00456     }
00457 
00458 #ifdef G4VERBOSE
00459     if (GetVerboseLevel()>1) {
00460       G4cout << "    weight: " << weight <<G4endl;
00461     }
00462 #endif
00463     
00464     // exit if number of Try exceeds 100
00465     if (numberOfTry++ >100) {
00466 #ifdef G4VERBOSE
00467       if (GetVerboseLevel()>0) {
00468         G4cout << "G4PhaseSpaceDecayChannel::ManyBodyDecayIt: ";
00469         G4cout << " can not determine Decay Kinematics " << G4endl;
00470         G4cout << "parent : " << *parent_name << G4endl;
00471         G4cout << "daughters : "; 
00472         for(index=0; index<numberOfDaughters; index++) {
00473           G4cout <<  *daughters_name[index] << " , ";
00474         }
00475         G4cout << G4endl;
00476       }
00477 #endif
00478       G4Exception("G4PhaseSpaceDecayChannel::ManyBodyDecayIt: ",
00479                   "PART113", JustWarning,
00480                   " Cannot decay :  Decay Kinematics cannot be calculated ");
00481       
00482       delete [] sm;
00483       delete [] daughtermass;
00484       delete [] daughtermomentum;
00485       return 0;  // Error detection
00486     }
00487   } while ( weight > G4UniformRand());
00488 
00489 #ifdef G4VERBOSE
00490   if (GetVerboseLevel()>1) {
00491       G4cout << "Start calulation of daughters momentum vector "<<G4endl;
00492   }
00493 #endif 
00494  
00495   G4double costheta, sintheta, phi;
00496   G4double beta;
00497   daughterparticle = new G4DynamicParticle*[numberOfDaughters];
00498 
00499   index = numberOfDaughters -2;
00500   costheta = 2.*G4UniformRand()-1.0;
00501   sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
00502   phi  = twopi*G4UniformRand()*rad;
00503   direction.setZ(costheta);
00504   direction.setY(sintheta*std::sin(phi));
00505   direction.setX(sintheta*std::cos(phi));
00506   daughterparticle[index] = new G4DynamicParticle( daughters[index], direction*daughtermomentum[index] );
00507   daughterparticle[index+1] = new G4DynamicParticle( daughters[index+1], direction*(-1.0*daughtermomentum[index]) );
00508 
00509   for (index = numberOfDaughters -3;  index >= 0; index--) {
00510     //calculate momentum direction
00511     costheta = 2.*G4UniformRand()-1.0;
00512     sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
00513     phi  = twopi*G4UniformRand()*rad;
00514     direction.setZ(costheta);
00515     direction.setY(sintheta*std::sin(phi));
00516     direction.setX(sintheta*std::cos(phi));
00517 
00518     // boost already created particles 
00519     beta = daughtermomentum[index];
00520     beta /= std::sqrt( daughtermomentum[index]*daughtermomentum[index] + sm[index+1]*sm[index+1] );
00521     for (index2 = index+1; index2<numberOfDaughters; index2++) {
00522       G4LorentzVector p4;
00523       // make G4LorentzVector for secondaries
00524       p4 = daughterparticle[index2]->Get4Momentum();
00525 
00526       // boost secondaries to  new frame 
00527       p4.boost( direction.x()*beta, direction.y()*beta, direction.z()*beta);
00528 
00529       // change energy/momentum
00530       daughterparticle[index2]->Set4Momentum(p4);
00531     }
00532     //create daughter G4DynamicParticle 
00533     daughterparticle[index]= new G4DynamicParticle( daughters[index], direction*(-1.0*daughtermomentum[index]));
00534   }
00535 
00536   //create G4Decayproducts
00537   G4DynamicParticle *parentparticle; 
00538   direction.setX(1.0);  direction.setY(0.0); direction.setZ(0.0);
00539   parentparticle = new G4DynamicParticle( parent, direction, 0.0);
00540   products = new G4DecayProducts(*parentparticle);
00541   delete parentparticle;
00542   for (index = 0; index<numberOfDaughters; index++) {
00543     products->PushProducts(daughterparticle[index]);
00544   }
00545 
00546 #ifdef G4VERBOSE
00547   if (GetVerboseLevel()>1) { 
00548     G4cout << "G4PhaseSpaceDecayChannel::ManyBodyDecayIt ";
00549     G4cout << "  create decay products in rest frame " << G4endl;
00550     products->DumpInfo();
00551   }
00552 #endif
00553 
00554   delete [] daughterparticle;
00555   delete [] daughtermomentum;
00556   delete [] daughtermass;
00557   delete [] sm;
00558   
00559   return products;
00560 }
00561 
00562 
00563 G4double G4PhaseSpaceDecayChannel::Pmx(G4double e, G4double p1, G4double p2)
00564 {
00565    // calcurate momentum of daughter particles in two-body decay
00566    G4double ppp = (e+p1+p2)*(e+p1-p2)*(e-p1+p2)*(e-p1-p2)/(4.0*e*e);
00567    if (ppp>0) return std::sqrt(ppp);
00568    else       return -1.;
00569 }
00570 
00571 
00572 

Generated on Mon May 27 17:49:18 2013 for Geant4 by  doxygen 1.4.7