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
00031 #include "G4RPGReaction.hh"
00032 #include "G4PhysicalConstants.hh"
00033 #include "G4SystemOfUnits.hh"
00034 #include "Randomize.hh"
00035
00036 G4bool G4RPGReaction::
00037 ReactionStage(const G4HadProjectile* ,
00038 G4ReactionProduct& ,
00039 G4bool& ,
00040 const G4DynamicParticle* ,
00041 G4ReactionProduct& ,
00042 G4bool& ,
00043 const G4Nucleus& ,
00044 G4ReactionProduct& ,
00045 G4FastVector<G4ReactionProduct,256>& ,
00046 G4int& ,
00047 G4bool ,
00048 G4ReactionProduct& )
00049 {
00050 G4cout << " G4RPGReactionStage must be overridden in a derived class "
00051 << G4endl;
00052 return false;
00053 }
00054
00055
00056 void G4RPGReaction::
00057 AddBlackTrackParticles(const G4double epnb,
00058 const G4int npnb,
00059 const G4double edta,
00060 const G4int ndta,
00061 const G4ReactionProduct& modifiedOriginal,
00062 G4int PinNucleus,
00063 G4int NinNucleus,
00064 const G4Nucleus& targetNucleus,
00065 G4FastVector<G4ReactionProduct,256>& vec,
00066 G4int& vecLen)
00067 {
00068
00069
00070
00071
00072
00073
00074
00075 G4ParticleDefinition* aProton = G4Proton::Proton();
00076 G4ParticleDefinition* aNeutron = G4Neutron::Neutron();
00077 G4ParticleDefinition* aDeuteron = G4Deuteron::Deuteron();
00078 G4ParticleDefinition* aTriton = G4Triton::Triton();
00079 G4ParticleDefinition* anAlpha = G4Alpha::Alpha();
00080
00081 const G4double ekOriginal = modifiedOriginal.GetKineticEnergy()/MeV;
00082 const G4double atomicWeight = targetNucleus.GetA_asInt();
00083 const G4double atomicNumber = targetNucleus.GetZ_asInt();
00084
00085 const G4double ika1 = 3.6;
00086 const G4double ika2 = 35.56;
00087 const G4double ika3 = 6.45;
00088
00089 G4int i;
00090 G4double pp;
00091 G4double kinetic = 0;
00092 G4double kinCreated = 0;
00093
00094 G4double remainingE = 0;
00095
00096
00097 if (npnb > 0) {
00098
00099 G4int local_npnb = npnb;
00100
00101 local_npnb = std::min(PinNucleus + NinNucleus , local_npnb);
00102 G4double local_epnb = epnb;
00103 if (ndta == 0) local_epnb += edta;
00104
00105
00106 remainingE = local_epnb;
00107 for (i = 0; i < local_npnb; ++i)
00108 {
00109 G4ReactionProduct* p1 = new G4ReactionProduct();
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122 if (G4UniformRand() > (1.0-atomicNumber/atomicWeight)) {
00123
00124
00125
00126 if (PinNucleus > 0) {
00127 p1->SetDefinition( aProton );
00128 PinNucleus--;
00129 } else {
00130 p1->SetDefinition( aNeutron );
00131 NinNucleus--;
00132
00133
00134
00135 }
00136 } else {
00137
00138
00139
00140 if (NinNucleus > 0) {
00141 p1->SetDefinition( aNeutron );
00142 NinNucleus--;
00143 } else {
00144 p1->SetDefinition( aProton );
00145 PinNucleus--;
00146
00147
00148
00149 }
00150 }
00151
00152 if (i < local_npnb - 1) {
00153 kinetic = remainingE*G4UniformRand();
00154 remainingE -= kinetic;
00155 } else {
00156 kinetic = remainingE;
00157 }
00158
00159 vec.SetElement( vecLen, p1 );
00160 G4double cost = G4UniformRand() * 2.0 - 1.0;
00161 G4double sint = std::sqrt(std::fabs(1.0-cost*cost));
00162 G4double phi = twopi * G4UniformRand();
00163 vec[vecLen]->SetNewlyAdded( true );
00164 vec[vecLen]->SetKineticEnergy( kinetic*GeV );
00165 kinCreated+=kinetic;
00166 pp = vec[vecLen]->GetTotalMomentum();
00167 vec[vecLen]->SetMomentum(pp*sint*std::sin(phi),
00168 pp*sint*std::cos(phi),
00169 pp*cost );
00170 vecLen++;
00171 }
00172
00173 if (NinNucleus > 0) {
00174 if( (atomicWeight >= 10.0) && (ekOriginal <= 2.0*GeV) )
00175 {
00176 G4double ekw = ekOriginal/GeV;
00177 G4int ika, kk = 0;
00178 if( ekw > 1.0 )ekw *= ekw;
00179 ekw = std::max( 0.1, ekw );
00180 ika = G4int(ika1*std::exp((atomicNumber*atomicNumber/
00181 atomicWeight-ika2)/ika3)/ekw);
00182 if( ika > 0 )
00183 {
00184 for( i=(vecLen-1); i>=0; --i )
00185 {
00186 if( (vec[i]->GetDefinition() == aProton) && vec[i]->GetNewlyAdded() )
00187 {
00188 vec[i]->SetDefinitionAndUpdateE( aNeutron );
00189 PinNucleus++;
00190 NinNucleus--;
00191 if( ++kk > ika )break;
00192 }
00193 }
00194 }
00195 }
00196 }
00197 }
00198
00199
00200
00201 G4double ran = 0;
00202 if (ndta > 0) {
00203
00204 G4int local_ndta=ndta;
00205
00206 G4double local_edta = edta;
00207 if (npnb == 0) local_edta += epnb;
00208
00209
00210 remainingE = local_edta;
00211 for (i = 0; i < local_ndta; ++i) {
00212 G4ReactionProduct* p2 = new G4ReactionProduct();
00213
00214
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224
00225 ran = G4UniformRand();
00226 if (ran < 0.60) {
00227 if (PinNucleus > 0 && NinNucleus > 0) {
00228 p2->SetDefinition( aDeuteron );
00229 PinNucleus--;
00230 NinNucleus--;
00231 } else if (NinNucleus > 0) {
00232 p2->SetDefinition( aNeutron );
00233 NinNucleus--;
00234 } else if (PinNucleus > 0) {
00235 p2->SetDefinition( aProton );
00236 PinNucleus--;
00237 } else {
00238 delete p2;
00239 break;
00240 }
00241 } else if (ran < 0.90) {
00242 if (PinNucleus > 0 && NinNucleus > 1) {
00243 p2->SetDefinition( aTriton );
00244 PinNucleus--;
00245 NinNucleus -= 2;
00246 } else if (PinNucleus > 0 && NinNucleus > 0) {
00247 p2->SetDefinition( aDeuteron );
00248 PinNucleus--;
00249 NinNucleus--;
00250 } else if (NinNucleus > 0) {
00251 p2->SetDefinition( aNeutron );
00252 NinNucleus--;
00253 } else if (PinNucleus > 0) {
00254 p2->SetDefinition( aProton );
00255 PinNucleus--;
00256 } else {
00257 delete p2;
00258 break;
00259 }
00260 } else {
00261 if (PinNucleus > 1 && NinNucleus > 1) {
00262 p2->SetDefinition( anAlpha );
00263 PinNucleus -= 2;
00264 NinNucleus -= 2;
00265 } else if (PinNucleus > 0 && NinNucleus > 1) {
00266 p2->SetDefinition( aTriton );
00267 PinNucleus--;
00268 NinNucleus -= 2;
00269 } else if (PinNucleus > 0 && NinNucleus > 0) {
00270 p2->SetDefinition( aDeuteron );
00271 PinNucleus--;
00272 NinNucleus--;
00273 } else if (NinNucleus > 0) {
00274 p2->SetDefinition( aNeutron );
00275 NinNucleus--;
00276 } else if (PinNucleus > 0) {
00277 p2->SetDefinition( aProton );
00278 PinNucleus--;
00279 } else {
00280 delete p2;
00281 break;
00282 }
00283 }
00284
00285 if (i < local_ndta - 1) {
00286 kinetic = remainingE*G4UniformRand();
00287 remainingE -= kinetic;
00288 } else {
00289 kinetic = remainingE;
00290 }
00291
00292 vec.SetElement( vecLen, p2 );
00293 G4double cost = 2.0*G4UniformRand() - 1.0;
00294 G4double sint = std::sqrt(std::max(0.0,(1.0-cost*cost)));
00295 G4double phi = twopi*G4UniformRand();
00296 vec[vecLen]->SetNewlyAdded( true );
00297 vec[vecLen]->SetKineticEnergy( kinetic*GeV );
00298 kinCreated+=kinetic;
00299
00300 pp = vec[vecLen]->GetTotalMomentum();
00301 vec[vecLen]->SetMomentum( pp*sint*std::sin(phi),
00302 pp*sint*std::cos(phi),
00303 pp*cost );
00304 vecLen++;
00305 }
00306 }
00307 }
00308
00309
00310 G4double
00311 G4RPGReaction::GenerateNBodyEvent(const G4double totalEnergy,
00312 const G4bool constantCrossSection,
00313 G4FastVector<G4ReactionProduct,256>& vec,
00314 G4int &vecLen)
00315 {
00316 G4int i;
00317 const G4double expxu = 82.;
00318 const G4double expxl = -expxu;
00319
00320 if (vecLen < 2) {
00321 G4cerr << "*** Error in G4RPGReaction::GenerateNBodyEvent" << G4endl;
00322 G4cerr << " number of particles < 2" << G4endl;
00323 G4cerr << "totalEnergy = " << totalEnergy << "MeV, vecLen = " << vecLen << G4endl;
00324 return -1.0;
00325 }
00326
00327 G4double mass[18];
00328 G4double energy[18];
00329 G4double pcm[3][18];
00330
00331 G4double totalMass = 0.0;
00332 G4double extraMass = 0;
00333 G4double sm[18];
00334
00335 for (i=0; i<vecLen; ++i) {
00336 mass[i] = vec[i]->GetMass()/GeV;
00337 if(vec[i]->GetSide() == -2) extraMass+=vec[i]->GetMass()/GeV;
00338 vec[i]->SetMomentum( 0.0, 0.0, 0.0 );
00339 pcm[0][i] = 0.0;
00340 pcm[1][i] = 0.0;
00341 pcm[2][i] = 0.0;
00342 energy[i] = mass[i];
00343 totalMass += mass[i];
00344 sm[i] = totalMass;
00345 }
00346
00347 G4double totalE = totalEnergy/GeV;
00348 if (totalMass > totalE) {
00349
00350
00351
00352 totalE = totalMass;
00353 return -1.0;
00354 }
00355
00356 G4double kineticEnergy = totalE - totalMass;
00357 G4double emm[18];
00358 emm[0] = mass[0];
00359 emm[vecLen-1] = totalE;
00360
00361 if (vecLen > 2) {
00362 G4double ran[18];
00363 for( i=0; i<vecLen; ++i )ran[i] = G4UniformRand();
00364 for (i=0; i<vecLen-2; ++i) {
00365 for (G4int j=vecLen-2; j>i; --j) {
00366 if (ran[i] > ran[j]) {
00367 G4double temp = ran[i];
00368 ran[i] = ran[j];
00369 ran[j] = temp;
00370 }
00371 }
00372 }
00373 for( i=1; i<vecLen-1; ++i )emm[i] = ran[i-1]*kineticEnergy + sm[i];
00374 }
00375
00376
00377
00378 G4bool lzero = true;
00379 G4double wtmax = 0.0;
00380 if (constantCrossSection) {
00381 G4double emmax = kineticEnergy + mass[0];
00382 G4double emmin = 0.0;
00383 for( i=1; i<vecLen; ++i )
00384 {
00385 emmin += mass[i-1];
00386 emmax += mass[i];
00387 G4double wtfc = 0.0;
00388 if( emmax*emmax > 0.0 )
00389 {
00390 G4double arg = emmax*emmax
00391 + (emmin*emmin-mass[i]*mass[i])*(emmin*emmin-mass[i]*mass[i])/(emmax*emmax)
00392 - 2.0*(emmin*emmin+mass[i]*mass[i]);
00393 if( arg > 0.0 )wtfc = 0.5*std::sqrt( arg );
00394 }
00395 if( wtfc == 0.0 )
00396 {
00397 lzero = false;
00398 break;
00399 }
00400 wtmax += std::log( wtfc );
00401 }
00402 if( lzero )
00403 wtmax = -wtmax;
00404 else
00405 wtmax = expxu;
00406 } else {
00407
00408 const G4double ffq[18] = { 0., 3.141592, 19.73921, 62.01255, 129.8788, 204.0131,
00409 256.3704, 268.4705, 240.9780, 189.2637,
00410 132.1308, 83.0202, 47.4210, 24.8295,
00411 12.0006, 5.3858, 2.2560, 0.8859 };
00412 wtmax = std::log( std::pow( kineticEnergy, vecLen-2 ) * ffq[vecLen-1] / totalE );
00413 }
00414
00415
00416
00417 lzero = true;
00418 G4double pd[50];
00419
00420 for( i=0; i<vecLen-1; ++i )
00421 {
00422 pd[i] = 0.0;
00423 if( emm[i+1]*emm[i+1] > 0.0 )
00424 {
00425 G4double arg = emm[i+1]*emm[i+1]
00426 + (emm[i]*emm[i]-mass[i+1]*mass[i+1])*(emm[i]*emm[i]-mass[i+1]*mass[i+1])
00427 /(emm[i+1]*emm[i+1])
00428 - 2.0*(emm[i]*emm[i]+mass[i+1]*mass[i+1]);
00429 if( arg > 0.0 )pd[i] = 0.5*std::sqrt( arg );
00430 }
00431 if( pd[i] <= 0.0 )
00432 lzero = false;
00433 else
00434 wtmax += std::log( pd[i] );
00435 }
00436 G4double weight = 0.0;
00437 if( lzero )weight = std::exp( std::max(std::min(wtmax,expxu),expxl) );
00438
00439 G4double bang, cb, sb, s0, s1, s2, c, esys, a, b, gama, beta;
00440 G4double ss;
00441 pcm[0][0] = 0.0;
00442 pcm[1][0] = pd[0];
00443 pcm[2][0] = 0.0;
00444 for( i=1; i<vecLen; ++i )
00445 {
00446 pcm[0][i] = 0.0;
00447 pcm[1][i] = -pd[i-1];
00448 pcm[2][i] = 0.0;
00449 bang = twopi*G4UniformRand();
00450 cb = std::cos(bang);
00451 sb = std::sin(bang);
00452 c = 2.0*G4UniformRand() - 1.0;
00453 ss = std::sqrt( std::fabs( 1.0-c*c ) );
00454 if( i < vecLen-1 )
00455 {
00456 esys = std::sqrt(pd[i]*pd[i] + emm[i]*emm[i]);
00457 beta = pd[i]/esys;
00458 gama = esys/emm[i];
00459 for( G4int j=0; j<=i; ++j )
00460 {
00461 s0 = pcm[0][j];
00462 s1 = pcm[1][j];
00463 s2 = pcm[2][j];
00464 energy[j] = std::sqrt( s0*s0 + s1*s1 + s2*s2 + mass[j]*mass[j] );
00465 a = s0*c - s1*ss;
00466 pcm[1][j] = s0*ss + s1*c;
00467 b = pcm[2][j];
00468 pcm[0][j] = a*cb - b*sb;
00469 pcm[2][j] = a*sb + b*cb;
00470 pcm[1][j] = gama*(pcm[1][j] + beta*energy[j]);
00471 }
00472 }
00473 else
00474 {
00475 for( G4int j=0; j<=i; ++j )
00476 {
00477 s0 = pcm[0][j];
00478 s1 = pcm[1][j];
00479 s2 = pcm[2][j];
00480 energy[j] = std::sqrt( s0*s0 + s1*s1 + s2*s2 + mass[j]*mass[j] );
00481 a = s0*c - s1*s;
00482 pcm[1][j] = s0*ss + s1*c;
00483 b = pcm[2][j];
00484 pcm[0][j] = a*cb - b*sb;
00485 pcm[2][j] = a*sb + b*cb;
00486 }
00487 }
00488 }
00489
00490 for (i=0; i<vecLen; ++i) {
00491 vec[i]->SetMomentum( pcm[0][i]*GeV, pcm[1][i]*GeV, pcm[2][i]*GeV );
00492 vec[i]->SetTotalEnergy( energy[i]*GeV );
00493 }
00494
00495 return weight;
00496 }
00497
00498
00499 G4double
00500 G4RPGReaction::GenerateNBodyEventT(const G4double totalEnergy,
00501 const G4bool constantCrossSection,
00502 std::vector<G4ReactionProduct*>& tempList)
00503 {
00504 G4int i;
00505 const G4double expxu = 82.;
00506 const G4double expxl = -expxu;
00507 G4int listLen = tempList.size();
00508
00509 if (listLen < 2) {
00510 G4cerr << "*** Error in G4RPGReaction::GenerateNBodyEvent" << G4endl;
00511 G4cerr << " number of particles < 2" << G4endl;
00512 G4cerr << "totalEnergy = " << totalEnergy << "MeV, listLen = " << listLen << G4endl;
00513 return -1.0;
00514 }
00515
00516 G4double mass[18];
00517 G4double energy[18];
00518 G4double pcm[3][18];
00519
00520 G4double totalMass = 0.0;
00521 G4double extraMass = 0;
00522 G4double sm[18];
00523
00524 for (i=0; i<listLen; ++i) {
00525 mass[i] = tempList[i]->GetMass()/GeV;
00526 if(tempList[i]->GetSide() == -2) extraMass+=tempList[i]->GetMass()/GeV;
00527 tempList[i]->SetMomentum( 0.0, 0.0, 0.0 );
00528 pcm[0][i] = 0.0;
00529 pcm[1][i] = 0.0;
00530 pcm[2][i] = 0.0;
00531 energy[i] = mass[i];
00532 totalMass += mass[i];
00533 sm[i] = totalMass;
00534 }
00535
00536 G4double totalE = totalEnergy/GeV;
00537 if (totalMass > totalE) {
00538 totalE = totalMass;
00539 return -1.0;
00540 }
00541
00542 G4double kineticEnergy = totalE - totalMass;
00543 G4double emm[18];
00544 emm[0] = mass[0];
00545 emm[listLen-1] = totalE;
00546
00547 if (listLen > 2) {
00548 G4double ran[18];
00549 for( i=0; i<listLen; ++i )ran[i] = G4UniformRand();
00550 for (i=0; i<listLen-2; ++i) {
00551 for (G4int j=listLen-2; j>i; --j) {
00552 if (ran[i] > ran[j]) {
00553 G4double temp = ran[i];
00554 ran[i] = ran[j];
00555 ran[j] = temp;
00556 }
00557 }
00558 }
00559 for( i=1; i<listLen-1; ++i )emm[i] = ran[i-1]*kineticEnergy + sm[i];
00560 }
00561
00562
00563
00564 G4bool lzero = true;
00565 G4double wtmax = 0.0;
00566 if (constantCrossSection) {
00567 G4double emmax = kineticEnergy + mass[0];
00568 G4double emmin = 0.0;
00569 for( i=1; i<listLen; ++i )
00570 {
00571 emmin += mass[i-1];
00572 emmax += mass[i];
00573 G4double wtfc = 0.0;
00574 if( emmax*emmax > 0.0 )
00575 {
00576 G4double arg = emmax*emmax
00577 + (emmin*emmin-mass[i]*mass[i])*(emmin*emmin-mass[i]*mass[i])/(emmax*emmax)
00578 - 2.0*(emmin*emmin+mass[i]*mass[i]);
00579 if( arg > 0.0 )wtfc = 0.5*std::sqrt( arg );
00580 }
00581 if( wtfc == 0.0 )
00582 {
00583 lzero = false;
00584 break;
00585 }
00586 wtmax += std::log( wtfc );
00587 }
00588 if( lzero )
00589 wtmax = -wtmax;
00590 else
00591 wtmax = expxu;
00592 } else {
00593
00594 const G4double ffq[18] = { 0., 3.141592, 19.73921, 62.01255, 129.8788, 204.0131,
00595 256.3704, 268.4705, 240.9780, 189.2637,
00596 132.1308, 83.0202, 47.4210, 24.8295,
00597 12.0006, 5.3858, 2.2560, 0.8859 };
00598 wtmax = std::log( std::pow( kineticEnergy, listLen-2 ) * ffq[listLen-1] / totalE );
00599 }
00600
00601
00602
00603 lzero = true;
00604 G4double pd[50];
00605
00606 for( i=0; i<listLen-1; ++i )
00607 {
00608 pd[i] = 0.0;
00609 if( emm[i+1]*emm[i+1] > 0.0 )
00610 {
00611 G4double arg = emm[i+1]*emm[i+1]
00612 + (emm[i]*emm[i]-mass[i+1]*mass[i+1])*(emm[i]*emm[i]-mass[i+1]*mass[i+1])
00613 /(emm[i+1]*emm[i+1])
00614 - 2.0*(emm[i]*emm[i]+mass[i+1]*mass[i+1]);
00615 if( arg > 0.0 )pd[i] = 0.5*std::sqrt( arg );
00616 }
00617 if( pd[i] <= 0.0 )
00618 lzero = false;
00619 else
00620 wtmax += std::log( pd[i] );
00621 }
00622 G4double weight = 0.0;
00623 if( lzero )weight = std::exp( std::max(std::min(wtmax,expxu),expxl) );
00624
00625 G4double bang, cb, sb, s0, s1, s2, c, esys, a, b, gama, beta;
00626 G4double ss;
00627 pcm[0][0] = 0.0;
00628 pcm[1][0] = pd[0];
00629 pcm[2][0] = 0.0;
00630 for( i=1; i<listLen; ++i )
00631 {
00632 pcm[0][i] = 0.0;
00633 pcm[1][i] = -pd[i-1];
00634 pcm[2][i] = 0.0;
00635 bang = twopi*G4UniformRand();
00636 cb = std::cos(bang);
00637 sb = std::sin(bang);
00638 c = 2.0*G4UniformRand() - 1.0;
00639 ss = std::sqrt( std::fabs( 1.0-c*c ) );
00640 if( i < listLen-1 )
00641 {
00642 esys = std::sqrt(pd[i]*pd[i] + emm[i]*emm[i]);
00643 beta = pd[i]/esys;
00644 gama = esys/emm[i];
00645 for( G4int j=0; j<=i; ++j )
00646 {
00647 s0 = pcm[0][j];
00648 s1 = pcm[1][j];
00649 s2 = pcm[2][j];
00650 energy[j] = std::sqrt( s0*s0 + s1*s1 + s2*s2 + mass[j]*mass[j] );
00651 a = s0*c - s1*ss;
00652 pcm[1][j] = s0*ss + s1*c;
00653 b = pcm[2][j];
00654 pcm[0][j] = a*cb - b*sb;
00655 pcm[2][j] = a*sb + b*cb;
00656 pcm[1][j] = gama*(pcm[1][j] + beta*energy[j]);
00657 }
00658 }
00659 else
00660 {
00661 for( G4int j=0; j<=i; ++j )
00662 {
00663 s0 = pcm[0][j];
00664 s1 = pcm[1][j];
00665 s2 = pcm[2][j];
00666 energy[j] = std::sqrt( s0*s0 + s1*s1 + s2*s2 + mass[j]*mass[j] );
00667 a = s0*c - s1*ss;
00668 pcm[1][j] = s0*ss + s1*c;
00669 b = pcm[2][j];
00670 pcm[0][j] = a*cb - b*sb;
00671 pcm[2][j] = a*sb + b*cb;
00672 }
00673 }
00674 }
00675
00676 for (i=0; i<listLen; ++i) {
00677 tempList[i]->SetMomentum(pcm[0][i]*GeV, pcm[1][i]*GeV, pcm[2][i]*GeV);
00678 tempList[i]->SetTotalEnergy(energy[i]*GeV);
00679 }
00680
00681 return weight;
00682 }
00683
00684
00685 G4double G4RPGReaction::normal()
00686 {
00687 G4double ran = -6.0;
00688 for( G4int i=0; i<12; ++i )ran += G4UniformRand();
00689 return ran;
00690 }
00691
00692
00693 void G4RPGReaction::Defs1(const G4ReactionProduct& modifiedOriginal,
00694 G4ReactionProduct& currentParticle,
00695 G4ReactionProduct& targetParticle,
00696 G4FastVector<G4ReactionProduct,256>& vec,
00697 G4int& vecLen)
00698 {
00699
00700
00701 const G4double pjx = modifiedOriginal.GetMomentum().x()/MeV;
00702 const G4double pjy = modifiedOriginal.GetMomentum().y()/MeV;
00703 const G4double pjz = modifiedOriginal.GetMomentum().z()/MeV;
00704 const G4double p = modifiedOriginal.GetMomentum().mag()/MeV;
00705
00706 if (pjx*pjx+pjy*pjy > 0.0) {
00707 G4double cost, sint, ph, cosp, sinp, pix, piy, piz;
00708 cost = pjz/p;
00709 sint = std::sqrt(std::abs((1.0-cost)*(1.0+cost)));
00710 if( pjy < 0.0 )
00711 ph = 3*halfpi;
00712 else
00713 ph = halfpi;
00714 if( std::abs( pjx ) > 0.001*MeV )ph = std::atan2(pjy,pjx);
00715 cosp = std::cos(ph);
00716 sinp = std::sin(ph);
00717 pix = currentParticle.GetMomentum().x()/MeV;
00718 piy = currentParticle.GetMomentum().y()/MeV;
00719 piz = currentParticle.GetMomentum().z()/MeV;
00720 currentParticle.SetMomentum((cost*cosp*pix - sinp*piy + sint*cosp*piz)*MeV,
00721 (cost*sinp*pix + cosp*piy + sint*sinp*piz)*MeV,
00722 (-sint*pix + cost*piz)*MeV);
00723 pix = targetParticle.GetMomentum().x()/MeV;
00724 piy = targetParticle.GetMomentum().y()/MeV;
00725 piz = targetParticle.GetMomentum().z()/MeV;
00726 targetParticle.SetMomentum((cost*cosp*pix - sinp*piy + sint*cosp*piz)*MeV,
00727 (cost*sinp*pix + cosp*piy + sint*sinp*piz)*MeV,
00728 (-sint*pix + cost*piz)*MeV);
00729
00730 for (G4int i=0; i<vecLen; ++i) {
00731 pix = vec[i]->GetMomentum().x()/MeV;
00732 piy = vec[i]->GetMomentum().y()/MeV;
00733 piz = vec[i]->GetMomentum().z()/MeV;
00734 vec[i]->SetMomentum((cost*cosp*pix - sinp*piy + sint*cosp*piz)*MeV,
00735 (cost*sinp*pix + cosp*piy + sint*sinp*piz)*MeV,
00736 (-sint*pix + cost*piz)*MeV);
00737 }
00738
00739 } else {
00740 if (pjz < 0.0) {
00741 currentParticle.SetMomentum( -currentParticle.GetMomentum().z() );
00742 targetParticle.SetMomentum( -targetParticle.GetMomentum().z() );
00743 for (G4int i=0; i<vecLen; ++i) vec[i]->SetMomentum( -vec[i]->GetMomentum().z() );
00744 }
00745 }
00746 }
00747
00748
00749 void G4RPGReaction::Rotate(
00750 const G4double numberofFinalStateNucleons,
00751 const G4ThreeVector &temp,
00752 const G4ReactionProduct &modifiedOriginal,
00753 const G4HadProjectile *originalIncident,
00754 const G4Nucleus &targetNucleus,
00755 G4ReactionProduct ¤tParticle,
00756 G4ReactionProduct &targetParticle,
00757 G4FastVector<G4ReactionProduct,256> &vec,
00758 G4int &vecLen )
00759 {
00760
00761
00762
00763
00764
00765 const G4double atomicWeight = targetNucleus.GetA_asInt();
00766 const G4double logWeight = std::log(atomicWeight);
00767
00768 G4ParticleDefinition *aPiMinus = G4PionMinus::PionMinus();
00769 G4ParticleDefinition *aPiPlus = G4PionPlus::PionPlus();
00770 G4ParticleDefinition *aPiZero = G4PionZero::PionZero();
00771
00772 G4int i;
00773 G4ThreeVector pseudoParticle[4];
00774 for( i=0; i<4; ++i )pseudoParticle[i].set(0,0,0);
00775 pseudoParticle[0] = currentParticle.GetMomentum()
00776 + targetParticle.GetMomentum();
00777 for( i=0; i<vecLen; ++i )
00778 pseudoParticle[0] = pseudoParticle[0] + (vec[i]->GetMomentum());
00779
00780
00781
00782 G4float pp, pp1;
00783 G4double alekw, p;
00784 G4double r1, r2, a1, ran1, ran2, xxh, exh, pxTemp, pyTemp, pzTemp;
00785
00786 r1 = twopi*G4UniformRand();
00787 r2 = G4UniformRand();
00788 a1 = std::sqrt(-2.0*std::log(r2));
00789 ran1 = a1*std::sin(r1)*0.020*numberofFinalStateNucleons*GeV;
00790 ran2 = a1*std::cos(r1)*0.020*numberofFinalStateNucleons*GeV;
00791 G4ThreeVector fermir(ran1, ran2, 0);
00792
00793 pseudoParticle[0] = pseudoParticle[0]+fermir;
00794 pseudoParticle[2] = temp;
00795 pseudoParticle[3] = pseudoParticle[0];
00796
00797 pseudoParticle[1] = pseudoParticle[2].cross(pseudoParticle[3]);
00798 G4double rotation = 2.*pi*G4UniformRand();
00799 pseudoParticle[1] = pseudoParticle[1].rotate(rotation, pseudoParticle[3]);
00800 pseudoParticle[2] = pseudoParticle[3].cross(pseudoParticle[1]);
00801 for(G4int ii=1; ii<=3; ii++)
00802 {
00803 p = pseudoParticle[ii].mag();
00804 if( p == 0.0 )
00805 pseudoParticle[ii]= G4ThreeVector( 0.0, 0.0, 0.0 );
00806 else
00807 pseudoParticle[ii]= pseudoParticle[ii] * (1./p);
00808 }
00809
00810 pxTemp = pseudoParticle[1].dot(currentParticle.GetMomentum());
00811 pyTemp = pseudoParticle[2].dot(currentParticle.GetMomentum());
00812 pzTemp = pseudoParticle[3].dot(currentParticle.GetMomentum());
00813 currentParticle.SetMomentum( pxTemp, pyTemp, pzTemp );
00814
00815 pxTemp = pseudoParticle[1].dot(targetParticle.GetMomentum());
00816 pyTemp = pseudoParticle[2].dot(targetParticle.GetMomentum());
00817 pzTemp = pseudoParticle[3].dot(targetParticle.GetMomentum());
00818 targetParticle.SetMomentum( pxTemp, pyTemp, pzTemp );
00819
00820 for( i=0; i<vecLen; ++i )
00821 {
00822 pxTemp = pseudoParticle[1].dot(vec[i]->GetMomentum());
00823 pyTemp = pseudoParticle[2].dot(vec[i]->GetMomentum());
00824 pzTemp = pseudoParticle[3].dot(vec[i]->GetMomentum());
00825 vec[i]->SetMomentum( pxTemp, pyTemp, pzTemp );
00826 }
00827
00828
00829
00830
00831 Defs1( modifiedOriginal, currentParticle, targetParticle, vec, vecLen );
00832 G4double ekin;
00833 G4double dekin = 0.0;
00834 G4double ek1 = 0.0;
00835 G4int npions = 0;
00836 if( atomicWeight >= 1.5 )
00837 {
00838
00839
00840 const G4double alem[] = { 1.40, 2.30, 2.70, 3.00, 3.40, 4.60, 7.00 };
00841 const G4double val0[] = { 0.00, 0.40, 0.48, 0.51, 0.54, 0.60, 0.65 };
00842 alekw = std::log( originalIncident->GetKineticEnergy()/GeV );
00843 exh = 1.0;
00844 if( alekw > alem[0] )
00845 {
00846 exh = val0[6];
00847 for( G4int j=1; j<7; ++j )
00848 {
00849 if( alekw < alem[j] )
00850 {
00851 G4double rcnve = (val0[j] - val0[j-1]) / (alem[j] - alem[j-1]);
00852 exh = rcnve * alekw + val0[j-1] - rcnve * alem[j-1];
00853 break;
00854 }
00855 }
00856 exh = 1.0 - exh;
00857 }
00858 const G4double cfa = 0.025*((atomicWeight-1.)/120.)*std::exp(-(atomicWeight-1.)/120.);
00859 ekin = currentParticle.GetKineticEnergy()/GeV - cfa*(1+normal()/2.0);
00860 ekin = std::max( 1.0e-6, ekin );
00861 xxh = 1.0;
00862 if( (modifiedOriginal.GetDefinition() == aPiPlus ||
00863 modifiedOriginal.GetDefinition() == aPiMinus) &&
00864 currentParticle.GetDefinition() == aPiZero &&
00865 G4UniformRand() <= logWeight) xxh = exh;
00866 dekin += ekin*(1.0-xxh);
00867 ekin *= xxh;
00868 if (currentParticle.GetDefinition()->GetParticleSubType() == "pi") {
00869 ++npions;
00870 ek1 += ekin;
00871 }
00872 currentParticle.SetKineticEnergy( ekin*GeV );
00873 pp = currentParticle.GetTotalMomentum()/MeV;
00874 pp1 = currentParticle.GetMomentum().mag()/MeV;
00875 if( pp1 < 0.001*MeV )
00876 {
00877 G4double costheta = 2.*G4UniformRand() - 1.;
00878 G4double sintheta = std::sqrt(1. - costheta*costheta);
00879 G4double phi = twopi*G4UniformRand();
00880 currentParticle.SetMomentum( pp*sintheta*std::cos(phi)*MeV,
00881 pp*sintheta*std::sin(phi)*MeV,
00882 pp*costheta*MeV ) ;
00883 }
00884 else
00885 currentParticle.SetMomentum( currentParticle.GetMomentum() * (pp/pp1) );
00886 ekin = targetParticle.GetKineticEnergy()/GeV - cfa*(1+normal()/2.0);
00887 ekin = std::max( 1.0e-6, ekin );
00888 xxh = 1.0;
00889 if( (modifiedOriginal.GetDefinition() == aPiPlus ||
00890 modifiedOriginal.GetDefinition() == aPiMinus) &&
00891 targetParticle.GetDefinition() == aPiZero &&
00892 G4UniformRand() < logWeight) xxh = exh;
00893 dekin += ekin*(1.0-xxh);
00894 ekin *= xxh;
00895 if (targetParticle.GetDefinition()->GetParticleSubType() == "pi") {
00896 ++npions;
00897 ek1 += ekin;
00898 }
00899 targetParticle.SetKineticEnergy( ekin*GeV );
00900 pp = targetParticle.GetTotalMomentum()/MeV;
00901 pp1 = targetParticle.GetMomentum().mag()/MeV;
00902 if( pp1 < 0.001*MeV )
00903 {
00904 G4double costheta = 2.*G4UniformRand() - 1.;
00905 G4double sintheta = std::sqrt(1. - costheta*costheta);
00906 G4double phi = twopi*G4UniformRand();
00907 targetParticle.SetMomentum( pp*sintheta*std::cos(phi)*MeV,
00908 pp*sintheta*std::sin(phi)*MeV,
00909 pp*costheta*MeV ) ;
00910 }
00911 else
00912 targetParticle.SetMomentum( targetParticle.GetMomentum() * (pp/pp1) );
00913 for( i=0; i<vecLen; ++i )
00914 {
00915 ekin = vec[i]->GetKineticEnergy()/GeV - cfa*(1+normal()/2.0);
00916 ekin = std::max( 1.0e-6, ekin );
00917 xxh = 1.0;
00918 if( (modifiedOriginal.GetDefinition() == aPiPlus ||
00919 modifiedOriginal.GetDefinition() == aPiMinus) &&
00920 vec[i]->GetDefinition() == aPiZero &&
00921 G4UniformRand() < logWeight) xxh = exh;
00922 dekin += ekin*(1.0-xxh);
00923 ekin *= xxh;
00924 if (vec[i]->GetDefinition()->GetParticleSubType() == "pi") {
00925 ++npions;
00926 ek1 += ekin;
00927 }
00928 vec[i]->SetKineticEnergy( ekin*GeV );
00929 pp = vec[i]->GetTotalMomentum()/MeV;
00930 pp1 = vec[i]->GetMomentum().mag()/MeV;
00931 if( pp1 < 0.001*MeV )
00932 {
00933 G4double costheta = 2.*G4UniformRand() - 1.;
00934 G4double sintheta = std::sqrt(1. - costheta*costheta);
00935 G4double phi = twopi*G4UniformRand();
00936 vec[i]->SetMomentum( pp*sintheta*std::cos(phi)*MeV,
00937 pp*sintheta*std::sin(phi)*MeV,
00938 pp*costheta*MeV ) ;
00939 }
00940 else
00941 vec[i]->SetMomentum( vec[i]->GetMomentum() * (pp/pp1) );
00942 }
00943 }
00944 if( (ek1 != 0.0) && (npions > 0) )
00945 {
00946 dekin = 1.0 + dekin/ek1;
00947
00948
00949
00950 if (currentParticle.GetDefinition()->GetParticleSubType() == "pi")
00951 {
00952 currentParticle.SetKineticEnergy(
00953 std::max( 0.001*MeV, dekin*currentParticle.GetKineticEnergy() ) );
00954 pp = currentParticle.GetTotalMomentum()/MeV;
00955 pp1 = currentParticle.GetMomentum().mag()/MeV;
00956 if( pp1 < 0.001 )
00957 {
00958 G4double costheta = 2.*G4UniformRand() - 1.;
00959 G4double sintheta = std::sqrt(1. - costheta*costheta);
00960 G4double phi = twopi*G4UniformRand();
00961 currentParticle.SetMomentum( pp*sintheta*std::cos(phi)*MeV,
00962 pp*sintheta*std::sin(phi)*MeV,
00963 pp*costheta*MeV ) ;
00964 } else {
00965 currentParticle.SetMomentum( currentParticle.GetMomentum() * (pp/pp1) );
00966 }
00967 }
00968
00969 if (targetParticle.GetDefinition()->GetParticleSubType() == "pi")
00970 {
00971 targetParticle.SetKineticEnergy(
00972 std::max( 0.001*MeV, dekin*targetParticle.GetKineticEnergy() ) );
00973 pp = targetParticle.GetTotalMomentum()/MeV;
00974 pp1 = targetParticle.GetMomentum().mag()/MeV;
00975 if( pp1 < 0.001 )
00976 {
00977 G4double costheta = 2.*G4UniformRand() - 1.;
00978 G4double sintheta = std::sqrt(1. - costheta*costheta);
00979 G4double phi = twopi*G4UniformRand();
00980 targetParticle.SetMomentum( pp*sintheta*std::cos(phi)*MeV,
00981 pp*sintheta*std::sin(phi)*MeV,
00982 pp*costheta*MeV ) ;
00983 } else {
00984 targetParticle.SetMomentum( targetParticle.GetMomentum() * (pp/pp1) );
00985 }
00986 }
00987
00988 for( i=0; i<vecLen; ++i )
00989 {
00990 if (vec[i]->GetDefinition()->GetParticleSubType() == "pi")
00991 {
00992 vec[i]->SetKineticEnergy( std::max( 0.001*MeV, dekin*vec[i]->GetKineticEnergy() ) );
00993 pp = vec[i]->GetTotalMomentum()/MeV;
00994 pp1 = vec[i]->GetMomentum().mag()/MeV;
00995 if( pp1 < 0.001 )
00996 {
00997 G4double costheta = 2.*G4UniformRand() - 1.;
00998 G4double sintheta = std::sqrt(1. - costheta*costheta);
00999 G4double phi = twopi*G4UniformRand();
01000 vec[i]->SetMomentum( pp*sintheta*std::cos(phi)*MeV,
01001 pp*sintheta*std::sin(phi)*MeV,
01002 pp*costheta*MeV ) ;
01003 } else {
01004 vec[i]->SetMomentum( vec[i]->GetMomentum() * (pp/pp1) );
01005 }
01006 }
01007 }
01008 }
01009 }
01010
01011 std::pair<G4int, G4int> G4RPGReaction::GetFinalStateNucleons(
01012 const G4DynamicParticle* originalTarget,
01013 const G4FastVector<G4ReactionProduct,256>& vec,
01014 const G4int& vecLen)
01015 {
01016
01017
01018 G4int protonsRemoved = 0;
01019 G4int neutronsRemoved = 0;
01020 if (originalTarget->GetDefinition()->GetParticleName() == "proton")
01021 protonsRemoved++;
01022 else
01023 neutronsRemoved++;
01024
01025 G4String secName;
01026 for (G4int i = 0; i < vecLen; i++) {
01027 secName = vec[i]->GetDefinition()->GetParticleName();
01028 if (secName == "proton") {
01029 protonsRemoved++;
01030 } else if (secName == "neutron") {
01031 neutronsRemoved++;
01032 } else if (secName == "anti_proton") {
01033 protonsRemoved--;
01034 } else if (secName == "anti_neutron") {
01035 neutronsRemoved--;
01036 }
01037 }
01038
01039 return std::pair<G4int, G4int>(protonsRemoved, neutronsRemoved);
01040 }
01041
01042
01043 G4ThreeVector G4RPGReaction::Isotropic(const G4double& pp)
01044 {
01045 G4double costheta = 2.*G4UniformRand() - 1.;
01046 G4double sintheta = std::sqrt(1. - costheta*costheta);
01047 G4double phi = twopi*G4UniformRand();
01048 return G4ThreeVector(pp*sintheta*std::cos(phi),
01049 pp*sintheta*std::sin(phi),
01050 pp*costheta);
01051 }
01052
01053
01054 void G4RPGReaction::MomentumCheck(
01055 const G4ReactionProduct &modifiedOriginal,
01056 G4ReactionProduct ¤tParticle,
01057 G4ReactionProduct &targetParticle,
01058 G4FastVector<G4ReactionProduct,256> &vec,
01059 G4int &vecLen )
01060 {
01061 const G4double pOriginal = modifiedOriginal.GetTotalMomentum()/MeV;
01062 G4double testMomentum = currentParticle.GetMomentum().mag()/MeV;
01063 G4double pMass;
01064 if( testMomentum >= pOriginal )
01065 {
01066 pMass = currentParticle.GetMass()/MeV;
01067 currentParticle.SetTotalEnergy(
01068 std::sqrt( pMass*pMass + pOriginal*pOriginal )*MeV );
01069 currentParticle.SetMomentum(
01070 currentParticle.GetMomentum() * (pOriginal/testMomentum) );
01071 }
01072 testMomentum = targetParticle.GetMomentum().mag()/MeV;
01073 if( testMomentum >= pOriginal )
01074 {
01075 pMass = targetParticle.GetMass()/MeV;
01076 targetParticle.SetTotalEnergy(
01077 std::sqrt( pMass*pMass + pOriginal*pOriginal )*MeV );
01078 targetParticle.SetMomentum(
01079 targetParticle.GetMomentum() * (pOriginal/testMomentum) );
01080 }
01081 for( G4int i=0; i<vecLen; ++i )
01082 {
01083 testMomentum = vec[i]->GetMomentum().mag()/MeV;
01084 if( testMomentum >= pOriginal )
01085 {
01086 pMass = vec[i]->GetMass()/MeV;
01087 vec[i]->SetTotalEnergy(
01088 std::sqrt( pMass*pMass + pOriginal*pOriginal )*MeV );
01089 vec[i]->SetMomentum( vec[i]->GetMomentum() * (pOriginal/testMomentum) );
01090 }
01091 }
01092 }
01093
01094 void G4RPGReaction::NuclearReaction(
01095 G4FastVector<G4ReactionProduct,4> &vec,
01096 G4int &vecLen,
01097 const G4HadProjectile *originalIncident,
01098 const G4Nucleus &targetNucleus,
01099 const G4double theAtomicMass,
01100 const G4double *mass )
01101 {
01102
01103
01104
01105
01106 G4ParticleDefinition *aGamma = G4Gamma::Gamma();
01107 G4ParticleDefinition *aProton = G4Proton::Proton();
01108 G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
01109 G4ParticleDefinition *aDeuteron = G4Deuteron::Deuteron();
01110 G4ParticleDefinition *aTriton = G4Triton::Triton();
01111 G4ParticleDefinition *anAlpha = G4Alpha::Alpha();
01112
01113 const G4double aProtonMass = aProton->GetPDGMass()/MeV;
01114 const G4double aNeutronMass = aNeutron->GetPDGMass()/MeV;
01115 const G4double aDeuteronMass = aDeuteron->GetPDGMass()/MeV;
01116 const G4double aTritonMass = aTriton->GetPDGMass()/MeV;
01117 const G4double anAlphaMass = anAlpha->GetPDGMass()/MeV;
01118
01119 G4ReactionProduct currentParticle;
01120 currentParticle = *originalIncident;
01121
01122
01123
01124
01125
01126 G4double p = currentParticle.GetTotalMomentum();
01127 G4double pp = currentParticle.GetMomentum().mag();
01128 if( pp <= 0.001*MeV )
01129 {
01130 G4double phinve = twopi*G4UniformRand();
01131 G4double rthnve = std::acos( std::max( -1.0, std::min( 1.0, -1.0 + 2.0*G4UniformRand() ) ) );
01132 currentParticle.SetMomentum( p*std::sin(rthnve)*std::cos(phinve),
01133 p*std::sin(rthnve)*std::sin(phinve),
01134 p*std::cos(rthnve) );
01135 }
01136 else
01137 currentParticle.SetMomentum( currentParticle.GetMomentum() * (p/pp) );
01138
01139
01140
01141 G4double currentKinetic = currentParticle.GetKineticEnergy()/MeV;
01142 G4double currentMass = currentParticle.GetDefinition()->GetPDGMass()/MeV;
01143 G4double qv = currentKinetic + theAtomicMass + currentMass;
01144
01145 G4double qval[9];
01146 qval[0] = qv - mass[0];
01147 qval[1] = qv - mass[1] - aNeutronMass;
01148 qval[2] = qv - mass[2] - aProtonMass;
01149 qval[3] = qv - mass[3] - aDeuteronMass;
01150 qval[4] = qv - mass[4] - aTritonMass;
01151 qval[5] = qv - mass[5] - anAlphaMass;
01152 qval[6] = qv - mass[6] - aNeutronMass - aNeutronMass;
01153 qval[7] = qv - mass[7] - aNeutronMass - aProtonMass;
01154 qval[8] = qv - mass[8] - aProtonMass - aProtonMass;
01155
01156 if( currentParticle.GetDefinition() == aNeutron )
01157 {
01158 const G4double A = targetNucleus.GetA_asInt();
01159 if( G4UniformRand() > ((A-1.0)/230.0)*((A-1.0)/230.0) )
01160 qval[0] = 0.0;
01161 if( G4UniformRand() >= currentKinetic/7.9254*A )
01162 qval[2] = qval[3] = qval[4] = qval[5] = qval[8] = 0.0;
01163 }
01164 else
01165 qval[0] = 0.0;
01166
01167 G4int i;
01168 qv = 0.0;
01169 for( i=0; i<9; ++i )
01170 {
01171 if( mass[i] < 500.0*MeV )qval[i] = 0.0;
01172 if( qval[i] < 0.0 )qval[i] = 0.0;
01173 qv += qval[i];
01174 }
01175 G4double qv1 = 0.0;
01176 G4double ran = G4UniformRand();
01177 G4int index;
01178 for( index=0; index<9; ++index )
01179 {
01180 if( qval[index] > 0.0 )
01181 {
01182 qv1 += qval[index]/qv;
01183 if( ran <= qv1 )break;
01184 }
01185 }
01186 if( index == 9 )
01187 {
01188 throw G4HadronicException(__FILE__, __LINE__,
01189 "G4RPGReaction::NuclearReaction: inelastic reaction kinematically not possible");
01190 }
01191 G4double ke = currentParticle.GetKineticEnergy()/GeV;
01192 G4int nt = 2;
01193 if( (index>=6) || (G4UniformRand()<std::min(0.5,ke*10.0)) )nt = 3;
01194
01195 G4ReactionProduct **v = new G4ReactionProduct * [3];
01196 v[0] = new G4ReactionProduct;
01197 v[1] = new G4ReactionProduct;
01198 v[2] = new G4ReactionProduct;
01199
01200 v[0]->SetMass( mass[index]*MeV );
01201 switch( index )
01202 {
01203 case 0:
01204 v[1]->SetDefinition( aGamma );
01205 v[2]->SetDefinition( aGamma );
01206 break;
01207 case 1:
01208 v[1]->SetDefinition( aNeutron );
01209 v[2]->SetDefinition( aGamma );
01210 break;
01211 case 2:
01212 v[1]->SetDefinition( aProton );
01213 v[2]->SetDefinition( aGamma );
01214 break;
01215 case 3:
01216 v[1]->SetDefinition( aDeuteron );
01217 v[2]->SetDefinition( aGamma );
01218 break;
01219 case 4:
01220 v[1]->SetDefinition( aTriton );
01221 v[2]->SetDefinition( aGamma );
01222 break;
01223 case 5:
01224 v[1]->SetDefinition( anAlpha );
01225 v[2]->SetDefinition( aGamma );
01226 break;
01227 case 6:
01228 v[1]->SetDefinition( aNeutron );
01229 v[2]->SetDefinition( aNeutron );
01230 break;
01231 case 7:
01232 v[1]->SetDefinition( aNeutron );
01233 v[2]->SetDefinition( aProton );
01234 break;
01235 case 8:
01236 v[1]->SetDefinition( aProton );
01237 v[2]->SetDefinition( aProton );
01238 break;
01239 }
01240
01241
01242
01243 G4ReactionProduct pseudo1;
01244 pseudo1.SetMass( theAtomicMass*MeV );
01245 pseudo1.SetTotalEnergy( theAtomicMass*MeV );
01246 G4ReactionProduct pseudo2 = currentParticle + pseudo1;
01247 pseudo2.SetMomentum( pseudo2.GetMomentum() * (-1.0) );
01248
01249
01250
01251 G4FastVector<G4ReactionProduct,256> tempV;
01252 tempV.Initialize( nt );
01253 G4int tempLen = 0;
01254 tempV.SetElement( tempLen++, v[0] );
01255 tempV.SetElement( tempLen++, v[1] );
01256 if( nt == 3 )tempV.SetElement( tempLen++, v[2] );
01257 G4bool constantCrossSection = true;
01258 GenerateNBodyEvent( pseudo2.GetMass()/MeV, constantCrossSection, tempV, tempLen );
01259 v[0]->Lorentz( *v[0], pseudo2 );
01260 v[1]->Lorentz( *v[1], pseudo2 );
01261 if( nt == 3 )v[2]->Lorentz( *v[2], pseudo2 );
01262
01263 G4bool particleIsDefined = false;
01264 if( v[0]->GetMass()/MeV - aProtonMass < 0.1 )
01265 {
01266 v[0]->SetDefinition( aProton );
01267 particleIsDefined = true;
01268 }
01269 else if( v[0]->GetMass()/MeV - aNeutronMass < 0.1 )
01270 {
01271 v[0]->SetDefinition( aNeutron );
01272 particleIsDefined = true;
01273 }
01274 else if( v[0]->GetMass()/MeV - aDeuteronMass < 0.1 )
01275 {
01276 v[0]->SetDefinition( aDeuteron );
01277 particleIsDefined = true;
01278 }
01279 else if( v[0]->GetMass()/MeV - aTritonMass < 0.1 )
01280 {
01281 v[0]->SetDefinition( aTriton );
01282 particleIsDefined = true;
01283 }
01284 else if( v[0]->GetMass()/MeV - anAlphaMass < 0.1 )
01285 {
01286 v[0]->SetDefinition( anAlpha );
01287 particleIsDefined = true;
01288 }
01289 currentParticle.SetKineticEnergy(
01290 std::max( 0.001, currentParticle.GetKineticEnergy()/MeV ) );
01291 p = currentParticle.GetTotalMomentum();
01292 pp = currentParticle.GetMomentum().mag();
01293 if( pp <= 0.001*MeV )
01294 {
01295 G4double phinve = twopi*G4UniformRand();
01296 G4double rthnve = std::acos( std::max( -1.0, std::min( 1.0, -1.0 + 2.0*G4UniformRand() ) ) );
01297 currentParticle.SetMomentum( p*std::sin(rthnve)*std::cos(phinve),
01298 p*std::sin(rthnve)*std::sin(phinve),
01299 p*std::cos(rthnve) );
01300 }
01301 else
01302 currentParticle.SetMomentum( currentParticle.GetMomentum() * (p/pp) );
01303
01304 if( particleIsDefined )
01305 {
01306 v[0]->SetKineticEnergy(
01307 std::max( 0.001, 0.5*G4UniformRand()*v[0]->GetKineticEnergy()/MeV ) );
01308 p = v[0]->GetTotalMomentum();
01309 pp = v[0]->GetMomentum().mag();
01310 if( pp <= 0.001*MeV )
01311 {
01312 G4double phinve = twopi*G4UniformRand();
01313 G4double rthnve = std::acos( std::max(-1.0,std::min(1.0,-1.0+2.0*G4UniformRand())) );
01314 v[0]->SetMomentum( p*std::sin(rthnve)*std::cos(phinve),
01315 p*std::sin(rthnve)*std::sin(phinve),
01316 p*std::cos(rthnve) );
01317 }
01318 else
01319 v[0]->SetMomentum( v[0]->GetMomentum() * (p/pp) );
01320 }
01321 if( (v[1]->GetDefinition() == aDeuteron) ||
01322 (v[1]->GetDefinition() == aTriton) ||
01323 (v[1]->GetDefinition() == anAlpha) )
01324 v[1]->SetKineticEnergy(
01325 std::max( 0.001, 0.5*G4UniformRand()*v[1]->GetKineticEnergy()/MeV ) );
01326 else
01327 v[1]->SetKineticEnergy( std::max( 0.001, v[1]->GetKineticEnergy()/MeV ) );
01328
01329 p = v[1]->GetTotalMomentum();
01330 pp = v[1]->GetMomentum().mag();
01331 if( pp <= 0.001*MeV )
01332 {
01333 G4double phinve = twopi*G4UniformRand();
01334 G4double rthnve = std::acos( std::max(-1.0,std::min(1.0,-1.0+2.0*G4UniformRand())) );
01335 v[1]->SetMomentum( p*std::sin(rthnve)*std::cos(phinve),
01336 p*std::sin(rthnve)*std::sin(phinve),
01337 p*std::cos(rthnve) );
01338 }
01339 else
01340 v[1]->SetMomentum( v[1]->GetMomentum() * (p/pp) );
01341
01342 if( nt == 3 )
01343 {
01344 if( (v[2]->GetDefinition() == aDeuteron) ||
01345 (v[2]->GetDefinition() == aTriton) ||
01346 (v[2]->GetDefinition() == anAlpha) )
01347 v[2]->SetKineticEnergy(
01348 std::max( 0.001, 0.5*G4UniformRand()*v[2]->GetKineticEnergy()/MeV ) );
01349 else
01350 v[2]->SetKineticEnergy( std::max( 0.001, v[2]->GetKineticEnergy()/MeV ) );
01351
01352 p = v[2]->GetTotalMomentum();
01353 pp = v[2]->GetMomentum().mag();
01354 if( pp <= 0.001*MeV )
01355 {
01356 G4double phinve = twopi*G4UniformRand();
01357 G4double rthnve = std::acos( std::max(-1.0,std::min(1.0,-1.0+2.0*G4UniformRand())) );
01358 v[2]->SetMomentum( p*std::sin(rthnve)*std::cos(phinve),
01359 p*std::sin(rthnve)*std::sin(phinve),
01360 p*std::cos(rthnve) );
01361 }
01362 else
01363 v[2]->SetMomentum( v[2]->GetMomentum() * (p/pp) );
01364 }
01365 G4int del;
01366 for(del=0; del<vecLen; del++) delete vec[del];
01367 vecLen = 0;
01368 if( particleIsDefined )
01369 {
01370 vec.SetElement( vecLen++, v[0] );
01371 }
01372 else
01373 {
01374 delete v[0];
01375 }
01376 vec.SetElement( vecLen++, v[1] );
01377 if( nt == 3 )
01378 {
01379 vec.SetElement( vecLen++, v[2] );
01380 }
01381 else
01382 {
01383 delete v[2];
01384 }
01385 delete [] v;
01386 return;
01387 }
01388
01389
01390