Changeset 1315 for trunk/source/processes/hadronic/models/high_energy/src
- Timestamp:
- Jun 18, 2010, 11:42:07 AM (14 years ago)
- Location:
- trunk/source/processes/hadronic/models/high_energy/src
- Files:
-
- 4 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/source/processes/hadronic/models/high_energy/src/G4HEAntiKaonZeroInelastic.cc
r1228 r1315 25 25 // 26 26 // 27 // $Id: G4HEAntiKaonZeroInelastic.cc,v 1.1 5 2008/03/17 20:49:17dennis Exp $28 // GEANT4 tag $Name: geant4-09-0 3$27 // $Id: G4HEAntiKaonZeroInelastic.cc,v 1.16 2010/02/09 22:02:04 dennis Exp $ 28 // GEANT4 tag $Name: geant4-09-04-beta-cand-01 $ 29 29 // 30 30 // … … 179 179 } 180 180 181 if (!successful) 182 {183 G4cout << "GHEInelasticInteraction::ApplyYourself fails to produce final state particles"<< G4endl;184 } 181 if (!successful) 182 G4cout << "GHEInelasticInteraction::ApplyYourself fails to produce final state particles" 183 << G4endl; 184 185 185 FillParticleChange(pv, vecLength); 186 186 delete [] pv; -
trunk/source/processes/hadronic/models/high_energy/src/G4HEKaonZeroInelastic.cc
r1228 r1315 25 25 // 26 26 // 27 // $Id: G4HEKaonZeroInelastic.cc,v 1.1 5 2008/03/17 20:49:17dennis Exp $28 // GEANT4 tag $Name: geant4-09-0 3$27 // $Id: G4HEKaonZeroInelastic.cc,v 1.16 2010/02/09 22:02:20 dennis Exp $ 28 // GEANT4 tag $Name: geant4-09-04-beta-cand-01 $ 29 29 // 30 30 // … … 182 182 } 183 183 184 if (!successful) 185 { 186 G4cout << "GHEInelasticInteraction::ApplyYourself fails to produce final state particles" << G4endl; 187 } 188 FillParticleChange(pv, vecLength); 189 delete [] pv; 190 theParticleChange.SetStatusChange(stopAndKill); 191 return & theParticleChange; 184 if (!successful) 185 G4cout << "GHEInelasticInteraction::ApplyYourself fails to produce final state particles" 186 << G4endl; 187 188 FillParticleChange(pv, vecLength); 189 190 delete [] pv; 191 theParticleChange.SetStatusChange(stopAndKill); 192 return & theParticleChange; 192 193 } 193 194 -
trunk/source/processes/hadronic/models/high_energy/src/G4HEKaonZeroLongInelastic.cc
r1228 r1315 24 24 // ******************************************************************** 25 25 // 26 // 27 // $Id: G4HEKaonZeroLongInelastic.cc,v 1.10 2006/06/29 20:30:22 gunter Exp $ 28 // GEANT4 tag $Name: geant4-09-03 $ 26 // $Id: G4HEKaonZeroLongInelastic.cc,v 1.11 2010/02/09 21:59:10 dennis Exp $ 27 // GEANT4 tag $Name: geant4-09-04-beta-cand-01 $ 29 28 // 30 29 // … … 38 37 // and the low energy two-body model. Not included are the low energy stuff like 39 38 // nuclear reactions, nuclear fission without any cascading and all processes for 40 // particles at rest. 41 // First work done by J.L.Chuma and F.W.Jones, TRIUMF, June 96. 42 // H. Fesefeldt, RWTH-Aachen, 23-October-1996 43 // Last modified: 29-July-1998 39 // particles at rest. 40 // 41 // New version by D.H. Wright (SLAC) to fix seg fault in old version 42 // 26 January 2010 43 44 44 45 45 #include "G4HEKaonZeroLongInelastic.hh" 46 46 47 G4HadFinalState * G4HEKaonZeroLongInelastic:: 48 ApplyYourself( const G4HadProjectile &aTrack, G4Nucleus &targetNucleus ) 49 { 50 G4HEKaonZeroInelastic theKaonZeroInelastic; 51 G4HEAntiKaonZeroInelastic theAntiKaonZeroInelastic; 52 theKaonZeroInelastic.SetVerboseLevel(verboseLevel); 53 theAntiKaonZeroInelastic.SetVerboseLevel(verboseLevel); 47 G4HadFinalState* G4HEKaonZeroLongInelastic:: 48 ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus &targetNucleus) 49 { 50 G4HEVector * pv = new G4HEVector[MAXPART]; 51 const G4HadProjectile *aParticle = &aTrack; 52 const G4double atomicWeight = targetNucleus.GetN(); 53 const G4double atomicNumber = targetNucleus.GetZ(); 54 G4HEVector incidentParticle(aParticle); 55 56 G4int incidentCode = incidentParticle.getCode(); 57 G4double incidentMass = incidentParticle.getMass(); 58 G4double incidentTotalEnergy = incidentParticle.getEnergy(); 59 G4double incidentTotalMomentum = incidentParticle.getTotalMomentum(); 60 G4double incidentKineticEnergy = incidentTotalEnergy - incidentMass; 61 62 if(incidentKineticEnergy < 1) 63 G4cout << "GHEKaonZeroLongInelastic: incident energy < 1 GeV " << G4endl; 64 65 if(verboseLevel > 1) { 66 G4cout << "G4HEKaonZeroLongInelastic::ApplyYourself" << G4endl; 67 G4cout << "incident particle " << incidentParticle.getName() 68 << "mass " << incidentMass 69 << "kinetic energy " << incidentKineticEnergy 70 << G4endl; 71 G4cout << "target material with (A,Z) = (" 72 << atomicWeight << "," << atomicNumber << ")" << G4endl; 73 } 54 74 55 if(G4UniformRand() < 0.50) 56 { 57 return theKaonZeroInelastic.ApplyYourself(aTrack, targetNucleus); 58 } 75 G4double inelasticity = NuclearInelasticity(incidentKineticEnergy, 76 atomicWeight, atomicNumber); 77 if(verboseLevel > 1) 78 G4cout << "nuclear inelasticity = " << inelasticity << G4endl; 79 80 incidentKineticEnergy -= inelasticity; 81 82 G4double excitationEnergyGNP = 0.; 83 G4double excitationEnergyDTA = 0.; 84 85 G4double excitation = NuclearExcitation(incidentKineticEnergy, 86 atomicWeight, atomicNumber, 87 excitationEnergyGNP, 88 excitationEnergyDTA); 89 if(verboseLevel > 1) 90 G4cout << "nuclear excitation = " << excitation << excitationEnergyGNP 91 << excitationEnergyDTA << G4endl; 92 93 incidentKineticEnergy -= excitation; 94 incidentTotalEnergy = incidentKineticEnergy + incidentMass; 95 incidentTotalMomentum = std::sqrt( (incidentTotalEnergy-incidentMass) 96 *(incidentTotalEnergy+incidentMass)); 97 98 99 G4HEVector targetParticle; 100 if(G4UniformRand() < atomicNumber/atomicWeight) { 101 targetParticle.setDefinition("Proton"); 102 } else { 103 targetParticle.setDefinition("Neutron"); 104 } 105 106 G4double targetMass = targetParticle.getMass(); 107 G4double centerOfMassEnergy = std::sqrt( incidentMass*incidentMass + targetMass*targetMass 108 + 2.0*targetMass*incidentTotalEnergy); 109 G4double availableEnergy = centerOfMassEnergy - targetMass - incidentMass; 110 // this was the meaning of inElastic in the 111 // original Gheisha stand-alone version. 112 // G4bool inElastic = InElasticCrossSectionInFirstInt 113 // (availableEnergy, incidentCode, incidentTotalMomentum); 114 // for unknown reasons, it has been replaced by the following code in Geant??? 115 116 G4bool inElastic = true; 117 // if (G4UniformRand() < elasticCrossSection/totalCrossSection) inElastic = false; 118 119 vecLength = 0; 120 121 if(verboseLevel > 1) 122 G4cout << "ApplyYourself: CallFirstIntInCascade for particle " 123 << incidentCode << G4endl; 124 125 G4bool successful = false; 126 127 if(inElastic || (!inElastic && atomicWeight < 1.5)) { 128 129 // Split K0L into K0 and K0bar 130 if (G4UniformRand() < 0.5) 131 FirstIntInCasAntiKaonZero(inElastic, availableEnergy, pv, vecLength, 132 incidentParticle, targetParticle ); 59 133 else 60 { 61 return theAntiKaonZeroInelastic.ApplyYourself(aTrack, targetNucleus); 62 } 134 FirstIntInCasKaonZero(inElastic, availableEnergy, pv, vecLength, 135 incidentParticle, targetParticle, atomicWeight ); 136 137 // Do nuclear interaction with either K0 or K0bar 138 if ((vecLength > 0) && (availableEnergy > 1.)) 139 StrangeParticlePairProduction(availableEnergy, centerOfMassEnergy, 140 pv, vecLength, 141 incidentParticle, targetParticle); 142 HighEnergyCascading(successful, pv, vecLength, 143 excitationEnergyGNP, excitationEnergyDTA, 144 incidentParticle, targetParticle, 145 atomicWeight, atomicNumber); 146 if (!successful) 147 HighEnergyClusterProduction(successful, pv, vecLength, 148 excitationEnergyGNP, excitationEnergyDTA, 149 incidentParticle, targetParticle, 150 atomicWeight, atomicNumber); 151 if (!successful) 152 MediumEnergyCascading(successful, pv, vecLength, 153 excitationEnergyGNP, excitationEnergyDTA, 154 incidentParticle, targetParticle, 155 atomicWeight, atomicNumber); 156 157 if (!successful) 158 MediumEnergyClusterProduction(successful, pv, vecLength, 159 excitationEnergyGNP, excitationEnergyDTA, 160 incidentParticle, targetParticle, 161 atomicWeight, atomicNumber); 162 if (!successful) 163 QuasiElasticScattering(successful, pv, vecLength, 164 excitationEnergyGNP, excitationEnergyDTA, 165 incidentParticle, targetParticle, 166 atomicWeight, atomicNumber); 167 } 168 169 if (!successful) 170 ElasticScattering(successful, pv, vecLength, 171 incidentParticle, 172 atomicWeight, atomicNumber); 173 174 if (!successful) 175 G4cout << "GHEInelasticInteraction::ApplyYourself fails to produce final state particles" 176 << G4endl; 177 178 // Check for K0, K0bar and change particle types to K0L, K0S if necessary 179 G4int kcode; 180 for (G4int i = 0; i < vecLength; i++) { 181 kcode = pv[i].getCode(); 182 if (kcode == KaonZero.getCode() || kcode == AntiKaonZero.getCode()) { 183 if (G4UniformRand() < 0.5) 184 pv[i] = KaonZeroShort; 185 else 186 pv[i] = KaonZeroLong; 187 } 63 188 } 64 65 66 67 68 69 70 71 72 189 190 // ................ 191 192 FillParticleChange(pv, vecLength); 193 delete [] pv; 194 theParticleChange.SetStatusChange(stopAndKill); 195 return & theParticleChange; 196 } 197 198 199 void 200 G4HEKaonZeroLongInelastic::FirstIntInCasKaonZero(G4bool &inElastic, 201 const G4double availableEnergy, 202 G4HEVector pv[], 203 G4int &vecLen, 204 G4HEVector incidentParticle, 205 G4HEVector targetParticle, 206 const G4double atomicWeight) 207 208 // Kaon0 undergoes interaction with nucleon within a nucleus. Check if it is 209 // energetically possible to produce pions/kaons. In not, assume nuclear excitation 210 // occurs and input particle is degraded in energy. No other particles are produced. 211 // If reaction is possible, find the correct number of pions/protons/neutrons 212 // produced using an interpolation to multiplicity data. Replace some pions or 213 // protons/neutrons by kaons or strange baryons according to the average 214 // multiplicity per inelastic reaction. 215 216 { 217 static const G4double expxu = std::log(MAXFLOAT); // upper bound for arg. of exp 218 static const G4double expxl = -expxu; // lower bound for arg. of exp 219 220 static const G4double protb = 0.7; 221 static const G4double neutb = 0.7; 222 static const G4double c = 1.25; 223 224 static const G4int numMul = 1200; 225 static const G4int numSec = 60; 226 227 G4int neutronCode = Neutron.getCode(); 228 G4int protonCode = Proton.getCode(); 229 230 G4int targetCode = targetParticle.getCode(); 231 G4double incidentTotalMomentum = incidentParticle.getTotalMomentum(); 232 233 static G4bool first = true; 234 static G4double protmul[numMul], protnorm[numSec]; // proton constants 235 static G4double neutmul[numMul], neutnorm[numSec]; // neutron constants 236 237 // misc. local variables 238 // np = number of pi+, nm = number of pi-, nz = number of pi0 239 240 G4int i, counter, nt, np, nm, nz; 241 242 if (first) { 243 // compute normalization constants, this will only be done once 244 first = false; 245 for( i=0; i<numMul; i++ )protmul[i] = 0.0; 246 for( i=0; i<numSec; i++ )protnorm[i] = 0.0; 247 counter = -1; 248 for (np=0; np<(numSec/3); np++) { 249 for (nm=std::max(0,np-1); nm<=(np+1); nm++) { 250 for (nz=0; nz<numSec/3; nz++) { 251 if (++counter < numMul) { 252 nt = np+nm+nz; 253 if( (nt>0) && (nt<=numSec) ) { 254 protmul[counter] = pmltpc(np,nm,nz,nt,protb,c) ; 255 protnorm[nt-1] += protmul[counter]; 256 } 257 } 258 } 259 } 260 } 261 262 for( i=0; i<numMul; i++ )neutmul[i] = 0.0; 263 for( i=0; i<numSec; i++ )neutnorm[i] = 0.0; 264 counter = -1; 265 for (np=0; np<numSec/3; np++) { 266 for (nm=np; nm<=(np+2); nm++) { 267 for (nz=0; nz<numSec/3; nz++) { 268 if (++counter < numMul) { 269 nt = np+nm+nz; 270 if( (nt>0) && (nt<=numSec) ) { 271 neutmul[counter] = pmltpc(np,nm,nz,nt,neutb,c); 272 neutnorm[nt-1] += neutmul[counter]; 273 } 274 } 275 } 276 } 277 } 278 279 for (i=0; i<numSec; i++) { 280 if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i]; 281 if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i]; 282 } 283 } // end of initialization 284 285 286 // Initialize the first two particles 287 // the same as beam and target 288 pv[0] = incidentParticle; 289 pv[1] = targetParticle; 290 vecLen = 2; 291 292 if( !inElastic ) { 293 // quasi-elastic scattering, no pions produced 294 if( targetCode == protonCode) { 295 G4double cech[] = {0.33,0.27,0.29,0.31,0.27,0.18,0.13,0.10,0.09,0.07}; 296 G4int iplab = G4int( std::min( 9.0, incidentTotalMomentum*5. ) ); 297 if( G4UniformRand() < cech[iplab]/std::pow(atomicWeight,0.42)) { 298 // charge exchange K+ n -> K0 p 299 pv[0] = KaonPlus; 300 pv[1] = Neutron; 301 } 302 } 303 return; 304 } else if (availableEnergy <= PionPlus.getMass()) { 305 return; 306 } 307 308 // Inelastic scattering 309 310 np = 0, nm = 0, nz = 0; 311 G4double eab = availableEnergy; 312 G4int ieab = G4int( eab*5.0 ); 313 314 G4double supp[] = {0., 0.4, 0.55, 0.65, 0.75, 0.82, 0.86, 0.90, 0.94, 0.98}; 315 if( (ieab <= 9) && (G4UniformRand() >= supp[ieab])) { 316 // Suppress high multiplicity events at low momentum 317 // only one additional pion will be produced 318 G4double w0, wp, wm, wt, ran; 319 if (targetCode == neutronCode) { 320 // target is a neutron 321 w0 = - sqr(1.+protb)/(2.*c*c); 322 w0 = std::exp(w0); 323 wm = - sqr(-1.+protb)/(2.*c*c); 324 wm = std::exp(wm); 325 w0 = w0/2.; 326 wm = wm*1.5; 327 if (G4UniformRand() < w0/(w0+wm) ) { 328 np = 0; 329 nm = 0; 330 nz = 1; 331 } else { 332 np = 0; 333 nm = 1; 334 nz = 0; 335 } 336 337 } else { 338 // target is a proton 339 w0 = -sqr(1.+neutb)/(2.*c*c); 340 wp = w0 = std::exp(w0); 341 wm = -sqr(-1.+neutb)/(2.*c*c); 342 wm = std::exp(wm); 343 wt = w0+wp+wm; 344 wp = w0+wp; 345 ran = G4UniformRand(); 346 if ( ran < w0/wt) { 347 np = 0; 348 nm = 0; 349 nz = 1; 350 } else if (ran < wp/wt) { 351 np = 1; 352 nm = 0; 353 nz = 0; 354 } else { 355 np = 0; 356 nm = 1; 357 nz = 0; 358 } 359 } 360 } else { 361 // number of total particles vs. centre of mass Energy - 2*proton mass 362 363 G4double aleab = std::log(availableEnergy); 364 G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514 365 + aleab*(0.117712+0.0136912*aleab))) - 2.0; 366 367 // Normalization constant for kno-distribution. 368 // Calculate first the sum of all constants, check for numerical problems. 369 G4double test, dum, anpn = 0.0; 370 371 for (nt=1; nt<=numSec; nt++) { 372 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) ); 373 dum = pi*nt/(2.0*n*n); 374 if (std::fabs(dum) < 1.0) { 375 if( test >= 1.0e-10 )anpn += dum*test; 376 } else { 377 anpn += dum*test; 378 } 379 } 380 381 G4double ran = G4UniformRand(); 382 G4double excs = 0.0; 383 if( targetCode == protonCode ) 384 { 385 counter = -1; 386 for( np=0; np<numSec/3; np++ ) 387 { 388 for( nm=std::max(0,np-1); nm<=(np+1); nm++ ) 389 { 390 for (nz=0; nz<numSec/3; nz++) { 391 if (++counter < numMul) { 392 nt = np+nm+nz; 393 if ( (nt>0) && (nt<=numSec) ) { 394 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) ); 395 dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n); 396 if (std::fabs(dum) < 1.0) { 397 if( test >= 1.0e-10 )excs += dum*test; 398 } else { 399 excs += dum*test; 400 } 401 if (ran < excs) goto outOfLoop; //-----------------------> 402 } 403 } 404 } 405 } 406 } 407 408 // 3 previous loops continued to the end 409 inElastic = false; // quasi-elastic scattering 410 return; 411 } 412 else 413 { // target must be a neutron 414 counter = -1; 415 for( np=0; np<numSec/3; np++ ) 416 { 417 for( nm=np; nm<=(np+2); nm++ ) 418 { 419 for (nz=0; nz<numSec/3; nz++) { 420 if (++counter < numMul) { 421 nt = np+nm+nz; 422 if ( (nt>=1) && (nt<=numSec) ) { 423 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) ); 424 dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n); 425 if (std::fabs(dum) < 1.0) { 426 if( test >= 1.0e-10 )excs += dum*test; 427 } else { 428 excs += dum*test; 429 } 430 if (ran < excs) goto outOfLoop; // --------------------------> 431 } 432 } 433 } 434 } 435 } 436 // 3 previous loops continued to the end 437 inElastic = false; // quasi-elastic scattering. 438 return; 439 } 440 } 441 outOfLoop: // <----------------------------------------------- 442 443 if( targetCode == neutronCode) 444 { 445 if( np == nm) 446 { 447 } 448 else if (np == (nm-1)) 449 { 450 if( G4UniformRand() < 0.5) 451 { 452 pv[0] = KaonPlus; 453 } 454 else 455 { 456 pv[1] = Proton; 457 } 458 } 459 else 460 { 461 pv[0] = KaonPlus; 462 pv[1] = Proton; 463 } 464 } 465 else 466 { 467 if( np == nm ) 468 { 469 if( G4UniformRand() < 0.25) 470 { 471 pv[0] = KaonPlus; 472 pv[1] = Neutron; 473 } 474 else 475 { 476 } 477 } 478 else if ( np == (nm+1)) 479 { 480 pv[1] = Neutron; 481 } 482 else 483 { 484 pv[0] = KaonPlus; 485 } 486 } 487 488 nt = np + nm + nz; 489 while (nt > 0) { 490 G4double ran = G4UniformRand(); 491 if (ran < (G4double)np/nt) { 492 if (np > 0) { 493 pv[vecLen++] = PionPlus; 494 np--; 495 } 496 } else if ( ran < (G4double)(np+nm)/nt) { 497 if (nm > 0) { 498 pv[vecLen++] = PionMinus; 499 nm--; 500 } 501 } else { 502 if (nz > 0) { 503 pv[vecLen++] = PionZero; 504 nz--; 505 } 506 } 507 nt = np + nm + nz; 508 } 509 510 if (verboseLevel > 1) { 511 G4cout << "Particles produced: " ; 512 G4cout << pv[0].getName() << " " ; 513 G4cout << pv[1].getName() << " " ; 514 for (i=2; i < vecLen; i++) G4cout << pv[i].getName() << " " ; 515 G4cout << G4endl; 516 } 517 518 return; 519 } 520 521 522 void 523 G4HEKaonZeroLongInelastic::FirstIntInCasAntiKaonZero(G4bool &inElastic, 524 const G4double availableEnergy, 525 G4HEVector pv[], 526 G4int &vecLen, 527 G4HEVector incidentParticle, 528 G4HEVector targetParticle ) 529 530 // AntiKaon0 undergoes interaction with nucleon within a nucleus. Check if it is 531 // energetically possible to produce pions/kaons. In not, assume nuclear excitation 532 // occurs and input particle is degraded in energy. No other particles are produced. 533 // If reaction is possible, find the correct number of pions/protons/neutrons 534 // produced using an interpolation to multiplicity data. Replace some pions or 535 // protons/neutrons by kaons or strange baryons according to the average 536 // multiplicity per inelastic reaction. 537 538 { 539 static const G4double expxu = std::log(MAXFLOAT); // upper bound for arg. of exp 540 static const G4double expxl = -expxu; // lower bound for arg. of exp 541 542 static const G4double protb = 0.7; 543 static const G4double neutb = 0.7; 544 static const G4double c = 1.25; 545 546 static const G4int numMul = 1200; 547 static const G4int numSec = 60; 548 549 G4int neutronCode = Neutron.getCode(); 550 G4int protonCode = Proton.getCode(); 551 G4int kaonMinusCode = KaonMinus.getCode(); 552 G4int kaonZeroCode = KaonZero.getCode(); 553 G4int antiKaonZeroCode = AntiKaonZero.getCode(); 554 555 G4int targetCode = targetParticle.getCode(); 556 G4double incidentTotalMomentum = incidentParticle.getTotalMomentum(); 557 558 static G4bool first = true; 559 static G4double protmul[numMul], protnorm[numSec]; // proton constants 560 static G4double neutmul[numMul], neutnorm[numSec]; // neutron constants 561 562 // misc. local variables 563 // np = number of pi+, nm = number of pi-, nz = number of pi0 564 565 G4int i, counter, nt, np, nm, nz; 566 567 if(first) { 568 // compute normalization constants, this will only be done once 569 first = false; 570 for( i=0; i<numMul; i++ )protmul[i] = 0.0; 571 for( i=0; i<numSec; i++ )protnorm[i] = 0.0; 572 counter = -1; 573 for(np=0; np<(numSec/3); np++) { 574 for(nm=std::max(0,np-2); nm<=np; nm++) { 575 for(nz=0; nz<numSec/3; nz++) { 576 if(++counter < numMul) { 577 nt = np+nm+nz; 578 if( (nt>0) && (nt<=numSec) ) { 579 protmul[counter] = pmltpc(np,nm,nz,nt,protb,c) ; 580 protnorm[nt-1] += protmul[counter]; 581 } 582 } 583 } 584 } 585 } 586 587 for( i=0; i<numMul; i++ )neutmul[i] = 0.0; 588 for( i=0; i<numSec; i++ )neutnorm[i] = 0.0; 589 counter = -1; 590 for(np=0; np<numSec/3; np++) { 591 for(nm=std::max(0,np-1); nm<=(np+1); nm++) { 592 for(nz=0; nz<numSec/3; nz++) { 593 if(++counter < numMul) { 594 nt = np+nm+nz; 595 if( (nt>0) && (nt<=numSec) ) { 596 neutmul[counter] = pmltpc(np,nm,nz,nt,neutb,c); 597 neutnorm[nt-1] += neutmul[counter]; 598 } 599 } 600 } 601 } 602 } 603 604 for(i=0; i<numSec; i++) { 605 if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i]; 606 if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i]; 607 } 608 } // end of initialization 609 610 // initialize the first two particles 611 // the same as beam and target 612 pv[0] = incidentParticle; 613 pv[1] = targetParticle; 614 vecLen = 2; 615 616 if (!inElastic || (availableEnergy <= PionPlus.getMass())) 617 return; 618 619 // Inelastic scattering 620 621 np = 0, nm = 0, nz = 0; 622 G4double cech[] = { 1., 1., 1., 0.70, 0.60, 0.55, 0.35, 0.25, 0.18, 0.15}; 623 G4int iplab = G4int( incidentTotalMomentum*5.); 624 if( (iplab < 10) && (G4UniformRand() < cech[iplab]) ) { 625 G4int iplab = std::min(19, G4int( incidentTotalMomentum*5.)); 626 G4double cnk0[] = {0.17, 0.18, 0.17, 0.24, 0.26, 0.20, 0.22, 0.21, 0.34, 0.45, 627 0.58, 0.55, 0.36, 0.29, 0.29, 0.32, 0.32, 0.33, 0.33, 0.33}; 628 if(G4UniformRand() < cnk0[iplab]) { 629 if(targetCode == protonCode) { 630 return; 631 } else { 632 pv[0] = KaonMinus; 633 pv[1] = Proton; 634 return; 635 } 636 } 637 638 G4double ran = G4UniformRand(); 639 if(targetCode == protonCode) { 640 641 // target is a proton 642 if( ran < 0.25 ) { 643 ; 644 } else if (ran < 0.50) { 645 pv[0] = PionPlus; 646 pv[1] = SigmaZero; 647 } else if (ran < 0.75) { 648 ; 649 } else { 650 pv[0] = PionPlus; 651 pv[1] = Lambda; 652 } 653 } else { 654 655 // target is a neutron 656 if( ran < 0.25 ) { 657 pv[0] = PionMinus; 658 pv[1] = SigmaPlus; 659 } else if (ran < 0.50) { 660 pv[0] = PionZero; 661 pv[1] = SigmaZero; 662 } else if (ran < 0.75) { 663 pv[0] = PionPlus; 664 pv[1] = SigmaMinus; 665 } else { 666 pv[0] = PionZero; 667 pv[1] = Lambda; 668 } 669 } 670 return; 671 672 } else { 673 // number of total particles vs. centre of mass Energy - 2*proton mass 674 675 G4double aleab = std::log(availableEnergy); 676 G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514 677 + aleab*(0.117712+0.0136912*aleab))) - 2.0; 678 679 // Normalization constant for kno-distribution. 680 // Calculate first the sum of all constants, check for numerical problems. 681 G4double test, dum, anpn = 0.0; 682 683 for (nt=1; nt<=numSec; nt++) { 684 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) ); 685 dum = pi*nt/(2.0*n*n); 686 if (std::fabs(dum) < 1.0) { 687 if( test >= 1.0e-10 )anpn += dum*test; 688 } else { 689 anpn += dum*test; 690 } 691 } 692 693 G4double ran = G4UniformRand(); 694 G4double excs = 0.0; 695 if (targetCode == protonCode) { 696 counter = -1; 697 for (np=0; np<numSec/3; np++) { 698 for (nm=std::max(0,np-2); nm<=np; nm++) { 699 for (nz=0; nz<numSec/3; nz++) { 700 if (++counter < numMul) { 701 nt = np+nm+nz; 702 if( (nt>0) && (nt<=numSec) ) { 703 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) ); 704 dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n); 705 706 if (std::fabs(dum) < 1.0) { 707 if( test >= 1.0e-10 )excs += dum*test; 708 } else { 709 excs += dum*test; 710 } 711 712 if (ran < excs) goto outOfLoop; //-----------------------> 713 } 714 } 715 } 716 } 717 } 718 // 3 previous loops continued to the end 719 inElastic = false; // quasi-elastic scattering 720 return; 721 722 } else { // target must be a neutron 723 counter = -1; 724 for (np=0; np<numSec/3; np++) { 725 for (nm=std::max(0,np-1); nm<=(np+1); nm++) { 726 for (nz=0; nz<numSec/3; nz++) { 727 if (++counter < numMul) { 728 nt = np+nm+nz; 729 if( (nt>=1) && (nt<=numSec) ) { 730 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) ); 731 dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n); 732 733 if (std::fabs(dum) < 1.0) { 734 if( test >= 1.0e-10 )excs += dum*test; 735 } else { 736 excs += dum*test; 737 } 738 739 if (ran < excs) goto outOfLoop; // --------------------------> 740 } 741 } 742 } 743 } 744 } 745 // 3 previous loops continued to the end 746 inElastic = false; // quasi-elastic scattering. 747 return; 748 } 749 } 750 outOfLoop: // <------------------------------------------------------------------------ 751 752 if( targetCode == protonCode) 753 { 754 if( np == nm) 755 { 756 } 757 else if (np == (1+nm)) 758 { 759 if( G4UniformRand() < 0.5) 760 { 761 pv[0] = KaonMinus; 762 } 763 else 764 { 765 pv[1] = Neutron; 766 } 767 } 768 else 769 { 770 pv[0] = KaonMinus; 771 pv[1] = Neutron; 772 } 773 } 774 else 775 { 776 if( np == nm) 777 { 778 if( G4UniformRand() < 0.75) 779 { 780 } 781 else 782 { 783 pv[0] = KaonMinus; 784 pv[1] = Proton; 785 } 786 } 787 else if ( np == (1+nm)) 788 { 789 pv[0] = KaonMinus; 790 } 791 else 792 { 793 pv[1] = Proton; 794 } 795 } 796 797 798 if( G4UniformRand() < 0.5 ) 799 { 800 if( ( (pv[0].getCode() == kaonMinusCode) 801 && (pv[1].getCode() == neutronCode) ) 802 || ( (pv[0].getCode() == kaonZeroCode) 803 && (pv[1].getCode() == protonCode) ) 804 || ( (pv[0].getCode() == antiKaonZeroCode) 805 && (pv[1].getCode() == protonCode) ) ) 806 { 807 G4double ran = G4UniformRand(); 808 if( pv[1].getCode() == protonCode) 809 { 810 if(ran < 0.68) 811 { 812 pv[0] = PionPlus; 813 pv[1] = Lambda; 814 } 815 else if (ran < 0.84) 816 { 817 pv[0] = PionZero; 818 pv[1] = SigmaPlus; 819 } 820 else 821 { 822 pv[0] = PionPlus; 823 pv[1] = SigmaZero; 824 } 825 } 826 else 827 { 828 if(ran < 0.68) 829 { 830 pv[0] = PionMinus; 831 pv[1] = Lambda; 832 } 833 else if (ran < 0.84) 834 { 835 pv[0] = PionMinus; 836 pv[1] = SigmaZero; 837 } 838 else 839 { 840 pv[0] = PionZero; 841 pv[1] = SigmaMinus; 842 } 843 } 844 } 845 else 846 { 847 G4double ran = G4UniformRand(); 848 if (ran < 0.67) 849 { 850 pv[0] = PionZero; 851 pv[1] = Lambda; 852 } 853 else if (ran < 0.78) 854 { 855 pv[0] = PionMinus; 856 pv[1] = SigmaPlus; 857 } 858 else if (ran < 0.89) 859 { 860 pv[0] = PionZero; 861 pv[1] = SigmaZero; 862 } 863 else 864 { 865 pv[0] = PionPlus; 866 pv[1] = SigmaMinus; 867 } 868 } 869 } 870 871 nt = np + nm + nz; 872 while ( nt > 0) { 873 G4double ran = G4UniformRand(); 874 if ( ran < (G4double)np/nt) { 875 if( np > 0 ) { 876 pv[vecLen++] = PionPlus; 877 np--; 878 } 879 } else if (ran < (G4double)(np+nm)/nt) { 880 if( nm > 0 ) { 881 pv[vecLen++] = PionMinus; 882 nm--; 883 } 884 } else { 885 if( nz > 0 ) { 886 pv[vecLen++] = PionZero; 887 nz--; 888 } 889 } 890 nt = np + nm + nz; 891 } 892 893 if (verboseLevel > 1) { 894 G4cout << "Particles produced: " ; 895 G4cout << pv[0].getName() << " " ; 896 G4cout << pv[1].getName() << " " ; 897 for (i=2; i < vecLen; i++) G4cout << pv[i].getName() << " " ; 898 G4cout << G4endl; 899 } 900 901 return; 902 } -
trunk/source/processes/hadronic/models/high_energy/src/G4HEKaonZeroShortInelastic.cc
r1228 r1315 24 24 // ******************************************************************** 25 25 // 26 // 27 // $Id: G4HEKaonZeroShortInelastic.cc,v 1.10 2006/06/29 20:30:24 gunter Exp $ 28 // GEANT4 tag $Name: geant4-09-03 $ 26 // $Id: G4HEKaonZeroShortInelastic.cc,v 1.11 2010/02/09 21:58:57 dennis Exp $ 27 // GEANT4 tag $Name: geant4-09-04-beta-cand-01 $ 29 28 // 30 29 // … … 38 37 // and the low energy two-body model. Not included are the low energy stuff like 39 38 // nuclear reactions, nuclear fission without any cascading and all processes for 40 // particles at rest. 41 // First work done by J.L.Chuma and F.W.Jones, TRIUMF, June 96. 42 // H. Fesefeldt, RWTH-Aachen, 23-October-1996 43 // Last modified: 29-July-1998 39 // particles at rest. 40 // 41 // New version by D.H. Wright (SLAC) to fix seg fault in old version 42 // 21 January 2010 43 44 44 45 45 #include "G4HEKaonZeroShortInelastic.hh" 46 46 47 G4HadFinalState * G4HEKaonZeroShortInelastic:: 48 ApplyYourself( const G4HadProjectile &aTrack, G4Nucleus &targetNucleus ) 49 { 50 G4HEKaonZeroInelastic theKaonZeroInelastic; 51 G4HEAntiKaonZeroInelastic theAntiKaonZeroInelastic; 52 theKaonZeroInelastic.SetVerboseLevel(verboseLevel); 53 theAntiKaonZeroInelastic.SetVerboseLevel(verboseLevel); 47 G4HadFinalState* G4HEKaonZeroShortInelastic:: 48 ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus &targetNucleus) 49 { 50 G4HEVector * pv = new G4HEVector[MAXPART]; 51 const G4HadProjectile *aParticle = &aTrack; 52 const G4double atomicWeight = targetNucleus.GetN(); 53 const G4double atomicNumber = targetNucleus.GetZ(); 54 G4HEVector incidentParticle(aParticle); 55 56 G4int incidentCode = incidentParticle.getCode(); 57 G4double incidentMass = incidentParticle.getMass(); 58 G4double incidentTotalEnergy = incidentParticle.getEnergy(); 59 G4double incidentTotalMomentum = incidentParticle.getTotalMomentum(); 60 G4double incidentKineticEnergy = incidentTotalEnergy - incidentMass; 61 62 if(incidentKineticEnergy < 1) 63 G4cout << "GHEKaonZeroShortInelastic: incident energy < 1 GeV" << G4endl; 64 65 if(verboseLevel > 1) { 66 G4cout << "G4HEKaonZeroShortInelastic::ApplyYourself" << G4endl; 67 G4cout << "incident particle " << incidentParticle.getName() 68 << "mass " << incidentMass 69 << "kinetic energy " << incidentKineticEnergy 70 << G4endl; 71 G4cout << "target material with (A,Z) = (" 72 << atomicWeight << "," << atomicNumber << ")" << G4endl; 73 } 54 74 55 if(G4UniformRand() < 0.50) 56 { 57 return theKaonZeroInelastic.ApplyYourself(aTrack, targetNucleus); 58 } 75 G4double inelasticity = NuclearInelasticity(incidentKineticEnergy, 76 atomicWeight, atomicNumber); 77 if(verboseLevel > 1) 78 G4cout << "nuclear inelasticity = " << inelasticity << G4endl; 79 80 incidentKineticEnergy -= inelasticity; 81 82 G4double excitationEnergyGNP = 0.; 83 G4double excitationEnergyDTA = 0.; 84 85 G4double excitation = NuclearExcitation(incidentKineticEnergy, 86 atomicWeight, atomicNumber, 87 excitationEnergyGNP, 88 excitationEnergyDTA); 89 if(verboseLevel > 1) 90 G4cout << "nuclear excitation = " << excitation << excitationEnergyGNP 91 << excitationEnergyDTA << G4endl; 92 93 incidentKineticEnergy -= excitation; 94 incidentTotalEnergy = incidentKineticEnergy + incidentMass; 95 incidentTotalMomentum = std::sqrt( (incidentTotalEnergy-incidentMass) 96 *(incidentTotalEnergy+incidentMass)); 97 98 99 G4HEVector targetParticle; 100 if(G4UniformRand() < atomicNumber/atomicWeight) { 101 targetParticle.setDefinition("Proton"); 102 } else { 103 targetParticle.setDefinition("Neutron"); 104 } 105 106 G4double targetMass = targetParticle.getMass(); 107 G4double centerOfMassEnergy = std::sqrt( incidentMass*incidentMass + targetMass*targetMass 108 + 2.0*targetMass*incidentTotalEnergy); 109 G4double availableEnergy = centerOfMassEnergy - targetMass - incidentMass; 110 // this was the meaning of inElastic in the 111 // original Gheisha stand-alone version. 112 // G4bool inElastic = InElasticCrossSectionInFirstInt 113 // (availableEnergy, incidentCode, incidentTotalMomentum); 114 // for unknown reasons, it has been replaced by the following code in Geant??? 115 116 G4bool inElastic = true; 117 // if (G4UniformRand() < elasticCrossSection/totalCrossSection) inElastic = false; 118 119 vecLength = 0; 120 121 if(verboseLevel > 1) 122 G4cout << "ApplyYourself: CallFirstIntInCascade for particle " 123 << incidentCode << G4endl; 124 125 G4bool successful = false; 126 127 if(inElastic || (!inElastic && atomicWeight < 1.5)) { 128 129 // Split K0L into K0 and K0bar 130 if (G4UniformRand() < 0.5) 131 FirstIntInCasAntiKaonZero(inElastic, availableEnergy, pv, vecLength, 132 incidentParticle, targetParticle ); 59 133 else 60 { 61 return theAntiKaonZeroInelastic.ApplyYourself(aTrack, targetNucleus); 62 } 134 FirstIntInCasKaonZero(inElastic, availableEnergy, pv, vecLength, 135 incidentParticle, targetParticle, atomicWeight ); 136 137 // Do nuclear interaction with either K0 or K0bar 138 if ((vecLength > 0) && (availableEnergy > 1.)) 139 StrangeParticlePairProduction(availableEnergy, centerOfMassEnergy, 140 pv, vecLength, 141 incidentParticle, targetParticle); 142 HighEnergyCascading(successful, pv, vecLength, 143 excitationEnergyGNP, excitationEnergyDTA, 144 incidentParticle, targetParticle, 145 atomicWeight, atomicNumber); 146 if (!successful) 147 HighEnergyClusterProduction(successful, pv, vecLength, 148 excitationEnergyGNP, excitationEnergyDTA, 149 incidentParticle, targetParticle, 150 atomicWeight, atomicNumber); 151 if (!successful) 152 MediumEnergyCascading(successful, pv, vecLength, 153 excitationEnergyGNP, excitationEnergyDTA, 154 incidentParticle, targetParticle, 155 atomicWeight, atomicNumber); 156 157 if (!successful) 158 MediumEnergyClusterProduction(successful, pv, vecLength, 159 excitationEnergyGNP, excitationEnergyDTA, 160 incidentParticle, targetParticle, 161 atomicWeight, atomicNumber); 162 if (!successful) 163 QuasiElasticScattering(successful, pv, vecLength, 164 excitationEnergyGNP, excitationEnergyDTA, 165 incidentParticle, targetParticle, 166 atomicWeight, atomicNumber); 167 } 168 169 if (!successful) 170 ElasticScattering(successful, pv, vecLength, 171 incidentParticle, 172 atomicWeight, atomicNumber); 173 174 if (!successful) 175 G4cout << "GHEInelasticInteraction::ApplyYourself fails to produce final state particles" 176 << G4endl; 177 178 // Check for K0, K0bar and change particle types to K0L, K0S if necessary 179 G4int kcode; 180 for (G4int i = 0; i < vecLength; i++) { 181 kcode = pv[i].getCode(); 182 if (kcode == KaonZero.getCode() || kcode == AntiKaonZero.getCode()) { 183 if (G4UniformRand() < 0.5) 184 pv[i] = KaonZeroShort; 185 else 186 pv[i] = KaonZeroLong; 187 } 63 188 } 64 65 66 67 68 69 70 71 72 189 190 // ................ 191 192 FillParticleChange(pv, vecLength); 193 delete [] pv; 194 theParticleChange.SetStatusChange(stopAndKill); 195 return & theParticleChange; 196 } 197 198 199 void 200 G4HEKaonZeroShortInelastic::FirstIntInCasKaonZero(G4bool &inElastic, 201 const G4double availableEnergy, 202 G4HEVector pv[], 203 G4int &vecLen, 204 G4HEVector incidentParticle, 205 G4HEVector targetParticle, 206 const G4double atomicWeight) 207 208 // Kaon0 undergoes interaction with nucleon within a nucleus. Check if it is 209 // energetically possible to produce pions/kaons. In not, assume nuclear excitation 210 // occurs and input particle is degraded in energy. No other particles are produced. 211 // If reaction is possible, find the correct number of pions/protons/neutrons 212 // produced using an interpolation to multiplicity data. Replace some pions or 213 // protons/neutrons by kaons or strange baryons according to the average 214 // multiplicity per inelastic reaction. 215 216 { 217 static const G4double expxu = std::log(MAXFLOAT); // upper bound for arg. of exp 218 static const G4double expxl = -expxu; // lower bound for arg. of exp 219 220 static const G4double protb = 0.7; 221 static const G4double neutb = 0.7; 222 static const G4double c = 1.25; 223 224 static const G4int numMul = 1200; 225 static const G4int numSec = 60; 226 227 G4int neutronCode = Neutron.getCode(); 228 G4int protonCode = Proton.getCode(); 229 230 G4int targetCode = targetParticle.getCode(); 231 G4double incidentTotalMomentum = incidentParticle.getTotalMomentum(); 232 233 static G4bool first = true; 234 static G4double protmul[numMul], protnorm[numSec]; // proton constants 235 static G4double neutmul[numMul], neutnorm[numSec]; // neutron constants 236 237 // misc. local variables 238 // np = number of pi+, nm = number of pi-, nz = number of pi0 239 240 G4int i, counter, nt, np, nm, nz; 241 242 if (first) { 243 // compute normalization constants, this will only be done once 244 first = false; 245 for( i=0; i<numMul; i++ )protmul[i] = 0.0; 246 for( i=0; i<numSec; i++ )protnorm[i] = 0.0; 247 counter = -1; 248 for (np=0; np<(numSec/3); np++) { 249 for (nm=std::max(0,np-1); nm<=(np+1); nm++) { 250 for (nz=0; nz<numSec/3; nz++) { 251 if (++counter < numMul) { 252 nt = np+nm+nz; 253 if( (nt>0) && (nt<=numSec) ) { 254 protmul[counter] = pmltpc(np,nm,nz,nt,protb,c) ; 255 protnorm[nt-1] += protmul[counter]; 256 } 257 } 258 } 259 } 260 } 261 262 for( i=0; i<numMul; i++ )neutmul[i] = 0.0; 263 for( i=0; i<numSec; i++ )neutnorm[i] = 0.0; 264 counter = -1; 265 for (np=0; np<numSec/3; np++) { 266 for (nm=np; nm<=(np+2); nm++) { 267 for (nz=0; nz<numSec/3; nz++) { 268 if (++counter < numMul) { 269 nt = np+nm+nz; 270 if( (nt>0) && (nt<=numSec) ) { 271 neutmul[counter] = pmltpc(np,nm,nz,nt,neutb,c); 272 neutnorm[nt-1] += neutmul[counter]; 273 } 274 } 275 } 276 } 277 } 278 279 for (i=0; i<numSec; i++) { 280 if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i]; 281 if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i]; 282 } 283 } // end of initialization 284 285 286 // Initialize the first two particles 287 // the same as beam and target 288 pv[0] = incidentParticle; 289 pv[1] = targetParticle; 290 vecLen = 2; 291 292 if( !inElastic ) { 293 // quasi-elastic scattering, no pions produced 294 if( targetCode == protonCode) { 295 G4double cech[] = {0.33,0.27,0.29,0.31,0.27,0.18,0.13,0.10,0.09,0.07}; 296 G4int iplab = G4int( std::min( 9.0, incidentTotalMomentum*5. ) ); 297 if( G4UniformRand() < cech[iplab]/std::pow(atomicWeight,0.42)) { 298 // charge exchange K+ n -> K0 p 299 pv[0] = KaonPlus; 300 pv[1] = Neutron; 301 } 302 } 303 return; 304 } else if (availableEnergy <= PionPlus.getMass()) { 305 return; 306 } 307 308 // Inelastic scattering 309 310 np = 0, nm = 0, nz = 0; 311 G4double eab = availableEnergy; 312 G4int ieab = G4int( eab*5.0 ); 313 314 G4double supp[] = {0., 0.4, 0.55, 0.65, 0.75, 0.82, 0.86, 0.90, 0.94, 0.98}; 315 if( (ieab <= 9) && (G4UniformRand() >= supp[ieab])) { 316 // Suppress high multiplicity events at low momentum 317 // only one additional pion will be produced 318 G4double w0, wp, wm, wt, ran; 319 if (targetCode == neutronCode) { 320 // target is a neutron 321 w0 = - sqr(1.+protb)/(2.*c*c); 322 w0 = std::exp(w0); 323 wm = - sqr(-1.+protb)/(2.*c*c); 324 wm = std::exp(wm); 325 w0 = w0/2.; 326 wm = wm*1.5; 327 if (G4UniformRand() < w0/(w0+wm) ) { 328 np = 0; 329 nm = 0; 330 nz = 1; 331 } else { 332 np = 0; 333 nm = 1; 334 nz = 0; 335 } 336 337 } else { 338 // target is a proton 339 w0 = -sqr(1.+neutb)/(2.*c*c); 340 wp = w0 = std::exp(w0); 341 wm = -sqr(-1.+neutb)/(2.*c*c); 342 wm = std::exp(wm); 343 wt = w0+wp+wm; 344 wp = w0+wp; 345 ran = G4UniformRand(); 346 if ( ran < w0/wt) { 347 np = 0; 348 nm = 0; 349 nz = 1; 350 } else if (ran < wp/wt) { 351 np = 1; 352 nm = 0; 353 nz = 0; 354 } else { 355 np = 0; 356 nm = 1; 357 nz = 0; 358 } 359 } 360 } else { 361 // number of total particles vs. centre of mass Energy - 2*proton mass 362 363 G4double aleab = std::log(availableEnergy); 364 G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514 365 + aleab*(0.117712+0.0136912*aleab))) - 2.0; 366 367 // Normalization constant for kno-distribution. 368 // Calculate first the sum of all constants, check for numerical problems. 369 G4double test, dum, anpn = 0.0; 370 371 for (nt=1; nt<=numSec; nt++) { 372 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) ); 373 dum = pi*nt/(2.0*n*n); 374 if (std::fabs(dum) < 1.0) { 375 if( test >= 1.0e-10 )anpn += dum*test; 376 } else { 377 anpn += dum*test; 378 } 379 } 380 381 G4double ran = G4UniformRand(); 382 G4double excs = 0.0; 383 if( targetCode == protonCode ) 384 { 385 counter = -1; 386 for( np=0; np<numSec/3; np++ ) 387 { 388 for( nm=std::max(0,np-1); nm<=(np+1); nm++ ) 389 { 390 for (nz=0; nz<numSec/3; nz++) { 391 if (++counter < numMul) { 392 nt = np+nm+nz; 393 if ( (nt>0) && (nt<=numSec) ) { 394 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) ); 395 dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n); 396 if (std::fabs(dum) < 1.0) { 397 if( test >= 1.0e-10 )excs += dum*test; 398 } else { 399 excs += dum*test; 400 } 401 if (ran < excs) goto outOfLoop; //-----------------------> 402 } 403 } 404 } 405 } 406 } 407 408 // 3 previous loops continued to the end 409 inElastic = false; // quasi-elastic scattering 410 return; 411 } 412 else 413 { // target must be a neutron 414 counter = -1; 415 for( np=0; np<numSec/3; np++ ) 416 { 417 for( nm=np; nm<=(np+2); nm++ ) 418 { 419 for (nz=0; nz<numSec/3; nz++) { 420 if (++counter < numMul) { 421 nt = np+nm+nz; 422 if ( (nt>=1) && (nt<=numSec) ) { 423 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) ); 424 dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n); 425 if (std::fabs(dum) < 1.0) { 426 if( test >= 1.0e-10 )excs += dum*test; 427 } else { 428 excs += dum*test; 429 } 430 if (ran < excs) goto outOfLoop; // --------------------------> 431 } 432 } 433 } 434 } 435 } 436 // 3 previous loops continued to the end 437 inElastic = false; // quasi-elastic scattering. 438 return; 439 } 440 } 441 outOfLoop: // <----------------------------------------------- 442 443 if( targetCode == neutronCode) 444 { 445 if( np == nm) 446 { 447 } 448 else if (np == (nm-1)) 449 { 450 if( G4UniformRand() < 0.5) 451 { 452 pv[0] = KaonPlus; 453 } 454 else 455 { 456 pv[1] = Proton; 457 } 458 } 459 else 460 { 461 pv[0] = KaonPlus; 462 pv[1] = Proton; 463 } 464 } 465 else 466 { 467 if( np == nm ) 468 { 469 if( G4UniformRand() < 0.25) 470 { 471 pv[0] = KaonPlus; 472 pv[1] = Neutron; 473 } 474 else 475 { 476 } 477 } 478 else if ( np == (nm+1)) 479 { 480 pv[1] = Neutron; 481 } 482 else 483 { 484 pv[0] = KaonPlus; 485 } 486 } 487 488 nt = np + nm + nz; 489 while (nt > 0) { 490 G4double ran = G4UniformRand(); 491 if (ran < (G4double)np/nt) { 492 if (np > 0) { 493 pv[vecLen++] = PionPlus; 494 np--; 495 } 496 } else if ( ran < (G4double)(np+nm)/nt) { 497 if (nm > 0) { 498 pv[vecLen++] = PionMinus; 499 nm--; 500 } 501 } else { 502 if (nz > 0) { 503 pv[vecLen++] = PionZero; 504 nz--; 505 } 506 } 507 nt = np + nm + nz; 508 } 509 510 if (verboseLevel > 1) { 511 G4cout << "Particles produced: " ; 512 G4cout << pv[0].getName() << " " ; 513 G4cout << pv[1].getName() << " " ; 514 for (i=2; i < vecLen; i++) G4cout << pv[i].getName() << " " ; 515 G4cout << G4endl; 516 } 517 518 return; 519 } 520 521 522 void 523 G4HEKaonZeroShortInelastic::FirstIntInCasAntiKaonZero(G4bool &inElastic, 524 const G4double availableEnergy, 525 G4HEVector pv[], 526 G4int &vecLen, 527 G4HEVector incidentParticle, 528 G4HEVector targetParticle ) 529 530 // AntiKaon0 undergoes interaction with nucleon within a nucleus. Check if it is 531 // energetically possible to produce pions/kaons. In not, assume nuclear excitation 532 // occurs and input particle is degraded in energy. No other particles are produced. 533 // If reaction is possible, find the correct number of pions/protons/neutrons 534 // produced using an interpolation to multiplicity data. Replace some pions or 535 // protons/neutrons by kaons or strange baryons according to the average 536 // multiplicity per inelastic reaction. 537 538 { 539 static const G4double expxu = std::log(MAXFLOAT); // upper bound for arg. of exp 540 static const G4double expxl = -expxu; // lower bound for arg. of exp 541 542 static const G4double protb = 0.7; 543 static const G4double neutb = 0.7; 544 static const G4double c = 1.25; 545 546 static const G4int numMul = 1200; 547 static const G4int numSec = 60; 548 549 G4int neutronCode = Neutron.getCode(); 550 G4int protonCode = Proton.getCode(); 551 G4int kaonMinusCode = KaonMinus.getCode(); 552 G4int kaonZeroCode = KaonZero.getCode(); 553 G4int antiKaonZeroCode = AntiKaonZero.getCode(); 554 555 G4int targetCode = targetParticle.getCode(); 556 G4double incidentTotalMomentum = incidentParticle.getTotalMomentum(); 557 558 static G4bool first = true; 559 static G4double protmul[numMul], protnorm[numSec]; // proton constants 560 static G4double neutmul[numMul], neutnorm[numSec]; // neutron constants 561 562 // misc. local variables 563 // np = number of pi+, nm = number of pi-, nz = number of pi0 564 565 G4int i, counter, nt, np, nm, nz; 566 567 if(first) { 568 // compute normalization constants, this will only be done once 569 first = false; 570 for( i=0; i<numMul; i++ )protmul[i] = 0.0; 571 for( i=0; i<numSec; i++ )protnorm[i] = 0.0; 572 counter = -1; 573 for(np=0; np<(numSec/3); np++) { 574 for(nm=std::max(0,np-2); nm<=np; nm++) { 575 for(nz=0; nz<numSec/3; nz++) { 576 if(++counter < numMul) { 577 nt = np+nm+nz; 578 if( (nt>0) && (nt<=numSec) ) { 579 protmul[counter] = pmltpc(np,nm,nz,nt,protb,c) ; 580 protnorm[nt-1] += protmul[counter]; 581 } 582 } 583 } 584 } 585 } 586 587 for( i=0; i<numMul; i++ )neutmul[i] = 0.0; 588 for( i=0; i<numSec; i++ )neutnorm[i] = 0.0; 589 counter = -1; 590 for(np=0; np<numSec/3; np++) { 591 for(nm=std::max(0,np-1); nm<=(np+1); nm++) { 592 for(nz=0; nz<numSec/3; nz++) { 593 if(++counter < numMul) { 594 nt = np+nm+nz; 595 if( (nt>0) && (nt<=numSec) ) { 596 neutmul[counter] = pmltpc(np,nm,nz,nt,neutb,c); 597 neutnorm[nt-1] += neutmul[counter]; 598 } 599 } 600 } 601 } 602 } 603 604 for(i=0; i<numSec; i++) { 605 if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i]; 606 if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i]; 607 } 608 } // end of initialization 609 610 // initialize the first two particles 611 // the same as beam and target 612 pv[0] = incidentParticle; 613 pv[1] = targetParticle; 614 vecLen = 2; 615 616 if (!inElastic || (availableEnergy <= PionPlus.getMass())) 617 return; 618 619 // Inelastic scattering 620 621 np = 0, nm = 0, nz = 0; 622 G4double cech[] = { 1., 1., 1., 0.70, 0.60, 0.55, 0.35, 0.25, 0.18, 0.15}; 623 G4int iplab = G4int( incidentTotalMomentum*5.); 624 if( (iplab < 10) && (G4UniformRand() < cech[iplab]) ) { 625 G4int iplab = std::min(19, G4int( incidentTotalMomentum*5.)); 626 G4double cnk0[] = {0.17, 0.18, 0.17, 0.24, 0.26, 0.20, 0.22, 0.21, 0.34, 0.45, 627 0.58, 0.55, 0.36, 0.29, 0.29, 0.32, 0.32, 0.33, 0.33, 0.33}; 628 if(G4UniformRand() < cnk0[iplab]) { 629 if(targetCode == protonCode) { 630 return; 631 } else { 632 pv[0] = KaonMinus; 633 pv[1] = Proton; 634 return; 635 } 636 } 637 638 G4double ran = G4UniformRand(); 639 if(targetCode == protonCode) { 640 641 // target is a proton 642 if( ran < 0.25 ) { 643 ; 644 } else if (ran < 0.50) { 645 pv[0] = PionPlus; 646 pv[1] = SigmaZero; 647 } else if (ran < 0.75) { 648 ; 649 } else { 650 pv[0] = PionPlus; 651 pv[1] = Lambda; 652 } 653 } else { 654 655 // target is a neutron 656 if( ran < 0.25 ) { 657 pv[0] = PionMinus; 658 pv[1] = SigmaPlus; 659 } else if (ran < 0.50) { 660 pv[0] = PionZero; 661 pv[1] = SigmaZero; 662 } else if (ran < 0.75) { 663 pv[0] = PionPlus; 664 pv[1] = SigmaMinus; 665 } else { 666 pv[0] = PionZero; 667 pv[1] = Lambda; 668 } 669 } 670 return; 671 672 } else { 673 // number of total particles vs. centre of mass Energy - 2*proton mass 674 675 G4double aleab = std::log(availableEnergy); 676 G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514 677 + aleab*(0.117712+0.0136912*aleab))) - 2.0; 678 679 // Normalization constant for kno-distribution. 680 // Calculate first the sum of all constants, check for numerical problems. 681 G4double test, dum, anpn = 0.0; 682 683 for (nt=1; nt<=numSec; nt++) { 684 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) ); 685 dum = pi*nt/(2.0*n*n); 686 if (std::fabs(dum) < 1.0) { 687 if( test >= 1.0e-10 )anpn += dum*test; 688 } else { 689 anpn += dum*test; 690 } 691 } 692 693 G4double ran = G4UniformRand(); 694 G4double excs = 0.0; 695 if (targetCode == protonCode) { 696 counter = -1; 697 for (np=0; np<numSec/3; np++) { 698 for (nm=std::max(0,np-2); nm<=np; nm++) { 699 for (nz=0; nz<numSec/3; nz++) { 700 if (++counter < numMul) { 701 nt = np+nm+nz; 702 if( (nt>0) && (nt<=numSec) ) { 703 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) ); 704 dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n); 705 706 if (std::fabs(dum) < 1.0) { 707 if( test >= 1.0e-10 )excs += dum*test; 708 } else { 709 excs += dum*test; 710 } 711 712 if (ran < excs) goto outOfLoop; //-----------------------> 713 } 714 } 715 } 716 } 717 } 718 // 3 previous loops continued to the end 719 inElastic = false; // quasi-elastic scattering 720 return; 721 722 } else { // target must be a neutron 723 counter = -1; 724 for (np=0; np<numSec/3; np++) { 725 for (nm=std::max(0,np-1); nm<=(np+1); nm++) { 726 for (nz=0; nz<numSec/3; nz++) { 727 if (++counter < numMul) { 728 nt = np+nm+nz; 729 if( (nt>=1) && (nt<=numSec) ) { 730 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) ); 731 dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n); 732 733 if (std::fabs(dum) < 1.0) { 734 if( test >= 1.0e-10 )excs += dum*test; 735 } else { 736 excs += dum*test; 737 } 738 739 if (ran < excs) goto outOfLoop; // --------------------------> 740 } 741 } 742 } 743 } 744 } 745 // 3 previous loops continued to the end 746 inElastic = false; // quasi-elastic scattering. 747 return; 748 } 749 } 750 outOfLoop: // <------------------------------------------------------------------------ 751 752 if( targetCode == protonCode) 753 { 754 if( np == nm) 755 { 756 } 757 else if (np == (1+nm)) 758 { 759 if( G4UniformRand() < 0.5) 760 { 761 pv[0] = KaonMinus; 762 } 763 else 764 { 765 pv[1] = Neutron; 766 } 767 } 768 else 769 { 770 pv[0] = KaonMinus; 771 pv[1] = Neutron; 772 } 773 } 774 else 775 { 776 if( np == nm) 777 { 778 if( G4UniformRand() < 0.75) 779 { 780 } 781 else 782 { 783 pv[0] = KaonMinus; 784 pv[1] = Proton; 785 } 786 } 787 else if ( np == (1+nm)) 788 { 789 pv[0] = KaonMinus; 790 } 791 else 792 { 793 pv[1] = Proton; 794 } 795 } 796 797 798 if( G4UniformRand() < 0.5 ) 799 { 800 if( ( (pv[0].getCode() == kaonMinusCode) 801 && (pv[1].getCode() == neutronCode) ) 802 || ( (pv[0].getCode() == kaonZeroCode) 803 && (pv[1].getCode() == protonCode) ) 804 || ( (pv[0].getCode() == antiKaonZeroCode) 805 && (pv[1].getCode() == protonCode) ) ) 806 { 807 G4double ran = G4UniformRand(); 808 if( pv[1].getCode() == protonCode) 809 { 810 if(ran < 0.68) 811 { 812 pv[0] = PionPlus; 813 pv[1] = Lambda; 814 } 815 else if (ran < 0.84) 816 { 817 pv[0] = PionZero; 818 pv[1] = SigmaPlus; 819 } 820 else 821 { 822 pv[0] = PionPlus; 823 pv[1] = SigmaZero; 824 } 825 } 826 else 827 { 828 if(ran < 0.68) 829 { 830 pv[0] = PionMinus; 831 pv[1] = Lambda; 832 } 833 else if (ran < 0.84) 834 { 835 pv[0] = PionMinus; 836 pv[1] = SigmaZero; 837 } 838 else 839 { 840 pv[0] = PionZero; 841 pv[1] = SigmaMinus; 842 } 843 } 844 } 845 else 846 { 847 G4double ran = G4UniformRand(); 848 if (ran < 0.67) 849 { 850 pv[0] = PionZero; 851 pv[1] = Lambda; 852 } 853 else if (ran < 0.78) 854 { 855 pv[0] = PionMinus; 856 pv[1] = SigmaPlus; 857 } 858 else if (ran < 0.89) 859 { 860 pv[0] = PionZero; 861 pv[1] = SigmaZero; 862 } 863 else 864 { 865 pv[0] = PionPlus; 866 pv[1] = SigmaMinus; 867 } 868 } 869 } 870 871 nt = np + nm + nz; 872 while ( nt > 0) { 873 G4double ran = G4UniformRand(); 874 if ( ran < (G4double)np/nt) { 875 if( np > 0 ) { 876 pv[vecLen++] = PionPlus; 877 np--; 878 } 879 } else if (ran < (G4double)(np+nm)/nt) { 880 if( nm > 0 ) { 881 pv[vecLen++] = PionMinus; 882 nm--; 883 } 884 } else { 885 if( nz > 0 ) { 886 pv[vecLen++] = PionZero; 887 nz--; 888 } 889 } 890 nt = np + nm + nz; 891 } 892 893 if (verboseLevel > 1) { 894 G4cout << "Particles produced: " ; 895 G4cout << pv[0].getName() << " " ; 896 G4cout << pv[1].getName() << " " ; 897 for (i=2; i < vecLen; i++) G4cout << pv[i].getName() << " " ; 898 G4cout << G4endl; 899 } 900 901 return; 902 }
Note: See TracChangeset
for help on using the changeset viewer.