47 G4bool& incidentHasChanged,
70 G4bool veryForward =
false;
78 G4double centerofmassEnergy = std::sqrt(mOriginal*mOriginal +
79 targetMass*targetMass +
80 2.0*targetMass*etOriginal);
84 if (currentMass == 0.0 && targetMass == 0.0) {
87 currentParticle = *vec[0];
88 targetParticle = *vec[1];
89 for (i = 0; i < (vecLen-2); ++i) *vec[i] = *vec[i+2];
91 for (
G4int j = 0; j < vecLen; j++)
delete vec[j];
94 "G4RPGTwoCluster::ReactionStage : Negative number of particles");
101 incidentHasChanged =
true;
102 targetHasChanged =
true;
116 G4int forwardCount = 1;
122 G4int backwardCount = 1;
128 for (i = 0; i < vecLen; ++i) {
129 if (vec[i]->GetSide() < 0) vec[i]->SetSide(-1);
132 if (vec[i]->GetSide() == -1) {
134 backwardMass += vec[i]->GetMass()/
GeV;
137 forwardMass += vec[i]->GetMass()/
GeV;
142 G4double term1 =
G4Log(centerofmassEnergy*centerofmassEnergy);
143 if(term1 < 0) term1 = 0.0001;
148 xtarg = afc * (a13-1.0) * (2*backwardCount+vecLen+2)/2.0;
150 xtarg = afc * (a13-1.0) * (2*backwardCount);
152 if( xtarg <= 0.0 )xtarg = 0.01;
155 if(atomicWeight<1.0001) nuclearExcitationCount = 0;
159 if( nuclearExcitationCount > 0 )
162 const G4double nucsup[] = { 1.0, 0.8, 0.6, 0.5, 0.4 };
167 for( i=0; i<nuclearExcitationCount; ++i )
185 else if( ran < 0.6819 )
205 G4double forwardEnergy = (centerofmassEnergy-cMass-bMass)/2.0 +cMass - forwardMass;
206 G4double backwardEnergy = (centerofmassEnergy-cMass-bMass)/2.0 +bMass - backwardMass;
207 G4double eAvailable = centerofmassEnergy - (forwardMass+backwardMass);
213 ed <<
" While count exceeded " <<
G4endl;
215 while( eAvailable <= 0.0 ) {
222 secondaryDeleted =
false;
223 for( i=(vecLen-1); i>=0; --i )
225 if( vec[i]->GetSide() == 1 && vec[i]->GetMayBeKilled())
227 pMass = vec[i]->GetMass()/
GeV;
228 for(
G4int j=i; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];
230 forwardEnergy += pMass;
231 forwardMass -= pMass;
232 secondaryDeleted =
true;
235 else if( vec[i]->GetSide() == -1 && vec[i]->GetMayBeKilled())
237 pMass = vec[i]->GetMass()/
GeV;
238 for(
G4int j=i; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];
240 backwardEnergy += pMass;
241 backwardMass -= pMass;
242 secondaryDeleted =
true;
247 if( secondaryDeleted )
249 delete vec[vecLen-1];
255 if( vecLen == 0 )
return false;
256 if( targetParticle.
GetSide() == -1 )
259 targetParticle = *vec[0];
260 for(
G4int j=0; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];
262 backwardEnergy += pMass;
263 backwardMass -= pMass;
264 secondaryDeleted =
true;
266 else if( targetParticle.
GetSide() == 1 )
269 targetParticle = *vec[0];
270 for(
G4int j=0; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];
272 forwardEnergy += pMass;
273 forwardMass -= pMass;
274 secondaryDeleted =
true;
277 if( secondaryDeleted )
279 delete vec[vecLen-1];
284 if( currentParticle.
GetSide() == -1 )
287 currentParticle = *vec[0];
288 for(
G4int j=0; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];
290 backwardEnergy += pMass;
291 backwardMass -= pMass;
292 secondaryDeleted =
true;
294 else if( currentParticle.
GetSide() == 1 )
297 currentParticle = *vec[0];
298 for(
G4int j=0; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];
300 forwardEnergy += pMass;
301 forwardMass -= pMass;
302 secondaryDeleted =
true;
304 if( secondaryDeleted )
306 delete vec[vecLen-1];
314 eAvailable = centerofmassEnergy - (forwardMass+backwardMass);
326 const G4double cpar[] = { 1.60, 1.35, 1.15, 1.10 };
327 const G4double gpar[] = { 2.60, 1.80, 1.30, 1.20 };
330 if (forwardCount < 1 || backwardCount < 1)
return false;
333 if (forwardCount > 1) {
338 if( backwardCount > 1 ) {
345 eda <<
" While count exceeded " <<
G4endl;
346 while( rmc+rmd > centerofmassEnergy ) {
353 if( (rmc <= forwardMass) && (rmd <= backwardMass) )
355 G4double temp = 0.999*centerofmassEnergy/(rmc+rmd);
361 rmc = 0.1*forwardMass + 0.9*rmc;
362 rmd = 0.1*backwardMass + 0.9*rmd;
367 for( i=0; i<8; ++i )pseudoParticle[i].SetZero();
371 pseudoParticle[1].
SetMomentum( 0.0, 0.0, pOriginal*GeV );
379 pseudoParticle[0] = pseudoParticle[1] + pseudoParticle[2];
380 pseudoParticle[1].
Lorentz( pseudoParticle[1], pseudoParticle[0] );
381 pseudoParticle[2].Lorentz( pseudoParticle[2], pseudoParticle[0] );
388 G4double pf = (centerofmassEnergy*centerofmassEnergy+rmd*rmd-rmc*rmc);
390 pf -= 4*centerofmassEnergy*centerofmassEnergy*rmd*rmd;
391 pf = std::sqrt(
std::max(pf,pfMin) )/(2.0*centerofmassEnergy);
395 pseudoParticle[3].SetMass( rmc*GeV );
396 pseudoParticle[3].SetTotalEnergy( std::sqrt(pf*pf+rmc*rmc)*GeV );
398 pseudoParticle[4].SetMass( rmd*GeV );
399 pseudoParticle[4].SetTotalEnergy( std::sqrt(pf*pf+rmd*rmd)*GeV );
407 G4double pin = pseudoParticle[1].GetMomentum().mag()/
GeV;
413 G4double sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
418 pseudoParticle[3].SetMomentum( pf*sintheta*std::cos(phi)*GeV,
419 pf*sintheta*std::sin(phi)*GeV,
421 pseudoParticle[4].SetMomentum( -pseudoParticle[3].GetMomentum());
427 if( nuclearExcitationCount > 0 )
432 if( ekOriginal <= 5.0 )
434 ekit1 *= ekOriginal*ekOriginal/25.0;
435 ekit2 *= ekOriginal*ekOriginal/25.0;
438 for( i=0; i<vecLen; ++i )
440 if( vec[i]->GetSide() == -2 )
443 vec[i]->SetKineticEnergy( kineticE*GeV );
445 G4double totalE = kineticE*GeV + vMass;
446 pp = std::sqrt(
std::abs(totalE*totalE-vMass*vMass) );
448 G4double sint = std::sqrt(1.0-cost*cost);
450 vec[i]->SetMomentum( pp*sint*std::cos(phi)*MeV,
451 pp*sint*std::sin(phi)*MeV,
453 vec[i]->Lorentz( *vec[i], pseudoParticle[0] );
462 currentParticle.
SetMomentum( pseudoParticle[3].GetMomentum() );
463 currentParticle.
SetTotalEnergy( pseudoParticle[3].GetTotalEnergy() );
465 targetParticle.
SetMomentum( pseudoParticle[4].GetMomentum() );
466 targetParticle.
SetTotalEnergy( pseudoParticle[4].GetTotalEnergy() );
468 pseudoParticle[5].SetMomentum( pseudoParticle[3].GetMomentum() * (-1.0) );
469 pseudoParticle[5].SetMass( pseudoParticle[3].GetMass() );
470 pseudoParticle[5].SetTotalEnergy( pseudoParticle[3].GetTotalEnergy() );
472 pseudoParticle[6].SetMomentum( pseudoParticle[4].GetMomentum() * (-1.0) );
473 pseudoParticle[6].SetMass( pseudoParticle[4].GetMass() );
474 pseudoParticle[6].SetTotalEnergy( pseudoParticle[4].GetTotalEnergy() );
478 if( forwardCount > 1 )
482 G4bool constantCrossSection =
true;
484 if( currentParticle.
GetSide() == 1 )
485 tempV.
SetElement( tempLen++, ¤tParticle );
486 if( targetParticle.
GetSide() == 1 )
487 tempV.
SetElement( tempLen++, &targetParticle );
488 for( i=0; i<vecLen; ++i )
490 if( vec[i]->GetSide() == 1 )
496 vec[i]->SetSide( -1 );
504 constantCrossSection, tempV, tempLen );
505 if( currentParticle.
GetSide() == 1 )
506 currentParticle.
Lorentz( currentParticle, pseudoParticle[5] );
507 if( targetParticle.
GetSide() == 1 )
508 targetParticle.
Lorentz( targetParticle, pseudoParticle[5] );
509 for( i=0; i<vecLen; ++i )
511 if( vec[i]->GetSide() == 1 )vec[i]->Lorentz( *vec[i], pseudoParticle[5] );
516 if( backwardCount > 1 )
520 G4bool constantCrossSection =
true;
522 if( currentParticle.
GetSide() == -1 )
523 tempV.
SetElement( tempLen++, ¤tParticle );
524 if( targetParticle.
GetSide() == -1 )
525 tempV.
SetElement( tempLen++, &targetParticle );
526 for( i=0; i<vecLen; ++i )
528 if( vec[i]->GetSide() == -1 )
534 vec[i]->SetSide( -2 );
535 vec[i]->SetKineticEnergy( 0.0 );
536 vec[i]->SetMomentum( 0.0, 0.0, 0.0 );
544 constantCrossSection, tempV, tempLen );
545 if( currentParticle.
GetSide() == -1 )
546 currentParticle.
Lorentz( currentParticle, pseudoParticle[6] );
547 if( targetParticle.
GetSide() == -1 )
548 targetParticle.
Lorentz( targetParticle, pseudoParticle[6] );
549 for( i=0; i<vecLen; ++i )
551 if( vec[i]->GetSide() == -1 )vec[i]->Lorentz( *vec[i], pseudoParticle[6] );
560 currentParticle.
Lorentz( currentParticle, pseudoParticle[2] );
561 targetParticle.
Lorentz( targetParticle, pseudoParticle[2] );
562 for( i=0; i<vecLen; ++i ) vec[i]->Lorentz( *vec[i], pseudoParticle[2] );
583 for( i=0; i<vecLen; ++i )
585 if( vec[i]->GetDefinition() == leadingStrangeParticle.
GetDefinition() )
596 if( ((leadMass < protonMass) && (targetParticle.
GetMass()/MeV < protonMass)) ||
597 ((leadMass >= protonMass) && (targetParticle.
GetMass()/MeV >= protonMass)) )
610 targetHasChanged =
true;
625 incidentHasChanged =
true;
633 std::pair<G4int, G4int> finalStateNucleons =
636 G4int protonsInFinalState = finalStateNucleons.first;
637 G4int neutronsInFinalState = finalStateNucleons.second;
639 G4int numberofFinalStateNucleons =
640 protonsInFinalState + neutronsInFinalState;
646 numberofFinalStateNucleons++;
648 numberofFinalStateNucleons =
std::max(1, numberofFinalStateNucleons);
658 pseudoParticle[4].SetMass( mOriginal*GeV );
659 pseudoParticle[4].SetTotalEnergy( etOriginal*GeV );
660 pseudoParticle[4].SetMomentum( 0.0, 0.0, pOriginal*GeV );
665 if(numberofFinalStateNucleons == 1) diff = 0;
666 pseudoParticle[5].SetMomentum( 0.0, 0.0, 0.0 );
667 pseudoParticle[5].SetMass( protonMass*(numberofFinalStateNucleons-diff)*MeV);
668 pseudoParticle[5].SetTotalEnergy( protonMass*(numberofFinalStateNucleons-diff)*MeV);
671 pseudoParticle[4].GetTotalEnergy()/GeV + pseudoParticle[5].GetTotalEnergy()/
GeV;
673 pseudoParticle[6] = pseudoParticle[4] + pseudoParticle[5];
674 pseudoParticle[4].Lorentz( pseudoParticle[4], pseudoParticle[6] );
675 pseudoParticle[5].Lorentz( pseudoParticle[5], pseudoParticle[6] );
680 tempR[0] = currentParticle;
681 tempR[1] = targetParticle;
682 for( i=0; i<vecLen; ++i )tempR[i+2] = *vec[i];
686 G4bool constantCrossSection =
true;
688 for( i=0; i<vecLen+2; ++i )tempV.
SetElement( tempLen++, &tempR[i] );
694 pseudoParticle[5].GetTotalEnergy()/MeV,
695 constantCrossSection, tempV, tempLen );
698 for (i = 0; i < tempLen; i++) Qvalue += tempV[i]->GetMass();
700 constantCrossSection, tempV, tempLen );
702 theoreticalKinetic = 0.0;
703 for( i=0; i<vecLen+2; ++i )
705 pseudoParticle[7].SetMomentum( tempV[i]->GetMomentum() );
706 pseudoParticle[7].SetMass( tempV[i]->GetMass() );
707 pseudoParticle[7].SetTotalEnergy( tempV[i]->GetTotalEnergy() );
708 pseudoParticle[7].Lorentz( pseudoParticle[7], pseudoParticle[5] );
709 theoreticalKinetic += pseudoParticle[7].GetKineticEnergy()/
GeV;
716 theoreticalKinetic -=
718 for( i=0; i<vecLen; ++i )theoreticalKinetic -= vec[i]->GetMass()/
GeV;
722 for( i=0; i<vecLen; ++i )simulatedKinetic += vec[i]->GetKineticEnergy()/
GeV;
727 if( simulatedKinetic != 0.0 )
729 wgt = (theoreticalKinetic)/simulatedKinetic;
733 if( pp1 < 0.001*MeV ) {
743 if( pp1 < 0.001*MeV ) {
750 for( i=0; i<vecLen; ++i )
752 vec[i]->SetKineticEnergy( wgt*vec[i]->GetKineticEnergy() );
753 pp = vec[i]->GetTotalMomentum()/
MeV;
754 pp1 = vec[i]->GetMomentum().mag()/
MeV;
757 vec[i]->SetMomentum( iso.
x(), iso.
y(), iso.
z() );
759 vec[i]->SetMomentum( vec[i]->GetMomentum() * (pp/pp1) );
765 Rotate( numberofFinalStateNucleons, pseudoParticle[4].GetMomentum(),
766 modifiedOriginal, originalIncident, targetNucleus,
767 currentParticle, targetParticle, vec, vecLen );
773 if( atomicWeight >= 1.5 )
804 if( epnb >= pnCutOff )
806 npnb =
G4Poisson((1.5+1.25*numberofFinalStateNucleons)*epnb/(epnb+edta));
807 if( numberofFinalStateNucleons + npnb > atomicWeight )
808 npnb =
G4int(atomicWeight - numberofFinalStateNucleons);
809 npnb =
std::min( npnb, 127-vecLen );
811 if( edta >= dtaCutOff )
813 ndta =
G4Poisson( (1.5+1.25*numberofFinalStateNucleons)*edta/(epnb+edta) );
814 ndta =
std::min( ndta, 127-vecLen );
816 if (npnb == 0 && ndta == 0) npnb = 1;
821 PinNucleus, NinNucleus, targetNucleus,
831 if( (atomicWeight >= 1.5) && (atomicWeight <= 230.0) && (ekOriginal <= 0.2) )
834 currentParticle.
SetTOF( 1.0 );