114 G4bool &incidentHasChanged,
132 if(vecLen == 0)
return false;
141 G4bool veryForward =
false;
148 G4double centerofmassEnergy = std::sqrt( mOriginal*mOriginal +
149 targetMass*targetMass +
150 2.0*targetMass*etOriginal );
152 targetMass = targetParticle.
GetMass()/GeV;
157 for( i=0; i<vecLen; ++i )
161 *vec[itemp] = *vec[i];
166 if( currentMass == 0.0 && targetMass == 0.0 )
173 currentParticle = *vec[0];
174 targetParticle = *vec[1];
175 for( i=0; i<(vecLen-2); ++i )*vec[i] = *vec[i+2];
178 temp = vec[vecLen-2];
181 currentMass = currentParticle.
GetMass()/GeV;
182 targetMass = targetParticle.
GetMass()/GeV;
183 incidentHasChanged =
true;
184 targetHasChanged =
true;
196 currentParticle = targetParticle;
197 targetParticle = temp;
198 incidentHasChanged =
true;
199 targetHasChanged =
true;
200 currentMass = currentParticle.
GetMass()/GeV;
201 targetMass = targetParticle.
GetMass()/GeV;
203 const G4double afc = std::min( 0.75,
204 0.312+0.200*std::log(std::log(centerofmassEnergy*centerofmassEnergy))+
205 std::pow(centerofmassEnergy*centerofmassEnergy,1.5)/6000.0 );
207 G4double freeEnergy = centerofmassEnergy-currentMass-targetMass;
208 G4double forwardEnergy = freeEnergy/2.;
209 G4int forwardCount = 1;
211 G4double backwardEnergy = freeEnergy/2.;
212 G4int backwardCount = 1;
216 if(currentParticle.
GetSide()==-1)
218 forwardEnergy += currentMass;
220 backwardEnergy -= currentMass;
223 if(targetParticle.
GetSide()!=-1)
225 backwardEnergy += targetMass;
227 forwardEnergy -= targetMass;
232 for( i=0; i<vecLen; ++i )
234 if( vec[i]->GetSide() == -1 )
237 backwardEnergy -= vec[i]->GetMass()/GeV;
240 forwardEnergy -= vec[i]->GetMass()/GeV;
250 xtarg = afc * (std::pow(atomicWeight,0.33)-1.0) * (2.0*backwardCount+vecLen+2)/2.0;
252 xtarg = afc * (std::pow(atomicWeight,0.33)-1.0) * (2.0*backwardCount);
253 if( xtarg <= 0.0 )xtarg = 0.01;
254 G4int nuclearExcitationCount = Poisson( xtarg );
255 if(atomicWeight<1.0001) nuclearExcitationCount = 0;
256 G4int extraNucleonCount = 0;
258 if( nuclearExcitationCount > 0 )
260 const G4double nucsup[] = { 1.00, 0.7, 0.5, 0.4, 0.35, 0.3 };
261 const G4double psup[] = { 3., 6., 20., 50., 100., 1000. };
262 G4int momentumBin = 0;
263 while( (momentumBin < 6) &&
266 momentumBin = std::min( 5, momentumBin );
272 for( i=0; i<nuclearExcitationCount; ++i )
283 backwardEnergy += pVec->
GetMass()/GeV;
284 extraNucleonMass += pVec->
GetMass()/GeV;
291 else if( ran < 0.6819 )
300 backwardEnergy -= pVec->
GetMass()/GeV;
308 while( forwardEnergy <= 0.0 )
313 G4int forwardParticlesLeft = 0;
314 for( i=(vecLen-1); i>=0; --i )
316 if( vec[i]->GetSide() == 1 && vec[i]->GetMayBeKilled())
318 forwardParticlesLeft = 1;
321 forwardEnergy += vec[i]->GetMass()/GeV;
322 for(
G4int j=i; j<(vecLen-1); j++ )*vec[j] = *vec[j+1];
324 delete vec[vecLen-1];
325 if( --vecLen == 0 )
return false;
331 if( forwardParticlesLeft == 0 )
334 for (i = 0; i < vecLen; i++) {
335 if (vec[i]->GetDefinition()->GetParticleSubType() ==
"pi") {
341 for (i = 0; i < vecLen; i++) {
342 if (vec[i]->GetDefinition()->GetParticleSubType() ==
"kaon") {
348 if (iremove == -1) iremove = 0;
350 forwardEnergy += vec[iremove]->GetMass()/GeV;
351 if (vec[iremove]->GetSide() > 0) --forwardCount;
353 for (i = iremove; i < vecLen-1; i++) *vec[i] = *vec[i+1];
354 delete vec[vecLen-1];
356 if (vecLen == 0)
return false;
362 while( backwardEnergy <= 0.0 )
367 G4int backwardParticlesLeft = 0;
368 for( i=(vecLen-1); i>=0; --i )
370 if( vec[i]->GetSide() < 0 && vec[i]->GetMayBeKilled())
372 backwardParticlesLeft = 1;
375 if( vec[i]->GetSide() == -2 )
378 extraNucleonMass -= vec[i]->GetMass()/GeV;
379 backwardEnergy -= vec[i]->GetTotalEnergy()/GeV;
381 backwardEnergy += vec[i]->GetTotalEnergy()/GeV;
382 for(
G4int j=i; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];
384 delete vec[vecLen-1];
385 if( --vecLen == 0 )
return false;
391 if( backwardParticlesLeft == 0 )
394 for (i = 0; i < vecLen; i++) {
395 if (vec[i]->GetDefinition()->GetParticleSubType() ==
"pi") {
401 for (i = 0; i < vecLen; i++) {
402 if (vec[i]->GetDefinition()->GetParticleSubType() ==
"kaon") {
408 if (iremove == -1) iremove = 0;
410 backwardEnergy += vec[iremove]->GetMass()/GeV;
411 if (vec[iremove]->GetSide() > 0) --backwardCount;
413 for (i = iremove; i < vecLen-1; i++) *vec[i] = *vec[i+1];
414 delete vec[vecLen-1];
416 if (vecLen == 0)
return false;
427 for( i=0; i<10; ++i )pseudoParticle[i].SetZero();
429 pseudoParticle[0].
SetMass( mOriginal*GeV );
430 pseudoParticle[0].
SetMomentum( 0.0, 0.0, pOriginal*GeV );
432 std::sqrt( pOriginal*pOriginal + mOriginal*mOriginal )*GeV );
434 pseudoParticle[1].
SetMass( protonMass*MeV );
437 pseudoParticle[3].
SetMass( protonMass*(1+extraNucleonCount)*MeV );
438 pseudoParticle[3].
SetTotalEnergy( protonMass*(1+extraNucleonCount)*MeV );
440 pseudoParticle[8].
SetMomentum( 1.0*GeV, 0.0, 0.0 );
442 pseudoParticle[2] = pseudoParticle[0] + pseudoParticle[1];
443 pseudoParticle[3] = pseudoParticle[3] + pseudoParticle[0];
445 pseudoParticle[0].
Lorentz( pseudoParticle[0], pseudoParticle[2] );
446 pseudoParticle[1].
Lorentz( pseudoParticle[1], pseudoParticle[2] );
453 G4double aspar, pt, et, x, pp, pp1, rmb, wgt;
454 G4int innerCounter, outerCounter;
455 G4bool eliminateThisParticle, resetEnergies, constantCrossSection;
457 G4double forwardKinetic = 0.0, backwardKinetic = 0.0;
463 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,
464 1.43,1.67,2.0,2.5,3.33,5.00,10.00};
465 G4int backwardNucleonCount = 0;
466 G4double totalEnergy, kineticEnergy, vecMass;
468 for( i=(vecLen-1); i>=0; --i )
471 if( vec[i]->GetNewlyAdded() )
473 if( vec[i]->GetSide() == -2 )
475 if( backwardNucleonCount < 18 )
477 if (vec[i]->GetDefinition()->GetParticleSubType() ==
"pi") {
478 for(
G4int j=0; j<vecLen; j++)
delete vec[j];
481 "G4ReactionDynamics::GenerateXandPt : a pion has been counted as a backward nucleon");
483 vec[i]->SetSide( -3 );
484 ++backwardNucleonCount;
493 vecMass = vec[i]->GetMass()/GeV;
495 if( vec[i]->GetSide() == -2 )
497 if (vec[i]->GetDefinition()->GetParticleSubType() ==
"kaon" ||
498 vec[i]->GetDefinition()->GetParticleSubType() ==
"pi") {
500 pt = std::sqrt( std::pow( ran, 1.7 ) );
503 pt = std::sqrt( std::pow( ran, 1.2 ) );
508 if (vec[i]->GetDefinition()->GetParticleSubType() ==
"pi") {
510 pt = std::sqrt( std::pow( ran, 1.7 ) );
511 }
else if (vec[i]->GetDefinition()->GetParticleSubType() ==
"kaon") {
513 pt = std::sqrt( std::pow( ran, 1.7 ) );
516 pt = std::sqrt( std::pow( ran, 1.5 ) );
519 pt = std::max( 0.001, pt );
520 vec[i]->SetMomentum( pt*std::cos(phi)*GeV, pt*std::sin(phi)*GeV );
521 for(
G4int j=0; j<20; ++j )binl[j] = j/(19.*pt);
522 if( vec[i]->GetSide() > 0 )
531 eliminateThisParticle =
true;
532 resetEnergies =
true;
533 while( ++outerCounter < 3 )
535 for( l=1; l<20; ++l )
537 x = (binl[l]+binl[l-1])/2.;
538 pt = std::max( 0.001, pt );
540 dndl[l] += dndl[l-1];
542 dndl[l] = et * aspar/std::sqrt( std::pow((1.+aspar*x*aspar*x),3) )
543 * (binl[l]-binl[l-1]) / std::sqrt( pt*x*et*pt*x*et + pt*pt + vecMass*vecMass )
547 vec[i]->SetMomentum( pt*std::cos(phi)*GeV, pt*std::sin(phi)*GeV );
551 while( ++innerCounter < 7 )
555 while( ( ran >= dndl[l] ) && ( l < 20 ) )l++;
556 l = std::min( 19, l );
557 x = std::min( 1.0, pt*(binl[l-1] +
G4UniformRand()*(binl[l]-binl[l-1]) ) );
558 if( vec[i]->GetSide() < 0 )x *= -1.;
559 vec[i]->SetMomentum( x*et*GeV );
560 totalEnergy = std::sqrt( x*et*x*et + pt*pt + vecMass*vecMass );
561 vec[i]->SetTotalEnergy( totalEnergy*GeV );
562 kineticEnergy = vec[i]->GetKineticEnergy()/GeV;
563 if( vec[i]->GetSide() > 0 )
565 if( (forwardKinetic+kineticEnergy) < 0.95*forwardEnergy )
567 pseudoParticle[4] = pseudoParticle[4] + (*vec[i]);
568 forwardKinetic += kineticEnergy;
569 pseudoParticle[6] = pseudoParticle[4] + pseudoParticle[5];
571 phi = pseudoParticle[6].Angle( pseudoParticle[8] );
572 if( pseudoParticle[6].GetMomentum().y()/MeV < 0.0 )phi = twopi - phi;
573 phi += pi + normal()*pi/12.0;
574 if( phi > twopi )phi -= twopi;
575 if( phi < 0.0 )phi = twopi - phi;
577 eliminateThisParticle =
false;
578 resetEnergies =
false;
581 if( innerCounter > 5 )
break;
582 if( backwardEnergy >= vecMass )
584 vec[i]->SetSide( -1 );
585 forwardEnergy += vecMass;
586 backwardEnergy -= vecMass;
590 if( extraNucleonCount > 19 )
592 G4double xxx = 0.95+0.05*extraNucleonCount/20.0;
593 if( (backwardKinetic+kineticEnergy) < xxx*backwardEnergy )
595 pseudoParticle[5] = pseudoParticle[5] + (*vec[i]);
596 backwardKinetic += kineticEnergy;
597 pseudoParticle[6] = pseudoParticle[4] + pseudoParticle[5];
599 phi = pseudoParticle[6].Angle( pseudoParticle[8] );
600 if( pseudoParticle[6].GetMomentum().y()/MeV < 0.0 )phi = twopi - phi;
601 phi += pi + normal() * pi / 12.0;
602 if( phi > twopi )phi -= twopi;
603 if( phi < 0.0 )phi = twopi - phi;
605 eliminateThisParticle =
false;
606 resetEnergies =
false;
609 if( innerCounter > 5 )
break;
610 if( forwardEnergy >= vecMass )
612 vec[i]->SetSide( 1 );
613 forwardEnergy -= vecMass;
614 backwardEnergy += vecMass;
619 vec[i]->SetMomentum( momentum.
x() * 0.9, momentum.
y() * 0.9 );
630 forwardKinetic = 0.0;
631 backwardKinetic = 0.0;
634 for( l=i+1; l<vecLen; ++l )
636 if (vec[l]->GetSide() > 0 ||
637 vec[l]->GetDefinition()->GetParticleSubType() ==
"kaon" ||
638 vec[l]->GetDefinition()->GetParticleSubType() ==
"pi") {
640 G4double tempMass = vec[l]->GetMass()/MeV;
641 totalEnergy = 0.95*vec[l]->GetTotalEnergy()/MeV + 0.05*tempMass;
642 totalEnergy = std::max( tempMass, totalEnergy );
643 vec[l]->SetTotalEnergy( totalEnergy*MeV );
644 pp = std::sqrt( std::abs( totalEnergy*totalEnergy - tempMass*tempMass ) );
645 pp1 = vec[l]->GetMomentum().mag()/MeV;
646 if( pp1 < 1.0e-6*GeV )
649 G4double sintheta = std::sqrt(1. - costheta*costheta);
651 vec[l]->SetMomentum( pp*sintheta*std::cos(phi2)*MeV,
652 pp*sintheta*std::sin(phi2)*MeV,
655 vec[l]->SetMomentum( vec[l]->GetMomentum() * (pp/pp1) );
657 G4double px = vec[l]->GetMomentum().x()/MeV;
658 G4double py = vec[l]->GetMomentum().y()/MeV;
659 pt = std::max( 1.0, std::sqrt( px*px + py*py ) )/GeV;
660 if( vec[l]->GetSide() > 0 )
662 forwardKinetic += vec[l]->GetKineticEnergy()/GeV;
663 pseudoParticle[4] = pseudoParticle[4] + (*vec[l]);
665 backwardKinetic += vec[l]->GetKineticEnergy()/GeV;
666 pseudoParticle[5] = pseudoParticle[5] + (*vec[l]);
673 if( eliminateThisParticle && vec[i]->GetMayBeKilled())
675 if( vec[i]->GetSide() > 0 )
678 forwardEnergy += vecMass;
680 if( vec[i]->GetSide() == -2 )
683 extraNucleonMass -= vecMass;
684 backwardEnergy -= vecMass;
687 backwardEnergy += vecMass;
689 for(
G4int j=i; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];
693 if( --vecLen == 0 )
return false;
694 pseudoParticle[6] = pseudoParticle[4] + pseudoParticle[5];
696 phi = pseudoParticle[6].Angle( pseudoParticle[8] );
697 if( pseudoParticle[6].GetMomentum().y()/MeV < 0.0 )phi = twopi - phi;
698 phi += pi + normal() * pi / 12.0;
699 if( phi > twopi )phi -= twopi;
700 if( phi < 0.0 )phi = twopi - phi;
713 pt = std::sqrt( std::pow( ran/6.0, 1.7 ) );
716 pt = std::sqrt( std::pow( ran/5.0, 1.4 ) );
719 pt = std::sqrt( std::pow( ran/4.0, 1.2 ) );
722 for(
G4int j=0; j<20; ++j )binl[j] = j/(19.*pt);
723 currentParticle.
SetMomentum( pt*std::cos(phi)*GeV, pt*std::sin(phi)*GeV );
726 vecMass = currentParticle.
GetMass()/GeV;
727 for( l=1; l<20; ++l )
729 x = (binl[l]+binl[l-1])/2.;
731 dndl[l] += dndl[l-1];
733 dndl[l] = aspar/std::sqrt( std::pow((1.+
sqr(aspar*x)),3) ) *
734 (binl[l]-binl[l-1]) * et / std::sqrt( pt*x*et*pt*x*et + pt*pt + vecMass*vecMass ) +
739 while( (ran>dndl[l]) && (l<20) )l++;
740 l = std::min( 19, l );
741 x = std::min( 1.0, pt*(binl[l-1] +
G4UniformRand()*(binl[l]-binl[l-1]) ) );
743 if( forwardEnergy < forwardKinetic )
744 totalEnergy = vecMass + 0.04*std::fabs(normal());
746 totalEnergy = vecMass + forwardEnergy - forwardKinetic;
748 pp = std::sqrt( std::abs( totalEnergy*totalEnergy - vecMass*vecMass ) )*GeV;
750 if( pp1 < 1.0e-6*GeV )
753 G4double sintheta = std::sqrt(1. - costheta*costheta);
755 currentParticle.
SetMomentum( pp*sintheta*std::cos(phi2)*MeV,
756 pp*sintheta*std::sin(phi2)*MeV,
761 pseudoParticle[4] = pseudoParticle[4] + currentParticle;
769 if( backwardNucleonCount < 18 )
772 ++backwardNucleonCount;
779 vecMass = targetParticle.
GetMass()/GeV;
782 pt = std::max( 0.001, std::sqrt( std::pow( ran/4.0, 1.2 ) ) );
783 targetParticle.
SetMomentum( pt*std::cos(phi)*GeV, pt*std::sin(phi)*GeV );
784 for(
G4int j=0; j<20; ++j )binl[j] = (j-1.)/(19.*pt);
788 eliminateThisParticle =
true;
789 resetEnergies =
true;
790 while( ++outerCounter < 3 )
792 for( l=1; l<20; ++l )
794 x = (binl[l]+binl[l-1])/2.;
796 dndl[l] += dndl[l-1];
798 dndl[l] = aspar/std::sqrt( std::pow((1.+aspar*x*aspar*x),3) ) *
799 (binl[l]-binl[l-1])*et / std::sqrt( pt*x*et*pt*x*et+pt*pt+vecMass*vecMass ) +
803 while( ++innerCounter < 7 )
807 while( ( ran >= dndl[l] ) && ( l < 20 ) )l++;
808 l = std::min( 19, l );
809 x = std::min( 1.0, pt*(binl[l-1] +
G4UniformRand()*(binl[l]-binl[l-1]) ) );
810 if( targetParticle.
GetSide() < 0 )x *= -1.;
812 totalEnergy = std::sqrt( x*et*x*et + pt*pt + vecMass*vecMass );
814 if( targetParticle.
GetSide() < 0 )
816 if( extraNucleonCount > 19 )x=0.999;
817 G4double xxx = 0.95+0.05*extraNucleonCount/20.0;
818 if( (backwardKinetic+totalEnergy-vecMass) < xxx*backwardEnergy )
820 pseudoParticle[5] = pseudoParticle[5] + targetParticle;
821 backwardKinetic += totalEnergy - vecMass;
822 pseudoParticle[6] = pseudoParticle[4] + pseudoParticle[5];
824 phi = pseudoParticle[6].Angle( pseudoParticle[8] );
825 if( pseudoParticle[6].GetMomentum().y()/MeV < 0.0 )phi = twopi - phi;
826 phi += pi + normal() * pi / 12.0;
827 if( phi > twopi )phi -= twopi;
828 if( phi < 0.0 )phi = twopi - phi;
830 eliminateThisParticle =
false;
831 resetEnergies =
false;
834 if( innerCounter > 5 )
break;
835 if( forwardEnergy >= vecMass )
838 forwardEnergy -= vecMass;
839 backwardEnergy += vecMass;
843 targetParticle.
SetMomentum( momentum.
x() * 0.9, momentum.
y() * 0.9 );
849 if( forwardEnergy < forwardKinetic )
850 totalEnergy = vecMass + 0.04*std::fabs(normal());
852 totalEnergy = vecMass + forwardEnergy - forwardKinetic;
854 pp = std::sqrt( std::abs( totalEnergy*totalEnergy - vecMass*vecMass ) )*GeV;
856 if( pp1 < 1.0e-6*GeV )
859 G4double sintheta = std::sqrt(1. - costheta*costheta);
861 targetParticle.
SetMomentum( pp*sintheta*std::cos(phi2)*MeV,
862 pp*sintheta*std::sin(phi2)*MeV,
868 pseudoParticle[4] = pseudoParticle[4] + targetParticle;
870 eliminateThisParticle =
false;
871 resetEnergies =
false;
882 forwardKinetic = backwardKinetic = 0.0;
885 for( l=0; l<vecLen; ++l )
887 if (vec[l]->GetSide() > 0 ||
888 vec[l]->GetDefinition()->GetParticleSubType() ==
"kaon" ||
889 vec[l]->GetDefinition()->GetParticleSubType() ==
"pi") {
890 G4double tempMass = vec[l]->GetMass()/GeV;
892 std::max( tempMass, 0.95*vec[l]->GetTotalEnergy()/GeV + 0.05*tempMass );
893 vec[l]->SetTotalEnergy( totalEnergy*GeV );
894 pp = std::sqrt( std::abs( totalEnergy*totalEnergy - tempMass*tempMass ) )*GeV;
895 pp1 = vec[l]->GetMomentum().mag()/MeV;
896 if( pp1 < 1.0e-6*GeV )
899 G4double sintheta = std::sqrt(1. - costheta*costheta);
901 vec[l]->SetMomentum( pp*sintheta*std::cos(phi2)*MeV,
902 pp*sintheta*std::sin(phi2)*MeV,
906 vec[l]->SetMomentum( vec[l]->GetMomentum() * (pp/pp1) );
908 pt = std::max( 0.001*GeV, std::sqrt(
sqr(vec[l]->GetMomentum().x()/MeV) +
909 sqr(vec[l]->GetMomentum().y()/MeV) ) )/GeV;
910 if( vec[l]->GetSide() > 0)
912 forwardKinetic += vec[l]->GetKineticEnergy()/GeV;
913 pseudoParticle[4] = pseudoParticle[4] + (*vec[l]);
915 backwardKinetic += vec[l]->GetKineticEnergy()/GeV;
916 pseudoParticle[5] = pseudoParticle[5] + (*vec[l]);
929 pseudoParticle[6].
Lorentz( pseudoParticle[3], pseudoParticle[2] );
930 pseudoParticle[6] = pseudoParticle[6] - pseudoParticle[4];
931 pseudoParticle[6] = pseudoParticle[6] - pseudoParticle[5];
932 if( backwardNucleonCount == 1 )
935 std::min( backwardEnergy-backwardKinetic, centerofmassEnergy/2.0-protonMass/GeV );
937 if( ekin < 0.04 )ekin = 0.04 * std::fabs( normal() );
938 vecMass = targetParticle.
GetMass()/GeV;
939 totalEnergy = ekin+vecMass;
941 pp = std::sqrt( std::abs( totalEnergy*totalEnergy - vecMass*vecMass ) )*GeV;
943 if( pp1 < 1.0e-6*GeV )
946 G4double sintheta = std::sqrt(1. - costheta*costheta);
948 targetParticle.
SetMomentum( pp*sintheta*std::cos(phi2)*MeV,
949 pp*sintheta*std::sin(phi2)*MeV,
952 targetParticle.
SetMomentum( pseudoParticle[6].GetMomentum() * (pp/pp1) );
954 pseudoParticle[5] = pseudoParticle[5] + targetParticle;
958 const G4double cpar[] = { 0.6, 0.6, 0.35, 0.15, 0.10 };
959 const G4double gpar[] = { 2.6, 2.6, 1.80, 1.30, 1.20 };
963 if (backwardNucleonCount < 5)
965 tempCount = backwardNucleonCount;
975 if( targetParticle.
GetSide() == -3 )
976 rmb0 += targetParticle.
GetMass()/GeV;
977 for( i=0; i<vecLen; ++i )
979 if( vec[i]->GetSide() == -3 )rmb0 += vec[i]->GetMass()/GeV;
981 rmb = rmb0 + std::pow(-std::log(1.0-
G4UniformRand()),cpar[tempCount]) / gpar[tempCount];
983 vecMass = std::min( rmb, totalEnergy );
984 pseudoParticle[6].
SetMass( vecMass*GeV );
985 pp = std::sqrt( std::abs( totalEnergy*totalEnergy - vecMass*vecMass ) )*GeV;
987 if( pp1 < 1.0e-6*GeV )
990 G4double sintheta = std::sqrt(1. - costheta*costheta);
992 pseudoParticle[6].
SetMomentum( -pp*sintheta*std::cos(phi2)*MeV,
993 -pp*sintheta*std::sin(phi2)*MeV,
997 pseudoParticle[6].
SetMomentum( pseudoParticle[6].GetMomentum() * (-pp/pp1) );
1002 if( targetParticle.
GetSide() == -3 )tempV.
SetElement( tempLen++, &targetParticle );
1003 for( i=0; i<vecLen; ++i )
1005 if( vec[i]->GetSide() == -3 )tempV.
SetElement( tempLen++, vec[i] );
1007 if( tempLen != backwardNucleonCount )
1009 G4cerr <<
"tempLen is not the same as backwardNucleonCount" <<
G4endl;
1010 G4cerr <<
"tempLen = " << tempLen;
1011 G4cerr <<
", backwardNucleonCount = " << backwardNucleonCount <<
G4endl;
1014 for( i=0; i<vecLen; ++i )
1015 G4cerr <<
"particle #" << i <<
" side = " << vec[i]->GetSide() <<
G4endl;
1016 G4Exception(
"G4ReactionDynamics::GenerateXandPt",
"601",
1019 constantCrossSection =
true;
1024 pseudoParticle[6].GetMass(), constantCrossSection, tempV, tempLen );
1026 if( targetParticle.
GetSide() == -3 )
1028 targetParticle.
Lorentz( targetParticle, pseudoParticle[6] );
1030 pseudoParticle[5] = pseudoParticle[5] + targetParticle;
1032 for( i=0; i<vecLen; ++i )
1034 if( vec[i]->GetSide() == -3 )
1036 vec[i]->Lorentz( *vec[i], pseudoParticle[6] );
1037 pseudoParticle[5] = pseudoParticle[5] + (*vec[i]);
1046 if( vecLen == 0 )
return false;
1049 currentParticle.
Lorentz( currentParticle, pseudoParticle[1] );
1050 targetParticle.
Lorentz( targetParticle, pseudoParticle[1] );
1051 for( i=0; i<vecLen; ++i ) vec[i]->Lorentz( *vec[i], pseudoParticle[1] );
1063 G4bool leadingStrangeParticleHasChanged =
true;
1067 leadingStrangeParticleHasChanged =
false;
1068 if( leadingStrangeParticleHasChanged &&
1070 leadingStrangeParticleHasChanged =
false;
1071 if( leadingStrangeParticleHasChanged )
1073 for( i=0; i<vecLen; i++ )
1075 if( vec[i]->GetDefinition() == leadingStrangeParticle.
GetDefinition() )
1077 leadingStrangeParticleHasChanged =
false;
1082 if( leadingStrangeParticleHasChanged )
1093 if( (leadTest&&targetTest) || !(leadTest||targetTest) )
1096 targetHasChanged =
true;
1102 incidentHasChanged =
false;
1111 std::pair<G4int, G4int> finalStateNucleons =
1112 GetFinalStateNucleons(originalTarget, vec, vecLen);
1114 G4int protonsInFinalState = finalStateNucleons.first;
1115 G4int neutronsInFinalState = finalStateNucleons.second;
1117 G4int numberofFinalStateNucleons =
1118 protonsInFinalState + neutronsInFinalState;
1124 numberofFinalStateNucleons++;
1126 numberofFinalStateNucleons = std::max(1, numberofFinalStateNucleons);
1128 G4int PinNucleus = std::max(0,
1129 targetNucleus.
GetZ_asInt() - protonsInFinalState);
1130 G4int NinNucleus = std::max(0,
1131 targetNucleus.
GetN_asInt() - neutronsInFinalState);
1133 pseudoParticle[3].
SetMomentum( 0.0, 0.0, pOriginal*GeV );
1134 pseudoParticle[3].
SetMass( mOriginal*GeV );
1136 std::sqrt( pOriginal*pOriginal + mOriginal*mOriginal )*GeV );
1141 if(numberofFinalStateNucleons == 1) diff = 0;
1143 pseudoParticle[4].
SetMass( protonMass*(numberofFinalStateNucleons-diff)*MeV );
1144 pseudoParticle[4].
SetTotalEnergy( protonMass*(numberofFinalStateNucleons-diff)*MeV );
1149 currentParticle.
GetMass()/MeV -
1155 pseudoParticle[5] = pseudoParticle[3] + pseudoParticle[4];
1156 pseudoParticle[3].
Lorentz( pseudoParticle[3], pseudoParticle[5] );
1157 pseudoParticle[4].
Lorentz( pseudoParticle[4], pseudoParticle[5] );
1160 pseudoParticle[7] = pseudoParticle[7] + currentParticle;
1161 pseudoParticle[7] = pseudoParticle[7] + targetParticle;
1163 for( i=0; i<vecLen; ++i )
1165 pseudoParticle[7] = pseudoParticle[7] + *vec[i];
1166 simulatedKinetic += vec[i]->GetKineticEnergy()/MeV;
1167 theoreticalKinetic -= vec[i]->GetMass()/MeV;
1170 if( vecLen <= 16 && vecLen > 0 )
1176 tempR[0] = currentParticle;
1177 tempR[1] = targetParticle;
1178 for( i=0; i<vecLen; ++i )tempR[i+2] = *vec[i];
1182 for( i=0; i<vecLen+2; ++i )tempV.
SetElement( tempLen++, &tempR[i] );
1183 constantCrossSection =
true;
1186 pseudoParticle[4].GetTotalEnergy()/MeV,
1187 constantCrossSection, tempV, tempLen );
1190 for (i = 0; i < tempLen; i++) Qvalue += tempV[i]->GetMass();
1192 constantCrossSection, tempV, tempLen );
1196 theoreticalKinetic = 0.0;
1197 for( i=0; i<tempLen; ++i )
1199 pseudoParticle[6].
Lorentz( *tempV[i], pseudoParticle[4] );
1208 if( simulatedKinetic != 0.0 )
1210 wgt = (theoreticalKinetic)/simulatedKinetic;
1212 simulatedKinetic = theoreticalKinetic;
1216 if( pp1 < 1.0e-6*GeV )
1219 G4double sintheta = std::sqrt(1. - costheta*costheta);
1221 currentParticle.
SetMomentum( pp*sintheta*std::cos(phi2)*MeV,
1222 pp*sintheta*std::sin(phi2)*MeV,
1231 simulatedKinetic += theoreticalKinetic;
1235 if( pp1 < 1.0e-6*GeV )
1238 G4double sintheta = std::sqrt(1. - costheta*costheta);
1240 targetParticle.
SetMomentum( pp*sintheta*std::cos(phi2)*MeV,
1241 pp*sintheta*std::sin(phi2)*MeV,
1246 for( i=0; i<vecLen; ++i )
1248 theoreticalKinetic = vec[i]->GetKineticEnergy()/MeV * wgt;
1249 simulatedKinetic += theoreticalKinetic;
1250 vec[i]->SetKineticEnergy( theoreticalKinetic*MeV );
1251 pp = vec[i]->GetTotalMomentum()/MeV;
1252 pp1 = vec[i]->GetMomentum().mag()/MeV;
1253 if( pp1 < 1.0e-6*GeV )
1256 G4double sintheta = std::sqrt(1. - costheta*costheta);
1258 vec[i]->SetMomentum( pp*sintheta*std::cos(phi2)*MeV,
1259 pp*sintheta*std::sin(phi2)*MeV,
1263 vec[i]->SetMomentum( vec[i]->GetMomentum() * (pp/pp1) );
1268 Rotate( numberofFinalStateNucleons, pseudoParticle[3].GetMomentum(),
1269 modifiedOriginal, originalIncident, targetNucleus,
1270 currentParticle, targetParticle, vec, vecLen );
1277 if( atomicWeight >= 1.5 )
1298 const G4double kineticMinimum = 1.e-6;
1299 const G4double kineticFactor = -0.010;
1302 if( ekIncident >= 5.0 )sprob = std::min( 1.0, 0.6*std::log(ekIncident-4.0) );
1303 if( epnb >= pnCutOff )
1305 npnb = Poisson((1.5+1.25*numberofFinalStateNucleons)*epnb/(epnb+edta));
1306 if( numberofFinalStateNucleons + npnb > atomicWeight )
1307 npnb =
G4int(atomicWeight+0.00001 - numberofFinalStateNucleons);
1308 npnb = std::min( npnb, 127-vecLen );
1310 if( edta >= dtaCutOff )
1312 ndta = Poisson( (1.5+1.25*numberofFinalStateNucleons)*edta/(epnb+edta) );
1313 ndta = std::min( ndta, 127-vecLen );
1315 if (npnb == 0 && ndta == 0) npnb = 1;
1319 AddBlackTrackParticles(epnb, npnb, edta, ndta, sprob, kineticMinimum,
1320 kineticFactor, modifiedOriginal,
1321 PinNucleus, NinNucleus, targetNucleus,
1331 if( (atomicWeight >= 1.5) && (atomicWeight <= 230.0) && (ekOriginal <= 0.2) )
1334 currentParticle.
SetTOF( 1.0 );
1346 G4bool &incidentHasChanged,
1347 G4bool &targetHasChanged )
1356 G4double cmEnergy = std::sqrt( mOriginal*mOriginal + targetMass*targetMass +
1357 2.0*targetMass*etOriginal );
1358 G4double eAvailable = cmEnergy - mOriginal - targetMass;
1359 for (
G4int i = 0; i < vecLen; i++) eAvailable -= vec[i]->GetMass()/GeV;
1390 if (eAvailable > nucleonMass - piMass) {
1395 incidentHasChanged =
true;
1405 if (eAvailable > nucleonMass - piMass) {
1410 targetHasChanged =
true;
1414 for(
G4int i=0; i<vecLen; ++i )
1417 vec[i]->GetDefinition() == aPiPlus ||
1418 vec[i]->GetDefinition() == aPiZero ||
1419 vec[i]->GetDefinition() == aPiMinus ) &&
1423 if (eAvailable > nucleonMass - piMass) {
1425 vec[i]->SetDefinitionAndUpdateE( aNeutron );
1427 vec[i]->SetDefinitionAndUpdateE( aProton );
1443 G4bool &incidentHasChanged,
1444 G4bool &targetHasChanged,
1467 G4bool veryForward =
false;
1475 G4double centerofmassEnergy = std::sqrt( mOriginal*mOriginal +
1476 targetMass*targetMass +
1477 2.0*targetMass*etOriginal );
1479 targetMass = targetParticle.
GetMass()/GeV;
1481 if( currentMass == 0.0 && targetMass == 0.0 )
1485 currentParticle = *vec[0];
1486 targetParticle = *vec[1];
1487 for( i=0; i<(vecLen-2); ++i )*vec[i] = *vec[i+2];
1490 for(
G4int j=0; j<vecLen; j++)
delete vec[j];
1493 "G4ReactionDynamics::TwoCluster: Negative number of particles");
1495 delete vec[vecLen-1];
1496 delete vec[vecLen-2];
1498 currentMass = currentParticle.
GetMass()/GeV;
1499 targetMass = targetParticle.
GetMass()/GeV;
1500 incidentHasChanged =
true;
1501 targetHasChanged =
true;
1514 G4int forwardCount = 1;
1520 G4int backwardCount = 1;
1521 G4int backwardNucleonCount = 1;
1526 for( i=0; i<vecLen; ++i )
1528 if( vec[i]->GetSide() < 0 )vec[i]->SetSide( -1 );
1531 if( vec[i]->GetSide() == -1 )
1534 backwardMass += vec[i]->GetMass()/GeV;
1539 forwardMass += vec[i]->GetMass()/GeV;
1545 G4double term1 = std::log(centerofmassEnergy*centerofmassEnergy);
1546 if(term1 < 0) term1 = 0.0001;
1547 const G4double afc = 0.312 + 0.2 * std::log(term1);
1550 xtarg = afc * (std::pow(atomicWeight,0.33)-1.0) * (2*backwardCount+vecLen+2)/2.0;
1552 xtarg = afc * (std::pow(atomicWeight,0.33)-1.0) * (2*backwardCount);
1553 if( xtarg <= 0.0 )xtarg = 0.01;
1554 G4int nuclearExcitationCount = Poisson( xtarg );
1555 if(atomicWeight<1.0001) nuclearExcitationCount = 0;
1556 G4int extraNucleonCount = 0;
1559 if( nuclearExcitationCount > 0 )
1561 G4int momentumBin = std::min( 4,
G4int(pOriginal/3.0) );
1562 const G4double nucsup[] = { 1.0, 0.8, 0.6, 0.5, 0.4 };
1568 for( i=0; i<nuclearExcitationCount; ++i )
1577 ++backwardNucleonCount;
1578 ++extraNucleonCount;
1579 extraNucleonMass += pVec->
GetMass()/GeV;
1586 else if( ran < 0.6819 )
1592 extraMass += pVec->
GetMass()/GeV;
1599 G4double forwardEnergy = (centerofmassEnergy-cMass-bMass)/2.0 +cMass - forwardMass;
1600 G4double backwardEnergy = (centerofmassEnergy-cMass-bMass)/2.0 +bMass - backwardMass;
1601 G4double eAvailable = centerofmassEnergy - (forwardMass+backwardMass);
1605 while( eAvailable <= 0.0 )
1607 secondaryDeleted =
false;
1608 for( i=(vecLen-1); i>=0; --i )
1610 if( vec[i]->GetSide() == 1 && vec[i]->GetMayBeKilled())
1612 pMass = vec[i]->GetMass()/GeV;
1613 for(
G4int j=i; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];
1615 forwardEnergy += pMass;
1616 forwardMass -= pMass;
1617 secondaryDeleted =
true;
1620 else if( vec[i]->GetSide() == -1 && vec[i]->GetMayBeKilled())
1622 pMass = vec[i]->GetMass()/GeV;
1623 for(
G4int j=i; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];
1625 backwardEnergy += pMass;
1626 backwardMass -= pMass;
1627 secondaryDeleted =
true;
1631 if( secondaryDeleted )
1633 delete vec[vecLen-1];
1643 if( targetParticle.
GetSide() == -1 )
1645 pMass = targetParticle.
GetMass()/GeV;
1646 targetParticle = *vec[0];
1647 for(
G4int j=0; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];
1649 backwardEnergy += pMass;
1650 backwardMass -= pMass;
1651 secondaryDeleted =
true;
1653 else if( targetParticle.
GetSide() == 1 )
1655 pMass = targetParticle.
GetMass()/GeV;
1656 targetParticle = *vec[0];
1657 for(
G4int j=0; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];
1659 forwardEnergy += pMass;
1660 forwardMass -= pMass;
1661 secondaryDeleted =
true;
1663 if( secondaryDeleted )
1665 delete vec[vecLen-1];
1671 if( currentParticle.
GetSide() == -1 )
1673 pMass = currentParticle.
GetMass()/GeV;
1674 currentParticle = *vec[0];
1675 for(
G4int j=0; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];
1677 backwardEnergy += pMass;
1678 backwardMass -= pMass;
1679 secondaryDeleted =
true;
1681 else if( currentParticle.
GetSide() == 1 )
1683 pMass = currentParticle.
GetMass()/GeV;
1684 currentParticle = *vec[0];
1685 for(
G4int j=0; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];
1687 forwardEnergy += pMass;
1688 forwardMass -= pMass;
1689 secondaryDeleted =
true;
1691 if( secondaryDeleted )
1693 delete vec[vecLen-1];
1700 eAvailable = centerofmassEnergy - (forwardMass+backwardMass);
1710 const G4double cpar[] = { 0.6, 0.6, 0.35, 0.15, 0.10 };
1711 const G4double gpar[] = { 2.6, 2.6, 1.8, 1.30, 1.20 };
1713 if (forwardCount <= 0 || backwardCount <= 0)
return false;
1715 if (forwardCount == 1) rmc = forwardMass;
1718 G4int ntc = std::max(1, std::min(5,forwardCount))-1;
1719 rmc = forwardMass + std::pow(-std::log(1.0-
G4UniformRand()),cpar[ntc-1])/gpar[ntc-1];
1722 if (backwardCount == 1) rmd = backwardMass;
1725 G4int ntc = std::max(1, std::min(5,backwardCount));
1726 rmd = backwardMass + std::pow(-std::log(1.0-
G4UniformRand()),cpar[ntc-1])/gpar[ntc-1];
1729 while( rmc+rmd > centerofmassEnergy )
1731 if( (rmc <= forwardMass) && (rmd <= backwardMass) )
1733 G4double temp = 0.999*centerofmassEnergy/(rmc+rmd);
1739 rmc = 0.1*forwardMass + 0.9*rmc;
1740 rmd = 0.1*backwardMass + 0.9*rmd;
1745 for( i=0; i<8; ++i )pseudoParticle[i].SetZero();
1747 pseudoParticle[1].
SetMass( mOriginal*GeV );
1749 pseudoParticle[1].
SetMomentum( 0.0, 0.0, pOriginal*GeV );
1751 pseudoParticle[2].
SetMass( protonMass*MeV );
1757 pseudoParticle[0] = pseudoParticle[1] + pseudoParticle[2];
1758 pseudoParticle[1].
Lorentz( pseudoParticle[1], pseudoParticle[0] );
1759 pseudoParticle[2].
Lorentz( pseudoParticle[2], pseudoParticle[0] );
1762 G4double pf = (centerofmassEnergy*centerofmassEnergy+rmd*rmd-rmc*rmc);
1764 pf -= 4*centerofmassEnergy*centerofmassEnergy*rmd*rmd;
1765 pf = std::sqrt( std::max(pf,pfMin) )/(2.0*centerofmassEnergy);
1769 pseudoParticle[3].
SetMass( rmc*GeV );
1770 pseudoParticle[3].
SetTotalEnergy( std::sqrt(pf*pf+rmc*rmc)*GeV );
1772 pseudoParticle[4].
SetMass( rmd*GeV );
1773 pseudoParticle[4].
SetTotalEnergy( std::sqrt(pf*pf+rmd*rmd)*GeV );
1782 G4double dtb = 4.0*pin*pf*std::max( bMin, b1+b2*std::log(pOriginal) );
1783 G4double factor = 1.0 - std::exp(-dtb);
1785 costheta = std::max(-1.0, std::min(1.0, costheta));
1786 G4double sintheta = std::sqrt(1.0-costheta*costheta);
1792 pseudoParticle[3].
SetMomentum( pf*sintheta*std::cos(phi)*GeV,
1793 pf*sintheta*std::sin(phi)*GeV,
1796 pseudoParticle[4].
SetMomentum( pseudoParticle[3].GetMomentum() * (-1.0) );
1801 if( nuclearExcitationCount > 0 )
1806 if( ekOriginal <= 5.0 )
1808 ekit1 *= ekOriginal*ekOriginal/25.0;
1809 ekit2 *= ekOriginal*ekOriginal/25.0;
1811 const G4double a = (1.0-ga)/(std::pow(ekit2,(1.0-ga)) - std::pow(ekit1,(1.0-ga)));
1812 for( i=0; i<vecLen; ++i )
1814 if( vec[i]->GetSide() == -2 )
1817 std::pow( (
G4UniformRand()*(1.0-ga)/a+std::pow(ekit1,(1.0-ga))), (1.0/(1.0-ga)) );
1818 vec[i]->SetKineticEnergy( kineticE*GeV );
1819 G4double vMass = vec[i]->GetMass()/MeV;
1820 G4double totalE = kineticE*GeV + vMass;
1821 pp = std::sqrt( std::abs(totalE*totalE-vMass*vMass) );
1823 G4double sint = std::sqrt( std::max( 0.0, (1.0-cost*cost) ) );
1825 vec[i]->SetMomentum( pp*sint*std::sin(phi)*MeV,
1826 pp*sint*std::cos(phi)*MeV,
1828 vec[i]->Lorentz( *vec[i], pseudoParticle[0] );
1836 currentParticle.
SetMomentum( pseudoParticle[3].GetMomentum() );
1837 currentParticle.
SetTotalEnergy( pseudoParticle[3].GetTotalEnergy() );
1839 targetParticle.
SetMomentum( pseudoParticle[4].GetMomentum() );
1840 targetParticle.
SetTotalEnergy( pseudoParticle[4].GetTotalEnergy() );
1842 pseudoParticle[5].
SetMomentum( pseudoParticle[3].GetMomentum() * (-1.0) );
1843 pseudoParticle[5].
SetMass( pseudoParticle[3].GetMass() );
1844 pseudoParticle[5].
SetTotalEnergy( pseudoParticle[3].GetTotalEnergy() );
1846 pseudoParticle[6].
SetMomentum( pseudoParticle[4].GetMomentum() * (-1.0) );
1847 pseudoParticle[6].
SetMass( pseudoParticle[4].GetMass() );
1848 pseudoParticle[6].
SetTotalEnergy( pseudoParticle[4].GetTotalEnergy() );
1852 if( forwardCount > 1 )
1856 G4bool constantCrossSection =
true;
1858 if( currentParticle.
GetSide() == 1 )
1859 tempV.
SetElement( tempLen++, ¤tParticle );
1860 if( targetParticle.
GetSide() == 1 )
1861 tempV.
SetElement( tempLen++, &targetParticle );
1862 for( i=0; i<vecLen; ++i )
1864 if( vec[i]->GetSide() == 1 )
1870 vec[i]->SetSide( -1 );
1878 constantCrossSection, tempV, tempLen );
1879 if( currentParticle.
GetSide() == 1 )
1880 currentParticle.
Lorentz( currentParticle, pseudoParticle[5] );
1881 if( targetParticle.
GetSide() == 1 )
1882 targetParticle.
Lorentz( targetParticle, pseudoParticle[5] );
1883 for( i=0; i<vecLen; ++i )
1885 if( vec[i]->GetSide() == 1 )vec[i]->Lorentz( *vec[i], pseudoParticle[5] );
1890 if( backwardCount > 1 )
1894 G4bool constantCrossSection =
true;
1896 if( currentParticle.
GetSide() == -1 )
1897 tempV.
SetElement( tempLen++, ¤tParticle );
1898 if( targetParticle.
GetSide() == -1 )
1899 tempV.
SetElement( tempLen++, &targetParticle );
1900 for( i=0; i<vecLen; ++i )
1902 if( vec[i]->GetSide() == -1 )
1908 vec[i]->SetSide( -2 );
1909 vec[i]->SetKineticEnergy( 0.0 );
1910 vec[i]->SetMomentum( 0.0, 0.0, 0.0 );
1918 constantCrossSection, tempV, tempLen );
1919 if( currentParticle.
GetSide() == -1 )
1920 currentParticle.
Lorentz( currentParticle, pseudoParticle[6] );
1921 if( targetParticle.
GetSide() == -1 )
1922 targetParticle.
Lorentz( targetParticle, pseudoParticle[6] );
1923 for( i=0; i<vecLen; ++i )
1925 if( vec[i]->GetSide() == -1 )vec[i]->Lorentz( *vec[i], pseudoParticle[6] );
1933 currentParticle.
Lorentz( currentParticle, pseudoParticle[2] );
1934 targetParticle.
Lorentz( targetParticle, pseudoParticle[2] );
1935 for( i=0; i<vecLen; ++i ) vec[i]->Lorentz( *vec[i], pseudoParticle[2] );
1958 for( i=0; i<vecLen; ++i )
1960 if( vec[i]->GetDefinition() == leadingStrangeParticle.
GetDefinition() )
1971 if( ((leadMass < protonMass) && (targetParticle.
GetMass()/MeV < protonMass)) ||
1972 ((leadMass >= protonMass) && (targetParticle.
GetMass()/MeV >= protonMass)) )
1982 G4double sintheta2 = std::sqrt(1. - costheta2*costheta2);
1984 targetParticle.
SetMomentum( pp*sintheta2*std::cos(phi2)*MeV,
1985 pp*sintheta2*std::sin(phi2)*MeV,
1986 pp*costheta2*MeV ) ;
1991 targetHasChanged =
true;
2003 G4double sintheta2 = std::sqrt(1. - costheta2*costheta2);
2005 currentParticle.
SetMomentum( pp*sintheta2*std::cos(phi2)*MeV,
2006 pp*sintheta2*std::sin(phi2)*MeV,
2007 pp*costheta2*MeV ) ;
2012 incidentHasChanged =
true;
2020 std::pair<G4int, G4int> finalStateNucleons =
2021 GetFinalStateNucleons(originalTarget, vec, vecLen);
2023 G4int protonsInFinalState = finalStateNucleons.first;
2024 G4int neutronsInFinalState = finalStateNucleons.second;
2026 G4int numberofFinalStateNucleons =
2027 protonsInFinalState + neutronsInFinalState;
2033 numberofFinalStateNucleons++;
2035 numberofFinalStateNucleons = std::max(1, numberofFinalStateNucleons);
2037 G4int PinNucleus = std::max(0,
2038 targetNucleus.
GetZ_asInt() - protonsInFinalState);
2039 G4int NinNucleus = std::max(0,
2040 targetNucleus.
GetN_asInt() - neutronsInFinalState);
2045 pseudoParticle[4].
SetMass( mOriginal*GeV );
2047 pseudoParticle[4].
SetMomentum( 0.0, 0.0, pOriginal*GeV );
2052 if(numberofFinalStateNucleons == 1) diff = 0;
2054 pseudoParticle[5].
SetMass( protonMass*(numberofFinalStateNucleons-diff)*MeV);
2055 pseudoParticle[5].
SetTotalEnergy( protonMass*(numberofFinalStateNucleons-diff)*MeV);
2060 pseudoParticle[6] = pseudoParticle[4] + pseudoParticle[5];
2061 pseudoParticle[4].
Lorentz( pseudoParticle[4], pseudoParticle[6] );
2062 pseudoParticle[5].
Lorentz( pseudoParticle[5], pseudoParticle[6] );
2067 tempR[0] = currentParticle;
2068 tempR[1] = targetParticle;
2069 for( i=0; i<vecLen; ++i )tempR[i+2] = *vec[i];
2073 G4bool constantCrossSection =
true;
2075 for( i=0; i<vecLen+2; ++i )tempV.
SetElement( tempLen++, &tempR[i] );
2081 pseudoParticle[5].GetTotalEnergy()/MeV,
2082 constantCrossSection, tempV, tempLen );
2085 for (i = 0; i < tempLen; i++) Qvalue += tempV[i]->GetMass();
2087 constantCrossSection, tempV, tempLen );
2089 theoreticalKinetic = 0.0;
2090 for( i=0; i<vecLen+2; ++i )
2092 pseudoParticle[7].
SetMomentum( tempV[i]->GetMomentum() );
2093 pseudoParticle[7].
SetMass( tempV[i]->GetMass() );
2095 pseudoParticle[7].
Lorentz( pseudoParticle[7], pseudoParticle[5] );
2103 theoreticalKinetic -=
2104 ( currentParticle.
GetMass()/GeV + targetParticle.
GetMass()/GeV );
2105 for( i=0; i<vecLen; ++i )theoreticalKinetic -= vec[i]->GetMass()/GeV;
2109 for( i=0; i<vecLen; ++i )simulatedKinetic += vec[i]->GetKineticEnergy()/GeV;
2115 if( simulatedKinetic != 0.0 )
2117 wgt = (theoreticalKinetic)/simulatedKinetic;
2121 if( pp1 < 0.001*MeV )
2124 G4double sintheta2 = std::sqrt(1. - costheta2*costheta2);
2126 currentParticle.
SetMomentum( pp*sintheta2*std::cos(phi2)*MeV,
2127 pp*sintheta2*std::sin(phi2)*MeV,
2128 pp*costheta2*MeV ) ;
2136 if( pp1 < 0.001*MeV )
2139 G4double sintheta2 = std::sqrt(1. - costheta2*costheta2);
2141 targetParticle.
SetMomentum( pp*sintheta2*std::cos(phi2)*MeV,
2142 pp*sintheta2*std::sin(phi2)*MeV,
2143 pp*costheta2*MeV ) ;
2148 for( i=0; i<vecLen; ++i )
2150 vec[i]->SetKineticEnergy( wgt*vec[i]->GetKineticEnergy() );
2151 pp = vec[i]->GetTotalMomentum()/MeV;
2152 pp1 = vec[i]->GetMomentum().mag()/MeV;
2156 G4double sintheta2 = std::sqrt(1. - costheta2*costheta2);
2158 vec[i]->SetMomentum( pp*sintheta2*std::cos(phi2)*MeV,
2159 pp*sintheta2*std::sin(phi2)*MeV,
2160 pp*costheta2*MeV ) ;
2163 vec[i]->SetMomentum( vec[i]->GetMomentum() * (pp/pp1) );
2168 Rotate( numberofFinalStateNucleons, pseudoParticle[4].GetMomentum(),
2169 modifiedOriginal, originalIncident, targetNucleus,
2170 currentParticle, targetParticle, vec, vecLen );
2176 if( atomicWeight >= 1.5 )
2199 const G4double kineticMinimum = 1.e-6;
2200 const G4double kineticFactor = -0.005;
2205 if( ekIncident >= 5.0 )sprob = std::min( 1.0, 0.6*std::log(ekIncident-4.0) );
2207 if( epnb >= pnCutOff )
2209 npnb = Poisson((1.5+1.25*numberofFinalStateNucleons)*epnb/(epnb+edta));
2210 if( numberofFinalStateNucleons + npnb > atomicWeight )
2211 npnb =
G4int(atomicWeight - numberofFinalStateNucleons);
2212 npnb = std::min( npnb, 127-vecLen );
2214 if( edta >= dtaCutOff )
2216 ndta = Poisson( (1.5+1.25*numberofFinalStateNucleons)*edta/(epnb+edta) );
2217 ndta = std::min( ndta, 127-vecLen );
2219 if (npnb == 0 && ndta == 0) npnb = 1;
2223 AddBlackTrackParticles(epnb, npnb, edta, ndta, sprob, kineticMinimum,
2224 kineticFactor, modifiedOriginal,
2225 PinNucleus, NinNucleus, targetNucleus,
2234 if( (atomicWeight >= 1.5) && (atomicWeight <= 230.0) && (ekOriginal <= 0.2) )
2237 currentParticle.
SetTOF( 1.0 );
2265 static const G4double expxl = -expxu;
2274 targetMass = targetParticle.
GetMass()/GeV;
2280 G4double cmEnergy = std::sqrt( currentMass*currentMass +
2281 targetMass*targetMass +
2282 2.0*targetMass*etCurrent );
2290 if( (pCurrent < 0.1) || (cmEnergy < 0.01) )
2292 targetParticle.
SetMass( 0.0 );
2310 G4double pf = cmEnergy*cmEnergy + targetMass*targetMass - currentMass*currentMass;
2311 pf = pf*pf - 4*cmEnergy*cmEnergy*targetMass*targetMass;
2315 for(
G4int i=0; i<vecLen; i++)
delete vec[i];
2317 throw G4HadronicException(__FILE__, __LINE__,
"G4ReactionDynamics::TwoBody: pf is too small ");
2320 pf = std::sqrt( pf ) / ( 2.0*cmEnergy );
2328 pseudoParticle[0].
SetMass( targetMass*GeV );
2330 pseudoParticle[0].
SetMomentum( 0.0, 0.0, pOriginal*GeV );
2333 pseudoParticle[1].
SetMass( mOriginal*GeV );
2337 pseudoParticle[0].
SetMass( currentMass*GeV );
2339 pseudoParticle[0].
SetMomentum( 0.0, 0.0, pCurrent*GeV );
2342 pseudoParticle[1].
SetMass( targetMass*GeV );
2348 pseudoParticle[2] = pseudoParticle[0] + pseudoParticle[1];
2349 pseudoParticle[0].
Lorentz( pseudoParticle[0], pseudoParticle[2] );
2350 pseudoParticle[1].
Lorentz( pseudoParticle[1], pseudoParticle[2] );
2354 currentParticle.
SetTotalEnergy( std::sqrt(pf*pf+currentMass*currentMass)*GeV );
2355 targetParticle.
SetTotalEnergy( std::sqrt(pf*pf+targetMass*targetMass)*GeV );
2365 G4double b = std::max( cb, b1+b2*std::log(pOriginal) );
2369 exindt += std::exp(std::max(-btrang,expxl));
2374 if( std::fabs(ctet) > 1.0 )ctet > 0.0 ? ctet = 1.0 : ctet = -1.0;
2375 G4double stet = std::sqrt( (1.0-ctet)*(1.0+ctet) );
2383 currentParticle.
SetMomentum( -pf*stet*std::sin(phi)*GeV,
2384 -pf*stet*std::cos(phi)*GeV,
2388 currentParticle.
SetMomentum( pf*stet*std::sin(phi)*GeV,
2389 pf*stet*std::cos(phi)*GeV,
2396 currentParticle.
Lorentz( currentParticle, pseudoParticle[1] );
2397 targetParticle.
Lorentz( targetParticle, pseudoParticle[1] );
2399 Defs1( modifiedOriginal, currentParticle, targetParticle, vec, vecLen );
2402 if( atomicWeight >= 1.5 )
2404 const G4double cfa = 0.025*((atomicWeight-1.)/120.)*std::exp(-(atomicWeight-1.)/120.);
2408 ekin = currentParticle.
GetKineticEnergy()/MeV - cfa*(1.0+0.5*normal())*GeV;
2409 ekin = std::max( 0.0001*GeV, ekin );
2418 ekin = std::max( 0.0001*GeV, ekin );
2429 std::pair<G4int, G4int> finalStateNucleons =
2430 GetFinalStateNucleons(originalTarget, vec, vecLen);
2431 G4int protonsInFinalState = finalStateNucleons.first;
2432 G4int neutronsInFinalState = finalStateNucleons.second;
2434 G4int PinNucleus = std::max(0,
2435 targetNucleus.
GetZ_asInt() - protonsInFinalState);
2436 G4int NinNucleus = std::max(0,
2437 targetNucleus.
GetN_asInt() - neutronsInFinalState);
2440 if( atomicWeight >= 1.5 )
2449 G4int npnb=0, ndta=0;
2455 const G4double kineticMinimum = 0.0001;
2456 const G4double kineticFactor = -0.010;
2458 if( epnb >= pnCutOff )
2460 npnb = Poisson( epnb/0.02 );
2461 if( npnb > atomicWeight )npnb =
G4int(atomicWeight);
2462 if( (epnb > pnCutOff) && (npnb <= 0) )npnb = 1;
2463 npnb = std::min( npnb, 127-vecLen );
2465 if( edta >= dtaCutOff )
2467 ndta =
G4int(2.0 * std::log(atomicWeight));
2468 ndta = std::min( ndta, 127-vecLen );
2471 if (npnb == 0 && ndta == 0) npnb = 1;
2475 AddBlackTrackParticles(epnb, npnb, edta, ndta, sprob, kineticMinimum,
2476 kineticFactor, modifiedOriginal,
2477 PinNucleus, NinNucleus, targetNucleus,
2484 if( (atomicWeight >= 1.5) && (atomicWeight <= 230.0) && (ekOriginal <= 0.2) )
2487 currentParticle.
SetTOF( 1.0 );
2493 const G4bool constantCrossSection,
2506 G4cerr <<
"*** Error in G4ReactionDynamics::GenerateNBodyEvent" <<
G4endl;
2508 G4cerr <<
"totalEnergy = " << totalEnergy <<
"MeV, vecLen = " << vecLen <<
G4endl;
2519 for( i=0; i<vecLen; ++i )
2521 mass[i] = vec[i]->GetMass()/GeV;
2522 if(vec[i]->GetSide() == -2) extraMass+=vec[i]->GetMass()/GeV;
2523 vec[i]->SetMomentum( 0.0, 0.0, 0.0 );
2527 energy[i] = mass[i];
2528 totalMass += mass[i];
2532 if( totalMass > totalE )
2540 G4double kineticEnergy = totalE - totalMass;
2544 emm[vecLen-1] = totalE;
2549 for( i=0; i<vecLen-2; ++i )
2551 for(
G4int j=vecLen-2; j>i; --j )
2553 if( ran[i] > ran[j] )
2561 for( i=1; i<vecLen-1; ++i )emm[i] = ran[i-1]*kineticEnergy + sm[i];
2566 if( constantCrossSection )
2568 G4double emmax = kineticEnergy + mass[0];
2570 for( i=1; i<vecLen; ++i )
2575 if( emmax*emmax > 0.0 )
2578 + (emmin*emmin-mass[i]*mass[i])*(emmin*emmin-mass[i]*mass[i])/(emmax*emmax)
2579 - 2.0*(emmin*emmin+mass[i]*mass[i]);
2580 if( arg > 0.0 )wtfc = 0.5*std::sqrt( arg );
2587 wtmax += std::log( wtfc );
2597 const G4double ffq[18] = { 0., 3.141592, 19.73921, 62.01255, 129.8788, 204.0131,
2598 256.3704, 268.4705, 240.9780, 189.2637,
2599 132.1308, 83.0202, 47.4210, 24.8295,
2600 12.0006, 5.3858, 2.2560, 0.8859 };
2601 wtmax = std::log( std::pow( kineticEnergy, vecLen-2 ) * ffq[vecLen-1] / totalE );
2606 for( i=0; i<vecLen-1; ++i )
2609 if( emm[i+1]*emm[i+1] > 0.0 )
2612 + (emm[i]*emm[i]-mass[i+1]*mass[i+1])*(emm[i]*emm[i]-mass[i+1]*mass[i+1])
2613 /(emm[i+1]*emm[i+1])
2614 - 2.0*(emm[i]*emm[i]+mass[i+1]*mass[i+1]);
2615 if( arg > 0.0 )pd[i] = 0.5*std::sqrt( arg );
2620 wtmax += std::log( pd[i] );
2623 if( lzero )weight = std::exp( std::max(std::min(wtmax,expxu),expxl) );
2625 G4double bang, cb, sb, s0, s1, s2, c, ss, esys, a, b, gama, beta;
2629 for( i=1; i<vecLen; ++i )
2632 pcm[1][i] = -pd[i-1];
2635 cb = std::cos(bang);
2636 sb = std::sin(bang);
2638 ss = std::sqrt( std::fabs( 1.0-c*c ) );
2641 esys = std::sqrt(pd[i]*pd[i] + emm[i]*emm[i]);
2644 for(
G4int j=0; j<=i; ++j )
2649 energy[j] = std::sqrt( s0*s0 + s1*s1 + s2*s2 + mass[j]*mass[j] );
2651 pcm[1][j] = s0*ss + s1*c;
2653 pcm[0][j] = a*cb - b*sb;
2654 pcm[2][j] = a*sb + b*cb;
2655 pcm[1][j] = gama*(pcm[1][j] + beta*energy[j]);
2660 for(
G4int j=0; j<=i; ++j )
2665 energy[j] = std::sqrt( s0*s0 + s1*s1 + s2*s2 + mass[j]*mass[j] );
2667 pcm[1][j] = s0*ss + s1*c;
2669 pcm[0][j] = a*cb - b*sb;
2670 pcm[2][j] = a*sb + b*cb;
2674 for( i=0; i<vecLen; ++i )
2676 vec[i]->SetMomentum( pcm[0][i]*GeV, pcm[1][i]*GeV, pcm[2][i]*GeV );
2677 vec[i]->SetTotalEnergy( energy[i]*GeV );
2685 G4ReactionDynamics::normal()
2693 G4ReactionDynamics::Poisson(
G4double x )
2699 iran =
static_cast<G4int>(std::max( 0.0, x+normal()*std::sqrt(x) ) );
2726 for(
G4int i=1; i<=mn; ++i )
2730 rrr = std::exp(i*std::log(x)-(i+0.5)*std::log((
G4double)i)+i-0.9189385);
2734 if( ran <= rr )
break;
2745 G4int mn = std::min(n,10);
2747 if( mn <= 1 )
return result;
2748 for(
G4int i=2; i<=mn; ++i )result *= i;
2752 void G4ReactionDynamics::Defs1(
2763 if( pjx*pjx+pjy*pjy > 0.0 )
2765 G4double cost, sint, ph, cosp, sinp, pix, piy, piz;
2767 sint = 0.5 * ( std::sqrt(std::abs((1.0-cost)*(1.0+cost))) + std::sqrt(pjx*pjx+pjy*pjy)/p );
2772 if( std::abs( pjx ) > 0.001*MeV )ph = std::atan2(pjy,pjx);
2773 cosp = std::cos(ph);
2774 sinp = std::sin(ph);
2778 currentParticle.
SetMomentum( cost*cosp*pix*MeV - sinp*piy+sint*cosp*piz*MeV,
2779 cost*sinp*pix*MeV + cosp*piy+sint*sinp*piz*MeV,
2780 -sint*pix*MeV + cost*piz*MeV );
2784 targetParticle.
SetMomentum( cost*cosp*pix*MeV - sinp*piy+sint*cosp*piz*MeV,
2785 cost*sinp*pix*MeV + cosp*piy+sint*sinp*piz*MeV,
2786 -sint*pix*MeV + cost*piz*MeV );
2787 for(
G4int i=0; i<vecLen; ++i )
2789 pix = vec[i]->GetMomentum().x()/MeV;
2790 piy = vec[i]->GetMomentum().y()/MeV;
2791 piz = vec[i]->GetMomentum().z()/MeV;
2792 vec[i]->SetMomentum( cost*cosp*pix*MeV - sinp*piy+sint*cosp*piz*MeV,
2793 cost*sinp*pix*MeV + cosp*piy+sint*sinp*piz*MeV,
2794 -sint*pix*MeV + cost*piz*MeV );
2803 for(
G4int i=0; i<vecLen; ++i )
2804 vec[i]->SetMomentum( -vec[i]->GetMomentum().z() );
2809 void G4ReactionDynamics::Rotate(
2810 const G4double numberofFinalStateNucleons,
2826 const G4double logWeight = std::log(atomicWeight);
2834 for( i=0; i<4; ++i )pseudoParticle[i].set(0,0,0);
2837 for( i=0; i<vecLen; ++i )
2838 pseudoParticle[0] = pseudoParticle[0] + (vec[i]->GetMomentum());
2844 G4double r1, r2, a1, ran1, ran2, xxh, exh, pxTemp, pyTemp, pzTemp;
2848 a1 = std::sqrt(-2.0*std::log(r2));
2849 ran1 = a1*std::sin(r1)*0.020*numberofFinalStateNucleons*GeV;
2850 ran2 = a1*std::cos(r1)*0.020*numberofFinalStateNucleons*GeV;
2853 pseudoParticle[0] = pseudoParticle[0]+frm;
2854 pseudoParticle[2] = temp;
2855 pseudoParticle[3] = pseudoParticle[0];
2857 pseudoParticle[1] = pseudoParticle[2].
cross(pseudoParticle[3]);
2859 pseudoParticle[1] = pseudoParticle[1].
rotate(rotation, pseudoParticle[3]);
2860 pseudoParticle[2] = pseudoParticle[3].
cross(pseudoParticle[1]);
2861 for(
G4int ii=1; ii<=3; ii++)
2863 p = pseudoParticle[ii].
mag();
2867 pseudoParticle[ii]= pseudoParticle[ii] * (1./p);
2873 currentParticle.
SetMomentum( pxTemp, pyTemp, pzTemp );
2878 targetParticle.
SetMomentum( pxTemp, pyTemp, pzTemp );
2880 for( i=0; i<vecLen; ++i )
2882 pxTemp = pseudoParticle[1].
dot(vec[i]->GetMomentum());
2883 pyTemp = pseudoParticle[2].
dot(vec[i]->GetMomentum());
2884 pzTemp = pseudoParticle[3].
dot(vec[i]->GetMomentum());
2885 vec[i]->SetMomentum( pxTemp, pyTemp, pzTemp );
2891 Defs1( modifiedOriginal, currentParticle, targetParticle, vec, vecLen );
2896 if( atomicWeight >= 1.5 )
2900 const G4double alem[] = { 1.40, 2.30, 2.70, 3.00, 3.40, 4.60, 7.00 };
2901 const G4double val0[] = { 0.00, 0.40, 0.48, 0.51, 0.54, 0.60, 0.65 };
2904 if( alekw > alem[0] )
2907 for(
G4int j=1; j<7; ++j )
2909 if( alekw < alem[j] )
2911 G4double rcnve = (val0[j] - val0[j-1]) / (alem[j] - alem[j-1]);
2912 exh = rcnve * alekw + val0[j-1] - rcnve * alem[j-1];
2918 const G4double cfa = 0.025*((atomicWeight-1.)/120.)*std::exp(-(atomicWeight-1.)/120.);
2920 ekin = std::max( 1.0e-6, ekin );
2926 dekin += ekin*(1.0-xxh);
2935 if( pp1 < 0.001*MeV )
2938 G4double sintheta = std::sqrt(1. - costheta*costheta);
2940 currentParticle.
SetMomentum( pp*sintheta*std::cos(phi)*MeV,
2941 pp*sintheta*std::sin(phi)*MeV,
2947 ekin = std::max( 1.0e-6, ekin );
2953 dekin += ekin*(1.0-xxh);
2962 if( pp1 < 0.001*MeV )
2965 G4double sintheta = std::sqrt(1. - costheta*costheta);
2967 targetParticle.
SetMomentum( pp*sintheta*std::cos(phi)*MeV,
2968 pp*sintheta*std::sin(phi)*MeV,
2973 for( i=0; i<vecLen; ++i )
2975 ekin = vec[i]->GetKineticEnergy()/GeV - cfa*(1+normal()/2.0);
2976 ekin = std::max( 1.0e-6, ekin );
2980 vec[i]->GetDefinition() == aPiZero &&
2982 dekin += ekin*(1.0-xxh);
2984 if (vec[i]->GetDefinition()->GetParticleSubType() ==
"pi") {
2988 vec[i]->SetKineticEnergy( ekin*GeV );
2989 pp = vec[i]->GetTotalMomentum()/MeV;
2990 pp1 = vec[i]->GetMomentum().mag()/MeV;
2991 if( pp1 < 0.001*MeV )
2994 G4double sintheta = std::sqrt(1. - costheta*costheta);
2996 vec[i]->SetMomentum( pp*sintheta*std::cos(phi)*MeV,
2997 pp*sintheta*std::sin(phi)*MeV,
3001 vec[i]->SetMomentum( vec[i]->GetMomentum() * (pp/pp1) );
3004 if( (ek1 != 0.0) && (npions > 0) )
3006 dekin = 1.0 + dekin/ek1;
3019 G4double sintheta = std::sqrt(1. - costheta*costheta);
3021 currentParticle.
SetMomentum( pp*sintheta*std::cos(phi)*MeV,
3022 pp*sintheta*std::sin(phi)*MeV,
3038 G4double sintheta = std::sqrt(1. - costheta*costheta);
3040 targetParticle.
SetMomentum( pp*sintheta*std::cos(phi)*MeV,
3041 pp*sintheta*std::sin(phi)*MeV,
3048 for( i=0; i<vecLen; ++i )
3050 if (vec[i]->GetDefinition()->GetParticleSubType() ==
"pi")
3052 vec[i]->SetKineticEnergy( std::max( 0.001*MeV, dekin*vec[i]->GetKineticEnergy() ) );
3053 pp = vec[i]->GetTotalMomentum()/MeV;
3054 pp1 = vec[i]->GetMomentum().mag()/MeV;
3058 G4double sintheta = std::sqrt(1. - costheta*costheta);
3060 vec[i]->SetMomentum( pp*sintheta*std::cos(phi)*MeV,
3061 pp*sintheta*std::sin(phi)*MeV,
3064 vec[i]->SetMomentum( vec[i]->GetMomentum() * (pp/pp1) );
3071 void G4ReactionDynamics::AddBlackTrackParticles(
3111 G4double cfa = 0.025*((atomicWeight-1.0)/120.0) * std::exp(-(atomicWeight-1.0)/120.0);
3118 G4int local_npnb = npnb;
3119 for( i=0; i<npnb; ++i )
if(
G4UniformRand() < sprob ) local_npnb--;
3121 if (ndta == 0) local_epnb += edta;
3122 G4double ekin = local_epnb/std::max(1,local_npnb);
3124 for( i=0; i<local_npnb; ++i )
3127 if( backwardKinetic > local_epnb )
3133 G4double kinetic = -ekin*std::log(ran) - cfa*(1.0+0.5*normal());
3134 if( kinetic < 0.0 )kinetic = -0.010*std::log(ran);
3135 backwardKinetic += kinetic;
3136 if( backwardKinetic > local_epnb )
3137 kinetic = std::max( kineticMinimum, local_epnb-(backwardKinetic-kinetic) );
3143 if (PinNucleus > 0) {
3146 }
else if (NinNucleus > 0) {
3157 if (NinNucleus > 0) {
3160 }
else if (PinNucleus > 0) {
3171 G4double sint = std::sqrt(std::fabs(1.0-cost*cost));
3173 vec[vecLen]->SetNewlyAdded(
true );
3174 vec[vecLen]->SetKineticEnergy( kinetic*GeV );
3175 kinCreated+=kinetic;
3176 pp = vec[vecLen]->GetTotalMomentum()/MeV;
3177 vec[vecLen]->SetMomentum( pp*sint*std::sin(phi)*MeV,
3178 pp*sint*std::cos(phi)*MeV,
3184 if (NinNucleus > 0) {
3185 if( (atomicWeight >= 10.0) && (ekOriginal <= 2.0*GeV) )
3189 if( ekw > 1.0 )ekw *= ekw;
3190 ekw = std::max( 0.1, ekw );
3191 ika =
G4int(ika1*std::exp((atomicNumber*atomicNumber/
3192 atomicWeight-ika2)/ika3)/ekw);
3195 for( i=(vecLen-1); i>=0; --i )
3197 if( (vec[i]->GetDefinition() == aProton) && vec[i]->GetNewlyAdded() )
3199 vec[i]->SetDefinitionAndUpdateE( aNeutron );
3202 if( ++kk > ika )
break;
3215 G4int local_ndta=ndta;
3216 for( i=0; i<ndta; ++i )
if(
G4UniformRand() < sprob )local_ndta--;
3218 if (npnb == 0) local_edta += epnb;
3219 G4double ekin = local_edta/std::max(1,local_ndta);
3221 for( i=0; i<local_ndta; ++i )
3224 if( backwardKinetic > local_edta )
3230 G4double kinetic = -ekin*std::log(ran)-cfa*(1.+0.5*normal());
3231 if( kinetic < 0.0 )kinetic = kineticFactor*std::log(ran);
3232 backwardKinetic += kinetic;
3233 if( backwardKinetic > local_edta )kinetic = local_edta-(backwardKinetic-kinetic);
3234 if( kinetic < 0.0 )kinetic = kineticMinimum;
3236 G4double sint = std::sqrt(std::max(0.0,(1.0-cost*cost)));
3240 if (PinNucleus > 0 && NinNucleus > 0) {
3244 }
else if (NinNucleus > 0) {
3247 }
else if (PinNucleus > 0) {
3254 }
else if (ran < 0.90) {
3255 if (PinNucleus > 0 && NinNucleus > 1) {
3259 }
else if (PinNucleus > 0 && NinNucleus > 0) {
3263 }
else if (NinNucleus > 0) {
3266 }
else if (PinNucleus > 0) {
3274 if (PinNucleus > 1 && NinNucleus > 1) {
3278 }
else if (PinNucleus > 0 && NinNucleus > 1) {
3282 }
else if (PinNucleus > 0 && NinNucleus > 0) {
3286 }
else if (NinNucleus > 0) {
3289 }
else if (PinNucleus > 0) {
3299 vec[vecLen]->SetNewlyAdded(
true );
3300 vec[vecLen]->SetKineticEnergy( kinetic*GeV );
3301 kinCreated+=kinetic;
3302 pp = vec[vecLen]->GetTotalMomentum()/MeV;
3303 vec[vecLen++]->SetMomentum( pp*sint*std::sin(phi)*MeV,
3304 pp*sint*std::cos(phi)*MeV,
3314 std::pair<G4int, G4int> G4ReactionDynamics::GetFinalStateNucleons(
3317 const G4int& vecLen)
3321 G4int protonsRemoved = 0;
3322 G4int neutronsRemoved = 0;
3329 for (
G4int i = 0; i < vecLen; i++) {
3330 secName = vec[i]->GetDefinition()->GetParticleName();
3331 if (secName ==
"proton") {
3333 }
else if (secName ==
"neutron") {
3335 }
else if (secName ==
"anti_proton") {
3337 }
else if (secName ==
"anti_neutron") {
3342 return std::pair<G4int, G4int>(protonsRemoved, neutronsRemoved);
3346 void G4ReactionDynamics::MomentumCheck(
3356 if( testMomentum >= pOriginal )
3358 pMass = currentParticle.
GetMass()/MeV;
3360 std::sqrt( pMass*pMass + pOriginal*pOriginal )*MeV );
3362 currentParticle.
GetMomentum() * (pOriginal/testMomentum) );
3365 if( testMomentum >= pOriginal )
3367 pMass = targetParticle.
GetMass()/MeV;
3369 std::sqrt( pMass*pMass + pOriginal*pOriginal )*MeV );
3371 targetParticle.
GetMomentum() * (pOriginal/testMomentum) );
3373 for(
G4int i=0; i<vecLen; ++i )
3375 testMomentum = vec[i]->GetMomentum().mag()/MeV;
3376 if( testMomentum >= pOriginal )
3378 pMass = vec[i]->GetMass()/MeV;
3379 vec[i]->SetTotalEnergy(
3380 std::sqrt( pMass*pMass + pOriginal*pOriginal )*MeV );
3381 vec[i]->SetMomentum( vec[i]->GetMomentum() * (pOriginal/testMomentum) );
3393 G4bool &incidentHasChanged,
3394 G4bool &targetHasChanged )
3405 if( vecLen == 0 )
return;
3409 if( currentParticle.
GetMass() == 0.0 || targetParticle.
GetMass() == 0.0 )
return;
3414 G4double centerofmassEnergy = std::sqrt( mOriginal*mOriginal +
3415 targetMass*targetMass +
3416 2.0*targetMass*etOriginal );
3418 G4double availableEnergy = centerofmassEnergy-(targetMass+currentMass);
3419 if( availableEnergy <= 1.0 )
return;
3443 const G4double avrs[] = {3.,4.,5.,6.,7.,8.,9.,10.,20.,30.,40.,50.};
3448 while( (i<12) && (centerofmassEnergy>avrs[i]) )++i;
3466 i4 = i3 =
G4int( vecLen*ran );
3471 i4 =
G4int( vecLen*ran );
3477 const G4double avkkb[] = { 0.0015, 0.005, 0.012, 0.0285, 0.0525, 0.075,
3478 0.0975, 0.123, 0.28, 0.398, 0.495, 0.573 };
3479 const G4double avky[] = { 0.005, 0.03, 0.064, 0.095, 0.115, 0.13,
3480 0.145, 0.155, 0.20, 0.205, 0.210, 0.212 };
3481 const G4double avnnb[] = { 0.00001, 0.0001, 0.0006, 0.0025, 0.01, 0.02,
3482 0.04, 0.05, 0.12, 0.15, 0.18, 0.20 };
3484 avk = (std::log(avkkb[ibin])-std::log(avkkb[ibin-1]))*(centerofmassEnergy-avrs[ibin-1])
3485 /(avrs[ibin]-avrs[ibin-1]) + std::log(avkkb[ibin-1]);
3486 avk = std::exp(avk);
3488 avy = (std::log(avky[ibin])-std::log(avky[ibin-1]))*(centerofmassEnergy-avrs[ibin-1])
3489 /(avrs[ibin]-avrs[ibin-1]) + std::log(avky[ibin-1]);
3490 avy = std::exp(avy);
3492 avn = (std::log(avnnb[ibin])-std::log(avnnb[ibin-1]))*(centerofmassEnergy-avrs[ibin-1])
3493 /(avrs[ibin]-avrs[ibin-1]) + std::log(avnnb[ibin-1]);
3494 avn = std::exp(avn);
3496 if( avk+avy+avn <= 0.0 )
return;
3498 if( currentMass < protonMass )avy /= 2.0;
3499 if( targetMass < protonMass )avy = 0.0;
3505 if( availableEnergy < 2.0 )
return;
3511 vec[0]->SetDefinition( aNeutron );
3514 vec[0]->SetMayBeKilled(
false);
3519 vec[0]->SetDefinition( aProton );
3522 vec[0]->SetMayBeKilled(
false);
3532 vec[i3]->SetDefinition( aNeutron );
3533 vec[i4]->SetDefinition( anAntiNeutron );
3534 vec[i3]->SetMayBeKilled(
false);
3535 vec[i4]->SetMayBeKilled(
false);
3539 vec[i3]->SetDefinition( aProton );
3540 vec[i4]->SetDefinition( anAntiProton );
3541 vec[i3]->SetMayBeKilled(
false);
3542 vec[i4]->SetMayBeKilled(
false);
3546 else if( ran < avk )
3548 if( availableEnergy < 1.0 )
return;
3550 const G4double kkb[] = { 0.2500, 0.3750, 0.5000, 0.5625, 0.6250,
3551 0.6875, 0.7500, 0.8750, 1.000 };
3552 const G4int ipakkb1[] = { 10, 10, 10, 11, 11, 12, 12, 11, 12 };
3553 const G4int ipakkb2[] = { 13, 11, 12, 11, 12, 11, 12, 13, 13 };
3556 while( (i<9) && (ran>=kkb[i]) )++i;
3562 switch( ipakkb1[i] )
3565 vec[i3]->SetDefinition( aKaonPlus );
3566 vec[i3]->SetMayBeKilled(
false);
3569 vec[i3]->SetDefinition( aKaonZS );
3570 vec[i3]->SetMayBeKilled(
false);
3573 vec[i3]->SetDefinition( aKaonZL );
3574 vec[i3]->SetMayBeKilled(
false);
3580 switch( ipakkb2[i] )
3601 switch( ipakkb2[i] )
3604 vec[i4]->SetDefinition( aKaonZS );
3605 vec[i4]->SetMayBeKilled(
false);
3608 vec[i4]->SetDefinition( aKaonZL );
3609 vec[i4]->SetMayBeKilled(
false);
3612 vec[i4]->SetDefinition( aKaonMinus );
3613 vec[i4]->SetMayBeKilled(
false);
3618 else if( ran < avy )
3620 if( availableEnergy < 1.6 )
return;
3622 const G4double ky[] = { 0.200, 0.300, 0.400, 0.550, 0.625, 0.700,
3623 0.800, 0.850, 0.900, 0.950, 0.975, 1.000 };
3624 const G4int ipaky1[] = { 18, 18, 18, 20, 20, 20, 21, 21, 21, 22, 22, 22 };
3625 const G4int ipaky2[] = { 10, 11, 12, 10, 11, 12, 10, 11, 12, 10, 11, 12 };
3626 const G4int ipakyb1[] = { 19, 19, 19, 23, 23, 23, 24, 24, 24, 25, 25, 25 };
3627 const G4int ipakyb2[] = { 13, 12, 11, 13, 12, 11, 13, 12, 11, 13, 12, 11 };
3630 while( (i<12) && (ran>ky[i]) )++i;
3631 if( i == 12 )
return;
3654 targetHasChanged =
true;
3658 vec[i3]->SetDefinition( aKaonPlus );
3659 vec[i3]->SetMayBeKilled(
false);
3662 vec[i3]->SetDefinition( aKaonZS );
3663 vec[i3]->SetMayBeKilled(
false);
3666 vec[i3]->SetDefinition( aKaonZL );
3667 vec[i3]->SetMayBeKilled(
false);
3678 (currentMass > sigmaMinusMass) )
3680 switch( ipakyb1[i] )
3695 incidentHasChanged =
true;
3696 switch( ipakyb2[i] )
3699 vec[i3]->SetDefinition( aKaonZS );
3700 vec[i3]->SetMayBeKilled(
false);
3703 vec[i3]->SetDefinition( aKaonZL );
3704 vec[i3]->SetMayBeKilled(
false);
3707 vec[i3]->SetDefinition( aKaonMinus );
3708 vec[i3]->SetMayBeKilled(
false);
3729 incidentHasChanged =
true;
3733 vec[i3]->SetDefinition( aKaonPlus );
3734 vec[i3]->SetMayBeKilled(
false);
3737 vec[i3]->SetDefinition( aKaonZS );
3738 vec[i3]->SetMayBeKilled(
false);
3741 vec[i3]->SetDefinition( aKaonZL );
3742 vec[i3]->SetMayBeKilled(
false);
3758 currentMass = currentParticle.
GetMass()/GeV;
3759 targetMass = targetParticle.
GetMass()/GeV;
3761 G4double energyCheck = centerofmassEnergy-(currentMass+targetMass);
3762 for( i=0; i<vecLen; ++i )
3764 energyCheck -= vec[i]->GetMass()/GeV;
3765 if( energyCheck < 0.0 )
3767 vecLen = std::max( 0, --i );
3769 for(j=i; j<vecLen; j++)
delete vec[j];
3803 currentParticle = *originalIncident;
3811 if( pp <= 0.001*MeV )
3815 currentParticle.
SetMomentum( p*std::sin(rthnve)*std::cos(phinve),
3816 p*std::sin(rthnve)*std::sin(phinve),
3817 p*std::cos(rthnve) );
3826 G4double qv = currentKinetic + theAtomicMass + currentMass;
3829 qval[0] = qv - mass[0];
3830 qval[1] = qv - mass[1] - aNeutronMass;
3831 qval[2] = qv - mass[2] - aProtonMass;
3832 qval[3] = qv - mass[3] - aDeuteronMass;
3833 qval[4] = qv - mass[4] - aTritonMass;
3834 qval[5] = qv - mass[5] - anAlphaMass;
3835 qval[6] = qv - mass[6] - aNeutronMass - aNeutronMass;
3836 qval[7] = qv - mass[7] - aNeutronMass - aProtonMass;
3837 qval[8] = qv - mass[8] - aProtonMass - aProtonMass;
3845 qval[2] = qval[3] = qval[4] = qval[5] = qval[8] = 0.0;
3852 for( i=0; i<9; ++i )
3854 if( mass[i] < 500.0*MeV )qval[i] = 0.0;
3855 if( qval[i] < 0.0 )qval[i] = 0.0;
3861 for( index=0; index<9; ++index )
3863 if( qval[index] > 0.0 )
3865 qv1 += qval[index]/qv;
3866 if( ran <= qv1 )
break;
3872 "G4ReactionDynamics::NuclearReaction: inelastic reaction kinematically not possible");
3876 if( (index>=6) || (
G4UniformRand()<std::min(0.5,ke*10.0)) )nt = 3;
3883 v[0]->
SetMass( mass[index]*MeV );
3927 pseudo1.
SetMass( theAtomicMass*MeV );
3939 if( nt == 3 )tempV.
SetElement( tempLen++, v[2] );
3940 G4bool constantCrossSection =
true;
3942 v[0]->
Lorentz( *v[0], pseudo2 );
3943 v[1]->
Lorentz( *v[1], pseudo2 );
3944 if( nt == 3 )v[2]->
Lorentz( *v[2], pseudo2 );
3946 G4bool particleIsDefined =
false;
3947 if( v[0]->GetMass()/MeV - aProtonMass < 0.1 )
3950 particleIsDefined =
true;
3952 else if( v[0]->GetMass()/MeV - aNeutronMass < 0.1 )
3955 particleIsDefined =
true;
3957 else if( v[0]->GetMass()/MeV - aDeuteronMass < 0.1 )
3960 particleIsDefined =
true;
3962 else if( v[0]->GetMass()/MeV - aTritonMass < 0.1 )
3965 particleIsDefined =
true;
3967 else if( v[0]->GetMass()/MeV - anAlphaMass < 0.1 )
3970 particleIsDefined =
true;
3976 if( pp <= 0.001*MeV )
3980 currentParticle.
SetMomentum( p*std::sin(rthnve)*std::cos(phinve),
3981 p*std::sin(rthnve)*std::sin(phinve),
3982 p*std::cos(rthnve) );
3987 if( particleIsDefined )
3990 std::max( 0.001, 0.5*
G4UniformRand()*v[0]->GetKineticEnergy()/MeV ) );
3993 if( pp <= 0.001*MeV )
3997 v[0]->
SetMomentum( p*std::sin(rthnve)*std::cos(phinve),
3998 p*std::sin(rthnve)*std::sin(phinve),
3999 p*std::cos(rthnve) );
4002 v[0]->
SetMomentum( v[0]->GetMomentum() * (p/pp) );
4004 if( (v[1]->GetDefinition() == aDeuteron) ||
4005 (v[1]->GetDefinition() == aTriton) ||
4006 (v[1]->GetDefinition() == anAlpha) )
4008 std::max( 0.001, 0.5*
G4UniformRand()*v[1]->GetKineticEnergy()/MeV ) );
4014 if( pp <= 0.001*MeV )
4018 v[1]->
SetMomentum( p*std::sin(rthnve)*std::cos(phinve),
4019 p*std::sin(rthnve)*std::sin(phinve),
4020 p*std::cos(rthnve) );
4023 v[1]->
SetMomentum( v[1]->GetMomentum() * (p/pp) );
4027 if( (v[2]->GetDefinition() == aDeuteron) ||
4028 (v[2]->GetDefinition() == aTriton) ||
4029 (v[2]->GetDefinition() == anAlpha) )
4031 std::max( 0.001, 0.5*
G4UniformRand()*v[2]->GetKineticEnergy()/MeV ) );
4037 if( pp <= 0.001*MeV )
4041 v[2]->
SetMomentum( p*std::sin(rthnve)*std::cos(phinve),
4042 p*std::sin(rthnve)*std::sin(phinve),
4043 p*std::cos(rthnve) );
4046 v[2]->
SetMomentum( v[2]->GetMomentum() * (p/pp) );
4049 for(del=0; del<vecLen; del++)
delete vec[del];
4051 if( particleIsDefined )
CLHEP::Hep3Vector G4ThreeVector
G4DLLIMPORT std::ostream G4cerr
Hep3Vector cross(const Hep3Vector &) const
double dot(const Hep3Vector &) const
Hep3Vector & rotate(double, const Hep3Vector &)
static G4AntiLambda * AntiLambda()
static G4AntiNeutron * AntiNeutron()
static G4AntiOmegaMinus * AntiOmegaMinus()
static G4AntiProton * AntiProton()
static G4AntiSigmaMinus * AntiSigmaMinus()
static G4AntiSigmaPlus * AntiSigmaPlus()
static G4AntiSigmaZero * AntiSigmaZero()
static G4AntiXiMinus * AntiXiMinus()
static G4AntiXiZero * AntiXiZero()
static G4Deuteron * Deuteron()
G4ParticleDefinition * GetDefinition() const
void SetElement(G4int anIndex, Type *anElement)
void Initialize(G4int items)
const G4ParticleDefinition * GetDefinition() const
G4double GetKineticEnergy() const
static G4KaonMinus * KaonMinus()
static G4KaonPlus * KaonPlus()
static G4KaonZeroLong * KaonZeroLong()
static G4KaonZeroShort * KaonZeroShort()
static G4Lambda * Lambda()
static G4Neutron * Neutron()
G4double GetAnnihilationPNBlackTrackEnergy() const
G4double GetPNBlackTrackEnergy() const
G4double GetAnnihilationDTABlackTrackEnergy() const
G4double GetDTABlackTrackEnergy() const
G4double GetPDGMass() const
G4int GetBaryonNumber() const
const G4String & GetParticleName() const
const G4String & GetParticleSubType() const
static G4PionMinus * PionMinus()
static G4PionPlus * PionPlus()
static G4PionZero * PionZero()
static G4Proton * Proton()
G4bool TwoCluster(G4FastVector< G4ReactionProduct, GHADLISTSIZE > &vec, G4int &vecLen, G4ReactionProduct &modifiedOriginal, const G4HadProjectile *originalIncident, G4ReactionProduct ¤tParticle, G4ReactionProduct &targetParticle, const G4DynamicParticle *originalTarget, const G4Nucleus &targetNucleus, G4bool &incidentHasChanged, G4bool &targetHasChanged, G4bool leadFlag, G4ReactionProduct &leadingStrangeParticle)
void SuppressChargedPions(G4FastVector< G4ReactionProduct, GHADLISTSIZE > &vec, G4int &vecLen, const G4ReactionProduct &modifiedOriginal, G4ReactionProduct ¤tParticle, G4ReactionProduct &targetParticle, const G4Nucleus &targetNucleus, G4bool &incidentHasChanged, G4bool &targetHasChanged)
void ProduceStrangeParticlePairs(G4FastVector< G4ReactionProduct, GHADLISTSIZE > &vec, G4int &vecLen, const G4ReactionProduct &modifiedOriginal, const G4DynamicParticle *originalTarget, G4ReactionProduct ¤tParticle, G4ReactionProduct &targetParticle, G4bool &incidentHasChanged, G4bool &targetHasChanged)
G4double GenerateNBodyEvent(const G4double totalEnergy, const G4bool constantCrossSection, G4FastVector< G4ReactionProduct, GHADLISTSIZE > &vec, G4int &vecLen)
G4bool GenerateXandPt(G4FastVector< G4ReactionProduct, GHADLISTSIZE > &vec, G4int &vecLen, G4ReactionProduct &modifiedOriginal, const G4HadProjectile *originalIncident, G4ReactionProduct ¤tParticle, G4ReactionProduct &targetParticle, const G4DynamicParticle *originalTarget, const G4Nucleus &targetNucleus, G4bool &incidentHasChanged, G4bool &targetHasChanged, G4bool leadFlag, G4ReactionProduct &leadingStrangeParticle)
void TwoBody(G4FastVector< G4ReactionProduct, GHADLISTSIZE > &vec, G4int &vecLen, G4ReactionProduct &modifiedOriginal, const G4DynamicParticle *originalTarget, G4ReactionProduct ¤tParticle, G4ReactionProduct &targetParticle, const G4Nucleus &targetNucleus, G4bool &targetHasChanged)
void NuclearReaction(G4FastVector< G4ReactionProduct, 4 > &vec, G4int &vecLen, const G4HadProjectile *originalIncident, const G4Nucleus &aNucleus, const G4double theAtomicMass, const G4double *massVec)
void SetMomentum(const G4double x, const G4double y, const G4double z)
void SetTotalEnergy(const G4double en)
G4double GetTotalMomentum() const
G4double GetKineticEnergy() const
G4double GetTotalEnergy() const
void SetMayBeKilled(const G4bool f)
void SetNewlyAdded(const G4bool f)
G4ThreeVector GetMomentum() const
void SetSide(const G4int sid)
void SetDefinitionAndUpdateE(G4ParticleDefinition *aParticleDefinition)
void Lorentz(const G4ReactionProduct &p1, const G4ReactionProduct &p2)
void SetKineticEnergy(const G4double en)
G4ParticleDefinition * GetDefinition() const
void SetTOF(const G4double t)
void SetDefinition(G4ParticleDefinition *aParticleDefinition)
void SetMass(const G4double mas)
static G4SigmaMinus * SigmaMinus()
static G4SigmaPlus * SigmaPlus()
static G4SigmaZero * SigmaZero()
static G4Triton * Triton()
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)