#include <G4RPGFragmentation.hh>
Inheritance diagram for G4RPGFragmentation:
Public Member Functions | |
G4RPGFragmentation () | |
void | FragmentationIntegral (G4double, G4double, G4double, G4double) |
G4bool | ReactionStage (const G4HadProjectile *, G4ReactionProduct &, G4bool &, const G4DynamicParticle *, G4ReactionProduct &, G4bool &, const G4Nucleus &, G4ReactionProduct &, G4FastVector< G4ReactionProduct, 256 > &, G4int &, G4bool, G4ReactionProduct &) |
Definition at line 49 of file G4RPGFragmentation.hh.
G4RPGFragmentation::G4RPGFragmentation | ( | ) |
Definition at line 41 of file G4RPGFragmentation.cc.
00042 : G4RPGReaction() 00043 { 00044 for (G4int i = 0; i < 20; i++) dndl[i] = 0.0; 00045 }
Definition at line 49 of file G4RPGFragmentation.cc.
Referenced by ReactionStage().
00050 { 00051 pt = std::max( 0.001, pt ); 00052 G4double dx = 1./(19.*pt); 00053 G4double x; 00054 G4double term1; 00055 G4double term2; 00056 00057 for (G4int i = 1; i < 20; i++) { 00058 x = (G4double(i) - 0.5)*dx; 00059 term1 = 1. + parMass*parMass*x*x; 00060 term2 = pt*x*et*pt*x*et + pt*pt + secMass*secMass; 00061 dndl[i] = dx / std::sqrt( term1*term1*term1*term2 ) 00062 + dndl[i-1]; 00063 } 00064 }
G4bool G4RPGFragmentation::ReactionStage | ( | const G4HadProjectile * | , | |
G4ReactionProduct & | , | |||
G4bool & | , | |||
const G4DynamicParticle * | , | |||
G4ReactionProduct & | , | |||
G4bool & | , | |||
const G4Nucleus & | , | |||
G4ReactionProduct & | , | |||
G4FastVector< G4ReactionProduct, 256 > & | , | |||
G4int & | , | |||
G4bool | , | |||
G4ReactionProduct & | ||||
) |
Reimplemented from G4RPGReaction.
Definition at line 68 of file G4RPGFragmentation.cc.
References G4RPGReaction::AddBlackTrackParticles(), FragmentationIntegral(), G4cout, G4endl, G4Poisson(), G4UniformRand, G4RPGReaction::GenerateNBodyEventT(), G4Nucleus::GetA_asInt(), G4Nucleus::GetAnnihilationDTABlackTrackEnergy(), G4Nucleus::GetAnnihilationPNBlackTrackEnergy(), G4ParticleDefinition::GetBaryonNumber(), G4HadProjectile::GetDefinition(), G4ReactionProduct::GetDefinition(), G4Nucleus::GetDTABlackTrackEnergy(), G4RPGReaction::GetFinalStateNucleons(), G4ReactionProduct::GetKineticEnergy(), G4ReactionProduct::GetMass(), G4ReactionProduct::GetMomentum(), G4ParticleDefinition::GetParticleSubType(), G4ParticleDefinition::GetPDGMass(), G4Nucleus::GetPNBlackTrackEnergy(), G4ReactionProduct::GetSide(), G4ReactionProduct::GetTotalEnergy(), G4ReactionProduct::GetTotalMomentum(), G4Nucleus::GetZ_asInt(), G4RPGReaction::Isotropic(), G4Lambda::Lambda(), G4ReactionProduct::Lorentz(), G4Neutron::Neutron(), G4RPGReaction::normal(), G4PionMinus::PionMinus(), G4PionPlus::PionPlus(), G4PionZero::PionZero(), G4InuclParticleNames::pp, G4Proton::Proton(), G4ReactionProduct::SetDefinitionAndUpdateE(), G4FastVector< Type, N >::SetElement(), G4ReactionProduct::SetKineticEnergy(), G4ReactionProduct::SetMass(), G4ReactionProduct::SetMomentum(), G4ReactionProduct::SetSide(), and G4ReactionProduct::SetTotalEnergy().
Referenced by G4RPGInelastic::CalculateMomenta().
00080 { 00081 // 00082 // Based on H. Fesefeldt's original FORTRAN code GENXPT 00083 // 00084 // Generation of x- and pT- values for incident, target, and all secondary 00085 // particles using a simple single variable description E D3S/DP3= F(Q) 00086 // with Q^2 = (M*X)^2 + PT^2. Final state kinematics are produced by an 00087 // FF-type iterative cascade method 00088 // 00089 // Internal units are GeV 00090 // 00091 00092 // Protection in case no secondary has been created. In that case use 00093 // two-body scattering 00094 // 00095 if (vecLen == 0) return false; 00096 00097 G4ParticleDefinition* aPiMinus = G4PionMinus::PionMinus(); 00098 G4ParticleDefinition* aProton = G4Proton::Proton(); 00099 G4ParticleDefinition* aNeutron = G4Neutron::Neutron(); 00100 G4ParticleDefinition* aPiPlus = G4PionPlus::PionPlus(); 00101 G4ParticleDefinition* aPiZero = G4PionZero::PionZero(); 00102 00103 G4int i, l; 00104 G4bool veryForward = false; 00105 00106 const G4double ekOriginal = modifiedOriginal.GetKineticEnergy()/GeV; 00107 const G4double etOriginal = modifiedOriginal.GetTotalEnergy()/GeV; 00108 const G4double mOriginal = modifiedOriginal.GetMass()/GeV; 00109 const G4double pOriginal = modifiedOriginal.GetMomentum().mag()/GeV; 00110 G4double targetMass = targetParticle.GetDefinition()->GetPDGMass()/GeV; 00111 G4double centerofmassEnergy = std::sqrt( mOriginal*mOriginal + 00112 targetMass*targetMass + 00113 2.0*targetMass*etOriginal ); // GeV 00114 G4double currentMass = currentParticle.GetMass()/GeV; 00115 targetMass = targetParticle.GetMass()/GeV; 00116 00117 // Randomize the order of the secondary particles. 00118 // Note that the current and target particles are not affected. 00119 00120 for (i=0; i<vecLen; ++i) { 00121 G4int itemp = G4int( G4UniformRand()*vecLen ); 00122 G4ReactionProduct pTemp = *vec[itemp]; 00123 *vec[itemp] = *vec[i]; 00124 *vec[i] = pTemp; 00125 } 00126 00127 if (currentMass == 0.0 && targetMass == 0.0) { 00128 // Target and projectile have annihilated. Replace them with the first 00129 // two secondaries in the list. Current particle KE is maintained. 00130 00131 G4double ek = currentParticle.GetKineticEnergy(); 00132 G4ThreeVector mom = currentParticle.GetMomentum(); 00133 currentParticle = *vec[0]; 00134 currentParticle.SetSide(1); 00135 targetParticle = *vec[1]; 00136 targetParticle.SetSide(-1); 00137 for( i=0; i<(vecLen-2); ++i )*vec[i] = *vec[i+2]; 00138 G4ReactionProduct *temp = vec[vecLen-1]; 00139 delete temp; 00140 temp = vec[vecLen-2]; 00141 delete temp; 00142 vecLen -= 2; 00143 currentMass = currentParticle.GetMass()/GeV; 00144 targetMass = targetParticle.GetMass()/GeV; 00145 incidentHasChanged = true; 00146 targetHasChanged = true; 00147 currentParticle.SetKineticEnergy( ek ); 00148 currentParticle.SetMomentum(mom); 00149 veryForward = true; 00150 } 00151 const G4double atomicWeight = targetNucleus.GetA_asInt(); 00152 const G4double atomicNumber = targetNucleus.GetZ_asInt(); 00153 const G4double protonMass = aProton->GetPDGMass(); 00154 00155 if (originalIncident->GetDefinition()->GetParticleSubType() == "kaon" 00156 && G4UniformRand() >= 0.7) { 00157 G4ReactionProduct temp = currentParticle; 00158 currentParticle = targetParticle; 00159 targetParticle = temp; 00160 incidentHasChanged = true; 00161 targetHasChanged = true; 00162 currentMass = currentParticle.GetMass()/GeV; 00163 targetMass = targetParticle.GetMass()/GeV; 00164 } 00165 const G4double afc = std::min( 0.75, 00166 0.312+0.200*std::log(std::log(centerofmassEnergy*centerofmassEnergy))+ 00167 std::pow(centerofmassEnergy*centerofmassEnergy,1.5)/6000.0 ); 00168 00169 G4double freeEnergy = centerofmassEnergy-currentMass-targetMass; 00170 G4double forwardEnergy = freeEnergy/2.; 00171 G4int forwardCount = 1; // number of particles in forward hemisphere 00172 00173 G4double backwardEnergy = freeEnergy/2.; 00174 G4int backwardCount = 1; // number of particles in backward hemisphere 00175 00176 if(veryForward) { 00177 if(currentParticle.GetSide()==-1) 00178 { 00179 forwardEnergy += currentMass; 00180 forwardCount --; 00181 backwardEnergy -= currentMass; 00182 backwardCount ++; 00183 } 00184 if(targetParticle.GetSide()!=-1) 00185 { 00186 backwardEnergy += targetMass; 00187 backwardCount --; 00188 forwardEnergy -= targetMass; 00189 forwardCount ++; 00190 } 00191 } 00192 00193 for (i=0; i<vecLen; ++i) { 00194 if( vec[i]->GetSide() == -1 ) 00195 { 00196 ++backwardCount; 00197 backwardEnergy -= vec[i]->GetMass()/GeV; 00198 } else { 00199 ++forwardCount; 00200 forwardEnergy -= vec[i]->GetMass()/GeV; 00201 } 00202 } 00203 00204 // Check that sum of forward particle masses does not exceed forwardEnergy, 00205 // and similarly for backward particles. If so, move particles from one 00206 // hemisphere to another. 00207 00208 if (backwardEnergy < 0.0) { 00209 for (i = 0; i < vecLen; ++i) { 00210 if (vec[i]->GetSide() == -1) { 00211 backwardEnergy += vec[i]->GetMass()/GeV; 00212 --backwardCount; 00213 vec[i]->SetSide(1); 00214 forwardEnergy -= vec[i]->GetMass()/GeV; 00215 ++forwardCount; 00216 if (backwardEnergy > 0.0) break; 00217 } 00218 } 00219 } 00220 00221 if (forwardEnergy < 0.0) { 00222 for (i = 0; i < vecLen; ++i) { 00223 if (vec[i]->GetSide() == 1) { 00224 forwardEnergy += vec[i]->GetMass()/GeV; 00225 --forwardCount; 00226 vec[i]->SetSide(-1); 00227 backwardEnergy -= vec[i]->GetMass()/GeV; 00228 ++backwardCount; 00229 if (forwardEnergy > 0.0) break; 00230 } 00231 } 00232 } 00233 00234 // Special cases for reactions near threshold 00235 00236 // 1. There is only one secondary 00237 if (forwardEnergy > 0.0 && backwardEnergy < 0.0) { 00238 forwardEnergy += backwardEnergy; 00239 backwardEnergy = 0; 00240 } 00241 00242 // 2. Nuclear binding energy is large 00243 if (forwardEnergy + backwardEnergy < 0.0) return false; 00244 00245 00246 // forwardEnergy now represents the total energy in the forward reaction 00247 // hemisphere which is available for kinetic energy and particle creation. 00248 // Similarly for backwardEnergy. 00249 00250 // Add particles from the intranuclear cascade. 00251 // nuclearExcitationCount = number of new secondaries produced by nuclear 00252 // excitation 00253 // extraCount = number of nucleons within these new secondaries 00254 // 00255 // Note: eventually have to make sure that enough nucleons are available 00256 // in the case of small target nuclei 00257 00258 G4double xtarg; 00259 if( centerofmassEnergy < (2.0+G4UniformRand()) ) 00260 xtarg = afc * (std::pow(atomicWeight,0.33)-1.0) * (2.0*backwardCount+vecLen+2)/2.0; 00261 else 00262 xtarg = afc * (std::pow(atomicWeight,0.33)-1.0) * (2.0*backwardCount); 00263 if( xtarg <= 0.0 )xtarg = 0.01; 00264 G4int nuclearExcitationCount = G4Poisson( xtarg ); 00265 // To do: try reducing nuclearExcitationCount with increasing energy 00266 // to simulate cut-off of cascade 00267 if(atomicWeight<1.0001) nuclearExcitationCount = 0; 00268 G4int extraNucleonCount = 0; 00269 G4double extraNucleonMass = 0.0; 00270 00271 if (nuclearExcitationCount > 0) { 00272 const G4double nucsup[] = { 1.00, 0.7, 0.5, 0.4, 0.35, 0.3 }; 00273 const G4double psup[] = { 3., 6., 20., 50., 100., 1000. }; 00274 G4int momentumBin = 0; 00275 while( (momentumBin < 6) && 00276 (modifiedOriginal.GetTotalMomentum()/GeV > psup[momentumBin]) ) 00277 ++momentumBin; 00278 momentumBin = std::min( 5, momentumBin ); 00279 00280 // NOTE: in GENXPT, these new particles were given negative codes 00281 // here I use NewlyAdded = true instead 00282 // 00283 00284 for (i = 0; i < nuclearExcitationCount; ++i) { 00285 G4ReactionProduct * pVec = new G4ReactionProduct(); 00286 if (G4UniformRand() < nucsup[momentumBin]) { 00287 00288 if (G4UniformRand() > 1.0-atomicNumber/atomicWeight) 00289 pVec->SetDefinition( aProton ); 00290 else 00291 pVec->SetDefinition( aNeutron ); 00292 00293 // nucleon comes from nucleus - 00294 // do not subtract its mass from backward energy 00295 pVec->SetSide( -2 ); // -2 means backside nucleon 00296 ++extraNucleonCount; 00297 extraNucleonMass += pVec->GetMass()/GeV; 00298 // To do: Remove chosen nucleon from target nucleus 00299 pVec->SetNewlyAdded( true ); 00300 vec.SetElement( vecLen++, pVec ); 00301 ++backwardCount; 00302 00303 } else { 00304 00305 G4double ran = G4UniformRand(); 00306 if( ran < 0.3181 ) 00307 pVec->SetDefinition( aPiPlus ); 00308 else if( ran < 0.6819 ) 00309 pVec->SetDefinition( aPiZero ); 00310 else 00311 pVec->SetDefinition( aPiMinus ); 00312 00313 if (backwardEnergy > pVec->GetMass()/GeV) { 00314 backwardEnergy -= pVec->GetMass()/GeV; // pion mass comes from free energy 00315 ++backwardCount; 00316 pVec->SetSide( -1 ); // backside particle, but not a nucleon 00317 pVec->SetNewlyAdded( true ); 00318 vec.SetElement( vecLen++, pVec ); 00319 } 00320 00321 // To do: Change proton to neutron (or vice versa) in target nucleus depending 00322 // on pion charge 00323 } 00324 } 00325 } 00326 00327 // Define initial state vectors for Lorentz transformations 00328 // The pseudoParticles have non-standard masses, hence the "pseudo" 00329 00330 G4ReactionProduct pseudoParticle[8]; 00331 for (i = 0; i < 8; ++i) pseudoParticle[i].SetZero(); 00332 00333 pseudoParticle[0].SetMass( mOriginal*GeV ); 00334 pseudoParticle[0].SetMomentum( 0.0, 0.0, pOriginal*GeV ); 00335 pseudoParticle[0].SetTotalEnergy( 00336 std::sqrt( pOriginal*pOriginal + mOriginal*mOriginal )*GeV ); 00337 00338 pseudoParticle[1].SetMass(protonMass); // this could be targetMass 00339 pseudoParticle[1].SetTotalEnergy(protonMass); 00340 00341 pseudoParticle[3].SetMass(protonMass*(1+extraNucleonCount) ); 00342 pseudoParticle[3].SetTotalEnergy(protonMass*(1+extraNucleonCount) ); 00343 00344 pseudoParticle[2] = pseudoParticle[0] + pseudoParticle[1]; 00345 pseudoParticle[3] = pseudoParticle[3] + pseudoParticle[0]; 00346 00347 pseudoParticle[0].Lorentz( pseudoParticle[0], pseudoParticle[2] ); 00348 pseudoParticle[1].Lorentz( pseudoParticle[1], pseudoParticle[2] ); 00349 00350 // Main loop for 4-momentum generation 00351 // See Pitha-report (Aachen) for a detailed description of the method 00352 00353 G4double aspar, pt, et, x, pp, pp1, wgt; 00354 G4int innerCounter, outerCounter; 00355 G4bool eliminateThisParticle, resetEnergies, constantCrossSection; 00356 00357 G4double forwardKinetic = 0.0; 00358 G4double backwardKinetic = 0.0; 00359 00360 // Process the secondary particles in reverse order 00361 // The incident particle is done after the secondaries 00362 // Nucleons, including the target, in the backward hemisphere are also 00363 // done later 00364 00365 G4int backwardNucleonCount = 0; // number of nucleons in backward hemisphere 00366 G4double totalEnergy, kineticEnergy, vecMass; 00367 G4double phi; 00368 00369 for (i = vecLen-1; i >= 0; --i) { 00370 00371 if (vec[i]->GetNewlyAdded()) { // added from intranuclear cascade 00372 if (vec[i]->GetSide() == -2) { // its a nucleon 00373 if (backwardNucleonCount < 18) { 00374 if (vec[i]->GetDefinition()->GetParticleSubType() == "pi") { 00375 for (G4int j = 0; j < vecLen; j++) delete vec[j]; 00376 vecLen = 0; 00377 throw G4HadReentrentException(__FILE__, __LINE__, 00378 "G4RPGFragmentation::ReactionStage : a pion has been counted as a backward nucleon"); 00379 } 00380 vec[i]->SetSide(-3); 00381 ++backwardNucleonCount; 00382 continue; // Don't generate momenta for the first 17 backward 00383 // cascade nucleons. This gets done by the cluster 00384 // model later on. 00385 } 00386 } 00387 } 00388 00389 // Set pt and phi values, they are changed somewhat in the iteration loop 00390 // Set mass parameter for lambda fragmentation model 00391 00392 vecMass = vec[i]->GetMass()/GeV; 00393 G4double ran = -std::log(1.0-G4UniformRand())/3.5; 00394 00395 if (vec[i]->GetSide() == -2) { // backward nucleon 00396 aspar = 0.20; 00397 pt = std::sqrt( std::pow( ran, 1.2 ) ); 00398 00399 } else { // not a backward nucleon 00400 if (vec[i]->GetDefinition()->GetParticleSubType() == "pi") { 00401 aspar = 0.75; 00402 pt = std::sqrt( std::pow( ran, 1.7 ) ); 00403 } else if (vec[i]->GetDefinition()->GetParticleSubType() == "kaon") { 00404 aspar = 0.70; 00405 pt = std::sqrt( std::pow( ran, 1.7 ) ); 00406 } else { // vec[i] must be a baryon or ion 00407 aspar = 0.65; 00408 pt = std::sqrt( std::pow( ran, 1.5 ) ); 00409 } 00410 } 00411 00412 pt = std::max( 0.001, pt ); 00413 phi = G4UniformRand()*twopi; 00414 vec[i]->SetMomentum( pt*std::cos(phi)*GeV, pt*std::sin(phi)*GeV ); 00415 if (vec[i]->GetSide() > 0) 00416 et = pseudoParticle[0].GetTotalEnergy()/GeV; 00417 else 00418 et = pseudoParticle[1].GetTotalEnergy()/GeV; 00419 00420 // 00421 // Start of outer iteration loop 00422 // 00423 outerCounter = 0; 00424 eliminateThisParticle = true; 00425 resetEnergies = true; 00426 dndl[0] = 0.0; 00427 00428 while (++outerCounter < 3) { 00429 00430 FragmentationIntegral(pt, et, aspar, vecMass); 00431 innerCounter = 0; 00432 vec[i]->SetMomentum( pt*std::cos(phi)*GeV, pt*std::sin(phi)*GeV ); 00433 00434 // Start of inner iteration loop 00435 00436 while (++innerCounter < 7) { 00437 00438 ran = G4UniformRand()*dndl[19]; 00439 l = 1; 00440 while( ( ran > dndl[l] ) && ( l < 19 ) ) l++; 00441 x = (G4double(l-1) + G4UniformRand())/19.; 00442 if (vec[i]->GetSide() < 0) x *= -1.; 00443 vec[i]->SetMomentum( x*et*GeV ); // set the z-momentum 00444 totalEnergy = std::sqrt( x*et*x*et + pt*pt + vecMass*vecMass ); 00445 vec[i]->SetTotalEnergy( totalEnergy*GeV ); 00446 kineticEnergy = vec[i]->GetKineticEnergy()/GeV; 00447 00448 if (vec[i]->GetSide() > 0) { // forward side 00449 if( (forwardKinetic+kineticEnergy) < 0.95*forwardEnergy ) { 00450 // Leave at least 5% of the forward free energy for the projectile 00451 00452 pseudoParticle[4] = pseudoParticle[4] + (*vec[i]); 00453 forwardKinetic += kineticEnergy; 00454 outerCounter = 2; // leave outer loop 00455 eliminateThisParticle = false; // don't eliminate this particle 00456 resetEnergies = false; 00457 break; // leave inner loop 00458 } 00459 if( innerCounter > 5 )break; // leave inner loop 00460 if( backwardEnergy >= vecMass ) // switch sides 00461 { 00462 vec[i]->SetSide(-1); 00463 forwardEnergy += vecMass; 00464 backwardEnergy -= vecMass; 00465 ++backwardCount; 00466 } 00467 } else { // backward side 00468 // if (extraNucleonCount > 19) x = 0.999; 00469 // G4double xxx = 0.95+0.05*extraNucleonCount/20.0; 00470 // DHW: I think above lines were meant to be as follows: 00471 G4double xxx = 0.999; 00472 if (extraNucleonCount < 20) xxx = 0.95+0.05*extraNucleonCount/20.0; 00473 00474 if ((backwardKinetic+kineticEnergy) < xxx*backwardEnergy) { 00475 pseudoParticle[5] = pseudoParticle[5] + (*vec[i]); 00476 backwardKinetic += kineticEnergy; 00477 outerCounter = 2; // leave outer loop 00478 eliminateThisParticle = false; // don't eliminate this particle 00479 resetEnergies = false; 00480 break; // leave inner loop 00481 } 00482 if (innerCounter > 5) break; // leave inner loop 00483 if (forwardEnergy >= vecMass) { // switch sides 00484 vec[i]->SetSide(1); 00485 forwardEnergy -= vecMass; 00486 backwardEnergy += vecMass; 00487 backwardCount--; 00488 } 00489 } 00490 G4ThreeVector momentum = vec[i]->GetMomentum(); 00491 vec[i]->SetMomentum( momentum.x() * 0.9, momentum.y() * 0.9 ); 00492 pt *= 0.9; 00493 dndl[19] *= 0.9; 00494 } // closes inner loop 00495 00496 // If we get here, the inner loop has been done 6 times. 00497 // If necessary, reduce energies of the previously done particles if 00498 // they are lighter than protons or are in the forward hemisphere. 00499 // Then continue with outer loop. 00500 00501 if (resetEnergies) 00502 ReduceEnergiesOfSecondaries(i+1, forwardKinetic, backwardKinetic, 00503 vec, vecLen, 00504 pseudoParticle[4], pseudoParticle[5], 00505 pt); 00506 00507 } // closes outer loop 00508 00509 if (eliminateThisParticle && vec[i]->GetMayBeKilled()) { 00510 // not enough energy, eliminate this particle 00511 00512 if (vec[i]->GetSide() > 0) { 00513 --forwardCount; 00514 forwardEnergy += vecMass; 00515 } else { 00516 --backwardCount; 00517 if (vec[i]->GetSide() == -2) { 00518 --extraNucleonCount; 00519 extraNucleonMass -= vecMass; 00520 } else { 00521 backwardEnergy += vecMass; 00522 } 00523 } 00524 00525 for( G4int j=i; j<(vecLen-1); ++j )*vec[j] = *vec[j+1]; // shift up 00526 G4ReactionProduct* temp = vec[vecLen-1]; 00527 delete temp; 00528 // To do: modify target nucleus according to particle eliminated 00529 00530 if( --vecLen == 0 ){ 00531 G4cout << " FALSE RETURN DUE TO ENERGY BALANCE " << G4endl; 00532 return false; 00533 } // all the secondaries have been eliminated 00534 } 00535 } // closes main loop over secondaries 00536 00537 // Re-balance forward and backward energy if possible and if necessary 00538 00539 G4double forwardKEDiff = forwardEnergy - forwardKinetic; 00540 G4double backwardKEDiff = backwardEnergy - backwardKinetic; 00541 00542 if (forwardKEDiff < 0.0 || backwardKEDiff < 0.0) { 00543 ReduceEnergiesOfSecondaries(0, forwardKinetic, backwardKinetic, 00544 vec, vecLen, 00545 pseudoParticle[4], pseudoParticle[5], 00546 pt); 00547 00548 forwardKEDiff = forwardEnergy - forwardKinetic; 00549 backwardKEDiff = backwardEnergy - backwardKinetic; 00550 if (backwardKEDiff < 0.0) { 00551 if (forwardKEDiff + backwardKEDiff > 0.0) { 00552 backwardEnergy = backwardKinetic; 00553 forwardEnergy += backwardKEDiff; 00554 forwardKEDiff = forwardEnergy - forwardKinetic; 00555 backwardKEDiff = 0.0; 00556 } else { 00557 G4cout << " False return due to insufficient backward energy " << G4endl; 00558 return false; 00559 } 00560 } 00561 00562 if (forwardKEDiff < 0.0) { 00563 if (forwardKEDiff + backwardKEDiff > 0.0) { 00564 forwardEnergy = forwardKinetic; 00565 backwardEnergy += forwardKEDiff; 00566 backwardKEDiff = backwardEnergy - backwardKinetic; 00567 forwardKEDiff = 0.0; 00568 } else { 00569 G4cout << " False return due to insufficient forward energy " << G4endl; 00570 return false; 00571 } 00572 } 00573 } 00574 00575 // Generate momentum for incident (current) particle, which was placed 00576 // in the forward hemisphere. 00577 // Set mass parameter for lambda fragmentation model. 00578 // Set pt and phi values, which are changed somewhat in the iteration loop 00579 00580 G4double ran = -std::log(1.0-G4UniformRand()); 00581 if (currentParticle.GetDefinition()->GetParticleSubType() == "pi") { 00582 aspar = 0.60; 00583 pt = std::sqrt( std::pow( ran/6.0, 1.7 ) ); 00584 } else if (currentParticle.GetDefinition()->GetParticleSubType() == "kaon") { 00585 aspar = 0.50; 00586 pt = std::sqrt( std::pow( ran/5.0, 1.4 ) ); 00587 } else { 00588 aspar = 0.40; 00589 pt = std::sqrt( std::pow( ran/4.0, 1.2 ) ); 00590 } 00591 00592 phi = G4UniformRand()*twopi; 00593 currentParticle.SetMomentum(pt*std::cos(phi)*GeV, pt*std::sin(phi)*GeV); 00594 et = pseudoParticle[0].GetTotalEnergy()/GeV; 00595 dndl[0] = 0.0; 00596 vecMass = currentParticle.GetMass()/GeV; 00597 00598 FragmentationIntegral(pt, et, aspar, vecMass); 00599 00600 ran = G4UniformRand()*dndl[19]; 00601 l = 1; 00602 while( ( ran > dndl[l] ) && ( l < 19 ) ) l++; 00603 x = (G4double(l-1) + G4UniformRand())/19.; 00604 currentParticle.SetMomentum( x*et*GeV ); // set the z-momentum 00605 00606 if (forwardEnergy < forwardKinetic) { 00607 totalEnergy = vecMass + 0.04*std::fabs(normal()); 00608 G4cout << " Not enough forward energy: forwardEnergy = " 00609 << forwardEnergy << " forwardKinetic = " 00610 << forwardKinetic << " total energy left = " 00611 << backwardKEDiff + forwardKEDiff << G4endl; 00612 } else { 00613 totalEnergy = vecMass + forwardEnergy - forwardKinetic; 00614 forwardKinetic = forwardEnergy; 00615 } 00616 currentParticle.SetTotalEnergy( totalEnergy*GeV ); 00617 pp = std::sqrt(std::abs( totalEnergy*totalEnergy - vecMass*vecMass) )*GeV; 00618 pp1 = currentParticle.GetMomentum().mag(); 00619 00620 if (pp1 < 1.0e-6*GeV) { 00621 G4ThreeVector iso = Isotropic(pp); 00622 currentParticle.SetMomentum( iso.x(), iso.y(), iso.z() ); 00623 } else { 00624 currentParticle.SetMomentum( currentParticle.GetMomentum() * (pp/pp1) ); 00625 } 00626 pseudoParticle[4] = pseudoParticle[4] + currentParticle; 00627 00628 // Current particle now finished 00629 00630 // Begin target particle 00631 00632 if (backwardNucleonCount < 18) { 00633 targetParticle.SetSide(-3); 00634 ++backwardNucleonCount; 00635 00636 } else { 00637 // Set pt and phi values, they are changed somewhat in the iteration loop 00638 // Set mass parameter for lambda fragmentation model 00639 00640 vecMass = targetParticle.GetMass()/GeV; 00641 ran = -std::log(1.0-G4UniformRand()); 00642 aspar = 0.40; 00643 pt = std::max( 0.001, std::sqrt( std::pow( ran/4.0, 1.2 ) ) ); 00644 phi = G4UniformRand()*twopi; 00645 targetParticle.SetMomentum(pt*std::cos(phi)*GeV, pt*std::sin(phi)*GeV); 00646 et = pseudoParticle[1].GetTotalEnergy()/GeV; 00647 outerCounter = 0; 00648 innerCounter = 0; 00649 G4bool marginalEnergy = true; 00650 dndl[0] = 0.0; 00651 G4double xxx = 0.999; 00652 if( extraNucleonCount < 20 ) xxx = 0.95+0.05*extraNucleonCount/20.0; 00653 G4ThreeVector momentum; 00654 00655 while (++outerCounter < 4) { 00656 FragmentationIntegral(pt, et, aspar, vecMass); 00657 00658 for (innerCounter = 0; innerCounter < 6; innerCounter++) { 00659 ran = G4UniformRand()*dndl[19]; 00660 l = 1; 00661 while( ( ran > dndl[l] ) && ( l < 19 ) ) l++; 00662 x = -(G4double(l-1) + G4UniformRand())/19.; 00663 targetParticle.SetMomentum( x*et*GeV ); // set the z-momentum 00664 totalEnergy = std::sqrt(x*et*x*et + pt*pt + vecMass*vecMass); 00665 targetParticle.SetTotalEnergy( totalEnergy*GeV ); 00666 00667 if ((backwardKinetic+totalEnergy-vecMass) < xxx*backwardEnergy) { 00668 pseudoParticle[5] = pseudoParticle[5] + targetParticle; 00669 backwardKinetic += totalEnergy - vecMass; 00670 outerCounter = 3; // leave outer loop 00671 marginalEnergy = false; 00672 break; // leave inner loop 00673 } 00674 momentum = targetParticle.GetMomentum(); 00675 targetParticle.SetMomentum(momentum.x() * 0.9, momentum.y() * 0.9); 00676 pt *= 0.9; 00677 dndl[19] *= 0.9; 00678 } 00679 } 00680 00681 if (marginalEnergy) { 00682 G4cout << " Extra backward kinetic energy = " 00683 << 0.999*backwardEnergy - backwardKinetic << G4endl; 00684 totalEnergy = vecMass + 0.999*backwardEnergy - backwardKinetic; 00685 targetParticle.SetTotalEnergy(totalEnergy*GeV); 00686 pp = std::sqrt(std::abs(totalEnergy*totalEnergy - vecMass*vecMass) )*GeV; 00687 targetParticle.SetMomentum(momentum.x()/0.9, momentum.y()/0.9); 00688 pp1 = targetParticle.GetMomentum().mag(); 00689 targetParticle.SetMomentum(targetParticle.GetMomentum() * pp/pp1 ); 00690 pseudoParticle[5] = pseudoParticle[5] + targetParticle; 00691 backwardKinetic = 0.999*backwardEnergy; 00692 } 00693 00694 } 00695 00696 if (backwardEnergy < backwardKinetic) 00697 G4cout << " Backward Edif = " << backwardEnergy - backwardKinetic << G4endl; 00698 if (forwardEnergy != forwardKinetic) 00699 G4cout << " Forward Edif = " << forwardEnergy - forwardKinetic << G4endl; 00700 00701 // Target particle finished. 00702 // Now produce backward nucleons with a cluster model 00703 // ps[2] = CM frame of projectile + target 00704 // ps[3] = sum of projectile + nucleon cluster in lab frame 00705 // ps[6] = proj + cluster 4-vector boosted into CM frame of proj + targ 00706 // with secondaries, current and target particles subtracted 00707 // = total 4-momentum of backward nucleon cluster 00708 00709 pseudoParticle[6].Lorentz( pseudoParticle[3], pseudoParticle[2] ); 00710 pseudoParticle[6] = pseudoParticle[6] - pseudoParticle[4]; 00711 pseudoParticle[6] = pseudoParticle[6] - pseudoParticle[5]; 00712 00713 if (backwardNucleonCount == 1) { 00714 // Target particle is the only backward nucleon. Give it the remainder 00715 // of the backward kinetic energy. 00716 00717 G4double ekin = 00718 std::min(backwardEnergy-backwardKinetic, centerofmassEnergy/2.0-protonMass/GeV); 00719 00720 if( ekin < 0.04 )ekin = 0.04 * std::fabs( normal() ); 00721 vecMass = targetParticle.GetMass()/GeV; 00722 totalEnergy = ekin + vecMass; 00723 targetParticle.SetTotalEnergy( totalEnergy*GeV ); 00724 pp = std::sqrt(std::abs(totalEnergy*totalEnergy - vecMass*vecMass) )*GeV; 00725 pp1 = pseudoParticle[6].GetMomentum().mag(); 00726 if (pp1 < 1.0e-6*GeV) { 00727 G4ThreeVector iso = Isotropic(pp); 00728 targetParticle.SetMomentum( iso.x(), iso.y(), iso.z() ); 00729 } else { 00730 targetParticle.SetMomentum( pseudoParticle[6].GetMomentum() * (pp/pp1)); 00731 } 00732 pseudoParticle[5] = pseudoParticle[5] + targetParticle; 00733 00734 } else if (backwardNucleonCount > 1) { 00735 // Share remaining energy with up to 17 backward nucleons 00736 00737 G4int tempCount = 5; 00738 if (backwardNucleonCount < 5) tempCount = backwardNucleonCount; 00739 tempCount -= 2; 00740 00741 G4double clusterMass = 0.; 00742 if (targetParticle.GetSide() == -3) 00743 clusterMass = targetParticle.GetMass()/GeV; 00744 for (i = 0; i < vecLen; ++i) 00745 if (vec[i]->GetSide() == -3) clusterMass += vec[i]->GetMass()/GeV; 00746 clusterMass += backwardEnergy - backwardKinetic; 00747 00748 totalEnergy = pseudoParticle[6].GetTotalEnergy()/GeV; 00749 pseudoParticle[6].SetMass(clusterMass*GeV); 00750 00751 pp = std::sqrt(std::abs(totalEnergy*totalEnergy - 00752 clusterMass*clusterMass) )*GeV; 00753 pp1 = pseudoParticle[6].GetMomentum().mag(); 00754 if (pp1 < 1.0e-6*GeV) { 00755 G4ThreeVector iso = Isotropic(pp); 00756 pseudoParticle[6].SetMomentum(iso.x(), iso.y(), iso.z()); 00757 } else { 00758 pseudoParticle[6].SetMomentum(pseudoParticle[6].GetMomentum() * (-pp/pp1)); 00759 } 00760 00761 std::vector<G4ReactionProduct*> tempList; // ptrs to backward nucleons 00762 if (targetParticle.GetSide() == -3) tempList.push_back(&targetParticle); 00763 for (i = 0; i < vecLen; ++i) 00764 if (vec[i]->GetSide() == -3) tempList.push_back(vec[i]); 00765 00766 constantCrossSection = true; 00767 00768 if (tempList.size() > 1) { 00769 G4int n_entry = 0; 00770 wgt = GenerateNBodyEventT(pseudoParticle[6].GetMass(), 00771 constantCrossSection, tempList); 00772 00773 if (targetParticle.GetSide() == -3) { 00774 targetParticle = *tempList[0]; 00775 targetParticle.Lorentz(targetParticle, pseudoParticle[6]); 00776 n_entry++; 00777 } 00778 00779 for (i = 0; i < vecLen; ++i) { 00780 if (vec[i]->GetSide() == -3) { 00781 *vec[i] = *tempList[n_entry]; 00782 vec[i]->Lorentz(*vec[i], pseudoParticle[6]); 00783 n_entry++; 00784 } 00785 } 00786 } 00787 } else return false; 00788 00789 if (vecLen == 0) return false; // all the secondaries have been eliminated 00790 00791 // Lorentz transformation to lab frame 00792 00793 currentParticle.Lorentz( currentParticle, pseudoParticle[1] ); 00794 targetParticle.Lorentz( targetParticle, pseudoParticle[1] ); 00795 for (i = 0; i < vecLen; ++i) vec[i]->Lorentz(*vec[i], pseudoParticle[1]); 00796 00797 // Set leading strange particle flag. 00798 // leadFlag will be true if original particle and incident particle are 00799 // both strange, in which case the incident particle becomes the leading 00800 // particle. 00801 // leadFlag will also be true if the target particle is strange, but the 00802 // incident particle is not, in which case the target particle becomes the 00803 // leading particle. 00804 00805 G4bool leadingStrangeParticleHasChanged = true; 00806 if (leadFlag) 00807 { 00808 if (currentParticle.GetDefinition() == leadingStrangeParticle.GetDefinition()) 00809 leadingStrangeParticleHasChanged = false; 00810 if (leadingStrangeParticleHasChanged && 00811 (targetParticle.GetDefinition() == leadingStrangeParticle.GetDefinition()) ) 00812 leadingStrangeParticleHasChanged = false; 00813 if( leadingStrangeParticleHasChanged ) 00814 { 00815 for( i=0; i<vecLen; i++ ) 00816 { 00817 if( vec[i]->GetDefinition() == leadingStrangeParticle.GetDefinition() ) 00818 { 00819 leadingStrangeParticleHasChanged = false; 00820 break; 00821 } 00822 } 00823 } 00824 if( leadingStrangeParticleHasChanged ) 00825 { 00826 G4bool leadTest = 00827 (leadingStrangeParticle.GetDefinition()->GetParticleSubType() == "kaon" || 00828 leadingStrangeParticle.GetDefinition()->GetParticleSubType() == "pi"); 00829 G4bool targetTest = 00830 (targetParticle.GetDefinition()->GetParticleSubType() == "kaon" || 00831 targetParticle.GetDefinition()->GetParticleSubType() == "pi"); 00832 00833 // following modified by JLC 22-Oct-97 00834 00835 if( (leadTest&&targetTest) || !(leadTest||targetTest) ) // both true or both false 00836 { 00837 targetParticle.SetDefinitionAndUpdateE( leadingStrangeParticle.GetDefinition() ); 00838 targetHasChanged = true; 00839 } 00840 else 00841 { 00842 currentParticle.SetDefinitionAndUpdateE( leadingStrangeParticle.GetDefinition() ); 00843 incidentHasChanged = false; 00844 } 00845 } 00846 } // end of if( leadFlag ) 00847 00848 // Get number of final state nucleons and nucleons remaining in 00849 // target nucleus 00850 00851 std::pair<G4int, G4int> finalStateNucleons = 00852 GetFinalStateNucleons(originalTarget, vec, vecLen); 00853 00854 G4int protonsInFinalState = finalStateNucleons.first; 00855 G4int neutronsInFinalState = finalStateNucleons.second; 00856 00857 G4int numberofFinalStateNucleons = 00858 protonsInFinalState + neutronsInFinalState; 00859 00860 if (currentParticle.GetDefinition()->GetBaryonNumber() == 1 && 00861 targetParticle.GetDefinition()->GetBaryonNumber() == 1 && 00862 originalIncident->GetDefinition()->GetPDGMass() < 00863 G4Lambda::Lambda()->GetPDGMass()) 00864 numberofFinalStateNucleons++; 00865 00866 numberofFinalStateNucleons = std::max(1, numberofFinalStateNucleons); 00867 00868 G4int PinNucleus = std::max(0, 00869 G4int(targetNucleus.GetZ_asInt()) - protonsInFinalState); 00870 G4int NinNucleus = std::max(0, 00871 G4int(targetNucleus.GetA_asInt()-targetNucleus.GetZ_asInt()) - neutronsInFinalState); 00872 00873 pseudoParticle[3].SetMomentum( 0.0, 0.0, pOriginal*GeV ); 00874 pseudoParticle[3].SetMass( mOriginal*GeV ); 00875 pseudoParticle[3].SetTotalEnergy( 00876 std::sqrt( pOriginal*pOriginal + mOriginal*mOriginal )*GeV ); 00877 00878 G4ParticleDefinition * aOrgDef = modifiedOriginal.GetDefinition(); 00879 G4int diff = 0; 00880 if(aOrgDef == G4Proton::Proton() || aOrgDef == G4Neutron::Neutron() ) diff = 1; 00881 if(numberofFinalStateNucleons == 1) diff = 0; 00882 pseudoParticle[4].SetMomentum( 0.0, 0.0, 0.0 ); 00883 pseudoParticle[4].SetMass( protonMass*(numberofFinalStateNucleons-diff) ); 00884 pseudoParticle[4].SetTotalEnergy( protonMass*(numberofFinalStateNucleons-diff) ); 00885 00886 G4double theoreticalKinetic = 00887 pseudoParticle[3].GetTotalEnergy() + pseudoParticle[4].GetTotalEnergy() - 00888 currentParticle.GetMass() - targetParticle.GetMass(); 00889 for (i = 0; i < vecLen; ++i) theoreticalKinetic -= vec[i]->GetMass(); 00890 00891 G4double simulatedKinetic = 00892 currentParticle.GetKineticEnergy() + targetParticle.GetKineticEnergy(); 00893 for (i = 0; i < vecLen; ++i) 00894 simulatedKinetic += vec[i]->GetKineticEnergy(); 00895 00896 pseudoParticle[5] = pseudoParticle[3] + pseudoParticle[4]; 00897 pseudoParticle[3].Lorentz( pseudoParticle[3], pseudoParticle[5] ); 00898 pseudoParticle[4].Lorentz( pseudoParticle[4], pseudoParticle[5] ); 00899 00900 pseudoParticle[7].SetZero(); 00901 pseudoParticle[7] = pseudoParticle[7] + currentParticle; 00902 pseudoParticle[7] = pseudoParticle[7] + targetParticle; 00903 for (i = 0; i < vecLen; ++i) 00904 pseudoParticle[7] = pseudoParticle[7] + *vec[i]; 00905 00906 /* 00907 // This code does not appear to do anything to the energy or angle spectra 00908 if( vecLen <= 16 && vecLen > 0 ) 00909 { 00910 // must create a new set of ReactionProducts here because GenerateNBody will 00911 // modify the momenta for the particles, and we don't want to do this 00912 // 00913 G4ReactionProduct tempR[130]; 00914 tempR[0] = currentParticle; 00915 tempR[1] = targetParticle; 00916 for( i=0; i<vecLen; ++i )tempR[i+2] = *vec[i]; 00917 G4FastVector<G4ReactionProduct,256> tempV1; 00918 tempV1.Initialize( vecLen+2 ); 00919 G4int tempLen1 = 0; 00920 for( i=0; i<vecLen+2; ++i )tempV1.SetElement( tempLen1++, &tempR[i] ); 00921 constantCrossSection = true; 00922 00923 wgt = GenerateNBodyEvent(pseudoParticle[3].GetTotalEnergy() + 00924 pseudoParticle[4].GetTotalEnergy(), 00925 constantCrossSection, tempV1, tempLen1); 00926 if (wgt == -1) { 00927 G4double Qvalue = 0; 00928 for (i = 0; i < tempLen1; i++) Qvalue += tempV1[i]->GetMass(); 00929 wgt = GenerateNBodyEvent(Qvalue, 00930 constantCrossSection, tempV1, tempLen1); 00931 } 00932 if(wgt>-.5) 00933 { 00934 theoreticalKinetic = 0.0; 00935 for( i=0; i<tempLen1; ++i ) 00936 { 00937 pseudoParticle[6].Lorentz( *tempV1[i], pseudoParticle[4] ); 00938 theoreticalKinetic += pseudoParticle[6].GetKineticEnergy(); 00939 } 00940 } 00941 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen); 00942 } 00943 */ 00944 00945 // 00946 // Make sure that the kinetic energies are correct 00947 // 00948 00949 if (simulatedKinetic != 0.0) { 00950 wgt = theoreticalKinetic/simulatedKinetic; 00951 theoreticalKinetic = currentParticle.GetKineticEnergy() * wgt; 00952 simulatedKinetic = theoreticalKinetic; 00953 currentParticle.SetKineticEnergy(theoreticalKinetic); 00954 pp = currentParticle.GetTotalMomentum(); 00955 pp1 = currentParticle.GetMomentum().mag(); 00956 if (pp1 < 1.0e-6*GeV) { 00957 G4ThreeVector iso = Isotropic(pp); 00958 currentParticle.SetMomentum( iso.x(), iso.y(), iso.z() ); 00959 } else { 00960 currentParticle.SetMomentum(currentParticle.GetMomentum() * (pp/pp1)); 00961 } 00962 00963 theoreticalKinetic = targetParticle.GetKineticEnergy() * wgt; 00964 targetParticle.SetKineticEnergy(theoreticalKinetic); 00965 simulatedKinetic += theoreticalKinetic; 00966 pp = targetParticle.GetTotalMomentum(); 00967 pp1 = targetParticle.GetMomentum().mag(); 00968 00969 if (pp1 < 1.0e-6*GeV) { 00970 G4ThreeVector iso = Isotropic(pp); 00971 targetParticle.SetMomentum(iso.x(), iso.y(), iso.z() ); 00972 } else { 00973 targetParticle.SetMomentum(targetParticle.GetMomentum() * (pp/pp1) ); 00974 } 00975 00976 for (i = 0; i < vecLen; ++i ) { 00977 theoreticalKinetic = vec[i]->GetKineticEnergy() * wgt; 00978 simulatedKinetic += theoreticalKinetic; 00979 vec[i]->SetKineticEnergy(theoreticalKinetic); 00980 pp = vec[i]->GetTotalMomentum(); 00981 pp1 = vec[i]->GetMomentum().mag(); 00982 if( pp1 < 1.0e-6*GeV ) { 00983 G4ThreeVector iso = Isotropic(pp); 00984 vec[i]->SetMomentum(iso.x(), iso.y(), iso.z() ); 00985 } else { 00986 vec[i]->SetMomentum(vec[i]->GetMomentum() * (pp/pp1) ); 00987 } 00988 } 00989 } 00990 00991 // Rotate(numberofFinalStateNucleons, pseudoParticle[3].GetMomentum(), 00992 // modifiedOriginal, originalIncident, targetNucleus, 00993 // currentParticle, targetParticle, vec, vecLen ); 00994 00995 // Add black track particles 00996 // the total number of particles produced is restricted to 198 00997 // this may have influence on very high energies 00998 00999 if( atomicWeight >= 1.5 ) 01000 { 01001 // npnb is number of proton/neutron black track particles 01002 // ndta is the number of deuterons, tritons, and alphas produced 01003 // epnb is the kinetic energy available for proton/neutron black track 01004 // particles 01005 // edta is the kinetic energy available for deuteron/triton/alpha particles 01006 01007 G4int npnb = 0; 01008 G4int ndta = 0; 01009 01010 G4double epnb, edta; 01011 if (veryForward) { 01012 epnb = targetNucleus.GetAnnihilationPNBlackTrackEnergy(); 01013 edta = targetNucleus.GetAnnihilationDTABlackTrackEnergy(); 01014 } else { 01015 epnb = targetNucleus.GetPNBlackTrackEnergy(); 01016 edta = targetNucleus.GetDTABlackTrackEnergy(); 01017 } 01018 /* 01019 G4ReactionProduct* fudge = new G4ReactionProduct(); 01020 fudge->SetDefinition( aProton ); 01021 G4double TT = epnb + edta; 01022 G4double MM = fudge->GetMass()/GeV; 01023 fudge->SetTotalEnergy(MM*GeV + TT*GeV); 01024 G4double pzz = std::sqrt(TT*(TT + 2.*MM)); 01025 fudge->SetMomentum(0.0, 0.0, pzz*GeV); 01026 vec.SetElement(vecLen++, fudge); 01027 // G4cout << " Fudge = " << vec[vecLen-1]->GetKineticEnergy()/GeV 01028 << G4endl; 01029 */ 01030 01031 const G4double pnCutOff = 0.001; 01032 const G4double dtaCutOff = 0.001; 01033 // const G4double kineticMinimum = 1.e-6; 01034 // const G4double kineticFactor = -0.010; 01035 // G4double sprob = 0.0; // sprob = probability of self-absorption in 01036 // heavy molecules 01037 // Not currently used (DHW 9 June 2008) const G4double ekIncident = originalIncident->GetKineticEnergy()/GeV; 01038 // if (ekIncident >= 5.0) sprob = std::min(1.0, 0.6*std::log(ekIncident-4.0)); 01039 if (epnb > pnCutOff) 01040 { 01041 npnb = G4Poisson((1.5+1.25*numberofFinalStateNucleons)*epnb/(epnb+edta)); 01042 if (numberofFinalStateNucleons + npnb > atomicWeight) 01043 npnb = G4int(atomicWeight+0.00001 - numberofFinalStateNucleons); 01044 npnb = std::min( npnb, 127-vecLen ); 01045 } 01046 if( edta >= dtaCutOff ) 01047 { 01048 ndta = G4Poisson((1.5+1.25*numberofFinalStateNucleons)*edta/(epnb+edta)); 01049 ndta = std::min( ndta, 127-vecLen ); 01050 } 01051 if (npnb == 0 && ndta == 0) npnb = 1; 01052 01053 AddBlackTrackParticles(epnb, npnb, edta, ndta, modifiedOriginal, 01054 PinNucleus, NinNucleus, targetNucleus, 01055 vec, vecLen); 01056 } 01057 01058 // if( centerofmassEnergy <= (4.0+G4UniformRand()) ) 01059 // MomentumCheck( modifiedOriginal, currentParticle, targetParticle, 01060 // vec, vecLen ); 01061 // 01062 // calculate time delay for nuclear reactions 01063 // 01064 01065 if( (atomicWeight >= 1.5) && (atomicWeight <= 230.0) && (ekOriginal <= 0.2) ) 01066 currentParticle.SetTOF( 01067 1.0-500.0*std::exp(-ekOriginal/0.04)*std::log(G4UniformRand()) ); 01068 else 01069 currentParticle.SetTOF( 1.0 ); 01070 return true; 01071 01072 }