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 "G4RPGAntiKZeroInelastic.hh"
00030 #include "Randomize.hh"
00031 #include "G4PhysicalConstants.hh"
00032 #include "G4SystemOfUnits.hh"
00033 #include "G4HadReentrentException.hh"
00034
00035 G4HadFinalState *
00036 G4RPGAntiKZeroInelastic::ApplyYourself( const G4HadProjectile &aTrack,
00037 G4Nucleus &targetNucleus )
00038 {
00039 const G4HadProjectile *originalIncident = &aTrack;
00040
00041
00042
00043 G4DynamicParticle *originalTarget = targetNucleus.ReturnTargetParticle();
00044
00045 if( verboseLevel > 1 )
00046 {
00047 const G4Material *targetMaterial = aTrack.GetMaterial();
00048 G4cout << "G4RPGAntiKZeroInelastic::ApplyYourself called" << G4endl;
00049 G4cout << "kinetic energy = " << originalIncident->GetKineticEnergy()/MeV << "MeV, ";
00050 G4cout << "target material = " << targetMaterial->GetName() << ", ";
00051 G4cout << "target particle = " << originalTarget->GetDefinition()->GetParticleName()
00052 << G4endl;
00053 }
00054
00055
00056
00057
00058 G4double ek = originalIncident->GetKineticEnergy()/MeV;
00059 G4double amas = originalIncident->GetDefinition()->GetPDGMass()/MeV;
00060 G4ReactionProduct modifiedOriginal;
00061 modifiedOriginal = *originalIncident;
00062
00063 G4double tkin = targetNucleus.Cinema( ek );
00064 ek += tkin;
00065 modifiedOriginal.SetKineticEnergy( ek*MeV );
00066 G4double et = ek + amas;
00067 G4double p = std::sqrt( std::abs((et-amas)*(et+amas)) );
00068 G4double pp = modifiedOriginal.GetMomentum().mag()/MeV;
00069 if( pp > 0.0 )
00070 {
00071 G4ThreeVector momentum = modifiedOriginal.GetMomentum();
00072 modifiedOriginal.SetMomentum( momentum * (p/pp) );
00073 }
00074
00075
00076
00077 tkin = targetNucleus.EvaporationEffects( ek );
00078 ek -= tkin;
00079 modifiedOriginal.SetKineticEnergy( ek*MeV );
00080 et = ek + amas;
00081 p = std::sqrt( std::abs((et-amas)*(et+amas)) );
00082 pp = modifiedOriginal.GetMomentum().mag()/MeV;
00083 if( pp > 0.0 )
00084 {
00085 G4ThreeVector momentum = modifiedOriginal.GetMomentum();
00086 modifiedOriginal.SetMomentum( momentum * (p/pp) );
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 if( currentParticle.GetKineticEnergy()/MeV > cutOff )
00102 Cascade( vec, vecLen,
00103 originalIncident, currentParticle, targetParticle,
00104 incidentHasChanged, targetHasChanged, quasiElastic );
00105
00106 try
00107 {
00108 CalculateMomenta( vec, vecLen,
00109 originalIncident, originalTarget, modifiedOriginal,
00110 targetNucleus, currentParticle, targetParticle,
00111 incidentHasChanged, targetHasChanged, quasiElastic );
00112 }
00113 catch(G4HadReentrentException aR)
00114 {
00115 aR.Report(G4cout);
00116 throw G4HadReentrentException(__FILE__, __LINE__, "Bailing out");
00117 }
00118 SetUpChange( vec, vecLen,
00119 currentParticle, targetParticle,
00120 incidentHasChanged );
00121
00122 delete originalTarget;
00123 return &theParticleChange;
00124 }
00125
00126 void G4RPGAntiKZeroInelastic::Cascade(
00127 G4FastVector<G4ReactionProduct,GHADLISTSIZE> &vec,
00128 G4int& vecLen,
00129 const G4HadProjectile *originalIncident,
00130 G4ReactionProduct ¤tParticle,
00131 G4ReactionProduct &targetParticle,
00132 G4bool &incidentHasChanged,
00133 G4bool &targetHasChanged,
00134 G4bool &quasiElastic )
00135 {
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146 const G4double mOriginal = originalIncident->GetDefinition()->GetPDGMass()/MeV;
00147 const G4double etOriginal = originalIncident->Get4Momentum().e()/MeV;
00148 const G4double pOriginal = originalIncident->Get4Momentum().vect().mag()/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 numSec = 60;
00158 static G4double protmul[numMul], protnorm[numSec];
00159 static G4double neutmul[numMul], neutnorm[numSec];
00160
00161
00162
00163 G4int counter, nt=0, np=0, nneg=0, nz=0;
00164 const G4double c = 1.25;
00165 const G4double b[] = { 0.7, 0.7 };
00166 if( first )
00167 {
00168 first = false;
00169 G4int i;
00170 for( i=0; i<numMul; ++i )protmul[i] = 0.0;
00171 for( i=0; i<numSec; ++i )protnorm[i] = 0.0;
00172 counter = -1;
00173 for( np=0; np<numSec/3; ++np )
00174 {
00175 for( nneg=std::max(0,np-2); nneg<=np; ++nneg )
00176 {
00177 for( nz=0; nz<numSec/3; ++nz )
00178 {
00179 if( ++counter < numMul )
00180 {
00181 nt = np+nneg+nz;
00182 if( nt>0 && nt<=numSec )
00183 {
00184 protmul[counter] = Pmltpc(np,nneg,nz,nt,b[0],c);
00185 protnorm[nt-1] += protmul[counter];
00186 }
00187 }
00188 }
00189 }
00190 }
00191 for( i=0; i<numMul; ++i )neutmul[i] = 0.0;
00192 for( i=0; i<numSec; ++i )neutnorm[i] = 0.0;
00193 counter = -1;
00194 for( np=0; np<(numSec/3); ++np )
00195 {
00196 for( nneg=std::max(0,np-1); nneg<=(np+1); ++nneg )
00197 {
00198 for( nz=0; nz<numSec/3; ++nz )
00199 {
00200 if( ++counter < numMul )
00201 {
00202 nt = np+nneg+nz;
00203 if( nt>0 && nt<=numSec )
00204 {
00205 neutmul[counter] = Pmltpc(np,nneg,nz,nt,b[1],c);
00206 neutnorm[nt-1] += neutmul[counter];
00207 }
00208 }
00209 }
00210 }
00211 }
00212 for( i=0; i<numSec; ++i )
00213 {
00214 if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
00215 if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
00216 }
00217 }
00218
00219 const G4double expxu = 82.;
00220 const G4double expxl = -expxu;
00221 G4ParticleDefinition *aKaonMinus = G4KaonMinus::KaonMinus();
00222 G4ParticleDefinition *aKaonZS = G4KaonZeroShort::KaonZeroShort();
00223 G4ParticleDefinition *aKaonZL = G4KaonZeroLong::KaonZeroLong();
00224 G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
00225 G4ParticleDefinition *aProton = G4Proton::Proton();
00226 G4ParticleDefinition *aPiPlus = G4PionPlus::PionPlus();
00227 G4ParticleDefinition *aPiMinus = G4PionMinus::PionMinus();
00228 G4ParticleDefinition *aPiZero = G4PionZero::PionZero();
00229 G4ParticleDefinition *aLambda = G4Lambda::Lambda();
00230 G4ParticleDefinition *aSigmaPlus = G4SigmaPlus::SigmaPlus();
00231 G4ParticleDefinition *aSigmaMinus = G4SigmaMinus::SigmaMinus();
00232 G4ParticleDefinition *aSigmaZero = G4SigmaZero::SigmaZero();
00233 const G4double cech[] = {1.,1.,1.,0.70,0.60,0.55,0.35,0.25,0.18,0.15};
00234 G4int iplab = G4int(std::min( 9.0, 5.0*pOriginal*MeV/GeV ));
00235
00236 if ((pOriginal*MeV/GeV <= 2.0) && (G4UniformRand() < cech[iplab]) ) {
00237 np = nneg = nz = nt = 0;
00238 iplab = G4int(std::min( 19.0, pOriginal*MeV/GeV*10.0 ));
00239 const G4double cnk0[] = {0.17,0.18,0.17,0.24,0.26,0.20,0.22,0.21,0.34,0.45,
00240 0.58,0.55,0.36,0.29,0.29,0.32,0.32,0.33,0.33,0.33};
00241 if (G4UniformRand() > cnk0[iplab] ) {
00242 G4double ran = G4UniformRand();
00243 if (ran < 0.25) {
00244 if (targetParticle.GetDefinition() == aNeutron) {
00245 currentParticle.SetDefinitionAndUpdateE( aPiMinus );
00246 targetParticle.SetDefinitionAndUpdateE( aSigmaPlus );
00247 incidentHasChanged = true;
00248 targetHasChanged = true;
00249 }
00250 } else if( ran < 0.50 ) {
00251 if( targetParticle.GetDefinition() == aNeutron )
00252 currentParticle.SetDefinitionAndUpdateE( aPiZero );
00253 else
00254 currentParticle.SetDefinitionAndUpdateE( aPiPlus );
00255 targetParticle.SetDefinitionAndUpdateE( aSigmaZero );
00256 incidentHasChanged = true;
00257 targetHasChanged = true;
00258 } else if( ran < 0.75 ) {
00259 if( targetParticle.GetDefinition() == aNeutron )
00260 {
00261 currentParticle.SetDefinitionAndUpdateE( aPiPlus );
00262 targetParticle.SetDefinitionAndUpdateE( aSigmaMinus );
00263 incidentHasChanged = true;
00264 targetHasChanged = true;
00265 }
00266 } else {
00267 if( targetParticle.GetDefinition() == aNeutron )
00268 currentParticle.SetDefinitionAndUpdateE( aPiZero );
00269 else
00270 currentParticle.SetDefinitionAndUpdateE( aPiPlus );
00271 targetParticle.SetDefinitionAndUpdateE( aLambda );
00272 incidentHasChanged = true;
00273 targetHasChanged = true;
00274 }
00275 } else {
00276 quasiElastic = true;
00277 if (targetParticle.GetDefinition() == aNeutron) {
00278 currentParticle.SetDefinitionAndUpdateE( aKaonMinus );
00279 targetParticle.SetDefinitionAndUpdateE( aProton );
00280 incidentHasChanged = true;
00281 targetHasChanged = true;
00282 }
00283 }
00284 } else {
00285 if (availableEnergy < aPiPlus->GetPDGMass()/MeV) {
00286 quasiElastic = true;
00287 return;
00288 }
00289 G4double n, anpn;
00290 GetNormalizationConstant( availableEnergy, n, anpn );
00291 G4double ran = G4UniformRand();
00292 G4double dum, test, excs = 0.0;
00293 if (targetParticle.GetDefinition() == aProton) {
00294 counter = -1;
00295 for( np=0; (np<numSec/3) && (ran>=excs); ++np )
00296 {
00297 for( nneg=std::max(0,np-2); nneg<=np && ran>=excs; ++nneg )
00298 {
00299 for( nz=0; nz<numSec/3 && ran>=excs; ++nz )
00300 {
00301 if( ++counter < numMul )
00302 {
00303 nt = np+nneg+nz;
00304 if( nt>0 && nt<=numSec )
00305 {
00306 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00307 dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
00308 if( std::fabs(dum) < 1.0 )
00309 {
00310 if( test >= 1.0e-10 )excs += dum*test;
00311 }
00312 else
00313 excs += dum*test;
00314 }
00315 }
00316 }
00317 }
00318 }
00319 if( ran >= excs )
00320 {
00321 quasiElastic = true;
00322 return;
00323 }
00324 np--; nneg--; nz--;
00325 switch( np-nneg )
00326 {
00327 case 1:
00328 if( G4UniformRand() < 0.5 )
00329 {
00330 currentParticle.SetDefinitionAndUpdateE( aKaonMinus );
00331 incidentHasChanged = true;
00332 }
00333 else
00334 {
00335 targetParticle.SetDefinitionAndUpdateE( aNeutron );
00336 targetHasChanged = true;
00337 }
00338 case 0:
00339 break;
00340 default:
00341 currentParticle.SetDefinitionAndUpdateE( aKaonMinus );
00342 targetParticle.SetDefinitionAndUpdateE( aNeutron );
00343 incidentHasChanged = true;
00344 targetHasChanged = true;
00345 break;
00346 }
00347 }
00348 else
00349 {
00350 counter = -1;
00351 for( np=0; np<numSec/3 && ran>=excs; ++np )
00352 {
00353 for( nneg=std::max(0,np-1); nneg<=(np+1) && ran>=excs; ++nneg )
00354 {
00355 for( nz=0; nz<numSec/3 && ran>=excs; ++nz )
00356 {
00357 if( ++counter < numMul )
00358 {
00359 nt = np+nneg+nz;
00360 if( nt>0 && nt<=numSec )
00361 {
00362 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
00363 dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
00364 if( std::fabs(dum) < 1.0 )
00365 {
00366 if( test >= 1.0e-10 )excs += dum*test;
00367 }
00368 else
00369 excs += dum*test;
00370 }
00371 }
00372 }
00373 }
00374 }
00375 if( ran >= excs )
00376 {
00377 quasiElastic = true;
00378 return;
00379 }
00380 np--; nneg--; nz--;
00381 switch( np-nneg )
00382 {
00383 case 0:
00384 currentParticle.SetDefinitionAndUpdateE( aKaonMinus );
00385 targetParticle.SetDefinitionAndUpdateE( aProton );
00386 incidentHasChanged = true;
00387 targetHasChanged = true;
00388 break;
00389 case 1:
00390 currentParticle.SetDefinitionAndUpdateE( aKaonMinus );
00391 incidentHasChanged = true;
00392 break;
00393 default:
00394 targetParticle.SetDefinitionAndUpdateE( aProton );
00395 targetHasChanged = true;
00396 break;
00397 }
00398 }
00399 if( G4UniformRand() >= 0.5 )
00400 {
00401 if( currentParticle.GetDefinition() == aKaonMinus &&
00402 targetParticle.GetDefinition() == aNeutron )
00403 {
00404 ran = G4UniformRand();
00405 if( ran < 0.68 )
00406 {
00407 currentParticle.SetDefinitionAndUpdateE( aPiMinus );
00408 targetParticle.SetDefinitionAndUpdateE( aLambda );
00409 }
00410 else if( ran < 0.84 )
00411 {
00412 currentParticle.SetDefinitionAndUpdateE( aPiMinus );
00413 targetParticle.SetDefinitionAndUpdateE( aSigmaZero );
00414 }
00415 else
00416 {
00417 currentParticle.SetDefinitionAndUpdateE( aPiZero );
00418 targetParticle.SetDefinitionAndUpdateE( aSigmaMinus );
00419 }
00420 }
00421 else if( (currentParticle.GetDefinition() == aKaonZS ||
00422 currentParticle.GetDefinition() == aKaonZL ) &&
00423 targetParticle.GetDefinition() == aProton )
00424 {
00425 ran = G4UniformRand();
00426 if( ran < 0.68 )
00427 {
00428 currentParticle.SetDefinitionAndUpdateE( aPiPlus );
00429 targetParticle.SetDefinitionAndUpdateE( aLambda );
00430 }
00431 else if( ran < 0.84 )
00432 {
00433 currentParticle.SetDefinitionAndUpdateE( aPiZero );
00434 targetParticle.SetDefinitionAndUpdateE( aSigmaPlus );
00435 }
00436 else
00437 {
00438 currentParticle.SetDefinitionAndUpdateE( aPiPlus );
00439 targetParticle.SetDefinitionAndUpdateE( aSigmaZero );
00440 }
00441 }
00442 else
00443 {
00444 ran = G4UniformRand();
00445 if( ran < 0.67 )
00446 {
00447 currentParticle.SetDefinitionAndUpdateE( aPiZero );
00448 targetParticle.SetDefinitionAndUpdateE( aLambda );
00449 }
00450 else if( ran < 0.78 )
00451 {
00452 currentParticle.SetDefinitionAndUpdateE( aPiMinus );
00453 targetParticle.SetDefinitionAndUpdateE( aSigmaPlus );
00454 }
00455 else if( ran < 0.89 )
00456 {
00457 currentParticle.SetDefinitionAndUpdateE( aPiZero );
00458 targetParticle.SetDefinitionAndUpdateE( aSigmaZero );
00459 }
00460 else
00461 {
00462 currentParticle.SetDefinitionAndUpdateE( aPiPlus );
00463 targetParticle.SetDefinitionAndUpdateE( aSigmaMinus );
00464 }
00465 }
00466 incidentHasChanged = true;
00467 targetHasChanged = true;
00468 }
00469 }
00470 if( currentParticle.GetDefinition() == aKaonZL )
00471 {
00472 if( G4UniformRand() >= 0.5 )
00473 {
00474 currentParticle.SetDefinitionAndUpdateE( aKaonZS );
00475 incidentHasChanged = true;
00476 }
00477 }
00478 if( targetParticle.GetDefinition() == aKaonZL )
00479 {
00480 if( G4UniformRand() >= 0.5 )
00481 {
00482 targetParticle.SetDefinitionAndUpdateE( aKaonZS );
00483 targetHasChanged = true;
00484 }
00485 }
00486 SetUpPions( np, nneg, nz, vec, vecLen );
00487 return;
00488 }
00489
00490
00491