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 #include <iostream>
00030 #include <signal.h>
00031
00032 #include "G4RPGTwoCluster.hh"
00033 #include "G4PhysicalConstants.hh"
00034 #include "G4SystemOfUnits.hh"
00035 #include "Randomize.hh"
00036 #include "G4Poisson.hh"
00037 #include "G4HadReentrentException.hh"
00038
00039 G4RPGTwoCluster::G4RPGTwoCluster()
00040 : G4RPGReaction() {}
00041
00042
00043 G4bool G4RPGTwoCluster::
00044 ReactionStage(const G4HadProjectile* originalIncident,
00045 G4ReactionProduct& modifiedOriginal,
00046 G4bool& incidentHasChanged,
00047 const G4DynamicParticle* originalTarget,
00048 G4ReactionProduct& targetParticle,
00049 G4bool& targetHasChanged,
00050 const G4Nucleus& targetNucleus,
00051 G4ReactionProduct& currentParticle,
00052 G4FastVector<G4ReactionProduct,256>& vec,
00053 G4int& vecLen,
00054 G4bool leadFlag,
00055 G4ReactionProduct& leadingStrangeParticle)
00056 {
00057
00058
00059
00060
00061
00062
00063 G4int i;
00064 G4ParticleDefinition* aProton = G4Proton::Proton();
00065 G4ParticleDefinition* aNeutron = G4Neutron::Neutron();
00066 G4ParticleDefinition* aPiPlus = G4PionPlus::PionPlus();
00067 G4ParticleDefinition* aPiMinus = G4PionMinus::PionMinus();
00068 G4ParticleDefinition* aPiZero = G4PionZero::PionZero();
00069 G4bool veryForward = false;
00070
00071 const G4double protonMass = aProton->GetPDGMass()/MeV;
00072 const G4double ekOriginal = modifiedOriginal.GetKineticEnergy()/GeV;
00073 const G4double etOriginal = modifiedOriginal.GetTotalEnergy()/GeV;
00074 const G4double mOriginal = modifiedOriginal.GetMass()/GeV;
00075 const G4double pOriginal = modifiedOriginal.GetMomentum().mag()/GeV;
00076 G4double targetMass = targetParticle.GetDefinition()->GetPDGMass()/GeV;
00077 G4double centerofmassEnergy = std::sqrt(mOriginal*mOriginal +
00078 targetMass*targetMass +
00079 2.0*targetMass*etOriginal);
00080 G4double currentMass = currentParticle.GetMass()/GeV;
00081 targetMass = targetParticle.GetMass()/GeV;
00082
00083 if (currentMass == 0.0 && targetMass == 0.0) {
00084 G4double ek = currentParticle.GetKineticEnergy();
00085 G4ThreeVector mom = currentParticle.GetMomentum();
00086 currentParticle = *vec[0];
00087 targetParticle = *vec[1];
00088 for (i = 0; i < (vecLen-2); ++i) *vec[i] = *vec[i+2];
00089 if (vecLen < 2) {
00090 for (G4int j = 0; j < vecLen; j++) delete vec[j];
00091 vecLen = 0;
00092 throw G4HadReentrentException(__FILE__, __LINE__,
00093 "G4RPGTwoCluster::ReactionStage : Negative number of particles");
00094 }
00095 delete vec[vecLen-1];
00096 delete vec[vecLen-2];
00097 vecLen -= 2;
00098 currentMass = currentParticle.GetMass()/GeV;
00099 targetMass = targetParticle.GetMass()/GeV;
00100 incidentHasChanged = true;
00101 targetHasChanged = true;
00102 currentParticle.SetKineticEnergy(ek);
00103 currentParticle.SetMomentum(mom);
00104 veryForward = true;
00105 }
00106
00107 const G4double atomicWeight = targetNucleus.GetA_asInt();
00108 const G4double atomicNumber = targetNucleus.GetZ_asInt();
00109
00110
00111
00112
00113
00114
00115 G4int forwardCount = 1;
00116 currentParticle.SetSide( 1 );
00117 G4double forwardMass = currentParticle.GetMass()/GeV;
00118 G4double cMass = forwardMass;
00119
00120
00121 G4int backwardCount = 1;
00122 targetParticle.SetSide( -1 );
00123 G4double backwardMass = targetParticle.GetMass()/GeV;
00124 G4double bMass = backwardMass;
00125
00126
00127 for (i = 0; i < vecLen; ++i) {
00128 if (vec[i]->GetSide() < 0) vec[i]->SetSide(-1);
00129
00130
00131 if (vec[i]->GetSide() == -1) {
00132 ++backwardCount;
00133 backwardMass += vec[i]->GetMass()/GeV;
00134 } else {
00135 ++forwardCount;
00136 forwardMass += vec[i]->GetMass()/GeV;
00137 }
00138 }
00139
00140
00141 G4double term1 = std::log(centerofmassEnergy*centerofmassEnergy);
00142 if(term1 < 0) term1 = 0.0001;
00143 const G4double afc = 0.312 + 0.2 * std::log(term1);
00144 G4double xtarg;
00145 if( centerofmassEnergy < 2.0+G4UniformRand() )
00146 xtarg = afc * (std::pow(atomicWeight,0.33)-1.0) * (2*backwardCount+vecLen+2)/2.0;
00147 else
00148 xtarg = afc * (std::pow(atomicWeight,0.33)-1.0) * (2*backwardCount);
00149 if( xtarg <= 0.0 )xtarg = 0.01;
00150 G4int nuclearExcitationCount = G4Poisson( xtarg );
00151
00152 if(atomicWeight<1.0001) nuclearExcitationCount = 0;
00153
00154
00155
00156 if( nuclearExcitationCount > 0 )
00157 {
00158 G4int momentumBin = std::min( 4, G4int(pOriginal/3.0) );
00159 const G4double nucsup[] = { 1.0, 0.8, 0.6, 0.5, 0.4 };
00160
00161
00162
00163
00164 for( i=0; i<nuclearExcitationCount; ++i )
00165 {
00166 G4ReactionProduct* pVec = new G4ReactionProduct();
00167 if( G4UniformRand() < nucsup[momentumBin] )
00168 {
00169 if( G4UniformRand() > 1.0-atomicNumber/atomicWeight )
00170 pVec->SetDefinition( aProton );
00171 else
00172 pVec->SetDefinition( aNeutron );
00173
00174
00175
00176 }
00177 else
00178 {
00179 G4double ran = G4UniformRand();
00180 if( ran < 0.3181 )
00181 pVec->SetDefinition( aPiPlus );
00182 else if( ran < 0.6819 )
00183 pVec->SetDefinition( aPiZero );
00184 else
00185 pVec->SetDefinition( aPiMinus );
00186
00187
00188
00189
00190 }
00191 pVec->SetSide( -2 );
00192
00193 pVec->SetNewlyAdded( true );
00194 vec.SetElement( vecLen++, pVec );
00195 }
00196 }
00197
00198
00199
00200
00201
00202 G4double forwardEnergy = (centerofmassEnergy-cMass-bMass)/2.0 +cMass - forwardMass;
00203 G4double backwardEnergy = (centerofmassEnergy-cMass-bMass)/2.0 +bMass - backwardMass;
00204 G4double eAvailable = centerofmassEnergy - (forwardMass+backwardMass);
00205 G4bool secondaryDeleted;
00206 G4double pMass;
00207
00208 while( eAvailable <= 0.0 )
00209 {
00210 secondaryDeleted = false;
00211 for( i=(vecLen-1); i>=0; --i )
00212 {
00213 if( vec[i]->GetSide() == 1 && vec[i]->GetMayBeKilled())
00214 {
00215 pMass = vec[i]->GetMass()/GeV;
00216 for( G4int j=i; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];
00217 --forwardCount;
00218 forwardEnergy += pMass;
00219 forwardMass -= pMass;
00220 secondaryDeleted = true;
00221 break;
00222 }
00223 else if( vec[i]->GetSide() == -1 && vec[i]->GetMayBeKilled())
00224 {
00225 pMass = vec[i]->GetMass()/GeV;
00226 for( G4int j=i; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];
00227 --backwardCount;
00228 backwardEnergy += pMass;
00229 backwardMass -= pMass;
00230 secondaryDeleted = true;
00231 break;
00232 }
00233 }
00234
00235 if( secondaryDeleted )
00236 {
00237 delete vec[vecLen-1];
00238 --vecLen;
00239
00240 }
00241 else
00242 {
00243 if( vecLen == 0 ) return false;
00244 if( targetParticle.GetSide() == -1 )
00245 {
00246 pMass = targetParticle.GetMass()/GeV;
00247 targetParticle = *vec[0];
00248 for( G4int j=0; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];
00249 --backwardCount;
00250 backwardEnergy += pMass;
00251 backwardMass -= pMass;
00252 secondaryDeleted = true;
00253 }
00254 else if( targetParticle.GetSide() == 1 )
00255 {
00256 pMass = targetParticle.GetMass()/GeV;
00257 targetParticle = *vec[0];
00258 for( G4int j=0; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];
00259 --forwardCount;
00260 forwardEnergy += pMass;
00261 forwardMass -= pMass;
00262 secondaryDeleted = true;
00263 }
00264
00265 if( secondaryDeleted )
00266 {
00267 delete vec[vecLen-1];
00268 --vecLen;
00269 }
00270 else
00271 {
00272 if( currentParticle.GetSide() == -1 )
00273 {
00274 pMass = currentParticle.GetMass()/GeV;
00275 currentParticle = *vec[0];
00276 for( G4int j=0; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];
00277 --backwardCount;
00278 backwardEnergy += pMass;
00279 backwardMass -= pMass;
00280 secondaryDeleted = true;
00281 }
00282 else if( currentParticle.GetSide() == 1 )
00283 {
00284 pMass = currentParticle.GetMass()/GeV;
00285 currentParticle = *vec[0];
00286 for( G4int j=0; j<(vecLen-1); ++j )*vec[j] = *vec[j+1];
00287 --forwardCount;
00288 forwardEnergy += pMass;
00289 forwardMass -= pMass;
00290 secondaryDeleted = true;
00291 }
00292 if( secondaryDeleted )
00293 {
00294 delete vec[vecLen-1];
00295 --vecLen;
00296 }
00297 else break;
00298
00299 }
00300 }
00301
00302 eAvailable = centerofmassEnergy - (forwardMass+backwardMass);
00303 }
00304
00305
00306
00307
00308
00309
00310
00311
00312
00313
00314 const G4double cpar[] = { 1.60, 1.35, 1.15, 1.10 };
00315 const G4double gpar[] = { 2.60, 1.80, 1.30, 1.20 };
00316 G4int ntc = 0;
00317
00318 if (forwardCount < 1 || backwardCount < 1) return false;
00319
00320 G4double rmc = forwardMass;
00321 if (forwardCount > 1) {
00322 ntc = std::min(3,forwardCount-2);
00323 rmc += std::pow(-std::log(1.0-G4UniformRand()),1./cpar[ntc])/gpar[ntc];
00324 }
00325 G4double rmd = backwardMass;
00326 if( backwardCount > 1 ) {
00327 ntc = std::min(3,backwardCount-2);
00328 rmd += std::pow(-std::log(1.0-G4UniformRand()),1./cpar[ntc])/gpar[ntc];
00329 }
00330
00331 while( rmc+rmd > centerofmassEnergy )
00332 {
00333 if( (rmc <= forwardMass) && (rmd <= backwardMass) )
00334 {
00335 G4double temp = 0.999*centerofmassEnergy/(rmc+rmd);
00336 rmc *= temp;
00337 rmd *= temp;
00338 }
00339 else
00340 {
00341 rmc = 0.1*forwardMass + 0.9*rmc;
00342 rmd = 0.1*backwardMass + 0.9*rmd;
00343 }
00344 }
00345
00346 G4ReactionProduct pseudoParticle[8];
00347 for( i=0; i<8; ++i )pseudoParticle[i].SetZero();
00348
00349 pseudoParticle[1].SetMass( mOriginal*GeV );
00350 pseudoParticle[1].SetTotalEnergy( etOriginal*GeV );
00351 pseudoParticle[1].SetMomentum( 0.0, 0.0, pOriginal*GeV );
00352
00353 pseudoParticle[2].SetMass( protonMass*MeV );
00354 pseudoParticle[2].SetTotalEnergy( protonMass*MeV );
00355 pseudoParticle[2].SetMomentum( 0.0, 0.0, 0.0 );
00356
00357
00358
00359 pseudoParticle[0] = pseudoParticle[1] + pseudoParticle[2];
00360 pseudoParticle[1].Lorentz( pseudoParticle[1], pseudoParticle[0] );
00361 pseudoParticle[2].Lorentz( pseudoParticle[2], pseudoParticle[0] );
00362
00363
00364
00365
00366
00367 const G4double pfMin = 0.0001;
00368 G4double pf = (centerofmassEnergy*centerofmassEnergy+rmd*rmd-rmc*rmc);
00369 pf *= pf;
00370 pf -= 4*centerofmassEnergy*centerofmassEnergy*rmd*rmd;
00371 pf = std::sqrt( std::max(pf,pfMin) )/(2.0*centerofmassEnergy);
00372
00373
00374
00375 pseudoParticle[3].SetMass( rmc*GeV );
00376 pseudoParticle[3].SetTotalEnergy( std::sqrt(pf*pf+rmc*rmc)*GeV );
00377
00378 pseudoParticle[4].SetMass( rmd*GeV );
00379 pseudoParticle[4].SetTotalEnergy( std::sqrt(pf*pf+rmd*rmd)*GeV );
00380
00381
00382
00383
00384 const G4double bMin = 0.01;
00385 const G4double b1 = 4.0;
00386 const G4double b2 = 1.6;
00387 G4double pin = pseudoParticle[1].GetMomentum().mag()/GeV;
00388 G4double dtb = 4.0*pin*pf*std::max( bMin, b1+b2*std::log(pOriginal) );
00389 G4double factor = 1.0 - std::exp(-dtb);
00390 G4double costheta = 1.0 + 2.0*std::log(1.0 - G4UniformRand()*factor) / dtb;
00391
00392 costheta = std::max(-1.0, std::min(1.0, costheta));
00393 G4double sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
00394 G4double phi = G4UniformRand() * twopi;
00395
00396
00397
00398 pseudoParticle[3].SetMomentum( pf*sintheta*std::cos(phi)*GeV,
00399 pf*sintheta*std::sin(phi)*GeV,
00400 pf*costheta*GeV );
00401 pseudoParticle[4].SetMomentum( -pseudoParticle[3].GetMomentum());
00402
00403
00404
00405
00406 G4double pp, pp1;
00407 if( nuclearExcitationCount > 0 )
00408 {
00409 const G4double ga = 1.2;
00410 G4double ekit1 = 0.04;
00411 G4double ekit2 = 0.6;
00412 if( ekOriginal <= 5.0 )
00413 {
00414 ekit1 *= ekOriginal*ekOriginal/25.0;
00415 ekit2 *= ekOriginal*ekOriginal/25.0;
00416 }
00417 G4double scale = std::pow(ekit2/ekit1, 1.0-ga) - 1.0;
00418 for( i=0; i<vecLen; ++i )
00419 {
00420 if( vec[i]->GetSide() == -2 )
00421 {
00422 G4double kineticE = ekit1*std::pow((1.0 + G4UniformRand()*scale), 1.0/(1.0-ga) );
00423 vec[i]->SetKineticEnergy( kineticE*GeV );
00424 G4double vMass = vec[i]->GetMass()/MeV;
00425 G4double totalE = kineticE*GeV + vMass;
00426 pp = std::sqrt( std::abs(totalE*totalE-vMass*vMass) );
00427 G4double cost = std::min( 1.0, std::max( -1.0, std::log(2.23*G4UniformRand()+0.383)/0.96 ) );
00428 G4double sint = std::sqrt(1.0-cost*cost);
00429 phi = twopi*G4UniformRand();
00430 vec[i]->SetMomentum( pp*sint*std::cos(phi)*MeV,
00431 pp*sint*std::sin(phi)*MeV,
00432 pp*cost*MeV );
00433 vec[i]->Lorentz( *vec[i], pseudoParticle[0] );
00434 }
00435 }
00436 }
00437
00438
00439
00440
00441
00442 currentParticle.SetMomentum( pseudoParticle[3].GetMomentum() );
00443 currentParticle.SetTotalEnergy( pseudoParticle[3].GetTotalEnergy() );
00444
00445 targetParticle.SetMomentum( pseudoParticle[4].GetMomentum() );
00446 targetParticle.SetTotalEnergy( pseudoParticle[4].GetTotalEnergy() );
00447
00448 pseudoParticle[5].SetMomentum( pseudoParticle[3].GetMomentum() * (-1.0) );
00449 pseudoParticle[5].SetMass( pseudoParticle[3].GetMass() );
00450 pseudoParticle[5].SetTotalEnergy( pseudoParticle[3].GetTotalEnergy() );
00451
00452 pseudoParticle[6].SetMomentum( pseudoParticle[4].GetMomentum() * (-1.0) );
00453 pseudoParticle[6].SetMass( pseudoParticle[4].GetMass() );
00454 pseudoParticle[6].SetTotalEnergy( pseudoParticle[4].GetTotalEnergy() );
00455
00456 G4double wgt;
00457
00458 if( forwardCount > 1 )
00459 {
00460 G4FastVector<G4ReactionProduct,256> tempV;
00461 tempV.Initialize( forwardCount );
00462 G4bool constantCrossSection = true;
00463 G4int tempLen = 0;
00464 if( currentParticle.GetSide() == 1 )
00465 tempV.SetElement( tempLen++, ¤tParticle );
00466 if( targetParticle.GetSide() == 1 )
00467 tempV.SetElement( tempLen++, &targetParticle );
00468 for( i=0; i<vecLen; ++i )
00469 {
00470 if( vec[i]->GetSide() == 1 )
00471 {
00472 if( tempLen < 18 )
00473 tempV.SetElement( tempLen++, vec[i] );
00474 else
00475 {
00476 vec[i]->SetSide( -1 );
00477 continue;
00478 }
00479 }
00480 }
00481 if( tempLen >= 2 )
00482 {
00483 wgt = GenerateNBodyEvent( pseudoParticle[3].GetMass()/MeV,
00484 constantCrossSection, tempV, tempLen );
00485 if( currentParticle.GetSide() == 1 )
00486 currentParticle.Lorentz( currentParticle, pseudoParticle[5] );
00487 if( targetParticle.GetSide() == 1 )
00488 targetParticle.Lorentz( targetParticle, pseudoParticle[5] );
00489 for( i=0; i<vecLen; ++i )
00490 {
00491 if( vec[i]->GetSide() == 1 )vec[i]->Lorentz( *vec[i], pseudoParticle[5] );
00492 }
00493 }
00494 }
00495
00496 if( backwardCount > 1 )
00497 {
00498 G4FastVector<G4ReactionProduct,256> tempV;
00499 tempV.Initialize( backwardCount );
00500 G4bool constantCrossSection = true;
00501 G4int tempLen = 0;
00502 if( currentParticle.GetSide() == -1 )
00503 tempV.SetElement( tempLen++, ¤tParticle );
00504 if( targetParticle.GetSide() == -1 )
00505 tempV.SetElement( tempLen++, &targetParticle );
00506 for( i=0; i<vecLen; ++i )
00507 {
00508 if( vec[i]->GetSide() == -1 )
00509 {
00510 if( tempLen < 18 )
00511 tempV.SetElement( tempLen++, vec[i] );
00512 else
00513 {
00514 vec[i]->SetSide( -2 );
00515 vec[i]->SetKineticEnergy( 0.0 );
00516 vec[i]->SetMomentum( 0.0, 0.0, 0.0 );
00517 continue;
00518 }
00519 }
00520 }
00521 if( tempLen >= 2 )
00522 {
00523 wgt = GenerateNBodyEvent( pseudoParticle[4].GetMass()/MeV,
00524 constantCrossSection, tempV, tempLen );
00525 if( currentParticle.GetSide() == -1 )
00526 currentParticle.Lorentz( currentParticle, pseudoParticle[6] );
00527 if( targetParticle.GetSide() == -1 )
00528 targetParticle.Lorentz( targetParticle, pseudoParticle[6] );
00529 for( i=0; i<vecLen; ++i )
00530 {
00531 if( vec[i]->GetSide() == -1 )vec[i]->Lorentz( *vec[i], pseudoParticle[6] );
00532 }
00533 }
00534 }
00535
00536
00537
00538
00539
00540 currentParticle.Lorentz( currentParticle, pseudoParticle[2] );
00541 targetParticle.Lorentz( targetParticle, pseudoParticle[2] );
00542 for( i=0; i<vecLen; ++i ) vec[i]->Lorentz( *vec[i], pseudoParticle[2] );
00543
00544
00545
00546
00547
00548 G4bool dum = true;
00549 if( leadFlag )
00550 {
00551
00552
00553
00554
00555
00556
00557 if( currentParticle.GetDefinition() == leadingStrangeParticle.GetDefinition() )
00558 dum = false;
00559 else if( targetParticle.GetDefinition() == leadingStrangeParticle.GetDefinition() )
00560 dum = false;
00561 else
00562 {
00563 for( i=0; i<vecLen; ++i )
00564 {
00565 if( vec[i]->GetDefinition() == leadingStrangeParticle.GetDefinition() )
00566 {
00567 dum = false;
00568 break;
00569 }
00570 }
00571 }
00572 if( dum )
00573 {
00574 G4double leadMass = leadingStrangeParticle.GetMass()/MeV;
00575 G4double ekin;
00576 if( ((leadMass < protonMass) && (targetParticle.GetMass()/MeV < protonMass)) ||
00577 ((leadMass >= protonMass) && (targetParticle.GetMass()/MeV >= protonMass)) )
00578 {
00579 ekin = targetParticle.GetKineticEnergy()/GeV;
00580 pp1 = targetParticle.GetMomentum().mag()/MeV;
00581 targetParticle.SetDefinition( leadingStrangeParticle.GetDefinition() );
00582 targetParticle.SetKineticEnergy( ekin*GeV );
00583 pp = targetParticle.GetTotalMomentum()/MeV;
00584 if( pp1 < 1.0e-3 ) {
00585 G4ThreeVector iso = Isotropic(pp);
00586 targetParticle.SetMomentum( iso.x(), iso.y(), iso.z() );
00587 } else {
00588 targetParticle.SetMomentum( targetParticle.GetMomentum() * (pp/pp1) );
00589 }
00590 targetHasChanged = true;
00591 }
00592 else
00593 {
00594 ekin = currentParticle.GetKineticEnergy()/GeV;
00595 pp1 = currentParticle.GetMomentum().mag()/MeV;
00596 currentParticle.SetDefinition( leadingStrangeParticle.GetDefinition() );
00597 currentParticle.SetKineticEnergy( ekin*GeV );
00598 pp = currentParticle.GetTotalMomentum()/MeV;
00599 if( pp1 < 1.0e-3 ) {
00600 G4ThreeVector iso = Isotropic(pp);
00601 currentParticle.SetMomentum( iso.x(), iso.y(), iso.z() );
00602 } else {
00603 currentParticle.SetMomentum( currentParticle.GetMomentum() * (pp/pp1) );
00604 }
00605 incidentHasChanged = true;
00606 }
00607 }
00608 }
00609
00610
00611
00612
00613 std::pair<G4int, G4int> finalStateNucleons =
00614 GetFinalStateNucleons(originalTarget, vec, vecLen);
00615
00616 G4int protonsInFinalState = finalStateNucleons.first;
00617 G4int neutronsInFinalState = finalStateNucleons.second;
00618
00619 G4int numberofFinalStateNucleons =
00620 protonsInFinalState + neutronsInFinalState;
00621
00622 if (currentParticle.GetDefinition()->GetBaryonNumber() == 1 &&
00623 targetParticle.GetDefinition()->GetBaryonNumber() == 1 &&
00624 originalIncident->GetDefinition()->GetPDGMass() <
00625 G4Lambda::Lambda()->GetPDGMass())
00626 numberofFinalStateNucleons++;
00627
00628 numberofFinalStateNucleons = std::max(1, numberofFinalStateNucleons);
00629
00630 G4int PinNucleus = std::max(0,
00631 G4int(targetNucleus.GetZ_asInt()) - protonsInFinalState);
00632 G4int NinNucleus = std::max(0,
00633 G4int(targetNucleus.GetA_asInt()-targetNucleus.GetZ_asInt()) - neutronsInFinalState);
00634
00635
00636
00637
00638 pseudoParticle[4].SetMass( mOriginal*GeV );
00639 pseudoParticle[4].SetTotalEnergy( etOriginal*GeV );
00640 pseudoParticle[4].SetMomentum( 0.0, 0.0, pOriginal*GeV );
00641
00642 G4ParticleDefinition * aOrgDef = modifiedOriginal.GetDefinition();
00643 G4int diff = 0;
00644 if(aOrgDef == G4Proton::Proton() || aOrgDef == G4Neutron::Neutron() ) diff = 1;
00645 if(numberofFinalStateNucleons == 1) diff = 0;
00646 pseudoParticle[5].SetMomentum( 0.0, 0.0, 0.0 );
00647 pseudoParticle[5].SetMass( protonMass*(numberofFinalStateNucleons-diff)*MeV);
00648 pseudoParticle[5].SetTotalEnergy( protonMass*(numberofFinalStateNucleons-diff)*MeV);
00649
00650 G4double theoreticalKinetic =
00651 pseudoParticle[4].GetTotalEnergy()/GeV + pseudoParticle[5].GetTotalEnergy()/GeV;
00652
00653 pseudoParticle[6] = pseudoParticle[4] + pseudoParticle[5];
00654 pseudoParticle[4].Lorentz( pseudoParticle[4], pseudoParticle[6] );
00655 pseudoParticle[5].Lorentz( pseudoParticle[5], pseudoParticle[6] );
00656
00657 if( vecLen < 16 )
00658 {
00659 G4ReactionProduct tempR[130];
00660 tempR[0] = currentParticle;
00661 tempR[1] = targetParticle;
00662 for( i=0; i<vecLen; ++i )tempR[i+2] = *vec[i];
00663
00664 G4FastVector<G4ReactionProduct,256> tempV;
00665 tempV.Initialize( vecLen+2 );
00666 G4bool constantCrossSection = true;
00667 G4int tempLen = 0;
00668 for( i=0; i<vecLen+2; ++i )tempV.SetElement( tempLen++, &tempR[i] );
00669
00670 if( tempLen >= 2 )
00671 {
00672
00673 wgt = GenerateNBodyEvent( pseudoParticle[4].GetTotalEnergy()/MeV +
00674 pseudoParticle[5].GetTotalEnergy()/MeV,
00675 constantCrossSection, tempV, tempLen );
00676 if (wgt == -1) {
00677 G4double Qvalue = 0;
00678 for (i = 0; i < tempLen; i++) Qvalue += tempV[i]->GetMass();
00679 wgt = GenerateNBodyEvent( Qvalue/MeV,
00680 constantCrossSection, tempV, tempLen );
00681 }
00682 theoreticalKinetic = 0.0;
00683 for( i=0; i<vecLen+2; ++i )
00684 {
00685 pseudoParticle[7].SetMomentum( tempV[i]->GetMomentum() );
00686 pseudoParticle[7].SetMass( tempV[i]->GetMass() );
00687 pseudoParticle[7].SetTotalEnergy( tempV[i]->GetTotalEnergy() );
00688 pseudoParticle[7].Lorentz( pseudoParticle[7], pseudoParticle[5] );
00689 theoreticalKinetic += pseudoParticle[7].GetKineticEnergy()/GeV;
00690 }
00691 }
00692
00693 }
00694 else
00695 {
00696 theoreticalKinetic -=
00697 ( currentParticle.GetMass()/GeV + targetParticle.GetMass()/GeV );
00698 for( i=0; i<vecLen; ++i )theoreticalKinetic -= vec[i]->GetMass()/GeV;
00699 }
00700 G4double simulatedKinetic =
00701 currentParticle.GetKineticEnergy()/GeV + targetParticle.GetKineticEnergy()/GeV;
00702 for( i=0; i<vecLen; ++i )simulatedKinetic += vec[i]->GetKineticEnergy()/GeV;
00703
00704
00705
00706
00707 if( simulatedKinetic != 0.0 )
00708 {
00709 wgt = (theoreticalKinetic)/simulatedKinetic;
00710 currentParticle.SetKineticEnergy( wgt*currentParticle.GetKineticEnergy() );
00711 pp = currentParticle.GetTotalMomentum()/MeV;
00712 pp1 = currentParticle.GetMomentum().mag()/MeV;
00713 if( pp1 < 0.001*MeV ) {
00714 G4ThreeVector iso = Isotropic(pp);
00715 currentParticle.SetMomentum( iso.x(), iso.y(), iso.z() );
00716 } else {
00717 currentParticle.SetMomentum( currentParticle.GetMomentum() * (pp/pp1) );
00718 }
00719
00720 targetParticle.SetKineticEnergy( wgt*targetParticle.GetKineticEnergy() );
00721 pp = targetParticle.GetTotalMomentum()/MeV;
00722 pp1 = targetParticle.GetMomentum().mag()/MeV;
00723 if( pp1 < 0.001*MeV ) {
00724 G4ThreeVector iso = Isotropic(pp);
00725 targetParticle.SetMomentum( iso.x(), iso.y(), iso.z() );
00726 } else {
00727 targetParticle.SetMomentum( targetParticle.GetMomentum() * (pp/pp1) );
00728 }
00729
00730 for( i=0; i<vecLen; ++i )
00731 {
00732 vec[i]->SetKineticEnergy( wgt*vec[i]->GetKineticEnergy() );
00733 pp = vec[i]->GetTotalMomentum()/MeV;
00734 pp1 = vec[i]->GetMomentum().mag()/MeV;
00735 if( pp1 < 0.001 ) {
00736 G4ThreeVector iso = Isotropic(pp);
00737 vec[i]->SetMomentum( iso.x(), iso.y(), iso.z() );
00738 } else {
00739 vec[i]->SetMomentum( vec[i]->GetMomentum() * (pp/pp1) );
00740 }
00741 }
00742 }
00743
00744
00745 Rotate( numberofFinalStateNucleons, pseudoParticle[4].GetMomentum(),
00746 modifiedOriginal, originalIncident, targetNucleus,
00747 currentParticle, targetParticle, vec, vecLen );
00748
00749
00750
00751
00752
00753 if( atomicWeight >= 1.5 )
00754 {
00755
00756
00757
00758
00759
00760
00761
00762 G4int npnb = 0;
00763 G4int ndta = 0;
00764
00765 G4double epnb, edta;
00766 if (veryForward) {
00767 epnb = targetNucleus.GetAnnihilationPNBlackTrackEnergy();
00768 edta = targetNucleus.GetAnnihilationDTABlackTrackEnergy();
00769 } else {
00770 epnb = targetNucleus.GetPNBlackTrackEnergy();
00771 edta = targetNucleus.GetDTABlackTrackEnergy();
00772 }
00773
00774 const G4double pnCutOff = 0.001;
00775 const G4double dtaCutOff = 0.001;
00776
00777
00778
00779
00780
00781
00782
00783
00784 if( epnb >= pnCutOff )
00785 {
00786 npnb = G4Poisson((1.5+1.25*numberofFinalStateNucleons)*epnb/(epnb+edta));
00787 if( numberofFinalStateNucleons + npnb > atomicWeight )
00788 npnb = G4int(atomicWeight - numberofFinalStateNucleons);
00789 npnb = std::min( npnb, 127-vecLen );
00790 }
00791 if( edta >= dtaCutOff )
00792 {
00793 ndta = G4Poisson( (1.5+1.25*numberofFinalStateNucleons)*edta/(epnb+edta) );
00794 ndta = std::min( ndta, 127-vecLen );
00795 }
00796 if (npnb == 0 && ndta == 0) npnb = 1;
00797
00798
00799
00800 AddBlackTrackParticles(epnb, npnb, edta, ndta, modifiedOriginal,
00801 PinNucleus, NinNucleus, targetNucleus,
00802 vec, vecLen );
00803
00804 }
00805
00806
00807
00808
00809
00810
00811 if( (atomicWeight >= 1.5) && (atomicWeight <= 230.0) && (ekOriginal <= 0.2) )
00812 currentParticle.SetTOF( 1.0-500.0*std::exp(-ekOriginal/0.04)*std::log(G4UniformRand()) );
00813 else
00814 currentParticle.SetTOF( 1.0 );
00815
00816 return true;
00817 }
00818
00819