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 "G4RPGAntiLambdaInelastic.hh"
00030 #include "G4PhysicalConstants.hh"
00031 #include "G4SystemOfUnits.hh"
00032 #include "Randomize.hh"
00033
00034 G4HadFinalState*
00035 G4RPGAntiLambdaInelastic::ApplyYourself( const G4HadProjectile &aTrack,
00036 G4Nucleus &targetNucleus )
00037 {
00038 const G4HadProjectile *originalIncident = &aTrack;
00039
00040
00041
00042 G4DynamicParticle *originalTarget = targetNucleus.ReturnTargetParticle();
00043
00044 if( verboseLevel > 1 )
00045 {
00046 const G4Material *targetMaterial = aTrack.GetMaterial();
00047 G4cout << "G4RPGAntiLambdaInelastic::ApplyYourself called" << G4endl;
00048 G4cout << "kinetic energy = " << originalIncident->GetKineticEnergy()/MeV << "MeV, ";
00049 G4cout << "target material = " << targetMaterial->GetName() << ", ";
00050 G4cout << "target particle = " << originalTarget->GetDefinition()->GetParticleName()
00051 << G4endl;
00052 }
00053
00054
00055
00056
00057 G4double ek = originalIncident->GetKineticEnergy()/MeV;
00058 G4double amas = originalIncident->GetDefinition()->GetPDGMass()/MeV;
00059 G4ReactionProduct modifiedOriginal;
00060 modifiedOriginal = *originalIncident;
00061
00062 G4double tkin = targetNucleus.Cinema( ek );
00063 ek += tkin;
00064 modifiedOriginal.SetKineticEnergy( ek*MeV );
00065 G4double et = ek + amas;
00066 G4double p = std::sqrt( std::abs((et-amas)*(et+amas)) );
00067 G4double pp = modifiedOriginal.GetMomentum().mag()/MeV;
00068 if( pp > 0.0 )
00069 {
00070 G4ThreeVector momentum = modifiedOriginal.GetMomentum();
00071 modifiedOriginal.SetMomentum( momentum * (p/pp) );
00072 }
00073
00074
00075
00076 tkin = targetNucleus.EvaporationEffects( ek );
00077 ek -= tkin;
00078 modifiedOriginal.SetKineticEnergy( ek*MeV );
00079 et = ek + amas;
00080 p = std::sqrt( std::abs((et-amas)*(et+amas)) );
00081 pp = modifiedOriginal.GetMomentum().mag()/MeV;
00082 if( pp > 0.0 )
00083 {
00084 G4ThreeVector momentum = modifiedOriginal.GetMomentum();
00085 modifiedOriginal.SetMomentum( momentum * (p/pp) );
00086 }
00087
00088 G4ReactionProduct currentParticle = modifiedOriginal;
00089 G4ReactionProduct targetParticle;
00090 targetParticle = *originalTarget;
00091 currentParticle.SetSide( 1 );
00092 targetParticle.SetSide( -1 );
00093 G4bool incidentHasChanged = false;
00094 G4bool targetHasChanged = false;
00095 G4bool quasiElastic = false;
00096 G4FastVector<G4ReactionProduct,GHADLISTSIZE> vec;
00097 G4int vecLen = 0;
00098 vec.Initialize( 0 );
00099
00100 const G4double cutOff = 0.1;
00101 const G4double anni = std::min( 1.3*currentParticle.GetTotalMomentum()/GeV, 0.4 );
00102 if( (originalIncident->GetKineticEnergy()/MeV > cutOff) || (G4UniformRand() > anni) )
00103 Cascade( vec, vecLen,
00104 originalIncident, currentParticle, targetParticle,
00105 incidentHasChanged, targetHasChanged, quasiElastic );
00106
00107 CalculateMomenta( vec, vecLen,
00108 originalIncident, originalTarget, modifiedOriginal,
00109 targetNucleus, currentParticle, targetParticle,
00110 incidentHasChanged, targetHasChanged, quasiElastic );
00111
00112 SetUpChange( vec, vecLen,
00113 currentParticle, targetParticle,
00114 incidentHasChanged );
00115
00116 delete originalTarget;
00117 return &theParticleChange;
00118 }
00119
00120
00121 void G4RPGAntiLambdaInelastic::Cascade(
00122 G4FastVector<G4ReactionProduct,GHADLISTSIZE> &vec,
00123 G4int &vecLen,
00124 const G4HadProjectile *originalIncident,
00125 G4ReactionProduct ¤tParticle,
00126 G4ReactionProduct &targetParticle,
00127 G4bool &incidentHasChanged,
00128 G4bool &targetHasChanged,
00129 G4bool &quasiElastic )
00130 {
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140 const G4double mOriginal = originalIncident->GetDefinition()->GetPDGMass()/MeV;
00141 const G4double etOriginal = originalIncident->GetTotalEnergy()/MeV;
00142 const G4double targetMass = targetParticle.GetMass()/MeV;
00143 const G4double pOriginal = originalIncident->GetTotalMomentum()/GeV;
00144 G4double centerofmassEnergy = std::sqrt( mOriginal*mOriginal +
00145 targetMass*targetMass +
00146 2.0*targetMass*etOriginal );
00147 G4double availableEnergy = centerofmassEnergy-(targetMass+mOriginal);
00148
00149 static G4bool first = true;
00150 const G4int numMul = 1200;
00151 const G4int numMulA = 400;
00152 const G4int numSec = 60;
00153 static G4double protmul[numMul], protnorm[numSec];
00154 static G4double neutmul[numMul], neutnorm[numSec];
00155 static G4double protmulA[numMulA], protnormA[numSec];
00156 static G4double neutmulA[numMulA], neutnormA[numSec];
00157
00158 G4int nt=0, np=0, nneg=0, nz=0;
00159 G4double test;
00160 const G4double c = 1.25;
00161 const G4double b[] = { 0.7, 0.7 };
00162 if( first )
00163 {
00164 first = false;
00165 G4int i;
00166 for( i=0; i<numMul; ++i )protmul[i] = 0.0;
00167 for( i=0; i<numSec; ++i )protnorm[i] = 0.0;
00168 G4int counter = -1;
00169 for( np=0; np<(numSec/3); ++np )
00170 {
00171 for( nneg=std::max(0,np-2); nneg<=(np+1); ++nneg )
00172 {
00173 for( nz=0; nz<numSec/3; ++nz )
00174 {
00175 if( ++counter < numMul )
00176 {
00177 nt = np+nneg+nz;
00178 if( nt>0 && nt<=numSec )
00179 {
00180 protmul[counter] = Pmltpc(np,nneg,nz,nt,b[0],c);
00181 protnorm[nt-1] += protmul[counter];
00182 }
00183 }
00184 }
00185 }
00186 }
00187 for( i=0; i<numMul; ++i )neutmul[i] = 0.0;
00188 for( i=0; i<numSec; ++i )neutnorm[i] = 0.0;
00189 counter = -1;
00190 for( np=0; np<numSec/3; ++np )
00191 {
00192 for( nneg=std::max(0,np-1); nneg<=(np+2); ++nneg )
00193 {
00194 for( nz=0; nz<numSec/3; ++nz )
00195 {
00196 if( ++counter < numMul )
00197 {
00198 nt = np+nneg+nz;
00199 if( nt>0 && nt<=numSec )
00200 {
00201 neutmul[counter] = Pmltpc(np,nneg,nz,nt,b[1],c);
00202 neutnorm[nt-1] += neutmul[counter];
00203 }
00204 }
00205 }
00206 }
00207 }
00208 for( i=0; i<numSec; ++i )
00209 {
00210 if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
00211 if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
00212 }
00213
00214
00215
00216 for( i=0; i<numMulA; ++i )protmulA[i] = 0.0;
00217 for( i=0; i<numSec; ++i )protnormA[i] = 0.0;
00218 counter = -1;
00219 for( np=1; np<(numSec/3); ++np )
00220 {
00221 nneg = np-1;
00222 for( nz=0; nz<numSec/3; ++nz )
00223 {
00224 if( ++counter < numMulA )
00225 {
00226 nt = np+nneg+nz;
00227 if( nt>1 && nt<=numSec )
00228 {
00229 protmulA[counter] = Pmltpc(np,nneg,nz,nt,b[0],c);
00230 protnormA[nt-1] += protmulA[counter];
00231 }
00232 }
00233 }
00234 }
00235 for( i=0; i<numMulA; ++i )neutmulA[i] = 0.0;
00236 for( i=0; i<numSec; ++i )neutnormA[i] = 0.0;
00237 counter = -1;
00238 for( np=0; np<numSec/3; ++np )
00239 {
00240 nneg = np;
00241 for( nz=0; nz<numSec/3; ++nz )
00242 {
00243 if( ++counter < numMulA )
00244 {
00245 nt = np+nneg+nz;
00246 if( nt>1 && nt<=numSec )
00247 {
00248 neutmulA[counter] = Pmltpc(np,nneg,nz,nt,b[1],c);
00249 neutnormA[nt-1] += neutmulA[counter];
00250 }
00251 }
00252 }
00253 }
00254 for( i=0; i<numSec; ++i )
00255 {
00256 if( protnormA[i] > 0.0 )protnormA[i] = 1.0/protnormA[i];
00257 if( neutnormA[i] > 0.0 )neutnormA[i] = 1.0/neutnormA[i];
00258 }
00259 }
00260 const G4double expxu = 82.;
00261 const G4double expxl = -expxu;
00262
00263 G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
00264 G4ParticleDefinition *aProton = G4Proton::Proton();
00265 G4ParticleDefinition *aPiPlus = G4PionPlus::PionPlus();
00266 G4ParticleDefinition *aPiMinus = G4PionMinus::PionMinus();
00267 G4ParticleDefinition *aPiZero = G4PionZero::PionZero();
00268 G4ParticleDefinition *aKaonPlus = G4KaonPlus::KaonPlus();
00269 G4ParticleDefinition *aKaonMinus = G4KaonMinus::KaonMinus();
00270 G4ParticleDefinition *aKaonZL = G4KaonZeroLong::KaonZeroLong();
00271 G4ParticleDefinition *anAntiSigmaZero = G4AntiSigmaZero::AntiSigmaZero();
00272 G4ParticleDefinition *anAntiSigmaPlus = G4AntiSigmaPlus::AntiSigmaPlus();
00273 G4ParticleDefinition *anAntiSigmaMinus = G4AntiSigmaMinus::AntiSigmaMinus();
00274 const G4double anhl[] = {1.00,1.00,1.00,1.00,1.00,1.00,1.00,1.00,0.97,0.88,
00275 0.85,0.81,0.75,0.64,0.64,0.55,0.55,0.45,0.47,0.40,
00276 0.39,0.36,0.33,0.10,0.01};
00277 G4int iplab = G4int( pOriginal*10.0 );
00278 if( iplab > 9 )iplab = G4int( (pOriginal- 1.0)*5.0 ) + 10;
00279 if( iplab > 14 )iplab = G4int( pOriginal- 2.0 ) + 15;
00280 if( iplab > 22 )iplab = G4int( (pOriginal-10.0)/10.0 ) + 23;
00281 if( iplab > 24 )iplab = 24;
00282 if( G4UniformRand() > anhl[iplab] )
00283 {
00284 if( availableEnergy <= aPiPlus->GetPDGMass()/MeV )
00285 {
00286 quasiElastic = true;
00287 return;
00288 }
00289 G4double n, anpn;
00290 GetNormalizationConstant( availableEnergy, n, anpn );
00291 G4double ran = G4UniformRand();
00292 G4double dum, excs = 0.0;
00293 if( targetParticle.GetDefinition() == aProton )
00294 {
00295 G4int counter = -1;
00296 for( np=0; np<numSec/3 && ran>=excs; ++np )
00297 {
00298 for( nneg=std::max(0,np-2); nneg<=(np+1) && ran>=excs; ++nneg )
00299 {
00300 for( nz=0; nz<numSec/3 && ran>=excs; ++nz )
00301 {
00302 if( ++counter < numMul )
00303 {
00304 nt = np+nneg+nz;
00305 if( nt>0 && nt<=numSec )
00306 {
00307 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00308 dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
00309 if( std::fabs(dum) < 1.0 )
00310 {
00311 if( test >= 1.0e-10 )excs += dum*test;
00312 }
00313 else
00314 excs += dum*test;
00315 }
00316 }
00317 }
00318 }
00319 }
00320 if( ran >= excs )
00321 {
00322 quasiElastic = true;
00323 return;
00324 }
00325 np--; nneg--; nz--;
00326 G4int ncht = std::min( 4, std::max( 1, np-nneg+2 ) );
00327 switch( ncht )
00328 {
00329 case 1:
00330 currentParticle.SetDefinitionAndUpdateE( anAntiSigmaMinus );
00331 incidentHasChanged = true;
00332 break;
00333 case 2:
00334 if( G4UniformRand() < 0.5 )
00335 {
00336 if( G4UniformRand() < 0.5 )
00337 {
00338 currentParticle.SetDefinitionAndUpdateE( anAntiSigmaZero );
00339 incidentHasChanged = true;
00340 }
00341 else
00342 {
00343 currentParticle.SetDefinitionAndUpdateE( anAntiSigmaMinus );
00344 incidentHasChanged = true;
00345 targetParticle.SetDefinitionAndUpdateE( aNeutron );
00346 targetHasChanged = true;
00347 }
00348 }
00349 else
00350 {
00351 if( G4UniformRand() >= 0.5 )
00352 {
00353 currentParticle.SetDefinitionAndUpdateE( anAntiSigmaMinus );
00354 incidentHasChanged = true;
00355 targetParticle.SetDefinitionAndUpdateE( aNeutron );
00356 targetHasChanged = true;
00357 }
00358 }
00359 break;
00360 case 3:
00361 if( G4UniformRand() < 0.5 )
00362 {
00363 if( G4UniformRand() < 0.5 )
00364 {
00365 currentParticle.SetDefinitionAndUpdateE( anAntiSigmaZero );
00366 incidentHasChanged = true;
00367 targetParticle.SetDefinitionAndUpdateE( aNeutron );
00368 targetHasChanged = true;
00369 }
00370 else
00371 {
00372 currentParticle.SetDefinitionAndUpdateE( anAntiSigmaPlus );
00373 incidentHasChanged = true;
00374 }
00375 }
00376 else
00377 {
00378 if( G4UniformRand() < 0.5 )
00379 {
00380 targetParticle.SetDefinitionAndUpdateE( aNeutron );
00381 targetHasChanged = true;
00382 }
00383 else
00384 {
00385 currentParticle.SetDefinitionAndUpdateE( anAntiSigmaPlus );
00386 incidentHasChanged = true;
00387 }
00388 }
00389 break;
00390 case 4:
00391 currentParticle.SetDefinitionAndUpdateE( anAntiSigmaPlus );
00392 incidentHasChanged = true;
00393 targetParticle.SetDefinitionAndUpdateE( aNeutron );
00394 targetHasChanged = true;
00395 break;
00396 }
00397 }
00398 else
00399 {
00400 G4int counter = -1;
00401 for( np=0; np<numSec/3 && ran>=excs; ++np )
00402 {
00403 for( nneg=std::max(0,np-1); nneg<=(np+2) && ran>=excs; ++nneg )
00404 {
00405 for( nz=0; nz<numSec/3 && ran>=excs; ++nz )
00406 {
00407 if( ++counter < numMul )
00408 {
00409 nt = np+nneg+nz;
00410 if( nt>0 && nt<=numSec )
00411 {
00412 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00413 dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
00414 if( std::fabs(dum) < 1.0 )
00415 {
00416 if( test >= 1.0e-10 )excs += dum*test;
00417 }
00418 else
00419 excs += dum*test;
00420 }
00421 }
00422 }
00423 }
00424 }
00425 if( ran >= excs )
00426 {
00427 quasiElastic = true;
00428 return;
00429 }
00430 np--; nneg--; nz--;
00431 G4int ncht = std::min( 4, std::max( 1, np-nneg+3 ) );
00432 switch( ncht )
00433 {
00434 case 1:
00435 currentParticle.SetDefinitionAndUpdateE( anAntiSigmaMinus );
00436 incidentHasChanged = true;
00437 targetParticle.SetDefinitionAndUpdateE( aProton );
00438 targetHasChanged = true;
00439 break;
00440 case 2:
00441 if( G4UniformRand() < 0.5 )
00442 {
00443 if( G4UniformRand() < 0.5 )
00444 {
00445 currentParticle.SetDefinitionAndUpdateE( anAntiSigmaZero );
00446 incidentHasChanged = true;
00447 targetParticle.SetDefinitionAndUpdateE( aProton );
00448 targetHasChanged = true;
00449 }
00450 else
00451 {
00452 currentParticle.SetDefinitionAndUpdateE( anAntiSigmaMinus );
00453 incidentHasChanged = true;
00454 }
00455 }
00456 else
00457 {
00458 if( G4UniformRand() < 0.5 )
00459 {
00460 targetParticle.SetDefinitionAndUpdateE( aProton );
00461 targetHasChanged = true;
00462 }
00463 else
00464 {
00465 currentParticle.SetDefinitionAndUpdateE( anAntiSigmaMinus );
00466 incidentHasChanged = true;
00467 }
00468 }
00469 break;
00470 case 3:
00471 if( G4UniformRand() < 0.5 )
00472 {
00473 if( G4UniformRand() < 0.5 )
00474 {
00475 currentParticle.SetDefinitionAndUpdateE( anAntiSigmaZero );
00476 incidentHasChanged = true;
00477 }
00478 else
00479 {
00480 currentParticle.SetDefinitionAndUpdateE( anAntiSigmaPlus );
00481 incidentHasChanged = true;
00482 targetParticle.SetDefinitionAndUpdateE( aProton );
00483 targetHasChanged = true;
00484 }
00485 }
00486 else
00487 {
00488 if( G4UniformRand() >= 0.5 )
00489 {
00490 currentParticle.SetDefinitionAndUpdateE( anAntiSigmaPlus );
00491 incidentHasChanged = true;
00492 targetParticle.SetDefinitionAndUpdateE( aProton );
00493 targetHasChanged = true;
00494 }
00495 }
00496 break;
00497 default:
00498 currentParticle.SetDefinitionAndUpdateE( anAntiSigmaPlus );
00499 incidentHasChanged = true;
00500 break;
00501 }
00502 }
00503 }
00504 else
00505 {
00506 if( centerofmassEnergy <= aPiPlus->GetPDGMass()/MeV+aKaonPlus->GetPDGMass()/MeV )
00507 {
00508 quasiElastic = true;
00509 return;
00510 }
00511
00512
00513
00514 G4double n, anpn;
00515 GetNormalizationConstant( -centerofmassEnergy, n, anpn );
00516 G4double ran = G4UniformRand();
00517 G4double dum, excs = 0.0;
00518 if( targetParticle.GetDefinition() == aProton )
00519 {
00520 G4int counter = -1;
00521 for( np=1; np<numSec/3 && ran>=excs; ++np )
00522 {
00523 nneg = np-1;
00524 for( nz=0; nz<numSec/3 && ran>=excs; ++nz )
00525 {
00526 if( ++counter < numMulA )
00527 {
00528 nt = np+nneg+nz;
00529 if( nt>1 && nt<=numSec )
00530 {
00531 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00532 dum = (pi/anpn)*nt*protmulA[counter]*protnormA[nt-1]/(2.0*n*n);
00533 if( std::fabs(dum) < 1.0 )
00534 {
00535 if( test >= 1.0e-10 )excs += dum*test;
00536 }
00537 else
00538 excs += dum*test;
00539 }
00540 }
00541 }
00542 }
00543 }
00544 else
00545 {
00546 G4int counter = -1;
00547 for( np=0; np<numSec/3 && ran>=excs; ++np )
00548 {
00549 nneg = np;
00550 for( nz=0; nz<numSec/3 && ran>=excs; ++nz )
00551 {
00552 if( ++counter < numMulA )
00553 {
00554 nt = np+nneg+nz;
00555 if( nt>1 && nt<=numSec )
00556 {
00557 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00558 dum = (pi/anpn)*nt*neutmulA[counter]*neutnormA[nt-1]/(2.0*n*n);
00559 if( std::fabs(dum) < 1.0 )
00560 {
00561 if( test >= 1.0e-10 )excs += dum*test;
00562 }
00563 else
00564 excs += dum*test;
00565 }
00566 }
00567 }
00568 }
00569 }
00570 if( ran >= excs )
00571 {
00572 quasiElastic = true;
00573 return;
00574 }
00575 np--; nz--;
00576 currentParticle.SetMass( 0.0 );
00577 targetParticle.SetMass( 0.0 );
00578 }
00579 SetUpPions( np, nneg, nz, vec, vecLen );
00580 if( currentParticle.GetMass() == 0.0 )
00581 {
00582 if( nz == 0 )
00583 {
00584 if( nneg > 0 )
00585 {
00586 for( G4int i=0; i<vecLen; ++i )
00587 {
00588 if( vec[i]->GetDefinition() == aPiMinus )
00589 {
00590 vec[i]->SetDefinitionAndUpdateE( aKaonMinus );
00591 break;
00592 }
00593 }
00594 }
00595 }
00596 else
00597 {
00598 if( nneg == 0 )
00599 {
00600 for( G4int i=0; i<vecLen; ++i )
00601 {
00602 if( vec[i]->GetDefinition() == aPiZero )
00603 {
00604 vec[i]->SetDefinitionAndUpdateE( aKaonZL );
00605 break;
00606 }
00607 }
00608 }
00609 else
00610 {
00611 if( G4UniformRand() < 0.5 )
00612 {
00613 if( nneg > 0 )
00614 {
00615 for( G4int i=0; i<vecLen; ++i )
00616 {
00617 if( vec[i]->GetDefinition() == aPiMinus )
00618 {
00619 vec[i]->SetDefinitionAndUpdateE( aKaonMinus );
00620 break;
00621 }
00622 }
00623 }
00624 }
00625 else
00626 {
00627 for( G4int i=0; i<vecLen; ++i )
00628 {
00629 if( vec[i]->GetDefinition() == aPiZero )
00630 {
00631 vec[i]->SetDefinitionAndUpdateE( aKaonZL );
00632 break;
00633 }
00634 }
00635 }
00636 }
00637 }
00638 }
00639 return;
00640 }
00641
00642
00643