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