G4QMDReaction Class Reference

#include <G4QMDReaction.hh>

Inheritance diagram for G4QMDReaction:

G4HadronicInteraction

Public Member Functions

 G4QMDReaction ()
 ~G4QMDReaction ()
std::vector< G4QMDSystem * > GetFinalStates ()
G4HadFinalState * ApplyYourself (const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)
void UnUseGEM ()
void UseFRAG ()
void EnableFermiG ()
void SetTMAX (G4int i)
void SetDT (G4double t)
void SetEF (G4double x)

Detailed Description

Definition at line 60 of file G4QMDReaction.hh.

Constructor & Destructor Documentation

G4QMDReaction::G4QMDReaction (  ) 

Definition at line 45 of file G4QMDReaction.cc.

References G4ExcitationHandler::SetEvaporation().

00046 : G4HadronicInteraction("QMDModel")
00047 , system ( NULL )
00048 , deltaT ( 1 ) // in fsec (c=1)
00049 , maxTime ( 100 ) // will have maxTime-th time step
00050 , envelopF ( 1.05 ) // 10% for Peripheral reactions
00051 , gem ( true )
00052 , frag ( false )
00053 {
00054 
00055    //090331
00056    shenXS = new G4IonsShenCrossSection();
00057    //genspaXS = new G4GeneralSpaceNNCrossSection();
00058    piNucXS = new G4PiNuclearCrossSection();
00059    meanField = new G4QMDMeanField();
00060    collision = new G4QMDCollision();
00061 
00062    excitationHandler = new G4ExcitationHandler;
00063    evaporation = new G4Evaporation;
00064    excitationHandler->SetEvaporation( evaporation );
00065    setEvaporationCh();
00066 }

G4QMDReaction::~G4QMDReaction (  ) 

Definition at line 70 of file G4QMDReaction.cc.

00071 {
00072    delete evaporation; 
00073    delete excitationHandler;
00074    delete collision;
00075    delete meanField;
00076 }

Member Function Documentation

G4HadFinalState * G4QMDReaction::ApplyYourself ( const G4HadProjectile &  aTrack,
G4Nucleus &  targetNucleus 
) [virtual]

Implements G4HadronicInteraction.

Definition at line 80 of file G4QMDReaction.cc.

References G4HadFinalState::AddSecondary(), G4Alpha::Alpha(), CLHEP::boostOf(), G4ExcitationHandler::BreakItUp(), G4QMDNucleus::CalEnergyAndAngularMomentumInCM(), G4QMDCollision::CalKinematicsOfBinaryCollisions(), G4QMDSystem::Clear(), G4HadFinalState::Clear(), G4QMDMeanField::DoClusterJudgment(), G4QMDMeanField::DoPropagation(), G4UniformRand, G4QMDParticipant::Get4Momentum(), G4DynamicParticle::Get4Momentum(), G4HadProjectile::Get4Momentum(), G4Nucleus::GetA_asInt(), G4ParticleDefinition::GetAtomicMass(), G4ParticleDefinition::GetAtomicNumber(), G4QMDParticipant::GetDefinition(), G4DynamicParticle::GetDefinition(), G4HadProjectile::GetDefinition(), G4ParticleTable::GetIon(), G4IonTable::GetIonMass(), G4ParticleTable::GetIonTable(), G4HadProjectile::GetKineticEnergy(), G4QMDParticipant::GetMomentum(), G4DynamicParticle::GetMomentum(), G4QMDSystem::GetNOCollision(), G4QMDSystem::GetParticipant(), G4ParticleDefinition::GetParticleName(), G4ParticleTable::GetParticleTable(), G4ParticleDefinition::GetParticleType(), G4ParticleDefinition::GetPDGCharge(), G4ParticleDefinition::GetPDGMass(), G4QMDParticipant::GetPosition(), G4HadProjectile::GetTotalMomentum(), G4QMDSystem::GetTotalNumberOfParticipant(), G4QMDMeanField::GetTotalPotential(), G4Nucleus::GetZ_asInt(), G4INCL::Math::pi, G4QMDCollision::SetMeanField(), G4QMDMeanField::SetNucleus(), G4QMDSystem::SetParticipant(), G4QMDParticipant::SetProjectile(), G4HadFinalState::SetStatusChange(), G4QMDMeanField::SetSystem(), G4QMDParticipant::SetTarget(), G4QMDNucleus::SetTotalPotential(), stopAndKill, and G4HadronicInteraction::theParticleChange.

00081 {
00082    //G4cout << "G4QMDReaction::ApplyYourself" << G4endl;
00083 
00084    theParticleChange.Clear();
00085 
00086    system = new G4QMDSystem;
00087 
00088    G4int proj_Z = 0;
00089    G4int proj_A = 0;
00090    G4ParticleDefinition* proj_pd = ( G4ParticleDefinition* ) projectile.GetDefinition();
00091    if ( proj_pd->GetParticleType() == "nucleus" )
00092    {
00093       proj_Z = proj_pd->GetAtomicNumber();
00094       proj_A = proj_pd->GetAtomicMass();
00095    }
00096    else
00097    {
00098       proj_Z = (int)( proj_pd->GetPDGCharge()/eplus );
00099       proj_A = 1;
00100    }
00101    //G4int targ_Z = int ( target.GetZ() + 0.5 );
00102    //G4int targ_A = int ( target.GetN() + 0.5 );
00103    //migrate to integer A and Z (GetN_asInt returns number of neutrons in the nucleus since this) 
00104    G4int targ_Z = target.GetZ_asInt();
00105    G4int targ_A = target.GetA_asInt();
00106    G4ParticleDefinition* targ_pd = G4ParticleTable::GetParticleTable()->GetIon( targ_Z , targ_A , 0.0 );
00107 
00108 
00109    //G4NistManager* nistMan = G4NistManager::Instance();
00110 //   G4Element* G4NistManager::FindOrBuildElement( targ_Z );
00111 
00112    const G4DynamicParticle* proj_dp = new G4DynamicParticle ( proj_pd , projectile.Get4Momentum() );
00113    //const G4Element* targ_ele =  nistMan->FindOrBuildElement( targ_Z ); 
00114    //G4double aTemp = projectile.GetMaterial()->GetTemperature();
00115 
00116      //090331
00117   
00118    G4VCrossSectionDataSet* theXS = shenXS;
00119 
00120    if ( proj_pd->GetParticleType() == "meson" ) theXS = piNucXS;
00121 
00122    //G4double xs_0 = genspaXS->GetCrossSection ( proj_dp , targ_ele , aTemp );
00123    //G4double xs_0 = theXS->GetCrossSection ( proj_dp , targ_ele , aTemp );
00124    //110822 
00125    G4double xs_0 = theXS->GetIsoCrossSection ( proj_dp , targ_Z , targ_A );
00126 
00127      G4double bmax_0 = std::sqrt( xs_0 / pi );
00128      //std::cout << "bmax_0 in fm (fermi) " <<  bmax_0/fermi << std::endl;
00129 
00130      //delete proj_dp; 
00131 
00132    G4bool elastic = true;
00133    
00134    std::vector< G4QMDNucleus* > nucleuses; // Secondary nuceluses 
00135    G4ThreeVector boostToReac; // ReactionSystem (CM or NN); 
00136    G4ThreeVector boostBackToLAB; // Reaction System to LAB; 
00137 
00138    G4LorentzVector targ4p( G4ThreeVector( 0.0 ) , targ_pd->GetPDGMass()/GeV );
00139    G4ThreeVector boostLABtoCM = targ4p.findBoostToCM( proj_dp->Get4Momentum()/GeV ); // CM of target and proj; 
00140 
00141    G4double p1 = proj_dp->GetMomentum().mag()/GeV/proj_A; 
00142    G4double m1 = proj_dp->GetDefinition()->GetPDGMass()/GeV/proj_A;
00143    G4double e1 = std::sqrt( p1*p1 + m1*m1 ); 
00144    G4double e2 = targ_pd->GetPDGMass()/GeV/targ_A;
00145    G4double beta_nn = -p1 / ( e1+e2 );
00146 
00147    G4ThreeVector boostLABtoNN ( 0. , 0. , beta_nn ); // CM of NN; 
00148 
00149    G4double beta_nncm = ( - boostLABtoCM.beta() + boostLABtoNN.beta() ) / ( 1 - boostLABtoCM.beta() * boostLABtoNN.beta() ) ;  
00150 
00151    //std::cout << targ4p << std::endl; 
00152    //std::cout << proj_dp->Get4Momentum()<< std::endl; 
00153    //std::cout << beta_nncm << std::endl; 
00154    G4ThreeVector boostNNtoCM( 0. , 0. , beta_nncm ); // 
00155    G4ThreeVector boostCMtoNN( 0. , 0. , -beta_nncm ); // 
00156 
00157    boostToReac = boostLABtoNN; 
00158    boostBackToLAB = -boostLABtoNN; 
00159 
00160    delete proj_dp; 
00161 
00162    while ( elastic ) 
00163    {
00164 
00165 // impact parameter 
00166       //G4double bmax = 1.05*(bmax_0/fermi);  // 10% for Peripheral reactions
00167       G4double bmax = envelopF*(bmax_0/fermi);
00168       G4double b = bmax * std::sqrt ( G4UniformRand() );
00169 //071112
00170       //G4double b = 0;
00171       //G4double b = bmax;
00172       //G4double b = bmax/1.05 * 0.7 * G4UniformRand();
00173 
00174       //G4cout << "G4QMDRESULT bmax_0 = " << bmax_0/fermi << " fm, bmax = " << bmax << " fm , b = " << b  << " fm " << G4endl; 
00175 
00176       G4double plab = projectile.GetTotalMomentum()/GeV;
00177       G4double elab = ( projectile.GetKineticEnergy() + proj_pd->GetPDGMass() + targ_pd->GetPDGMass() )/GeV;
00178 
00179       calcOffSetOfCollision( b , proj_pd , targ_pd , plab , elab , bmax , boostCMtoNN );
00180 
00181 // Projectile
00182       G4LorentzVector proj4pLAB = projectile.Get4Momentum()/GeV;
00183 
00184       G4QMDGroundStateNucleus* proj(NULL); 
00185       if ( projectile.GetDefinition()->GetParticleType() == "nucleus" 
00186         || projectile.GetDefinition()->GetParticleName() == "proton"
00187         || projectile.GetDefinition()->GetParticleName() == "neutron" )
00188       {
00189 
00190          proj_Z = proj_pd->GetAtomicNumber();
00191          proj_A = proj_pd->GetAtomicMass();
00192 
00193          proj = new G4QMDGroundStateNucleus( proj_Z , proj_A );
00194          //proj->ShowParticipants();
00195 
00196 
00197          meanField->SetSystem ( proj );
00198          proj->SetTotalPotential( meanField->GetTotalPotential() );
00199          proj->CalEnergyAndAngularMomentumInCM();
00200 
00201       }
00202 
00203 // Target
00204       //G4int iz = int ( target.GetZ() );
00205       //G4int ia = int ( target.GetN() );
00206       //migrate to integer A and Z (GetN_asInt returns number of neutrons in the nucleus since this) 
00207       G4int iz = int ( target.GetZ_asInt() );
00208       G4int ia = int ( target.GetA_asInt() );
00209 
00210       G4QMDGroundStateNucleus* targ = new G4QMDGroundStateNucleus( iz , ia );
00211 
00212       meanField->SetSystem (targ );
00213       targ->SetTotalPotential( meanField->GetTotalPotential() );
00214       targ->CalEnergyAndAngularMomentumInCM();
00215    
00216       //G4LorentzVector targ4p( G4ThreeVector( 0.0 ) , targ->GetNuclearMass()/GeV );
00217 // Boost Vector to CM
00218       //boostToCM = targ4p.findBoostToCM( proj4pLAB );
00219 
00220 //    Target 
00221       for ( G4int i = 0 ; i < targ->GetTotalNumberOfParticipant() ; i++ )
00222       {
00223 
00224          G4ThreeVector p0 = targ->GetParticipant( i )->GetMomentum();
00225          G4ThreeVector r0 = targ->GetParticipant( i )->GetPosition();
00226 
00227          G4ThreeVector p ( p0.x() + coulomb_collision_px_targ 
00228                          , p0.y() 
00229                          , p0.z() * coulomb_collision_gamma_targ + coulomb_collision_pz_targ ); 
00230 
00231          G4ThreeVector r ( r0.x() + coulomb_collision_rx_targ 
00232                          , r0.y() 
00233                          , r0.z() / coulomb_collision_gamma_targ + coulomb_collision_rz_targ ); 
00234      
00235          system->SetParticipant( new G4QMDParticipant( targ->GetParticipant( i )->GetDefinition() , p , r ) );
00236          system->GetParticipant( i )->SetTarget();
00237 
00238       }
00239 
00240       G4LorentzVector proj4pCM = CLHEP::boostOf ( proj4pLAB , boostToReac );
00241       G4LorentzVector targ4pCM = CLHEP::boostOf ( targ4p , boostToReac );
00242 
00243 //    Projectile
00244       if ( proj != NULL )
00245       {
00246 
00247 //    projectile is nucleus
00248 
00249          for ( G4int i = 0 ; i < proj->GetTotalNumberOfParticipant() ; i++ )
00250          {
00251 
00252             G4ThreeVector p0 = proj->GetParticipant( i )->GetMomentum();
00253             G4ThreeVector r0 = proj->GetParticipant( i )->GetPosition();
00254 
00255             G4ThreeVector p ( p0.x() + coulomb_collision_px_proj 
00256                             , p0.y() 
00257                             , p0.z() * coulomb_collision_gamma_proj + coulomb_collision_pz_proj ); 
00258 
00259             G4ThreeVector r ( r0.x() + coulomb_collision_rx_proj 
00260                             , r0.y() 
00261                             , r0.z() / coulomb_collision_gamma_proj + coulomb_collision_rz_proj ); 
00262      
00263             system->SetParticipant( new G4QMDParticipant( proj->GetParticipant( i )->GetDefinition() , p  , r ) );
00264             system->GetParticipant ( i + targ->GetTotalNumberOfParticipant() )->SetProjectile();
00265          }
00266 
00267       }
00268       else
00269       {
00270 
00271 //       projectile is particle
00272 
00273          // avoid multiple set in "elastic" loop
00274          if ( system->GetTotalNumberOfParticipant() == targ->GetTotalNumberOfParticipant() )
00275          {
00276 
00277             G4int i = targ->GetTotalNumberOfParticipant(); 
00278       
00279             G4ThreeVector p0( 0 ); 
00280             G4ThreeVector r0( 0 );
00281 
00282             G4ThreeVector p ( p0.x() + coulomb_collision_px_proj 
00283                             , p0.y() 
00284                             , p0.z() * coulomb_collision_gamma_proj + coulomb_collision_pz_proj ); 
00285 
00286             G4ThreeVector r ( r0.x() + coulomb_collision_rx_proj 
00287                             , r0.y() 
00288                             , r0.z() / coulomb_collision_gamma_proj + coulomb_collision_rz_proj ); 
00289 
00290             system->SetParticipant( new G4QMDParticipant( (G4ParticleDefinition*)projectile.GetDefinition() , p , r ) );
00291             // This is not important becase only 1 projectile particle.
00292             system->GetParticipant ( i )->SetProjectile();
00293          }
00294 
00295       }
00296       //system->ShowParticipants();
00297 
00298       delete targ;
00299       delete proj;
00300 
00301    meanField->SetSystem ( system );
00302    collision->SetMeanField ( meanField );
00303 
00304 // Time Evolution 
00305    //std::cout << "Start time evolution " << std::endl;
00306    //system->ShowParticipants();
00307    for ( G4int i = 0 ; i < maxTime ; i++ )
00308    {
00309       //G4cout << " do Paropagate " << i << " th time step. " << G4endl;
00310       meanField->DoPropagation( deltaT );
00311       //system->ShowParticipants();
00312       collision->CalKinematicsOfBinaryCollisions( deltaT );
00313 
00314       if ( i / 10 * 10 == i ) 
00315       {
00316          //G4cout << i << " th time step. " << G4endl;
00317          //system->ShowParticipants();
00318       } 
00319       //system->ShowParticipants();
00320    }
00321    //system->ShowParticipants();
00322 
00323 
00324    //std::cout << "Doing Cluster Judgment " << std::endl;
00325 
00326    nucleuses = meanField->DoClusterJudgment();
00327 
00328 // Elastic Judgment  
00329 
00330    G4int numberOfSecondary = int ( nucleuses.size() ) + system->GetTotalNumberOfParticipant(); 
00331 
00332    G4int sec_a_Z = 0;
00333    G4int sec_a_A = 0;
00334    G4ParticleDefinition* sec_a_pd = NULL;
00335    G4int sec_b_Z = 0;
00336    G4int sec_b_A = 0;
00337    G4ParticleDefinition* sec_b_pd = NULL;
00338 
00339    if ( numberOfSecondary == 2 )
00340    {
00341 
00342       G4bool elasticLike_system = false;
00343       if ( nucleuses.size() == 2 ) 
00344       {
00345 
00346          sec_a_Z = nucleuses[0]->GetAtomicNumber();
00347          sec_a_A = nucleuses[0]->GetMassNumber();
00348          sec_b_Z = nucleuses[1]->GetAtomicNumber();
00349          sec_b_A = nucleuses[1]->GetMassNumber();
00350 
00351          if ( ( sec_a_Z == proj_Z && sec_a_A == proj_A && sec_b_Z == targ_Z && sec_b_A == targ_A )
00352            || ( sec_a_Z == targ_Z && sec_a_A == targ_A && sec_b_Z == proj_Z && sec_b_A == proj_A ) )
00353          {
00354             elasticLike_system = true;
00355          } 
00356 
00357       }
00358       else if ( nucleuses.size() == 1 ) 
00359       {
00360 
00361          sec_a_Z = nucleuses[0]->GetAtomicNumber();
00362          sec_a_A = nucleuses[0]->GetMassNumber();
00363          sec_b_pd = system->GetParticipant( 0 )->GetDefinition();
00364 
00365          if ( ( sec_a_Z == proj_Z && sec_a_A == proj_A && sec_b_pd == targ_pd )
00366            || ( sec_a_Z == targ_Z && sec_a_A == targ_A && sec_b_pd == proj_pd ) )
00367          {
00368             elasticLike_system = true;
00369          } 
00370 
00371       }  
00372       else
00373       {
00374 
00375          sec_a_pd = system->GetParticipant( 0 )->GetDefinition();
00376          sec_b_pd = system->GetParticipant( 1 )->GetDefinition();
00377  
00378          if ( ( sec_a_pd == proj_pd && sec_b_pd == targ_pd ) 
00379            || ( sec_a_pd == targ_pd && sec_b_pd == proj_pd ) ) 
00380          {
00381             elasticLike_system = true;
00382          } 
00383 
00384       } 
00385 
00386       if ( elasticLike_system == true )
00387       {
00388 
00389          G4bool elasticLike_energy = true;
00390 //    Cal ExcitationEnergy 
00391          for ( G4int i = 0 ; i < int ( nucleuses.size() ) ; i++ )
00392          { 
00393 
00394             //meanField->SetSystem( nucleuses[i] );
00395             meanField->SetNucleus( nucleuses[i] );
00396             //nucleuses[i]->SetTotalPotential( meanField->GetTotalPotential() );
00397             //nucleuses[i]->CalEnergyAndAngularMomentumInCM();
00398 
00399             if ( nucleuses[i]->GetExcitationEnergy()*GeV > 1.0*MeV ) elasticLike_energy = false;  
00400 
00401          } 
00402 
00403 //    Check Collision 
00404          G4bool withCollision = true;
00405          if ( system->GetNOCollision() == 0 ) withCollision = false;
00406 
00407 //    Final judegement for Inelasitc or Elastic;
00408 //
00409 //       ElasticLike without Collision 
00410          //if ( elasticLike_energy == true && withCollision == false ) elastic = true;  // ielst = 0
00411 //       ElasticLike with Collision 
00412          //if ( elasticLike_energy == true && withCollision == true ) elastic = true;   // ielst = 1 
00413 //       InelasticLike without Collision 
00414          //if ( elasticLike_energy == false ) elastic = false;                          // ielst = 2                
00415          if ( frag == true )
00416             if ( elasticLike_energy == false ) elastic = false;
00417 //       InelasticLike with Collision 
00418          if ( elasticLike_energy == false && withCollision == true ) elastic = false; // ielst = 3
00419 
00420       }
00421 
00422       }
00423       else
00424       {
00425 
00426 //       numberOfSecondary != 2 
00427          elastic = false;
00428 
00429       }
00430 
00431 //071115
00432       //G4cout << elastic << G4endl;
00433       // if elastic is true try again from sampling of impact parameter 
00434 
00435       if ( elastic == true )
00436       {
00437          // delete this nucleues
00438          for ( std::vector< G4QMDNucleus* >::iterator
00439                it = nucleuses.begin() ; it != nucleuses.end() ; it++ )
00440          {
00441             delete *it;
00442          }
00443          nucleuses.clear();
00444       }
00445    } 
00446 
00447 
00448 // Statical Decay Phase
00449 
00450    for ( std::vector< G4QMDNucleus* >::iterator it
00451        = nucleuses.begin() ; it != nucleuses.end() ; it++ )
00452    {
00453 
00454 /*
00455       std::cout << "G4QMDRESULT "
00456                 << (*it)->GetAtomicNumber() 
00457                 << " " 
00458                 << (*it)->GetMassNumber() 
00459                 << " " 
00460                 << (*it)->Get4Momentum() 
00461                 << " " 
00462                 << (*it)->Get4Momentum().vect() 
00463                 << " " 
00464                 << (*it)->Get4Momentum().restMass() 
00465                 << " " 
00466                 << (*it)->GetNuclearMass()/GeV 
00467                 << std::endl;
00468 */
00469 
00470       meanField->SetNucleus ( *it );
00471 
00472       if ( (*it)->GetAtomicNumber() == 0  // neutron cluster
00473         || (*it)->GetAtomicNumber() == (*it)->GetMassNumber() ) // proton cluster
00474       {
00475          // push back system 
00476          for ( G4int i = 0 ; i < (*it)->GetTotalNumberOfParticipant() ; i++ )
00477          {
00478             G4QMDParticipant* aP = new G4QMDParticipant( ( (*it)->GetParticipant( i ) )->GetDefinition() , ( (*it)->GetParticipant( i ) )->GetMomentum() , ( (*it)->GetParticipant( i ) )->GetPosition() );  
00479             system->SetParticipant ( aP );  
00480          } 
00481          continue;  
00482       }
00483 
00484       G4double nucleus_e = std::sqrt ( std::pow ( (*it)->GetNuclearMass()/GeV , 2 ) + std::pow ( (*it)->Get4Momentum().vect().mag() , 2 ) );
00485       G4LorentzVector nucleus_p4CM ( (*it)->Get4Momentum().vect() , nucleus_e ); 
00486 
00487 //      std::cout << "G4QMDRESULT nucleus deltaQ " << deltaQ << std::endl;
00488 
00489       G4int ia = (*it)->GetMassNumber();
00490       G4int iz = (*it)->GetAtomicNumber();
00491 
00492       G4LorentzVector lv ( G4ThreeVector( 0.0 ) , (*it)->GetExcitationEnergy()*GeV + G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass( iz , ia ) );
00493 
00494       G4Fragment* aFragment = new G4Fragment( ia , iz , lv );
00495 
00496       G4ReactionProductVector* rv;
00497       rv = excitationHandler->BreakItUp( *aFragment );
00498       G4bool notBreak = true;
00499       for ( G4ReactionProductVector::iterator itt
00500           = rv->begin() ; itt != rv->end() ; itt++ )
00501       {
00502 
00503           notBreak = false;
00504           // Secondary from this nucleus (*it) 
00505           G4ParticleDefinition* pd = (*itt)->GetDefinition();
00506 
00507           G4LorentzVector p4 ( (*itt)->GetMomentum()/GeV , (*itt)->GetTotalEnergy()/GeV );  //in nucleus(*it) rest system
00508           G4LorentzVector p4_CM = CLHEP::boostOf( p4 , -nucleus_p4CM.findBoostToCM() );  // Back to CM
00509           G4LorentzVector p4_LAB = CLHEP::boostOf( p4_CM , boostBackToLAB ); // Back to LAB  
00510 
00511 
00512 //090122
00513           //theParticleChange.AddSecondary( dp ); 
00514           if ( !( pd->GetAtomicNumber() == 4 && pd->GetAtomicMass() == 8 ) )
00515           {
00516              G4DynamicParticle* dp = new G4DynamicParticle( pd , p4_LAB*GeV );  
00517              theParticleChange.AddSecondary( dp ); 
00518           }
00519           else
00520           {
00521              //Be8 -> Alpha + Alpha + Q
00522              G4ThreeVector randomized_direction( G4UniformRand() , G4UniformRand() , G4UniformRand() );
00523              randomized_direction = randomized_direction.unit();
00524              G4double q_decay = (*itt)->GetMass() - 2*G4Alpha::Alpha()->GetPDGMass();
00525              G4double p_decay = std::sqrt ( std::pow(G4Alpha::Alpha()->GetPDGMass()+q_decay/2,2) - std::pow(G4Alpha::Alpha()->GetPDGMass() , 2 ) ); 
00526              G4LorentzVector p4_a1 ( p_decay*randomized_direction , G4Alpha::Alpha()->GetPDGMass()+q_decay/2 );  //in Be8 rest system
00527              
00528              G4LorentzVector p4_a1_Be8 = CLHEP::boostOf ( p4_a1/GeV , -p4.findBoostToCM() );
00529              G4LorentzVector p4_a1_CM = CLHEP::boostOf ( p4_a1_Be8 , -nucleus_p4CM.findBoostToCM() );
00530              G4LorentzVector p4_a1_LAB = CLHEP::boostOf ( p4_a1_CM , boostBackToLAB );
00531 
00532              G4LorentzVector p4_a2 ( -p_decay*randomized_direction , G4Alpha::Alpha()->GetPDGMass()+q_decay/2 );  //in Be8 rest system
00533              
00534              G4LorentzVector p4_a2_Be8 = CLHEP::boostOf ( p4_a2/GeV , -p4.findBoostToCM() );
00535              G4LorentzVector p4_a2_CM = CLHEP::boostOf ( p4_a2_Be8 , -nucleus_p4CM.findBoostToCM() );
00536              G4LorentzVector p4_a2_LAB = CLHEP::boostOf ( p4_a2_CM , boostBackToLAB );
00537              
00538              G4DynamicParticle* dp1 = new G4DynamicParticle( G4Alpha::Alpha() , p4_a1_LAB*GeV );  
00539              G4DynamicParticle* dp2 = new G4DynamicParticle( G4Alpha::Alpha() , p4_a2_LAB*GeV );  
00540              theParticleChange.AddSecondary( dp1 ); 
00541              theParticleChange.AddSecondary( dp2 ); 
00542           }
00543 //090122
00544 
00545 /*
00546           std::cout
00547                 << "Regist Secondary "
00548                 << (*itt)->GetDefinition()->GetParticleName()
00549                 << " "
00550                 << (*itt)->GetMomentum()/GeV
00551                 << " "
00552                 << (*itt)->GetKineticEnergy()/GeV
00553                 << " "
00554                 << (*itt)->GetMass()/GeV
00555                 << " "
00556                 << (*itt)->GetTotalEnergy()/GeV
00557                 << " "
00558                 << (*itt)->GetTotalEnergy()/GeV * (*itt)->GetTotalEnergy()/GeV
00559                  - (*itt)->GetMomentum()/GeV * (*itt)->GetMomentum()/GeV
00560                 << " "
00561                 << nucleus_p4CM.findBoostToCM() 
00562                 << " "
00563                 << p4
00564                 << " "
00565                 << p4_CM
00566                 << " "
00567                 << p4_LAB
00568                 << std::endl;
00569 */
00570 
00571       }
00572       if ( notBreak == true )
00573       {
00574 
00575          G4ParticleDefinition* pd = G4ParticleTable::GetParticleTable()->GetIon( (*it)->GetAtomicNumber() , (*it)->GetMassNumber(), (*it)->GetExcitationEnergy()*GeV );
00576          G4LorentzVector p4_CM = nucleus_p4CM;
00577          G4LorentzVector p4_LAB = CLHEP::boostOf( p4_CM , boostBackToLAB ); // Back to LAB  
00578          G4DynamicParticle* dp = new G4DynamicParticle( pd , p4_LAB*GeV );  
00579          theParticleChange.AddSecondary( dp ); 
00580 
00581       }
00582 
00583       for ( G4ReactionProductVector::iterator itt
00584           = rv->begin() ; itt != rv->end() ; itt++ )
00585       {
00586           delete *itt;
00587       }
00588       delete rv;
00589 
00590       delete aFragment;
00591 
00592    }
00593 
00594 
00595 
00596    for ( G4int i = 0 ; i < system->GetTotalNumberOfParticipant() ; i++ )
00597    {
00598 
00599       // Secondary particles 
00600 
00601       G4ParticleDefinition* pd = system->GetParticipant( i )->GetDefinition();
00602       G4LorentzVector p4_CM = system->GetParticipant( i )->Get4Momentum();
00603       G4LorentzVector p4_LAB = CLHEP::boostOf( p4_CM , boostBackToLAB );
00604       G4DynamicParticle* dp = new G4DynamicParticle( pd , p4_LAB*GeV );  
00605       theParticleChange.AddSecondary( dp ); 
00606 
00607 /*
00608       G4cout << "G4QMDRESULT "
00609       << "r" << i << " " << system->GetParticipant ( i ) -> GetPosition() << " "
00610       << "p" << i << " " << system->GetParticipant ( i ) -> Get4Momentum()
00611       << G4endl;
00612 */
00613 
00614    }
00615 
00616    for ( std::vector< G4QMDNucleus* >::iterator it
00617        = nucleuses.begin() ; it != nucleuses.end() ; it++ )
00618    {
00619       delete *it;  // delete nulceuse 
00620    }
00621    nucleuses.clear();
00622 
00623    system->Clear();
00624    delete system; 
00625 
00626    theParticleChange.SetStatusChange( stopAndKill );
00627 
00628    return &theParticleChange;
00629 
00630 }

void G4QMDReaction::EnableFermiG (  )  [inline]

Definition at line 74 of file G4QMDReaction.hh.

00074 { heg = true; setHighEnergyModel(); };

std::vector< G4QMDSystem* > G4QMDReaction::GetFinalStates (  ) 

void G4QMDReaction::SetDT ( G4double  t  )  [inline]

Definition at line 77 of file G4QMDReaction.hh.

00077 { deltaT = t; };

void G4QMDReaction::SetEF ( G4double  x  )  [inline]

Definition at line 78 of file G4QMDReaction.hh.

00078 { envelopF = x; };

void G4QMDReaction::SetTMAX ( G4int  i  )  [inline]

Definition at line 76 of file G4QMDReaction.hh.

00076 { maxTime = i; };

void G4QMDReaction::UnUseGEM (  )  [inline]

Definition at line 71 of file G4QMDReaction.hh.

00071 { gem = false; setEvaporationCh(); };

void G4QMDReaction::UseFRAG (  )  [inline]

Definition at line 72 of file G4QMDReaction.hh.

00072 { frag = true; };

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