00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054 #include <iostream>
00055 #include <signal.h>
00056
00057 #include "G4ReactionDynamics.hh"
00058 #include "G4AntiProton.hh"
00059 #include "G4AntiNeutron.hh"
00060 #include "G4PhysicalConstants.hh"
00061 #include "G4SystemOfUnits.hh"
00062 #include "Randomize.hh"
00063 #include "G4HadReentrentException.hh"
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105 G4bool G4ReactionDynamics::GenerateXandPt(
00106 G4FastVector<G4ReactionProduct,GHADLISTSIZE> &vec,
00107 G4int &vecLen,
00108 G4ReactionProduct &modifiedOriginal,
00109 const G4HadProjectile *originalIncident,
00110 G4ReactionProduct ¤tParticle,
00111 G4ReactionProduct &targetParticle,
00112 const G4DynamicParticle* originalTarget,
00113 const G4Nucleus &targetNucleus,
00114 G4bool &incidentHasChanged,
00115 G4bool &targetHasChanged,
00116 G4bool leadFlag,
00117 G4ReactionProduct &leadingStrangeParticle )
00118 {
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132 if(vecLen == 0) return false;
00133
00134 G4ParticleDefinition *aPiMinus = G4PionMinus::PionMinus();
00135 G4ParticleDefinition *aProton = G4Proton::Proton();
00136 G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
00137 G4ParticleDefinition *aPiPlus = G4PionPlus::PionPlus();
00138 G4ParticleDefinition *aPiZero = G4PionZero::PionZero();
00139
00140 G4int i, l;
00141 G4bool veryForward = false;
00142
00143 const G4double ekOriginal = modifiedOriginal.GetKineticEnergy()/GeV;
00144 const G4double etOriginal = modifiedOriginal.GetTotalEnergy()/GeV;
00145 const G4double mOriginal = modifiedOriginal.GetMass()/GeV;
00146 const G4double pOriginal = modifiedOriginal.GetMomentum().mag()/GeV;
00147 G4double targetMass = targetParticle.GetDefinition()->GetPDGMass()/GeV;
00148 G4double centerofmassEnergy = std::sqrt( mOriginal*mOriginal +
00149 targetMass*targetMass +
00150 2.0*targetMass*etOriginal );
00151 G4double currentMass = currentParticle.GetMass()/GeV;
00152 targetMass = targetParticle.GetMass()/GeV;
00153
00154
00155
00156
00157 for( i=0; i<vecLen; ++i )
00158 {
00159 G4int itemp = G4int( G4UniformRand()*vecLen );
00160 G4ReactionProduct pTemp = *vec[itemp];
00161 *vec[itemp] = *vec[i];
00162 *vec[i] = pTemp;
00163
00164 }
00165
00166 if( currentMass == 0.0 && targetMass == 0.0 )
00167 {
00168
00169
00170
00171 G4double ek = currentParticle.GetKineticEnergy();
00172 G4ThreeVector mom = currentParticle.GetMomentum();
00173 currentParticle = *vec[0];
00174 targetParticle = *vec[1];
00175 for( i=0; i<(vecLen-2); ++i )*vec[i] = *vec[i+2];
00176 G4ReactionProduct *temp = vec[vecLen-1];
00177 delete temp;
00178 temp = vec[vecLen-2];
00179 delete temp;
00180 vecLen -= 2;
00181 currentMass = currentParticle.GetMass()/GeV;
00182 targetMass = targetParticle.GetMass()/GeV;
00183 incidentHasChanged = true;
00184 targetHasChanged = true;
00185 currentParticle.SetKineticEnergy( ek );
00186 currentParticle.SetMomentum( mom );
00187 veryForward = true;
00188 }
00189 const G4double atomicWeight = G4double(targetNucleus.GetA_asInt());
00190 const G4double atomicNumber = G4double(targetNucleus.GetZ_asInt());
00191 const G4double protonMass = aProton->GetPDGMass()/MeV;
00192
00193 if (originalIncident->GetDefinition()->GetParticleSubType() == "kaon"
00194 && G4UniformRand() >= 0.7) {
00195 G4ReactionProduct temp = currentParticle;
00196 currentParticle = targetParticle;
00197 targetParticle = temp;
00198 incidentHasChanged = true;
00199 targetHasChanged = true;
00200 currentMass = currentParticle.GetMass()/GeV;
00201 targetMass = targetParticle.GetMass()/GeV;
00202 }
00203 const G4double afc = std::min( 0.75,
00204 0.312+0.200*std::log(std::log(centerofmassEnergy*centerofmassEnergy))+
00205 std::pow(centerofmassEnergy*centerofmassEnergy,1.5)/6000.0 );
00206
00207 G4double freeEnergy = centerofmassEnergy-currentMass-targetMass;
00208 G4double forwardEnergy = freeEnergy/2.;
00209 G4int forwardCount = 1;
00210
00211 G4double backwardEnergy = freeEnergy/2.;
00212 G4int backwardCount = 1;
00213
00214 if(veryForward)
00215 {
00216 if(currentParticle.GetSide()==-1)
00217 {
00218 forwardEnergy += currentMass;
00219 forwardCount --;
00220 backwardEnergy -= currentMass;
00221 backwardCount ++;
00222 }
00223 if(targetParticle.GetSide()!=-1)
00224 {
00225 backwardEnergy += targetMass;
00226 backwardCount --;
00227 forwardEnergy -= targetMass;
00228 forwardCount ++;
00229 }
00230 }
00231
00232 for( i=0; i<vecLen; ++i )
00233 {
00234 if( vec[i]->GetSide() == -1 )
00235 {
00236 ++backwardCount;
00237 backwardEnergy -= vec[i]->GetMass()/GeV;
00238 } else {
00239 ++forwardCount;
00240 forwardEnergy -= vec[i]->GetMass()/GeV;
00241 }
00242 }
00243
00244
00245
00246
00247
00248 G4double xtarg;
00249 if( centerofmassEnergy < (2.0+G4UniformRand()) )
00250 xtarg = afc * (std::pow(atomicWeight,0.33)-1.0) * (2.0*backwardCount+vecLen+2)/2.0;
00251 else
00252 xtarg = afc * (std::pow(atomicWeight,0.33)-1.0) * (2.0*backwardCount);
00253 if( xtarg <= 0.0 )xtarg = 0.01;
00254 G4int nuclearExcitationCount = Poisson( xtarg );
00255 if(atomicWeight<1.0001) nuclearExcitationCount = 0;
00256 G4int extraNucleonCount = 0;
00257 G4double extraNucleonMass = 0.0;
00258 if( nuclearExcitationCount > 0 )
00259 {
00260 const G4double nucsup[] = { 1.00, 0.7, 0.5, 0.4, 0.35, 0.3 };
00261 const G4double psup[] = { 3., 6., 20., 50., 100., 1000. };
00262 G4int momentumBin = 0;
00263 while( (momentumBin < 6) &&
00264 (modifiedOriginal.GetTotalMomentum()/GeV > psup[momentumBin]) )
00265 ++momentumBin;
00266 momentumBin = std::min( 5, momentumBin );
00267
00268
00269
00270
00271
00272 for( i=0; i<nuclearExcitationCount; ++i )
00273 {
00274 G4ReactionProduct * pVec = new G4ReactionProduct();
00275 if( G4UniformRand() < nucsup[momentumBin] )
00276 {
00277 if( G4UniformRand() > 1.0-atomicNumber/atomicWeight )
00278 pVec->SetDefinition( aProton );
00279 else
00280 pVec->SetDefinition( aNeutron );
00281 pVec->SetSide( -2 );
00282 ++extraNucleonCount;
00283 backwardEnergy += pVec->GetMass()/GeV;
00284 extraNucleonMass += pVec->GetMass()/GeV;
00285 }
00286 else
00287 {
00288 G4double ran = G4UniformRand();
00289 if( ran < 0.3181 )
00290 pVec->SetDefinition( aPiPlus );
00291 else if( ran < 0.6819 )
00292 pVec->SetDefinition( aPiZero );
00293 else
00294 pVec->SetDefinition( aPiMinus );
00295 pVec->SetSide( -1 );
00296 }
00297 pVec->SetNewlyAdded( true );
00298 vec.SetElement( vecLen++, pVec );
00299
00300 backwardEnergy -= pVec->GetMass()/GeV;
00301 ++backwardCount;
00302 }
00303 }
00304
00305
00306
00307 G4int is, iskip;
00308 while( forwardEnergy <= 0.0 )
00309 {
00310
00311 iskip = G4int(G4UniformRand()*forwardCount) + 1;
00312 is = 0;
00313 G4int forwardParticlesLeft = 0;
00314 for( i=(vecLen-1); i>=0; --i )
00315 {
00316 if( vec[i]->GetSide() == 1 && vec[i]->GetMayBeKilled())
00317 {
00318 forwardParticlesLeft = 1;
00319 if( ++is == iskip )
00320 {
00321 forwardEnergy += vec[i]->GetMass()/GeV;
00322 for( G4int j=i; j<(vecLen-1); j++ )*vec[j] = *vec[j+1];
00323 --forwardCount;
00324 delete vec[vecLen-1];
00325 if( --vecLen == 0 )return false;
00326 break;
00327 }
00328 }
00329 }
00330
00331 if( forwardParticlesLeft == 0 )
00332 {
00333 G4int iremove = -1;
00334 for (i = 0; i < vecLen; i++) {
00335 if (vec[i]->GetDefinition()->GetParticleSubType() == "pi") {
00336 iremove = i;
00337 break;
00338 }
00339 }
00340 if (iremove == -1) {
00341 for (i = 0; i < vecLen; i++) {
00342 if (vec[i]->GetDefinition()->GetParticleSubType() == "kaon") {
00343 iremove = i;
00344 break;
00345 }
00346 }
00347 }
00348 if (iremove == -1) iremove = 0;
00349
00350 forwardEnergy += vec[iremove]->GetMass()/GeV;
00351 if (vec[iremove]->GetSide() > 0) --forwardCount;
00352
00353 for (i = iremove; i < vecLen-1; i++) *vec[i] = *vec[i+1];
00354 delete vec[vecLen-1];
00355 vecLen--;
00356 if (vecLen == 0) return false;
00357 break;
00358 }
00359 }
00360
00361
00362 while( backwardEnergy <= 0.0 )
00363 {
00364
00365 iskip = G4int(G4UniformRand()*backwardCount) + 1;
00366 is = 0;
00367 G4int backwardParticlesLeft = 0;
00368 for( i=(vecLen-1); i>=0; --i )
00369 {
00370 if( vec[i]->GetSide() < 0 && vec[i]->GetMayBeKilled())
00371 {
00372 backwardParticlesLeft = 1;
00373 if( ++is == iskip )
00374 {
00375 if( vec[i]->GetSide() == -2 )
00376 {
00377 --extraNucleonCount;
00378 extraNucleonMass -= vec[i]->GetMass()/GeV;
00379 backwardEnergy -= vec[i]->GetTotalEnergy()/GeV;
00380 }
00381 backwardEnergy += vec[i]->GetTotalEnergy()/GeV;
00382 for( G4int j=i; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];
00383 --backwardCount;
00384 delete vec[vecLen-1];
00385 if( --vecLen == 0 )return false;
00386 break;
00387 }
00388 }
00389 }
00390
00391 if( backwardParticlesLeft == 0 )
00392 {
00393 G4int iremove = -1;
00394 for (i = 0; i < vecLen; i++) {
00395 if (vec[i]->GetDefinition()->GetParticleSubType() == "pi") {
00396 iremove = i;
00397 break;
00398 }
00399 }
00400 if (iremove == -1) {
00401 for (i = 0; i < vecLen; i++) {
00402 if (vec[i]->GetDefinition()->GetParticleSubType() == "kaon") {
00403 iremove = i;
00404 break;
00405 }
00406 }
00407 }
00408 if (iremove == -1) iremove = 0;
00409
00410 backwardEnergy += vec[iremove]->GetMass()/GeV;
00411 if (vec[iremove]->GetSide() > 0) --backwardCount;
00412
00413 for (i = iremove; i < vecLen-1; i++) *vec[i] = *vec[i+1];
00414 delete vec[vecLen-1];
00415 vecLen--;
00416 if (vecLen == 0) return false;
00417 break;
00418 }
00419 }
00420
00421
00422
00423
00424
00425
00426 G4ReactionProduct pseudoParticle[10];
00427 for( i=0; i<10; ++i )pseudoParticle[i].SetZero();
00428
00429 pseudoParticle[0].SetMass( mOriginal*GeV );
00430 pseudoParticle[0].SetMomentum( 0.0, 0.0, pOriginal*GeV );
00431 pseudoParticle[0].SetTotalEnergy(
00432 std::sqrt( pOriginal*pOriginal + mOriginal*mOriginal )*GeV );
00433
00434 pseudoParticle[1].SetMass( protonMass*MeV );
00435 pseudoParticle[1].SetTotalEnergy( protonMass*MeV );
00436
00437 pseudoParticle[3].SetMass( protonMass*(1+extraNucleonCount)*MeV );
00438 pseudoParticle[3].SetTotalEnergy( protonMass*(1+extraNucleonCount)*MeV );
00439
00440 pseudoParticle[8].SetMomentum( 1.0*GeV, 0.0, 0.0 );
00441
00442 pseudoParticle[2] = pseudoParticle[0] + pseudoParticle[1];
00443 pseudoParticle[3] = pseudoParticle[3] + pseudoParticle[0];
00444
00445 pseudoParticle[0].Lorentz( pseudoParticle[0], pseudoParticle[2] );
00446 pseudoParticle[1].Lorentz( pseudoParticle[1], pseudoParticle[2] );
00447
00448 G4double dndl[20];
00449
00450
00451
00452
00453 G4double aspar, pt, et, x, pp, pp1, rmb, wgt;
00454 G4int innerCounter, outerCounter;
00455 G4bool eliminateThisParticle, resetEnergies, constantCrossSection;
00456
00457 G4double forwardKinetic = 0.0, backwardKinetic = 0.0;
00458
00459
00460
00461
00462
00463 G4double binl[20] = {0.,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.11,1.25,
00464 1.43,1.67,2.0,2.5,3.33,5.00,10.00};
00465 G4int backwardNucleonCount = 0;
00466 G4double totalEnergy, kineticEnergy, vecMass;
00467
00468 for( i=(vecLen-1); i>=0; --i )
00469 {
00470 G4double phi = G4UniformRand()*twopi;
00471 if( vec[i]->GetNewlyAdded() )
00472 {
00473 if( vec[i]->GetSide() == -2 )
00474 {
00475 if( backwardNucleonCount < 18 )
00476 {
00477 if (vec[i]->GetDefinition()->GetParticleSubType() == "pi") {
00478 for(G4int j=0; j<vecLen; j++) delete vec[j];
00479 vecLen = 0;
00480 throw G4HadReentrentException(__FILE__, __LINE__,
00481 "G4ReactionDynamics::GenerateXandPt : a pion has been counted as a backward nucleon");
00482 }
00483 vec[i]->SetSide( -3 );
00484 ++backwardNucleonCount;
00485 continue;
00486 }
00487 }
00488 }
00489
00490
00491
00492
00493 vecMass = vec[i]->GetMass()/GeV;
00494 G4double ran = -std::log(1.0-G4UniformRand())/3.5;
00495 if( vec[i]->GetSide() == -2 )
00496 {
00497 if (vec[i]->GetDefinition()->GetParticleSubType() == "kaon" ||
00498 vec[i]->GetDefinition()->GetParticleSubType() == "pi") {
00499 aspar = 0.75;
00500 pt = std::sqrt( std::pow( ran, 1.7 ) );
00501 } else {
00502 aspar = 0.20;
00503 pt = std::sqrt( std::pow( ran, 1.2 ) );
00504 }
00505
00506 } else {
00507
00508 if (vec[i]->GetDefinition()->GetParticleSubType() == "pi") {
00509 aspar = 0.75;
00510 pt = std::sqrt( std::pow( ran, 1.7 ) );
00511 } else if (vec[i]->GetDefinition()->GetParticleSubType() == "kaon") {
00512 aspar = 0.70;
00513 pt = std::sqrt( std::pow( ran, 1.7 ) );
00514 } else {
00515 aspar = 0.65;
00516 pt = std::sqrt( std::pow( ran, 1.5 ) );
00517 }
00518 }
00519 pt = std::max( 0.001, pt );
00520 vec[i]->SetMomentum( pt*std::cos(phi)*GeV, pt*std::sin(phi)*GeV );
00521 for( G4int j=0; j<20; ++j )binl[j] = j/(19.*pt);
00522 if( vec[i]->GetSide() > 0 )
00523 et = pseudoParticle[0].GetTotalEnergy()/GeV;
00524 else
00525 et = pseudoParticle[1].GetTotalEnergy()/GeV;
00526 dndl[0] = 0.0;
00527
00528
00529
00530 outerCounter = 0;
00531 eliminateThisParticle = true;
00532 resetEnergies = true;
00533 while( ++outerCounter < 3 )
00534 {
00535 for( l=1; l<20; ++l )
00536 {
00537 x = (binl[l]+binl[l-1])/2.;
00538 pt = std::max( 0.001, pt );
00539 if( x > 1.0/pt )
00540 dndl[l] += dndl[l-1];
00541 else
00542 dndl[l] = et * aspar/std::sqrt( std::pow((1.+aspar*x*aspar*x),3) )
00543 * (binl[l]-binl[l-1]) / std::sqrt( pt*x*et*pt*x*et + pt*pt + vecMass*vecMass )
00544 + dndl[l-1];
00545 }
00546 innerCounter = 0;
00547 vec[i]->SetMomentum( pt*std::cos(phi)*GeV, pt*std::sin(phi)*GeV );
00548
00549
00550
00551 while( ++innerCounter < 7 )
00552 {
00553 ran = G4UniformRand()*dndl[19];
00554 l = 1;
00555 while( ( ran >= dndl[l] ) && ( l < 20 ) )l++;
00556 l = std::min( 19, l );
00557 x = std::min( 1.0, pt*(binl[l-1] + G4UniformRand()*(binl[l]-binl[l-1]) ) );
00558 if( vec[i]->GetSide() < 0 )x *= -1.;
00559 vec[i]->SetMomentum( x*et*GeV );
00560 totalEnergy = std::sqrt( x*et*x*et + pt*pt + vecMass*vecMass );
00561 vec[i]->SetTotalEnergy( totalEnergy*GeV );
00562 kineticEnergy = vec[i]->GetKineticEnergy()/GeV;
00563 if( vec[i]->GetSide() > 0 )
00564 {
00565 if( (forwardKinetic+kineticEnergy) < 0.95*forwardEnergy )
00566 {
00567 pseudoParticle[4] = pseudoParticle[4] + (*vec[i]);
00568 forwardKinetic += kineticEnergy;
00569 pseudoParticle[6] = pseudoParticle[4] + pseudoParticle[5];
00570 pseudoParticle[6].SetMomentum( 0.0 );
00571 phi = pseudoParticle[6].Angle( pseudoParticle[8] );
00572 if( pseudoParticle[6].GetMomentum().y()/MeV < 0.0 )phi = twopi - phi;
00573 phi += pi + normal()*pi/12.0;
00574 if( phi > twopi )phi -= twopi;
00575 if( phi < 0.0 )phi = twopi - phi;
00576 outerCounter = 2;
00577 eliminateThisParticle = false;
00578 resetEnergies = false;
00579 break;
00580 }
00581 if( innerCounter > 5 )break;
00582 if( backwardEnergy >= vecMass )
00583 {
00584 vec[i]->SetSide( -1 );
00585 forwardEnergy += vecMass;
00586 backwardEnergy -= vecMass;
00587 ++backwardCount;
00588 }
00589 } else {
00590 if( extraNucleonCount > 19 )
00591 x = 0.999;
00592 G4double xxx = 0.95+0.05*extraNucleonCount/20.0;
00593 if( (backwardKinetic+kineticEnergy) < xxx*backwardEnergy )
00594 {
00595 pseudoParticle[5] = pseudoParticle[5] + (*vec[i]);
00596 backwardKinetic += kineticEnergy;
00597 pseudoParticle[6] = pseudoParticle[4] + pseudoParticle[5];
00598 pseudoParticle[6].SetMomentum( 0.0 );
00599 phi = pseudoParticle[6].Angle( pseudoParticle[8] );
00600 if( pseudoParticle[6].GetMomentum().y()/MeV < 0.0 )phi = twopi - phi;
00601 phi += pi + normal() * pi / 12.0;
00602 if( phi > twopi )phi -= twopi;
00603 if( phi < 0.0 )phi = twopi - phi;
00604 outerCounter = 2;
00605 eliminateThisParticle = false;
00606 resetEnergies = false;
00607 break;
00608 }
00609 if( innerCounter > 5 )break;
00610 if( forwardEnergy >= vecMass )
00611 {
00612 vec[i]->SetSide( 1 );
00613 forwardEnergy -= vecMass;
00614 backwardEnergy += vecMass;
00615 backwardCount--;
00616 }
00617 }
00618 G4ThreeVector momentum = vec[i]->GetMomentum();
00619 vec[i]->SetMomentum( momentum.x() * 0.9, momentum.y() * 0.9 );
00620 pt *= 0.9;
00621 dndl[19] *= 0.9;
00622 }
00623 if( resetEnergies )
00624 {
00625
00626
00627
00628
00629
00630 forwardKinetic = 0.0;
00631 backwardKinetic = 0.0;
00632 pseudoParticle[4].SetZero();
00633 pseudoParticle[5].SetZero();
00634 for( l=i+1; l<vecLen; ++l )
00635 {
00636 if (vec[l]->GetSide() > 0 ||
00637 vec[l]->GetDefinition()->GetParticleSubType() == "kaon" ||
00638 vec[l]->GetDefinition()->GetParticleSubType() == "pi") {
00639
00640 G4double tempMass = vec[l]->GetMass()/MeV;
00641 totalEnergy = 0.95*vec[l]->GetTotalEnergy()/MeV + 0.05*tempMass;
00642 totalEnergy = std::max( tempMass, totalEnergy );
00643 vec[l]->SetTotalEnergy( totalEnergy*MeV );
00644 pp = std::sqrt( std::abs( totalEnergy*totalEnergy - tempMass*tempMass ) );
00645 pp1 = vec[l]->GetMomentum().mag()/MeV;
00646 if( pp1 < 1.0e-6*GeV )
00647 {
00648 G4double costheta = 2.*G4UniformRand() - 1.;
00649 G4double sintheta = std::sqrt(1. - costheta*costheta);
00650 G4double phi2 = twopi*G4UniformRand();
00651 vec[l]->SetMomentum( pp*sintheta*std::cos(phi2)*MeV,
00652 pp*sintheta*std::sin(phi2)*MeV,
00653 pp*costheta*MeV ) ;
00654 } else {
00655 vec[l]->SetMomentum( vec[l]->GetMomentum() * (pp/pp1) );
00656 }
00657 G4double px = vec[l]->GetMomentum().x()/MeV;
00658 G4double py = vec[l]->GetMomentum().y()/MeV;
00659 pt = std::max( 1.0, std::sqrt( px*px + py*py ) )/GeV;
00660 if( vec[l]->GetSide() > 0 )
00661 {
00662 forwardKinetic += vec[l]->GetKineticEnergy()/GeV;
00663 pseudoParticle[4] = pseudoParticle[4] + (*vec[l]);
00664 } else {
00665 backwardKinetic += vec[l]->GetKineticEnergy()/GeV;
00666 pseudoParticle[5] = pseudoParticle[5] + (*vec[l]);
00667 }
00668 }
00669 }
00670 }
00671 }
00672
00673 if( eliminateThisParticle && vec[i]->GetMayBeKilled())
00674 {
00675 if( vec[i]->GetSide() > 0 )
00676 {
00677 --forwardCount;
00678 forwardEnergy += vecMass;
00679 } else {
00680 if( vec[i]->GetSide() == -2 )
00681 {
00682 --extraNucleonCount;
00683 extraNucleonMass -= vecMass;
00684 backwardEnergy -= vecMass;
00685 }
00686 --backwardCount;
00687 backwardEnergy += vecMass;
00688 }
00689 for( G4int j=i; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];
00690 G4ReactionProduct *temp = vec[vecLen-1];
00691 delete temp;
00692
00693 if( --vecLen == 0 )return false;
00694 pseudoParticle[6] = pseudoParticle[4] + pseudoParticle[5];
00695 pseudoParticle[6].SetMomentum( 0.0 );
00696 phi = pseudoParticle[6].Angle( pseudoParticle[8] );
00697 if( pseudoParticle[6].GetMomentum().y()/MeV < 0.0 )phi = twopi - phi;
00698 phi += pi + normal() * pi / 12.0;
00699 if( phi > twopi )phi -= twopi;
00700 if( phi < 0.0 )phi = twopi - phi;
00701 }
00702 }
00703
00704
00705
00706
00707
00708
00709 G4double phi = G4UniformRand()*twopi;
00710 G4double ran = -std::log(1.0-G4UniformRand());
00711 if (currentParticle.GetDefinition()->GetParticleSubType() == "pi") {
00712 aspar = 0.60;
00713 pt = std::sqrt( std::pow( ran/6.0, 1.7 ) );
00714 } else if (currentParticle.GetDefinition()->GetParticleSubType() == "kaon") {
00715 aspar = 0.50;
00716 pt = std::sqrt( std::pow( ran/5.0, 1.4 ) );
00717 } else {
00718 aspar = 0.40;
00719 pt = std::sqrt( std::pow( ran/4.0, 1.2 ) );
00720 }
00721
00722 for( G4int j=0; j<20; ++j )binl[j] = j/(19.*pt);
00723 currentParticle.SetMomentum( pt*std::cos(phi)*GeV, pt*std::sin(phi)*GeV );
00724 et = pseudoParticle[0].GetTotalEnergy()/GeV;
00725 dndl[0] = 0.0;
00726 vecMass = currentParticle.GetMass()/GeV;
00727 for( l=1; l<20; ++l )
00728 {
00729 x = (binl[l]+binl[l-1])/2.;
00730 if( x > 1.0/pt )
00731 dndl[l] += dndl[l-1];
00732 else
00733 dndl[l] = aspar/std::sqrt( std::pow((1.+sqr(aspar*x)),3) ) *
00734 (binl[l]-binl[l-1]) * et / std::sqrt( pt*x*et*pt*x*et + pt*pt + vecMass*vecMass ) +
00735 dndl[l-1];
00736 }
00737 ran = G4UniformRand()*dndl[19];
00738 l = 1;
00739 while( (ran>dndl[l]) && (l<20) )l++;
00740 l = std::min( 19, l );
00741 x = std::min( 1.0, pt*(binl[l-1] + G4UniformRand()*(binl[l]-binl[l-1]) ) );
00742 currentParticle.SetMomentum( x*et*GeV );
00743 if( forwardEnergy < forwardKinetic )
00744 totalEnergy = vecMass + 0.04*std::fabs(normal());
00745 else
00746 totalEnergy = vecMass + forwardEnergy - forwardKinetic;
00747 currentParticle.SetTotalEnergy( totalEnergy*GeV );
00748 pp = std::sqrt( std::abs( totalEnergy*totalEnergy - vecMass*vecMass ) )*GeV;
00749 pp1 = currentParticle.GetMomentum().mag()/MeV;
00750 if( pp1 < 1.0e-6*GeV )
00751 {
00752 G4double costheta = 2.*G4UniformRand() - 1.;
00753 G4double sintheta = std::sqrt(1. - costheta*costheta);
00754 G4double phi2 = twopi*G4UniformRand();
00755 currentParticle.SetMomentum( pp*sintheta*std::cos(phi2)*MeV,
00756 pp*sintheta*std::sin(phi2)*MeV,
00757 pp*costheta*MeV ) ;
00758 } else {
00759 currentParticle.SetMomentum( currentParticle.GetMomentum() * (pp/pp1) );
00760 }
00761 pseudoParticle[4] = pseudoParticle[4] + currentParticle;
00762
00763
00764
00765
00766
00767
00768
00769 if( backwardNucleonCount < 18 )
00770 {
00771 targetParticle.SetSide( -3 );
00772 ++backwardNucleonCount;
00773 }
00774 else
00775 {
00776
00777
00778
00779 vecMass = targetParticle.GetMass()/GeV;
00780 ran = -std::log(1.0-G4UniformRand());
00781 aspar = 0.40;
00782 pt = std::max( 0.001, std::sqrt( std::pow( ran/4.0, 1.2 ) ) );
00783 targetParticle.SetMomentum( pt*std::cos(phi)*GeV, pt*std::sin(phi)*GeV );
00784 for( G4int j=0; j<20; ++j )binl[j] = (j-1.)/(19.*pt);
00785 et = pseudoParticle[1].GetTotalEnergy()/GeV;
00786 dndl[0] = 0.0;
00787 outerCounter = 0;
00788 eliminateThisParticle = true;
00789 resetEnergies = true;
00790 while( ++outerCounter < 3 )
00791 {
00792 for( l=1; l<20; ++l )
00793 {
00794 x = (binl[l]+binl[l-1])/2.;
00795 if( x > 1.0/pt )
00796 dndl[l] += dndl[l-1];
00797 else
00798 dndl[l] = aspar/std::sqrt( std::pow((1.+aspar*x*aspar*x),3) ) *
00799 (binl[l]-binl[l-1])*et / std::sqrt( pt*x*et*pt*x*et+pt*pt+vecMass*vecMass ) +
00800 dndl[l-1];
00801 }
00802 innerCounter = 0;
00803 while( ++innerCounter < 7 )
00804 {
00805 l = 1;
00806 ran = G4UniformRand()*dndl[19];
00807 while( ( ran >= dndl[l] ) && ( l < 20 ) )l++;
00808 l = std::min( 19, l );
00809 x = std::min( 1.0, pt*(binl[l-1] + G4UniformRand()*(binl[l]-binl[l-1]) ) );
00810 if( targetParticle.GetSide() < 0 )x *= -1.;
00811 targetParticle.SetMomentum( x*et*GeV );
00812 totalEnergy = std::sqrt( x*et*x*et + pt*pt + vecMass*vecMass );
00813 targetParticle.SetTotalEnergy( totalEnergy*GeV );
00814 if( targetParticle.GetSide() < 0 )
00815 {
00816 if( extraNucleonCount > 19 )x=0.999;
00817 G4double xxx = 0.95+0.05*extraNucleonCount/20.0;
00818 if( (backwardKinetic+totalEnergy-vecMass) < xxx*backwardEnergy )
00819 {
00820 pseudoParticle[5] = pseudoParticle[5] + targetParticle;
00821 backwardKinetic += totalEnergy - vecMass;
00822 pseudoParticle[6] = pseudoParticle[4] + pseudoParticle[5];
00823 pseudoParticle[6].SetMomentum( 0.0 );
00824 phi = pseudoParticle[6].Angle( pseudoParticle[8] );
00825 if( pseudoParticle[6].GetMomentum().y()/MeV < 0.0 )phi = twopi - phi;
00826 phi += pi + normal() * pi / 12.0;
00827 if( phi > twopi )phi -= twopi;
00828 if( phi < 0.0 )phi = twopi - phi;
00829 outerCounter = 2;
00830 eliminateThisParticle = false;
00831 resetEnergies = false;
00832 break;
00833 }
00834 if( innerCounter > 5 )break;
00835 if( forwardEnergy >= vecMass )
00836 {
00837 targetParticle.SetSide( 1 );
00838 forwardEnergy -= vecMass;
00839 backwardEnergy += vecMass;
00840 --backwardCount;
00841 }
00842 G4ThreeVector momentum = targetParticle.GetMomentum();
00843 targetParticle.SetMomentum( momentum.x() * 0.9, momentum.y() * 0.9 );
00844 pt *= 0.9;
00845 dndl[19] *= 0.9;
00846 }
00847 else
00848 {
00849 if( forwardEnergy < forwardKinetic )
00850 totalEnergy = vecMass + 0.04*std::fabs(normal());
00851 else
00852 totalEnergy = vecMass + forwardEnergy - forwardKinetic;
00853 targetParticle.SetTotalEnergy( totalEnergy*GeV );
00854 pp = std::sqrt( std::abs( totalEnergy*totalEnergy - vecMass*vecMass ) )*GeV;
00855 pp1 = targetParticle.GetMomentum().mag()/MeV;
00856 if( pp1 < 1.0e-6*GeV )
00857 {
00858 G4double costheta = 2.*G4UniformRand() - 1.;
00859 G4double sintheta = std::sqrt(1. - costheta*costheta);
00860 G4double phi2 = twopi*G4UniformRand();
00861 targetParticle.SetMomentum( pp*sintheta*std::cos(phi2)*MeV,
00862 pp*sintheta*std::sin(phi2)*MeV,
00863 pp*costheta*MeV ) ;
00864 }
00865 else
00866 targetParticle.SetMomentum( targetParticle.GetMomentum() * (pp/pp1) );
00867
00868 pseudoParticle[4] = pseudoParticle[4] + targetParticle;
00869 outerCounter = 2;
00870 eliminateThisParticle = false;
00871 resetEnergies = false;
00872 break;
00873 }
00874 }
00875 if( resetEnergies )
00876 {
00877
00878
00879
00880
00881
00882 forwardKinetic = backwardKinetic = 0.0;
00883 pseudoParticle[4].SetZero();
00884 pseudoParticle[5].SetZero();
00885 for( l=0; l<vecLen; ++l )
00886 {
00887 if (vec[l]->GetSide() > 0 ||
00888 vec[l]->GetDefinition()->GetParticleSubType() == "kaon" ||
00889 vec[l]->GetDefinition()->GetParticleSubType() == "pi") {
00890 G4double tempMass = vec[l]->GetMass()/GeV;
00891 totalEnergy =
00892 std::max( tempMass, 0.95*vec[l]->GetTotalEnergy()/GeV + 0.05*tempMass );
00893 vec[l]->SetTotalEnergy( totalEnergy*GeV );
00894 pp = std::sqrt( std::abs( totalEnergy*totalEnergy - tempMass*tempMass ) )*GeV;
00895 pp1 = vec[l]->GetMomentum().mag()/MeV;
00896 if( pp1 < 1.0e-6*GeV )
00897 {
00898 G4double costheta = 2.*G4UniformRand() - 1.;
00899 G4double sintheta = std::sqrt(1. - costheta*costheta);
00900 G4double phi2 = twopi*G4UniformRand();
00901 vec[l]->SetMomentum( pp*sintheta*std::cos(phi2)*MeV,
00902 pp*sintheta*std::sin(phi2)*MeV,
00903 pp*costheta*MeV ) ;
00904 }
00905 else
00906 vec[l]->SetMomentum( vec[l]->GetMomentum() * (pp/pp1) );
00907
00908 pt = std::max( 0.001*GeV, std::sqrt( sqr(vec[l]->GetMomentum().x()/MeV) +
00909 sqr(vec[l]->GetMomentum().y()/MeV) ) )/GeV;
00910 if( vec[l]->GetSide() > 0)
00911 {
00912 forwardKinetic += vec[l]->GetKineticEnergy()/GeV;
00913 pseudoParticle[4] = pseudoParticle[4] + (*vec[l]);
00914 } else {
00915 backwardKinetic += vec[l]->GetKineticEnergy()/GeV;
00916 pseudoParticle[5] = pseudoParticle[5] + (*vec[l]);
00917 }
00918 }
00919 }
00920 }
00921 }
00922 }
00923
00924
00925
00926
00927
00928
00929 pseudoParticle[6].Lorentz( pseudoParticle[3], pseudoParticle[2] );
00930 pseudoParticle[6] = pseudoParticle[6] - pseudoParticle[4];
00931 pseudoParticle[6] = pseudoParticle[6] - pseudoParticle[5];
00932 if( backwardNucleonCount == 1 )
00933 {
00934 G4double ekin =
00935 std::min( backwardEnergy-backwardKinetic, centerofmassEnergy/2.0-protonMass/GeV );
00936
00937 if( ekin < 0.04 )ekin = 0.04 * std::fabs( normal() );
00938 vecMass = targetParticle.GetMass()/GeV;
00939 totalEnergy = ekin+vecMass;
00940 targetParticle.SetTotalEnergy( totalEnergy*GeV );
00941 pp = std::sqrt( std::abs( totalEnergy*totalEnergy - vecMass*vecMass ) )*GeV;
00942 pp1 = pseudoParticle[6].GetMomentum().mag()/MeV;
00943 if( pp1 < 1.0e-6*GeV )
00944 {
00945 G4double costheta = 2.*G4UniformRand() - 1.;
00946 G4double sintheta = std::sqrt(1. - costheta*costheta);
00947 G4double phi2 = twopi*G4UniformRand();
00948 targetParticle.SetMomentum( pp*sintheta*std::cos(phi2)*MeV,
00949 pp*sintheta*std::sin(phi2)*MeV,
00950 pp*costheta*MeV ) ;
00951 } else {
00952 targetParticle.SetMomentum( pseudoParticle[6].GetMomentum() * (pp/pp1) );
00953 }
00954 pseudoParticle[5] = pseudoParticle[5] + targetParticle;
00955 }
00956 else
00957 {
00958 const G4double cpar[] = { 0.6, 0.6, 0.35, 0.15, 0.10 };
00959 const G4double gpar[] = { 2.6, 2.6, 1.80, 1.30, 1.20 };
00960
00961
00962 G4int tempCount;
00963 if (backwardNucleonCount < 5)
00964 {
00965 tempCount = backwardNucleonCount;
00966 }
00967 else
00968 {
00969 tempCount = 5;
00970 }
00971 tempCount--;
00972
00973
00974 G4double rmb0 = 0.0;
00975 if( targetParticle.GetSide() == -3 )
00976 rmb0 += targetParticle.GetMass()/GeV;
00977 for( i=0; i<vecLen; ++i )
00978 {
00979 if( vec[i]->GetSide() == -3 )rmb0 += vec[i]->GetMass()/GeV;
00980 }
00981 rmb = rmb0 + std::pow(-std::log(1.0-G4UniformRand()),cpar[tempCount]) / gpar[tempCount];
00982 totalEnergy = pseudoParticle[6].GetTotalEnergy()/GeV;
00983 vecMass = std::min( rmb, totalEnergy );
00984 pseudoParticle[6].SetMass( vecMass*GeV );
00985 pp = std::sqrt( std::abs( totalEnergy*totalEnergy - vecMass*vecMass ) )*GeV;
00986 pp1 = pseudoParticle[6].GetMomentum().mag()/MeV;
00987 if( pp1 < 1.0e-6*GeV )
00988 {
00989 G4double costheta = 2.*G4UniformRand() - 1.;
00990 G4double sintheta = std::sqrt(1. - costheta*costheta);
00991 G4double phi2 = twopi*G4UniformRand();
00992 pseudoParticle[6].SetMomentum( -pp*sintheta*std::cos(phi2)*MeV,
00993 -pp*sintheta*std::sin(phi2)*MeV,
00994 -pp*costheta*MeV ) ;
00995 }
00996 else
00997 pseudoParticle[6].SetMomentum( pseudoParticle[6].GetMomentum() * (-pp/pp1) );
00998
00999 G4FastVector<G4ReactionProduct,GHADLISTSIZE> tempV;
01000 tempV.Initialize( backwardNucleonCount );
01001 G4int tempLen = 0;
01002 if( targetParticle.GetSide() == -3 )tempV.SetElement( tempLen++, &targetParticle );
01003 for( i=0; i<vecLen; ++i )
01004 {
01005 if( vec[i]->GetSide() == -3 )tempV.SetElement( tempLen++, vec[i] );
01006 }
01007 if( tempLen != backwardNucleonCount )
01008 {
01009 G4cerr << "tempLen is not the same as backwardNucleonCount" << G4endl;
01010 G4cerr << "tempLen = " << tempLen;
01011 G4cerr << ", backwardNucleonCount = " << backwardNucleonCount << G4endl;
01012 G4cerr << "targetParticle side = " << targetParticle.GetSide() << G4endl;
01013 G4cerr << "currentParticle side = " << currentParticle.GetSide() << G4endl;
01014 for( i=0; i<vecLen; ++i )
01015 G4cerr << "particle #" << i << " side = " << vec[i]->GetSide() << G4endl;
01016 G4Exception("G4ReactionDynamics::GenerateXandPt", "601",
01017 FatalException, "Mismatch in nucleon count");
01018 }
01019 constantCrossSection = true;
01020
01021 if( tempLen >= 2 )
01022 {
01023 wgt = GenerateNBodyEvent(
01024 pseudoParticle[6].GetMass(), constantCrossSection, tempV, tempLen );
01025
01026 if( targetParticle.GetSide() == -3 )
01027 {
01028 targetParticle.Lorentz( targetParticle, pseudoParticle[6] );
01029
01030 pseudoParticle[5] = pseudoParticle[5] + targetParticle;
01031 }
01032 for( i=0; i<vecLen; ++i )
01033 {
01034 if( vec[i]->GetSide() == -3 )
01035 {
01036 vec[i]->Lorentz( *vec[i], pseudoParticle[6] );
01037 pseudoParticle[5] = pseudoParticle[5] + (*vec[i]);
01038 }
01039 }
01040
01041 }
01042 }
01043
01044
01045
01046 if( vecLen == 0 )return false;
01047
01048
01049 currentParticle.Lorentz( currentParticle, pseudoParticle[1] );
01050 targetParticle.Lorentz( targetParticle, pseudoParticle[1] );
01051 for( i=0; i<vecLen; ++i ) vec[i]->Lorentz( *vec[i], pseudoParticle[1] );
01052
01053
01054
01055
01056
01057
01058
01059
01060
01061
01062
01063 G4bool leadingStrangeParticleHasChanged = true;
01064 if( leadFlag )
01065 {
01066 if( currentParticle.GetDefinition() == leadingStrangeParticle.GetDefinition() )
01067 leadingStrangeParticleHasChanged = false;
01068 if( leadingStrangeParticleHasChanged &&
01069 ( targetParticle.GetDefinition() == leadingStrangeParticle.GetDefinition() ) )
01070 leadingStrangeParticleHasChanged = false;
01071 if( leadingStrangeParticleHasChanged )
01072 {
01073 for( i=0; i<vecLen; i++ )
01074 {
01075 if( vec[i]->GetDefinition() == leadingStrangeParticle.GetDefinition() )
01076 {
01077 leadingStrangeParticleHasChanged = false;
01078 break;
01079 }
01080 }
01081 }
01082 if( leadingStrangeParticleHasChanged )
01083 {
01084 G4bool leadTest =
01085 (leadingStrangeParticle.GetDefinition()->GetParticleSubType() == "kaon" ||
01086 leadingStrangeParticle.GetDefinition()->GetParticleSubType() == "pi");
01087 G4bool targetTest =
01088 (targetParticle.GetDefinition()->GetParticleSubType() == "kaon" ||
01089 targetParticle.GetDefinition()->GetParticleSubType() == "pi");
01090
01091
01092
01093 if( (leadTest&&targetTest) || !(leadTest||targetTest) )
01094 {
01095 targetParticle.SetDefinitionAndUpdateE( leadingStrangeParticle.GetDefinition() );
01096 targetHasChanged = true;
01097
01098 }
01099 else
01100 {
01101 currentParticle.SetDefinitionAndUpdateE( leadingStrangeParticle.GetDefinition() );
01102 incidentHasChanged = false;
01103
01104 }
01105 }
01106 }
01107
01108
01109
01110
01111 std::pair<G4int, G4int> finalStateNucleons =
01112 GetFinalStateNucleons(originalTarget, vec, vecLen);
01113
01114 G4int protonsInFinalState = finalStateNucleons.first;
01115 G4int neutronsInFinalState = finalStateNucleons.second;
01116
01117 G4int numberofFinalStateNucleons =
01118 protonsInFinalState + neutronsInFinalState;
01119
01120 if (currentParticle.GetDefinition()->GetBaryonNumber() == 1 &&
01121 targetParticle.GetDefinition()->GetBaryonNumber() == 1 &&
01122 originalIncident->GetDefinition()->GetPDGMass() <
01123 G4Lambda::Lambda()->GetPDGMass())
01124 numberofFinalStateNucleons++;
01125
01126 numberofFinalStateNucleons = std::max(1, numberofFinalStateNucleons);
01127
01128 G4int PinNucleus = std::max(0,
01129 targetNucleus.GetZ_asInt() - protonsInFinalState);
01130 G4int NinNucleus = std::max(0,
01131 targetNucleus.GetN_asInt() - neutronsInFinalState);
01132
01133 pseudoParticle[3].SetMomentum( 0.0, 0.0, pOriginal*GeV );
01134 pseudoParticle[3].SetMass( mOriginal*GeV );
01135 pseudoParticle[3].SetTotalEnergy(
01136 std::sqrt( pOriginal*pOriginal + mOriginal*mOriginal )*GeV );
01137
01138 G4ParticleDefinition * aOrgDef = modifiedOriginal.GetDefinition();
01139 G4int diff = 0;
01140 if(aOrgDef == G4Proton::Proton() || aOrgDef == G4Neutron::Neutron() ) diff = 1;
01141 if(numberofFinalStateNucleons == 1) diff = 0;
01142 pseudoParticle[4].SetMomentum( 0.0, 0.0, 0.0 );
01143 pseudoParticle[4].SetMass( protonMass*(numberofFinalStateNucleons-diff)*MeV );
01144 pseudoParticle[4].SetTotalEnergy( protonMass*(numberofFinalStateNucleons-diff)*MeV );
01145
01146 G4double theoreticalKinetic =
01147 pseudoParticle[3].GetTotalEnergy()/MeV +
01148 pseudoParticle[4].GetTotalEnergy()/MeV -
01149 currentParticle.GetMass()/MeV -
01150 targetParticle.GetMass()/MeV;
01151
01152 G4double simulatedKinetic =
01153 currentParticle.GetKineticEnergy()/MeV + targetParticle.GetKineticEnergy()/MeV;
01154
01155 pseudoParticle[5] = pseudoParticle[3] + pseudoParticle[4];
01156 pseudoParticle[3].Lorentz( pseudoParticle[3], pseudoParticle[5] );
01157 pseudoParticle[4].Lorentz( pseudoParticle[4], pseudoParticle[5] );
01158
01159 pseudoParticle[7].SetZero();
01160 pseudoParticle[7] = pseudoParticle[7] + currentParticle;
01161 pseudoParticle[7] = pseudoParticle[7] + targetParticle;
01162
01163 for( i=0; i<vecLen; ++i )
01164 {
01165 pseudoParticle[7] = pseudoParticle[7] + *vec[i];
01166 simulatedKinetic += vec[i]->GetKineticEnergy()/MeV;
01167 theoreticalKinetic -= vec[i]->GetMass()/MeV;
01168 }
01169
01170 if( vecLen <= 16 && vecLen > 0 )
01171 {
01172
01173
01174
01175 G4ReactionProduct tempR[130];
01176 tempR[0] = currentParticle;
01177 tempR[1] = targetParticle;
01178 for( i=0; i<vecLen; ++i )tempR[i+2] = *vec[i];
01179 G4FastVector<G4ReactionProduct,GHADLISTSIZE> tempV;
01180 tempV.Initialize( vecLen+2 );
01181 G4int tempLen = 0;
01182 for( i=0; i<vecLen+2; ++i )tempV.SetElement( tempLen++, &tempR[i] );
01183 constantCrossSection = true;
01184
01185 wgt = GenerateNBodyEvent( pseudoParticle[3].GetTotalEnergy()/MeV+
01186 pseudoParticle[4].GetTotalEnergy()/MeV,
01187 constantCrossSection, tempV, tempLen );
01188 if (wgt == -1) {
01189 G4double Qvalue = 0;
01190 for (i = 0; i < tempLen; i++) Qvalue += tempV[i]->GetMass();
01191 wgt = GenerateNBodyEvent( Qvalue/MeV,
01192 constantCrossSection, tempV, tempLen );
01193 }
01194 if(wgt>-.5)
01195 {
01196 theoreticalKinetic = 0.0;
01197 for( i=0; i<tempLen; ++i )
01198 {
01199 pseudoParticle[6].Lorentz( *tempV[i], pseudoParticle[4] );
01200 theoreticalKinetic += pseudoParticle[6].GetKineticEnergy()/MeV;
01201 }
01202 }
01203
01204 }
01205
01206
01207
01208 if( simulatedKinetic != 0.0 )
01209 {
01210 wgt = (theoreticalKinetic)/simulatedKinetic;
01211 theoreticalKinetic = currentParticle.GetKineticEnergy()/MeV * wgt;
01212 simulatedKinetic = theoreticalKinetic;
01213 currentParticle.SetKineticEnergy( theoreticalKinetic*MeV );
01214 pp = currentParticle.GetTotalMomentum()/MeV;
01215 pp1 = currentParticle.GetMomentum().mag()/MeV;
01216 if( pp1 < 1.0e-6*GeV )
01217 {
01218 G4double costheta = 2.*G4UniformRand() - 1.;
01219 G4double sintheta = std::sqrt(1. - costheta*costheta);
01220 G4double phi2 = twopi*G4UniformRand();
01221 currentParticle.SetMomentum( pp*sintheta*std::cos(phi2)*MeV,
01222 pp*sintheta*std::sin(phi2)*MeV,
01223 pp*costheta*MeV ) ;
01224 }
01225 else
01226 {
01227 currentParticle.SetMomentum( currentParticle.GetMomentum() * (pp/pp1) );
01228 }
01229 theoreticalKinetic = targetParticle.GetKineticEnergy()/MeV * wgt;
01230 targetParticle.SetKineticEnergy( theoreticalKinetic*MeV );
01231 simulatedKinetic += theoreticalKinetic;
01232 pp = targetParticle.GetTotalMomentum()/MeV;
01233 pp1 = targetParticle.GetMomentum().mag()/MeV;
01234
01235 if( pp1 < 1.0e-6*GeV )
01236 {
01237 G4double costheta = 2.*G4UniformRand() - 1.;
01238 G4double sintheta = std::sqrt(1. - costheta*costheta);
01239 G4double phi2 = twopi*G4UniformRand();
01240 targetParticle.SetMomentum( pp*sintheta*std::cos(phi2)*MeV,
01241 pp*sintheta*std::sin(phi2)*MeV,
01242 pp*costheta*MeV ) ;
01243 } else {
01244 targetParticle.SetMomentum( targetParticle.GetMomentum() * (pp/pp1) );
01245 }
01246 for( i=0; i<vecLen; ++i )
01247 {
01248 theoreticalKinetic = vec[i]->GetKineticEnergy()/MeV * wgt;
01249 simulatedKinetic += theoreticalKinetic;
01250 vec[i]->SetKineticEnergy( theoreticalKinetic*MeV );
01251 pp = vec[i]->GetTotalMomentum()/MeV;
01252 pp1 = vec[i]->GetMomentum().mag()/MeV;
01253 if( pp1 < 1.0e-6*GeV )
01254 {
01255 G4double costheta = 2.*G4UniformRand() - 1.;
01256 G4double sintheta = std::sqrt(1. - costheta*costheta);
01257 G4double phi2 = twopi*G4UniformRand();
01258 vec[i]->SetMomentum( pp*sintheta*std::cos(phi2)*MeV,
01259 pp*sintheta*std::sin(phi2)*MeV,
01260 pp*costheta*MeV ) ;
01261 }
01262 else
01263 vec[i]->SetMomentum( vec[i]->GetMomentum() * (pp/pp1) );
01264 }
01265 }
01266
01267
01268 Rotate( numberofFinalStateNucleons, pseudoParticle[3].GetMomentum(),
01269 modifiedOriginal, originalIncident, targetNucleus,
01270 currentParticle, targetParticle, vec, vecLen );
01271
01272
01273
01274
01275
01276
01277 if( atomicWeight >= 1.5 )
01278 {
01279
01280
01281
01282
01283
01284 G4int npnb = 0;
01285 G4int ndta = 0;
01286
01287 G4double epnb, edta;
01288 if (veryForward) {
01289 epnb = targetNucleus.GetAnnihilationPNBlackTrackEnergy();
01290 edta = targetNucleus.GetAnnihilationDTABlackTrackEnergy();
01291 } else {
01292 epnb = targetNucleus.GetPNBlackTrackEnergy();
01293 edta = targetNucleus.GetDTABlackTrackEnergy();
01294 }
01295
01296 const G4double pnCutOff = 0.001;
01297 const G4double dtaCutOff = 0.001;
01298 const G4double kineticMinimum = 1.e-6;
01299 const G4double kineticFactor = -0.010;
01300 G4double sprob = 0.0;
01301 const G4double ekIncident = originalIncident->GetKineticEnergy()/GeV;
01302 if( ekIncident >= 5.0 )sprob = std::min( 1.0, 0.6*std::log(ekIncident-4.0) );
01303 if( epnb >= pnCutOff )
01304 {
01305 npnb = Poisson((1.5+1.25*numberofFinalStateNucleons)*epnb/(epnb+edta));
01306 if( numberofFinalStateNucleons + npnb > atomicWeight )
01307 npnb = G4int(atomicWeight+0.00001 - numberofFinalStateNucleons);
01308 npnb = std::min( npnb, 127-vecLen );
01309 }
01310 if( edta >= dtaCutOff )
01311 {
01312 ndta = Poisson( (1.5+1.25*numberofFinalStateNucleons)*edta/(epnb+edta) );
01313 ndta = std::min( ndta, 127-vecLen );
01314 }
01315 if (npnb == 0 && ndta == 0) npnb = 1;
01316
01317
01318
01319 AddBlackTrackParticles(epnb, npnb, edta, ndta, sprob, kineticMinimum,
01320 kineticFactor, modifiedOriginal,
01321 PinNucleus, NinNucleus, targetNucleus,
01322 vec, vecLen);
01323
01324
01325 }
01326
01327
01328
01329
01330
01331 if( (atomicWeight >= 1.5) && (atomicWeight <= 230.0) && (ekOriginal <= 0.2) )
01332 currentParticle.SetTOF( 1.0-500.0*std::exp(-ekOriginal/0.04)*std::log(G4UniformRand()) );
01333 else
01334 currentParticle.SetTOF( 1.0 );
01335 return true;
01336
01337 }
01338
01339 void G4ReactionDynamics::SuppressChargedPions(
01340 G4FastVector<G4ReactionProduct,GHADLISTSIZE> &vec,
01341 G4int &vecLen,
01342 const G4ReactionProduct &modifiedOriginal,
01343 G4ReactionProduct ¤tParticle,
01344 G4ReactionProduct &targetParticle,
01345 const G4Nucleus &targetNucleus,
01346 G4bool &incidentHasChanged,
01347 G4bool &targetHasChanged )
01348 {
01349
01350
01351
01352
01353 G4double mOriginal = modifiedOriginal.GetMass()/GeV;
01354 G4double etOriginal = modifiedOriginal.GetTotalEnergy()/GeV;
01355 G4double targetMass = targetParticle.GetDefinition()->GetPDGMass()/GeV;
01356 G4double cmEnergy = std::sqrt( mOriginal*mOriginal + targetMass*targetMass +
01357 2.0*targetMass*etOriginal );
01358 G4double eAvailable = cmEnergy - mOriginal - targetMass;
01359 for (G4int i = 0; i < vecLen; i++) eAvailable -= vec[i]->GetMass()/GeV;
01360
01361 const G4double atomicWeight = G4double(targetNucleus.GetA_asInt());
01362 const G4double atomicNumber = G4double(targetNucleus.GetZ_asInt());
01363 const G4double pOriginal = modifiedOriginal.GetTotalMomentum()/GeV;
01364
01365 G4ParticleDefinition *aPiMinus = G4PionMinus::PionMinus();
01366 G4ParticleDefinition *aPiPlus = G4PionPlus::PionPlus();
01367 G4ParticleDefinition* aPiZero = G4PionZero::PionZero();
01368 G4ParticleDefinition *aProton = G4Proton::Proton();
01369 G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
01370 G4double piMass = aPiPlus->GetPDGMass()/GeV;
01371 G4double nucleonMass = aNeutron->GetPDGMass()/GeV;
01372
01373 const G4bool antiTest =
01374 modifiedOriginal.GetDefinition() != G4AntiProton::AntiProton() &&
01375 modifiedOriginal.GetDefinition() != G4AntiNeutron::AntiNeutron() &&
01376 modifiedOriginal.GetDefinition() != G4AntiLambda::AntiLambda() &&
01377 modifiedOriginal.GetDefinition() != G4AntiSigmaPlus::AntiSigmaPlus() &&
01378 modifiedOriginal.GetDefinition() != G4AntiSigmaMinus::AntiSigmaMinus() &&
01379 modifiedOriginal.GetDefinition() != G4AntiXiZero::AntiXiZero() &&
01380 modifiedOriginal.GetDefinition() != G4AntiXiMinus::AntiXiMinus() &&
01381 modifiedOriginal.GetDefinition() != G4AntiOmegaMinus::AntiOmegaMinus();
01382
01383 if( antiTest && (
01384 currentParticle.GetDefinition() == aPiPlus ||
01385 currentParticle.GetDefinition() == aPiZero ||
01386 currentParticle.GetDefinition() == aPiMinus ) &&
01387 ( G4UniformRand() <= (10.0-pOriginal)/6.0 ) &&
01388 ( G4UniformRand() <= atomicWeight/300.0 ) )
01389 {
01390 if (eAvailable > nucleonMass - piMass) {
01391 if( G4UniformRand() > atomicNumber/atomicWeight )
01392 currentParticle.SetDefinitionAndUpdateE( aNeutron );
01393 else
01394 currentParticle.SetDefinitionAndUpdateE( aProton );
01395 incidentHasChanged = true;
01396 }
01397 }
01398 if( antiTest && (
01399 targetParticle.GetDefinition() == aPiPlus ||
01400 targetParticle.GetDefinition() == aPiZero ||
01401 targetParticle.GetDefinition() == aPiMinus ) &&
01402 ( G4UniformRand() <= (10.0-pOriginal)/6.0 ) &&
01403 ( G4UniformRand() <= atomicWeight/300.0 ) )
01404 {
01405 if (eAvailable > nucleonMass - piMass) {
01406 if( G4UniformRand() > atomicNumber/atomicWeight )
01407 targetParticle.SetDefinitionAndUpdateE( aNeutron );
01408 else
01409 targetParticle.SetDefinitionAndUpdateE( aProton );
01410 targetHasChanged = true;
01411 }
01412 }
01413
01414 for( G4int i=0; i<vecLen; ++i )
01415 {
01416 if( antiTest && (
01417 vec[i]->GetDefinition() == aPiPlus ||
01418 vec[i]->GetDefinition() == aPiZero ||
01419 vec[i]->GetDefinition() == aPiMinus ) &&
01420 ( G4UniformRand() <= (10.0-pOriginal)/6.0 ) &&
01421 ( G4UniformRand() <= atomicWeight/300.0 ) )
01422 {
01423 if (eAvailable > nucleonMass - piMass) {
01424 if( G4UniformRand() > atomicNumber/atomicWeight )
01425 vec[i]->SetDefinitionAndUpdateE( aNeutron );
01426 else
01427 vec[i]->SetDefinitionAndUpdateE( aProton );
01428 }
01429 }
01430 }
01431
01432 }
01433
01434 G4bool G4ReactionDynamics::TwoCluster(
01435 G4FastVector<G4ReactionProduct,GHADLISTSIZE> &vec,
01436 G4int &vecLen,
01437 G4ReactionProduct &modifiedOriginal,
01438 const G4HadProjectile *originalIncident,
01439 G4ReactionProduct ¤tParticle,
01440 G4ReactionProduct &targetParticle,
01441 const G4DynamicParticle* originalTarget,
01442 const G4Nucleus &targetNucleus,
01443 G4bool &incidentHasChanged,
01444 G4bool &targetHasChanged,
01445 G4bool leadFlag,
01446 G4ReactionProduct &leadingStrangeParticle )
01447 {
01448
01449
01450
01451
01452
01453
01454
01455
01456
01457
01458
01459
01460
01461 G4int i;
01462 G4ParticleDefinition *aProton = G4Proton::Proton();
01463 G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
01464 G4ParticleDefinition *aPiPlus = G4PionPlus::PionPlus();
01465 G4ParticleDefinition *aPiMinus = G4PionMinus::PionMinus();
01466 G4ParticleDefinition *aPiZero = G4PionZero::PionZero();
01467 G4bool veryForward = false;
01468
01469 const G4double protonMass = aProton->GetPDGMass()/MeV;
01470 const G4double ekOriginal = modifiedOriginal.GetKineticEnergy()/GeV;
01471 const G4double etOriginal = modifiedOriginal.GetTotalEnergy()/GeV;
01472 const G4double mOriginal = modifiedOriginal.GetMass()/GeV;
01473 const G4double pOriginal = modifiedOriginal.GetMomentum().mag()/GeV;
01474 G4double targetMass = targetParticle.GetDefinition()->GetPDGMass()/GeV;
01475 G4double centerofmassEnergy = std::sqrt( mOriginal*mOriginal +
01476 targetMass*targetMass +
01477 2.0*targetMass*etOriginal );
01478 G4double currentMass = currentParticle.GetMass()/GeV;
01479 targetMass = targetParticle.GetMass()/GeV;
01480
01481 if( currentMass == 0.0 && targetMass == 0.0 )
01482 {
01483 G4double ek = currentParticle.GetKineticEnergy();
01484 G4ThreeVector mom = currentParticle.GetMomentum();
01485 currentParticle = *vec[0];
01486 targetParticle = *vec[1];
01487 for( i=0; i<(vecLen-2); ++i )*vec[i] = *vec[i+2];
01488 if(vecLen<2)
01489 {
01490 for(G4int j=0; j<vecLen; j++) delete vec[j];
01491 vecLen = 0;
01492 throw G4HadReentrentException(__FILE__, __LINE__,
01493 "G4ReactionDynamics::TwoCluster: Negative number of particles");
01494 }
01495 delete vec[vecLen-1];
01496 delete vec[vecLen-2];
01497 vecLen -= 2;
01498 currentMass = currentParticle.GetMass()/GeV;
01499 targetMass = targetParticle.GetMass()/GeV;
01500 incidentHasChanged = true;
01501 targetHasChanged = true;
01502 currentParticle.SetKineticEnergy( ek );
01503 currentParticle.SetMomentum( mom );
01504 veryForward = true;
01505 }
01506
01507 const G4double atomicWeight = G4double(targetNucleus.GetA_asInt());
01508 const G4double atomicNumber = G4double(targetNucleus.GetZ_asInt());
01509
01510
01511
01512
01513
01514 G4int forwardCount = 1;
01515 currentParticle.SetSide( 1 );
01516 G4double forwardMass = currentParticle.GetMass()/GeV;
01517 G4double cMass = forwardMass;
01518
01519
01520 G4int backwardCount = 1;
01521 G4int backwardNucleonCount = 1;
01522 targetParticle.SetSide( -1 );
01523 G4double backwardMass = targetParticle.GetMass()/GeV;
01524 G4double bMass = backwardMass;
01525
01526 for( i=0; i<vecLen; ++i )
01527 {
01528 if( vec[i]->GetSide() < 0 )vec[i]->SetSide( -1 );
01529
01530
01531 if( vec[i]->GetSide() == -1 )
01532 {
01533 ++backwardCount;
01534 backwardMass += vec[i]->GetMass()/GeV;
01535 }
01536 else
01537 {
01538 ++forwardCount;
01539 forwardMass += vec[i]->GetMass()/GeV;
01540 }
01541 }
01542
01543
01544
01545 G4double term1 = std::log(centerofmassEnergy*centerofmassEnergy);
01546 if(term1 < 0) term1 = 0.0001;
01547 const G4double afc = 0.312 + 0.2 * std::log(term1);
01548 G4double xtarg;
01549 if( centerofmassEnergy < 2.0+G4UniformRand() )
01550 xtarg = afc * (std::pow(atomicWeight,0.33)-1.0) * (2*backwardCount+vecLen+2)/2.0;
01551 else
01552 xtarg = afc * (std::pow(atomicWeight,0.33)-1.0) * (2*backwardCount);
01553 if( xtarg <= 0.0 )xtarg = 0.01;
01554 G4int nuclearExcitationCount = Poisson( xtarg );
01555 if(atomicWeight<1.0001) nuclearExcitationCount = 0;
01556 G4int extraNucleonCount = 0;
01557 G4double extraMass = 0.0;
01558 G4double extraNucleonMass = 0.0;
01559 if( nuclearExcitationCount > 0 )
01560 {
01561 G4int momentumBin = std::min( 4, G4int(pOriginal/3.0) );
01562 const G4double nucsup[] = { 1.0, 0.8, 0.6, 0.5, 0.4 };
01563
01564
01565
01566
01567
01568 for( i=0; i<nuclearExcitationCount; ++i )
01569 {
01570 G4ReactionProduct* pVec = new G4ReactionProduct();
01571 if( G4UniformRand() < nucsup[momentumBin] )
01572 {
01573 if( G4UniformRand() > 1.0-atomicNumber/atomicWeight )
01574 pVec->SetDefinition( aProton );
01575 else
01576 pVec->SetDefinition( aNeutron );
01577 ++backwardNucleonCount;
01578 ++extraNucleonCount;
01579 extraNucleonMass += pVec->GetMass()/GeV;
01580 }
01581 else
01582 {
01583 G4double ran = G4UniformRand();
01584 if( ran < 0.3181 )
01585 pVec->SetDefinition( aPiPlus );
01586 else if( ran < 0.6819 )
01587 pVec->SetDefinition( aPiZero );
01588 else
01589 pVec->SetDefinition( aPiMinus );
01590 }
01591 pVec->SetSide( -2 );
01592 extraMass += pVec->GetMass()/GeV;
01593 pVec->SetNewlyAdded( true );
01594 vec.SetElement( vecLen++, pVec );
01595
01596 }
01597 }
01598
01599 G4double forwardEnergy = (centerofmassEnergy-cMass-bMass)/2.0 +cMass - forwardMass;
01600 G4double backwardEnergy = (centerofmassEnergy-cMass-bMass)/2.0 +bMass - backwardMass;
01601 G4double eAvailable = centerofmassEnergy - (forwardMass+backwardMass);
01602 G4bool secondaryDeleted;
01603 G4double pMass;
01604
01605 while( eAvailable <= 0.0 )
01606 {
01607 secondaryDeleted = false;
01608 for( i=(vecLen-1); i>=0; --i )
01609 {
01610 if( vec[i]->GetSide() == 1 && vec[i]->GetMayBeKilled())
01611 {
01612 pMass = vec[i]->GetMass()/GeV;
01613 for( G4int j=i; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];
01614 --forwardCount;
01615 forwardEnergy += pMass;
01616 forwardMass -= pMass;
01617 secondaryDeleted = true;
01618 break;
01619 }
01620 else if( vec[i]->GetSide() == -1 && vec[i]->GetMayBeKilled())
01621 {
01622 pMass = vec[i]->GetMass()/GeV;
01623 for( G4int j=i; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];
01624 --backwardCount;
01625 backwardEnergy += pMass;
01626 backwardMass -= pMass;
01627 secondaryDeleted = true;
01628 break;
01629 }
01630 }
01631 if( secondaryDeleted )
01632 {
01633 delete vec[vecLen-1];
01634 --vecLen;
01635
01636 }
01637 else
01638 {
01639 if( vecLen == 0 )
01640 {
01641 return false;
01642 }
01643 if( targetParticle.GetSide() == -1 )
01644 {
01645 pMass = targetParticle.GetMass()/GeV;
01646 targetParticle = *vec[0];
01647 for( G4int j=0; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];
01648 --backwardCount;
01649 backwardEnergy += pMass;
01650 backwardMass -= pMass;
01651 secondaryDeleted = true;
01652 }
01653 else if( targetParticle.GetSide() == 1 )
01654 {
01655 pMass = targetParticle.GetMass()/GeV;
01656 targetParticle = *vec[0];
01657 for( G4int j=0; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];
01658 --forwardCount;
01659 forwardEnergy += pMass;
01660 forwardMass -= pMass;
01661 secondaryDeleted = true;
01662 }
01663 if( secondaryDeleted )
01664 {
01665 delete vec[vecLen-1];
01666 --vecLen;
01667
01668 }
01669 else
01670 {
01671 if( currentParticle.GetSide() == -1 )
01672 {
01673 pMass = currentParticle.GetMass()/GeV;
01674 currentParticle = *vec[0];
01675 for( G4int j=0; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];
01676 --backwardCount;
01677 backwardEnergy += pMass;
01678 backwardMass -= pMass;
01679 secondaryDeleted = true;
01680 }
01681 else if( currentParticle.GetSide() == 1 )
01682 {
01683 pMass = currentParticle.GetMass()/GeV;
01684 currentParticle = *vec[0];
01685 for( G4int j=0; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];
01686 --forwardCount;
01687 forwardEnergy += pMass;
01688 forwardMass -= pMass;
01689 secondaryDeleted = true;
01690 }
01691 if( secondaryDeleted )
01692 {
01693 delete vec[vecLen-1];
01694 --vecLen;
01695
01696 }
01697 else break;
01698 }
01699 }
01700 eAvailable = centerofmassEnergy - (forwardMass+backwardMass);
01701 }
01702
01703
01704
01705
01706
01707
01708
01709 G4double rmc = 0.0, rmd = 0.0;
01710 const G4double cpar[] = { 0.6, 0.6, 0.35, 0.15, 0.10 };
01711 const G4double gpar[] = { 2.6, 2.6, 1.8, 1.30, 1.20 };
01712
01713 if (forwardCount <= 0 || backwardCount <= 0) return false;
01714
01715 if (forwardCount == 1) rmc = forwardMass;
01716 else
01717 {
01718 G4int ntc = std::max(1, std::min(5,forwardCount))-1;
01719 rmc = forwardMass + std::pow(-std::log(1.0-G4UniformRand()),cpar[ntc-1])/gpar[ntc-1];
01720 }
01721
01722 if (backwardCount == 1) rmd = backwardMass;
01723 else
01724 {
01725 G4int ntc = std::max(1, std::min(5,backwardCount));
01726 rmd = backwardMass + std::pow(-std::log(1.0-G4UniformRand()),cpar[ntc-1])/gpar[ntc-1];
01727 }
01728
01729 while( rmc+rmd > centerofmassEnergy )
01730 {
01731 if( (rmc <= forwardMass) && (rmd <= backwardMass) )
01732 {
01733 G4double temp = 0.999*centerofmassEnergy/(rmc+rmd);
01734 rmc *= temp;
01735 rmd *= temp;
01736 }
01737 else
01738 {
01739 rmc = 0.1*forwardMass + 0.9*rmc;
01740 rmd = 0.1*backwardMass + 0.9*rmd;
01741 }
01742 }
01743
01744 G4ReactionProduct pseudoParticle[8];
01745 for( i=0; i<8; ++i )pseudoParticle[i].SetZero();
01746
01747 pseudoParticle[1].SetMass( mOriginal*GeV );
01748 pseudoParticle[1].SetTotalEnergy( etOriginal*GeV );
01749 pseudoParticle[1].SetMomentum( 0.0, 0.0, pOriginal*GeV );
01750
01751 pseudoParticle[2].SetMass( protonMass*MeV );
01752 pseudoParticle[2].SetTotalEnergy( protonMass*MeV );
01753 pseudoParticle[2].SetMomentum( 0.0, 0.0, 0.0 );
01754
01755
01756
01757 pseudoParticle[0] = pseudoParticle[1] + pseudoParticle[2];
01758 pseudoParticle[1].Lorentz( pseudoParticle[1], pseudoParticle[0] );
01759 pseudoParticle[2].Lorentz( pseudoParticle[2], pseudoParticle[0] );
01760
01761 const G4double pfMin = 0.0001;
01762 G4double pf = (centerofmassEnergy*centerofmassEnergy+rmd*rmd-rmc*rmc);
01763 pf *= pf;
01764 pf -= 4*centerofmassEnergy*centerofmassEnergy*rmd*rmd;
01765 pf = std::sqrt( std::max(pf,pfMin) )/(2.0*centerofmassEnergy);
01766
01767
01768
01769 pseudoParticle[3].SetMass( rmc*GeV );
01770 pseudoParticle[3].SetTotalEnergy( std::sqrt(pf*pf+rmc*rmc)*GeV );
01771
01772 pseudoParticle[4].SetMass( rmd*GeV );
01773 pseudoParticle[4].SetTotalEnergy( std::sqrt(pf*pf+rmd*rmd)*GeV );
01774
01775
01776
01777 const G4double bMin = 0.01;
01778 const G4double b1 = 4.0;
01779 const G4double b2 = 1.6;
01780
01781 G4double pin = pseudoParticle[1].GetMomentum().mag()/GeV;
01782 G4double dtb = 4.0*pin*pf*std::max( bMin, b1+b2*std::log(pOriginal) );
01783 G4double factor = 1.0 - std::exp(-dtb);
01784 G4double costheta = 1.0 + 2.0*std::log(1.0 - G4UniformRand()*factor) / dtb;
01785 costheta = std::max(-1.0, std::min(1.0, costheta));
01786 G4double sintheta = std::sqrt(1.0-costheta*costheta);
01787 G4double phi = G4UniformRand() * twopi;
01788
01789
01790
01791
01792 pseudoParticle[3].SetMomentum( pf*sintheta*std::cos(phi)*GeV,
01793 pf*sintheta*std::sin(phi)*GeV,
01794 pf*costheta*GeV );
01795
01796 pseudoParticle[4].SetMomentum( pseudoParticle[3].GetMomentum() * (-1.0) );
01797
01798
01799
01800 G4double pp, pp1;
01801 if( nuclearExcitationCount > 0 )
01802 {
01803 const G4double ga = 1.2;
01804 G4double ekit1 = 0.04;
01805 G4double ekit2 = 0.6;
01806 if( ekOriginal <= 5.0 )
01807 {
01808 ekit1 *= ekOriginal*ekOriginal/25.0;
01809 ekit2 *= ekOriginal*ekOriginal/25.0;
01810 }
01811 const G4double a = (1.0-ga)/(std::pow(ekit2,(1.0-ga)) - std::pow(ekit1,(1.0-ga)));
01812 for( i=0; i<vecLen; ++i )
01813 {
01814 if( vec[i]->GetSide() == -2 )
01815 {
01816 G4double kineticE =
01817 std::pow( (G4UniformRand()*(1.0-ga)/a+std::pow(ekit1,(1.0-ga))), (1.0/(1.0-ga)) );
01818 vec[i]->SetKineticEnergy( kineticE*GeV );
01819 G4double vMass = vec[i]->GetMass()/MeV;
01820 G4double totalE = kineticE*GeV + vMass;
01821 pp = std::sqrt( std::abs(totalE*totalE-vMass*vMass) );
01822 G4double cost = std::min( 1.0, std::max( -1.0, std::log(2.23*G4UniformRand()+0.383)/0.96 ) );
01823 G4double sint = std::sqrt( std::max( 0.0, (1.0-cost*cost) ) );
01824 phi = twopi*G4UniformRand();
01825 vec[i]->SetMomentum( pp*sint*std::sin(phi)*MeV,
01826 pp*sint*std::cos(phi)*MeV,
01827 pp*cost*MeV );
01828 vec[i]->Lorentz( *vec[i], pseudoParticle[0] );
01829 }
01830 }
01831
01832 }
01833
01834
01835
01836 currentParticle.SetMomentum( pseudoParticle[3].GetMomentum() );
01837 currentParticle.SetTotalEnergy( pseudoParticle[3].GetTotalEnergy() );
01838
01839 targetParticle.SetMomentum( pseudoParticle[4].GetMomentum() );
01840 targetParticle.SetTotalEnergy( pseudoParticle[4].GetTotalEnergy() );
01841
01842 pseudoParticle[5].SetMomentum( pseudoParticle[3].GetMomentum() * (-1.0) );
01843 pseudoParticle[5].SetMass( pseudoParticle[3].GetMass() );
01844 pseudoParticle[5].SetTotalEnergy( pseudoParticle[3].GetTotalEnergy() );
01845
01846 pseudoParticle[6].SetMomentum( pseudoParticle[4].GetMomentum() * (-1.0) );
01847 pseudoParticle[6].SetMass( pseudoParticle[4].GetMass() );
01848 pseudoParticle[6].SetTotalEnergy( pseudoParticle[4].GetTotalEnergy() );
01849
01850 G4double wgt;
01851
01852 if( forwardCount > 1 )
01853 {
01854 G4FastVector<G4ReactionProduct,GHADLISTSIZE> tempV;
01855 tempV.Initialize( forwardCount );
01856 G4bool constantCrossSection = true;
01857 G4int tempLen = 0;
01858 if( currentParticle.GetSide() == 1 )
01859 tempV.SetElement( tempLen++, ¤tParticle );
01860 if( targetParticle.GetSide() == 1 )
01861 tempV.SetElement( tempLen++, &targetParticle );
01862 for( i=0; i<vecLen; ++i )
01863 {
01864 if( vec[i]->GetSide() == 1 )
01865 {
01866 if( tempLen < 18 )
01867 tempV.SetElement( tempLen++, vec[i] );
01868 else
01869 {
01870 vec[i]->SetSide( -1 );
01871 continue;
01872 }
01873 }
01874 }
01875 if( tempLen >= 2 )
01876 {
01877 wgt = GenerateNBodyEvent( pseudoParticle[3].GetMass()/MeV,
01878 constantCrossSection, tempV, tempLen );
01879 if( currentParticle.GetSide() == 1 )
01880 currentParticle.Lorentz( currentParticle, pseudoParticle[5] );
01881 if( targetParticle.GetSide() == 1 )
01882 targetParticle.Lorentz( targetParticle, pseudoParticle[5] );
01883 for( i=0; i<vecLen; ++i )
01884 {
01885 if( vec[i]->GetSide() == 1 )vec[i]->Lorentz( *vec[i], pseudoParticle[5] );
01886 }
01887 }
01888 }
01889
01890 if( backwardCount > 1 )
01891 {
01892 G4FastVector<G4ReactionProduct,GHADLISTSIZE> tempV;
01893 tempV.Initialize( backwardCount );
01894 G4bool constantCrossSection = true;
01895 G4int tempLen = 0;
01896 if( currentParticle.GetSide() == -1 )
01897 tempV.SetElement( tempLen++, ¤tParticle );
01898 if( targetParticle.GetSide() == -1 )
01899 tempV.SetElement( tempLen++, &targetParticle );
01900 for( i=0; i<vecLen; ++i )
01901 {
01902 if( vec[i]->GetSide() == -1 )
01903 {
01904 if( tempLen < 18 )
01905 tempV.SetElement( tempLen++, vec[i] );
01906 else
01907 {
01908 vec[i]->SetSide( -2 );
01909 vec[i]->SetKineticEnergy( 0.0 );
01910 vec[i]->SetMomentum( 0.0, 0.0, 0.0 );
01911 continue;
01912 }
01913 }
01914 }
01915 if( tempLen >= 2 )
01916 {
01917 wgt = GenerateNBodyEvent( pseudoParticle[4].GetMass()/MeV,
01918 constantCrossSection, tempV, tempLen );
01919 if( currentParticle.GetSide() == -1 )
01920 currentParticle.Lorentz( currentParticle, pseudoParticle[6] );
01921 if( targetParticle.GetSide() == -1 )
01922 targetParticle.Lorentz( targetParticle, pseudoParticle[6] );
01923 for( i=0; i<vecLen; ++i )
01924 {
01925 if( vec[i]->GetSide() == -1 )vec[i]->Lorentz( *vec[i], pseudoParticle[6] );
01926 }
01927 }
01928 }
01929
01930
01931
01932
01933 currentParticle.Lorentz( currentParticle, pseudoParticle[2] );
01934 targetParticle.Lorentz( targetParticle, pseudoParticle[2] );
01935 for( i=0; i<vecLen; ++i ) vec[i]->Lorentz( *vec[i], pseudoParticle[2] );
01936
01937
01938
01939
01940
01941 G4bool dum = true;
01942 if( leadFlag )
01943 {
01944
01945
01946
01947
01948
01949
01950
01951
01952 if( currentParticle.GetDefinition() == leadingStrangeParticle.GetDefinition() )
01953 dum = false;
01954 else if( targetParticle.GetDefinition() == leadingStrangeParticle.GetDefinition() )
01955 dum = false;
01956 else
01957 {
01958 for( i=0; i<vecLen; ++i )
01959 {
01960 if( vec[i]->GetDefinition() == leadingStrangeParticle.GetDefinition() )
01961 {
01962 dum = false;
01963 break;
01964 }
01965 }
01966 }
01967 if( dum )
01968 {
01969 G4double leadMass = leadingStrangeParticle.GetMass()/MeV;
01970 G4double ekin;
01971 if( ((leadMass < protonMass) && (targetParticle.GetMass()/MeV < protonMass)) ||
01972 ((leadMass >= protonMass) && (targetParticle.GetMass()/MeV >= protonMass)) )
01973 {
01974 ekin = targetParticle.GetKineticEnergy()/GeV;
01975 pp1 = targetParticle.GetMomentum().mag()/MeV;
01976 targetParticle.SetDefinition( leadingStrangeParticle.GetDefinition() );
01977 targetParticle.SetKineticEnergy( ekin*GeV );
01978 pp = targetParticle.GetTotalMomentum()/MeV;
01979 if( pp1 < 1.0e-3 )
01980 {
01981 G4double costheta2 = 2.*G4UniformRand() - 1.;
01982 G4double sintheta2 = std::sqrt(1. - costheta2*costheta2);
01983 G4double phi2 = twopi*G4UniformRand();
01984 targetParticle.SetMomentum( pp*sintheta2*std::cos(phi2)*MeV,
01985 pp*sintheta2*std::sin(phi2)*MeV,
01986 pp*costheta2*MeV ) ;
01987 }
01988 else
01989 targetParticle.SetMomentum( targetParticle.GetMomentum() * (pp/pp1) );
01990
01991 targetHasChanged = true;
01992 }
01993 else
01994 {
01995 ekin = currentParticle.GetKineticEnergy()/GeV;
01996 pp1 = currentParticle.GetMomentum().mag()/MeV;
01997 currentParticle.SetDefinition( leadingStrangeParticle.GetDefinition() );
01998 currentParticle.SetKineticEnergy( ekin*GeV );
01999 pp = currentParticle.GetTotalMomentum()/MeV;
02000 if( pp1 < 1.0e-3 )
02001 {
02002 G4double costheta2 = 2.*G4UniformRand() - 1.;
02003 G4double sintheta2 = std::sqrt(1. - costheta2*costheta2);
02004 G4double phi2 = twopi*G4UniformRand();
02005 currentParticle.SetMomentum( pp*sintheta2*std::cos(phi2)*MeV,
02006 pp*sintheta2*std::sin(phi2)*MeV,
02007 pp*costheta2*MeV ) ;
02008 }
02009 else
02010 currentParticle.SetMomentum( currentParticle.GetMomentum() * (pp/pp1) );
02011
02012 incidentHasChanged = true;
02013 }
02014 }
02015 }
02016
02017
02018
02019
02020 std::pair<G4int, G4int> finalStateNucleons =
02021 GetFinalStateNucleons(originalTarget, vec, vecLen);
02022
02023 G4int protonsInFinalState = finalStateNucleons.first;
02024 G4int neutronsInFinalState = finalStateNucleons.second;
02025
02026 G4int numberofFinalStateNucleons =
02027 protonsInFinalState + neutronsInFinalState;
02028
02029 if (currentParticle.GetDefinition()->GetBaryonNumber() == 1 &&
02030 targetParticle.GetDefinition()->GetBaryonNumber() == 1 &&
02031 originalIncident->GetDefinition()->GetPDGMass() <
02032 G4Lambda::Lambda()->GetPDGMass())
02033 numberofFinalStateNucleons++;
02034
02035 numberofFinalStateNucleons = std::max(1, numberofFinalStateNucleons);
02036
02037 G4int PinNucleus = std::max(0,
02038 targetNucleus.GetZ_asInt() - protonsInFinalState);
02039 G4int NinNucleus = std::max(0,
02040 targetNucleus.GetN_asInt() - neutronsInFinalState);
02041
02042
02043
02044
02045 pseudoParticle[4].SetMass( mOriginal*GeV );
02046 pseudoParticle[4].SetTotalEnergy( etOriginal*GeV );
02047 pseudoParticle[4].SetMomentum( 0.0, 0.0, pOriginal*GeV );
02048
02049 G4ParticleDefinition * aOrgDef = modifiedOriginal.GetDefinition();
02050 G4int diff = 0;
02051 if(aOrgDef == G4Proton::Proton() || aOrgDef == G4Neutron::Neutron() ) diff = 1;
02052 if(numberofFinalStateNucleons == 1) diff = 0;
02053 pseudoParticle[5].SetMomentum( 0.0, 0.0, 0.0 );
02054 pseudoParticle[5].SetMass( protonMass*(numberofFinalStateNucleons-diff)*MeV);
02055 pseudoParticle[5].SetTotalEnergy( protonMass*(numberofFinalStateNucleons-diff)*MeV);
02056
02057 G4double theoreticalKinetic =
02058 pseudoParticle[4].GetTotalEnergy()/GeV + pseudoParticle[5].GetTotalEnergy()/GeV;
02059
02060 pseudoParticle[6] = pseudoParticle[4] + pseudoParticle[5];
02061 pseudoParticle[4].Lorentz( pseudoParticle[4], pseudoParticle[6] );
02062 pseudoParticle[5].Lorentz( pseudoParticle[5], pseudoParticle[6] );
02063
02064 if( vecLen < 16 )
02065 {
02066 G4ReactionProduct tempR[130];
02067 tempR[0] = currentParticle;
02068 tempR[1] = targetParticle;
02069 for( i=0; i<vecLen; ++i )tempR[i+2] = *vec[i];
02070
02071 G4FastVector<G4ReactionProduct,GHADLISTSIZE> tempV;
02072 tempV.Initialize( vecLen+2 );
02073 G4bool constantCrossSection = true;
02074 G4int tempLen = 0;
02075 for( i=0; i<vecLen+2; ++i )tempV.SetElement( tempLen++, &tempR[i] );
02076
02077 if( tempLen >= 2 )
02078 {
02079
02080 wgt = GenerateNBodyEvent( pseudoParticle[4].GetTotalEnergy()/MeV +
02081 pseudoParticle[5].GetTotalEnergy()/MeV,
02082 constantCrossSection, tempV, tempLen );
02083 if (wgt == -1) {
02084 G4double Qvalue = 0;
02085 for (i = 0; i < tempLen; i++) Qvalue += tempV[i]->GetMass();
02086 wgt = GenerateNBodyEvent( Qvalue/MeV,
02087 constantCrossSection, tempV, tempLen );
02088 }
02089 theoreticalKinetic = 0.0;
02090 for( i=0; i<vecLen+2; ++i )
02091 {
02092 pseudoParticle[7].SetMomentum( tempV[i]->GetMomentum() );
02093 pseudoParticle[7].SetMass( tempV[i]->GetMass() );
02094 pseudoParticle[7].SetTotalEnergy( tempV[i]->GetTotalEnergy() );
02095 pseudoParticle[7].Lorentz( pseudoParticle[7], pseudoParticle[5] );
02096 theoreticalKinetic += pseudoParticle[7].GetKineticEnergy()/GeV;
02097 }
02098 }
02099
02100 }
02101 else
02102 {
02103 theoreticalKinetic -=
02104 ( currentParticle.GetMass()/GeV + targetParticle.GetMass()/GeV );
02105 for( i=0; i<vecLen; ++i )theoreticalKinetic -= vec[i]->GetMass()/GeV;
02106 }
02107 G4double simulatedKinetic =
02108 currentParticle.GetKineticEnergy()/GeV + targetParticle.GetKineticEnergy()/GeV;
02109 for( i=0; i<vecLen; ++i )simulatedKinetic += vec[i]->GetKineticEnergy()/GeV;
02110
02111
02112
02113
02114
02115 if( simulatedKinetic != 0.0 )
02116 {
02117 wgt = (theoreticalKinetic)/simulatedKinetic;
02118 currentParticle.SetKineticEnergy( wgt*currentParticle.GetKineticEnergy() );
02119 pp = currentParticle.GetTotalMomentum()/MeV;
02120 pp1 = currentParticle.GetMomentum().mag()/MeV;
02121 if( pp1 < 0.001*MeV )
02122 {
02123 G4double costheta2 = 2.*G4UniformRand() - 1.;
02124 G4double sintheta2 = std::sqrt(1. - costheta2*costheta2);
02125 G4double phi2 = twopi*G4UniformRand();
02126 currentParticle.SetMomentum( pp*sintheta2*std::cos(phi2)*MeV,
02127 pp*sintheta2*std::sin(phi2)*MeV,
02128 pp*costheta2*MeV ) ;
02129 }
02130 else
02131 currentParticle.SetMomentum( currentParticle.GetMomentum() * (pp/pp1) );
02132
02133 targetParticle.SetKineticEnergy( wgt*targetParticle.GetKineticEnergy() );
02134 pp = targetParticle.GetTotalMomentum()/MeV;
02135 pp1 = targetParticle.GetMomentum().mag()/MeV;
02136 if( pp1 < 0.001*MeV )
02137 {
02138 G4double costheta2 = 2.*G4UniformRand() - 1.;
02139 G4double sintheta2 = std::sqrt(1. - costheta2*costheta2);
02140 G4double phi2 = twopi*G4UniformRand();
02141 targetParticle.SetMomentum( pp*sintheta2*std::cos(phi2)*MeV,
02142 pp*sintheta2*std::sin(phi2)*MeV,
02143 pp*costheta2*MeV ) ;
02144 }
02145 else
02146 targetParticle.SetMomentum( targetParticle.GetMomentum() * (pp/pp1) );
02147
02148 for( i=0; i<vecLen; ++i )
02149 {
02150 vec[i]->SetKineticEnergy( wgt*vec[i]->GetKineticEnergy() );
02151 pp = vec[i]->GetTotalMomentum()/MeV;
02152 pp1 = vec[i]->GetMomentum().mag()/MeV;
02153 if( pp1 < 0.001 )
02154 {
02155 G4double costheta2 = 2.*G4UniformRand() - 1.;
02156 G4double sintheta2 = std::sqrt(1. - costheta2*costheta2);
02157 G4double phi2 = twopi*G4UniformRand();
02158 vec[i]->SetMomentum( pp*sintheta2*std::cos(phi2)*MeV,
02159 pp*sintheta2*std::sin(phi2)*MeV,
02160 pp*costheta2*MeV ) ;
02161 }
02162 else
02163 vec[i]->SetMomentum( vec[i]->GetMomentum() * (pp/pp1) );
02164 }
02165 }
02166
02167
02168 Rotate( numberofFinalStateNucleons, pseudoParticle[4].GetMomentum(),
02169 modifiedOriginal, originalIncident, targetNucleus,
02170 currentParticle, targetParticle, vec, vecLen );
02171
02172
02173
02174
02175
02176 if( atomicWeight >= 1.5 )
02177 {
02178
02179
02180
02181
02182
02183
02184
02185 G4int npnb = 0;
02186 G4int ndta = 0;
02187
02188 G4double epnb, edta;
02189 if (veryForward) {
02190 epnb = targetNucleus.GetAnnihilationPNBlackTrackEnergy();
02191 edta = targetNucleus.GetAnnihilationDTABlackTrackEnergy();
02192 } else {
02193 epnb = targetNucleus.GetPNBlackTrackEnergy();
02194 edta = targetNucleus.GetDTABlackTrackEnergy();
02195 }
02196
02197 const G4double pnCutOff = 0.001;
02198 const G4double dtaCutOff = 0.001;
02199 const G4double kineticMinimum = 1.e-6;
02200 const G4double kineticFactor = -0.005;
02201
02202 G4double sprob = 0.0;
02203
02204 const G4double ekIncident = originalIncident->GetKineticEnergy()/GeV;
02205 if( ekIncident >= 5.0 )sprob = std::min( 1.0, 0.6*std::log(ekIncident-4.0) );
02206
02207 if( epnb >= pnCutOff )
02208 {
02209 npnb = Poisson((1.5+1.25*numberofFinalStateNucleons)*epnb/(epnb+edta));
02210 if( numberofFinalStateNucleons + npnb > atomicWeight )
02211 npnb = G4int(atomicWeight - numberofFinalStateNucleons);
02212 npnb = std::min( npnb, 127-vecLen );
02213 }
02214 if( edta >= dtaCutOff )
02215 {
02216 ndta = Poisson( (1.5+1.25*numberofFinalStateNucleons)*edta/(epnb+edta) );
02217 ndta = std::min( ndta, 127-vecLen );
02218 }
02219 if (npnb == 0 && ndta == 0) npnb = 1;
02220
02221
02222
02223 AddBlackTrackParticles(epnb, npnb, edta, ndta, sprob, kineticMinimum,
02224 kineticFactor, modifiedOriginal,
02225 PinNucleus, NinNucleus, targetNucleus,
02226 vec, vecLen );
02227
02228 }
02229
02230
02231
02232
02233
02234 if( (atomicWeight >= 1.5) && (atomicWeight <= 230.0) && (ekOriginal <= 0.2) )
02235 currentParticle.SetTOF( 1.0-500.0*std::exp(-ekOriginal/0.04)*std::log(G4UniformRand()) );
02236 else
02237 currentParticle.SetTOF( 1.0 );
02238
02239
02240 return true;
02241 }
02242
02243 void G4ReactionDynamics::TwoBody(
02244 G4FastVector<G4ReactionProduct,GHADLISTSIZE> &vec,
02245 G4int &vecLen,
02246 G4ReactionProduct &modifiedOriginal,
02247 const G4DynamicParticle* originalTarget,
02248 G4ReactionProduct ¤tParticle,
02249 G4ReactionProduct &targetParticle,
02250 const G4Nucleus &targetNucleus,
02251 G4bool & )
02252 {
02253
02254
02255
02256
02257
02258
02259
02260
02261
02262
02263
02264 static const G4double expxu = 82.;
02265 static const G4double expxl = -expxu;
02266
02267 const G4double ekOriginal = modifiedOriginal.GetKineticEnergy()/GeV;
02268 const G4double etOriginal = modifiedOriginal.GetTotalEnergy()/GeV;
02269 const G4double mOriginal = modifiedOriginal.GetMass()/GeV;
02270 const G4double pOriginal = modifiedOriginal.GetMomentum().mag()/GeV;
02271 G4double currentMass = currentParticle.GetMass()/GeV;
02272 G4double targetMass = targetParticle.GetDefinition()->GetPDGMass()/GeV;
02273
02274 targetMass = targetParticle.GetMass()/GeV;
02275 const G4double atomicWeight = G4double(targetNucleus.GetA_asInt());
02276
02277 G4double etCurrent = currentParticle.GetTotalEnergy()/GeV;
02278 G4double pCurrent = currentParticle.GetTotalMomentum()/GeV;
02279
02280 G4double cmEnergy = std::sqrt( currentMass*currentMass +
02281 targetMass*targetMass +
02282 2.0*targetMass*etCurrent );
02283
02284
02285
02286
02287
02288
02289
02290 if( (pCurrent < 0.1) || (cmEnergy < 0.01) )
02291 {
02292 targetParticle.SetMass( 0.0 );
02293 }
02294 else
02295 {
02296
02297
02298
02299
02300
02301
02302
02303
02304
02305
02306
02307
02308
02309
02310 G4double pf = cmEnergy*cmEnergy + targetMass*targetMass - currentMass*currentMass;
02311 pf = pf*pf - 4*cmEnergy*cmEnergy*targetMass*targetMass;
02312
02313 if( pf < 0.001 )
02314 {
02315 for(G4int i=0; i<vecLen; i++) delete vec[i];
02316 vecLen = 0;
02317 throw G4HadronicException(__FILE__, __LINE__, "G4ReactionDynamics::TwoBody: pf is too small ");
02318 }
02319
02320 pf = std::sqrt( pf ) / ( 2.0*cmEnergy );
02321
02322
02323
02324 G4ReactionProduct pseudoParticle[3];
02325
02326 if (targetParticle.GetDefinition()->GetParticleSubType() == "kaon" ||
02327 targetParticle.GetDefinition()->GetParticleSubType() == "pi") {
02328 pseudoParticle[0].SetMass( targetMass*GeV );
02329 pseudoParticle[0].SetTotalEnergy( etOriginal*GeV );
02330 pseudoParticle[0].SetMomentum( 0.0, 0.0, pOriginal*GeV );
02331
02332 pseudoParticle[1].SetMomentum( 0.0, 0.0, 0.0 );
02333 pseudoParticle[1].SetMass( mOriginal*GeV );
02334 pseudoParticle[1].SetKineticEnergy( 0.0 );
02335
02336 } else {
02337 pseudoParticle[0].SetMass( currentMass*GeV );
02338 pseudoParticle[0].SetTotalEnergy( etCurrent*GeV );
02339 pseudoParticle[0].SetMomentum( 0.0, 0.0, pCurrent*GeV );
02340
02341 pseudoParticle[1].SetMomentum( 0.0, 0.0, 0.0 );
02342 pseudoParticle[1].SetMass( targetMass*GeV );
02343 pseudoParticle[1].SetKineticEnergy( 0.0 );
02344 }
02345
02346
02347
02348 pseudoParticle[2] = pseudoParticle[0] + pseudoParticle[1];
02349 pseudoParticle[0].Lorentz( pseudoParticle[0], pseudoParticle[2] );
02350 pseudoParticle[1].Lorentz( pseudoParticle[1], pseudoParticle[2] );
02351
02352
02353
02354 currentParticle.SetTotalEnergy( std::sqrt(pf*pf+currentMass*currentMass)*GeV );
02355 targetParticle.SetTotalEnergy( std::sqrt(pf*pf+targetMass*targetMass)*GeV );
02356
02357
02358
02359 const G4double cb = 0.01;
02360 const G4double b1 = 4.225;
02361 const G4double b2 = 1.795;
02362
02363
02364
02365 G4double b = std::max( cb, b1+b2*std::log(pOriginal) );
02366 G4double btrang = b * 4.0 * pf * pseudoParticle[0].GetMomentum().mag()/GeV;
02367
02368 G4double exindt = -1.0;
02369 exindt += std::exp(std::max(-btrang,expxl));
02370
02371
02372
02373 G4double ctet = 1.0 + 2*std::log( 1.0+G4UniformRand()*exindt ) / btrang;
02374 if( std::fabs(ctet) > 1.0 )ctet > 0.0 ? ctet = 1.0 : ctet = -1.0;
02375 G4double stet = std::sqrt( (1.0-ctet)*(1.0+ctet) );
02376 G4double phi = twopi * G4UniformRand();
02377
02378
02379
02380 if (targetParticle.GetDefinition()->GetParticleSubType() == "kaon" ||
02381 targetParticle.GetDefinition()->GetParticleSubType() == "pi") {
02382
02383 currentParticle.SetMomentum( -pf*stet*std::sin(phi)*GeV,
02384 -pf*stet*std::cos(phi)*GeV,
02385 -pf*ctet*GeV );
02386 } else {
02387
02388 currentParticle.SetMomentum( pf*stet*std::sin(phi)*GeV,
02389 pf*stet*std::cos(phi)*GeV,
02390 pf*ctet*GeV );
02391 }
02392 targetParticle.SetMomentum( currentParticle.GetMomentum() * (-1.0) );
02393
02394
02395
02396 currentParticle.Lorentz( currentParticle, pseudoParticle[1] );
02397 targetParticle.Lorentz( targetParticle, pseudoParticle[1] );
02398
02399 Defs1( modifiedOriginal, currentParticle, targetParticle, vec, vecLen );
02400
02401 G4double pp, pp1, ekin;
02402 if( atomicWeight >= 1.5 )
02403 {
02404 const G4double cfa = 0.025*((atomicWeight-1.)/120.)*std::exp(-(atomicWeight-1.)/120.);
02405 pp1 = currentParticle.GetMomentum().mag()/MeV;
02406 if( pp1 >= 1.0 )
02407 {
02408 ekin = currentParticle.GetKineticEnergy()/MeV - cfa*(1.0+0.5*normal())*GeV;
02409 ekin = std::max( 0.0001*GeV, ekin );
02410 currentParticle.SetKineticEnergy( ekin*MeV );
02411 pp = currentParticle.GetTotalMomentum()/MeV;
02412 currentParticle.SetMomentum( currentParticle.GetMomentum() * (pp/pp1) );
02413 }
02414 pp1 = targetParticle.GetMomentum().mag()/MeV;
02415 if( pp1 >= 1.0 )
02416 {
02417 ekin = targetParticle.GetKineticEnergy()/MeV - cfa*(1.0+normal()/2.)*GeV;
02418 ekin = std::max( 0.0001*GeV, ekin );
02419 targetParticle.SetKineticEnergy( ekin*MeV );
02420 pp = targetParticle.GetTotalMomentum()/MeV;
02421 targetParticle.SetMomentum( targetParticle.GetMomentum() * (pp/pp1) );
02422 }
02423 }
02424 }
02425
02426
02427
02428
02429 std::pair<G4int, G4int> finalStateNucleons =
02430 GetFinalStateNucleons(originalTarget, vec, vecLen);
02431 G4int protonsInFinalState = finalStateNucleons.first;
02432 G4int neutronsInFinalState = finalStateNucleons.second;
02433
02434 G4int PinNucleus = std::max(0,
02435 targetNucleus.GetZ_asInt() - protonsInFinalState);
02436 G4int NinNucleus = std::max(0,
02437 targetNucleus.GetN_asInt() - neutronsInFinalState);
02438
02439
02440 if( atomicWeight >= 1.5 )
02441 {
02442
02443
02444
02445
02446
02447
02448 G4double epnb, edta;
02449 G4int npnb=0, ndta=0;
02450
02451 epnb = targetNucleus.GetPNBlackTrackEnergy();
02452 edta = targetNucleus.GetDTABlackTrackEnergy();
02453 const G4double pnCutOff = 0.0001;
02454 const G4double dtaCutOff = 0.0001;
02455 const G4double kineticMinimum = 0.0001;
02456 const G4double kineticFactor = -0.010;
02457 G4double sprob = 0.0;
02458 if( epnb >= pnCutOff )
02459 {
02460 npnb = Poisson( epnb/0.02 );
02461 if( npnb > atomicWeight )npnb = G4int(atomicWeight);
02462 if( (epnb > pnCutOff) && (npnb <= 0) )npnb = 1;
02463 npnb = std::min( npnb, 127-vecLen );
02464 }
02465 if( edta >= dtaCutOff )
02466 {
02467 ndta = G4int(2.0 * std::log(atomicWeight));
02468 ndta = std::min( ndta, 127-vecLen );
02469 }
02470
02471 if (npnb == 0 && ndta == 0) npnb = 1;
02472
02473
02474
02475 AddBlackTrackParticles(epnb, npnb, edta, ndta, sprob, kineticMinimum,
02476 kineticFactor, modifiedOriginal,
02477 PinNucleus, NinNucleus, targetNucleus,
02478 vec, vecLen);
02479
02480 }
02481
02482
02483
02484 if( (atomicWeight >= 1.5) && (atomicWeight <= 230.0) && (ekOriginal <= 0.2) )
02485 currentParticle.SetTOF( 1.0-500.0*std::exp(-ekOriginal/0.04)*std::log(G4UniformRand()) );
02486 else
02487 currentParticle.SetTOF( 1.0 );
02488 return;
02489 }
02490
02491 G4double G4ReactionDynamics::GenerateNBodyEvent(
02492 const G4double totalEnergy,
02493 const G4bool constantCrossSection,
02494 G4FastVector<G4ReactionProduct,GHADLISTSIZE> &vec,
02495 G4int &vecLen )
02496 {
02497
02498
02499
02500
02501 G4int i;
02502 const G4double expxu = 82.;
02503 const G4double expxl = -expxu;
02504 if( vecLen < 2 )
02505 {
02506 G4cerr << "*** Error in G4ReactionDynamics::GenerateNBodyEvent" << G4endl;
02507 G4cerr << " number of particles < 2" << G4endl;
02508 G4cerr << "totalEnergy = " << totalEnergy << "MeV, vecLen = " << vecLen << G4endl;
02509 return -1.0;
02510 }
02511 G4double mass[18];
02512 G4double energy[18];
02513 G4double pcm[3][18];
02514
02515 G4double totalMass = 0.0;
02516 G4double extraMass = 0;
02517 G4double sm[18];
02518
02519 for( i=0; i<vecLen; ++i )
02520 {
02521 mass[i] = vec[i]->GetMass()/GeV;
02522 if(vec[i]->GetSide() == -2) extraMass+=vec[i]->GetMass()/GeV;
02523 vec[i]->SetMomentum( 0.0, 0.0, 0.0 );
02524 pcm[0][i] = 0.0;
02525 pcm[1][i] = 0.0;
02526 pcm[2][i] = 0.0;
02527 energy[i] = mass[i];
02528 totalMass += mass[i];
02529 sm[i] = totalMass;
02530 }
02531 G4double totalE = totalEnergy/GeV;
02532 if( totalMass > totalE )
02533 {
02534
02535
02536
02537 totalE = totalMass;
02538 return -1.0;
02539 }
02540 G4double kineticEnergy = totalE - totalMass;
02541 G4double emm[18];
02542
02543 emm[0] = mass[0];
02544 emm[vecLen-1] = totalE;
02545 if( vecLen > 2 )
02546 {
02547 G4double ran[18];
02548 for( i=0; i<vecLen; ++i )ran[i] = G4UniformRand();
02549 for( i=0; i<vecLen-2; ++i )
02550 {
02551 for( G4int j=vecLen-2; j>i; --j )
02552 {
02553 if( ran[i] > ran[j] )
02554 {
02555 G4double temp = ran[i];
02556 ran[i] = ran[j];
02557 ran[j] = temp;
02558 }
02559 }
02560 }
02561 for( i=1; i<vecLen-1; ++i )emm[i] = ran[i-1]*kineticEnergy + sm[i];
02562 }
02563
02564 G4bool lzero = true;
02565 G4double wtmax = 0.0;
02566 if( constantCrossSection )
02567 {
02568 G4double emmax = kineticEnergy + mass[0];
02569 G4double emmin = 0.0;
02570 for( i=1; i<vecLen; ++i )
02571 {
02572 emmin += mass[i-1];
02573 emmax += mass[i];
02574 G4double wtfc = 0.0;
02575 if( emmax*emmax > 0.0 )
02576 {
02577 G4double arg = emmax*emmax
02578 + (emmin*emmin-mass[i]*mass[i])*(emmin*emmin-mass[i]*mass[i])/(emmax*emmax)
02579 - 2.0*(emmin*emmin+mass[i]*mass[i]);
02580 if( arg > 0.0 )wtfc = 0.5*std::sqrt( arg );
02581 }
02582 if( wtfc == 0.0 )
02583 {
02584 lzero = false;
02585 break;
02586 }
02587 wtmax += std::log( wtfc );
02588 }
02589 if( lzero )
02590 wtmax = -wtmax;
02591 else
02592 wtmax = expxu;
02593 }
02594 else
02595 {
02596
02597 const G4double ffq[18] = { 0., 3.141592, 19.73921, 62.01255, 129.8788, 204.0131,
02598 256.3704, 268.4705, 240.9780, 189.2637,
02599 132.1308, 83.0202, 47.4210, 24.8295,
02600 12.0006, 5.3858, 2.2560, 0.8859 };
02601 wtmax = std::log( std::pow( kineticEnergy, vecLen-2 ) * ffq[vecLen-1] / totalE );
02602 }
02603 lzero = true;
02604 G4double pd[50];
02605
02606 for( i=0; i<vecLen-1; ++i )
02607 {
02608 pd[i] = 0.0;
02609 if( emm[i+1]*emm[i+1] > 0.0 )
02610 {
02611 G4double arg = emm[i+1]*emm[i+1]
02612 + (emm[i]*emm[i]-mass[i+1]*mass[i+1])*(emm[i]*emm[i]-mass[i+1]*mass[i+1])
02613 /(emm[i+1]*emm[i+1])
02614 - 2.0*(emm[i]*emm[i]+mass[i+1]*mass[i+1]);
02615 if( arg > 0.0 )pd[i] = 0.5*std::sqrt( arg );
02616 }
02617 if( pd[i] <= 0.0 )
02618 lzero = false;
02619 else
02620 wtmax += std::log( pd[i] );
02621 }
02622 G4double weight = 0.0;
02623 if( lzero )weight = std::exp( std::max(std::min(wtmax,expxu),expxl) );
02624
02625 G4double bang, cb, sb, s0, s1, s2, c, ss, esys, a, b, gama, beta;
02626 pcm[0][0] = 0.0;
02627 pcm[1][0] = pd[0];
02628 pcm[2][0] = 0.0;
02629 for( i=1; i<vecLen; ++i )
02630 {
02631 pcm[0][i] = 0.0;
02632 pcm[1][i] = -pd[i-1];
02633 pcm[2][i] = 0.0;
02634 bang = twopi*G4UniformRand();
02635 cb = std::cos(bang);
02636 sb = std::sin(bang);
02637 c = 2.0*G4UniformRand() - 1.0;
02638 ss = std::sqrt( std::fabs( 1.0-c*c ) );
02639 if( i < vecLen-1 )
02640 {
02641 esys = std::sqrt(pd[i]*pd[i] + emm[i]*emm[i]);
02642 beta = pd[i]/esys;
02643 gama = esys/emm[i];
02644 for( G4int j=0; j<=i; ++j )
02645 {
02646 s0 = pcm[0][j];
02647 s1 = pcm[1][j];
02648 s2 = pcm[2][j];
02649 energy[j] = std::sqrt( s0*s0 + s1*s1 + s2*s2 + mass[j]*mass[j] );
02650 a = s0*c - s1*ss;
02651 pcm[1][j] = s0*ss + s1*c;
02652 b = pcm[2][j];
02653 pcm[0][j] = a*cb - b*sb;
02654 pcm[2][j] = a*sb + b*cb;
02655 pcm[1][j] = gama*(pcm[1][j] + beta*energy[j]);
02656 }
02657 }
02658 else
02659 {
02660 for( G4int j=0; j<=i; ++j )
02661 {
02662 s0 = pcm[0][j];
02663 s1 = pcm[1][j];
02664 s2 = pcm[2][j];
02665 energy[j] = std::sqrt( s0*s0 + s1*s1 + s2*s2 + mass[j]*mass[j] );
02666 a = s0*c - s1*ss;
02667 pcm[1][j] = s0*ss + s1*c;
02668 b = pcm[2][j];
02669 pcm[0][j] = a*cb - b*sb;
02670 pcm[2][j] = a*sb + b*cb;
02671 }
02672 }
02673 }
02674 for( i=0; i<vecLen; ++i )
02675 {
02676 vec[i]->SetMomentum( pcm[0][i]*GeV, pcm[1][i]*GeV, pcm[2][i]*GeV );
02677 vec[i]->SetTotalEnergy( energy[i]*GeV );
02678 }
02679
02680
02681 return weight;
02682 }
02683
02684 G4double
02685 G4ReactionDynamics::normal()
02686 {
02687 G4double ran = -6.0;
02688 for( G4int i=0; i<12; ++i )ran += G4UniformRand();
02689 return ran;
02690 }
02691
02692 G4int
02693 G4ReactionDynamics::Poisson( G4double x )
02694 {
02695 G4int iran;
02696 G4double ran;
02697
02698 if( x > 9.9 )
02699 iran = static_cast<G4int>(std::max( 0.0, x+normal()*std::sqrt(x) ) );
02700 else {
02701 G4int mn = G4int(5.0*x);
02702 if( mn <= 0 )
02703 {
02704 G4double p1 = x*std::exp(-x);
02705 G4double p2 = x*p1/2.0;
02706 G4double p3 = x*p2/3.0;
02707 ran = G4UniformRand();
02708 if( ran < p3 )
02709 iran = 3;
02710 else if( ran < p2 )
02711 iran = 2;
02712 else if( ran < p1 )
02713 iran = 1;
02714 else
02715 iran = 0;
02716 }
02717 else
02718 {
02719 iran = 0;
02720 G4double r = std::exp(-x);
02721 ran = G4UniformRand();
02722 if( ran > r )
02723 {
02724 G4double rrr;
02725 G4double rr = r;
02726 for( G4int i=1; i<=mn; ++i )
02727 {
02728 iran++;
02729 if( i > 5 )
02730 rrr = std::exp(i*std::log(x)-(i+0.5)*std::log((G4double)i)+i-0.9189385);
02731 else
02732 rrr = std::pow(x,i)/Factorial(i);
02733 rr += r*rrr;
02734 if( ran <= rr )break;
02735 }
02736 }
02737 }
02738 }
02739 return iran;
02740 }
02741
02742 G4int
02743 G4ReactionDynamics::Factorial( G4int n )
02744 {
02745 G4int mn = std::min(n,10);
02746 G4int result = 1;
02747 if( mn <= 1 )return result;
02748 for( G4int i=2; i<=mn; ++i )result *= i;
02749 return result;
02750 }
02751
02752 void G4ReactionDynamics::Defs1(
02753 const G4ReactionProduct &modifiedOriginal,
02754 G4ReactionProduct ¤tParticle,
02755 G4ReactionProduct &targetParticle,
02756 G4FastVector<G4ReactionProduct,GHADLISTSIZE> &vec,
02757 G4int &vecLen )
02758 {
02759 const G4double pjx = modifiedOriginal.GetMomentum().x()/MeV;
02760 const G4double pjy = modifiedOriginal.GetMomentum().y()/MeV;
02761 const G4double pjz = modifiedOriginal.GetMomentum().z()/MeV;
02762 const G4double p = modifiedOriginal.GetMomentum().mag()/MeV;
02763 if( pjx*pjx+pjy*pjy > 0.0 )
02764 {
02765 G4double cost, sint, ph, cosp, sinp, pix, piy, piz;
02766 cost = pjz/p;
02767 sint = 0.5 * ( std::sqrt(std::abs((1.0-cost)*(1.0+cost))) + std::sqrt(pjx*pjx+pjy*pjy)/p );
02768 if( pjy < 0.0 )
02769 ph = 3*halfpi;
02770 else
02771 ph = halfpi;
02772 if( std::abs( pjx ) > 0.001*MeV )ph = std::atan2(pjy,pjx);
02773 cosp = std::cos(ph);
02774 sinp = std::sin(ph);
02775 pix = currentParticle.GetMomentum().x()/MeV;
02776 piy = currentParticle.GetMomentum().y()/MeV;
02777 piz = currentParticle.GetMomentum().z()/MeV;
02778 currentParticle.SetMomentum( cost*cosp*pix*MeV - sinp*piy+sint*cosp*piz*MeV,
02779 cost*sinp*pix*MeV + cosp*piy+sint*sinp*piz*MeV,
02780 -sint*pix*MeV + cost*piz*MeV );
02781 pix = targetParticle.GetMomentum().x()/MeV;
02782 piy = targetParticle.GetMomentum().y()/MeV;
02783 piz = targetParticle.GetMomentum().z()/MeV;
02784 targetParticle.SetMomentum( cost*cosp*pix*MeV - sinp*piy+sint*cosp*piz*MeV,
02785 cost*sinp*pix*MeV + cosp*piy+sint*sinp*piz*MeV,
02786 -sint*pix*MeV + cost*piz*MeV );
02787 for( G4int i=0; i<vecLen; ++i )
02788 {
02789 pix = vec[i]->GetMomentum().x()/MeV;
02790 piy = vec[i]->GetMomentum().y()/MeV;
02791 piz = vec[i]->GetMomentum().z()/MeV;
02792 vec[i]->SetMomentum( cost*cosp*pix*MeV - sinp*piy+sint*cosp*piz*MeV,
02793 cost*sinp*pix*MeV + cosp*piy+sint*sinp*piz*MeV,
02794 -sint*pix*MeV + cost*piz*MeV );
02795 }
02796 }
02797 else
02798 {
02799 if( pjz < 0.0 )
02800 {
02801 currentParticle.SetMomentum( -currentParticle.GetMomentum().z() );
02802 targetParticle.SetMomentum( -targetParticle.GetMomentum().z() );
02803 for( G4int i=0; i<vecLen; ++i )
02804 vec[i]->SetMomentum( -vec[i]->GetMomentum().z() );
02805 }
02806 }
02807 }
02808
02809 void G4ReactionDynamics::Rotate(
02810 const G4double numberofFinalStateNucleons,
02811 const G4ThreeVector &temp,
02812 const G4ReactionProduct &modifiedOriginal,
02813 const G4HadProjectile *originalIncident,
02814 const G4Nucleus &targetNucleus,
02815 G4ReactionProduct ¤tParticle,
02816 G4ReactionProduct &targetParticle,
02817 G4FastVector<G4ReactionProduct,GHADLISTSIZE> &vec,
02818 G4int &vecLen )
02819 {
02820
02821
02822
02823
02824
02825 const G4double atomicWeight = G4double(targetNucleus.GetA_asInt());
02826 const G4double logWeight = std::log(atomicWeight);
02827
02828 G4ParticleDefinition *aPiMinus = G4PionMinus::PionMinus();
02829 G4ParticleDefinition *aPiPlus = G4PionPlus::PionPlus();
02830 G4ParticleDefinition *aPiZero = G4PionZero::PionZero();
02831
02832 G4int i;
02833 G4ThreeVector pseudoParticle[4];
02834 for( i=0; i<4; ++i )pseudoParticle[i].set(0,0,0);
02835 pseudoParticle[0] = currentParticle.GetMomentum()
02836 + targetParticle.GetMomentum();
02837 for( i=0; i<vecLen; ++i )
02838 pseudoParticle[0] = pseudoParticle[0] + (vec[i]->GetMomentum());
02839
02840
02841
02842 G4float pp, pp1;
02843 G4double alekw, p;
02844 G4double r1, r2, a1, ran1, ran2, xxh, exh, pxTemp, pyTemp, pzTemp;
02845
02846 r1 = twopi*G4UniformRand();
02847 r2 = G4UniformRand();
02848 a1 = std::sqrt(-2.0*std::log(r2));
02849 ran1 = a1*std::sin(r1)*0.020*numberofFinalStateNucleons*GeV;
02850 ran2 = a1*std::cos(r1)*0.020*numberofFinalStateNucleons*GeV;
02851 G4ThreeVector frm(ran1, ran2, 0);
02852
02853 pseudoParticle[0] = pseudoParticle[0]+frm;
02854 pseudoParticle[2] = temp;
02855 pseudoParticle[3] = pseudoParticle[0];
02856
02857 pseudoParticle[1] = pseudoParticle[2].cross(pseudoParticle[3]);
02858 G4double rotation = 2.*pi*G4UniformRand();
02859 pseudoParticle[1] = pseudoParticle[1].rotate(rotation, pseudoParticle[3]);
02860 pseudoParticle[2] = pseudoParticle[3].cross(pseudoParticle[1]);
02861 for(G4int ii=1; ii<=3; ii++)
02862 {
02863 p = pseudoParticle[ii].mag();
02864 if( p == 0.0 )
02865 pseudoParticle[ii]= G4ThreeVector( 0.0, 0.0, 0.0 );
02866 else
02867 pseudoParticle[ii]= pseudoParticle[ii] * (1./p);
02868 }
02869
02870 pxTemp = pseudoParticle[1].dot(currentParticle.GetMomentum());
02871 pyTemp = pseudoParticle[2].dot(currentParticle.GetMomentum());
02872 pzTemp = pseudoParticle[3].dot(currentParticle.GetMomentum());
02873 currentParticle.SetMomentum( pxTemp, pyTemp, pzTemp );
02874
02875 pxTemp = pseudoParticle[1].dot(targetParticle.GetMomentum());
02876 pyTemp = pseudoParticle[2].dot(targetParticle.GetMomentum());
02877 pzTemp = pseudoParticle[3].dot(targetParticle.GetMomentum());
02878 targetParticle.SetMomentum( pxTemp, pyTemp, pzTemp );
02879
02880 for( i=0; i<vecLen; ++i )
02881 {
02882 pxTemp = pseudoParticle[1].dot(vec[i]->GetMomentum());
02883 pyTemp = pseudoParticle[2].dot(vec[i]->GetMomentum());
02884 pzTemp = pseudoParticle[3].dot(vec[i]->GetMomentum());
02885 vec[i]->SetMomentum( pxTemp, pyTemp, pzTemp );
02886 }
02887
02888
02889
02890
02891 Defs1( modifiedOriginal, currentParticle, targetParticle, vec, vecLen );
02892 G4double ekin;
02893 G4double dekin = 0.0;
02894 G4double ek1 = 0.0;
02895 G4int npions = 0;
02896 if( atomicWeight >= 1.5 )
02897 {
02898
02899
02900 const G4double alem[] = { 1.40, 2.30, 2.70, 3.00, 3.40, 4.60, 7.00 };
02901 const G4double val0[] = { 0.00, 0.40, 0.48, 0.51, 0.54, 0.60, 0.65 };
02902 alekw = std::log( originalIncident->GetKineticEnergy()/GeV );
02903 exh = 1.0;
02904 if( alekw > alem[0] )
02905 {
02906 exh = val0[6];
02907 for( G4int j=1; j<7; ++j )
02908 {
02909 if( alekw < alem[j] )
02910 {
02911 G4double rcnve = (val0[j] - val0[j-1]) / (alem[j] - alem[j-1]);
02912 exh = rcnve * alekw + val0[j-1] - rcnve * alem[j-1];
02913 break;
02914 }
02915 }
02916 exh = 1.0 - exh;
02917 }
02918 const G4double cfa = 0.025*((atomicWeight-1.)/120.)*std::exp(-(atomicWeight-1.)/120.);
02919 ekin = currentParticle.GetKineticEnergy()/GeV - cfa*(1+normal()/2.0);
02920 ekin = std::max( 1.0e-6, ekin );
02921 xxh = 1.0;
02922 if( (modifiedOriginal.GetDefinition() == aPiPlus ||
02923 modifiedOriginal.GetDefinition() == aPiMinus) &&
02924 currentParticle.GetDefinition() == aPiZero &&
02925 G4UniformRand() <= logWeight) xxh = exh;
02926 dekin += ekin*(1.0-xxh);
02927 ekin *= xxh;
02928 if (currentParticle.GetDefinition()->GetParticleSubType() == "pi") {
02929 ++npions;
02930 ek1 += ekin;
02931 }
02932 currentParticle.SetKineticEnergy( ekin*GeV );
02933 pp = currentParticle.GetTotalMomentum()/MeV;
02934 pp1 = currentParticle.GetMomentum().mag()/MeV;
02935 if( pp1 < 0.001*MeV )
02936 {
02937 G4double costheta = 2.*G4UniformRand() - 1.;
02938 G4double sintheta = std::sqrt(1. - costheta*costheta);
02939 G4double phi = twopi*G4UniformRand();
02940 currentParticle.SetMomentum( pp*sintheta*std::cos(phi)*MeV,
02941 pp*sintheta*std::sin(phi)*MeV,
02942 pp*costheta*MeV ) ;
02943 }
02944 else
02945 currentParticle.SetMomentum( currentParticle.GetMomentum() * (pp/pp1) );
02946 ekin = targetParticle.GetKineticEnergy()/GeV - cfa*(1+normal()/2.0);
02947 ekin = std::max( 1.0e-6, ekin );
02948 xxh = 1.0;
02949 if( (modifiedOriginal.GetDefinition() == aPiPlus ||
02950 modifiedOriginal.GetDefinition() == aPiMinus) &&
02951 targetParticle.GetDefinition() == aPiZero &&
02952 G4UniformRand() < logWeight) xxh = exh;
02953 dekin += ekin*(1.0-xxh);
02954 ekin *= xxh;
02955 if (targetParticle.GetDefinition()->GetParticleSubType() == "pi") {
02956 ++npions;
02957 ek1 += ekin;
02958 }
02959 targetParticle.SetKineticEnergy( ekin*GeV );
02960 pp = targetParticle.GetTotalMomentum()/MeV;
02961 pp1 = targetParticle.GetMomentum().mag()/MeV;
02962 if( pp1 < 0.001*MeV )
02963 {
02964 G4double costheta = 2.*G4UniformRand() - 1.;
02965 G4double sintheta = std::sqrt(1. - costheta*costheta);
02966 G4double phi = twopi*G4UniformRand();
02967 targetParticle.SetMomentum( pp*sintheta*std::cos(phi)*MeV,
02968 pp*sintheta*std::sin(phi)*MeV,
02969 pp*costheta*MeV ) ;
02970 }
02971 else
02972 targetParticle.SetMomentum( targetParticle.GetMomentum() * (pp/pp1) );
02973 for( i=0; i<vecLen; ++i )
02974 {
02975 ekin = vec[i]->GetKineticEnergy()/GeV - cfa*(1+normal()/2.0);
02976 ekin = std::max( 1.0e-6, ekin );
02977 xxh = 1.0;
02978 if( (modifiedOriginal.GetDefinition() == aPiPlus ||
02979 modifiedOriginal.GetDefinition() == aPiMinus) &&
02980 vec[i]->GetDefinition() == aPiZero &&
02981 G4UniformRand() < logWeight) xxh = exh;
02982 dekin += ekin*(1.0-xxh);
02983 ekin *= xxh;
02984 if (vec[i]->GetDefinition()->GetParticleSubType() == "pi") {
02985 ++npions;
02986 ek1 += ekin;
02987 }
02988 vec[i]->SetKineticEnergy( ekin*GeV );
02989 pp = vec[i]->GetTotalMomentum()/MeV;
02990 pp1 = vec[i]->GetMomentum().mag()/MeV;
02991 if( pp1 < 0.001*MeV )
02992 {
02993 G4double costheta = 2.*G4UniformRand() - 1.;
02994 G4double sintheta = std::sqrt(1. - costheta*costheta);
02995 G4double phi = twopi*G4UniformRand();
02996 vec[i]->SetMomentum( pp*sintheta*std::cos(phi)*MeV,
02997 pp*sintheta*std::sin(phi)*MeV,
02998 pp*costheta*MeV ) ;
02999 }
03000 else
03001 vec[i]->SetMomentum( vec[i]->GetMomentum() * (pp/pp1) );
03002 }
03003 }
03004 if( (ek1 != 0.0) && (npions > 0) )
03005 {
03006 dekin = 1.0 + dekin/ek1;
03007
03008
03009
03010 if (currentParticle.GetDefinition()->GetParticleSubType() == "pi")
03011 {
03012 currentParticle.SetKineticEnergy(
03013 std::max( 0.001*MeV, dekin*currentParticle.GetKineticEnergy() ) );
03014 pp = currentParticle.GetTotalMomentum()/MeV;
03015 pp1 = currentParticle.GetMomentum().mag()/MeV;
03016 if( pp1 < 0.001 )
03017 {
03018 G4double costheta = 2.*G4UniformRand() - 1.;
03019 G4double sintheta = std::sqrt(1. - costheta*costheta);
03020 G4double phi = twopi*G4UniformRand();
03021 currentParticle.SetMomentum( pp*sintheta*std::cos(phi)*MeV,
03022 pp*sintheta*std::sin(phi)*MeV,
03023 pp*costheta*MeV ) ;
03024 } else {
03025 currentParticle.SetMomentum( currentParticle.GetMomentum() * (pp/pp1) );
03026 }
03027 }
03028
03029 if (targetParticle.GetDefinition()->GetParticleSubType() == "pi")
03030 {
03031 targetParticle.SetKineticEnergy(
03032 std::max( 0.001*MeV, dekin*targetParticle.GetKineticEnergy() ) );
03033 pp = targetParticle.GetTotalMomentum()/MeV;
03034 pp1 = targetParticle.GetMomentum().mag()/MeV;
03035 if( pp1 < 0.001 )
03036 {
03037 G4double costheta = 2.*G4UniformRand() - 1.;
03038 G4double sintheta = std::sqrt(1. - costheta*costheta);
03039 G4double phi = twopi*G4UniformRand();
03040 targetParticle.SetMomentum( pp*sintheta*std::cos(phi)*MeV,
03041 pp*sintheta*std::sin(phi)*MeV,
03042 pp*costheta*MeV ) ;
03043 } else {
03044 targetParticle.SetMomentum( targetParticle.GetMomentum() * (pp/pp1) );
03045 }
03046 }
03047
03048 for( i=0; i<vecLen; ++i )
03049 {
03050 if (vec[i]->GetDefinition()->GetParticleSubType() == "pi")
03051 {
03052 vec[i]->SetKineticEnergy( std::max( 0.001*MeV, dekin*vec[i]->GetKineticEnergy() ) );
03053 pp = vec[i]->GetTotalMomentum()/MeV;
03054 pp1 = vec[i]->GetMomentum().mag()/MeV;
03055 if( pp1 < 0.001 )
03056 {
03057 G4double costheta = 2.*G4UniformRand() - 1.;
03058 G4double sintheta = std::sqrt(1. - costheta*costheta);
03059 G4double phi = twopi*G4UniformRand();
03060 vec[i]->SetMomentum( pp*sintheta*std::cos(phi)*MeV,
03061 pp*sintheta*std::sin(phi)*MeV,
03062 pp*costheta*MeV ) ;
03063 } else {
03064 vec[i]->SetMomentum( vec[i]->GetMomentum() * (pp/pp1) );
03065 }
03066 }
03067 }
03068 }
03069 }
03070
03071 void G4ReactionDynamics::AddBlackTrackParticles(
03072 const G4double epnb,
03073 const G4int npnb,
03074 const G4double edta,
03075 const G4int ndta,
03076 const G4double sprob,
03077 const G4double kineticMinimum,
03078 const G4double kineticFactor,
03079 const G4ReactionProduct &modifiedOriginal,
03080 G4int PinNucleus,
03081 G4int NinNucleus,
03082 const G4Nucleus &targetNucleus,
03083 G4FastVector<G4ReactionProduct,GHADLISTSIZE> &vec,
03084 G4int &vecLen )
03085 {
03086
03087
03088
03089
03090
03091
03092
03093
03094 G4ParticleDefinition *aProton = G4Proton::Proton();
03095 G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
03096 G4ParticleDefinition *aDeuteron = G4Deuteron::Deuteron();
03097 G4ParticleDefinition *aTriton = G4Triton::Triton();
03098 G4ParticleDefinition *anAlpha = G4Alpha::Alpha();
03099
03100 const G4double ekOriginal = modifiedOriginal.GetKineticEnergy()/MeV;
03101 const G4double atomicWeight = G4double(targetNucleus.GetA_asInt());
03102 const G4double atomicNumber = G4double(targetNucleus.GetZ_asInt());
03103
03104 const G4double ika1 = 3.6;
03105 const G4double ika2 = 35.56;
03106 const G4double ika3 = 6.45;
03107
03108 G4int i;
03109 G4double pp;
03110 G4double kinCreated = 0;
03111 G4double cfa = 0.025*((atomicWeight-1.0)/120.0) * std::exp(-(atomicWeight-1.0)/120.0);
03112
03113
03114
03115 if (npnb > 0)
03116 {
03117 G4double backwardKinetic = 0.0;
03118 G4int local_npnb = npnb;
03119 for( i=0; i<npnb; ++i ) if( G4UniformRand() < sprob ) local_npnb--;
03120 G4double local_epnb = epnb;
03121 if (ndta == 0) local_epnb += edta;
03122 G4double ekin = local_epnb/std::max(1,local_npnb);
03123
03124 for( i=0; i<local_npnb; ++i )
03125 {
03126 G4ReactionProduct * p1 = new G4ReactionProduct();
03127 if( backwardKinetic > local_epnb )
03128 {
03129 delete p1;
03130 break;
03131 }
03132 G4double ran = G4UniformRand();
03133 G4double kinetic = -ekin*std::log(ran) - cfa*(1.0+0.5*normal());
03134 if( kinetic < 0.0 )kinetic = -0.010*std::log(ran);
03135 backwardKinetic += kinetic;
03136 if( backwardKinetic > local_epnb )
03137 kinetic = std::max( kineticMinimum, local_epnb-(backwardKinetic-kinetic) );
03138
03139 if (G4UniformRand() > (1.0-atomicNumber/atomicWeight)) {
03140
03141
03142
03143 if (PinNucleus > 0) {
03144 p1->SetDefinition( aProton );
03145 PinNucleus--;
03146 } else if (NinNucleus > 0) {
03147 p1->SetDefinition( aNeutron );
03148 NinNucleus--;
03149 } else {
03150 delete p1;
03151 break;
03152 }
03153 } else {
03154
03155
03156
03157 if (NinNucleus > 0) {
03158 p1->SetDefinition( aNeutron );
03159 NinNucleus--;
03160 } else if (PinNucleus > 0) {
03161 p1->SetDefinition( aProton );
03162 PinNucleus--;
03163 } else {
03164 delete p1;
03165 break;
03166 }
03167 }
03168
03169 vec.SetElement( vecLen, p1 );
03170 G4double cost = G4UniformRand() * 2.0 - 1.0;
03171 G4double sint = std::sqrt(std::fabs(1.0-cost*cost));
03172 G4double phi = twopi * G4UniformRand();
03173 vec[vecLen]->SetNewlyAdded( true );
03174 vec[vecLen]->SetKineticEnergy( kinetic*GeV );
03175 kinCreated+=kinetic;
03176 pp = vec[vecLen]->GetTotalMomentum()/MeV;
03177 vec[vecLen]->SetMomentum( pp*sint*std::sin(phi)*MeV,
03178 pp*sint*std::cos(phi)*MeV,
03179 pp*cost*MeV );
03180 vecLen++;
03181
03182 }
03183
03184 if (NinNucleus > 0) {
03185 if( (atomicWeight >= 10.0) && (ekOriginal <= 2.0*GeV) )
03186 {
03187 G4double ekw = ekOriginal/GeV;
03188 G4int ika, kk = 0;
03189 if( ekw > 1.0 )ekw *= ekw;
03190 ekw = std::max( 0.1, ekw );
03191 ika = G4int(ika1*std::exp((atomicNumber*atomicNumber/
03192 atomicWeight-ika2)/ika3)/ekw);
03193 if( ika > 0 )
03194 {
03195 for( i=(vecLen-1); i>=0; --i )
03196 {
03197 if( (vec[i]->GetDefinition() == aProton) && vec[i]->GetNewlyAdded() )
03198 {
03199 vec[i]->SetDefinitionAndUpdateE( aNeutron );
03200 PinNucleus++;
03201 NinNucleus--;
03202 if( ++kk > ika )break;
03203 }
03204 }
03205 }
03206 }
03207 }
03208 }
03209
03210
03211
03212 if (ndta > 0)
03213 {
03214 G4double backwardKinetic = 0.0;
03215 G4int local_ndta=ndta;
03216 for( i=0; i<ndta; ++i )if( G4UniformRand() < sprob )local_ndta--;
03217 G4double local_edta = edta;
03218 if (npnb == 0) local_edta += epnb;
03219 G4double ekin = local_edta/std::max(1,local_ndta);
03220
03221 for( i=0; i<local_ndta; ++i )
03222 {
03223 G4ReactionProduct *p2 = new G4ReactionProduct();
03224 if( backwardKinetic > local_edta )
03225 {
03226 delete p2;
03227 break;
03228 }
03229 G4double ran = G4UniformRand();
03230 G4double kinetic = -ekin*std::log(ran)-cfa*(1.+0.5*normal());
03231 if( kinetic < 0.0 )kinetic = kineticFactor*std::log(ran);
03232 backwardKinetic += kinetic;
03233 if( backwardKinetic > local_edta )kinetic = local_edta-(backwardKinetic-kinetic);
03234 if( kinetic < 0.0 )kinetic = kineticMinimum;
03235 G4double cost = 2.0*G4UniformRand() - 1.0;
03236 G4double sint = std::sqrt(std::max(0.0,(1.0-cost*cost)));
03237 G4double phi = twopi*G4UniformRand();
03238 ran = G4UniformRand();
03239 if (ran < 0.60) {
03240 if (PinNucleus > 0 && NinNucleus > 0) {
03241 p2->SetDefinition( aDeuteron );
03242 PinNucleus--;
03243 NinNucleus--;
03244 } else if (NinNucleus > 0) {
03245 p2->SetDefinition( aNeutron );
03246 NinNucleus--;
03247 } else if (PinNucleus > 0) {
03248 p2->SetDefinition( aProton );
03249 PinNucleus--;
03250 } else {
03251 delete p2;
03252 break;
03253 }
03254 } else if (ran < 0.90) {
03255 if (PinNucleus > 0 && NinNucleus > 1) {
03256 p2->SetDefinition( aTriton );
03257 PinNucleus--;
03258 NinNucleus -= 2;
03259 } else if (PinNucleus > 0 && NinNucleus > 0) {
03260 p2->SetDefinition( aDeuteron );
03261 PinNucleus--;
03262 NinNucleus--;
03263 } else if (NinNucleus > 0) {
03264 p2->SetDefinition( aNeutron );
03265 NinNucleus--;
03266 } else if (PinNucleus > 0) {
03267 p2->SetDefinition( aProton );
03268 PinNucleus--;
03269 } else {
03270 delete p2;
03271 break;
03272 }
03273 } else {
03274 if (PinNucleus > 1 && NinNucleus > 1) {
03275 p2->SetDefinition( anAlpha );
03276 PinNucleus -= 2;
03277 NinNucleus -= 2;
03278 } else if (PinNucleus > 0 && NinNucleus > 1) {
03279 p2->SetDefinition( aTriton );
03280 PinNucleus--;
03281 NinNucleus -= 2;
03282 } else if (PinNucleus > 0 && NinNucleus > 0) {
03283 p2->SetDefinition( aDeuteron );
03284 PinNucleus--;
03285 NinNucleus--;
03286 } else if (NinNucleus > 0) {
03287 p2->SetDefinition( aNeutron );
03288 NinNucleus--;
03289 } else if (PinNucleus > 0) {
03290 p2->SetDefinition( aProton );
03291 PinNucleus--;
03292 } else {
03293 delete p2;
03294 break;
03295 }
03296 }
03297
03298 vec.SetElement( vecLen, p2 );
03299 vec[vecLen]->SetNewlyAdded( true );
03300 vec[vecLen]->SetKineticEnergy( kinetic*GeV );
03301 kinCreated+=kinetic;
03302 pp = vec[vecLen]->GetTotalMomentum()/MeV;
03303 vec[vecLen++]->SetMomentum( pp*sint*std::sin(phi)*MeV,
03304 pp*sint*std::cos(phi)*MeV,
03305 pp*cost*MeV );
03306
03307 }
03308 }
03309
03310
03311 }
03312
03313
03314 std::pair<G4int, G4int> G4ReactionDynamics::GetFinalStateNucleons(
03315 const G4DynamicParticle* originalTarget,
03316 const G4FastVector<G4ReactionProduct,GHADLISTSIZE>& vec,
03317 const G4int& vecLen)
03318 {
03319
03320
03321 G4int protonsRemoved = 0;
03322 G4int neutronsRemoved = 0;
03323 if (originalTarget->GetDefinition()->GetParticleName() == "proton")
03324 protonsRemoved++;
03325 else
03326 neutronsRemoved++;
03327
03328 G4String secName;
03329 for (G4int i = 0; i < vecLen; i++) {
03330 secName = vec[i]->GetDefinition()->GetParticleName();
03331 if (secName == "proton") {
03332 protonsRemoved++;
03333 } else if (secName == "neutron") {
03334 neutronsRemoved++;
03335 } else if (secName == "anti_proton") {
03336 protonsRemoved--;
03337 } else if (secName == "anti_neutron") {
03338 neutronsRemoved--;
03339 }
03340 }
03341
03342 return std::pair<G4int, G4int>(protonsRemoved, neutronsRemoved);
03343 }
03344
03345
03346 void G4ReactionDynamics::MomentumCheck(
03347 const G4ReactionProduct &modifiedOriginal,
03348 G4ReactionProduct ¤tParticle,
03349 G4ReactionProduct &targetParticle,
03350 G4FastVector<G4ReactionProduct,GHADLISTSIZE> &vec,
03351 G4int &vecLen )
03352 {
03353 const G4double pOriginal = modifiedOriginal.GetTotalMomentum()/MeV;
03354 G4double testMomentum = currentParticle.GetMomentum().mag()/MeV;
03355 G4double pMass;
03356 if( testMomentum >= pOriginal )
03357 {
03358 pMass = currentParticle.GetMass()/MeV;
03359 currentParticle.SetTotalEnergy(
03360 std::sqrt( pMass*pMass + pOriginal*pOriginal )*MeV );
03361 currentParticle.SetMomentum(
03362 currentParticle.GetMomentum() * (pOriginal/testMomentum) );
03363 }
03364 testMomentum = targetParticle.GetMomentum().mag()/MeV;
03365 if( testMomentum >= pOriginal )
03366 {
03367 pMass = targetParticle.GetMass()/MeV;
03368 targetParticle.SetTotalEnergy(
03369 std::sqrt( pMass*pMass + pOriginal*pOriginal )*MeV );
03370 targetParticle.SetMomentum(
03371 targetParticle.GetMomentum() * (pOriginal/testMomentum) );
03372 }
03373 for( G4int i=0; i<vecLen; ++i )
03374 {
03375 testMomentum = vec[i]->GetMomentum().mag()/MeV;
03376 if( testMomentum >= pOriginal )
03377 {
03378 pMass = vec[i]->GetMass()/MeV;
03379 vec[i]->SetTotalEnergy(
03380 std::sqrt( pMass*pMass + pOriginal*pOriginal )*MeV );
03381 vec[i]->SetMomentum( vec[i]->GetMomentum() * (pOriginal/testMomentum) );
03382 }
03383 }
03384 }
03385
03386 void G4ReactionDynamics::ProduceStrangeParticlePairs(
03387 G4FastVector<G4ReactionProduct,GHADLISTSIZE> &vec,
03388 G4int &vecLen,
03389 const G4ReactionProduct &modifiedOriginal,
03390 const G4DynamicParticle *originalTarget,
03391 G4ReactionProduct ¤tParticle,
03392 G4ReactionProduct &targetParticle,
03393 G4bool &incidentHasChanged,
03394 G4bool &targetHasChanged )
03395 {
03396
03397
03398
03399
03400
03401
03402
03403
03404
03405 if( vecLen == 0 )return;
03406
03407
03408
03409 if( currentParticle.GetMass() == 0.0 || targetParticle.GetMass() == 0.0 )return;
03410
03411 const G4double etOriginal = modifiedOriginal.GetTotalEnergy()/GeV;
03412 const G4double mOriginal = modifiedOriginal.GetDefinition()->GetPDGMass()/GeV;
03413 G4double targetMass = originalTarget->GetDefinition()->GetPDGMass()/GeV;
03414 G4double centerofmassEnergy = std::sqrt( mOriginal*mOriginal +
03415 targetMass*targetMass +
03416 2.0*targetMass*etOriginal );
03417 G4double currentMass = currentParticle.GetMass()/GeV;
03418 G4double availableEnergy = centerofmassEnergy-(targetMass+currentMass);
03419 if( availableEnergy <= 1.0 )return;
03420
03421 G4ParticleDefinition *aProton = G4Proton::Proton();
03422 G4ParticleDefinition *anAntiProton = G4AntiProton::AntiProton();
03423 G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
03424 G4ParticleDefinition *anAntiNeutron = G4AntiNeutron::AntiNeutron();
03425 G4ParticleDefinition *aSigmaMinus = G4SigmaMinus::SigmaMinus();
03426 G4ParticleDefinition *aSigmaPlus = G4SigmaPlus::SigmaPlus();
03427 G4ParticleDefinition *aSigmaZero = G4SigmaZero::SigmaZero();
03428 G4ParticleDefinition *anAntiSigmaMinus = G4AntiSigmaMinus::AntiSigmaMinus();
03429 G4ParticleDefinition *anAntiSigmaPlus = G4AntiSigmaPlus::AntiSigmaPlus();
03430 G4ParticleDefinition *anAntiSigmaZero = G4AntiSigmaZero::AntiSigmaZero();
03431 G4ParticleDefinition *aKaonMinus = G4KaonMinus::KaonMinus();
03432 G4ParticleDefinition *aKaonPlus = G4KaonPlus::KaonPlus();
03433 G4ParticleDefinition *aKaonZL = G4KaonZeroLong::KaonZeroLong();
03434 G4ParticleDefinition *aKaonZS = G4KaonZeroShort::KaonZeroShort();
03435 G4ParticleDefinition *aLambda = G4Lambda::Lambda();
03436 G4ParticleDefinition *anAntiLambda = G4AntiLambda::AntiLambda();
03437
03438 const G4double protonMass = aProton->GetPDGMass()/GeV;
03439 const G4double sigmaMinusMass = aSigmaMinus->GetPDGMass()/GeV;
03440
03441
03442
03443 const G4double avrs[] = {3.,4.,5.,6.,7.,8.,9.,10.,20.,30.,40.,50.};
03444
03445 G4int ibin, i3, i4;
03446 G4double avk, avy, avn, ran;
03447 G4int i = 1;
03448 while( (i<12) && (centerofmassEnergy>avrs[i]) )++i;
03449 if( i == 12 )
03450 ibin = 11;
03451 else
03452 ibin = i;
03453
03454
03455
03456
03457 if( vecLen == 1 )
03458 {
03459 i3 = 0;
03460 i4 = 1;
03461 }
03462 else
03463 {
03464 ran = G4UniformRand();
03465 while( ran == 1.0 )ran = G4UniformRand();
03466 i4 = i3 = G4int( vecLen*ran );
03467 while( i3 == i4 )
03468 {
03469 ran = G4UniformRand();
03470 while( ran == 1.0 )ran = G4UniformRand();
03471 i4 = G4int( vecLen*ran );
03472 }
03473 }
03474
03475
03476
03477 const G4double avkkb[] = { 0.0015, 0.005, 0.012, 0.0285, 0.0525, 0.075,
03478 0.0975, 0.123, 0.28, 0.398, 0.495, 0.573 };
03479 const G4double avky[] = { 0.005, 0.03, 0.064, 0.095, 0.115, 0.13,
03480 0.145, 0.155, 0.20, 0.205, 0.210, 0.212 };
03481 const G4double avnnb[] = { 0.00001, 0.0001, 0.0006, 0.0025, 0.01, 0.02,
03482 0.04, 0.05, 0.12, 0.15, 0.18, 0.20 };
03483
03484 avk = (std::log(avkkb[ibin])-std::log(avkkb[ibin-1]))*(centerofmassEnergy-avrs[ibin-1])
03485 /(avrs[ibin]-avrs[ibin-1]) + std::log(avkkb[ibin-1]);
03486 avk = std::exp(avk);
03487
03488 avy = (std::log(avky[ibin])-std::log(avky[ibin-1]))*(centerofmassEnergy-avrs[ibin-1])
03489 /(avrs[ibin]-avrs[ibin-1]) + std::log(avky[ibin-1]);
03490 avy = std::exp(avy);
03491
03492 avn = (std::log(avnnb[ibin])-std::log(avnnb[ibin-1]))*(centerofmassEnergy-avrs[ibin-1])
03493 /(avrs[ibin]-avrs[ibin-1]) + std::log(avnnb[ibin-1]);
03494 avn = std::exp(avn);
03495
03496 if( avk+avy+avn <= 0.0 )return;
03497
03498 if( currentMass < protonMass )avy /= 2.0;
03499 if( targetMass < protonMass )avy = 0.0;
03500 avy += avk+avn;
03501 avk += avn;
03502 ran = G4UniformRand();
03503 if( ran < avn )
03504 {
03505 if( availableEnergy < 2.0 )return;
03506 if( vecLen == 1 )
03507 {
03508 G4ReactionProduct *p1 = new G4ReactionProduct;
03509 if( G4UniformRand() < 0.5 )
03510 {
03511 vec[0]->SetDefinition( aNeutron );
03512 p1->SetDefinition( anAntiNeutron );
03513 (G4UniformRand() < 0.5) ? p1->SetSide( -1 ) : p1->SetSide( 1 );
03514 vec[0]->SetMayBeKilled(false);
03515 p1->SetMayBeKilled(false);
03516 }
03517 else
03518 {
03519 vec[0]->SetDefinition( aProton );
03520 p1->SetDefinition( anAntiProton );
03521 (G4UniformRand() < 0.5) ? p1->SetSide( -1 ) : p1->SetSide( 1 );
03522 vec[0]->SetMayBeKilled(false);
03523 p1->SetMayBeKilled(false);
03524 }
03525 vec.SetElement( vecLen++, p1 );
03526
03527 }
03528 else
03529 {
03530 if( G4UniformRand() < 0.5 )
03531 {
03532 vec[i3]->SetDefinition( aNeutron );
03533 vec[i4]->SetDefinition( anAntiNeutron );
03534 vec[i3]->SetMayBeKilled(false);
03535 vec[i4]->SetMayBeKilled(false);
03536 }
03537 else
03538 {
03539 vec[i3]->SetDefinition( aProton );
03540 vec[i4]->SetDefinition( anAntiProton );
03541 vec[i3]->SetMayBeKilled(false);
03542 vec[i4]->SetMayBeKilled(false);
03543 }
03544 }
03545 }
03546 else if( ran < avk )
03547 {
03548 if( availableEnergy < 1.0 )return;
03549
03550 const G4double kkb[] = { 0.2500, 0.3750, 0.5000, 0.5625, 0.6250,
03551 0.6875, 0.7500, 0.8750, 1.000 };
03552 const G4int ipakkb1[] = { 10, 10, 10, 11, 11, 12, 12, 11, 12 };
03553 const G4int ipakkb2[] = { 13, 11, 12, 11, 12, 11, 12, 13, 13 };
03554 ran = G4UniformRand();
03555 i = 0;
03556 while( (i<9) && (ran>=kkb[i]) )++i;
03557 if( i == 9 )return;
03558
03559
03560
03561
03562 switch( ipakkb1[i] )
03563 {
03564 case 10:
03565 vec[i3]->SetDefinition( aKaonPlus );
03566 vec[i3]->SetMayBeKilled(false);
03567 break;
03568 case 11:
03569 vec[i3]->SetDefinition( aKaonZS );
03570 vec[i3]->SetMayBeKilled(false);
03571 break;
03572 case 12:
03573 vec[i3]->SetDefinition( aKaonZL );
03574 vec[i3]->SetMayBeKilled(false);
03575 break;
03576 }
03577 if( vecLen == 1 )
03578 {
03579 G4ReactionProduct *p1 = new G4ReactionProduct;
03580 switch( ipakkb2[i] )
03581 {
03582 case 11:
03583 p1->SetDefinition( aKaonZS );
03584 p1->SetMayBeKilled(false);
03585 break;
03586 case 12:
03587 p1->SetDefinition( aKaonZL );
03588 p1->SetMayBeKilled(false);
03589 break;
03590 case 13:
03591 p1->SetDefinition( aKaonMinus );
03592 p1->SetMayBeKilled(false);
03593 break;
03594 }
03595 (G4UniformRand() < 0.5) ? p1->SetSide( -1 ) : p1->SetSide( 1 );
03596 vec.SetElement( vecLen++, p1 );
03597
03598 }
03599 else
03600 {
03601 switch( ipakkb2[i] )
03602 {
03603 case 11:
03604 vec[i4]->SetDefinition( aKaonZS );
03605 vec[i4]->SetMayBeKilled(false);
03606 break;
03607 case 12:
03608 vec[i4]->SetDefinition( aKaonZL );
03609 vec[i4]->SetMayBeKilled(false);
03610 break;
03611 case 13:
03612 vec[i4]->SetDefinition( aKaonMinus );
03613 vec[i4]->SetMayBeKilled(false);
03614 break;
03615 }
03616 }
03617 }
03618 else if( ran < avy )
03619 {
03620 if( availableEnergy < 1.6 )return;
03621
03622 const G4double ky[] = { 0.200, 0.300, 0.400, 0.550, 0.625, 0.700,
03623 0.800, 0.850, 0.900, 0.950, 0.975, 1.000 };
03624 const G4int ipaky1[] = { 18, 18, 18, 20, 20, 20, 21, 21, 21, 22, 22, 22 };
03625 const G4int ipaky2[] = { 10, 11, 12, 10, 11, 12, 10, 11, 12, 10, 11, 12 };
03626 const G4int ipakyb1[] = { 19, 19, 19, 23, 23, 23, 24, 24, 24, 25, 25, 25 };
03627 const G4int ipakyb2[] = { 13, 12, 11, 13, 12, 11, 13, 12, 11, 13, 12, 11 };
03628 ran = G4UniformRand();
03629 i = 0;
03630 while( (i<12) && (ran>ky[i]) )++i;
03631 if( i == 12 )return;
03632 if( (currentMass<protonMass) || (G4UniformRand()<0.5) )
03633 {
03634
03635
03636
03637
03638
03639 switch( ipaky1[i] )
03640 {
03641 case 18:
03642 targetParticle.SetDefinition( aLambda );
03643 break;
03644 case 20:
03645 targetParticle.SetDefinition( aSigmaPlus );
03646 break;
03647 case 21:
03648 targetParticle.SetDefinition( aSigmaZero );
03649 break;
03650 case 22:
03651 targetParticle.SetDefinition( aSigmaMinus );
03652 break;
03653 }
03654 targetHasChanged = true;
03655 switch( ipaky2[i] )
03656 {
03657 case 10:
03658 vec[i3]->SetDefinition( aKaonPlus );
03659 vec[i3]->SetMayBeKilled(false);
03660 break;
03661 case 11:
03662 vec[i3]->SetDefinition( aKaonZS );
03663 vec[i3]->SetMayBeKilled(false);
03664 break;
03665 case 12:
03666 vec[i3]->SetDefinition( aKaonZL );
03667 vec[i3]->SetMayBeKilled(false);
03668 break;
03669 }
03670 }
03671 else
03672 {
03673
03674
03675 if( (currentParticle.GetDefinition() == anAntiProton) ||
03676 (currentParticle.GetDefinition() == anAntiNeutron) ||
03677 (currentParticle.GetDefinition() == anAntiLambda) ||
03678 (currentMass > sigmaMinusMass) )
03679 {
03680 switch( ipakyb1[i] )
03681 {
03682 case 19:
03683 currentParticle.SetDefinitionAndUpdateE( anAntiLambda );
03684 break;
03685 case 23:
03686 currentParticle.SetDefinitionAndUpdateE( anAntiSigmaPlus );
03687 break;
03688 case 24:
03689 currentParticle.SetDefinitionAndUpdateE( anAntiSigmaZero );
03690 break;
03691 case 25:
03692 currentParticle.SetDefinitionAndUpdateE( anAntiSigmaMinus );
03693 break;
03694 }
03695 incidentHasChanged = true;
03696 switch( ipakyb2[i] )
03697 {
03698 case 11:
03699 vec[i3]->SetDefinition( aKaonZS );
03700 vec[i3]->SetMayBeKilled(false);
03701 break;
03702 case 12:
03703 vec[i3]->SetDefinition( aKaonZL );
03704 vec[i3]->SetMayBeKilled(false);
03705 break;
03706 case 13:
03707 vec[i3]->SetDefinition( aKaonMinus );
03708 vec[i3]->SetMayBeKilled(false);
03709 break;
03710 }
03711 }
03712 else
03713 {
03714 switch( ipaky1[i] )
03715 {
03716 case 18:
03717 currentParticle.SetDefinitionAndUpdateE( aLambda );
03718 break;
03719 case 20:
03720 currentParticle.SetDefinitionAndUpdateE( aSigmaPlus );
03721 break;
03722 case 21:
03723 currentParticle.SetDefinitionAndUpdateE( aSigmaZero );
03724 break;
03725 case 22:
03726 currentParticle.SetDefinitionAndUpdateE( aSigmaMinus );
03727 break;
03728 }
03729 incidentHasChanged = true;
03730 switch( ipaky2[i] )
03731 {
03732 case 10:
03733 vec[i3]->SetDefinition( aKaonPlus );
03734 vec[i3]->SetMayBeKilled(false);
03735 break;
03736 case 11:
03737 vec[i3]->SetDefinition( aKaonZS );
03738 vec[i3]->SetMayBeKilled(false);
03739 break;
03740 case 12:
03741 vec[i3]->SetDefinition( aKaonZL );
03742 vec[i3]->SetMayBeKilled(false);
03743 break;
03744 }
03745 }
03746 }
03747 }
03748 else return;
03749
03750
03751
03752
03753
03754
03755
03756
03757
03758 currentMass = currentParticle.GetMass()/GeV;
03759 targetMass = targetParticle.GetMass()/GeV;
03760
03761 G4double energyCheck = centerofmassEnergy-(currentMass+targetMass);
03762 for( i=0; i<vecLen; ++i )
03763 {
03764 energyCheck -= vec[i]->GetMass()/GeV;
03765 if( energyCheck < 0.0 )
03766 {
03767 vecLen = std::max( 0, --i );
03768 G4int j;
03769 for(j=i; j<vecLen; j++) delete vec[j];
03770 break;
03771 }
03772 }
03773 return;
03774 }
03775
03776 void
03777 G4ReactionDynamics::NuclearReaction(
03778 G4FastVector<G4ReactionProduct,4> &vec,
03779 G4int &vecLen,
03780 const G4HadProjectile *originalIncident,
03781 const G4Nucleus &targetNucleus,
03782 const G4double theAtomicMass,
03783 const G4double *mass )
03784 {
03785
03786
03787
03788
03789 G4ParticleDefinition *aGamma = G4Gamma::Gamma();
03790 G4ParticleDefinition *aProton = G4Proton::Proton();
03791 G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
03792 G4ParticleDefinition *aDeuteron = G4Deuteron::Deuteron();
03793 G4ParticleDefinition *aTriton = G4Triton::Triton();
03794 G4ParticleDefinition *anAlpha = G4Alpha::Alpha();
03795
03796 const G4double aProtonMass = aProton->GetPDGMass()/MeV;
03797 const G4double aNeutronMass = aNeutron->GetPDGMass()/MeV;
03798 const G4double aDeuteronMass = aDeuteron->GetPDGMass()/MeV;
03799 const G4double aTritonMass = aTriton->GetPDGMass()/MeV;
03800 const G4double anAlphaMass = anAlpha->GetPDGMass()/MeV;
03801
03802 G4ReactionProduct currentParticle;
03803 currentParticle = *originalIncident;
03804
03805
03806
03807
03808
03809 G4double p = currentParticle.GetTotalMomentum();
03810 G4double pp = currentParticle.GetMomentum().mag();
03811 if( pp <= 0.001*MeV )
03812 {
03813 G4double phinve = twopi*G4UniformRand();
03814 G4double rthnve = std::acos( std::max( -1.0, std::min( 1.0, -1.0 + 2.0*G4UniformRand() ) ) );
03815 currentParticle.SetMomentum( p*std::sin(rthnve)*std::cos(phinve),
03816 p*std::sin(rthnve)*std::sin(phinve),
03817 p*std::cos(rthnve) );
03818 }
03819 else
03820 currentParticle.SetMomentum( currentParticle.GetMomentum() * (p/pp) );
03821
03822
03823
03824 G4double currentKinetic = currentParticle.GetKineticEnergy()/MeV;
03825 G4double currentMass = currentParticle.GetDefinition()->GetPDGMass()/MeV;
03826 G4double qv = currentKinetic + theAtomicMass + currentMass;
03827
03828 G4double qval[9];
03829 qval[0] = qv - mass[0];
03830 qval[1] = qv - mass[1] - aNeutronMass;
03831 qval[2] = qv - mass[2] - aProtonMass;
03832 qval[3] = qv - mass[3] - aDeuteronMass;
03833 qval[4] = qv - mass[4] - aTritonMass;
03834 qval[5] = qv - mass[5] - anAlphaMass;
03835 qval[6] = qv - mass[6] - aNeutronMass - aNeutronMass;
03836 qval[7] = qv - mass[7] - aNeutronMass - aProtonMass;
03837 qval[8] = qv - mass[8] - aProtonMass - aProtonMass;
03838
03839 if( currentParticle.GetDefinition() == aNeutron )
03840 {
03841 const G4double A = G4double(targetNucleus.GetA_asInt());
03842 if( G4UniformRand() > ((A-1.0)/230.0)*((A-1.0)/230.0) )
03843 qval[0] = 0.0;
03844 if( G4UniformRand() >= currentKinetic/7.9254*A )
03845 qval[2] = qval[3] = qval[4] = qval[5] = qval[8] = 0.0;
03846 }
03847 else
03848 qval[0] = 0.0;
03849
03850 G4int i;
03851 qv = 0.0;
03852 for( i=0; i<9; ++i )
03853 {
03854 if( mass[i] < 500.0*MeV )qval[i] = 0.0;
03855 if( qval[i] < 0.0 )qval[i] = 0.0;
03856 qv += qval[i];
03857 }
03858 G4double qv1 = 0.0;
03859 G4double ran = G4UniformRand();
03860 G4int index;
03861 for( index=0; index<9; ++index )
03862 {
03863 if( qval[index] > 0.0 )
03864 {
03865 qv1 += qval[index]/qv;
03866 if( ran <= qv1 )break;
03867 }
03868 }
03869 if( index == 9 )
03870 {
03871 throw G4HadronicException(__FILE__, __LINE__,
03872 "G4ReactionDynamics::NuclearReaction: inelastic reaction kinematically not possible");
03873 }
03874 G4double ke = currentParticle.GetKineticEnergy()/GeV;
03875 G4int nt = 2;
03876 if( (index>=6) || (G4UniformRand()<std::min(0.5,ke*10.0)) )nt = 3;
03877
03878 G4ReactionProduct **v = new G4ReactionProduct * [3];
03879 v[0] = new G4ReactionProduct;
03880 v[1] = new G4ReactionProduct;
03881 v[2] = new G4ReactionProduct;
03882
03883 v[0]->SetMass( mass[index]*MeV );
03884 switch( index )
03885 {
03886 case 0:
03887 v[1]->SetDefinition( aGamma );
03888 v[2]->SetDefinition( aGamma );
03889 break;
03890 case 1:
03891 v[1]->SetDefinition( aNeutron );
03892 v[2]->SetDefinition( aGamma );
03893 break;
03894 case 2:
03895 v[1]->SetDefinition( aProton );
03896 v[2]->SetDefinition( aGamma );
03897 break;
03898 case 3:
03899 v[1]->SetDefinition( aDeuteron );
03900 v[2]->SetDefinition( aGamma );
03901 break;
03902 case 4:
03903 v[1]->SetDefinition( aTriton );
03904 v[2]->SetDefinition( aGamma );
03905 break;
03906 case 5:
03907 v[1]->SetDefinition( anAlpha );
03908 v[2]->SetDefinition( aGamma );
03909 break;
03910 case 6:
03911 v[1]->SetDefinition( aNeutron );
03912 v[2]->SetDefinition( aNeutron );
03913 break;
03914 case 7:
03915 v[1]->SetDefinition( aNeutron );
03916 v[2]->SetDefinition( aProton );
03917 break;
03918 case 8:
03919 v[1]->SetDefinition( aProton );
03920 v[2]->SetDefinition( aProton );
03921 break;
03922 }
03923
03924
03925
03926 G4ReactionProduct pseudo1;
03927 pseudo1.SetMass( theAtomicMass*MeV );
03928 pseudo1.SetTotalEnergy( theAtomicMass*MeV );
03929 G4ReactionProduct pseudo2 = currentParticle + pseudo1;
03930 pseudo2.SetMomentum( pseudo2.GetMomentum() * (-1.0) );
03931
03932
03933
03934 G4FastVector<G4ReactionProduct,GHADLISTSIZE> tempV;
03935 tempV.Initialize( nt );
03936 G4int tempLen = 0;
03937 tempV.SetElement( tempLen++, v[0] );
03938 tempV.SetElement( tempLen++, v[1] );
03939 if( nt == 3 )tempV.SetElement( tempLen++, v[2] );
03940 G4bool constantCrossSection = true;
03941 GenerateNBodyEvent( pseudo2.GetMass()/MeV, constantCrossSection, tempV, tempLen );
03942 v[0]->Lorentz( *v[0], pseudo2 );
03943 v[1]->Lorentz( *v[1], pseudo2 );
03944 if( nt == 3 )v[2]->Lorentz( *v[2], pseudo2 );
03945
03946 G4bool particleIsDefined = false;
03947 if( v[0]->GetMass()/MeV - aProtonMass < 0.1 )
03948 {
03949 v[0]->SetDefinition( aProton );
03950 particleIsDefined = true;
03951 }
03952 else if( v[0]->GetMass()/MeV - aNeutronMass < 0.1 )
03953 {
03954 v[0]->SetDefinition( aNeutron );
03955 particleIsDefined = true;
03956 }
03957 else if( v[0]->GetMass()/MeV - aDeuteronMass < 0.1 )
03958 {
03959 v[0]->SetDefinition( aDeuteron );
03960 particleIsDefined = true;
03961 }
03962 else if( v[0]->GetMass()/MeV - aTritonMass < 0.1 )
03963 {
03964 v[0]->SetDefinition( aTriton );
03965 particleIsDefined = true;
03966 }
03967 else if( v[0]->GetMass()/MeV - anAlphaMass < 0.1 )
03968 {
03969 v[0]->SetDefinition( anAlpha );
03970 particleIsDefined = true;
03971 }
03972 currentParticle.SetKineticEnergy(
03973 std::max( 0.001, currentParticle.GetKineticEnergy()/MeV ) );
03974 p = currentParticle.GetTotalMomentum();
03975 pp = currentParticle.GetMomentum().mag();
03976 if( pp <= 0.001*MeV )
03977 {
03978 G4double phinve = twopi*G4UniformRand();
03979 G4double rthnve = std::acos( std::max( -1.0, std::min( 1.0, -1.0 + 2.0*G4UniformRand() ) ) );
03980 currentParticle.SetMomentum( p*std::sin(rthnve)*std::cos(phinve),
03981 p*std::sin(rthnve)*std::sin(phinve),
03982 p*std::cos(rthnve) );
03983 }
03984 else
03985 currentParticle.SetMomentum( currentParticle.GetMomentum() * (p/pp) );
03986
03987 if( particleIsDefined )
03988 {
03989 v[0]->SetKineticEnergy(
03990 std::max( 0.001, 0.5*G4UniformRand()*v[0]->GetKineticEnergy()/MeV ) );
03991 p = v[0]->GetTotalMomentum();
03992 pp = v[0]->GetMomentum().mag();
03993 if( pp <= 0.001*MeV )
03994 {
03995 G4double phinve = twopi*G4UniformRand();
03996 G4double rthnve = std::acos( std::max(-1.0,std::min(1.0,-1.0+2.0*G4UniformRand())) );
03997 v[0]->SetMomentum( p*std::sin(rthnve)*std::cos(phinve),
03998 p*std::sin(rthnve)*std::sin(phinve),
03999 p*std::cos(rthnve) );
04000 }
04001 else
04002 v[0]->SetMomentum( v[0]->GetMomentum() * (p/pp) );
04003 }
04004 if( (v[1]->GetDefinition() == aDeuteron) ||
04005 (v[1]->GetDefinition() == aTriton) ||
04006 (v[1]->GetDefinition() == anAlpha) )
04007 v[1]->SetKineticEnergy(
04008 std::max( 0.001, 0.5*G4UniformRand()*v[1]->GetKineticEnergy()/MeV ) );
04009 else
04010 v[1]->SetKineticEnergy( std::max( 0.001, v[1]->GetKineticEnergy()/MeV ) );
04011
04012 p = v[1]->GetTotalMomentum();
04013 pp = v[1]->GetMomentum().mag();
04014 if( pp <= 0.001*MeV )
04015 {
04016 G4double phinve = twopi*G4UniformRand();
04017 G4double rthnve = std::acos( std::max(-1.0,std::min(1.0,-1.0+2.0*G4UniformRand())) );
04018 v[1]->SetMomentum( p*std::sin(rthnve)*std::cos(phinve),
04019 p*std::sin(rthnve)*std::sin(phinve),
04020 p*std::cos(rthnve) );
04021 }
04022 else
04023 v[1]->SetMomentum( v[1]->GetMomentum() * (p/pp) );
04024
04025 if( nt == 3 )
04026 {
04027 if( (v[2]->GetDefinition() == aDeuteron) ||
04028 (v[2]->GetDefinition() == aTriton) ||
04029 (v[2]->GetDefinition() == anAlpha) )
04030 v[2]->SetKineticEnergy(
04031 std::max( 0.001, 0.5*G4UniformRand()*v[2]->GetKineticEnergy()/MeV ) );
04032 else
04033 v[2]->SetKineticEnergy( std::max( 0.001, v[2]->GetKineticEnergy()/MeV ) );
04034
04035 p = v[2]->GetTotalMomentum();
04036 pp = v[2]->GetMomentum().mag();
04037 if( pp <= 0.001*MeV )
04038 {
04039 G4double phinve = twopi*G4UniformRand();
04040 G4double rthnve = std::acos( std::max(-1.0,std::min(1.0,-1.0+2.0*G4UniformRand())) );
04041 v[2]->SetMomentum( p*std::sin(rthnve)*std::cos(phinve),
04042 p*std::sin(rthnve)*std::sin(phinve),
04043 p*std::cos(rthnve) );
04044 }
04045 else
04046 v[2]->SetMomentum( v[2]->GetMomentum() * (p/pp) );
04047 }
04048 G4int del;
04049 for(del=0; del<vecLen; del++) delete vec[del];
04050 vecLen = 0;
04051 if( particleIsDefined )
04052 {
04053 vec.SetElement( vecLen++, v[0] );
04054 }
04055 else
04056 {
04057 delete v[0];
04058 }
04059 vec.SetElement( vecLen++, v[1] );
04060 if( nt == 3 )
04061 {
04062 vec.SetElement( vecLen++, v[2] );
04063 }
04064 else
04065 {
04066 delete v[2];
04067 }
04068 delete [] v;
04069 return;
04070 }
04071
04072
04073